refineable_axisym_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 namespace oomph
34 {
35 
36 
37 //========================================================================
38 /// Add element's contribution to the elemental
39 /// residual vector and/or Jacobian matrix.
40 /// flag=1: compute both
41 /// flag=0: compute only residual vector
42 //=======================================================================
45  DenseMatrix<double> &jacobian,
46  DenseMatrix<double> &mass_matrix,
47  unsigned flag)
48 {
49  //The dimension is actually two
50  unsigned DIM=2;
51 
52 //Find out how many nodes there are
53 unsigned n_node = nnode();
54 
55 // Get continuous time from timestepper of first node
56 double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
57 
58 //Find out how many pressure dofs there are
59 unsigned n_pres = npres_axi_nst();
60 
61 // Get the local indices of the nodal coordinates
62  unsigned u_nodal_index[3];
63  for(unsigned i=0;i<3;++i) {u_nodal_index[i] = u_index_axi_nst(i);}
64 
65 // Which nodal value represents the pressure? (Negative if pressure
66 // is not based on nodal interpolation).
67 int p_index = this->p_nodal_index_axi_nst();
68 
69 // Local array of booleans that are true if the l-th pressure value is
70 // hanging (avoid repeated virtual function calls)
71  bool pressure_dof_is_hanging[n_pres];
72  //If the pressure is stored at a node
73  if(p_index >= 0)
74  {
75  //Read out whether the pressure is hanging
76  for(unsigned l=0;l<n_pres;++l)
77  {
78  pressure_dof_is_hanging[l] =
79  pressure_node_pt(l)->is_hanging(p_index);
80  }
81  }
82  //Otherwise the pressure is not stored at a node and so cannot hang
83  else
84  {
85  for(unsigned l=0;l<n_pres;++l)
86  {pressure_dof_is_hanging[l] = false;}
87  }
88 
89 
90 //Set up memory for the shape and test functions
91 Shape psif(n_node), testf(n_node);
92 DShape dpsifdx(n_node,DIM), dtestfdx(n_node,DIM);
93 
94 
95 //Set up memory for pressure shape and test functions
96 Shape psip(n_pres), testp(n_pres);
97 
98 //Set the value of Nintpt
99 unsigned Nintpt = integral_pt()->nweight();
100 
101 //Set the Vector to hold local coordinates
102 Vector<double> s(DIM);
103 
104 //Get Physical Variables from Element
105 //Reynolds number must be multiplied by the density ratio
106 double scaled_re = re()*density_ratio();
107 double scaled_re_st = re_st()*density_ratio();
108 double scaled_re_inv_fr = re_invfr()*density_ratio();
109 double scaled_re_inv_ro = re_invro()*density_ratio();
110 double visc_ratio = viscosity_ratio(); // hierher -- rewrite and
111  // make consistent with
112  //non-refineable version
113 Vector<double> G = g();
114 
115 //Integers to store the local equation and unknown numbers
116 int local_eqn=0, local_unknown=0;
117 
118 //Local storage for pointers to hang info objects
119  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
120 
121 //Loop over the integration points
122 for(unsigned ipt=0;ipt<Nintpt;ipt++)
123 {
124 
125  //Assign values of s
126  for(unsigned i=0;i<DIM;i++) {s[i] = integral_pt()->knot(ipt,i);}
127 
128  //Get the integral weight
129  double w = integral_pt()->weight(ipt);
130 
131  //Call the derivatives of the shape and test functions
132  double J =
133  dshape_and_dtest_eulerian_at_knot_axi_nst(ipt,psif,dpsifdx,testf,dtestfdx);
134 
135  //Call the pressure shape and test functions
136  pshape_axi_nst(s,psip,testp);
137 
138  //Premultiply the weights and the Jacobian
139  double W = w*J;
140 
141  //Calculate local values of the pressure and velocity components
142  //--------------------------------------------------------------
143  double interpolated_p=0.0;
144  Vector<double> interpolated_u(DIM+1,0.0);
146  Vector<double> mesh_veloc(DIM,0.0);
147  Vector<double> dudt(DIM+1,0.0);
148  DenseMatrix<double> interpolated_dudx(DIM+1,DIM,0.0);
149 
150  //Calculate pressure
151  for(unsigned l=0;l<n_pres;l++) {interpolated_p += p_axi_nst(l)*psip[l];}
152 
153 
154  //Calculate velocities and derivatives
155 
156  // Loop over nodes
157  for(unsigned l=0;l<n_node;l++)
158  {
159  //Cache the shape function
160  const double psif_ = psif(l);
161  //Loop over directions
162  for(unsigned i=0;i<DIM;i++)
163  {
164  interpolated_x[i] += nodal_position(l,i)*psif_;
165  //mesh_veloc[i] +=dnodal_position_dt(l,i)*psif(l);
166  }
167 
168  for(unsigned i=0;i<DIM+1;i++)
169  {
170  const double u_value = nodal_value(l,u_nodal_index[i]);
171  interpolated_u[i] += u_value*psif_;
172  dudt[i]+=du_dt_axi_nst(l,i)*psif_;
173  //Loop over derivative directions for gradients
174  for(unsigned j=0;j<DIM;j++)
175  {
176  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
177  }
178  }
179  }
180 
181  //Get the mesh velocity if ALE is enabled
182  if(!ALE_is_disabled)
183  {
184  // Loop over nodes
185  for(unsigned l=0;l<n_node;l++)
186  {
187  //Loop over the two coordinate directions
188  for(unsigned i=0;i<2;i++)
189  {
190  mesh_veloc[i] += this->dnodal_position_dt(l,i)*psif(l);
191  }
192  }
193  }
194 
195 
196  //Get the user-defined body force terms
197  Vector<double> body_force(DIM+1);
198  get_body_force_axi_nst(time,ipt,s,interpolated_x,body_force);
199 
200  //Get the user-defined source function
201  double source=get_source_fct(time,ipt,interpolated_x);
202 
203  // r is the first postition component
204  double r = interpolated_x[0];
205 
206  //MOMENTUM EQUATIONS
207  //==================
208  //Number of master nodes and storage for the weight of the shape function
209  unsigned n_master=1; double hang_weight=1.0;
210 
211  //Loop over the nodes for the test functions/equations
212  //----------------------------------------------------
213  for(unsigned l=0;l<n_node;l++)
214  {
215  //Local boolean that indicates the hanging status of the node
216  bool is_node_hanging = node_pt(l)->is_hanging();
217 
218  //If the node is hanging
219  if(is_node_hanging)
220  {
221  //Get the hanging pointer
222  hang_info_pt = node_pt(l)->hanging_pt();
223  //Read out number of master nodes from hanging data
224  n_master = hang_info_pt->nmaster();
225  }
226  //Otherwise the node is its own master
227  else
228  {
229  n_master = 1;
230  }
231 
232  //Loop over the master nodes
233  for(unsigned m=0;m<n_master;m++)
234  {
235  // Loop over velocity components for equations
236  for(unsigned i=0;i<DIM+1;i++)
237  {
238  //Get the equation number
239  //If the node is hanging
240  if(is_node_hanging)
241  {
242  //Get the equation number from the master node
243  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
244  u_nodal_index[i]);
245  //Get the hang weight from the master node
246  hang_weight = hang_info_pt->master_weight(m);
247  }
248  //If the node is not hanging
249  else
250  {
251  // Local equation number
252  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
253  // Node contributes with full weight
254  hang_weight = 1.0;
255  }
256 
257  //If it's not a boundary condition...
258  if(local_eqn>= 0)
259  {
260  // initialise for residual calculation
261  double sum=0.0;
262 
263  switch (i)
264  {
265  // RADIAL MOMENTUM EQUATION
266  case 0:
267  //Add the user-defined body force terms
268  sum +=
269  r*body_force[0]*testf[l]*W*hang_weight;
270 
271  //Add the gravitational body force term
272  sum += r*scaled_re_inv_fr*testf[l]*G[0]*W*hang_weight;
273 
274  //Add the pressure gradient term
275  sum +=
276  interpolated_p*(testf[l] + r*dtestfdx(l,0))*W*hang_weight;
277 
278  //Add in the stress tensor terms
279  //The viscosity ratio needs to go in here to ensure
280  //continuity of normal stress is satisfied even in flows
281  //with zero pressure gradient!
282  sum -= visc_ratio*
283  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W
284  *hang_weight;
285 
286  sum -= visc_ratio*r*
287  (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
288  dtestfdx(l,1)*W*hang_weight;
289 
290  sum -=
291  visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]
292  *testf[l]*W*hang_weight/r;
293 
294  //Add in the inertial terms
295  //du/dt term
296  sum -= scaled_re_st*r*dudt[0]*testf[l]*W*hang_weight;
297 
298  //Convective terms
299  sum -=
300  scaled_re*(r*interpolated_u[0]*interpolated_dudx(0,0)
301  - interpolated_u[2]*interpolated_u[2]
302  + r*interpolated_u[1]*interpolated_dudx(0,1))*testf[l]*W*
303  hang_weight;
304 
305  //Mesh velocity terms
306  if(!ALE_is_disabled)
307  {
308  for(unsigned k=0;k<2;k++)
309  {
310  sum +=
311  scaled_re_st*r*mesh_veloc[k]*interpolated_dudx(0,k)*testf[l]*W*
312  hang_weight;
313  }
314  }
315 
316  //Coriolis term
317  sum +=
318  2.0*r*scaled_re_inv_ro*interpolated_u[2]*testf[l]*W*hang_weight;
319 
320  break;
321 
322  // AXIAL MOMENTUM EQUATION
323  case 1:
324  //If it's not a boundary condition
325  //Add the user-defined body force terms
326  //Remember to multiply by the density ratio!
327  sum += r*body_force[1]*testf[l]*W*hang_weight;
328 
329  //Add the gravitational body force term
330  sum += r*scaled_re_inv_fr*testf[l]*G[1]*W*hang_weight;
331 
332  //Add the pressure gradient term
333  sum += r*interpolated_p*dtestfdx(l,1)*W*hang_weight;
334 
335  //Add in the stress tensor terms
336  //The viscosity ratio needs to go in here to ensure
337  //continuity of normal stress is satisfied even in flows
338  //with zero pressure gradient!
339  sum -= visc_ratio*
340  r*(interpolated_dudx(1,0) +
341  Gamma[1]*interpolated_dudx(0,1))*dtestfdx(l,0)*W*
342  hang_weight;
343 
344  sum -= visc_ratio*r*
345  (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W*
346  hang_weight;
347 
348  //Add in the inertial terms
349  //du/dt term
350  sum -= scaled_re_st*r*dudt[1]*testf[l]*W*hang_weight;
351 
352  //Convective terms
353  sum -=
354  scaled_re*(r*interpolated_u[0]*interpolated_dudx(1,0)
355  + r*interpolated_u[1]*interpolated_dudx(1,1))*
356  testf[l]*W*hang_weight;
357 
358  //Mesh velocity terms
359  if(!ALE_is_disabled)
360  {
361  for(unsigned k=0;k<2;k++)
362  {
363  sum +=
364  scaled_re_st*r*mesh_veloc[k]*interpolated_dudx(1,k)*testf[l]*W
365  *hang_weight;
366  }
367  }
368  break;
369 
370  // AZIMUTHAL MOMENTUM EQUATION
371  case 2:
372  //Add the user-defined body force terms
373  //Remember to multiply by the density ratio!
374  sum += r*body_force[2]*testf[l]*W*hang_weight;
375 
376  //Add the gravitational body force term
377  sum += r*scaled_re_inv_fr*testf[l]*G[2]*W*hang_weight;
378 
379  //There is NO pressure gradient term
380 
381  //Add in the stress tensor terms
382  //The viscosity ratio needs to go in here to ensure
383  //continuity of normal stress is satisfied even in flows
384  //with zero pressure gradient!
385  sum -= visc_ratio*
386  (r*interpolated_dudx(2,0) - Gamma[0]*interpolated_u[2])
387  *dtestfdx(l,0)*W*hang_weight;
388 
389  sum -= visc_ratio*r*interpolated_dudx(2,1)*
390  dtestfdx(l,1)*W*hang_weight;
391 
392  sum -= visc_ratio*
393  ((interpolated_u[2]/r) - Gamma[0]*interpolated_dudx(2,0))*
394  testf[l]*W*hang_weight;
395 
396 
397  //Add in the inertial terms
398  //du/dt term
399  sum -= scaled_re_st*r*dudt[2]*testf[l]*W*hang_weight;
400 
401  //Convective terms
402  sum -=
403  scaled_re*(r*interpolated_u[0]*interpolated_dudx(2,0)
404  + interpolated_u[0]*interpolated_u[2]
405  + r*interpolated_u[1]*interpolated_dudx(2,1))*testf[l]*W
406  *hang_weight;
407 
408  //Mesh velocity terms
409  if(!ALE_is_disabled)
410  {
411  for(unsigned k=0;k<2;k++)
412  {
413  sum += scaled_re_st*r*mesh_veloc[k]*
414  interpolated_dudx(2,k)*testf[l]*W*hang_weight;
415  }
416  }
417 
418  //Coriolis term
419  sum -=
420  2.0*r*scaled_re_inv_ro*interpolated_u[0]*testf[l]*W*hang_weight;
421 
422  break;
423  }
424 
425  // Add contribution
426  residuals[local_eqn] += sum;
427 
428  //CALCULATE THE JACOBIAN
429  if (flag)
430  {
431  //Number of master nodes and weights
432  unsigned n_master2=1; double hang_weight2=1.0;
433  //Loop over the velocity nodes for columns
434  for(unsigned l2=0;l2<n_node;l2++)
435  {
436  //Local boolean for hanging status
437  bool is_node2_hanging = node_pt(l2)->is_hanging();
438 
439  //If the node is hanging
440  if(is_node2_hanging)
441  {
442  hang_info2_pt = node_pt(l2)->hanging_pt();
443  //Read out number of master nodes from hanging data
444  n_master2 = hang_info2_pt->nmaster();
445  }
446  //Otherwise the node is its own master
447  else
448  {
449  n_master2 = 1;
450  }
451 
452  //Loop over the master nodes
453  for(unsigned m2=0;m2<n_master2;m2++)
454  {
455  //Loop over the velocity components
456  for(unsigned i2=0;i2<DIM+1;i2++)
457  {
458  //Get the number of the unknown
459  //If the node is hanging
460  if(is_node2_hanging)
461  {
462  //Get the equation number from the master node
463  local_unknown =
464  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
465  u_nodal_index[i2]);
466  hang_weight2 = hang_info2_pt->master_weight(m2);
467  }
468  else
469  {
470  local_unknown = this->nodal_local_eqn(l2,u_nodal_index[i2]);
471  hang_weight2 = 1.0;
472  }
473 
474  // If the unknown is non-zero
475  if(local_unknown >= 0)
476  {
477  //Different results for i and i2
478  switch (i)
479  {
480  // RADIAL MOMENTUM EQUATION
481  case 0:
482  switch (i2)
483  {
484  // radial component
485  case 0:
486 
487  //Add the mass matrix entries
488  if(flag==2)
489  {
490  mass_matrix(local_eqn,local_unknown) +=
491  scaled_re_st*r*psif[l2]*testf[l]*W
492  *hang_weight*hang_weight2;
493  }
494 
495  //Add contribution to the Jacobian matrix
496  jacobian(local_eqn,local_unknown)
497  -= visc_ratio*r*(1.0+Gamma[0])
498  *dpsifdx(l2,0)*dtestfdx(l,0)*W*
499  hang_weight*hang_weight2;
500 
501  jacobian(local_eqn,local_unknown)
502  -= visc_ratio*r*dpsifdx(l2,1)*
503  dtestfdx(l,1)*W*hang_weight*hang_weight2;
504 
505  jacobian(local_eqn,local_unknown)
506  -= visc_ratio*(1.0 + Gamma[0])*psif[l2]*
507  testf[l]*W*hang_weight*hang_weight2/r;
508 
509  //Add in the inertial terms
510  //du/dt term
511  jacobian(local_eqn,local_unknown)
512  -= scaled_re_st*r*node_pt(l2)->time_stepper_pt()->
513  weight(1,0)*psif[l2]*testf[l]*W
514  *hang_weight*hang_weight2;
515 
516  //Convective terms
517  jacobian(local_eqn,local_unknown) -=
518  scaled_re*(r*psif[l2]*interpolated_dudx(0,0)
519  + r*interpolated_u[0]*dpsifdx(l2,0)
520  + r*interpolated_u[1]*dpsifdx(l2,1))
521  *testf[l]*W*hang_weight*hang_weight2;
522 
523  //Mesh velocity terms
524  if(!ALE_is_disabled)
525  {
526  for(unsigned k=0;k<2;k++)
527  {
528  jacobian(local_eqn,local_unknown)
529  += scaled_re_st*r*mesh_veloc[k]*dpsifdx(l2,k)
530  *testf[l]*W*hang_weight*hang_weight2;
531  }
532  }
533  break;
534 
535  // axial component
536  case 1:
537  jacobian(local_eqn,local_unknown) -=
538  visc_ratio*r*Gamma[0]*dpsifdx(l2,0)
539  *dtestfdx(l,1)*W*hang_weight*hang_weight2;
540 
541  //Convective terms
542  jacobian(local_eqn,local_unknown) -=
543  scaled_re*r*psif[l2]*interpolated_dudx(0,1)*
544  testf[l]*W*hang_weight*hang_weight2;
545  break;
546 
547  // azimuthal component
548  case 2:
549  //Convective terms
550  jacobian(local_eqn,local_unknown) -=
551  - scaled_re*2.0*interpolated_u[2]*psif[l2]
552  *testf[l]*W*hang_weight*hang_weight2;
553 
554  //Coriolis terms
555  jacobian(local_eqn,local_unknown) +=
556  2.0*r*scaled_re_inv_ro*psif[l2]*testf[l]*W
557  *hang_weight*hang_weight2;
558 
559  break;
560  } /*End of contribution radial momentum eqn*/
561  break;
562 
563  // AXIAL MOMENTUM EQUATION
564  case 1:
565  switch (i2)
566  {
567  // radial component
568  case 0:
569  //Add in the stress tensor terms
570  //The viscosity ratio needs to go in here to ensure
571  //continuity of normal stress is satisfied even
572  //in flows with zero pressure gradient!
573  jacobian(local_eqn,local_unknown) -=
574  visc_ratio*r*Gamma[1]*dpsifdx(l2,1)*
575  dtestfdx(l,0)*W*hang_weight*hang_weight2;
576 
577  //Convective terms
578  jacobian(local_eqn,local_unknown) -=
579  scaled_re*r*psif[l2]*interpolated_dudx(1,0)*testf[l]*W
580  *hang_weight*hang_weight2;
581  break;
582 
583  // axial component
584  case 1:
585 
586  //Add the mass matrix terms
587  if(flag==2)
588  {
589  mass_matrix(local_eqn,local_unknown) +=
590  scaled_re_st*r*psif[l2]*testf[l]*W
591  *hang_weight*hang_weight2;
592  }
593 
594  //Add in the stress tensor terms
595  //The viscosity ratio needs to go in here to ensure
596  //continuity of normal stress is satisfied even in
597  //flows with zero pressure gradient!
598  jacobian(local_eqn,local_unknown) -=
599  visc_ratio*r*dpsifdx(l2,0)*dtestfdx(l,0)*W
600  *hang_weight*hang_weight2;
601 
602  jacobian(local_eqn,local_unknown) -=
603  visc_ratio*r*(1.0 + Gamma[1])*dpsifdx(l2,1)*
604  dtestfdx(l,1)*W*hang_weight*hang_weight2;
605 
606  //Add in the inertial terms
607  //du/dt term
608  jacobian(local_eqn,local_unknown) -=
609  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->
610  weight(1,0)*psif[l2]*testf[l]*W
611  *hang_weight*hang_weight2;
612 
613  //Convective terms
614  jacobian(local_eqn,local_unknown) -=
615  scaled_re*(r*interpolated_u[0]*dpsifdx(l2,0)
616  + r*psif[l2]*interpolated_dudx(1,1)
617  + r*interpolated_u[1]*dpsifdx(l2,1))
618  *testf[l]*W*hang_weight*hang_weight2;
619 
620  //Mesh velocity terms
621  if(!ALE_is_disabled)
622  {
623  for(unsigned k=0;k<2;k++)
624  {
625  jacobian(local_eqn,local_unknown)
626  += scaled_re_st*r*mesh_veloc[k]*dpsifdx(l2,k)
627  *testf[l]*W*hang_weight*hang_weight2;
628  }
629  }
630  break;
631 
632  // azimuthal component
633  case 2:
634  //There are no azimithal terms in the axial
635  //momentum equation
636  break;
637  }
638  break;
639 
640  // AZIMUTHAL MOMENTUM EQUATION
641  case 2:
642  switch (i2)
643  {
644  // radial component
645  case 0:
646  //Convective terms
647  jacobian(local_eqn,local_unknown) -=
648  scaled_re*(r*psif[l2]*interpolated_dudx(2,0)
649  + psif[l2]*interpolated_u[2])*testf[l]
650  *W*hang_weight*hang_weight2;
651 
652  //Coriolis term
653  jacobian(local_eqn,local_unknown) -=
654  2.0*r*scaled_re_inv_ro*psif[l2]*testf[l]*W
655  *hang_weight*hang_weight2;
656 
657  break;
658 
659  // axial component
660  case 1:
661  //Convective terms
662  jacobian(local_eqn,local_unknown) -=
663  scaled_re*r*psif[l2]*interpolated_dudx(2,1)*testf[l]
664  *W*hang_weight*hang_weight2;
665  break;
666 
667  // azimuthal component
668  case 2:
669 
670  //Add the mass matrix terms
671  if(flag==2)
672  {
673  mass_matrix(local_eqn,local_unknown) +=
674  scaled_re_st*r*psif[l2]*testf[l]*W
675  *hang_weight*hang_weight2;
676  }
677 
678  //Add in the stress tensor terms
679  //The viscosity ratio needs to go in here to ensure
680  //continuity of normal stress is satisfied even in
681  //flows with zero pressure gradient!
682  jacobian(local_eqn,local_unknown) -=
683  visc_ratio*
684  (r*dpsifdx(l2,0) - Gamma[0]*psif[l2])
685  *dtestfdx(l,0)*W*hang_weight*hang_weight2;
686 
687  jacobian(local_eqn,local_unknown) -=
688  visc_ratio*r*dpsifdx(l2,1)*dtestfdx(l,1)
689  *W*hang_weight*hang_weight2;
690 
691  jacobian(local_eqn,local_unknown) -=
692  visc_ratio*
693  ((psif[l2]/r) - Gamma[0]*dpsifdx(l2,0))
694  *testf[l]*W*hang_weight*hang_weight2;
695 
696  //Add in the inertial terms
697  //du/dt term
698  jacobian(local_eqn,local_unknown) -=
699  scaled_re_st*r*node_pt(l2)
700  ->time_stepper_pt()->weight(1,0)*
701  psif[l2]*testf[l]*W*hang_weight*hang_weight2;
702 
703  //Convective terms
704  jacobian(local_eqn,local_unknown) -=
705  scaled_re*(r*interpolated_u[0]*dpsifdx(l2,0)
706  + interpolated_u[0]*psif[l2]
707  + r*interpolated_u[1]*dpsifdx(l2,1))
708  *testf[l]*W*hang_weight*hang_weight2;
709 
710  //Mesh velocity terms
711  if(!ALE_is_disabled)
712  {
713  for(unsigned k=0;k<2;k++)
714  {
715  jacobian(local_eqn,local_unknown)
716  += scaled_re_st*r*mesh_veloc[k]*dpsifdx(l2,k)
717  *testf[l]*W*hang_weight*hang_weight2;
718  }
719  }
720  break;
721  }
722  break;
723  }
724  }
725  }
726  }
727  } //End of loop over the nodes
728 
729 
730  //Loop over the pressure shape functions
731  for(unsigned l2=0;l2<n_pres;l2++)
732  {
733  //If the pressure dof is hanging
734  if(pressure_dof_is_hanging[l2])
735  {
736  // Pressure dof is hanging so it must be nodal-based
737  hang_info2_pt =
738  pressure_node_pt(l2)->hanging_pt(p_index);
739 
740  //Get the number of master nodes from the pressure node
741  n_master2 = hang_info2_pt->nmaster();
742  }
743  //Otherwise the node is its own master
744  else
745  {
746  n_master2 = 1;
747  }
748 
749  //Loop over the master nodes
750  for(unsigned m2=0;m2<n_master2;m2++)
751  {
752  //Get the number of the unknown
753  //If the pressure dof is hanging
754  if(pressure_dof_is_hanging[l2])
755  {
756  //Get the unknown from the node
757  local_unknown =
758  local_hang_eqn(hang_info2_pt->master_node_pt(m2),
759  p_index);
760  //Get the weight from the hanging object
761  hang_weight2 = hang_info2_pt->master_weight(m2);
762  }
763  else
764  {
765  local_unknown = p_local_eqn(l2);
766  hang_weight2 = 1.0;
767  }
768 
769  //If the unknown is not pinned
770  if(local_unknown >= 0)
771  {
772  //Add in contributions to different equations
773  switch (i)
774  {
775  // RADIAL MOMENTUM EQUATION
776  case 0:
777  jacobian(local_eqn,local_unknown)
778  += psip[l2]*(testf[l] + r*dtestfdx(l,0))*W*
779  hang_weight*hang_weight2;
780  break;
781 
782  // AXIAL MOMENTUM EQUATION
783  case 1:
784  jacobian(local_eqn,local_unknown)
785  += r*psip[l2]*dtestfdx(l,1)*W*hang_weight*hang_weight2;
786  break;
787 
788  // AZIMUTHAL MOMENTUM EQUATION
789  case 2:
790  break;
791  }
792  }
793  }
794  } //End of loop over pressure dofs
795  }// End of Jacobian calculation
796  } //End of if not boundary condition statement
797  } //End of loop over velocity components
798  } //End of loop over master nodes
799  } //End of loop over nodes
800 
801 
802  //CONTINUITY EQUATION
803  //===================
804 
805  //Loop over the pressure shape functions
806  for(unsigned l=0;l<n_pres;l++)
807  {
808  //If the pressure dof is hanging
809  if(pressure_dof_is_hanging[l])
810  {
811  // Pressure dof is hanging so it must be nodal-based
812  hang_info_pt = pressure_node_pt(l)->hanging_pt(p_index);
813  //Get the number of master nodes from the pressure node
814  n_master = hang_info_pt->nmaster();
815  }
816  //Otherwise the node is its own master
817  else
818  {
819  n_master = 1;
820  }
821 
822  //Loop over the master nodes
823  for(unsigned m=0;m<n_master;m++)
824  {
825  //Get the number of the unknown
826  //If the pressure dof is hanging
827  if(pressure_dof_is_hanging[l])
828  {
829  local_eqn = local_hang_eqn(hang_info_pt->master_node_pt(m),p_index);
830  hang_weight = hang_info_pt->master_weight(m);
831  }
832  else
833  {
834  local_eqn = p_local_eqn(l);
835  hang_weight = 1.0;
836  }
837 
838  //If the equation is not pinned
839  if(local_eqn >= 0)
840  {
841  // Source term
842  residuals[local_eqn] -= r*source*testp[l]*W*hang_weight;
843 
844  //Gradient terms
845  residuals[local_eqn] +=
846  (interpolated_u[0] + r*interpolated_dudx(0,0)
847  + r*interpolated_dudx(1,1))*testp[l]*W*hang_weight;
848 
849  //CALCULATE THE JACOBIAN
850  //======================
851  if (flag)
852  {
853  unsigned n_master2=1; double hang_weight2=1.0;
854  //Loop over the velocity nodes for columns
855  for(unsigned l2=0;l2<n_node;l2++)
856  {
857  //Local boolean to indicate whether the node is hanging
858  bool is_node2_hanging = node_pt(l2)->is_hanging();
859 
860  //If the node is hanging
861  if(is_node2_hanging)
862  {
863  hang_info2_pt = node_pt(l2)->hanging_pt();
864  //Read out number of master nodes from hanging data
865  n_master2 = hang_info2_pt->nmaster();
866  }
867  //Otherwise the node is its own master
868  else
869  {
870  n_master2 = 1;
871  }
872 
873  //Loop over the master nodes
874  for(unsigned m2=0;m2<n_master2;m2++)
875  {
876  //Loop over the velocity components
877  for(unsigned i2=0;i2<DIM+1;i2++)
878  {
879  //Get the number of the unknown
880  //If the node is hanging
881  if(is_node2_hanging)
882  {
883  //Get the equation number from the master node
884  local_unknown =
885  local_hang_eqn(hang_info2_pt->master_node_pt(m2),
886  u_nodal_index[i2]);
887  hang_weight2 = hang_info2_pt->master_weight(m2);
888  }
889  else
890  {
891  local_unknown = this->nodal_local_eqn(l2,u_nodal_index[i2]);
892  hang_weight2 = 1.0;
893  }
894 
895  //If the unknown is not pinned
896  if(local_unknown >= 0)
897  {
898  switch (i2)
899  {
900  // radial component
901  case 0:
902  jacobian(local_eqn,local_unknown) +=
903  (psif[l2] + r*dpsifdx(l2,0))*testp[l]*W*
904  hang_weight*hang_weight2;
905  break;
906 
907  // axial component
908  case 1:
909  jacobian(local_eqn,local_unknown) +=
910  r*dpsifdx(l2,1)*testp[l]*W*hang_weight*hang_weight2;
911  break;
912 
913  // azimuthal component
914  case 2:
915  break;
916  }
917  }
918  }
919  }
920  } //End of loop over nodes
921 
922  //NO PRESSURE CONTRIBUTIONS TO CONTINUITY EQUATION
923  } //End of Jacobian calculation
924  }
925  }
926  } // End of loop over pressure nodes
927 
928 } // End of loop over integration points
929 
930 
931 } // End of fill_in_generic_residual_contribution_axi_nst(...)
932 
933 
934 //======================================================================
935 /// Compute derivatives of elemental residual vector with respect
936 /// to nodal coordinates.
937 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
938 /// Overloads the FD-based version in the FiniteElement base class.
939 //======================================================================
942  RankThreeTensor<double>& dresidual_dnodal_coordinates)
943 {
944  // Create an Oomph Lib warning
945  std::string warning_message = "This function has not been tested.\n";
946  std::string function = "RefineableAxisymmetricNavierStokesEquations::\n";
947  function += "get_dresidual_dnodal_coordinates(...)";
948  OomphLibWarning(warning_message,function,OOMPH_EXCEPTION_LOCATION);
949 
950  // Return immediately if there are no dofs
951  if(ndof()==0) { return; }
952 
953  // Determine number of nodes in element
954  const unsigned n_node = nnode();
955 
956  // Get continuous time from timestepper of first node
957  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
958 
959  // Determine number of pressure dofs in element
960  const unsigned n_pres = this->npres_axi_nst();
961 
962  // Find the indices at which the local velocities are stored
963  unsigned u_nodal_index[3];
964  for(unsigned i=0;i<3;i++) { u_nodal_index[i] = this->u_index_axi_nst(i); }
965 
966  // Which nodal value represents the pressure? (Negative if pressure
967  // is not based on nodal interpolation).
968  const int p_index = this->p_nodal_index_axi_nst();
969 
970  // Local array of booleans that are true if the l-th pressure value is
971  // hanging (avoid repeated virtual function calls)
972  bool pressure_dof_is_hanging[n_pres];
973 
974  // If the pressure is stored at a node
975  if(p_index >= 0)
976  {
977  // Read out whether the pressure is hanging
978  for(unsigned l=0;l<n_pres;++l)
979  {
980  pressure_dof_is_hanging[l] =
981  pressure_node_pt(l)->is_hanging(p_index);
982  }
983  }
984  // Otherwise the pressure is not stored at a node and so cannot hang
985  else
986  {
987  for(unsigned l=0;l<n_pres;++l) { pressure_dof_is_hanging[l] = false; }
988  }
989 
990  // Set up memory for the shape and test functions
991  // Note that there are only two dimensions, r and z, in this problem
992  Shape psif(n_node), testf(n_node);
993  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
994 
995  // Set up memory for pressure shape and test functions
996  Shape psip(n_pres), testp(n_pres);
997 
998  // Determine number of shape controlling nodes
999  const unsigned n_shape_controlling_node = nshape_controlling_nodes();
1000 
1001  // Deriatives of shape fct derivatives w.r.t. nodal coords
1002  RankFourTensor<double> d_dpsifdx_dX(2,n_shape_controlling_node,n_node,2);
1003  RankFourTensor<double> d_dtestfdx_dX(2,n_shape_controlling_node,n_node,2);
1004 
1005  // Derivative of Jacobian of mapping w.r.t. to nodal coords
1006  DenseMatrix<double> dJ_dX(2,n_shape_controlling_node);
1007 
1008  // Derivatives of derivative of u w.r.t. nodal coords
1009  // Note that the entry d_dudx_dX(p,q,i,k) contains the derivative w.r.t.
1010  // nodal coordinate X_pq of du_i/dx_k. Since there are three velocity
1011  // components, i can take on the values 0, 1 and 2, while k and p only
1012  // take on the values 0 and 1 since there are only two spatial dimensions.
1013  RankFourTensor<double> d_dudx_dX(2,n_shape_controlling_node,3,2);
1014 
1015  // Derivatives of nodal velocities w.r.t. nodal coords:
1016  // Assumption: Interaction only local via no-slip so
1017  // X_pq only affects U_iq.
1018  // Note that the entry d_U_dX(p,q,i) contains the derivative w.r.t. nodal
1019  // coordinate X_pq of nodal value U_iq. In other words this entry is the
1020  // derivative of the i-th velocity component w.r.t. the p-th spatial
1021  // coordinate, taken at the q-th local node.
1022  RankThreeTensor<double> d_U_dX(2,n_shape_controlling_node,3,0.0);
1023 
1024  // Determine the number of integration points
1025  const unsigned n_intpt = integral_pt()->nweight();
1026 
1027  // Vector to hold local coordinates (two dimensions)
1028  Vector<double> s(2);
1029 
1030  // Get physical variables from element
1031  // (Reynolds number must be multiplied by the density ratio)
1032  const double scaled_re = this->re()*this->density_ratio();
1033  const double scaled_re_st = this->re_st()*this->density_ratio();
1034  const double scaled_re_inv_fr = this->re_invfr()*this->density_ratio();
1035  const double scaled_re_inv_ro = this->re_invro()*this->density_ratio();
1036  const double visc_ratio = this->viscosity_ratio();
1037  const Vector<double> G = this->g();
1038 
1039  // FD step
1041 
1042  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1043  // Assumption: Interaction only local via no-slip so
1044  // X_ij only affects U_ij.
1045  bool element_has_node_with_aux_node_update_fct=false;
1046 
1047  std::map<Node*,unsigned> local_shape_controlling_node_lookup=
1049 
1050  // FD loop over shape-controlling nodes
1051  for(std::map<Node*,unsigned>::iterator it=
1052  local_shape_controlling_node_lookup.begin();
1053  it!=local_shape_controlling_node_lookup.end();
1054  it++)
1055  {
1056  // Get pointer to q-th local shape-controlling node
1057  Node* nod_pt=it->first;
1058 
1059  // Get its number
1060  unsigned q=it->second;
1061 
1062  // Only compute if there's a node-update fct involved
1063  if(nod_pt->has_auxiliary_node_update_fct_pt())
1064  {
1065  element_has_node_with_aux_node_update_fct=true;
1066 
1067  // This functionality has not been tested so issue a warning
1068  // to this effect
1069  std::ostringstream warning_stream;
1070  warning_stream << "\nThe functionality to evaluate the additional"
1071  << "\ncontribution to the deriv of the residual eqn"
1072  << "\nw.r.t. the nodal coordinates which comes about"
1073  << "\nif a node's values are updated using an auxiliary"
1074  << "\nnode update function has NOT been tested for"
1075  << "\nrefineable axisymmetric Navier-Stokes elements."
1076  << "\nUse at your own risk"
1077  << std::endl;
1079  warning_stream.str(),
1080  "RefineableAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates",
1081  OOMPH_EXCEPTION_LOCATION);
1082 
1083  // Current nodal velocity
1084  Vector<double> u_ref(3);
1085  for(unsigned i=0;i<3;i++)
1086  {
1087  u_ref[i]=*(nod_pt->value_pt(u_nodal_index[i]));
1088  }
1089 
1090  // FD
1091  for(unsigned p=0;p<2;p++)
1092  {
1093  // Make backup
1094  double backup=nod_pt->x(p);
1095 
1096  // Do FD step. No node update required as we're
1097  // attacking the coordinate directly...
1098  nod_pt->x(p)+=eps_fd;
1099 
1100  // Do auxiliary node update (to apply no slip)
1102 
1103  // Loop over velocity components
1104  for(unsigned i=0;i<3;i++)
1105  {
1106  // Evaluate
1107  d_U_dX(p,q,i)=(*(nod_pt->value_pt(u_nodal_index[p]))-u_ref[p])/eps_fd;
1108  }
1109 
1110  // Reset
1111  nod_pt->x(p)=backup;
1112 
1113  // Do auxiliary node update (to apply no slip)
1115  }
1116  }
1117  }
1118 
1119  // Integer to store the local equation number
1120  int local_eqn=0;
1121 
1122  // Pointers to hang info object
1123  HangInfo *hang_info_pt=0;
1124 
1125  // Loop over the integration points
1126  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1127  {
1128  // Assign values of s
1129  for(unsigned i=0;i<2;i++) { s[i] = integral_pt()->knot(ipt,i); }
1130 
1131  // Get the integral weight
1132  const double w = integral_pt()->weight(ipt);
1133 
1134  // Call the derivatives of the shape and test functions
1135  const double J = this->dshape_and_dtest_eulerian_at_knot_axi_nst(
1136  ipt,psif,dpsifdx,d_dpsifdx_dX,testf,dtestfdx,d_dtestfdx_dX,dJ_dX);
1137 
1138  // Call the pressure shape and test functions
1139  this->pshape_axi_nst(s,psip,testp);
1140 
1141  // Allocate storage for the position and the derivative of the
1142  // mesh positions w.r.t. time
1144  Vector<double> mesh_velocity(2,0.0);
1145 
1146  // Allocate storage for the pressure, velocity components and their
1147  // time and spatial derivatives
1148  double interpolated_p=0.0;
1149  Vector<double> interpolated_u(3,0.0);
1150  Vector<double> dudt(3,0.0);
1151  DenseMatrix<double> interpolated_dudx(3,2,0.0);
1152 
1153  // Calculate pressure at integration point
1154  for(unsigned l=0;l<n_pres;l++)
1155  {
1156  interpolated_p += this->p_axi_nst(l)*psip[l];
1157  }
1158 
1159  // Calculate velocities and derivatives at integration point
1160  // ---------------------------------------------------------
1161 
1162  // Loop over nodes
1163  for(unsigned l=0;l<n_node;l++)
1164  {
1165  // Cache the shape function
1166  const double psif_ = psif(l);
1167 
1168  // Loop over the two coordinate directions
1169  for(unsigned i=0;i<2;i++)
1170  {
1171  interpolated_x[i] += nodal_position(l,i)*psif_;
1172  }
1173 
1174  // Loop over the three velocity directions
1175  for(unsigned i=0;i<3;i++)
1176  {
1177  // Get the nodal value
1178  const double u_value = nodal_value(l,u_nodal_index[i]);
1179  interpolated_u[i] += u_value*psif_;
1180  dudt[i] += this->du_dt_axi_nst(l,i)*psif_;
1181 
1182  // Loop over derivative directions
1183  for(unsigned j=0;j<2;j++)
1184  {
1185  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
1186  }
1187  }
1188  }
1189 
1190  // Get the mesh velocity if ALE is enabled
1191  if(!this->ALE_is_disabled)
1192  {
1193  // Loop over nodes
1194  for(unsigned l=0;l<n_node;l++)
1195  {
1196  // Loop over the two coordinate directions
1197  for(unsigned i=0;i<2;i++)
1198  {
1199  mesh_velocity[i] += this->dnodal_position_dt(l,i)*psif[l];
1200  }
1201  }
1202  }
1203 
1204  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1205 
1206  // Loop over shape-controlling nodes
1207  for(unsigned q=0;q<n_shape_controlling_node;q++)
1208  {
1209  // Loop over the two coordinate directions
1210  for(unsigned p=0;p<2;p++)
1211  {
1212  // Loop over the three velocity components
1213  for(unsigned i=0;i<3;i++)
1214  {
1215  // Loop over the two coordinate directions
1216  for(unsigned k=0;k<2;k++)
1217  {
1218  double aux = 0.0;
1219 
1220  // Loop over nodes and add contribution
1221  for(unsigned j=0;j<n_node;j++)
1222  {
1223  aux += nodal_value(j,u_nodal_index[i])*d_dpsifdx_dX(p,q,j,k);
1224  }
1225  d_dudx_dX(p,q,i,k) = aux;
1226  }
1227  }
1228  }
1229  }
1230 
1231  // Get weight of actual nodal position/value in computation of mesh
1232  // velocity from positional/value time stepper
1233  const double pos_time_weight
1234  = node_pt(0)->position_time_stepper_pt()->weight(1,0);
1235  const double val_time_weight = node_pt(0)->time_stepper_pt()->weight(1,0);
1236 
1237  // Get the user-defined body force terms
1238  Vector<double> body_force(3);
1239  this->get_body_force_axi_nst(time,ipt,s,interpolated_x,body_force);
1240 
1241  // Get the user-defined source function
1242  const double source = this->get_source_fct(time,ipt,interpolated_x);
1243 
1244  // Get gradient of body force function
1245  DenseMatrix<double> d_body_force_dx(3,2,0.0);
1246  this->get_body_force_gradient_axi_nst(time,ipt,s,interpolated_x,
1247  d_body_force_dx);
1248 
1249  // Get gradient of source function
1250  Vector<double> source_gradient(2,0.0);
1251  this->get_source_fct_gradient(time,ipt,interpolated_x,
1252  source_gradient);
1253 
1254  // Cache r (the first position component)
1255  const double r = interpolated_x[0];
1256 
1257  // Assemble shape derivatives
1258  // --------------------------
1259 
1260  // ==================
1261  // MOMENTUM EQUATIONS
1262  // ==================
1263 
1264  // Number of master nodes and storage for the weight of the shape function
1265  unsigned n_master=1; double hang_weight=1.0;
1266 
1267  // Loop over the test functions
1268  for(unsigned l=0;l<n_node;l++)
1269  {
1270  // Cache the test function
1271  const double testf_ = testf[l];
1272 
1273  // Local boolean to indicate whether the node is hanging
1274  bool is_node_hanging = node_pt(l)->is_hanging();
1275 
1276  // If the node is hanging
1277  if(is_node_hanging)
1278  {
1279  hang_info_pt = node_pt(l)->hanging_pt();
1280 
1281  // Read out number of master nodes from hanging data
1282  n_master = hang_info_pt->nmaster();
1283  }
1284  // Otherwise the node is its own master
1285  else
1286  {
1287  n_master = 1;
1288  }
1289 
1290  // Loop over the master nodes
1291  for(unsigned m=0;m<n_master;m++)
1292  {
1293 
1294  // --------------------------------
1295  // FIRST (RADIAL) MOMENTUM EQUATION
1296  // --------------------------------
1297 
1298  // Get the equation number
1299  // If the node is hanging
1300  if(is_node_hanging)
1301  {
1302  // Get the equation number from the master node
1303  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1304  u_nodal_index[0]);
1305  // Get the hang weight from the master node
1306  hang_weight = hang_info_pt->master_weight(m);
1307  }
1308  // If the node is not hanging
1309  else
1310  {
1311  // Local equation number
1312  local_eqn = this->nodal_local_eqn(l,u_nodal_index[0]);
1313 
1314  // Node contributes with full weight
1315  hang_weight = 1.0;
1316  }
1317 
1318  // If it's not a boundary condition
1319  if(local_eqn >= 0)
1320  {
1321  // Loop over the two coordinate directions
1322  for(unsigned p=0;p<2;p++)
1323  {
1324  // Loop over shape controlling nodes
1325  for(unsigned q=0;q<n_shape_controlling_node;q++)
1326  {
1327 
1328  // Residual x deriv of Jacobian
1329  // ----------------------------
1330 
1331  // Add the user-defined body force terms
1332  double sum = r*body_force[0]*testf_;
1333 
1334  // Add the gravitational body force term
1335  sum += r*scaled_re_inv_fr*testf_*G[0];
1336 
1337  // Add the pressure gradient term
1338  sum += interpolated_p*(testf_ + r*dtestfdx(l,0));
1339 
1340  // Add the stress tensor terms
1341  // The viscosity ratio needs to go in here to ensure
1342  // continuity of normal stress is satisfied even in flows
1343  // with zero pressure gradient!
1344  sum -= visc_ratio*
1345  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0);
1346 
1347  sum -= visc_ratio*r*
1348  (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
1349  dtestfdx(l,1);
1350 
1351  sum -= visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]*testf_/r;
1352 
1353  // Add the du/dt term
1354  sum -= scaled_re_st*r*dudt[0]*testf_;
1355 
1356  // Add the convective terms
1357  sum -=
1358  scaled_re*(r*interpolated_u[0]*interpolated_dudx(0,0)
1359  - interpolated_u[2]*interpolated_u[2]
1360  + r*interpolated_u[1]*interpolated_dudx(0,1))*testf_;
1361 
1362  // Add the mesh velocity terms
1363  if(!ALE_is_disabled)
1364  {
1365  for(unsigned k=0;k<2;k++)
1366  {
1367  sum += scaled_re_st*r*mesh_velocity[k]
1368  *interpolated_dudx(0,k)*testf_;
1369  }
1370  }
1371 
1372  // Add the Coriolis term
1373  sum += 2.0*r*scaled_re_inv_ro*interpolated_u[2]*testf_;
1374 
1375  // Multiply through by deriv of Jacobian and integration weight
1376  dresidual_dnodal_coordinates(local_eqn,p,q)
1377  += sum*dJ_dX(p,q)*w*hang_weight;
1378 
1379  // Derivative of residual x Jacobian
1380  // ---------------------------------
1381 
1382  // Body force
1383  sum = r*d_body_force_dx(0,p)*psif[q]*testf_;
1384  if(p==0) { sum += body_force[0]*psif[q]*testf_; }
1385 
1386  // Gravitational body force
1387  if(p==0) { sum += scaled_re_inv_fr*G[0]*psif[q]*testf_; }
1388 
1389  // Pressure gradient term
1390  sum += r*interpolated_p*d_dtestfdx_dX(p,q,l,0);
1391  if(p==0) { sum += interpolated_p*psif[q]*dtestfdx(l,0); }
1392 
1393  // Viscous terms
1394  sum -= r*visc_ratio*(
1395  (1.0+Gamma[0])*(
1396  d_dudx_dX(p,q,0,0)*dtestfdx(l,0)
1397  + interpolated_dudx(0,0)*d_dtestfdx_dX(p,q,l,0))
1398  + (d_dudx_dX(p,q,0,1)
1399  + Gamma[0]*d_dudx_dX(p,q,1,0))*dtestfdx(l,1)
1400  + (interpolated_dudx(0,1)
1401  + Gamma[0]*interpolated_dudx(1,0))
1402  *d_dtestfdx_dX(p,q,l,1));
1403  if(p==0)
1404  {
1405  sum -= visc_ratio*(
1406  (1.0+Gamma[0])*(
1407  interpolated_dudx(0,0)*psif[q]*dtestfdx(l,0)
1408  - interpolated_u[0]*psif[q]*testf_/(r*r))
1409  + (interpolated_dudx(0,1)
1410  + Gamma[0]*interpolated_dudx(1,0))*psif[q]*dtestfdx(l,1));
1411  }
1412 
1413  // Convective terms, including mesh velocity
1414  for(unsigned k=0;k<2;k++)
1415  {
1416  double tmp = scaled_re*interpolated_u[k];
1417  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
1418  sum -= r*tmp*d_dudx_dX(p,q,0,k)*testf_;
1419  if(p==0) { sum -= tmp*interpolated_dudx(0,k)*psif[q]*testf_; }
1420  }
1421  if(!ALE_is_disabled)
1422  {
1423  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(0,p)
1424  *psif[q]*testf_;
1425  }
1426 
1427  // du/dt term
1428  if(p==0) { sum -= scaled_re_st*dudt[0]*psif[q]*testf_; }
1429 
1430  // Coriolis term
1431  if(p==0)
1432  {
1433  sum += 2.0*scaled_re_inv_ro*interpolated_u[2]*psif[q]*testf_;
1434  }
1435 
1436  // Multiply through by Jacobian and integration weight
1437  dresidual_dnodal_coordinates(local_eqn,p,q)
1438  += sum*J*w*hang_weight;
1439 
1440  } // End of loop over shape controlling nodes q
1441  } // End of loop over coordinate directions p
1442 
1443  // Derivs w.r.t. to nodal velocities
1444  // ---------------------------------
1445  if(element_has_node_with_aux_node_update_fct)
1446  {
1447  // Loop over local nodes
1448  for (unsigned q_local=0;q_local<n_node;q_local++)
1449  {
1450  // Number of master nodes and storage for the weight of
1451  // the shape function
1452  unsigned n_master2=1;
1453  double hang_weight2=1.0;
1454  HangInfo* hang_info2_pt=0;
1455 
1456  // Local boolean to indicate whether the node is hanging
1457  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1458 
1459  Node* actual_shape_controlling_node_pt=node_pt(q_local);
1460 
1461  // If the node is hanging
1462  if(is_node_hanging2)
1463  {
1464  hang_info2_pt = node_pt(q_local)->hanging_pt();
1465 
1466  // Read out number of master nodes from hanging data
1467  n_master2 = hang_info2_pt->nmaster();
1468  }
1469  // Otherwise the node is its own master
1470  else
1471  {
1472  n_master2 = 1;
1473  }
1474 
1475  // Loop over the master nodes
1476  for(unsigned mm=0;mm<n_master2;mm++)
1477  {
1478  if(is_node_hanging2)
1479  {
1480  actual_shape_controlling_node_pt=
1481  hang_info2_pt->master_node_pt(mm);
1482  hang_weight2 = hang_info2_pt->master_weight(mm);
1483  }
1484 
1485  // Find the corresponding number
1486  unsigned q=local_shape_controlling_node_lookup[
1487  actual_shape_controlling_node_pt];
1488 
1489  // Loop over the two coordinate directions
1490  for(unsigned p=0;p<2;p++)
1491  {
1492  // Contribution from deriv of first velocity component
1493  double tmp = scaled_re_st*r*val_time_weight
1494  *psif[q_local]*testf_;
1495  tmp += r*scaled_re*interpolated_dudx(0,0)
1496  *psif[q_local]*testf_;
1497 
1498  for(unsigned k=0;k<2;k++)
1499  {
1500  double tmp2 = scaled_re*interpolated_u[k];
1501  if(!ALE_is_disabled)
1502  {
1503  tmp2 -= scaled_re_st*mesh_velocity[k];
1504  }
1505  tmp += r*tmp2*dpsifdx(q_local,k)*testf_;
1506  }
1507 
1508  tmp += r*visc_ratio*(1.0+Gamma[0])*dpsifdx(q_local,0)
1509  *dtestfdx(l,0);
1510  tmp += r*visc_ratio*dpsifdx(q_local,1)*dtestfdx(l,1);
1511  tmp += visc_ratio*(1.0+Gamma[0])*psif[q_local]*testf_/r;
1512 
1513  // Multiply through by dU_0q/dX_pq
1514  double sum = -d_U_dX(p,q_local,0)*tmp;
1515 
1516  // Contribution from deriv of second velocity component
1517  tmp = scaled_re*r*interpolated_dudx(0,1)*psif[q_local]*testf_;
1518  tmp += r*visc_ratio*Gamma[0]*dpsifdx(q_local,0)*dtestfdx(l,1);
1519 
1520  // Multiply through by dU_1q/dX_pq
1521  sum -= d_U_dX(p,q,1)*tmp;
1522 
1523  // Contribution from deriv of third velocity component
1524  tmp = 2.0*(r*scaled_re_inv_ro + scaled_re*interpolated_u[2])
1525  *psif[q_local]*testf_;
1526 
1527  // Multiply through by dU_2q/dX_pq
1528  sum += d_U_dX(p,q,2)*tmp;
1529 
1530  // Multiply through by Jacobian and integration weight
1531  dresidual_dnodal_coordinates(local_eqn,p,q) +=
1532  sum*J*w*hang_weight*hang_weight2;
1533  }
1534  } // End of loop over master nodes
1535  } // End of loop over local nodes
1536  } // End of if(element_has_node_with_aux_node_update_fct)
1537  } // End of if it's not a boundary condition
1538 
1539  // --------------------------------
1540  // SECOND (AXIAL) MOMENTUM EQUATION
1541  // --------------------------------
1542 
1543  // Get the equation number
1544  // If the node is hanging
1545  if(is_node_hanging)
1546  {
1547  // Get the equation number from the master node
1548  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1549  u_nodal_index[1]);
1550  // Get the hang weight from the master node
1551  hang_weight = hang_info_pt->master_weight(m);
1552  }
1553  // If the node is not hanging
1554  else
1555  {
1556  // Local equation number
1557  local_eqn = this->nodal_local_eqn(l,u_nodal_index[1]);
1558 
1559  // Node contributes with full weight
1560  hang_weight = 1.0;
1561  }
1562 
1563  // If it's not a boundary condition
1564  if(local_eqn >= 0)
1565  {
1566  // Loop over the two coordinate directions
1567  for(unsigned p=0;p<2;p++)
1568  {
1569  // Loop over shape controlling nodes
1570  for(unsigned q=0;q<n_shape_controlling_node;q++)
1571  {
1572 
1573  // Residual x deriv of Jacobian
1574  // ----------------------------
1575 
1576  // Add the user-defined body force terms
1577  double sum = r*body_force[1]*testf_;
1578 
1579  // Add the gravitational body force term
1580  sum += r*scaled_re_inv_fr*testf_*G[1];
1581 
1582  // Add the pressure gradient term
1583  sum += r*interpolated_p*dtestfdx(l,1);
1584 
1585  // Add the stress tensor terms
1586  // The viscosity ratio needs to go in here to ensure
1587  // continuity of normal stress is satisfied even in flows
1588  // with zero pressure gradient!
1589  sum -= visc_ratio*
1590  r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
1591  *dtestfdx(l,0);
1592 
1593  sum -= visc_ratio*r*
1594  (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1);
1595 
1596  // Add the du/dt term
1597  sum -= scaled_re_st*r*dudt[1]*testf_;
1598 
1599  // Add the convective terms
1600  sum -=
1601  scaled_re*(r*interpolated_u[0]*interpolated_dudx(1,0)
1602  + r*interpolated_u[1]*interpolated_dudx(1,1))*testf_;
1603 
1604  // Add the mesh velocity terms
1605  if(!ALE_is_disabled)
1606  {
1607  for(unsigned k=0;k<2;k++)
1608  {
1609  sum += scaled_re_st*r*mesh_velocity[k]
1610  *interpolated_dudx(1,k)*testf_;
1611  }
1612  }
1613 
1614  // Multiply through by deriv of Jacobian and integration weight
1615  dresidual_dnodal_coordinates(local_eqn,p,q)
1616  += sum*dJ_dX(p,q)*w*hang_weight;
1617 
1618  // Derivative of residual x Jacobian
1619  // ---------------------------------
1620 
1621  // Body force
1622  sum = r*d_body_force_dx(1,p)*psif[q]*testf_;
1623  if(p==0) { sum += body_force[1]*psif[q]*testf_; }
1624 
1625  // Gravitational body force
1626  if(p==0) { sum += scaled_re_inv_fr*G[1]*psif[q]*testf_; }
1627 
1628  // Pressure gradient term
1629  sum += r*interpolated_p*d_dtestfdx_dX(p,q,l,1);
1630  if(p==0) { sum += interpolated_p*psif[q]*dtestfdx(l,1); }
1631 
1632  // Viscous terms
1633  sum -= r*visc_ratio*(
1634  (d_dudx_dX(p,q,1,0) + Gamma[1]*d_dudx_dX(p,q,0,1))*dtestfdx(l,0)
1635  + (interpolated_dudx(1,0)
1636  + Gamma[1]*interpolated_dudx(0,1))*d_dtestfdx_dX(p,q,l,0)
1637  + (1.0+Gamma[1])*d_dudx_dX(p,q,1,1)*dtestfdx(l,1)
1638  + (1.0+Gamma[1])*interpolated_dudx(1,1)*d_dtestfdx_dX(p,q,l,1));
1639  if(p==0)
1640  {
1641  sum -= visc_ratio*(
1642  (interpolated_dudx(1,0)
1643  + Gamma[1]*interpolated_dudx(0,1))*psif[q]*dtestfdx(l,0)
1644  + (1.0+Gamma[1])*interpolated_dudx(1,1)*psif[q]*dtestfdx(l,1));
1645  }
1646 
1647  // Convective terms, including mesh velocity
1648  for(unsigned k=0;k<2;k++)
1649  {
1650  double tmp = scaled_re*interpolated_u[k];
1651  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
1652  sum -= r*tmp*d_dudx_dX(p,q,1,k)*testf_;
1653  if(p==0) { sum -= tmp*interpolated_dudx(1,k)*psif[q]*testf_; }
1654  }
1655  if(!ALE_is_disabled)
1656  {
1657  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(1,p)
1658  *psif[q]*testf_;
1659  }
1660 
1661  // du/dt term
1662  if(p==0) { sum -= scaled_re_st*dudt[1]*psif[q]*testf_; }
1663 
1664  // Multiply through by Jacobian and integration weight
1665  dresidual_dnodal_coordinates(local_eqn,p,q)
1666  += sum*J*w*hang_weight;
1667 
1668  } // End of loop over shape controlling nodes q
1669  } // End of loop over coordinate directions p
1670 
1671  // Derivs w.r.t. to nodal velocities
1672  // ---------------------------------
1673  if(element_has_node_with_aux_node_update_fct)
1674  {
1675  // Loop over local nodes
1676  for (unsigned q_local=0;q_local<n_node;q_local++)
1677  {
1678  // Number of master nodes and storage for the weight of
1679  // the shape function
1680  unsigned n_master2=1;
1681  double hang_weight2=1.0;
1682  HangInfo* hang_info2_pt=0;
1683 
1684  // Local boolean to indicate whether the node is hanging
1685  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1686 
1687  Node* actual_shape_controlling_node_pt=node_pt(q_local);
1688 
1689  // If the node is hanging
1690  if(is_node_hanging2)
1691  {
1692  hang_info2_pt = node_pt(q_local)->hanging_pt();
1693 
1694  // Read out number of master nodes from hanging data
1695  n_master2 = hang_info2_pt->nmaster();
1696  }
1697  // Otherwise the node is its own master
1698  else
1699  {
1700  n_master2 = 1;
1701  }
1702 
1703  // Loop over the master nodes
1704  for(unsigned mm=0;mm<n_master2;mm++)
1705  {
1706  if(is_node_hanging2)
1707  {
1708  actual_shape_controlling_node_pt=
1709  hang_info2_pt->master_node_pt(mm);
1710  hang_weight2 = hang_info2_pt->master_weight(mm);
1711  }
1712 
1713  // Find the corresponding number
1714  unsigned q=local_shape_controlling_node_lookup[
1715  actual_shape_controlling_node_pt];
1716 
1717  // Loop over the two coordinate directions
1718  for(unsigned p=0;p<2;p++)
1719  {
1720  // Contribution from deriv of first velocity component
1721  double tmp = scaled_re*r*interpolated_dudx(1,0)
1722  *psif[q_local]*testf_;
1723  tmp += r*visc_ratio*Gamma[1]*dpsifdx(q_local,1)*dtestfdx(l,0);
1724 
1725  // Multiply through by dU_0q/dX_pq
1726  double sum = -d_U_dX(p,q,0)*tmp;
1727 
1728  // Contribution from deriv of second velocity component
1729  tmp = scaled_re_st*r*val_time_weight*psif[q_local]*testf_;
1730  tmp += r*scaled_re*interpolated_dudx(1,1)
1731  *psif[q_local]*testf_;
1732 
1733  for(unsigned k=0;k<2;k++)
1734  {
1735  double tmp2 = scaled_re*interpolated_u[k];
1736  if(!ALE_is_disabled)
1737  {
1738  tmp2 -= scaled_re_st*mesh_velocity[k];
1739  }
1740  tmp += r*tmp2*dpsifdx(q_local,k)*testf_;
1741  }
1742  tmp += r*visc_ratio*(dpsifdx(q_local,0)*dtestfdx(l,0)
1743  + (1.0+Gamma[1])*dpsifdx(q_local,1)
1744  *dtestfdx(l,1));
1745 
1746  // Multiply through by dU_1q/dX_pq
1747  sum -= d_U_dX(p,q,1)*tmp;
1748 
1749  // Multiply through by Jacobian and integration weight
1750  dresidual_dnodal_coordinates(local_eqn,p,q) +=
1751  sum*J*w*hang_weight*hang_weight2;
1752  }
1753  } // End of loop over master nodes
1754  } // End of loop over local nodes
1755  } // End of if(element_has_node_with_aux_node_update_fct)
1756  } // End of if it's not a boundary condition
1757 
1758  // -----------------------------------
1759  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
1760  // -----------------------------------
1761 
1762  // Get the equation number
1763  // If the node is hanging
1764  if(is_node_hanging)
1765  {
1766  // Get the equation number from the master node
1767  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1768  u_nodal_index[2]);
1769  // Get the hang weight from the master node
1770  hang_weight = hang_info_pt->master_weight(m);
1771  }
1772  // If the node is not hanging
1773  else
1774  {
1775  // Local equation number
1776  local_eqn = this->nodal_local_eqn(l,u_nodal_index[2]);
1777 
1778  // Node contributes with full weight
1779  hang_weight = 1.0;
1780  }
1781 
1782  // If it's not a boundary condition
1783  if(local_eqn >= 0)
1784  {
1785  // Loop over the two coordinate directions
1786  for(unsigned p=0;p<2;p++)
1787  {
1788  // Loop over shape controlling nodes
1789  for(unsigned q=0;q<n_shape_controlling_node;q++)
1790  {
1791 
1792  // Residual x deriv of Jacobian
1793  // ----------------------------
1794 
1795  // Add the user-defined body force terms
1796  double sum = r*body_force[2]*testf_;
1797 
1798  // Add the gravitational body force term
1799  sum += r*scaled_re_inv_fr*testf_*G[2];
1800 
1801  // There is NO pressure gradient term
1802 
1803  // Add in the stress tensor terms
1804  // The viscosity ratio needs to go in here to ensure
1805  // continuity of normal stress is satisfied even in flows
1806  // with zero pressure gradient!
1807  sum -= visc_ratio*(r*interpolated_dudx(2,0) -
1808  Gamma[0]*interpolated_u[2])*dtestfdx(l,0);
1809 
1810  sum -= visc_ratio*r*interpolated_dudx(2,1)*dtestfdx(l,1);
1811 
1812  sum -= visc_ratio*((interpolated_u[2]/r)
1813  - Gamma[0]*interpolated_dudx(2,0))*testf_;
1814 
1815  // Add the du/dt term
1816  sum -= scaled_re_st*r*dudt[2]*testf_;
1817 
1818  // Add the convective terms
1819  sum -=
1820  scaled_re*(r*interpolated_u[0]*interpolated_dudx(2,0)
1821  + interpolated_u[0]*interpolated_u[2]
1822  + r*interpolated_u[1]*interpolated_dudx(2,1))*testf_;
1823 
1824  // Add the mesh velocity terms
1825  if(!ALE_is_disabled)
1826  {
1827  for(unsigned k=0;k<2;k++)
1828  {
1829  sum += scaled_re_st*r*mesh_velocity[k]
1830  *interpolated_dudx(2,k)*testf_;
1831  }
1832  }
1833 
1834  // Add the Coriolis term
1835  sum -= 2.0*r*scaled_re_inv_ro*interpolated_u[0]*testf_;
1836 
1837  // Multiply through by deriv of Jacobian and integration weight
1838  dresidual_dnodal_coordinates(local_eqn,p,q)
1839  += sum*dJ_dX(p,q)*w*hang_weight;
1840 
1841  // Derivative of residual x Jacobian
1842  // ---------------------------------
1843 
1844  // Body force
1845  sum = r*d_body_force_dx(2,p)*psif[q]*testf_;
1846  if(p==0) { sum += body_force[2]*psif[q]*testf_; }
1847 
1848  // Gravitational body force
1849  if(p==0) { sum += scaled_re_inv_fr*G[2]*psif[q]*testf_; }
1850 
1851  // There is NO pressure gradient term
1852 
1853  // Viscous terms
1854  sum -= visc_ratio*r*(
1855  d_dudx_dX(p,q,2,0)*dtestfdx(l,0)
1856  + interpolated_dudx(2,0)*d_dtestfdx_dX(p,q,l,0)
1857  + d_dudx_dX(p,q,2,1)*dtestfdx(l,1)
1858  + interpolated_dudx(2,1)*d_dtestfdx_dX(p,q,l,1));
1859 
1860  sum += visc_ratio*Gamma[0]*d_dudx_dX(p,q,2,0)*testf_;
1861 
1862  if(p==0)
1863  {
1864  sum -=
1865  visc_ratio*(interpolated_dudx(2,0)*psif[q]*dtestfdx(l,0)
1866  + interpolated_dudx(2,1)*psif[q]*dtestfdx(l,1)
1867  + interpolated_u[2]*psif[q]*testf_/(r*r));
1868  }
1869 
1870  // Convective terms, including mesh velocity
1871  for(unsigned k=0;k<2;k++)
1872  {
1873  double tmp = scaled_re*interpolated_u[k];
1874  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
1875  sum -= r*tmp*d_dudx_dX(p,q,2,k)*testf_;
1876  if(p==0) { sum -= tmp*interpolated_dudx(2,k)*psif[q]*testf_; }
1877  }
1878  if(!ALE_is_disabled)
1879  {
1880  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(2,p)
1881  *psif[q]*testf_;
1882  }
1883 
1884  // du/dt term
1885  if(p==0) { sum -= scaled_re_st*dudt[2]*psif[q]*testf_; }
1886 
1887  // Coriolis term
1888  if(p==0)
1889  {
1890  sum -= 2.0*scaled_re_inv_ro*interpolated_u[0]*psif[q]*testf_;
1891  }
1892 
1893  // Multiply through by Jacobian and integration weight
1894  dresidual_dnodal_coordinates(local_eqn,p,q)
1895  += sum*J*w*hang_weight;
1896 
1897  } // End of loop over shape controlling nodes q
1898  } // End of loop over coordinate directions p
1899 
1900  // Derivs w.r.t. to nodal velocities
1901  // ---------------------------------
1902  if(element_has_node_with_aux_node_update_fct)
1903  {
1904  // Loop over local nodes
1905  for (unsigned q_local=0;q_local<n_node;q_local++)
1906  {
1907  // Number of master nodes and storage for the weight of
1908  // the shape function
1909  unsigned n_master2=1;
1910  double hang_weight2=1.0;
1911  HangInfo* hang_info2_pt=0;
1912 
1913  // Local boolean to indicate whether the node is hanging
1914  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1915 
1916  Node* actual_shape_controlling_node_pt=node_pt(q_local);
1917 
1918  // If the node is hanging
1919  if(is_node_hanging2)
1920  {
1921  hang_info2_pt = node_pt(q_local)->hanging_pt();
1922 
1923  // Read out number of master nodes from hanging data
1924  n_master2 = hang_info2_pt->nmaster();
1925  }
1926  // Otherwise the node is its own master
1927  else
1928  {
1929  n_master2 = 1;
1930  }
1931 
1932  // Loop over the master nodes
1933  for(unsigned mm=0;mm<n_master2;mm++)
1934  {
1935  if(is_node_hanging2)
1936  {
1937  actual_shape_controlling_node_pt=
1938  hang_info2_pt->master_node_pt(mm);
1939  hang_weight2 = hang_info2_pt->master_weight(mm);
1940  }
1941 
1942  // Find the corresponding number
1943  unsigned q=local_shape_controlling_node_lookup[
1944  actual_shape_controlling_node_pt];
1945 
1946  // Loop over the two coordinate directions
1947  for(unsigned p=0;p<2;p++)
1948  {
1949  // Contribution from deriv of first velocity component
1950  double tmp = (2.0*r*scaled_re_inv_ro
1951  + r*scaled_re*interpolated_dudx(2,0)
1952  - scaled_re*interpolated_u[2])
1953  *psif[q_local]*testf_;
1954 
1955  // Multiply through by dU_0q/dX_pq
1956  double sum = -d_U_dX(p,q,0)*tmp;
1957 
1958  // Contribution from deriv of second velocity component
1959  // Multiply through by dU_1q/dX_pq
1960  sum -= d_U_dX(p,q,1)*scaled_re*r*interpolated_dudx(2,1)
1961  *psif[q_local]*testf_;
1962 
1963  // Contribution from deriv of third velocity component
1964  tmp = scaled_re_st*r*val_time_weight*psif[q_local]*testf_;
1965  tmp -= scaled_re*interpolated_u[0]*psif[q_local]*testf_;
1966 
1967  for(unsigned k=0;k<2;k++)
1968  {
1969  double tmp2 = scaled_re*interpolated_u[k];
1970  if(!ALE_is_disabled)
1971  {
1972  tmp2 -= scaled_re_st*mesh_velocity[k];
1973  }
1974  tmp += r*tmp2*dpsifdx(q_local,k)*testf_;
1975  }
1976  tmp += visc_ratio*(r*dpsifdx(q_local,0)
1977  - Gamma[0]*psif[q_local])*dtestfdx(l,0);
1978  tmp += r*visc_ratio*dpsifdx(q_local,1)*dtestfdx(l,1);
1979  tmp += visc_ratio*((psif[q_local]/r)
1980  - Gamma[0]*dpsifdx(q_local,0))*testf_;
1981 
1982  // Multiply through by dU_2q/dX_pq
1983  sum -= d_U_dX(p,q,2)*tmp;
1984 
1985  // Multiply through by Jacobian and integration weight
1986  dresidual_dnodal_coordinates(local_eqn,p,q) +=
1987  sum*J*w*hang_weight*hang_weight2;
1988  }
1989  } // End of loop over master nodes
1990  } // End of loop over local nodes
1991  } // End of if(element_has_node_with_aux_node_update_fct)
1992  } // End of if it's not a boundary condition
1993  } // End of loop over master nodes
1994  } // End of loop over test functions
1995 
1996 
1997  // ===================
1998  // CONTINUITY EQUATION
1999  // ===================
2000 
2001  // Loop over the shape functions
2002  for(unsigned l=0;l<n_pres;l++)
2003  {
2004 
2005  // If the pressure dof is hanging
2006  if(pressure_dof_is_hanging[l])
2007  {
2008  // Pressure dof is hanging so it must be nodal-based
2009  // Get the hang info object
2010  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
2011 
2012  // Get the number of master nodes from the pressure node
2013  n_master = hang_info_pt->nmaster();
2014  }
2015  // Otherwise the node is its own master
2016  else
2017  {
2018  n_master = 1;
2019  }
2020 
2021  // Loop over the master nodes
2022  for(unsigned m=0;m<n_master;m++)
2023  {
2024  // Get the number of the unknown
2025  // If the pressure dof is hanging
2026  if(pressure_dof_is_hanging[l])
2027  {
2028  // Get the local equation from the master node
2029  local_eqn =
2030  this->local_hang_eqn(hang_info_pt->master_node_pt(m),p_index);
2031  // Get the weight for the node
2032  hang_weight = hang_info_pt->master_weight(m);
2033  }
2034  else
2035  {
2036  local_eqn = this->p_local_eqn(l);
2037  hang_weight = 1.0;
2038  }
2039 
2040  // Cache the test function
2041  const double testp_ = testp[l];
2042 
2043  // If not a boundary conditions
2044  if(local_eqn >= 0)
2045  {
2046 
2047  // Loop over the two coordinate directions
2048  for(unsigned p=0;p<2;p++)
2049  {
2050  // Loop over shape controlling nodes
2051  for(unsigned q=0;q<n_shape_controlling_node;q++)
2052  {
2053 
2054  // Residual x deriv of Jacobian
2055  //-----------------------------
2056 
2057  // Source term
2058  double aux = -r*source;
2059 
2060  // Gradient terms
2061  aux += (interpolated_u[0]
2062  + r*interpolated_dudx(0,0)
2063  + r*interpolated_dudx(1,1));
2064 
2065  // Multiply through by deriv of Jacobian and integration weight
2066  dresidual_dnodal_coordinates(local_eqn,p,q)
2067  += aux*dJ_dX(p,q)*testp_*w*hang_weight;
2068 
2069  // Derivative of residual x Jacobian
2070  // ---------------------------------
2071 
2072  // Gradient of source function
2073  aux = -r*source_gradient[p]*psif[q];
2074  if(p==0) { aux -= source*psif[q]; }
2075 
2076  // Gradient terms
2077  aux += r*(d_dudx_dX(p,q,0,0) + d_dudx_dX(p,q,1,1));
2078  if(p==0)
2079  {
2080  aux += (interpolated_dudx(0,0)
2081  + interpolated_dudx(1,1))*psif[q];
2082  }
2083 
2084  // Multiply through by Jacobian and integration weight
2085  dresidual_dnodal_coordinates(local_eqn,p,q)
2086  += aux*testp_*J*w*hang_weight;
2087  }
2088  }
2089 
2090  // Derivs w.r.t. to nodal velocities
2091  // ---------------------------------
2092  if(element_has_node_with_aux_node_update_fct)
2093  {
2094  // Loop over local nodes
2095  for(unsigned q_local=0;q_local<n_node;q_local++)
2096  {
2097  // Number of master nodes and storage for the weight of
2098  // the shape function
2099  unsigned n_master2=1;
2100  double hang_weight2=1.0;
2101  HangInfo* hang_info2_pt=0;
2102 
2103  // Local boolean to indicate whether the node is hanging
2104  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
2105 
2106  Node* actual_shape_controlling_node_pt=node_pt(q_local);
2107 
2108  // If the node is hanging
2109  if(is_node_hanging2)
2110  {
2111  hang_info2_pt = node_pt(q_local)->hanging_pt();
2112 
2113  // Read out number of master nodes from hanging data
2114  n_master2 = hang_info2_pt->nmaster();
2115  }
2116  // Otherwise the node is its own master
2117  else
2118  {
2119  n_master2 = 1;
2120  }
2121 
2122  // Loop over the master nodes
2123  for(unsigned mm=0;mm<n_master2;mm++)
2124  {
2125 
2126  if(is_node_hanging2)
2127  {
2128  actual_shape_controlling_node_pt=
2129  hang_info2_pt->master_node_pt(mm);
2130  hang_weight2 = hang_info2_pt->master_weight(mm);
2131  }
2132 
2133  // Find the corresponding number
2134  unsigned q=local_shape_controlling_node_lookup[
2135  actual_shape_controlling_node_pt];
2136 
2137  // Loop over the two coordinate directions
2138  for(unsigned p=0;p<2;p++)
2139  {
2140  double aux =
2141  d_U_dX(p,q,0)*(psif[q_local] + r*dpsifdx(q_local,0))
2142  + d_U_dX(p,q,1)*r*dpsifdx(q_local,1);
2143 
2144  // Multiply through by Jacobian, test function and
2145  // integration weight
2146  dresidual_dnodal_coordinates(local_eqn,p,q) +=
2147  aux*testp_*J*w*hang_weight*hang_weight2;
2148  }
2149  } // End of loop over (mm) master nodes
2150  } // End of loop over local nodes q_local
2151  } // End of if(element_has_node_with_aux_node_update_fct)
2152  } // End of if it's not a boundary condition
2153  } // End of loop over (m) master nodes
2154  } // End of loop over shape functions for continuity eqn
2155 
2156  } // End of loop over integration points
2157 
2158 } // End of get_dresidual_dnodal_coordinates(...)
2159 
2160 
2161 
2162 } // End of namespace
virtual void pshape_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
virtual unsigned u_index_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.
virtual double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and derivatives w.r.t. global coords at ipt-th integration point Return J...
double get_source_fct(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
const double & re_invfr() const
Global inverse Froude number.
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
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 void get_source_fct_gradient(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2458
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition: elements.h:2226
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed...
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:557
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
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
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
virtual Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interp...
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_body_force_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Calculate the body force fct at a given time and Eulerian position.
virtual void get_body_force_gradient_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
Get gradient of body force term at (Eulerian) position x. Computed via function pointer (if set) or b...
virtual int p_local_eqn(const unsigned &n) const =0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. This function computes these terms analytically and overwrites the default implementation in the FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}.
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
virtual unsigned npres_axi_nst() const =0
Function to return number of pressure degrees of freedom.
unsigned nshape_controlling_nodes()
Number of shape-controlling nodes = the number of non-hanging nodes plus the number of master nodes a...
A Rank 3 Tensor class.
Definition: matrices.h:1337
double du_dt_axi_nst(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes...
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:540
static char t char * s
Definition: cfortran.h:572
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:130
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
Class that contains data for hanging nodes.
Definition: nodes.h:684
const double & re() const
Reynolds number.
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2215
std::map< Node *, unsigned > shape_controlling_node_lookup()
Return lookup scheme for unique number associated with any of the nodes that actively control the sha...
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:834
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
const Vector< double > & g() const
Vector of gravitational components.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
virtual int p_nodal_index_axi_nst() const
Which nodal value represents the pressure?
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
virtual double p_axi_nst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes...
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164
void fill_in_generic_residual_contribution_axi_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add element's contribution to the elemental residual vector and/or Jacobian matrix and mass matrix fl...