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