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