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