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