refineable_linearised_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-2025 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Non-inline functions for the refineable linearised axisymmetric
27// Navier-Stokes elements
28
29// oomph-lib includes
30#include "refineable_linearised_axisym_navier_stokes_elements.h"
31
32namespace oomph
33{
34 //=======================================================================
35 /// Compute the residuals for the refineable linearised axisymmetric
36 /// Navier--Stokes equations; flag=1(or 0): do (or don't) compute the
37 /// Jacobian as well.
38 //=======================================================================
39 void RefineableLinearisedAxisymmetricNavierStokesEquations::
40 fill_in_generic_residual_contribution_linearised_axi_nst(
41 Vector<double>& residuals,
42 DenseMatrix<double>& jacobian,
43 DenseMatrix<double>& mass_matrix,
44 unsigned flag)
45 {
46 // Get the time from the first node in the element
47 const double time = this->node_pt(0)->time_stepper_pt()->time();
48
49 // Determine number of nodes in the element
50 const unsigned n_node = nnode();
51
52 // Determine how many pressure values there are associated with
53 // a single pressure component
54 const unsigned n_pres = npres_linearised_axi_nst();
55
56 // Get the nodal indices at which the velocity is stored
57 unsigned u_nodal_index[6];
58 for (unsigned i = 0; i < 6; ++i)
59 {
60 u_nodal_index[i] = u_index_linearised_axi_nst(i);
61 }
62
63 // Which nodal values represent the two pressure components?
64 // (Negative if pressure is not based on nodal interpolation).
65 Vector<int> p_index(2);
66 for (unsigned i = 0; i < 2; i++)
67 {
68 p_index[i] = this->p_index_linearised_axi_nst(i);
69 }
70
71 // Local array of booleans that are true if the l-th pressure value is
72 // hanging (avoid repeated virtual function calls)
73 bool pressure_dof_is_hanging[n_pres];
74
75 // If the pressure is stored at a node
76 if (p_index[0] >= 0)
77 {
78 // Read out whether the pressure is hanging
79 for (unsigned l = 0; l < n_pres; ++l)
80 {
81 pressure_dof_is_hanging[l] =
82 pressure_node_pt(l)->is_hanging(p_index[0]);
83 }
84 }
85 // Otherwise the pressure is not stored at a node and so cannot hang
86 else
87 {
88 for (unsigned l = 0; l < n_pres; ++l)
89 {
90 pressure_dof_is_hanging[l] = false;
91 }
92 }
93
94 // Set up memory for the fluid shape and test functions
95 // Note that there are two dimensions, r and z, in this problem
96 Shape psif(n_node), testf(n_node);
97 DShape dpsifdx(n_node, 2), dtestfdx(n_node, 2);
98
99 // Set up memory for the pressure shape and test functions
100 Shape psip(n_pres), testp(n_pres);
101
102 // Determine number of integration points
103 const unsigned n_intpt = integral_pt()->nweight();
104
105 // Set up memory for the vector to hold local coordinates (two dimensions)
106 Vector<double> s(2);
107
108 // Get physical variables from the element
109 // (Reynolds number must be multiplied by the density ratio)
110 const double scaled_re = re() * density_ratio();
111 const double scaled_re_st = re_st() * density_ratio();
112 const double visc_ratio = viscosity_ratio();
113 const int k = azimuthal_mode_number();
114
115 // Integers used to store the local equation and unknown numbers
116 int local_eqn = 0, local_unknown = 0;
117
118 // Local storage for pointers to hang info objects
119 HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
120
121 // Loop over the integration points
122 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
123 {
124 // Assign values of the local coordinates s
125 for (unsigned i = 0; i < 2; i++)
126 {
127 s[i] = integral_pt()->knot(ipt, i);
128 }
129
130 // Get the integral weight
131 const double w = integral_pt()->weight(ipt);
132
133 // Calculate the fluid shape and test functions, and their derivatives
134 // w.r.t. the global coordinates
135 const double J = dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(
136 ipt, psif, dpsifdx, testf, dtestfdx);
137
138 // Calculate the pressure shape and test functions
139 pshape_linearised_axi_nst(s, psip, testp);
140
141 // Premultiply the weights and the Jacobian of the mapping between
142 // local and global coordinates
143 const double W = w * J;
144
145 // Allocate storage for the position and the derivative of the
146 // mesh positions w.r.t. time
147 Vector<double> interpolated_x(2, 0.0);
148 Vector<double> mesh_velocity(2, 0.0);
149
150 // Allocate storage for the velocity components (six of these)
151 // and their derivatives w.r.t. time
152 Vector<double> interpolated_u(6, 0.0);
153 Vector<double> dudt(6, 0.0);
154
155 // Allocate storage for the pressure components (two of these)
156 Vector<double> interpolated_p(2, 0.0);
157
158 // Allocate storage for the derivatives of the velocity components
159 // w.r.t. global coordinates (r and z)
160 DenseMatrix<double> interpolated_dudx(6, 2, 0.0);
161
162 // Calculate pressure at the integration point
163 // -------------------------------------------
164
165 // Loop over pressure degrees of freedom (associated with a single
166 // pressure component) in the element
167 for (unsigned l = 0; l < n_pres; l++)
168 {
169 // Cache the shape function
170 const double psip_ = psip(l);
171
172 // Loop over the two pressure components
173 for (unsigned i = 0; i < 2; i++)
174 {
175 // Get the value
176 const double p_value = this->p_linearised_axi_nst(l, i);
177
178 // Add contribution
179 interpolated_p[i] += p_value * psip_;
180 }
181 } // End of loop over the pressure degrees of freedom in the element
182
183 // Calculate velocities and their derivatives at the integration point
184 // -------------------------------------------------------------------
185
186 // Loop over the element's nodes
187 for (unsigned l = 0; l < n_node; l++)
188 {
189 // Cache the shape function
190 const double psif_ = psif(l);
191
192 // Loop over the two coordinate directions
193 for (unsigned i = 0; i < 2; i++)
194 {
195 interpolated_x[i] += this->nodal_position(l, i) * psif_;
196 }
197
198 // Loop over the six velocity components
199 for (unsigned i = 0; i < 6; i++)
200 {
201 // Get the value
202 const double u_value = this->nodal_value(l, u_nodal_index[i]);
203
204 // Add contribution
205 interpolated_u[i] += u_value * psif_;
206
207 // Add contribution to dudt
208 dudt[i] += du_dt_linearised_axi_nst(l, i) * psif_;
209
210 // Loop over two coordinate directions (for derivatives)
211 for (unsigned j = 0; j < 2; j++)
212 {
213 interpolated_dudx(i, j) += u_value * dpsifdx(l, j);
214 }
215 }
216 } // End of loop over the element's nodes
217
218 // Get the mesh velocity if ALE is enabled
219 if (!ALE_is_disabled)
220 {
221 // Loop over the element's nodes
222 for (unsigned l = 0; l < n_node; l++)
223 {
224 // Loop over the two coordinate directions
225 for (unsigned i = 0; i < 2; i++)
226 {
227 mesh_velocity[i] += this->dnodal_position_dt(l, i) * psif(l);
228 }
229 }
230 }
231
232 // Get velocities and their derivatives from base flow problem
233 // -----------------------------------------------------------
234
235 // Allocate storage for the velocity components of the base state
236 // solution (initialise to zero)
237 Vector<double> base_flow_u(3, 0.0);
238
239 // Get the user-defined base state solution velocity components
240 get_base_flow_u(time, ipt, interpolated_x, base_flow_u);
241
242 // Allocate storage for the derivatives of the base state solution's
243 // velocity components w.r.t. global coordinate (r and z)
244 // N.B. the derivatives of the base flow components w.r.t. the
245 // azimuthal coordinate direction (theta) are always zero since the
246 // base flow is axisymmetric
247 DenseMatrix<double> base_flow_dudx(3, 2, 0.0);
248
249 // Get the derivatives of the user-defined base state solution
250 // velocity components w.r.t. global coordinates
251 get_base_flow_dudx(time, ipt, interpolated_x, base_flow_dudx);
252
253 // Cache base flow velocities and their derivatives
254 const double interpolated_ur = base_flow_u[0];
255 const double interpolated_uz = base_flow_u[1];
256 const double interpolated_utheta = base_flow_u[2];
257 const double interpolated_durdr = base_flow_dudx(0, 0);
258 const double interpolated_durdz = base_flow_dudx(0, 1);
259 const double interpolated_duzdr = base_flow_dudx(1, 0);
260 const double interpolated_duzdz = base_flow_dudx(1, 1);
261 const double interpolated_duthetadr = base_flow_dudx(2, 0);
262 const double interpolated_duthetadz = base_flow_dudx(2, 1);
263
264 // Cache r-component of position
265 const double r = interpolated_x[0];
266
267 // Cache unknowns
268 const double interpolated_UC = interpolated_u[0];
269 const double interpolated_US = interpolated_u[1];
270 const double interpolated_WC = interpolated_u[2];
271 const double interpolated_WS = interpolated_u[3];
272 const double interpolated_VC = interpolated_u[4];
273 const double interpolated_VS = interpolated_u[5];
274 const double interpolated_PC = interpolated_p[0];
275 const double interpolated_PS = interpolated_p[1];
276
277 // Cache derivatives of the unknowns
278 const double interpolated_dUCdr = interpolated_dudx(0, 0);
279 const double interpolated_dUCdz = interpolated_dudx(0, 1);
280 const double interpolated_dUSdr = interpolated_dudx(1, 0);
281 const double interpolated_dUSdz = interpolated_dudx(1, 1);
282 const double interpolated_dWCdr = interpolated_dudx(2, 0);
283 const double interpolated_dWCdz = interpolated_dudx(2, 1);
284 const double interpolated_dWSdr = interpolated_dudx(3, 0);
285 const double interpolated_dWSdz = interpolated_dudx(3, 1);
286 const double interpolated_dVCdr = interpolated_dudx(4, 0);
287 const double interpolated_dVCdz = interpolated_dudx(4, 1);
288 const double interpolated_dVSdr = interpolated_dudx(5, 0);
289 const double interpolated_dVSdz = interpolated_dudx(5, 1);
290
291 // ==================
292 // MOMENTUM EQUATIONS
293 // ==================
294
295 // Number of master nodes
296 unsigned n_master = 1;
297
298 // Storage for the weight of the shape function
299 double hang_weight = 1.0;
300
301 // Loop over the nodes for the test functions/equations
302 for (unsigned l = 0; l < n_node; l++)
303 {
304 // Cache test functions and their derivatives
305 const double testf_ = testf(l);
306 const double dtestfdr = dtestfdx(l, 0);
307 const double dtestfdz = dtestfdx(l, 1);
308
309 // Local boolean that indicates the hanging status of the node
310 bool is_node_hanging = node_pt(l)->is_hanging();
311
312 // If the node is hanging
313 if (is_node_hanging)
314 {
315 // Get the hanging pointer
316 hang_info_pt = node_pt(l)->hanging_pt();
317
318 // Read out number of master nodes from hanging data
319 n_master = hang_info_pt->nmaster();
320 }
321 // Otherwise the node is its own master
322 else
323 {
324 n_master = 1;
325 }
326
327 // Loop over the master nodes
328 for (unsigned m = 0; m < n_master; m++)
329 {
330 // Loop over velocity components for equations
331 for (unsigned i = 0; i < 6; i++)
332 {
333 // Get the equation number
334 // -----------------------
335
336 // If the node is hanging
337 if (is_node_hanging)
338 {
339 // Get the equation number from the master node
340 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
341 u_nodal_index[i]);
342
343 // Get the hang weight from the master node
344 hang_weight = hang_info_pt->master_weight(m);
345 }
346 // If the node is not hanging
347 else
348 {
349 // Local equation number
350 local_eqn = this->nodal_local_eqn(l, u_nodal_index[i]);
351
352 // Node contributes with full weight
353 hang_weight = 1.0;
354 }
355
356 // If it's not a boundary condition...
357 if (local_eqn >= 0)
358 {
359 // initialise for residual calculation
360 double sum = 0.0;
361
362 switch (i)
363 {
364 // ---------------------------------------------
365 // FIRST (RADIAL) MOMENTUM EQUATION: COSINE PART
366 // ---------------------------------------------
367
368 case 0:
369
370 // Pressure gradient term
371 sum +=
372 interpolated_PC * (testf_ + r * dtestfdr) * W * hang_weight;
373
374 // Stress tensor terms
375 sum -= visc_ratio * r * (1.0 + Gamma[0]) *
376 interpolated_dUCdr * dtestfdr * W * hang_weight;
377
378 sum -= visc_ratio * r *
379 (interpolated_dUCdz + Gamma[0] * interpolated_dWCdr) *
380 dtestfdz * W * hang_weight;
381
382 sum += visc_ratio *
383 ((k * Gamma[0] * interpolated_dVSdr) -
384 (k * (2.0 + Gamma[0]) * interpolated_VS / r) -
385 ((1.0 + Gamma[0] + (k * k)) * interpolated_UC / r)) *
386 testf_ * W * hang_weight;
387
388 // Inertial terms (du/dt)
389 sum -= scaled_re_st * r * dudt[0] * testf_ * W * hang_weight;
390
391 // Inertial terms (convective)
392 sum -= scaled_re *
393 (r * interpolated_ur * interpolated_dUCdr +
394 r * interpolated_durdr * interpolated_UC +
395 k * interpolated_utheta * interpolated_US -
396 2 * interpolated_utheta * interpolated_VC +
397 r * interpolated_uz * interpolated_dUCdz +
398 r * interpolated_durdz * interpolated_WC) *
399 testf_ * W * hang_weight;
400
401 // Mesh velocity terms
402 if (!ALE_is_disabled)
403 {
404 for (unsigned j = 0; j < 2; j++)
405 {
406 sum += scaled_re_st * r * mesh_velocity[j] *
407 interpolated_dudx(0, j) * testf_ * W * hang_weight;
408 }
409 }
410
411 break;
412
413 // --------------------------------------------
414 // SECOND (RADIAL) MOMENTUM EQUATION: SINE PART
415 // --------------------------------------------
416
417 case 1:
418
419 // Pressure gradient term
420 sum +=
421 interpolated_PS * (testf_ + r * dtestfdr) * W * hang_weight;
422
423 // Stress tensor terms
424 sum -= visc_ratio * r * (1.0 + Gamma[0]) *
425 interpolated_dUSdr * dtestfdr * W * hang_weight;
426
427 sum -= visc_ratio * r *
428 (interpolated_dUSdz + Gamma[0] * interpolated_dWSdr) *
429 dtestfdz * W * hang_weight;
430
431 sum -= visc_ratio *
432 ((k * Gamma[0] * interpolated_dVCdr) -
433 (k * (2.0 + Gamma[0]) * interpolated_VC / r) +
434 ((1.0 + Gamma[0] + (k * k)) * interpolated_US / r)) *
435 testf_ * W * hang_weight;
436
437 // Inertial terms (du/dt)
438 sum -= scaled_re_st * r * dudt[1] * testf_ * W * hang_weight;
439
440 // Inertial terms (convective)
441 sum -= scaled_re *
442 (r * interpolated_ur * interpolated_dUSdr +
443 r * interpolated_durdr * interpolated_US -
444 k * interpolated_utheta * interpolated_UC -
445 2 * interpolated_utheta * interpolated_VS +
446 r * interpolated_uz * interpolated_dUSdz +
447 r * interpolated_durdz * interpolated_WS) *
448 testf_ * W * hang_weight;
449
450 // Mesh velocity terms
451 if (!ALE_is_disabled)
452 {
453 for (unsigned j = 0; j < 2; j++)
454 {
455 sum += scaled_re_st * r * mesh_velocity[j] *
456 interpolated_dudx(1, j) * testf_ * W * hang_weight;
457 }
458 }
459
460 break;
461
462 // --------------------------------------------
463 // THIRD (AXIAL) MOMENTUM EQUATION: COSINE PART
464 // --------------------------------------------
465
466 case 2:
467
468 // Pressure gradient term
469 sum += r * interpolated_PC * dtestfdz * W * hang_weight;
470
471 // Stress tensor terms
472 sum -= visc_ratio * r *
473 (interpolated_dWCdr + Gamma[1] * interpolated_dUCdz) *
474 dtestfdr * W * hang_weight;
475
476 sum -= visc_ratio * r * (1.0 + Gamma[1]) *
477 interpolated_dWCdz * dtestfdz * W * hang_weight;
478
479 sum +=
480 visc_ratio * k *
481 (Gamma[1] * interpolated_dVSdz - k * interpolated_WC / r) *
482 testf_ * W * hang_weight;
483
484 // Inertial terms (du/dt)
485 sum -= scaled_re_st * r * dudt[2] * testf_ * W * hang_weight;
486
487 // Inertial terms (convective)
488 sum -= scaled_re *
489 (r * interpolated_ur * interpolated_dWCdr +
490 r * interpolated_duzdr * interpolated_UC +
491 k * interpolated_utheta * interpolated_WS +
492 r * interpolated_uz * interpolated_dWCdz +
493 r * interpolated_duzdz * interpolated_WC) *
494 testf_ * W * hang_weight;
495
496 // Mesh velocity terms
497 if (!ALE_is_disabled)
498 {
499 for (unsigned j = 0; j < 2; j++)
500 {
501 sum += scaled_re_st * r * mesh_velocity[j] *
502 interpolated_dudx(2, j) * testf_ * W * hang_weight;
503 }
504 }
505
506 break;
507
508 // -------------------------------------------
509 // FOURTH (AXIAL) MOMENTUM EQUATION: SINE PART
510 // -------------------------------------------
511
512 case 3:
513
514 // Pressure gradient term
515 sum += r * interpolated_PS * dtestfdz * W * hang_weight;
516
517 // Stress tensor terms
518 sum -= visc_ratio * r *
519 (interpolated_dWSdr + Gamma[1] * interpolated_dUSdz) *
520 dtestfdr * W * hang_weight;
521
522 sum -= visc_ratio * r * (1.0 + Gamma[1]) *
523 interpolated_dWSdz * dtestfdz * W * hang_weight;
524
525 sum -=
526 visc_ratio * k *
527 (Gamma[1] * interpolated_dVCdz + k * interpolated_WS / r) *
528 testf_ * W * hang_weight;
529
530 // Inertial terms (du/dt)
531 sum -= scaled_re_st * r * dudt[3] * testf_ * W * hang_weight;
532
533 // Inertial terms (convective)
534 sum -= scaled_re *
535 (r * interpolated_ur * interpolated_dWSdr +
536 r * interpolated_duzdr * interpolated_US -
537 k * interpolated_utheta * interpolated_WC +
538 r * interpolated_uz * interpolated_dWSdz +
539 r * interpolated_duzdz * interpolated_WS) *
540 testf_ * W * hang_weight;
541
542 // Mesh velocity terms
543 if (!ALE_is_disabled)
544 {
545 for (unsigned j = 0; j < 2; j++)
546 {
547 sum += scaled_re_st * r * mesh_velocity[j] *
548 interpolated_dudx(3, j) * testf_ * W * hang_weight;
549 }
550 }
551
552 break;
553
554 // ------------------------------------------------
555 // FIFTH (AZIMUTHAL) MOMENTUM EQUATION: COSINE PART
556 // ------------------------------------------------
557
558 case 4:
559
560 // Pressure gradient term
561 sum -= k * interpolated_PS * testf_ * W * hang_weight;
562
563 // Stress tensor terms
564 sum +=
565 visc_ratio *
566 (-r * interpolated_dVCdr - k * Gamma[0] * interpolated_US +
567 Gamma[0] * interpolated_VC) *
568 dtestfdr * W * hang_weight;
569
570 sum -=
571 visc_ratio *
572 (k * Gamma[0] * interpolated_WS + r * interpolated_dVCdz) *
573 dtestfdz * W * hang_weight;
574
575 sum += visc_ratio *
576 (Gamma[0] * interpolated_dVCdr +
577 k * (2.0 + Gamma[0]) * interpolated_US / r -
578 (1.0 + (k * k) + (k * k * Gamma[0])) *
579 interpolated_VC / r) *
580 testf_ * W * hang_weight;
581
582 // Inertial terms (du/dt)
583 sum -= scaled_re_st * r * dudt[4] * testf_ * W * hang_weight;
584
585 // Inertial terms (convective)
586 sum -= scaled_re *
587 (r * interpolated_ur * interpolated_dVCdr +
588 r * interpolated_duthetadr * interpolated_UC +
589 k * interpolated_utheta * interpolated_VS +
590 interpolated_utheta * interpolated_UC +
591 interpolated_ur * interpolated_VC +
592 r * interpolated_uz * interpolated_dVCdz +
593 r * interpolated_duthetadz * interpolated_WC) *
594 testf_ * W * hang_weight;
595
596 // Mesh velocity terms
597 if (!ALE_is_disabled)
598 {
599 for (unsigned j = 0; j < 2; j++)
600 {
601 sum += scaled_re_st * r * mesh_velocity[j] *
602 interpolated_dudx(4, j) * testf_ * W * hang_weight;
603 }
604 }
605
606 break;
607
608 // ----------------------------------------------
609 // SIXTH (AZIMUTHAL) MOMENTUM EQUATION: SINE PART
610 // ----------------------------------------------
611
612 case 5:
613
614 // Pressure gradient term
615 sum += k * interpolated_PC * testf_ * W * hang_weight;
616
617 // Stress tensor terms
618 sum +=
619 visc_ratio *
620 (-r * interpolated_dVSdr + k * Gamma[0] * interpolated_UC +
621 Gamma[0] * interpolated_VS) *
622 dtestfdr * W * hang_weight;
623
624 sum +=
625 visc_ratio *
626 (k * Gamma[0] * interpolated_WC - r * interpolated_dVSdz) *
627 dtestfdz * W * hang_weight;
628
629 sum += visc_ratio *
630 (Gamma[0] * interpolated_dVSdr -
631 k * (2.0 + Gamma[0]) * interpolated_UC / r -
632 (1.0 + (k * k) + (k * k * Gamma[0])) *
633 interpolated_VS / r) *
634 testf_ * W * hang_weight;
635
636 // Inertial terms (du/dt)
637 sum -= scaled_re_st * r * dudt[5] * testf_ * W * hang_weight;
638
639 // Inertial terms (convective)
640 sum -= scaled_re *
641 (r * interpolated_ur * interpolated_dVSdr +
642 r * interpolated_duthetadr * interpolated_US -
643 k * interpolated_utheta * interpolated_VC +
644 interpolated_utheta * interpolated_US +
645 interpolated_ur * interpolated_VS +
646 r * interpolated_uz * interpolated_dVSdz +
647 r * interpolated_duthetadz * interpolated_WS) *
648 testf_ * W * hang_weight;
649
650 // Mesh velocity terms
651 if (!ALE_is_disabled)
652 {
653 for (unsigned j = 0; j < 2; j++)
654 {
655 sum += scaled_re_st * r * mesh_velocity[j] *
656 interpolated_dudx(5, j) * testf_ * W * hang_weight;
657 }
658 }
659
660 break;
661
662 } // End of switch statement for momentum equations
663
664 // Add contribution to elemental residual vector
665 residuals[local_eqn] += sum;
666
667 // ======================
668 // Calculate the Jacobian
669 // ======================
670 if (flag)
671 {
672 // Number of master nodes
673 unsigned n_master2 = 1;
674
675 // Storage for the weight of the shape function
676 double hang_weight2 = 1.0;
677
678 // Loop over the velocity nodes for columns
679 for (unsigned l2 = 0; l2 < n_node; l2++)
680 {
681 // Cache velocity shape functions and their derivatives
682 const double psif_ = psif[l2];
683 const double dpsifdr = dpsifdx(l2, 0);
684 const double dpsifdz = dpsifdx(l2, 1);
685
686 // Local boolean that indicates the hanging status of the node
687 bool is_node2_hanging = node_pt(l2)->is_hanging();
688
689 // If the node is hanging
690 if (is_node2_hanging)
691 {
692 // Get the hanging pointer
693 hang_info2_pt = node_pt(l2)->hanging_pt();
694
695 // Read out number of master nodes from hanging data
696 n_master2 = hang_info2_pt->nmaster();
697 }
698 // Otherwise the node is its own master
699 else
700 {
701 n_master2 = 1;
702 }
703
704 // Loop over the master nodes
705 for (unsigned m2 = 0; m2 < n_master2; m2++)
706 {
707 // Loop over the velocity components
708 for (unsigned i2 = 0; i2 < 6; i2++)
709 {
710 // Get the number of the unknown
711 // -----------------------------
712
713 // If the node is hanging
714 if (is_node2_hanging)
715 {
716 // Get the equation number from the master node
717 local_unknown = this->local_hang_eqn(
718 hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
719
720 // Get the hang weight from the master node
721 hang_weight2 = hang_info2_pt->master_weight(m2);
722 }
723 // If the node is not hanging
724 else
725 {
726 // Local equation number
727 local_unknown =
728 this->nodal_local_eqn(l2, u_nodal_index[i2]);
729
730 // Node contributes with full weight
731 hang_weight2 = 1.0;
732 }
733
734 // If the unknown is not pinned
735 if (local_unknown >= 0)
736 {
737 // Different results for i and i2
738 switch (i)
739 {
740 // ---------------------------------------------
741 // FIRST (RADIAL) MOMENTUM EQUATION: COSINE PART
742 // ---------------------------------------------
743
744 case 0:
745
746 switch (i2)
747 {
748 // Radial velocity component (cosine part) U_k^C
749
750 case 0:
751
752 // Add the mass matrix entries
753 if (flag == 2)
754 {
755 mass_matrix(local_eqn, local_unknown) +=
756 scaled_re_st * r * psif_ * testf_ * W *
757 hang_weight * hang_weight2 * hang_weight *
758 hang_weight2;
759 }
760
761 // Add contributions to the Jacobian matrix
762
763 // Stress tensor terms
764 jacobian(local_eqn, local_unknown) -=
765 visc_ratio * r * (1.0 + Gamma[0]) * dpsifdr *
766 dtestfdr * W * hang_weight * hang_weight2;
767
768 jacobian(local_eqn, local_unknown) -=
769 visc_ratio * r * dpsifdz * dtestfdz * W *
770 hang_weight * hang_weight2;
771
772 jacobian(local_eqn, local_unknown) -=
773 visc_ratio * (1.0 + Gamma[0] + (k * k)) *
774 psif_ * testf_ * W * hang_weight *
775 hang_weight2 / r;
776
777 // Inertial terms (du/dt)
778 jacobian(local_eqn, local_unknown) -=
779 scaled_re_st * r *
780 node_pt(l2)->time_stepper_pt()->weight(1, 0) *
781 psif_ * testf_ * W * hang_weight *
782 hang_weight2;
783
784 // Inertial terms (convective)
785 jacobian(local_eqn, local_unknown) -=
786 scaled_re * r *
787 (psif_ * interpolated_durdr +
788 interpolated_ur * dpsifdr +
789 interpolated_uz * dpsifdz) *
790 testf_ * W * hang_weight * hang_weight2;
791
792 // Mesh velocity terms
793 if (!ALE_is_disabled)
794 {
795 for (unsigned j = 0; j < 2; j++)
796 {
797 jacobian(local_eqn, local_unknown) +=
798 scaled_re_st * r * mesh_velocity[j] *
799 dpsifdx(l2, j) * testf_ * W *
800 hang_weight * hang_weight2;
801 }
802 }
803
804 break;
805
806 // Radial velocity component (sine part) U_k^S
807
808 case 1:
809
810 // Convective term
811 jacobian(local_eqn, local_unknown) -=
812 scaled_re * k * interpolated_utheta * psif_ *
813 testf_ * W * hang_weight * hang_weight2;
814
815 break;
816
817 // Axial velocity component (cosine part) W_k^C
818
819 case 2:
820
821 // Stress tensor term
822 jacobian(local_eqn, local_unknown) -=
823 visc_ratio * Gamma[0] * r * dpsifdr *
824 dtestfdz * W * hang_weight * hang_weight2;
825
826 // Convective term
827 jacobian(local_eqn, local_unknown) -=
828 scaled_re * r * interpolated_durdz * psif_ *
829 testf_ * W * hang_weight * hang_weight2;
830
831 break;
832
833 // Axial velocity component (sine part) W_k^S
834 // has no contribution
835
836 case 3:
837 break;
838
839 // Azimuthal velocity component (cosine part)
840 // V_k^C
841
842 case 4:
843
844 // Convective term
845 jacobian(local_eqn, local_unknown) +=
846 scaled_re * 2 * interpolated_utheta * psif_ *
847 testf_ * W * hang_weight * hang_weight2;
848
849 break;
850
851 // Azimuthal velocity component (sine part)
852 // V_k^S
853
854 case 5:
855
856 // Stress tensor term
857 jacobian(local_eqn, local_unknown) +=
858 visc_ratio *
859 ((Gamma[0] * k * dpsifdr) -
860 (k * (2.0 + Gamma[0]) * psif_ / r)) *
861 testf_ * W * hang_weight * hang_weight2;
862
863 break;
864
865 } // End of first (radial) momentum equation
866 break;
867
868 // --------------------------------------------
869 // SECOND (RADIAL) MOMENTUM EQUATION: SINE PART
870 // --------------------------------------------
871
872 case 1:
873
874 switch (i2)
875 {
876 // Radial velocity component (cosine part) U_k^C
877
878 case 0:
879
880 // Convective term
881 jacobian(local_eqn, local_unknown) +=
882 scaled_re * k * interpolated_utheta * psif_ *
883 testf_ * W * hang_weight * hang_weight2;
884
885 break;
886
887 // Radial velocity component (sine part) U_k^S
888
889 case 1:
890
891 // Add the mass matrix entries
892 if (flag == 2)
893 {
894 mass_matrix(local_eqn, local_unknown) +=
895 scaled_re_st * r * psif_ * testf_ * W *
896 hang_weight * hang_weight2;
897 }
898
899 // Add contributions to the Jacobian matrix
900
901 // Stress tensor terms
902 jacobian(local_eqn, local_unknown) -=
903 visc_ratio * r * (1.0 + Gamma[0]) * dpsifdr *
904 dtestfdr * W * hang_weight * hang_weight2;
905
906 jacobian(local_eqn, local_unknown) -=
907 visc_ratio * r * dpsifdz * dtestfdz * W *
908 hang_weight * hang_weight2;
909
910 jacobian(local_eqn, local_unknown) -=
911 visc_ratio * (1.0 + Gamma[0] + (k * k)) *
912 psif_ * testf_ * W * hang_weight *
913 hang_weight2 / r;
914
915 // Inertial terms (du/dt)
916 jacobian(local_eqn, local_unknown) -=
917 scaled_re_st * r *
918 node_pt(l2)->time_stepper_pt()->weight(1, 0) *
919 psif_ * testf_ * W * hang_weight *
920 hang_weight2;
921
922 // Inertial terms (convective)
923 jacobian(local_eqn, local_unknown) -=
924 scaled_re * r *
925 (psif_ * interpolated_durdr +
926 interpolated_ur * dpsifdr +
927 interpolated_uz * dpsifdz) *
928 testf_ * W * hang_weight * hang_weight2;
929
930 // Mesh velocity terms
931 if (!ALE_is_disabled)
932 {
933 for (unsigned j = 0; j < 2; j++)
934 {
935 jacobian(local_eqn, local_unknown) +=
936 scaled_re_st * r * mesh_velocity[j] *
937 dpsifdx(l2, j) * testf_ * W *
938 hang_weight * hang_weight2;
939 }
940 }
941
942 break;
943
944 // Axial velocity component (cosine part) W_k^C
945 // has no contribution
946
947 case 2:
948 break;
949
950 // Axial velocity component (sine part) W_k^S
951
952 case 3:
953
954 // Stress tensor term
955 jacobian(local_eqn, local_unknown) -=
956 visc_ratio * Gamma[0] * r * dpsifdr *
957 dtestfdz * W * hang_weight * hang_weight2;
958
959 // Convective term
960 jacobian(local_eqn, local_unknown) -=
961 scaled_re * r * interpolated_durdz * psif_ *
962 testf_ * W * hang_weight * hang_weight2;
963
964 break;
965
966 // Azimuthal velocity component (cosine part)
967 // V_k^C
968
969 case 4:
970
971 // Stress tensor terms
972 jacobian(local_eqn, local_unknown) -=
973 visc_ratio *
974 ((Gamma[0] * k * dpsifdr) -
975 (k * (2.0 + Gamma[0]) * psif_ / r)) *
976 testf_ * W * hang_weight * hang_weight2;
977
978 break;
979
980 // Azimuthal velocity component (sine part)
981 // V_k^S
982
983 case 5:
984
985 // Convective term
986 jacobian(local_eqn, local_unknown) +=
987 scaled_re * 2 * interpolated_utheta * psif_ *
988 testf_ * W * hang_weight * hang_weight2;
989
990 break;
991
992 } // End of second (radial) momentum equation
993 break;
994
995 // --------------------------------------------
996 // THIRD (AXIAL) MOMENTUM EQUATION: COSINE PART
997 // --------------------------------------------
998
999 case 2:
1000
1001 switch (i2)
1002 {
1003 // Radial velocity component (cosine part) U_k^C
1004
1005 case 0:
1006
1007 // Stress tensor term
1008 jacobian(local_eqn, local_unknown) -=
1009 visc_ratio * r * Gamma[1] * dpsifdz *
1010 dtestfdr * W * hang_weight * hang_weight2;
1011
1012 // Convective term
1013 jacobian(local_eqn, local_unknown) -=
1014 scaled_re * r * psif_ * interpolated_duzdr *
1015 testf_ * W * hang_weight * hang_weight2;
1016
1017 break;
1018
1019 // Radial velocity component (sine part) U_k^S
1020 // has no contribution
1021
1022 case 1:
1023 break;
1024
1025 // Axial velocity component (cosine part) W_k^C
1026
1027 case 2:
1028
1029 // Add the mass matrix entries
1030 if (flag == 2)
1031 {
1032 mass_matrix(local_eqn, local_unknown) +=
1033 scaled_re_st * r * psif_ * testf_ * W *
1034 hang_weight * hang_weight2;
1035 }
1036
1037 // Add contributions to the Jacobian matrix
1038
1039 // Stress tensor terms
1040 jacobian(local_eqn, local_unknown) -=
1041 visc_ratio * r * dpsifdr * dtestfdr * W *
1042 hang_weight * hang_weight2;
1043
1044 jacobian(local_eqn, local_unknown) -=
1045 visc_ratio * r * (1.0 + Gamma[1]) * dpsifdz *
1046 dtestfdz * W * hang_weight * hang_weight2;
1047
1048 jacobian(local_eqn, local_unknown) -=
1049 visc_ratio * k * k * psif_ * testf_ * W *
1050 hang_weight * hang_weight2 / r;
1051
1052 // Inertial terms (du/dt)
1053 jacobian(local_eqn, local_unknown) -=
1054 scaled_re_st * r *
1055 node_pt(l2)->time_stepper_pt()->weight(1, 0) *
1056 psif_ * testf_ * W * hang_weight *
1057 hang_weight2;
1058
1059 // Inertial terms (convective)
1060 jacobian(local_eqn, local_unknown) -=
1061 scaled_re * r *
1062 (interpolated_ur * dpsifdr +
1063 psif_ * interpolated_duzdz +
1064 interpolated_uz * dpsifdz) *
1065 testf_ * W * hang_weight * hang_weight2;
1066
1067
1068 // Mesh velocity terms
1069 if (!ALE_is_disabled)
1070 {
1071 for (unsigned j = 0; j < 2; j++)
1072 {
1073 jacobian(local_eqn, local_unknown) +=
1074 scaled_re_st * r * mesh_velocity[j] *
1075 dpsifdx(l2, j) * testf_ * W *
1076 hang_weight * hang_weight2;
1077 }
1078 }
1079
1080 break;
1081
1082 // Axial velocity component (sine part) W_k^S
1083
1084 case 3:
1085
1086 // Convective term
1087 jacobian(local_eqn, local_unknown) -=
1088 scaled_re * k * interpolated_utheta * psif_ *
1089 testf_ * W * hang_weight * hang_weight2;
1090
1091 break;
1092
1093 // Azimuthal velocity component (cosine part)
1094 // V_k^C has no contribution
1095
1096 case 4:
1097 break;
1098
1099 // Azimuthal velocity component (sine part)
1100 // V_k^S
1101
1102 case 5:
1103
1104 // Stress tensor term
1105 jacobian(local_eqn, local_unknown) +=
1106 visc_ratio * Gamma[1] * k * dpsifdz * testf_ *
1107 W * hang_weight * hang_weight2;
1108
1109 break;
1110
1111 } // End of third (axial) momentum equation
1112 break;
1113
1114 // -------------------------------------------
1115 // FOURTH (AXIAL) MOMENTUM EQUATION: SINE PART
1116 // -------------------------------------------
1117
1118 case 3:
1119
1120 switch (i2)
1121 {
1122 // Radial velocity component (cosine part) U_k^C
1123 // has no contribution
1124
1125 case 0:
1126 break;
1127
1128 // Radial velocity component (sine part) U_k^S
1129
1130 case 1:
1131
1132 // Stress tensor term
1133 jacobian(local_eqn, local_unknown) -=
1134 visc_ratio * r * Gamma[1] * dpsifdz *
1135 dtestfdr * W * hang_weight * hang_weight2;
1136
1137 // Convective term
1138 jacobian(local_eqn, local_unknown) -=
1139 scaled_re * r * psif_ * interpolated_duzdr *
1140 testf_ * W * hang_weight * hang_weight2;
1141
1142 break;
1143
1144 // Axial velocity component (cosine part) W_k^S
1145
1146 case 2:
1147
1148 // Convective term
1149 jacobian(local_eqn, local_unknown) +=
1150 scaled_re * k * interpolated_utheta * psif_ *
1151 testf_ * W * hang_weight * hang_weight2;
1152
1153 break;
1154
1155 // Axial velocity component (sine part) W_k^S
1156
1157 case 3:
1158
1159 // Add the mass matrix entries
1160 if (flag == 2)
1161 {
1162 mass_matrix(local_eqn, local_unknown) +=
1163 scaled_re_st * r * psif_ * testf_ * W *
1164 hang_weight * hang_weight2;
1165 }
1166
1167 // Add contributions to the Jacobian matrix
1168
1169 // Stress tensor terms
1170 jacobian(local_eqn, local_unknown) -=
1171 visc_ratio * r * dpsifdr * dtestfdr * W *
1172 hang_weight * hang_weight2;
1173
1174 jacobian(local_eqn, local_unknown) -=
1175 visc_ratio * r * (1.0 + Gamma[1]) * dpsifdz *
1176 dtestfdz * W * hang_weight * hang_weight2;
1177
1178 jacobian(local_eqn, local_unknown) -=
1179 visc_ratio * k * k * psif_ * testf_ * W *
1180 hang_weight * hang_weight2 / r;
1181
1182 // Inertial terms (du/dt)
1183 jacobian(local_eqn, local_unknown) -=
1184 scaled_re_st * r *
1185 node_pt(l2)->time_stepper_pt()->weight(1, 0) *
1186 psif_ * testf_ * W * hang_weight *
1187 hang_weight2;
1188
1189 // Inertial terms (convective)
1190 jacobian(local_eqn, local_unknown) -=
1191 scaled_re * r *
1192 (interpolated_ur * dpsifdr +
1193 psif_ * interpolated_duzdz +
1194 interpolated_uz * dpsifdz) *
1195 testf_ * W * hang_weight * hang_weight2;
1196
1197
1198 // Mesh velocity terms
1199 if (!ALE_is_disabled)
1200 {
1201 for (unsigned j = 0; j < 2; j++)
1202 {
1203 jacobian(local_eqn, local_unknown) +=
1204 scaled_re_st * r * mesh_velocity[j] *
1205 dpsifdx(l2, j) * testf_ * W *
1206 hang_weight * hang_weight2;
1207 }
1208 }
1209
1210 break;
1211
1212 // Azimuthal velocity component (cosine part)
1213 // V_k^C
1214
1215 case 4:
1216
1217 // Stress tensor term
1218 jacobian(local_eqn, local_unknown) -=
1219 visc_ratio * Gamma[1] * k * dpsifdz * testf_ *
1220 W * hang_weight * hang_weight2;
1221
1222 break;
1223
1224 // Azimuthal velocity component (sine part)
1225 // V_k^S has no contribution
1226
1227 case 5:
1228 break;
1229
1230 } // End of fourth (axial) momentum equation
1231 break;
1232
1233 // ------------------------------------------------
1234 // FIFTH (AZIMUTHAL) MOMENTUM EQUATION: COSINE PART
1235 // ------------------------------------------------
1236
1237 case 4:
1238
1239 switch (i2)
1240 {
1241 // Radial velocity component (cosine part) U_k^C
1242
1243 case 0:
1244
1245 // Convective terms
1246 jacobian(local_eqn, local_unknown) -=
1247 scaled_re *
1248 (r * interpolated_duthetadr +
1249 interpolated_utheta) *
1250 psif_ * testf_ * W * hang_weight *
1251 hang_weight2;
1252
1253 break;
1254
1255 // Radial velocity component (sine part) U_k^S
1256
1257 case 1:
1258
1259 // Stress tensor terms
1260 jacobian(local_eqn, local_unknown) +=
1261 visc_ratio * k * psif_ *
1262 (((2.0 + Gamma[0]) * testf_ / r) -
1263 (Gamma[0] * dtestfdr)) *
1264 W * hang_weight * hang_weight2;
1265
1266 break;
1267
1268 // Axial velocity component (cosine part) W_k^C
1269
1270 case 2:
1271
1272 // Convective term
1273 jacobian(local_eqn, local_unknown) -=
1274 scaled_re * r * psif_ *
1275 interpolated_duthetadz * testf_ * W *
1276 hang_weight * hang_weight2;
1277
1278 break;
1279
1280 // Axial velocity component (sine part) W_k^S
1281
1282 case 3:
1283
1284 // Stress tensor term
1285 jacobian(local_eqn, local_unknown) -=
1286 visc_ratio * k * Gamma[0] * psif_ * dtestfdz *
1287 W * hang_weight * hang_weight2;
1288
1289 break;
1290
1291 // Azimuthal velocity component (cosine part)
1292 // V_k^C
1293
1294 case 4:
1295
1296 // Add the mass matrix entries
1297 if (flag == 2)
1298 {
1299 mass_matrix(local_eqn, local_unknown) +=
1300 scaled_re_st * r * psif_ * testf_ * W *
1301 hang_weight * hang_weight2;
1302 }
1303
1304 // Add contributions to the Jacobian matrix
1305
1306 // Stress tensor terms
1307 jacobian(local_eqn, local_unknown) +=
1308 visc_ratio *
1309 (Gamma[0] * psif_ - r * dpsifdr) * dtestfdr *
1310 W * hang_weight * hang_weight2;
1311
1312 jacobian(local_eqn, local_unknown) -=
1313 visc_ratio * r * dpsifdz * dtestfdz * W *
1314 hang_weight * hang_weight2;
1315
1316 jacobian(local_eqn, local_unknown) +=
1317 visc_ratio *
1318 (Gamma[0] * dpsifdr -
1319 (1.0 + (k * k) + (k * k * Gamma[0])) *
1320 psif_ / r) *
1321 testf_ * W * hang_weight * hang_weight2;
1322
1323 // Inertial terms (du/dt)
1324 jacobian(local_eqn, local_unknown) -=
1325 scaled_re_st * r *
1326 node_pt(l2)->time_stepper_pt()->weight(1, 0) *
1327 psif_ * testf_ * W * hang_weight *
1328 hang_weight2;
1329
1330 // Inertial terms (convective)
1331 jacobian(local_eqn, local_unknown) -=
1332 scaled_re *
1333 (r * interpolated_ur * dpsifdr +
1334 interpolated_ur * psif_ +
1335 r * interpolated_uz * dpsifdz) *
1336 testf_ * W * hang_weight * hang_weight2;
1337
1338 // Mesh velocity terms
1339 if (!ALE_is_disabled)
1340 {
1341 for (unsigned j = 0; j < 2; j++)
1342 {
1343 jacobian(local_eqn, local_unknown) +=
1344 scaled_re_st * r * mesh_velocity[j] *
1345 dpsifdx(l2, j) * testf_ * W *
1346 hang_weight * hang_weight2;
1347 }
1348 }
1349
1350 break;
1351
1352 // Azimuthal velocity component (sine part)
1353 // V_k^S
1354
1355 case 5:
1356
1357 // Convective term
1358 jacobian(local_eqn, local_unknown) -=
1359 scaled_re * k * interpolated_utheta * psif_ *
1360 testf_ * W * hang_weight * hang_weight2;
1361
1362 break;
1363
1364 } // End of fifth (azimuthal) momentum equation
1365 break;
1366
1367 // ----------------------------------------------
1368 // SIXTH (AZIMUTHAL) MOMENTUM EQUATION: SINE PART
1369 // ----------------------------------------------
1370
1371 case 5:
1372
1373 switch (i2)
1374 {
1375 // Radial velocity component (cosine part) U_k^C
1376
1377 case 0:
1378
1379 // Stress tensor terms
1380 jacobian(local_eqn, local_unknown) +=
1381 visc_ratio * k * psif_ *
1382 ((Gamma[0] * dtestfdr) -
1383 ((2.0 + Gamma[0]) * testf_ / r)) *
1384 W * hang_weight * hang_weight2;
1385
1386 break;
1387
1388 // Radial velocity component (sine part) U_k^S
1389
1390 case 1:
1391
1392 // Convective terms
1393 jacobian(local_eqn, local_unknown) -=
1394 scaled_re *
1395 (r * interpolated_duthetadr +
1396 interpolated_utheta) *
1397 psif_ * testf_ * W * hang_weight *
1398 hang_weight2;
1399
1400 break;
1401
1402 // Axial velocity component (cosine part) W_k^C
1403
1404 case 2:
1405
1406 // Stress tensor term
1407 jacobian(local_eqn, local_unknown) +=
1408 visc_ratio * k * Gamma[0] * psif_ * dtestfdz *
1409 W * hang_weight * hang_weight2;
1410
1411 break;
1412
1413 // Axial velocity component (sine part) W_k^S
1414
1415 case 3:
1416
1417 // Convective term
1418 jacobian(local_eqn, local_unknown) -=
1419 scaled_re * r * psif_ *
1420 interpolated_duthetadz * testf_ * W *
1421 hang_weight * hang_weight2;
1422
1423 break;
1424
1425 // Azimuthal velocity component (cosine part)
1426 // V_k^C
1427
1428 case 4:
1429
1430 // Convective term
1431 jacobian(local_eqn, local_unknown) +=
1432 scaled_re * k * interpolated_utheta * psif_ *
1433 testf_ * W * hang_weight * hang_weight2;
1434
1435 break;
1436
1437 // Azimuthal velocity component (sine part)
1438 // V_k^S
1439
1440 case 5:
1441
1442 // Add the mass matrix entries
1443 if (flag == 2)
1444 {
1445 mass_matrix(local_eqn, local_unknown) +=
1446 scaled_re_st * r * psif_ * testf_ * W *
1447 hang_weight * hang_weight2;
1448 }
1449
1450 // Add contributions to the Jacobian matrix
1451
1452 // Stress tensor terms
1453 jacobian(local_eqn, local_unknown) +=
1454 visc_ratio *
1455 (Gamma[0] * psif_ - r * dpsifdr) * dtestfdr *
1456 W * hang_weight * hang_weight2;
1457
1458 jacobian(local_eqn, local_unknown) -=
1459 visc_ratio * r * dpsifdz * dtestfdz * W *
1460 hang_weight * hang_weight2;
1461
1462 jacobian(local_eqn, local_unknown) +=
1463 visc_ratio *
1464 (Gamma[0] * dpsifdr -
1465 (1.0 + (k * k) + (k * k * Gamma[0])) *
1466 psif_ / r) *
1467 testf_ * W * hang_weight * hang_weight2;
1468
1469 // Inertial terms (du/dt)
1470 jacobian(local_eqn, local_unknown) -=
1471 scaled_re_st * r *
1472 node_pt(l2)->time_stepper_pt()->weight(1, 0) *
1473 psif_ * testf_ * W * hang_weight *
1474 hang_weight2;
1475
1476 // Convective terms
1477 jacobian(local_eqn, local_unknown) -=
1478 scaled_re *
1479 (r * interpolated_ur * dpsifdr +
1480 interpolated_ur * psif_ +
1481 r * interpolated_uz * dpsifdz) *
1482 testf_ * W * hang_weight * hang_weight2;
1483
1484 // Mesh velocity terms
1485 if (!ALE_is_disabled)
1486 {
1487 for (unsigned j = 0; j < 2; j++)
1488 {
1489 jacobian(local_eqn, local_unknown) +=
1490 scaled_re_st * r * mesh_velocity[j] *
1491 dpsifdx(l2, j) * testf_ * W *
1492 hang_weight * hang_weight2;
1493 }
1494 }
1495
1496 break;
1497
1498 } // End of sixth (azimuthal) momentum equation
1499 break;
1500 }
1501 } // End of if not boundary condition statement
1502 } // End of loop over velocity components
1503 } // End of loop over master (m2) nodes
1504 } // End of loop over the velocity nodes (l2)
1505
1506
1507 // Loop over the pressure shape functions
1508 for (unsigned l2 = 0; l2 < n_pres; l2++)
1509 {
1510 // If the pressure dof is hanging
1511 if (pressure_dof_is_hanging[l2])
1512 {
1513 // Pressure dof is hanging so it must be nodal-based
1514 hang_info2_pt =
1515 pressure_node_pt(l2)->hanging_pt(p_index[0]);
1516
1517 // Get the number of master nodes from the pressure node
1518 n_master2 = hang_info2_pt->nmaster();
1519 }
1520 // Otherwise the node is its own master
1521 else
1522 {
1523 n_master2 = 1;
1524 }
1525
1526 // Loop over the master nodes
1527 for (unsigned m2 = 0; m2 < n_master2; m2++)
1528 {
1529 // Loop over the two pressure components
1530 for (unsigned i2 = 0; i2 < 2; i2++)
1531 {
1532 // Get the number of the unknown
1533 // -----------------------------
1534
1535 // If the pressure dof is hanging
1536 if (pressure_dof_is_hanging[l2])
1537 {
1538 // Get the unknown from the node
1539 local_unknown = local_hang_eqn(
1540 hang_info2_pt->master_node_pt(m2), p_index[i2]);
1541
1542 // Get the weight from the hanging object
1543 hang_weight2 = hang_info2_pt->master_weight(m2);
1544 }
1545 else
1546 {
1547 local_unknown = p_local_eqn(l2, i2);
1548 hang_weight2 = 1.0;
1549 }
1550
1551 // If the unknown is not pinned
1552 if (local_unknown >= 0)
1553 {
1554 // Add in contributions to momentum equations
1555 switch (i)
1556 {
1557 // ---------------------------------------------
1558 // FIRST (RADIAL) MOMENTUM EQUATION: COSINE PART
1559 // ---------------------------------------------
1560
1561 case 0:
1562
1563 switch (i2)
1564 {
1565 // Cosine part P_k^C
1566 case 0:
1567
1568 jacobian(local_eqn, local_unknown) +=
1569 psip[l2] * (testf_ + r * dtestfdr) * W *
1570 hang_weight * hang_weight2;
1571
1572 break;
1573
1574 // Sine part P_k^S has no contribution
1575 case 1:
1576 break;
1577
1578 } // End of first (radial) momentum equation
1579 break;
1580
1581 // --------------------------------------------
1582 // SECOND (RADIAL) MOMENTUM EQUATION: SINE PART
1583 // --------------------------------------------
1584
1585 case 1:
1586
1587 switch (i2)
1588 {
1589 // Cosine part P_k^C has no contribution
1590 case 0:
1591 break;
1592
1593 // Sine part P_k^S
1594 case 1:
1595
1596 jacobian(local_eqn, local_unknown) +=
1597 psip[l2] * (testf_ + r * dtestfdr) * W *
1598 hang_weight * hang_weight2;
1599
1600 break;
1601
1602 } // End of second (radial) momentum equation
1603 break;
1604
1605 // --------------------------------------------
1606 // THIRD (AXIAL) MOMENTUM EQUATION: COSINE PART
1607 // --------------------------------------------
1608
1609 case 2:
1610
1611 switch (i2)
1612 {
1613 // Cosine part P_k^C
1614 case 0:
1615
1616 jacobian(local_eqn, local_unknown) +=
1617 r * psip[l2] * dtestfdz * W * hang_weight *
1618 hang_weight2;
1619
1620 break;
1621
1622 // Sine part P_k^S has no contribution
1623 case 1:
1624 break;
1625
1626 } // End of third (axial) momentum equation
1627 break;
1628
1629 // -------------------------------------------
1630 // FOURTH (AXIAL) MOMENTUM EQUATION: SINE PART
1631 // -------------------------------------------
1632
1633 case 3:
1634
1635 switch (i2)
1636 {
1637 // Cosine part P_k^C has no contribution
1638 case 0:
1639 break;
1640
1641 // Sine part P_k^S
1642 case 1:
1643
1644 jacobian(local_eqn, local_unknown) +=
1645 r * psip[l2] * dtestfdz * W * hang_weight *
1646 hang_weight2;
1647
1648 break;
1649
1650 } // End of fourth (axial) momentum equation
1651 break;
1652
1653 // ------------------------------------------------
1654 // FIFTH (AZIMUTHAL) MOMENTUM EQUATION: COSINE PART
1655 // ------------------------------------------------
1656
1657 case 4:
1658
1659 switch (i2)
1660 {
1661 // Cosine part P_k^C has no contribution
1662 case 0:
1663 break;
1664
1665 // Sine part P_k^S
1666 case 1:
1667
1668 jacobian(local_eqn, local_unknown) -=
1669 k * psip[l2] * testf_ * W * hang_weight *
1670 hang_weight2;
1671
1672 break;
1673
1674 } // End of fifth (azimuthal) momentum equation
1675 break;
1676
1677 // ----------------------------------------------
1678 // SIXTH (AZIMUTHAL) MOMENTUM EQUATION: SINE PART
1679 // ----------------------------------------------
1680
1681 case 5:
1682
1683 switch (i2)
1684 {
1685 // Cosine part P_k^C
1686 case 0:
1687
1688 jacobian(local_eqn, local_unknown) +=
1689 k * psip[l2] * testf_ * W * hang_weight *
1690 hang_weight2;
1691
1692 break;
1693
1694 // Sine part P_k^S has no contribution
1695 case 1:
1696 break;
1697 } // End of sixth (azimuthal) momentum equation
1698 } // End of add in contributions to different equations
1699 } // End of if unknown is not pinned statement
1700 } // End of loop over pressure components
1701 } // End of loop over master (m2) nodes
1702 } // End of loop over pressure "nodes"
1703 } // End of Jacobian calculation
1704 } // End of if not boundary condition statement
1705 } // End of loop over velocity components
1706 } // End of loop over master (m) nodes
1707 } // End of loop over nodes
1708
1709 // ====================
1710 // CONTINUITY EQUATIONS
1711 // ====================
1712
1713 // Loop over the pressure shape functions
1714 for (unsigned l = 0; l < n_pres; l++)
1715 {
1716 // Cache test function
1717 const double testp_ = testp[l];
1718
1719 // If the pressure dof is hanging
1720 if (pressure_dof_is_hanging[l])
1721 {
1722 // Pressure dof is hanging so it must be nodal-based
1723 hang_info_pt = pressure_node_pt(l)->hanging_pt(p_index[0]);
1724
1725 // Get the number of master nodes from the pressure node
1726 n_master = hang_info_pt->nmaster();
1727 }
1728 // Otherwise the node is its own master
1729 else
1730 {
1731 n_master = 1;
1732 }
1733
1734 // Loop over the master nodes
1735 for (unsigned m = 0; m < n_master; m++)
1736 {
1737 // Loop over the two pressure components
1738 for (unsigned i = 0; i < 2; i++)
1739 {
1740 // Get the equation number
1741 // -----------------------
1742
1743 // If the pressure dof is hanging
1744 if (pressure_dof_is_hanging[l])
1745 {
1746 // Get the equation number from the master node
1747 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1748 p_index[i]);
1749
1750 // Get the hang weight from the master node
1751 hang_weight = hang_info_pt->master_weight(m);
1752 }
1753 else
1754 {
1755 // Local equation number
1756 local_eqn = this->p_local_eqn(l, i);
1757
1758 // Node contributes with full weight
1759 hang_weight = 1.0;
1760 }
1761
1762 // If it's not a boundary condition...
1763 if (local_eqn >= 0)
1764 {
1765 switch (i)
1766 {
1767 // --------------------------------------
1768 // FIRST CONTINUITY EQUATION: COSINE PART
1769 // --------------------------------------
1770
1771 case 0:
1772
1773 // Gradient terms
1774 residuals[local_eqn] +=
1775 (interpolated_UC + r * interpolated_dUCdr +
1776 k * interpolated_VS + r * interpolated_dWCdz) *
1777 testp_ * W * hang_weight;
1778
1779 break;
1780
1781 // -------------------------------------
1782 // SECOND CONTINUITY EQUATION: SINE PART
1783 // -------------------------------------
1784
1785 case 1:
1786
1787 // Gradient terms
1788 residuals[local_eqn] +=
1789 (interpolated_US + r * interpolated_dUSdr -
1790 k * interpolated_VC + r * interpolated_dWSdz) *
1791 testp_ * W * hang_weight;
1792
1793 break;
1794 }
1795
1796 // ======================
1797 // Calculate the Jacobian
1798 // ======================
1799 if (flag)
1800 {
1801 // Number of master nodes
1802 unsigned n_master2 = 1;
1803
1804 // Storage for the weight of the shape function
1805 double hang_weight2 = 1.0;
1806
1807 // Loop over the velocity nodes for columns
1808 for (unsigned l2 = 0; l2 < n_node; l2++)
1809 {
1810 // Cache velocity shape functions and their derivatives
1811 const double psif_ = psif[l2];
1812 const double dpsifdr = dpsifdx(l2, 0);
1813 const double dpsifdz = dpsifdx(l2, 1);
1814
1815 // Local boolean to indicate whether the node is hanging
1816 bool is_node2_hanging = node_pt(l2)->is_hanging();
1817
1818 // If the node is hanging
1819 if (is_node2_hanging)
1820 {
1821 // Get the hanging pointer
1822 hang_info2_pt = node_pt(l2)->hanging_pt();
1823
1824 // Read out number of master nodes from hanging data
1825 n_master2 = hang_info2_pt->nmaster();
1826 }
1827 // Otherwise the node is its own master
1828 else
1829 {
1830 n_master2 = 1;
1831 }
1832
1833 // Loop over the master nodes
1834 for (unsigned m2 = 0; m2 < n_master2; m2++)
1835 {
1836 // Loop over the velocity components
1837 for (unsigned i2 = 0; i2 < 6; i2++)
1838 {
1839 // Get the number of the unknown
1840 // -----------------------------
1841
1842 // If the node is hanging
1843 if (is_node2_hanging)
1844 {
1845 // Get the equation number from the master node
1846 local_unknown = local_hang_eqn(
1847 hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
1848
1849 // Get the hang weight from the master node
1850 hang_weight2 = hang_info2_pt->master_weight(m2);
1851 }
1852 // If the node is not hanging
1853 else
1854 {
1855 // Local equation number
1856 local_unknown =
1857 this->nodal_local_eqn(l2, u_nodal_index[i2]);
1858
1859 // Node contributes with full weight
1860 hang_weight2 = 1.0;
1861 }
1862
1863 // If the unknown is not pinned
1864 if (local_unknown >= 0)
1865 {
1866 // Different results for i and i2
1867 switch (i)
1868 {
1869 // --------------------------------------
1870 // FIRST CONTINUITY EQUATION: COSINE PART
1871 // --------------------------------------
1872
1873 case 0:
1874
1875 switch (i2)
1876 {
1877 // Radial velocity component (cosine part) U_k^C
1878
1879 case 0:
1880
1881 jacobian(local_eqn, local_unknown) +=
1882 (psif_ + r * dpsifdr) * testp_ * W *
1883 hang_weight * hang_weight2;
1884
1885 break;
1886
1887 // Radial velocity component (sine part) U_k^S
1888 // has no contribution
1889
1890 case 1:
1891 break;
1892
1893 // Axial velocity component (cosine part) W_k^C
1894
1895 case 2:
1896
1897 jacobian(local_eqn, local_unknown) +=
1898 r * dpsifdz * testp_ * W * hang_weight *
1899 hang_weight2;
1900
1901 break;
1902
1903 // Axial velocity component (sine part) W_k^S
1904 // has no contribution
1905
1906 case 3:
1907 break;
1908
1909 // Azimuthal velocity component (cosine part)
1910 // V_k^C has no contribution
1911
1912 case 4:
1913 break;
1914
1915 // Azimuthal velocity component (sine part)
1916 // V_k^S
1917
1918 case 5:
1919
1920 jacobian(local_eqn, local_unknown) +=
1921 k * psif_ * testp_ * W * hang_weight *
1922 hang_weight2;
1923
1924 break;
1925
1926 } // End of first continuity equation
1927 break;
1928
1929 // -------------------------------------
1930 // SECOND CONTINUITY EQUATION: SINE PART
1931 // -------------------------------------
1932
1933 case 1:
1934
1935 switch (i2)
1936 {
1937 // Radial velocity component (cosine part) U_k^C
1938 // has no contribution
1939
1940 case 0:
1941 break;
1942
1943 // Radial velocity component (sine part) U_k^S
1944
1945 case 1:
1946
1947 jacobian(local_eqn, local_unknown) +=
1948 (psif_ + r * dpsifdr) * testp_ * W *
1949 hang_weight * hang_weight2;
1950
1951 break;
1952
1953 // Axial velocity component (cosine part) W_k^C
1954 // has no contribution
1955
1956 case 2:
1957 break;
1958
1959 // Axial velocity component (sine part) W_k^S
1960
1961 case 3:
1962
1963 jacobian(local_eqn, local_unknown) +=
1964 r * dpsifdz * testp_ * W * hang_weight *
1965 hang_weight2;
1966
1967 break;
1968
1969 // Azimuthal velocity component (cosine part)
1970 // V_k^C
1971
1972 case 4:
1973
1974 jacobian(local_eqn, local_unknown) -=
1975 k * psif_ * testp_ * W * hang_weight *
1976 hang_weight2;
1977
1978 break;
1979
1980 // Azimuthal velocity component (sine part)
1981 // V_k^S has no contribution
1982
1983 case 5:
1984 break;
1985
1986 } // End of second continuity equation
1987 break;
1988 }
1989 } // End of if unknown is not pinned statement
1990 } // End of loop over velocity components
1991 } // End of loop over master (m2) nodes
1992 } // End of loop over velocity nodes
1993
1994 // Real and imaginary pressure components, P_k^C and P_k^S,
1995 // have no contribution to Jacobian
1996
1997 } // End of Jacobian calculation
1998 } // End of if not boundary condition statement
1999 } // End of loop over the two pressure components
2000 } // End of loop over master nodes
2001 } // End of loop over pressure nodes
2002
2003 } // End of loop over integration points
2004
2005 } // End of fill_in_generic_residual_contribution_linearised_axi_nst
2006
2007} // End of namespace oomph
s
static char t char * s
Definition cfortran.h:568
i
cstr elem_len * i
Definition cfortran.h:603
oomph::DShape
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition shape.h:278
oomph::Data::time_stepper_pt
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition nodes.h:238
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition elements.h:1967
oomph::FiniteElement::nodal_value
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition elements.h:2597
oomph::FiniteElement::interpolated_x
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition elements.cc:3992
oomph::FiniteElement::nodal_local_eqn
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition elements.h:1436
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition elements.h:2214
oomph::FiniteElement::nodal_position
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition elements.h:2321
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
oomph::FiniteElement::dnodal_position_dt
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition elements.h:2337
oomph::HangInfo
Class that contains data for hanging nodes.
Definition nodes.h:742
oomph::HangInfo::master_node_pt
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition nodes.h:791
oomph::Integral::knot
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
oomph::Integral::nweight
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
oomph::Integral::weight
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
oomph::LinearisedAxisymmetricNavierStokesEquations::azimuthal_mode_number
const int & azimuthal_mode_number() const
Azimuthal mode number k in e^ik(theta) decomposition.
Definition linearised_axisym_navier_stokes_elements.h:236
oomph::LinearisedAxisymmetricNavierStokesEquations::re_st
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Definition linearised_axisym_navier_stokes_elements.h:230
oomph::LinearisedAxisymmetricNavierStokesEquations::Gamma
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
Definition linearised_axisym_navier_stokes_elements.h:218
oomph::LinearisedAxisymmetricNavierStokesEquations::du_dt_linearised_axi_nst
double du_dt_linearised_axi_nst(const unsigned &n, const unsigned &i) const
Return the i-th component of du/dt at local node n. Uses suitably interpolated value for hanging node...
Definition linearised_axisym_navier_stokes_elements.h:319
oomph::LinearisedAxisymmetricNavierStokesEquations::get_base_flow_u
virtual void get_base_flow_u(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result) const
Calculate the velocity components of the base flow solution at a given time and Eulerian position.
Definition linearised_axisym_navier_stokes_elements.h:139
oomph::LinearisedAxisymmetricNavierStokesEquations::p_local_eqn
virtual int p_local_eqn(const unsigned &n, const unsigned &i)=0
Access function for the local equation number information for the i-th component of the pressure....
oomph::LinearisedAxisymmetricNavierStokesEquations::ALE_is_disabled
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed....
Definition linearised_axisym_navier_stokes_elements.h:101
oomph::LinearisedAxisymmetricNavierStokesEquations::pshape_linearised_axi_nst
virtual void pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
oomph::LinearisedAxisymmetricNavierStokesEquations::get_base_flow_dudx
virtual void get_base_flow_dudx(const double &time, const unsigned &ipt, const Vector< double > &x, DenseMatrix< double > &result) const
Calculate the derivatives of the velocity components of the base flow solution w.r....
Definition linearised_axisym_navier_stokes_elements.h:163
oomph::LinearisedAxisymmetricNavierStokesEquations::dshape_and_dtest_eulerian_at_knot_linearised_axi_nst
virtual double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and their derivatives w.r.t. global coordinates at the ipt-th integration...
oomph::LinearisedAxisymmetricNavierStokesEquations::npres_linearised_axi_nst
virtual unsigned npres_linearised_axi_nst() const =0
Return the number of pressure degrees of freedom associated with a single pressure component in the e...
oomph::LinearisedAxisymmetricNavierStokesEquations::viscosity_ratio
const double & viscosity_ratio() const
Viscosity ratio for element: element's viscosity relative to the viscosity used in the definition of ...
Definition linearised_axisym_navier_stokes_elements.h:261
oomph::LinearisedAxisymmetricNavierStokesEquations::u_index_linearised_axi_nst
virtual unsigned u_index_linearised_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
Definition linearised_axisym_navier_stokes_elements.h:312
oomph::LinearisedAxisymmetricNavierStokesEquations::p_linearised_axi_nst
virtual double p_linearised_axi_nst(const unsigned &n_p, const unsigned &i) const =0
Return the i-th pressure value at local pressure "node" n_p. Uses suitably interpolated value for han...
oomph::LinearisedAxisymmetricNavierStokesEquations::p_index_linearised_axi_nst
virtual int p_index_linearised_axi_nst(const unsigned &i) const
Which nodal value represents the pressure?
Definition linearised_axisym_navier_stokes_elements.h:372
oomph::LinearisedAxisymmetricNavierStokesEquations::density_ratio
const double & density_ratio() const
Density ratio for element: element's density relative to the viscosity used in the definition of the ...
Definition linearised_axisym_navier_stokes_elements.h:274
oomph::Node::is_hanging
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition nodes.h:1285
oomph::Node::hanging_pt
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition nodes.h:1228
oomph::RefineableElement::local_hang_eqn
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
Definition refineable_elements.h:278
oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::pressure_node_pt
virtual Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interp...
Definition refineable_linearised_axisym_navier_stokes_elements.h:55
oomph::RefineableLinearisedAxisymmetricNavierStokesEquations::fill_in_generic_residual_contribution_linearised_axi_nst
void fill_in_generic_residual_contribution_linearised_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 flag=1: compute bo...
Definition refineable_linearised_axisym_navier_stokes_elements.cc:40
oomph::Shape
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition shape.h:76
oomph::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Definition Tadvection_diffusion_reaction_elements.h:66
oomph::TimeStepper::weight
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition timesteppers.h:594
oomph::TimeStepper::time
double & time()
Return current value of continous time.
Definition timesteppers.h:332
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition advection_diffusion_elements.cc:30
refineable_linearised_axisym_navier_stokes_elements.h