linearised_axisym_navier_stokes_elements.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Header file for linearised axisymmetric Navier-Stokes elements
27 
28 #ifndef OOMPH_LINEARISED_AXISYM_NAVIER_STOKES_ELEMENTS_HEADER
29 #define OOMPH_LINEARISED_AXISYM_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 includes
37 #include "../generic/Qelements.h"
38 #include "../generic/fsi.h"
39 
40 namespace oomph
41 {
42  //=======================================================================
43  /// A class for elements that solve the linearised version of the
44  /// unsteady Navier--Stokes equations in cylindrical polar coordinates,
45  /// where we have Fourier-decomposed in the azimuthal direction so that
46  /// the theta-dependance is replaced by an azimuthal mode number.
47  //=======================================================================
49  : public virtual FiniteElement
50  {
51  private:
52  /// Static "magic" number that indicates that the pressure is not
53  /// stored at a node
55 
56  /// Static default value for the physical constants
57  /// (all initialised to zero)
59 
60  /// Static default value for the azimuthal mode number (zero)
62 
63  /// Static default value for the physical ratios (all initialised to one)
65 
66  protected:
67  // Physical constants
68  // ------------------
69 
70  /// Pointer to the viscosity ratio (relative to the
71  /// viscosity used in the definition of the Reynolds number)
73 
74  /// Pointer to the density ratio (relative to the
75  /// density used in the definition of the Reynolds number)
77 
78  /// Pointer to global Reynolds number
79  double* Re_pt;
80 
81  /// Pointer to global Reynolds number x Strouhal number (=Womersley)
82  double* ReSt_pt;
83 
84  /// Pointer to azimuthal mode number k in e^ik(theta) decomposition
86 
87  /// Pointer to base flow solution (velocity components) function
88  void (*Base_flow_u_fct_pt)(const double& time,
89  const Vector<double>& x,
90  Vector<double>& result);
91 
92  /// Pointer to derivatives of base flow solution velocity
93  /// components w.r.t. global coordinates (r and z) function
94  void (*Base_flow_dudx_fct_pt)(const double& time,
95  const Vector<double>& x,
96  DenseMatrix<double>& result);
97 
98  /// Boolean flag to indicate if ALE formulation is disabled when
99  /// the time-derivatives are computed. Only set to true if you're sure
100  /// that the mesh is stationary.
102 
103  /// Access function for the local equation number
104  /// information for the i-th component of the pressure.
105  /// p_local_eqn[n,i] = local equation number or < 0 if pinned.
106  virtual int p_local_eqn(const unsigned& n, const unsigned& i) = 0;
107 
108  /// Compute the shape functions and their derivatives
109  /// w.r.t. global coordinates at local coordinate s.
110  /// Return Jacobian of mapping between local and global coordinates.
112  const Vector<double>& s,
113  Shape& psi,
114  DShape& dpsidx,
115  Shape& test,
116  DShape& dtestdx) const = 0;
117 
118  /// Compute the shape functions and their derivatives
119  /// w.r.t. global coordinates at the ipt-th integration point.
120  /// Return Jacobian of mapping between local and global coordinates.
122  const unsigned& ipt,
123  Shape& psi,
124  DShape& dpsidx,
125  Shape& test,
126  DShape& dtestdx) const = 0;
127 
128  /// Compute the pressure shape functions at local coordinate s
130  Shape& psi) const = 0;
131 
132  /// Compute the pressure shape and test functions at local coordinate s
134  Shape& psi,
135  Shape& test) const = 0;
136 
137  /// Calculate the velocity components of the base flow solution
138  /// at a given time and Eulerian position
139  virtual void get_base_flow_u(const double& time,
140  const unsigned& ipt,
141  const Vector<double>& x,
142  Vector<double>& result) const
143  {
144  // If the function pointer is zero return zero
145  if (Base_flow_u_fct_pt == 0)
146  {
147  // Loop over velocity components and set base flow solution to zero
148  for (unsigned i = 0; i < 3; i++)
149  {
150  result[i] = 0.0;
151  }
152  }
153  // Otherwise call the function
154  else
155  {
156  (*Base_flow_u_fct_pt)(time, x, result);
157  }
158  }
159 
160  /// Calculate the derivatives of the velocity components of the
161  /// base flow solution w.r.t. global coordinates (r and z) at a given
162  /// time and Eulerian position
163  virtual void get_base_flow_dudx(const double& time,
164  const unsigned& ipt,
165  const Vector<double>& x,
166  DenseMatrix<double>& result) const
167  {
168  // If the function pointer is zero return zero
169  if (Base_flow_dudx_fct_pt == 0)
170  {
171  // Loop over velocity components
172  for (unsigned i = 0; i < 3; i++)
173  {
174  // Loop over coordinate directions and set to zero
175  for (unsigned j = 0; j < 2; j++)
176  {
177  result(i, j) = 0.0;
178  }
179  }
180  }
181  // Otherwise call the function
182  else
183  {
184  (*Base_flow_dudx_fct_pt)(time, x, result);
185  }
186  }
187 
188  /// Compute the residuals for the Navier-Stokes equations;
189  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
191  Vector<double>& residuals,
192  DenseMatrix<double>& jacobian,
193  DenseMatrix<double>& mass_matrix,
194  unsigned flag);
195 
196  public:
197  /// Constructor: NULL the base flow solution and the
198  /// derivatives of the base flow function
201  {
202  // Set all the physical parameter pointers to the default value of zero
205 
206  // Set the azimuthal mode number to default (zero)
208 
209  // Set the physical ratios to the default value of one
212  }
213 
214  /// Vector to decide whether the stress-divergence form is used or not.
215  // N.B. This needs to be public so that the intel compiler gets things
216  // correct. Somehow the access function messes things up when going to
217  // refineable navier--stokes
219 
220  // Access functions for the physical constants
221  // -------------------------------------------
222 
223  /// Reynolds number
224  const double& re() const
225  {
226  return *Re_pt;
227  }
228 
229  /// Product of Reynolds and Strouhal number (=Womersley number)
230  const double& re_st() const
231  {
232  return *ReSt_pt;
233  }
234 
235  /// Azimuthal mode number k in e^ik(theta) decomposition
236  const int& azimuthal_mode_number() const
237  {
238  return *Azimuthal_Mode_Number_pt;
239  }
240 
241  /// Pointer to Reynolds number
242  double*& re_pt()
243  {
244  return Re_pt;
245  }
246 
247  /// Pointer to product of Reynolds and Strouhal number (=Womersley number)
248  double*& re_st_pt()
249  {
250  return ReSt_pt;
251  }
252 
253  /// Pointer to azimuthal mode number k in e^ik(theta) decomposition
255  {
257  }
258 
259  /// Viscosity ratio for element: element's viscosity relative
260  /// to the viscosity used in the definition of the Reynolds number
261  const double& viscosity_ratio() const
262  {
263  return *Viscosity_Ratio_pt;
264  }
265 
266  /// Pointer to the viscosity ratio
268  {
269  return Viscosity_Ratio_pt;
270  }
271 
272  /// Density ratio for element: element's density relative
273  /// to the viscosity used in the definition of the Reynolds number
274  const double& density_ratio() const
275  {
276  return *Density_Ratio_pt;
277  }
278 
279  /// Pointer to the density ratio
280  double*& density_ratio_pt()
281  {
282  return Density_Ratio_pt;
283  }
284 
285  /// Access function for the base flow solution pointer
286  void (*&base_flow_u_fct_pt())(const double& time,
287  const Vector<double>& x,
288  Vector<double>& f)
289  {
290  return Base_flow_u_fct_pt;
291  }
292 
293  /// Access function for the derivatives of the base flow
294  /// w.r.t. global coordinates solution pointer
295  void (*&base_flow_dudx_fct_pt())(const double& time,
296  const Vector<double>& x,
298  {
299  return Base_flow_dudx_fct_pt;
300  }
301 
302  /// Return the number of pressure degrees of freedom
303  /// associated with a single pressure component in the element
304  virtual unsigned npres_linearised_axi_nst() const = 0;
305 
306  /// Return the index at which the i-th unknown velocity
307  /// component is stored. The default value, i, is appropriate for
308  /// single-physics problems. In derived multi-physics elements, this
309  /// function should be overloaded to reflect the chosen storage scheme.
310  /// Note that these equations require that the unknowns are always
311  /// stored at the same indices at each node.
312  virtual inline unsigned u_index_linearised_axi_nst(const unsigned& i) const
313  {
314  return i;
315  }
316 
317  /// Return the i-th component of du/dt at local node n.
318  /// Uses suitably interpolated value for hanging nodes.
319  double du_dt_linearised_axi_nst(const unsigned& n, const unsigned& i) const
320  {
321  // Get the data's timestepper
323 
324  // Initialise dudt
325  double dudt = 0.0;
326 
327  // Loop over the timesteps, if there is a non-steady timestepper
328  if (!time_stepper_pt->is_steady())
329  {
330  // Get the index at which the velocity is stored
331  const unsigned u_nodal_index = u_index_linearised_axi_nst(i);
332 
333  // Determine number of timsteps (past & present)
334  const unsigned n_time = time_stepper_pt->ntstorage();
335 
336  // Add the contributions to the time derivative
337  for (unsigned t = 0; t < n_time; t++)
338  {
339  dudt +=
340  time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
341  }
342  }
343 
344  return dudt;
345  }
346 
347  /// Disable ALE, i.e. assert the mesh is not moving -- you do this
348  /// at your own risk!
349  void disable_ALE()
350  {
351  ALE_is_disabled = true;
352  }
353 
354  /// (Re-)enable ALE, i.e. take possible mesh motion into account
355  /// when evaluating the time-derivative. Note: By default, ALE is
356  /// enabled, at the expense of possibly creating unnecessary work
357  /// in problems where the mesh is, in fact, stationary.
358  void enable_ALE()
359  {
360  ALE_is_disabled = false;
361  }
362 
363  /// Return the i-th pressure value at local pressure "node" n_p.
364  /// Uses suitably interpolated value for hanging nodes.
365  virtual double p_linearised_axi_nst(const unsigned& n_p,
366  const unsigned& i) const = 0;
367 
368  /// Which nodal value represents the pressure?
369  // N.B. This function has return type "int" (rather than "unsigned"
370  // as in the u_index case) so that we can return the "magic" number
371  // "Pressure_not_stored_at_node" ( = -100 )
372  virtual inline int p_index_linearised_axi_nst(const unsigned& i) const
373  {
375  }
376 
377  /// Strain-rate tensor: \f$ e_{ij} \f$
378  /// where \f$ i,j = r,z,\theta \f$ (in that order)
379  void strain_rate(const Vector<double>& s,
381 
382  /// Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
383  /// in tecplot format. Default number of plot points
384  void output(std::ostream& outfile)
385  {
386  const unsigned nplot = 5;
387  output(outfile, nplot);
388  }
389 
390  /// Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
391  /// in tecplot format. nplot points in each coordinate direction
392  void output(std::ostream& outfile, const unsigned& nplot);
393 
394  /// Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
395  /// in tecplot format. Default number of plot points
396  void output(FILE* file_pt)
397  {
398  const unsigned nplot = 5;
399  output(file_pt, nplot);
400  }
401 
402  /// Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
403  /// in tecplot format. nplot points in each coordinate direction
404  void output(FILE* file_pt, const unsigned& nplot);
405 
406  /// Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S,
407  /// in tecplot format. nplot points in each coordinate direction
408  /// at timestep t (t=0: present; t>0: previous timestep)
409  void output_veloc(std::ostream& outfile,
410  const unsigned& nplot,
411  const unsigned& t);
412 
413  /// Compute the element's residual Vector
415  {
416  // Call the generic residuals function with flag set to 0
417  // and using a dummy matrix argument
419  residuals,
422  0);
423  }
424 
425  /// Compute the element's residual Vector and the jacobian matrix.
426  /// Virtual function can be overloaded by hanging-node version.
428  DenseMatrix<double>& jacobian)
429  {
430  // Call the generic routine with the flag set to 1
432  residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
433  }
434 
435  /// Add the element's contribution to its residuals vector,
436  /// jacobian matrix and mass matrix
438  Vector<double>& residuals,
439  DenseMatrix<double>& jacobian,
440  DenseMatrix<double>& mass_matrix)
441  {
442  // Call the generic routine with the flag set to 2
444  residuals, jacobian, mass_matrix, 2);
445  }
446 
447  /// Return the i-th component of the FE interpolated velocity
448  /// u[i] at local coordinate s
450  const unsigned& i) const
451  {
452  // Determine number of nodes in the element
453  const unsigned n_node = nnode();
454 
455  // Provide storage for local shape functions
456  Shape psi(n_node);
457 
458  // Find values of shape functions
459  shape(s, psi);
460 
461  // Get the index at which the velocity is stored
462  const unsigned u_nodal_index = u_index_linearised_axi_nst(i);
463 
464  // Initialise value of u
465  double interpolated_u = 0.0;
466 
467  // Loop over the local nodes and sum
468  for (unsigned l = 0; l < n_node; l++)
469  {
470  interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
471  }
472 
473  return (interpolated_u);
474  }
475 
476  /// Return the i-th component of the FE interpolated pressure
477  /// p[i] at local coordinate s
479  const unsigned& i) const
480  {
481  // Determine number of pressure nodes in the element
482  const unsigned n_pressure_nodes = npres_linearised_axi_nst();
483 
484  // Provide storage for local shape functions
485  Shape psi(n_pressure_nodes);
486 
487  // Find values of shape functions
489 
490  // Initialise value of p
491  double interpolated_p = 0.0;
492 
493  // Loop over the local nodes and sum
494  for (unsigned l = 0; l < n_pressure_nodes; l++)
495  {
496  // N.B. The pure virtual function p_linearised_axi_nst(...)
497  // automatically calculates the index at which the pressure value
498  // is stored, so we don't need to worry about this here
499  interpolated_p += p_linearised_axi_nst(l, i) * psi[l];
500  }
501 
502  return (interpolated_p);
503  }
504 
505  }; // End of LinearisedAxisymmetricNavierStokesEquations class definition
506 
507 
508  /// ///////////////////////////////////////////////////////////////////////////
509  /// ///////////////////////////////////////////////////////////////////////////
510  /// ///////////////////////////////////////////////////////////////////////////
511 
512 
513  //=======================================================================
514  /// Crouzeix-Raviart elements are Navier-Stokes elements with quadratic
515  /// interpolation for velocities and positions, but a discontinuous
516  /// linear pressure interpolation
517  //=======================================================================
519  : public virtual QElement<2, 3>,
521  {
522  private:
523  /// Static array of ints to hold required number of variables at nodes
524  static const unsigned Initial_Nvalue[];
525 
526  protected:
527  /// Internal indices that indicate at which internal data the
528  /// pressure values are stored. We note that there are two pressure
529  /// values, corresponding to the functions P^C(r,z,t) and P^S(r,z,t)
530  /// which multiply the cosine and sine terms respectively.
532 
533  /// Velocity shape and test functions and their derivatives
534  /// w.r.t. global coordinates at local coordinate s (taken from geometry).
535  /// Return Jacobian of mapping between local and global coordinates.
537  const Vector<double>& s,
538  Shape& psi,
539  DShape& dpsidx,
540  Shape& test,
541  DShape& dtestdx) const;
542 
543  /// Velocity shape and test functions and their derivatives
544  /// w.r.t. global coordinates at the ipt-th integation point
545  /// (taken from geometry).
546  /// Return Jacobian of mapping between local and global coordinates.
548  const unsigned& ipt,
549  Shape& psi,
550  DShape& dpsidx,
551  Shape& test,
552  DShape& dtestdx) const;
553 
554  /// Compute the pressure shape functions at local coordinate s
555  inline void pshape_linearised_axi_nst(const Vector<double>& s,
556  Shape& psi) const;
557 
558  /// Compute the pressure shape and test functions at local coordinate s
559  inline void pshape_linearised_axi_nst(const Vector<double>& s,
560  Shape& psi,
561  Shape& test) const;
562 
563  public:
564  /// Constructor: there are three internal values for each
565  /// of the two pressure components
567  : QElement<2, 3>(),
570  {
571  // Loop over the two pressure components
572  for (unsigned i = 0; i < 2; i++)
573  {
574  // Allocate and add one internal data object for each of the two
575  // pressure components that store the three pressure values
577  this->add_internal_data(new Data(3));
578  }
579  }
580 
581  /// Return number of values (pinned or dofs) required at local node n
582  virtual unsigned required_nvalue(const unsigned& n) const;
583 
584  /// Return the pressure value i at internal dof i_internal
585  /// (Discontinous pressure interpolation -- no need to cater for
586  /// hanging nodes)
587  double p_linearised_axi_nst(const unsigned& i_internal,
588  const unsigned& i) const
589  {
591  ->value(i_internal);
592  }
593 
594  /// Return number of pressure values corresponding to a
595  /// single pressure component
596  unsigned npres_linearised_axi_nst() const
597  {
598  return 3;
599  }
600 
601  /// Fix both components of the internal pressure degrees
602  /// of freedom p_dof to pvalue
603  void fix_pressure(const unsigned& p_dof, const double& pvalue)
604  {
605  // Loop over the two pressure components
606  for (unsigned i = 0; i < 2; i++)
607  {
608  this->internal_data_pt(P_linearised_axi_nst_internal_index[i])
609  ->pin(p_dof);
611  ->set_value(p_dof, pvalue);
612  }
613  }
614 
615  /// Overload the access function for the i-th component of the
616  /// pressure's local equation numbers
617  inline int p_local_eqn(const unsigned& n, const unsigned& i)
618  {
620  }
621 
622  /// Redirect output to NavierStokesEquations output
623  void output(std::ostream& outfile)
624  {
626  }
627 
628  /// Redirect output to NavierStokesEquations output
629  void output(std::ostream& outfile, const unsigned& n_plot)
630  {
632  }
633 
634  /// Redirect output to NavierStokesEquations output
635  void output(FILE* file_pt)
636  {
638  }
639 
640  /// Redirect output to NavierStokesEquations output
641  void output(FILE* file_pt, const unsigned& n_plot)
642  {
644  }
645 
646  /// The number of "dof-blocks" that degrees of freedom in this
647  /// element are sub-divided into: Velocity and pressure.
648  unsigned ndof_types() const
649  {
650  return 8;
651  }
652 
653  }; // End of LinearisedAxisymmetricQCrouzeixRaviartElement class definition
654 
655 
656  // Inline functions
657  // ----------------
658 
659  //=======================================================================
660  /// Derivatives of the shape functions and test functions w.r.t.
661  /// global (Eulerian) coordinates at local coordinate s.
662  /// Return Jacobian of mapping between local and global coordinates.
663  //=======================================================================
666  Shape& psi,
667  DShape& dpsidx,
668  Shape& test,
669  DShape& dtestdx) const
670  {
671  // Call the geometrical shape functions and derivatives
672  const double J = this->dshape_eulerian(s, psi, dpsidx);
673 
674  // Loop over the test functions and derivatives and set them
675  // equal to the shape functions
676  for (unsigned i = 0; i < 9; i++)
677  {
678  test[i] = psi[i];
679  dtestdx(i, 0) = dpsidx(i, 0);
680  dtestdx(i, 1) = dpsidx(i, 1);
681  }
682 
683  // Return the Jacobian
684  return J;
685  }
686 
687  //=======================================================================
688  /// Derivatives of the shape functions and test functions w.r.t.
689  /// global (Eulerian) coordinates at the ipt-th integration point.
690  /// Return Jacobian of mapping between local and global coordinates.
691  //=======================================================================
694  Shape& psi,
695  DShape& dpsidx,
696  Shape& test,
697  DShape& dtestdx) const
698  {
699  // Call the geometrical shape functions and derivatives
700  const double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
701 
702  // Loop over the test functions and derivatives and set them
703  // equal to the shape functions
704  for (unsigned i = 0; i < 9; i++)
705  {
706  test[i] = psi[i];
707  dtestdx(i, 0) = dpsidx(i, 0);
708  dtestdx(i, 1) = dpsidx(i, 1);
709  }
710 
711  // Return the Jacobian
712  return J;
713  }
714 
715  //=======================================================================
716  /// Pressure shape functions
717  //=======================================================================
720  {
721  psi[0] = 1.0;
722  psi[1] = s[0];
723  psi[2] = s[1];
724  }
725 
726  //=======================================================================
727  /// Define the pressure shape and test functions
728  //=======================================================================
731  Shape& psi,
732  Shape& test) const
733  {
734  // Call the pressure shape functions
736 
737  // Loop over the test functions and set them equal to the shape functions
738  for (unsigned i = 0; i < 3; i++)
739  {
740  test[i] = psi[i];
741  }
742  }
743 
744  //=======================================================================
745  /// Face geometry of the linearised axisym Crouzeix-Raviart elements
746  //=======================================================================
747  template<>
749  : public virtual QElement<1, 3>
750  {
751  public:
752  FaceGeometry() : QElement<1, 3>() {}
753  };
754 
755  //=======================================================================
756  /// Face geometry of face geometry of the linearised axisymmetric
757  /// Crouzeix Raviart elements
758  //=======================================================================
759  template<>
762  : public virtual PointElement
763  {
764  public:
766  };
767 
768 
769  /// ///////////////////////////////////////////////////////////////////////////
770  /// ///////////////////////////////////////////////////////////////////////////
771  /// ///////////////////////////////////////////////////////////////////////////
772 
773 
774  //=======================================================================
775  /// Taylor--Hood elements are Navier--Stokes elements with quadratic
776  /// interpolation for velocities and positions and continuous linear
777  /// pressure interpolation
778  //=======================================================================
780  : public virtual QElement<2, 3>,
782  {
783  private:
784  /// Static array of ints to hold number of variables at node
785  static const unsigned Initial_Nvalue[];
786 
787  protected:
788  /// Static array of ints to hold conversion from pressure
789  /// node numbers to actual node numbers
790  static const unsigned Pconv[];
791 
792  /// Velocity shape and test functions and their derivatives
793  /// w.r.t. global coordinates at local coordinate s (taken from geometry).
794  /// Return Jacobian of mapping between local and global coordinates.
796  const Vector<double>& s,
797  Shape& psi,
798  DShape& dpsidx,
799  Shape& test,
800  DShape& dtestdx) const;
801 
802  /// Velocity shape and test functions and their derivatives
803  /// w.r.t. global coordinates the ipt-th integation point
804  /// (taken from geometry).
805  /// Return Jacobian of mapping between local and global coordinates.
807  const unsigned& ipt,
808  Shape& psi,
809  DShape& dpsidx,
810  Shape& test,
811  DShape& dtestdx) const;
812 
813  /// Compute the pressure shape functions at local coordinate s
814  inline void pshape_linearised_axi_nst(const Vector<double>& s,
815  Shape& psi) const;
816 
817  /// Compute the pressure shape and test functions at local coordinte s
818  inline void pshape_linearised_axi_nst(const Vector<double>& s,
819  Shape& psi,
820  Shape& test) const;
821 
822  public:
823  /// Constructor, no internal data points
826  {
827  }
828 
829  /// Number of values (pinned or dofs) required at node n. Can
830  /// be overwritten for hanging node version
831  inline virtual unsigned required_nvalue(const unsigned& n) const
832  {
833  return Initial_Nvalue[n];
834  }
835 
836  /// Which nodal value represents the pressure? Overload version
837  /// in base class which returns static int "Pressure_not_stored_at_node"
838  virtual int p_index_linearised_axi_nst(const unsigned& i) const
839  {
840  return (6 + i);
841  }
842 
843  /// Access function for the i-th component of pressure
844  /// at local pressure node n_p (const version).
845  double p_linearised_axi_nst(const unsigned& n_p, const unsigned& i) const
846  {
848  }
849 
850  /// Return number of pressure values corresponding to a
851  /// single pressure component
852  unsigned npres_linearised_axi_nst() const
853  {
854  return 4;
855  }
856 
857  /// Fix both components of the pressure at local pressure
858  /// node n_p to pvalue
859  void fix_pressure(const unsigned& n_p, const double& pvalue)
860  {
861  // Loop over the two pressure components
862  for (unsigned i = 0; i < 2; i++)
863  {
864  this->node_pt(Pconv[n_p])->pin(p_index_linearised_axi_nst(i));
865  this->node_pt(Pconv[n_p])
867  }
868  }
869 
870  /// Overload the access function for the i-th component of the
871  /// pressure's local equation numbers
872  inline int p_local_eqn(const unsigned& n, const unsigned& i)
873  {
875  }
876 
877  /// Redirect output to NavierStokesEquations output
878  void output(std::ostream& outfile)
879  {
881  }
882 
883  /// Redirect output to NavierStokesEquations output
884  void output(std::ostream& outfile, const unsigned& n_plot)
885  {
887  }
888 
889  /// Redirect output to NavierStokesEquations output
890  void output(FILE* file_pt)
891  {
893  }
894 
895  /// Redirect output to NavierStokesEquations output
896  void output(FILE* file_pt, const unsigned& n_plot)
897  {
899  }
900 
901  /// Returns the number of "dof-blocks" that degrees of freedom
902  /// in this element are sub-divided into: Velocity and pressure.
903  unsigned ndof_types() const
904  {
905  return 8;
906  }
907 
908  }; // End of LinearisedAxisymmetricQTaylorHoodElement class definition
909 
910 
911  // Inline functions
912  // ----------------
913 
914  //=======================================================================
915  /// Derivatives of the shape functions and test functions w.r.t
916  /// global (Eulerian) coordinates at local coordinate s.
917  /// Return Jacobian of mapping between local and global coordinates.
918  //=======================================================================
921  Shape& psi,
922  DShape& dpsidx,
923  Shape& test,
924  DShape& dtestdx) const
925  {
926  // Call the geometrical shape functions and derivatives
927  const double J = this->dshape_eulerian(s, psi, dpsidx);
928 
929  // Loop over the test functions and derivatives and set them
930  // equal to the shape functions
931  for (unsigned i = 0; i < 9; i++)
932  {
933  test[i] = psi[i];
934  dtestdx(i, 0) = dpsidx(i, 0);
935  dtestdx(i, 1) = dpsidx(i, 1);
936  }
937 
938  // Return the Jacobian
939  return J;
940  }
941 
942  //=======================================================================
943  /// Derivatives of the shape functions and test functions w.r.t
944  /// global (Eulerian) coordinates at the ipt-th integration point.
945  /// Return Jacobian of mapping between local and global coordinates.
946  //=======================================================================
949  Shape& psi,
950  DShape& dpsidx,
951  Shape& test,
952  DShape& dtestdx) const
953  {
954  // Call the geometrical shape functions and derivatives
955  const double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
956 
957  // Loop over the test functions and derivatives and set them
958  // equal to the shape functions
959  for (unsigned i = 0; i < 9; i++)
960  {
961  test[i] = psi[i];
962  dtestdx(i, 0) = dpsidx(i, 0);
963  dtestdx(i, 1) = dpsidx(i, 1);
964  }
965 
966  // Return the Jacobian
967  return J;
968  }
969 
970  //=======================================================================
971  /// Pressure shape functions
972  //=======================================================================
975  {
976  // Allocate local storage for the pressure shape functions
977  double psi1[2], psi2[2];
978 
979  // Call the one-dimensional shape functions
980  OneDimLagrange::shape<2>(s[0], psi1);
981  OneDimLagrange::shape<2>(s[1], psi2);
982 
983  // Now let's loop over the nodal points in the element
984  // s1 is the "r" coordinate, s2 the "z"
985  for (unsigned i = 0; i < 2; i++)
986  {
987  for (unsigned j = 0; j < 2; j++)
988  {
989  // Multiply the two 1D functions together to get the 2D function
990  psi[2 * i + j] = psi2[i] * psi1[j];
991  }
992  }
993  }
994 
995  //=======================================================================
996  /// Pressure shape and test functions
997  //=======================================================================
1000  Shape& psi,
1001  Shape& test) const
1002  {
1003  // Call the pressure shape functions
1005 
1006  // Loop over the test functions and set them equal to the shape functions
1007  for (unsigned i = 0; i < 4; i++)
1008  {
1009  test[i] = psi[i];
1010  }
1011  }
1012 
1013  //=======================================================================
1014  /// Face geometry of the linearised axisymmetric Taylor Hood elements
1015  //=======================================================================
1016  template<>
1018  : public virtual QElement<1, 3>
1019  {
1020  public:
1021  FaceGeometry() : QElement<1, 3>() {}
1022  };
1023 
1024  //=======================================================================
1025  /// Face geometry of the face geometry of the linearised
1026  /// axisymmetric Taylor Hood elements
1027  //=======================================================================
1028  template<>
1030  : public virtual PointElement
1031  {
1032  public:
1034  };
1035 
1036 
1037 } // namespace oomph
1038 
1039 #endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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...
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 nnode() const
Return the number of nodes.
Definition: elements.h:2210
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
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
A class for elements that solve the linearised version of the unsteady Navier–Stokes equations in cyl...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element's residual Vector.
virtual double dshape_and_dtest_eulerian_linearised_axi_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and their derivatives w.r.t. global coordinates at local coordinate s....
void output(FILE *file_pt)
Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S in tecplot format. Default number of pl...
int *& azimuthal_mode_number_pt()
Pointer to azimuthal mode number k in e^ik(theta) decomposition.
void(*&)(const double &time, const Vector< double > &x, DenseMatrix< double > &f) base_flow_dudx_fct_pt()
Access function for the derivatives of the base flow w.r.t. global coordinates solution pointer.
void(* Base_flow_dudx_fct_pt)(const double &time, const Vector< double > &x, DenseMatrix< double > &result)
Pointer to derivatives of base flow solution velocity components w.r.t. global coordinates (r and z) ...
const double & viscosity_ratio() const
Viscosity ratio for element: element's viscosity relative to the viscosity used in the definition of ...
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 * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
static int Default_Azimuthal_Mode_Number_Value
Static default value for the azimuthal mode number (zero)
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
int * Azimuthal_Mode_Number_pt
Pointer to azimuthal mode number k in e^ik(theta) decomposition.
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
double du_dt_linearised_axi_nst(const unsigned &n, const unsigned &i) const
Return the i-th component of du/dt at local node n. Uses suitably interpolated value for hanging node...
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 ...
void output(std::ostream &outfile)
Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S in tecplot format. Default number of pl...
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all initialised to one)
virtual void get_base_flow_u(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result) const
Calculate the velocity components of the base flow solution at a given time and Eulerian position.
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
virtual int p_local_eqn(const unsigned &n, const unsigned &i)=0
Access function for the local equation number information for the i-th component of the pressure....
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed....
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative....
virtual void pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
double interpolated_u_linearised_axi_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated velocity u[i] at local coordinate s.
void(*&)(const double &time, const Vector< double > &x, Vector< double > &f) base_flow_u_fct_pt()
Access function for the base flow solution pointer.
void output_veloc(std::ostream &outfile, const unsigned &nplot, const unsigned &t)
Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, in tecplot format. nplot points in each coordina...
virtual void get_base_flow_dudx(const double &time, const unsigned &ipt, const Vector< double > &x, DenseMatrix< double > &result) const
Calculate the derivatives of the velocity components of the base flow solution w.r....
virtual double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and their derivatives w.r.t. global coordinates at the ipt-th integration...
virtual unsigned npres_linearised_axi_nst() const =0
Return the number of pressure degrees of freedom associated with a single pressure component in the e...
double interpolated_p_linearised_axi_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated pressure p[i] at local coordinate s.
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
const double & density_ratio() const
Density ratio for element: element's density relative to the viscosity used in the definition of the ...
virtual unsigned u_index_linearised_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
const int & azimuthal_mode_number() const
Azimuthal mode number k in e^ik(theta) decomposition.
void(* Base_flow_u_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to base flow solution (velocity components) function.
virtual void pshape_linearised_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_residual_contribution_linearised_axi_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier-Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
virtual double p_linearised_axi_nst(const unsigned &n_p, const unsigned &i) const =0
Return the i-th pressure value at local pressure "node" n_p. Uses suitably interpolated value for han...
virtual int p_index_linearised_axi_nst(const unsigned &i) const
Which nodal value represents the pressure?
LinearisedAxisymmetricNavierStokesEquations()
Constructor: NULL the base flow solution and the derivatives of the base flow function.
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void output(std::ostream &outfile, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
unsigned ndof_types() const
The number of "dof-blocks" that degrees of freedom in this element are sub-divided into: Velocity and...
double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivatives w.r.t. global coordinates at the ipt-th integ...
double p_linearised_axi_nst(const unsigned &i_internal, const unsigned &i) const
Return the pressure value i at internal dof i_internal (Discontinous pressure interpolation – no need...
Vector< unsigned > P_linearised_axi_nst_internal_index
Internal indices that indicate at which internal data the pressure values are stored....
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
LinearisedAxisymmetricQCrouzeixRaviartElement()
Constructor: there are three internal values for each of the two pressure components.
double dshape_and_dtest_eulerian_linearised_axi_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivatives w.r.t. global coordinates at local coordinate...
virtual unsigned required_nvalue(const unsigned &n) const
Return number of values (pinned or dofs) required at local node n.
void pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi) const
Compute the pressure shape functions at local coordinate s.
void fix_pressure(const unsigned &p_dof, const double &pvalue)
Fix both components of the internal pressure degrees of freedom p_dof to pvalue.
unsigned npres_linearised_axi_nst() const
Return number of pressure values corresponding to a single pressure component.
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
int p_local_eqn(const unsigned &n, const unsigned &i)
Overload the access function for the i-th component of the pressure's local equation numbers.
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
unsigned ndof_types() const
Returns the number of "dof-blocks" that degrees of freedom in this element are sub-divided into: Velo...
double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivatives w.r.t. global coordinates the ipt-th integati...
void fix_pressure(const unsigned &n_p, const double &pvalue)
Fix both components of the pressure at local pressure node n_p to pvalue.
void pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi) const
Compute the pressure shape functions at local coordinate s.
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.
double p_linearised_axi_nst(const unsigned &n_p, const unsigned &i) const
Access function for the i-th component of pressure at local pressure node n_p (const version).
void output(std::ostream &outfile, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
unsigned npres_linearised_axi_nst() const
Return number of pressure values corresponding to a single pressure component.
virtual int p_index_linearised_axi_nst(const unsigned &i) const
Which nodal value represents the pressure? Overload version in base class which returns static int "P...
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
double dshape_and_dtest_eulerian_linearised_axi_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivatives w.r.t. global coordinates at local coordinate...
LinearisedAxisymmetricQTaylorHoodElement()
Constructor, no internal data points.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
int p_local_eqn(const unsigned &n, const unsigned &i)
Overload the access function for the i-th component of the pressure's local equation numbers.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: elements.h:3439
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
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
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:608
//////////////////////////////////////////////////////////////////// ////////////////////////////////...