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