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