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: 1297 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-08-21 13:56:24 +0100 (Mon, 21 Aug 2017) $
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_NAVIER_STOKES_ELEMENTS_HEADER
33 #define OOMPH_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/projection.h"
45 
46 namespace oomph
47 {
48 
49 
50 //======================================================================
51 /// Helper class for elements that impose Robin boundary conditions
52 /// on pressure advection diffusion problem required by Fp preconditioner
53 /// (class used to get around some templating issues)
54 //======================================================================
56 {
57  public:
58 
59  /// Constructor
61 
62  /// Empty virtual destructor
64 
65  /// \short This function returns the residuals for the
66  /// traction function. flag=1 (or 0): do (or don't) compute the
67  /// Jacobian as well.
69  Vector<double> &residuals,
70  DenseMatrix<double> &jacobian,
71  unsigned flag)=0;
72 
73 };
74 
75 
76 
77 ///////////////////////////////////////////////////////////////////////
78 ///////////////////////////////////////////////////////////////////////
79 ///////////////////////////////////////////////////////////////////////
80 
81 
82 
83 //======================================================================
84 /// A class for elements that allow the imposition of Robin boundary
85 /// conditions for the pressure advection diffusion problem in the
86 /// Fp preconditioner.
87 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
88 /// class and and thus, we can be generic enough without the need to have
89 /// a separate equations class
90 //======================================================================
91 template <class ELEMENT>
92 class FpPressureAdvDiffRobinBCElement : public virtual FaceGeometry<ELEMENT>,
93  public virtual FaceElement,
95 {
96 
97  public:
98 
99  ///Constructor, which takes a "bulk" element and the value of the index
100  ///and its limit. Optional boolean flag indicates if it's called
101  // refineable constructor.
103  FiniteElement* const &element_pt,
104  const int &face_index,
105  const bool& called_from_refineable_constructor=false)
106  : FaceGeometry<ELEMENT>(), FaceElement()
107  {
108 
109  //Attach the geometrical information to the element. N.B. This function
110  //also assigns nbulk_value from the required_nvalue of the bulk element
111  element_pt->build_face_element(face_index,this);
112 
113 #ifdef PARANOID
114  {
115  //Check that the element is not a refineable 3d element
116  if (!called_from_refineable_constructor)
117  {
118  //If it's three-d
119  if(element_pt->dim()==3)
120  {
121  //Is it refineable
122  RefineableElement* ref_el_pt=
123  dynamic_cast<RefineableElement*>(element_pt);
124  if(ref_el_pt!=0)
125  {
126  if (this->has_hanging_nodes())
127  {
128  throw OomphLibError(
129  "This flux element will not work correctly if nodes are hanging\n",
130  OOMPH_CURRENT_FUNCTION,
131  OOMPH_EXCEPTION_LOCATION);
132  }
133  }
134  }
135  }
136  }
137 #endif
138 
139  }
140 
141  /// Empty destructor
143 
144  /// \short This function returns the residuals for the
145  /// traction function. flag=1 (or 0): do (or don't) compute the
146  /// Jacobian as well.
148  Vector<double> &residuals,
149  DenseMatrix<double> &jacobian,
150  unsigned flag);
151 
152 
153  ///This function returns just the residuals
155  {
156  std::ostringstream error_message;
157  error_message
158  << "fill_in_contribution_to_residuals() must not be called directly.\n"
159  << "since it uses the local equation numbering of the bulk element\n"
160  << "which calls the relevant helper function directly.\n";
161  throw OomphLibError(
162  error_message.str(),
163  OOMPH_CURRENT_FUNCTION,
164  OOMPH_EXCEPTION_LOCATION);
165  }
166 
167  ///This function returns the residuals and the jacobian
169  DenseMatrix<double> &jacobian)
170  {
171  std::ostringstream error_message;
172  error_message
173  << "fill_in_contribution_to_jacobian() must not be called directly.\n"
174  << "since it uses the local equation numbering of the bulk element\n"
175  << "which calls the relevant helper function directly.\n";
176  throw OomphLibError(
177  error_message.str(),
178  OOMPH_CURRENT_FUNCTION,
179  OOMPH_EXCEPTION_LOCATION);
180  }
181 
182  ///Overload the output function
183  void output(std::ostream &outfile) {FiniteElement::output(outfile);}
184 
185  ///Output function: x,y,[z],u,v,[w],p in tecplot format
186  void output(std::ostream &outfile, const unsigned &nplot)
187  {FiniteElement::output(outfile,nplot);}
188 
189 };
190 
191 ///////////////////////////////////////////////////////////////////////
192 ///////////////////////////////////////////////////////////////////////
193 ///////////////////////////////////////////////////////////////////////
194 
195 
196 
197 
198 //============================================================================
199 /// Get residuals and Jacobian of Robin boundary conditions in pressure
200 /// advection diffusion problem in Fp preconditoner
201 //============================================================================
202 template<class ELEMENT>
205  Vector<double> &residuals,
206  DenseMatrix<double> &jacobian,
207  unsigned flag)
208 {
209  //Storage for local coordinates in FaceElement and associted bulk element
210  unsigned my_dim=this->dim();
211  Vector<double> s(my_dim);
212  Vector<double> s_bulk(my_dim+1);
213 
214  // Storage for outer unit normal
215  Vector<double> unit_normal(my_dim+1);
216 
217  //Storage for veloc in bulk element
218  Vector<double> veloc(my_dim+1);
219 
220  //Set the value of n_intpt
221  unsigned n_intpt = this->integral_pt()->nweight();
222 
223  //Integers to store local equation numbers
224  int local_eqn=0;
225  int local_unknown=0;
226 
227  // Get cast bulk element
228  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(this->bulk_element_pt());
229 
230  //Find out how many pressure dofs there are in the bulk element
231  unsigned n_pres = bulk_el_pt->npres_nst();
232 
233  // Get the Reynolds number from the bulk element
234  double re = bulk_el_pt->re();
235 
236  //Set up memory for pressure shape and test functions
237  Shape psip(n_pres), testp(n_pres);
238 
239  //Loop over the integration points
240  for(unsigned ipt=0;ipt<n_intpt;ipt++)
241  {
242  //Get the integral weight
243  double w = this->integral_pt()->weight(ipt);
244 
245  //Assign values of local coordinate in FaceElement
246  for(unsigned i=0;i<my_dim;i++) s[i] = this->integral_pt()->knot(ipt,i);
247 
248  // Find corresponding coordinate in the the bulk element
249  s_bulk=this->local_coordinate_in_bulk(s);
250 
251  /// Get outer unit normal
252  this->outer_unit_normal(ipt,unit_normal);
253 
254  // Get velocity in bulk element
255  bulk_el_pt->interpolated_u_nst(s_bulk,veloc);
256 
257  // Get normal component of veloc
258  double flux=0.0;
259  for (unsigned i=0;i<my_dim+1;i++)
260  {
261  flux+=veloc[i]*unit_normal[i];
262  }
263 
264  // Modify bc: If outflow (flux>0) apply Neumann condition instead
265  if (flux>0.0) flux=0.0;
266 
267  // Get pressure
268  double interpolated_press=bulk_el_pt->interpolated_p_nst(s_bulk);
269 
270  //Call the pressure shape and test functions in bulk element
271  bulk_el_pt->pshape_nst(s_bulk,psip,testp);
272 
273  //Find the Jacobian of the mapping within the FaceElement
274  double J = this->J_eulerian(s);
275 
276  //Premultiply the weights and the Jacobian
277  double W = w*J;
278 
279  //Loop over the pressure shape functions in bulk
280  //(wasteful but they'll be zero on the boundary)
281  for(unsigned l=0;l<n_pres;l++)
282  {
283  local_eqn=bulk_el_pt->p_local_eqn(l);
284 
285  //If not a boundary conditions
286  if(local_eqn >= 0)
287  {
288  residuals[local_eqn] -=
289  re*flux*interpolated_press*testp[l]*W;
290 
291  // Jacobian too?
292  if(flag)
293  {
294  //Loop over the shape functions in bulk
295  for(unsigned l2=0;l2<n_pres;l2++)
296  {
297  local_unknown = bulk_el_pt->p_local_eqn(l2);
298 
299  //If not a boundary conditions
300  if(local_unknown >= 0)
301  {
302  jacobian(local_eqn,local_unknown)-=
303  re*flux*psip[l2]*testp[l]*W;
304  }
305  }
306  } /*End of Jacobian calculation*/
307  } //End of if not boundary condition
308  }//End of loop over l
309  }
310 
311 }
312 
313 ///////////////////////////////////////////////////////////////////////
314 ///////////////////////////////////////////////////////////////////////
315 ///////////////////////////////////////////////////////////////////////
316 
317 
318 //======================================================================
319 /// Template-free base class for Navier-Stokes equations to avoid
320 /// casting problems
321 //======================================================================
324  public virtual FiniteElement
325 {
326 
327  public:
328 
329  /// Constructor (empty)
331 
332  /// Virtual destructor (empty)
334 
335  /// \short Compute the residuals for the associated pressure advection
336  /// diffusion problem. Used by the Fp preconditioner.
338  residuals)=0;
339 
340  /// \short Compute the residuals and Jacobian for the associated
341  /// pressure advection diffusion problem. Used by the Fp preconditioner.
343  Vector<double>& residuals, DenseMatrix<double> &jacobian)=0;
344 
345  /// \short Return the index at which the pressure is stored if it is
346  /// stored at the nodes. If not stored at the nodes this will return
347  /// a negative number.
348  virtual int p_nodal_index_nst() const =0;
349 
350  /// \short Access function for the local equation number information for
351  /// the pressure.
352  /// p_local_eqn[n] = local equation number or < 0 if pinned
353  virtual int p_local_eqn(const unsigned &n)const=0;
354 
355  /// \short Global eqn number of pressure dof that's pinned in pressure
356  /// adv diff problem
357  virtual int& pinned_fp_pressure_eqn()=0;
358 
359 
360  /// \short Pin all non-pressure dofs and backup eqn numbers of all Data
361  virtual void pin_all_non_pressure_dofs(std::map<Data*,std::vector<int> >&
362  eqn_number_backup)=0;
363 
364  /// \short Build FaceElements that apply the Robin boundary condition
365  /// to the pressure advection diffusion problem required by
366  /// Fp preconditioner
367  virtual void build_fp_press_adv_diff_robin_bc_element(const unsigned&
368  face_index)=0;
369 
370  /// \short Delete the FaceElements that apply the Robin boundary condition
371  /// to the pressure advection diffusion problem required by
372  /// Fp preconditioner
374 
375 
376  /// \short Compute the diagonal of the velocity/pressure mass matrices.
377  /// If which one=0, both are computed, otherwise only the pressure
378  /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
379  /// LSC version of the preconditioner only needs that one)
381  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
382  const unsigned& which_one=0)=0;
383 
384 };
385 
386 
387 ///////////////////////////////////////////////////////////////////////
388 ///////////////////////////////////////////////////////////////////////
389 ///////////////////////////////////////////////////////////////////////
390 
391 
392 
393 
394 
395 //======================================================================
396 /// A class for elements that solve the cartesian Navier--Stokes equations,
397 /// templated by the dimension DIM.
398 /// This contains the generic maths -- any concrete implementation must
399 /// be derived from this.
400 ///
401 /// We're solving:
402 ///
403 /// \f$ { Re \left( St \frac{\partial u_i}{\partial t} +
404 /// (u_j - u_j^{M}) \frac{\partial u_i}{\partial x_j} \right) =
405 /// - \frac{\partial p}{\partial x_i} - R_\rho B_i(x_j) -
406 /// \frac{Re}{Fr} G_i +
407 /// \frac{\partial }{\partial x_j} \left[ R_\mu \left(
408 /// \frac{\partial u_i}{\partial x_j} +
409 /// \frac{\partial u_j}{\partial x_i} \right) \right] } \f$
410 ///
411 /// and
412 ///
413 /// \f$ { \frac{\partial u_i}{\partial x_i} = Q } \f$
414 ///
415 /// We also provide all functions required to use this element
416 /// in FSI problems, by deriving it from the FSIFluidElement base
417 /// class.
418 //======================================================================
419 template <unsigned DIM>
420  class NavierStokesEquations : public virtual FSIFluidElement,
422 {
423 
424  public:
425 
426  /// \short Function pointer to body force function fct(t,x,f(x))
427  /// x is a Vector!
428  typedef void (*NavierStokesBodyForceFctPt)(const double& time,
429  const Vector<double>& x,
430  Vector<double>& body_force);
431 
432  /// \short Function pointer to source function fct(t,x)
433  /// x is a Vector!
434  typedef double (*NavierStokesSourceFctPt)(const double& time,
435  const Vector<double>& x);
436 
437 
438  /// \short Function pointer to source function fct(x) for the
439  /// pressure advection diffusion equation (only used during
440  /// validation!). x is a Vector!
442  const Vector<double>& x);
443 
444  private:
445 
446  /// \short Static "magic" number that indicates that the pressure is
447  /// not stored at a node
449 
450  /// Static default value for the physical constants (all initialised to zero)
452 
453  /// Static default value for the physical ratios (all are initialised to one)
455 
456  /// Static default value for the gravity vector
458 
459  protected:
460 
461  //Physical constants
462 
463  /// \short Pointer to the viscosity ratio (relative to the
464  /// viscosity used in the definition of the Reynolds number)
466 
467  /// \short Pointer to the density ratio (relative to the density used in the
468  /// definition of the Reynolds number)
470 
471  // Pointers to global physical constants
472 
473  /// Pointer to global Reynolds number
474  double *Re_pt;
475 
476  /// Pointer to global Reynolds number x Strouhal number (=Womersley)
477  double *ReSt_pt;
478 
479  /// \short Pointer to global Reynolds number x inverse Froude number
480  /// (= Bond number / Capillary number)
481  double *ReInvFr_pt;
482 
483  /// Pointer to global gravity Vector
485 
486  /// Pointer to body force function
488 
489  /// Pointer to volumetric source function
491 
492  /// \short Pointer to source function pressure advection diffusion equation
493  /// (only to be used during validation)
495 
496  /// \short Boolean flag to indicate if ALE formulation is disabled when
497  /// time-derivatives are computed. Only set to true if you're sure
498  /// that the mesh is stationary.
500 
501  /// \short Storage for FaceElements that apply Robin BC for pressure adv diff
502  /// equation used in Fp preconditioner.
505 
506  /// \short Global eqn number of pressure dof that's pinned in
507  /// pressure advection diffusion problem (defaults to -1)
509 
510  /// \short Compute the shape functions and derivatives
511  /// w.r.t. global coords at local coordinate s.
512  /// Return Jacobian of mapping between local and global coordinates.
513  virtual double dshape_and_dtest_eulerian_nst(const Vector<double> &s,
514  Shape &psi,
515  DShape &dpsidx, Shape &test,
516  DShape &dtestdx) const=0;
517 
518  /// \short Compute the shape functions and derivatives
519  /// w.r.t. global coords at ipt-th integration point
520  /// Return Jacobian of mapping between local and global coordinates.
521  virtual double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt,
522  Shape &psi,
523  DShape &dpsidx,
524  Shape &test,
525  DShape &dtestdx) const=0;
526 
527  /// \short Shape/test functions and derivs w.r.t. to global coords at
528  /// integration point ipt; return Jacobian of mapping (J). Also compute
529  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
531  const unsigned &ipt,
532  Shape &psi,
533  DShape &dpsidx,
534  RankFourTensor<double> &d_dpsidx_dX,
535  Shape &test,
536  DShape &dtestdx,
537  RankFourTensor<double> &d_dtestdx_dX,
538  DenseMatrix<double> &djacobian_dX) const=0;
539 
540  /// \short Calculate the body force at a given time and local and/or
541  /// Eulerian position. This function is virtual so that it can be
542  /// overloaded in multi-physics elements where the body force might
543  /// depend on another variable.
544  virtual void get_body_force_nst(const double& time,
545  const unsigned& ipt,
546  const Vector<double> &s,
547  const Vector<double> &x,
548  Vector<double> &result)
549  {
550  //If the function pointer is zero return zero
551  if(Body_force_fct_pt == 0)
552  {
553  //Loop over dimensions and set body forces to zero
554  for(unsigned i=0;i<DIM;i++) {result[i] = 0.0;}
555  }
556  //Otherwise call the function
557  else
558  {
559  (*Body_force_fct_pt)(time,x,result);
560  }
561  }
562 
563  /// Get gradient of body force term at (Eulerian) position x. This function is
564  /// virtual to allow overloading in multi-physics problems where
565  /// the strength of the source function might be determined by
566  /// another system of equations. Computed via function pointer
567  /// (if set) or by finite differencing (default)
568  inline virtual void get_body_force_gradient_nst(
569  const double& time,
570  const unsigned& ipt,
571  const Vector<double>& s,
572  const Vector<double>& x,
573  DenseMatrix<double>& d_body_force_dx)
574  {
575 // hierher: Implement function pointer version
576 /* //If no gradient function has been set, FD it */
577 /* if(Body_force_fct_gradient_pt==0) */
578  {
579  // Reference value
580  Vector<double> body_force(DIM,0.0);
581  get_body_force_nst(time,ipt,s,x,body_force);
582 
583  // FD it
585  Vector<double> body_force_pls(DIM,0.0);
586  Vector<double> x_pls(x);
587  for (unsigned i=0;i<DIM;i++)
588  {
589  x_pls[i]+=eps_fd;
590  get_body_force_nst(time,ipt,s,x_pls,body_force_pls);
591  for (unsigned j=0;j<DIM;j++)
592  {
593  d_body_force_dx(j,i)=(body_force_pls[j]-body_force[j])/eps_fd;
594  }
595  x_pls[i]=x[i];
596  }
597  }
598 /* else */
599 /* { */
600 /* // Get gradient */
601 /* (*Source_fct_gradient_pt)(time,x,gradient); */
602 /* } */
603  }
604 
605 
606 
607  /// \short Calculate the source fct at given time and
608  /// Eulerian position
609  virtual double get_source_nst(const double& time, const unsigned& ipt,
610  const Vector<double> &x)
611  {
612  //If the function pointer is zero return zero
613  if (Source_fct_pt == 0) {return 0;}
614  //Otherwise call the function
615  else {return (*Source_fct_pt)(time,x);}
616  }
617 
618 
619  /// Get gradient of source term at (Eulerian) position x. This function is
620  /// virtual to allow overloading in multi-physics problems where
621  /// the strength of the source function might be determined by
622  /// another system of equations. Computed via function pointer
623  /// (if set) or by finite differencing (default)
624  inline virtual void get_source_gradient_nst(
625  const double& time,
626  const unsigned& ipt,
627  const Vector<double>& x,
628  Vector<double>& gradient)
629  {
630 // hierher: Implement function pointer version
631 /* //If no gradient function has been set, FD it */
632 /* if(Source_fct_gradient_pt==0) */
633  {
634  // Reference value
635  double source=get_source_nst(time,ipt,x);
636 
637  // FD it
639  double source_pls=0.0;
640  Vector<double> x_pls(x);
641  for (unsigned i=0;i<DIM;i++)
642  {
643  x_pls[i]+=eps_fd;
644  source_pls=get_source_nst(time,ipt,x_pls);
645  gradient[i]=(source_pls-source)/eps_fd;
646  x_pls[i]=x[i];
647  }
648  }
649 /* else */
650 /* { */
651 /* // Get gradient */
652 /* (*Source_fct_gradient_pt)(time,x,gradient); */
653 /* } */
654  }
655 
656 
657  ///\short Compute the residuals for the Navier--Stokes equations.
658  /// Flag=1 (or 0): do (or don't) compute the Jacobian as well.
659  /// Flag=2: Fill in mass matrix too.
661  Vector<double> &residuals, DenseMatrix<double> &jacobian,
662  DenseMatrix<double> &mass_matrix, unsigned flag);
663 
664 
665  /// \short Compute the residuals for the associated pressure advection
666  /// diffusion problem. Used by the Fp preconditioner.
667  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
669  Vector<double> &residuals, DenseMatrix<double> &jacobian, unsigned flag);
670 
671  ///\short Compute the derivatives of the
672  /// residuals for the Navier--Stokes equations with respect to a parameter
673  /// Flag=1 (or 0): do (or don't) compute the Jacobian as well.
674  /// Flag=2: Fill in mass matrix too.
676  double* const &parameter_pt,
677  Vector<double> &dres_dparam, DenseMatrix<double> &djac_dparam,
678  DenseMatrix<double> &dmass_matrix_dparam, unsigned flag);
679 
680  /// \short Compute the hessian tensor vector products required to
681  /// perform continuation of bifurcations analytically
683  Vector<double> const &Y,
684  DenseMatrix<double> const &C,
685  DenseMatrix<double> &product);
686 
687 
688 public:
689 
690  /// \short Constructor: NULL the body force and source function
691  /// and make sure the ALE terms are included by default.
695  {
696  //Set all the Physical parameter pointers to the default value zero
701  //Set the Physical ratios to the default value of 1
704  }
705 
706  /// Vector to decide whether the stress-divergence form is used or not
707  // N.B. This needs to be public so that the intel compiler gets things correct
708  // somehow the access function messes things up when going to refineable
709  // navier--stokes
711 
712  //Access functions for the physical constants
713 
714  /// Reynolds number
715  const double &re() const {return *Re_pt;}
716 
717  /// Product of Reynolds and Strouhal number (=Womersley number)
718  const double &re_st() const {return *ReSt_pt;}
719 
720  /// Pointer to Reynolds number
721  double* &re_pt() {return Re_pt;}
722 
723  /// Pointer to product of Reynolds and Strouhal number (=Womersley number)
724  double* &re_st_pt() {return ReSt_pt;}
725 
726  /// \short Viscosity ratio for element: Element's viscosity relative to the
727  /// viscosity used in the definition of the Reynolds number
728  const double &viscosity_ratio() const {return *Viscosity_Ratio_pt;}
729 
730  /// Pointer to Viscosity Ratio
732 
733  /// \short Density ratio for element: Element's density relative to the
734  /// viscosity used in the definition of the Reynolds number
735  const double &density_ratio() const {return *Density_Ratio_pt;}
736 
737  /// Pointer to Density ratio
738  double* &density_ratio_pt() {return Density_Ratio_pt;}
739 
740  /// Global inverse Froude number
741  const double &re_invfr() const {return *ReInvFr_pt;}
742 
743  /// Pointer to global inverse Froude number
744  double* &re_invfr_pt() {return ReInvFr_pt;}
745 
746  /// Vector of gravitational components
747  const Vector<double> &g() const {return *G_pt;}
748 
749  /// Pointer to Vector of gravitational components
750  Vector<double>* &g_pt() {return G_pt;}
751 
752  /// Access function for the body-force pointer
754  {return Body_force_fct_pt;}
755 
756  /// Access function for the body-force pointer. Const version
758  {return Body_force_fct_pt;}
759 
760  ///Access function for the source-function pointer
762 
763  ///Access function for the source-function pointer. Const version
765 
766  /// \short Access function for the source-function pointer for pressure
767  /// advection diffusion (used for validation only).
770 
771  /// \short Access function for the source-function pointer for pressure
772  /// advection diffusion (used for validation only). Const version.
774  const
776 
777  /// \short Global eqn number of pressure dof that's pinned in pressure
778  /// adv diff problem
780 
781  ///Function to return number of pressure degrees of freedom
782  virtual unsigned npres_nst() const=0;
783 
784  /// Compute the pressure shape functions at local coordinate s
785  virtual void pshape_nst(const Vector<double> &s, Shape &psi) const=0;
786 
787  /// \short Compute the pressure shape and test functions
788  /// at local coordinate s
789  virtual void pshape_nst(const Vector<double> &s, Shape &psi,
790  Shape &test) const=0;
791 
792  /// \short Compute the pressure shape and test functions and derivatives
793  /// w.r.t. global coords at local coordinate s.
794  /// Return Jacobian of mapping between local and global coordinates.
795  virtual double dpshape_and_dptest_eulerian_nst(const Vector<double> &s,
796  Shape &ppsi,
797  DShape &dppsidx,
798  Shape &ptest,
799  DShape &dptestdx) const=0;
800 
801  /// \short Velocity i at local node n. Uses suitably interpolated value
802  /// for hanging nodes. The use of u_index_nst() permits the use of this
803  /// element as the basis for multi-physics elements. The default
804  /// is to assume that the i-th velocity component is stored at the
805  /// i-th location of the node
806  double u_nst(const unsigned &n, const unsigned &i) const
807  {return nodal_value(n,u_index_nst(i));}
808 
809  /// \short Velocity i at local node n at timestep t (t=0: present;
810  /// t>0: previous). Uses suitably interpolated value for hanging nodes.
811  double u_nst(const unsigned &t, const unsigned &n,
812  const unsigned &i) const
813  {return nodal_value(t,n,u_index_nst(i));}
814 
815  /// \short Return the index at which the i-th unknown velocity component
816  /// is stored. The default value, i, is appropriate for single-physics
817  /// problems.
818  /// In derived multi-physics elements, this function should be overloaded
819  /// to reflect the chosen storage scheme. Note that these equations require
820  /// that the unknowns are always stored at the same indices at each node.
821  virtual inline unsigned u_index_nst(const unsigned &i) const {return i;}
822 
823  /// \short Return the number of velocity components
824  /// Used in FluidInterfaceElements
825  inline unsigned n_u_nst() const {return DIM;}
826 
827  /// \short i-th component of du/dt at local node n.
828  /// Uses suitably interpolated value for hanging nodes.
829  double du_dt_nst(const unsigned &n, const unsigned &i) const
830  {
831  // Get the data's timestepper
833 
834  //Initialise dudt
835  double dudt=0.0;
836 
837  //Loop over the timesteps, if there is a non Steady timestepper
838  if (!time_stepper_pt->is_steady())
839  {
840  //Find the index at which the dof is stored
841  const unsigned u_nodal_index = this->u_index_nst(i);
842 
843  // Number of timsteps (past & present)
844  const unsigned n_time = time_stepper_pt->ntstorage();
845  // Loop over the timesteps
846  for(unsigned t=0;t<n_time;t++)
847  {
848  dudt+=time_stepper_pt->weight(1,t)*nodal_value(t,n,u_nodal_index);
849  }
850  }
851 
852  return dudt;
853  }
854 
855  /// \short Disable ALE, i.e. assert the mesh is not moving -- you do this
856  /// at your own risk!
857  void disable_ALE()
858  {
859  ALE_is_disabled=true;
860  }
861 
862  /// \short (Re-)enable ALE, i.e. take possible mesh motion into account
863  /// when evaluating the time-derivative. Note: By default, ALE is
864  /// enabled, at the expense of possibly creating unnecessary work
865  /// in problems where the mesh is, in fact, stationary.
866  void enable_ALE()
867  {
868  ALE_is_disabled=false;
869  }
870 
871  /// \short Pressure at local pressure "node" n_p
872  /// Uses suitably interpolated value for hanging nodes.
873  virtual double p_nst(const unsigned &n_p)const=0;
874 
875  /// \short Pressure at local pressure "node" n_p at time level t
876  virtual double p_nst(const unsigned &t, const unsigned &n_p)const=0;
877 
878  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
879  virtual void fix_pressure(const unsigned &p_dof, const double &p_value)=0;
880 
881  /// \short Return the index at which the pressure is stored if it is
882  /// stored at the nodes. If not stored at the nodes this will return
883  /// a negative number.
884  virtual int p_nodal_index_nst() const {return Pressure_not_stored_at_node;}
885 
886  /// Integral of pressure over element
887  double pressure_integral() const;
888 
889  /// \short Return integral of dissipation over element
890  double dissipation() const;
891 
892  /// \short Return dissipation at local coordinate s
893  double dissipation(const Vector<double>& s) const;
894 
895  /// \short Compute the vorticity vector at local coordinate s
896  void get_vorticity(const Vector<double>& s, Vector<double>& vorticity) const;
897 
898  /// \short Get integral of kinetic energy over element
899  double kin_energy() const;
900 
901  /// \short Get integral of time derivative of kinetic energy over element
902  double d_kin_energy_dt() const;
903 
904  /// Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)
905  void strain_rate(const Vector<double>& s,
907 
908  /// \short Compute traction (on the viscous scale) exerted onto
909  /// the fluid at local coordinate s. N has to be outer unit normal
910  /// to the fluid.
911  void get_traction(const Vector<double>& s, const Vector<double>& N,
912  Vector<double>& traction);
913 
914  /// \short Compute traction (on the viscous scale) exerted onto
915  /// the fluid at local coordinate s, decomposed into pressure and
916  /// normal and tangential viscous components.
917  /// N has to be outer unit normal to the fluid.
918  void get_traction(const Vector<double>& s, const Vector<double>& N,
919  Vector<double>& traction_p,
920  Vector<double>& traction_visc_n,
921  Vector<double>& traction_visc_t);
922 
923  /// \short This implements a pure virtual function defined
924  /// in the FSIFluidElement class. The function computes
925  /// the traction (on the viscous scale), at the element's local
926  /// coordinate s, that the fluid element exerts onto an adjacent
927  /// solid element. The number of arguments is imposed by
928  /// the interface defined in the FSIFluidElement -- only the
929  /// unit normal N (pointing into the fluid!) is actually used
930  /// in the computation.
931  void get_load(const Vector<double> &s,
932  const Vector<double> &N,
933  Vector<double> &load)
934  {
935  // Note: get_traction() computes the traction onto the fluid
936  // if N is the outer unit normal onto the fluid; here we're
937  // exepcting N to point into the fluid so we're getting the
938  // traction onto the adjacent wall instead!
939  get_traction(s,N,load);
940  }
941 
942  /// \short Compute the diagonal of the velocity/pressure mass matrices.
943  /// If which one=0, both are computed, otherwise only the pressure
944  /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
945  /// LSC version of the preconditioner only needs that one)
947  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
948  const unsigned& which_one=0);
949 
950  /// \short Number of scalars/fields output by this element. Reimplements
951  /// broken virtual function in base class.
952  unsigned nscalar_paraview() const
953  {
954  return DIM+1;
955  }
956 
957  /// \short Write values of the i-th scalar field at the plot points. Needs
958  /// to be implemented for each new specific element type.
959  void scalar_value_paraview(std::ofstream& file_out,
960  const unsigned& i,
961  const unsigned& nplot) const
962  {
963  // Vector of local coordinates
964  Vector<double> s(DIM);
965 
966 
967  // Loop over plot points
968  unsigned num_plot_points=nplot_points_paraview(nplot);
969  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
970  {
971 
972  // Get local coordinates of plot point
973  get_s_plot(iplot,nplot,s);
974 
975  // Velocities
976  if(i<DIM) {file_out << interpolated_u_nst(s,i) << std::endl;}
977 
978  // Pressure
979  else if(i==DIM) {file_out << interpolated_p_nst(s) << std::endl;}
980 
981  // Never get here
982  else
983  {
984 #ifdef PARANOID
985  std::stringstream error_stream;
986  error_stream
987  << "These Navier Stokes elements only store " << DIM+1 << " fields, "
988  << "but i is currently " << i << std::endl;
989  throw OomphLibError(
990  error_stream.str(),
991  OOMPH_CURRENT_FUNCTION,
992  OOMPH_EXCEPTION_LOCATION);
993 #endif
994  }
995  }
996  }
997 
998  /// \short Name of the i-th scalar field. Default implementation
999  /// returns V1 for the first one, V2 for the second etc. Can (should!) be
1000  /// overloaded with more meaningful names in specific elements.
1001  std::string scalar_name_paraview(const unsigned& i) const
1002  {
1003  // Velocities
1004  if(i<DIM)
1005  {
1006  return "Velocity "+StringConversion::to_string(i);
1007  }
1008  // Preussre
1009  else if(i==DIM)
1010  {
1011  return "Pressure";
1012  }
1013  // Never get here
1014  else
1015  {
1016  std::stringstream error_stream;
1017  error_stream
1018  << "These Navier Stokes elements only store " << DIM+1 << " fields,\n"
1019  << "but i is currently " << i << std::endl;
1020  throw OomphLibError(
1021  error_stream.str(),
1022  OOMPH_CURRENT_FUNCTION,
1023  OOMPH_EXCEPTION_LOCATION);
1024  // Dummy return
1025  return " ";
1026  }
1027  }
1028 
1029  /// \short Output function: x,y,[z],u,v,[w],p
1030  /// in tecplot format. Default number of plot points
1031  void output(std::ostream &outfile)
1032  {
1033  unsigned nplot=5;
1034  output(outfile,nplot);
1035  }
1036 
1037  /// \short Output function: x,y,[z],u,v,[w],p
1038  /// in tecplot format. nplot points in each coordinate direction
1039  void output(std::ostream &outfile, const unsigned &nplot);
1040 
1041  /// \short C-style output function: x,y,[z],u,v,[w],p
1042  /// in tecplot format. Default number of plot points
1043  void output(FILE* file_pt)
1044  {
1045  unsigned nplot=5;
1046  output(file_pt,nplot);
1047  }
1048 
1049  /// \short C-style output function: x,y,[z],u,v,[w],p
1050  /// in tecplot format. nplot points in each coordinate direction
1051  void output(FILE* file_pt, const unsigned &nplot);
1052 
1053  /// \short Full output function:
1054  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
1055  /// in tecplot format. Default number of plot points
1056  void full_output(std::ostream &outfile)
1057  {
1058  unsigned nplot=5;
1059  full_output(outfile,nplot);
1060  }
1061 
1062  /// \short Full output function:
1063  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
1064  /// in tecplot format. nplot points in each coordinate direction
1065  void full_output(std::ostream &outfile, const unsigned &nplot);
1066 
1067 
1068  /// \short Output function: x,y,[z],u,v,[w] in tecplot format.
1069  /// nplot points in each coordinate direction at timestep t
1070  /// (t=0: present; t>0: previous timestep)
1071  void output_veloc(std::ostream &outfile, const unsigned &nplot,
1072  const unsigned& t);
1073 
1074 
1075  /// \short Output function: x,y,[z], [omega_x,omega_y,[and/or omega_z]]
1076  /// in tecplot format. nplot points in each coordinate direction
1077  void output_vorticity(std::ostream &outfile,
1078  const unsigned &nplot);
1079 
1080  /// \short Output exact solution specified via function pointer
1081  /// at a given number of plot points. Function prints as
1082  /// many components as are returned in solution Vector
1083  void output_fct(std::ostream &outfile, const unsigned &nplot,
1085 
1086  /// \short Output exact solution specified via function pointer
1087  /// at a given time and at a given number of plot points.
1088  /// Function prints as many components as are returned in solution Vector.
1089  void output_fct(std::ostream &outfile, const unsigned &nplot,
1090  const double& time,
1092 
1093  /// \short Validate against exact solution at given time
1094  /// Solution is provided via function pointer.
1095  /// Plot at a given number of plot points and compute L2 error
1096  /// and L2 norm of velocity solution over element
1097  void compute_error(std::ostream &outfile,
1099  const double& time,
1100  double& error, double& norm);
1101 
1102  /// \short Validate against exact solution.
1103  /// Solution is provided via function pointer.
1104  /// Plot at a given number of plot points and compute L2 error
1105  /// and L2 norm of velocity solution over element
1106  void compute_error(std::ostream &outfile,
1108  double& error, double& norm);
1109 
1110  /// Compute the element's residual Vector
1112  {
1113  //Call the generic residuals function with flag set to 0
1114  //and using a dummy matrix argument
1118  }
1119 
1120  ///\short Compute the element's residual Vector and the jacobian matrix
1121  /// Virtual function can be overloaded by hanging-node version
1123  DenseMatrix<double> &jacobian)
1124  {
1125  //Call the generic routine with the flag set to 1
1126  fill_in_generic_residual_contribution_nst(residuals,jacobian,
1128  }
1129 
1130  /// \short Add the element's contribution to its residuals vector,
1131  /// jacobian matrix and mass matrix
1133  Vector<double> &residuals, DenseMatrix<double> &jacobian,
1134  DenseMatrix<double> &mass_matrix)
1135  {
1136  //Call the generic routine with the flag set to 2
1137  fill_in_generic_residual_contribution_nst(residuals,jacobian,mass_matrix,2);
1138  }
1139 
1140  /// Compute the element's residual Vector
1142  double* const &parameter_pt,Vector<double> &dres_dparam)
1143  {
1144  //Call the generic residuals function with flag set to 0
1145  //and using a dummy matrix argument
1147  parameter_pt,
1150  }
1151 
1152  ///\short Compute the element's residual Vector and the jacobian matrix
1153  /// Virtual function can be overloaded by hanging-node version
1155  double* const &parameter_pt,Vector<double> &dres_dparam,
1156  DenseMatrix<double> &djac_dparam)
1157  {
1158  //Call the generic routine with the flag set to 1
1160  parameter_pt,
1161  dres_dparam,djac_dparam,GeneralisedElement::Dummy_matrix,1);
1162  }
1163 
1164  /// Add the element's contribution to its residuals vector,
1165  /// jacobian matrix and mass matrix
1167  double* const &parameter_pt,
1168  Vector<double> &dres_dparam,
1169  DenseMatrix<double> &djac_dparam,
1170  DenseMatrix<double> &dmass_matrix_dparam)
1171  {
1172  //Call the generic routine with the flag set to 2
1174  parameter_pt,dres_dparam,djac_dparam,dmass_matrix_dparam,2);
1175  }
1176 
1177 
1178  /// \short Compute the residuals for the associated pressure advection
1179  /// diffusion problem. Used by the Fp preconditioner.
1181  {
1183  residuals,GeneralisedElement::Dummy_matrix,0);
1184  }
1185 
1186  /// \short Compute the residuals and Jacobian for the associated
1187  /// pressure advection diffusion problem. Used by the Fp preconditioner.
1189  Vector<double>& residuals, DenseMatrix<double> &jacobian)
1190  {
1192  residuals,jacobian,1);
1193  }
1194 
1195 
1196  /// \short Pin all non-pressure dofs and backup eqn numbers
1197  void pin_all_non_pressure_dofs(std::map<Data*,std::vector<int> >&
1198  eqn_number_backup)
1199  {
1200  // Loop over internal data and pin the values (having established that
1201  // pressure dofs aren't amongst those)
1202  unsigned nint=this->ninternal_data();
1203  for (unsigned j=0;j<nint;j++)
1204  {
1205  Data* data_pt=this->internal_data_pt(j);
1206  if (eqn_number_backup[data_pt].size()==0)
1207  {
1208  unsigned nvalue=data_pt->nvalue();
1209  eqn_number_backup[data_pt].resize(nvalue);
1210  for (unsigned i=0;i<nvalue;i++)
1211  {
1212  // Backup
1213  eqn_number_backup[data_pt][i]=data_pt->eqn_number(i);
1214 
1215  // Pin everything
1216  data_pt->pin(i);
1217  }
1218  }
1219  }
1220 
1221  // Now deal with nodal values
1222  unsigned nnod=this->nnode();
1223  for (unsigned j=0;j<nnod;j++)
1224  {
1225 
1226  Node* nod_pt=this->node_pt(j);
1227  if (eqn_number_backup[nod_pt].size()==0)
1228  {
1229 
1230  unsigned nvalue=nod_pt->nvalue();
1231  eqn_number_backup[nod_pt].resize(nvalue);
1232  for (unsigned i=0;i<nvalue;i++)
1233  {
1234  // Pin everything apart from the nodal pressure
1235  // value
1236  if (int(i)!=this->p_nodal_index_nst())
1237  {
1238  // Backup
1239  eqn_number_backup[nod_pt][i]=nod_pt->eqn_number(i);
1240 
1241  // Pin
1242  nod_pt->pin(i);
1243  }
1244  // Else it's a pressure value
1245  else
1246  {
1247  // Exclude non-nodal pressure based elements
1248  if (this->p_nodal_index_nst()>=0)
1249  {
1250  // Backup
1251  eqn_number_backup[nod_pt][i]=nod_pt->eqn_number(i);
1252  }
1253  }
1254  }
1255 
1256 
1257  // If it's a solid node deal with its positional data too
1258  SolidNode* solid_nod_pt=dynamic_cast<SolidNode*>(nod_pt);
1259  if (solid_nod_pt!=0)
1260  {
1261  Data* solid_posn_data_pt=solid_nod_pt->variable_position_pt();
1262  if (eqn_number_backup[solid_posn_data_pt].size()==0)
1263  {
1264  unsigned nvalue=solid_posn_data_pt->nvalue();
1265  eqn_number_backup[solid_posn_data_pt].resize(nvalue);
1266  for (unsigned i=0;i<nvalue;i++)
1267  {
1268  // Backup
1269  eqn_number_backup[solid_posn_data_pt][i]=
1270  solid_posn_data_pt->eqn_number(i);
1271 
1272  // Pin
1273  solid_posn_data_pt->pin(i);
1274  }
1275  }
1276  }
1277  }
1278  }
1279  }
1280 
1281 
1282  /// \short Build FaceElements that apply the Robin boundary condition
1283  /// to the pressure advection diffusion problem required by
1284  /// Fp preconditioner
1285  virtual void build_fp_press_adv_diff_robin_bc_element(const unsigned&
1286  face_index)=0;
1287 
1288  /// \short Output the FaceElements that apply the Robin boundary condition
1289  /// to the pressure advection diffusion problem required by
1290  /// Fp preconditioner
1292  {
1294  for (unsigned e=0;e<nel;e++)
1295  {
1297  outfile << "ZONE" << std::endl;
1298  Vector<double> unit_normal(DIM);
1299  Vector<double> x(DIM);
1300  Vector<double> s(DIM-1);
1301  unsigned n=face_el_pt->integral_pt()->nweight();
1302  for (unsigned ipt=0;ipt<n;ipt++)
1303  {
1304  for(unsigned i=0;i<DIM-1;i++)
1305  {
1306  s[i]=face_el_pt->integral_pt()->knot(ipt,i);
1307  }
1308  face_el_pt->interpolated_x(s,x);
1309  face_el_pt->outer_unit_normal(ipt,unit_normal);
1310  for (unsigned i=0;i<DIM;i++)
1311  {
1312  outfile << x[i] << " ";
1313  }
1314  for (unsigned i=0;i<DIM;i++)
1315  {
1316  outfile << unit_normal[i] << " ";
1317  }
1318  outfile << std::endl;
1319  }
1320  }
1321  }
1322 
1323  /// \short Delete the FaceElements that apply the Robin boundary condition
1324  /// to the pressure advection diffusion problem required by
1325  /// Fp preconditioner
1327  {
1329  for (unsigned e=0;e<nel;e++)
1330  {
1332  }
1334  }
1335 
1336  /// \short Compute derivatives of elemental residual vector with respect
1337  /// to nodal coordinates. Overwrites default implementation in
1338  /// FiniteElement base class.
1339  /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
1341  dresidual_dnodal_coordinates);
1342 
1343 
1344 
1345  /// Compute vector of FE interpolated velocity u at local coordinate s
1347  {
1348  //Find number of nodes
1349  unsigned n_node = nnode();
1350  //Local shape function
1351  Shape psi(n_node);
1352  //Find values of shape function
1353  shape(s,psi);
1354 
1355  for (unsigned i=0;i<DIM;i++)
1356  {
1357  //Index at which the nodal value is stored
1358  unsigned u_nodal_index = u_index_nst(i);
1359  //Initialise value of u
1360  veloc[i] = 0.0;
1361  //Loop over the local nodes and sum
1362  for(unsigned l=0;l<n_node;l++)
1363  {
1364  veloc[i] += nodal_value(l,u_nodal_index)*psi[l];
1365  }
1366  }
1367  }
1368 
1369  /// Return FE interpolated velocity u[i] at local coordinate s
1370  double interpolated_u_nst(const Vector<double> &s, const unsigned &i) const
1371  {
1372  //Find number of nodes
1373  unsigned n_node = nnode();
1374  //Local shape function
1375  Shape psi(n_node);
1376  //Find values of shape function
1377  shape(s,psi);
1378 
1379  //Get nodal index at which i-th velocity is stored
1380  unsigned u_nodal_index = u_index_nst(i);
1381 
1382  //Initialise value of u
1383  double interpolated_u = 0.0;
1384  //Loop over the local nodes and sum
1385  for(unsigned l=0;l<n_node;l++)
1386  {
1387  interpolated_u += nodal_value(l,u_nodal_index)*psi[l];
1388  }
1389 
1390  return(interpolated_u);
1391  }
1392 
1393  /// \short Return FE interpolated velocity u[i] at local coordinate s
1394  /// at time level t (t=0: present; t>0: history)
1395  double interpolated_u_nst(const unsigned& t,
1396  const Vector<double> &s,
1397  const unsigned &i) const
1398  {
1399  //Find number of nodes
1400  unsigned n_node = nnode();
1401 
1402  //Local shape function
1403  Shape psi(n_node);
1404 
1405  //Find values of shape function
1406  shape(s,psi);
1407 
1408  //Get nodal index at which i-th velocity is stored
1409  unsigned u_nodal_index = u_index_nst(i);
1410 
1411  //Initialise value of u
1412  double interpolated_u = 0.0;
1413  //Loop over the local nodes and sum
1414  for(unsigned l=0;l<n_node;l++)
1415  {
1416  interpolated_u += nodal_value(t,l,u_nodal_index)*psi[l];
1417  }
1418 
1419  return(interpolated_u);
1420  }
1421 
1422  /// \short Compute the derivatives of the i-th component of
1423  /// velocity at point s with respect
1424  /// to all data that can affect its value. In addition, return the global
1425  /// equation numbers corresponding to the data. The function is virtual
1426  /// so that it can be overloaded in the refineable version
1428  const unsigned &i,
1429  Vector<double> &du_ddata,
1430  Vector<unsigned> &global_eqn_number)
1431  {
1432  //Find number of nodes
1433  unsigned n_node = nnode();
1434  //Local shape function
1435  Shape psi(n_node);
1436  //Find values of shape function
1437  shape(s,psi);
1438 
1439  //Find the index at which the velocity component is stored
1440  const unsigned u_nodal_index = u_index_nst(i);
1441 
1442  //Find the number of dofs associated with interpolated u
1443  unsigned n_u_dof=0;
1444  for(unsigned l=0;l<n_node;l++)
1445  {
1446  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
1447  //If it's positive add to the count
1448  if(global_eqn >= 0) {++n_u_dof;}
1449  }
1450 
1451  //Now resize the storage schemes
1452  du_ddata.resize(n_u_dof,0.0);
1453  global_eqn_number.resize(n_u_dof,0);
1454 
1455  //Loop over th nodes again and set the derivatives
1456  unsigned count=0;
1457  //Loop over the local nodes and sum
1458  for(unsigned l=0;l<n_node;l++)
1459  {
1460  //Get the global equation number
1461  int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
1462  if (global_eqn >= 0)
1463  {
1464  //Set the global equation number
1465  global_eqn_number[count] = global_eqn;
1466  //Set the derivative with respect to the unknown
1467  du_ddata[count] = psi[l];
1468  //Increase the counter
1469  ++count;
1470  }
1471  }
1472  }
1473 
1474 
1475  /// Return FE interpolated pressure at local coordinate s
1476  virtual double interpolated_p_nst(const Vector<double> &s) const
1477  {
1478  //Find number of nodes
1479  unsigned n_pres = npres_nst();
1480  //Local shape function
1481  Shape psi(n_pres);
1482  //Find values of shape function
1483  pshape_nst(s,psi);
1484 
1485  //Initialise value of p
1486  double interpolated_p = 0.0;
1487  //Loop over the local nodes and sum
1488  for(unsigned l=0;l<n_pres;l++)
1489  {
1490  interpolated_p += p_nst(l)*psi[l];
1491  }
1492 
1493  return(interpolated_p);
1494  }
1495 
1496 
1497  /// Return FE interpolated pressure at local coordinate s at time level t
1498  double interpolated_p_nst(const unsigned &t, const Vector<double> &s) const
1499  {
1500  //Find number of nodes
1501  unsigned n_pres = npres_nst();
1502  //Local shape function
1503  Shape psi(n_pres);
1504  //Find values of shape function
1505  pshape_nst(s,psi);
1506 
1507  //Initialise value of p
1508  double interpolated_p = 0.0;
1509  //Loop over the local nodes and sum
1510  for(unsigned l=0;l<n_pres;l++)
1511  {
1512  interpolated_p += p_nst(t,l)*psi[l];
1513  }
1514 
1515  return(interpolated_p);
1516  }
1517 
1518 
1519  /// \short Output solution in data vector at local cordinates s:
1520  /// x,y [,z], u,v,[w], p
1522  {
1523  // Dimension
1524  unsigned dim=s.size();
1525 
1526  // Resize data for values
1527  data.resize(2*dim+1);
1528 
1529  // Write values in the vector
1530  for (unsigned i=0; i<dim; i++)
1531  {
1532  data[i]=interpolated_x(s,i);
1533  data[i+dim]=this->interpolated_u_nst(s,i);
1534  }
1535  data[2*dim]=this->interpolated_p_nst(s);
1536  }
1537 
1538 };
1539 
1540 //////////////////////////////////////////////////////////////////////////////
1541 //////////////////////////////////////////////////////////////////////////////
1542 //////////////////////////////////////////////////////////////////////////////
1543 
1544 
1545 //==========================================================================
1546 /// Crouzeix_Raviart elements are Navier--Stokes elements with quadratic
1547 /// interpolation for velocities and positions, but a discontinuous linear
1548 /// pressure interpolation. They can be used within oomph-lib's
1549 /// block preconditioning framework.
1550 //==========================================================================
1551 template <unsigned DIM>
1552 class QCrouzeixRaviartElement : public virtual QElement<DIM,3>,
1553  public virtual NavierStokesEquations<DIM>
1554 {
1555  private:
1556 
1557  /// Static array of ints to hold required number of variables at nodes
1558  static const unsigned Initial_Nvalue[];
1559 
1560  protected:
1561 
1562  /// Internal index that indicates at which internal data the pressure
1563  /// is stored
1565 
1566 
1567  /// \short Velocity shape and test functions and their derivs
1568  /// w.r.t. to global coords at local coordinate s (taken from geometry)
1569  ///Return Jacobian of mapping between local and global coordinates.
1570  inline double dshape_and_dtest_eulerian_nst(const Vector<double> &s,
1571  Shape &psi,
1572  DShape &dpsidx,
1573  Shape &test,
1574  DShape &dtestdx) const;
1575 
1576  /// \short Velocity shape and test functions and their derivs
1577  /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
1578  ///Return Jacobian of mapping between local and global coordinates.
1579  inline double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt,
1580  Shape &psi,
1581  DShape &dpsidx,
1582  Shape &test,
1583  DShape &dtestdx) const;
1584 
1585  /// \short Shape/test functions and derivs w.r.t. to global coords at
1586  /// integration point ipt; return Jacobian of mapping (J). Also compute
1587  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
1589  const unsigned &ipt,
1590  Shape &psi,
1591  DShape &dpsidx,
1592  RankFourTensor<double> &d_dpsidx_dX,
1593  Shape &test,
1594  DShape &dtestdx,
1595  RankFourTensor<double> &d_dtestdx_dX,
1596  DenseMatrix<double> &djacobian_dX) const;
1597 
1598 
1599 public:
1600 
1601  /// Constructor, there are DIM+1 internal values (for the pressure)
1603  {
1604  //Allocate and add one Internal data object that stored DIM+1 pressure
1605  //values;
1606  P_nst_internal_index = this->add_internal_data(new Data(DIM+1));
1607  }
1608 
1609  /// \short Number of values (pinned or dofs) required at local node n.
1610  virtual unsigned required_nvalue(const unsigned &n) const;
1611 
1612 
1613  /// Pressure shape functions at local coordinate s
1614  inline void pshape_nst(const Vector<double> &s, Shape &psi) const;
1615 
1616  /// Pressure shape and test functions at local coordinte s
1617  inline void pshape_nst(const Vector<double> &s, Shape &psi,
1618  Shape &test) const;
1619 
1620  /// \short Return the i-th pressure value
1621  /// (Discontinous pressure interpolation -- no need to cater for hanging
1622  /// nodes).
1623  double p_nst(const unsigned &i) const
1624  {return this->internal_data_pt(P_nst_internal_index)->value(i);}
1625 
1626  /// \short Return the i-th pressure value
1627  /// (Discontinous pressure interpolation -- no need to cater for hanging
1628  /// nodes).
1629  double p_nst(const unsigned &t, const unsigned &i) const
1630  {return this->internal_data_pt(P_nst_internal_index)->value(t,i);}
1631 
1632  /// Return number of pressure values
1633  unsigned npres_nst() const {return DIM+1;}
1634 
1635  /// \short Pressure shape and test functions and their derivs
1636  /// w.r.t. to global coords at local coordinate s (taken from geometry)
1637  /// Return Jacobian of mapping between local and global coordinates.
1638  inline double dpshape_and_dptest_eulerian_nst(const Vector<double> &s,
1639  Shape &ppsi,
1640  DShape &dppsidx,
1641  Shape &ptest,
1642  DShape &dptestdx) const;
1643 
1644  /// Return the local equation numbers for the pressure values.
1645  inline int p_local_eqn(const unsigned &n) const
1646  {return this->internal_local_eqn(P_nst_internal_index,n);}
1647 
1648  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
1649  void fix_pressure(const unsigned &p_dof, const double &p_value)
1650  {
1651  this->internal_data_pt(P_nst_internal_index)->pin(p_dof);
1652  this->internal_data_pt(P_nst_internal_index)->set_value(p_dof,p_value);
1653  }
1654 
1655 
1656  /// \short Build FaceElements that apply the Robin boundary condition
1657  /// to the pressure advection diffusion problem required by
1658  /// Fp preconditioner
1660  face_index)
1661  {
1664  this, face_index));
1665  }
1666 
1667  /// \short Add to the set \c paired_load_data pairs containing
1668  /// - the pointer to a Data object
1669  /// and
1670  /// - the index of the value in that Data object
1671  /// .
1672  /// for all values (pressures, velocities) that affect the
1673  /// load computed in the \c get_load(...) function.
1674  void identify_load_data(
1675  std::set<std::pair<Data*,unsigned> > &paired_load_data);
1676 
1677  /// \short Add to the set \c paired_pressure_data pairs
1678  /// containing
1679  /// - the pointer to a Data object
1680  /// and
1681  /// - the index of the value in that Data object
1682  /// .
1683  /// for all pressure values that affect the
1684  /// load computed in the \c get_load(...) function.
1686  std::set<std::pair<Data*,unsigned> > &paired_pressure_data);
1687 
1688 
1689  /// Redirect output to NavierStokesEquations output
1690  void output(std::ostream &outfile)
1692 
1693  /// Redirect output to NavierStokesEquations output
1694  void output(std::ostream &outfile, const unsigned &nplot)
1695  {NavierStokesEquations<DIM>::output(outfile,nplot);}
1696 
1697 
1698  /// Redirect output to NavierStokesEquations output
1699  void output(FILE* file_pt) {NavierStokesEquations<DIM>::output(file_pt);}
1700 
1701  /// Redirect output to NavierStokesEquations output
1702  void output(FILE* file_pt, const unsigned &nplot)
1703  {NavierStokesEquations<DIM>::output(file_pt,nplot);}
1704 
1705 
1706  /// \short Full output function:
1707  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
1708  /// in tecplot format. Default number of plot points
1709  void full_output(std::ostream &outfile)
1711 
1712  /// \short Full output function:
1713  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
1714  /// in tecplot format. nplot points in each coordinate direction
1715  void full_output(std::ostream &outfile, const unsigned &nplot)
1716  {NavierStokesEquations<DIM>::full_output(outfile,nplot);}
1717 
1718 
1719  /// \short The number of "DOF types" that degrees of freedom in this element
1720  /// are sub-divided into: Velocity and pressure.
1721  unsigned ndof_types() const
1722  {
1723  return DIM+1;
1724  }
1725 
1726  /// \short Create a list of pairs for all unknowns in this element,
1727  /// so that the first entry in each pair contains the global equation
1728  /// number of the unknown, while the second one contains the number
1729  /// of the "DOF type" that this unknown is associated with.
1730  /// (Function can obviously only be called if the equation numbering
1731  /// scheme has been set up.) Velocity=0; Pressure=1
1733  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const;
1734 
1735 };
1736 
1737 //Inline functions
1738 
1739 //=======================================================================
1740 /// Derivatives of the shape functions and test functions w.r.t. to global
1741 /// (Eulerian) coordinates. Return Jacobian of mapping between
1742 /// local and global coordinates.
1743 //=======================================================================
1744 template<unsigned DIM>
1746  const Vector<double> &s, Shape &psi,
1747  DShape &dpsidx, Shape &test,
1748  DShape &dtestdx) const
1749 {
1750  //Call the geometrical shape functions and derivatives
1751  double J = this->dshape_eulerian(s,psi,dpsidx);
1752  //The test functions are equal to the shape functions
1753  test = psi;
1754  dtestdx = dpsidx;
1755  //Return the jacobian
1756  return J;
1757 }
1758 
1759 //=======================================================================
1760 /// Derivatives of the shape functions and test functions w.r.t. to global
1761 /// (Eulerian) coordinates. Return Jacobian of mapping between
1762 /// local and global coordinates.
1763 //=======================================================================
1764 template<unsigned DIM>
1765 inline double QCrouzeixRaviartElement<DIM>::
1767  const unsigned &ipt, Shape &psi,
1768  DShape &dpsidx, Shape &test,
1769  DShape &dtestdx) const
1770 {
1771  //Call the geometrical shape functions and derivatives
1772  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
1773  //The test functions are equal to the shape functions
1774  test = psi;
1775  dtestdx = dpsidx;
1776  //Return the jacobian
1777  return J;
1778 }
1779 
1780 
1781 //=======================================================================
1782 /// 2D
1783 /// Define the shape functions (psi) and test functions (test) and
1784 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
1785 /// and return Jacobian of mapping (J). Additionally compute the
1786 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
1787 ///
1788 /// Galerkin: Test functions = shape functions
1789 //=======================================================================
1790 template<>
1791 inline double QCrouzeixRaviartElement<2>::
1793  const unsigned &ipt, Shape &psi, DShape &dpsidx,
1794  RankFourTensor<double> &d_dpsidx_dX,
1795  Shape &test, DShape &dtestdx,
1796  RankFourTensor<double> &d_dtestdx_dX,
1797  DenseMatrix<double> &djacobian_dX) const
1798  {
1799  // Call the geometrical shape functions and derivatives
1800  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
1801  djacobian_dX,
1802  d_dpsidx_dX);
1803 
1804  // Loop over the test functions and derivatives and set them equal to the
1805  // shape functions
1806  for(unsigned i=0;i<9;i++)
1807  {
1808  test[i] = psi[i];
1809 
1810  for(unsigned k=0;k<2;k++)
1811  {
1812  dtestdx(i,k) = dpsidx(i,k);
1813 
1814  for(unsigned p=0;p<2;p++)
1815  {
1816  for(unsigned q=0;q<9;q++)
1817  {
1818  d_dtestdx_dX(p,q,i,k) = d_dpsidx_dX(p,q,i,k);
1819  }
1820  }
1821  }
1822  }
1823 
1824  // Return the jacobian
1825  return J;
1826  }
1827 
1828 
1829 //=======================================================================
1830 /// 3D
1831 /// Define the shape functions (psi) and test functions (test) and
1832 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
1833 /// and return Jacobian of mapping (J). Additionally compute the
1834 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
1835 ///
1836 /// Galerkin: Test functions = shape functions
1837 //=======================================================================
1838 template<>
1839 inline double QCrouzeixRaviartElement<3>::
1841  const unsigned &ipt, Shape &psi, DShape &dpsidx,
1842  RankFourTensor<double> &d_dpsidx_dX,
1843  Shape &test, DShape &dtestdx,
1844  RankFourTensor<double> &d_dtestdx_dX,
1845  DenseMatrix<double> &djacobian_dX) const
1846  {
1847  // Call the geometrical shape functions and derivatives
1848  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
1849  djacobian_dX,
1850  d_dpsidx_dX);
1851 
1852  // Loop over the test functions and derivatives and set them equal to the
1853  // shape functions
1854  for(unsigned i=0;i<27;i++)
1855  {
1856  test[i] = psi[i];
1857 
1858  for(unsigned k=0;k<3;k++)
1859  {
1860  dtestdx(i,k) = dpsidx(i,k);
1861 
1862  for(unsigned p=0;p<3;p++)
1863  {
1864  for(unsigned q=0;q<27;q++)
1865  {
1866  d_dtestdx_dX(p,q,i,k) = d_dpsidx_dX(p,q,i,k);
1867  }
1868  }
1869  }
1870  }
1871 
1872  // Return the jacobian
1873  return J;
1874  }
1875 
1876 
1877 
1878 //=======================================================================
1879 /// 2D :
1880 /// Pressure shape functions
1881 //=======================================================================
1882 template<>
1884  Shape &psi)
1885 const
1886 {
1887  psi[0] = 1.0;
1888  psi[1] = s[0];
1889  psi[2] = s[1];
1890 }
1891 
1892 
1893 
1894 //==========================================================================
1895 /// 2D :
1896 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
1897 /// Return Jacobian of mapping between local and global coordinates.
1898 //==========================================================================
1899 template<>
1901  const Vector<double> &s,
1902  Shape &ppsi,
1903  DShape &dppsidx,
1904  Shape &ptest,
1905  DShape &dptestdx) const
1906  {
1907 
1908  // Initalise with shape fcts and derivs. w.r.t. to local coordinates
1909  ppsi[0] = 1.0;
1910  ppsi[1] = s[0];
1911  ppsi[2] = s[1];
1912 
1913  dppsidx(0,0) = 0.0;
1914  dppsidx(1,0) = 1.0;
1915  dppsidx(2,0) = 0.0;
1916 
1917  dppsidx(0,1) = 0.0;
1918  dppsidx(1,1) = 0.0;
1919  dppsidx(2,1) = 1.0;
1920 
1921 
1922  //Get the values of the shape functions and their local derivatives
1923  Shape psi(9);
1924  DShape dpsi(9,2);
1925  dshape_local(s,psi,dpsi);
1926 
1927  //Allocate memory for the inverse 2x2 jacobian
1928  DenseMatrix<double> inverse_jacobian(2);
1929 
1930  //Now calculate the inverse jacobian
1931  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
1932 
1933  //Now set the values of the derivatives to be derivs w.r.t. to the
1934  // Eulerian coordinates
1935  transform_derivatives(inverse_jacobian,dppsidx);
1936 
1937  //The test functions are equal to the shape functions
1938  ptest = ppsi;
1939  dptestdx = dppsidx;
1940 
1941  //Return the determinant of the jacobian
1942  return det;
1943 
1944  }
1945 
1946 
1947 //=======================================================================
1948 /// Ppressure shape and test functions
1949 //=======================================================================
1950 template<unsigned DIM>
1951  inline void QCrouzeixRaviartElement<DIM>::
1953  Shape &psi,
1954  Shape &test) const
1955 {
1956  //Call the pressure shape functions
1957  this->pshape_nst(s,psi);
1958  //Test functions are equal to shape functions
1959  test = psi;
1960 }
1961 
1962 
1963 //=======================================================================
1964 /// 3D :
1965 /// Pressure shape functions
1966 //=======================================================================
1967 template<>
1969  Shape &psi)
1970 const
1971 {
1972  psi[0] = 1.0;
1973  psi[1] = s[0];
1974  psi[2] = s[1];
1975  psi[3] = s[2];
1976 }
1977 
1978 
1979 //==========================================================================
1980 /// 3D :
1981 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
1982 /// Return Jacobian of mapping between local and global coordinates.
1983 //==========================================================================
1984 template<>
1986  const Vector<double> &s,
1987  Shape &ppsi,
1988  DShape &dppsidx,
1989  Shape &ptest,
1990  DShape &dptestdx) const
1991  {
1992 
1993  // Initalise with shape fcts and derivs. w.r.t. to local coordinates
1994  ppsi[0] = 1.0;
1995  ppsi[1] = s[0];
1996  ppsi[2] = s[1];
1997  ppsi[3] = s[2];
1998 
1999  dppsidx(0,0) = 0.0;
2000  dppsidx(1,0) = 1.0;
2001  dppsidx(2,0) = 0.0;
2002  dppsidx(3,0) = 0.0;
2003 
2004  dppsidx(0,1) = 0.0;
2005  dppsidx(1,1) = 0.0;
2006  dppsidx(2,1) = 1.0;
2007  dppsidx(3,1) = 0.0;
2008 
2009  dppsidx(0,2) = 0.0;
2010  dppsidx(1,2) = 0.0;
2011  dppsidx(2,2) = 0.0;
2012  dppsidx(3,2) = 1.0;
2013 
2014 
2015  //Get the values of the shape functions and their local derivatives
2016  Shape psi(27);
2017  DShape dpsi(27,3);
2018  dshape_local(s,psi,dpsi);
2019 
2020  // Allocate memory for the inverse 3x3 jacobian
2021  DenseMatrix<double> inverse_jacobian(3);
2022 
2023  // Now calculate the inverse jacobian
2024  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
2025 
2026  // Now set the values of the derivatives to be derivs w.r.t. to the
2027  // Eulerian coordinates
2028  transform_derivatives(inverse_jacobian,dppsidx);
2029 
2030  //The test functions are equal to the shape functions
2031  ptest = ppsi;
2032  dptestdx = dppsidx;
2033 
2034  // Return the determinant of the jacobian
2035  return det;
2036 
2037  }
2038 
2039 
2040 //=======================================================================
2041 /// Face geometry of the 2D Crouzeix_Raviart elements
2042 //=======================================================================
2043 template<>
2044 class FaceGeometry<QCrouzeixRaviartElement<2> >: public virtual QElement<1,3>
2045 {
2046  public:
2047  FaceGeometry() : QElement<1,3>() {}
2048 };
2049 
2050 //=======================================================================
2051 /// Face geometry of the 3D Crouzeix_Raviart elements
2052 //=======================================================================
2053 template<>
2054 class FaceGeometry<QCrouzeixRaviartElement<3> >: public virtual QElement<2,3>
2055 {
2056 
2057  public:
2058  FaceGeometry() : QElement<2,3>() {}
2059 };
2060 
2061 //=======================================================================
2062 /// Face geometry of the FaceGeometry of the 2D Crouzeix_Raviart elements
2063 //=======================================================================
2064 template<>
2066 public virtual PointElement
2067 {
2068  public:
2070 };
2071 
2072 
2073 //=======================================================================
2074 /// Face geometry of the FaceGeometry of the 3D Crouzeix_Raviart elements
2075 //=======================================================================
2076 template<>
2078 public virtual QElement<1,3>
2079 {
2080  public:
2081  FaceGeometry() : QElement<1,3>() {}
2082 };
2083 
2084 
2085 
2086 ////////////////////////////////////////////////////////////////////////////
2087 ////////////////////////////////////////////////////////////////////////////
2088 
2089 
2090 //=======================================================================
2091 /// Taylor--Hood elements are Navier--Stokes elements
2092 /// with quadratic interpolation for velocities and positions and
2093 /// continuous linear pressure interpolation. They can be used
2094 /// within oomph-lib's block-preconditioning framework.
2095 //=======================================================================
2096 template <unsigned DIM>
2097 class QTaylorHoodElement : public virtual QElement<DIM,3>,
2098  public virtual NavierStokesEquations<DIM>
2099 {
2100  private:
2101 
2102  /// Static array of ints to hold number of variables at node
2103  static const unsigned Initial_Nvalue[];
2104 
2105  protected:
2106 
2107  /// \short Static array of ints to hold conversion from pressure
2108  /// node numbers to actual node numbers
2109  static const unsigned Pconv[];
2110 
2111  /// \short Velocity shape and test functions and their derivs
2112  /// w.r.t. to global coords at local coordinate s (taken from geometry)
2113  /// Return Jacobian of mapping between local and global coordinates.
2114  inline double dshape_and_dtest_eulerian_nst(const Vector<double> &s,
2115  Shape &psi,
2116  DShape &dpsidx, Shape &test,
2117  DShape &dtestdx) const;
2118 
2119  /// \short Velocity shape and test functions and their derivs
2120  /// w.r.t. to global coords at local coordinate s (taken from geometry)
2121  /// Return Jacobian of mapping between local and global coordinates.
2122  inline double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt,
2123  Shape &psi,
2124  DShape &dpsidx,
2125  Shape &test,
2126  DShape &dtestdx) const;
2127 
2128  /// \short Shape/test functions and derivs w.r.t. to global coords at
2129  /// integration point ipt; return Jacobian of mapping (J). Also compute
2130  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
2132  const unsigned &ipt,
2133  Shape &psi,
2134  DShape &dpsidx,
2135  RankFourTensor<double> &d_dpsidx_dX,
2136  Shape &test,
2137  DShape &dtestdx,
2138  RankFourTensor<double> &d_dtestdx_dX,
2139  DenseMatrix<double> &djacobian_dX) const;
2140 
2141  public:
2142 
2143  /// Constructor, no internal data points
2145 
2146  /// \short Number of values (pinned or dofs) required at node n. Can
2147  /// be overwritten for hanging node version
2148  inline virtual unsigned required_nvalue(const unsigned &n) const
2149  {return Initial_Nvalue[n];}
2150 
2151 
2152  /// Pressure shape functions at local coordinate s
2153  inline void pshape_nst(const Vector<double> &s, Shape &psi) const;
2154 
2155  /// Pressure shape and test functions at local coordinte s
2156  inline void pshape_nst(const Vector<double> &s, Shape &psi,
2157  Shape &test) const;
2158 
2159  /// \short Set the value at which the pressure is stored in the nodes
2160  virtual int p_nodal_index_nst() const {return static_cast<int>(DIM);}
2161 
2162  /// Return the local equation numbers for the pressure values.
2163  inline int p_local_eqn(const unsigned &n) const
2164  {return this->nodal_local_eqn(Pconv[n],p_nodal_index_nst());}
2165 
2166  /// \short Access function for the pressure values at local pressure
2167  /// node n_p (const version)
2168  double p_nst(const unsigned &n_p) const
2169  {return this->nodal_value(Pconv[n_p],this->p_nodal_index_nst());}
2170 
2171  /// \short Access function for the pressure values at local pressure
2172  /// node n_p (const version)
2173  double p_nst(const unsigned &t, const unsigned &n_p) const
2174  {return this->nodal_value(t,Pconv[n_p],this->p_nodal_index_nst());}
2175 
2176  /// \short Pressure shape and test functions and their derivs
2177  /// w.r.t. to global coords at local coordinate s (taken from geometry).
2178  /// Return Jacobian of mapping between local and global coordinates.
2179  inline double dpshape_and_dptest_eulerian_nst(const Vector<double> &s,
2180  Shape &ppsi,
2181  DShape &dppsidx,
2182  Shape &ptest,
2183  DShape &dptestdx) const;
2184 
2185  /// Return number of pressure values
2186  unsigned npres_nst() const
2187  {return static_cast<unsigned>(pow(2.0,static_cast<int>(DIM)));}
2188 
2189  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
2190  void fix_pressure(const unsigned &p_dof, const double &p_value)
2191  {
2192  this->node_pt(Pconv[p_dof])->pin(this->p_nodal_index_nst());
2193  this->node_pt(Pconv[p_dof])->set_value(this->p_nodal_index_nst(),p_value);
2194  }
2195 
2196 
2197 
2198  /// \short Build FaceElements that apply the Robin boundary condition
2199  /// to the pressure advection diffusion problem required by
2200  /// Fp preconditioner
2202  face_index)
2203  {
2206  this, face_index));
2207  }
2208 
2209 
2210  /// \short Add to the set \c paired_load_data pairs containing
2211  /// - the pointer to a Data object
2212  /// and
2213  /// - the index of the value in that Data object
2214  /// .
2215  /// for all values (pressures, velocities) that affect the
2216  /// load computed in the \c get_load(...) function.
2217  void identify_load_data(
2218  std::set<std::pair<Data*,unsigned> > &paired_load_data);
2219 
2220 
2221  /// \short Add to the set \c paired_pressure_data pairs
2222  /// containing
2223  /// - the pointer to a Data object
2224  /// and
2225  /// - the index of the value in that Data object
2226  /// .
2227  /// for all pressure values that affect the
2228  /// load computed in the \c get_load(...) function.
2230  std::set<std::pair<Data*,unsigned> > &paired_pressure_data);
2231 
2232 
2233  /// Redirect output to NavierStokesEquations output
2234  void output(std::ostream &outfile)
2236 
2237  /// Redirect output to NavierStokesEquations output
2238  void output(std::ostream &outfile, const unsigned &nplot)
2239  {NavierStokesEquations<DIM>::output(outfile,nplot);}
2240 
2241  /// Redirect output to NavierStokesEquations output
2242  void output(FILE* file_pt) {NavierStokesEquations<DIM>::output(file_pt);}
2243 
2244  /// Redirect output to NavierStokesEquations output
2245  void output(FILE* file_pt, const unsigned &nplot)
2246  {NavierStokesEquations<DIM>::output(file_pt,nplot);}
2247 
2248 
2249  /// \short Returns the number of "DOF types" that degrees of freedom
2250  /// in this element are sub-divided into: Velocity and pressure.
2251  unsigned ndof_types() const
2252  {
2253  return DIM+1;
2254  }
2255 
2256  /// \short Create a list of pairs for all unknowns in this element,
2257  /// so that the first entry in each pair contains the global equation
2258  /// number of the unknown, while the second one contains the number
2259  /// of the "DOF type" that this unknown is associated with.
2260  /// (Function can obviously only be called if the equation numbering
2261  /// scheme has been set up.) Velocity=0; Pressure=1
2263  std::list<std::pair<unsigned long, unsigned> >& dof_lookup_list) const;
2264 
2265 
2266 };
2267 
2268 //Inline functions
2269 
2270 //==========================================================================
2271 /// Derivatives of the shape functions and test functions w.r.t to
2272 /// global (Eulerian) coordinates. Return Jacobian of mapping between
2273 /// local and global coordinates.
2274 //==========================================================================
2275 template<unsigned DIM>
2277  const Vector<double> &s,
2278  Shape &psi,
2279  DShape &dpsidx, Shape &test,
2280  DShape &dtestdx) const
2281 {
2282  //Call the geometrical shape functions and derivatives
2283  double J = this->dshape_eulerian(s,psi,dpsidx);
2284 
2285  //The test functions are equal to the shape functions
2286  test = psi;
2287  dtestdx = dpsidx;
2288 
2289  //Return the jacobian
2290  return J;
2291 }
2292 
2293 
2294 //==========================================================================
2295 /// Derivatives of the shape functions and test functions w.r.t to
2296 /// global (Eulerian) coordinates. Return Jacobian of mapping between
2297 /// local and global coordinates.
2298 //==========================================================================
2299 template<unsigned DIM>
2301  const unsigned &ipt,Shape &psi, DShape &dpsidx, Shape &test,
2302  DShape &dtestdx) const
2303 {
2304  //Call the geometrical shape functions and derivatives
2305  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
2306 
2307  //The test functions are equal to the shape functions
2308  test = psi;
2309  dtestdx = dpsidx;
2310 
2311  //Return the jacobian
2312  return J;
2313 }
2314 
2315 
2316 //==========================================================================
2317 /// 2D :
2318 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
2319 /// Return Jacobian of mapping between local and global coordinates.
2320 //==========================================================================
2321 template<>
2323  const Vector<double> &s,
2324  Shape &ppsi,
2325  DShape &dppsidx,
2326  Shape &ptest,
2327  DShape &dptestdx) const
2328  {
2329 
2330  //Local storage
2331  double psi1[2], psi2[2];
2332  double dpsi1[2],dpsi2[2];
2333 
2334  //Call the OneDimensional Shape functions
2335  OneDimLagrange::shape<2>(s[0],psi1);
2336  OneDimLagrange::shape<2>(s[1],psi2);
2337  OneDimLagrange::dshape<2>(s[0],dpsi1);
2338  OneDimLagrange::dshape<2>(s[1],dpsi2);
2339 
2340  //Now let's loop over the nodal points in the element
2341  //s1 is the "x" coordinate, s2 the "y"
2342  for(unsigned i=0;i<2;i++)
2343  {
2344  for(unsigned j=0;j<2;j++)
2345  {
2346  /*Multiply the two 1D functions together to get the 2D function*/
2347  ppsi[2*i+j] = psi2[i]*psi1[j];
2348  dppsidx(2*i+j,0)= psi2[i]*dpsi1[j];
2349  dppsidx(2*i+j,1)=dpsi2[i]* psi1[j];
2350  }
2351  }
2352 
2353 
2354  //Get the values of the shape functions and their local derivatives
2355  Shape psi(9);
2356  DShape dpsi(9,2);
2357  dshape_local(s,psi,dpsi);
2358 
2359  // Allocate memory for the inverse 2x2 jacobian
2360  DenseMatrix<double> inverse_jacobian(2);
2361 
2362  // Now calculate the inverse jacobian
2363  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
2364 
2365  // Now set the values of the derivatives to be derivs w.r.t. to the
2366  // Eulerian coordinates
2367  transform_derivatives(inverse_jacobian,dppsidx);
2368 
2369  //The test functions are equal to the shape functions
2370  ptest = ppsi;
2371  dptestdx = dppsidx;
2372 
2373  // Return the determinant of the jacobian
2374  return det;
2375 
2376  }
2377 
2378 
2379 
2380 //==========================================================================
2381 /// 2D :
2382 /// Define the shape functions (psi) and test functions (test) and
2383 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
2384 /// and return Jacobian of mapping (J). Additionally compute the
2385 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
2386 ///
2387 /// Galerkin: Test functions = shape functions
2388 //==========================================================================
2389 template<>
2390 inline double QTaylorHoodElement<2>::
2392  const unsigned &ipt, Shape &psi, DShape &dpsidx,
2393  RankFourTensor<double> &d_dpsidx_dX,
2394  Shape &test, DShape &dtestdx,
2395  RankFourTensor<double> &d_dtestdx_dX,
2396  DenseMatrix<double> &djacobian_dX) const
2397  {
2398  // Call the geometrical shape functions and derivatives
2399  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
2400  djacobian_dX,
2401  d_dpsidx_dX);
2402 
2403  // Loop over the test functions and derivatives and set them equal to the
2404  // shape functions
2405  for(unsigned i=0;i<9;i++)
2406  {
2407  test[i] = psi[i];
2408 
2409  for(unsigned k=0;k<2;k++)
2410  {
2411  dtestdx(i,k) = dpsidx(i,k);
2412 
2413  for(unsigned p=0;p<2;p++)
2414  {
2415  for(unsigned q=0;q<9;q++)
2416  {
2417  d_dtestdx_dX(p,q,i,k) = d_dpsidx_dX(p,q,i,k);
2418  }
2419  }
2420  }
2421  }
2422 
2423  // Return the jacobian
2424  return J;
2425  }
2426 
2427 
2428 //==========================================================================
2429 /// 3D :
2430 /// Define the shape functions (psi) and test functions (test) and
2431 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
2432 /// and return Jacobian of mapping (J). Additionally compute the
2433 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
2434 ///
2435 /// Galerkin: Test functions = shape functions
2436 //==========================================================================
2437 template<>
2438  inline double QTaylorHoodElement<3>::
2440  const unsigned &ipt, Shape &psi, DShape &dpsidx,
2441  RankFourTensor<double> &d_dpsidx_dX,
2442  Shape &test, DShape &dtestdx,
2443  RankFourTensor<double> &d_dtestdx_dX,
2444  DenseMatrix<double> &djacobian_dX) const
2445  {
2446  // Call the geometrical shape functions and derivatives
2447  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
2448  djacobian_dX,
2449  d_dpsidx_dX);
2450 
2451  // Loop over the test functions and derivatives and set them equal to the
2452  // shape functions
2453  for(unsigned i=0;i<27;i++)
2454  {
2455  test[i] = psi[i];
2456 
2457  for(unsigned k=0;k<3;k++)
2458  {
2459  dtestdx(i,k) = dpsidx(i,k);
2460 
2461  for(unsigned p=0;p<3;p++)
2462  {
2463  for(unsigned q=0;q<27;q++)
2464  {
2465  d_dtestdx_dX(p,q,i,k) = d_dpsidx_dX(p,q,i,k);
2466  }
2467  }
2468  }
2469  }
2470 
2471  // Return the jacobian
2472  return J;
2473  }
2474 
2475 
2476 
2477 //==========================================================================
2478 /// 2D :
2479 /// Pressure shape functions
2480 //==========================================================================
2481 template<>
2483  Shape &psi)
2484 const
2485 {
2486  //Local storage
2487  double psi1[2], psi2[2];
2488  //Call the OneDimensional Shape functions
2489  OneDimLagrange::shape<2>(s[0],psi1);
2490  OneDimLagrange::shape<2>(s[1],psi2);
2491 
2492  //Now let's loop over the nodal points in the element
2493  //s1 is the "x" coordinate, s2 the "y"
2494  for(unsigned i=0;i<2;i++)
2495  {
2496  for(unsigned j=0;j<2;j++)
2497  {
2498  /*Multiply the two 1D functions together to get the 2D function*/
2499  psi[2*i + j] = psi2[i]*psi1[j];
2500  }
2501  }
2502 }
2503 
2504 
2505 //==========================================================================
2506 /// 3D :
2507 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
2508 /// Return Jacobian of mapping between local and global coordinates.
2509 //==========================================================================
2510 template<>
2512  const Vector<double> &s,
2513  Shape &ppsi,
2514  DShape &dppsidx,
2515  Shape &ptest,
2516  DShape &dptestdx) const
2517  {
2518  //Local storage
2519  double psi1[2], psi2[2], psi3[2];
2520  double dpsi1[2],dpsi2[2],dpsi3[2];
2521 
2522  //Call the OneDimensional Shape functions
2523  OneDimLagrange::shape<2>(s[0],psi1);
2524  OneDimLagrange::shape<2>(s[1],psi2);
2525  OneDimLagrange::shape<2>(s[2],psi3);
2526  OneDimLagrange::dshape<2>(s[0],dpsi1);
2527  OneDimLagrange::dshape<2>(s[1],dpsi2);
2528  OneDimLagrange::dshape<2>(s[2],dpsi3);
2529 
2530  //Now let's loop over the nodal points in the element
2531  //s0 is the "x" coordinate, s1 the "y", s2 is the "z"
2532  for(unsigned i=0;i<2;i++)
2533  {
2534  for(unsigned j=0;j<2;j++)
2535  {
2536  for(unsigned k=0;k<2;k++)
2537  {
2538  /*Multiply the three 1D functions together to get the 3D function*/
2539  ppsi[4*i + 2*j + k] = psi3[i]*psi2[j]*psi1[k];
2540  dppsidx(4*i + 2*j + k,0) = psi3[i] * psi2[j] *dpsi1[k];
2541  dppsidx(4*i + 2*j + k,1) = psi3[i] *dpsi2[j] * psi1[k];
2542  dppsidx(4*i + 2*j + k,2) = dpsi3[i] * psi2[j] * psi1[k];
2543  }
2544  }
2545  }
2546 
2547 
2548  //Get the values of the shape functions and their local derivatives
2549  Shape psi(27);
2550  DShape dpsi(27,3);
2551  dshape_local(s,psi,dpsi);
2552 
2553  // Allocate memory for the inverse 3x3 jacobian
2554  DenseMatrix<double> inverse_jacobian(3);
2555 
2556  // Now calculate the inverse jacobian
2557  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
2558 
2559  // Now set the values of the derivatives to be derivs w.r.t. to the
2560  // Eulerian coordinates
2561  transform_derivatives(inverse_jacobian,dppsidx);
2562 
2563  //The test functions are equal to the shape functions
2564  ptest = ppsi;
2565  dptestdx = dppsidx;
2566 
2567  // Return the determinant of the jacobian
2568  return det;
2569 
2570 }
2571 
2572 //==========================================================================
2573 /// 3D :
2574 /// Pressure shape functions
2575 //==========================================================================
2576 template<>
2578  Shape &psi)
2579 const
2580 {
2581  //Local storage
2582  double psi1[2], psi2[2], psi3[2];
2583 
2584  //Call the OneDimensional Shape functions
2585  OneDimLagrange::shape<2>(s[0],psi1);
2586  OneDimLagrange::shape<2>(s[1],psi2);
2587  OneDimLagrange::shape<2>(s[2],psi3);
2588 
2589  //Now let's loop over the nodal points in the element
2590  //s0 is the "x" coordinate, s1 the "y", s2 is the "z"
2591  for(unsigned i=0;i<2;i++)
2592  {
2593  for(unsigned j=0;j<2;j++)
2594  {
2595  for(unsigned k=0;k<2;k++)
2596  {
2597  /*Multiply the three 1D functions together to get the 3D function*/
2598  psi[4*i + 2*j + k] = psi3[i]*psi2[j]*psi1[k];
2599  }
2600  }
2601  }
2602 }
2603 
2604 
2605 //==========================================================================
2606 /// Pressure shape and test functions
2607 //==========================================================================
2608 template<unsigned DIM>
2610  Shape &psi,
2611  Shape &test) const
2612 {
2613  //Call the pressure shape functions
2614  this->pshape_nst(s,psi);
2615  //Test functions are shape functions
2616  test = psi;
2617 }
2618 
2619 
2620 //=======================================================================
2621 /// Face geometry of the 2D Taylor_Hood elements
2622 //=======================================================================
2623 template<>
2624 class FaceGeometry<QTaylorHoodElement<2> >: public virtual QElement<1,3>
2625 {
2626  public:
2627  FaceGeometry() : QElement<1,3>() {}
2628 };
2629 
2630 //=======================================================================
2631 /// Face geometry of the 3D Taylor_Hood elements
2632 //=======================================================================
2633 template<>
2634 class FaceGeometry<QTaylorHoodElement<3> >: public virtual QElement<2,3>
2635 {
2636 
2637  public:
2638  FaceGeometry() : QElement<2,3>() {}
2639 };
2640 
2641 
2642 //=======================================================================
2643 /// Face geometry of the FaceGeometry of the 2D Taylor Hoodelements
2644 //=======================================================================
2645 template<>
2647 public virtual PointElement
2648 {
2649  public:
2651 };
2652 
2653 
2654 //=======================================================================
2655 /// Face geometry of the FaceGeometry of the 3D Taylor_Hood elements
2656 //=======================================================================
2657 template<>
2659 public virtual QElement<1,3>
2660 {
2661  public:
2662  FaceGeometry() : QElement<1,3>() {}
2663 };
2664 
2665 
2666 
2667 
2668 
2669 
2670 
2671 ////////////////////////////////////////////////////////////////////
2672 ////////////////////////////////////////////////////////////////////
2673 ////////////////////////////////////////////////////////////////////
2674 
2675 
2676 //==========================================================
2677 /// Taylor Hood upgraded to become projectable
2678 //==========================================================
2679  template<class TAYLOR_HOOD_ELEMENT>
2681  public virtual ProjectableElement<TAYLOR_HOOD_ELEMENT>
2682  {
2683 
2684  public:
2685 
2686  /// \short Constructor [this was only required explicitly
2687  /// from gcc 4.5.2 onwards...]
2689 
2690 
2691  /// \short Specify the values associated with field fld.
2692  /// The information is returned in a vector of pairs which comprise
2693  /// the Data object and the value within it, that correspond to field fld.
2694  /// In the underlying Taylor Hood elements the fld-th velocities are stored
2695  /// at the fld-th value of the nodes; the pressures (the dim-th
2696  /// field) are the dim-th values at the vertex nodes etc.
2698  {
2699  // Create the vector
2700  Vector<std::pair<Data*,unsigned> > data_values;
2701 
2702  // Velocities dofs
2703  if (fld<this->dim())
2704  {
2705  // Loop over all nodes
2706  unsigned nnod=this->nnode();
2707  for (unsigned j=0;j<nnod;j++)
2708  {
2709  // Add the data value associated with the velocity components
2710  data_values.push_back(std::make_pair(this->node_pt(j),fld));
2711  }
2712  }
2713  // Pressure
2714  else
2715  {
2716  // Loop over all vertex nodes
2717  unsigned Pconv_size=this->dim()+1;
2718  for (unsigned j=0;j<Pconv_size;j++)
2719  {
2720  // Add the data value associated with the pressure components
2721  unsigned vertex_index=this->Pconv[j];
2722  data_values.push_back(std::make_pair(this->node_pt(vertex_index),fld));
2723  }
2724  }
2725 
2726  // Return the vector
2727  return data_values;
2728 
2729  }
2730 
2731  /// \short Number of fields to be projected: dim+1, corresponding to
2732  /// velocity components and pressure
2734  {
2735  return this->dim()+1;
2736  }
2737 
2738  /// \short Number of history values to be stored for fld-th field. Whatever
2739  /// the timestepper has set up for the velocity components and
2740  /// none for the pressure field.
2741  /// (Note: count includes current value!)
2742  unsigned nhistory_values_for_projection(const unsigned &fld)
2743  {
2744  if (fld==this->dim())
2745  {
2746  //pressure doesn't have history values
2747  return this->node_pt(0)->ntstorage();//1;
2748  }
2749  else
2750  {
2751  return this->node_pt(0)->ntstorage();
2752  }
2753  }
2754 
2755  /// \short Number of positional history values
2756  /// (Note: count includes current value!)
2758  {
2759  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
2760  }
2761 
2762  /// \short Return Jacobian of mapping and shape functions of field fld
2763  /// at local coordinate s
2764  double jacobian_and_shape_of_field(const unsigned &fld,
2765  const Vector<double> &s,
2766  Shape &psi)
2767  {
2768  unsigned n_dim=this->dim();
2769  unsigned n_node=this->nnode();
2770 
2771  if (fld==n_dim)
2772  {
2773  //We are dealing with the pressure
2774  this->pshape_nst(s,psi);
2775 
2776  Shape psif(n_node),testf(n_node);
2777  DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
2778 
2779  //Domain Shape
2780  double J=this->dshape_and_dtest_eulerian_nst(s,psif,dpsifdx,
2781  testf,dtestfdx);
2782  return J;
2783  }
2784  else
2785  {
2786  Shape testf(n_node);
2787  DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
2788 
2789  //Domain Shape
2790  double J=this->dshape_and_dtest_eulerian_nst(s,psi,dpsifdx,
2791  testf,dtestfdx);
2792  return J;
2793  }
2794  }
2795 
2796 
2797 
2798  /// \short Return interpolated field fld at local coordinate s, at time level
2799  /// t (t=0: present; t>0: history values)
2800  double get_field(const unsigned &t,
2801  const unsigned &fld,
2802  const Vector<double>& s)
2803  {
2804  unsigned n_dim =this->dim();
2805  unsigned n_node=this->nnode();
2806 
2807  //If fld=n_dim, we deal with the pressure
2808  if (fld==n_dim)
2809  {
2810  return this->interpolated_p_nst(t,s);
2811  }
2812  // Velocity
2813  else
2814  {
2815  //Find the index at which the variable is stored
2816  unsigned u_nodal_index = this->u_index_nst(fld);
2817 
2818  //Local shape function
2819  Shape psi(n_node);
2820 
2821  //Find values of shape function
2822  this->shape(s,psi);
2823 
2824  //Initialise value of u
2825  double interpolated_u = 0.0;
2826 
2827  //Sum over the local nodes at that time
2828  for(unsigned l=0;l<n_node;l++)
2829  {
2830  interpolated_u += this->nodal_value(t,l,u_nodal_index)*psi[l];
2831  }
2832  return interpolated_u;
2833  }
2834  }
2835 
2836 
2837 
2838  ///Return number of values in field fld
2839  unsigned nvalue_of_field(const unsigned &fld)
2840  {
2841  if (fld==this->dim())
2842  {
2843  return this->npres_nst();
2844  }
2845  else
2846  {
2847  return this->nnode();
2848  }
2849  }
2850 
2851 
2852  ///Return local equation number of value j in field fld.
2853  int local_equation(const unsigned &fld,
2854  const unsigned &j)
2855  {
2856  if (fld==this->dim())
2857  {
2858  return this->p_local_eqn(j);
2859  }
2860  else
2861  {
2862  const unsigned u_nodal_index = this->u_index_nst(fld);
2863  return this->nodal_local_eqn(j,u_nodal_index);
2864  }
2865  }
2866 
2867  };
2868 
2869 
2870 //=======================================================================
2871 /// Face geometry for element is the same as that for the underlying
2872 /// wrapped element
2873 //=======================================================================
2874  template<class ELEMENT>
2876  : public virtual FaceGeometry<ELEMENT>
2877  {
2878  public:
2879  FaceGeometry() : FaceGeometry<ELEMENT>() {}
2880  };
2881 
2882 
2883 //=======================================================================
2884 /// Face geometry of the Face Geometry for element is the same as
2885 /// that for the underlying wrapped element
2886 //=======================================================================
2887  template<class ELEMENT>
2889  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
2890  {
2891  public:
2893  };
2894 
2895 
2896 //==========================================================
2897 /// Crouzeix Raviart upgraded to become projectable
2898 //==========================================================
2899  template<class CROUZEIX_RAVIART_ELEMENT>
2901  public virtual ProjectableElement<CROUZEIX_RAVIART_ELEMENT>
2902  {
2903 
2904  public:
2905 
2906  /// \short Constructor [this was only required explicitly
2907  /// from gcc 4.5.2 onwards...]
2909 
2910  /// \short Specify the values associated with field fld.
2911  /// The information is returned in a vector of pairs which comprise
2912  /// the Data object and the value within it, that correspond to field fld.
2913  /// In the underlying Crouzeix Raviart elements the
2914  /// fld-th velocities are stored
2915  /// at the fld-th value of the nodes; the pressures are stored internally
2917  {
2918  // Create the vector
2919  Vector<std::pair<Data*,unsigned> > data_values;
2920 
2921  // Velocities dofs
2922  if (fld < this->dim())
2923  {
2924  // Loop over all nodes
2925  const unsigned n_node=this->nnode();
2926  for (unsigned n=0;n<n_node;n++)
2927  {
2928  // Add the data value associated with the velocity components
2929  data_values.push_back(std::make_pair(this->node_pt(n),fld));
2930  }
2931  }
2932  // Pressure
2933  else
2934  {
2935  //Need to push back the internal data
2936  const unsigned n_press = this->npres_nst();
2937  //Loop over all pressure values
2938  for(unsigned j=0;j<n_press;j++)
2939  {
2940  data_values.push_back(
2941  std::make_pair(
2942  this->internal_data_pt(this->P_nst_internal_index),j));
2943  }
2944  }
2945 
2946  // Return the vector
2947  return data_values;
2948  }
2949 
2950  /// \short Number of fields to be projected: dim+1, corresponding to
2951  /// velocity components and pressure
2953  {
2954  return this->dim()+1;
2955  }
2956 
2957  /// \short Number of history values to be stored for fld-th field. Whatever
2958  /// the timestepper has set up for the velocity components and
2959  /// none for the pressure field.
2960  /// (Note: count includes current value!)
2961  unsigned nhistory_values_for_projection(const unsigned &fld)
2962  {
2963  if (fld==this->dim())
2964  {
2965  //pressure doesn't have history values
2966  return 1;
2967  }
2968  else
2969  {
2970  return this->node_pt(0)->ntstorage();
2971  }
2972  }
2973 
2974  ///\short Number of positional history values.
2975  /// (Note: count includes current value!)
2977  {
2978  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
2979  }
2980 
2981  /// \short Return Jacobian of mapping and shape functions of field fld
2982  /// at local coordinate s
2983  double jacobian_and_shape_of_field(const unsigned &fld,
2984  const Vector<double> &s,
2985  Shape &psi)
2986  {
2987  unsigned n_dim=this->dim();
2988  unsigned n_node=this->nnode();
2989 
2990  if (fld==n_dim)
2991  {
2992  //We are dealing with the pressure
2993  this->pshape_nst(s,psi);
2994 
2995  Shape psif(n_node),testf(n_node);
2996  DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
2997 
2998  //Domain Shape
2999  double J=this->dshape_and_dtest_eulerian_nst(s,psif,dpsifdx,
3000  testf,dtestfdx);
3001  return J;
3002  }
3003  else
3004  {
3005  Shape testf(n_node);
3006  DShape dpsifdx(n_node,n_dim), dtestfdx(n_node,n_dim);
3007 
3008  //Domain Shape
3009  double J=this->dshape_and_dtest_eulerian_nst(s,psi,dpsifdx,
3010  testf,dtestfdx);
3011  return J;
3012  }
3013  }
3014 
3015 
3016 
3017  /// \short Return interpolated field fld at local coordinate s, at time level
3018  /// t (t=0: present; t>0: history values)
3019  double get_field(const unsigned &t,
3020  const unsigned &fld,
3021  const Vector<double>& s)
3022  {
3023  unsigned n_dim =this->dim();
3024 
3025  //If fld=n_dim, we deal with the pressure
3026  if (fld==n_dim)
3027  {
3028  return this->interpolated_p_nst(s);
3029  }
3030  // Velocity
3031  else
3032  {
3033  return this->interpolated_u_nst(t,s,fld);
3034  }
3035  }
3036 
3037 
3038  ///Return number of values in field fld
3039  unsigned nvalue_of_field(const unsigned &fld)
3040  {
3041  if (fld==this->dim())
3042  {
3043  return this->npres_nst();
3044  }
3045  else
3046  {
3047  return this->nnode();
3048  }
3049  }
3050 
3051 
3052  ///Return local equation number of value j in field fld.
3053  int local_equation(const unsigned &fld,
3054  const unsigned &j)
3055  {
3056  if (fld==this->dim())
3057  {
3058  return this->p_local_eqn(j);
3059  }
3060  else
3061  {
3062  const unsigned u_nodal_index = this->u_index_nst(fld);
3063  return this->nodal_local_eqn(j,u_nodal_index);
3064  }
3065  }
3066 
3067  };
3068 
3069 
3070 //=======================================================================
3071 /// Face geometry for element is the same as that for the underlying
3072 /// wrapped element
3073 //=======================================================================
3074  template<class ELEMENT>
3076  : public virtual FaceGeometry<ELEMENT>
3077  {
3078  public:
3079  FaceGeometry() : FaceGeometry<ELEMENT>() {}
3080  };
3081 
3082 
3083 //=======================================================================
3084 /// Face geometry of the Face Geometry for element is the same as
3085 /// that for the underlying wrapped element
3086 //=======================================================================
3087  template<class ELEMENT>
3089  : public virtual FaceGeometry<FaceGeometry<ELEMENT> >
3090  {
3091  public:
3093  };
3094 
3095 
3096 }
3097 
3098 #endif
double dshape_and_dtest_eulerian_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...
virtual void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)=0
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed...
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
void output_pressure_advection_diffusion_robin_elements(std::ostream &outfile)
Output the FaceElements that apply the Robin boundary condition to the pressure advection diffusion p...
unsigned n_u_nst() const
Return the number of velocity components Used in FluidInterfaceElements.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:2938
FpPressureAdvDiffRobinBCElement(FiniteElement *const &element_pt, const int &face_index, const bool &called_from_refineable_constructor=false)
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 nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. Whatever the timestepper has set up for the v...
double p_nst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version) ...
Vector< FpPressureAdvDiffRobinBCElementBase * > Pressure_advection_diffusion_robin_element_pt
Storage for FaceElements that apply Robin BC for pressure adv diff equation used in Fp preconditioner...
double interpolated_p_nst(const unsigned &t, const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s at time level t.
virtual unsigned npres_nst() const =0
Function to return number of pressure degrees of freedom.
virtual void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)=0
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
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 nscalar_paraview() const
Number of scalars/fields output by this element. Reimplements broken virtual function in base class...
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...
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 ...
void dshape< 2 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:604
void output(FILE *file_pt, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
Taylor Hood upgraded to become projectable.
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing...
Definition: elements.h:2872
QTaylorHoodElement()
Constructor, no internal data points.
void fill_in_contribution_to_djacobian_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam)
Compute the element's residual Vector and the jacobian matrix Virtual function can be overloaded by h...
const double & re_invfr() const
Global inverse Froude number.
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:847
virtual void dinterpolated_u_nst_ddata(const Vector< double > &s, const unsigned &i, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
Compute the derivatives of the i-th component of velocity at point s with respect to all data that ca...
cstr elem_len * i
Definition: cfortran.h:607
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at ipt-th integation point...
unsigned ndof_types() const
Returns the number of "DOF types" that degrees of freedom in this element are sub-divided into: Veloc...
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
double interpolated_u_nst(const unsigned &t, const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s at time level t (t=0: present; t>0: histor...
void pin_all_non_pressure_dofs(std::map< Data *, std::vector< int > > &eqn_number_backup)
Pin all non-pressure dofs and backup eqn numbers.
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns just the residuals.
double(* NavierStokesSourceFctPt)(const double &time, const Vector< double > &x)
Function pointer to source function fct(t,x) x is a Vector!
NavierStokesPressureAdvDiffSourceFctPt & source_fct_for_pressure_adv_diff()
Access function for the source-function pointer for pressure advection diffusion (used for validation...
double * Re_pt
Pointer to global Reynolds number.
NavierStokesSourceFctPt & source_fct_pt()
Access function for the source-function pointer.
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
virtual void pshape_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
unsigned nfields_for_projection()
Number of fields to be projected: dim+1, corresponding to velocity components and pressure...
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
A general Finite Element class.
Definition: elements.h:1271
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
virtual void get_body_force_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(* NavierStokesPressureAdvDiffSourceFctPt)(const Vector< double > &x)
Function pointer to source function fct(x) for the pressure advection diffusion equation (only used d...
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
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:969
double d_kin_energy_dt() const
Get integral of time derivative of kinetic energy over element.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
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 fill_in_pressure_advection_diffusion_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the residuals and Jacobian for the associated pressure advection diffusion problem...
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number) ...
virtual void get_body_force_gradient_nst(const double &time, const unsigned &ipt, const Vector< double > &s, const Vector< double > &x, DenseMatrix< double > &d_body_force_dx)
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
ProjectableTaylorHoodElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
virtual double p_nst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes...
ProjectableCrouzeixRaviartElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
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 *& re_pt()
Pointer to Reynolds number.
virtual double get_source_nst(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:375
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:562
const Vector< double > & g() const
Vector of gravitational components.
double p_nst(const unsigned &i) const
Return the i-th pressure value (Discontinous pressure interpolation – no need to cater for hanging no...
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:5695
e
Definition: cfortran.h:575
virtual unsigned nplot_points_paraview(const unsigned &nplot) const
Return the number of actual plot points for paraview plot with parameter nplot. Broken virtual; can b...
Definition: elements.h:2700
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...
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
A Rank 4 Tensor class.
Definition: matrices.h:1625
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
int p_local_eqn(const unsigned &n) const
Return the local equation numbers for the pressure values.
QCrouzeixRaviartElement()
Constructor, there are DIM+1 internal values (for the pressure)
virtual void fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)=0
This function returns the residuals for the traction function. flag=1 (or 0): do (or don't) compute t...
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
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...
virtual int p_nodal_index_nst() const
Return the index at which the pressure is stored if it is stored at the nodes. If not stored at the n...
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 u_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...
virtual int p_nodal_index_nst() const
Set the value at which the pressure is stored in the nodes.
virtual void fill_in_pressure_advection_diffusion_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)=0
Compute the residuals and Jacobian for the associated pressure advection diffusion problem...
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
double p_nst(const unsigned &t, const unsigned &i) const
Return the i-th pressure value (Discontinous pressure interpolation – no need to cater for hanging no...
virtual ~FpPressureAdvDiffRobinBCElementBase()
Empty virtual destructor.
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
Data *const & variable_position_pt() const
Pointer to variable_position data (const version)
Definition: nodes.h:1654
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_dresidual_contribution_nst(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, unsigned flag)
Compute the derivatives of the residuals for the Navier–Stokes equations with respect to a parameter ...
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
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...
NavierStokesEquations()
Constructor: NULL the body force and source function and make sure the ALE terms are included by defa...
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)
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
NavierStokesPressureAdvDiffSourceFctPt source_fct_for_pressure_adv_diff() const
Access function for the source-function pointer for pressure advection diffusion (used for validation...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3841
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 size() const
Definition: elements.cc:4121
void output(FILE *file_pt, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
int Pinned_fp_pressure_eqn
Global eqn number of pressure dof that's pinned in pressure advection diffusion problem (defaults to ...
void fill_in_contribution_to_hessian_vector_products(Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Compute the hessian tensor vector products required to perform continuation of bifurcations analytica...
void output(FILE *file_pt)
C-style output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points...
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocity and ...
virtual ~TemplateFreeNavierStokesEquationsBase()
Virtual destructor (empty)
virtual void pin_all_non_pressure_dofs(std::map< Data *, std::vector< int > > &eqn_number_backup)=0
Pin all non-pressure dofs and backup eqn numbers of all Data.
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. Overwrites default implementation in FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}.
A Rank 3 Tensor class.
Definition: matrices.h:1337
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...
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values. (Note: count includes current value!)
unsigned nfields_for_projection()
Number of fields to be projected: dim+1, corresponding to velocity components and pressure...
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...
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
void get_vorticity(const Vector< double > &s, Vector< double > &vorticity) const
Compute the vorticity vector at local coordinate s.
double dpshape_and_dptest_eulerian_nst(const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const
Pressure shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
virtual void fill_in_generic_pressure_advection_diffusion_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
void delete_pressure_advection_diffusion_robin_elements()
Delete the FaceElements that apply the Robin boundary condition to the pressure advection diffusion p...
void point_output_data(const Vector< double > &s, Vector< double > &data)
Output solution in data vector at local cordinates s: x,y [,z], u,v,[w], p.
virtual void delete_pressure_advection_diffusion_robin_elements()=0
Delete the FaceElements that apply the Robin boundary condition to the pressure advection diffusion p...
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
double interpolated_u_nst(const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s.
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.
static char t char * s
Definition: cfortran.h:572
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
int & pinned_fp_pressure_eqn()
Global eqn number of pressure dof that's pinned in pressure adv diff problem.
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
const double & re() const
Reynolds number.
virtual void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)=0
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
double *& re_invfr_pt()
Pointer to global inverse Froude number.
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
virtual int & pinned_fp_pressure_eqn()=0
Global eqn number of pressure dof that's pinned in pressure adv diff problem.
double p_nst(const unsigned &t, const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version) ...
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
Vector< double > * G_pt
Pointer to global gravity Vector.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void scalar_value_paraview(std::ofstream &file_out, const unsigned &i, const unsigned &nplot) const
Write values of the i-th scalar field at the plot points. Needs to be implemented for each new specif...
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...
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
Definition: elements.h:1720
double *& density_ratio_pt()
Pointer to Density ratio.
void fill_in_contribution_to_dresiduals_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam)
Compute the element's residual Vector.
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double pressure_integral() const
Integral of pressure over element.
double dshape_and_dtest_eulerian_at_knot_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...
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
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
double kin_energy() const
Get integral of kinetic energy over element.
Crouzeix Raviart upgraded to become projectable.
void output(std::ostream &outfile, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
double du_dt_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 fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element's residual Vector.
unsigned ninternal_data() const
Return the number of internal data objects.
Definition: elements.h:828
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
Definition: nodes.h:1568
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 interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4277
double dshape_and_dtest_eulerian_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(std::ostream &outfile, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
double dissipation() const
Return integral of dissipation over element.
virtual double dpshape_and_dptest_eulerian_nst(const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const =0
Compute the pressure shape and test functions and derivatives w.r.t. global coords at local coordinat...
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
NavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2470
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:246
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
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 * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
virtual void fill_in_pressure_advection_diffusion_residuals(Vector< double > &residuals)=0
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
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 npres_nst() const
Return number of pressure values.
virtual double dshape_and_dtest_eulerian_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...
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
NavierStokesPressureAdvDiffSourceFctPt Press_adv_diff_source_fct_pt
Pointer to source function pressure advection diffusion equation (only to be used during validation) ...
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2344
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
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
void output(std::ostream &outfile)
Overload the output function.
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
void fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam)
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
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 * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number) ...
virtual void get_source_gradient_nst(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &gradient)
virtual void fill_in_generic_residual_contribution_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 Jacob...
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. Whatever the timestepper has set up for the v...
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Add the element's contribution to its residuals vector, jacobian matrix and mass matrix.
double 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
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
std::string scalar_name_paraview(const unsigned &i) const
Name of the i-th scalar field. Default implementation returns V1 for the first one, V2 for the second etc. Can (should!) be overloaded with more meaningful names in specific elements.
void fill_in_pressure_advection_diffusion_residuals(Vector< double > &residuals)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
unsigned npres_nst() const
Return number of pressure values.
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
virtual int p_nodal_index_nst() const =0
Return the index at which the pressure is stored if it is stored at the nodes. If not stored at the n...
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1164
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...
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.
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
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.
int p_local_eqn(const unsigned &n) const
Return the local equation numbers for the pressure values.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the jacobian.
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types)...
virtual double dshape_and_dtest_eulerian_at_knot_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...
NavierStokesSourceFctPt source_fct_pt() const
Access function for the source-function pointer. Const version.
double u_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_ns...
NavierStokesBodyForceFctPt body_force_fct_pt() const
Access function for the body-force pointer. Const version.
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:596
void interpolated_u_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
double dpshape_and_dptest_eulerian_nst(const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const
Pressure shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
virtual void fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
This function returns the residuals for the traction function. flag=1 (or 0): do (or don't) compute t...
virtual unsigned u_index_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.