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