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