generalised_newtonian_refineable_navier_stokes_elements.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1194 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-27 08:44:43 +0100 (Fri, 27 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
31 
32 
33 
34 namespace oomph
35 {
36 
37 
38 
39 
40 
41 
42 
43 //===================================================================
44 /// Compute the diagonal of the velocity/pressure mass matrices.
45 /// If which one=0, both are computed, otherwise only the pressure
46 /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
47 /// LSC version of the preconditioner only needs that one)
48 //===================================================================
49  template<unsigned DIM>
52  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
53  const unsigned& which_one)
54  {
55 
56  // Resize and initialise
57  unsigned n_dof=ndof();
58 
59  if ((which_one==0)||(which_one==1)) press_mass_diag.assign(n_dof,0.0);
60  if ((which_one==0)||(which_one==2)) veloc_mass_diag.assign(n_dof,0.0);
61 
62  //Pointers to hang info object
63  HangInfo *hang_info_pt=0;
64 
65  //Number of master nodes and weight for shape fcts
66  unsigned n_master=1;
67  double hang_weight=1.0;
68 
69  // find out how many nodes there are
70  unsigned n_node = nnode();
71 
72  //Set up memory for veloc shape functions
73  Shape psi(n_node);
74 
75  //Find number of pressure dofs
76  unsigned n_pres = this->npres_nst();
77 
78  // Pressure shape function
79  Shape psi_p(n_pres);
80 
81  // Local coordinates
82  Vector<double> s(DIM);
83 
84  // find the indices at which the local velocities are stored
85  Vector<unsigned> u_nodal_index(DIM);
86  for(unsigned i=0; i<DIM; i++)
87  {
88  u_nodal_index[i] = this->u_index_nst(i);
89  }
90 
91  // Which nodal value represents the pressure? (Negative if pressure
92  // is not based on nodal interpolation).
93  int p_index = this->p_nodal_index_nst();
94 
95  // Local array of booleans that are true if the l-th pressure value is
96  // hanging (avoid repeated virtual function calls)
97  bool pressure_dof_is_hanging[n_pres];
98 
99  //If the pressure is stored at a node
100  if(p_index >= 0)
101  {
102  //Read out whether the pressure is hanging
103  for(unsigned l=0;l<n_pres;++l)
104  {
105  pressure_dof_is_hanging[l] =
106  pressure_node_pt(l)->is_hanging(p_index);
107  }
108  }
109  //Otherwise the pressure is not stored at a node and so cannot hang
110  else
111  {
112  for(unsigned l=0;l<n_pres;++l)
113  {pressure_dof_is_hanging[l] = false;}
114  }
115 
116 
117  //Number of integration points
118  unsigned n_intpt = integral_pt()->nweight();
119 
120  //Integer to store the local equations no
121  int local_eqn=0;
122 
123  //Loop over the integration points
124  for(unsigned ipt=0; ipt<n_intpt; ipt++)
125  {
126 
127  //Get the integral weight
128  double w = integral_pt()->weight(ipt);
129 
130  //Get determinant of Jacobian of the mapping
131  double J = J_eulerian_at_knot(ipt);
132 
133  //Assign values of s
134  for(unsigned i=0;i<DIM;i++)
135  {
136  s[i] = integral_pt()->knot(ipt,i);
137  }
138 
139  //Premultiply weights and Jacobian
140  double W = w*J;
141 
142 
143 
144  // Do we want the velocity one?
145  if ((which_one==0)||(which_one==2))
146  {
147 
148  //Get the velocity shape functions
149  shape_at_knot(ipt,psi);
150 
151 
152  //Number of master nodes and storage for the weight of the shape function
153  unsigned n_master=1; double hang_weight=1.0;
154 
155  //Loop over the nodes for the test functions/equations
156  //----------------------------------------------------
157  for(unsigned l=0;l<n_node;l++)
158  {
159  //Local boolean to indicate whether the node is hanging
160  bool is_node_hanging = node_pt(l)->is_hanging();
161 
162  //If the node is hanging
163  if(is_node_hanging)
164  {
165  hang_info_pt = node_pt(l)->hanging_pt();
166 
167  //Read out number of master nodes from hanging data
168  n_master = hang_info_pt->nmaster();
169  }
170  //Otherwise the node is its own master
171  else
172  {
173  n_master = 1;
174  }
175 
176  //Loop over the master nodes
177  for(unsigned m=0;m<n_master;m++)
178  {
179  // Loop over velocity components for equations
180  for(unsigned i=0;i<DIM;i++)
181  {
182  //Get the equation number
183  //If the node is hanging
184  if(is_node_hanging)
185  {
186  //Get the equation number from the master node
187  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
188  u_nodal_index[i]);
189  //Get the hang weight from the master node
190  hang_weight = hang_info_pt->master_weight(m);
191  }
192  //If the node is not hanging
193  else
194  {
195  // Local equation number
196  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
197 
198  // Node contributes with full weight
199  hang_weight = 1.0;
200  }
201 
202  //If it's not a boundary condition...
203  if(local_eqn>= 0)
204  {
205 
206 // //Loop over the veclocity shape functions
207 // for(unsigned l=0; l<n_node; l++)
208 // {
209 // //Loop over the velocity components
210 // for(unsigned i=0; i<DIM; i++)
211 // {
212 // local_eqn = nodal_local_eqn(l,u_nodal_index[i]);
213 
214 // //If not a boundary condition
215 // if(local_eqn >= 0)
216 // {
217 
218 
219 
220  //Add the contribution
221  veloc_mass_diag[local_eqn] += pow(psi[l]*hang_weight,2)*W;
222  }
223  }
224  }
225  }
226  }
227 
228  // Do we want the pressure one?
229  if ((which_one==0)||(which_one==1))
230  {
231  //Get the pressure shape functions
232  this->pshape_nst(s,psi_p);
233 
234  //Loop over the pressure shape functions
235  for(unsigned l=0;l<n_pres;l++)
236  {
237  //If the pressure dof is hanging
238  if(pressure_dof_is_hanging[l])
239  {
240  // Pressure dof is hanging so it must be nodal-based
241  // Get the hang info object
242  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
243 
244  //Get the number of master nodes from the pressure node
245  n_master = hang_info_pt->nmaster();
246  }
247  //Otherwise the node is its own master
248  else
249  {
250  n_master = 1;
251  }
252 
253  //Loop over the master nodes
254  for(unsigned m=0;m<n_master;m++)
255  {
256  //Get the number of the unknown
257  //If the pressure dof is hanging
258  if(pressure_dof_is_hanging[l])
259  {
260  //Get the local equation from the master node
261  local_eqn =
262  this->local_hang_eqn(hang_info_pt->master_node_pt(m),p_index);
263  //Get the weight for the node
264  hang_weight = hang_info_pt->master_weight(m);
265  }
266  else
267  {
268  local_eqn = this->p_local_eqn(l);
269  hang_weight = 1.0;
270  }
271 
272  //If the equation is not pinned
273  if(local_eqn >= 0)
274  {
275 
276 // //Loop over the veclocity shape functions
277 // for(unsigned l=0; l<n_pres; l++)
278 // {
279 // // Get equation number
280 // local_eqn = p_local_eqn(l);
281 
282 // //If not a boundary condition
283 // if(local_eqn >= 0)
284 // {
285 
286 
287  //Add the contribution
288  press_mass_diag[local_eqn] += pow(psi_p[l]*hang_weight,2) * W;
289  }
290  }
291  }
292 
293  }
294  }
295  }
296 
297 
298 
299 //========================================================================
300 /// Add element's contribution to the elemental
301 /// residual vector and/or Jacobian matrix.
302 /// flag=1: compute both
303 /// flag=0: compute only residual vector
304 //========================================================================
305 template<unsigned DIM>
308  DenseMatrix<double> &jacobian,
309  DenseMatrix<double> &mass_matrix,
310  unsigned flag)
311 {
312 //Find out how many nodes there are
313 unsigned n_node = nnode();
314 
315 // Get continuous time from timestepper of first node
316 double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
317 
318 //Find out how many pressure dofs there are
319 unsigned n_pres = this->npres_nst();
320 
321 // Get the indices at which the velocity components are stored
322  unsigned u_nodal_index[DIM];
323  for(unsigned i=0;i<DIM;i++) {u_nodal_index[i] = this->u_index_nst(i);}
324 
325 // Which nodal value represents the pressure? (Negative if pressure
326 // is not based on nodal interpolation).
327 int p_index = this->p_nodal_index_nst();
328 
329 // Local array of booleans that are true if the l-th pressure value is
330 // hanging (avoid repeated virtual function calls)
331  bool pressure_dof_is_hanging[n_pres];
332  //If the pressure is stored at a node
333  if(p_index >= 0)
334  {
335  //Read out whether the pressure is hanging
336  for(unsigned l=0;l<n_pres;++l)
337  {
338  pressure_dof_is_hanging[l] =
339  pressure_node_pt(l)->is_hanging(p_index);
340  }
341  }
342  //Otherwise the pressure is not stored at a node and so cannot hang
343  else
344  {
345  for(unsigned l=0;l<n_pres;++l)
346  {pressure_dof_is_hanging[l] = false;}
347  }
348 
349 //Set up memory for the shape and test functions
350 Shape psif(n_node), testf(n_node);
351 DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
352 
353 
354 //Set up memory for pressure shape and test functions
355 Shape psip(n_pres), testp(n_pres);
356 
357 //Set the value of n_intpt
358 unsigned n_intpt = integral_pt()->nweight();
359 
360 //Set the Vector to hold local coordinates
361 Vector<double> s(DIM);
362 
363 //Get Physical Variables from Element
364 //Reynolds number must be multiplied by the density ratio
365 double scaled_re = this->re()*this->density_ratio();
366 double scaled_re_st = this->re_st()*this->density_ratio();
367 double scaled_re_inv_fr = this->re_invfr()*this->density_ratio();
368 double visc_ratio = this->viscosity_ratio();
369 Vector<double> G = this->g();
370 
371 //Integers that store the local equations and unknowns
372 int local_eqn=0, local_unknown=0;
373 
374 //Pointers to hang info objects
375  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
376 
377 //Local boolean for ALE (or not)
378  bool ALE_is_disabled_flag = this->ALE_is_disabled;
379 
380 //Loop over the integration points
381 for(unsigned ipt=0;ipt<n_intpt;ipt++)
382 {
383 
384  //Assign values of s
385  for(unsigned i=0;i<DIM;i++) {s[i] = integral_pt()->knot(ipt,i);}
386 
387  //Get the integral weight
388  double w = integral_pt()->weight(ipt);
389 
390  //Call the derivatives of the shape and test functions
391  double J =
392  this->dshape_and_dtest_eulerian_at_knot_nst(ipt,psif,dpsifdx,testf,dtestfdx);
393 
394  //Call the pressure shape and test functions
395  this->pshape_nst(s,psip,testp);
396 
397  //Premultiply the weights and the Jacobian
398  double W = w*J;
399 
400  //Calculate local values of the pressure and velocity components
401  //--------------------------------------------------------------
402  double interpolated_p=0.0;
403  Vector<double> interpolated_u(DIM,0.0);
404  Vector<double> interpolated_x(DIM,0.0);
405  Vector<double> mesh_veloc(DIM,0.0);
406  Vector<double> dudt(DIM,0.0);
407  DenseMatrix<double> interpolated_dudx(DIM,DIM,0.0);
408 
409  //Calculate pressure
410  for(unsigned l=0;l<n_pres;l++) {interpolated_p += this->p_nst(l)*psip[l];}
411 
412 
413  //Calculate velocities and derivatives
414 
415  // Loop over nodes
416  for(unsigned l=0;l<n_node;l++)
417  {
418  //Loop over directions
419  for(unsigned i=0;i<DIM;i++)
420  {
421  //Get the nodal value
422  double u_value = this->nodal_value(l,u_nodal_index[i]);
423  interpolated_u[i] += u_value*psif[l];
424  interpolated_x[i] += this->nodal_position(l,i)*psif[l];
425  dudt[i] += this->du_dt_nst(l,i)*psif[l];
426 
427  //Loop over derivative directions for velocity gradients
428  for(unsigned j=0;j<DIM;j++)
429  {
430  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
431  }
432  }
433  }
434 
435  if (!ALE_is_disabled_flag)
436  {
437  // Loop over nodes
438  for(unsigned l=0;l<n_node;l++)
439  {
440  //Loop over directions
441  for(unsigned i=0;i<DIM;i++)
442  {
443  mesh_veloc[i] += this->dnodal_position_dt(l,i)*psif[l];
444  }
445  }
446  }
447 
448  // The strainrate used to compute the second invariant
449  DenseMatrix<double> strainrate_to_compute_second_invariant(DIM,DIM,0.0);
450 
451  // the strainrate used to calculate the second invariant
452  // can be either the current one or the one extrapolated from
453  // previous velocity values
454  if(!this->Use_extrapolated_strainrate_to_compute_second_invariant)
455  {
456  this->strain_rate(s,strainrate_to_compute_second_invariant);
457  }
458  else
459  {
460  this->extrapolated_strain_rate(ipt,strainrate_to_compute_second_invariant);
461  }
462 
463  // Calculate the second invariant
465  strainrate_to_compute_second_invariant);
466 
467  // Get the viscosity according to the constitutive equation
468  double viscosity=this->Constitutive_eqn_pt->viscosity(second_invariant);
469 
470  //Get the user-defined body force terms
471  Vector<double> body_force(DIM);
472  this->get_body_force_nst(time,ipt,s,interpolated_x,body_force);
473 
474  //Get the user-defined source function
475  double source = this->get_source_nst(time,ipt,interpolated_x);
476 
477  // The generalised Newtonian viscosity differentiated w.r.t.
478  // the unknown velocities
479  DenseMatrix<double> dviscosity_dunknown(DIM,n_node,0.0);
480 
481  if(!this->Use_extrapolated_strainrate_to_compute_second_invariant)
482  {
483  // Calculate the derivate of the viscosity w.r.t. the second invariant
484  double dviscosity_dsecond_invariant=
485  this->Constitutive_eqn_pt->dviscosity_dinvariant(second_invariant);
486 
487  // calculate the strainrate
488  DenseMatrix<double> strainrate(DIM,DIM,0.0);
489  this->strain_rate(s,strainrate);
490 
491  // pre-compute the derivative of the second invariant w.r.t. the
492  // entries in the rate of strain tensor
493  DenseMatrix<double> dinvariant_dstrainrate(DIM,DIM,0.0);
494 
495  // setting up a Kronecker Delta matrix, which has ones at the diagonal
496  // and zeros off-diagonally
497  DenseMatrix<double> kroneckerdelta(DIM,DIM,0.0);
498 
499  for(unsigned i=0;i<DIM;i++)
500  {
501  for(unsigned j=0;j<DIM;j++)
502  {
503  if(i==j)
504  {
505  // Set the diagonal entries of the Kronecker delta
506  kroneckerdelta(i,i)=1.0;
507 
508  double tmp=0.0;
509  for(unsigned k=0;k<DIM;k++)
510  {
511  if(k!=i)
512  {
513  tmp+=strainrate(k,k);
514  }
515  }
516  dinvariant_dstrainrate(i,i)=tmp;
517  }
518  else
519  {
520  dinvariant_dstrainrate(i,j)=-strainrate(j,i);
521  }
522  }
523  }
524 
525  // a rank four tensor storing the derivative of the strainrate
526  // w.r.t. the unknowns
527  RankFourTensor<double> dstrainrate_dunknown(DIM,DIM,DIM,n_node);
528 
529  // loop over the nodes
530  for(unsigned l=0;l<n_node;l++)
531  {
532  // loop over the velocity components
533  for(unsigned k=0;k<DIM;k++)
534  {
535  // loop over the entries of the rate of strain tensor
536  for(unsigned i=0;i<DIM;i++)
537  {
538  for(unsigned j=0;j<DIM;j++)
539  {
540  // initialise the entry to zero
541  dstrainrate_dunknown(i,j,k,l)=0.0;
542 
543  // loop over velocities and directions
544  for(unsigned m=0;m<DIM;m++)
545  {
546  for(unsigned n=0;n<DIM;n++)
547  {
548  dstrainrate_dunknown(i,j,k,l)+=
549  0.5*(kroneckerdelta(i,m)*kroneckerdelta(j,n)+
550  kroneckerdelta(i,n)*kroneckerdelta(j,m))*
551  kroneckerdelta(m,k)*dpsifdx(l,n);
552  }
553  }
554  }
555  }
556  }
557  }
558 
559  // Now calculate the derivative of the viscosity w.r.t. the unknowns
560  // loop over the nodes
561  for(unsigned l=0;l<n_node;l++)
562  {
563  // loop over the velocities
564  for(unsigned k=0;k<DIM;k++)
565  {
566  // loop over the entries in the rate of strain tensor
567  for(unsigned i=0;i<DIM;i++)
568  {
569  for(unsigned j=0;j<DIM;j++)
570  {
571  dviscosity_dunknown(k,l)+=dviscosity_dsecond_invariant*
572  dinvariant_dstrainrate(i,j)*dstrainrate_dunknown(i,j,k,l);
573  }
574  }
575  }
576  }
577 
578  }
579 
580 
581  //MOMENTUM EQUATIONS
582  //==================
583 
584  //Number of master nodes and storage for the weight of the shape function
585  unsigned n_master=1; double hang_weight=1.0;
586 
587  //Loop over the nodes for the test functions/equations
588  //----------------------------------------------------
589  for(unsigned l=0;l<n_node;l++)
590  {
591  //Local boolean to indicate whether the node is hanging
592  bool is_node_hanging = node_pt(l)->is_hanging();
593 
594  //If the node is hanging
595  if(is_node_hanging)
596  {
597  hang_info_pt = node_pt(l)->hanging_pt();
598 
599  //Read out number of master nodes from hanging data
600  n_master = hang_info_pt->nmaster();
601  }
602  //Otherwise the node is its own master
603  else
604  {
605  n_master = 1;
606  }
607 
608  //Loop over the master nodes
609  for(unsigned m=0;m<n_master;m++)
610  {
611  // Loop over velocity components for equations
612  for(unsigned i=0;i<DIM;i++)
613  {
614  //Get the equation number
615  //If the node is hanging
616  if(is_node_hanging)
617  {
618  //Get the equation number from the master node
619  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
620  u_nodal_index[i]);
621  //Get the hang weight from the master node
622  hang_weight = hang_info_pt->master_weight(m);
623  }
624  //If the node is not hanging
625  else
626  {
627  // Local equation number
628  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
629 
630  // Node contributes with full weight
631  hang_weight = 1.0;
632  }
633 
634  //If it's not a boundary condition...
635  if(local_eqn>= 0)
636  {
637  //Temporary variable to hold the residuals
638  double sum=0.0;
639 
640  //Add the user-defined body force terms
641  sum += body_force[i];
642 
643  //Add the gravitational body force term
644  sum += scaled_re_inv_fr*G[i];
645 
646  //Add in the inertial term
647  sum -= scaled_re_st*dudt[i];
648 
649  //Convective terms, including mesh velocity
650  for(unsigned k=0;k<DIM;k++)
651  {
652  double tmp=scaled_re*interpolated_u[k];
653  if (!ALE_is_disabled_flag)
654  {tmp -= scaled_re_st*mesh_veloc[k];}
655  sum -= tmp*interpolated_dudx(i,k);
656  }
657 
658  //Add the pressure gradient term
659  // Potentially pre-multiply by viscosity ratio
660  if(this->Pre_multiply_by_viscosity_ratio)
661  {
662  sum = visc_ratio*viscosity*
663  (sum*testf[l] + interpolated_p*dtestfdx(l,i))*W*hang_weight;
664  }
665  else
666  {
667  sum = (sum*testf[l] + interpolated_p*dtestfdx(l,i))*W*hang_weight;
668  }
669 
670  //Add in the stress tensor terms
671  //The viscosity ratio needs to go in here to ensure
672  //continuity of normal stress is satisfied even in flows
673  //with zero pressure gradient!
674  for(unsigned k=0;k<DIM;k++)
675  {
676  sum -= visc_ratio*viscosity*
677  (interpolated_dudx(i,k) + this->Gamma[i]*interpolated_dudx(k,i))
678  *dtestfdx(l,k)*W*hang_weight;
679  }
680 
681  // Add contribution
682  residuals[local_eqn] += sum;
683 
684  //CALCULATE THE JACOBIAN
685  if(flag)
686  {
687  //Number of master nodes and weights
688  unsigned n_master2=1; double hang_weight2=1.0;
689  //Loop over the velocity nodes for columns
690  for(unsigned l2=0;l2<n_node;l2++)
691  {
692  //Local boolean to indicate whether the node is hanging
693  bool is_node2_hanging = node_pt(l2)->is_hanging();
694 
695  //If the node is hanging
696  if(is_node2_hanging)
697  {
698  hang_info2_pt = node_pt(l2)->hanging_pt();
699  //Read out number of master nodes from hanging data
700  n_master2 = hang_info2_pt->nmaster();
701  }
702  //Otherwise the node is its own master
703  else
704  {
705  n_master2 = 1;
706  }
707 
708  //Loop over the master nodes
709  for(unsigned m2=0;m2<n_master2;m2++)
710  {
711  //Loop over the velocity components
712  for(unsigned i2=0;i2<DIM;i2++)
713  {
714  //Get the number of the unknown
715  //If the node is hanging
716  if(is_node2_hanging)
717  {
718  //Get the equation number from the master node
719  local_unknown =
720  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
721  u_nodal_index[i2]);
722  //Get the hang weights
723  hang_weight2 = hang_info2_pt->master_weight(m2);
724  }
725  else
726  {
727  local_unknown = this->nodal_local_eqn(l2,u_nodal_index[i2]);
728  hang_weight2 = 1.0;
729  }
730 
731  // If the unknown is non-zero
732  if(local_unknown >= 0)
733  {
734  //Add contribution to Elemental Matrix
735  jacobian(local_eqn,local_unknown)
736  -= visc_ratio*viscosity*this->Gamma[i]*dpsifdx(l2,i)*
737  dtestfdx(l,i2)*W*hang_weight*hang_weight2;
738 
739  //Now add in the inertial terms
740  jacobian(local_eqn,local_unknown)
741  -= scaled_re*psif[l2]*interpolated_dudx(i,i2)*testf[l]*W*
742  hang_weight*hang_weight2;
743 
744  if(!this->
745  Use_extrapolated_strainrate_to_compute_second_invariant)
746  {
747  for(unsigned k=0;k<DIM;k++)
748  {
749  jacobian(local_eqn,local_unknown) -=
750  visc_ratio*dviscosity_dunknown(i2,l2)*
751  (interpolated_dudx(i,k) +
752  this->Gamma[i]*interpolated_dudx(k,i))
753  *dtestfdx(l,k)*W*hang_weight*hang_weight2;
754  }
755  }
756 
757  //Extra diagonal components if i2=i
758  if(i2 == i)
759  {
760  //Mass matrix entries
761  //Again note the positive sign because the mass
762  //matrix is taken on the other side of the equation
763  if(flag==2)
764  {
765  mass_matrix(local_eqn,local_unknown) +=
766  scaled_re_st*psif[l2]*testf[l]*W*
767  hang_weight*hang_weight2;
768  }
769 
770  // du/dt term
771  jacobian(local_eqn,local_unknown)
772  -= scaled_re_st*
773  node_pt(l2)->time_stepper_pt()->weight(1,0)*
774  psif[l2]*testf[l]*W*hang_weight*hang_weight2;
775 
776  //Extra advective terms
777  for(unsigned k=0;k<DIM;k++)
778  {
779  double tmp=scaled_re*interpolated_u[k];
780  if (!ALE_is_disabled_flag)
781  {tmp -= scaled_re_st*mesh_veloc[k];}
782 
783  jacobian(local_eqn,local_unknown) -=
784  tmp*dpsifdx(l2,k)*testf[l]*W*hang_weight*hang_weight2;
785  }
786 
787  // Extra viscous terms
788  for(unsigned k=0;k<DIM;k++)
789  {
790  jacobian(local_eqn,local_unknown)
791  -= visc_ratio*viscosity*dpsifdx(l2,k)*
792  dtestfdx(l,k)*W*hang_weight*hang_weight2;
793  }
794 
795  }
796  }
797  }
798  }
799  }
800 
801  //Loop over the pressure shape functions
802  for(unsigned l2=0;l2<n_pres;l2++)
803  {
804  //If the pressure dof is hanging
805  if(pressure_dof_is_hanging[l2])
806  {
807  hang_info2_pt =
808  this->pressure_node_pt(l2)->hanging_pt(p_index);
809  // Pressure dof is hanging so it must be nodal-based
810  //Get the number of master nodes from the pressure node
811  n_master2 = hang_info2_pt->nmaster();
812  }
813  //Otherwise the node is its own master
814  else
815  {
816  n_master2 = 1;
817  }
818 
819  //Loop over the master nodes
820  for(unsigned m2=0;m2<n_master2;m2++)
821  {
822  //Get the number of the unknown
823  //If the pressure dof is hanging
824  if(pressure_dof_is_hanging[l2])
825  {
826  //Get the unknown from the master node
827  local_unknown =
828  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
829  p_index);
830  //Get the weight from the hanging object
831  hang_weight2 = hang_info2_pt->master_weight(m2);
832  }
833  else
834  {
835  local_unknown = this->p_local_eqn(l2);
836  hang_weight2 = 1.0;
837  }
838 
839  //If the unknown is not pinned
840  if(local_unknown >= 0)
841  {
842  if(this->Pre_multiply_by_viscosity_ratio)
843  {
844  jacobian(local_eqn,local_unknown)
845  += visc_ratio*viscosity*
846  psip[l2]*dtestfdx(l,i)*W*hang_weight*hang_weight2;
847  }
848  else
849  {
850  jacobian(local_eqn,local_unknown)
851  += psip[l2]*dtestfdx(l,i)*W*hang_weight*hang_weight2;
852  }
853  }
854  }
855  }
856 
857  }// End of Jacobian calculation
858 
859  } //End of if not boundary condition statement
860 
861  } //End of loop over components of non-hanging node
862 
863  } //End of loop over master nodes
864 
865  } // End of loop over nodes for equations
866 
867 
868 
869  //CONTINUITY EQUATION
870  //===================
871 
872  //Loop over the pressure shape functions
873  for(unsigned l=0;l<n_pres;l++)
874  {
875  //If the pressure dof is hanging
876  if(pressure_dof_is_hanging[l])
877  {
878  // Pressure dof is hanging so it must be nodal-based
879  // Get the hang info object
880  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
881 
882  //Get the number of master nodes from the pressure node
883  n_master = hang_info_pt->nmaster();
884  }
885  //Otherwise the node is its own master
886  else
887  {
888  n_master = 1;
889  }
890 
891  //Loop over the master nodes
892  for(unsigned m=0;m<n_master;m++)
893  {
894  //Get the number of the unknown
895  //If the pressure dof is hanging
896  if(pressure_dof_is_hanging[l])
897  {
898  //Get the local equation from the master node
899  local_eqn =
900  this->local_hang_eqn(hang_info_pt->master_node_pt(m),p_index);
901  //Get the weight for the node
902  hang_weight = hang_info_pt->master_weight(m);
903  }
904  else
905  {
906  local_eqn = this->p_local_eqn(l);
907  hang_weight = 1.0;
908  }
909 
910  //If the equation is not pinned
911  if(local_eqn >= 0)
912  {
913  // Source term
914  residuals[local_eqn] -= source*testp[l]*W*hang_weight;
915 
916  // Loop over velocity components
917  for(unsigned k=0;k<DIM;k++)
918  {
919  // Potentially pre-multiply by viscosity ratio
920  if(this->Pre_multiply_by_viscosity_ratio)
921  {
922  residuals[local_eqn] += visc_ratio*viscosity*
923  interpolated_dudx(k,k)*testp[l]*W*hang_weight;
924  }
925  else
926  {
927  residuals[local_eqn] +=
928  interpolated_dudx(k,k)*testp[l]*W*hang_weight;
929  }
930  }
931 
932  //CALCULATE THE JACOBIAN
933  if(flag)
934  {
935  unsigned n_master2=1; double hang_weight2=1.0;
936  //Loop over the velocity nodes for columns
937  for(unsigned l2=0;l2<n_node;l2++)
938  {
939  //Local boolean to indicate whether the node is hanging
940  bool is_node2_hanging = node_pt(l2)->is_hanging();
941 
942  //If the node is hanging
943  if(is_node2_hanging)
944  {
945  hang_info2_pt = node_pt(l2)->hanging_pt();
946  //Read out number of master nodes from hanging data
947  n_master2 = hang_info2_pt->nmaster();
948  }
949  //Otherwise the node is its own master
950  else
951  {
952  n_master2 = 1;
953  }
954 
955  //Loop over the master nodes
956  for(unsigned m2=0;m2<n_master2;m2++)
957  {
958  //Loop over the velocity components
959  for(unsigned i2=0;i2<DIM;i2++)
960  {
961  //Get the number of the unknown
962  //If the node is hanging
963  if(is_node2_hanging)
964  {
965  //Get the equation number from the master node
966  local_unknown =
967  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
968  u_nodal_index[i2]);
969  hang_weight2 = hang_info2_pt->master_weight(m2);
970  }
971  else
972  {
973  local_unknown = this->nodal_local_eqn(l2,u_nodal_index[i2]);
974  hang_weight2 = 1.0;
975  }
976 
977  //If the unknown is not pinned
978  if(local_unknown >= 0)
979  {
980  if(this->Pre_multiply_by_viscosity_ratio)
981  {
982  jacobian(local_eqn,local_unknown)+= visc_ratio*viscosity*
983  dpsifdx(l2,i2)*testp[l]*W*hang_weight*hang_weight2;
984  }
985  else
986  {
987  jacobian(local_eqn,local_unknown)
988  += dpsifdx(l2,i2)*testp[l]*W*hang_weight*hang_weight2;
989  }
990  }
991  }
992  }
993  }
994 
995  // NO PRESSURE CONTRIBUTION TO THE JACOBIAN
996 
997  } //End of jacobian calculation
998  }
999  }
1000  } //End of loop over pressure variables
1001 
1002 } //End of loop over integration points
1003 }
1004 
1005 
1006 
1007 
1008 //======================================================================
1009 /// Compute derivatives of elemental residual vector with respect
1010 /// to nodal coordinates.
1011 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
1012 /// Overloads the FD-based version in the FE base class.
1013 //======================================================================
1014 template <unsigned DIM>
1017  dresidual_dnodal_coordinates)
1018 {
1019  // Return immediately if there are no dofs
1020  if(ndof()==0) { return; }
1021 
1022  // Determine number of nodes in element
1023  const unsigned n_node = nnode();
1024 
1025  // Get continuous time from timestepper of first node
1026  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
1027 
1028  // Determine number of pressure dofs in element
1029  const unsigned n_pres = this->npres_nst();
1030 
1031  // Find the indices at which the local velocities are stored
1032  unsigned u_nodal_index[DIM];
1033  for(unsigned i=0;i<DIM;i++) { u_nodal_index[i] = this->u_index_nst(i); }
1034 
1035  // Which nodal value represents the pressure? (Negative if pressure
1036  // is not based on nodal interpolation).
1037  const int p_index = this->p_nodal_index_nst();
1038 
1039  // Local array of booleans that are true if the l-th pressure value is
1040  // hanging (avoid repeated virtual function calls)
1041  bool pressure_dof_is_hanging[n_pres];
1042 
1043  // If the pressure is stored at a node
1044  if(p_index >= 0)
1045  {
1046  // Read out whether the pressure is hanging
1047  for(unsigned l=0;l<n_pres;++l)
1048  {
1049  pressure_dof_is_hanging[l] =
1050  pressure_node_pt(l)->is_hanging(p_index);
1051  }
1052  }
1053  // Otherwise the pressure is not stored at a node and so cannot hang
1054  else
1055  {
1056  for(unsigned l=0;l<n_pres;++l) { pressure_dof_is_hanging[l] = false; }
1057  }
1058 
1059  // Set up memory for the shape and test functions
1060  Shape psif(n_node), testf(n_node);
1061  DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
1062 
1063  // Set up memory for pressure shape and test functions
1064  Shape psip(n_pres), testp(n_pres);
1065 
1066  // Determine number of shape controlling nodes
1067  const unsigned n_shape_controlling_node = nshape_controlling_nodes();
1068 
1069  // Deriatives of shape fct derivatives w.r.t. nodal coords
1070  RankFourTensor<double> d_dpsifdx_dX(DIM,n_shape_controlling_node,n_node,DIM);
1071  RankFourTensor<double> d_dtestfdx_dX(DIM,n_shape_controlling_node,n_node,DIM);
1072 
1073  // Derivative of Jacobian of mapping w.r.t. to nodal coords
1074  DenseMatrix<double> dJ_dX(DIM,n_shape_controlling_node);
1075 
1076  // Derivatives of derivative of u w.r.t. nodal coords
1077  RankFourTensor<double> d_dudx_dX(DIM,n_shape_controlling_node,DIM,DIM);
1078 
1079  // Derivatives of nodal velocities w.r.t. nodal coords:
1080  // Assumption: Interaction only local via no-slip so
1081  // X_ij only affects U_ij.
1082  DenseMatrix<double> d_U_dX(DIM,n_shape_controlling_node,0.0);
1083 
1084  // Determine the number of integration points
1085  const unsigned n_intpt = integral_pt()->nweight();
1086 
1087  // Vector to hold local coordinates
1088  Vector<double> s(DIM);
1089 
1090  // Get physical variables from element
1091  // (Reynolds number must be multiplied by the density ratio)
1092  double scaled_re = this->re()*this->density_ratio();
1093  double scaled_re_st = this->re_st()*this->density_ratio();
1094  double scaled_re_inv_fr = this->re_invfr()*this->density_ratio();
1095  double visc_ratio = this->viscosity_ratio();
1096  Vector<double> G = this->g();
1097 
1098  // FD step
1100 
1101  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1102  // Assumption: Interaction only local via no-slip so
1103  // X_ij only affects U_ij.
1104  bool element_has_node_with_aux_node_update_fct=false;
1105 
1106  std::map<Node*,unsigned> local_shape_controlling_node_lookup=
1107  shape_controlling_node_lookup();
1108 
1109  // FD loop over shape-controlling nodes
1110  for(std::map<Node*,unsigned>::iterator it=
1111  local_shape_controlling_node_lookup.begin();
1112  it!=local_shape_controlling_node_lookup.end();
1113  it++)
1114  {
1115  // Get node
1116  Node* nod_pt=it->first;
1117 
1118  // Get its number
1119  unsigned q=it->second;
1120 
1121  // Only compute if there's a node-update fct involved
1122  if(nod_pt->has_auxiliary_node_update_fct_pt())
1123  {
1124  element_has_node_with_aux_node_update_fct=true;
1125 
1126  // Current nodal velocity
1127  Vector<double> u_ref(DIM);
1128  for(unsigned i=0;i<DIM;i++)
1129  {
1130  u_ref[i]=*(nod_pt->value_pt(u_nodal_index[i]));
1131  }
1132 
1133  // FD
1134  for(unsigned p=0;p<DIM;p++)
1135  {
1136  // Make backup
1137  double backup=nod_pt->x(p);
1138 
1139  // Do FD step. No node update required as we're
1140  // attacking the coordinate directly...
1141  nod_pt->x(p)+=eps_fd;
1142 
1143  // Do auxiliary node update (to apply no slip)
1145 
1146  // Evaluate
1147  d_U_dX(p,q)=(*(nod_pt->value_pt(u_nodal_index[p]))-u_ref[p])/eps_fd;
1148 
1149  // Reset
1150  nod_pt->x(p)=backup;
1151 
1152  // Do auxiliary node update (to apply no slip)
1154  }
1155  }
1156  }
1157 
1158  // Integer to store the local equation number
1159  int local_eqn=0;
1160 
1161  // Pointers to hang info object
1162  HangInfo *hang_info_pt=0;
1163 
1164  // Loop over the integration points
1165  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1166  {
1167  // Assign values of s
1168  for(unsigned i=0;i<DIM;i++) { s[i] = integral_pt()->knot(ipt,i); }
1169 
1170  // Get the integral weight
1171  const double w = integral_pt()->weight(ipt);
1172 
1173  // Call the derivatives of the shape and test functions
1174  const double J = this->dshape_and_dtest_eulerian_at_knot_nst(
1175  ipt,psif,dpsifdx,d_dpsifdx_dX,testf,dtestfdx,d_dtestfdx_dX,dJ_dX);
1176 
1177  // Call the pressure shape and test functions
1178  this->pshape_nst(s,psip,testp);
1179 
1180  // Calculate local values of the pressure and velocity components
1181  // Allocate
1182  double interpolated_p=0.0;
1183  Vector<double> interpolated_u(DIM,0.0);
1184  Vector<double> interpolated_x(DIM,0.0);
1185  Vector<double> mesh_velocity(DIM,0.0);
1186  Vector<double> dudt(DIM,0.0);
1187  DenseMatrix<double> interpolated_dudx(DIM,DIM,0.0);
1188 
1189  // Calculate pressure
1190  for(unsigned l=0;l<n_pres;l++) { interpolated_p += this->p_nst(l)*psip[l]; }
1191 
1192  // Calculate velocities and derivatives:
1193 
1194  // Loop over nodes
1195  for(unsigned l=0;l<n_node;l++)
1196  {
1197  // Loop over directions
1198  for(unsigned i=0;i<DIM;i++)
1199  {
1200  // Get the nodal value
1201  const double u_value = nodal_value(l,u_nodal_index[i]);
1202  interpolated_u[i] += u_value*psif[l];
1203  interpolated_x[i] += nodal_position(l,i)*psif[l];
1204  dudt[i] += this->du_dt_nst(l,i)*psif[l];
1205 
1206  // Loop over derivative directions
1207  for(unsigned j=0;j<DIM;j++)
1208  {
1209  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
1210  }
1211  }
1212  }
1213 
1214  if(!this->ALE_is_disabled)
1215  {
1216  // Loop over nodes
1217  for(unsigned l=0;l<n_node;l++)
1218  {
1219  // Loop over directions
1220  for(unsigned i=0;i<DIM;i++)
1221  {
1222  mesh_velocity[i] += this->dnodal_position_dt(l,i)*psif[l];
1223  }
1224  }
1225  }
1226 
1227  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1228 
1229  // Loop over shape-controlling nodes
1230  for(unsigned q=0;q<n_shape_controlling_node;q++)
1231  {
1232  // Loop over coordinate directions
1233  for(unsigned p=0;p<DIM;p++)
1234  {
1235  for(unsigned i=0;i<DIM;i++)
1236  {
1237  for(unsigned k=0;k<DIM;k++)
1238  {
1239  double aux=0.0;
1240  for(unsigned j=0;j<n_node;j++)
1241  {
1242  aux += nodal_value(j,u_nodal_index[i])*d_dpsifdx_dX(p,q,j,k);
1243  }
1244  d_dudx_dX(p,q,i,k) = aux;
1245  }
1246  }
1247  }
1248  }
1249 
1250  // Get weight of actual nodal position/value in computation of mesh
1251  // velocity from positional/value time stepper
1252  const double pos_time_weight
1253  = node_pt(0)->position_time_stepper_pt()->weight(1,0);
1254  const double val_time_weight = node_pt(0)->time_stepper_pt()->weight(1,0);
1255 
1256  // Get the user-defined body force terms
1257  Vector<double> body_force(DIM);
1258  this->get_body_force_nst(time,ipt,s,interpolated_x,body_force);
1259 
1260  // Get the user-defined source function
1261  const double source = this->get_source_nst(time,ipt,interpolated_x);
1262 
1263  // Get gradient of body force function
1264  DenseMatrix<double> d_body_force_dx(DIM,DIM,0.0);
1265  this->get_body_force_gradient_nst(time,ipt,s,
1266  interpolated_x, d_body_force_dx);
1267 
1268  // Get gradient of source function
1269  Vector<double> source_gradient(DIM,0.0);
1270  this->get_source_gradient_nst(time,ipt,interpolated_x, source_gradient);
1271 
1272 
1273  // Assemble shape derivatives
1274  //---------------------------
1275 
1276  // MOMENTUM EQUATIONS
1277  // ------------------
1278 
1279  // Number of master nodes and storage for the weight of the shape function
1280  unsigned n_master=1; double hang_weight=1.0;
1281 
1282  // Loop over the test functions
1283  for(unsigned l=0;l<n_node;l++)
1284  {
1285 
1286  // Local boolean to indicate whether the node is hanging
1287  bool is_node_hanging = node_pt(l)->is_hanging();
1288 
1289  // If the node is hanging
1290  if(is_node_hanging)
1291  {
1292  hang_info_pt = node_pt(l)->hanging_pt();
1293 
1294  // Read out number of master nodes from hanging data
1295  n_master = hang_info_pt->nmaster();
1296  }
1297  // Otherwise the node is its own master
1298  else
1299  {
1300  n_master = 1;
1301  }
1302 
1303  // Loop over the master nodes
1304  for(unsigned m=0;m<n_master;m++)
1305  {
1306 
1307  // Loop over coordinate directions
1308  for(unsigned i=0;i<DIM;i++)
1309  {
1310 
1311  // Get the equation number
1312  // If the node is hanging
1313  if(is_node_hanging)
1314  {
1315  // Get the equation number from the master node
1316  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1317  u_nodal_index[i]);
1318  // Get the hang weight from the master node
1319  hang_weight = hang_info_pt->master_weight(m);
1320  }
1321  // If the node is not hanging
1322  else
1323  {
1324  // Local equation number
1325  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
1326 
1327  // Node contributes with full weight
1328  hang_weight = 1.0;
1329  }
1330 
1331  // IF it's not a boundary condition
1332  if(local_eqn >= 0)
1333  {
1334  // Loop over coordinate directions
1335  for (unsigned p=0;p<DIM;p++)
1336  {
1337  // Loop over shape controlling nodes
1338  for (unsigned q=0;q<n_shape_controlling_node;q++)
1339  {
1340  // Residual x deriv of Jacobian
1341  // ----------------------------
1342 
1343  //Add the user-defined body force terms
1344  double sum = body_force[i]*testf[l];
1345 
1346  // Add the gravitational body force term
1347  sum += scaled_re_inv_fr*testf[l]*G[i];
1348 
1349  // Add the pressure gradient term
1350  sum += interpolated_p*dtestfdx(l,i);
1351 
1352  // Add in the stress tensor terms
1353  // The viscosity ratio needs to go in here to ensure
1354  // continuity of normal stress is satisfied even in flows
1355  // with zero pressure gradient!
1356  for(unsigned k=0;k<DIM;k++)
1357  {
1358  sum -= visc_ratio*
1359  (interpolated_dudx(i,k) +
1360  this->Gamma[i]*interpolated_dudx(k,i))*dtestfdx(l,k);
1361  }
1362 
1363  // Add in the inertial terms
1364 
1365  // du/dt term
1366  sum -= scaled_re_st*dudt[i]*testf[l];
1367 
1368  // Convective terms, including mesh velocity
1369  for(unsigned k=0;k<DIM;k++)
1370  {
1371  double tmp=scaled_re*interpolated_u[k];
1372  if (!this->ALE_is_disabled)
1373  {
1374  tmp-=scaled_re_st*mesh_velocity[k];
1375  }
1376  sum -= tmp*interpolated_dudx(i,k)*testf[l];
1377  }
1378 
1379  // Multiply throsugh by deriv of Jacobian and integration weight
1380  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1381  sum*dJ_dX(p,q)*w*hang_weight;
1382 
1383  // Derivative of residual x Jacobian
1384  // ---------------------------------
1385 
1386  // Body force
1387  sum=d_body_force_dx(i,p)*psif(q)*testf(l);
1388 
1389  // Pressure gradient term
1390  sum += interpolated_p*d_dtestfdx_dX(p,q,l,i);
1391 
1392  // Viscous term
1393  for (unsigned k=0;k<DIM;k++)
1394  {
1395  sum -= visc_ratio*(
1396  (interpolated_dudx(i,k)+
1397  this->Gamma[i]*interpolated_dudx(k,i))
1398  *d_dtestfdx_dX(p,q,l,k)+
1399  (d_dudx_dX(p,q,i,k) +
1400  this->Gamma[i]*d_dudx_dX(p,q,k,i))
1401  *dtestfdx(l,k));
1402  }
1403 
1404  // Convective terms, including mesh velocity
1405  for(unsigned k=0;k<DIM;k++)
1406  {
1407  double tmp=scaled_re*interpolated_u[k];
1408  if (!this->ALE_is_disabled)
1409  {
1410  tmp-=scaled_re_st*mesh_velocity[k];
1411  }
1412  sum -= tmp*d_dudx_dX(p,q,i,k)*testf(l);
1413  }
1414  if(!this->ALE_is_disabled)
1415  {
1416  sum+=scaled_re_st*pos_time_weight*
1417  psif(q)*interpolated_dudx(i,p)*testf(l);
1418  }
1419 
1420  // Multiply through by Jacobian and integration weight
1421  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1422  sum*J*w*hang_weight;
1423 
1424  } // End of loop over shape controlling nodes q
1425  } // End of loop over coordinate directions p
1426 
1427 
1428 
1429  // Derivs w.r.t. to nodal velocities
1430  // ---------------------------------
1431  if(element_has_node_with_aux_node_update_fct)
1432  {
1433  // Loop over local nodes
1434  for (unsigned q_local=0;q_local<n_node;q_local++)
1435  {
1436 
1437  // Number of master nodes and storage for the weight of
1438  // the shape function
1439  unsigned n_master2=1;
1440  double hang_weight2=1.0;
1441  HangInfo* hang_info2_pt=0;
1442 
1443  // Local boolean to indicate whether the node is hanging
1444  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1445 
1446  Node* actual_shape_controlling_node_pt=node_pt(q_local);
1447 
1448  // If the node is hanging
1449  if(is_node_hanging2)
1450  {
1451  hang_info2_pt = node_pt(q_local)->hanging_pt();
1452 
1453  // Read out number of master nodes from hanging data
1454  n_master2 = hang_info2_pt->nmaster();
1455  }
1456  // Otherwise the node is its own master
1457  else
1458  {
1459  n_master2 = 1;
1460  }
1461 
1462  // Loop over the master nodes
1463  for(unsigned mm=0;mm<n_master2;mm++)
1464  {
1465 
1466  if(is_node_hanging2)
1467  {
1468  actual_shape_controlling_node_pt=
1469  hang_info2_pt->master_node_pt(mm);
1470  hang_weight2 = hang_info2_pt->master_weight(mm);
1471  }
1472 
1473  // Find the corresponding number
1474  unsigned q=local_shape_controlling_node_lookup[
1475  actual_shape_controlling_node_pt];
1476 
1477  // Loop over coordinate directions
1478  for(unsigned p=0;p<DIM;p++)
1479  {
1480  double sum=
1481  -visc_ratio*this->Gamma[i]*dpsifdx(q_local,i)*
1482  dtestfdx(l,p)
1483  -scaled_re*psif(q_local)*interpolated_dudx(i,p)*testf(l);
1484  if (i==p)
1485  {
1486  sum-=scaled_re_st*val_time_weight*psif(q_local)*testf(l);
1487  for (unsigned k=0;k<DIM;k++)
1488  {
1489  sum-=visc_ratio*dpsifdx(q_local,k)*dtestfdx(l,k);
1490  double tmp=scaled_re*interpolated_u[k];
1491  if (!this->ALE_is_disabled)
1492  {
1493  tmp-=scaled_re_st*mesh_velocity[k];
1494  }
1495  sum-=tmp*dpsifdx(q_local,k)*testf(l);
1496  }
1497  }
1498 
1499  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1500  sum*d_U_dX(p,q)*J*w*hang_weight*hang_weight2;
1501  }
1502  } // End of loop over master nodes
1503  } // End of loop over local nodes
1504  } // End of if(element_has_node_with_aux_node_update_fct)
1505 
1506 
1507  } // local_eqn>=0
1508  }
1509  }
1510  } // End of loop over test functions
1511 
1512 
1513  // CONTINUITY EQUATION
1514  // -------------------
1515 
1516  // Loop over the shape functions
1517  for(unsigned l=0;l<n_pres;l++)
1518  {
1519 
1520  // If the pressure dof is hanging
1521  if(pressure_dof_is_hanging[l])
1522  {
1523  // Pressure dof is hanging so it must be nodal-based
1524  // Get the hang info object
1525  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1526 
1527  // Get the number of master nodes from the pressure node
1528  n_master = hang_info_pt->nmaster();
1529  }
1530  // Otherwise the node is its own master
1531  else
1532  {
1533  n_master = 1;
1534  }
1535 
1536  // Loop over the master nodes
1537  for(unsigned m=0;m<n_master;m++)
1538  {
1539  // Get the number of the unknown
1540  // If the pressure dof is hanging
1541  if(pressure_dof_is_hanging[l])
1542  {
1543  // Get the local equation from the master node
1544  local_eqn =
1545  this->local_hang_eqn(hang_info_pt->master_node_pt(m),p_index);
1546  // Get the weight for the node
1547  hang_weight = hang_info_pt->master_weight(m);
1548  }
1549  else
1550  {
1551  local_eqn = this->p_local_eqn(l);
1552  hang_weight = 1.0;
1553  }
1554 
1555  // If not a boundary conditions
1556  if(local_eqn >= 0)
1557  {
1558  // Loop over coordinate directions
1559  for (unsigned p=0;p<DIM;p++)
1560  {
1561  // Loop over nodes
1562  for (unsigned q=0;q<n_shape_controlling_node;q++)
1563  {
1564 
1565  // Residual x deriv of Jacobian
1566  //-----------------------------
1567 
1568  // Source term
1569  double aux=-source;
1570 
1571  // Loop over velocity components
1572  for(unsigned k=0;k<DIM;k++)
1573  {
1574  aux += interpolated_dudx(k,k);
1575  }
1576 
1577  // Multiply through by deriv of Jacobian and integration weight
1578  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1579  aux*dJ_dX(p,q)*testp[l]*w*hang_weight;
1580 
1581 
1582  // Derivative of residual x Jacobian
1583  // ---------------------------------
1584 
1585  // Loop over velocity components
1586  aux=-source_gradient[p]*psif(q);
1587  for(unsigned k=0;k<DIM;k++)
1588  {
1589  aux += d_dudx_dX(p,q,k,k);
1590  }
1591  // Multiply through by Jacobian and integration weight
1592  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1593  aux*testp[l]*J*w*hang_weight;
1594  }
1595  }
1596 
1597 
1598 
1599  // Derivs w.r.t. to nodal velocities
1600  // ---------------------------------
1601  if(element_has_node_with_aux_node_update_fct)
1602  {
1603  // Loop over local nodes
1604  for(unsigned q_local=0;q_local<n_node;q_local++)
1605  {
1606 
1607  // Number of master nodes and storage for the weight of
1608  // the shape function
1609  unsigned n_master2=1;
1610  double hang_weight2=1.0;
1611  HangInfo* hang_info2_pt=0;
1612 
1613  // Local boolean to indicate whether the node is hanging
1614  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1615 
1616  Node* actual_shape_controlling_node_pt=node_pt(q_local);
1617 
1618  // If the node is hanging
1619  if(is_node_hanging2)
1620  {
1621  hang_info2_pt = node_pt(q_local)->hanging_pt();
1622 
1623  // Read out number of master nodes from hanging data
1624  n_master2 = hang_info2_pt->nmaster();
1625  }
1626  // Otherwise the node is its own master
1627  else
1628  {
1629  n_master2 = 1;
1630  }
1631 
1632  // Loop over the master nodes
1633  for(unsigned mm=0;mm<n_master2;mm++)
1634  {
1635 
1636  if(is_node_hanging2)
1637  {
1638  actual_shape_controlling_node_pt=
1639  hang_info2_pt->master_node_pt(mm);
1640  hang_weight2 = hang_info2_pt->master_weight(mm);
1641  }
1642 
1643  // Find the corresponding number
1644  unsigned q=local_shape_controlling_node_lookup[
1645  actual_shape_controlling_node_pt];
1646 
1647  // Loop over coordinate directions
1648  for(unsigned p=0;p<DIM;p++)
1649  {
1650  double aux=dpsifdx(q_local,p)*testp(l);
1651  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1652  aux*d_U_dX(p,q)*J*w*hang_weight*hang_weight2;
1653  }
1654  } // End of loop over (mm) master nodes
1655  } // End of loop over local nodes q_local
1656  } // End of if(element_has_node_with_aux_node_update_fct)
1657  } // End of if it's not a boundary condition
1658  } // End of loop over (m) master nodes
1659  } // End of loop over shape functions for continuity eqn
1660 
1661  } // End of loop over integration points
1662 }
1663 
1664 //======================================================================
1665 /// 2D Further build for Crouzeix_Raviart interpolates the internal
1666 /// pressure dofs from father element: Make sure pressure values and
1667 /// dp/ds agree between fathers and sons at the midpoints of the son
1668 /// elements.
1669 //======================================================================
1670 template<>
1672 {
1673  if (this->tree_pt()->father_pt()!=0)
1674  {
1675  //Call the generic further build (if there is a father)
1677  }
1678  // Now do the PRefineableQElement further_build()
1680 
1681  // Resize internal pressure storage
1682  if (this->internal_data_pt(this->P_nst_internal_index)->nvalue()
1683  <= this->npres_nst())
1684  {
1685  this->internal_data_pt(this->P_nst_internal_index)
1686  ->resize(this->npres_nst());
1687  }
1688  else
1689  {
1690  Data* new_data_pt = new Data(this->npres_nst());
1691  delete internal_data_pt(this->P_nst_internal_index);
1692  internal_data_pt(this->P_nst_internal_index) = new_data_pt;
1693  }
1694 
1695  if(this->tree_pt()->father_pt()!=0)
1696  {
1697  // Pointer to my father (in C-R element impersonation)
1700  (quadtree_pt()->father_pt()->object_pt());
1701 
1702  // If element has same p-order as father then do the projection problem
1703  // (called after h-refinement)
1704  if(father_element_pt->p_order()==this->p_order())
1705  {
1706  using namespace QuadTreeNames;
1707 
1708  // What type of son am I? Ask my quadtree representation...
1709  int son_type=quadtree_pt()->son_type();
1710 
1711  Vector<double> s_father(2);
1712 
1713  // Son midpoint is located at the following coordinates in father element:
1714  switch(son_type)
1715  {
1716  case SW:
1717  // South west son
1718  s_father[0]=-0.5;
1719  s_father[1]=-0.5;
1720  break;
1721  case SE:
1722  // South east son
1723  s_father[0]= 0.5;
1724  s_father[1]=-0.5;
1725  break;
1726  case NE:
1727  // North east son
1728  s_father[0]= 0.5;
1729  s_father[1]= 0.5;
1730  break;
1731  case NW:
1732  // North west son
1733  s_father[0]=-0.5;
1734  s_father[1]= 0.5;
1735  break;
1736  default:
1737  throw OomphLibError(
1738  "Invalid son type in",
1739  OOMPH_CURRENT_FUNCTION,
1740  OOMPH_EXCEPTION_LOCATION);
1741  break;
1742  }
1743 
1744  // Get pressure value in father element
1745  double press=father_element_pt->interpolated_p_nst(s_father);
1746 
1747  // Reset all pressures to zero
1748  for(unsigned p=0; p<this->npres_nst(); p++)
1749  {
1750  internal_data_pt(this->P_nst_internal_index)->set_value(p,0.0);
1751  }
1752 
1753  // Set pressure values from father (BENFLAG: projection problem hack)
1754  if(this->npres_nst()==1)
1755  {
1756  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1757  }
1758  else if(this->npres_nst()==3)
1759  {
1760  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1761  internal_data_pt(this->P_nst_internal_index)->set_value(1,press);
1762  internal_data_pt(this->P_nst_internal_index)->set_value(2,press);
1763  }
1764  else
1765  {
1766  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1767  internal_data_pt(this->P_nst_internal_index)->set_value(1,press);
1768  internal_data_pt(this->P_nst_internal_index)->set_value(2,press);
1769  internal_data_pt(this->P_nst_internal_index)->set_value(3,press);
1770  }
1771  }// Otherwise this is called after p-refinement
1772  }
1773 }
1774 
1775 //======================================================================
1776 /// 2D Further build for Crouzeix_Raviart interpolates the internal
1777 /// pressure dofs from father element: Make sure pressure values and
1778 /// dp/ds agree between fathers and sons at the midpoints of the son
1779 /// elements.
1780 //======================================================================
1781 template<>
1783 {
1784  if (this->tree_pt()->father_pt()!=0)
1785  {
1786  //Call the generic further build (if there is a father)
1788  }
1789  // Now do the PRefineableQElement further_build()
1791 
1792  // Resize internal pressure storage
1793  if (this->internal_data_pt(this->P_nst_internal_index)->nvalue()
1794  <= this->npres_nst())
1795  {
1796  this->internal_data_pt(this->P_nst_internal_index)
1797  ->resize(this->npres_nst());
1798  }
1799  else
1800  {
1801  Data* new_data_pt = new Data(this->npres_nst());
1802  delete internal_data_pt(this->P_nst_internal_index);
1803  internal_data_pt(this->P_nst_internal_index) = new_data_pt;
1804  }
1805 
1806  if(this->tree_pt()->father_pt()!=0)
1807  {
1808  // Pointer to my father (in C-R element impersonation)
1811  (octree_pt()->father_pt()->object_pt());
1812 
1813  // If element has same p-order as father then do the projection problem
1814  // (called after h-refinement)
1815  if(father_element_pt->p_order()==this->p_order())
1816  {
1817  using namespace OcTreeNames;
1818 
1819  // What type of son am I? Ask my quadtree representation...
1820  int son_type=octree_pt()->son_type();
1821 
1822  Vector<double> s_father(3);
1823 
1824 
1825  // Son midpoint is located at the following coordinates in father element:
1826  for(unsigned i=0;i<3;i++)
1827  {
1828  s_father[i]=0.5*OcTree::Direction_to_vector[son_type][i];
1829  }
1830 
1831  // Get pressure value in father element
1832  double press=father_element_pt->interpolated_p_nst(s_father);
1833 
1834  // Reset all pressures to zero
1835  for(unsigned p=0; p<this->npres_nst(); p++)
1836  {
1837  internal_data_pt(this->P_nst_internal_index)->set_value(p,0.0);
1838  }
1839 
1840  // Set pressure values from father (BENFLAG: projection problem hack)
1841  if(this->npres_nst()==1)
1842  {
1843  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1844  }
1845  else
1846  {
1847  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1848  internal_data_pt(this->P_nst_internal_index)->set_value(1,press);
1849  internal_data_pt(this->P_nst_internal_index)->set_value(2,press);
1850  internal_data_pt(this->P_nst_internal_index)->set_value(3,press);
1851  internal_data_pt(this->P_nst_internal_index)->set_value(4,press);
1852  internal_data_pt(this->P_nst_internal_index)->set_value(5,press);
1853  internal_data_pt(this->P_nst_internal_index)->set_value(6,press);
1854  internal_data_pt(this->P_nst_internal_index)->set_value(7,press);
1855  }
1856  }// Otherwise this is called after p-refinement
1857  }
1858 }
1859 
1860 
1861 
1862 //====================================================================
1863 //// Force build of templates
1864 //====================================================================
1873 }
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
cstr elem_len * i
Definition: cfortran.h:607
void perform_auxiliary_node_update_fct()
Execute auxiliary update function (if any) – this can be used to update any nodal values following th...
Definition: nodes.h:1509
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
Definition: nodes.h:322
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
A Rank 4 Tensor class.
Definition: matrices.h:1625
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
bool has_auxiliary_node_update_fct_pt()
Boolean to indicate if node has a pointer to and auxiliary update function.
Definition: nodes.h:1500
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. Overwrites default implementation in FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
A Rank 3 Tensor class.
Definition: matrices.h:1337
static Vector< Vector< int > > Direction_to_vector
For each direction, i.e. a son_type (vertex), a face or an edge, this defines a vector that indicates...
Definition: octree.h:342
static char t char * s
Definition: cfortran.h:572
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void fill_in_generic_residual_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
Class that contains data for hanging nodes.
Definition: nodes.h:684
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed...
double second_invariant(const DenseMatrix< double > &tensor)
Compute second invariant of tensor.
Refineable version of Crouzeix Raviart elements. Generic class definitions.
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164