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