spherical_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 Navier Stokes elements
31 
32 #ifndef OOMPH_SPHERICAL_NAVIER_STOKES_ELEMENTS_HEADER
33 #define OOMPH_SPHERICAL_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 
41 //OOMPH-LIB headers
42 #include "../generic/Qelements.h"
43 #include "../generic/fsi.h"
44 //#include "generic/block_preconditioner.h"
45 
46 namespace oomph
47 {
48 
49 
50 //======================================================================
51 /// A class for elements that solve the Navier--Stokes equations,
52 /// in axisymmetric spherical polar coordinates.
53 /// This contains the generic maths -- any concrete implementation must
54 /// be derived from this.
55 ///
56 /// We also provide all functions required to use this element
57 /// in FSI problems, by deriving it from the FSIFluidElement base
58 /// class.
59 //======================================================================
62 {
63 
64  public:
65 
66  /// \short Function pointer to body force function fct(t,x,f(x))
67  /// x is a Vector!
69  const double& time,const Vector<double>& x,Vector<double>& body_force);
70 
71  /// \short Function pointer to source function fct(t,x)
72  /// x is a Vector!
73  typedef double (*SphericalNavierStokesSourceFctPt)(const double& time,
74  const Vector<double>& x);
75 
76  private:
77 
78  /// \short Static "magic" number that indicates that the pressure is
79  /// not stored at a node
81 
82  /// Static default value for the physical constants (all initialised to zero)
84 
85  /// Static default value for the physical ratios (all are initialised to one)
87 
88  /// Static default value for the gravity vector
90 
91  protected:
92 
93  //Physical constants
94 
95  /// \short Pointer to the viscosity ratio (relative to the
96  /// viscosity used in the definition of the Reynolds number)
98 
99  /// \short Pointer to the density ratio (relative to the density used in the
100  /// definition of the Reynolds number)
102 
103  // Pointers to global physical constants
104 
105  /// Pointer to global Reynolds number
106  double *Re_pt;
107 
108  /// Pointer to global Reynolds number x Strouhal number (=Womersley)
109  double *ReSt_pt;
110 
111  /// \short Pointer to global Reynolds number x inverse Froude number
112  /// (= Bond number / Capillary number)
113  double *ReInvFr_pt;
114 
115  /// \short Pointer to global Reynolds number x inverse Rossby number
116  /// (used when in a rotating frame)
117  double *ReInvRo_pt;
118 
119  /// Pointer to global gravity Vector
121 
122  /// Pointer to body force function
124 
125  /// Pointer to volumetric source function
127 
128  /// \short Boolean flag to indicate if ALE formulation is disabled when
129  /// time-derivatives are computed. Only set to true if you're sure
130  /// that the mesh is stationary.
132 
133  /// \short Access function for the local equation number information for
134  /// the pressure.
135  /// p_local_eqn[n] = local equation number or < 0 if pinned
136  virtual int p_local_eqn(const unsigned &n)const=0;
137 
138  /// \short Compute the shape functions and derivatives
139  /// w.r.t. global coords at local coordinate s.
140  /// Return Jacobian of mapping between local and global coordinates.
142  const Vector<double> &s,
143  Shape &psi,
144  DShape &dpsidx, Shape &test,
145  DShape &dtestdx) const=0;
146 
147  /// \short Compute the shape functions and derivatives
148  /// w.r.t. global coords at ipt-th integration point
149  /// Return Jacobian of mapping between local and global coordinates.
151  const unsigned &ipt,
152  Shape &psi,
153  DShape &dpsidx,
154  Shape &test,
155  DShape &dtestdx) const=0;
156 
157  /// Compute the pressure shape functions at local coordinate s
158  virtual void pshape_spherical_nst(const Vector<double> &s, Shape &psi)
159  const=0;
160 
161  /// \short Compute the pressure shape and test functions
162  /// at local coordinate s
163  virtual void pshape_spherical_nst(const Vector<double> &s, Shape &psi,
164  Shape &test) const=0;
165 
166 
167  /// \short Calculate the body force at a given time and local and/or
168  /// Eulerian position. This function is virtual so that it can be
169  /// overloaded in multi-physics elements where the body force might
170  /// depend on another variable.
171  virtual void get_body_force_spherical_nst(const double &time,
172  const unsigned& ipt,
173  const Vector<double> &s,
174  const Vector<double> &x,
175  Vector<double> &result)
176  {
177  //If the function pointer is zero return zero
178  if(Body_force_fct_pt == 0)
179  {
180  //Loop over three spatial dimensions and set body forces to zero
181  for(unsigned i=0;i<3;i++) {result[i] = 0.0;}
182  }
183  //Otherwise call the function
184  else
185  {
186  (*Body_force_fct_pt)(time,x,result);
187  }
188  }
189 
190  /// \short Calculate the source fct at given time and
191  /// Eulerian position
192  virtual double get_source_spherical_nst(double time, const unsigned& ipt,
193  const Vector<double> &x)
194  {
195  //If the function pointer is zero return zero
196  if (Source_fct_pt == 0) {return 0;}
197  //Otherwise call the function
198  else {return (*Source_fct_pt)(time,x);}
199  }
200 
201  ///\short Compute the residuals for the Navier--Stokes equations;
202  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
204  Vector<double> &residuals, DenseMatrix<double> &jacobian,
205  DenseMatrix<double> &mass_matrix, unsigned flag);
206 
207 
208 public:
209 
210  /// Include a cot function to simplify the
211  /// NS momentum and jacobian expressions
212  inline double cot(const double &th) const {return(1/tan(th));}
213 
214 
215  /// \short Constructor: NULL the body force and source function
216  /// and make sure the ALE terms are included by default.
218  ALE_is_disabled(false)
219  {
220  //Set all the Physical parameter pointers to the default value zero
226  //Set the Physical ratios to the default value of 1
229  }
230 
231  /// Vector to decide whether the stress-divergence form is used or not
232  // N.B. This needs to be public so that the intel compiler gets things correct
233  // somehow the access function messes things up when going to refineable
234  // navier--stokes
236 
237  //Access functions for the physical constants
238 
239  /// Reynolds number
240  const double &re() const {return *Re_pt;}
241 
242  /// Product of Reynolds and Strouhal number (=Womersley number)
243  const double &re_st() const {return *ReSt_pt;}
244 
245  /// Pointer to Reynolds number
246  double* &re_pt() {return Re_pt;}
247 
248  /// Pointer to product of Reynolds and Strouhal number (=Womersley number)
249  double* &re_st_pt() {return ReSt_pt;}
250 
251  /// \short Viscosity ratio for element: Element's viscosity relative to the
252  /// viscosity used in the definition of the Reynolds number
253  const double &viscosity_ratio() const {return *Viscosity_Ratio_pt;}
254 
255  /// Pointer to Viscosity Ratio
257 
258  /// \short Density ratio for element: Element's density relative to the
259  /// viscosity used in the definition of the Reynolds number
260  const double &density_ratio() const {return *Density_Ratio_pt;}
261 
262  /// Pointer to Density ratio
263  double* &density_ratio_pt() {return Density_Ratio_pt;}
264 
265  /// Global inverse Froude number
266  const double &re_invfr() const {return *ReInvFr_pt;}
267 
268  /// Pointer to global inverse Froude number
269  double* &re_invfr_pt() {return ReInvFr_pt;}
270 
271  /// Global Reynolds number multiplied by inverse Rossby number
272  const double &re_invro() const {return *ReInvRo_pt;}
273 
274  /// Pointer to global inverse Froude number
275  double* &re_invro_pt() {return ReInvRo_pt;}
276 
277  /// Vector of gravitational components
278  const Vector<double> &g() const {return *G_pt;}
279 
280  /// Pointer to Vector of gravitational components
281  Vector<double>* &g_pt() {return G_pt;}
282 
283  /// Access function for the body-force pointer
285  {return Body_force_fct_pt;}
286 
287  /// Access function for the body-force pointer. Const version
289  {return Body_force_fct_pt;}
290 
291  ///Access function for the source-function pointer
293 
294  ///Access function for the source-function pointer. Const version
296 
297  ///Function to return number of pressure degrees of freedom
298  virtual unsigned npres_spherical_nst() const=0;
299 
300  /// \short Velocity i at local node n. Uses suitably interpolated value
301  /// for hanging nodes.
302  /// The use of u_index_spherical_nst() permits the use of this
303  /// element as the basis for multi-physics elements. The default
304  /// is to assume that the i-th velocity component is stored at the
305  /// i-th location of the node
306  double u_spherical_nst(const unsigned &n, const unsigned &i) const
307  {return nodal_value(n,u_index_spherical_nst(i));}
308 
309  /// \short Velocity i at local node n at timestep t (t=0: present;
310  /// t>0: previous). Uses suitably interpolated value for hanging nodes.
311  double u_spherical_nst(const unsigned &t, const unsigned &n,
312  const unsigned &i) const
313  {return nodal_value(t,n,u_index_spherical_nst(i));}
314 
315  /// \short Return the index at which the i-th unknown velocity component
316  /// is stored. The default value, i, is appropriate for single-physics
317  /// problems.
318  /// In derived multi-physics elements, this function should be overloaded
319  /// to reflect the chosen storage scheme. Note that these equations require
320  /// that the unknowns are always stored at the same indices at each node.
321  virtual inline unsigned u_index_spherical_nst(const unsigned &i)
322  const {return i;}
323 
324 
325  /// \short i-th component of du/dt at local node n.
326  /// Uses suitably interpolated value for hanging nodes.
327  double du_dt_spherical_nst(const unsigned &n, const unsigned &i) const
328  {
329  // Get the data's timestepper
331 
332  //Initialise dudt
333  double dudt=0.0;
334 
335  //Loop over the timesteps, if there is a non Steady timestepper
336  if (time_stepper_pt->type()!="Steady")
337  {
338  //Find the index at which the dof is stored
339  const unsigned u_nodal_index = this->u_index_spherical_nst(i);
340 
341  // Number of timsteps (past & present)
342  const unsigned n_time = time_stepper_pt->ntstorage();
343  // Loop over the timesteps
344  for(unsigned t=0;t<n_time;t++)
345  {
346  dudt+=time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
347  }
348  }
349 
350  return dudt;
351  }
352 
353  /// \short Disable ALE, i.e. assert the mesh is not moving -- you do this
354  /// at your own risk!
355  void disable_ALE()
356  {
357  ALE_is_disabled=true;
358  }
359 
360  /// \short (Re-)enable ALE, i.e. take possible mesh motion into account
361  /// when evaluating the time-derivative. Note: By default, ALE is
362  /// enabled, at the expense of possibly creating unnecessary work
363  /// in problems where the mesh is, in fact, stationary.
364  void enable_ALE()
365  {
366  ALE_is_disabled=false;
367  }
368 
369  /// \short Pressure at local pressure "node" n_p
370  /// Uses suitably interpolated value for hanging nodes.
371  virtual double p_spherical_nst(const unsigned &n_p)const=0;
372 
373  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
374  virtual void fix_pressure(const unsigned &p_dof, const double &p_value)=0;
375 
376  /// \short Return the index at which the pressure is stored if it is
377  /// stored at the nodes. If not stored at the nodes this will return
378  /// a negative number.
379  virtual int p_nodal_index_spherical_nst() const
381 
382  /// Integral of pressure over element
383  double pressure_integral() const;
384 
385  /// \short Return integral of dissipation over element
386  double dissipation() const;
387 
388  /// \short Return dissipation at local coordinate s
389  double dissipation(const Vector<double>& s) const;
390 
391  /// \short Compute the vorticity vector at local coordinate s
392  void get_vorticity(const Vector<double>& s, Vector<double>& vorticity) const;
393 
394  /// \short Get integral of kinetic energy over element
395  double kin_energy() const;
396 
397  /// \short Get integral of time derivative of kinetic energy over element
398  double d_kin_energy_dt() const;
399 
400  /// Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)
401  void strain_rate(const Vector<double>& s,
403 
404  /// \short Compute traction (on the viscous scale) exerted onto
405  /// the fluid at local coordinate s. N has to be outer unit normal
406  /// to the fluid.
407  void get_traction(const Vector<double>& s, const Vector<double>& N,
408  Vector<double>& traction);
409 
410 
411  /// \short This implements a pure virtual function defined
412  /// in the FSIFluidElement class. The function computes
413  /// the traction (on the viscous scale), at the element's local
414  /// coordinate s, that the fluid element exerts onto an adjacent
415  /// solid element. The number of arguments is imposed by
416  /// the interface defined in the FSIFluidElement -- only the
417  /// unit normal N (pointing into the fluid!) is actually used
418  /// in the computation.
419  void get_load(const Vector<double> &s,
420  const Vector<double> &N,
421  Vector<double> &load)
422  {
423  // Note: get_traction() computes the traction onto the fluid
424  // if N is the outer unit normal onto the fluid; here we're
425  // exepcting N to point into the fluid so we're getting the
426  // traction onto the adjacent wall instead!
427  get_traction(s,N,load);
428  }
429 
430  /// \short Compute the diagonal of the velocity/pressure mass matrices.
431  /// If which one=0, both are computed, otherwise only the pressure
432  /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
433  /// LSC version of the preconditioner only needs that one)
434  /// NOTE: pressure versions isn't implemented yet because this
435  /// element has never been tried with Fp preconditoner.
437  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
438  const unsigned& which_one=0);
439 
440 
441  /// \short Output function: x,y,[z],u,v,[w],p
442  /// in tecplot format. Default number of plot points
443  void output(std::ostream &outfile)
444  {
445  unsigned nplot=5;
446  output(outfile,nplot);
447  }
448 
449  /// \short Output function: x,y,[z],u,v,[w],p
450  /// in tecplot format. nplot points in each coordinate direction
451  void output(std::ostream &outfile, const unsigned &nplot);
452 
453  /// \short C-style output function: x,y,[z],u,v,[w],p
454  /// in tecplot format. Default number of plot points
455  void output(FILE* file_pt)
456  {
457  unsigned nplot=5;
458  output(file_pt,nplot);
459  }
460 
461  /// \short C-style output function: x,y,[z],u,v,[w],p
462  /// in tecplot format. nplot points in each coordinate direction
463  void output(FILE* file_pt, const unsigned &nplot);
464 
465  /// \short Full output function:
466  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
467  /// in tecplot format. Default number of plot points
468  void full_output(std::ostream &outfile)
469  {
470  unsigned nplot=5;
471  full_output(outfile,nplot);
472  }
473 
474  /// \short Full output function:
475  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
476  /// in tecplot format. nplot points in each coordinate direction
477  void full_output(std::ostream &outfile, const unsigned &nplot);
478 
479 
480  /// \short Output function: x,y,[z],u,v,[w] in tecplot format.
481  /// nplot points in each coordinate direction at timestep t
482  /// (t=0: present; t>0: previous timestep)
483  void output_veloc(std::ostream &outfile, const unsigned &nplot,
484  const unsigned& t);
485 
486 
487  /// \short Output function: x,y,[z], [omega_x,omega_y,[and/or omega_z]]
488  /// in tecplot format. nplot points in each coordinate direction
489  void output_vorticity(std::ostream &outfile,
490  const unsigned &nplot);
491 
492  /// \short Output exact solution specified via function pointer
493  /// at a given number of plot points. Function prints as
494  /// many components as are returned in solution Vector
495  void output_fct(std::ostream &outfile, const unsigned &nplot,
497 
498  /// \short Output exact solution specified via function pointer
499  /// at a given time and at a given number of plot points.
500  /// Function prints as many components as are returned in solution Vector.
501  void output_fct(std::ostream &outfile, const unsigned &nplot,
502  const double& time,
504 
505  /// \short Validate against exact solution at given time
506  /// Solution is provided via function pointer.
507  /// Plot at a given number of plot points and compute L2 error
508  /// and L2 norm of velocity solution over element
509  void compute_error(std::ostream &outfile,
511  const double& time,
512  double& error, double& norm);
513 
514  /// \short Validate against exact solution.
515  /// Solution is provided via function pointer.
516  /// Plot at a given number of plot points and compute L2 error
517  /// and L2 norm of velocity solution over element
518  void compute_error(std::ostream &outfile,
520  double& error, double& norm);
521 
522  /// Validate against exact solution.
523  /// Solution is provided direct from exact_soln function.
524  /// Plot at a given number of plot points and compute the energy error
525  /// and energy norm of the velocity solution over the element.
526  void compute_error_e(std::ostream &outfile,
530  exact_soln_dtheta_pt,
531  double& u_error, double& u_norm,
532  double& p_error, double& p_norm);
533 
534  // Listing for the shear stress function
535  void compute_shear_stress(std::ostream &outfile);
536 
537  // Listing for the velocity extraction function
538  void extract_velocity(std::ostream &outfile);
539 
540  // Calculate the analytic solution at the point x and output a vector
542  {
543  const double Re = this->re();
544 
545  double r = x[0];
546  double theta = x[1];
547  Vector<double> ans(4,0.0);
548 
549  ans[2] = r*sin(theta);
550  ans[3] = 0.5*Re*r*r*sin(theta)*sin(theta);
551 
552  return(ans);
553  }
554 
555  // Calculate the r-derivatives of the analytic solution at the point x and output a vector
557  {
558  const double Re = this->re();
559 
560  double r = x[0];
561  double theta = x[1];
562  Vector<double> ans(4,0.0);
563 
564  ans[2] = sin(theta);
565  ans[3] = Re*r*sin(theta)*sin(theta);
566 
567  return(ans);
568  }
569 
570  // Calculate the theta-derivatives of the analytic solution at the point x and output a vector
572  {
573  const double Re = this->re();
574 
575  double r = x[0];
576  double theta = x[1];
577  Vector<double> ans(4,0.0);
578 
579  ans[2] = r*cos(theta);
580  ans[3] = Re*r*r*sin(theta)*cos(theta);
581 
582  return(ans);
583  }
584 
585  /// Compute the element's residual Vector
587  {
588  //Call the generic residuals function with flag set to 0
589  //and using a dummy matrix argument
593  }
594 
595  //\short Compute the element's residual Vector and the jacobian matrix
596  // Virtual function can be overloaded by hanging-node version
598  DenseMatrix<double> &jacobian)
599  {
600  //Call the generic routine with the flag set to 1
602  residuals, jacobian,
604  }
605 
606  /// Add the element's contribution to its residuals vector,
607  /// jacobian matrix and mass matrix
609  Vector<double> &residuals, DenseMatrix<double> &jacobian,
610  DenseMatrix<double> &mass_matrix)
611  {
612  //Call the generic routine with the flag set to 2
614  jacobian,mass_matrix,2);
615  }
616 
617  /// Compute vector of FE interpolated velocity u at local coordinate s
619  {
620  //Find number of nodes
621  unsigned n_node = nnode();
622  //Local shape function
623  Shape psi(n_node);
624  //Find values of shape function
625  shape(s,psi);
626 
627  for (unsigned i=0;i<3;i++)
628  {
629  //Index at which the nodal value is stored
630  unsigned u_nodal_index = u_index_spherical_nst(i);
631  //Initialise value of u
632  veloc[i] = 0.0;
633  //Loop over the local nodes and sum
634  for(unsigned l=0;l<n_node;l++)
635  {
636  veloc[i] += nodal_value(l,u_nodal_index)*psi[l];
637  }
638  }
639  }
640 
641  /// Return FE interpolated velocity u[i] at local coordinate s
642  double interpolated_u_spherical_nst(const Vector<double> &s, const unsigned &i) const
643  {
644  //Find number of nodes
645  unsigned n_node = nnode();
646  //Local shape function
647  Shape psi(n_node);
648  //Find values of shape function
649  shape(s,psi);
650 
651  //Get nodal index at which i-th velocity is stored
652  unsigned u_nodal_index = u_index_spherical_nst(i);
653 
654  //Initialise value of u
655  double interpolated_u = 0.0;
656  //Loop over the local nodes and sum
657  for(unsigned l=0;l<n_node;l++)
658  {
659  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
660  }
661 
662  return(interpolated_u);
663  }
664 
665  /// Return matrix entry dudx(i,j) of the FE interpolated velocity derivative at local coordinate s
667  const unsigned &i, const unsigned &j) const
668  {
669  //Find number of nodes
670  unsigned n_node = nnode();
671  //Local shape function
672  Shape psi(n_node);
673  DShape dpsidx(n_node,2);
674  //Find values of shape function
675  (void) dshape_eulerian(s,psi,dpsidx);
676 
677  //Get nodal index at which i-th velocity is stored
678  unsigned u_nodal_index = u_index_spherical_nst(i);
679 
680  //Initialise value of u
681  double interpolated_dudx = 0.0;
682  //Loop over the local nodes and sum
683  for(unsigned l=0;l<n_node;l++)
684  {
685  interpolated_dudx += nodal_value(l,u_nodal_index)*dpsidx(l,j);
686  }
687 
688  return(interpolated_dudx);
689  }
690 
691  /// Return FE interpolated pressure at local coordinate s
693  {
694  //Find number of nodes
695  unsigned n_pres = npres_spherical_nst();
696  //Local shape function
697  Shape psi(n_pres);
698  //Find values of shape function
699  pshape_spherical_nst(s,psi);
700 
701  //Initialise value of p
702  double interpolated_p = 0.0;
703  //Loop over the local nodes and sum
704  for(unsigned l=0;l<n_pres;l++)
705  {
706  interpolated_p += p_spherical_nst(l)*psi[l];
707  }
708 
709  return(interpolated_p);
710  }
711 
712 
713 };
714 
715 //////////////////////////////////////////////////////////////////////////////
716 //////////////////////////////////////////////////////////////////////////////
717 //////////////////////////////////////////////////////////////////////////////
718 
719 
720 //==========================================================================
721 /// Crouzeix_Raviart elements are Navier--Stokes elements with quadratic
722 /// interpolation for velocities and positions, but a discontinuous linear
723 /// pressure interpolation. They can be used within oomph-lib's
724 /// block preconditioning framework.
725 //==========================================================================
726 class QSphericalCrouzeixRaviartElement : public virtual QElement<2,3>,
727  public virtual SphericalNavierStokesEquations
728 {
729  private:
730 
731  /// Static array of ints to hold required number of variables at nodes
732  static const unsigned Initial_Nvalue[];
733 
734  protected:
735 
736  /// Internal index that indicates at which internal data the pressure
737  /// is stored
739 
740 
741  /// \short Velocity shape and test functions and their derivs
742  /// w.r.t. to global coords at local coordinate s (taken from geometry)
743  ///Return Jacobian of mapping between local and global coordinates.
745  Shape &psi,
746  DShape &dpsidx, Shape &test,
747  DShape &dtestdx) const;
748 
749  /// \short Velocity shape and test functions and their derivs
750  /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
751  ///Return Jacobian of mapping between local and global coordinates.
752  inline double dshape_and_dtest_eulerian_at_knot_spherical_nst(const unsigned &ipt,
753  Shape &psi,
754  DShape &dpsidx,
755  Shape &test,
756  DShape &dtestdx) const;
757 
758  /// Pressure shape functions at local coordinate s
759  inline void pshape_spherical_nst(const Vector<double> &s, Shape &psi) const;
760 
761  /// Pressure shape and test functions at local coordinte s
762  inline void pshape_spherical_nst(const Vector<double> &s, Shape &psi,
763  Shape &test) const;
764 
765  /// Return the local equation numbers for the pressure values.
766  inline int p_local_eqn(const unsigned &n)const
768 
769 public:
770 
771  /// Constructor, there are 3 internal values (for the pressure)
774  {
775  //Allocate and add one Internal data object that stored 2+1 pressure
776  //values;
778  }
779 
780  /// \short Number of values (pinned or dofs) required at local node n.
781  inline unsigned required_nvalue(const unsigned &n) const {return 3;}
782 
783 
784  /// \short Return the i-th pressure value
785  /// (Discontinous pressure interpolation -- no need to cater for hanging
786  /// nodes).
787  double p_spherical_nst(const unsigned &i) const
789 
790  /// Return number of pressure values
791  unsigned npres_spherical_nst() const {return 3;}
792 
793  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
794  void fix_pressure(const unsigned &p_dof, const double &p_value)
795  {
798  }
799 
800  /// \short Add to the set \c paired_load_data pairs containing
801  /// - the pointer to a Data object
802  /// and
803  /// - the index of the value in that Data object
804  /// .
805  /// for all values (pressures, velocities) that affect the
806  /// load computed in the \c get_load(...) function.
807  void identify_load_data(std::set<std::pair<Data*,unsigned> >
808  &paired_load_data);
809 
810  /// \short Add to the set \c paired_pressure_data pairs containing
811  /// - the pointer to a Data object
812  /// and
813  /// - the index of the value in that Data object
814  /// .
815  /// for pressure values that affect the
816  /// load computed in the \c get_load(...) function.
817  void identify_pressure_data(std::set<std::pair<Data*,unsigned> >
818  &paired_pressure_data);
819 
820  /// Redirect output to SphericalNavierStokesEquations output
821  void output(std::ostream &outfile)
823 
824  /// Redirect output to SphericalNavierStokesEquations output
825  void output(std::ostream &outfile, const unsigned &nplot)
827 
828 
829  /// Redirect output to SphericalNavierStokesEquations output
830  void output(FILE* file_pt)
832 
833  /// Redirect output to SphericalNavierStokesEquations output
834  void output(FILE* file_pt, const unsigned &nplot)
836 
837 
838  /// \short Full output function:
839  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
840  /// in tecplot format. Default number of plot points
841  void full_output(std::ostream &outfile)
843 
844  /// \short Full output function:
845  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
846  /// in tecplot format. nplot points in each coordinate direction
847  void full_output(std::ostream &outfile, const unsigned &nplot)
849 
850 
851  /// \short The number of "DOF types" that degrees of freedom in this element
852  /// are sub-divided into: Velocity (three comp) and pressure.
853  unsigned ndof_types() const
854  {
855  return 4;
856  }
857 
858  /// \short Create a list of pairs for all unknowns in this element,
859  /// so that the first entry in each pair contains the global equation
860  /// number of the unknown, while the second one contains the number
861  /// of the "DOF type" that this unknown is associated with.
862  /// (Function can obviously only be called if the equation numbering
863  /// scheme has been set up.)
865  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const;
866 
867 };
868 
869 //Inline functions
870 
871 //=======================================================================
872 /// Derivatives of the shape functions and test functions w.r.t. to global
873 /// (Eulerian) coordinates. Return Jacobian of mapping between
874 /// local and global coordinates.
875 //=======================================================================
878  const Vector<double> &s, Shape &psi,
879  DShape &dpsidx, Shape &test,
880  DShape &dtestdx) const
881 {
882  //Call the geometrical shape functions and derivatives
883  double J = this->dshape_eulerian(s,psi,dpsidx);
884  //Loop over the test functions and derivatives and set them equal to the
885  //shape functions
886  for(unsigned i=0;i<9;i++)
887  {
888  test[i] = psi[i];
889  dtestdx(i,0) = dpsidx(i,0);
890  dtestdx(i,1) = dpsidx(i,1);
891  }
892  //Return the jacobian
893  return J;
894 }
895 
896 //=======================================================================
897 /// Derivatives of the shape functions and test functions w.r.t. to global
898 /// (Eulerian) coordinates. Return Jacobian of mapping between
899 /// local and global coordinates.
900 //=======================================================================
903  const unsigned &ipt, Shape &psi,
904  DShape &dpsidx, Shape &test,
905  DShape &dtestdx) const
906 {
907  //Call the geometrical shape functions and derivatives
908  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
909  //Loop over the test functions and derivatives and set them equal to the
910  //shape functions
911  for(unsigned i=0;i<9;i++)
912  {
913  test[i] = psi[i];
914  dtestdx(i,0) = dpsidx(i,0);
915  dtestdx(i,1) = dpsidx(i,1);
916  }
917  //Return the jacobian
918  return J;
919 }
920 
921 
922 //=======================================================================
923 /// Pressure shape functions
924 //=======================================================================
927 const
928 {
929  psi[0] = 1.0;
930  psi[1] = s[0];
931  psi[2] = s[1];
932 }
933 
934 ///Define the pressure shape and test functions
936  const Vector<double> &s,
937  Shape &psi,
938  Shape &test) const
939 {
940  //Call the pressure shape functions
941  pshape_spherical_nst(s,psi);
942  //Loop over the test functions and set them equal to the shape functions
943  for(unsigned i=0;i<3;i++) test[i] = psi[i];
944 }
945 
946 //=======================================================================
947 /// Face geometry of the Spherical Crouzeix_Raviart elements
948 //=======================================================================
949 template<>
951 public virtual QElement<1,3>
952 {
953  public:
954  FaceGeometry() : QElement<1,3>() {}
955 };
956 
957 //=======================================================================
958 /// Face geometry of the FaceGeometry of the Spherical
959 /// Crouzeix_Raviart elements
960 //=======================================================================
961 template<>
963 public virtual PointElement
964 {
965  public:
967 };
968 
969 
970 ////////////////////////////////////////////////////////////////////////////
971 ////////////////////////////////////////////////////////////////////////////
972 
973 
974 //=======================================================================
975 /// Taylor--Hood elements are Navier--Stokes elements
976 /// with quadratic interpolation for velocities and positions and
977 /// continous linear pressure interpolation. They can be used
978 /// within oomph-lib's block-preconditioning framework.
979 //=======================================================================
980 class QSphericalTaylorHoodElement : public virtual QElement<2,3>,
981  public virtual SphericalNavierStokesEquations
982 {
983  private:
984 
985  /// Static array of ints to hold number of variables at node
986  static const unsigned Initial_Nvalue[];
987 
988  protected:
989 
990  /// \short Static array of ints to hold conversion from pressure
991  /// node numbers to actual node numbers
992  static const unsigned Pconv[];
993 
994  /// \short Velocity shape and test functions and their derivs
995  /// w.r.t. to global coords at local coordinate s (taken from geometry)
996  /// Return Jacobian of mapping between local and global coordinates.
998  Shape &psi,
999  DShape &dpsidx, Shape &test,
1000  DShape &dtestdx) const;
1001 
1002  /// \short Velocity shape and test functions and their derivs
1003  /// w.r.t. to global coords at local coordinate s (taken from geometry)
1004  /// Return Jacobian of mapping between local and global coordinates.
1005  inline double dshape_and_dtest_eulerian_at_knot_spherical_nst(const unsigned &ipt,
1006  Shape &psi,
1007  DShape &dpsidx,
1008  Shape &test,
1009  DShape &dtestdx) const;
1010 
1011  /// Pressure shape functions at local coordinate s
1012  inline void pshape_spherical_nst(const Vector<double> &s, Shape &psi) const;
1013 
1014  /// Pressure shape and test functions at local coordinte s
1015  inline void pshape_spherical_nst(const Vector<double> &s, Shape &psi,
1016  Shape &test) const;
1017 
1018  /// Return the local equation numbers for the pressure values.
1019  inline int p_local_eqn(const unsigned &n)const
1020  {return this->nodal_local_eqn(Pconv[n],p_nodal_index_spherical_nst());}
1021 
1022 public:
1023 
1024  /// Constructor, no internal data points
1027 
1028  /// \short Number of values (pinned or dofs) required at node n. Can
1029  /// be overwritten for hanging node version
1030  inline virtual unsigned required_nvalue(const unsigned &n) const
1031  {return Initial_Nvalue[n];}
1032 
1033  /// \short Set the value at which the pressure is stored in the nodes
1034  /// In this case the third index because there are three velocity components
1035  int p_nodal_index_spherical_nst() const {return 3;}
1036 
1037  /// \short Access function for the pressure values at local pressure
1038  /// node n_p (const version)
1039  double p_spherical_nst(const unsigned &n_p) const
1040  {return this->nodal_value(Pconv[n_p],this->p_nodal_index_spherical_nst());}
1041 
1042  /// Return number of pressure values
1043  unsigned npres_spherical_nst() const {return 4;}
1044 
1045  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
1046  void fix_pressure(const unsigned &p_dof, const double &p_value)
1047  {
1048  this->node_pt(Pconv[p_dof])->pin(this->p_nodal_index_spherical_nst());
1049  this->node_pt(Pconv[p_dof])->set_value(this->p_nodal_index_spherical_nst(),p_value);
1050  }
1051 
1052 
1053  /// \short Add to the set \c paired_load_data pairs containing
1054  /// - the pointer to a Data object
1055  /// and
1056  /// - the index of the value in that Data object
1057  /// .
1058  /// for all values (pressures, velocities) that affect the
1059  /// load computed in the \c get_load(...) function.
1060  void identify_load_data(std::set<std::pair<Data*,unsigned> > &
1061  paired_load_data);
1062 
1063 
1064  /// \short Add to the set \c paired_pressure_data pairs containing
1065  /// - the pointer to a Data object
1066  /// and
1067  /// - the index of the value in that Data object
1068  /// .
1069  /// for pressure values that affect the
1070  /// load computed in the \c get_load(...) function.
1072  std::set<std::pair<Data*,unsigned> > &paired_pressure_data);
1073 
1074  /// Redirect output to SphericalNavierStokesEquations output
1075  void output(std::ostream &outfile)
1077 
1078  /// Redirect output to SphericalNavierStokesEquations output
1079  void output(std::ostream &outfile, const unsigned &nplot)
1080  {SphericalNavierStokesEquations::output(outfile,nplot);}
1081 
1082  /// Redirect output to SphericalNavierStokesEquations output
1083  void output(FILE* file_pt) {SphericalNavierStokesEquations::output(file_pt);}
1084 
1085  /// Redirect output to SphericalNavierStokesEquations output
1086  void output(FILE* file_pt, const unsigned &nplot)
1087  {SphericalNavierStokesEquations::output(file_pt,nplot);}
1088 
1089 
1090  /// \short The number of "DOF types" that degrees of freedom in this element
1091  /// are sub-divided into: Velocity (3 components) and pressure.
1092  unsigned ndof_types() const
1093  {
1094  return 4;
1095  }
1096 
1097  /// \short Create a list of pairs for all unknowns in this element,
1098  /// so that the first entry in each pair contains the global equation
1099  /// number of the unknown, while the second one contains the number
1100  /// of the "Dof type" that this unknown is associated with.
1101  /// (Function can obviously only be called if the equation numbering
1102  /// scheme has been set up.)
1104  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list)const;
1105 
1106 
1107 };
1108 
1109 //Inline functions
1110 
1111 //==========================================================================
1112 /// 2D :
1113 /// Derivatives of the shape functions and test functions w.r.t to
1114 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1115 /// local and global coordinates.
1116 //==========================================================================
1117 inline double QSphericalTaylorHoodElement::
1119  const Vector<double> &s,
1120  Shape &psi,
1121  DShape &dpsidx, Shape &test,
1122  DShape &dtestdx) const
1123 {
1124  //Call the geometrical shape functions and derivatives
1125  double J = this->dshape_eulerian(s,psi,dpsidx);
1126  //Loop over the test functions and derivatives and set them equal to the
1127  //shape functions
1128  for(unsigned i=0;i<9;i++)
1129  {
1130  test[i] = psi[i];
1131  dtestdx(i,0) = dpsidx(i,0);
1132  dtestdx(i,1) = dpsidx(i,1);
1133  }
1134  //Return the jacobian
1135  return J;
1136 }
1137 
1138 
1139 //==========================================================================
1140 /// Derivatives of the shape functions and test functions w.r.t to
1141 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1142 /// local and global coordinates.
1143 //==========================================================================
1144 inline double QSphericalTaylorHoodElement::
1146  const unsigned &ipt,Shape &psi, DShape &dpsidx, Shape &test,
1147  DShape &dtestdx) const
1148 {
1149  //Call the geometrical shape functions and derivatives
1150  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
1151  //Loop over the test functions and derivatives and set them equal to the
1152  //shape functions
1153  for(unsigned i=0;i<9;i++)
1154  {
1155  test[i] = psi[i];
1156  dtestdx(i,0) = dpsidx(i,0);
1157  dtestdx(i,1) = dpsidx(i,1);
1158  }
1159  //Return the jacobian
1160  return J;
1161 }
1162 
1163 //==========================================================================
1164 /// Pressure shape functions
1165 //==========================================================================
1168  const
1169 {
1170  //Local storage
1171  double psi1[2], psi2[2];
1172  //Call the OneDimensional Shape functions
1173  OneDimLagrange::shape<2>(s[0],psi1);
1174  OneDimLagrange::shape<2>(s[1],psi2);
1175 
1176  //Now let's loop over the nodal points in the element
1177  //s1 is the "x" coordinate, s2 the "y"
1178  for(unsigned i=0;i<2;i++)
1179  {
1180  for(unsigned j=0;j<2;j++)
1181  {
1182  /*Multiply the two 1D functions together to get the 2D function*/
1183  psi[2*i + j] = psi2[i]*psi1[j];
1184  }
1185  }
1186 }
1187 
1188 
1189 //==========================================================================
1190 /// Pressure shape and test functions
1191 //==========================================================================
1193  Shape &psi,
1194  Shape &test) const
1195 {
1196  //Call the pressure shape functions
1197  pshape_spherical_nst(s,psi);
1198  //Loop over the test functions and set them equal to the shape functions
1199  for(unsigned i=0;i<4;i++) test[i] = psi[i];
1200 }
1201 
1202 
1203 //=======================================================================
1204 /// Face geometry of the Spherical Taylor_Hood elements
1205 //=======================================================================
1206 template<>
1208 {
1209  public:
1210  FaceGeometry() : QElement<1,3>() {}
1211 };
1212 
1213 //=======================================================================
1214 /// Face geometry of the FaceGeometry of the 2D Taylor Hoodelements
1215 //=======================================================================
1216 template<>
1218 public virtual PointElement
1219 {
1220  public:
1222 };
1223 
1224 
1225 }
1226 
1227 #endif
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
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
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
SphericalNavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocity (3 c...
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
void output_vorticity(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z], [omega_x,omega_y,[and/or omega_z]] in tecplot format. nplot points in each ...
virtual double p_spherical_nst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes...
virtual void get_body_force_spherical_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, Vector< double > &result)
Calculate the body force at a given time and local and/or Eulerian position. This function is virtual...
double *& re_invfr_pt()
Pointer to global inverse Froude number.
double pressure_integral() const
Integral of pressure over element.
void fix_pressure(const unsigned &p_dof, const double &p_value)
Pin p_dof-th pressure dof and set it to value specified by p_value.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element's residual Vector.
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers...
void output(FILE *file_pt)
C-style output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points...
virtual double get_source_spherical_nst(double time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
SphericalNavierStokesSourceFctPt & source_fct_pt()
Access function for the source-function pointer.
void output(FILE *file_pt)
Redirect output to SphericalNavierStokesEquations output.
const double & re() const
Reynolds number.
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
double u_spherical_nst(const unsigned &t, const unsigned &n, const unsigned &i) const
Velocity i at local node n at timestep t (t=0: present; t>0: previous). Uses suitably interpolated va...
void output(std::ostream &outfile, const unsigned &nplot)
Redirect output to SphericalNavierStokesEquations output.
cstr elem_len * i
Definition: cfortran.h:607
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
double dissipation() const
Return integral of dissipation over element.
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
void full_output(std::ostream &outfile)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format...
void get_load(const Vector< double > &s, const Vector< double > &N, Vector< double > &load)
This implements a pure virtual function defined in the FSIFluidElement class. The function computes t...
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
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
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
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 & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
std::string type() const
Return string that indicates the type of the timestepper (e.g. "BDF", "Newmark", etc.)
Definition: timesteppers.h:465
double interpolated_dudx_spherical_nst(const Vector< double > &s, const unsigned &i, const unsigned &j) const
Return matrix entry dudx(i,j) of the FE interpolated velocity derivative at local coordinate s...
const Vector< double > & g() const
Vector of gravitational components.
unsigned npres_spherical_nst() const
Return number of pressure values.
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
virtual double dshape_and_dtest_eulerian_at_knot_spherical_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and derivatives w.r.t. global coords at ipt-th integration point Return J...
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Validate against exact solution at given time Solution is provided via function pointer. Plot at a given number of plot points and compute L2 error and L2 norm of velocity solution over element.
unsigned npres_spherical_nst() const
Return number of pressure values.
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
QSphericalTaylorHoodElement()
Constructor, no internal data points.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
double *& density_ratio_pt()
Pointer to Density ratio.
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)
Compute traction (on the viscous scale) exerted onto the fluid at local coordinate s...
void compute_error_e(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, FiniteElement::SteadyExactSolutionFctPt exact_soln_dr_pt, FiniteElement::SteadyExactSolutionFctPt exact_soln_dtheta_pt, double &u_error, double &u_norm, double &p_error, double &p_norm)
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
virtual double dshape_and_dtest_eulerian_spherical_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and derivatives w.r.t. global coords at local coordinate s...
virtual unsigned npres_spherical_nst() const =0
Function to return number of pressure degrees of freedom.
Vector< double > actual_dth(const Vector< double > &x)
Vector< double > actual_dr(const Vector< double > &x)
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
jacobian matrix and mass matrix
Vector< double > * G_pt
Pointer to global gravity Vector.
double(* SphericalNavierStokesSourceFctPt)(const double &time, const Vector< double > &x)
Function pointer to source function fct(t,x) x is a Vector!
QSphericalCrouzeixRaviartElement()
Constructor, there are 3 internal values (for the pressure)
void output(std::ostream &outfile)
Redirect output to SphericalNavierStokesEquations output.
double d_kin_energy_dt() const
Get integral of time derivative of kinetic energy over element.
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 char t char * s
Definition: cfortran.h:572
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
void full_output(std::ostream &outfile, const unsigned &nplot)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format...
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
virtual void fill_in_generic_residual_contribution_spherical_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...
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
double * ReInvRo_pt
Pointer to global Reynolds number x inverse Rossby number (used when in a rotating frame) ...
const double & re_invfr() const
Global inverse Froude number.
int p_local_eqn(const unsigned &n) const
Return the local equation numbers for the pressure values.
SphericalNavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed...
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
SphericalNavierStokesEquations()
Constructor: NULL the body force and source function and make sure the ALE terms are included by defa...
void(* SphericalNavierStokesBodyForceFctPt)(const double &time, const Vector< double > &x, Vector< double > &body_force)
Function pointer to body force function fct(t,x,f(x)) x is a Vector!
double dshape_and_dtest_eulerian_spherical_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
void output_veloc(std::ostream &outfile, const unsigned &nplot, const unsigned &t)
Output function: x,y,[z],u,v,[w] in tecplot format. nplot points in each coordinate direction at time...
double p_spherical_nst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version) ...
virtual void pshape_spherical_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
void output(FILE *file_pt, const unsigned &nplot)
Redirect output to SphericalNavierStokesEquations output.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double *& re_invro_pt()
Pointer to global inverse Froude number.
double p_spherical_nst(const unsigned &i) const
Return the i-th pressure value (Discontinous pressure interpolation – no need to cater for hanging no...
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number) ...
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
void output(std::ostream &outfile, const unsigned &nplot)
Redirect output to SphericalNavierStokesEquations output.
virtual unsigned u_index_spherical_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
double dshape_and_dtest_eulerian_spherical_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
SphericalNavierStokesBodyForceFctPt body_force_fct_pt() const
Access function for the body-force pointer. Const version.
virtual void fix_pressure(const unsigned &p_dof, const double &p_value)=0
Pin p_dof-th pressure dof and set it to value specified by p_value.
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number) ...
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
double dshape_and_dtest_eulerian_at_knot_spherical_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at ipt-th integation point...
double * Re_pt
Pointer to global Reynolds number.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
void get_vorticity(const Vector< double > &s, Vector< double > &vorticity) const
Compute the vorticity vector at local coordinate s.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as ...
Definition: elements.h:1726
void output(FILE *file_pt, const unsigned &nplot)
Redirect output to SphericalNavierStokesEquations output.
virtual int p_local_eqn(const unsigned &n) const =0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
double dshape_and_dtest_eulerian_at_knot_spherical_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
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
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the elemental contribution to the jacobian matrix. and the residuals vector. Note that this funct...
int p_local_eqn(const unsigned &n) const
Return the local equation numbers for the pressure values.
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.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
SphericalNavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
double *& re_pt()
Pointer to Reynolds number.
virtual int p_nodal_index_spherical_nst() const
Return the index at which the pressure is stored if it is stored at the nodes. If not stored at the n...
void fix_pressure(const unsigned &p_dof, const double &p_value)
Pin p_dof-th pressure dof and set it to value specified by p_value.
void pshape_spherical_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
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
SphericalNavierStokesSourceFctPt source_fct_pt() const
Access function for the source-function pointer. Const version.
double interpolated_u_spherical_nst(const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s.
void output(std::ostream &outfile)
Redirect output to SphericalNavierStokesEquations output.
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
double u_spherical_nst(const unsigned &n, const unsigned &i) const
Velocity i at local node n. Uses suitably interpolated value for hanging nodes. The use of u_index_sp...
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
void output(FILE *file_pt)
Redirect output to SphericalNavierStokesEquations output.
void pshape_spherical_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
double kin_energy() const
Get integral of kinetic energy over element.
int p_nodal_index_spherical_nst() const
Set the value at which the pressure is stored in the nodes In this case the third index because there...
void interpolated_u_spherical_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
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...
Vector< double > actual(const Vector< double > &x)
The FSIFluidElement class is a base class for all fluid finite elements that apply a load (traction) ...
Definition: fsi.h:67
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
double interpolated_p_spherical_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution specified via function pointer at a given number of plot points. Function prints as many components as are returned in solution Vector.
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
double du_dt_spherical_nst(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes...
void full_output(std::ostream &outfile)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format...
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocity (thr...
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:596
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)
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