axisym_navier_stokes_elements.h
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: 1197 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-06-07 17:30:57 +0100 (Tue, 07 Jun 2016) $
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 //Header file for Navier Stokes elements
31 
32 #ifndef OOMPH_AXISYMMETRIC_NAVIER_STOKES_ELEMENTS_HEADER
33 #define OOMPH_AXISYMMETRIC_NAVIER_STOKES_ELEMENTS_HEADER
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37 #include <oomph-lib-config.h>
38 #endif
39 
40 //OOMPH-LIB headers
41 #include "../generic/Qelements.h"
42 #include "../generic/fsi.h"
43 #include "../generic/projection.h"
44 
45 namespace oomph
46 {
47 
48 
49 //======================================================================
50 /// A class for elements that solve the unsteady
51 /// axisymmetric Navier--Stokes equations in
52 /// cylindrical polar coordinates, \f$ x_0^* = r^*\f$ and \f$ x_1^* = z^* \f$
53 /// with \f$ \partial / \partial \theta = 0 \f$. We're solving for the
54 /// radial, axial and azimuthal (swirl) velocities,
55 /// \f$ u_0^* = u_r^*(r^*,z^*,t^*) = u^*(r^*,z^*,t^*),
56 /// \ u_1^* = u_z^*(r^*,z^*,t^*) = w^*(r^*,z^*,t^*)\f$ and
57 /// \f$ u_2^* = u_\theta^*(r^*,z^*,t^*) = v^*(r^*,z^*,t^*) \f$,
58 /// respectively, and the pressure \f$ p(r^*,z^*,t^*) \f$.
59 /// This class contains the generic maths -- any concrete
60 /// implementation must be derived from this.
61 ///
62 /// In dimensional form the axisymmetric Navier-Stokes equations are given
63 /// by the momentum equations (for the \f$ r^* \f$, \f$ z^* \f$ and \f$ \theta
64 /// \f$
65 /// directions, respectively)
66 /// \f[
67 /// \rho\left(\frac{\partial u^*}{\partial t^*} + {u^*}\frac{\partial
68 /// u^*}{\partial r^*} - \frac{{v^*}^2}{r^*}
69 /// + {w^*}\frac{\partial u^*}{\partial z^*} \right) =
70 /// B_r^*\left(r^*,z^*,t^*\right)+ \rho G_r^*+
71 /// \frac{1}{r^*}
72 /// \frac{\partial\left({r^*}\sigma_{rr}^*\right)}{\partial r^*}
73 /// - \frac{\sigma_{\theta\theta}^*}{r^*} +
74 /// \frac{\partial\sigma_{rz}^*}{\partial z^*},
75 /// \f]
76 /// \f[
77 /// \rho\left(\frac{\partial w^*}{\partial t^*} + {u^*}\frac{\partial
78 /// w^*}{\partial r^*} + {w^*}\frac{\partial
79 /// w^*}{\partial z^*} \right) =
80 /// B_z^*\left(r^*,z^*,t^*\right)+\rho G_z^*+
81 /// \frac{1}{r^*}\frac{\partial\left({r^*}\sigma_{zr}^*\right)}{\partial
82 /// r^*} + \frac{\partial\sigma_{zz}^*}{\partial z^*},
83 /// \f]
84 /// \f[
85 /// \rho\left(\frac{\partial v^*}{\partial t^*} +
86 /// {u^*}\frac{\partial v^*}{\partial r^*} +
87 /// \frac{u^* v^*}{r^*}
88 /// +{w^*}\frac{\partial v^*}{\partial z^*} \right)=
89 /// B_\theta^*\left(r^*,z^*,t^*\right)+ \rho G_\theta^*+
90 /// \frac{1}{r^*}\frac{\partial\left({r^*}\sigma_{\theta
91 /// r}^*\right)}{\partial r^*} + \frac{\sigma_{r\theta}^*}{r^*} +
92 /// \frac{\partial\sigma_{\theta z}^*}{\partial z^*},
93 /// \f]
94 /// and
95 /// \f[
96 /// \frac{1}{r^*}\frac{\partial\left(r^*u^*\right)}{\partial r^*} +
97 /// \frac{\partial w^*}{\partial z^*} = Q^*.
98 /// \f]
99 /// The dimensional, symmetric stress tensor is defined as:
100 /// \f[
101 /// \sigma_{rr}^* = -p^* + 2\mu\frac{\partial u^*}{\partial r^*},
102 /// \qquad
103 /// \sigma_{\theta\theta}^* = -p^* +2\mu\frac{u^*}{r^*},
104 /// \f]
105 /// \f[
106 /// \sigma_{zz}^* = -p^* + 2\mu\frac{\partial w^*}{\partial z^*},
107 /// \qquad
108 /// \sigma_{rz}^* = \mu\left(\frac{\partial u^*}{\partial z^*} +
109 /// \frac{\partial w^*}{\partial r^*}\right),
110 /// \f]
111 /// \f[
112 /// \sigma_{\theta r}^* = \mu r^*\frac{\partial}{\partial r^*}
113 /// \left(\frac{v^*}{r^*}\right),
114 /// \qquad
115 /// \sigma_{\theta z}^* = \mu\frac{\partial v^*}{\partial z^*}.
116 /// \f]
117 /// Here, the (dimensional) velocity components are denoted
118 /// by \f$ u^* \f$, \f$ w^* \f$
119 /// and \f$ v^* \f$ for the radial, axial and azimuthal velocities,
120 /// respectively, and we
121 /// have split the body force into two components: A constant
122 /// vector \f$ \rho \ G_i^* \f$ which typically represents gravitational
123 /// forces; and a variable body force, \f$ B_i^*(r^*,z^*,t^*) \f$.
124 /// \f$ Q^*(r^*,z^*,t^*) \f$ is a volumetric source term for the
125 /// continuity equation and is typically equal to zero.
126 /// \n\n
127 /// We non-dimensionalise the equations, using problem-specific reference
128 /// quantities for the velocity, \f$ U \f$, length, \f$ L \f$, and time,
129 /// \f$ T \f$, and scale the constant body force vector on the
130 /// gravitational acceleration, \f$ g \f$, so that
131 /// \f[
132 /// u^* = U\, u, \qquad
133 /// w^* = U\, w, \qquad
134 /// v^* = U\, v,
135 /// \f]
136 /// \f[
137 /// r^* = L\, r, \qquad
138 /// z^* = L\, z, \qquad
139 /// t^* = T\, t,
140 /// \f]
141 /// \f[
142 /// G_i^* = g\, G_i, \qquad
143 /// B_i^* = \frac{U\mu_{ref}}{L^2}\, B_i, \qquad
144 /// p^* = \frac{\mu_{ref} U}{L}\, p, \qquad
145 /// Q^* = \frac{U}{L}\, Q.
146 /// \f]
147 /// where we note that the pressure and the variable body force have
148 /// been non-dimensionalised on the viscous scale. \f$ \mu_{ref} \f$
149 /// and \f$ \rho_{ref} \f$ (used below) are reference values
150 /// for the fluid viscosity and density, respectively. In single-fluid
151 /// problems, they are identical to the viscosity \f$ \mu \f$ and
152 /// density \f$ \rho \f$ of the (one and only) fluid in the problem.
153 /// \n\n
154 /// The non-dimensional form of the axisymmetric Navier-Stokes equations
155 /// is then given by
156 /// \f[
157 /// R_{\rho} Re\left(St\frac{\partial u}{\partial t} + {u}\frac{\partial
158 /// u}{\partial r} - \frac{{v}^2}{r}
159 /// + {w}\frac{\partial u}{\partial z} \right) =
160 /// B_r\left(r,z,t\right)+ R_\rho \frac{Re}{Fr} G_r +
161 /// \frac{1}{r}
162 /// \frac{\partial\left({r}\sigma_{rr}\right)}{\partial r}
163 /// - \frac{\sigma_{\theta\theta}}{r} +
164 /// \frac{\partial\sigma_{rz}}{\partial z},
165 /// \f]
166 /// \f[
167 /// R_{\rho} Re\left(St\frac{\partial w}{\partial t} + {u}\frac{\partial
168 /// w}{\partial r} + {w}\frac{\partial
169 /// w}{\partial z} \right) =
170 /// B_z\left(r,z,t\right)+ R_\rho \frac{Re}{Fr} G_z+
171 /// \frac{1}{r}\frac{\partial\left({r}\sigma_{zr}\right)}{\partial
172 /// r} + \frac{\partial\sigma_{zz}}{\partial z},
173 /// \f]
174 /// \f[
175 /// R_{\rho} Re\left(St\frac{\partial v}{\partial t} +
176 /// {u}\frac{\partial v}{\partial r} +
177 /// \frac{u v}{r}
178 /// +{w}\frac{\partial v}{\partial z} \right)=
179 /// B_\theta\left(r,z,t\right)+ R_\rho \frac{Re}{Fr} G_\theta+
180 /// \frac{1}{r}\frac{\partial\left({r}\sigma_{\theta
181 /// r}\right)}{\partial r} + \frac{\sigma_{r\theta}}{r} +
182 /// \frac{\partial\sigma_{\theta z}}{\partial z},
183 /// \f]
184 /// and
185 /// \f[
186 /// \frac{1}{r}\frac{\partial\left(ru\right)}{\partial r} +
187 /// \frac{\partial w}{\partial z} = Q.
188 /// \f]
189 /// Here the non-dimensional, symmetric stress tensor is defined as:
190 /// \f[
191 /// \sigma_{rr} = -p + 2R_\mu \frac{\partial u}{\partial r},
192 /// \qquad
193 /// \sigma_{\theta\theta} = -p +2R_\mu \frac{u}{r},
194 /// \f]
195 /// \f[
196 /// \sigma_{zz} = -p + 2R_\mu \frac{\partial w}{\partial z},
197 /// \qquad
198 /// \sigma_{rz} = R_\mu \left(\frac{\partial u}{\partial z} +
199 /// \frac{\partial w}{\partial r}\right),
200 /// \f]
201 /// \f[
202 /// \sigma_{\theta r} = R_\mu r
203 /// \frac{\partial}{\partial r}\left(\frac{v}{r}\right),
204 /// \qquad
205 /// \sigma_{\theta z} = R_\mu \frac{\partial v}{\partial z}.
206 /// \f]
207 /// and the dimensionless parameters
208 /// \f[
209 /// Re = \frac{UL\rho_{ref}}{\mu_{ref}}, \qquad
210 /// St = \frac{L}{UT}, \qquad
211 /// Fr = \frac{U^2}{gL},
212 /// \f]
213 /// are the Reynolds number, Strouhal number and Froude number
214 /// respectively. \f$ R_\rho=\rho/\rho_{ref} \f$ and
215 /// \f$ R_\mu(T) =\mu(T)/\mu_{ref}\f$ represent the ratios
216 /// of the fluid's density and its dynamic viscosity, relative to the
217 /// density and viscosity values used to form the non-dimensional
218 /// parameters (By default, \f$ R_\rho = R_\mu = 1 \f$; other values
219 /// tend to be used in problems involving multiple fluids).
220 //======================================================================
223 {
224  private:
225 
226  /// \short Static "magic" number that indicates that the pressure is
227  /// not stored at a node
229 
230  /// Static default value for the physical constants (all initialised to zero)
232 
233  /// Static default value for the physical ratios (all are initialised to one)
235 
236  /// Static default value for the gravity vector
238 
239 protected:
240 
241  //Physical constants
242 
243  /// \short Pointer to the viscosity ratio (relative to the
244  /// viscosity used in the definition of the Reynolds number)
246 
247  /// \short Pointer to the density ratio (relative to the density used in the
248  /// definition of the Reynolds number)
250 
251  // Pointers to global physical constants
252 
253  /// Pointer to global Reynolds number
254  double *Re_pt;
255 
256  /// Pointer to global Reynolds number x Strouhal number (=Womersley)
257  double *ReSt_pt;
258 
259  /// \short Pointer to global Reynolds number x inverse Froude number
260  /// (= Bond number / Capillary number)
261  double *ReInvFr_pt;
262 
263  /// \short Pointer to global Reynolds number x inverse Rossby number
264  /// (used when in a rotating frame)
265  double *ReInvRo_pt;
266 
267  /// Pointer to global gravity Vector
269 
270  /// Pointer to body force function
271  void (*Body_force_fct_pt)(const double& time,
272  const Vector<double> &x,
273  Vector<double> &result);
274 
275  /// Pointer to volumetric source function
276  double (*Source_fct_pt)(const double& time,
277  const Vector<double> &x);
278 
279  /// \short Boolean flag to indicate if ALE formulation is disabled when
280  /// the time-derivatives are computed. Only set to true if you're sure
281  /// that the mesh is stationary
283 
284  /// \short Access function for the local equation number information for
285  /// the pressure.
286  /// p_local_eqn[n] = local equation number or < 0 if pinned
287  virtual int p_local_eqn(const unsigned &n)const=0;
288 
289  /// \short Compute the shape functions and derivatives
290  /// w.r.t. global coords at local coordinate s.
291  /// Return Jacobian of mapping between local and global coordinates.
292  virtual double dshape_and_dtest_eulerian_axi_nst(const Vector<double> &s,
293  Shape &psi,
294  DShape &dpsidx, Shape &test,
295  DShape &dtestdx) const=0;
296 
297  /// \short Compute the shape functions and derivatives
298  /// w.r.t. global coords at ipt-th integration point
299  /// Return Jacobian of mapping between local and global coordinates.
300  virtual double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt,
301  Shape &psi,
302  DShape &dpsidx,
303  Shape &test,
304  DShape &dtestdx)
305  const=0;
306 
307  /// \short Shape/test functions and derivs w.r.t. to global coords at
308  /// integration point ipt; return Jacobian of mapping (J). Also compute
309  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
311  const unsigned &ipt,
312  Shape &psi,
313  DShape &dpsidx,
314  RankFourTensor<double> &d_dpsidx_dX,
315  Shape &test,
316  DShape &dtestdx,
317  RankFourTensor<double> &d_dtestdx_dX,
318  DenseMatrix<double> &djacobian_dX) const=0;
319 
320  /// Compute the pressure shape functions at local coordinate s
321  virtual void pshape_axi_nst(const Vector<double> &s, Shape &psi) const=0;
322 
323  /// \short Compute the pressure shape and test functions
324  /// at local coordinate s
325  virtual void pshape_axi_nst(const Vector<double> &s, Shape &psi,
326  Shape &test) const=0;
327 
328  /// Calculate the body force fct at a given time and Eulerian position
329  virtual void get_body_force_axi_nst(const double& time,
330  const unsigned& ipt,
331  const Vector<double> &s,
332  const Vector<double> &x,
333  Vector<double> &result)
334  {
335  //If the function pointer is zero return zero
336  if(Body_force_fct_pt == 0)
337  {
338  //Loop over dimensions and set body forces to zero
339  for(unsigned i=0;i<3;i++) {result[i] = 0.0;}
340  }
341  //Otherwise call the function
342  else
343  {
344  (*Body_force_fct_pt)(time,x,result);
345  }
346  }
347 
348  /// \short Get gradient of body force term at (Eulerian) position x.
349  /// Computed via function pointer (if set) or by finite differencing
350  /// (default)
351  inline virtual void get_body_force_gradient_axi_nst(
352  const double& time,
353  const unsigned& ipt,
354  const Vector<double> &s,
355  const Vector<double> &x,
356  DenseMatrix<double> &d_body_force_dx)
357  {
358 // hierher: Implement function pointer version
359 /* //If no gradient function has been set, FD it */
360 /* if(Body_force_fct_gradient_pt==0) */
361  {
362  // Reference value
363  Vector<double> body_force(3,0.0);
364  get_body_force_axi_nst(time,ipt,s,x,body_force);
365 
366  // FD it
368  Vector<double> body_force_pls(3,0.0);
369  Vector<double> x_pls(x);
370  for(unsigned i=0;i<2;i++)
371  {
372  x_pls[i] += eps_fd;
373  get_body_force_axi_nst(time,ipt,s,x_pls,body_force_pls);
374  for(unsigned j=0;j<3;j++)
375  {
376  d_body_force_dx(j,i)=(body_force_pls[j]-body_force[j])/eps_fd;
377  }
378  x_pls[i]=x[i];
379  }
380  }
381 /* else */
382 /* { */
383 /* // Get gradient */
384 /* (*Source_fct_gradient_pt)(time,x,gradient); */
385 /* } */
386  }
387 
388  /// Calculate the source fct at given time and Eulerian position
389  double get_source_fct(const double& time,
390  const unsigned& ipt,
391  const Vector<double> &x)
392  {
393  // If the function pointer is zero return zero
394  if(Source_fct_pt == 0) { return 0; }
395 
396  // Otherwise call the function
397  else { return (*Source_fct_pt)(time,x); }
398  }
399 
400  /// Get gradient of source term at (Eulerian) position x. Computed
401  /// via function pointer (if set) or by finite differencing (default)
402  inline virtual void get_source_fct_gradient(
403  const double& time,
404  const unsigned& ipt,
405  const Vector<double>& x,
406  Vector<double>& gradient)
407  {
408 // hierher: Implement function pointer version
409 /* //If no gradient function has been set, FD it */
410 /* if(Source_fct_gradient_pt==0) */
411  {
412  // Reference value
413  const double source = get_source_fct(time,ipt,x);
414 
415  // FD it
416  const double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
417  double source_pls = 0.0;
418  Vector<double> x_pls(x);
419  for(unsigned i=0;i<2;i++)
420  {
421  x_pls[i] += eps_fd;
422  source_pls = get_source_fct(time,ipt,x_pls);
423  gradient[i] = (source_pls-source)/eps_fd;
424  x_pls[i] = x[i];
425  }
426  }
427 /* else */
428 /* { */
429 /* // Get gradient */
430 /* (*Source_fct_gradient_pt)(time,x,gradient); */
431 /* } */
432  }
433 
434  /// \short Calculate the viscosity ratio relative to the
435  /// viscosity used in the definition of the Reynolds number
436  /// at given time and Eulerian position
437  virtual void get_viscosity_ratio_axisym_nst(const unsigned& ipt,
438  const Vector<double> &s,
439  const Vector<double> &x,
440  double &visc_ratio)
441  {
442  visc_ratio = *Viscosity_Ratio_pt;
443  }
444 
445  ///\short Compute the residuals for the Navier--Stokes equations;
446  /// flag=2 or 1 or 0: compute the Jacobian and/or mass matrix
447  /// as well.
449  Vector<double> &residuals,
450  DenseMatrix<double> &jacobian,
451  DenseMatrix<double> &mass_matrix,
452  unsigned flag);
453 
454  ///\short Compute the derivative of residuals for the
455  /// Navier--Stokes equations; with respect to a parameeter
456  /// flag=2 or 1 or 0: compute the Jacobian and/or mass matrix as well
458  double* const &parameter_pt,
459  Vector<double> &dres_dparam,
460  DenseMatrix<double> &djac_dparam,
461  DenseMatrix<double> &dmass_matrix_dparam,
462  unsigned flag);
463 
464  /// \short Compute the hessian tensor vector products required to
465  /// perform continuation of bifurcations analytically
467  Vector<double> const &Y,
468  DenseMatrix<double> const &C,
469  DenseMatrix<double> &product);
470 
471 public:
472 
473  /// \short Constructor: NULL the body force and source function
475  Body_force_fct_pt(0),
476  Source_fct_pt(0),
477  ALE_is_disabled(false)
478  {
479  //Set all the Physical parameter pointers to the default value zero
485  //Set the Physical ratios to the default value of 1
488  }
489 
490  /// Vector to decide whether the stress-divergence form is used or not
491  // N.B. This needs to be public so that the intel compiler gets things correct
492  // somehow the access function messes things up when going to refineable
493  // navier--stokes
495 
496  //Access functions for the physical constants
497 
498  /// Reynolds number
499  const double &re() const {return *Re_pt;}
500 
501  /// Product of Reynolds and Strouhal number (=Womersley number)
502  const double &re_st() const {return *ReSt_pt;}
503 
504  /// Pointer to Reynolds number
505  double* &re_pt() {return Re_pt;}
506 
507  /// Pointer to product of Reynolds and Strouhal number (=Womersley number)
508  double* &re_st_pt() {return ReSt_pt;}
509 
510  /// Global inverse Froude number
511  const double &re_invfr() const {return *ReInvFr_pt;}
512 
513  /// Pointer to global inverse Froude number
514  double* &re_invfr_pt() {return ReInvFr_pt;}
515 
516  /// Global Reynolds number multiplied by inverse Rossby number
517  const double &re_invro() const {return *ReInvRo_pt;}
518 
519  /// Pointer to global inverse Froude number
520  double* &re_invro_pt() {return ReInvRo_pt;}
521 
522  /// Vector of gravitational components
523  const Vector<double> &g() const {return *G_pt;}
524 
525  /// Pointer to Vector of gravitational components
526  Vector<double>* &g_pt() {return G_pt;}
527 
528  /// \short Density ratio for element: Element's density relative to the
529  /// viscosity used in the definition of the Reynolds number
530  const double &density_ratio() const {return *Density_Ratio_pt;}
531 
532  /// Pointer to Density ratio
533  double* &density_ratio_pt() {return Density_Ratio_pt;}
534 
535  /// \short Viscosity ratio for element: Element's viscosity relative to the
536  /// viscosity used in the definition of the Reynolds number
537  const double &viscosity_ratio() const {return *Viscosity_Ratio_pt;}
538 
539  /// Pointer to Viscosity Ratio
541 
542  /// Access function for the body-force pointer
543  void (* &axi_nst_body_force_fct_pt())(const double& time,
544  const Vector<double>& x,
545  Vector<double> & f)
546  {return Body_force_fct_pt;}
547 
548  ///Access function for the source-function pointer
549  double (* &source_fct_pt())(const double& time,
550  const Vector<double>& x)
551  {return Source_fct_pt;}
552 
553  ///Function to return number of pressure degrees of freedom
554  virtual unsigned npres_axi_nst() const=0;
555 
556  /// \short Return the index at which the i-th unknown velocity component
557  /// is stored. The default value, i, is appropriate for single-physics
558  /// problems.
559  /// In derived multi-physics elements, this function should be overloaded
560  /// to reflect the chosen storage scheme. Note that these equations require
561  /// that the unknowns are always stored at the same indices at each node.
562  virtual inline unsigned u_index_axi_nst(const unsigned &i) const {return i;}
563 
564  /// \short Return the index at which the i-th unknown velocity component
565  /// is stored with a common interface for use in general
566  /// FluidInterface and similar elements.
567  /// To do: Merge all common storage etc to a common base class for
568  /// Navier--Stokes elements in all coordinate systems.
569  inline unsigned u_index_nst(const unsigned &i) const
570  {return this->u_index_axi_nst(i);}
571 
572  /// \short Return the number of velocity components for use in
573  /// general FluidInterface clas
574  inline unsigned n_u_nst() const {return 3;}
575 
576 
577  /// \short i-th component of du/dt at local node n.
578  /// Uses suitably interpolated value for hanging nodes.
579  double du_dt_axi_nst(const unsigned &n, const unsigned &i) const
580  {
581  // Get the data's timestepper
583 
584  //Initialise dudt
585  double dudt=0.0;
586  //Loop over the timesteps, if there is a non Steady timestepper
587  if (!time_stepper_pt->is_steady())
588  {
589  //Get the index at which the velocity is stored
590  const unsigned u_nodal_index = u_index_axi_nst(i);
591 
592  // Number of timsteps (past & present)
593  const unsigned n_time = time_stepper_pt->ntstorage();
594 
595  //Add the contributions to the time derivative
596  for(unsigned t=0;t<n_time;t++)
597  {
598  dudt+=time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
599  }
600  }
601 
602  return dudt;
603  }
604 
605  /// \short Disable ALE, i.e. assert the mesh is not moving -- you do this
606  /// at your own risk!
607  void disable_ALE()
608  {
609  ALE_is_disabled=true;
610  }
611 
612  /// \short (Re-)enable ALE, i.e. take possible mesh motion into account
613  /// when evaluating the time-derivative. Note: By default, ALE is
614  /// enabled, at the expense of possibly creating unnecessary work
615  /// in problems where the mesh is, in fact, stationary.
616  void enable_ALE()
617  {
618  ALE_is_disabled=false;
619  }
620 
621  /// \short Pressure at local pressure "node" n_p
622  /// Uses suitably interpolated value for hanging nodes.
623  virtual double p_axi_nst(const unsigned &n_p)const=0;
624 
625  /// \short Which nodal value represents the pressure?
626  virtual int p_nodal_index_axi_nst() const
628 
629  /// Integral of pressure over element
630  double pressure_integral() const;
631 
632  /// \short Return integral of dissipation over element
633  double dissipation() const;
634 
635  /// \short Return dissipation at local coordinate s
636  double dissipation(const Vector<double>& s) const;
637 
638 /// \short Get integral of kinetic energy over element
639  double kin_energy() const;
640 
641  /// \short Strain-rate tensor: \f$ e_{ij} \f$ where \f$ i,j = r,z,\theta \f$
642  /// (in that order)
643  void strain_rate(const Vector<double>& s,
645 
646  /// \short Compute traction (on the viscous scale) at local coordinate s
647  /// for outer unit normal N
648  void traction(const Vector<double>& s,
649  const Vector<double>& N,
650  Vector<double>& traction) const;
651 
652  /// \short Compute the diagonal of the velocity/pressure mass matrices.
653  /// If which one=0, both are computed, otherwise only the pressure
654  /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
655  /// LSC version of the preconditioner only needs that one)
656  /// NOTE: pressure versions isn't implemented yet because this
657  /// element has never been tried with Fp preconditoner.
659  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
660  const unsigned& which_one=0);
661 
662  /// \short Number of scalars/fields output by this element. Reimplements
663  /// broken virtual function in base class.
664  unsigned nscalar_paraview() const
665  {
666  return 4;
667  }
668 
669  /// \short Write values of the i-th scalar field at the plot points. Needs
670  /// to be implemented for each new specific element type.
671  void scalar_value_paraview(std::ofstream& file_out,
672  const unsigned& i,
673  const unsigned& nplot) const
674  {
675  // Vector of local coordinates
676  Vector<double> s(2);
677 
678  // Loop over plot points
679  unsigned num_plot_points=nplot_points_paraview(nplot);
680  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
681  {
682 
683  // Get local coordinates of plot point
684  get_s_plot(iplot,nplot,s);
685 
686  // Velocities
687  if(i<3)
688  {
689  file_out << interpolated_u_axi_nst(s,i) << std::endl;
690  }
691  // Pressure
692  else if(i==3)
693  {
694  file_out << interpolated_p_axi_nst(s) << std::endl;
695  }
696  // Never get here
697  else
698  {
699  std::stringstream error_stream;
700  error_stream
701  << "Axisymmetric Navier-Stokes Elements only store 4 fields "
702  << std::endl;
703  throw OomphLibError(
704  error_stream.str(),
705  OOMPH_CURRENT_FUNCTION,
706  OOMPH_EXCEPTION_LOCATION);
707  }
708  }
709  }
710 
711  /// \short Name of the i-th scalar field. Default implementation
712  /// returns V1 for the first one, V2 for the second etc. Can (should!) be
713  /// overloaded with more meaningful names in specific elements.
714  std::string scalar_name_paraview(const unsigned& i) const
715  {
716  // Veloc
717  if(i<3)
718  {
719  return "Velocity "+StringConversion::to_string(i);
720  }
721  // Pressure field
722  else if(i==3)
723  {
724  return "Pressure";
725  }
726  // Never get here
727  else
728  {
729  std::stringstream error_stream;
730  error_stream
731  << "Axisymmetric Navier-Stokes Elements only store 4 fields "
732  << std::endl;
733  throw OomphLibError(
734  error_stream.str(),
735  OOMPH_CURRENT_FUNCTION,
736  OOMPH_EXCEPTION_LOCATION);
737  return " ";
738  }
739  }
740 
741  /// \short Output solution in data vector at local cordinates s:
742  /// r,z,u_r,u_z,u_phi,p
744  {
745  // Output the components of the position
746  for(unsigned i=0;i<2;i++)
747  {
748  data.push_back(interpolated_x(s,i));
749  }
750 
751  // Output the components of the FE representation of u at s
752  for(unsigned i=0;i<3;i++)
753  {
754  data.push_back(interpolated_u_axi_nst(s,i));
755  }
756 
757  // Output FE representation of p at s
758  data.push_back(interpolated_p_axi_nst(s));
759  }
760 
761 
762  /// \short Output function: x,y,[z],u,v,[w],p
763  /// in tecplot format. Default number of plot points
764  void output(std::ostream &outfile)
765  {
766  unsigned nplot=5;
767  output(outfile,nplot);
768  }
769 
770  /// \short Output function: x,y,[z],u,v,[w],p
771  /// in tecplot format. nplot points in each coordinate direction
772  void output(std::ostream &outfile, const unsigned &nplot);
773 
774 
775  /// \short Output function: x,y,[z],u,v,[w],p
776  /// in tecplot format. Default number of plot points
777  void output(FILE* file_pt)
778  {
779  unsigned nplot=5;
780  output(file_pt,nplot);
781  }
782 
783  /// \short Output function: x,y,[z],u,v,[w],p
784  /// in tecplot format. nplot points in each coordinate direction
785  void output(FILE* file_pt, const unsigned &nplot);
786 
787  /// \short Output function: x,y,[z],u,v,[w] in tecplot format.
788  /// nplot points in each coordinate direction at timestep t
789  /// (t=0: present; t>0: previous timestep)
790  void output_veloc(std::ostream &outfile,
791  const unsigned &nplot,
792  const unsigned& t);
793 
794  /// \short Output exact solution specified via function pointer
795  /// at a given number of plot points. Function prints as
796  /// many components as are returned in solution Vector
797  void output_fct(std::ostream &outfile,
798  const unsigned &nplot,
800 
801  /// \short Output exact solution specified via function pointer
802  /// at a given time and at a given number of plot points.
803  /// Function prints as many components as are returned in solution Vector.
804  void output_fct(std::ostream &outfile,
805  const unsigned &nplot,
806  const double& time,
808 
809  /// \short Validate against exact solution at given time
810  /// Solution is provided via function pointer.
811  /// Plot at a given number of plot points and compute L2 error
812  /// and L2 norm of velocity solution over element
813  void compute_error(std::ostream &outfile,
815  const double& time,
816  double& error, double& norm);
817 
818  /// \short Validate against exact solution.
819  /// Solution is provided via function pointer.
820  /// Plot at a given number of plot points and compute L2 error
821  /// and L2 norm of velocity solution over element
822  void compute_error(std::ostream &outfile,
824  double& error,
825  double& norm);
826 
827  /// Compute the element's residual Vector
829  {
830  //Call the generic residuals function with flag set to 0
831  //and using a dummy matrix argument
835  }
836 
837  ///\short Compute the element's residual Vector and the jacobian matrix
838  /// Virtual function can be overloaded by hanging-node version
840  DenseMatrix<double> &jacobian)
841  {
842  //Call the generic routine with the flag set to 1
844  residuals,jacobian,GeneralisedElement::Dummy_matrix,1);
845  }
846 
847  /// Add the element's contribution to its residuals vector,
848  /// jacobian matrix and mass matrix
850  Vector<double> &residuals,
851  DenseMatrix<double> &jacobian,
852  DenseMatrix<double> &mass_matrix)
853  {
854  //Call the generic routine with the flag set to 2
856  residuals,jacobian,mass_matrix,2);
857  }
858 
859  /// \short Compute derivatives of elemental residual vector with respect to
860  /// nodal coordinates. This function computes these terms analytically and
861  /// overwrites the default implementation in the FiniteElement base class.
862  /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
864  dresidual_dnodal_coordinates);
865 
866  /// Compute the element's residual Vector
868  double* const &parameter_pt,Vector<double> &dres_dparam)
869  {
870  //Call the generic residuals function with flag set to 0
871  //and using a dummy matrix argument
875  }
876 
877  /// \short Compute the element's residual Vector and the jacobian matrix
878  /// Virtual function can be overloaded by hanging-node version
880  double* const &parameter_pt,Vector<double> &dres_dparam,
881  DenseMatrix<double> &djac_dparam)
882  {
883  //Call the generic routine with the flag set to 1
885  parameter_pt,
886  dres_dparam,djac_dparam,GeneralisedElement::Dummy_matrix,1);
887  }
888 
889  /// Add the element's contribution to its residuals vector,
890  /// jacobian matrix and mass matrix
892  double* const &parameter_pt,
893  Vector<double> &dres_dparam,
894  DenseMatrix<double> &djac_dparam,
895  DenseMatrix<double> &dmass_matrix_dparam)
896  {
897  //Call the generic routine with the flag set to 2
899  parameter_pt,dres_dparam,djac_dparam,dmass_matrix_dparam,2);
900  }
901 
902 
903  /// Compute vector of FE interpolated velocity u at local coordinate s
905  Vector<double>& veloc) const
906  {
907  //Find number of nodes
908  unsigned n_node = nnode();
909  //Local shape function
910  Shape psi(n_node);
911  //Find values of shape function
912  shape(s,psi);
913 
914  for (unsigned i=0;i<3;i++)
915  {
916  //Index at which the nodal value is stored
917  unsigned u_nodal_index = u_index_axi_nst(i);
918  //Initialise value of u
919  veloc[i] = 0.0;
920  //Loop over the local nodes and sum
921  for(unsigned l=0;l<n_node;l++)
922  {
923  veloc[i] += nodal_value(l,u_nodal_index)*psi[l];
924  }
925  }
926  }
927 
928  /// Return FE interpolated velocity u[i] at local coordinate s
930  const unsigned &i) const
931  {
932  //Find number of nodes
933  unsigned n_node = nnode();
934  //Local shape function
935  Shape psi(n_node);
936  //Find values of shape function
937  shape(s,psi);
938 
939  //Get the index at which the velocity is stored
940  unsigned u_nodal_index = u_index_axi_nst(i);
941 
942  //Initialise value of u
943  double interpolated_u = 0.0;
944  //Loop over the local nodes and sum
945  for(unsigned l=0;l<n_node;l++)
946  {
947  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
948  }
949 
950  return(interpolated_u);
951  }
952 
953 
954  /// Return FE interpolated velocity u[i] at local coordinate s
955  // at time level t (t=0: present, t>0: history)
956  double interpolated_u_axi_nst(const unsigned &t,
957  const Vector<double> &s,
958  const unsigned &i) const
959  {
960  //Find number of nodes
961  unsigned n_node = nnode();
962  //Local shape function
963  Shape psi(n_node);
964  //Find values of shape function
965  shape(s,psi);
966 
967  //Get the index at which the velocity is stored
968  unsigned u_nodal_index = u_index_axi_nst(i);
969 
970  //Initialise value of u
971  double interpolated_u = 0.0;
972  //Loop over the local nodes and sum
973  for(unsigned l=0;l<n_node;l++)
974  {
975  interpolated_u += nodal_value(t,l,u_nodal_index)*psi[l];
976  }
977 
978  return(interpolated_u);
979  }
980 
981 
982  /// \short Compute the derivatives of the i-th component of
983  /// velocity at point s with respect
984  /// to all data that can affect its value. In addition, return the global
985  /// equation numbers corresponding to the data. The function is virtual
986  /// so that it can be overloaded in the refineable version
988  const Vector<double> &s,
989  const unsigned &i,
990  Vector<double> &du_ddata,
991  Vector<unsigned> &global_eqn_number)
992  {
993  //Find number of nodes
994  unsigned n_node = nnode();
995  //Local shape function
996  Shape psi(n_node);
997  //Find values of shape function
998  shape(s,psi);
999 
1000  //Find the index at which the velocity component is stored
1001  const unsigned u_nodal_index = u_index_axi_nst(i);
1002 
1003  //Find the number of dofs associated with interpolated u
1004  unsigned n_u_dof=0;
1005  for(unsigned l=0;l<n_node;l++)
1006  {
1007  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
1008  //If it's positive add to the count
1009  if(global_eqn >= 0) {++n_u_dof;}
1010  }
1011 
1012  //Now resize the storage schemes
1013  du_ddata.resize(n_u_dof,0.0);
1014  global_eqn_number.resize(n_u_dof,0);
1015 
1016  //Loop over th nodes again and set the derivatives
1017  unsigned count=0;
1018  //Loop over the local nodes and sum
1019  for(unsigned l=0;l<n_node;l++)
1020  {
1021  //Get the global equation number
1022  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
1023  if (global_eqn >= 0)
1024  {
1025  //Set the global equation number
1026  global_eqn_number[count] = global_eqn;
1027  //Set the derivative with respect to the unknown
1028  du_ddata[count] = psi[l];
1029  //Increase the counter
1030  ++count;
1031  }
1032  }
1033  }
1034 
1035 
1036  /// Return FE interpolated pressure at local coordinate s
1038  {
1039  //Find number of nodes
1040  unsigned n_pres = npres_axi_nst();
1041  //Local shape function
1042  Shape psi(n_pres);
1043  //Find values of shape function
1044  pshape_axi_nst(s,psi);
1045 
1046  //Initialise value of p
1047  double interpolated_p = 0.0;
1048  //Loop over the local nodes and sum
1049  for(unsigned l=0;l<n_pres;l++)
1050  {
1051  interpolated_p += p_axi_nst(l)*psi[l];
1052  }
1053 
1054  return(interpolated_p);
1055  }
1056 
1057  /// Return FE interpolated derivatives of velocity component u[i]
1058  /// w.r.t spatial local coordinate direction s[j] at local coordinate s
1060  const unsigned &i,
1061  const unsigned &j) const
1062  {
1063  // Determine number of nodes
1064  const unsigned n_node = nnode();
1065 
1066  // Allocate storage for local shape function and its derivatives
1067  // with respect to space
1068  Shape psif(n_node);
1069  DShape dpsifds(n_node,2);
1070 
1071  // Find values of shape function (ignore jacobian)
1072  (void)this->dshape_local(s,psif,dpsifds);
1073 
1074  // Get the index at which the velocity is stored
1075  const unsigned u_nodal_index = u_index_axi_nst(i);
1076 
1077  // Initialise value of duds
1078  double interpolated_duds = 0.0;
1079 
1080  // Loop over the local nodes and sum
1081  for(unsigned l=0;l<n_node;l++)
1082  {
1083  interpolated_duds += nodal_value(l,u_nodal_index)*dpsifds(l,j);
1084  }
1085 
1086  return(interpolated_duds);
1087  }
1088 
1089  /// Return FE interpolated derivatives of velocity component u[i]
1090  /// w.r.t spatial global coordinate direction x[j] at local coordinate s
1092  const unsigned &i,
1093  const unsigned &j) const
1094  {
1095  // Determine number of nodes
1096  const unsigned n_node = nnode();
1097 
1098  // Allocate storage for local shape function and its derivatives
1099  // with respect to space
1100  Shape psif(n_node);
1101  DShape dpsifdx(n_node,2);
1102 
1103  // Find values of shape function (ignore jacobian)
1104  (void)this->dshape_eulerian(s,psif,dpsifdx);
1105 
1106  // Get the index at which the velocity is stored
1107  const unsigned u_nodal_index = u_index_axi_nst(i);
1108 
1109  // Initialise value of dudx
1110  double interpolated_dudx = 0.0;
1111 
1112  // Loop over the local nodes and sum
1113  for(unsigned l=0;l<n_node;l++)
1114  {
1115  interpolated_dudx += nodal_value(l,u_nodal_index)*dpsifdx(l,j);
1116  }
1117 
1118  return(interpolated_dudx);
1119  }
1120 
1121  /// Return FE interpolated derivatives of velocity component u[i]
1122  /// w.r.t time at local coordinate s
1124  const unsigned &i) const
1125  {
1126  // Determine number of nodes
1127  const unsigned n_node = nnode();
1128 
1129  // Allocate storage for local shape function
1130  Shape psif(n_node);
1131 
1132  // Find values of shape function
1133  shape(s,psif);
1134 
1135  // Initialise value of dudt
1136  double interpolated_dudt = 0.0;
1137 
1138  // Loop over the local nodes and sum
1139  for(unsigned l=0;l<n_node;l++)
1140  {
1141  interpolated_dudt += du_dt_axi_nst(l,i)*psif[l];
1142  }
1143 
1144  return(interpolated_dudt);
1145  }
1146 
1147  /// \short Return FE interpolated derivatives w.r.t. nodal coordinates
1148  /// X_{pq} of the spatial derivatives of the velocity components
1149  /// du_i/dx_k at local coordinate s
1151  const unsigned &p,
1152  const unsigned &q,
1153  const unsigned &i,
1154  const unsigned &k) const
1155  {
1156  // Determine number of nodes
1157  const unsigned n_node = nnode();
1158 
1159  // Allocate storage for local shape function and its derivatives
1160  // with respect to space
1161  Shape psif(n_node);
1162  DShape dpsifds(n_node,2);
1163 
1164  // Find values of shape function (ignore jacobian)
1165  (void)this->dshape_local(s,psif,dpsifds);
1166 
1167  // Allocate memory for the jacobian and the inverse of the jacobian
1168  DenseMatrix<double> jacobian(2), inverse_jacobian(2);
1169 
1170  // Allocate memory for the derivative of the jacobian w.r.t. nodal coords
1171  DenseMatrix<double> djacobian_dX(2,n_node);
1172 
1173  // Allocate memory for the derivative w.r.t. nodal coords of the
1174  // spatial derivatives of the shape functions
1175  RankFourTensor<double> d_dpsifdx_dX(2,n_node,3,2);
1176 
1177  // Now calculate the inverse jacobian
1178  const double det =
1179  local_to_eulerian_mapping(dpsifds,jacobian,inverse_jacobian);
1180 
1181  // Calculate the derivative of the jacobian w.r.t. nodal coordinates
1182  // Note: must call this before "transform_derivatives(...)" since this
1183  // function requires dpsids rather than dpsidx
1184  dJ_eulerian_dnodal_coordinates(jacobian,dpsifds,djacobian_dX);
1185 
1186  // Calculate the derivative of dpsidx w.r.t. nodal coordinates
1187  // Note: this function also requires dpsids rather than dpsidx
1188  d_dshape_eulerian_dnodal_coordinates(det,jacobian,djacobian_dX,
1189  inverse_jacobian,dpsifds,d_dpsifdx_dX);
1190 
1191  // Get the index at which the velocity is stored
1192  const unsigned u_nodal_index = u_index_axi_nst(i);
1193 
1194  // Initialise value of dudx
1195  double interpolated_d_dudx_dX = 0.0;
1196 
1197  // Loop over the local nodes and sum
1198  for(unsigned l=0;l<n_node;l++)
1199  {
1200  interpolated_d_dudx_dX +=
1201  nodal_value(l,u_nodal_index)*d_dpsifdx_dX(p,q,l,k);
1202  }
1203 
1204  return(interpolated_d_dudx_dX);
1205  }
1206 
1207  /// Compute the volume of the element
1208  double compute_physical_size() const
1209  {
1210  // Initialise result
1211  double result = 0.0;
1212 
1213  // Figure out the global (Eulerian) spatial dimension of the
1214  // element by checking the Eulerian dimension of the nodes
1215  const unsigned n_dim_eulerian = nodal_dimension();
1216 
1217  // Allocate Vector of global Eulerian coordinates
1218  Vector<double> x(n_dim_eulerian);
1219 
1220  // Set the value of n_intpt
1221  const unsigned n_intpt = integral_pt()->nweight();
1222 
1223  // Vector of local coordinates
1224  const unsigned n_dim = dim();
1225  Vector<double> s(n_dim);
1226 
1227  // Loop over the integration points
1228  for(unsigned ipt=0;ipt<n_intpt;ipt++)
1229  {
1230  // Assign the values of s
1231  for(unsigned i=0;i<n_dim;i++) {s[i] = integral_pt()->knot(ipt,i);}
1232 
1233  // Assign the values of the global Eulerian coordinates
1234  for(unsigned i=0;i<n_dim_eulerian;i++) {x[i] = interpolated_x(s,i);}
1235 
1236  // Get the integral weight
1237  const double w = integral_pt()->weight(ipt);
1238 
1239  // Get Jacobian of mapping
1240  const double J = J_eulerian(s);
1241 
1242  // The integrand at the current integration point is r
1243  result += x[0]*w*J;
1244  }
1245 
1246  // Multiply by 2pi (integrating in azimuthal direction)
1247  return (2.0*MathematicalConstants::Pi*result);
1248  }
1249 
1250 };
1251 
1252 //////////////////////////////////////////////////////////////////////////////
1253 //////////////////////////////////////////////////////////////////////////////
1254 //////////////////////////////////////////////////////////////////////////////
1255 
1256 
1257 //==========================================================================
1258 ///Crouzeix_Raviart elements are Navier--Stokes elements with quadratic
1259 ///interpolation for velocities and positions, but a discontinuous linear
1260 ///pressure interpolation
1261 //==========================================================================
1262 class AxisymmetricQCrouzeixRaviartElement : public virtual QElement<2,3>,
1263  public virtual AxisymmetricNavierStokesEquations
1264 {
1265  private:
1266 
1267  /// Static array of ints to hold required number of variables at nodes
1268  static const unsigned Initial_Nvalue[];
1269 
1270  protected:
1271 
1272  /// Internal index that indicates at which internal data the pressure is
1273  /// stored
1275 
1276  /// \short Velocity shape and test functions and their derivs
1277  /// w.r.t. to global coords at local coordinate s (taken from geometry)
1278  ///Return Jacobian of mapping between local and global coordinates.
1279  inline double dshape_and_dtest_eulerian_axi_nst(const Vector<double> &s,
1280  Shape &psi,
1281  DShape &dpsidx, Shape &test,
1282  DShape &dtestdx) const;
1283 
1284  /// \short Velocity shape and test functions and their derivs
1285  /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
1286  ///Return Jacobian of mapping between local and global coordinates.
1287  inline double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt,
1288  Shape &psi,
1289  DShape &dpsidx,
1290  Shape &test,
1291  DShape &dtestdx) const;
1292 
1293  /// \short Shape/test functions and derivs w.r.t. to global coords at
1294  /// integration point ipt; return Jacobian of mapping (J). Also compute
1295  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
1297  const unsigned &ipt,
1298  Shape &psi,
1299  DShape &dpsidx,
1300  RankFourTensor<double> &d_dpsidx_dX,
1301  Shape &test,
1302  DShape &dtestdx,
1303  RankFourTensor<double> &d_dtestdx_dX,
1304  DenseMatrix<double> &djacobian_dX) const;
1305 
1306 
1307  /// Pressure shape functions at local coordinate s
1308  inline void pshape_axi_nst(const Vector<double> &s, Shape &psi) const;
1309 
1310  /// Pressure shape and test functions at local coordinte s
1311  inline void pshape_axi_nst(const Vector<double> &s, Shape &psi,
1312  Shape &test) const;
1313 
1314 
1315 public:
1316 
1317  /// Constructor, there are three internal values (for the pressure)
1320  {
1321  //Allocate and add one Internal data object that stores the three
1322  //pressure values
1324  }
1325 
1326  /// \short Number of values (pinned or dofs) required at local node n.
1327  virtual unsigned required_nvalue(const unsigned &n) const;
1328 
1329  /// \short Return the pressure values at internal dof i_internal
1330  /// (Discontinous pressure interpolation -- no need to cater for hanging
1331  /// nodes).
1332  double p_axi_nst(const unsigned &i) const
1334 
1335  /// Return number of pressure values
1336  unsigned npres_axi_nst() const {return 3;}
1337 
1338  ///Function to fix the internal pressure dof idof_internal
1339  void fix_pressure(const unsigned &p_dof, const double &pvalue)
1340  {
1343  }
1344 
1345  /// \short Compute traction at local coordinate s for outer unit normal N
1346  void get_traction(const Vector<double>& s, const Vector<double>& N,
1347  Vector<double>& traction) const;
1348 
1349  /// \short Overload the access function for the pressure's local
1350  /// equation numbers
1351  inline int p_local_eqn(const unsigned &n)const
1353 
1354  /// Redirect output to NavierStokesEquations output
1355  void output(std::ostream &outfile)
1357 
1358  /// Redirect output to NavierStokesEquations output
1359  void output(std::ostream &outfile, const unsigned &n_plot)
1361 
1362 
1363  /// Redirect output to NavierStokesEquations output
1364  void output(FILE* file_pt)
1366 
1367  /// Redirect output to NavierStokesEquations output
1368  void output(FILE* file_pt, const unsigned &n_plot)
1370 
1371  /// \short The number of "DOF types" that degrees of freedom in this element
1372  /// are sub-divided into: Velocity (3 components) and pressure.
1373  unsigned ndof_types() const
1374  {
1375  return 4;
1376  }
1377 
1378  /// \short Create a list of pairs for all unknowns in this element,
1379  /// so that the first entry in each pair contains the global equation
1380  /// number of the unknown, while the second one contains the number
1381  /// of the "DOF type" that this unknown is associated with.
1382  /// (Function can obviously only be called if the equation numbering
1383  /// scheme has been set up.)
1385  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const;
1386 
1387 
1388 };
1389 
1390 //Inline functions
1391 
1392 //=======================================================================
1393 /// Derivatives of the shape functions and test functions w.r.t. to global
1394 /// (Eulerian) coordinates. Return Jacobian of mapping between
1395 /// local and global coordinates.
1396 //=======================================================================
1399  DShape &dpsidx, Shape &test,
1400  DShape &dtestdx) const
1401 {
1402  //Call the geometrical shape functions and derivatives
1403  double J = this->dshape_eulerian(s,psi,dpsidx);
1404  //Loop over the test functions and derivatives and set them equal to the
1405  //shape functions
1406  for(unsigned i=0;i<9;i++)
1407  {
1408  test[i] = psi[i];
1409  dtestdx(i,0) = dpsidx(i,0);
1410  dtestdx(i,1) = dpsidx(i,1);
1411  }
1412  //Return the jacobian
1413  return J;
1414 }
1415 
1416 //=======================================================================
1417 /// Derivatives of the shape functions and test functions w.r.t. to global
1418 /// (Eulerian) coordinates. Return Jacobian of mapping between
1419 /// local and global coordinates.
1420 //=======================================================================
1423  DShape &dpsidx, Shape &test,
1424  DShape &dtestdx) const
1425 {
1426  //Call the geometrical shape functions and derivatives
1427  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
1428  //Loop over the test functions and derivatives and set them equal to the
1429  //shape functions
1430  for(unsigned i=0;i<9;i++)
1431  {
1432  test[i] = psi[i];
1433  dtestdx(i,0) = dpsidx(i,0);
1434  dtestdx(i,1) = dpsidx(i,1);
1435  }
1436  //Return the jacobian
1437  return J;
1438 }
1439 
1440 //=======================================================================
1441 /// Define the shape functions (psi) and test functions (test) and
1442 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
1443 /// and return Jacobian of mapping (J). Additionally compute the
1444 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
1445 ///
1446 /// Galerkin: Test functions = shape functions
1447 //=======================================================================
1450  const unsigned &ipt, Shape &psi, DShape &dpsidx,
1451  RankFourTensor<double> &d_dpsidx_dX,
1452  Shape &test, DShape &dtestdx,
1453  RankFourTensor<double> &d_dtestdx_dX,
1454  DenseMatrix<double> &djacobian_dX) const
1455  {
1456  // Call the geometrical shape functions and derivatives
1457  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
1458  djacobian_dX,d_dpsidx_dX);
1459 
1460  // Loop over the test functions and derivatives and set them equal to the
1461  // shape functions
1462  for(unsigned i=0;i<9;i++)
1463  {
1464  test[i] = psi[i];
1465 
1466  for(unsigned k=0;k<2;k++)
1467  {
1468  dtestdx(i,k) = dpsidx(i,k);
1469 
1470  for(unsigned p=0;p<2;p++)
1471  {
1472  for(unsigned q=0;q<9;q++)
1473  {
1474  d_dtestdx_dX(p,q,i,k) = d_dpsidx_dX(p,q,i,k);
1475  }
1476  }
1477  }
1478  }
1479 
1480  // Return the jacobian
1481  return J;
1482  }
1483 
1484 //=======================================================================
1485 /// Pressure shape functions
1486 //=======================================================================
1489  const
1490 {
1491  psi[0] = 1.0;
1492  psi[1] = s[0];
1493  psi[2] = s[1];
1494 }
1495 
1496 ///Define the pressure shape and test functions
1498 pshape_axi_nst(const Vector<double> &s, Shape &psi, Shape &test) const
1499 {
1500  //Call the pressure shape functions
1501  pshape_axi_nst(s,psi);
1502  //Loop over the test functions and set them equal to the shape functions
1503  for(unsigned i=0;i<3;i++) test[i] = psi[i];
1504 }
1505 
1506 
1507 //=======================================================================
1508 /// Face geometry of the Axisymmetric Crouzeix_Raviart elements
1509 //=======================================================================
1510 template<>
1512 public virtual QElement<1,3>
1513 {
1514  public:
1515  FaceGeometry() : QElement<1,3>() {}
1516 };
1517 
1518 //=======================================================================
1519 /// Face geometry of face geometry of the Axisymmetric Crouzeix_Raviart
1520 /// elements
1521 //=======================================================================
1522 template<>
1524 public virtual PointElement
1525 {
1526  public:
1528 };
1529 
1530 
1531 ////////////////////////////////////////////////////////////////////////////
1532 ////////////////////////////////////////////////////////////////////////////
1533 
1534 //=======================================================================
1535 /// Taylor--Hood elements are Navier--Stokes elements
1536 /// with quadratic interpolation for velocities and positions and
1537 /// continous linear pressure interpolation
1538 //=======================================================================
1540 public virtual QElement<2,3>,
1541 public virtual AxisymmetricNavierStokesEquations
1542 {
1543  private:
1544 
1545  /// Static array of ints to hold number of variables at node
1546  static const unsigned Initial_Nvalue[];
1547 
1548  protected:
1549 
1550  /// \short Static array of ints to hold conversion from pressure
1551  /// node numbers to actual node numbers
1552  static const unsigned Pconv[];
1553 
1554  /// \short Velocity shape and test functions and their derivs
1555  /// w.r.t. to global coords at local coordinate s (taken from geometry)
1556  /// Return Jacobian of mapping between local and global coordinates.
1557  inline double dshape_and_dtest_eulerian_axi_nst(const Vector<double> &s,
1558  Shape &psi,
1559  DShape &dpsidx, Shape &test,
1560  DShape &dtestdx) const;
1561 
1562  /// \short Velocity shape and test functions and their derivs
1563  /// w.r.t. to global coords at local coordinate s (taken from geometry)
1564  /// Return Jacobian of mapping between local and global coordinates.
1565  inline double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt,
1566  Shape &psi,
1567  DShape &dpsidx,
1568  Shape &test,
1569  DShape &dtestdx) const;
1570 
1571  /// \short Shape/test functions and derivs w.r.t. to global coords at
1572  /// integration point ipt; return Jacobian of mapping (J). Also compute
1573  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
1575  const unsigned &ipt,
1576  Shape &psi,
1577  DShape &dpsidx,
1578  RankFourTensor<double> &d_dpsidx_dX,
1579  Shape &test,
1580  DShape &dtestdx,
1581  RankFourTensor<double> &d_dtestdx_dX,
1582  DenseMatrix<double> &djacobian_dX) const;
1583 
1584  /// Pressure shape functions at local coordinate s
1585  inline void pshape_axi_nst(const Vector<double> &s, Shape &psi) const;
1586 
1587  /// Pressure shape and test functions at local coordinte s
1588  inline void pshape_axi_nst(const Vector<double> &s, Shape &psi,
1589  Shape &test) const;
1590 
1591  public:
1592 
1593  /// Constructor, no internal data points
1596 
1597  /// \short Number of values (pinned or dofs) required at node n. Can
1598  /// be overwritten for hanging node version
1599  inline virtual unsigned required_nvalue(const unsigned &n) const
1600  {return Initial_Nvalue[n];}
1601 
1602  /// \short Which nodal value represents the pressure?
1603  virtual int p_nodal_index_axi_nst() const {return 3;}
1604 
1605  /// \short Access function for the pressure values at local pressure
1606  /// node n_p (const version)
1607  double p_axi_nst(const unsigned &n_p) const
1608  {return nodal_value(Pconv[n_p],p_nodal_index_axi_nst());}
1609 
1610  /// Return number of pressure values
1611  unsigned npres_axi_nst() const {return 4;}
1612 
1613  /// Fix the pressure at local pressure node n_p
1614  void fix_pressure(const unsigned &n_p, const double &pvalue)
1615  {
1616  this->node_pt(Pconv[n_p])->pin(p_nodal_index_axi_nst());
1617  this->node_pt(Pconv[n_p])->set_value(p_nodal_index_axi_nst(),pvalue);
1618  }
1619 
1620  /// \short Compute traction at local coordinate s for outer unit normal N
1621  void get_traction(const Vector<double>& s, const Vector<double>& N,
1622  Vector<double>& traction) const;
1623 
1624  /// \short Overload the access function for the pressure's local
1625  /// equation numbers
1626  inline int p_local_eqn(const unsigned &n)const
1628 
1629  /// Redirect output to NavierStokesEquations output
1630  void output(std::ostream &outfile)
1632 
1633  /// Redirect output to NavierStokesEquations output
1634  void output(std::ostream &outfile, const unsigned &n_plot)
1636 
1637  /// Redirect output to NavierStokesEquations output
1638  void output(FILE* file_pt)
1640 
1641  /// Redirect output to NavierStokesEquations output
1642  void output(FILE* file_pt, const unsigned &n_plot)
1644 
1645  /// \short Returns the number of "DOF types" that degrees of freedom
1646  /// in this element are sub-divided into: Velocity (3 components) and pressure.
1647  unsigned ndof_types() const
1648  {
1649  return 4;
1650  }
1651 
1652  /// \short Create a list of pairs for all unknowns in this element,
1653  /// so that the first entry in each pair contains the global equation
1654  /// number of the unknown, while the second one contains the number
1655  /// of the "DOF type" that this unknown is associated with.
1656  /// (Function can obviously only be called if the equation numbering
1657  /// scheme has been set up.)
1659  std::list<std::pair<unsigned long, unsigned> >& dof_lookup_list) const;
1660 
1661 
1662 };
1663 
1664 //Inline functions
1665 
1666 //==========================================================================
1667 /// Derivatives of the shape functions and test functions w.r.t to
1668 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1669 /// local and global coordinates.
1670 //==========================================================================
1673  Shape &psi,
1674  DShape &dpsidx, Shape &test,
1675  DShape &dtestdx) const
1676 {
1677  //Call the geometrical shape functions and derivatives
1678  double J = this->dshape_eulerian(s,psi,dpsidx);
1679  //Loop over the test functions and derivatives and set them equal to the
1680  //shape functions
1681  for(unsigned i=0;i<9;i++)
1682  {
1683  test[i] = psi[i];
1684  dtestdx(i,0) = dpsidx(i,0);
1685  dtestdx(i,1) = dpsidx(i,1);
1686  }
1687  //Return the jacobian
1688  return J;
1689 }
1690 
1691 //==========================================================================
1692 /// Derivatives of the shape functions and test functions w.r.t to
1693 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1694 /// local and global coordinates.
1695 //==========================================================================
1698  Shape &psi, DShape &dpsidx,
1699  Shape &test,
1700  DShape &dtestdx) const
1701 {
1702  //Call the geometrical shape functions and derivatives
1703  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
1704  //Loop over the test functions and derivatives and set them equal to the
1705  //shape functions
1706  for(unsigned i=0;i<9;i++)
1707  {
1708  test[i] = psi[i];
1709  dtestdx(i,0) = dpsidx(i,0);
1710  dtestdx(i,1) = dpsidx(i,1);
1711  }
1712  //Return the jacobian
1713  return J;
1714 }
1715 
1716 //==========================================================================
1717 /// Define the shape functions (psi) and test functions (test) and
1718 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
1719 /// and return Jacobian of mapping (J). Additionally compute the
1720 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
1721 ///
1722 /// Galerkin: Test functions = shape functions
1723 //==========================================================================
1726  const unsigned &ipt, Shape &psi, DShape &dpsidx,
1727  RankFourTensor<double> &d_dpsidx_dX,
1728  Shape &test, DShape &dtestdx,
1729  RankFourTensor<double> &d_dtestdx_dX,
1730  DenseMatrix<double> &djacobian_dX) const
1731  {
1732  // Call the geometrical shape functions and derivatives
1733  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
1734  djacobian_dX,d_dpsidx_dX);
1735 
1736  // Loop over the test functions and derivatives and set them equal to the
1737  // shape functions
1738  for(unsigned i=0;i<9;i++)
1739  {
1740  test[i] = psi[i];
1741 
1742  for(unsigned k=0;k<2;k++)
1743  {
1744  dtestdx(i,k) = dpsidx(i,k);
1745 
1746  for(unsigned p=0;p<2;p++)
1747  {
1748  for(unsigned q=0;q<9;q++)
1749  {
1750  d_dtestdx_dX(p,q,i,k) = d_dpsidx_dX(p,q,i,k);
1751  }
1752  }
1753  }
1754  }
1755 
1756  // Return the jacobian
1757  return J;
1758  }
1759 
1760 //==========================================================================
1761 /// Pressure shape functions
1762 //==========================================================================
1765  const
1766 {
1767  //Local storage
1768  double psi1[2], psi2[2];
1769  //Call the OneDimensional Shape functions
1770  OneDimLagrange::shape<2>(s[0],psi1);
1771  OneDimLagrange::shape<2>(s[1],psi2);
1772 
1773  //Now let's loop over the nodal points in the element
1774  //s1 is the "x" coordinate, s2 the "y"
1775  for(unsigned i=0;i<2;i++)
1776  {
1777  for(unsigned j=0;j<2;j++)
1778  {
1779  /*Multiply the two 1D functions together to get the 2D function*/
1780  psi[2*i + j] = psi2[i]*psi1[j];
1781  }
1782  }
1783 }
1784 
1785 //==========================================================================
1786 /// Pressure shape and test functions
1787 //==========================================================================
1789 pshape_axi_nst(const Vector<double> &s, Shape &psi, Shape &test) const
1790 {
1791  //Call the pressure shape functions
1792  pshape_axi_nst(s,psi);
1793  //Loop over the test functions and set them equal to the shape functions
1794  for(unsigned i=0;i<4;i++) test[i] = psi[i];
1795 }
1796 
1797 //=======================================================================
1798 /// Face geometry of the Axisymmetric Taylor_Hood elements
1799 //=======================================================================
1800 template<>
1802 {
1803  public:
1804  FaceGeometry() : QElement<1,3>() {}
1805 };
1806 
1807 //=======================================================================
1808 /// Face geometry of the face geometry of the Axisymmetric Taylor_Hood elements
1809 //=======================================================================
1810 template<>
1812 public virtual PointElement
1813 {
1814  public:
1816 };
1817 
1818 
1819 //==========================================================
1820 /// Axisymmetric Taylor Hood upgraded to become projectable
1821 //==========================================================
1822 template<class TAYLOR_HOOD_ELEMENT>
1824 public virtual ProjectableElement<TAYLOR_HOOD_ELEMENT>
1825  {
1826 
1827  public:
1828 
1829  /// \short Specify the values associated with field fld.
1830  /// The information is returned in a vector of pairs which comprise
1831  /// the Data object and the value within it, that correspond to field fld.
1832  /// In the underlying Taylor Hood elements the fld-th velocities are stored
1833  /// at the fld-th value of the nodes; the pressures (the dim-th
1834  /// field) are the dim-th values at the vertex nodes etc.
1836  {
1837  // Create the vector
1838  Vector<std::pair<Data*,unsigned> > data_values;
1839 
1840  // Velocities dofs
1841  if (fld<3)
1842  {
1843  // Loop over all nodes
1844  unsigned nnod=this->nnode();
1845  for (unsigned j=0;j<nnod;j++)
1846  {
1847  // Add the data value associated with the velocity components
1848  data_values.push_back(std::make_pair(this->node_pt(j),fld));
1849  }
1850  }
1851  // Pressure
1852  else
1853  {
1854  // Loop over all vertex nodes
1855  unsigned Pconv_size=3;
1856  for (unsigned j=0;j<Pconv_size;j++)
1857  {
1858  // Add the data value associated with the pressure components
1859  unsigned vertex_index=this->Pconv[j];
1860  data_values.push_back(std::make_pair(this->node_pt(vertex_index),fld));
1861  }
1862  }
1863 
1864  // Return the vector
1865  return data_values;
1866 
1867  }
1868 
1869  /// \short Number of fields to be projected: dim+1, corresponding to
1870  /// velocity components and pressure
1872  {
1873  return 4;
1874  }
1875 
1876  /// \short Number of history values to be stored for fld-th field. Whatever
1877  /// the timestepper has set up for the velocity components and
1878  /// none for the pressure field (includes current value!)
1879  unsigned nhistory_values_for_projection(const unsigned &fld)
1880  {
1881  if (fld==3)
1882  {
1883  //pressure doesn't have history values
1884  return 1;
1885  }
1886  else
1887  {
1888  return this->node_pt(0)->ntstorage();
1889  }
1890  }
1891 
1892  ///\short Number of positional history values (includes current value!)
1894  {
1895  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
1896  }
1897 
1898  /// \short Return Jacobian of mapping and shape functions of field fld
1899  /// at local coordinate s
1900  double jacobian_and_shape_of_field(const unsigned &fld,
1901  const Vector<double> &s,
1902  Shape &psi)
1903  {
1904  unsigned n_dim=this->dim();
1905  unsigned n_node=this->nnode();
1906 
1907  if (fld==3)
1908  {
1909  //We are dealing with the pressure
1910  this->pshape_axi_nst(s,psi);
1911 
1912  Shape psif(n_node),testf(n_node);
1913  DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
1914 
1915  //Domain Shape
1916  double J=this->dshape_and_dtest_eulerian_axi_nst(s,psif,dpsifdx,
1917  testf,dtestfdx);
1918  return J;
1919  }
1920  else
1921  {
1922  Shape testf(n_node);
1923  DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
1924 
1925  //Domain Shape
1926  double J=this->dshape_and_dtest_eulerian_axi_nst(s,psi,dpsifdx,
1927  testf,dtestfdx);
1928  return J;
1929  }
1930  }
1931 
1932 
1933 
1934  /// \short Return interpolated field fld at local coordinate s, at time level
1935  /// t (t=0: present; t>0: history values)
1936  double get_field(const unsigned &t,
1937  const unsigned &fld,
1938  const Vector<double>& s)
1939  {
1940  unsigned n_node=this->nnode();
1941 
1942  //If fld=3, we deal with the pressure
1943  if (fld==3)
1944  {
1945  return this->interpolated_p_axi_nst(s);
1946  }
1947  // Velocity
1948  else
1949  {
1950  //Find the index at which the variable is stored
1951  unsigned u_nodal_index = this->u_index_axi_nst(fld);
1952 
1953  //Local shape function
1954  Shape psi(n_node);
1955 
1956  //Find values of shape function
1957  this->shape(s,psi);
1958 
1959  //Initialise value of u
1960  double interpolated_u = 0.0;
1961 
1962  //Sum over the local nodes at that time
1963  for(unsigned l=0;l<n_node;l++)
1964  {
1965  interpolated_u += this->nodal_value(t,l,u_nodal_index)*psi[l];
1966  }
1967  return interpolated_u;
1968  }
1969  }
1970 
1971 
1972 
1973  ///Return number of values in field fld
1974  unsigned nvalue_of_field(const unsigned &fld)
1975  {
1976  if (fld==3)
1977  {
1978  return this->npres_axi_nst();
1979  }
1980  else
1981  {
1982  return this->nnode();
1983  }
1984  }
1985 
1986 
1987  ///Return local equation number of value j in field fld.
1988  int local_equation(const unsigned &fld,
1989  const unsigned &j)
1990  {
1991  if (fld==3)
1992  {
1993  return this->p_local_eqn(j);
1994  }
1995  else
1996  {
1997  const unsigned u_nodal_index = this->u_index_axi_nst(fld);
1998  return this->nodal_local_eqn(j,u_nodal_index);
1999  }
2000  }
2001 
2002  };
2003 
2004 
2005 //=======================================================================
2006 /// Face geometry for element is the same as that for the underlying
2007 /// wrapped element
2008 //=======================================================================
2009  template<class ELEMENT>
2011  : public virtual FaceGeometry<ELEMENT>
2012  {
2013  public:
2014  FaceGeometry() : FaceGeometry<ELEMENT>() {}
2015  };
2016 
2017 
2018 //=======================================================================
2019 /// Face geometry of the Face Geometry for element is the same as
2020 /// that for the underlying wrapped element
2021 //=======================================================================
2022  template<class ELEMENT>
2024  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
2025  {
2026  public:
2028  };
2029 
2030 
2031 //==========================================================
2032 /// Crouzeix Raviart upgraded to become projectable
2033 //==========================================================
2034  template<class CROUZEIX_RAVIART_ELEMENT>
2036  public virtual ProjectableElement<CROUZEIX_RAVIART_ELEMENT>
2037  {
2038 
2039  public:
2040 
2041  /// \short Specify the values associated with field fld.
2042  /// The information is returned in a vector of pairs which comprise
2043  /// the Data object and the value within it, that correspond to field fld.
2044  /// In the underlying Crouzeix Raviart elements the
2045  /// fld-th velocities are stored
2046  /// at the fld-th value of the nodes; the pressures are stored internally
2048  {
2049  // Create the vector
2050  Vector<std::pair<Data*,unsigned> > data_values;
2051 
2052  // Velocities dofs
2053  if (fld < 3)
2054  {
2055  // Loop over all nodes
2056  const unsigned n_node=this->nnode();
2057  for (unsigned n=0;n<n_node;n++)
2058  {
2059  // Add the data value associated with the velocity components
2060  data_values.push_back(std::make_pair(this->node_pt(n),fld));
2061  }
2062  }
2063  // Pressure
2064  else
2065  {
2066  //Need to push back the internal data
2067  const unsigned n_press = this->npres_axi_nst();
2068  //Loop over all pressure values
2069  for(unsigned j=0;j<n_press;j++)
2070  {
2071  data_values.push_back(
2072  std::make_pair(
2073  this->internal_data_pt(this->P_axi_nst_internal_index),j));
2074  }
2075  }
2076 
2077  // Return the vector
2078  return data_values;
2079  }
2080 
2081  /// \short Number of fields to be projected: dim+1, corresponding to
2082  /// velocity components and pressure
2084  {
2085  return 4;
2086  }
2087 
2088  /// \short Number of history values to be stored for fld-th field. Whatever
2089  /// the timestepper has set up for the velocity components and
2090  /// none for the pressure field (includes current value!)
2091  unsigned nhistory_values_for_projection(const unsigned &fld)
2092  {
2093  if (fld==3)
2094  {
2095  //pressure doesn't have history values
2096  return 1;
2097  }
2098  else
2099  {
2100  return this->node_pt(0)->ntstorage();
2101  }
2102  }
2103 
2104  ///\short Number of positional history values (includes current value!)
2106  {
2107  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
2108  }
2109 
2110  /// \short Return Jacobian of mapping and shape functions of field fld
2111  /// at local coordinate s
2112  double jacobian_and_shape_of_field(const unsigned &fld,
2113  const Vector<double> &s,
2114  Shape &psi)
2115  {
2116  unsigned n_dim=this->dim();
2117  unsigned n_node=this->nnode();
2118 
2119  if (fld==3)
2120  {
2121  //We are dealing with the pressure
2122  this->pshape_axi_nst(s,psi);
2123 
2124  Shape psif(n_node),testf(n_node);
2125  DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
2126 
2127  //Domain Shape
2128  double J=this->dshape_and_dtest_eulerian_axi_nst(s,psif,dpsifdx,
2129  testf,dtestfdx);
2130  return J;
2131  }
2132  else
2133  {
2134  Shape testf(n_node);
2135  DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
2136 
2137  //Domain Shape
2138  double J=this->dshape_and_dtest_eulerian_axi_nst(s,psi,dpsifdx,
2139  testf,dtestfdx);
2140  return J;
2141  }
2142  }
2143 
2144 
2145 
2146  /// \short Return interpolated field fld at local coordinate s, at time level
2147  /// t (t=0: present; t>0: history values)
2148  double get_field(const unsigned &t,
2149  const unsigned &fld,
2150  const Vector<double>& s)
2151  {
2152  //unsigned n_dim =this->dim();
2153  //unsigned n_node=this->nnode();
2154 
2155  //If fld=n_dim, we deal with the pressure
2156  if (fld==3)
2157  {
2158  return this->interpolated_p_axi_nst(s);
2159  }
2160  // Velocity
2161  else
2162  {
2163  return this->interpolated_u_axi_nst(t,s,fld);
2164  }
2165  }
2166 
2167 
2168  ///Return number of values in field fld
2169  unsigned nvalue_of_field(const unsigned &fld)
2170  {
2171  if (fld==3)
2172  {
2173  return this->npres_axi_nst();
2174  }
2175  else
2176  {
2177  return this->nnode();
2178  }
2179  }
2180 
2181 
2182  ///Return local equation number of value j in field fld.
2183  int local_equation(const unsigned &fld,
2184  const unsigned &j)
2185  {
2186  if (fld==3)
2187  {
2188  return this->p_local_eqn(j);
2189  }
2190  else
2191  {
2192  const unsigned u_nodal_index = this->u_index_axi_nst(fld);
2193  return this->nodal_local_eqn(j,u_nodal_index);
2194  }
2195  }
2196 
2197  };
2198 
2199 
2200 //=======================================================================
2201 /// Face geometry for element is the same as that for the underlying
2202 /// wrapped element
2203 //=======================================================================
2204  template<class ELEMENT>
2206  : public virtual FaceGeometry<ELEMENT>
2207  {
2208  public:
2209  FaceGeometry() : FaceGeometry<ELEMENT>() {}
2210  };
2211 
2212 
2213 //=======================================================================
2214 /// Face geometry of the Face Geometry for element is the same as
2215 /// that for the underlying wrapped element
2216 //=======================================================================
2217  template<class ELEMENT>
2220  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
2221  {
2222  public:
2224  };
2225 
2226 
2227 
2228 
2229 }
2230 
2231 #endif
virtual void pshape_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
double *& re_invro_pt()
Pointer to global inverse Froude number.
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
virtual unsigned u_index_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data...
Definition: elements.h:268
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
double *& re_invfr_pt()
Pointer to global inverse Froude number.
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
Compute traction at local coordinate s for outer unit normal N.
double interpolated_duds_axi_nst(const Vector< double > &s, const unsigned &i, const unsigned &j) const
virtual double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and derivatives w.r.t. global coords at ipt-th integration point Return J...
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
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
void(*&)(const double &time, const Vector< double > &x, Vector< double > &f) axi_nst_body_force_fct_pt()
Access function for the body-force pointer.
virtual int p_nodal_index_axi_nst() const
Which nodal value represents the pressure?
virtual void dinterpolated_u_axi_nst_ddata(const Vector< double > &s, const unsigned &i, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
Compute the derivatives of the i-th component of velocity at point s with respect to all data that ca...
unsigned nfields_for_projection()
Number of fields to be projected: dim+1, corresponding to velocity components and pressure...
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. Whatever the timestepper has set up for the v...
double get_source_fct(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
double interpolated_u_axi_nst(const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s.
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
const double & re_invfr() const
Global inverse Froude number.
virtual void fill_in_generic_residual_contribution_axi_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=2 or 1 or 0: compute the Jacobian and/or ...
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
virtual void fill_in_generic_dresidual_contribution_axi_nst(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, unsigned flag)
Compute the derivative of residuals for the Navier–Stokes equations; with respect to a parameeter fla...
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:847
virtual void get_viscosity_ratio_axisym_nst(const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, double &visc_ratio)
Calculate the viscosity ratio relative to the viscosity used in the definition of the Reynolds number...
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
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
void fill_in_contribution_to_djacobian_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Compute the element's residual Vector and the jacobian matrix Virtual function can be overloaded by h...
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
cstr elem_len * i
Definition: cfortran.h:607
unsigned npres_axi_nst() const
Return number of pressure values.
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
const double Pi
50 digits from maple
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at ipt-th integation point...
double kin_energy() const
Get integral of kinetic energy over element.
Crouzeix Raviart upgraded to become projectable.
void interpolated_u_axi_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
virtual void get_source_fct_gradient(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution specified via function pointer at a given number of plot points. Function prints as many components as are returned in solution Vector.
double interpolated_p_axi_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
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
A general Finite Element class.
Definition: elements.h:1271
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s...
Definition: elements.cc:3981
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (includes current value!)
virtual double dshape_and_dtest_eulerian_axi_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and derivatives w.r.t. global coords at local coordinate s...
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed...
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
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
double interpolated_d_dudx_dX_axi_nst(const Vector< double > &s, const unsigned &p, const unsigned &q, const unsigned &i, const unsigned &k) const
Return FE interpolated derivatives w.r.t. nodal coordinates X_{pq} of the spatial derivatives of the ...
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers...
void output_veloc(std::ostream &outfile, const unsigned &nplot, const unsigned &t)
Output function: x,y,[z],u,v,[w] in tecplot format. nplot points in each coordinate direction at time...
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
void output(std::ostream &outfile, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2358
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element's residual Vector.
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:375
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
unsigned u_index_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored with a common interface for u...
void fill_in_contribution_to_dresiduals_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam)
Compute the element's residual Vector.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocity (3 c...
void traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
Compute traction (on the viscous scale) at local coordinate s for outer unit normal N...
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2700
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class...
int p_local_eqn(const unsigned &n) const
Overload the access function for the pressure's local equation numbers.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
A Rank 4 Tensor class.
Definition: matrices.h:1625
virtual void dJ_eulerian_dnodal_coordinates(const DenseMatrix< double > &jacobian, const DShape &dpsids, DenseMatrix< double > &djacobian_dX) const
A template-free interface that calculates the derivative of the jacobian of a mapping with respect to...
Definition: elements.cc:2609
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
void(* Body_force_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to body force function.
virtual double local_to_eulerian_mapping(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates, given the derivatives of the shape function...
Definition: elements.h:1461
double * ReInvRo_pt
Pointer to global Reynolds number x inverse Rossby number (used when in a rotating frame) ...
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. Whatever the timestepper has set up for the v...
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (includes current value!)
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
AxisymmetricQCrouzeixRaviartElement()
Constructor, there are three internal values (for the pressure)
virtual void get_body_force_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Calculate the body force fct at a given time and Eulerian position.
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
virtual void get_body_force_gradient_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
Get gradient of body force term at (Eulerian) position x. Computed via function pointer (if set) or b...
void fix_pressure(const unsigned &n_p, const double &pvalue)
Fix the pressure at local pressure node n_p.
virtual int p_local_eqn(const unsigned &n) const =0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
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
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at node n. Can be overwritten for hanging node version...
virtual unsigned npres_axi_nst() const =0
Function to return number of pressure degrees of freedom.
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative. Note: By default, ALE is enabled, at the expense of possibly creating unnecessary work in problems where the mesh is, in fact, stationary.
void pshape_axi_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
A Rank 3 Tensor class.
Definition: matrices.h:1337
double du_dt_axi_nst(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes...
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Validate against exact solution at given time Solution is provided via function pointer. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element.
double pressure_integral() const
Integral of pressure over element.
static char t char * s
Definition: cfortran.h:572
void output(std::ostream &outfile, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
double * Re_pt
Pointer to global Reynolds number.
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
unsigned nfields_for_projection()
Number of fields to be projected: dim+1, corresponding to velocity components and pressure...
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. This function computes these terms analytically and overwrites the default implementation in the FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
double(*&)(const double &time, const Vector< double > &x) source_fct_pt()
Access function for the source-function pointer.
unsigned ndof_types() const
Returns the number of "DOF types" that degrees of freedom in this element are sub-divided into: Veloc...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double p_axi_nst(const unsigned &i) const
Return the pressure values at internal dof i_internal (Discontinous pressure interpolation – no need ...
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
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
const double & re() const
Reynolds number.
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
double *& density_ratio_pt()
Pointer to Density ratio.
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
Compute traction at local coordinate s for outer unit normal N.
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
double *& re_pt()
Pointer to Reynolds number.
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
unsigned n_u_nst() const
Return the number of velocity components for use in general FluidInterface clas.
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number) ...
double p_axi_nst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version) ...
void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed...
Axisymmetric Taylor Hood upgraded to become projectable.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
void fix_pressure(const unsigned &p_dof, const double &pvalue)
Function to fix the internal pressure dof idof_internal.
double interpolated_dudx_axi_nst(const Vector< double > &s, const unsigned &i, const unsigned &j) const
void pshape_axi_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
int p_local_eqn(const unsigned &n) const
Overload the access function for the pressure's local.
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
const Vector< double > & g() const
Vector of gravitational components.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
double dissipation() const
Return integral of dissipation over element.
double interpolated_dudt_axi_nst(const Vector< double > &s, const unsigned &i) const
void point_output_data(const Vector< double > &s, Vector< double > &data)
Output solution in data vector at local cordinates s: r,z,u_r,u_z,u_phi,p.
double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
Definition: elements.cc:66
Vector< double > * G_pt
Pointer to global gravity Vector.
virtual int p_nodal_index_axi_nst() const
Which nodal value represents the pressure?
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
double compute_physical_size() const
Compute the volume of the element.
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:291
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
void fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
virtual double p_axi_nst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes...
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
unsigned npres_axi_nst() const
Return number of pressure values.
void fill_in_contribution_to_hessian_vector_products(Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Compute the hessian tensor vector products required to perform continuation of bifurcations analytica...
double dshape_and_dtest_eulerian_axi_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual Vector and the jacobian matrix Virtual function can be overloaded by h...
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...
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number) ...
virtual void d_dshape_eulerian_dnodal_coordinates(const double &det_jacobian, const DenseMatrix< double > &jacobian, const DenseMatrix< double > &djacobian_dX, const DenseMatrix< double > &inverse_jacobian, const DShape &dpsids, RankFourTensor< double > &d_dpsidx_dX) const
A template-free interface that calculates the derivative w.r.t. the nodal coordinates of the derivat...
Definition: elements.cc:2695
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.
double(* Source_fct_pt)(const double &time, const Vector< double > &x)
Pointer to volumetric source function.
AxisymmetricQTaylorHoodElement()
Constructor, no internal data points.
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:228
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
Definition: elements.h:1914
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
double dshape_and_dtest_eulerian_axi_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
double interpolated_u_axi_nst(const unsigned &t, const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s.
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types)...
void output(FILE *file_pt)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
AxisymmetricNavierStokesEquations()
Constructor: NULL the body force and source function.
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:596
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3252