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-2022 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//====================================================================
27
28
29namespace 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 //=======================================================================
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
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);
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
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
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
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 (
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 //======================================================================
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
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 =
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
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)
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)
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 =
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
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 =
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);
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
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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
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
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
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1198
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
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...
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,...
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...
virtual void extrapolated_strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Extrapolated strain-rate tensor: where (in that order) based on the previous time steps evaluated a...
virtual void pshape_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
virtual unsigned npres_axi_nst() const =0
Function to return number of pressure degrees of freedom.
GeneralisedNewtonianConstitutiveEquation< 3 > * Constitutive_eqn_pt
Pointer to the generalised Newtonian constitutive equation.
double du_dt_axi_nst(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes.
double get_source_fct(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
virtual void get_body_force_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Calculate the body force fct at a given time and Eulerian position.
virtual int p_nodal_index_axi_nst() const
Which nodal value represents the pressure?
const Vector< double > & g() const
Vector of gravitational components.
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
virtual double p_axi_nst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes.
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
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...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed....
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
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...
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
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
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
bool has_auxiliary_node_update_fct_pt()
Boolean to indicate if node has a pointer to and auxiliary update function.
Definition: nodes.h:1606
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
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
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1701
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
std::map< Node *, unsigned > shape_controlling_node_lookup()
Return lookup scheme for unique number associated with any of the nodes that actively control the sha...
unsigned nshape_controlling_nodes()
Number of shape-controlling nodes = the number of non-hanging nodes plus the number of master nodes a...
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 ...
void fill_in_generic_residual_contribution_axi_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add element's contribution to the elemental residual vector and/or Jacobian matrix and mass matrix fl...
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
virtual Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interp...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
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
Time *const & time_pt() const
Access function for the pointer to time (const version)
Definition: timesteppers.h:572
double & time()
Return the current value of the continuous time.
Definition: timesteppers.h:123
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
double second_invariant(const DenseMatrix< double > &tensor)
Compute second invariant of tensor.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...