polar_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 // Created (18/03/08) by recopying from my_oomph-codes/polar_navier_stokes.h
31 // Now I know what I'm doing it needs updating to current oomph conventions.
32 
33 // 17/06/08 - Weakform adjusted to "correct" traction version
34 
35 // "_pnst" added to: - npres()
36 // - pshape()
37 // - u()
38 // - p()
39 // - du_dt()
40 // - interpolated_u()
41 // - interpolated_p()
42 // - interpolated_dudx()
43 // - dshape_and_dtest_eulerian()
44 // - dshape_and_dtest_eulerian_at_knot()
45 
46 // The following generality is necessary for elements with multiple physics.
47 // That is where we can't just assume that values one and two at the nodes
48 // are velocities and the third a pressure.
49 
50 // Then renamed "index_of_nodal_pressure_value to p_nodal_index_pnst.
51 
52 // Added u_index_pnst (by default just returns the value you give it)
53 
54 // The internal pressure data in Crouzeix Raviart elements is now stored
55 // as one data object with three values:
56 // Added P_pnst_internal_index which stores the index of that internal
57 // data, specified in the constructor.
58 // Adjusted: - p_pnst
59 // - fix_pressure
60 // - get_load_data
61 // To reflect this change in internal data storage.
62 
63 // Plus added Pressure_not_stored_at_node, default = -100.
64 
65 // assign_additional_local_eqn_numbers removed in favour of
66 // - p_local_eqn
67 // And we have:
68 // - local_eqn = nodal_local_eqn(l,u_nodal_index[i]);
69 // - local_eqn = p_local_eqn(l);
70 // in fill_in_generic_residual_contribution
71 
72 // Also removed the following functions for pinning/unpinning pressures
73 // as they're only needed in the refineable case:
74 //----------------------------------------------------------------------
75 // - unpin_all_nodal_pressure_dofs()
76 // - unpin_all_internal_pressure_dofs()
77 // - pin_all_nodal_pressure_dofs()
78 // - unpin_proper_nodal_pressure_dofs()
79 // - pin_redundant_nodal_pressures()
80 // - unpin_all_pressure_dofs()
81 // - pressure_dof_is_hanging()
82 //----------------------------------------------------------------------
83 
84 // strain_rate adapted to return the polar strain components.
85 
86 // interpolated positions / velocities now assembled using nodal_position
87 // and nodal_value - both of which should cope with hanging nodes.
88 
89 // Also stokes_output() added to compute convergence rates when exact
90 // (Stokes) solution is known.
91 
92 #ifndef OOMPH_POLAR_NAVIER_STOKES
93 #define OOMPH_POLAR_NAVIER_STOKES
94 
95 // Config header generated by autoconfig
96 #ifdef HAVE_CONFIG_H
97 #include <oomph-lib-config.h>
98 #endif
99 
100 //OOMPH-LIB headers
101 #include "../generic/Qelements.h"
102 
103 
104 namespace oomph
105 {
106 
107 //=====================================================================
108 /// A class for elements that solve the polar Navier--Stokes equations,
109 /// This contains the generic maths -- any concrete implementation must
110 /// be derived from this.
111 //======================================================================
113 {
114 
115  public:
116 
117  /// \short Function pointer to body force function fct(t,x,f(x))
118  /// x is a Vector!
119  typedef void (*NavierStokesBodyForceFctPt)(const double& time,
120  const Vector<double>& x,
121  Vector<double>& body_force);
122 
123  /// \short Function pointer to source function fct(t,x)
124  /// x is a Vector!
125  typedef double (*NavierStokesSourceFctPt)(const double& time,
126  const Vector<double>& x);
127 
128  private:
129 
130  /// \short Static "magic" number that indicates that the pressure is
131  /// not stored at a node
133 
134  /// Static default value for the physical constants (all initialised to zero)
136 
137  /// Static default value for the physical ratios (all are initialised to one)
139 
140  /// Static default value for the gravity vector
142 
143  protected:
144 
145  //Physical constants
146 
147  /// \short Pointer to the viscosity ratio (relative to the
148  /// viscosity used in the definition of the Reynolds number)
150 
151  /// \short Pointer to the density ratio (relative to the density used in the
152  /// definition of the Reynolds number)
154 
155  // Pointers to global physical constants
156 
157  /// Pointer to the angle alpha
158  double *Alpha_pt;
159 
160  /// Pointer to global Reynolds number
161  double *Re_pt;
162 
163  /// Pointer to global Reynolds number x Strouhal number (=Womersley)
164  double *ReSt_pt;
165 
166  /// \short Pointer to global Reynolds number x inverse Froude number
167  /// (= Bond number / Capillary number)
168  double *ReInvFr_pt;
169 
170  /// Pointer to global gravity Vector
172 
173  /// Pointer to body force function
175 
176  /// Pointer to volumetric source function
178 
179  /// \short Access function for the local equation number information for
180  /// the pressure.
181  /// p_local_eqn[n] = local equation number or < 0 if pinned
182  virtual int p_local_eqn(const unsigned &n)=0;
183 
184  /// \short Compute the shape functions and derivatives
185  /// w.r.t. global coords at local coordinate s.
186  /// Return Jacobian of mapping between local and global coordinates.
187  virtual double dshape_and_dtest_eulerian_pnst(const Vector<double> &s, Shape &psi,
188  DShape &dpsidx, Shape &test,
189  DShape &dtestdx) const=0;
190 
191  /// \short Compute the shape functions and derivatives
192  /// w.r.t. global coords at ipt-th integration point
193  /// Return Jacobian of mapping between local and global coordinates.
194  virtual double dshape_and_dtest_eulerian_at_knot_pnst(const unsigned &ipt,
195  Shape &psi,
196  DShape &dpsidx, Shape &test,
197  DShape &dtestdx) const=0;
198 
199  /// Compute the pressure shape functions at local coordinate s
200  virtual void pshape_pnst(const Vector<double> &s, Shape &psi) const=0;
201 
202  /// \short Compute the pressure shape and test functions
203  /// at local coordinate s
204  virtual void pshape_pnst(const Vector<double> &s, Shape &psi, Shape &test) const=0;
205 
206 
207  /// Calculate the body force at a given time and Eulerian position
208  void get_body_force(double time, const Vector<double> &x,
209  Vector<double> &result)
210  {
211  //If the function pointer is zero return zero
212  if(Body_force_fct_pt == 0)
213  {
214  //Loop over dimensions and set body forces to zero
215  for(unsigned i=0;i<2;i++) {result[i] = 0.0;}
216  }
217  //Otherwise call the function
218  else
219  {
220  (*Body_force_fct_pt)(time,x,result);
221  }
222  }
223 
224  /// \short Calculate the source fct at given time and
225  /// Eulerian position
226  double get_source_fct(double time, const Vector<double> &x)
227  {
228  //If the function pointer is zero return zero
229  if (Source_fct_pt == 0) {return 0;}
230  //Otherwise call the function
231  else {return (*Source_fct_pt)(time,x);}
232  }
233 
234 public:
235 
236  /// \short Constructor: NULL the body force and source function
238  {
239  //Set all the Physical parameter pointers to the default value zero
244  //Set the Physical ratios to the default value of 1
247  }
248 
249  /// Vector to decide whether the stress-divergence form is used or not
250  // N.B. This needs to be public so that the intel compiler gets things correct
251  // somehow the access function messes things up when going to refineable
252  // navier--stokes
254 
255  //Access functions for the physical constants
256 
257  /// Reynolds number
258  const double &re() const {return *Re_pt;}
259 
260  /// Alpha
261  const double &alpha() const {return *Alpha_pt;}
262 
263  /// Product of Reynolds and Strouhal number (=Womersley number)
264  const double &re_st() const {return *ReSt_pt;}
265 
266  /// Pointer to Reynolds number
267  double* &re_pt() {return Re_pt;}
268 
269  /// Pointer to Alpha
270  double* &alpha_pt() {return Alpha_pt;}
271 
272  /// Pointer to product of Reynolds and Strouhal number (=Womersley number)
273  double* &re_st_pt() {return ReSt_pt;}
274 
275  /// \short Viscosity ratio for element: Element's viscosity relative to the
276  /// viscosity used in the definition of the Reynolds number
277  const double &viscosity_ratio() const {return *Viscosity_Ratio_pt;}
278 
279  /// Pointer to Viscosity Ratio
281 
282  /// \short Density ratio for element: Element's density relative to the
283  /// viscosity used in the definition of the Reynolds number
284  const double &density_ratio() const {return *Density_Ratio_pt;}
285 
286  /// Pointer to Density ratio
287  double* &density_ratio_pt() {return Density_Ratio_pt;}
288 
289  /// Global inverse Froude number
290  const double &re_invfr() const {return *ReInvFr_pt;}
291 
292  /// Pointer to global inverse Froude number
293  double* &re_invfr_pt() {return ReInvFr_pt;}
294 
295  /// Vector of gravitational components
296  const Vector<double> &g() const {return *G_pt;}
297 
298  /// Pointer to Vector of gravitational components
299  Vector<double>* &g_pt() {return G_pt;}
300 
301  /// Access function for the body-force pointer
303  {return Body_force_fct_pt;}
304 
305  /// Access function for the body-force pointer. Const version
307  {return Body_force_fct_pt;}
308 
309  ///Access function for the source-function pointer
311 
312  ///Access function for the source-function pointer. Const version
314 
315  ///Function to return number of pressure degrees of freedom
316  virtual unsigned npres_pnst() const=0;
317 
318  /// \short Velocity i at local node n. Uses suitably interpolated value
319  /// for hanging nodes.
320  virtual double u_pnst(const unsigned &n, const unsigned &i) const=0;
321 
322  /// \short Velocity i at local node n at timestep t (t=0: present;
323  /// t>0: previous). Uses suitably interpolated value for hanging nodes.
324  virtual double u_pnst(const unsigned &t, const unsigned &n,
325  const unsigned &i) const=0;
326 
327  /// \short Return the index at which the i-th unknown velocity component
328  /// is stored. The default value, i, is appropriate for single-physics
329  /// problems.
330  /// In derived multi-physics elements, this function should be overloaded
331  /// to reflect the chosen storage scheme. Note that these equations require
332  /// that the unknowns are always stored at the same indices at each node.
333  virtual inline unsigned u_index_pnst(const unsigned &i) const {return i;}
334 
335  /// \short i-th component of du/dt at local node n.
336  /// Uses suitably interpolated value for hanging nodes.
337  double du_dt_pnst(const unsigned &n, const unsigned &i) const
338  {
339  // Get the data's timestepper
341 
342  // Number of timsteps (past & present)
343  unsigned n_time = time_stepper_pt->ntstorage();
344 
345  //Initialise dudt
346  double dudt=0.0;
347  //Loop over the timesteps, if there is a non Steady timestepper
348  if (time_stepper_pt->type()!="Steady")
349  {
350  for(unsigned t=0;t<n_time;t++)
351  {
352  dudt+=time_stepper_pt->weight(1,t)*u_pnst(t,n,i);
353  }
354  }
355 
356  return dudt;
357  }
358 
359  /// \short Pressure at local pressure "node" n_p
360  /// Uses suitably interpolated value for hanging nodes.
361  virtual double p_pnst(const unsigned &n_p)const=0;
362 
363  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
364  virtual void fix_pressure(const unsigned &p_dof, const double &p_value)=0;
365 
366  /// \short Which nodal value represents the pressure? (Default: negative,
367  /// indicating that pressure is not based on nodal interpolation).
369 
370  /// \short Pointer to n_p-th pressure node (Default: NULL,
371  /// indicating that pressure is not based on nodal interpolation).
372  virtual Node* pressure_node_pt(const unsigned& n_p) {return NULL;}
373 
374  /// Integral of pressure over element
375  double pressure_integral() const;
376 
377  /// \short Return integral of dissipation over element
378  double dissipation() const;
379 
380  /// \short Return dissipation at local coordinate s
381  double dissipation(const Vector<double>& s) const;
382 
383 /// \short Get integral of kinetic energy over element
384  double kin_energy() const;
385 
386  /// Strain-rate tensor: Now returns polar strain
387  void strain_rate(const Vector<double>& s,
389 
390  /// Function to return polar strain multiplied by r
391  void strain_rate_by_r(const Vector<double>& s,
393 
394  /// \short Compute traction (on the viscous scale) at local coordinate s
395  /// for outer unit normal N
396  void get_traction(const Vector<double>& s, const Vector<double>& N,
397  Vector<double>& traction);
398 
399  /// \short The potential loading on an external SolidElement is always
400  /// provided by the traction function
401  void get_load(const Vector<double> &s, const Vector<double> &xi,
402  const Vector<double> &x, const Vector<double> &N,
403  Vector<double> &load)
404  {get_traction(s,N,load);}
405 
406 
407  /// \short Output functionget_vels(const Vector<double>& x_to_get, Vector<double>& vels): x,y,[z],u,v,[w],p
408  /// in tecplot format. Default number of plot points
409  void output(std::ostream &outfile)
410  {
411  unsigned nplot=5;
412  output(outfile,nplot);
413  }
414 
415  /// \short Output function: x,y,[z],u,v,[w],p
416  /// in tecplot format. nplot points in each coordinate direction
417  void output(std::ostream &outfile, const unsigned &nplot);
418 
419  /// \short C-style output function: x,y,[z],u,v,[w],p
420  /// in tecplot format. Default number of plot points
421  void output(FILE* file_pt)
422  {
423  unsigned nplot=5;
424  output(file_pt,nplot);
425  }
426 
427  /// \short C-style output function: x,y,[z],u,v,[w],p
428  /// in tecplot format. nplot points in each coordinate direction
429  void output(FILE* file_pt, const unsigned &nplot);
430 
431  /// \short Full output function:
432  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
433  /// in tecplot format. Default number of plot points
434  void full_output(std::ostream &outfile)
435  {
436  unsigned nplot=5;
437  full_output(outfile,nplot);
438  }
439 
440  /// \short Full output function:
441  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
442  /// in tecplot format. nplot points in each coordinate direction
443  void full_output(std::ostream &outfile, const unsigned &nplot);
444 
445 
446  /// \short Output function: x,y,[z],u,v,[w] in tecplot format.
447  /// nplot points in each coordinate direction at timestep t
448  /// (t=0: present; t>0: previous timestep)
449  void output_veloc(std::ostream &outfile, const unsigned &nplot,
450  const unsigned& t);
451 
452  /// \short Output exact solution specified via function pointer
453  /// at a given number of plot points. Function prints as
454  /// many components as are returned in solution Vector
455  void output_fct(std::ostream &outfile, const unsigned &nplot,
457 
458  /// \short Output exact solution specified via function pointer
459  /// at a given time and at a given number of plot points.
460  /// Function prints as many components as are returned in solution Vector.
461  void output_fct(std::ostream &outfile, const unsigned &nplot,
462  const double& time,
464 
465  /// \short Validate against exact solution at given time
466  /// Solution is provided via function pointer.
467  /// Plot at a given number of plot points and compute L2 error
468  /// and L2 norm of velocity solution over element
469  void compute_error(std::ostream &outfile,
471  const double& time,
472  double& error, double& norm);
473 
474  /// \short Validate against exact solution.
475  /// Solution is provided via function pointer.
476  /// Plot at a given number of plot points and compute L2 error
477  /// and L2 norm of velocity solution over element
478  void compute_error(std::ostream &outfile,
480  double& error, double& norm);
481 
482  /// Compute the element's residual Vector
484  {
485  //Call the generic residuals function with flag set to 0
489  }
490 
491  ///\short Compute the element's residual Vector and the jacobian matrix
492  /// Virtual function can be overloaded by hanging-node version
494  DenseMatrix<double> &jacobian)
495  {
496  //Call the generic routine with the flag set to 1
497  fill_in_generic_residual_contribution(residuals,jacobian,
499 
500  /*
501  // Only needed for test purposes
502  //Create a new matrix
503  unsigned n_dof = ndof();
504  DenseMatrix<double> jacobian_fd(n_dof,n_dof,0.0);
505  fill_in_jacobian_from_internal_by_fd(residuals,jacobian_fd);
506  fill_in_jacobian_from_nodal_by_fd(residuals,jacobian_fd);
507 
508  if(n_dof==21)
509  {
510  for(unsigned i=0;i<n_dof;i++)
511  {
512  for(unsigned j=0;j<n_dof;j++)
513  {
514  if(std::fabs(jacobian(i,j) - jacobian_fd(i,j)) > 1.0e-3)
515  {
516  std::cout << i << " " << j << " " << std::fabs(jacobian(i,j) - jacobian_fd(i,j)) << std::endl;
517  }
518  }
519  }
520  }
521  // Only needed for test purposes
522 
523  // Only needed for test purposes
524  //Create a new matrix
525  unsigned n_dof = ndof();
526  DenseMatrix<double> jacobian_new(n_dof,n_dof,0.0),jacobian_old(n_dof,n_dof,0.0);
527  Vector<double> residuals_new(n_dof,0.0),residuals_old(n_dof,0.0);
528 
529  //Call the generic routine with the flag set to 1
530  PolarNavierStokesEquations::fill_in_generic_residual_contribution(residuals_old,jacobian_old,
531  GeneralisedElement::Dummy_matrix,1);
532  //Call the generic routine with the flag set to 1
533  fill_in_generic_residual_contribution(residuals_new,jacobian_new,
534  GeneralisedElement::Dummy_matrix,1);
535 
536  if(n_dof==21)
537  {
538  for(unsigned i=0;i<n_dof;i++)
539  {
540  if(std::fabs(residuals_new[i] - residuals_old[i])>1.0e-8) std::cout << i << " " << std::fabs(residuals_new[i] - residuals_old[i]) << std::endl;
541  //std::cout << residuals_new[i] << std::endl;
542  for(unsigned j=0;j<n_dof;j++)
543  {
544  if(std::fabs(jacobian_new(i,j) - jacobian_old(i,j)) > 1.0e-8)
545  {
546  std::cout << i << " " << j << " " << std::fabs(jacobian_new(i,j) - jacobian_old(i,j)) << std::endl;
547  }
548  }
549  }
550  }
551  // Only needed for test purposes
552  */
553  } //End of fill_in_contribution_to_jacobian
554 
555  ///\short Compute the element's residual Vector and the jacobian matrix
556  /// Plus the mass matrix especially for eigenvalue problems
558  Vector<double> &residuals,
559  DenseMatrix<double> &jacobian,DenseMatrix<double> &mass_matrix)
560  {
561  //Call the generic routine with the flag set to 2
562  fill_in_generic_residual_contribution(residuals,jacobian,mass_matrix,2);
563  }
564 
565  ///\short Compute the residuals for the Navier--Stokes equations;
566  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
568  Vector<double> &residuals, DenseMatrix<double> &jacobian,
569  DenseMatrix<double> &mass_matrix, unsigned flag);
570 
571  /// Compute vector of FE interpolated velocity u at local coordinate s
573  {
574  //Find number of nodes
575  unsigned n_node = nnode();
576  //Local shape function
577  Shape psi(n_node);
578  //Find values of shape function
579  shape(s,psi);
580 
581  for (unsigned i=0;i<2;i++)
582  {
583  //Initialise value of u
584  veloc[i] = 0.0;
585  //Loop over the local nodes and sum
586  for(unsigned l=0;l<n_node;l++)
587  {
588  veloc[i] += u_pnst(l,i)*psi[l];
589  }
590  }
591  }
592 
593  /// Return FE interpolated velocity u[i] at local coordinate s
594  double interpolated_u_pnst(const Vector<double> &s, const unsigned &i) const
595  {
596  //Find number of nodes
597  unsigned n_node = nnode();
598  //Local shape function
599  Shape psi(n_node);
600  //Find values of shape function
601  shape(s,psi);
602 
603  //Initialise value of u
604  double interpolated_u = 0.0;
605  //Loop over the local nodes and sum
606  for(unsigned l=0;l<n_node;l++)
607  {
608  interpolated_u += u_pnst(l,i)*psi[l];
609  }
610 
611  return(interpolated_u);
612  }
613 
614  /// Return FE interpolated pressure at local coordinate s
615  double interpolated_p_pnst(const Vector<double> &s) const
616  {
617  //Find number of nodes
618  unsigned n_pres = npres_pnst();
619  //Local shape function
620  Shape psi(n_pres);
621  //Find values of shape function
622  pshape_pnst(s,psi);
623 
624  //Initialise value of p
625  double interpolated_p = 0.0;
626  //Loop over the local nodes and sum
627  for(unsigned l=0;l<n_pres;l++)
628  {
629  interpolated_p += p_pnst(l)*psi[l];
630  }
631 
632  return(interpolated_p);
633  }
634 
635  ///////////////////////////////////////////////////////////////////
636  /// My own work: ///
637  /// Return FE interpolated velocity derivative du[i]/dx[j] ///
638  /// at local coordinate s ///
639  ///////////////////////////////////////////////////////////////////
640  double interpolated_dudx_pnst(const Vector<double> &s, const unsigned &i,const unsigned &j) const
641  {
642  //Find number of nodes
643  unsigned n_node = nnode();
644 
645  //Set up memory for the shape and test functions
646  Shape psif(n_node);
647  DShape dpsifdx(n_node,2);
648 
649  //double J =
650  dshape_eulerian(s,psif,dpsifdx);
651 
652  //Initialise to zero
653  double interpolated_dudx = 0.0;
654 
655  //Calculate velocity derivative:
656 
657  // Loop over nodes
658  for(unsigned l=0;l<n_node;l++)
659  {
660  interpolated_dudx += u_pnst(l,i)*dpsifdx(l,j);
661  }
662 
663  return(interpolated_dudx);
664  }
665 
666 };
667 
668 //////////////////////////////////////////////////////////////////////////////
669 //////////////////////////////////////////////////////////////////////////////
670 //////////////////////////////////////////////////////////////////////////////
671 
672 
673 //==========================================================================
674 ///Crouzeix_Raviart elements are Navier--Stokes elements with quadratic
675 ///interpolation for velocities and positions, but a discontinuous linear
676 ///pressure interpolation
677 //==========================================================================
678 class PolarCrouzeixRaviartElement : public virtual QElement<2,3>,
679  public virtual PolarNavierStokesEquations
680 {
681  private:
682 
683  /// Static array of ints to hold required number of variables at nodes
684  static const unsigned Initial_Nvalue[];
685 
686  protected:
687 
688  /// Internal index that indicates at which internal data the pressure
689  /// is stored
691 
692  /// \short Velocity shape and test functions and their derivs
693  /// w.r.t. to global coords at local coordinate s (taken from geometry)
694  ///Return Jacobian of mapping between local and global coordinates.
695  inline double dshape_and_dtest_eulerian_pnst(const Vector<double> &s, Shape &psi,
696  DShape &dpsidx, Shape &test,
697  DShape &dtestdx) const;
698 
699  /// \short Velocity shape and test functions and their derivs
700  /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
701  ///Return Jacobian of mapping between local and global coordinates.
702  inline double dshape_and_dtest_eulerian_at_knot_pnst(const unsigned &ipt,
703  Shape &psi,
704  DShape &dpsidx, Shape &test,
705  DShape &dtestdx) const;
706 
707  /// Pressure shape functions at local coordinate s
708  inline void pshape_pnst(const Vector<double> &s, Shape &psi) const;
709 
710  /// Pressure shape and test functions at local coordinte s
711  inline void pshape_pnst(const Vector<double> &s, Shape &psi, Shape &test) const;
712 
713  /// Return the local equation numbers for the pressure values.
714  inline int p_local_eqn(const unsigned &n)
715  {return this->internal_local_eqn(P_pnst_internal_index,n);}
716 
717 public:
718 
719  /// Constructor, there are DIM+1 internal values (for the pressure)
721  {
722  //Allocate and add one Internal data object that stores 3 pressure
723  //values;
725  }
726 
727  /// \short Number of values (pinned or dofs) required at local node n.
728  virtual unsigned required_nvalue(const unsigned &n) const;
729 
730  /// \short Return the velocity component u[i] at local node
731  /// n. Uses suitably interpolated value for hanging nodes.
732  double u_pnst(const unsigned &n, const unsigned &i) const
733  {return this->nodal_value(n,i);}
734 
735  /// \short Return the velocity component u[i] at local node
736  /// n at timestep t (t=0: present; t>0: previous timestep).
737  /// Uses suitably interpolated value for hanging nodes.
738  double u_pnst(const unsigned &t, const unsigned &n, const unsigned &i)
739  const
740  {return this->nodal_value(t,n,i);}
741 
742  /// \short Return the pressure values at internal dof i_internal
743  /// (Discontinous pressure interpolation -- no need to cater for hanging
744  /// nodes).
745  double p_pnst(const unsigned &i_internal) const
746  {return *this->internal_data_pt(P_pnst_internal_index)->value_pt(i_internal);}
747 
748  /// Return number of pressure values
749  unsigned npres_pnst() const {return 3;}
750 
751  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
752  void fix_pressure(const unsigned &p_dof, const double &p_value)
753  {
755  *this->internal_data_pt(P_pnst_internal_index)->value_pt(p_dof) = p_value;
756  }
757 
758  /// \short Add to the set paired_load_data
759  /// pairs of pointers to data objects and unsigned integers that
760  /// index the values in the data object that affect the load (traction),
761  /// as specified in the get_load() function.
762  void get_load_data(std::set<std::pair<Data*,unsigned> > &paired_load_data);
763 
764  /// Redirect output to NavierStokesEquations output
765  void output(std::ostream &outfile) {PolarNavierStokesEquations::output(outfile);}
766 
767  /// Redirect output to NavierStokesEquations output
768  void output(std::ostream &outfile, const unsigned &Nplot)
769  {PolarNavierStokesEquations::output(outfile,Nplot);}
770 
771 
772  /// Redirect output to NavierStokesEquations output
773  void output(FILE* file_pt) {PolarNavierStokesEquations::output(file_pt);}
774 
775  /// Redirect output to NavierStokesEquations output
776  void output(FILE* file_pt, const unsigned &Nplot)
777  {PolarNavierStokesEquations::output(file_pt,Nplot);}
778 
779 
780  /// \short Full output function:
781  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
782  /// in tecplot format. Default number of plot points
783  void full_output(std::ostream &outfile)
785 
786  /// \short Full output function:
787  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
788  /// in tecplot format. nplot points in each coordinate direction
789  void full_output(std::ostream &outfile, const unsigned &nplot)
791 
792 };
793 
794 //Inline functions
795 
796 //=======================================================================
797 /// 2D
798 /// Derivatives of the shape functions and test functions w.r.t. to global
799 /// (Eulerian) coordinates. Return Jacobian of mapping between
800 /// local and global coordinates.
801 //=======================================================================
803  const Vector<double> &s, Shape &psi,
804  DShape &dpsidx, Shape &test,
805  DShape &dtestdx) const
806 {
807  //Call the geometrical shape functions and derivatives
808  double J = QElement<2,3>::dshape_eulerian(s,psi,dpsidx);
809  //Loop over the test functions and derivatives and set them equal to the
810  //shape functions
811  for(unsigned i=0;i<9;i++)
812  {
813  test[i] = psi[i];
814  dtestdx(i,0) = dpsidx(i,0);
815  dtestdx(i,1) = dpsidx(i,1);
816  }
817  //Return the jacobian
818  return J;
819 }
820 
821 
822 //=======================================================================
823 /// 2D
824 /// Derivatives of the shape functions and test functions w.r.t. to global
825 /// (Eulerian) coordinates. Return Jacobian of mapping between
826 /// local and global coordinates.
827 //=======================================================================
829  const unsigned &ipt, Shape &psi,
830  DShape &dpsidx, Shape &test,
831  DShape &dtestdx) const
832 {
833  //Call the geometrical shape functions and derivatives
834  double J = QElement<2,3>::dshape_eulerian_at_knot(ipt,psi,dpsidx);
835  //Loop over the test functions and derivatives and set them equal to the
836  //shape functions
837  for(unsigned i=0;i<9;i++)
838  {
839  test[i] = psi[i];
840  dtestdx(i,0) = dpsidx(i,0);
841  dtestdx(i,1) = dpsidx(i,1);
842  }
843  //Return the jacobian
844  return J;
845 }
846 
847 //=======================================================================
848 /// 2D :
849 /// Pressure shape functions
850 //=======================================================================
852 const
853 {
854  psi[0] = 1.0;
855  psi[1] = s[0];
856  psi[2] = s[1];
857 }
858 
859 ///Define the pressure shape and test functions
861 Shape &test) const
862 {
863  //Call the pressure shape functions
864  pshape_pnst(s,psi);
865  //Loop over the test functions and set them equal to the shape functions
866  for(unsigned i=0;i<3;i++) test[i] = psi[i];
867 }
868 
869 //=======================================================================
870 /// Face geometry of the 2D Crouzeix_Raviart elements
871 //=======================================================================
872 template<>
874 {
875  public:
876  FaceGeometry() : QElement<1,3>() {}
877 };
878 
879 ////////////////////////////////////////////////////////////////////////////
880 ////////////////////////////////////////////////////////////////////////////
881 
882 //=======================================================================
883 /// Taylor--Hood elements are Navier--Stokes elements
884 /// with quadratic interpolation for velocities and positions and
885 /// continous linear pressure interpolation
886 //=======================================================================
887 class PolarTaylorHoodElement : public virtual QElement<2,3>,
888  public virtual PolarNavierStokesEquations
889 {
890  private:
891 
892  /// Static array of ints to hold number of variables at node
893  static const unsigned Initial_Nvalue[];
894 
895  protected:
896 
897  /// \short Static array of ints to hold conversion from pressure
898  /// node numbers to actual node numbers
899  static const unsigned Pconv[];
900 
901  /// \short Velocity shape and test functions and their derivs
902  /// w.r.t. to global coords at local coordinate s (taken from geometry)
903  /// Return Jacobian of mapping between local and global coordinates.
904  inline double dshape_and_dtest_eulerian_pnst(const Vector<double> &s, Shape &psi,
905  DShape &dpsidx, Shape &test,
906  DShape &dtestdx) const;
907 
908  /// \short Velocity shape and test functions and their derivs
909  /// w.r.t. to global coords at local coordinate s (taken from geometry)
910  /// Return Jacobian of mapping between local and global coordinates.
911  inline double dshape_and_dtest_eulerian_at_knot_pnst(const unsigned &ipt,
912  Shape &psi,
913  DShape &dpsidx, Shape &test,
914  DShape &dtestdx) const;
915 
916  /// Pressure shape functions at local coordinate s
917  inline void pshape_pnst(const Vector<double> &s, Shape &psi) const;
918 
919  /// Pressure shape and test functions at local coordinte s
920  inline void pshape_pnst(const Vector<double> &s, Shape &psi, Shape &test) const;
921 
922  /// Return the local equation numbers for the pressure values.
923  inline int p_local_eqn(const unsigned &n)
924  {return this->nodal_local_eqn(Pconv[n],p_nodal_index_pnst());}
925 
926 public:
927 
928  /// Constructor, no internal data points
930 
931  /// \short Number of values (pinned or dofs) required at node n. Can
932  /// be overwritten for hanging node version
933  inline virtual unsigned required_nvalue(const unsigned &n) const
934  {return Initial_Nvalue[n];}
935 
936  /// \short Which nodal value represents the pressure?
937  virtual int p_nodal_index_pnst() {return 2;}
938 
939  /// \short Pointer to n_p-th pressure node
940  Node* pressure_node_pt(const unsigned &n_p)
941  {return this->node_pt(Pconv[n_p]);}
942 
943  ///\short Return the velocity component u[i] at local node
944  /// n. Uses suitably interpolated value for hanging nodes.
945  double u_pnst(const unsigned &n, const unsigned &i) const
946  {return this->nodal_value(n,i);}
947 
948  /// \short Return the velocity component u[i] at local node
949  /// n at timestep t (t=0: present; t>0: previous timestep).
950  /// Uses suitably interpolated value for hanging nodes.
951  double u_pnst(const unsigned &t, const unsigned &n,
952  const unsigned &i) const {return this->nodal_value(t,n,i);}
953 
954  /// \short Access function for the pressure values at local pressure
955  /// node n_p (const version)
956  double p_pnst(const unsigned &n_p) const
957  {return this->nodal_value(Pconv[n_p],2);}
958  // {return this->nodal_value(Pconv[n_p],this->p_nodal_index_pnst());}
959 
960  /// Return number of pressure values
961  unsigned npres_pnst() const {return 4;}
962 
963  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
964  void fix_pressure(const unsigned &p_dof, const double &p_value)
965  {
966  this->node_pt(Pconv[p_dof])->pin(this->p_nodal_index_pnst());
967  *this->node_pt(Pconv[p_dof])->value_pt(this->p_nodal_index_pnst()) = p_value;
968  }
969 
970 
971  /// \short Add to the set paired_load_data
972  /// pairs of pointers to data objects and unsigned integers that
973  /// index the values in the data object that affect the load (traction),
974  /// as specified in the get_load() function.
975  void get_load_data(std::set<std::pair<Data*,unsigned> > &paired_load_data);
976 
977  /// Redirect output to NavierStokesEquations output
978  void output(std::ostream &outfile) {PolarNavierStokesEquations::output(outfile);}
979 
980  /// Redirect output to NavierStokesEquations output
981  void output(std::ostream &outfile, const unsigned &Nplot)
982  {PolarNavierStokesEquations::output(outfile,Nplot);}
983 
984  /// Redirect output to NavierStokesEquations output
985  void output(FILE* file_pt) {PolarNavierStokesEquations::output(file_pt);}
986 
987  /// Redirect output to NavierStokesEquations output
988  void output(FILE* file_pt, const unsigned &Nplot)
989  {PolarNavierStokesEquations::output(file_pt,Nplot);}
990 
991 
992 };
993 
994 //Inline functions
995 
996 //==========================================================================
997 /// 2D :
998 /// Derivatives of the shape functions and test functions w.r.t to
999 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1000 /// local and global coordinates.
1001 //==========================================================================
1003  const Vector<double> &s,
1004  Shape &psi,
1005  DShape &dpsidx, Shape &test,
1006  DShape &dtestdx) const
1007 {
1008  //Call the geometrical shape functions and derivatives
1009  double J = QElement<2,3>::dshape_eulerian(s,psi,dpsidx);
1010  //Loop over the test functions and derivatives and set them equal to the
1011  //shape functions
1012  for(unsigned i=0;i<9;i++)
1013  {
1014  test[i] = psi[i];
1015  dtestdx(i,0) = dpsidx(i,0);
1016  dtestdx(i,1) = dpsidx(i,1);
1017  }
1018  //Return the jacobian
1019  return J;
1020 }
1021 
1022 
1023 //==========================================================================
1024 /// 2D :
1025 /// Derivatives of the shape functions and test functions w.r.t to
1026 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1027 /// local and global coordinates.
1028 //==========================================================================
1030  const unsigned &ipt,Shape &psi, DShape &dpsidx, Shape &test,
1031  DShape &dtestdx) const
1032 {
1033  //Call the geometrical shape functions and derivatives
1034  double J = QElement<2,3>::dshape_eulerian_at_knot(ipt,psi,dpsidx);
1035  //Loop over the test functions and derivatives and set them equal to the
1036  //shape functions
1037  for(unsigned i=0;i<9;i++)
1038  {
1039  test[i] = psi[i];
1040  dtestdx(i,0) = dpsidx(i,0);
1041  dtestdx(i,1) = dpsidx(i,1);
1042  }
1043  //Return the jacobian
1044  return J;
1045 }
1046 
1047 
1048 //==========================================================================
1049 /// 2D :
1050 /// Pressure shape functions
1051 //==========================================================================
1053 const
1054 {
1055  //Call the OneDimensional Shape functions
1056  //Local storage
1057  double psi1[2], psi2[2];
1058  //Call the OneDimensional Shape functions
1059  OneDimLagrange::shape<2>(s[0],psi1);
1060  OneDimLagrange::shape<2>(s[1],psi2);
1061 
1062  //Now let's loop over the nodal points in the element
1063  //s1 is the "x" coordinate, s2 the "y"
1064  for(unsigned i=0;i<2;i++)
1065  {
1066  for(unsigned j=0;j<2;j++)
1067  {
1068  /*Multiply the two 1D functions together to get the 2D function*/
1069  psi[2*i + j] = psi2[i]*psi1[j];
1070  }
1071  }
1072 }
1073 
1074 
1075 //==========================================================================
1076 /// 2D :
1077 /// Pressure shape and test functions
1078 //==========================================================================
1080  Shape &test) const
1081 {
1082  //Call the pressure shape functions
1083  pshape_pnst(s,psi);
1084  //Loop over the test functions and set them equal to the shape functions
1085  for(unsigned i=0;i<4;i++) test[i] = psi[i];
1086 }
1087 
1088 //=======================================================================
1089 /// Face geometry of the 2D Taylor_Hood elements
1090 //=======================================================================
1091 template<>
1092 class FaceGeometry<PolarTaylorHoodElement>: public virtual QElement<1,3>
1093 {
1094  public:
1095  FaceGeometry() : QElement<1,3>() {}
1096 };
1097 
1098 }
1099 #endif
PolarCrouzeixRaviartElement()
Constructor, there are DIM+1 internal values (for the pressure)
void get_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs of pointers to data objects and unsigned integers that index th...
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 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.
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number) ...
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
const double & re() const
Reynolds number.
virtual double dshape_and_dtest_eulerian_pnst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and derivatives w.r.t. global coords at local coordinate s...
void output(std::ostream &outfile, const unsigned &Nplot)
Redirect output to NavierStokesEquations output.
double u_pnst(const unsigned &t, const unsigned &n, const unsigned &i) const
Return the velocity component u[i] at local node n at timestep t (t=0: present; t>0: previous timeste...
void output(FILE *file_pt, const unsigned &Nplot)
Redirect output to NavierStokesEquations output.
Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node.
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...
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Compute the element's residual Vector and the jacobian matrix Plus the mass matrix especially for eig...
double interpolated_u_pnst(const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s.
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
virtual unsigned npres_pnst() const =0
Function to return number of pressure degrees of freedom.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
void pshape_pnst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
unsigned npres_pnst() const
Return number of pressure values.
void(* NavierStokesBodyForceFctPt)(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_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
cstr elem_len * i
Definition: cfortran.h:607
void output(FILE *file_pt, const unsigned &Nplot)
Redirect output to NavierStokesEquations output.
double dshape_and_dtest_eulerian_pnst(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 double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
void get_body_force(double time, const Vector< double > &x, Vector< double > &result)
Calculate the body force at a given time and Eulerian position.
void strain_rate_by_r(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Function to return polar strain multiplied by r.
unsigned npres_pnst() const
Return number of pressure values.
double interpolated_p_pnst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
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_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...
A general Finite Element class.
Definition: elements.h:1271
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
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
double kin_energy() const
Get integral of kinetic energy over element.
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
Definition: nodes.h:322
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
virtual void pshape_pnst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
std::string type() const
Return string that indicates the type of the timestepper (e.g. "BDF", "Newmark", etc.)
Definition: timesteppers.h:465
int p_local_eqn(const unsigned &n)
Return the local equation numbers for the pressure values.
virtual double p_pnst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes...
double p_pnst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version) ...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
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...
double du_dt_pnst(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes...
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
void get_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs of pointers to data objects and unsigned integers that index th...
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at node n. Can be overwritten for hanging node version...
double *& re_pt()
Pointer to Reynolds number.
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
double dshape_and_dtest_eulerian_at_knot_pnst(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...
double dshape_and_dtest_eulerian_at_knot_pnst(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...
void interpolated_u_pnst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
const double & re_invfr() const
Global inverse Froude number.
double dshape_and_dtest_eulerian_pnst(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 pshape_pnst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:197
virtual void fill_in_generic_residual_contribution(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 * Alpha_pt
Pointer to the angle alpha.
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)
Compute traction (on the viscous scale) at local coordinate s for outer unit normal N...
virtual int p_nodal_index_pnst()
Which nodal value represents the pressure? (Default: negative, indicating that pressure is not based ...
NavierStokesSourceFctPt & source_fct_pt()
Access function for the source-function pointer.
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 output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number) ...
void output(FILE *file_pt)
C-style output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points...
NavierStokesSourceFctPt source_fct_pt() const
Access function for the source-function pointer. Const version.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: Now returns polar strain.
virtual Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interp...
static char t char * s
Definition: cfortran.h:572
double dissipation() const
Return integral of dissipation over element.
void get_load(const Vector< double > &s, const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
The potential loading on an external SolidElement is always provided by the traction function...
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers...
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
virtual int p_nodal_index_pnst()
Which nodal value represents the pressure?
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
const Vector< double > & g() const
Vector of gravitational components.
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
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 * Re_pt
Pointer to global Reynolds number.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
Vector< double > * G_pt
Pointer to global gravity Vector.
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
double interpolated_dudx_pnst(const Vector< double > &s, const unsigned &i, const unsigned &j) const
void output(std::ostream &outfile, const unsigned &Nplot)
Redirect output to NavierStokesEquations output.
double u_pnst(const unsigned &t, const unsigned &n, const unsigned &i) const
Return the velocity component u[i] at local node n at timestep t (t=0: present; t>0: previous timeste...
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
NavierStokesBodyForceFctPt body_force_fct_pt() const
Access function for the body-force pointer. Const version.
PolarNavierStokesEquations()
Constructor: NULL the body force and source function.
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
PolarTaylorHoodElement()
Constructor, no internal data points.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual Vector and the jacobian matrix Virtual function can be overloaded by h...
double get_source_fct(double time, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
double p_pnst(const unsigned &i_internal) const
Return the pressure values at internal dof i_internal (Discontinous pressure interpolation – no need ...
virtual double dshape_and_dtest_eulerian_at_knot_pnst(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 output(std::ostream &outfile)
Output functionget_vels(const Vector<double>& x_to_get, Vector<double>& vels): x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
virtual int p_local_eqn(const unsigned &n)=0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
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.
virtual unsigned u_index_pnst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.
double(* NavierStokesSourceFctPt)(const double &time, const Vector< double > &x)
Function pointer to source function fct(t,x) x is a Vector!
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element's residual Vector.
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
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 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 u_pnst(const unsigned &n, const unsigned &i) const
Return the velocity component u[i] at local node n. Uses suitably interpolated value for hanging node...
double *& re_invfr_pt()
Pointer to global inverse Froude number.
virtual double u_pnst(const unsigned &n, const unsigned &i) const =0
Velocity i at local node n. Uses suitably interpolated value for hanging nodes.
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
NavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
double u_pnst(const unsigned &n, const unsigned &i) const
Return the velocity component u[i] at local node n. Uses suitably interpolated value for hanging node...
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
int p_local_eqn(const unsigned &n)
Return the local equation numbers for the pressure values.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
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...
double pressure_integral() const
Integral of pressure over element.
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 *& density_ratio_pt()
Pointer to Density ratio.
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:596