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