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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 /// Compute the residuals for the associated pressure advection
301 /// diffusion problem. Used by the Fp preconditioner.
302 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
303 //==============================================================
304 template<unsigned DIM>
307  Vector<double> &residuals,
308  DenseMatrix<double> &jacobian,
309  unsigned flag)
310 {
311  // Return immediately if there are no dofs
312  if (ndof()==0) return;
313 
314  //Find out how many nodes there are
315  unsigned n_node = nnode();
316 
317  //Find out how many pressure dofs there are
318  unsigned n_pres = this->npres_nst();
319 
320  //Find the indices at which the local velocities are stored
321  unsigned u_nodal_index[DIM];
322  for(unsigned i=0;i<DIM;i++) {u_nodal_index[i] = this->u_index_nst(i);}
323 
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  // pressure advection diffusion doesn't work for this one!
346  throw OomphLibError(
347  "Pressure advection diffusion does not work in this case\n",
348  OOMPH_CURRENT_FUNCTION,
349  OOMPH_EXCEPTION_LOCATION);
350 
351  for(unsigned l=0;l<n_pres;++l)
352  {pressure_dof_is_hanging[l] = false;}
353  }
354 
355  //Set up memory for the velocity shape fcts
356  Shape psif(n_node);
357  DShape dpsidx(n_node,DIM);
358 
359  //Set up memory for pressure shape and test functions
360  Shape psip(n_pres), testp(n_pres);
361  DShape dpsip(n_pres,DIM);
362  DShape dtestp(n_pres,DIM);
363 
364  //Number of integration points
365  unsigned n_intpt = integral_pt()->nweight();
366 
367  //Set the Vector to hold local coordinates
368  Vector<double> s(DIM);
369 
370  //Get Physical Variables from Element
371  //Reynolds number must be multiplied by the density ratio
372  double scaled_re = this->re()*this->density_ratio();
373 
374  //Integers to store the local equations and unknowns
375  int local_eqn=0, local_unknown=0;
376 
377 //Pointers to hang info objects
378  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
379 
380  //Loop over the integration points
381  for(unsigned ipt=0;ipt<n_intpt;ipt++)
382  {
383  //Assign values of s
384  for(unsigned i=0;i<DIM;i++) s[i] = integral_pt()->knot(ipt,i);
385 
386  //Get the integral weight
387  double w = integral_pt()->weight(ipt);
388 
389  // Call the derivatives of the veloc shape functions
390  // (Derivs not needed but they are free)
391  double J = this->dshape_eulerian_at_knot(ipt,psif,dpsidx);
392 
393  //Call the pressure shape and test functions
394  this->dpshape_and_dptest_eulerian_nst(s,psip,dpsip,testp,dtestp);
395 
396  //Premultiply the weights and the Jacobian
397  double W = w*J;
398 
399  //Calculate local values of the pressure and velocity components
400  //Allocate
401  Vector<double> interpolated_u(DIM,0.0);
402  Vector<double> interpolated_x(DIM,0.0);
403  Vector<double> interpolated_dpdx(DIM,0.0);
404 
405  //Calculate pressure gradient
406  for(unsigned l=0;l<n_pres;l++)
407  {
408  for (unsigned i=0;i<DIM;i++)
409  {
410  interpolated_dpdx[i] += this->p_nst(l)*dpsip(l,i);
411  }
412  }
413 
414  //Calculate velocities
415 
416  // Loop over nodes
417  for(unsigned l=0;l<n_node;l++)
418  {
419  //Loop over directions
420  for(unsigned i=0;i<DIM;i++)
421  {
422  //Get the nodal value
423  double u_value = nodal_value(l,u_nodal_index[i]);
424  interpolated_u[i] += u_value*psif[l];
425  interpolated_x[i] += nodal_position(l,i)*psif[l];
426  }
427  }
428 
429  // Source function (for validaton only)
430  double source=0.0;
431  if (this->Press_adv_diff_source_fct_pt!=0)
432  {
433  source=this->Press_adv_diff_source_fct_pt(interpolated_x);
434  }
435 
436 
437 
438  //Number of master nodes and storage for the weight of the shape function
439  unsigned n_master=1; double hang_weight=1.0;
440 
441 
442  //Loop over the pressure shape functions
443  for(unsigned l=0;l<n_pres;l++)
444  {
445  //If the pressure dof is hanging
446  if(pressure_dof_is_hanging[l])
447  {
448  // Pressure dof is hanging so it must be nodal-based
449  // Get the hang info object
450  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
451 
452  //Get the number of master nodes from the pressure node
453  n_master = hang_info_pt->nmaster();
454  }
455  //Otherwise the node is its own master
456  else
457  {
458  n_master = 1;
459  }
460 
461  //Loop over the master nodes
462  for(unsigned m=0;m<n_master;m++)
463  {
464  //Get the number of the unknown
465  //If the pressure dof is hanging
466  if(pressure_dof_is_hanging[l])
467  {
468  //Get the local equation from the master node
469  local_eqn =
470  this->local_hang_eqn(hang_info_pt->master_node_pt(m),p_index);
471  //Get the weight for the node
472  hang_weight = hang_info_pt->master_weight(m);
473  }
474  else
475  {
476  local_eqn = this->p_local_eqn(l);
477  hang_weight = 1.0;
478  }
479 
480  //If the equation is not pinned
481  if(local_eqn >= 0)
482  {
483  residuals[local_eqn] -= source*testp[l]*W*hang_weight;
484  for(unsigned k=0;k<DIM;k++)
485  {
486  residuals[local_eqn] += interpolated_dpdx[k]*
487  (scaled_re*interpolated_u[k]*testp[l]+dtestp(l,k))*W*hang_weight;
488  }
489 
490  // Jacobian too?
491  if(flag)
492  {
493  //Number of master nodes and weights
494  unsigned n_master2=1; double hang_weight2=1.0;
495 
496  //Loop over the pressure shape functions
497  for(unsigned l2=0;l2<n_pres;l2++)
498  {
499  //If the pressure dof is hanging
500  if(pressure_dof_is_hanging[l2])
501  {
502  hang_info2_pt =
503  this->pressure_node_pt(l2)->hanging_pt(p_index);
504  // Pressure dof is hanging so it must be nodal-based
505  //Get the number of master nodes from the pressure node
506  n_master2 = hang_info2_pt->nmaster();
507  }
508  //Otherwise the node is its own master
509  else
510  {
511  n_master2 = 1;
512  }
513 
514  //Loop over the master nodes
515  for(unsigned m2=0;m2<n_master2;m2++)
516  {
517  //Get the number of the unknown
518  //If the pressure dof is hanging
519  if(pressure_dof_is_hanging[l2])
520  {
521  //Get the unknown from the master node
522  local_unknown =
523  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
524  p_index);
525  //Get the weight from the hanging object
526  hang_weight2 = hang_info2_pt->master_weight(m2);
527  }
528  else
529  {
530  local_unknown = this->p_local_eqn(l2);
531  hang_weight2 = 1.0;
532  }
533 
534  //If the unknown is not pinned
535  if(local_unknown >= 0)
536  {
537 
538  if ((int(eqn_number(local_eqn))!=
539  this->Pinned_fp_pressure_eqn)&&
540  (int(eqn_number(local_unknown))!=
541  this->Pinned_fp_pressure_eqn))
542  {
543  for(unsigned k=0;k<DIM;k++)
544  {
545  jacobian(local_eqn,local_unknown)+=dtestp(l2,k)*
546  (scaled_re*interpolated_u[k]*testp[l]+dtestp(l,k))*
547  W*hang_weight*hang_weight2;
548  }
549  }
550  else
551  {
552  if ((int(eqn_number(local_eqn))==
553  this->Pinned_fp_pressure_eqn)&&
554  (int(eqn_number(local_unknown))
555  ==this->Pinned_fp_pressure_eqn))
556  {
557  jacobian(local_eqn,local_unknown)=1.0;
558  }
559  }
560  }
561  }
562  }
563  } /*End of Jacobian calculation*/
564  } //End of if not boundary condition
565  }//End of loop over master nodes
566  }//End of loop over l
567  }//end of integration loop
568 
569  // Now add boundary contributions from Robin BCs
570  unsigned nrobin=this->Pressure_advection_diffusion_robin_element_pt.size();
571  for (unsigned e=0;e<nrobin;e++)
572  {
573  this->Pressure_advection_diffusion_robin_element_pt[e]->
574  fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(
575  residuals,jacobian,flag);
576  }
577 }
578 
579 
580 
581 
582 //========================================================================
583 /// Add element's contribution to the elemental
584 /// residual vector and/or Jacobian matrix.
585 /// flag=1: compute both
586 /// flag=0: compute only residual vector
587 //========================================================================
588 template<unsigned DIM>
591  DenseMatrix<double> &jacobian,
592  DenseMatrix<double> &mass_matrix,
593  unsigned flag)
594 {
595 //Find out how many nodes there are
596 unsigned n_node = nnode();
597 
598 // Get continuous time from timestepper of first node
599 double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
600 
601 //Find out how many pressure dofs there are
602 unsigned n_pres = this->npres_nst();
603 
604 // Get the indices at which the velocity components are stored
605  unsigned u_nodal_index[DIM];
606  for(unsigned i=0;i<DIM;i++) {u_nodal_index[i] = this->u_index_nst(i);}
607 
608 // Which nodal value represents the pressure? (Negative if pressure
609 // is not based on nodal interpolation).
610 int p_index = this->p_nodal_index_nst();
611 
612 // Local array of booleans that are true if the l-th pressure value is
613 // hanging (avoid repeated virtual function calls)
614  bool pressure_dof_is_hanging[n_pres];
615  //If the pressure is stored at a node
616  if(p_index >= 0)
617  {
618  //Read out whether the pressure is hanging
619  for(unsigned l=0;l<n_pres;++l)
620  {
621  pressure_dof_is_hanging[l] =
622  pressure_node_pt(l)->is_hanging(p_index);
623  }
624  }
625  //Otherwise the pressure is not stored at a node and so cannot hang
626  else
627  {
628  for(unsigned l=0;l<n_pres;++l)
629  {pressure_dof_is_hanging[l] = false;}
630  }
631 
632 //Set up memory for the shape and test functions
633 Shape psif(n_node), testf(n_node);
634 DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
635 
636 
637 //Set up memory for pressure shape and test functions
638 Shape psip(n_pres), testp(n_pres);
639 
640 //Set the value of n_intpt
641 unsigned n_intpt = integral_pt()->nweight();
642 
643 //Set the Vector to hold local coordinates
644 Vector<double> s(DIM);
645 
646 //Get Physical Variables from Element
647 //Reynolds number must be multiplied by the density ratio
648 double scaled_re = this->re()*this->density_ratio();
649 double scaled_re_st = this->re_st()*this->density_ratio();
650 double scaled_re_inv_fr = this->re_invfr()*this->density_ratio();
651 double visc_ratio = this->viscosity_ratio();
652 Vector<double> G = this->g();
653 
654 //Integers that store the local equations and unknowns
655 int local_eqn=0, local_unknown=0;
656 
657 //Pointers to hang info objects
658  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
659 
660 //Local boolean for ALE (or not)
661  bool ALE_is_disabled_flag = this->ALE_is_disabled;
662 
663 //Loop over the integration points
664 for(unsigned ipt=0;ipt<n_intpt;ipt++)
665 {
666 
667  //Assign values of s
668  for(unsigned i=0;i<DIM;i++) {s[i] = integral_pt()->knot(ipt,i);}
669 
670  //Get the integral weight
671  double w = integral_pt()->weight(ipt);
672 
673  //Call the derivatives of the shape and test functions
674  double J =
675  this->dshape_and_dtest_eulerian_at_knot_nst(ipt,psif,dpsifdx,testf,dtestfdx);
676 
677  //Call the pressure shape and test functions
678  this->pshape_nst(s,psip,testp);
679 
680  //Premultiply the weights and the Jacobian
681  double W = w*J;
682 
683  //Calculate local values of the pressure and velocity components
684  //--------------------------------------------------------------
685  double interpolated_p=0.0;
686  Vector<double> interpolated_u(DIM,0.0);
687  Vector<double> interpolated_x(DIM,0.0);
688  Vector<double> mesh_veloc(DIM,0.0);
689  Vector<double> dudt(DIM,0.0);
690  DenseMatrix<double> interpolated_dudx(DIM,DIM,0.0);
691 
692  //Calculate pressure
693  for(unsigned l=0;l<n_pres;l++) {interpolated_p += this->p_nst(l)*psip[l];}
694 
695 
696  //Calculate velocities and derivatives
697 
698  // Loop over nodes
699  for(unsigned l=0;l<n_node;l++)
700  {
701  //Loop over directions
702  for(unsigned i=0;i<DIM;i++)
703  {
704  //Get the nodal value
705  double u_value = this->nodal_value(l,u_nodal_index[i]);
706  interpolated_u[i] += u_value*psif[l];
707  interpolated_x[i] += this->nodal_position(l,i)*psif[l];
708  dudt[i] += this->du_dt_nst(l,i)*psif[l];
709 
710  //Loop over derivative directions for velocity gradients
711  for(unsigned j=0;j<DIM;j++)
712  {
713  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
714  }
715  }
716  }
717 
718  if (!ALE_is_disabled_flag)
719  {
720  // Loop over nodes
721  for(unsigned l=0;l<n_node;l++)
722  {
723  //Loop over directions
724  for(unsigned i=0;i<DIM;i++)
725  {
726  mesh_veloc[i] += this->dnodal_position_dt(l,i)*psif[l];
727  }
728  }
729  }
730 
731  //Get the user-defined body force terms
732  Vector<double> body_force(DIM);
733  this->get_body_force_nst(time,ipt,s,interpolated_x,body_force);
734 
735  //Get the user-defined source function
736  double source = this->get_source_nst(time,ipt,interpolated_x);
737 
738  //MOMENTUM EQUATIONS
739  //==================
740 
741  //Number of master nodes and storage for the weight of the shape function
742  unsigned n_master=1; double hang_weight=1.0;
743 
744  //Loop over the nodes for the test functions/equations
745  //----------------------------------------------------
746  for(unsigned l=0;l<n_node;l++)
747  {
748  //Local boolean to indicate whether the node is hanging
749  bool is_node_hanging = node_pt(l)->is_hanging();
750 
751  //If the node is hanging
752  if(is_node_hanging)
753  {
754  hang_info_pt = node_pt(l)->hanging_pt();
755 
756  //Read out number of master nodes from hanging data
757  n_master = hang_info_pt->nmaster();
758  }
759  //Otherwise the node is its own master
760  else
761  {
762  n_master = 1;
763  }
764 
765  //Loop over the master nodes
766  for(unsigned m=0;m<n_master;m++)
767  {
768  // Loop over velocity components for equations
769  for(unsigned i=0;i<DIM;i++)
770  {
771  //Get the equation number
772  //If the node is hanging
773  if(is_node_hanging)
774  {
775  //Get the equation number from the master node
776  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
777  u_nodal_index[i]);
778  //Get the hang weight from the master node
779  hang_weight = hang_info_pt->master_weight(m);
780  }
781  //If the node is not hanging
782  else
783  {
784  // Local equation number
785  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
786 
787  // Node contributes with full weight
788  hang_weight = 1.0;
789  }
790 
791  //If it's not a boundary condition...
792  if(local_eqn>= 0)
793  {
794  //Temporary variable to hold the residuals
795  double sum=0.0;
796 
797  //Add the user-defined body force terms
798  sum += body_force[i];
799 
800  //Add the gravitational body force term
801  sum += scaled_re_inv_fr*G[i];
802 
803  //Add in the inertial term
804  sum -= scaled_re_st*dudt[i];
805 
806  //Convective terms, including mesh velocity
807  for(unsigned k=0;k<DIM;k++)
808  {
809  double tmp=scaled_re*interpolated_u[k];
810  if (!ALE_is_disabled_flag)
811  {tmp -= scaled_re_st*mesh_veloc[k];}
812  sum -= tmp*interpolated_dudx(i,k);
813  }
814 
815  //Add the pressure gradient term
816  sum = (sum*testf[l] + interpolated_p*dtestfdx(l,i))*W*hang_weight;
817 
818  //Add in the stress tensor terms
819  //The viscosity ratio needs to go in here to ensure
820  //continuity of normal stress is satisfied even in flows
821  //with zero pressure gradient!
822  for(unsigned k=0;k<DIM;k++)
823  {
824  sum -= visc_ratio*
825  (interpolated_dudx(i,k) + this->Gamma[i]*interpolated_dudx(k,i))
826  *dtestfdx(l,k)*W*hang_weight;
827  }
828 
829  // Add contribution
830  residuals[local_eqn] += sum;
831 
832  //CALCULATE THE JACOBIAN
833  if(flag)
834  {
835  //Number of master nodes and weights
836  unsigned n_master2=1; double hang_weight2=1.0;
837  //Loop over the velocity nodes for columns
838  for(unsigned l2=0;l2<n_node;l2++)
839  {
840  //Local boolean to indicate whether the node is hanging
841  bool is_node2_hanging = node_pt(l2)->is_hanging();
842 
843  //If the node is hanging
844  if(is_node2_hanging)
845  {
846  hang_info2_pt = node_pt(l2)->hanging_pt();
847  //Read out number of master nodes from hanging data
848  n_master2 = hang_info2_pt->nmaster();
849  }
850  //Otherwise the node is its own master
851  else
852  {
853  n_master2 = 1;
854  }
855 
856  //Loop over the master nodes
857  for(unsigned m2=0;m2<n_master2;m2++)
858  {
859  //Loop over the velocity components
860  for(unsigned i2=0;i2<DIM;i2++)
861  {
862  //Get the number of the unknown
863  //If the node is hanging
864  if(is_node2_hanging)
865  {
866  //Get the equation number from the master node
867  local_unknown =
868  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
869  u_nodal_index[i2]);
870  //Get the hang weights
871  hang_weight2 = hang_info2_pt->master_weight(m2);
872  }
873  else
874  {
875  local_unknown = this->nodal_local_eqn(l2,u_nodal_index[i2]);
876  hang_weight2 = 1.0;
877  }
878 
879  // If the unknown is non-zero
880  if(local_unknown >= 0)
881  {
882  //Add contribution to Elemental Matrix
883  jacobian(local_eqn,local_unknown)
884  -= visc_ratio*this->Gamma[i]*dpsifdx(l2,i)*
885  dtestfdx(l,i2)*W*hang_weight*hang_weight2;
886 
887  //Now add in the inertial terms
888  jacobian(local_eqn,local_unknown)
889  -= scaled_re*psif[l2]*interpolated_dudx(i,i2)*testf[l]*W*
890  hang_weight*hang_weight2;
891 
892  //Extra diagonal components if i2=i
893  if(i2 == i)
894  {
895  //Mass matrix entries
896  //Again note the positive sign because the mass
897  //matrix is taken on the other side of the equation
898  if(flag==2)
899  {
900  mass_matrix(local_eqn,local_unknown) +=
901  scaled_re_st*psif[l2]*testf[l]*W*
902  hang_weight*hang_weight2;
903  }
904 
905  // du/dt term
906  jacobian(local_eqn,local_unknown)
907  -= scaled_re_st*
908  node_pt(l2)->time_stepper_pt()->weight(1,0)*
909  psif[l2]*testf[l]*W*hang_weight*hang_weight2;
910 
911  //Extra advective terms
912  for(unsigned k=0;k<DIM;k++)
913  {
914  double tmp=scaled_re*interpolated_u[k];
915  if (!ALE_is_disabled_flag)
916  {tmp -= scaled_re_st*mesh_veloc[k];}
917 
918  jacobian(local_eqn,local_unknown) -=
919  tmp*dpsifdx(l2,k)*testf[l]*W*hang_weight*hang_weight2;
920  }
921 
922  // Extra viscous terms
923  for(unsigned k=0;k<DIM;k++)
924  {
925  jacobian(local_eqn,local_unknown)
926  -= visc_ratio*dpsifdx(l2,k)*
927  dtestfdx(l,k)*W*hang_weight*hang_weight2;
928  }
929 
930  }
931  }
932  }
933  }
934  }
935 
936  //Loop over the pressure shape functions
937  for(unsigned l2=0;l2<n_pres;l2++)
938  {
939  //If the pressure dof is hanging
940  if(pressure_dof_is_hanging[l2])
941  {
942  hang_info2_pt =
943  this->pressure_node_pt(l2)->hanging_pt(p_index);
944  // Pressure dof is hanging so it must be nodal-based
945  //Get the number of master nodes from the pressure node
946  n_master2 = hang_info2_pt->nmaster();
947  }
948  //Otherwise the node is its own master
949  else
950  {
951  n_master2 = 1;
952  }
953 
954  //Loop over the master nodes
955  for(unsigned m2=0;m2<n_master2;m2++)
956  {
957  //Get the number of the unknown
958  //If the pressure dof is hanging
959  if(pressure_dof_is_hanging[l2])
960  {
961  //Get the unknown from the master node
962  local_unknown =
963  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
964  p_index);
965  //Get the weight from the hanging object
966  hang_weight2 = hang_info2_pt->master_weight(m2);
967  }
968  else
969  {
970  local_unknown = this->p_local_eqn(l2);
971  hang_weight2 = 1.0;
972  }
973 
974  //If the unknown is not pinned
975  if(local_unknown >= 0)
976  {
977  jacobian(local_eqn,local_unknown)
978  += psip[l2]*dtestfdx(l,i)*W*hang_weight*hang_weight2;
979  }
980  }
981  }
982 
983  }// End of Jacobian calculation
984 
985  } //End of if not boundary condition statement
986 
987  } //End of loop over components of non-hanging node
988 
989  } //End of loop over master nodes
990 
991  } // End of loop over nodes for equations
992 
993 
994 
995  //CONTINUITY EQUATION
996  //===================
997 
998  //Loop over the pressure shape functions
999  for(unsigned l=0;l<n_pres;l++)
1000  {
1001  //If the pressure dof is hanging
1002  if(pressure_dof_is_hanging[l])
1003  {
1004  // Pressure dof is hanging so it must be nodal-based
1005  // Get the hang info object
1006  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1007 
1008  //Get the number of master nodes from the pressure node
1009  n_master = hang_info_pt->nmaster();
1010  }
1011  //Otherwise the node is its own master
1012  else
1013  {
1014  n_master = 1;
1015  }
1016 
1017  //Loop over the master nodes
1018  for(unsigned m=0;m<n_master;m++)
1019  {
1020  //Get the number of the unknown
1021  //If the pressure dof is hanging
1022  if(pressure_dof_is_hanging[l])
1023  {
1024  //Get the local equation from the master node
1025  local_eqn =
1026  this->local_hang_eqn(hang_info_pt->master_node_pt(m),p_index);
1027  //Get the weight for the node
1028  hang_weight = hang_info_pt->master_weight(m);
1029  }
1030  else
1031  {
1032  local_eqn = this->p_local_eqn(l);
1033  hang_weight = 1.0;
1034  }
1035 
1036  //If the equation is not pinned
1037  if(local_eqn >= 0)
1038  {
1039  // Source term
1040  residuals[local_eqn] -= source*testp[l]*W*hang_weight;
1041 
1042  // Loop over velocity components
1043  for(unsigned k=0;k<DIM;k++)
1044  {
1045  residuals[local_eqn] += interpolated_dudx(k,k)*testp[l]*W*hang_weight;
1046  }
1047 
1048  //CALCULATE THE JACOBIAN
1049  if(flag)
1050  {
1051  unsigned n_master2=1; double hang_weight2=1.0;
1052  //Loop over the velocity nodes for columns
1053  for(unsigned l2=0;l2<n_node;l2++)
1054  {
1055  //Local boolean to indicate whether the node is hanging
1056  bool is_node2_hanging = node_pt(l2)->is_hanging();
1057 
1058  //If the node is hanging
1059  if(is_node2_hanging)
1060  {
1061  hang_info2_pt = node_pt(l2)->hanging_pt();
1062  //Read out number of master nodes from hanging data
1063  n_master2 = hang_info2_pt->nmaster();
1064  }
1065  //Otherwise the node is its own master
1066  else
1067  {
1068  n_master2 = 1;
1069  }
1070 
1071  //Loop over the master nodes
1072  for(unsigned m2=0;m2<n_master2;m2++)
1073  {
1074  //Loop over the velocity components
1075  for(unsigned i2=0;i2<DIM;i2++)
1076  {
1077  //Get the number of the unknown
1078  //If the node is hanging
1079  if(is_node2_hanging)
1080  {
1081  //Get the equation number from the master node
1082  local_unknown =
1083  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
1084  u_nodal_index[i2]);
1085  hang_weight2 = hang_info2_pt->master_weight(m2);
1086  }
1087  else
1088  {
1089  local_unknown = this->nodal_local_eqn(l2,u_nodal_index[i2]);
1090  hang_weight2 = 1.0;
1091  }
1092 
1093  //If the unknown is not pinned
1094  if(local_unknown >= 0)
1095  {
1096  jacobian(local_eqn,local_unknown)
1097  += dpsifdx(l2,i2)*testp[l]*W*hang_weight*hang_weight2;
1098  }
1099  }
1100  }
1101  }
1102 
1103  // NO PRESSURE CONTRIBUTION TO THE JACOBIAN
1104 
1105  } //End of jacobian calculation
1106  }
1107  }
1108  } //End of loop over pressure variables
1109 
1110 } //End of loop over integration points
1111 }
1112 
1113 
1114 
1115 
1116 //======================================================================
1117 /// Compute derivatives of elemental residual vector with respect
1118 /// to nodal coordinates.
1119 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
1120 /// Overloads the FD-based version in the FE base class.
1121 //======================================================================
1122 template <unsigned DIM>
1125  dresidual_dnodal_coordinates)
1126 {
1127  // Return immediately if there are no dofs
1128  if(ndof()==0) { return; }
1129 
1130  // Determine number of nodes in element
1131  const unsigned n_node = nnode();
1132 
1133  // Get continuous time from timestepper of first node
1134  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
1135 
1136  // Determine number of pressure dofs in element
1137  const unsigned n_pres = this->npres_nst();
1138 
1139  // Find the indices at which the local velocities are stored
1140  unsigned u_nodal_index[DIM];
1141  for(unsigned i=0;i<DIM;i++) { u_nodal_index[i] = this->u_index_nst(i); }
1142 
1143  // Which nodal value represents the pressure? (Negative if pressure
1144  // is not based on nodal interpolation).
1145  const int p_index = this->p_nodal_index_nst();
1146 
1147  // Local array of booleans that are true if the l-th pressure value is
1148  // hanging (avoid repeated virtual function calls)
1149  bool pressure_dof_is_hanging[n_pres];
1150 
1151  // If the pressure is stored at a node
1152  if(p_index >= 0)
1153  {
1154  // Read out whether the pressure is hanging
1155  for(unsigned l=0;l<n_pres;++l)
1156  {
1157  pressure_dof_is_hanging[l] =
1158  pressure_node_pt(l)->is_hanging(p_index);
1159  }
1160  }
1161  // Otherwise the pressure is not stored at a node and so cannot hang
1162  else
1163  {
1164  for(unsigned l=0;l<n_pres;++l) { pressure_dof_is_hanging[l] = false; }
1165  }
1166 
1167  // Set up memory for the shape and test functions
1168  Shape psif(n_node), testf(n_node);
1169  DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
1170 
1171  // Set up memory for pressure shape and test functions
1172  Shape psip(n_pres), testp(n_pres);
1173 
1174  // Determine number of shape controlling nodes
1175  const unsigned n_shape_controlling_node = nshape_controlling_nodes();
1176 
1177  // Deriatives of shape fct derivatives w.r.t. nodal coords
1178  RankFourTensor<double> d_dpsifdx_dX(DIM,n_shape_controlling_node,n_node,DIM);
1179  RankFourTensor<double> d_dtestfdx_dX(DIM,n_shape_controlling_node,n_node,DIM);
1180 
1181  // Derivative of Jacobian of mapping w.r.t. to nodal coords
1182  DenseMatrix<double> dJ_dX(DIM,n_shape_controlling_node);
1183 
1184  // Derivatives of derivative of u w.r.t. nodal coords
1185  RankFourTensor<double> d_dudx_dX(DIM,n_shape_controlling_node,DIM,DIM);
1186 
1187  // Derivatives of nodal velocities w.r.t. nodal coords:
1188  // Assumption: Interaction only local via no-slip so
1189  // X_ij only affects U_ij.
1190  DenseMatrix<double> d_U_dX(DIM,n_shape_controlling_node,0.0);
1191 
1192  // Determine the number of integration points
1193  const unsigned n_intpt = integral_pt()->nweight();
1194 
1195  // Vector to hold local coordinates
1196  Vector<double> s(DIM);
1197 
1198  // Get physical variables from element
1199  // (Reynolds number must be multiplied by the density ratio)
1200  double scaled_re = this->re()*this->density_ratio();
1201  double scaled_re_st = this->re_st()*this->density_ratio();
1202  double scaled_re_inv_fr = this->re_invfr()*this->density_ratio();
1203  double visc_ratio = this->viscosity_ratio();
1204  Vector<double> G = this->g();
1205 
1206  // FD step
1208 
1209  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1210  // Assumption: Interaction only local via no-slip so
1211  // X_ij only affects U_ij.
1212  bool element_has_node_with_aux_node_update_fct=false;
1213 
1214  std::map<Node*,unsigned> local_shape_controlling_node_lookup=
1215  shape_controlling_node_lookup();
1216 
1217  // FD loop over shape-controlling nodes
1218  for(std::map<Node*,unsigned>::iterator it=
1219  local_shape_controlling_node_lookup.begin();
1220  it!=local_shape_controlling_node_lookup.end();
1221  it++)
1222  {
1223  // Get node
1224  Node* nod_pt=it->first;
1225 
1226  // Get its number
1227  unsigned q=it->second;
1228 
1229  // Only compute if there's a node-update fct involved
1230  if(nod_pt->has_auxiliary_node_update_fct_pt())
1231  {
1232  element_has_node_with_aux_node_update_fct=true;
1233 
1234  // Current nodal velocity
1235  Vector<double> u_ref(DIM);
1236  for(unsigned i=0;i<DIM;i++)
1237  {
1238  u_ref[i]=*(nod_pt->value_pt(u_nodal_index[i]));
1239  }
1240 
1241  // FD
1242  for(unsigned p=0;p<DIM;p++)
1243  {
1244  // Make backup
1245  double backup=nod_pt->x(p);
1246 
1247  // Do FD step. No node update required as we're
1248  // attacking the coordinate directly...
1249  nod_pt->x(p)+=eps_fd;
1250 
1251  // Do auxiliary node update (to apply no slip)
1253 
1254  // Evaluate
1255  d_U_dX(p,q)=(*(nod_pt->value_pt(u_nodal_index[p]))-u_ref[p])/eps_fd;
1256 
1257  // Reset
1258  nod_pt->x(p)=backup;
1259 
1260  // Do auxiliary node update (to apply no slip)
1262  }
1263  }
1264  }
1265 
1266  // Integer to store the local equation number
1267  int local_eqn=0;
1268 
1269  // Pointers to hang info object
1270  HangInfo *hang_info_pt=0;
1271 
1272  // Loop over the integration points
1273  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1274  {
1275  // Assign values of s
1276  for(unsigned i=0;i<DIM;i++) { s[i] = integral_pt()->knot(ipt,i); }
1277 
1278  // Get the integral weight
1279  const double w = integral_pt()->weight(ipt);
1280 
1281  // Call the derivatives of the shape and test functions
1282  const double J = this->dshape_and_dtest_eulerian_at_knot_nst(
1283  ipt,psif,dpsifdx,d_dpsifdx_dX,testf,dtestfdx,d_dtestfdx_dX,dJ_dX);
1284 
1285  // Call the pressure shape and test functions
1286  this->pshape_nst(s,psip,testp);
1287 
1288  // Calculate local values of the pressure and velocity components
1289  // Allocate
1290  double interpolated_p=0.0;
1291  Vector<double> interpolated_u(DIM,0.0);
1292  Vector<double> interpolated_x(DIM,0.0);
1293  Vector<double> mesh_velocity(DIM,0.0);
1294  Vector<double> dudt(DIM,0.0);
1295  DenseMatrix<double> interpolated_dudx(DIM,DIM,0.0);
1296 
1297  // Calculate pressure
1298  for(unsigned l=0;l<n_pres;l++) { interpolated_p += this->p_nst(l)*psip[l]; }
1299 
1300  // Calculate velocities and derivatives:
1301 
1302  // Loop over nodes
1303  for(unsigned l=0;l<n_node;l++)
1304  {
1305  // Loop over directions
1306  for(unsigned i=0;i<DIM;i++)
1307  {
1308  // Get the nodal value
1309  const double u_value = nodal_value(l,u_nodal_index[i]);
1310  interpolated_u[i] += u_value*psif[l];
1311  interpolated_x[i] += nodal_position(l,i)*psif[l];
1312  dudt[i] += this->du_dt_nst(l,i)*psif[l];
1313 
1314  // Loop over derivative directions
1315  for(unsigned j=0;j<DIM;j++)
1316  {
1317  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
1318  }
1319  }
1320  }
1321 
1322  if(!this->ALE_is_disabled)
1323  {
1324  // Loop over nodes
1325  for(unsigned l=0;l<n_node;l++)
1326  {
1327  // Loop over directions
1328  for(unsigned i=0;i<DIM;i++)
1329  {
1330  mesh_velocity[i] += this->dnodal_position_dt(l,i)*psif[l];
1331  }
1332  }
1333  }
1334 
1335  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1336 
1337  // Loop over shape-controlling nodes
1338  for(unsigned q=0;q<n_shape_controlling_node;q++)
1339  {
1340  // Loop over coordinate directions
1341  for(unsigned p=0;p<DIM;p++)
1342  {
1343  for(unsigned i=0;i<DIM;i++)
1344  {
1345  for(unsigned k=0;k<DIM;k++)
1346  {
1347  double aux=0.0;
1348  for(unsigned j=0;j<n_node;j++)
1349  {
1350  aux += nodal_value(j,u_nodal_index[i])*d_dpsifdx_dX(p,q,j,k);
1351  }
1352  d_dudx_dX(p,q,i,k) = aux;
1353  }
1354  }
1355  }
1356  }
1357 
1358  // Get weight of actual nodal position/value in computation of mesh
1359  // velocity from positional/value time stepper
1360  const double pos_time_weight
1361  = node_pt(0)->position_time_stepper_pt()->weight(1,0);
1362  const double val_time_weight = node_pt(0)->time_stepper_pt()->weight(1,0);
1363 
1364  // Get the user-defined body force terms
1365  Vector<double> body_force(DIM);
1366  this->get_body_force_nst(time,ipt,s,interpolated_x,body_force);
1367 
1368  // Get the user-defined source function
1369  const double source = this->get_source_nst(time,ipt,interpolated_x);
1370 
1371  // Get gradient of body force function
1372  DenseMatrix<double> d_body_force_dx(DIM,DIM,0.0);
1373  this->get_body_force_gradient_nst(time,ipt,s,
1374  interpolated_x, d_body_force_dx);
1375 
1376  // Get gradient of source function
1377  Vector<double> source_gradient(DIM,0.0);
1378  this->get_source_gradient_nst(time,ipt,interpolated_x, source_gradient);
1379 
1380 
1381  // Assemble shape derivatives
1382  //---------------------------
1383 
1384  // MOMENTUM EQUATIONS
1385  // ------------------
1386 
1387  // Number of master nodes and storage for the weight of the shape function
1388  unsigned n_master=1; double hang_weight=1.0;
1389 
1390  // Loop over the test functions
1391  for(unsigned l=0;l<n_node;l++)
1392  {
1393 
1394  // Local boolean to indicate whether the node is hanging
1395  bool is_node_hanging = node_pt(l)->is_hanging();
1396 
1397  // If the node is hanging
1398  if(is_node_hanging)
1399  {
1400  hang_info_pt = node_pt(l)->hanging_pt();
1401 
1402  // Read out number of master nodes from hanging data
1403  n_master = hang_info_pt->nmaster();
1404  }
1405  // Otherwise the node is its own master
1406  else
1407  {
1408  n_master = 1;
1409  }
1410 
1411  // Loop over the master nodes
1412  for(unsigned m=0;m<n_master;m++)
1413  {
1414 
1415  // Loop over coordinate directions
1416  for(unsigned i=0;i<DIM;i++)
1417  {
1418 
1419  // Get the equation number
1420  // If the node is hanging
1421  if(is_node_hanging)
1422  {
1423  // Get the equation number from the master node
1424  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1425  u_nodal_index[i]);
1426  // Get the hang weight from the master node
1427  hang_weight = hang_info_pt->master_weight(m);
1428  }
1429  // If the node is not hanging
1430  else
1431  {
1432  // Local equation number
1433  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
1434 
1435  // Node contributes with full weight
1436  hang_weight = 1.0;
1437  }
1438 
1439  // IF it's not a boundary condition
1440  if(local_eqn >= 0)
1441  {
1442  // Loop over coordinate directions
1443  for (unsigned p=0;p<DIM;p++)
1444  {
1445  // Loop over shape controlling nodes
1446  for (unsigned q=0;q<n_shape_controlling_node;q++)
1447  {
1448  // Residual x deriv of Jacobian
1449  // ----------------------------
1450 
1451  //Add the user-defined body force terms
1452  double sum = body_force[i]*testf[l];
1453 
1454  // Add the gravitational body force term
1455  sum += scaled_re_inv_fr*testf[l]*G[i];
1456 
1457  // Add the pressure gradient term
1458  sum += interpolated_p*dtestfdx(l,i);
1459 
1460  // Add in the stress tensor terms
1461  // The viscosity ratio needs to go in here to ensure
1462  // continuity of normal stress is satisfied even in flows
1463  // with zero pressure gradient!
1464  for(unsigned k=0;k<DIM;k++)
1465  {
1466  sum -= visc_ratio*
1467  (interpolated_dudx(i,k) +
1468  this->Gamma[i]*interpolated_dudx(k,i))*dtestfdx(l,k);
1469  }
1470 
1471  // Add in the inertial terms
1472 
1473  // du/dt term
1474  sum -= scaled_re_st*dudt[i]*testf[l];
1475 
1476  // Convective terms, including mesh velocity
1477  for(unsigned k=0;k<DIM;k++)
1478  {
1479  double tmp=scaled_re*interpolated_u[k];
1480  if (!this->ALE_is_disabled)
1481  {
1482  tmp-=scaled_re_st*mesh_velocity[k];
1483  }
1484  sum -= tmp*interpolated_dudx(i,k)*testf[l];
1485  }
1486 
1487  // Multiply throsugh by deriv of Jacobian and integration weight
1488  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1489  sum*dJ_dX(p,q)*w*hang_weight;
1490 
1491  // Derivative of residual x Jacobian
1492  // ---------------------------------
1493 
1494  // Body force
1495  sum=d_body_force_dx(i,p)*psif(q)*testf(l);
1496 
1497  // Pressure gradient term
1498  sum += interpolated_p*d_dtestfdx_dX(p,q,l,i);
1499 
1500  // Viscous term
1501  for (unsigned k=0;k<DIM;k++)
1502  {
1503  sum -= visc_ratio*(
1504  (interpolated_dudx(i,k)+
1505  this->Gamma[i]*interpolated_dudx(k,i))
1506  *d_dtestfdx_dX(p,q,l,k)+
1507  (d_dudx_dX(p,q,i,k) +
1508  this->Gamma[i]*d_dudx_dX(p,q,k,i))
1509  *dtestfdx(l,k));
1510  }
1511 
1512  // Convective terms, including mesh velocity
1513  for(unsigned k=0;k<DIM;k++)
1514  {
1515  double tmp=scaled_re*interpolated_u[k];
1516  if (!this->ALE_is_disabled)
1517  {
1518  tmp-=scaled_re_st*mesh_velocity[k];
1519  }
1520  sum -= tmp*d_dudx_dX(p,q,i,k)*testf(l);
1521  }
1522  if(!this->ALE_is_disabled)
1523  {
1524  sum+=scaled_re_st*pos_time_weight*
1525  psif(q)*interpolated_dudx(i,p)*testf(l);
1526  }
1527 
1528  // Multiply through by Jacobian and integration weight
1529  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1530  sum*J*w*hang_weight;
1531 
1532  } // End of loop over shape controlling nodes q
1533  } // End of loop over coordinate directions p
1534 
1535 
1536 
1537  // Derivs w.r.t. to nodal velocities
1538  // ---------------------------------
1539  if(element_has_node_with_aux_node_update_fct)
1540  {
1541  // Loop over local nodes
1542  for (unsigned q_local=0;q_local<n_node;q_local++)
1543  {
1544 
1545  // Number of master nodes and storage for the weight of
1546  // the shape function
1547  unsigned n_master2=1;
1548  double hang_weight2=1.0;
1549  HangInfo* hang_info2_pt=0;
1550 
1551  // Local boolean to indicate whether the node is hanging
1552  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1553 
1554  Node* actual_shape_controlling_node_pt=node_pt(q_local);
1555 
1556  // If the node is hanging
1557  if(is_node_hanging2)
1558  {
1559  hang_info2_pt = node_pt(q_local)->hanging_pt();
1560 
1561  // Read out number of master nodes from hanging data
1562  n_master2 = hang_info2_pt->nmaster();
1563  }
1564  // Otherwise the node is its own master
1565  else
1566  {
1567  n_master2 = 1;
1568  }
1569 
1570  // Loop over the master nodes
1571  for(unsigned mm=0;mm<n_master2;mm++)
1572  {
1573 
1574  if(is_node_hanging2)
1575  {
1576  actual_shape_controlling_node_pt=
1577  hang_info2_pt->master_node_pt(mm);
1578  hang_weight2 = hang_info2_pt->master_weight(mm);
1579  }
1580 
1581  // Find the corresponding number
1582  unsigned q=local_shape_controlling_node_lookup[
1583  actual_shape_controlling_node_pt];
1584 
1585  // Loop over coordinate directions
1586  for(unsigned p=0;p<DIM;p++)
1587  {
1588  double sum=
1589  -visc_ratio*this->Gamma[i]*dpsifdx(q_local,i)*
1590  dtestfdx(l,p)
1591  -scaled_re*psif(q_local)*interpolated_dudx(i,p)*testf(l);
1592  if (i==p)
1593  {
1594  sum-=scaled_re_st*val_time_weight*psif(q_local)*testf(l);
1595  for (unsigned k=0;k<DIM;k++)
1596  {
1597  sum-=visc_ratio*dpsifdx(q_local,k)*dtestfdx(l,k);
1598  double tmp=scaled_re*interpolated_u[k];
1599  if (!this->ALE_is_disabled)
1600  {
1601  tmp-=scaled_re_st*mesh_velocity[k];
1602  }
1603  sum-=tmp*dpsifdx(q_local,k)*testf(l);
1604  }
1605  }
1606 
1607  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1608  sum*d_U_dX(p,q)*J*w*hang_weight*hang_weight2;
1609  }
1610  } // End of loop over master nodes
1611  } // End of loop over local nodes
1612  } // End of if(element_has_node_with_aux_node_update_fct)
1613 
1614 
1615  } // local_eqn>=0
1616  }
1617  }
1618  } // End of loop over test functions
1619 
1620 
1621  // CONTINUITY EQUATION
1622  // -------------------
1623 
1624  // Loop over the shape functions
1625  for(unsigned l=0;l<n_pres;l++)
1626  {
1627 
1628  // If the pressure dof is hanging
1629  if(pressure_dof_is_hanging[l])
1630  {
1631  // Pressure dof is hanging so it must be nodal-based
1632  // Get the hang info object
1633  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1634 
1635  // Get the number of master nodes from the pressure node
1636  n_master = hang_info_pt->nmaster();
1637  }
1638  // Otherwise the node is its own master
1639  else
1640  {
1641  n_master = 1;
1642  }
1643 
1644  // Loop over the master nodes
1645  for(unsigned m=0;m<n_master;m++)
1646  {
1647  // Get the number of the unknown
1648  // If the pressure dof is hanging
1649  if(pressure_dof_is_hanging[l])
1650  {
1651  // Get the local equation from the master node
1652  local_eqn =
1653  this->local_hang_eqn(hang_info_pt->master_node_pt(m),p_index);
1654  // Get the weight for the node
1655  hang_weight = hang_info_pt->master_weight(m);
1656  }
1657  else
1658  {
1659  local_eqn = this->p_local_eqn(l);
1660  hang_weight = 1.0;
1661  }
1662 
1663  // If not a boundary conditions
1664  if(local_eqn >= 0)
1665  {
1666  // Loop over coordinate directions
1667  for (unsigned p=0;p<DIM;p++)
1668  {
1669  // Loop over nodes
1670  for (unsigned q=0;q<n_shape_controlling_node;q++)
1671  {
1672 
1673  // Residual x deriv of Jacobian
1674  //-----------------------------
1675 
1676  // Source term
1677  double aux=-source;
1678 
1679  // Loop over velocity components
1680  for(unsigned k=0;k<DIM;k++)
1681  {
1682  aux += interpolated_dudx(k,k);
1683  }
1684 
1685  // Multiply through by deriv of Jacobian and integration weight
1686  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1687  aux*dJ_dX(p,q)*testp[l]*w*hang_weight;
1688 
1689 
1690  // Derivative of residual x Jacobian
1691  // ---------------------------------
1692 
1693  // Loop over velocity components
1694  aux=-source_gradient[p]*psif(q);
1695  for(unsigned k=0;k<DIM;k++)
1696  {
1697  aux += d_dudx_dX(p,q,k,k);
1698  }
1699  // Multiply through by Jacobian and integration weight
1700  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1701  aux*testp[l]*J*w*hang_weight;
1702  }
1703  }
1704 
1705 
1706 
1707  // Derivs w.r.t. to nodal velocities
1708  // ---------------------------------
1709  if(element_has_node_with_aux_node_update_fct)
1710  {
1711  // Loop over local nodes
1712  for(unsigned q_local=0;q_local<n_node;q_local++)
1713  {
1714 
1715  // Number of master nodes and storage for the weight of
1716  // the shape function
1717  unsigned n_master2=1;
1718  double hang_weight2=1.0;
1719  HangInfo* hang_info2_pt=0;
1720 
1721  // Local boolean to indicate whether the node is hanging
1722  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1723 
1724  Node* actual_shape_controlling_node_pt=node_pt(q_local);
1725 
1726  // If the node is hanging
1727  if(is_node_hanging2)
1728  {
1729  hang_info2_pt = node_pt(q_local)->hanging_pt();
1730 
1731  // Read out number of master nodes from hanging data
1732  n_master2 = hang_info2_pt->nmaster();
1733  }
1734  // Otherwise the node is its own master
1735  else
1736  {
1737  n_master2 = 1;
1738  }
1739 
1740  // Loop over the master nodes
1741  for(unsigned mm=0;mm<n_master2;mm++)
1742  {
1743 
1744  if(is_node_hanging2)
1745  {
1746  actual_shape_controlling_node_pt=
1747  hang_info2_pt->master_node_pt(mm);
1748  hang_weight2 = hang_info2_pt->master_weight(mm);
1749  }
1750 
1751  // Find the corresponding number
1752  unsigned q=local_shape_controlling_node_lookup[
1753  actual_shape_controlling_node_pt];
1754 
1755  // Loop over coordinate directions
1756  for(unsigned p=0;p<DIM;p++)
1757  {
1758  double aux=dpsifdx(q_local,p)*testp(l);
1759  dresidual_dnodal_coordinates(local_eqn,p,q)+=
1760  aux*d_U_dX(p,q)*J*w*hang_weight*hang_weight2;
1761  }
1762  } // End of loop over (mm) master nodes
1763  } // End of loop over local nodes q_local
1764  } // End of if(element_has_node_with_aux_node_update_fct)
1765  } // End of if it's not a boundary condition
1766  } // End of loop over (m) master nodes
1767  } // End of loop over shape functions for continuity eqn
1768 
1769  } // End of loop over integration points
1770 }
1771 
1772 //======================================================================
1773 /// 2D Further build for Crouzeix_Raviart interpolates the internal
1774 /// pressure dofs from father element: Make sure pressure values and
1775 /// dp/ds agree between fathers and sons at the midpoints of the son
1776 /// elements.
1777 //======================================================================
1778 template<>
1780 {
1781  if (this->tree_pt()->father_pt()!=0)
1782  {
1783  //Call the generic further build (if there is a father)
1785  }
1786  // Now do the PRefineableQElement further_build()
1788 
1789  // Resize internal pressure storage
1790  if (this->internal_data_pt(this->P_nst_internal_index)->nvalue()
1791  <= this->npres_nst())
1792  {
1793  this->internal_data_pt(this->P_nst_internal_index)
1794  ->resize(this->npres_nst());
1795  }
1796  else
1797  {
1798  Data* new_data_pt = new Data(this->npres_nst());
1799  delete internal_data_pt(this->P_nst_internal_index);
1800  internal_data_pt(this->P_nst_internal_index) = new_data_pt;
1801  }
1802 
1803  if(this->tree_pt()->father_pt()!=0)
1804  {
1805  // Pointer to my father (in C-R element impersonation)
1806  PRefineableQCrouzeixRaviartElement<2>* father_element_pt=
1808  (quadtree_pt()->father_pt()->object_pt());
1809 
1810  // If element has same p-order as father then do the projection problem
1811  // (called after h-refinement)
1812  if(father_element_pt->p_order()==this->p_order())
1813  {
1814  using namespace QuadTreeNames;
1815 
1816  // What type of son am I? Ask my quadtree representation...
1817  int son_type=quadtree_pt()->son_type();
1818 
1819  Vector<double> s_father(2);
1820 
1821  // Son midpoint is located at the following coordinates in father element:
1822  switch(son_type)
1823  {
1824  case SW:
1825  // South west son
1826  s_father[0]=-0.5;
1827  s_father[1]=-0.5;
1828  break;
1829  case SE:
1830  // South east son
1831  s_father[0]= 0.5;
1832  s_father[1]=-0.5;
1833  break;
1834  case NE:
1835  // North east son
1836  s_father[0]= 0.5;
1837  s_father[1]= 0.5;
1838  break;
1839  case NW:
1840  // North west son
1841  s_father[0]=-0.5;
1842  s_father[1]= 0.5;
1843  break;
1844  default:
1845  throw OomphLibError(
1846  "Invalid son type in",
1847  OOMPH_CURRENT_FUNCTION,
1848  OOMPH_EXCEPTION_LOCATION);
1849  break;
1850  }
1851 
1852  // Get pressure value in father element
1853  double press=father_element_pt->interpolated_p_nst(s_father);
1854 
1855  // Reset all pressures to zero
1856  for(unsigned p=0; p<this->npres_nst(); p++)
1857  {
1858  internal_data_pt(this->P_nst_internal_index)->set_value(p,0.0);
1859  }
1860 
1861  // Set pressure values from father (BENFLAG: projection problem hack)
1862  if(this->npres_nst()==1)
1863  {
1864  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1865  }
1866  else if(this->npres_nst()==3)
1867  {
1868  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1869  internal_data_pt(this->P_nst_internal_index)->set_value(1,press);
1870  internal_data_pt(this->P_nst_internal_index)->set_value(2,press);
1871  }
1872  else
1873  {
1874  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1875  internal_data_pt(this->P_nst_internal_index)->set_value(1,press);
1876  internal_data_pt(this->P_nst_internal_index)->set_value(2,press);
1877  internal_data_pt(this->P_nst_internal_index)->set_value(3,press);
1878  }
1879  }// Otherwise this is called after p-refinement
1880  }
1881 }
1882 
1883 //======================================================================
1884 /// 2D Further build for Crouzeix_Raviart interpolates the internal
1885 /// pressure dofs from father element: Make sure pressure values and
1886 /// dp/ds agree between fathers and sons at the midpoints of the son
1887 /// elements.
1888 //======================================================================
1889 template<>
1891 {
1892  if (this->tree_pt()->father_pt()!=0)
1893  {
1894  //Call the generic further build (if there is a father)
1896  }
1897  // Now do the PRefineableQElement further_build()
1899 
1900  // Resize internal pressure storage
1901  if (this->internal_data_pt(this->P_nst_internal_index)->nvalue()
1902  <= this->npres_nst())
1903  {
1904  this->internal_data_pt(this->P_nst_internal_index)
1905  ->resize(this->npres_nst());
1906  }
1907  else
1908  {
1909  Data* new_data_pt = new Data(this->npres_nst());
1910  delete internal_data_pt(this->P_nst_internal_index);
1911  internal_data_pt(this->P_nst_internal_index) = new_data_pt;
1912  }
1913 
1914  if(this->tree_pt()->father_pt()!=0)
1915  {
1916  // Pointer to my father (in C-R element impersonation)
1917  PRefineableQCrouzeixRaviartElement<3>* father_element_pt=
1919  (octree_pt()->father_pt()->object_pt());
1920 
1921  // If element has same p-order as father then do the projection problem
1922  // (called after h-refinement)
1923  if(father_element_pt->p_order()==this->p_order())
1924  {
1925  using namespace OcTreeNames;
1926 
1927  // What type of son am I? Ask my quadtree representation...
1928  int son_type=octree_pt()->son_type();
1929 
1930  Vector<double> s_father(3);
1931 
1932 
1933  // Son midpoint is located at the following coordinates in father element:
1934  for(unsigned i=0;i<3;i++)
1935  {
1936  s_father[i]=0.5*OcTree::Direction_to_vector[son_type][i];
1937  }
1938 
1939  // Get pressure value in father element
1940  double press=father_element_pt->interpolated_p_nst(s_father);
1941 
1942  // Reset all pressures to zero
1943  for(unsigned p=0; p<this->npres_nst(); p++)
1944  {
1945  internal_data_pt(this->P_nst_internal_index)->set_value(p,0.0);
1946  }
1947 
1948  // Set pressure values from father (BENFLAG: projection problem hack)
1949  if(this->npres_nst()==1)
1950  {
1951  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1952  }
1953  else
1954  {
1955  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1956  internal_data_pt(this->P_nst_internal_index)->set_value(1,press);
1957  internal_data_pt(this->P_nst_internal_index)->set_value(2,press);
1958  internal_data_pt(this->P_nst_internal_index)->set_value(3,press);
1959  internal_data_pt(this->P_nst_internal_index)->set_value(4,press);
1960  internal_data_pt(this->P_nst_internal_index)->set_value(5,press);
1961  internal_data_pt(this->P_nst_internal_index)->set_value(6,press);
1962  internal_data_pt(this->P_nst_internal_index)->set_value(7,press);
1963  }
1964  }// Otherwise this is called after p-refinement
1965  }
1966 }
1967 
1968 
1969 
1970 //====================================================================
1971 //// Force build of templates
1972 //====================================================================
1973 template class RefineableNavierStokesEquations<2>;
1974 template class RefineableNavierStokesEquations<3>;
1975 template class RefineableQTaylorHoodElement<2>;
1976 template class RefineableQTaylorHoodElement<3>;
1981 }
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
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}.
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
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
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
e
Definition: cfortran.h:575
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
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
A Rank 3 Tensor class.
Definition: matrices.h:1337
void further_build()
Further build, pass the pointers down to the sons.
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
Refineable version of Crouzeix Raviart elements. Generic class definitions.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
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
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...
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...
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
void fill_in_generic_pressure_advection_diffusion_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164