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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 linearised axisymmetric Navier-Stokes elements
31 
32 // oomph-lib includes
34 
35 namespace oomph
36 {
37 
38  //=======================================================================
39  /// Linearised axisymmetric Navier--Stokes equations static data
40  //=======================================================================
41 
42  // Use the stress-divergence form by default (Gamma=1)
44 
45  // "Magic" number to indicate pressure is not stored at node
48 
49  // Physical constants default to zero
52 
53  // Azimuthal mode number defaults to zero
56 
57  // Density/viscosity ratios default to one
60 
61 
62 
63  //=======================================================================
64  /// Output function in tecplot format: Velocities only
65  /// r, z, U^C, U^S, V^C, V^S, W^C, W^S
66  /// at specified previous timestep (t=0 present; t>0 previous timestep).
67  /// Specified number of plot points in each coordinate direction.
68  //=======================================================================
70  std::ostream &outfile, const unsigned &nplot, const unsigned &t)
71  {
72  // Determine number of nodes in element
73  const unsigned n_node = nnode();
74 
75  // Provide storage for local shape functions
76  Shape psi(n_node);
77 
78  // Provide storage for vectors of local coordinates and
79  // global coordinates and velocities
80  Vector<double> s(2);
82  Vector<double> interpolated_u(6);
83 
84  // Tecplot header info
85  outfile << tecplot_zone_string(nplot);
86 
87  // Determine number of plot points
88  const unsigned n_plot_points = nplot_points(nplot);
89 
90  // Loop over plot points
91  for(unsigned iplot=0;iplot<n_plot_points;iplot++)
92  {
93  // Get local coordinates of plot point
94  get_s_plot(iplot,nplot,s);
95 
96  // Get shape functions
97  shape(s,psi);
98 
99  // Loop over coordinate directions
100  for(unsigned i=0;i<2;i++)
101  {
102  // Initialise global coordinate
103  interpolated_x[i]=0.0;
104 
105  // Loop over the local nodes and sum
106  for(unsigned l=0;l<n_node;l++)
107  {
108  interpolated_x[i] += nodal_position(t,l,i)*psi[l];
109  }
110  }
111 
112  // Loop over the velocity components
113  for(unsigned i=0;i<6;i++)
114  {
115  // Get the index at which the velocity is stored
116  const unsigned u_nodal_index = u_index_linearised_axi_nst(i);
117 
118  // Initialise velocity
119  interpolated_u[i]=0.0;
120 
121  // Loop over the local nodes and sum
122  for(unsigned l=0;l<n_node;l++)
123  {
124  interpolated_u[i] += nodal_value(t,l,u_nodal_index)*psi[l];
125  }
126  }
127 
128  // Output global coordinates to file
129  for(unsigned i=0;i<2;i++) { outfile << interpolated_x[i] << " "; }
130 
131  // Output velocities to file
132  for(unsigned i=0;i<6;i++) { outfile << interpolated_u[i] << " "; }
133 
134  outfile << std::endl;
135  }
136 
137  // Write tecplot footer (e.g. FE connectivity lists)
138  write_tecplot_zone_footer(outfile,nplot);
139 
140  } // End of output_veloc
141 
142 
143 
144  //=======================================================================
145  /// Output function in tecplot format:
146  /// r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
147  /// in tecplot format. Specified number of plot points in each
148  /// coordinate direction.
149  //=======================================================================
151  std::ostream &outfile, const unsigned &nplot)
152  {
153  // Provide storage for vector of local coordinates
154  Vector<double> s(2);
155 
156  // Tecplot header info
157  outfile << tecplot_zone_string(nplot);
158 
159  // Determine number of plot points
160  const unsigned n_plot_points = nplot_points(nplot);
161 
162  // Loop over plot points
163  for(unsigned iplot=0;iplot<n_plot_points;iplot++)
164  {
165  // Get local coordinates of plot point
166  get_s_plot(iplot,nplot,s);
167 
168  // Output global coordinates to file
169  for(unsigned i=0;i<2;i++) { outfile << interpolated_x(s,i) << " "; }
170 
171  // Output velocities to file
172  for(unsigned i=0;i<6;i++)
173  {
174  outfile << interpolated_u_linearised_axi_nst(s,i) << " ";
175  }
176 
177  // Output pressure to file
178  for(unsigned i=0;i<2;i++)
179  {
180  outfile << interpolated_p_linearised_axi_nst(s,i) << " ";
181  }
182 
183  outfile << std::endl;
184  }
185  outfile << std::endl;
186 
187  // Write tecplot footer (e.g. FE connectivity lists)
188  write_tecplot_zone_footer(outfile,nplot);
189 
190  } // End of output
191 
192 
193 
194 
195  //=======================================================================
196  /// Output function in tecplot format:
197  /// r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
198  /// Specified number of plot points in each coordinate direction.
199  //=======================================================================
201  FILE* file_pt, const unsigned &nplot)
202  {
203  // Provide storage for vector of local coordinates
204  Vector<double> s(2);
205 
206  // Tecplot header info
207  fprintf(file_pt,"%s ",tecplot_zone_string(nplot).c_str());
208 
209  // Determine number of plot points
210  const unsigned n_plot_points = nplot_points(nplot);
211 
212  // Loop over plot points
213  for(unsigned iplot=0;iplot<n_plot_points;iplot++)
214  {
215  // Get local coordinates of plot point
216  get_s_plot(iplot,nplot,s);
217 
218  // Output global coordinates to file
219  for(unsigned i=0;i<2;i++) { fprintf(file_pt,"%g ",interpolated_x(s,i)); }
220 
221  // Output velocities to file
222  for(unsigned i=0;i<6;i++)
223  {
224  fprintf(file_pt,"%g ",interpolated_u_linearised_axi_nst(s,i));
225  }
226 
227  // Output pressure to file
228  for(unsigned i=0;i<2;i++)
229  {
230  fprintf(file_pt,"%g ",interpolated_p_linearised_axi_nst(s,i));
231  }
232 
233  fprintf(file_pt,"\n");
234  }
235 
236  fprintf(file_pt,"\n");
237 
238  // Write tecplot footer (e.g. FE connectivity lists)
239  write_tecplot_zone_footer(file_pt,nplot);
240 
241  } // End of output
242 
243 
244 
245  //=======================================================================
246  /// Get strain-rate tensor: \f$ e_{ij} \f$ where
247  /// \f$ i,j = r,z,\theta \f$ (in that order). We evaluate this tensor
248  /// at a value of theta such that the product of theta and the azimuthal
249  /// mode number (k) gives \f$ \pi/4 \f$. Therefore
250  /// \f$ \cos(k \theta) = \sin(k \theta) = 1/\sqrt{2} \f$.
251  //=======================================================================
254  {
255 #ifdef PARANOID
256  if ((strainrate.ncol()!=3)||(strainrate.nrow()!=3))
257  {
258  std::ostringstream error_message;
259  error_message << "The strain rate has incorrect dimensions "
260  << strainrate.ncol() << " x "
261  << strainrate.nrow() << " Not 3" << std::endl;
262 
263  throw OomphLibError(
264  error_message.str(),
265  "LinearisedAxisymmetricNavierStokeEquations::strain_rate()",
266  OOMPH_EXCEPTION_LOCATION);
267  }
268 #endif
269 
270  // Determine number of nodes in element
271  const unsigned n_node = nnode();
272 
273  // Set up memory for the shape and test functions
274  Shape psi(n_node);
275  DShape dpsidx(n_node,2);
276 
277  // Call the derivatives of the shape functions
278  dshape_eulerian(s,psi,dpsidx);
279 
280  // Radius
281  double interpolated_r = 0.0;
282 
283  // Velocity components and their derivatives
284  double UC = 0.0, US = 0.0;
285  double dUCdr = 0.0, dUSdr = 0.0;
286  double dUCdz = 0.0, dUSdz = 0.0;
287  double WC = 0.0, WS = 0.0;
288  double dWCdr = 0.0, dWSdr = 0.0;
289  double dWCdz = 0.0, dWSdz = 0.0;
290  double VC = 0.0, VS = 0.0;
291  double dVCdr = 0.0, dVSdr = 0.0;
292  double dVCdz = 0.0, dVSdz = 0.0;
293 
294  // Get the local storage for the indices
295  unsigned u_nodal_index[6];
296  for(unsigned i=0;i<6;++i)
297  {
298  u_nodal_index[i] = u_index_linearised_axi_nst(i);
299  }
300 
301  // Loop over nodes to assemble velocities and their derivatives
302  // w.r.t. r and z (x_0 and x_1)
303  for(unsigned l=0;l<n_node;l++)
304  {
305  interpolated_r += nodal_position(l,0)*psi[l];
306 
307  UC += nodal_value(l,u_nodal_index[0])*psi[l];
308  US += nodal_value(l,u_nodal_index[1])*psi[l];
309  WC += nodal_value(l,u_nodal_index[2])*psi[l];
310  WS += nodal_value(l,u_nodal_index[3])*psi[l];
311  VC += nodal_value(l,u_nodal_index[4])*psi[l];
312  VS += nodal_value(l,u_nodal_index[4])*psi[l];
313 
314  dUCdr += nodal_value(l,u_nodal_index[0])*dpsidx(l,0);
315  dUSdr += nodal_value(l,u_nodal_index[1])*dpsidx(l,0);
316  dWCdr += nodal_value(l,u_nodal_index[2])*dpsidx(l,0);
317  dWSdr += nodal_value(l,u_nodal_index[3])*dpsidx(l,0);
318  dVCdr += nodal_value(l,u_nodal_index[4])*dpsidx(l,0);
319  dVSdr += nodal_value(l,u_nodal_index[5])*dpsidx(l,0);
320 
321  dUCdz += nodal_value(l,u_nodal_index[0])*dpsidx(l,1);
322  dUSdz += nodal_value(l,u_nodal_index[1])*dpsidx(l,1);
323  dWCdz += nodal_value(l,u_nodal_index[2])*dpsidx(l,1);
324  dWSdz += nodal_value(l,u_nodal_index[3])*dpsidx(l,1);
325  dVCdz += nodal_value(l,u_nodal_index[4])*dpsidx(l,1);
326  dVSdz += nodal_value(l,u_nodal_index[5])*dpsidx(l,1);
327  }
328 
329  // Cache azimuthal mode number
330  const int k = this->azimuthal_mode_number();
331 
332  // We wish to evaluate the strain-rate tensor at a value of theta
333  // such that k*theta = pi/4 radians. That way we pick up equal
334  // contributions from the real and imaginary parts of the velocities.
335  // sin(pi/4) = cos(pi/4) = 1/sqrt(2)
336  const double cosktheta = 1.0/sqrt(2);
337  const double sinktheta = cosktheta;
338 
339  // Assemble velocities and their derivatives w.r.t. r, z and theta
340  // from real and imaginary parts
341  const double ur = UC*cosktheta + US*sinktheta;
342  const double utheta = VC*cosktheta + VS*sinktheta;
343 
344  const double durdr = dUCdr*cosktheta + dUSdr*sinktheta;
345  const double durdz = dUCdz*cosktheta + dUSdz*sinktheta;
346  const double durdtheta = k*US*cosktheta - k*UC*sinktheta;
347 
348  const double duzdr = dWCdr*cosktheta + dWSdr*sinktheta;
349  const double duzdz = dWCdz*cosktheta + dWSdz*sinktheta;
350  const double duzdtheta = k*WS*cosktheta - k*WC*sinktheta;
351 
352  const double duthetadr = dVCdr*cosktheta + dVSdr*sinktheta;
353  const double duthetadz = dVCdz*cosktheta + dVSdz*sinktheta;
354  const double duthetadtheta = k*VS*cosktheta - k*VC*sinktheta;
355 
356  // Assign strain rates without negative powers of the radius
357  // and zero those with:
358  strainrate(0,0)=durdr;
359  strainrate(0,1)=0.5*(durdz+duzdr);
360  strainrate(1,0)=strainrate(0,1);
361  strainrate(0,2)=0.5*duthetadr;
362  strainrate(2,0)=strainrate(0,2);
363  strainrate(1,1)=duzdz;
364  strainrate(1,2)=0.5*duthetadz;
365  strainrate(2,1)=strainrate(1,2);
366  strainrate(2,2)=0.0;
367 
368  // Overwrite the strain rates with negative powers of the radius
369  // unless we're at the origin
370  if (std::abs(interpolated_r)>1.0e-16)
371  {
372  double inverse_radius=1.0/interpolated_r;
373  strainrate(0,2)=0.5*(duthetadr + inverse_radius*(durdtheta - utheta));
374  strainrate(2,0)=strainrate(0,2);
375  strainrate(2,2)=inverse_radius*(ur + duthetadtheta);
376  strainrate(1,2)=0.5*(duthetadz + inverse_radius*duzdtheta);
377  strainrate(2,1)=strainrate(1,2);
378  }
379 
380  } // End of strain_rate
381 
382 
383 
384  //=======================================================================
385  /// Compute the residuals for the linearised axisymmetric Navier--Stokes
386  /// equations; flag=1(or 0): do (or don't) compute the Jacobian as well.
387  //=======================================================================
390  Vector<double> &residuals,
391  DenseMatrix<double> &jacobian,
392  DenseMatrix<double> &mass_matrix,
393  unsigned flag)
394  {
395  // Get the time from the first node in the element
396  const double time = this->node_pt(0)->time_stepper_pt()->time();
397 
398  // Determine number of nodes in the element
399  const unsigned n_node = nnode();
400 
401  // Determine how many pressure values there are associated with
402  // a single pressure component
403  const unsigned n_pres = npres_linearised_axi_nst();
404 
405  // Get the nodal indices at which the velocity is stored
406  unsigned u_nodal_index[6];
407  for(unsigned i=0;i<6;++i)
408  {
409  u_nodal_index[i] = u_index_linearised_axi_nst(i);
410  }
411 
412  // Set up memory for the fluid shape and test functions
413  // Note that there are two spatial dimensions, r and z, in this problem
414  Shape psif(n_node), testf(n_node);
415  DShape dpsifdx(n_node,2), dtestfdx(n_node,2);
416 
417  // Set up memory for the pressure shape and test functions
418  Shape psip(n_pres), testp(n_pres);
419 
420  // Determine number of integration points
421  const unsigned n_intpt = integral_pt()->nweight();
422 
423  // Set up memory for the vector to hold local coordinates (two dimensions)
424  Vector<double> s(2);
425 
426  // Get physical variables from the element
427  // (Reynolds number must be multiplied by the density ratio)
428  const double scaled_re = re()*density_ratio();
429  const double scaled_re_st = re_st()*density_ratio();
430  const double visc_ratio = viscosity_ratio();
431  const int k = azimuthal_mode_number();
432 
433  // Integers used to store the local equation and unknown numbers
434  int local_eqn = 0, local_unknown = 0;
435 
436  // Loop over the integration points
437  for(unsigned ipt=0;ipt<n_intpt;ipt++)
438  {
439  // Assign values of the local coordinates s
440  for(unsigned i=0;i<2;i++) { s[i] = integral_pt()->knot(ipt,i); }
441 
442  // Get the integral weight
443  const double w = integral_pt()->weight(ipt);
444 
445  // Calculate the fluid shape and test functions, and their derivatives
446  // w.r.t. the global coordinates
447  const double J =
449  testf,dtestfdx);
450 
451  // Calculate the pressure shape and test functions
452  pshape_linearised_axi_nst(s,psip,testp);
453 
454  // Premultiply the weights and the Jacobian of the mapping between
455  // local and global coordinates
456  const double W = w*J;
457 
458  // Allocate storage for the position and the derivative of the
459  // mesh positions w.r.t. time
461  Vector<double> mesh_velocity(2,0.0);
462 
463  // Allocate storage for the velocity components (six of these)
464  // and their derivatives w.r.t. time
465  Vector<double> interpolated_u(6,0.0);
466  Vector<double> dudt(6,0.0);
467 
468  // Allocate storage for the pressure components (two of these)
469  Vector<double> interpolated_p(2,0.0);
470 
471  // Allocate storage for the derivatives of the velocity components
472  // w.r.t. global coordinates (r and z)
473  DenseMatrix<double> interpolated_dudx(6,2,0.0);
474 
475  // Calculate pressure at the integration point
476  // -------------------------------------------
477 
478  // Loop over pressure degrees of freedom (associated with a single
479  // pressure component) in the element
480  for(unsigned l=0;l<n_pres;l++)
481  {
482  // Cache the shape function
483  const double psip_ = psip(l);
484 
485  // Loop over the two pressure components
486  for(unsigned i=0;i<2;i++)
487  {
488  // Get the value
489  const double p_value = this->p_linearised_axi_nst(l,i);
490 
491  // Add contribution
492  interpolated_p[i] += p_value*psip_;
493  }
494  } // End of loop over the pressure degrees of freedom in the element
495 
496  // Calculate velocities and their derivatives at the integration point
497  // -------------------------------------------------------------------
498 
499  // Loop over the element's nodes
500  for(unsigned l=0;l<n_node;l++)
501  {
502  // Cache the shape function
503  const double psif_ = psif(l);
504 
505  // Loop over the two coordinate directions
506  for(unsigned i=0;i<2;i++)
507  {
508  interpolated_x[i] += this->raw_nodal_position(l,i)*psif_;
509  }
510 
511  // Loop over the six velocity components
512  for(unsigned i=0;i<6;i++)
513  {
514  // Get the value
515  const double u_value = this->raw_nodal_value(l,u_nodal_index[i]);
516 
517  // Add contribution
518  interpolated_u[i] += u_value*psif_;
519 
520  // Add contribution to dudt
521  dudt[i] += du_dt_linearised_axi_nst(l,i)*psif_;
522 
523  // Loop over two coordinate directions (for derivatives)
524  for(unsigned j=0;j<2;j++)
525  {
526  interpolated_dudx(i,j) += u_value*dpsifdx(l,j);
527  }
528  }
529  } // End of loop over the element's nodes
530 
531  // Get the mesh velocity if ALE is enabled
532  if(!ALE_is_disabled)
533  {
534  // Loop over the element's nodes
535  for(unsigned l=0;l<n_node;l++)
536  {
537  // Loop over the two coordinate directions
538  for(unsigned i=0;i<2;i++)
539  {
540  mesh_velocity[i] += this->raw_dnodal_position_dt(l,i)*psif(l);
541  }
542  }
543  }
544 
545  // Get velocities and their derivatives from base flow problem
546  // -----------------------------------------------------------
547 
548  // Allocate storage for the velocity components of the base state
549  // solution (initialise to zero)
550  Vector<double> base_flow_u(3,0.0);
551 
552  // Get the user-defined base state solution velocity components
553  get_base_flow_u(time,ipt,interpolated_x,base_flow_u);
554 
555  // Allocate storage for the derivatives of the base state solution's
556  // velocity components w.r.t. global coordinate (r and z)
557  // N.B. the derivatives of the base flow components w.r.t. the
558  // azimuthal coordinate direction (theta) are always zero since the
559  // base flow is axisymmetric
560  DenseMatrix<double> base_flow_dudx(3,2,0.0);
561 
562  // Get the derivatives of the user-defined base state solution
563  // velocity components w.r.t. global coordinates
564  get_base_flow_dudx(time,ipt,interpolated_x,base_flow_dudx);
565 
566  // Cache base flow velocities and their derivatives
567  const double interpolated_ur = base_flow_u[0];
568  const double interpolated_uz = base_flow_u[1];
569  const double interpolated_utheta = base_flow_u[2];
570  const double interpolated_durdr = base_flow_dudx(0,0);
571  const double interpolated_durdz = base_flow_dudx(0,1);
572  const double interpolated_duzdr = base_flow_dudx(1,0);
573  const double interpolated_duzdz = base_flow_dudx(1,1);
574  const double interpolated_duthetadr = base_flow_dudx(2,0);
575  const double interpolated_duthetadz = base_flow_dudx(2,1);
576 
577  // Cache r-component of position
578  const double r = interpolated_x[0];
579 
580  // Cache unknowns
581  const double interpolated_UC = interpolated_u[0];
582  const double interpolated_US = interpolated_u[1];
583  const double interpolated_WC = interpolated_u[2];
584  const double interpolated_WS = interpolated_u[3];
585  const double interpolated_VC = interpolated_u[4];
586  const double interpolated_VS = interpolated_u[5];
587  const double interpolated_PC = interpolated_p[0];
588  const double interpolated_PS = interpolated_p[1];
589 
590  // Cache derivatives of the unknowns
591  const double interpolated_dUCdr = interpolated_dudx(0,0);
592  const double interpolated_dUCdz = interpolated_dudx(0,1);
593  const double interpolated_dUSdr = interpolated_dudx(1,0);
594  const double interpolated_dUSdz = interpolated_dudx(1,1);
595  const double interpolated_dWCdr = interpolated_dudx(2,0);
596  const double interpolated_dWCdz = interpolated_dudx(2,1);
597  const double interpolated_dWSdr = interpolated_dudx(3,0);
598  const double interpolated_dWSdz = interpolated_dudx(3,1);
599  const double interpolated_dVCdr = interpolated_dudx(4,0);
600  const double interpolated_dVCdz = interpolated_dudx(4,1);
601  const double interpolated_dVSdr = interpolated_dudx(5,0);
602  const double interpolated_dVSdz = interpolated_dudx(5,1);
603 
604  // ==================
605  // MOMENTUM EQUATIONS
606  // ==================
607 
608  // Loop over the test functions
609  for(unsigned l=0;l<n_node;l++)
610  {
611  // Cache test functions and their derivatives
612  const double testf_ = testf(l);
613  const double dtestfdr = dtestfdx(l,0);
614  const double dtestfdz = dtestfdx(l,1);
615 
616  // ---------------------------------------------
617  // FIRST (RADIAL) MOMENTUM EQUATION: COSINE PART
618  // ---------------------------------------------
619 
620  // Get local equation number of first velocity value at this node
621  local_eqn = nodal_local_eqn(l,u_nodal_index[0]);
622 
623  // If it's not a boundary condition
624  if(local_eqn >= 0)
625  {
626  // Pressure gradient term
627  residuals[local_eqn] += interpolated_PC*(testf_ + r*dtestfdr)*W;
628 
629  // Stress tensor terms
630  residuals[local_eqn] -= visc_ratio*
631  r*(1.0+Gamma[0])*interpolated_dUCdr*dtestfdr*W;
632 
633  residuals[local_eqn] -= visc_ratio*r*
634  (interpolated_dUCdz + Gamma[0]*interpolated_dWCdr)*
635  dtestfdz*W;
636 
637  residuals[local_eqn] += visc_ratio*(
638  (k*Gamma[0]*interpolated_dVSdr)
639  - (k*(2.0+Gamma[0])*interpolated_VS/r)
640  - ((1.0+Gamma[0]+(k*k))*interpolated_UC/r))*testf_*W;
641 
642  // Inertial terms (du/dt)
643  residuals[local_eqn] -= scaled_re_st*r*dudt[0]*testf_*W;
644 
645  // Inertial terms (convective)
646  residuals[local_eqn] -=
647  scaled_re*(r*interpolated_ur*interpolated_dUCdr
648  + r*interpolated_durdr*interpolated_UC
649  + k*interpolated_utheta*interpolated_US
650  - 2*interpolated_utheta*interpolated_VC
651  + r*interpolated_uz*interpolated_dUCdz
652  + r*interpolated_durdz*interpolated_WC)*testf_*W;
653 
654  // Mesh velocity terms
655  if(!ALE_is_disabled)
656  {
657  for(unsigned j=0;j<2;j++)
658  {
659  residuals[local_eqn] +=
660  scaled_re_st*r*mesh_velocity[j]
661  *interpolated_dudx(0,j)*testf_*W;
662  }
663  }
664 
665  // Calculate the Jacobian
666  // ----------------------
667 
668  if(flag)
669  {
670  // Loop over the velocity shape functions again
671  for(unsigned l2=0;l2<n_node;l2++)
672  {
673  // Cache velocity shape functions and their derivatives
674  const double psif_ = psif[l2];
675  const double dpsifdr = dpsifdx(l2,0);
676  const double dpsifdz = dpsifdx(l2,1);
677 
678  // Radial velocity component (cosine part) U_k^C
679  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
680  if(local_unknown >= 0)
681  {
682  if(flag==2)
683  {
684  // Add the mass matrix
685  mass_matrix(local_eqn,local_unknown) +=
686  scaled_re_st*r*psif_*testf_*W;
687  }
688 
689  // Add contributions to the Jacobian matrix
690 
691  // Stress tensor terms
692  jacobian(local_eqn,local_unknown)
693  -= visc_ratio*r*(1.0+Gamma[0])*dpsifdr*dtestfdr*W;
694 
695  jacobian(local_eqn,local_unknown)
696  -= visc_ratio*r*dpsifdz*dtestfdz*W;
697 
698  jacobian(local_eqn,local_unknown)
699  -= visc_ratio*(1.0+Gamma[0]+(k*k))*psif_*testf_*W/r;
700 
701  // Inertial terms (du/dt)
702  jacobian(local_eqn,local_unknown) -=
703  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
704  psif_*testf_*W;
705 
706  // Inertial terms (convective)
707  jacobian(local_eqn,local_unknown) -=
708  scaled_re*r*(psif_*interpolated_durdr
709  + interpolated_ur*dpsifdr
710  + interpolated_uz*dpsifdz)*testf_*W;
711 
712  // Mesh velocity terms
713  if(!ALE_is_disabled)
714  {
715  for(unsigned j=0;j<2;j++)
716  {
717  jacobian(local_eqn,local_unknown) +=
718  scaled_re_st*r*mesh_velocity[j]
719  *dpsifdx(l2,j)*testf_*W;
720  }
721  }
722  }
723 
724  // Radial velocity component (sine part) U_k^S
725  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
726  if(local_unknown >= 0)
727  {
728  // Convective term
729  jacobian(local_eqn,local_unknown) -=
730  scaled_re*k*interpolated_utheta*psif_*testf_*W;
731  }
732 
733  // Axial velocity component (cosine part) W_k^C
734  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
735  if(local_unknown >= 0)
736  {
737  // Stress tensor term
738  jacobian(local_eqn,local_unknown) -=
739  visc_ratio*Gamma[0]*r*dpsifdr*dtestfdz*W;
740 
741  // Convective term
742  jacobian(local_eqn,local_unknown) -=
743  scaled_re*r*interpolated_durdz*psif_*testf_*W;
744  }
745 
746  // Axial velocity component (sine part) W_k^S
747  // has no contribution
748 
749  // Azimuthal velocity component (cosine part) V_k^C
750  local_unknown = nodal_local_eqn(l2,u_nodal_index[4]);
751  if(local_unknown >= 0)
752  {
753  // Convective term
754  jacobian(local_eqn,local_unknown) +=
755  scaled_re*2*interpolated_utheta*psif_*testf_*W;
756  }
757 
758  // Azimuthal velocity component (sine part) V_k^S
759  local_unknown = nodal_local_eqn(l2,u_nodal_index[5]);
760  if(local_unknown >= 0)
761  {
762  // Stress tensor term
763  jacobian(local_eqn,local_unknown) += visc_ratio*(
764  (Gamma[0]*k*dpsifdr) - (k*(2.0+Gamma[0])*psif_/r))*testf_*W;
765  }
766  } // End of loop over velocity shape functions
767 
768  // Now loop over pressure shape functions
769  // (This is the contribution from pressure gradient)
770  for(unsigned l2=0;l2<n_pres;l2++)
771  {
772  // Cosine part P_k^C
773  local_unknown = p_local_eqn(l2,0);
774  if(local_unknown >= 0)
775  {
776  jacobian(local_eqn,local_unknown) +=
777  psip[l2]*(testf_ + r*dtestfdr)*W;
778  }
779 
780  // Sine part P_k^S has no contribution
781 
782  } // End of loop over pressure shape functions
783  } // End of Jacobian calculation
784 
785  } // End of if not boundary condition statement
786 
787  // --------------------------------------------
788  // SECOND (RADIAL) MOMENTUM EQUATION: SINE PART
789  // --------------------------------------------
790 
791  // Get local equation number of second velocity value at this node
792  local_eqn = nodal_local_eqn(l,u_nodal_index[1]);
793 
794  // If it's not a boundary condition
795  if(local_eqn >= 0)
796  {
797  // Pressure gradient term
798  residuals[local_eqn] += interpolated_PS*(testf_ + r*dtestfdr)*W;
799 
800  // Stress tensor terms
801  residuals[local_eqn] -= visc_ratio*
802  r*(1.0+Gamma[0])*interpolated_dUSdr*dtestfdr*W;
803 
804  residuals[local_eqn] -= visc_ratio*r*
805  (interpolated_dUSdz + Gamma[0]*interpolated_dWSdr)*
806  dtestfdz*W;
807 
808  residuals[local_eqn] -= visc_ratio*(
809  (k*Gamma[0]*interpolated_dVCdr)
810  - (k*(2.0+Gamma[0])*interpolated_VC/r)
811  + ((1.0+Gamma[0]+(k*k))*interpolated_US/r))*testf_*W;
812 
813  // Inertial terms (du/dt)
814  residuals[local_eqn] -= scaled_re_st*r*dudt[1]*testf_*W;
815 
816  // Inertial terms (convective)
817  residuals[local_eqn] -=
818  scaled_re*(r*interpolated_ur*interpolated_dUSdr
819  + r*interpolated_durdr*interpolated_US
820  - k*interpolated_utheta*interpolated_UC
821  - 2*interpolated_utheta*interpolated_VS
822  + r*interpolated_uz*interpolated_dUSdz
823  + r*interpolated_durdz*interpolated_WS)*testf_*W;
824 
825  // Mesh velocity terms
826  if(!ALE_is_disabled)
827  {
828  for(unsigned j=0;j<2;j++)
829  {
830  residuals[local_eqn] +=
831  scaled_re_st*r*mesh_velocity[j]
832  *interpolated_dudx(1,j)*testf_*W;
833  }
834  }
835 
836  // Calculate the Jacobian
837  // ----------------------
838 
839  if(flag)
840  {
841  // Loop over the velocity shape functions again
842  for(unsigned l2=0;l2<n_node;l2++)
843  {
844  // Cache velocity shape functions and their derivatives
845  const double psif_ = psif[l2];
846  const double dpsifdr = dpsifdx(l2,0);
847  const double dpsifdz = dpsifdx(l2,1);
848 
849  // Radial velocity component (cosine part) U_k^C
850  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
851  if(local_unknown >= 0)
852  {
853  // Convective term
854  jacobian(local_eqn,local_unknown) +=
855  scaled_re*k*interpolated_utheta*psif_*testf_*W;
856  }
857 
858  // Radial velocity component (sine part) U_k^S
859  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
860  if(local_unknown >= 0)
861  {
862  if(flag==2)
863  {
864  // Add the mass matrix
865  mass_matrix(local_eqn,local_unknown) +=
866  scaled_re_st*r*psif_*testf_*W;
867  }
868 
869  // Add contributions to the Jacobian matrix
870 
871  // Stress tensor terms
872  jacobian(local_eqn,local_unknown)
873  -= visc_ratio*r*(1.0+Gamma[0])*dpsifdr*dtestfdr*W;
874 
875  jacobian(local_eqn,local_unknown)
876  -= visc_ratio*r*dpsifdz*dtestfdz*W;
877 
878  jacobian(local_eqn,local_unknown)
879  -= visc_ratio*(1.0+Gamma[0]+(k*k))*psif_*testf_*W/r;
880 
881  // Inertial terms (du/dt)
882  jacobian(local_eqn,local_unknown)
883  -= scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
884  psif_*testf_*W;
885 
886  // Inertial terms (convective)
887  jacobian(local_eqn,local_unknown) -=
888  scaled_re*r*(psif_*interpolated_durdr
889  + interpolated_ur*dpsifdr
890  + interpolated_uz*dpsifdz)*testf_*W;
891 
892  // Mesh velocity terms
893  if(!ALE_is_disabled)
894  {
895  for(unsigned j=0;j<2;j++)
896  {
897  jacobian(local_eqn,local_unknown) +=
898  scaled_re_st*r*mesh_velocity[j]
899  *dpsifdx(l2,j)*testf_*W;
900  }
901  }
902  }
903 
904  // Axial velocity component (cosine part) W_k^C
905  // has no contribution
906 
907  // Axial velocity component (sine part) W_k^S
908  local_unknown = nodal_local_eqn(l2,u_nodal_index[3]);
909  if(local_unknown >= 0)
910  {
911  // Stress tensor term
912  jacobian(local_eqn,local_unknown) -=
913  visc_ratio*Gamma[0]*r*dpsifdr*dtestfdz*W;
914 
915  // Convective term
916  jacobian(local_eqn,local_unknown) -=
917  scaled_re*r*interpolated_durdz*psif_*testf_*W;
918  }
919 
920  // Azimuthal velocity component (cosine part) V_k^C
921  local_unknown = nodal_local_eqn(l2,u_nodal_index[4]);
922  if(local_unknown >= 0)
923  {
924  // Stress tensor terms
925  jacobian(local_eqn,local_unknown) -= visc_ratio*(
926  (Gamma[0]*k*dpsifdr)
927  - (k*(2.0+Gamma[0])*psif_/r))*testf_*W;
928  }
929 
930  // Azimuthal velocity component (sine part) V_k^S
931  local_unknown = nodal_local_eqn(l2,u_nodal_index[5]);
932  if(local_unknown >= 0)
933  {
934  // Convective term
935  jacobian(local_eqn,local_unknown) +=
936  scaled_re*2*interpolated_utheta*psif_*testf_*W;
937  }
938  } // End of loop over velocity shape functions
939 
940  // Now loop over pressure shape functions
941  // (This is the contribution from pressure gradient)
942  for(unsigned l2=0;l2<n_pres;l2++)
943  {
944  // Cosine part P_k^C has no contribution
945 
946  // Sine part P_k^S
947  local_unknown = p_local_eqn(l2,1);
948  if(local_unknown >= 0)
949  {
950  jacobian(local_eqn,local_unknown)
951  += psip[l2]*(testf_ + r*dtestfdr)*W;
952  }
953  } // End of loop over pressure shape functions
954  } // End of Jacobian calculation
955 
956  } // End of if not boundary condition statement
957 
958  // --------------------------------------------
959  // THIRD (AXIAL) MOMENTUM EQUATION: COSINE PART
960  // --------------------------------------------
961 
962  // Get local equation number of third velocity value at this node
963  local_eqn = nodal_local_eqn(l,u_nodal_index[2]);
964 
965  // If it's not a boundary condition
966  if(local_eqn >= 0)
967  {
968  // Pressure gradient term
969  residuals[local_eqn] += r*interpolated_PC*dtestfdz*W;
970 
971  // Stress tensor terms
972  residuals[local_eqn] -= visc_ratio*r*
973  (interpolated_dWCdr + Gamma[1]*interpolated_dUCdz)
974  *dtestfdr*W;
975 
976  residuals[local_eqn] -= visc_ratio*r*
977  (1.0 + Gamma[1])*interpolated_dWCdz*dtestfdz*W;
978 
979  residuals[local_eqn] += visc_ratio*k*
980  (Gamma[1]*interpolated_dVSdz - k*interpolated_WC/r)*testf_*W;
981 
982  // Inertial terms (du/dt)
983  residuals[local_eqn] -= scaled_re_st*r*dudt[2]*testf_*W;
984 
985  // Inertial terms (convective)
986  residuals[local_eqn] -=
987  scaled_re*(r*interpolated_ur*interpolated_dWCdr
988  + r*interpolated_duzdr*interpolated_UC
989  + k*interpolated_utheta*interpolated_WS
990  + r*interpolated_uz*interpolated_dWCdz
991  + r*interpolated_duzdz*interpolated_WC)*testf_*W;
992 
993  // Mesh velocity terms
994  if(!ALE_is_disabled)
995  {
996  for(unsigned j=0;j<2;j++)
997  {
998  residuals[local_eqn] +=
999  scaled_re_st*r*mesh_velocity[j]
1000  *interpolated_dudx(2,j)*testf_*W;
1001  }
1002  }
1003 
1004  // Calculate the Jacobian
1005  // ----------------------
1006 
1007  if(flag)
1008  {
1009  // Loop over the velocity shape functions again
1010  for(unsigned l2=0;l2<n_node;l2++)
1011  {
1012  // Cache velocity shape functions and their derivatives
1013  const double psif_ = psif[l2];
1014  const double dpsifdr = dpsifdx(l2,0);
1015  const double dpsifdz = dpsifdx(l2,1);
1016 
1017  // Radial velocity component (cosine part) U_k^C
1018  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
1019  if(local_unknown >= 0)
1020  {
1021  // Stress tensor term
1022  jacobian(local_eqn,local_unknown) -=
1023  visc_ratio*r*Gamma[1]*dpsifdz*dtestfdr*W;
1024 
1025  // Convective term
1026  jacobian(local_eqn,local_unknown) -=
1027  scaled_re*r*psif_*interpolated_duzdr*testf_*W;
1028  }
1029 
1030  // Radial velocity component (sine part) U_k^S
1031  // has no contribution
1032 
1033  // Axial velocity component (cosine part) W_k^C
1034  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
1035  if(local_unknown >= 0)
1036  {
1037  if(flag==2)
1038  {
1039  // Add the mass matrix
1040  mass_matrix(local_eqn,local_unknown) +=
1041  scaled_re_st*r*psif_*testf_*W;
1042  }
1043 
1044  // Add contributions to the Jacobian matrix
1045 
1046  // Stress tensor terms
1047  jacobian(local_eqn,local_unknown) -=
1048  visc_ratio*r*dpsifdr*dtestfdr*W;
1049 
1050  jacobian(local_eqn,local_unknown) -=
1051  visc_ratio*r*(1.0+Gamma[1])*dpsifdz*dtestfdz*W;
1052 
1053  jacobian(local_eqn,local_unknown) -=
1054  visc_ratio*k*k*psif_*testf_*W/r;
1055 
1056  // Inertial terms (du/dt)
1057  jacobian(local_eqn,local_unknown) -=
1058  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1059  psif_*testf_*W;
1060 
1061  // Inertial terms (convective)
1062  jacobian(local_eqn,local_unknown) -=
1063  scaled_re*r*(interpolated_ur*dpsifdr
1064  + psif_*interpolated_duzdz
1065  + interpolated_uz*dpsifdz)*testf_*W;
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  }
1077  }
1078  }
1079 
1080  // Axial velocity component (sine part) W_k^S
1081  local_unknown = nodal_local_eqn(l2,u_nodal_index[3]);
1082  if(local_unknown >= 0)
1083  {
1084  // Convective term
1085  jacobian(local_eqn,local_unknown) -=
1086  scaled_re*k*interpolated_utheta*psif_*testf_*W;
1087  }
1088 
1089  // Azimuthal velocity component (cosine part) V_k^C
1090  // has no contribution
1091 
1092  // Azimuthal velocity component (sine part) V_k^S
1093  local_unknown = nodal_local_eqn(l2,u_nodal_index[5]);
1094  if(local_unknown >= 0)
1095  {
1096  // Stress tensor term
1097  jacobian(local_eqn,local_unknown) +=
1098  visc_ratio*Gamma[1]*k*dpsifdz*testf_*W;
1099  }
1100  } // End of loop over velocity shape functions
1101 
1102  // Now loop over pressure shape functions
1103  // (This is the contribution from pressure gradient)
1104  for(unsigned l2=0;l2<n_pres;l2++)
1105  {
1106  // Cosine part P_k^C
1107  local_unknown = p_local_eqn(l2,0);
1108  if(local_unknown >= 0)
1109  {
1110  jacobian(local_eqn,local_unknown) +=
1111  r*psip[l2]*dtestfdz*W;
1112  }
1113 
1114  // Sine part P_k^S has no contribution
1115 
1116  } // End of loop over pressure shape functions
1117  } // End of Jacobian calculation
1118 
1119  } // End of if not boundary condition statement
1120 
1121  // -------------------------------------------
1122  // FOURTH (AXIAL) MOMENTUM EQUATION: SINE PART
1123  // -------------------------------------------
1124 
1125  // Get local equation number of fourth velocity value at this node
1126  local_eqn = nodal_local_eqn(l,u_nodal_index[3]);
1127 
1128  // If it's not a boundary condition
1129  if(local_eqn >= 0)
1130  {
1131  // Pressure gradient term
1132  residuals[local_eqn] += r*interpolated_PS*dtestfdz*W;
1133 
1134  // Stress tensor terms
1135  residuals[local_eqn] -= visc_ratio*r*
1136  (interpolated_dWSdr + Gamma[1]*interpolated_dUSdz)
1137  *dtestfdr*W;
1138 
1139  residuals[local_eqn] -= visc_ratio*r*
1140  (1.0+Gamma[1])*interpolated_dWSdz*dtestfdz*W;
1141 
1142  residuals[local_eqn] -= visc_ratio*k*
1143  (Gamma[1]*interpolated_dVCdz + k*interpolated_WS/r)*testf_*W;
1144 
1145  // Inertial terms (du/dt)
1146  residuals[local_eqn] -= scaled_re_st*r*dudt[3]*testf_*W;
1147 
1148  // Inertial terms (convective)
1149  residuals[local_eqn] -=
1150  scaled_re*(r*interpolated_ur*interpolated_dWSdr
1151  + r*interpolated_duzdr*interpolated_US
1152  - k*interpolated_utheta*interpolated_WC
1153  + r*interpolated_uz*interpolated_dWSdz
1154  + r*interpolated_duzdz*interpolated_WS)*testf_*W;
1155 
1156  // Mesh velocity terms
1157  if(!ALE_is_disabled)
1158  {
1159  for(unsigned j=0;j<2;j++)
1160  {
1161  residuals[local_eqn] +=
1162  scaled_re_st*r*mesh_velocity[j]
1163  *interpolated_dudx(3,j)*testf_*W;
1164  }
1165  }
1166 
1167  // Calculate the Jacobian
1168  // ----------------------
1169 
1170  if(flag)
1171  {
1172  // Loop over the velocity shape functions again
1173  for(unsigned l2=0;l2<n_node;l2++)
1174  {
1175  // Cache velocity shape functions and their derivatives
1176  const double psif_ = psif[l2];
1177  const double dpsifdr = dpsifdx(l2,0);
1178  const double dpsifdz = dpsifdx(l2,1);
1179 
1180  // Radial velocity component (cosine part) U_k^C
1181  // has no contribution
1182 
1183  // Radial velocity component (sine part) U_k^S
1184  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
1185  if(local_unknown >= 0)
1186  {
1187  // Stress tensor term
1188  jacobian(local_eqn,local_unknown) -=
1189  visc_ratio*r*Gamma[1]*dpsifdz*dtestfdr*W;
1190 
1191  // Convective term
1192  jacobian(local_eqn,local_unknown) -=
1193  scaled_re*r*psif_*interpolated_duzdr*testf_*W;
1194  }
1195 
1196  // Axial velocity component (cosine part) W_k^S
1197  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
1198  if(local_unknown >= 0)
1199  {
1200  // Convective term
1201  jacobian(local_eqn,local_unknown) +=
1202  scaled_re*k*interpolated_utheta*psif_*testf_*W;
1203  }
1204 
1205  // Axial velocity component (sine part) W_k^S
1206  local_unknown = nodal_local_eqn(l2,u_nodal_index[3]);
1207  if(local_unknown >= 0)
1208  {
1209  if(flag==2)
1210  {
1211  // Add the mass matrix
1212  mass_matrix(local_eqn,local_unknown) +=
1213  scaled_re_st*r*psif_*testf_*W;
1214  }
1215 
1216  // Add contributions to the Jacobian matrix
1217 
1218  // Stress tensor terms
1219  jacobian(local_eqn,local_unknown) -=
1220  visc_ratio*r*dpsifdr*dtestfdr*W;
1221 
1222  jacobian(local_eqn,local_unknown) -=
1223  visc_ratio*r*(1.0+Gamma[1])*dpsifdz*dtestfdz*W;
1224 
1225  jacobian(local_eqn,local_unknown) -=
1226  visc_ratio*k*k*psif_*testf_*W/r;
1227 
1228  // Inertial terms (du/dt)
1229  jacobian(local_eqn,local_unknown) -=
1230  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1231  psif_*testf_*W;
1232 
1233  // Inertial terms (convective)
1234  jacobian(local_eqn,local_unknown) -=
1235  scaled_re*r*(interpolated_ur*dpsifdr
1236  + psif_*interpolated_duzdz
1237  + interpolated_uz*dpsifdz)*testf_*W;
1238 
1239 
1240  // Mesh velocity terms
1241  if(!ALE_is_disabled)
1242  {
1243  for(unsigned j=0;j<2;j++)
1244  {
1245  jacobian(local_eqn,local_unknown) +=
1246  scaled_re_st*r*mesh_velocity[j]
1247  *dpsifdx(l2,j)*testf_*W;
1248  }
1249  }
1250  }
1251 
1252  // Azimuthal velocity component (cosine part) V_k^C
1253  local_unknown = nodal_local_eqn(l2,u_nodal_index[4]);
1254  if(local_unknown >= 0)
1255  {
1256  // Stress tensor term
1257  jacobian(local_eqn,local_unknown) -=
1258  visc_ratio*Gamma[1]*k*dpsifdz*testf_*W;
1259  }
1260 
1261  // Azimuthal velocity component (sine part) V_k^S
1262  // has no contribution
1263 
1264  } // End of loop over velocity shape functions
1265 
1266  // Now loop over pressure shape functions
1267  // (This is the contribution from pressure gradient)
1268  for(unsigned l2=0;l2<n_pres;l2++)
1269  {
1270  // Cosine part P_k^C has no contribution
1271 
1272  // Sine part P_k^S
1273  local_unknown = p_local_eqn(l2,1);
1274  if(local_unknown >= 0)
1275  {
1276  jacobian(local_eqn,local_unknown) +=
1277  r*psip[l2]*dtestfdz*W;
1278  }
1279  } // End of loop over pressure shape functions
1280  } // End of Jacobian calculation
1281 
1282  } // End of if not boundary condition statement
1283 
1284  // ------------------------------------------------
1285  // FIFTH (AZIMUTHAL) MOMENTUM EQUATION: COSINE PART
1286  // ------------------------------------------------
1287 
1288  // Get local equation number of fifth velocity value at this node
1289  local_eqn = nodal_local_eqn(l,u_nodal_index[4]);
1290 
1291  // If it's not a boundary condition
1292  if(local_eqn >= 0)
1293  {
1294  // Pressure gradient term
1295  residuals[local_eqn] -= k*interpolated_PS*testf_*W;
1296 
1297  // Stress tensor terms
1298  residuals[local_eqn] += visc_ratio*
1299  (-r*interpolated_dVCdr
1300  - k*Gamma[0]*interpolated_US
1301  + Gamma[0]*interpolated_VC)*dtestfdr*W;
1302 
1303  residuals[local_eqn] -= visc_ratio*
1304  (k*Gamma[0]*interpolated_WS + r*interpolated_dVCdz)*dtestfdz*W;
1305 
1306  residuals[local_eqn] += visc_ratio*
1307  (Gamma[0]*interpolated_dVCdr
1308  + k*(2.0+Gamma[0])*interpolated_US/r
1309  - (1.0+(k*k)+(k*k*Gamma[0]))*interpolated_VC/r)*testf_*W;
1310 
1311  // Inertial terms (du/dt)
1312  residuals[local_eqn] -= scaled_re_st*r*dudt[4]*testf_*W;
1313 
1314  // Inertial terms (convective)
1315  residuals[local_eqn] -=
1316  scaled_re*(r*interpolated_ur*interpolated_dVCdr
1317  + r*interpolated_duthetadr*interpolated_UC
1318  + k*interpolated_utheta*interpolated_VS
1319  + interpolated_utheta*interpolated_UC
1320  + interpolated_ur*interpolated_VC
1321  + r*interpolated_uz*interpolated_dVCdz
1322  + r*interpolated_duthetadz*interpolated_WC)*testf_*W;
1323 
1324  // Mesh velocity terms
1325  if(!ALE_is_disabled)
1326  {
1327  for(unsigned j=0;j<2;j++)
1328  {
1329  residuals[local_eqn] +=
1330  scaled_re_st*r*mesh_velocity[j]
1331  *interpolated_dudx(4,j)*testf_*W;
1332  }
1333  }
1334 
1335  // Calculate the Jacobian
1336  // ----------------------
1337 
1338  if(flag)
1339  {
1340  // Loop over the velocity shape functions again
1341  for(unsigned l2=0;l2<n_node;l2++)
1342  {
1343  // Cache velocity shape functions and their derivatives
1344  const double psif_ = psif[l2];
1345  const double dpsifdr = dpsifdx(l2,0);
1346  const double dpsifdz = dpsifdx(l2,1);
1347 
1348  // Radial velocity component (cosine part) U_k^C
1349  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
1350  if(local_unknown >= 0)
1351  {
1352  // Convective terms
1353  jacobian(local_eqn,local_unknown) -=
1354  scaled_re*(r*interpolated_duthetadr
1355  + interpolated_utheta)*psif_*testf_*W;
1356  }
1357 
1358  // Radial velocity component (sine part) U_k^S
1359  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
1360  if(local_unknown >= 0)
1361  {
1362  // Stress tensor terms
1363  jacobian(local_eqn,local_unknown) += visc_ratio*k*psif_*(
1364  ((2.0+Gamma[0])*testf_/r) - (Gamma[0]*dtestfdr))*W;
1365  }
1366 
1367  // Axial velocity component (cosine part) W_k^C
1368  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
1369  if(local_unknown >= 0)
1370  {
1371  // Convective term
1372  jacobian(local_eqn,local_unknown) -=
1373  scaled_re*r*psif_*interpolated_duthetadz*testf_*W;
1374  }
1375 
1376  // Axial velocity component (sine part) W_k^S
1377  local_unknown = nodal_local_eqn(l2,u_nodal_index[3]);
1378  if(local_unknown >= 0)
1379  {
1380  // Stress tensor term
1381  jacobian(local_eqn,local_unknown) -=
1382  visc_ratio*k*Gamma[0]*psif_*dtestfdz*W;
1383  }
1384 
1385  // Azimuthal velocity component (cosine part) V_k^C
1386  local_unknown = nodal_local_eqn(l2,u_nodal_index[4]);
1387  if(local_unknown >= 0)
1388  {
1389  if(flag==2)
1390  {
1391  // Add the mass matrix
1392  mass_matrix(local_eqn,local_unknown) +=
1393  scaled_re_st*r*psif_*testf_*W;
1394  }
1395 
1396  // Add contributions to the Jacobian matrix
1397 
1398  // Stress tensor terms
1399  jacobian(local_eqn,local_unknown) +=
1400  visc_ratio*(Gamma[0]*psif_ - r*dpsifdr)*dtestfdr*W;
1401 
1402  jacobian(local_eqn,local_unknown) -=
1403  visc_ratio*r*dpsifdz*dtestfdz*W;
1404 
1405  jacobian(local_eqn,local_unknown) +=
1406  visc_ratio*(Gamma[0]*dpsifdr
1407  - (1.0+(k*k)+(k*k*Gamma[0]))*psif_/r)*testf_*W;
1408 
1409  // Inertial terms (du/dt)
1410  jacobian(local_eqn,local_unknown) -=
1411  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1412  psif_*testf_*W;
1413 
1414  // Inertial terms (convective)
1415  jacobian(local_eqn,local_unknown) -=
1416  scaled_re*(r*interpolated_ur*dpsifdr
1417  + interpolated_ur*psif_
1418  + r*interpolated_uz*dpsifdz)*testf_*W;
1419 
1420  // Mesh velocity terms
1421  if(!ALE_is_disabled)
1422  {
1423  for(unsigned j=0;j<2;j++)
1424  {
1425  jacobian(local_eqn,local_unknown) +=
1426  scaled_re_st*r*mesh_velocity[j]
1427  *dpsifdx(l2,j)*testf_*W;
1428  }
1429  }
1430  }
1431 
1432  // Azimuthal velocity component (sine part) V_k^S
1433  local_unknown = nodal_local_eqn(l2,u_nodal_index[5]);
1434  if(local_unknown >= 0)
1435  {
1436  // Convective term
1437  jacobian(local_eqn,local_unknown) -=
1438  scaled_re*k*interpolated_utheta*psif_*testf_*W;
1439  }
1440  } // End of loop over velocity shape functions
1441 
1442  // Now loop over pressure shape functions
1443  // (This is the contribution from pressure gradient)
1444  for(unsigned l2=0;l2<n_pres;l2++)
1445  {
1446  // Cosine part P_k^C has no contribution
1447 
1448  // Sine part P_k^S
1449  local_unknown = p_local_eqn(l2,1);
1450  if(local_unknown >= 0)
1451  {
1452  jacobian(local_eqn,local_unknown) -=
1453  k*psip[l2]*testf_*W;
1454  }
1455  } // End of loop over pressure shape functions
1456  } // End of Jacobian calculation
1457 
1458  } // End of if not boundary condition statement
1459 
1460  // ----------------------------------------------
1461  // SIXTH (AZIMUTHAL) MOMENTUM EQUATION: SINE PART
1462  // ----------------------------------------------
1463 
1464  // Get local equation number of sixth velocity value at this node
1465  local_eqn = nodal_local_eqn(l,u_nodal_index[5]);
1466 
1467  // If it's not a boundary condition
1468  if(local_eqn >= 0)
1469  {
1470  // Pressure gradient term
1471  residuals[local_eqn] += k*interpolated_PC*testf_*W;
1472 
1473  // Stress tensor terms
1474  residuals[local_eqn] += visc_ratio*
1475  (-r*interpolated_dVSdr
1476  + k*Gamma[0]*interpolated_UC
1477  + Gamma[0]*interpolated_VS)*dtestfdr*W;
1478 
1479  residuals[local_eqn] += visc_ratio*
1480  (k*Gamma[0]*interpolated_WC - r*interpolated_dVSdz)*dtestfdz*W;
1481 
1482  residuals[local_eqn] += visc_ratio*
1483  (Gamma[0]*interpolated_dVSdr
1484  - k*(2.0+Gamma[0])*interpolated_UC/r
1485  - (1.0+(k*k)+(k*k*Gamma[0]))*interpolated_VS/r)*testf_*W;
1486 
1487  // Inertial terms (du/dt)
1488  residuals[local_eqn] -= scaled_re_st*r*dudt[5]*testf_*W;
1489 
1490  // Inertial terms (convective)
1491  residuals[local_eqn] -=
1492  scaled_re*(r*interpolated_ur*interpolated_dVSdr
1493  + r*interpolated_duthetadr*interpolated_US
1494  - k*interpolated_utheta*interpolated_VC
1495  + interpolated_utheta*interpolated_US
1496  + interpolated_ur*interpolated_VS
1497  + r*interpolated_uz*interpolated_dVSdz
1498  + r*interpolated_duthetadz*interpolated_WS)*testf_*W;
1499 
1500  // Mesh velocity terms
1501  if(!ALE_is_disabled)
1502  {
1503  for(unsigned j=0;j<2;j++)
1504  {
1505  residuals[local_eqn] +=
1506  scaled_re_st*r*mesh_velocity[j]
1507  *interpolated_dudx(5,j)*testf_*W;
1508  }
1509  }
1510 
1511  // Calculate the Jacobian
1512  // ----------------------
1513 
1514  if(flag)
1515  {
1516  // Loop over the velocity shape functions again
1517  for(unsigned l2=0;l2<n_node;l2++)
1518  {
1519  // Cache velocity shape functions and their derivatives
1520  const double psif_ = psif[l2];
1521  const double dpsifdr = dpsifdx(l2,0);
1522  const double dpsifdz = dpsifdx(l2,1);
1523 
1524  // Radial velocity component (cosine part) U_k^C
1525  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
1526  if(local_unknown >= 0)
1527  {
1528  // Stress tensor terms
1529  jacobian(local_eqn,local_unknown) += visc_ratio*k*psif_*(
1530  (Gamma[0]*dtestfdr) - ((2.0+Gamma[0])*testf_/r))*W;
1531  }
1532 
1533  // Radial velocity component (sine part) U_k^S
1534  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
1535  if(local_unknown >= 0)
1536  {
1537  // Convective terms
1538  jacobian(local_eqn,local_unknown) -=
1539  scaled_re*(r*interpolated_duthetadr
1540  + interpolated_utheta)*psif_*testf_*W;
1541  }
1542 
1543  // Axial velocity component (cosine part) W_k^C
1544  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
1545  if(local_unknown >= 0)
1546  {
1547  // Stress tensor term
1548  jacobian(local_eqn,local_unknown) +=
1549  visc_ratio*k*Gamma[0]*psif_*dtestfdz*W;
1550  }
1551 
1552  // Axial velocity component (sine part) W_k^S
1553  local_unknown = nodal_local_eqn(l2,u_nodal_index[3]);
1554  if(local_unknown >= 0)
1555  {
1556  // Convective term
1557  jacobian(local_eqn,local_unknown) -=
1558  scaled_re*r*psif_*interpolated_duthetadz*testf_*W;
1559  }
1560 
1561  // Azimuthal velocity component (cosine part) V_k^C
1562  local_unknown = nodal_local_eqn(l2,u_nodal_index[4]);
1563  if(local_unknown >= 0)
1564  {
1565  // Convective term
1566  jacobian(local_eqn,local_unknown) +=
1567  scaled_re*k*interpolated_utheta*psif_*testf_*W;
1568  }
1569 
1570  // Azimuthal velocity component (sine part) V_k^S
1571  local_unknown = nodal_local_eqn(l2,u_nodal_index[5]);
1572  if(local_unknown >= 0)
1573  {
1574  if(flag==2)
1575  {
1576  // Add the mass matrix
1577  mass_matrix(local_eqn,local_unknown) +=
1578  scaled_re_st*r*psif_*testf_*W;
1579  }
1580 
1581  // Add contributions to the Jacobian matrix
1582 
1583  // Stress tensor terms
1584  jacobian(local_eqn,local_unknown) +=
1585  visc_ratio*(Gamma[0]*psif_ - r*dpsifdr)*dtestfdr*W;
1586 
1587  jacobian(local_eqn,local_unknown) -=
1588  visc_ratio*r*dpsifdz*dtestfdz*W;
1589 
1590  jacobian(local_eqn,local_unknown) +=
1591  visc_ratio*(Gamma[0]*dpsifdr
1592  - (1.0+(k*k)+(k*k*Gamma[0]))*psif_/r)*testf_*W;
1593 
1594  // Inertial terms (du/dt)
1595  jacobian(local_eqn,local_unknown) -=
1596  scaled_re_st*r*node_pt(l2)->time_stepper_pt()->weight(1,0)*
1597  psif_*testf_*W;
1598 
1599  // Convective terms
1600  jacobian(local_eqn,local_unknown) -=
1601  scaled_re*(r*interpolated_ur*dpsifdr
1602  + interpolated_ur*psif_
1603  + r*interpolated_uz*dpsifdz)*testf_*W;
1604 
1605  // Mesh velocity terms
1606  if(!ALE_is_disabled)
1607  {
1608  for(unsigned j=0;j<2;j++)
1609  {
1610  jacobian(local_eqn,local_unknown) +=
1611  scaled_re_st*r*mesh_velocity[j]
1612  *dpsifdx(l2,j)*testf_*W;
1613  }
1614  }
1615  }
1616  } // End of loop over velocity shape functions
1617 
1618  // Now loop over pressure shape functions
1619  // (This is the contribution from pressure gradient)
1620  for(unsigned l2=0;l2<n_pres;l2++)
1621  {
1622  // Cosine part P_k^C
1623  local_unknown = p_local_eqn(l2,0);
1624  if(local_unknown >= 0)
1625  {
1626  jacobian(local_eqn,local_unknown) +=
1627  k*psip[l2]*testf_*W;
1628  }
1629 
1630  // Sine part P_k^S has no contribution
1631 
1632  } // End of loop over pressure shape functions
1633  } // End of Jacobian calculation
1634 
1635  } // End of if not boundary condition statement
1636 
1637  } // End of loop over shape functions
1638 
1639 
1640  // ====================
1641  // CONTINUITY EQUATIONS
1642  // ====================
1643 
1644  // Loop over the pressure shape functions
1645  for(unsigned l=0;l<n_pres;l++)
1646  {
1647  // Cache test function
1648  const double testp_ = testp[l];
1649 
1650  // --------------------------------------
1651  // FIRST CONTINUITY EQUATION: COSINE PART
1652  // --------------------------------------
1653 
1654  // Get local equation number of first pressure value at this node
1655  local_eqn = p_local_eqn(l,0);
1656 
1657  // If it's not a boundary condition
1658  if(local_eqn >= 0)
1659  {
1660  // Gradient terms
1661  residuals[local_eqn] +=
1662  (interpolated_UC + r*interpolated_dUCdr
1663  + k*interpolated_VS + r*interpolated_dWCdz)*testp_*W;
1664 
1665  // Calculate the Jacobian
1666  // ----------------------
1667 
1668  if(flag)
1669  {
1670  // Loop over the velocity shape functions
1671  for(unsigned l2=0;l2<n_node;l2++)
1672  {
1673  // Cache velocity shape functions and their derivatives
1674  const double psif_ = psif[l2];
1675  const double dpsifdr = dpsifdx(l2,0);
1676  const double dpsifdz = dpsifdx(l2,1);
1677 
1678  // Radial velocity component (cosine part) U_k^C
1679  local_unknown = nodal_local_eqn(l2,u_nodal_index[0]);
1680  if(local_unknown >= 0)
1681  {
1682  jacobian(local_eqn,local_unknown) +=
1683  (psif_ + r*dpsifdr)*testp_*W;
1684  }
1685 
1686  // Radial velocity component (sine part) U_k^S
1687  // has no contribution
1688 
1689  // Axial velocity component (cosine part) W_k^C
1690  local_unknown = nodal_local_eqn(l2,u_nodal_index[2]);
1691  if(local_unknown >= 0)
1692  {
1693  jacobian(local_eqn,local_unknown) +=
1694  r*dpsifdz*testp_*W;
1695  }
1696 
1697  // Axial velocity component (sine part) W_k^S
1698  // has no contribution
1699 
1700  // Azimuthal velocity component (cosine part) V_k^C
1701  // has no contribution
1702 
1703  // Azimuthal velocity component (sine part) V_k^S
1704  local_unknown = nodal_local_eqn(l2,u_nodal_index[5]);
1705  if(local_unknown >= 0)
1706  {
1707  jacobian(local_eqn,local_unknown) +=
1708  k*psif_*testp_*W;
1709  }
1710  } // End of loop over velocity shape functions
1711 
1712  // Real and imaginary pressure components, P_k^C and P_k^S,
1713  // have no contribution
1714 
1715  } // End of Jacobian calculation
1716 
1717  } // End of if not boundary condition statement
1718 
1719  // -------------------------------------
1720  // SECOND CONTINUITY EQUATION: SINE PART
1721  // -------------------------------------
1722 
1723  // Get local equation number of second pressure value at this node
1724  local_eqn = p_local_eqn(l,1);
1725 
1726  // If it's not a boundary condition
1727  if(local_eqn >= 0)
1728  {
1729  // Gradient terms
1730  residuals[local_eqn] +=
1731  (interpolated_US + r*interpolated_dUSdr
1732  - k*interpolated_VC + r*interpolated_dWSdz)*testp_*W;
1733 
1734  // Calculate the Jacobian
1735  // ----------------------
1736 
1737  if(flag)
1738  {
1739  // Loop over the velocity shape functions
1740  for(unsigned l2=0;l2<n_node;l2++)
1741  {
1742  // Cache velocity shape functions and their derivatives
1743  const double psif_ = psif[l2];
1744  const double dpsifdr = dpsifdx(l2,0);
1745  const double dpsifdz = dpsifdx(l2,1);
1746 
1747  // Radial velocity component (cosine part) U_k^C
1748  // has no contribution
1749 
1750  // Radial velocity component (sine part) U_k^S
1751  local_unknown = nodal_local_eqn(l2,u_nodal_index[1]);
1752  if(local_unknown >= 0)
1753  {
1754  jacobian(local_eqn,local_unknown) +=
1755  (psif_ + r*dpsifdr)*testp_*W;
1756  }
1757 
1758  // Axial velocity component (cosine part) W_k^C
1759  // has no contribution
1760 
1761  // Axial velocity component (sine part) W_k^S
1762  local_unknown = nodal_local_eqn(l2,u_nodal_index[3]);
1763  if(local_unknown >= 0)
1764  {
1765  jacobian(local_eqn,local_unknown) +=
1766  r*dpsifdz*testp_*W;
1767  }
1768 
1769  // Azimuthal velocity component (cosine part) V_k^C
1770  local_unknown = nodal_local_eqn(l2,u_nodal_index[4]);
1771  if(local_unknown >= 0)
1772  {
1773  jacobian(local_eqn,local_unknown) -=
1774  k*psif_*testp_*W;
1775  }
1776 
1777  // Azimuthal velocity component (sine part) V_k^S
1778  // has no contribution
1779 
1780  } // End of loop over velocity shape functions
1781 
1782  // Real and imaginary pressure components, P_k^C and P_k^S,
1783  // have no contribution
1784 
1785  } // End of Jacobian calculation
1786 
1787  } // End of if not boundary condition statement
1788 
1789  } // End of loop over pressure shape functions
1790 
1791  } // End of loop over the integration points
1792 
1793  } // End of fill_in_generic_residual_contribution_linearised_axi_nst
1794 
1795 
1796 
1797 ///////////////////////////////////////////////////////////////////////////
1798 ///////////////////////////////////////////////////////////////////////////
1799 ///////////////////////////////////////////////////////////////////////////
1800 
1801 
1802  /// Linearised axisymmetric Crouzeix-Raviart elements
1803  /// -------------------------------------------------
1804 
1805 
1806  //=======================================================================
1807  /// Set the data for the number of variables at each node
1808  //=======================================================================
1810  Initial_Nvalue[9]={6,6,6,6,6,6,6,6,6};
1811 
1812 
1813 
1814  //========================================================================
1815  /// Number of values (pinned or dofs) required at node n
1816  //========================================================================
1818  required_nvalue(const unsigned &n) const { return Initial_Nvalue[n]; }
1819 
1820 
1821 
1822 ///////////////////////////////////////////////////////////////////////////
1823 ///////////////////////////////////////////////////////////////////////////
1824 ///////////////////////////////////////////////////////////////////////////
1825 
1826 
1827  /// Linearised axisymmetric Taylor-Hood elements
1828  /// --------------------------------------------
1829 
1830 
1831  //=======================================================================
1832  /// Set the data for the number of variables at each node
1833  //=======================================================================
1835  Initial_Nvalue[9]={8,6,8,6,6,6,8,6,8};
1836 
1837 
1838 
1839  //=======================================================================
1840  /// Set the data for the pressure conversion array
1841  //=======================================================================
1843  Pconv[4]={0,2,6,8};
1844 
1845 
1846 
1847 } // End of oomph namespace
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers...
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
virtual void fill_in_generic_residual_contribution_linearised_axi_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier-Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
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...
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
cstr elem_len * i
Definition: cfortran.h:607
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:2962
double & time()
Return current value of continous time.
Definition: timesteppers.h:324
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...
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2976
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...
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:2458
const int & azimuthal_mode_number() const
Azimuthal mode number k in e^ik(theta) decomposition.
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
char t
Definition: cfortran.h:572
const double & viscosity_ratio() const
Viscosity ratio for element: element's viscosity relative to the viscosity used in the definition of ...
void output_veloc(std::ostream &outfile, const unsigned &nplot, const unsigned &t)
Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, in tecplot format. nplot points in each coordina...
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...
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation. NOTE: Moved to cc file because of a possible compiler bug in gcc (yes, really!). The move to the cc file avoids inlining which appears to cause problems (only) when compiled with gcc and -O3; offensive "illegal read" is in optimised-out section of code and data that is allegedly illegal is readily readable (by other means) just before this function is called so I can't really see how we could possibly be responsible for this...
Definition: elements.cc:1657
e
Definition: cfortran.h:575
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:504
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
virtual void pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
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:3841
double interpolated_p_linearised_axi_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated pressure p[i] at local coordinate s...
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
static int Default_Azimuthal_Mode_Number_Value
Static default value for the azimuthal mode number (zero)
double interpolated_u_linearised_axi_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated velocity u[i] at local coordinate s...
static char t char * s
Definition: cfortran.h:572
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
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...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:507
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed...
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
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...
virtual unsigned required_nvalue(const unsigned &n) const
Return number of values (pinned or dofs) required at local node n.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all initialised to one)
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:1383
double raw_dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n. Do not use the hanging node repre...
Definition: elements.h:2170
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:2215
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
Definition: elements.h:2446
void output(std::ostream &outfile)
Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S in tecplot format. Default number of pl...
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
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.
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
const double & density_ratio() const
Density ratio for element: element's density relative to the viscosity used in the definition of the ...
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...