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