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-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // 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 
32 namespace 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::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
oomph::FiniteElement::nodal_value
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
oomph::FiniteElement::interpolated_x
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
oomph::FiniteElement::nodal_local_eqn
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
oomph::FiniteElement::nnode
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
oomph::FiniteElement::integral_pt
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
oomph::FiniteElement::nodal_position
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
oomph::FiniteElement::dnodal_position_dt
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition: elements.h:2333
oomph::HangInfo
Class that contains data for hanging nodes.
Definition: nodes.h:742
oomph::HangInfo::master_weight
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
oomph::HangInfo::master_node_pt
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
oomph::HangInfo::nmaster
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
oomph::Integral::knot
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
oomph::Integral::nweight
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
oomph::Integral::weight
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
oomph::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::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::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::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::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::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::Node::hanging_pt
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
oomph::Node::is_hanging
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
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::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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: advection_diffusion_elements.cc:30
refineable_linearised_axisym_navier_stokes_elements.h