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-2022 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
44namespace 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.
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
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
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!
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.
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.
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}
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>,
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>,
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
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
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
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
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
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
/////////////////////////////////////////////////////////////////////////
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
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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
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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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.
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
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...
NavierStokesSourceFctPt & source_fct_pt()
Access function for the source-function pointer.
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)
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.
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
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.
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.
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 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.
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
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...
const double & re_invfr() const
Global inverse Froude number.
NavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
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...
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.
GeneralisedNewtonianConstitutiveEquation< DIM > *& constitutive_eqn_pt()
Access function for the constitutive equation pointer.
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
Pointer to body force function.
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.
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
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...
double *& re_invfr_pt()
Pointer to global inverse Froude number.
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
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)
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....
const Vector< double > & g() const
Vector of gravitational components.
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 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 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 identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
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.
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
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 ...
/////////////////////////////////////////////////////////////////////////
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.
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
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 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.
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
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.
double p_nst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version)
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(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
void output(std::ostream &outfile, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
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 int p_local_eqn(const unsigned &n) const =0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
virtual void 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.
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.
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...
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
double 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.
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...
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 ...
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).
//////////////////////////////////////////////////////////////////// ////////////////////////////////...