generalised_newtonian_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: 1194 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-27 08:44:43 +0100 (Fri, 27 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
31 
32 
33 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  // The strainrate used to compute the second invariant
196  DenseMatrix<double> strainrate_to_compute_second_invariant(3,3,0.0);
197 
198  // the strainrate used to calculate the second invariant
199  // can be either the current one or the one extrapolated from
200  // previous velocity values
202  {
203  strain_rate(s,strainrate_to_compute_second_invariant);
204  }
205  else
206  {
207  extrapolated_strain_rate(ipt,strainrate_to_compute_second_invariant);
208  }
209 
210  // Calculate the second invariant
212  strainrate_to_compute_second_invariant);
213 
214  // Get the viscosity according to the constitutive equation
215  double viscosity=Constitutive_eqn_pt->viscosity(second_invariant);
216 
217  //Get the user-defined body force terms
218  Vector<double> body_force(DIM+1);
219  get_body_force_axi_nst(time,ipt,s,interpolated_x,body_force);
220 
221  //Get the user-defined source function
222  double source=get_source_fct(time,ipt,interpolated_x);
223 
224  // r is the first postition component
225  double r = interpolated_x[0];
226 
227  // obacht set up vectors of the viscosity differentiated w.r.t.
228  // the velocity components (radial, axial, azimuthal)
229  Vector<double> dviscosity_dUr(n_node,0.0);
230  Vector<double> dviscosity_dUz(n_node,0.0);
231  Vector<double> dviscosity_dUphi(n_node,0.0);
232 
234  {
235  // Calculate the derivate of the viscosity w.r.t. the second invariant
236  double dviscosity_dsecond_invariant=
237  Constitutive_eqn_pt->dviscosity_dinvariant(second_invariant);
238 
239  // calculate a reference strainrate
240  DenseMatrix<double> strainrate_ref(3,3,0.0);
241  strain_rate(s,strainrate_ref);
242 
243  // pre-compute the derivative of the second invariant w.r.t. the
244  // entries in the rate of strain tensor
245  DenseMatrix<double> dinvariant_dstrainrate(3,3,0.0);
246 
247  // d I_2 / d epsilon_{r,r}
248  dinvariant_dstrainrate(0,0)=strainrate_ref(1,1)+strainrate_ref(2,2);
249  // d I_2 / d epsilon_{z,z}
250  dinvariant_dstrainrate(1,1)=strainrate_ref(0,0)+strainrate_ref(2,2);
251  // d I_2 / d epsilon_{phi,phi}
252  dinvariant_dstrainrate(2,2)=strainrate_ref(0,0)+strainrate_ref(1,1);
253  // d I_2 / d epsilon_{r,z}
254  dinvariant_dstrainrate(0,1)=-strainrate_ref(1,0);
255  // d I_2 / d epsilon_{z,r}
256  dinvariant_dstrainrate(1,0)=-strainrate_ref(0,1);
257  // d I_2 / d epsilon_{r,phi}
258  dinvariant_dstrainrate(0,2)=-strainrate_ref(2,0);
259  // d I_2 / d epsilon_{phi,r}
260  dinvariant_dstrainrate(2,0)=-strainrate_ref(0,2);
261  // d I_2 / d epsilon_{phi,z}
262  dinvariant_dstrainrate(2,1)=-strainrate_ref(1,2);
263  // d I_2 / d epsilon_{z,phi}
264  dinvariant_dstrainrate(1,2)=-strainrate_ref(2,1);
265 
266  // loop over the nodes
267  for(unsigned l=0;l<n_node;l++)
268  {
269  // Get pointer to l-th local node
270  //Node* nod_pt = node_pt(l);
271 
272  // loop over the three velocity components
273  for(unsigned i=0;i<3;i++)
274  {
275  // initialise the derivative of the second invariant w.r.t. the
276  // unknown velocity U_{i,l}
277  double dinvariant_dunknown=0.0;
278 
279  // loop over the entries of the rate of strain tensor
280  for(unsigned m=0;m<3;m++)
281  {
282  for(unsigned n=0;n<3;n++)
283  {
284 
285  // initialise the derivative of the strainrate w.r.t. the
286  // unknown velocity U_{i,l}
287  double dstrainrate_dunknown=0.0;
288 
289  // switch based on first index
290  switch(m)
291  {
292  // epsilon_{r ...}
293  case 0:
294 
295  // switch for second index
296  switch(n)
297  {
298  // epsilon_{r r}
299  case 0:
300  if(i==0)
301  {
302  dstrainrate_dunknown=dpsifdx(l,0);
303  }
304  break;
305 
306  // epsilon_{r z}
307  case 1:
308  if(i==0)
309  {
310  dstrainrate_dunknown=0.5*dpsifdx(l,1);
311  }
312  else if(i==1)
313  {
314  dstrainrate_dunknown=0.5*dpsifdx(l,0);
315  }
316  break;
317 
318  // epsilon_{r phi}
319  case 2:
320  if(i==2)
321  {
322  dstrainrate_dunknown=0.5*dpsifdx(l,0)-0.5/r*psif[l];
323  }
324  break;
325 
326  default:
327  std::ostringstream error_stream;
328  error_stream << "Should never get here...";
329  throw OomphLibError(
330  error_stream.str(),
331  OOMPH_CURRENT_FUNCTION,
332  OOMPH_EXCEPTION_LOCATION);
333 
334  }
335 
336  break;
337 
338  // epsilon_{z ...}
339  case 1:
340 
341  // switch for second index
342  switch(n)
343  {
344  // epsilon_{z r}
345  case 0:
346  if(i==0)
347  {
348  dstrainrate_dunknown=0.5*dpsifdx(l,1);
349  }
350  else if(i==1)
351  {
352  dstrainrate_dunknown=0.5*dpsifdx(l,0);
353  }
354  break;
355 
356  // epsilon_{z z}
357  case 1:
358  if(i==1)
359  {
360  dstrainrate_dunknown=dpsifdx(l,1);
361  }
362  else
363  {
364  //dstrainrate_dunknown=0.0;
365  }
366  break;
367 
368  // epsilon_{z phi}
369  case 2:
370  if(i==2)
371  {
372  dstrainrate_dunknown=0.5*dpsifdx(l,1);
373  }
374  break;
375 
376  default:
377  std::ostringstream error_stream;
378  error_stream << "Should never get here...";
379  throw OomphLibError(
380  error_stream.str(),
381  OOMPH_CURRENT_FUNCTION,
382  OOMPH_EXCEPTION_LOCATION);
383 
384  }
385 
386  break;
387 
388  // epsilon_{phi ...}
389  case 2:
390 
391  // switch for second index
392  switch(n)
393  {
394  // epsilon_{phi r}
395  case 0:
396  if(i==2)
397  {
398  dstrainrate_dunknown=0.5*dpsifdx(l,0)-0.5/r*psif[l];
399  }
400  break;
401 
402  // epsilon_{phi z}
403  case 1:
404  if(i==2)
405  {
406  dstrainrate_dunknown=0.5*dpsifdx(l,1);
407  }
408  break;
409 
410  // epsilon_{phi phi}
411  case 2:
412  if(i==0)
413  {
414  dstrainrate_dunknown=1.0/r*psif[l];
415  }
416  break;
417 
418  default:
419  std::ostringstream error_stream;
420  error_stream << "Should never get here...";
421  throw OomphLibError(
422  error_stream.str(),
423  OOMPH_CURRENT_FUNCTION,
424  OOMPH_EXCEPTION_LOCATION);
425 
426  }
427 
428  break;
429 
430  default:
431  std::ostringstream error_stream;
432  error_stream << "Should never get here...";
433  throw OomphLibError(
434  error_stream.str(),
435  OOMPH_CURRENT_FUNCTION,
436  OOMPH_EXCEPTION_LOCATION);
437 
438  }
439 
440  dinvariant_dunknown += dinvariant_dstrainrate(m,n)*
441  dstrainrate_dunknown;
442  }
443  }
444 
445  switch(i)
446  {
447  case 0:
448  dviscosity_dUr[l] = dviscosity_dsecond_invariant*dinvariant_dunknown;
449  break;
450 
451  case 1:
452  dviscosity_dUz[l] = dviscosity_dsecond_invariant*dinvariant_dunknown;
453  break;
454 
455  case 2:
456  dviscosity_dUphi[l] = dviscosity_dsecond_invariant*
457  dinvariant_dunknown;
458  break;
459 
460  default:
461  std::ostringstream error_stream;
462  error_stream << "Should never get here...";
463  throw OomphLibError(
464  error_stream.str(),
465  OOMPH_CURRENT_FUNCTION,
466  OOMPH_EXCEPTION_LOCATION);
467  }
468  }
469  }
470  }
471 
472 
473  //MOMENTUM EQUATIONS
474  //==================
475  //Number of master nodes and storage for the weight of the shape function
476  unsigned n_master=1; double hang_weight=1.0;
477 
478  //Loop over the nodes for the test functions/equations
479  //----------------------------------------------------
480  for(unsigned l=0;l<n_node;l++)
481  {
482  //Local boolean that indicates the hanging status of the node
483  bool is_node_hanging = node_pt(l)->is_hanging();
484 
485  //If the node is hanging
486  if(is_node_hanging)
487  {
488  //Get the hanging pointer
489  hang_info_pt = node_pt(l)->hanging_pt();
490  //Read out number of master nodes from hanging data
491  n_master = hang_info_pt->nmaster();
492  }
493  //Otherwise the node is its own master
494  else
495  {
496  n_master = 1;
497  }
498 
499  //Loop over the master nodes
500  for(unsigned m=0;m<n_master;m++)
501  {
502  // Loop over velocity components for equations
503  for(unsigned i=0;i<DIM+1;i++)
504  {
505  //Get the equation number
506  //If the node is hanging
507  if(is_node_hanging)
508  {
509  //Get the equation number from the master node
510  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
511  u_nodal_index[i]);
512  //Get the hang weight from the master node
513  hang_weight = hang_info_pt->master_weight(m);
514  }
515  //If the node is not hanging
516  else
517  {
518  // Local equation number
519  local_eqn = this->nodal_local_eqn(l,u_nodal_index[i]);
520  // Node contributes with full weight
521  hang_weight = 1.0;
522  }
523 
524  //If it's not a boundary condition...
525  if(local_eqn>= 0)
526  {
527  // initialise for residual calculation
528  double sum=0.0;
529 
530  switch (i)
531  {
532  // RADIAL MOMENTUM EQUATION
533  case 0:
534  //Add the user-defined body force terms
535  sum +=
536  r*body_force[0]*testf[l]*W*hang_weight;
537 
538  //Add the gravitational body force term
539  sum += r*scaled_re_inv_fr*testf[l]*G[0]*W*hang_weight;
540 
541  //Add the pressure gradient term
542  sum +=
543  interpolated_p*(testf[l] + r*dtestfdx(l,0))*W*hang_weight;
544 
545  //Add in the stress tensor terms
546  //The viscosity ratio needs to go in here to ensure
547  //continuity of normal stress is satisfied even in flows
548  //with zero pressure gradient!
549  // stress term 1
550  sum -= visc_ratio*viscosity*
551  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0)*W
552  *hang_weight;
553 
554  // stress term 2
555  sum -= visc_ratio*viscosity*r*
556  (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
557  dtestfdx(l,1)*W*hang_weight;
558 
559  // stress term 3
560  sum -=
561  visc_ratio*viscosity*(1.0 + Gamma[0])*interpolated_u[0]
562  *testf[l]*W*hang_weight/r;
563 
564  //Add in the inertial terms
565  //du/dt term
566  sum -= scaled_re_st*r*dudt[0]*testf[l]*W*hang_weight;
567 
568  //Convective terms
569  sum -=
570  scaled_re*(r*interpolated_u[0]*interpolated_dudx(0,0)
571  - interpolated_u[2]*interpolated_u[2]
572  + r*interpolated_u[1]*interpolated_dudx(0,1))*testf[l]*W*
573  hang_weight;
574 
575  //Mesh velocity terms
576  if(!ALE_is_disabled)
577  {
578  for(unsigned k=0;k<2;k++)
579  {
580  sum +=
581  scaled_re_st*r*mesh_veloc[k]*interpolated_dudx(0,k)*testf[l]*W*
582  hang_weight;
583  }
584  }
585 
586  //Coriolis term
587  sum +=
588  2.0*r*scaled_re_inv_ro*interpolated_u[2]*testf[l]*W*hang_weight;
589 
590  break;
591 
592  // AXIAL MOMENTUM EQUATION
593  case 1:
594  //If it's not a boundary condition
595  //Add the user-defined body force terms
596  //Remember to multiply by the density ratio!
597  sum += r*body_force[1]*testf[l]*W*hang_weight;
598 
599  //Add the gravitational body force term
600  sum += r*scaled_re_inv_fr*testf[l]*G[1]*W*hang_weight;
601 
602  //Add the pressure gradient term
603  sum += r*interpolated_p*dtestfdx(l,1)*W*hang_weight;
604 
605  //Add in the stress tensor terms
606  //The viscosity ratio needs to go in here to ensure
607  //continuity of normal stress is satisfied even in flows
608  //with zero pressure gradient!
609  // stress term 1
610  sum -= visc_ratio*viscosity*
611  r*(interpolated_dudx(1,0) +
612  Gamma[1]*interpolated_dudx(0,1))*dtestfdx(l,0)*W*
613  hang_weight;
614 
615  // stress term 2
616  sum -= visc_ratio*viscosity*r*
617  (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1)*W*
618  hang_weight;
619 
620  //Add in the inertial terms
621  //du/dt term
622  sum -= scaled_re_st*r*dudt[1]*testf[l]*W*hang_weight;
623 
624  //Convective terms
625  sum -=
626  scaled_re*(r*interpolated_u[0]*interpolated_dudx(1,0)
627  + r*interpolated_u[1]*interpolated_dudx(1,1))*
628  testf[l]*W*hang_weight;
629 
630  //Mesh velocity terms
631  if(!ALE_is_disabled)
632  {
633  for(unsigned k=0;k<2;k++)
634  {
635  sum +=
636  scaled_re_st*r*mesh_veloc[k]*interpolated_dudx(1,k)*testf[l]*W
637  *hang_weight;
638  }
639  }
640  break;
641 
642  // AZIMUTHAL MOMENTUM EQUATION
643  case 2:
644  //Add the user-defined body force terms
645  //Remember to multiply by the density ratio!
646  sum += r*body_force[2]*testf[l]*W*hang_weight;
647 
648  //Add the gravitational body force term
649  sum += r*scaled_re_inv_fr*testf[l]*G[2]*W*hang_weight;
650 
651  //There is NO pressure gradient term
652 
653  //Add in the stress tensor terms
654  //The viscosity ratio needs to go in here to ensure
655  //continuity of normal stress is satisfied even in flows
656  //with zero pressure gradient!
657  // stress term 1
658  sum -= visc_ratio*viscosity*
659  (r*interpolated_dudx(2,0) - Gamma[0]*interpolated_u[2])
660  *dtestfdx(l,0)*W*hang_weight;
661 
662  // stress term 2
663  sum -= visc_ratio*viscosity*r*interpolated_dudx(2,1)*
664  dtestfdx(l,1)*W*hang_weight;
665 
666  // stress term 3
667  sum -= visc_ratio*viscosity*
668  ((interpolated_u[2]/r) - Gamma[0]*interpolated_dudx(2,0))*
669  testf[l]*W*hang_weight;
670 
671 
672  //Add in the inertial terms
673  //du/dt term
674  sum -= scaled_re_st*r*dudt[2]*testf[l]*W*hang_weight;
675 
676  //Convective terms
677  sum -=
678  scaled_re*(r*interpolated_u[0]*interpolated_dudx(2,0)
679  + interpolated_u[0]*interpolated_u[2]
680  + r*interpolated_u[1]*interpolated_dudx(2,1))*testf[l]*W
681  *hang_weight;
682 
683  //Mesh velocity terms
684  if(!ALE_is_disabled)
685  {
686  for(unsigned k=0;k<2;k++)
687  {
688  sum += scaled_re_st*r*mesh_veloc[k]*
689  interpolated_dudx(2,k)*testf[l]*W*hang_weight;
690  }
691  }
692 
693  //Coriolis term
694  sum -=
695  2.0*r*scaled_re_inv_ro*interpolated_u[0]*testf[l]*W*hang_weight;
696 
697  break;
698  }
699 
700  // Add contribution
701  residuals[local_eqn] += sum;
702 
703  //CALCULATE THE JACOBIAN
704  if (flag)
705  {
706  //Number of master nodes and weights
707  unsigned n_master2=1; double hang_weight2=1.0;
708  //Loop over the velocity nodes for columns
709  for(unsigned l2=0;l2<n_node;l2++)
710  {
711  //Local boolean for hanging status
712  bool is_node2_hanging = node_pt(l2)->is_hanging();
713 
714  //If the node is hanging
715  if(is_node2_hanging)
716  {
717  hang_info2_pt = node_pt(l2)->hanging_pt();
718  //Read out number of master nodes from hanging data
719  n_master2 = hang_info2_pt->nmaster();
720  }
721  //Otherwise the node is its own master
722  else
723  {
724  n_master2 = 1;
725  }
726 
727  //Loop over the master nodes
728  for(unsigned m2=0;m2<n_master2;m2++)
729  {
730  //Loop over the velocity components
731  for(unsigned i2=0;i2<DIM+1;i2++)
732  {
733  //Get the number of the unknown
734  //If the node is hanging
735  if(is_node2_hanging)
736  {
737  //Get the equation number from the master node
738  local_unknown =
739  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
740  u_nodal_index[i2]);
741  hang_weight2 = hang_info2_pt->master_weight(m2);
742  }
743  else
744  {
745  local_unknown = this->nodal_local_eqn(l2,u_nodal_index[i2]);
746  hang_weight2 = 1.0;
747  }
748 
749  // If the unknown is non-zero
750  if(local_unknown >= 0)
751  {
752  //Different results for i and i2
753  switch (i)
754  {
755  // RADIAL MOMENTUM EQUATION
756  case 0:
757  switch (i2)
758  {
759  // radial component
760  case 0:
761 
762  //Add the mass matrix entries
763  if(flag==2)
764  {
765  mass_matrix(local_eqn,local_unknown) +=
766  scaled_re_st*r*psif[l2]*testf[l]*W
767  *hang_weight*hang_weight2;
768  }
769 
770  //Add contribution to the Jacobian matrix
771  // stress term 1
772  jacobian(local_eqn,local_unknown)
773  -= visc_ratio*viscosity*r*(1.0+Gamma[0])
774  *dpsifdx(l2,0)*dtestfdx(l,0)*W*
775  hang_weight*hang_weight2;
776 
778  {
779  // stress term 1 contribution from
780  // generalised Newtonian model
781  jacobian(local_eqn,local_unknown)
782  -= visc_ratio*dviscosity_dUr[l2]*
783  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*
784  dtestfdx(l,0)*W*hang_weight*hang_weight2;
785  }
786 
787  // stress term 2
788  jacobian(local_eqn,local_unknown)
789  -= visc_ratio*viscosity*r*dpsifdx(l2,1)*
790  dtestfdx(l,1)*W*hang_weight*hang_weight2;
791 
793  {
794  // stress term 2 contribution from
795  // generalised Newtonian model
796  jacobian(local_eqn,local_unknown)
797  -= visc_ratio*dviscosity_dUr[l2]*
798  r*(interpolated_dudx(0,1) + Gamma[0]*
799  interpolated_dudx(1,0))*
800  dtestfdx(l,1)*W*hang_weight*hang_weight2;
801  }
802 
803  // stress term 3
804  jacobian(local_eqn,local_unknown)
805  -= visc_ratio*viscosity*(1.0 + Gamma[0])*psif[l2]*
806  testf[l]*W*hang_weight*hang_weight2/r;
807 
809  {
810  // stress term 3 contribution from
811  // generalised Newtonian model
812  jacobian(local_eqn,local_unknown)
813  -= visc_ratio*dviscosity_dUr[l2]*
814  (1.0 + Gamma[0])*interpolated_u[0]*testf[l]*W*
815  hang_weight*hang_weight2/r;
816  }
817 
818  //Add in the inertial terms
819  //du/dt term
820  jacobian(local_eqn,local_unknown)
821  -= scaled_re_st*r*node_pt(l2)->time_stepper_pt()->
822  weight(1,0)*psif[l2]*testf[l]*W
823  *hang_weight*hang_weight2;
824 
825  //Convective terms
826  jacobian(local_eqn,local_unknown) -=
827  scaled_re*(r*psif[l2]*interpolated_dudx(0,0)
828  + r*interpolated_u[0]*dpsifdx(l2,0)
829  + r*interpolated_u[1]*dpsifdx(l2,1))
830  *testf[l]*W*hang_weight*hang_weight2;
831 
832  //Mesh velocity terms
833  if(!ALE_is_disabled)
834  {
835  for(unsigned k=0;k<2;k++)
836  {
837  jacobian(local_eqn,local_unknown)
838  += scaled_re_st*r*mesh_veloc[k]*dpsifdx(l2,k)
839  *testf[l]*W*hang_weight*hang_weight2;
840  }
841  }
842  break;
843 
844  // axial component
845  case 1:
846 
847  // no stress term 1
848 
850  {
851  // stress term 1 contribution from
852  // generalised Newtonian model
853  jacobian(local_eqn,local_unknown)
854  -= visc_ratio*dviscosity_dUz[l2]*
855  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*
856  dtestfdx(l,0)*W*hang_weight*hang_weight2;
857  }
858 
859  // stress term 2
860  jacobian(local_eqn,local_unknown) -=
861  visc_ratio*viscosity*r*Gamma[0]*dpsifdx(l2,0)
862  *dtestfdx(l,1)*W*hang_weight*hang_weight2;
863 
865  {
866  // stress term 2 contribution from
867  // generalised Newtonian model
868  jacobian(local_eqn,local_unknown) -=
869  visc_ratio*dviscosity_dUz[l2]*
870  r*(interpolated_dudx(0,1) + Gamma[0]*
871  interpolated_dudx(1,0))*
872  dtestfdx(l,1)*W*hang_weight*hang_weight2;
873  }
874 
875  // no stress term 3
876 
878  {
879  // stress term 3 contribution from
880  // generalised Newtonian model
881  jacobian(local_eqn,local_unknown)
882  -= visc_ratio*dviscosity_dUz[l2]*
883  (1.0 + Gamma[0])*interpolated_u[0]*testf[l]*
884  W*hang_weight*hang_weight2/r;
885  }
886 
887  //Convective terms
888  jacobian(local_eqn,local_unknown) -=
889  scaled_re*r*psif[l2]*interpolated_dudx(0,1)*
890  testf[l]*W*hang_weight*hang_weight2;
891  break;
892 
893  // azimuthal component
894  case 2:
895 
896  // no stress term 1
897 
899  {
900  // stress term 1 contribution from
901  // generalised Newtonian model
902  jacobian(local_eqn,local_unknown)
903  -= visc_ratio*dviscosity_dUphi[l2]*
904  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*
905  dtestfdx(l,0)*W*hang_weight*hang_weight2;
906  }
907 
908  // no stress term 2
909 
911  {
912  // stress term 2 contribution from
913  // generalised Newtonian model
914  jacobian(local_eqn,local_unknown) -=
915  visc_ratio*dviscosity_dUphi[l2]*
916  r*(interpolated_dudx(0,1) + Gamma[0]*
917  interpolated_dudx(1,0))*
918  dtestfdx(l,1)*W*hang_weight*hang_weight2;
919  }
920 
921  // no stress term 3
922 
924  {
925  // stress term 3 contribution from
926  // generalised Newtonian model
927  jacobian(local_eqn,local_unknown)
928  -= visc_ratio*dviscosity_dUphi[l2]*
929  (1.0 + Gamma[0])*interpolated_u[0]*testf[l]*
930  W*hang_weight*hang_weight2/r;
931  }
932 
933  //Convective terms
934  jacobian(local_eqn,local_unknown) -=
935  - scaled_re*2.0*interpolated_u[2]*psif[l2]
936  *testf[l]*W*hang_weight*hang_weight2;
937 
938  //Coriolis terms
939  jacobian(local_eqn,local_unknown) +=
940  2.0*r*scaled_re_inv_ro*psif[l2]*testf[l]*W
941  *hang_weight*hang_weight2;
942 
943  break;
944  } /*End of contribution radial momentum eqn*/
945  break;
946 
947  // AXIAL MOMENTUM EQUATION
948  case 1:
949  switch (i2)
950  {
951  // radial component
952  case 0:
953  //Add in the stress tensor terms
954  //The viscosity ratio needs to go in here to ensure
955  //continuity of normal stress is satisfied even
956  //in flows with zero pressure gradient!
957  // stress term 1
958  jacobian(local_eqn,local_unknown) -=
959  visc_ratio*viscosity*r*Gamma[1]*dpsifdx(l2,1)*
960  dtestfdx(l,0)*W*hang_weight*hang_weight2;
961 
963  {
964  // stress term 1 contribution from
965  // generalised Newtonian model
966  jacobian(local_eqn,local_unknown) -=
967  visc_ratio*dviscosity_dUr[l2]*
968  r*(interpolated_dudx(1,0) + Gamma[1]*
969  interpolated_dudx(0,1))
970  *dtestfdx(l,0)*W*hang_weight*hang_weight2;
971  }
972 
973  // no stress term 2
974 
976  {
977  // stress term 2 contribution from
978  // generalised Newtonian model
979  jacobian(local_eqn,local_unknown) -=
980  visc_ratio*dviscosity_dUr[l2]*
981  r*(1.0 + Gamma[1])*interpolated_dudx(1,1)*
982  dtestfdx(l,1)*W*hang_weight*hang_weight2;
983  }
984 
985  //Convective terms
986  jacobian(local_eqn,local_unknown) -=
987  scaled_re*r*psif[l2]*interpolated_dudx(1,0)*testf[l]*W
988  *hang_weight*hang_weight2;
989  break;
990 
991  // axial component
992  case 1:
993 
994  //Add the mass matrix terms
995  if(flag==2)
996  {
997  mass_matrix(local_eqn,local_unknown) +=
998  scaled_re_st*r*psif[l2]*testf[l]*W
999  *hang_weight*hang_weight2;
1000  }
1001 
1002  //Add in the stress tensor terms
1003  //The viscosity ratio needs to go in here to ensure
1004  //continuity of normal stress is satisfied even in
1005  //flows with zero pressure gradient!
1006  // stress term 1
1007  jacobian(local_eqn,local_unknown) -=
1008  visc_ratio*viscosity*r*dpsifdx(l2,0)*dtestfdx(l,0)*W
1009  *hang_weight*hang_weight2;
1010 
1012  {
1013  // stress term 1 contribution from
1014  // generalised Newtonian model
1015  jacobian(local_eqn,local_unknown) -=
1016  visc_ratio*dviscosity_dUz[l2]*
1017  r*(interpolated_dudx(1,0) + Gamma[1]*
1018  interpolated_dudx(0,1))
1019  *dtestfdx(l,0)*W*hang_weight*hang_weight2;
1020  }
1021 
1022  // stress term 2
1023  jacobian(local_eqn,local_unknown) -=
1024  visc_ratio*viscosity*r*(1.0 + Gamma[1])*dpsifdx(l2,1)*
1025  dtestfdx(l,1)*W*hang_weight*hang_weight2;
1026 
1028  {
1029  // stress term 2 contribution from
1030  // generalised Newtonian model
1031  jacobian(local_eqn,local_unknown) -=
1032  visc_ratio*dviscosity_dUz[l2]*
1033  r*(1.0 + Gamma[1])*interpolated_dudx(1,1)*
1034  dtestfdx(l,1)*W*hang_weight*hang_weight2;
1035  }
1036 
1037  //Add in the inertial terms
1038  //du/dt term
1039  jacobian(local_eqn,local_unknown) -=
1040  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->
1041  weight(1,0)*psif[l2]*testf[l]*W
1042  *hang_weight*hang_weight2;
1043 
1044  //Convective terms
1045  jacobian(local_eqn,local_unknown) -=
1046  scaled_re*(r*interpolated_u[0]*dpsifdx(l2,0)
1047  + r*psif[l2]*interpolated_dudx(1,1)
1048  + r*interpolated_u[1]*dpsifdx(l2,1))
1049  *testf[l]*W*hang_weight*hang_weight2;
1050 
1051  //Mesh velocity terms
1052  if(!ALE_is_disabled)
1053  {
1054  for(unsigned k=0;k<2;k++)
1055  {
1056  jacobian(local_eqn,local_unknown)
1057  += scaled_re_st*r*mesh_veloc[k]*dpsifdx(l2,k)
1058  *testf[l]*W*hang_weight*hang_weight2;
1059  }
1060  }
1061  break;
1062 
1063  // azimuthal component
1064  case 2:
1065  //There are no azimithal terms in the axial
1066  //momentum equation
1067  // edit: for the generalised Newtonian elements
1068  // there are...
1069 
1071  {
1072  // stress term 1 contribution from
1073  // generalised Newtonian model
1074  jacobian(local_eqn,local_unknown) -=
1075  visc_ratio*dviscosity_dUphi[l2]*
1076  r*(interpolated_dudx(1,0) + Gamma[1]*
1077  interpolated_dudx(0,1))
1078  *dtestfdx(l,0)*W*hang_weight*hang_weight2;
1079 
1080  // stress term 2 contribution from
1081  // generalised Newtonian model
1082  jacobian(local_eqn,local_unknown) -=
1083  visc_ratio*dviscosity_dUphi[l2]*
1084  r*(1.0 + Gamma[1])*interpolated_dudx(1,1)*
1085  dtestfdx(l,1)*W*hang_weight*hang_weight2;
1086  }
1087  break;
1088  }
1089  break;
1090 
1091  // AZIMUTHAL MOMENTUM EQUATION
1092  case 2:
1093  switch (i2)
1094  {
1095  // radial component
1096  case 0:
1098  {
1099  // stress term 1 contribution from
1100  // generalised Newtonian model
1101  jacobian(local_eqn,local_unknown) -=
1102  visc_ratio*dviscosity_dUr[l2]*
1103  (r*interpolated_dudx(2,0) -
1104  Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*
1105  W*hang_weight*hang_weight2;
1106 
1107  // stress term 2 contribution from
1108  // generalised Newtonian model
1109  jacobian(local_eqn,local_unknown) -=
1110  visc_ratio*dviscosity_dUr[l2]*
1111  r*interpolated_dudx(2,1)*dtestfdx(l,1)*
1112  W*hang_weight*hang_weight2;
1113 
1114  // stress term 3 contribution from
1115  // generalised Newtonian model
1116  jacobian(local_eqn,local_unknown) -=
1117  visc_ratio*dviscosity_dUr[l2]*
1118  ((interpolated_u[2]/r) -Gamma[0]*
1119  interpolated_dudx(2,0))*
1120  testf[l]*W*hang_weight*hang_weight2;
1121  }
1122 
1123  //Convective terms
1124  jacobian(local_eqn,local_unknown) -=
1125  scaled_re*(r*psif[l2]*interpolated_dudx(2,0)
1126  + psif[l2]*interpolated_u[2])*testf[l]
1127  *W*hang_weight*hang_weight2;
1128 
1129  //Coriolis term
1130  jacobian(local_eqn,local_unknown) -=
1131  2.0*r*scaled_re_inv_ro*psif[l2]*testf[l]*W
1132  *hang_weight*hang_weight2;
1133 
1134  break;
1135 
1136  // axial component
1137  case 1:
1139  {
1140  // stress term 1 contribution from
1141  // generalised Newtonian model
1142  jacobian(local_eqn,local_unknown) -=
1143  visc_ratio*dviscosity_dUz[l2]*
1144  (r*interpolated_dudx(2,0) -
1145  Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*W*
1146  hang_weight*hang_weight2;
1147 
1148  // stress term 2 contribution from
1149  // generalised Newtonian model
1150  jacobian(local_eqn,local_unknown) -=
1151  visc_ratio*dviscosity_dUz[l2]*
1152  r*interpolated_dudx(2,1)*dtestfdx(l,1)*
1153  W*hang_weight*hang_weight2;
1154 
1155  // stress term 3 contribution from
1156  // generalised Newtonian model
1157  jacobian(local_eqn,local_unknown) -=
1158  visc_ratio*dviscosity_dUz[l2]*
1159  ((interpolated_u[2]/r) -Gamma[0]*
1160  interpolated_dudx(2,0))*
1161  testf[l]*W*hang_weight*hang_weight2;
1162  }
1163 
1164  //Convective terms
1165  jacobian(local_eqn,local_unknown) -=
1166  scaled_re*r*psif[l2]*interpolated_dudx(2,1)*testf[l]
1167  *W*hang_weight*hang_weight2;
1168  break;
1169 
1170  // azimuthal component
1171  case 2:
1172 
1173  //Add the mass matrix terms
1174  if(flag==2)
1175  {
1176  mass_matrix(local_eqn,local_unknown) +=
1177  scaled_re_st*r*psif[l2]*testf[l]*W
1178  *hang_weight*hang_weight2;
1179  }
1180 
1181  //Add in the stress tensor terms
1182  //The viscosity ratio needs to go in here to ensure
1183  //continuity of normal stress is satisfied even in
1184  //flows with zero pressure gradient!
1185  // stress term 1
1186  jacobian(local_eqn,local_unknown) -=
1187  visc_ratio*viscosity*
1188  (r*dpsifdx(l2,0) - Gamma[0]*psif[l2])
1189  *dtestfdx(l,0)*W*hang_weight*hang_weight2;
1190 
1192  {
1193  // stress term 1 contribution from
1194  // generalised Newtonian model
1195  jacobian(local_eqn,local_unknown) -=
1196  visc_ratio*dviscosity_dUphi[l2]*
1197  (r*interpolated_dudx(2,0) -
1198  Gamma[0]*interpolated_u[2])*dtestfdx(l,0)*
1199  W*hang_weight*hang_weight2;
1200  }
1201 
1202  // stress term 2
1203  jacobian(local_eqn,local_unknown) -=
1204  visc_ratio*viscosity*r*dpsifdx(l2,1)*dtestfdx(l,1)
1205  *W*hang_weight*hang_weight2;
1206 
1208  {
1209  // stress term 2 contribution from
1210  // generalised Newtonian model
1211  jacobian(local_eqn,local_unknown) -=
1212  visc_ratio*dviscosity_dUphi[l2]*
1213  r*interpolated_dudx(2,1)*dtestfdx(l,1)*
1214  W*hang_weight*hang_weight2;
1215  }
1216 
1217  // stress term 3
1218  jacobian(local_eqn,local_unknown) -=
1219  visc_ratio*viscosity*
1220  ((psif[l2]/r) - Gamma[0]*dpsifdx(l2,0))
1221  *testf[l]*W*hang_weight*hang_weight2;
1222 
1224  {
1225  // stress term 3 contribution from
1226  // generalised Newtonian model
1227  jacobian(local_eqn,local_unknown) -=
1228  visc_ratio*dviscosity_dUphi[l2]*
1229  ((interpolated_u[2]/r) -Gamma[0]*
1230  interpolated_dudx(2,0))*
1231  testf[l]*W*hang_weight*hang_weight2;
1232  }
1233 
1234  //Add in the inertial terms
1235  //du/dt term
1236  jacobian(local_eqn,local_unknown) -=
1237  scaled_re_st*r*node_pt(l2)
1238  ->time_stepper_pt()->weight(1,0)*
1239  psif[l2]*testf[l]*W*hang_weight*hang_weight2;
1240 
1241  //Convective terms
1242  jacobian(local_eqn,local_unknown) -=
1243  scaled_re*(r*interpolated_u[0]*dpsifdx(l2,0)
1244  + interpolated_u[0]*psif[l2]
1245  + r*interpolated_u[1]*dpsifdx(l2,1))
1246  *testf[l]*W*hang_weight*hang_weight2;
1247 
1248  //Mesh velocity terms
1249  if(!ALE_is_disabled)
1250  {
1251  for(unsigned k=0;k<2;k++)
1252  {
1253  jacobian(local_eqn,local_unknown)
1254  += scaled_re_st*r*mesh_veloc[k]*dpsifdx(l2,k)
1255  *testf[l]*W*hang_weight*hang_weight2;
1256  }
1257  }
1258  break;
1259  }
1260  break;
1261  }
1262  }
1263  }
1264  }
1265  } //End of loop over the nodes
1266 
1267 
1268  //Loop over the pressure shape functions
1269  for(unsigned l2=0;l2<n_pres;l2++)
1270  {
1271  //If the pressure dof is hanging
1272  if(pressure_dof_is_hanging[l2])
1273  {
1274  // Pressure dof is hanging so it must be nodal-based
1275  hang_info2_pt =
1276  pressure_node_pt(l2)->hanging_pt(p_index);
1277 
1278  //Get the number of master nodes from the pressure node
1279  n_master2 = hang_info2_pt->nmaster();
1280  }
1281  //Otherwise the node is its own master
1282  else
1283  {
1284  n_master2 = 1;
1285  }
1286 
1287  //Loop over the master nodes
1288  for(unsigned m2=0;m2<n_master2;m2++)
1289  {
1290  //Get the number of the unknown
1291  //If the pressure dof is hanging
1292  if(pressure_dof_is_hanging[l2])
1293  {
1294  //Get the unknown from the node
1295  local_unknown =
1296  local_hang_eqn(hang_info2_pt->master_node_pt(m2),
1297  p_index);
1298  //Get the weight from the hanging object
1299  hang_weight2 = hang_info2_pt->master_weight(m2);
1300  }
1301  else
1302  {
1303  local_unknown = p_local_eqn(l2);
1304  hang_weight2 = 1.0;
1305  }
1306 
1307  //If the unknown is not pinned
1308  if(local_unknown >= 0)
1309  {
1310  //Add in contributions to different equations
1311  switch (i)
1312  {
1313  // RADIAL MOMENTUM EQUATION
1314  case 0:
1315  jacobian(local_eqn,local_unknown)
1316  += psip[l2]*(testf[l] + r*dtestfdx(l,0))*W*
1317  hang_weight*hang_weight2;
1318  break;
1319 
1320  // AXIAL MOMENTUM EQUATION
1321  case 1:
1322  jacobian(local_eqn,local_unknown)
1323  += r*psip[l2]*dtestfdx(l,1)*W*hang_weight*hang_weight2;
1324  break;
1325 
1326  // AZIMUTHAL MOMENTUM EQUATION
1327  case 2:
1328  break;
1329  }
1330  }
1331  }
1332  } //End of loop over pressure dofs
1333  }// End of Jacobian calculation
1334  } //End of if not boundary condition statement
1335  } //End of loop over velocity components
1336  } //End of loop over master nodes
1337  } //End of loop over nodes
1338 
1339 
1340  //CONTINUITY EQUATION
1341  //===================
1342 
1343  //Loop over the pressure shape functions
1344  for(unsigned l=0;l<n_pres;l++)
1345  {
1346  //If the pressure dof is hanging
1347  if(pressure_dof_is_hanging[l])
1348  {
1349  // Pressure dof is hanging so it must be nodal-based
1350  hang_info_pt = pressure_node_pt(l)->hanging_pt(p_index);
1351  //Get the number of master nodes from the pressure node
1352  n_master = hang_info_pt->nmaster();
1353  }
1354  //Otherwise the node is its own master
1355  else
1356  {
1357  n_master = 1;
1358  }
1359 
1360  //Loop over the master nodes
1361  for(unsigned m=0;m<n_master;m++)
1362  {
1363  //Get the number of the unknown
1364  //If the pressure dof is hanging
1365  if(pressure_dof_is_hanging[l])
1366  {
1367  local_eqn = local_hang_eqn(hang_info_pt->master_node_pt(m),p_index);
1368  hang_weight = hang_info_pt->master_weight(m);
1369  }
1370  else
1371  {
1372  local_eqn = p_local_eqn(l);
1373  hang_weight = 1.0;
1374  }
1375 
1376  //If the equation is not pinned
1377  if(local_eqn >= 0)
1378  {
1379  // Source term
1380  residuals[local_eqn] -= r*source*testp[l]*W*hang_weight;
1381 
1382  //Gradient terms
1383  residuals[local_eqn] +=
1384  (interpolated_u[0] + r*interpolated_dudx(0,0)
1385  + r*interpolated_dudx(1,1))*testp[l]*W*hang_weight;
1386 
1387  //CALCULATE THE JACOBIAN
1388  //======================
1389  if (flag)
1390  {
1391  unsigned n_master2=1; double hang_weight2=1.0;
1392  //Loop over the velocity nodes for columns
1393  for(unsigned l2=0;l2<n_node;l2++)
1394  {
1395  //Local boolean to indicate whether the node is hanging
1396  bool is_node2_hanging = node_pt(l2)->is_hanging();
1397 
1398  //If the node is hanging
1399  if(is_node2_hanging)
1400  {
1401  hang_info2_pt = node_pt(l2)->hanging_pt();
1402  //Read out number of master nodes from hanging data
1403  n_master2 = hang_info2_pt->nmaster();
1404  }
1405  //Otherwise the node is its own master
1406  else
1407  {
1408  n_master2 = 1;
1409  }
1410 
1411  //Loop over the master nodes
1412  for(unsigned m2=0;m2<n_master2;m2++)
1413  {
1414  //Loop over the velocity components
1415  for(unsigned i2=0;i2<DIM+1;i2++)
1416  {
1417  //Get the number of the unknown
1418  //If the node is hanging
1419  if(is_node2_hanging)
1420  {
1421  //Get the equation number from the master node
1422  local_unknown =
1423  local_hang_eqn(hang_info2_pt->master_node_pt(m2),
1424  u_nodal_index[i2]);
1425  hang_weight2 = hang_info2_pt->master_weight(m2);
1426  }
1427  else
1428  {
1429  local_unknown = this->nodal_local_eqn(l2,u_nodal_index[i2]);
1430  hang_weight2 = 1.0;
1431  }
1432 
1433  //If the unknown is not pinned
1434  if(local_unknown >= 0)
1435  {
1436  switch (i2)
1437  {
1438  // radial component
1439  case 0:
1440  jacobian(local_eqn,local_unknown) +=
1441  (psif[l2] + r*dpsifdx(l2,0))*testp[l]*W*
1442  hang_weight*hang_weight2;
1443  break;
1444 
1445  // axial component
1446  case 1:
1447  jacobian(local_eqn,local_unknown) +=
1448  r*dpsifdx(l2,1)*testp[l]*W*hang_weight*hang_weight2;
1449  break;
1450 
1451  // azimuthal component
1452  case 2:
1453  break;
1454  }
1455  }
1456  }
1457  }
1458  } //End of loop over nodes
1459 
1460  //NO PRESSURE CONTRIBUTIONS TO CONTINUITY EQUATION
1461  } //End of Jacobian calculation
1462  }
1463  }
1464  } // End of loop over pressure nodes
1465 
1466 } // End of loop over integration points
1467 
1468 
1469 } // End of fill_in_generic_residual_contribution_axi_nst(...)
1470 
1471 
1472 //======================================================================
1473 /// Compute derivatives of elemental residual vector with respect
1474 /// to nodal coordinates.
1475 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
1476 /// Overloads the FD-based version in the FiniteElement base class.
1477 //======================================================================
1480  RankThreeTensor<double>& dresidual_dnodal_coordinates)
1481 {
1482  // Create an Oomph Lib warning
1483  std::string warning_message = "This function has not been tested.\n";
1484  std::string function = "RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::\n";
1485  function += "get_dresidual_dnodal_coordinates(...)";
1486  OomphLibWarning(warning_message,function,OOMPH_EXCEPTION_LOCATION);
1487 
1488  // Return immediately if there are no dofs
1489  if(ndof()==0) { return; }
1490 
1491  // Determine number of nodes in element
1492  const unsigned n_node = nnode();
1493 
1494  // Get continuous time from timestepper of first node
1495  double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
1496 
1497  // Determine number of pressure dofs in element
1498  const unsigned n_pres = this->npres_axi_nst();
1499 
1500  // Find the indices at which the local velocities are stored
1501  unsigned u_nodal_index[3];
1502  for(unsigned i=0;i<3;i++) { u_nodal_index[i] = this->u_index_axi_nst(i); }
1503 
1504  // Which nodal value represents the pressure? (Negative if pressure
1505  // is not based on nodal interpolation).
1506  const int p_index = this->p_nodal_index_axi_nst();
1507 
1508  // Local array of booleans that are true if the l-th pressure value is
1509  // hanging (avoid repeated virtual function calls)
1510  bool pressure_dof_is_hanging[n_pres];
1511 
1512  // If the pressure is stored at a node
1513  if(p_index >= 0)
1514  {
1515  // Read out whether the pressure is hanging
1516  for(unsigned l=0;l<n_pres;++l)
1517  {
1518  pressure_dof_is_hanging[l] =
1519  pressure_node_pt(l)->is_hanging(p_index);
1520  }
1521  }
1522  // Otherwise the pressure is not stored at a node and so cannot hang
1523  else
1524  {
1525  for(unsigned l=0;l<n_pres;++l) { pressure_dof_is_hanging[l] = false; }
1526  }
1527 
1528  // Set up memory for the shape and test functions
1529  // Note that there are only two dimensions, r and z, in this problem
1530  Shape psif(n_node), testf(n_node);
1531  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
1532 
1533  // Set up memory for pressure shape and test functions
1534  Shape psip(n_pres), testp(n_pres);
1535 
1536  // Determine number of shape controlling nodes
1537  const unsigned n_shape_controlling_node = nshape_controlling_nodes();
1538 
1539  // Deriatives of shape fct derivatives w.r.t. nodal coords
1540  RankFourTensor<double> d_dpsifdx_dX(2,n_shape_controlling_node,n_node,2);
1541  RankFourTensor<double> d_dtestfdx_dX(2,n_shape_controlling_node,n_node,2);
1542 
1543  // Derivative of Jacobian of mapping w.r.t. to nodal coords
1544  DenseMatrix<double> dJ_dX(2,n_shape_controlling_node);
1545 
1546  // Derivatives of derivative of u w.r.t. nodal coords
1547  // Note that the entry d_dudx_dX(p,q,i,k) contains the derivative w.r.t.
1548  // nodal coordinate X_pq of du_i/dx_k. Since there are three velocity
1549  // components, i can take on the values 0, 1 and 2, while k and p only
1550  // take on the values 0 and 1 since there are only two spatial dimensions.
1551  RankFourTensor<double> d_dudx_dX(2,n_shape_controlling_node,3,2);
1552 
1553  // Derivatives of nodal velocities w.r.t. nodal coords:
1554  // Assumption: Interaction only local via no-slip so
1555  // X_pq only affects U_iq.
1556  // Note that the entry d_U_dX(p,q,i) contains the derivative w.r.t. nodal
1557  // coordinate X_pq of nodal value U_iq. In other words this entry is the
1558  // derivative of the i-th velocity component w.r.t. the p-th spatial
1559  // coordinate, taken at the q-th local node.
1560  RankThreeTensor<double> d_U_dX(2,n_shape_controlling_node,3,0.0);
1561 
1562  // Determine the number of integration points
1563  const unsigned n_intpt = integral_pt()->nweight();
1564 
1565  // Vector to hold local coordinates (two dimensions)
1566  Vector<double> s(2);
1567 
1568  // Get physical variables from element
1569  // (Reynolds number must be multiplied by the density ratio)
1570  const double scaled_re = this->re()*this->density_ratio();
1571  const double scaled_re_st = this->re_st()*this->density_ratio();
1572  const double scaled_re_inv_fr = this->re_invfr()*this->density_ratio();
1573  const double scaled_re_inv_ro = this->re_invro()*this->density_ratio();
1574  const double visc_ratio = this->viscosity_ratio();
1575  const Vector<double> G = this->g();
1576 
1577  // FD step
1579 
1580  // Pre-compute derivatives of nodal velocities w.r.t. nodal coords:
1581  // Assumption: Interaction only local via no-slip so
1582  // X_ij only affects U_ij.
1583  bool element_has_node_with_aux_node_update_fct=false;
1584 
1585  std::map<Node*,unsigned> local_shape_controlling_node_lookup=
1587 
1588  // FD loop over shape-controlling nodes
1589  for(std::map<Node*,unsigned>::iterator it=
1590  local_shape_controlling_node_lookup.begin();
1591  it!=local_shape_controlling_node_lookup.end();
1592  it++)
1593  {
1594  // Get pointer to q-th local shape-controlling node
1595  Node* nod_pt=it->first;
1596 
1597  // Get its number
1598  unsigned q=it->second;
1599 
1600  // Only compute if there's a node-update fct involved
1601  if(nod_pt->has_auxiliary_node_update_fct_pt())
1602  {
1603  element_has_node_with_aux_node_update_fct=true;
1604 
1605  // This functionality has not been tested so issue a warning
1606  // to this effect
1607  std::ostringstream warning_stream;
1608  warning_stream << "\nThe functionality to evaluate the additional"
1609  << "\ncontribution to the deriv of the residual eqn"
1610  << "\nw.r.t. the nodal coordinates which comes about"
1611  << "\nif a node's values are updated using an auxiliary"
1612  << "\nnode update function has NOT been tested for"
1613  << "\nrefineable axisymmetric Navier-Stokes elements."
1614  << "\nUse at your own risk"
1615  << std::endl;
1617  warning_stream.str(),
1618  "RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations::get_dresidual_dnodal_coordinates",
1619  OOMPH_EXCEPTION_LOCATION);
1620 
1621  // Current nodal velocity
1622  Vector<double> u_ref(3);
1623  for(unsigned i=0;i<3;i++)
1624  {
1625  u_ref[i]=*(nod_pt->value_pt(u_nodal_index[i]));
1626  }
1627 
1628  // FD
1629  for(unsigned p=0;p<2;p++)
1630  {
1631  // Make backup
1632  double backup=nod_pt->x(p);
1633 
1634  // Do FD step. No node update required as we're
1635  // attacking the coordinate directly...
1636  nod_pt->x(p)+=eps_fd;
1637 
1638  // Do auxiliary node update (to apply no slip)
1640 
1641  // Loop over velocity components
1642  for(unsigned i=0;i<3;i++)
1643  {
1644  // Evaluate
1645  d_U_dX(p,q,i)=(*(nod_pt->value_pt(u_nodal_index[p]))-u_ref[p])/eps_fd;
1646  }
1647 
1648  // Reset
1649  nod_pt->x(p)=backup;
1650 
1651  // Do auxiliary node update (to apply no slip)
1653  }
1654  }
1655  }
1656 
1657  // Integer to store the local equation number
1658  int local_eqn=0;
1659 
1660  // Pointers to hang info object
1661  HangInfo *hang_info_pt=0;
1662 
1663  // Loop over the integration points
1664  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1665  {
1666  // Assign values of s
1667  for(unsigned i=0;i<2;i++) { s[i] = integral_pt()->knot(ipt,i); }
1668 
1669  // Get the integral weight
1670  const double w = integral_pt()->weight(ipt);
1671 
1672  // Call the derivatives of the shape and test functions
1673  const double J = this->dshape_and_dtest_eulerian_at_knot_axi_nst(
1674  ipt,psif,dpsifdx,d_dpsifdx_dX,testf,dtestfdx,d_dtestfdx_dX,dJ_dX);
1675 
1676  // Call the pressure shape and test functions
1677  this->pshape_axi_nst(s,psip,testp);
1678 
1679  // Allocate storage for the position and the derivative of the
1680  // mesh positions w.r.t. time
1682  Vector<double> mesh_velocity(2,0.0);
1683 
1684  // Allocate storage for the pressure, velocity components and their
1685  // time and spatial derivatives
1686  double interpolated_p=0.0;
1687  Vector<double> interpolated_u(3,0.0);
1688  Vector<double> dudt(3,0.0);
1689  DenseMatrix<double> interpolated_dudx(3,2,0.0);
1690 
1691  // Calculate pressure at integration point
1692  for(unsigned l=0;l<n_pres;l++)
1693  {
1694  interpolated_p += this->p_axi_nst(l)*psip[l];
1695  }
1696 
1697  // Calculate velocities and derivatives at integration point
1698  // ---------------------------------------------------------
1699 
1700  // Loop over nodes
1701  for(unsigned l=0;l<n_node;l++)
1702  {
1703  // Cache the shape function
1704  const double psif_ = psif(l);
1705 
1706  // Loop over the two coordinate directions
1707  for(unsigned i=0;i<2;i++)
1708  {
1709  interpolated_x[i] += nodal_position(l,i)*psif_;
1710  }
1711 
1712  // Loop over the three velocity directions
1713  for(unsigned i=0;i<3;i++)
1714  {
1715  // Get the nodal value
1716  const double u_value = nodal_value(l,u_nodal_index[i]);
1717  interpolated_u[i] += u_value*psif_;
1718  dudt[i] += this->du_dt_axi_nst(l,i)*psif_;
1719 
1720  // Loop over derivative directions
1721  for(unsigned j=0;j<2;j++)
1722  {
1723  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
1724  }
1725  }
1726  }
1727 
1728  // Get the mesh velocity if ALE is enabled
1729  if(!this->ALE_is_disabled)
1730  {
1731  // Loop over nodes
1732  for(unsigned l=0;l<n_node;l++)
1733  {
1734  // Loop over the two coordinate directions
1735  for(unsigned i=0;i<2;i++)
1736  {
1737  mesh_velocity[i] += this->dnodal_position_dt(l,i)*psif[l];
1738  }
1739  }
1740  }
1741 
1742  // Calculate derivative of du_i/dx_k w.r.t. nodal positions X_{pq}
1743 
1744  // Loop over shape-controlling nodes
1745  for(unsigned q=0;q<n_shape_controlling_node;q++)
1746  {
1747  // Loop over the two coordinate directions
1748  for(unsigned p=0;p<2;p++)
1749  {
1750  // Loop over the three velocity components
1751  for(unsigned i=0;i<3;i++)
1752  {
1753  // Loop over the two coordinate directions
1754  for(unsigned k=0;k<2;k++)
1755  {
1756  double aux = 0.0;
1757 
1758  // Loop over nodes and add contribution
1759  for(unsigned j=0;j<n_node;j++)
1760  {
1761  aux += nodal_value(j,u_nodal_index[i])*d_dpsifdx_dX(p,q,j,k);
1762  }
1763  d_dudx_dX(p,q,i,k) = aux;
1764  }
1765  }
1766  }
1767  }
1768 
1769  // Get weight of actual nodal position/value in computation of mesh
1770  // velocity from positional/value time stepper
1771  const double pos_time_weight
1772  = node_pt(0)->position_time_stepper_pt()->weight(1,0);
1773  const double val_time_weight = node_pt(0)->time_stepper_pt()->weight(1,0);
1774 
1775  // Get the user-defined body force terms
1776  Vector<double> body_force(3);
1777  this->get_body_force_axi_nst(time,ipt,s,interpolated_x,body_force);
1778 
1779  // Get the user-defined source function
1780  const double source = this->get_source_fct(time,ipt,interpolated_x);
1781 
1782  // Get gradient of body force function
1783  DenseMatrix<double> d_body_force_dx(3,2,0.0);
1784  this->get_body_force_gradient_axi_nst(time,ipt,s,interpolated_x,
1785  d_body_force_dx);
1786 
1787  // Get gradient of source function
1788  Vector<double> source_gradient(2,0.0);
1789  this->get_source_fct_gradient(time,ipt,interpolated_x,
1790  source_gradient);
1791 
1792  // Cache r (the first position component)
1793  const double r = interpolated_x[0];
1794 
1795  // Assemble shape derivatives
1796  // --------------------------
1797 
1798  // ==================
1799  // MOMENTUM EQUATIONS
1800  // ==================
1801 
1802  // Number of master nodes and storage for the weight of the shape function
1803  unsigned n_master=1; double hang_weight=1.0;
1804 
1805  // Loop over the test functions
1806  for(unsigned l=0;l<n_node;l++)
1807  {
1808  // Cache the test function
1809  const double testf_ = testf[l];
1810 
1811  // Local boolean to indicate whether the node is hanging
1812  bool is_node_hanging = node_pt(l)->is_hanging();
1813 
1814  // If the node is hanging
1815  if(is_node_hanging)
1816  {
1817  hang_info_pt = node_pt(l)->hanging_pt();
1818 
1819  // Read out number of master nodes from hanging data
1820  n_master = hang_info_pt->nmaster();
1821  }
1822  // Otherwise the node is its own master
1823  else
1824  {
1825  n_master = 1;
1826  }
1827 
1828  // Loop over the master nodes
1829  for(unsigned m=0;m<n_master;m++)
1830  {
1831 
1832  // --------------------------------
1833  // FIRST (RADIAL) MOMENTUM EQUATION
1834  // --------------------------------
1835 
1836  // Get the equation number
1837  // If the node is hanging
1838  if(is_node_hanging)
1839  {
1840  // Get the equation number from the master node
1841  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1842  u_nodal_index[0]);
1843  // Get the hang weight from the master node
1844  hang_weight = hang_info_pt->master_weight(m);
1845  }
1846  // If the node is not hanging
1847  else
1848  {
1849  // Local equation number
1850  local_eqn = this->nodal_local_eqn(l,u_nodal_index[0]);
1851 
1852  // Node contributes with full weight
1853  hang_weight = 1.0;
1854  }
1855 
1856  // If it's not a boundary condition
1857  if(local_eqn >= 0)
1858  {
1859  // Loop over the two coordinate directions
1860  for(unsigned p=0;p<2;p++)
1861  {
1862  // Loop over shape controlling nodes
1863  for(unsigned q=0;q<n_shape_controlling_node;q++)
1864  {
1865 
1866  // Residual x deriv of Jacobian
1867  // ----------------------------
1868 
1869  // Add the user-defined body force terms
1870  double sum = r*body_force[0]*testf_;
1871 
1872  // Add the gravitational body force term
1873  sum += r*scaled_re_inv_fr*testf_*G[0];
1874 
1875  // Add the pressure gradient term
1876  sum += interpolated_p*(testf_ + r*dtestfdx(l,0));
1877 
1878  // Add the stress tensor terms
1879  // The viscosity ratio needs to go in here to ensure
1880  // continuity of normal stress is satisfied even in flows
1881  // with zero pressure gradient!
1882  sum -= visc_ratio*
1883  r*(1.0+Gamma[0])*interpolated_dudx(0,0)*dtestfdx(l,0);
1884 
1885  sum -= visc_ratio*r*
1886  (interpolated_dudx(0,1) + Gamma[0]*interpolated_dudx(1,0))*
1887  dtestfdx(l,1);
1888 
1889  sum -= visc_ratio*(1.0 + Gamma[0])*interpolated_u[0]*testf_/r;
1890 
1891  // Add the du/dt term
1892  sum -= scaled_re_st*r*dudt[0]*testf_;
1893 
1894  // Add the convective terms
1895  sum -=
1896  scaled_re*(r*interpolated_u[0]*interpolated_dudx(0,0)
1897  - interpolated_u[2]*interpolated_u[2]
1898  + r*interpolated_u[1]*interpolated_dudx(0,1))*testf_;
1899 
1900  // Add the mesh velocity terms
1901  if(!ALE_is_disabled)
1902  {
1903  for(unsigned k=0;k<2;k++)
1904  {
1905  sum += scaled_re_st*r*mesh_velocity[k]
1906  *interpolated_dudx(0,k)*testf_;
1907  }
1908  }
1909 
1910  // Add the Coriolis term
1911  sum += 2.0*r*scaled_re_inv_ro*interpolated_u[2]*testf_;
1912 
1913  // Multiply through by deriv of Jacobian and integration weight
1914  dresidual_dnodal_coordinates(local_eqn,p,q)
1915  += sum*dJ_dX(p,q)*w*hang_weight;
1916 
1917  // Derivative of residual x Jacobian
1918  // ---------------------------------
1919 
1920  // Body force
1921  sum = r*d_body_force_dx(0,p)*psif[q]*testf_;
1922  if(p==0) { sum += body_force[0]*psif[q]*testf_; }
1923 
1924  // Gravitational body force
1925  if(p==0) { sum += scaled_re_inv_fr*G[0]*psif[q]*testf_; }
1926 
1927  // Pressure gradient term
1928  sum += r*interpolated_p*d_dtestfdx_dX(p,q,l,0);
1929  if(p==0) { sum += interpolated_p*psif[q]*dtestfdx(l,0); }
1930 
1931  // Viscous terms
1932  sum -= r*visc_ratio*(
1933  (1.0+Gamma[0])*(
1934  d_dudx_dX(p,q,0,0)*dtestfdx(l,0)
1935  + interpolated_dudx(0,0)*d_dtestfdx_dX(p,q,l,0))
1936  + (d_dudx_dX(p,q,0,1)
1937  + Gamma[0]*d_dudx_dX(p,q,1,0))*dtestfdx(l,1)
1938  + (interpolated_dudx(0,1)
1939  + Gamma[0]*interpolated_dudx(1,0))
1940  *d_dtestfdx_dX(p,q,l,1));
1941  if(p==0)
1942  {
1943  sum -= visc_ratio*(
1944  (1.0+Gamma[0])*(
1945  interpolated_dudx(0,0)*psif[q]*dtestfdx(l,0)
1946  - interpolated_u[0]*psif[q]*testf_/(r*r))
1947  + (interpolated_dudx(0,1)
1948  + Gamma[0]*interpolated_dudx(1,0))*psif[q]*dtestfdx(l,1));
1949  }
1950 
1951  // Convective terms, including mesh velocity
1952  for(unsigned k=0;k<2;k++)
1953  {
1954  double tmp = scaled_re*interpolated_u[k];
1955  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
1956  sum -= r*tmp*d_dudx_dX(p,q,0,k)*testf_;
1957  if(p==0) { sum -= tmp*interpolated_dudx(0,k)*psif[q]*testf_; }
1958  }
1959  if(!ALE_is_disabled)
1960  {
1961  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(0,p)
1962  *psif[q]*testf_;
1963  }
1964 
1965  // du/dt term
1966  if(p==0) { sum -= scaled_re_st*dudt[0]*psif[q]*testf_; }
1967 
1968  // Coriolis term
1969  if(p==0)
1970  {
1971  sum += 2.0*scaled_re_inv_ro*interpolated_u[2]*psif[q]*testf_;
1972  }
1973 
1974  // Multiply through by Jacobian and integration weight
1975  dresidual_dnodal_coordinates(local_eqn,p,q)
1976  += sum*J*w*hang_weight;
1977 
1978  } // End of loop over shape controlling nodes q
1979  } // End of loop over coordinate directions p
1980 
1981  // Derivs w.r.t. to nodal velocities
1982  // ---------------------------------
1983  if(element_has_node_with_aux_node_update_fct)
1984  {
1985  // Loop over local nodes
1986  for (unsigned q_local=0;q_local<n_node;q_local++)
1987  {
1988  // Number of master nodes and storage for the weight of
1989  // the shape function
1990  unsigned n_master2=1;
1991  double hang_weight2=1.0;
1992  HangInfo* hang_info2_pt=0;
1993 
1994  // Local boolean to indicate whether the node is hanging
1995  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1996 
1997  Node* actual_shape_controlling_node_pt=node_pt(q_local);
1998 
1999  // If the node is hanging
2000  if(is_node_hanging2)
2001  {
2002  hang_info2_pt = node_pt(q_local)->hanging_pt();
2003 
2004  // Read out number of master nodes from hanging data
2005  n_master2 = hang_info2_pt->nmaster();
2006  }
2007  // Otherwise the node is its own master
2008  else
2009  {
2010  n_master2 = 1;
2011  }
2012 
2013  // Loop over the master nodes
2014  for(unsigned mm=0;mm<n_master2;mm++)
2015  {
2016  if(is_node_hanging2)
2017  {
2018  actual_shape_controlling_node_pt=
2019  hang_info2_pt->master_node_pt(mm);
2020  hang_weight2 = hang_info2_pt->master_weight(mm);
2021  }
2022 
2023  // Find the corresponding number
2024  unsigned q=local_shape_controlling_node_lookup[
2025  actual_shape_controlling_node_pt];
2026 
2027  // Loop over the two coordinate directions
2028  for(unsigned p=0;p<2;p++)
2029  {
2030  // Contribution from deriv of first velocity component
2031  double tmp = scaled_re_st*r*val_time_weight
2032  *psif[q_local]*testf_;
2033  tmp += r*scaled_re*interpolated_dudx(0,0)
2034  *psif[q_local]*testf_;
2035 
2036  for(unsigned k=0;k<2;k++)
2037  {
2038  double tmp2 = scaled_re*interpolated_u[k];
2039  if(!ALE_is_disabled)
2040  {
2041  tmp2 -= scaled_re_st*mesh_velocity[k];
2042  }
2043  tmp += r*tmp2*dpsifdx(q_local,k)*testf_;
2044  }
2045 
2046  tmp += r*visc_ratio*(1.0+Gamma[0])*dpsifdx(q_local,0)
2047  *dtestfdx(l,0);
2048  tmp += r*visc_ratio*dpsifdx(q_local,1)*dtestfdx(l,1);
2049  tmp += visc_ratio*(1.0+Gamma[0])*psif[q_local]*testf_/r;
2050 
2051  // Multiply through by dU_0q/dX_pq
2052  double sum = -d_U_dX(p,q_local,0)*tmp;
2053 
2054  // Contribution from deriv of second velocity component
2055  tmp = scaled_re*r*interpolated_dudx(0,1)*psif[q_local]*testf_;
2056  tmp += r*visc_ratio*Gamma[0]*dpsifdx(q_local,0)*dtestfdx(l,1);
2057 
2058  // Multiply through by dU_1q/dX_pq
2059  sum -= d_U_dX(p,q,1)*tmp;
2060 
2061  // Contribution from deriv of third velocity component
2062  tmp = 2.0*(r*scaled_re_inv_ro + scaled_re*interpolated_u[2])
2063  *psif[q_local]*testf_;
2064 
2065  // Multiply through by dU_2q/dX_pq
2066  sum += d_U_dX(p,q,2)*tmp;
2067 
2068  // Multiply through by Jacobian and integration weight
2069  dresidual_dnodal_coordinates(local_eqn,p,q) +=
2070  sum*J*w*hang_weight*hang_weight2;
2071  }
2072  } // End of loop over master nodes
2073  } // End of loop over local nodes
2074  } // End of if(element_has_node_with_aux_node_update_fct)
2075  } // End of if it's not a boundary condition
2076 
2077  // --------------------------------
2078  // SECOND (AXIAL) MOMENTUM EQUATION
2079  // --------------------------------
2080 
2081  // Get the equation number
2082  // If the node is hanging
2083  if(is_node_hanging)
2084  {
2085  // Get the equation number from the master node
2086  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
2087  u_nodal_index[1]);
2088  // Get the hang weight from the master node
2089  hang_weight = hang_info_pt->master_weight(m);
2090  }
2091  // If the node is not hanging
2092  else
2093  {
2094  // Local equation number
2095  local_eqn = this->nodal_local_eqn(l,u_nodal_index[1]);
2096 
2097  // Node contributes with full weight
2098  hang_weight = 1.0;
2099  }
2100 
2101  // If it's not a boundary condition
2102  if(local_eqn >= 0)
2103  {
2104  // Loop over the two coordinate directions
2105  for(unsigned p=0;p<2;p++)
2106  {
2107  // Loop over shape controlling nodes
2108  for(unsigned q=0;q<n_shape_controlling_node;q++)
2109  {
2110 
2111  // Residual x deriv of Jacobian
2112  // ----------------------------
2113 
2114  // Add the user-defined body force terms
2115  double sum = r*body_force[1]*testf_;
2116 
2117  // Add the gravitational body force term
2118  sum += r*scaled_re_inv_fr*testf_*G[1];
2119 
2120  // Add the pressure gradient term
2121  sum += r*interpolated_p*dtestfdx(l,1);
2122 
2123  // Add the stress tensor terms
2124  // The viscosity ratio needs to go in here to ensure
2125  // continuity of normal stress is satisfied even in flows
2126  // with zero pressure gradient!
2127  sum -= visc_ratio*
2128  r*(interpolated_dudx(1,0) + Gamma[1]*interpolated_dudx(0,1))
2129  *dtestfdx(l,0);
2130 
2131  sum -= visc_ratio*r*
2132  (1.0 + Gamma[1])*interpolated_dudx(1,1)*dtestfdx(l,1);
2133 
2134  // Add the du/dt term
2135  sum -= scaled_re_st*r*dudt[1]*testf_;
2136 
2137  // Add the convective terms
2138  sum -=
2139  scaled_re*(r*interpolated_u[0]*interpolated_dudx(1,0)
2140  + r*interpolated_u[1]*interpolated_dudx(1,1))*testf_;
2141 
2142  // Add the mesh velocity terms
2143  if(!ALE_is_disabled)
2144  {
2145  for(unsigned k=0;k<2;k++)
2146  {
2147  sum += scaled_re_st*r*mesh_velocity[k]
2148  *interpolated_dudx(1,k)*testf_;
2149  }
2150  }
2151 
2152  // Multiply through by deriv of Jacobian and integration weight
2153  dresidual_dnodal_coordinates(local_eqn,p,q)
2154  += sum*dJ_dX(p,q)*w*hang_weight;
2155 
2156  // Derivative of residual x Jacobian
2157  // ---------------------------------
2158 
2159  // Body force
2160  sum = r*d_body_force_dx(1,p)*psif[q]*testf_;
2161  if(p==0) { sum += body_force[1]*psif[q]*testf_; }
2162 
2163  // Gravitational body force
2164  if(p==0) { sum += scaled_re_inv_fr*G[1]*psif[q]*testf_; }
2165 
2166  // Pressure gradient term
2167  sum += r*interpolated_p*d_dtestfdx_dX(p,q,l,1);
2168  if(p==0) { sum += interpolated_p*psif[q]*dtestfdx(l,1); }
2169 
2170  // Viscous terms
2171  sum -= r*visc_ratio*(
2172  (d_dudx_dX(p,q,1,0) + Gamma[1]*d_dudx_dX(p,q,0,1))*dtestfdx(l,0)
2173  + (interpolated_dudx(1,0)
2174  + Gamma[1]*interpolated_dudx(0,1))*d_dtestfdx_dX(p,q,l,0)
2175  + (1.0+Gamma[1])*d_dudx_dX(p,q,1,1)*dtestfdx(l,1)
2176  + (1.0+Gamma[1])*interpolated_dudx(1,1)*d_dtestfdx_dX(p,q,l,1));
2177  if(p==0)
2178  {
2179  sum -= visc_ratio*(
2180  (interpolated_dudx(1,0)
2181  + Gamma[1]*interpolated_dudx(0,1))*psif[q]*dtestfdx(l,0)
2182  + (1.0+Gamma[1])*interpolated_dudx(1,1)*psif[q]*dtestfdx(l,1));
2183  }
2184 
2185  // Convective terms, including mesh velocity
2186  for(unsigned k=0;k<2;k++)
2187  {
2188  double tmp = scaled_re*interpolated_u[k];
2189  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
2190  sum -= r*tmp*d_dudx_dX(p,q,1,k)*testf_;
2191  if(p==0) { sum -= tmp*interpolated_dudx(1,k)*psif[q]*testf_; }
2192  }
2193  if(!ALE_is_disabled)
2194  {
2195  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(1,p)
2196  *psif[q]*testf_;
2197  }
2198 
2199  // du/dt term
2200  if(p==0) { sum -= scaled_re_st*dudt[1]*psif[q]*testf_; }
2201 
2202  // Multiply through by Jacobian and integration weight
2203  dresidual_dnodal_coordinates(local_eqn,p,q)
2204  += sum*J*w*hang_weight;
2205 
2206  } // End of loop over shape controlling nodes q
2207  } // End of loop over coordinate directions p
2208 
2209  // Derivs w.r.t. to nodal velocities
2210  // ---------------------------------
2211  if(element_has_node_with_aux_node_update_fct)
2212  {
2213  // Loop over local nodes
2214  for (unsigned q_local=0;q_local<n_node;q_local++)
2215  {
2216  // Number of master nodes and storage for the weight of
2217  // the shape function
2218  unsigned n_master2=1;
2219  double hang_weight2=1.0;
2220  HangInfo* hang_info2_pt=0;
2221 
2222  // Local boolean to indicate whether the node is hanging
2223  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
2224 
2225  Node* actual_shape_controlling_node_pt=node_pt(q_local);
2226 
2227  // If the node is hanging
2228  if(is_node_hanging2)
2229  {
2230  hang_info2_pt = node_pt(q_local)->hanging_pt();
2231 
2232  // Read out number of master nodes from hanging data
2233  n_master2 = hang_info2_pt->nmaster();
2234  }
2235  // Otherwise the node is its own master
2236  else
2237  {
2238  n_master2 = 1;
2239  }
2240 
2241  // Loop over the master nodes
2242  for(unsigned mm=0;mm<n_master2;mm++)
2243  {
2244  if(is_node_hanging2)
2245  {
2246  actual_shape_controlling_node_pt=
2247  hang_info2_pt->master_node_pt(mm);
2248  hang_weight2 = hang_info2_pt->master_weight(mm);
2249  }
2250 
2251  // Find the corresponding number
2252  unsigned q=local_shape_controlling_node_lookup[
2253  actual_shape_controlling_node_pt];
2254 
2255  // Loop over the two coordinate directions
2256  for(unsigned p=0;p<2;p++)
2257  {
2258  // Contribution from deriv of first velocity component
2259  double tmp = scaled_re*r*interpolated_dudx(1,0)
2260  *psif[q_local]*testf_;
2261  tmp += r*visc_ratio*Gamma[1]*dpsifdx(q_local,1)*dtestfdx(l,0);
2262 
2263  // Multiply through by dU_0q/dX_pq
2264  double sum = -d_U_dX(p,q,0)*tmp;
2265 
2266  // Contribution from deriv of second velocity component
2267  tmp = scaled_re_st*r*val_time_weight*psif[q_local]*testf_;
2268  tmp += r*scaled_re*interpolated_dudx(1,1)
2269  *psif[q_local]*testf_;
2270 
2271  for(unsigned k=0;k<2;k++)
2272  {
2273  double tmp2 = scaled_re*interpolated_u[k];
2274  if(!ALE_is_disabled)
2275  {
2276  tmp2 -= scaled_re_st*mesh_velocity[k];
2277  }
2278  tmp += r*tmp2*dpsifdx(q_local,k)*testf_;
2279  }
2280  tmp += r*visc_ratio*(dpsifdx(q_local,0)*dtestfdx(l,0)
2281  + (1.0+Gamma[1])*dpsifdx(q_local,1)
2282  *dtestfdx(l,1));
2283 
2284  // Multiply through by dU_1q/dX_pq
2285  sum -= d_U_dX(p,q,1)*tmp;
2286 
2287  // Multiply through by Jacobian and integration weight
2288  dresidual_dnodal_coordinates(local_eqn,p,q) +=
2289  sum*J*w*hang_weight*hang_weight2;
2290  }
2291  } // End of loop over master nodes
2292  } // End of loop over local nodes
2293  } // End of if(element_has_node_with_aux_node_update_fct)
2294  } // End of if it's not a boundary condition
2295 
2296  // -----------------------------------
2297  // THIRD (AZIMUTHAL) MOMENTUM EQUATION
2298  // -----------------------------------
2299 
2300  // Get the equation number
2301  // If the node is hanging
2302  if(is_node_hanging)
2303  {
2304  // Get the equation number from the master node
2305  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
2306  u_nodal_index[2]);
2307  // Get the hang weight from the master node
2308  hang_weight = hang_info_pt->master_weight(m);
2309  }
2310  // If the node is not hanging
2311  else
2312  {
2313  // Local equation number
2314  local_eqn = this->nodal_local_eqn(l,u_nodal_index[2]);
2315 
2316  // Node contributes with full weight
2317  hang_weight = 1.0;
2318  }
2319 
2320  // If it's not a boundary condition
2321  if(local_eqn >= 0)
2322  {
2323  // Loop over the two coordinate directions
2324  for(unsigned p=0;p<2;p++)
2325  {
2326  // Loop over shape controlling nodes
2327  for(unsigned q=0;q<n_shape_controlling_node;q++)
2328  {
2329 
2330  // Residual x deriv of Jacobian
2331  // ----------------------------
2332 
2333  // Add the user-defined body force terms
2334  double sum = r*body_force[2]*testf_;
2335 
2336  // Add the gravitational body force term
2337  sum += r*scaled_re_inv_fr*testf_*G[2];
2338 
2339  // There is NO pressure gradient term
2340 
2341  // Add in the stress tensor terms
2342  // The viscosity ratio needs to go in here to ensure
2343  // continuity of normal stress is satisfied even in flows
2344  // with zero pressure gradient!
2345  sum -= visc_ratio*(r*interpolated_dudx(2,0) -
2346  Gamma[0]*interpolated_u[2])*dtestfdx(l,0);
2347 
2348  sum -= visc_ratio*r*interpolated_dudx(2,1)*dtestfdx(l,1);
2349 
2350  sum -= visc_ratio*((interpolated_u[2]/r)
2351  - Gamma[0]*interpolated_dudx(2,0))*testf_;
2352 
2353  // Add the du/dt term
2354  sum -= scaled_re_st*r*dudt[2]*testf_;
2355 
2356  // Add the convective terms
2357  sum -=
2358  scaled_re*(r*interpolated_u[0]*interpolated_dudx(2,0)
2359  + interpolated_u[0]*interpolated_u[2]
2360  + r*interpolated_u[1]*interpolated_dudx(2,1))*testf_;
2361 
2362  // Add the mesh velocity terms
2363  if(!ALE_is_disabled)
2364  {
2365  for(unsigned k=0;k<2;k++)
2366  {
2367  sum += scaled_re_st*r*mesh_velocity[k]
2368  *interpolated_dudx(2,k)*testf_;
2369  }
2370  }
2371 
2372  // Add the Coriolis term
2373  sum -= 2.0*r*scaled_re_inv_ro*interpolated_u[0]*testf_;
2374 
2375  // Multiply through by deriv of Jacobian and integration weight
2376  dresidual_dnodal_coordinates(local_eqn,p,q)
2377  += sum*dJ_dX(p,q)*w*hang_weight;
2378 
2379  // Derivative of residual x Jacobian
2380  // ---------------------------------
2381 
2382  // Body force
2383  sum = r*d_body_force_dx(2,p)*psif[q]*testf_;
2384  if(p==0) { sum += body_force[2]*psif[q]*testf_; }
2385 
2386  // Gravitational body force
2387  if(p==0) { sum += scaled_re_inv_fr*G[2]*psif[q]*testf_; }
2388 
2389  // There is NO pressure gradient term
2390 
2391  // Viscous terms
2392  sum -= visc_ratio*r*(
2393  d_dudx_dX(p,q,2,0)*dtestfdx(l,0)
2394  + interpolated_dudx(2,0)*d_dtestfdx_dX(p,q,l,0)
2395  + d_dudx_dX(p,q,2,1)*dtestfdx(l,1)
2396  + interpolated_dudx(2,1)*d_dtestfdx_dX(p,q,l,1));
2397 
2398  sum += visc_ratio*Gamma[0]*d_dudx_dX(p,q,2,0)*testf_;
2399 
2400  if(p==0)
2401  {
2402  sum -=
2403  visc_ratio*(interpolated_dudx(2,0)*psif[q]*dtestfdx(l,0)
2404  + interpolated_dudx(2,1)*psif[q]*dtestfdx(l,1)
2405  + interpolated_u[2]*psif[q]*testf_/(r*r));
2406  }
2407 
2408  // Convective terms, including mesh velocity
2409  for(unsigned k=0;k<2;k++)
2410  {
2411  double tmp = scaled_re*interpolated_u[k];
2412  if(!ALE_is_disabled) { tmp -= scaled_re_st*mesh_velocity[k]; }
2413  sum -= r*tmp*d_dudx_dX(p,q,2,k)*testf_;
2414  if(p==0) { sum -= tmp*interpolated_dudx(2,k)*psif[q]*testf_; }
2415  }
2416  if(!ALE_is_disabled)
2417  {
2418  sum += scaled_re_st*r*pos_time_weight*interpolated_dudx(2,p)
2419  *psif[q]*testf_;
2420  }
2421 
2422  // du/dt term
2423  if(p==0) { sum -= scaled_re_st*dudt[2]*psif[q]*testf_; }
2424 
2425  // Coriolis term
2426  if(p==0)
2427  {
2428  sum -= 2.0*scaled_re_inv_ro*interpolated_u[0]*psif[q]*testf_;
2429  }
2430 
2431  // Multiply through by Jacobian and integration weight
2432  dresidual_dnodal_coordinates(local_eqn,p,q)
2433  += sum*J*w*hang_weight;
2434 
2435  } // End of loop over shape controlling nodes q
2436  } // End of loop over coordinate directions p
2437 
2438  // Derivs w.r.t. to nodal velocities
2439  // ---------------------------------
2440  if(element_has_node_with_aux_node_update_fct)
2441  {
2442  // Loop over local nodes
2443  for (unsigned q_local=0;q_local<n_node;q_local++)
2444  {
2445  // Number of master nodes and storage for the weight of
2446  // the shape function
2447  unsigned n_master2=1;
2448  double hang_weight2=1.0;
2449  HangInfo* hang_info2_pt=0;
2450 
2451  // Local boolean to indicate whether the node is hanging
2452  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
2453 
2454  Node* actual_shape_controlling_node_pt=node_pt(q_local);
2455 
2456  // If the node is hanging
2457  if(is_node_hanging2)
2458  {
2459  hang_info2_pt = node_pt(q_local)->hanging_pt();
2460 
2461  // Read out number of master nodes from hanging data
2462  n_master2 = hang_info2_pt->nmaster();
2463  }
2464  // Otherwise the node is its own master
2465  else
2466  {
2467  n_master2 = 1;
2468  }
2469 
2470  // Loop over the master nodes
2471  for(unsigned mm=0;mm<n_master2;mm++)
2472  {
2473  if(is_node_hanging2)
2474  {
2475  actual_shape_controlling_node_pt=
2476  hang_info2_pt->master_node_pt(mm);
2477  hang_weight2 = hang_info2_pt->master_weight(mm);
2478  }
2479 
2480  // Find the corresponding number
2481  unsigned q=local_shape_controlling_node_lookup[
2482  actual_shape_controlling_node_pt];
2483 
2484  // Loop over the two coordinate directions
2485  for(unsigned p=0;p<2;p++)
2486  {
2487  // Contribution from deriv of first velocity component
2488  double tmp = (2.0*r*scaled_re_inv_ro
2489  + r*scaled_re*interpolated_dudx(2,0)
2490  - scaled_re*interpolated_u[2])
2491  *psif[q_local]*testf_;
2492 
2493  // Multiply through by dU_0q/dX_pq
2494  double sum = -d_U_dX(p,q,0)*tmp;
2495 
2496  // Contribution from deriv of second velocity component
2497  // Multiply through by dU_1q/dX_pq
2498  sum -= d_U_dX(p,q,1)*scaled_re*r*interpolated_dudx(2,1)
2499  *psif[q_local]*testf_;
2500 
2501  // Contribution from deriv of third velocity component
2502  tmp = scaled_re_st*r*val_time_weight*psif[q_local]*testf_;
2503  tmp -= scaled_re*interpolated_u[0]*psif[q_local]*testf_;
2504 
2505  for(unsigned k=0;k<2;k++)
2506  {
2507  double tmp2 = scaled_re*interpolated_u[k];
2508  if(!ALE_is_disabled)
2509  {
2510  tmp2 -= scaled_re_st*mesh_velocity[k];
2511  }
2512  tmp += r*tmp2*dpsifdx(q_local,k)*testf_;
2513  }
2514  tmp += visc_ratio*(r*dpsifdx(q_local,0)
2515  - Gamma[0]*psif[q_local])*dtestfdx(l,0);
2516  tmp += r*visc_ratio*dpsifdx(q_local,1)*dtestfdx(l,1);
2517  tmp += visc_ratio*((psif[q_local]/r)
2518  - Gamma[0]*dpsifdx(q_local,0))*testf_;
2519 
2520  // Multiply through by dU_2q/dX_pq
2521  sum -= d_U_dX(p,q,2)*tmp;
2522 
2523  // Multiply through by Jacobian and integration weight
2524  dresidual_dnodal_coordinates(local_eqn,p,q) +=
2525  sum*J*w*hang_weight*hang_weight2;
2526  }
2527  } // End of loop over master nodes
2528  } // End of loop over local nodes
2529  } // End of if(element_has_node_with_aux_node_update_fct)
2530  } // End of if it's not a boundary condition
2531  } // End of loop over master nodes
2532  } // End of loop over test functions
2533 
2534 
2535  // ===================
2536  // CONTINUITY EQUATION
2537  // ===================
2538 
2539  // Loop over the shape functions
2540  for(unsigned l=0;l<n_pres;l++)
2541  {
2542 
2543  // If the pressure dof is hanging
2544  if(pressure_dof_is_hanging[l])
2545  {
2546  // Pressure dof is hanging so it must be nodal-based
2547  // Get the hang info object
2548  hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
2549 
2550  // Get the number of master nodes from the pressure node
2551  n_master = hang_info_pt->nmaster();
2552  }
2553  // Otherwise the node is its own master
2554  else
2555  {
2556  n_master = 1;
2557  }
2558 
2559  // Loop over the master nodes
2560  for(unsigned m=0;m<n_master;m++)
2561  {
2562  // Get the number of the unknown
2563  // If the pressure dof is hanging
2564  if(pressure_dof_is_hanging[l])
2565  {
2566  // Get the local equation from the master node
2567  local_eqn =
2568  this->local_hang_eqn(hang_info_pt->master_node_pt(m),p_index);
2569  // Get the weight for the node
2570  hang_weight = hang_info_pt->master_weight(m);
2571  }
2572  else
2573  {
2574  local_eqn = this->p_local_eqn(l);
2575  hang_weight = 1.0;
2576  }
2577 
2578  // Cache the test function
2579  const double testp_ = testp[l];
2580 
2581  // If not a boundary conditions
2582  if(local_eqn >= 0)
2583  {
2584 
2585  // Loop over the two coordinate directions
2586  for(unsigned p=0;p<2;p++)
2587  {
2588  // Loop over shape controlling nodes
2589  for(unsigned q=0;q<n_shape_controlling_node;q++)
2590  {
2591 
2592  // Residual x deriv of Jacobian
2593  //-----------------------------
2594 
2595  // Source term
2596  double aux = -r*source;
2597 
2598  // Gradient terms
2599  aux += (interpolated_u[0]
2600  + r*interpolated_dudx(0,0)
2601  + r*interpolated_dudx(1,1));
2602 
2603  // Multiply through by deriv of Jacobian and integration weight
2604  dresidual_dnodal_coordinates(local_eqn,p,q)
2605  += aux*dJ_dX(p,q)*testp_*w*hang_weight;
2606 
2607  // Derivative of residual x Jacobian
2608  // ---------------------------------
2609 
2610  // Gradient of source function
2611  aux = -r*source_gradient[p]*psif[q];
2612  if(p==0) { aux -= source*psif[q]; }
2613 
2614  // Gradient terms
2615  aux += r*(d_dudx_dX(p,q,0,0) + d_dudx_dX(p,q,1,1));
2616  if(p==0)
2617  {
2618  aux += (interpolated_dudx(0,0)
2619  + interpolated_dudx(1,1))*psif[q];
2620  }
2621 
2622  // Multiply through by Jacobian and integration weight
2623  dresidual_dnodal_coordinates(local_eqn,p,q)
2624  += aux*testp_*J*w*hang_weight;
2625  }
2626  }
2627 
2628  // Derivs w.r.t. to nodal velocities
2629  // ---------------------------------
2630  if(element_has_node_with_aux_node_update_fct)
2631  {
2632  // Loop over local nodes
2633  for(unsigned q_local=0;q_local<n_node;q_local++)
2634  {
2635  // Number of master nodes and storage for the weight of
2636  // the shape function
2637  unsigned n_master2=1;
2638  double hang_weight2=1.0;
2639  HangInfo* hang_info2_pt=0;
2640 
2641  // Local boolean to indicate whether the node is hanging
2642  bool is_node_hanging2 = node_pt(q_local)->is_hanging();
2643 
2644  Node* actual_shape_controlling_node_pt=node_pt(q_local);
2645 
2646  // If the node is hanging
2647  if(is_node_hanging2)
2648  {
2649  hang_info2_pt = node_pt(q_local)->hanging_pt();
2650 
2651  // Read out number of master nodes from hanging data
2652  n_master2 = hang_info2_pt->nmaster();
2653  }
2654  // Otherwise the node is its own master
2655  else
2656  {
2657  n_master2 = 1;
2658  }
2659 
2660  // Loop over the master nodes
2661  for(unsigned mm=0;mm<n_master2;mm++)
2662  {
2663 
2664  if(is_node_hanging2)
2665  {
2666  actual_shape_controlling_node_pt=
2667  hang_info2_pt->master_node_pt(mm);
2668  hang_weight2 = hang_info2_pt->master_weight(mm);
2669  }
2670 
2671  // Find the corresponding number
2672  unsigned q=local_shape_controlling_node_lookup[
2673  actual_shape_controlling_node_pt];
2674 
2675  // Loop over the two coordinate directions
2676  for(unsigned p=0;p<2;p++)
2677  {
2678  double aux =
2679  d_U_dX(p,q,0)*(psif[q_local] + r*dpsifdx(q_local,0))
2680  + d_U_dX(p,q,1)*r*dpsifdx(q_local,1);
2681 
2682  // Multiply through by Jacobian, test function and
2683  // integration weight
2684  dresidual_dnodal_coordinates(local_eqn,p,q) +=
2685  aux*testp_*J*w*hang_weight*hang_weight2;
2686  }
2687  } // End of loop over (mm) master nodes
2688  } // End of loop over local nodes q_local
2689  } // End of if(element_has_node_with_aux_node_update_fct)
2690  } // End of if it's not a boundary condition
2691  } // End of loop over (m) master nodes
2692  } // End of loop over shape functions for continuity eqn
2693 
2694  } // End of loop over integration points
2695 
2696 } // End of get_dresidual_dnodal_coordinates(...)
2697 
2698 
2699 
2700 } // End of namespace
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 void extrapolated_strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Extrapolated strain-rate tensor: where (in that order) based on the previous time steps evaluated a...
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}.
GeneralisedNewtonianConstitutiveEquation< 3 > * Constitutive_eqn_pt
Pointer to the generalised Newtonian constitutive equation.
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
virtual double dviscosity_dinvariant(const double &second_invariant_of_rate_of_strain_tensor)=0
cstr elem_len * i
Definition: cfortran.h:607
virtual unsigned u_index_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.
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
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
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
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
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
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...
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...
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
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...
A Rank 4 Tensor class.
Definition: matrices.h:1625
const Vector< double > & g() const
Vector of gravitational components.
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
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.
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
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
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
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.
virtual int p_nodal_index_axi_nst() const
Which nodal value represents the pressure?
double get_source_fct(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
virtual void get_source_fct_gradient(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
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
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...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
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...
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
virtual void pshape_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
Class that contains data for hanging nodes.
Definition: nodes.h:684
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
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
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.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed...
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.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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...
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
virtual unsigned npres_axi_nst() const =0
Function to return number of pressure degrees of freedom.
virtual double viscosity(const double &second_invariant_of_rate_of_strain_tensor)=0
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
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...
double second_invariant(const DenseMatrix< double > &tensor)
Compute second invariant of tensor.
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164