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