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-2022 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
44namespace 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 //======================================================================
117 : public virtual FiniteElement,
119 {
120 private:
121 /// Static "magic" number that indicates that the pressure is
122 /// not stored at a node
124
125 /// Static default value for the physical constants (all initialised to
126 /// zero)
128
129 /// Static default value for the physical ratios (all are initialised to
130 /// one)
132
133 /// Static default value for the gravity vector
135
136 protected:
137 /// Static boolean telling us whether we pre-multiply the pressure and
138 /// continuity by the 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)
146
147 /// Pointer to the density ratio (relative to the density used in the
148 /// definition of the Reynolds number)
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
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
180
181 // Boolean that indicates whether we use the extrapolated strain rate
182 // for calculation of viscosity or not
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
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.
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.
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.
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)
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
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
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.
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
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
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
395 Source_fct_pt(0),
398 ALE_is_disabled(false)
399 {
400 // Set all the Physical parameter pointers to the default value zero
406 // Set the Physical ratios to the default value of 1
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
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
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
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
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,
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
521 {
522 return Constitutive_eqn_pt;
523 }
524
525 /// Switch to fully implict evaluation of non-Newtonian const eqn
527 {
529 }
530
531 /// Use extrapolation for non-Newtonian const eqn
533 {
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
574
575 // Initialise dudt
576 double dudt = 0.0;
577 // Loop over the timesteps, if there is a non Steady timestepper
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!
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.
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 {
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,
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,
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(
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
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);
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.
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
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,
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,
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,
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,
850 double& error,
851 double& norm);
852
853 /// Compute the element's residual Vector
855 {
856 // Call the generic residuals function with flag set to 0
857 // and using a dummy matrix argument
859 residuals,
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
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
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
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
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}
897 RankThreeTensor<double>& dresidual_dnodal_coordinates);
898
899 /// Compute the element's residual Vector
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
906 parameter_pt,
907 dres_dparam,
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
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
922 parameter_pt,
923 dres_dparam,
924 djac_dparam,
926 1);
927 }
928
929 /// Add the element's contribution to its residuals vector,
930 /// jacobian matrix and mass matrix
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
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
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
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
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
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
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
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
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
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
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
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 //==========================================================================
1311 : public virtual QElement<2, 3>,
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
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.
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.
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.
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)
1368 : QElement<2, 3>(),
1370 {
1371 // Allocate and add one Internal data object that stores the three
1372 // pressure values
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 {
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);
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 {
1410 }
1411
1412 /// Redirect output to NavierStokesEquations output
1413 void output(std::ostream& outfile)
1414 {
1416 }
1417
1418 /// Redirect output to NavierStokesEquations output
1419 void output(std::ostream& outfile, const unsigned& n_plot)
1420 {
1422 n_plot);
1423 }
1424
1425
1426 /// Redirect output to NavierStokesEquations output
1427 void output(FILE* file_pt)
1428 {
1430 }
1431
1432 /// Redirect output to NavierStokesEquations output
1433 void output(FILE* file_pt, const unsigned& n_plot)
1434 {
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
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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 //=======================================================================
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
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<>
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<>
1598 : public virtual PointElement
1599 {
1600 public:
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 //=======================================================================
1614 : public virtual QElement<2, 3>,
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.
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.
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.
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
1669 : QElement<2, 3>(),
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 {
1717 }
1718
1719 /// Redirect output to NavierStokesEquations output
1720 void output(std::ostream& outfile)
1721 {
1723 }
1724
1725 /// Redirect output to NavierStokesEquations output
1726 void output(std::ostream& outfile, const unsigned& n_plot)
1727 {
1729 n_plot);
1730 }
1731
1732 /// Redirect output to NavierStokesEquations output
1733 void output(FILE* file_pt)
1734 {
1736 }
1737
1738 /// Redirect output to NavierStokesEquations output
1739 void output(FILE* file_pt, const unsigned& n_plot)
1740 {
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
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 //==========================================================================
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 //==========================================================================
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 //==========================================================================
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 //==========================================================================
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 //==========================================================================
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<>
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<>
1917 : public virtual PointElement
1918 {
1919 public:
1921 };
1922
1923
1924 //==========================================================
1925 /// GeneralisedNewtonianAxisymmetric Taylor Hood upgraded to become
1926 /// projectable
1927 //==========================================================
1928 template<class TAYLOR_HOOD_ELEMENT>
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.
1940 {
1941 // Create the vector
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
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!)
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>
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>
2126 : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
2127 {
2128 public:
2130 };
2131
2132
2133 //==========================================================
2134 /// Crouzeix Raviart upgraded to become projectable
2135 //==========================================================
2136 template<class CROUZEIX_RAVIART_ELEMENT>
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
2148 {
2149 // Create the vector
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
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!)
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>
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>
2317 : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
2318 {
2319 public:
2321 };
2322
2323
2324} // namespace oomph
2325
2326#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
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
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
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
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
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
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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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
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
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
virtual 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
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
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
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
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
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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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
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
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
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
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
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1198
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
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
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
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
A class for elements that solve the unsteady axisymmetric Navier–Stokes equations in cylindrical pola...
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...
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...
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...
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,...
static bool Pre_multiply_by_viscosity_ratio
Static boolean telling us whether we pre-multiply the pressure and continuity by the viscosity ratio.
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...
double(*&)(const double &time, const Vector< double > &x) source_fct_pt()
Access function for the source-function pointer.
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...
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,...
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...
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.
virtual void fill_in_generic_dresidual_contribution_axi_nst(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, unsigned flag)
Compute the derivative of residuals for the Navier–Stokes equations; with respect to a parameeter fla...
double interpolated_u_axi_nst(const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s.
void fill_in_contribution_to_hessian_vector_products(Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Compute the hessian tensor vector products required to perform continuation of bifurcations analytica...
double 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 ...
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
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.
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....
unsigned n_u_nst() const
Return the number of velocity components for use in general FluidInterface class.
double * ReInvRo_pt
Pointer to global Reynolds number x inverse Rossby number (used when in a rotating frame)
void output(FILE *file_pt)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
virtual void extrapolated_strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Extrapolated strain-rate tensor: where (in that order) based on the previous time steps evaluated a...
virtual void pshape_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
virtual unsigned npres_axi_nst() const =0
Function to return number of pressure degrees of freedom.
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative....
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.
void use_current_strainrate_to_compute_second_invariant()
Switch to fully implict evaluation of non-Newtonian const eqn.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one,...
GeneralisedNewtonianConstitutiveEquation< 3 > *& constitutive_eqn_pt()
Access function for the constitutive equation pointer.
GeneralisedNewtonianConstitutiveEquation< 3 > * Constitutive_eqn_pt
Pointer to the generalised Newtonian constitutive equation.
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number)
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
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...
GeneralisedNewtonianAxisymmetricNavierStokesEquations()
Constructor: NULL the body force and source function.
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
void(*&)(const double &time, const Vector< double > &x, Vector< double > &f) axi_nst_body_force_fct_pt()
Access function for the body-force pointer.
virtual 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 ...
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....
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
void fill_in_contribution_to_dresiduals_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam)
Compute the element's residual Vector.
double du_dt_axi_nst(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes.
void point_output_data(const Vector< double > &s, Vector< double > &data)
Output solution in data vector at local cordinates s: r,z,u_r,u_z,u_phi,p.
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
double get_source_fct(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
void traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
Compute traction (on the viscous scale) at local coordinate s for outer unit normal N.
virtual double dshape_and_dtest_eulerian_axi_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and derivatives w.r.t. global coords at local coordinate s....
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element's residual Vector.
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...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual Vector and the jacobian matrix Virtual function can be overloaded by h...
virtual void get_body_force_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Calculate the body force fct at a given time and Eulerian position.
void(* Body_force_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to body force function.
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...
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
virtual int p_nodal_index_axi_nst() const
Which nodal value represents the pressure?
unsigned nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class.
double(* Source_fct_pt)(const double &time, const Vector< double > &x)
Pointer to volumetric source function.
const Vector< double > & g() const
Vector of gravitational components.
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
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.
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
double interpolated_p_axi_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
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...
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.
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.
virtual void fill_in_generic_residual_contribution_axi_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=2 or 1 or 0: compute the Jacobian and/or ...
void 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...
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.
void use_extrapolated_strainrate_to_compute_second_invariant()
Use extrapolation for non-Newtonian const eqn.
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
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...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed....
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
virtual void get_body_force_gradient_axi_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
Get gradient of body force term at (Eulerian) position x. Computed via function pointer (if set) or b...
void interpolated_u_axi_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void pshape_axi_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
double p_axi_nst(const unsigned &i) const
Return the pressure values at internal dof i_internal (Discontinous pressure interpolation – no need ...
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
Compute traction at local coordinate s for outer unit normal N.
double 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...
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
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...
int p_local_eqn(const unsigned &n) const
Overload the access function for the pressure's local equation numbers.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void output(std::ostream &outfile, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement()
Constructor, there are three internal values (for the pressure)
void fix_pressure(const unsigned &p_dof, const double &pvalue)
Function to fix the internal pressure dof idof_internal.
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocity and ...
unsigned P_axi_nst_internal_index
Internal index that indicates at which internal data the pressure is stored.
/////////////////////////////////////////////////////////////////////////
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.
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
unsigned ndof_types() const
Returns the number of "DOF types" that degrees of freedom in this element are sub-divided into: Veloc...
virtual int p_nodal_index_axi_nst() const
Which nodal value represents the pressure?
void output(std::ostream &outfile, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
void pshape_axi_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
int p_local_eqn(const unsigned &n) const
Overload the access function for the pressure's local equation numbers.
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...
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
double p_axi_nst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version)
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction) const
Compute traction at local coordinate s for outer unit normal N.
double dshape_and_dtest_eulerian_axi_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
void fix_pressure(const unsigned &n_p, const double &pvalue)
Fix the pressure at local pressure node n_p.
A Base class defining the generalise Newtonian constitutive relation.
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...
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (includes current value!)
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...
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 ...
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
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.
unsigned nfields_for_projection()
Number of fields to be projected: dim+1, corresponding to velocity components and pressure.
GeneralisedNewtonianAxisymmetric Taylor Hood upgraded to become projectable.
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.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
unsigned nfields_for_projection()
Number of fields to be projected: dim+1, corresponding to velocity components and pressure.
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. Whatever the timestepper has set up for the v...
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 ...
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (includes current value!)
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...
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:192
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:5231
A GeneralisedNewtonianConstitutiveEquation class defining a Newtonian fluid.
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: elements.h:3439
Wrapper class for projectable elements. Adds "projectability" to the underlying ELEMENT.
Definition: projection.h:183
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1701
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
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
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
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
const double Pi
50 digits from maple
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:608
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).
//////////////////////////////////////////////////////////////////// ////////////////////////////////...