discontinuous_galerkin_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_DISCONTINUOUS_GALERKIN_SPACETIME_NAVIER_STOKES_ELEMENTS_HEADER
28#define OOMPH_DISCONTINUOUS_GALERKIN_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 just 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* St_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
555 /// Compute the shape functions and derivatives
556 /// w.r.t. global coords at local coordinate s.
557 /// Return Jacobian of mapping between local and global coordinates.
559 Shape& psi,
560 DShape& dpsidx,
561 Shape& test,
562 DShape& dtestdx) const = 0;
563
564
565 /// Compute the shape functions and derivatives
566 /// w.r.t. global coords at ipt-th integration point
567 /// Return Jacobian of mapping between local and global coordinates.
569 const unsigned& ipt,
570 Shape& psi,
571 DShape& dpsidx,
572 Shape& test,
573 DShape& dtestdx) const = 0;
574
575
576 /// Shape/test functions and derivs w.r.t. to global coords at
577 /// integration point ipt; return Jacobian of mapping (J). Also compute
578 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
580 const unsigned& ipt,
581 Shape& psi,
582 DShape& dpsidx,
583 RankFourTensor<double>& d_dpsidx_dX,
584 Shape& test,
585 DShape& dtestdx,
586 RankFourTensor<double>& d_dtestdx_dX,
587 DenseMatrix<double>& djacobian_dX) const = 0;
588
589
590 /// Pressure shape functions and their derivs w.r.t. to global coords
591 /// at local coordinate s (taken from geometry). Return Jacobian of mapping
592 /// between local and global coordinates.
593 virtual double dpshape_eulerian(const Vector<double>& s,
594 Shape& ppsi,
595 DShape& dppsidx) const = 0;
596
597
598 /// Pressure test functions and their derivs w.r.t. to global coords
599 /// at local coordinate s (taken from geometry). Return Jacobian of mapping
600 /// between local and global coordinates.
601 virtual double dptest_eulerian(const Vector<double>& s,
602 Shape& ptest,
603 DShape& dptestdx) const = 0;
604
605
606 /// Compute the pressure shape and test functions and derivatives
607 /// w.r.t. global coords at local coordinate s.
608 /// Return Jacobian of mapping between local and global coordinates.
610 Shape& ppsi,
611 DShape& dppsidx,
612 Shape& ptest,
613 DShape& dptestdx) const = 0;
614
615
616 /// Calculate the body force at a given time and local and/or
617 /// Eulerian position. This function is virtual so that it can be
618 /// overloaded in multi-physics elements where the body force might
619 /// depend on another variable.
620 virtual void get_body_force_nst(const double& time,
621 const unsigned& ipt,
622 const Vector<double>& s,
623 const Vector<double>& x,
624 Vector<double>& result)
625 {
626 // If a function has not been provided
627 if (Body_force_fct_pt == 0)
628 {
629 // Set the body forces to zero
630 result.initialise(0.0);
631 }
632 // If the function pointer is non-zero
633 else
634 {
635 // Call the function
636 (*Body_force_fct_pt)(time, x, result);
637 } // if (Body_force_fct_pt!=0)
638 } // End of get_body_force_nst
639
640
641 /// Get gradient of body force term at (Eulerian) position x. This function
642 /// is virtual to allow overloading in multi-physics problems where
643 /// the strength of the source function might be determined by another
644 /// system of equations. Computed via function pointer (if set) or by
645 /// finite differencing (default)
646 inline virtual void get_body_force_gradient_nst(
647 const double& time,
648 const unsigned& ipt,
649 const Vector<double>& s,
650 const Vector<double>& x,
651 DenseMatrix<double>& d_body_force_dx)
652 {
653 // Reference value
654 Vector<double> body_force(DIM, 0.0);
655
656 // Get the body force vector
657 get_body_force_nst(time, ipt, s, x, body_force);
658
659 // Get the finite-difference step size
661
662 // The body force computed at the perturbed coordinates
663 Vector<double> body_force_pls(DIM, 0.0);
664
665 // Copy the (Eulerian) coordinate vector x
666 Vector<double> x_pls(x);
667
668 // Loop over the coordinate directions
669 for (unsigned i = 0; i < DIM; i++)
670 {
671 // Update the i-th entry
672 x_pls[i] += eps_fd;
673
674 // Get the body force at the current time
675 get_body_force_nst(time, ipt, s, x_pls, body_force_pls);
676
677 // Loop over the coordinate directions
678 for (unsigned j = 0; j < DIM; j++)
679 {
680 // Finite-difference the body force derivative
681 d_body_force_dx(j, i) = (body_force_pls[j] - body_force[j]) / eps_fd;
682 }
683
684 // Reset the i-th entry
685 x_pls[i] = x[i];
686 }
687 } // End of get_body_force_gradient_nst
688
689
690 /// Calculate the source fct at a given time and Eulerian position
691 virtual double get_source_nst(const double& time,
692 const unsigned& ipt,
693 const Vector<double>& x)
694 {
695 // If the function pointer is null
696 if (Source_fct_pt == 0)
697 {
698 // Set the source function value to zero
699 return 0;
700 }
701 // Otherwise call the function pointer
702 else
703 {
704 // Return the appropriate value
705 return (*Source_fct_pt)(time, x);
706 }
707 } // End of get_source_nst
708
709
710 /// Get gradient of source term at (Eulerian) position x. This function
711 /// is virtual to allow overloading in multi-physics problems where the
712 /// strength of the source function might be determined by another system
713 /// of equations. Computed via function pointer (if set) or by finite
714 /// differencing (default)
715 inline virtual void get_source_gradient_nst(const double& time,
716 const unsigned& ipt,
717 const Vector<double>& x,
718 Vector<double>& gradient)
719 {
720 // Reference value
721 double source = get_source_nst(time, ipt, x);
722
723 // Get the finite-difference step size
725
726 // The source function value computed at the perturbed coordinates
727 double source_pls = 0.0;
728
729 // Copy the (Eulerian) coordinate vector, x
730 Vector<double> x_pls(x);
731
732 // Loop over the coordinate directions
733 for (unsigned i = 0; i < DIM; i++)
734 {
735 // Update the i-th entry
736 x_pls[i] += eps_fd;
737
738 // Get the body force at the current time
739 source_pls = get_source_nst(time, ipt, x_pls);
740
741 // Loop over the coordinate directions
742 for (unsigned j = 0; j < DIM; j++)
743 {
744 // Finite-difference the source function gradient
745 gradient[i] = (source_pls - source) / eps_fd;
746 }
747
748 // Reset the i-th entry
749 x_pls[i] = x[i];
750 }
751 } // End of get_source_gradient_nst
752
753
754 /// Compute the residuals for the Navier-Stokes equations.
755 /// Flag=1 (or 0): do (or don't) compute the Jacobian as well.
756 /// Flag=2: Fill in mass matrix too.
758 Vector<double>& residuals,
759 DenseMatrix<double>& jacobian,
760 DenseMatrix<double>& mass_matrix,
761 const unsigned& flag);
762
763
764 /// Compute the residuals for the associated pressure advection
765 /// diffusion problem. Used by the Fp preconditioner.
766 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
768 Vector<double>& residuals,
769 DenseMatrix<double>& jacobian,
770 const unsigned& flag);
771
772
773 /// Compute the derivatives of the
774 /// residuals for the Navier-Stokes equations with respect to a parameter
775 /// Flag=1 (or 0): do (or don't) compute the Jacobian as well.
776 /// Flag=2: Fill in mass matrix too.
778 double* const& parameter_pt,
779 Vector<double>& dres_dparam,
780 DenseMatrix<double>& djac_dparam,
781 DenseMatrix<double>& dmass_matrix_dparam,
782 const unsigned& flag);
783
784
785 /// Compute the hessian tensor vector products required to
786 /// perform continuation of bifurcations analytically
788 Vector<double> const& Y,
789 DenseMatrix<double> const& C,
790 DenseMatrix<double>& product);
791
792 public:
793 /// Constructor: NULL the body force and source function
794 /// and make sure the ALE terms are included by default.
797 Source_fct_pt(0),
799 ALE_is_disabled(false),
801 {
802 // Set all the Physical parameter pointers:
803
804 // Set the Reynolds number to the value zero
806
807 // Set the Strouhal number to the value zero
809
810 // The Strouhal number (which is a proxy for the period here) is not
811 // stored as external data
813
814 // Set the Reynolds / Froude value to zero
816
817 // Set the gravity vector to be the zero vector
819
820 // Set the Physical ratios:
821
822 // Set the viscosity ratio to the value one
824
825 // Set the density ratio to the value one
827 }
828
829 /// Vector to decide whether the stress-divergence form is used or not
830 /// N.B. This needs to be public so that the intel compiler gets things
831 /// correct somehow the access function messes things up when going to
832 /// refineable Navier-Stokes
833 static Vector<double> Gamma;
834
835
836 /// Function that tells us whether the period is stored as external data
837 void store_strouhal_as_external_data(Data* strouhal_data_pt)
838 {
839#ifdef PARANOID
840 // Sanity check; make sure it's not a null pointer
841 if (strouhal_data_pt == 0)
842 {
843 // Create an output stream
844 std::ostringstream error_message_stream;
845
846 // Create an error message
847 error_message_stream
848 << "User supplied Strouhal number as external data\n"
849 << "but the pointer provided is a null pointer!" << std::endl;
850
851 // Throw an error
852 throw OomphLibError(error_message_stream.str(),
853 OOMPH_CURRENT_FUNCTION,
854 OOMPH_EXCEPTION_LOCATION);
855 }
856#endif
857
858 // Set the Strouhal value pointer (the Strouhal number is stored as the
859 // only piece of internal data in the phase condition element)
860 this->add_external_data(strouhal_data_pt);
861
862 // Indicate that the Strouhal number is store as external data
864 } // End of store_strouhal_as_external_data
865
866
867 // Access functions for the physical constants:
868
869 /// Reynolds number
870 const double& re() const
871 {
872 // Use the Reynolds number pointer to get the Reynolds value
873 return *Re_pt;
874 } // End of re
875
876
877 /// Pointer to Reynolds number
878 double*& re_pt()
879 {
880 // Return the Reynolds number pointer
881 return Re_pt;
882 } // End of re_pt
883
884
885 /// Are we storing the Strouhal number as external data?
887 {
888 // Return the flags value
890 } // End of is_strouhal_stored_as_external_data
891
892
893 /// Strouhal parameter (const. version)
894 const double& st() const
895 {
896 // If the st is stored as external data
898 {
899 // The index of the external data (which contains the st)
900 unsigned data_index = 0;
901
902 // The index of the value at which the Strouhal is stored
903 unsigned strouhal_index = 0;
904
905 // Return the value of the st in the external data
906 return *(this->external_data_pt(data_index)->value_pt(strouhal_index));
907 }
908 // Otherwise the st is just stored as a pointer
909 else
910 {
911 // Return the value of St
912 return *St_pt;
913 }
914 } // End of st
915
916
917 /// Pointer to Strouhal parameter (const. version)
918 double* st_pt() const
919 {
920 // If the strouhal is stored as external data
922 {
923 // The index of the external data (which contains the strouhal)
924 unsigned data_index = 0;
925
926 // The index of the value at which the strouhal is stored (the only
927 // value)
928 unsigned strouhal_index = 0;
929
930 // Return the value of the strouhal in the external data
931 return this->external_data_pt(data_index)->value_pt(strouhal_index);
932 }
933 // Otherwise the strouhal is just stored as a pointer
934 else
935 {
936 // Return the value of Strouhal
937 return St_pt;
938 }
939 } // End of st_pt
940
941
942 /// Pointer to Strouhal number (can only assign to private member data)
943 double*& st_pt()
944 {
945 // Return the Strouhal number pointer
946 return St_pt;
947 } // End of st_pt
948
949
950 /// Product of Reynolds and Strouhal number (=Womersley number)
951 double re_st() const
952 {
953 // Return the product of the appropriate physical constants
954 return (*Re_pt) * (*st_pt());
955 } // End of re_st
956
957
958 /// Viscosity ratio for element: Element's viscosity relative to the
959 /// viscosity used in the definition of the Reynolds number
960 const double& viscosity_ratio() const
961 {
962 return *Viscosity_Ratio_pt;
963 }
964
965 /// Pointer to Viscosity Ratio
967 {
968 return Viscosity_Ratio_pt;
969 }
970
971 /// Density ratio for element: Element's density relative to the
972 /// viscosity used in the definition of the Reynolds number
973 const double& density_ratio() const
974 {
975 return *Density_Ratio_pt;
976 }
977
978 /// Pointer to Density ratio
980 {
981 return Density_Ratio_pt;
982 }
983
984 /// Global inverse Froude number
985 const double& re_invfr() const
986 {
987 return *ReInvFr_pt;
988 }
989
990 /// Pointer to global inverse Froude number
991 double*& re_invfr_pt()
992 {
993 return ReInvFr_pt;
994 }
995
996 /// Vector of gravitational components
997 const Vector<double>& g() const
998 {
999 return *G_pt;
1000 }
1001
1002 /// Pointer to Vector of gravitational components
1004 {
1005 return G_pt;
1006 }
1007
1008 /// Access function for the body-force pointer
1010 {
1011 return Body_force_fct_pt;
1012 }
1013
1014 /// Access function for the body-force pointer. Const version
1016 {
1017 return Body_force_fct_pt;
1018 }
1019
1020 /// Access function for the source-function pointer
1022 {
1023 return Source_fct_pt;
1024 }
1025
1026 /// Access function for the source-function pointer. Const version
1028 {
1029 return Source_fct_pt;
1030 }
1031
1032 /// Access function for the source-function pointer for pressure
1033 /// advection diffusion (used for validation only).
1035 {
1037 }
1038
1039 /// Access function for the source-function pointer for pressure
1040 /// advection diffusion (used for validation only). Const version.
1042 const
1043 {
1045 }
1046
1047 /// Global eqn number of pressure dof that's pinned in pressure
1048 /// adv diff problem
1050 {
1051 // Return the appropriate equation number
1053 }
1054
1055 /// Function to return number of pressure degrees of freedom
1056 virtual unsigned npres_nst() const = 0;
1057
1058 /// Compute the pressure shape functions at local coordinate s
1059 virtual void pshape_nst(const Vector<double>& s, Shape& psi) const = 0;
1060
1061 /// Compute the pressure shape and test functions
1062 /// at local coordinate s
1063 virtual void pshape_nst(const Vector<double>& s,
1064 Shape& psi,
1065 Shape& test) const = 0;
1066
1067 /// Velocity i at local node n. Uses suitably interpolated value
1068 /// for hanging nodes. The use of u_index_nst() permits the use of this
1069 /// element as the basis for multi-physics elements. The default
1070 /// is to assume that the i-th velocity component is stored at the
1071 /// i-th location of the node
1072 double u_nst(const unsigned& n, const unsigned& i) const
1073 {
1074 return nodal_value(n, u_index_nst(i));
1075 }
1076
1077 /// Velocity i at local node n at timestep t (t=0: present;
1078 /// t>0: previous). Uses suitably interpolated value for hanging nodes.
1079 double u_nst(const unsigned& t, const unsigned& n, const unsigned& i) const
1080 {
1081#ifdef PARANOID
1082 // Since we're using space-time elements, this only makes sense if t=0
1083 if (t != 0)
1084 {
1085 // Throw an error
1086 throw OomphLibError("Space-time elements cannot have history values!",
1087 OOMPH_CURRENT_FUNCTION,
1088 OOMPH_EXCEPTION_LOCATION);
1089 }
1090#endif
1091
1092 // Return the appropriate nodal value (noting that t=0 here)
1093 return nodal_value(t, n, u_index_nst(i));
1094 } // End of u_nst
1095
1096 /// Return the index at which the i-th unknown velocity component
1097 /// is stored. The default value, i, is appropriate for single-physics
1098 /// problems.
1099 /// In derived multi-physics elements, this function should be overloaded
1100 /// to reflect the chosen storage scheme. Note that these equations require
1101 /// that the unknowns are always stored at the same indices at each node.
1102 virtual inline unsigned u_index_nst(const unsigned& i) const
1103 {
1104#ifdef PARANOID
1105 if (i > DIM - 1)
1106 {
1107 // Create an output stream
1108 std::ostringstream error_message_stream;
1109
1110 // Create an error message
1111 error_message_stream << "Input index " << i << " does not correspond "
1112 << "to a velocity component when DIM=" << DIM
1113 << "!" << std::endl;
1114
1115 // Throw an error
1116 throw OomphLibError(error_message_stream.str(),
1117 OOMPH_CURRENT_FUNCTION,
1118 OOMPH_EXCEPTION_LOCATION);
1119 }
1120#endif
1121
1122 // Return the appropriate entry
1123 return i;
1124 } // End of u_index_nst
1125
1126
1127 /// Return the number of velocity components (used in
1128 /// FluidInterfaceElements)
1129 inline unsigned n_u_nst() const
1130 {
1131 // Return the number of equations to solve
1132 return DIM;
1133 } // End of n_u_nst
1134
1135
1136 /// i-th component of du/dt at local node n.
1137 /// Uses suitably interpolated value for hanging nodes.
1138 /// NOTE: This is essentially a wrapper for du_dt_nst()
1139 /// so it can be called externally.
1140 double get_du_dt(const unsigned& n, const unsigned& i) const
1141 {
1142 // Return the value calculated by du_dt_vdp
1143 return du_dt_nst(n, i);
1144 } // End of get_du_dt
1145
1146
1147 /// i-th component of du/dt at local node n.
1148 /// Uses suitably interpolated value for hanging nodes.
1149 double du_dt_nst(const unsigned& n, const unsigned& i) const
1150 {
1151 // Storage for the local coordinates
1152 Vector<double> s(DIM + 1, 0.0);
1153
1154 // Get the local coordinate at the n-th node
1156
1157 // Return the interpolated du_i/dt value
1158 return interpolated_du_dt_nst(s, i);
1159 } // End of du_dt_nst
1160
1161
1162 /// Return FE representation of function value du_i/dt(s) at local
1163 /// coordinate s
1165 const unsigned& i) const
1166 {
1167 // Find number of nodes
1168 unsigned n_node = nnode();
1169
1170 // Local shape function
1171 Shape psi(n_node);
1172
1173 // Allocate space for the derivatives of the shape functions
1174 DShape dpsidx(n_node, DIM + 1);
1175
1176 // Compute the geometric shape functions and also first derivatives
1177 // w.r.t. global coordinates at local coordinate s
1178 dshape_eulerian(s, psi, dpsidx);
1179
1180 // Initialise value of du_i/dt
1181 double interpolated_dudt = 0.0;
1182
1183 // Find the index at which the variable is stored
1184 unsigned u_nodal_index = u_index_nst(i);
1185
1186 // Loop over the local nodes and sum
1187 for (unsigned l = 0; l < n_node; l++)
1188 {
1189 // Update the interpolated du_i/dt value
1190 interpolated_dudt += nodal_value(l, u_nodal_index) * dpsidx(l, DIM);
1191 }
1192
1193 // Return the interpolated du_i/dt value
1194 return interpolated_dudt;
1195 } // End of interpolated_du_dt_nst
1196
1197
1198 /// Disable ALE, i.e. assert the mesh is not moving -- you do this
1199 /// at your own risk!
1201 {
1202 ALE_is_disabled = true;
1203 } // End of disable_ALE
1204
1205
1206 /// (Re-)enable ALE, i.e. take possible mesh motion into account
1207 /// when evaluating the time-derivative. Note: By default, ALE is
1208 /// enabled, at the expense of possibly creating unnecessary work
1209 /// in problems where the mesh is, in fact, stationary.
1211 {
1212 ALE_is_disabled = false;
1213 } // End of enable_ALE
1214
1215
1216 /// Pressure at local pressure "node" n_p
1217 /// Uses suitably interpolated value for hanging nodes.
1218 virtual double p_nst(const unsigned& n_p) const = 0;
1219
1220 /// Pressure at local pressure "node" n_p at time level t
1221 virtual double p_nst(const unsigned& t, const unsigned& n_p) const = 0;
1222
1223 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
1224 virtual void fix_pressure(const unsigned& p_dof, const double& p_value) = 0;
1225
1226 /// Return the index at which the pressure is stored if it is
1227 /// stored at the nodes. If not stored at the nodes this will return
1228 /// a negative number.
1229 virtual int p_nodal_index_nst() const
1230 {
1232 }
1233
1234 /// Integral of pressure over element
1235 double pressure_integral() const;
1236
1237 /// Return integral of dissipation over element
1238 double dissipation() const;
1239
1240 /// Return dissipation at local coordinate s
1241 double dissipation(const Vector<double>& s) const;
1242
1243 /// Compute the vorticity vector at local coordinate s
1245 Vector<double>& vorticity) const;
1246
1247 /// Compute the vorticity vector at local coordinate s
1248 void get_vorticity(const Vector<double>& s, double& vorticity) const;
1249
1250 /// Get integral of kinetic energy over element
1251 double kin_energy() const;
1252
1253 /// Get integral of time derivative of kinetic energy over element
1254 double d_kin_energy_dt() const;
1255
1256 /// Strain-rate tensor: 1/2 (du_i/dx_j+du_j/dx_i)
1259
1260 /// Compute traction (on the viscous scale) exerted onto
1261 /// the fluid at local coordinate s. N has to be outer unit normal
1262 /// to the fluid.
1264 const Vector<double>& N,
1265 Vector<double>& traction);
1266
1267 /// Compute traction (on the viscous scale) exerted onto
1268 /// the fluid at local coordinate s, decomposed into pressure and
1269 /// normal and tangential viscous components.
1270 /// N has to be outer unit normal to the fluid.
1272 const Vector<double>& N,
1273 Vector<double>& traction_p,
1274 Vector<double>& traction_visc_n,
1275 Vector<double>& traction_visc_t);
1276
1277 /// This implements a pure virtual function defined
1278 /// in the FSIFluidElement class. The function computes
1279 /// the traction (on the viscous scale), at the element's local
1280 /// coordinate s, that the fluid element exerts onto an adjacent
1281 /// solid element. The number of arguments is imposed by
1282 /// the interface defined in the FSIFluidElement -- only the
1283 /// unit normal N (pointing into the fluid!) is actually used
1284 /// in the computation.
1286 const Vector<double>& N,
1287 Vector<double>& load)
1288 {
1289 // Note: get_traction() computes the traction onto the fluid
1290 // if N is the outer unit normal onto the fluid; here we're
1291 // expecting N to point into the fluid so we're getting the
1292 // traction onto the adjacent wall instead!
1293 get_traction(s, N, load);
1294 }
1295
1296 /// Compute the diagonal of the velocity/pressure mass matrices.
1297 /// If which one=0, both are computed, otherwise only the pressure
1298 /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
1299 /// LSC version of the preconditioner only needs that one)
1301 Vector<double>& press_mass_diag,
1302 Vector<double>& veloc_mass_diag,
1303 const unsigned& which_one = 0);
1304
1305 /// Number of scalars/fields output by this element. Reimplements
1306 /// broken virtual function in base class.
1307 unsigned nscalar_paraview() const
1308 {
1309 // The number of velocity components plus the pressure field
1310 return DIM + 1;
1311 }
1312
1313 /// Write values of the i-th scalar field at the plot points. Needs
1314 /// to be implemented for each new specific element type.
1315 void scalar_value_paraview(std::ofstream& file_out,
1316 const unsigned& i,
1317 const unsigned& nplot) const
1318 {
1319 // Vector of local coordinates
1320 Vector<double> s(DIM + 1, 0.0);
1321
1322 // How many plot points do we have in total?
1323 unsigned num_plot_points = nplot_points_paraview(nplot);
1324
1325 // Loop over plot points
1326 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1327 {
1328 // Get the local coordinates of the iplot-th plot point
1329 get_s_plot(iplot, nplot, s);
1330
1331 // Velocities
1332 if (i < DIM)
1333 {
1334 // Output the i-th velocity component
1335 file_out << interpolated_u_nst(s, i) << std::endl;
1336 }
1337 // Pressure
1338 else if (i == DIM)
1339 {
1340 // Output the pressure at this point
1341 file_out << interpolated_p_nst(s) << std::endl;
1342 }
1343 // Never get here
1344 else
1345 {
1346#ifdef PARANOID
1347 // Create an output stream
1348 std::stringstream error_stream;
1349
1350 // Create the error message
1351 error_stream << "These Navier Stokes elements only store " << DIM + 1
1352 << " fields, "
1353 << "but i is currently " << i << std::endl;
1354
1355 // Throw the error message
1356 throw OomphLibError(error_stream.str(),
1357 OOMPH_CURRENT_FUNCTION,
1358 OOMPH_EXCEPTION_LOCATION);
1359#endif
1360 }
1361 } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
1362 } // End of scalar_value_paraview
1363
1364
1365 /// Write values of the i-th scalar field at the plot points. Needs
1366 /// to be implemented for each new specific element type.
1368 std::ofstream& file_out,
1369 const unsigned& i,
1370 const unsigned& nplot,
1371 const double& time,
1373 {
1374#ifdef PARANOID
1375 if (i > DIM)
1376 {
1377 // Create an output stream
1378 std::stringstream error_stream;
1379
1380 // Create the error message
1381 error_stream << "These Navier Stokes elements only store " << DIM + 1
1382 << " fields, but i is currently " << i << std::endl;
1383
1384 // Throw the error message
1385 throw OomphLibError(
1386 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1387 }
1388#endif
1389
1390 // Vector of local coordinates
1391 Vector<double> s(DIM + 1, 0.0);
1392
1393 // Storage for the time value
1394 double interpolated_t = 0.0;
1395
1396 // Storage for the spatial coordinates
1397 Vector<double> spatial_coordinates(DIM, 0.0);
1398
1399 // How many plot points do we have in total?
1400 unsigned num_plot_points = nplot_points_paraview(nplot);
1401
1402 // Loop over plot points
1403 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1404 {
1405 // Get the local coordinates of the iplot-th plot point
1406 get_s_plot(iplot, nplot, s);
1407
1408 // Loop over the spatial coordinates
1409 for (unsigned i = 0; i < DIM; i++)
1410 {
1411 // Assign the i-th spatial coordinate
1412 spatial_coordinates[i] = interpolated_x(s, i);
1413 }
1414
1415 // Get the time value
1416 interpolated_t = interpolated_x(s, DIM);
1417
1418 // ------ DRAIG: REMOVE ----------------------------------------
1419 // Exact solution vector (here it's simply a scalar)
1420 Vector<double> exact_soln(DIM + 1, 0.0);
1421 // DRAIG: Needs to be changed to a Vector of size DIM
1422 // ------ DRAIG: REMOVE ----------------------------------------
1423
1424 // Get the exact solution at this point
1425 (*exact_soln_pt)(interpolated_t, spatial_coordinates, exact_soln);
1426
1427 // Output the interpolated solution value
1428 file_out << exact_soln[i] << std::endl;
1429 } // for (unsigned iplot=0;iplot<num_plot_points;iplot++)
1430 } // End of scalar_value_fct_paraview
1431
1432
1433 /// Name of the i-th scalar field. Default implementation
1434 /// returns V1 for the first one, V2 for the second etc. Can (should!) be
1435 /// overloaded with more meaningful names in specific elements.
1436 std::string scalar_name_paraview(const unsigned& i) const
1437 {
1438 // Velocities
1439 if (i < DIM)
1440 {
1441 // Return the appropriate string for the i-th velocity component
1442 return "Velocity " + StringConversion::to_string(i);
1443 }
1444 // Pressure
1445 else if (i == DIM)
1446 {
1447 // Return the name for the pressure
1448 return "Pressure";
1449 }
1450 // Never get here
1451 else
1452 {
1453 // Create an output stream
1454 std::stringstream error_stream;
1455
1456 // Create the error message
1457 error_stream << "These Navier Stokes elements only store " << DIM + 1
1458 << " fields,\nbut i is currently " << i << std::endl;
1459
1460 // Throw the error message
1461 throw OomphLibError(
1462 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1463
1464 // Dummy return so the compiler doesn't complain
1465 return " ";
1466 }
1467 } // End of scalar_name_paraview
1468
1469
1470 /// Output function: x,y,t,u,v,p in tecplot format. The
1471 /// default number of plot points is five
1472 void output(std::ostream& outfile)
1473 {
1474 // Set the number of plot points in each direction
1475 unsigned n_plot = 5;
1476
1477 // Forward the call to the appropriate output function
1478 output(outfile, n_plot);
1479 } // End of output
1480
1481
1482 /// Output function: x,y,[z],u,v,[w],p in tecplot format. Here,
1483 /// we use n_plot plot points in each coordinate direction
1484 void output(std::ostream& outfile, const unsigned& n_plot);
1485
1486
1487 /// C-style output function: x,y,[z],u,v,[w],p in tecplot format. The
1488 /// default number of plot points is five
1489 void output(FILE* file_pt)
1490 {
1491 // Set the number of plot points in each direction
1492 unsigned n_plot = 5;
1493
1494 // Forward the call to the appropriate output function
1495 output(file_pt, n_plot);
1496 } // End of output
1497
1498
1499 /// C-style output function: x,y,[z],u,v,[w],p in tecplot format. Use
1500 /// n_plot points in each coordinate direction
1501 void output(FILE* file_pt, const unsigned& n_plot);
1502
1503
1504 /// Full output function:
1505 /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format. The
1506 /// default number of plot points is five
1507 void full_output(std::ostream& outfile)
1508 {
1509 // Set the number of plot points in each direction
1510 unsigned n_plot = 5;
1511
1512 // Forward the call to the appropriate output function
1513 full_output(outfile, n_plot);
1514 } // End of full_output
1515
1516
1517 /// Full output function:
1518 /// x,y,t,u,v,p,du/dt,dv/dt,dissipation in tecplot format. Use
1519 /// n_plot plot points in each coordinate direction
1520 void full_output(std::ostream& outfile, const unsigned& n_plot);
1521
1522
1523 /// Output function: x,y,t,u,v in tecplot format. Use
1524 /// n_plot points in each coordinate direction at timestep t (t=0:
1525 /// present; t>0: previous timestep)
1526 void output_veloc(std::ostream& outfile,
1527 const unsigned& nplot,
1528 const unsigned& t);
1529
1530
1531 /// Output function: x,y,t,omega
1532 /// in tecplot format. nplot points in each coordinate direction
1533 void output_vorticity(std::ostream& outfile, const unsigned& nplot);
1534
1535
1536 /// Output exact solution specified via function pointer
1537 /// at a given number of plot points. Function prints as
1538 /// many components as are returned in solution Vector
1539 void output_fct(std::ostream& outfile,
1540 const unsigned& nplot,
1542
1543
1544 /// Output exact solution specified via function pointer
1545 /// at a given time and at a given number of plot points.
1546 /// Function prints as many components as are returned in solution Vector.
1547 void output_fct(std::ostream& outfile,
1548 const unsigned& nplot,
1549 const double& time,
1551
1552
1553 /// Compute the vector norm of the FEM solution
1555
1556
1557 /// Validate against exact solution at given time
1558 /// Solution is provided via function pointer.
1559 /// Plot at a given number of plot points and compute L2 error
1560 /// and L2 norm of velocity solution over element
1561 void compute_error(std::ostream& outfile,
1563 const double& time,
1564 double& error,
1565 double& norm);
1566
1567
1568 /// Validate against exact solution at given time
1569 /// Solution is provided via function pointer.
1570 /// Plot at a given number of plot points and compute L2 error
1571 /// and L2 norm of velocity solution over element
1572 void compute_error(std::ostream& outfile,
1574 const double& time,
1575 Vector<double>& error,
1576 Vector<double>& norm);
1577
1578
1579 /// Validate against exact solution.
1580 /// Solution is provided via function pointer.
1581 /// Plot at a given number of plot points and compute L2 error
1582 /// and L2 norm of velocity solution over element
1583 void compute_error(std::ostream& outfile,
1585 double& error,
1586 double& norm);
1587
1588
1589 /// Validate against exact solution. Solution is provided via
1590 /// function pointer. Compute L2 error and L2 norm of velocity solution
1591 /// over element.
1593 const double& time,
1594 double& error,
1595 double& norm);
1596
1597
1598 /// Validate against exact solution. Solution is provided via
1599 /// function pointer. Compute L2 error and L2 norm of velocity solution
1600 /// over element.
1602 double& error,
1603 double& norm);
1604
1605
1606 /// Compute the element's residual Vector
1608 {
1609 // Do we want to compute the Jacobian? ANSWER: No!
1610 unsigned compute_jacobian_flag = 0;
1611
1612 // Call the generic residuals function using a dummy matrix argument
1614 residuals,
1617 compute_jacobian_flag);
1618 } // End of fill_in_contribution_to_residuals
1619
1620
1621 /// Compute the element's residual Vector and the jacobian matrix
1622 /// Virtual function can be overloaded by hanging-node version
1624 DenseMatrix<double>& jacobian)
1625 {
1626 // Do we want to compute the Jacobian? ANSWER: Yes!
1627 unsigned compute_jacobian_flag = 1;
1628
1629 // Call the generic routine with the flag set to 1
1631 residuals,
1632 jacobian,
1634 compute_jacobian_flag);
1635 } // End of fill_in_contribution_to_residuals
1636
1637
1638 /// Add the element's contribution to its residuals vector,
1639 /// jacobian matrix and mass matrix
1641 Vector<double>& residuals,
1642 DenseMatrix<double>& jacobian,
1643 DenseMatrix<double>& mass_matrix)
1644 {
1645 // Do we want to compute the Jacobian AND mass matrix? ANSWER: Yes!
1646 unsigned compute_matrices_flag = 2;
1647
1648 // Call the generic routine with the flag set to 2
1650 residuals, jacobian, mass_matrix, compute_matrices_flag);
1651 } // End of fill_in_contribution_to_jacobian_and_mass_matrix
1652
1653
1654 /// Compute the element's residual Vector (differentiated w.r.t. a
1655 /// parameter)
1657 double* const& parameter_pt, Vector<double>& dres_dparam)
1658 {
1659 // Do we want to compute the Jacobian? ANSWER: No!
1660 unsigned compute_jacobian_flag = 0;
1661
1662 // Call the generic residuals function using a dummy matrix argument
1664 parameter_pt,
1665 dres_dparam,
1668 compute_jacobian_flag);
1669 } // End of fill_in_contribution_to_dresiduals_dparameter
1670
1671
1672 /// Compute the element's residual Vector and the jacobian matrix
1673 /// Virtual function can be overloaded by hanging-node version
1675 double* const& parameter_pt,
1676 Vector<double>& dres_dparam,
1677 DenseMatrix<double>& djac_dparam)
1678 {
1679 // Do we want to compute the Jacobian? ANSWER: Yes!
1680 unsigned compute_jacobian_flag = 1;
1681
1682 // Call the generic routine with the flag set to 1
1684 parameter_pt,
1685 dres_dparam,
1686 djac_dparam,
1688 compute_jacobian_flag);
1689 } // End of fill_in_contribution_to_djacobian_dparameter
1690
1691
1692 /// Add the element's contribution to its residuals vector,
1693 /// jacobian matrix and mass matrix
1695 double* const& parameter_pt,
1696 Vector<double>& dres_dparam,
1697 DenseMatrix<double>& djac_dparam,
1698 DenseMatrix<double>& dmass_matrix_dparam)
1699 {
1700 // Do we want to compute the Jacobian AND mass matrix? ANSWER: Yes!
1701 unsigned compute_matrices_flag = 2;
1702
1703 // Call the generic routine with the appropriate flag
1705 dres_dparam,
1706 djac_dparam,
1707 dmass_matrix_dparam,
1708 compute_matrices_flag);
1709 } // End of fill_in_contribution_to_djacobian_and_dmass_matrix_dparameter
1710
1711
1712 /// Compute the residuals for the associated pressure advection
1713 /// diffusion problem. Used by the Fp preconditioner.
1715 Vector<double>& residuals)
1716 {
1717 // Do we want to compute the Jacobian? ANSWER: No!
1718 unsigned compute_jacobian_flag = 0;
1719
1720 // Call the generic routine with the appropriate flag
1722 residuals, GeneralisedElement::Dummy_matrix, compute_jacobian_flag);
1723 } // End of fill_in_pressure_advection_diffusion_residuals
1724
1725
1726 /// Compute the residuals and Jacobian for the associated
1727 /// pressure advection diffusion problem. Used by the Fp preconditioner.
1729 Vector<double>& residuals, DenseMatrix<double>& jacobian)
1730 {
1731 // Do we want to compute the Jacobian? ANSWER: Yes!
1732 unsigned compute_jacobian_flag = 1;
1733
1734 // Call the generic routine with the appropriate flag
1736 residuals, jacobian, compute_jacobian_flag);
1737 } // End of fill_in_pressure_advection_diffusion_jacobian
1738
1739
1740 /// Pin all non-pressure dofs and backup eqn numbers
1742 std::map<Data*, std::vector<int>>& eqn_number_backup)
1743 {
1744 // Loop over internal data and pin the values (having established that
1745 // pressure dofs aren't amongst those)
1746 unsigned nint = this->ninternal_data();
1747 for (unsigned j = 0; j < nint; j++)
1748 {
1749 Data* data_pt = this->internal_data_pt(j);
1750 if (eqn_number_backup[data_pt].size() == 0)
1751 {
1752 unsigned nvalue = data_pt->nvalue();
1753 eqn_number_backup[data_pt].resize(nvalue);
1754 for (unsigned i = 0; i < nvalue; i++)
1755 {
1756 // Backup
1757 eqn_number_backup[data_pt][i] = data_pt->eqn_number(i);
1758
1759 // Pin everything
1760 data_pt->pin(i);
1761 }
1762 }
1763 }
1764
1765 // Now deal with nodal values
1766 unsigned nnod = this->nnode();
1767 for (unsigned j = 0; j < nnod; j++)
1768 {
1769 Node* nod_pt = this->node_pt(j);
1770 if (eqn_number_backup[nod_pt].size() == 0)
1771 {
1772 unsigned nvalue = nod_pt->nvalue();
1773 eqn_number_backup[nod_pt].resize(nvalue);
1774 for (unsigned i = 0; i < nvalue; i++)
1775 {
1776 // Pin everything apart from the nodal pressure
1777 // value
1778 if (int(i) != this->p_nodal_index_nst())
1779 {
1780 // Backup
1781 eqn_number_backup[nod_pt][i] = nod_pt->eqn_number(i);
1782
1783 // Pin
1784 nod_pt->pin(i);
1785 }
1786 // Else it's a pressure value
1787 else
1788 {
1789 // Exclude non-nodal pressure based elements
1790 if (this->p_nodal_index_nst() >= 0)
1791 {
1792 // Backup
1793 eqn_number_backup[nod_pt][i] = nod_pt->eqn_number(i);
1794 }
1795 }
1796 }
1797
1798
1799 // If it's a solid node deal with its positional data too
1800 SolidNode* solid_nod_pt = dynamic_cast<SolidNode*>(nod_pt);
1801 if (solid_nod_pt != 0)
1802 {
1803 Data* solid_posn_data_pt = solid_nod_pt->variable_position_pt();
1804 if (eqn_number_backup[solid_posn_data_pt].size() == 0)
1805 {
1806 unsigned nvalue = solid_posn_data_pt->nvalue();
1807 eqn_number_backup[solid_posn_data_pt].resize(nvalue);
1808 for (unsigned i = 0; i < nvalue; i++)
1809 {
1810 // Backup
1811 eqn_number_backup[solid_posn_data_pt][i] =
1812 solid_posn_data_pt->eqn_number(i);
1813
1814 // Pin
1815 solid_posn_data_pt->pin(i);
1816 }
1817 }
1818 }
1819 }
1820 }
1821 } // End of pin_all_non_pressure_dofs
1822
1823
1824 /// Build FaceElements that apply the Robin boundary condition
1825 /// to the pressure advection diffusion problem required by
1826 /// Fp preconditioner
1828 const unsigned& face_index) = 0;
1829
1830
1831 /// Output the FaceElements that apply the Robin boundary condition
1832 /// to the pressure advection diffusion problem required by
1833 /// Fp preconditioner
1835 std::ostream& outfile)
1836 {
1838 for (unsigned e = 0; e < nel; e++)
1839 {
1840 FaceElement* face_el_pt =
1842 outfile << "ZONE" << std::endl;
1843 Vector<double> unit_normal(DIM);
1844 Vector<double> x(DIM + 1);
1845 Vector<double> s(DIM);
1846 unsigned n = face_el_pt->integral_pt()->nweight();
1847 for (unsigned ipt = 0; ipt < n; ipt++)
1848 {
1849 for (unsigned i = 0; i < DIM; i++)
1850 {
1851 s[i] = face_el_pt->integral_pt()->knot(ipt, i);
1852 }
1853 face_el_pt->interpolated_x(s, x);
1854 face_el_pt->outer_unit_normal(ipt, unit_normal);
1855 for (unsigned i = 0; i < DIM + 1; i++)
1856 {
1857 outfile << x[i] << " ";
1858 }
1859 for (unsigned i = 0; i < DIM; i++)
1860 {
1861 outfile << unit_normal[i] << " ";
1862 }
1863 outfile << std::endl;
1864 }
1865 }
1866 } // End of output_pressure_advection_diffusion_robin_elements
1867
1868
1869 /// Delete the FaceElements that apply the Robin boundary condition
1870 /// to the pressure advection diffusion problem required by
1871 /// Fp preconditioner
1873 {
1875 for (unsigned e = 0; e < nel; e++)
1876 {
1878 }
1880 } // End of delete_pressure_advection_diffusion_robin_elements
1881
1882
1883 /// Compute derivatives of elemental residual vector with respect
1884 /// to nodal coordinates. Overwrites default implementation in
1885 /// FiniteElement base class.
1886 /// dresidual_dnodal_coordinates(l,i,j)=d res(l) / dX_{ij}
1888 RankThreeTensor<double>& dresidual_dnodal_coordinates);
1889
1890
1891 /// Compute vector of FE interpolated velocity u at local coordinate s
1893 Vector<double>& velocity) const
1894 {
1895 // Find the number of nodes
1896 unsigned n_node = nnode();
1897
1898 // Local shape function
1899 Shape psi(n_node);
1900
1901 // Find values of shape function
1902 shape(s, psi);
1903
1904 // Loop over the velocity components
1905 for (unsigned i = 0; i < DIM; i++)
1906 {
1907 // Index at which the nodal value is stored
1908 unsigned u_nodal_index = u_index_nst(i);
1909
1910 // Initialise the i-th value of the velocity vector
1911 velocity[i] = 0.0;
1912
1913 // Loop over the local nodes and sum
1914 for (unsigned l = 0; l < n_node; l++)
1915 {
1916 // Update the i-th velocity component approximation
1917 velocity[i] += nodal_value(l, u_nodal_index) * psi[l];
1918 }
1919 } // for (unsigned i=0;i<DIM;i++)
1920 } // End of interpolated_u_nst
1921
1922
1923 /// Return FE interpolated velocity u[i] at local coordinate s
1924 double interpolated_u_nst(const Vector<double>& s, const unsigned& i) const
1925 {
1926 // Find the number of nodes
1927 unsigned n_node = nnode();
1928
1929 // Local shape function
1930 Shape psi(n_node);
1931
1932 // Find values of shape function
1933 shape(s, psi);
1934
1935 // Get nodal index at which i-th velocity is stored
1936 unsigned u_nodal_index = u_index_nst(i);
1937
1938 // Initialise value of u
1939 double interpolated_u = 0.0;
1940
1941 // Loop over the local nodes and sum
1942 for (unsigned l = 0; l < n_node; l++)
1943 {
1944 // Update the i-th velocity component approximation
1945 interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
1946 }
1947
1948 // Return the calculated approximation to the i-th velocity component
1949 return interpolated_u;
1950 } // End of interpolated_u_nst
1951
1952
1953 /// Return FE interpolated velocity u[i] at local coordinate s
1954 /// time level, t. Purposely broken for space-time elements.
1955 double interpolated_u_nst(const unsigned& t,
1956 const Vector<double>& s,
1957 const unsigned& i) const
1958 {
1959 // Create an output stream
1960 std::ostringstream error_message_stream;
1961
1962 // Create an error message
1963 error_message_stream << "This interface doesn't make sense in "
1964 << "space-time elements!" << std::endl;
1965
1966 // Throw an error
1967 throw OomphLibError(error_message_stream.str(),
1968 OOMPH_CURRENT_FUNCTION,
1969 OOMPH_EXCEPTION_LOCATION);
1970 } // End of interpolated_u_nst
1971
1972
1973 /// Compute the derivatives of the i-th component of
1974 /// velocity at point s with respect
1975 /// to all data that can affect its value. In addition, return the global
1976 /// equation numbers corresponding to the data. The function is virtual
1977 /// so that it can be overloaded in the refineable version
1979 const unsigned& i,
1980 Vector<double>& du_ddata,
1981 Vector<unsigned>& global_eqn_number)
1982 {
1983 // Find the number of nodes
1984 unsigned n_node = nnode();
1985
1986 // Local shape function
1987 Shape psi(n_node);
1988
1989 // Find values of shape function
1990 shape(s, psi);
1991
1992 // Get nodal index at which i-th velocity is stored
1993 unsigned u_nodal_index = u_index_nst(i);
1994
1995 // Find the number of dofs associated with interpolated u
1996 unsigned n_u_dof = 0;
1997
1998 // Loop over the nodes
1999 for (unsigned l = 0; l < n_node; l++)
2000 {
2001 // Get the global equation number of the dof
2002 int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
2003
2004 // If it's positive add to the count
2005 if (global_eqn >= 0)
2006 {
2007 // Increment the counter
2008 n_u_dof++;
2009 }
2010 } // End of dinterpolated_u_nst_ddata
2011
2012 // Now resize the storage schemes
2013 du_ddata.resize(n_u_dof, 0.0);
2014 global_eqn_number.resize(n_u_dof, 0);
2015
2016 // Loop over th nodes again and set the derivatives
2017 unsigned count = 0;
2018 // Loop over the local nodes and sum
2019 for (unsigned l = 0; l < n_node; l++)
2020 {
2021 // Get the global equation number
2022 int global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
2023 if (global_eqn >= 0)
2024 {
2025 // Set the global equation number
2026 global_eqn_number[count] = global_eqn;
2027 // Set the derivative with respect to the unknown
2028 du_ddata[count] = psi[l];
2029 // Increase the counter
2030 ++count;
2031 }
2032 }
2033 } // End of dinterpolated_u_nst_ddata
2034
2035
2036 /// Return FE interpolated pressure at local coordinate s
2037 virtual double interpolated_p_nst(const Vector<double>& s) const
2038 {
2039 // Find number of nodes
2040 unsigned n_pres = npres_nst();
2041 // Local shape function
2042 Shape psi(n_pres);
2043 // Find values of shape function
2044 pshape_nst(s, psi);
2045
2046 // Initialise value of p
2047 double interpolated_p = 0.0;
2048 // Loop over the local nodes and sum
2049 for (unsigned l = 0; l < n_pres; l++)
2050 {
2051 interpolated_p += p_nst(l) * psi[l];
2052 }
2053
2054 return (interpolated_p);
2055 }
2056
2057
2058 /// Return FE interpolated pressure at local coordinate s at time level t
2059 double interpolated_p_nst(const unsigned& t, const Vector<double>& s) const
2060 {
2061 // Find number of nodes
2062 unsigned n_pres = npres_nst();
2063 // Local shape function
2064 Shape psi(n_pres);
2065 // Find values of shape function
2066 pshape_nst(s, psi);
2067
2068 // Initialise value of p
2069 double interpolated_p = 0.0;
2070 // Loop over the local nodes and sum
2071 for (unsigned l = 0; l < n_pres; l++)
2072 {
2073 interpolated_p += p_nst(t, l) * psi[l];
2074 }
2075
2076 return (interpolated_p);
2077 }
2078
2079
2080 /// Output solution in data vector at local cordinates s:
2081 /// x,y,z,u,v,p
2083 {
2084 // Resize data for values
2085 data.resize(2 * DIM + 2);
2086
2087 // Write values in the vector
2088 for (unsigned i = 0; i < DIM + 1; i++)
2089 {
2090 // Output the i-th (Eulerian) coordinate at these local coordinates
2091 data[i] = interpolated_x(s, i);
2092 }
2093
2094 // Write values in the vector
2095 for (unsigned i = 0; i < DIM; i++)
2096 {
2097 // Output the i-th velocity component at these local coordinates
2098 data[i + (DIM + 1)] = this->interpolated_u_nst(s, i);
2099 }
2100
2101 // Output the pressure field value at these local coordinates
2102 data[2 * DIM + 1] = this->interpolated_p_nst(s);
2103 } // End of point_output_data
2104 };
2105
2106 /// ///////////////////////////////////////////////////////////////////////////
2107 /// ///////////////////////////////////////////////////////////////////////////
2108 /// ///////////////////////////////////////////////////////////////////////////
2109
2110 //=======================================================================
2111 /// Taylor-Hood elements are Navier-Stokes elements with quadratic
2112 /// interpolation for velocities and positions and continuous linear
2113 /// pressure interpolation. They can be used within oomph-lib's
2114 /// block-preconditioning framework.
2115 //=======================================================================
2116 template<unsigned DIM>
2117 class QTaylorHoodSpaceTimeElement
2118 : public virtual QElement<DIM + 1, 3>,
2119 public virtual SpaceTimeNavierStokesEquations<DIM>
2120 {
2121 private:
2122 /// Static array of ints to hold number of variables at node
2123 static const unsigned Initial_Nvalue[];
2124
2125 protected:
2126 /// Static array of ints to hold conversion from pressure
2127 /// node numbers to actual node numbers
2128 static const unsigned Pconv[];
2129
2130
2131 /// Velocity shape and test functions and their derivs
2132 /// w.r.t. to global coords at local coordinate s (taken from geometry)
2133 /// Return Jacobian of mapping between local and global coordinates.
2135 Shape& psi,
2136 DShape& dpsidx,
2137 Shape& test,
2138 DShape& dtestdx) const;
2139
2140
2141 /// Velocity shape and test functions and their derivs
2142 /// w.r.t. to global coords at local coordinate s (taken from geometry)
2143 /// Return Jacobian of mapping between local and global coordinates.
2144 inline double dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
2145 Shape& psi,
2146 DShape& dpsidx,
2147 Shape& test,
2148 DShape& dtestdx) const;
2149
2150
2151 /// Shape/test functions and derivs w.r.t. to global coords at
2152 /// integration point ipt; return Jacobian of mapping (J). Also compute
2153 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
2155 const unsigned& ipt,
2156 Shape& psi,
2157 DShape& dpsidx,
2158 RankFourTensor<double>& d_dpsidx_dX,
2159 Shape& test,
2160 DShape& dtestdx,
2161 RankFourTensor<double>& d_dtestdx_dX,
2162 DenseMatrix<double>& djacobian_dX) const;
2163
2164
2165 /// Pressure shape functions and their derivs w.r.t. to global coords
2166 /// at local coordinate s (taken from geometry). Return Jacobian of mapping
2167 /// between local and global coordinates.
2168 inline double dpshape_eulerian(const Vector<double>& s,
2169 Shape& ppsi,
2170 DShape& dppsidx) const;
2171
2172
2173 /// Pressure test functions and their derivs w.r.t. to global coords
2174 /// at local coordinate s (taken from geometry). Return Jacobian of mapping
2175 /// between local and global coordinates.
2176 inline double dptest_eulerian(const Vector<double>& s,
2177 Shape& ptest,
2178 DShape& dptestdx) const;
2179
2180
2181 /// Pressure shape and test functions and their derivs
2182 /// w.r.t. to global coords at local coordinate s (taken from geometry).
2183 /// Return Jacobian of mapping between local and global coordinates.
2185 Shape& ppsi,
2186 DShape& dppsidx,
2187 Shape& ptest,
2188 DShape& dptestdx) const;
2189
2190 public:
2191 /// Constructor, no internal data points
2193 : QElement<DIM + 1, 3>(), SpaceTimeNavierStokesEquations<DIM>()
2194 {
2195 }
2196
2197 // Destructor: delete the integration scheme
2199 {
2200 // Delete the integral pointer
2201 delete this->integral_pt();
2202 } // End of ~QTaylorHoodSpaceTimeElement
2203
2204 /// Number of values (pinned or dofs) required at node n. Can
2205 /// be overwritten for hanging node version
2206 inline virtual unsigned required_nvalue(const unsigned& n) const
2207 {
2208 // Return the appropriate entry from Initial_Nvalue
2209 return Initial_Nvalue[n];
2210 } // End of required_nvalue
2211
2212
2213 /// Pressure shape functions at local coordinate s
2214 inline void pshape_nst(const Vector<double>& s, Shape& psi) const;
2215
2216
2217 /// Pressure test functions at local coordinate s
2218 inline void ptest_nst(const Vector<double>& s, Shape& psi) const;
2219
2220
2221 /// Pressure shape and test functions at local coordinte s
2222 inline void pshape_nst(const Vector<double>& s,
2223 Shape& psi,
2224 Shape& test) const;
2225
2226
2227 /// Set the value at which the pressure is stored in the nodes
2228 virtual int p_nodal_index_nst() const
2229 {
2230 // The pressure is stored straight after the velocity components
2231 return static_cast<int>(DIM);
2232 } // End of p_nodal_index_nst
2233
2234
2235 /// Return the local equation numbers for the pressure values.
2236 inline int p_local_eqn(const unsigned& n) const
2237 {
2238 // Get the local equation number associated with the n-th pressure dof
2239 return this->nodal_local_eqn(Pconv[n], p_nodal_index_nst());
2240 } // End of p_local_eqn
2241
2242
2243 /// Access function for the pressure values at local pressure
2244 /// node n_p (const version)
2245 double p_nst(const unsigned& n_p) const
2246 {
2247 // Get the nodal value associated with the n_p-th pressure dof
2248 return this->nodal_value(Pconv[n_p], this->p_nodal_index_nst());
2249 } // End of p_nst
2250
2251
2252 /// Access function for the pressure values at local pressure
2253 /// node n_p (const version)
2254 double p_nst(const unsigned& t, const unsigned& n_p) const
2255 {
2256 // Return the pressure value of the n_p-th pressure dof at time-level t
2257 return this->nodal_value(t, Pconv[n_p], this->p_nodal_index_nst());
2258 } // End of p_nst
2259
2260
2261 /// Return number of pressure values
2262 unsigned npres_nst() const
2263 {
2264 // There are 2^{DIM+1} pressure dofs where DIM is the spatial dimension
2265 // (rememebering that these are space-time elements)
2266 return static_cast<unsigned>(pow(2.0, static_cast<int>(DIM + 1)));
2267 } // End of npres_nst
2268
2269
2270 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
2271 void fix_pressure(const unsigned& p_dof, const double& p_value)
2272 {
2273 // Pin the pressure dof
2274 this->node_pt(Pconv[p_dof])->pin(this->p_nodal_index_nst());
2275
2276 // Now set its value
2277 this->node_pt(Pconv[p_dof])
2278 ->set_value(this->p_nodal_index_nst(), p_value);
2279 } // End of fix_pressure
2280
2281
2282 /// Build FaceElements that apply the Robin boundary condition
2283 /// to the pressure advection diffusion problem required by
2284 /// Fp preconditioner
2285 void build_fp_press_adv_diff_robin_bc_element(const unsigned& face_index)
2286 {
2287 // Create a new Robic BC element and add it to the storage
2290 QTaylorHoodSpaceTimeElement<DIM>>(this, face_index));
2291 } // End of build_fp_press_adv_diff_robin_bc_element
2292
2293
2294 /// Add to the set \c paired_load_data pairs containing
2295 /// - the pointer to a Data object
2296 /// and
2297 /// - the index of the value in that Data object
2298 /// .
2299 /// for all values (pressures, velocities) that affect the
2300 /// load computed in the \c get_load(...) function.
2302 std::set<std::pair<Data*, unsigned>>& paired_load_data);
2303
2304
2305 /// Add to the set \c paired_pressure_data pairs
2306 /// containing
2307 /// - the pointer to a Data object
2308 /// and
2309 /// - the index of the value in that Data object
2310 /// .
2311 /// for all pressure values that affect the
2312 /// load computed in the \c get_load(...) function.
2314 std::set<std::pair<Data*, unsigned>>& paired_pressure_data);
2315
2316
2317 /// Redirect output to NavierStokesEquations output
2318 void output(std::ostream& outfile)
2319 {
2320 // Call the base class implementation
2322 } // End of output
2323
2324
2325 /// Redirect output to NavierStokesEquations output
2326 void output(std::ostream& outfile, const unsigned& nplot)
2327 {
2328 // Call the base class implementation
2330 } // End of output
2331
2332
2333 /// Redirect output to NavierStokesEquations output
2334 void output(FILE* file_pt)
2335 {
2336 // Call the base class implementation
2338 } // End of output
2339
2340
2341 /// Redirect output to NavierStokesEquations output
2342 void output(FILE* file_pt, const unsigned& nplot)
2343 {
2344 // Call the base class implementation
2346 } // End of output
2347
2348 /*
2349 /// DRAIG: Delete...
2350 /// The number of "DOF types" that degrees of freedom in this element
2351 /// are sub-divided into (DIM velocities + 1 pressure)
2352 unsigned ndof_types() const
2353 {
2354 // Return the number of dof types being used in the mesh
2355 return DIM+1;
2356 } // End of ndof_types
2357
2358 /// Create a list of pairs for all unknowns in this element,
2359 /// so the first entry in each pair contains the global equation
2360 /// number of the unknown, while the second one contains the number
2361 /// of the "DOF type" that this unknown is associated with.
2362 /// (Function can obviously only be called if the equation numbering
2363 /// scheme has been set up.)
2364 ///
2365 /// The dof type enumeration (in 3D) is as follows:
2366 /// S_x = 0
2367 /// S_y = 1
2368 /// S_z = 2
2369 void get_dof_numbers_for_unknowns(
2370 std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const;
2371 */
2372 };
2373
2374 // Inline functions:
2375
2376 //==========================================================================
2377 /// Derivatives of the shape functions and test functions w.r.t to
2378 /// global (Eulerian) coordinates. Return Jacobian of mapping between
2379 /// local and global coordinates.
2380 //==========================================================================
2381 template<unsigned DIM>
2383 const Vector<double>& s,
2384 Shape& psi,
2385 DShape& dpsidx,
2386 Shape& test,
2387 DShape& dtestdx) const
2388 {
2389 //--------------------------
2390 // Call the shape functions:
2391 //--------------------------
2392 // Find the element dimension
2393 const unsigned el_dim = this->dim();
2394
2395 // Get the values of the shape functions and their local derivatives,
2396 // temporarily stored in dpsi
2397 this->dshape_local(s, psi, dpsidx);
2398
2399 // Allocate memory for the inverse jacobian
2400 DenseMatrix<double> inverse_jacobian(el_dim);
2401
2402 // Now calculate the inverse jacobian
2403 const double det =
2404 this->local_to_eulerian_mapping(dpsidx, inverse_jacobian);
2405
2406 // Now set the values of the derivatives to be dpsidx
2407 this->transform_derivatives(inverse_jacobian, dpsidx);
2408
2409 //-------------------------
2410 // Call the test functions:
2411 //-------------------------
2412 // Make sure we're using 3D elements
2413 if (el_dim != 3)
2414 {
2415 // Create an output stream
2416 std::ostringstream error_message_stream;
2417
2418 // Create an error message
2419 error_message_stream << "Need 3D space-time elements for this to work!"
2420 << std::endl;
2421
2422 // Throw the error message
2423 throw OomphLibError(error_message_stream.str(),
2424 OOMPH_CURRENT_FUNCTION,
2425 OOMPH_EXCEPTION_LOCATION);
2426 }
2427
2428 //--------start_of_dshape_local--------------------------------------
2429 // Local storage
2430 double test_values[3][3];
2431 double dtest_values[3][3];
2432
2433 // Index of the total shape function
2434 unsigned index = 0;
2435
2436 // Call the 1D shape functions and derivatives
2437 OneDimLagrange::shape<3>(s[0], test_values[0]);
2438 OneDimLagrange::shape<3>(s[1], test_values[1]);
2439 OneDimLagrange::dshape<3>(s[0], dtest_values[0]);
2440 OneDimLagrange::dshape<3>(s[1], dtest_values[1]);
2441
2442 // Set the time discretisation
2443 // OneDimLagrange::shape<3>(s[2],test_values[2]);
2444 // OneDimLagrange::dshape<3>(s[2],dtest_values[2]);
2445 OneDimDiscontinuousGalerkin::shape<3>(s[2], test_values[2]);
2446 OneDimDiscontinuousGalerkin::dshape<3>(s[2], dtest_values[2]);
2447
2448 // Loop over the nodes in the third local coordinate direction
2449 for (unsigned k = 0; k < 3; k++)
2450 {
2451 // Loop over the nodes in the second local coordinate direction
2452 for (unsigned j = 0; j < 3; j++)
2453 {
2454 // Loop over the nodes in the first local coordinate direction
2455 for (unsigned i = 0; i < 3; i++)
2456 {
2457 // Calculate dtest/ds_0
2458 dtestdx(index, 0) =
2459 dtest_values[0][i] * test_values[1][j] * test_values[2][k];
2460
2461 // Calculate dtest/ds_1
2462 dtestdx(index, 1) =
2463 test_values[0][i] * dtest_values[1][j] * test_values[2][k];
2464
2465 // Calculate dtest/ds_2
2466 dtestdx(index, 2) =
2467 test_values[0][i] * test_values[1][j] * dtest_values[2][k];
2468
2469 // Calculate the index-th entry of test
2470 test[index] =
2471 test_values[0][i] * test_values[1][j] * test_values[2][k];
2472
2473 // Increment the index
2474 index++;
2475 }
2476 } // for (unsigned j=0;j<3;j++)
2477 } // for (unsigned k=0;k<3;k++)
2478 //--------end_of_dshape_local----------------------------------------
2479
2480 // Transform derivatives from dtest/ds to dtest/dx
2481 this->transform_derivatives(inverse_jacobian, dtestdx);
2482
2483 // Return the determinant value
2484 return det;
2485 } // End of dshape_and_dtest_eulerian_nst
2486
2487 //==========================================================================
2488 /// Derivatives of the shape functions and test functions w.r.t to
2489 /// global (Eulerian) coordinates. Return Jacobian of mapping between
2490 /// local and global coordinates.
2491 //==========================================================================
2492 template<unsigned DIM>
2493 inline double QTaylorHoodSpaceTimeElement<
2494 DIM>::dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
2495 Shape& psi,
2496 DShape& dpsidx,
2497 Shape& test,
2498 DShape& dtestdx) const
2499 {
2500 // Calculate the element dimension
2501 const unsigned el_dim = DIM + 1;
2502
2503 // Storage for the local coordinates of the integration point
2504 Vector<double> s(el_dim, 0.0);
2505
2506 // Set the local coordinate
2507 for (unsigned i = 0; i < el_dim; i++)
2508 {
2509 // Calculate the i-th local coordinate at the ipt-th knot point
2510 s[i] = this->integral_pt()->knot(ipt, i);
2511 }
2512
2513 // Return the Jacobian of the geometrical shape functions and derivatives
2514 return dshape_and_dtest_eulerian_nst(s, psi, dpsidx, test, dtestdx);
2515 } // End of dshape_and_dtest_eulerian_at_knot_nst
2516
2517
2518 //==========================================================================
2519 /// 2D (in space): Pressure shape and test functions and derivs w.r.t. to
2520 /// Eulerian coordinates. Return Jacobian of mapping between local and
2521 /// global coordinates.
2522 //==========================================================================
2523 template<>
2525 const Vector<double>& s, Shape& ppsi, DShape& dppsidx) const
2526 {
2527 // Local storage for the shape function (x-direction)
2528 double psi1[2];
2529
2530 // Local storage for the shape function (y-direction)
2531 double psi2[2];
2532
2533 // Local storage for the shape function (z-direction)
2534 double psi3[2];
2535
2536 // Local storage for the shape function derivatives (x-direction)
2537 double dpsi1[2];
2538
2539 // Local storage for the test function derivatives (y-direction)
2540 double dpsi2[2];
2541
2542 // Local storage for the test function derivatives (z-direction)
2543 double dpsi3[2];
2544
2545 // Call the OneDimensional Shape functions
2546 OneDimLagrange::shape<2>(s[0], psi1);
2547
2548 // Call the OneDimensional Shape functions
2549 OneDimLagrange::shape<2>(s[1], psi2);
2550
2551 // Call the OneDimensional Shape functions
2552 OneDimLagrange::shape<2>(s[2], psi3);
2553
2554 // Call the OneDimensional Shape functions
2555 OneDimLagrange::dshape<2>(s[0], dpsi1);
2556
2557 // Call the OneDimensional Shape functions
2558 OneDimLagrange::dshape<2>(s[1], dpsi2);
2559
2560 // Call the OneDimensional Shape functions
2561 OneDimLagrange::dshape<2>(s[2], dpsi3);
2562
2563 //--------------------------------------------------------------------
2564 // Now let's loop over the nodal points in the element with s1 being
2565 // the "x" coordinate, s2 being the "y" coordinate and s3 being the
2566 // "z" coordinate:
2567 //--------------------------------------------------------------------
2568 // Loop over the points in the z-direction
2569 for (unsigned i = 0; i < 2; i++)
2570 {
2571 // Loop over the points in the y-direction
2572 for (unsigned j = 0; j < 2; j++)
2573 {
2574 // Loop over the points in the x-direction
2575 for (unsigned k = 0; k < 2; k++)
2576 {
2577 // Multiply the three 1D functions together to get the 3D function
2578 ppsi[4 * i + 2 * j + k] = psi3[i] * psi2[j] * psi1[k];
2579
2580 // Multiply the appropriate shape and shape function derivatives
2581 // together
2582 dppsidx(4 * i + 2 * j + k, 0) = psi3[i] * psi2[j] * dpsi1[k];
2583
2584 // Multiply the appropriate shape and shape function derivatives
2585 // together
2586 dppsidx(4 * i + 2 * j + k, 1) = psi3[i] * dpsi2[j] * psi1[k];
2587
2588 // Multiply the appropriate shape and shape function derivatives
2589 // together
2590 dppsidx(4 * i + 2 * j + k, 2) = dpsi3[i] * psi2[j] * psi1[k];
2591 }
2592 } // for (unsigned j=0;j<2;j++)
2593 } // for (unsigned i=0;i<2;i++)
2594
2595 // Allocate space for the geometrical shape functions
2596 Shape psi(27);
2597
2598 // Allocate space for the geometrical shape function derivatives
2599 DShape dpsi(27, 3);
2600
2601 // Get the values of the shape functions and their derivatives
2602 dshape_local(s, psi, dpsi);
2603
2604 // Allocate memory for the 3x3 inverse jacobian
2605 DenseMatrix<double> inverse_jacobian(3);
2606
2607 // Now calculate the inverse jacobian
2608 const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
2609
2610 // Now set the values of the derivatives to be derivatives w.r.t. to
2611 // the Eulerian coordinates
2612 transform_derivatives(inverse_jacobian, dppsidx);
2613
2614 // Return the determinant of the jacobian
2615 return det;
2616 } // End of dpshape_eulerian
2617
2618
2619 //==========================================================================
2620 /// 2D (in space): Pressure shape and test functions and derivs w.r.t. to
2621 /// Eulerian coordinates. Return Jacobian of mapping between local and
2622 /// global coordinates.
2623 //==========================================================================
2624 template<>
2626 const Vector<double>& s, Shape& ptest, DShape& dptestdx) const
2627 {
2628 // Local storage for the shape function (x-direction)
2629 double test1[2];
2630
2631 // Local storage for the shape function (y-direction)
2632 double test2[2];
2633
2634 // Local storage for the shape function (z-direction)
2635 double test3[2];
2636
2637 // Local storage for the shape function derivatives (x-direction)
2638 double dtest1[2];
2639
2640 // Local storage for the test function derivatives (y-direction)
2641 double dtest2[2];
2642
2643 // Local storage for the test function derivatives (z-direction)
2644 double dtest3[2];
2645
2646 // Call the OneDimensional Shape functions
2647 OneDimLagrange::shape<2>(s[0], test1);
2648
2649 // Call the OneDimensional Shape functions
2650 OneDimLagrange::shape<2>(s[1], test2);
2651
2652 // Call the OneDimensional Shape functions
2653 // OneDimLagrange::shape<2>(s[2],test3);
2655
2656 // Call the OneDimensional Shape functions
2657 OneDimLagrange::dshape<2>(s[0], dtest1);
2658
2659 // Call the OneDimensional Shape functions
2660 OneDimLagrange::dshape<2>(s[1], dtest2);
2661
2662 // Call the OneDimensional Shape functions
2663 // OneDimLagrange::dshape<2>(s[2],dtest3);
2665
2666 //--------------------------------------------------------------------
2667 // Now let's loop over the nodal points in the element with s1 being
2668 // the "x" coordinate, s2 being the "y" coordinate and s3 being the
2669 // "z" coordinate:
2670 //--------------------------------------------------------------------
2671 // Loop over the points in the z-direction
2672 for (unsigned i = 0; i < 2; i++)
2673 {
2674 // Loop over the points in the y-direction
2675 for (unsigned j = 0; j < 2; j++)
2676 {
2677 // Loop over the points in the x-direction
2678 for (unsigned k = 0; k < 2; k++)
2679 {
2680 // Multiply the three 1D functions together to get the 3D function
2681 ptest[4 * i + 2 * j + k] = test3[i] * test2[j] * test1[k];
2682
2683 // Multiply the appropriate shape and shape function derivatives
2684 // together
2685 dptestdx(4 * i + 2 * j + k, 0) = test3[i] * test2[j] * dtest1[k];
2686
2687 // Multiply the appropriate shape and shape function derivatives
2688 // together
2689 dptestdx(4 * i + 2 * j + k, 1) = test3[i] * dtest2[j] * test1[k];
2690
2691 // Multiply the appropriate shape and shape function derivatives
2692 // together
2693 dptestdx(4 * i + 2 * j + k, 2) = dtest3[i] * test2[j] * test1[k];
2694 }
2695 } // for (unsigned j=0;j<2;j++)
2696 } // for (unsigned i=0;i<2;i++)
2697
2698 // Allocate space for the geometrical shape functions
2699 Shape psi(27);
2700
2701 // Allocate space for the geometrical shape function derivatives
2702 DShape dpsi(27, 3);
2703
2704 // Get the values of the shape functions and their derivatives
2705 dshape_local(s, psi, dpsi);
2706
2707 // Allocate memory for the 3x3 inverse jacobian
2708 DenseMatrix<double> inverse_jacobian(3);
2709
2710 // Now calculate the inverse jacobian
2711 const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
2712
2713 // Now set the values of the derivatives to be derivatives w.r.t. to
2714 // the Eulerian coordinates
2715 transform_derivatives(inverse_jacobian, dptestdx);
2716
2717 // Return the determinant of the jacobian
2718 return det;
2719 } // End of dptest_eulerian
2720
2721
2722 //==========================================================================
2723 /// 2D (in space): Pressure shape and test functions and derivs w.r.t. to
2724 /// Eulerian coordinates. Return Jacobian of mapping between local and
2725 /// global coordinates.
2726 //==========================================================================
2727 template<>
2729 const Vector<double>& s,
2730 Shape& ppsi,
2731 DShape& dppsidx,
2732 Shape& ptest,
2733 DShape& dptestdx) const
2734 {
2735 // Call the test functions and derivatives
2736 dptest_eulerian(s, ptest, dptestdx);
2737
2738 // Call the shape functions and derivatives and the Jacobian of the mapping
2739 return this->dpshape_eulerian(s, ppsi, dppsidx);
2740 } // End of dpshape_and_dptest_eulerian_nst
2741
2742
2743 //==========================================================================
2744 /// 2D (in space):
2745 /// Define the shape functions (psi) and test functions (test) and
2746 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
2747 /// and return Jacobian of mapping (J). Additionally compute the
2748 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
2749 ///
2750 /// Galerkin: Test functions = shape functions
2751 //==========================================================================
2752 template<>
2755 const unsigned& ipt,
2756 Shape& psi,
2757 DShape& dpsidx,
2758 RankFourTensor<double>& d_dpsidx_dX,
2759 Shape& test,
2760 DShape& dtestdx,
2761 RankFourTensor<double>& d_dtestdx_dX,
2762 DenseMatrix<double>& djacobian_dX) const
2763 {
2764 // Call the geometrical shape functions and derivatives
2765 const double J = this->dshape_eulerian_at_knot(
2766 ipt, psi, dpsidx, djacobian_dX, d_dpsidx_dX);
2767
2768 // DRAIG: Delete!
2769 throw OomphLibError("Hasn't been implemented properly yet!",
2770 OOMPH_CURRENT_FUNCTION,
2771 OOMPH_EXCEPTION_LOCATION);
2772
2773 // Loop over the test functions and derivatives
2774 for (unsigned i = 0; i < 27; i++)
2775 {
2776 // The test functions are the same as the shape functions
2777 test[i] = psi[i];
2778
2779 // Loop over the spatial derivatives
2780 for (unsigned k = 0; k < 3; k++)
2781 {
2782 // Set the test function derivatives to the shape function derivatives
2783 dtestdx(i, k) = dpsidx(i, k);
2784
2785 // Loop over the dimensions
2786 for (unsigned p = 0; p < 3; p++)
2787 {
2788 // Loop over test functions
2789 for (unsigned q = 0; q < 27; q++)
2790 {
2791 // Set the test function derivatives to the shape function
2792 // derivatives
2793 d_dtestdx_dX(p, q, i, k) = d_dpsidx_dX(p, q, i, k);
2794 }
2795 } // for (unsigned p=0;p<3;p++)
2796 } // for (unsigned k=0;k<3;k++)
2797 } // for (unsigned i=0;i<27;i++)
2798
2799 // Return the jacobian
2800 return J;
2801 } // End of dshape_and_dtest_eulerian_at_knot_nst
2802
2803
2804 //==========================================================================
2805 /// 2D (in space): Pressure shape functions
2806 //==========================================================================
2807 template<>
2809 const Vector<double>& s, Shape& psi) const
2810 {
2811 // Local storage for the shape function value (in the x-direction)
2812 double psi1[2];
2813
2814 // Local storage for the shape function value (in the y-direction)
2815 double psi2[2];
2816
2817 // Local storage for the shape function value (in the z-direction)
2818 double psi3[2];
2819
2820 // Call the OneDimensional Shape functions
2821 OneDimLagrange::shape<2>(s[0], psi1);
2822
2823 // Call the OneDimensional Shape functions
2824 OneDimLagrange::shape<2>(s[1], psi2);
2825
2826 // Call the OneDimensional Shape functions
2827 OneDimLagrange::shape<2>(s[2], psi3);
2828
2829 //--------------------------------------------------------------------
2830 // Now let's loop over the nodal points in the element with s1 being
2831 // the "x" coordinate, s2 being the "y" coordinate and s3 being the
2832 // "z" coordinate:
2833 //--------------------------------------------------------------------
2834 // Loop over the points in the z-direction
2835 for (unsigned i = 0; i < 2; i++)
2836 {
2837 // Loop over the points in the y-direction
2838 for (unsigned j = 0; j < 2; j++)
2839 {
2840 // Loop over the points in the x-direction
2841 for (unsigned k = 0; k < 2; k++)
2842 {
2843 // Multiply the three 1D functions together to get the 3D function
2844 psi[4 * i + 2 * j + k] = psi3[i] * psi2[j] * psi1[k];
2845 }
2846 } // for (unsigned j=0;j<2;j++)
2847 } // for (unsigned i=0;i<2;i++)
2848 } // End of pshape_nst
2849
2850
2851 //==========================================================================
2852 /// 2D (in space): Pressure shape functions
2853 //==========================================================================
2854 template<>
2856 Shape& test) const
2857 {
2858 // Local storage for the shape function value (in the x-direction)
2859 double test1[2];
2860
2861 // Local storage for the shape function value (in the y-direction)
2862 double test2[2];
2863
2864 // Local storage for the shape function value (in the z-direction)
2865 double test3[2];
2866
2867 // Call the OneDimensional Shape functions
2868 OneDimLagrange::shape<2>(s[0], test1);
2869
2870 // Call the OneDimensional Shape functions
2871 OneDimLagrange::shape<2>(s[1], test2);
2872
2873 // Call the OneDimensional Shape functions
2875
2876 //--------------------------------------------------------------------
2877 // Now let's loop over the nodal points in the element with s1 being
2878 // the "x" coordinate, s2 being the "y" coordinate and s3 being the
2879 // "z" coordinate:
2880 //--------------------------------------------------------------------
2881 // Loop over the points in the z-direction
2882 for (unsigned i = 0; i < 2; i++)
2883 {
2884 // Loop over the points in the y-direction
2885 for (unsigned j = 0; j < 2; j++)
2886 {
2887 // Loop over the points in the x-direction
2888 for (unsigned k = 0; k < 2; k++)
2889 {
2890 // Multiply the three 1D functions together to get the 3D function
2891 test[4 * i + 2 * j + k] = test3[i] * test2[j] * test1[k];
2892 }
2893 } // for (unsigned j=0;j<2;j++)
2894 } // for (unsigned i=0;i<2;i++)
2895 } // End of ptest_nst
2896
2897
2898 //==========================================================================
2899 /// Pressure shape and test functions
2900 //==========================================================================
2901 template<unsigned DIM>
2903 const Vector<double>& s, Shape& psi, Shape& test) const
2904 {
2905 // Call the pressure shape functions
2906 pshape_nst(s, psi);
2907
2908 // Call the pressure test functions
2909 ptest_nst(s, test);
2910 } // End of pshape_nst
2911
2912
2913 //=======================================================================
2914 /// Face geometry of the 2D Taylor_Hood elements
2915 //=======================================================================
2916 template<>
2917 class FaceGeometry<QTaylorHoodSpaceTimeElement<2>>
2918 : public virtual QElement<2, 3>
2919 {
2920 public:
2921 /// Constructor; empty
2922 FaceGeometry() : QElement<2, 3>() {}
2923 };
2924
2925 //=======================================================================
2926 /// Face geometry of the FaceGeometry of the 2D Taylor Hoodelements
2927 //=======================================================================
2928 template<>
2929 class FaceGeometry<FaceGeometry<QTaylorHoodSpaceTimeElement<2>>>
2930 : public virtual QElement<1, 3>
2931 {
2932 public:
2933 /// Constructor; empty
2934 FaceGeometry() : QElement<1, 3>() {}
2935 };
2936
2937 /// /////////////////////////////////////////////////////////////////
2938 /// /////////////////////////////////////////////////////////////////
2939 /// /////////////////////////////////////////////////////////////////
2940
2941 //==========================================================
2942 /// Taylor Hood upgraded to become projectable
2943 //==========================================================
2944 template<class TAYLOR_HOOD_ELEMENT>
2945 class ProjectableTaylorHoodSpaceTimeElement
2946 : public virtual ProjectableElement<TAYLOR_HOOD_ELEMENT>
2947 {
2948 public:
2949 /// Constructor [this was only required explicitly
2950 /// from gcc 4.5.2 onwards...]
2952
2953
2954 /// Specify the values associated with field fld.
2955 /// The information is returned in a vector of pairs which comprise
2956 /// the Data object and the value within it, that correspond to field fld.
2957 /// In the underlying Taylor Hood elements the fld-th velocities are stored
2958 /// at the fld-th value of the nodes; the pressures (the dim-th
2959 /// field) are the dim-th values at the vertex nodes etc.
2961 {
2962 // Create the vector
2964
2965 // If we're dealing with the velocity dofs
2966 if (fld < this->dim() - 1)
2967 {
2968 // How many nodes in the element?
2969 unsigned nnod = this->nnode();
2970
2971 // Loop over all nodes
2972 for (unsigned j = 0; j < nnod; j++)
2973 {
2974 // Add the data value associated with the velocity components
2975 data_values.push_back(std::make_pair(this->node_pt(j), fld));
2976 }
2977 }
2978 // If we're dealing with the pressure dof
2979 else
2980 {
2981 // How many pressure dofs are there?
2982 // DRAIG: Shouldn't there be more?
2983 unsigned Pconv_size = this->dim();
2984
2985 // Loop over all vertex nodes
2986 for (unsigned j = 0; j < Pconv_size; j++)
2987 {
2988 // Get the vertex index associated with the j-th pressure dof
2989 unsigned vertex_index = this->Pconv[j];
2990
2991 // Add the data value associated with the pressure components
2992 data_values.push_back(
2993 std::make_pair(this->node_pt(vertex_index), fld));
2994 }
2995 } // if (fld<this->dim())
2996
2997 // Return the vector
2998 return data_values;
2999 } // End of data_values_of_field
3000
3001
3002 /// Number of fields to be projected: dim+1, corresponding to
3003 /// velocity components and pressure
3005 {
3006 // There are dim velocity dofs and 1 pressure dof
3007 return this->dim();
3008 } // End of nfields_for_projection
3009
3010
3011 /// Number of history values to be stored for fld-th field. Whatever
3012 /// the timestepper has set up for the velocity components and none for
3013 /// the pressure field. NOTE: The count includes the current value!
3014 unsigned nhistory_values_for_projection(const unsigned& fld)
3015 {
3016 // If we're dealing with the pressure dof
3017 if (fld == this->dim())
3018 {
3019 // Pressure doesn't have history values
3020 return this->node_pt(0)->ntstorage();
3021 }
3022 // If we're dealing with the velocity dofs
3023 else
3024 {
3025 // The velocity dofs have ntstorage() history values
3026 return this->node_pt(0)->ntstorage();
3027 }
3028 } // End of nhistory_values_for_projection
3029
3030
3031 /// Number of positional history values
3032 /// (Note: count includes current value!)
3034 {
3035 // Return the number of positional history values
3036 return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
3037 } // End of nhistory_values_for_coordinate_projection
3038
3039
3040 /// Return the Jacobian of the mapping and the shape functions
3041 /// of field fld at local coordinate s
3042 double jacobian_and_shape_of_field(const unsigned& fld,
3043 const Vector<double>& s,
3044 Shape& psi)
3045 {
3046 // How many dimensions in this element?
3047 unsigned n_dim = this->dim();
3048
3049 // How many nodes in this element?
3050 unsigned n_node = this->nnode();
3051
3052 // If we're on the pressure dof
3053 if (fld == n_dim)
3054 {
3055 // Call the pressure interpolation function
3056 this->pshape_nst(s, psi);
3057
3058 // Allocate space for the pressure shape function
3059 Shape psif(n_node);
3060
3061 // Allocate space for the pressure test function
3062 Shape testf(n_node);
3063
3064 // Allocate space for the pressure shape function derivatives
3065 DShape dpsifdx(n_node, n_dim);
3066
3067 // Allocate space for the pressure test function derivatives
3068 DShape dtestfdx(n_node, n_dim);
3069
3070 // Calculate the Jacobian of the mapping
3071 double J = this->dshape_and_dtest_eulerian_nst(
3072 s, psif, dpsifdx, testf, dtestfdx);
3073
3074 // Return the Jacobian
3075 return J;
3076 }
3077 // If we're on the velocity components
3078 else
3079 {
3080 // Allocate space for the test functions
3081 Shape testf(n_node);
3082
3083 // Allocate space for the shape function derivatives
3084 DShape dpsifdx(n_node, n_dim);
3085
3086 // Allocate space for the test function derivatives
3087 DShape dtestfdx(n_node, n_dim);
3088
3089 // Calculate the Jacobian of the mapping
3090 double J =
3091 this->dshape_and_dtest_eulerian_nst(s, psi, dpsifdx, testf, dtestfdx);
3092
3093 // Return the Jacobian
3094 return J;
3095 }
3096 } // End of jacobian_and_shape_of_field
3097
3098
3099 /// Return interpolated field fld at local coordinate s, at time
3100 /// level t (t=0: present; t>0: history values)
3101 double get_field(const unsigned& t,
3102 const unsigned& fld,
3103 const Vector<double>& s)
3104 {
3105 unsigned n_dim = this->dim();
3106 unsigned n_node = this->nnode();
3107
3108 // If fld=n_dim, we deal with the pressure
3109 if (fld == n_dim)
3110 {
3111 return this->interpolated_p_nst(t, s);
3112 }
3113 // Velocity
3114 else
3115 {
3116 // Find the index at which the variable is stored
3117 unsigned u_nodal_index = this->u_index_nst(fld);
3118
3119 // Local shape function
3120 Shape psi(n_node);
3121
3122 // Find values of shape function
3123 this->shape(s, psi);
3124
3125 // Initialise value of u
3126 double interpolated_u = 0.0;
3127
3128 // Sum over the local nodes at that time
3129 for (unsigned l = 0; l < n_node; l++)
3130 {
3131 interpolated_u += this->nodal_value(t, l, u_nodal_index) * psi[l];
3132 }
3133 return interpolated_u;
3134 }
3135 }
3136
3137
3138 /// Return number of values in field fld
3139 unsigned nvalue_of_field(const unsigned& fld)
3140 {
3141 if (fld == this->dim())
3142 {
3143 return this->npres_nst();
3144 }
3145 else
3146 {
3147 return this->nnode();
3148 }
3149 }
3150
3151
3152 /// Return local equation number of value j in field fld.
3153 int local_equation(const unsigned& fld, const unsigned& j)
3154 {
3155 if (fld == this->dim())
3156 {
3157 return this->p_local_eqn(j);
3158 }
3159 else
3160 {
3161 const unsigned u_nodal_index = this->u_index_nst(fld);
3162 return this->nodal_local_eqn(j, u_nodal_index);
3163 }
3164 }
3165 };
3166
3167
3168 //=======================================================================
3169 /// Face geometry for element is the same as that for the underlying
3170 /// wrapped element
3171 //=======================================================================
3172 template<class ELEMENT>
3173 class FaceGeometry<ProjectableTaylorHoodSpaceTimeElement<ELEMENT>>
3174 : public virtual FaceGeometry<ELEMENT>
3175 {
3176 public:
3177 FaceGeometry() : FaceGeometry<ELEMENT>() {}
3178 };
3179
3180 //=======================================================================
3181 /// Face geometry of the Face Geometry for element is the same as
3182 /// that for the underlying wrapped element
3183 //=======================================================================
3184 template<class ELEMENT>
3185 class FaceGeometry<
3186 FaceGeometry<ProjectableTaylorHoodSpaceTimeElement<ELEMENT>>>
3187 : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
3188 {
3189 public:
3191 };
3192
3193} // End of namespace oomph
3194
3195#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 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 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.
double dptest_eulerian(const Vector< double > &s, Shape &ptest, DShape &dptestdx) const
Pressure test functions and their derivs w.r.t. to global coords at local coordinate s (taken from ge...
double p_nst(const unsigned &t, const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version)
double dpshape_eulerian(const Vector< double > &s, Shape &ppsi, DShape &dppsidx) const
Pressure shape functions and their derivs w.r.t. to global coords at local coordinate s (taken from g...
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...
void ptest_nst(const Vector< double > &s, Shape &psi) const
Pressure test functions at local coordinate s.
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.
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...
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....
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.
double * st_pt() const
Pointer to Strouhal parameter (const. version)
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 ...
unsigned n_u_nst() const
Return the number of velocity components (used in FluidInterfaceElements)
double * St_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
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.
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.
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.
bool is_strouhal_stored_as_external_data() const
Are we storing the Strouhal number as external data?
NavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
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...
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,...
double re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
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 ...
virtual double dptest_eulerian(const Vector< double > &s, Shape &ptest, DShape &dptestdx) const =0
Pressure test functions and their derivs w.r.t. to global coords at local coordinate s (taken from ge...
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....
const double & st() const
Strouhal parameter (const. version)
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...
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.
virtual double dpshape_eulerian(const Vector< double > &s, Shape &ppsi, DShape &dppsidx) const =0
Pressure shape functions and their derivs w.r.t. to global coords at local coordinate s (taken from g...
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 Strouhal_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 *& st_pt()
Pointer to Strouhal number (can only assign to private member data)
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 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 dshape< 3 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:820
void shape< 3 >(const double &s, double *Psi)
1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:810
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< 3 >(const double &s, double *DPsi)
Derivatives of 1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:645
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< 3 >(const double &s, double *Psi)
1D shape functions specialised to quadratic order (3 Nodes)
Definition: shape.h:635
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).
//////////////////////////////////////////////////////////////////// ////////////////////////////////...