discontinuous_galerkin_equal_order_pressure_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_EQUAL_ORDER_PRESSURE_SPACETIME_NAVIER_STOKES_ELEMENTS_HEADER
28#define OOMPH_DISCONTINUOUS_GALERKIN_EQUAL_ORDER_PRESSURE_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 //======================================================================
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 //======================================================================
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)
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)
506
507 /// Pointer to the density ratio (relative to the density used in the
508 /// definition of the Reynolds number)
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.
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
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)
1257 void strain_rate(const Vector<double>& s,
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.
1263 void get_traction(const Vector<double>& s,
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.
1271 void get_traction(const Vector<double>& s,
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
1554 void compute_norm(Vector<double>& norm);
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>
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 /// Number of values (pinned or dofs) required at node n. Can
2198 /// be overwritten for hanging node version
2199 inline virtual unsigned required_nvalue(const unsigned& n) const
2200 {
2201 // Return the appropriate entry from Initial_Nvalue
2202 return Initial_Nvalue[n];
2203 } // End of required_nvalue
2204
2205
2206 /// Pressure shape functions at local coordinate s
2207 inline void pshape_nst(const Vector<double>& s, Shape& psi) const;
2208
2209
2210 /// Pressure test functions at local coordinate s
2211 inline void ptest_nst(const Vector<double>& s, Shape& psi) const;
2212
2213
2214 /// Pressure shape and test functions at local coordinte s
2215 inline void pshape_nst(const Vector<double>& s,
2216 Shape& psi,
2217 Shape& test) const;
2218
2219
2220 /// Set the value at which the pressure is stored in the nodes
2221 virtual int p_nodal_index_nst() const
2222 {
2223 // The pressure is stored straight after the velocity components
2224 return static_cast<int>(DIM);
2225 } // End of p_nodal_index_nst
2226
2227
2228 /// Return the local equation numbers for the pressure values.
2229 inline int p_local_eqn(const unsigned& n) const
2230 {
2231 // Get the local equation number associated with the n-th pressure dof
2232 return this->nodal_local_eqn(Pconv[n], p_nodal_index_nst());
2233 } // End of p_local_eqn
2234
2235
2236 /// Access function for the pressure values at local pressure
2237 /// node n_p (const version)
2238 double p_nst(const unsigned& n_p) const
2239 {
2240 // Get the nodal value associated with the n_p-th pressure dof
2241 return this->nodal_value(Pconv[n_p], this->p_nodal_index_nst());
2242 } // End of p_nst
2243
2244
2245 /// Access function for the pressure values at local pressure
2246 /// node n_p (const version)
2247 double p_nst(const unsigned& t, const unsigned& n_p) const
2248 {
2249 // Return the pressure value of the n_p-th pressure dof at time-level t
2250 return this->nodal_value(t, Pconv[n_p], this->p_nodal_index_nst());
2251 } // End of p_nst
2252
2253
2254 /// Return number of pressure values
2255 unsigned npres_nst() const
2256 {
2257 // There are 3*(2^{DIM}) pressure dofs where DIM is the spatial dimension
2258 // (remembering that these are space-time elements)
2259 return 3.0 * static_cast<unsigned>(pow(2.0, static_cast<int>(DIM)));
2260 } // End of npres_nst
2261
2262
2263 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
2264 void fix_pressure(const unsigned& p_dof, const double& p_value)
2265 {
2266 // Pin the pressure dof
2267 this->node_pt(Pconv[p_dof])->pin(this->p_nodal_index_nst());
2268
2269 // Now set its value
2270 this->node_pt(Pconv[p_dof])
2271 ->set_value(this->p_nodal_index_nst(), p_value);
2272 } // End of fix_pressure
2273
2274
2275 /// Build FaceElements that apply the Robin boundary condition
2276 /// to the pressure advection diffusion problem required by
2277 /// Fp preconditioner
2278 void build_fp_press_adv_diff_robin_bc_element(const unsigned& face_index)
2279 {
2280 // Create a new Robic BC element and add it to the storage
2283 QTaylorHoodSpaceTimeElement<DIM>>(this, face_index));
2284 } // End of build_fp_press_adv_diff_robin_bc_element
2285
2286
2287 /// Add to the set \c paired_load_data pairs containing
2288 /// - the pointer to a Data object
2289 /// and
2290 /// - the index of the value in that Data object
2291 /// .
2292 /// for all values (pressures, velocities) that affect the
2293 /// load computed in the \c get_load(...) function.
2294 void identify_load_data(
2295 std::set<std::pair<Data*, unsigned>>& paired_load_data);
2296
2297
2298 /// Add to the set \c paired_pressure_data pairs
2299 /// containing
2300 /// - the pointer to a Data object
2301 /// and
2302 /// - the index of the value in that Data object
2303 /// .
2304 /// for all pressure values that affect the
2305 /// load computed in the \c get_load(...) function.
2307 std::set<std::pair<Data*, unsigned>>& paired_pressure_data);
2308
2309
2310 /// Redirect output to NavierStokesEquations output
2311 void output(std::ostream& outfile)
2312 {
2313 // Call the base class implementation
2315 } // End of output
2316
2317
2318 /// Redirect output to NavierStokesEquations output
2319 void output(std::ostream& outfile, const unsigned& nplot)
2320 {
2321 // Call the base class implementation
2323 } // End of output
2324
2325
2326 /// Redirect output to NavierStokesEquations output
2327 void output(FILE* file_pt)
2328 {
2329 // Call the base class implementation
2331 } // End of output
2332
2333
2334 /// Redirect output to NavierStokesEquations output
2335 void output(FILE* file_pt, const unsigned& nplot)
2336 {
2337 // Call the base class implementation
2339 } // End of output
2340
2341
2342 /// Returns the number of "DOF types" that degrees of freedom
2343 /// in this element are sub-divided into: Velocity and pressure.
2344 unsigned ndof_types() const
2345 {
2346 // There are DIM velocity components and 1 pressure component
2347 return DIM + 1;
2348 } // End of ndof_types
2349
2350
2351 /// Create a list of pairs for all unknowns in this element,
2352 /// so that the first entry in each pair contains the global equation
2353 /// number of the unknown, while the second one contains the number
2354 /// of the "DOF type" that this unknown is associated with.
2355 /// (Function can obviously only be called if the equation numbering
2356 /// scheme has been set up.) Velocity=0; Pressure=1
2358 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const;
2359 };
2360
2361 // Inline functions:
2362
2363 //==========================================================================
2364 /// Derivatives of the shape functions and test functions w.r.t to
2365 /// global (Eulerian) coordinates. Return Jacobian of mapping between
2366 /// local and global coordinates.
2367 //==========================================================================
2368 template<unsigned DIM>
2370 const Vector<double>& s,
2371 Shape& psi,
2372 DShape& dpsidx,
2373 Shape& test,
2374 DShape& dtestdx) const
2375 {
2376 //--------------------------
2377 // Call the shape functions:
2378 //--------------------------
2379 // Find the element dimension
2380 const unsigned el_dim = this->dim();
2381
2382 // Get the values of the shape functions and their local derivatives,
2383 // temporarily stored in dpsi
2384 this->dshape_local(s, psi, dpsidx);
2385
2386 // Allocate memory for the inverse jacobian
2387 DenseMatrix<double> inverse_jacobian(el_dim);
2388
2389 // Now calculate the inverse jacobian
2390 const double det =
2391 this->local_to_eulerian_mapping(dpsidx, inverse_jacobian);
2392
2393 // Now set the values of the derivatives to be dpsidx
2394 this->transform_derivatives(inverse_jacobian, dpsidx);
2395
2396 //-------------------------
2397 // Call the test functions:
2398 //-------------------------
2399 // Make sure we're using 3D elements
2400 if (el_dim != 3)
2401 {
2402 // Create an output stream
2403 std::ostringstream error_message_stream;
2404
2405 // Create an error message
2406 error_message_stream << "Need 3D space-time elements for this to work!"
2407 << std::endl;
2408
2409 // Throw the error message
2410 throw OomphLibError(error_message_stream.str(),
2411 OOMPH_CURRENT_FUNCTION,
2412 OOMPH_EXCEPTION_LOCATION);
2413 }
2414
2415 //--------start_of_dshape_local--------------------------------------
2416 // Local storage
2417 double test_values[3][3];
2418 double dtest_values[3][3];
2419
2420 // Index of the total shape function
2421 unsigned index = 0;
2422
2423 // Call the 1D shape functions and derivatives
2424 OneDimLagrange::shape<3>(s[0], test_values[0]);
2425 OneDimLagrange::shape<3>(s[1], test_values[1]);
2426 OneDimLagrange::dshape<3>(s[0], dtest_values[0]);
2427 OneDimLagrange::dshape<3>(s[1], dtest_values[1]);
2428
2429 // Set the time discretisation
2430 OneDimDiscontinuousGalerkin::shape<3>(s[2], test_values[2]);
2431 OneDimDiscontinuousGalerkin::dshape<3>(s[2], dtest_values[2]);
2432
2433 // Loop over the nodes in the third local coordinate direction
2434 for (unsigned k = 0; k < 3; k++)
2435 {
2436 // Loop over the nodes in the second local coordinate direction
2437 for (unsigned j = 0; j < 3; j++)
2438 {
2439 // Loop over the nodes in the first local coordinate direction
2440 for (unsigned i = 0; i < 3; i++)
2441 {
2442 // Calculate dtest/ds_0
2443 dtestdx(index, 0) =
2444 dtest_values[0][i] * test_values[1][j] * test_values[2][k];
2445
2446 // Calculate dtest/ds_1
2447 dtestdx(index, 1) =
2448 test_values[0][i] * dtest_values[1][j] * test_values[2][k];
2449
2450 // Calculate dtest/ds_2
2451 dtestdx(index, 2) =
2452 test_values[0][i] * test_values[1][j] * dtest_values[2][k];
2453
2454 // Calculate the index-th entry of test
2455 test[index] =
2456 test_values[0][i] * test_values[1][j] * test_values[2][k];
2457
2458 // Increment the index
2459 index++;
2460 }
2461 } // for (unsigned j=0;j<3;j++)
2462 } // for (unsigned k=0;k<3;k++)
2463 //--------end_of_dshape_local----------------------------------------
2464
2465 // Transform derivatives from dtest/ds to dtest/dx
2466 this->transform_derivatives(inverse_jacobian, dtestdx);
2467
2468 // Return the determinant value
2469 return det;
2470 } // End of dshape_and_dtest_eulerian_nst
2471
2472 //==========================================================================
2473 /// Derivatives of the shape functions and test functions w.r.t to
2474 /// global (Eulerian) coordinates. Return Jacobian of mapping between
2475 /// local and global coordinates.
2476 //==========================================================================
2477 template<unsigned DIM>
2478 inline double QTaylorHoodSpaceTimeElement<
2479 DIM>::dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
2480 Shape& psi,
2481 DShape& dpsidx,
2482 Shape& test,
2483 DShape& dtestdx) const
2484 {
2485 // Calculate the element dimension
2486 const unsigned el_dim = DIM + 1;
2487
2488 // Storage for the local coordinates of the integration point
2489 Vector<double> s(el_dim, 0.0);
2490
2491 // Set the local coordinate
2492 for (unsigned i = 0; i < el_dim; i++)
2493 {
2494 // Calculate the i-th local coordinate at the ipt-th knot point
2495 s[i] = this->integral_pt()->knot(ipt, i);
2496 }
2497
2498 // Return the Jacobian of the geometrical shape functions and derivatives
2499 return dshape_and_dtest_eulerian_nst(s, psi, dpsidx, test, dtestdx);
2500 } // End of dshape_and_dtest_eulerian_at_knot_nst
2501
2502
2503 //==========================================================================
2504 /// 2D (in space): Pressure shape and test functions and derivs w.r.t. to
2505 /// Eulerian coordinates. Return Jacobian of mapping between local and
2506 /// global coordinates.
2507 //==========================================================================
2508 template<>
2510 const Vector<double>& s, Shape& ppsi, DShape& dppsidx) const
2511 {
2512 // Local storage for the shape function (x-direction)
2513 double psi1[2];
2514
2515 // Local storage for the shape function (y-direction)
2516 double psi2[2];
2517
2518 // Local storage for the shape function (t-direction)
2519 double psi3[3];
2520
2521 // Local storage for the shape function derivatives (x-direction)
2522 double dpsi1[2];
2523
2524 // Local storage for the test function derivatives (y-direction)
2525 double dpsi2[2];
2526
2527 // Local storage for the test function derivatives (t-direction)
2528 double dpsi3[3];
2529
2530 // Call the OneDimensional Shape functions
2531 OneDimLagrange::shape<2>(s[0], psi1);
2532
2533 // Call the OneDimensional Shape functions
2534 OneDimLagrange::shape<2>(s[1], psi2);
2535
2536 // Call the OneDimensional Shape functions
2537 OneDimLagrange::shape<3>(s[2], psi3);
2538
2539 // Call the OneDimensional Shape functions
2540 OneDimLagrange::dshape<2>(s[0], dpsi1);
2541
2542 // Call the OneDimensional Shape functions
2543 OneDimLagrange::dshape<2>(s[1], dpsi2);
2544
2545 // Call the OneDimensional Shape functions
2546 OneDimLagrange::dshape<3>(s[2], dpsi3);
2547
2548 //--------------------------------------------------------------------
2549 // Now let's loop over the nodal points in the element with s1 being
2550 // the "x" coordinate, s2 being the "y" coordinate and s3 being the
2551 // "z" coordinate:
2552 //--------------------------------------------------------------------
2553 // Loop over the points in the t-direction
2554 for (unsigned i = 0; i < 3; i++)
2555 {
2556 // Loop over the points in the y-direction
2557 for (unsigned j = 0; j < 2; j++)
2558 {
2559 // Loop over the points in the x-direction
2560 for (unsigned k = 0; k < 2; k++)
2561 {
2562 // Multiply the three 1D functions together to get the 3D function
2563 ppsi[4 * i + 2 * j + k] = psi3[i] * psi2[j] * psi1[k];
2564
2565 // Multiply the appropriate shape and shape function derivatives
2566 // together
2567 dppsidx(4 * i + 2 * j + k, 0) = psi3[i] * psi2[j] * dpsi1[k];
2568
2569 // Multiply the appropriate shape and shape function derivatives
2570 // together
2571 dppsidx(4 * i + 2 * j + k, 1) = psi3[i] * dpsi2[j] * psi1[k];
2572
2573 // Multiply the appropriate shape and shape function derivatives
2574 // together
2575 dppsidx(4 * i + 2 * j + k, 2) = dpsi3[i] * psi2[j] * psi1[k];
2576 }
2577 } // for (unsigned j=0;j<2;j++)
2578 } // for (unsigned i=0;i<3;i++)
2579
2580 // Allocate space for the geometrical shape functions
2581 Shape psi(27);
2582
2583 // Allocate space for the geometrical shape function derivatives
2584 DShape dpsi(27, 3);
2585
2586 // Get the values of the shape functions and their derivatives
2587 dshape_local(s, psi, dpsi);
2588
2589 // Allocate memory for the 3x3 inverse jacobian
2590 DenseMatrix<double> inverse_jacobian(3);
2591
2592 // Now calculate the inverse jacobian
2593 const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
2594
2595 // Now set the values of the derivatives to be derivatives w.r.t. to
2596 // the Eulerian coordinates
2597 transform_derivatives(inverse_jacobian, dppsidx);
2598
2599 // Return the determinant of the jacobian
2600 return det;
2601 } // End of dpshape_eulerian
2602
2603
2604 //==========================================================================
2605 /// 2D (in space): Pressure shape and test functions and derivs w.r.t. to
2606 /// Eulerian coordinates. Return Jacobian of mapping between local and
2607 /// global coordinates.
2608 //==========================================================================
2609 template<>
2611 const Vector<double>& s, Shape& ptest, DShape& dptestdx) const
2612 {
2613 // Local storage for the shape function (x-direction)
2614 double test1[2];
2615
2616 // Local storage for the shape function (y-direction)
2617 double test2[2];
2618
2619 // Local storage for the shape function (t-direction)
2620 double test3[3];
2621
2622 // Local storage for the shape function derivatives (x-direction)
2623 double dtest1[2];
2624
2625 // Local storage for the test function derivatives (y-direction)
2626 double dtest2[2];
2627
2628 // Local storage for the test function derivatives (t-direction)
2629 double dtest3[3];
2630
2631 // Call the OneDimensional Shape functions
2632 OneDimLagrange::shape<2>(s[0], test1);
2633
2634 // Call the OneDimensional Shape functions
2635 OneDimLagrange::shape<2>(s[1], test2);
2636
2637 // Call the OneDimensional Shape functions
2639
2640 // Call the OneDimensional Shape functions
2641 OneDimLagrange::dshape<2>(s[0], dtest1);
2642
2643 // Call the OneDimensional Shape functions
2644 OneDimLagrange::dshape<2>(s[1], dtest2);
2645
2646 // Call the OneDimensional Shape functions
2648
2649 //--------------------------------------------------------------------
2650 // Now let's loop over the nodal points in the element with s1 being
2651 // the "x" coordinate, s2 being the "y" coordinate and s3 being the
2652 // "z" coordinate:
2653 //--------------------------------------------------------------------
2654 // Loop over the points in the t-direction
2655 for (unsigned i = 0; i < 3; i++)
2656 {
2657 // Loop over the points in the y-direction
2658 for (unsigned j = 0; j < 2; j++)
2659 {
2660 // Loop over the points in the x-direction
2661 for (unsigned k = 0; k < 2; k++)
2662 {
2663 // Multiply the three 1D functions together to get the 3D function
2664 ptest[4 * i + 2 * j + k] = test3[i] * test2[j] * test1[k];
2665
2666 // Multiply the appropriate shape and shape function derivatives
2667 // together
2668 dptestdx(4 * i + 2 * j + k, 0) = test3[i] * test2[j] * dtest1[k];
2669
2670 // Multiply the appropriate shape and shape function derivatives
2671 // together
2672 dptestdx(4 * i + 2 * j + k, 1) = test3[i] * dtest2[j] * test1[k];
2673
2674 // Multiply the appropriate shape and shape function derivatives
2675 // together
2676 dptestdx(4 * i + 2 * j + k, 2) = dtest3[i] * test2[j] * test1[k];
2677 }
2678 } // for (unsigned j=0;j<2;j++)
2679 } // for (unsigned i=0;i<3;i++)
2680
2681 // Allocate space for the geometrical shape functions
2682 Shape psi(27);
2683
2684 // Allocate space for the geometrical shape function derivatives
2685 DShape dpsi(27, 3);
2686
2687 // Get the values of the shape functions and their derivatives
2688 dshape_local(s, psi, dpsi);
2689
2690 // Allocate memory for the 3x3 inverse jacobian
2691 DenseMatrix<double> inverse_jacobian(3);
2692
2693 // Now calculate the inverse jacobian
2694 const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
2695
2696 // Now set the values of the derivatives to be derivatives w.r.t. to
2697 // the Eulerian coordinates
2698 transform_derivatives(inverse_jacobian, dptestdx);
2699
2700 // Return the determinant of the jacobian
2701 return det;
2702 } // End of dptest_eulerian
2703
2704
2705 //==========================================================================
2706 /// 2D (in space): Pressure shape and test functions and derivs w.r.t. to
2707 /// Eulerian coordinates. Return Jacobian of mapping between local and
2708 /// global coordinates.
2709 //==========================================================================
2710 template<>
2712 const Vector<double>& s,
2713 Shape& ppsi,
2714 DShape& dppsidx,
2715 Shape& ptest,
2716 DShape& dptestdx) const
2717 {
2718 // Call the test functions and derivatives
2719 dptest_eulerian(s, ptest, dptestdx);
2720
2721 // Call the shape functions and derivatives and the Jacobian of the mapping
2722 return this->dpshape_eulerian(s, ppsi, dppsidx);
2723 } // End of dpshape_and_dptest_eulerian_nst
2724
2725
2726 //==========================================================================
2727 /// 2D (in space):
2728 /// Define the shape functions (psi) and test functions (test) and
2729 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
2730 /// and return Jacobian of mapping (J). Additionally compute the
2731 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
2732 ///
2733 /// Galerkin: Test functions = shape functions
2734 //==========================================================================
2735 template<>
2738 const unsigned& ipt,
2739 Shape& psi,
2740 DShape& dpsidx,
2741 RankFourTensor<double>& d_dpsidx_dX,
2742 Shape& test,
2743 DShape& dtestdx,
2744 RankFourTensor<double>& d_dtestdx_dX,
2745 DenseMatrix<double>& djacobian_dX) const
2746 {
2747 // Call the geometrical shape functions and derivatives
2748 const double J = this->dshape_eulerian_at_knot(
2749 ipt, psi, dpsidx, djacobian_dX, d_dpsidx_dX);
2750
2751 // DRAIG: Delete!
2752 throw OomphLibError("Hasn't been implemented properly yet!",
2753 OOMPH_CURRENT_FUNCTION,
2754 OOMPH_EXCEPTION_LOCATION);
2755
2756 // Loop over the test functions and derivatives
2757 for (unsigned i = 0; i < 27; i++)
2758 {
2759 // The test functions are the same as the shape functions
2760 test[i] = psi[i];
2761
2762 // Loop over the spatial derivatives
2763 for (unsigned k = 0; k < 3; k++)
2764 {
2765 // Set the test function derivatives to the shape function derivatives
2766 dtestdx(i, k) = dpsidx(i, k);
2767
2768 // Loop over the dimensions
2769 for (unsigned p = 0; p < 3; p++)
2770 {
2771 // Loop over test functions
2772 for (unsigned q = 0; q < 27; q++)
2773 {
2774 // Set the test function derivatives to the shape function
2775 // derivatives
2776 d_dtestdx_dX(p, q, i, k) = d_dpsidx_dX(p, q, i, k);
2777 }
2778 } // for (unsigned p=0;p<3;p++)
2779 } // for (unsigned k=0;k<3;k++)
2780 } // for (unsigned i=0;i<27;i++)
2781
2782 // Return the jacobian
2783 return J;
2784 } // End of dshape_and_dtest_eulerian_at_knot_nst
2785
2786
2787 //==========================================================================
2788 /// 2D (in space): Pressure shape functions
2789 //==========================================================================
2790 template<>
2792 const Vector<double>& s, Shape& psi) const
2793 {
2794 // Local storage for the shape function value (in the x-direction)
2795 double psi1[2];
2796
2797 // Local storage for the shape function value (in the y-direction)
2798 double psi2[2];
2799
2800 // Local storage for the shape function value (in the z-direction)
2801 double psi3[3];
2802
2803 // Call the OneDimensional Shape functions
2804 OneDimLagrange::shape<2>(s[0], psi1);
2805
2806 // Call the OneDimensional Shape functions
2807 OneDimLagrange::shape<2>(s[1], psi2);
2808
2809 // Call the OneDimensional Shape functions
2810 OneDimLagrange::shape<3>(s[2], psi3);
2811
2812 //--------------------------------------------------------------------
2813 // Now let's loop over the nodal points in the element with s1 being
2814 // the "x" coordinate, s2 being the "y" coordinate and s3 being the
2815 // "z" coordinate:
2816 //--------------------------------------------------------------------
2817 // Loop over the points in the t-direction
2818 for (unsigned i = 0; i < 3; i++)
2819 {
2820 // Loop over the points in the y-direction
2821 for (unsigned j = 0; j < 2; j++)
2822 {
2823 // Loop over the points in the x-direction
2824 for (unsigned k = 0; k < 2; k++)
2825 {
2826 // Multiply the three 1D functions together to get the 3D function
2827 psi[4 * i + 2 * j + k] = psi3[i] * psi2[j] * psi1[k];
2828 }
2829 } // for (unsigned j=0;j<2;j++)
2830 } // for (unsigned i=0;i<3;i++)
2831 } // End of pshape_nst
2832
2833
2834 //==========================================================================
2835 /// 2D (in space): Pressure shape functions
2836 //==========================================================================
2837 template<>
2839 Shape& test) const
2840 {
2841 // Local storage for the shape function value (in the x-direction)
2842 double test1[2];
2843
2844 // Local storage for the shape function value (in the y-direction)
2845 double test2[2];
2846
2847 // Local storage for the shape function value (in the z-direction)
2848 double test3[3];
2849
2850 // Call the OneDimensional Shape functions
2851 OneDimLagrange::shape<2>(s[0], test1);
2852
2853 // Call the OneDimensional Shape functions
2854 OneDimLagrange::shape<2>(s[1], test2);
2855
2856 // Call the OneDimensional Shape functions
2858
2859 //--------------------------------------------------------------------
2860 // Now let's loop over the nodal points in the element with s1 being
2861 // the "x" coordinate, s2 being the "y" coordinate and s3 being the
2862 // "z" coordinate:
2863 //--------------------------------------------------------------------
2864 // Loop over the points in the z-direction
2865 for (unsigned i = 0; i < 3; i++)
2866 {
2867 // Loop over the points in the y-direction
2868 for (unsigned j = 0; j < 2; j++)
2869 {
2870 // Loop over the points in the x-direction
2871 for (unsigned k = 0; k < 2; k++)
2872 {
2873 // Multiply the three 1D functions together to get the 3D function
2874 test[4 * i + 2 * j + k] = test3[i] * test2[j] * test1[k];
2875 }
2876 } // for (unsigned j=0;j<2;j++)
2877 } // for (unsigned i=0;i<3;i++)
2878 } // End of ptest_nst
2879
2880
2881 //==========================================================================
2882 /// Pressure shape and test functions
2883 //==========================================================================
2884 template<unsigned DIM>
2886 const Vector<double>& s, Shape& psi, Shape& test) const
2887 {
2888 // Call the pressure shape functions
2889 pshape_nst(s, psi);
2890
2891 // Call the pressure test functions
2892 ptest_nst(s, test);
2893 } // End of pshape_nst
2894
2895
2896 //=======================================================================
2897 /// Face geometry of the 2D Taylor_Hood elements
2898 //=======================================================================
2899 template<>
2901 : public virtual QElement<2, 3>
2902 {
2903 public:
2904 /// Constructor; empty
2905 FaceGeometry() : QElement<2, 3>() {}
2906 };
2907
2908 //=======================================================================
2909 /// Face geometry of the FaceGeometry of the 2D Taylor Hoodelements
2910 //=======================================================================
2911 template<>
2913 : public virtual QElement<1, 3>
2914 {
2915 public:
2916 /// Constructor; empty
2917 FaceGeometry() : QElement<1, 3>() {}
2918 };
2919
2920 /// /////////////////////////////////////////////////////////////////
2921 /// /////////////////////////////////////////////////////////////////
2922 /// /////////////////////////////////////////////////////////////////
2923
2924 //==========================================================
2925 /// Taylor Hood upgraded to become projectable
2926 //==========================================================
2927 template<class TAYLOR_HOOD_ELEMENT>
2929 : public virtual ProjectableElement<TAYLOR_HOOD_ELEMENT>
2930 {
2931 public:
2932 /// Constructor [this was only required explicitly
2933 /// from gcc 4.5.2 onwards...]
2935
2936
2937 /// Specify the values associated with field fld.
2938 /// The information is returned in a vector of pairs which comprise
2939 /// the Data object and the value within it, that correspond to field fld.
2940 /// In the underlying Taylor Hood elements the fld-th velocities are stored
2941 /// at the fld-th value of the nodes; the pressures (the dim-th
2942 /// field) are the dim-th values at the vertex nodes etc.
2944 {
2945 // Create the vector
2947
2948 // If we're dealing with the velocity dofs
2949 if (fld < this->dim() - 1)
2950 {
2951 // How many nodes in the element?
2952 unsigned nnod = this->nnode();
2953
2954 // Loop over all nodes
2955 for (unsigned j = 0; j < nnod; j++)
2956 {
2957 // Add the data value associated with the velocity components
2958 data_values.push_back(std::make_pair(this->node_pt(j), fld));
2959 }
2960 }
2961 // If we're dealing with the pressure dof
2962 else
2963 {
2964 // How many pressure dofs are there?
2965 // DRAIG: Shouldn't there be more?
2966 unsigned Pconv_size = this->dim();
2967
2968 // Loop over all vertex nodes
2969 for (unsigned j = 0; j < Pconv_size; j++)
2970 {
2971 // Get the vertex index associated with the j-th pressure dof
2972 unsigned vertex_index = this->Pconv[j];
2973
2974 // Add the data value associated with the pressure components
2975 data_values.push_back(
2976 std::make_pair(this->node_pt(vertex_index), fld));
2977 }
2978 } // if (fld<this->dim())
2979
2980 // Return the vector
2981 return data_values;
2982 } // End of data_values_of_field
2983
2984
2985 /// Number of fields to be projected: dim+1, corresponding to
2986 /// velocity components and pressure
2988 {
2989 // There are dim velocity dofs and 1 pressure dof
2990 return this->dim();
2991 } // End of nfields_for_projection
2992
2993
2994 /// Number of history values to be stored for fld-th field. Whatever
2995 /// the timestepper has set up for the velocity components and none for
2996 /// the pressure field. NOTE: The count includes the current value!
2997 unsigned nhistory_values_for_projection(const unsigned& fld)
2998 {
2999 // If we're dealing with the pressure dof
3000 if (fld == this->dim())
3001 {
3002 // Pressure doesn't have history values
3003 return this->node_pt(0)->ntstorage();
3004 }
3005 // If we're dealing with the velocity dofs
3006 else
3007 {
3008 // The velocity dofs have ntstorage() history values
3009 return this->node_pt(0)->ntstorage();
3010 }
3011 } // End of nhistory_values_for_projection
3012
3013
3014 /// Number of positional history values
3015 /// (Note: count includes current value!)
3017 {
3018 // Return the number of positional history values
3019 return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
3020 } // End of nhistory_values_for_coordinate_projection
3021
3022
3023 /// Return the Jacobian of the mapping and the shape functions
3024 /// of field fld at local coordinate s
3025 double jacobian_and_shape_of_field(const unsigned& fld,
3026 const Vector<double>& s,
3027 Shape& psi)
3028 {
3029 // How many dimensions in this element?
3030 unsigned n_dim = this->dim();
3031
3032 // How many nodes in this element?
3033 unsigned n_node = this->nnode();
3034
3035 // If we're on the pressure dof
3036 if (fld == n_dim)
3037 {
3038 // Call the pressure interpolation function
3039 this->pshape_nst(s, psi);
3040
3041 // Allocate space for the pressure shape function
3042 Shape psif(n_node);
3043
3044 // Allocate space for the pressure test function
3045 Shape testf(n_node);
3046
3047 // Allocate space for the pressure shape function derivatives
3048 DShape dpsifdx(n_node, n_dim);
3049
3050 // Allocate space for the pressure test function derivatives
3051 DShape dtestfdx(n_node, n_dim);
3052
3053 // Calculate the Jacobian of the mapping
3054 double J = this->dshape_and_dtest_eulerian_nst(
3055 s, psif, dpsifdx, testf, dtestfdx);
3056
3057 // Return the Jacobian
3058 return J;
3059 }
3060 // If we're on the velocity components
3061 else
3062 {
3063 // Allocate space for the test functions
3064 Shape testf(n_node);
3065
3066 // Allocate space for the shape function derivatives
3067 DShape dpsifdx(n_node, n_dim);
3068
3069 // Allocate space for the test function derivatives
3070 DShape dtestfdx(n_node, n_dim);
3071
3072 // Calculate the Jacobian of the mapping
3073 double J =
3074 this->dshape_and_dtest_eulerian_nst(s, psi, dpsifdx, testf, dtestfdx);
3075
3076 // Return the Jacobian
3077 return J;
3078 }
3079 } // End of jacobian_and_shape_of_field
3080
3081
3082 /// Return interpolated field fld at local coordinate s, at time
3083 /// level t (t=0: present; t>0: history values)
3084 double get_field(const unsigned& t,
3085 const unsigned& fld,
3086 const Vector<double>& s)
3087 {
3088 unsigned n_dim = this->dim();
3089 unsigned n_node = this->nnode();
3090
3091 // If fld=n_dim, we deal with the pressure
3092 if (fld == n_dim)
3093 {
3094 return this->interpolated_p_nst(t, s);
3095 }
3096 // Velocity
3097 else
3098 {
3099 // Find the index at which the variable is stored
3100 unsigned u_nodal_index = this->u_index_nst(fld);
3101
3102 // Local shape function
3103 Shape psi(n_node);
3104
3105 // Find values of shape function
3106 this->shape(s, psi);
3107
3108 // Initialise value of u
3109 double interpolated_u = 0.0;
3110
3111 // Sum over the local nodes at that time
3112 for (unsigned l = 0; l < n_node; l++)
3113 {
3114 interpolated_u += this->nodal_value(t, l, u_nodal_index) * psi[l];
3115 }
3116 return interpolated_u;
3117 }
3118 }
3119
3120
3121 /// Return number of values in field fld
3122 unsigned nvalue_of_field(const unsigned& fld)
3123 {
3124 if (fld == this->dim())
3125 {
3126 return this->npres_nst();
3127 }
3128 else
3129 {
3130 return this->nnode();
3131 }
3132 }
3133
3134
3135 /// Return local equation number of value j in field fld.
3136 int local_equation(const unsigned& fld, const unsigned& j)
3137 {
3138 if (fld == this->dim())
3139 {
3140 return this->p_local_eqn(j);
3141 }
3142 else
3143 {
3144 const unsigned u_nodal_index = this->u_index_nst(fld);
3145 return this->nodal_local_eqn(j, u_nodal_index);
3146 }
3147 }
3148 };
3149
3150
3151 //=======================================================================
3152 /// Face geometry for element is the same as that for the underlying
3153 /// wrapped element
3154 //=======================================================================
3155 template<class ELEMENT>
3157 : public virtual FaceGeometry<ELEMENT>
3158 {
3159 public:
3160 FaceGeometry() : FaceGeometry<ELEMENT>() {}
3161 };
3162
3163 //=======================================================================
3164 /// Face geometry of the Face Geometry for element is the same as
3165 /// that for the underlying wrapped element
3166 //=======================================================================
3167 template<class ELEMENT>
3170 : public virtual FaceGeometry<FaceGeometry<ELEMENT>>
3171 {
3172 public:
3174 };
3175
3176} // End of namespace oomph
3177
3178#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.
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:5231
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....
Wrapper class for projectable elements. Adds "projectability" to the underlying ELEMENT.
Definition: projection.h:183
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
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 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.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
double 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...
unsigned ndof_types() const
Returns the number of "DOF types" that degrees of freedom in this element are sub-divided into: Veloc...
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 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 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.
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.
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 get_vorticity(const Vector< double > &s, double &vorticity) const
Compute the vorticity vector at local coordinate s.
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...
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.
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.
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.
double d_kin_energy_dt() const
Get integral of time derivative of kinetic energy 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.
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.
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...
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...
NavierStokesPressureAdvDiffSourceFctPt & source_fct_for_pressure_adv_diff()
Access function for the source-function pointer for pressure advection diffusion (used for validation...
void pin_all_non_pressure_dofs(std::map< Data *, std::vector< int > > &eqn_number_backup)
Pin all non-pressure dofs and backup eqn numbers.
virtual void pshape_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
virtual unsigned u_index_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
virtual double get_source_nst(const double &time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at a given time and Eulerian position.
void output_veloc(std::ostream &outfile, const unsigned &nplot, const unsigned &t)
Output function: x,y,t,u,v in tecplot format. Use n_plot points in each coordinate direction at times...
double interpolated_u_nst(const unsigned &t, const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s time level, t. Purposely broken for space-...
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Validate against exact solution at given time Solution is provided via function pointer....
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
void get_vorticity(const Vector< double > &s, Vector< double > &vorticity) const
Compute the vorticity vector at local coordinate s.
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 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).
//////////////////////////////////////////////////////////////////// ////////////////////////////////...