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