spherical_navier_stokes_elements.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for Navier Stokes elements
27
28#ifndef OOMPH_SPHERICAL_NAVIER_STOKES_ELEMENTS_HEADER
29#define OOMPH_SPHERICAL_NAVIER_STOKES_ELEMENTS_HEADER
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36
37// OOMPH-LIB headers
38#include "../generic/Qelements.h"
39#include "../generic/fsi.h"
40//#include "generic/block_preconditioner.h"
41
42namespace oomph
43{
44 //======================================================================
45 /// A class for elements that solve the Navier--Stokes equations,
46 /// in axisymmetric spherical polar coordinates.
47 /// This contains the generic maths -- any concrete implementation must
48 /// be derived from this.
49 ///
50 /// We also provide all functions required to use this element
51 /// in FSI problems, by deriving it from the FSIFluidElement base
52 /// class.
53 //======================================================================
55 : public virtual FSIFluidElement,
57 {
58 public:
59 /// Function pointer to body force function fct(t,x,f(x))
60 /// x is a Vector!
62 const double& time, const Vector<double>& x, Vector<double>& body_force);
63
64 /// Function pointer to source function fct(t,x)
65 /// x is a Vector!
66 typedef double (*SphericalNavierStokesSourceFctPt)(const double& time,
67 const Vector<double>& x);
68
69 private:
70 /// Static "magic" number that indicates that the pressure is
71 /// not stored at a node
73
74 /// Static default value for the physical constants (all initialised to
75 /// zero)
77
78 /// Static default value for the physical ratios (all are initialised to
79 /// one)
81
82 /// Static default value for the gravity vector
84
85 protected:
86 // Physical constants
87
88 /// Pointer to the viscosity ratio (relative to the
89 /// viscosity used in the definition of the Reynolds number)
91
92 /// Pointer to the density ratio (relative to the density used in the
93 /// definition of the Reynolds number)
95
96 // Pointers to global physical constants
97
98 /// Pointer to global Reynolds number
99 double* Re_pt;
100
101 /// Pointer to global Reynolds number x Strouhal number (=Womersley)
102 double* ReSt_pt;
103
104 /// Pointer to global Reynolds number x inverse Froude number
105 /// (= Bond number / Capillary number)
106 double* ReInvFr_pt;
107
108 /// Pointer to global Reynolds number x inverse Rossby number
109 /// (used when in a rotating frame)
110 double* ReInvRo_pt;
111
112 /// Pointer to global gravity Vector
114
115 /// Pointer to body force function
117
118 /// Pointer to volumetric source function
120
121 /// Boolean flag to indicate if ALE formulation is disabled when
122 /// time-derivatives are computed. Only set to true if you're sure
123 /// that the mesh is stationary.
125
126 /// Access function for the local equation number information for
127 /// the pressure.
128 /// p_local_eqn[n] = local equation number or < 0 if pinned
129 virtual int p_local_eqn(const unsigned& n) const = 0;
130
131 /// Compute the shape functions and derivatives
132 /// w.r.t. global coords at local coordinate s.
133 /// Return Jacobian of mapping between local and global coordinates.
135 const Vector<double>& s,
136 Shape& psi,
137 DShape& dpsidx,
138 Shape& test,
139 DShape& dtestdx) const = 0;
140
141 /// Compute the shape functions and derivatives
142 /// w.r.t. global coords at ipt-th integration point
143 /// Return Jacobian of mapping between local and global coordinates.
145 const unsigned& ipt,
146 Shape& psi,
147 DShape& dpsidx,
148 Shape& test,
149 DShape& dtestdx) const = 0;
150
151 /// Compute the pressure shape functions at local coordinate s
153 Shape& psi) const = 0;
154
155 /// Compute the pressure shape and test functions
156 /// at local coordinate s
158 Shape& psi,
159 Shape& test) const = 0;
160
161
162 /// Calculate the body force at a given time and local and/or
163 /// Eulerian position. This function is virtual so that it can be
164 /// overloaded in multi-physics elements where the body force might
165 /// depend on another variable.
166 virtual void get_body_force_spherical_nst(const double& time,
167 const unsigned& ipt,
168 const Vector<double>& s,
169 const Vector<double>& x,
170 Vector<double>& result)
171 {
172 // If the function pointer is zero return zero
173 if (Body_force_fct_pt == 0)
174 {
175 // Loop over three spatial dimensions and set body forces to zero
176 for (unsigned i = 0; i < 3; i++)
177 {
178 result[i] = 0.0;
179 }
180 }
181 // Otherwise call the function
182 else
183 {
184 (*Body_force_fct_pt)(time, x, result);
185 }
186 }
187
188 /// Calculate the source fct at given time and
189 /// Eulerian position
190 virtual double get_source_spherical_nst(double time,
191 const unsigned& ipt,
192 const Vector<double>& x)
193 {
194 // If the function pointer is zero return zero
195 if (Source_fct_pt == 0)
196 {
197 return 0;
198 }
199 // Otherwise call the function
200 else
201 {
202 return (*Source_fct_pt)(time, x);
203 }
204 }
205
206 /// Compute the residuals for the Navier--Stokes equations;
207 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
209 Vector<double>& residuals,
210 DenseMatrix<double>& jacobian,
211 DenseMatrix<double>& mass_matrix,
212 unsigned flag);
213
214
215 public:
216 /// Include a cot function to simplify the
217 /// NS momentum and jacobian expressions
218 inline double cot(const double& th) const
219 {
220 return (1 / tan(th));
221 }
222
223
224 /// Constructor: NULL the body force and source function
225 /// and make sure the ALE terms are included by default.
228 {
229 // Set all the Physical parameter pointers to the default value zero
235 // Set the Physical ratios to the default value of 1
238 }
239
240 /// Vector to decide whether the stress-divergence form is used or not
241 // N.B. This needs to be public so that the intel compiler gets things
242 // correct somehow the access function messes things up when going to
243 // refineable navier--stokes
245
246 // Access functions for the physical constants
247
248 /// Reynolds number
249 const double& re() const
250 {
251 return *Re_pt;
252 }
253
254 /// Product of Reynolds and Strouhal number (=Womersley number)
255 const double& re_st() const
256 {
257 return *ReSt_pt;
258 }
259
260 /// Pointer to Reynolds number
261 double*& re_pt()
262 {
263 return Re_pt;
264 }
265
266 /// Pointer to product of Reynolds and Strouhal number (=Womersley number)
267 double*& re_st_pt()
268 {
269 return ReSt_pt;
270 }
271
272 /// Viscosity ratio for element: Element's viscosity relative to the
273 /// viscosity used in the definition of the Reynolds number
274 const double& viscosity_ratio() const
275 {
276 return *Viscosity_Ratio_pt;
277 }
278
279 /// Pointer to Viscosity Ratio
281 {
282 return Viscosity_Ratio_pt;
283 }
284
285 /// Density ratio for element: Element's density relative to the
286 /// viscosity used in the definition of the Reynolds number
287 const double& density_ratio() const
288 {
289 return *Density_Ratio_pt;
290 }
291
292 /// Pointer to Density ratio
294 {
295 return Density_Ratio_pt;
296 }
297
298 /// Global inverse Froude number
299 const double& re_invfr() const
300 {
301 return *ReInvFr_pt;
302 }
303
304 /// Pointer to global inverse Froude number
305 double*& re_invfr_pt()
306 {
307 return ReInvFr_pt;
308 }
309
310 /// Global Reynolds number multiplied by inverse Rossby number
311 const double& re_invro() const
312 {
313 return *ReInvRo_pt;
314 }
315
316 /// Pointer to global inverse Froude number
317 double*& re_invro_pt()
318 {
319 return ReInvRo_pt;
320 }
321
322 /// Vector of gravitational components
323 const Vector<double>& g() const
324 {
325 return *G_pt;
326 }
327
328 /// Pointer to Vector of gravitational components
330 {
331 return G_pt;
332 }
333
334 /// Access function for the body-force pointer
336 {
337 return Body_force_fct_pt;
338 }
339
340 /// Access function for the body-force pointer. Const version
342 {
343 return Body_force_fct_pt;
344 }
345
346 /// Access function for the source-function pointer
348 {
349 return Source_fct_pt;
350 }
351
352 /// Access function for the source-function pointer. Const version
354 {
355 return Source_fct_pt;
356 }
357
358 /// Function to return number of pressure degrees of freedom
359 virtual unsigned npres_spherical_nst() const = 0;
360
361 /// Velocity i at local node n. Uses suitably interpolated value
362 /// for hanging nodes.
363 /// The use of u_index_spherical_nst() permits the use of this
364 /// element as the basis for multi-physics elements. The default
365 /// is to assume that the i-th velocity component is stored at the
366 /// i-th location of the node
367 double u_spherical_nst(const unsigned& n, const unsigned& i) const
368 {
370 }
371
372 /// Velocity i at local node n at timestep t (t=0: present;
373 /// t>0: previous). Uses suitably interpolated value for hanging nodes.
374 double u_spherical_nst(const unsigned& t,
375 const unsigned& n,
376 const unsigned& i) const
377 {
379 }
380
381 /// Return the index at which the i-th unknown velocity component
382 /// is stored. The default value, i, is appropriate for single-physics
383 /// problems.
384 /// In derived multi-physics elements, this function should be overloaded
385 /// to reflect the chosen storage scheme. Note that these equations require
386 /// that the unknowns are always stored at the same indices at each node.
387 virtual inline unsigned u_index_spherical_nst(const unsigned& i) const
388 {
389 return i;
390 }
391
392
393 /// i-th component of du/dt at local node n.
394 /// Uses suitably interpolated value for hanging nodes.
395 double du_dt_spherical_nst(const unsigned& n, const unsigned& i) const
396 {
397 // Get the data's timestepper
399
400 // Initialise dudt
401 double dudt = 0.0;
402
403 // Loop over the timesteps, if there is a non Steady timestepper
404 if (time_stepper_pt->type() != "Steady")
405 {
406 // Find the index at which the dof is stored
407 const unsigned u_nodal_index = this->u_index_spherical_nst(i);
408
409 // Number of timsteps (past & present)
410 const unsigned n_time = time_stepper_pt->ntstorage();
411 // Loop over the timesteps
412 for (unsigned t = 0; t < n_time; t++)
413 {
414 dudt +=
415 time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
416 }
417 }
418
419 return dudt;
420 }
421
422 /// Disable ALE, i.e. assert the mesh is not moving -- you do this
423 /// at your own risk!
425 {
426 ALE_is_disabled = true;
427 }
428
429 /// (Re-)enable ALE, i.e. take possible mesh motion into account
430 /// when evaluating the time-derivative. Note: By default, ALE is
431 /// enabled, at the expense of possibly creating unnecessary work
432 /// in problems where the mesh is, in fact, stationary.
434 {
435 ALE_is_disabled = false;
436 }
437
438 /// Pressure at local pressure "node" n_p
439 /// Uses suitably interpolated value for hanging nodes.
440 virtual double p_spherical_nst(const unsigned& n_p) const = 0;
441
442 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
443 virtual void fix_pressure(const unsigned& p_dof, const double& p_value) = 0;
444
445 /// Return the index at which the pressure is stored if it is
446 /// stored at the nodes. If not stored at the nodes this will return
447 /// a negative number.
448 virtual int p_nodal_index_spherical_nst() const
449 {
451 }
452
453 /// Integral of pressure over element
454 double pressure_integral() const;
455
456 /// Return integral of dissipation over element
457 double dissipation() const;
458
459 /// Return dissipation at local coordinate s
460 double dissipation(const Vector<double>& s) const;
461
462 /// Compute the vorticity vector at local coordinate s
463 void get_vorticity(const Vector<double>& s,
464 Vector<double>& vorticity) const;
465
466 /// Get integral of kinetic energy over element
467 double kin_energy() const;
468
469 /// Get integral of time derivative of kinetic energy over element
470 double d_kin_energy_dt() const;
471
472 /// Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)
473 void strain_rate(const Vector<double>& s,
475
476 /// Compute traction (on the viscous scale) exerted onto
477 /// the fluid at local coordinate s. N has to be outer unit normal
478 /// to the fluid.
479 void get_traction(const Vector<double>& s,
480 const Vector<double>& N,
481 Vector<double>& traction);
482
483
484 /// This implements a pure virtual function defined
485 /// in the FSIFluidElement class. The function computes
486 /// the traction (on the viscous scale), at the element's local
487 /// coordinate s, that the fluid element exerts onto an adjacent
488 /// solid element. The number of arguments is imposed by
489 /// the interface defined in the FSIFluidElement -- only the
490 /// unit normal N (pointing into the fluid!) is actually used
491 /// in the computation.
493 const Vector<double>& N,
494 Vector<double>& load)
495 {
496 // Note: get_traction() computes the traction onto the fluid
497 // if N is the outer unit normal onto the fluid; here we're
498 // exepcting N to point into the fluid so we're getting the
499 // traction onto the adjacent wall instead!
500 get_traction(s, N, load);
501 }
502
503 /// Compute the diagonal of the velocity/pressure mass matrices.
504 /// If which one=0, both are computed, otherwise only the pressure
505 /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
506 /// LSC version of the preconditioner only needs that one)
507 /// NOTE: pressure versions isn't implemented yet because this
508 /// element has never been tried with Fp preconditoner.
510 Vector<double>& press_mass_diag,
511 Vector<double>& veloc_mass_diag,
512 const unsigned& which_one = 0);
513
514
515 /// Output function: x,y,[z],u,v,[w],p
516 /// in tecplot format. Default number of plot points
517 void output(std::ostream& outfile)
518 {
519 unsigned nplot = 5;
520 output(outfile, nplot);
521 }
522
523 /// Output function: x,y,[z],u,v,[w],p
524 /// in tecplot format. nplot points in each coordinate direction
525 void output(std::ostream& outfile, const unsigned& nplot);
526
527 /// C-style output function: x,y,[z],u,v,[w],p
528 /// in tecplot format. Default number of plot points
529 void output(FILE* file_pt)
530 {
531 unsigned nplot = 5;
532 output(file_pt, nplot);
533 }
534
535 /// C-style output function: x,y,[z],u,v,[w],p
536 /// in tecplot format. nplot points in each coordinate direction
537 void output(FILE* file_pt, const unsigned& nplot);
538
539 /// Full output function:
540 /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
541 /// in tecplot format. Default number of plot points
542 void full_output(std::ostream& outfile)
543 {
544 unsigned nplot = 5;
545 full_output(outfile, nplot);
546 }
547
548 /// Full output function:
549 /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
550 /// in tecplot format. nplot points in each coordinate direction
551 void full_output(std::ostream& outfile, const unsigned& nplot);
552
553
554 /// Output function: x,y,[z],u,v,[w] in tecplot format.
555 /// nplot points in each coordinate direction at timestep t
556 /// (t=0: present; t>0: previous timestep)
557 void output_veloc(std::ostream& outfile,
558 const unsigned& nplot,
559 const unsigned& t);
560
561
562 /// Output function: x,y,[z], [omega_x,omega_y,[and/or omega_z]]
563 /// in tecplot format. nplot points in each coordinate direction
564 void output_vorticity(std::ostream& outfile, const unsigned& nplot);
565
566 /// Output exact solution specified via function pointer
567 /// at a given number of plot points. Function prints as
568 /// many components as are returned in solution Vector
569 void output_fct(std::ostream& outfile,
570 const unsigned& nplot,
572
573 /// Output exact solution specified via function pointer
574 /// at a given time and at a given number of plot points.
575 /// Function prints as many components as are returned in solution Vector.
576 void output_fct(std::ostream& outfile,
577 const unsigned& nplot,
578 const double& time,
580
581 /// Validate against exact solution at given time
582 /// Solution is provided via function pointer.
583 /// Plot at a given number of plot points and compute L2 error
584 /// and L2 norm of velocity solution over element
585 void compute_error(std::ostream& outfile,
587 const double& time,
588 double& error,
589 double& norm);
590
591 /// Validate against exact solution.
592 /// Solution is provided via function pointer.
593 /// Plot at a given number of plot points and compute L2 error
594 /// and L2 norm of velocity solution over element
595 void compute_error(std::ostream& outfile,
597 double& error,
598 double& norm);
599
600 /// Validate against exact solution.
601 /// Solution is provided direct from exact_soln function.
602 /// Plot at a given number of plot points and compute the energy error
603 /// and energy norm of the velocity solution over the element.
604 void compute_error_e(
605 std::ostream& outfile,
608 FiniteElement::SteadyExactSolutionFctPt exact_soln_dtheta_pt,
609 double& u_error,
610 double& u_norm,
611 double& p_error,
612 double& p_norm);
613
614 // Listing for the shear stress function
615 void compute_shear_stress(std::ostream& outfile);
616
617 // Listing for the velocity extraction function
618 void extract_velocity(std::ostream& outfile);
619
620 // Calculate the analytic solution at the point x and output a vector
622 {
623 const double Re = this->re();
624
625 double r = x[0];
626 double theta = x[1];
627 Vector<double> ans(4, 0.0);
628
629 ans[2] = r * sin(theta);
630 ans[3] = 0.5 * Re * r * r * sin(theta) * sin(theta);
631
632 return (ans);
633 }
634
635 // Calculate the r-derivatives of the analytic solution at the point x and
636 // output a vector
638 {
639 const double Re = this->re();
640
641 double r = x[0];
642 double theta = x[1];
643 Vector<double> ans(4, 0.0);
644
645 ans[2] = sin(theta);
646 ans[3] = Re * r * sin(theta) * sin(theta);
647
648 return (ans);
649 }
650
651 // Calculate the theta-derivatives of the analytic solution at the point x
652 // and output a vector
654 {
655 const double Re = this->re();
656
657 double r = x[0];
658 double theta = x[1];
659 Vector<double> ans(4, 0.0);
660
661 ans[2] = r * cos(theta);
662 ans[3] = Re * r * r * sin(theta) * cos(theta);
663
664 return (ans);
665 }
666
667 /// Compute the element's residual Vector
669 {
670 // Call the generic residuals function with flag set to 0
671 // and using a dummy matrix argument
673 residuals,
676 0);
677 }
678
679 // Compute the element's residual Vector and the jacobian matrix
680 // Virtual function can be overloaded by hanging-node version
682 DenseMatrix<double>& jacobian)
683 {
684 // Call the generic routine with the flag set to 1
686 residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
687 }
688
689 /// Add the element's contribution to its residuals vector,
690 /// jacobian matrix and mass matrix
692 Vector<double>& residuals,
693 DenseMatrix<double>& jacobian,
694 DenseMatrix<double>& mass_matrix)
695 {
696 // Call the generic routine with the flag set to 2
698 residuals, jacobian, mass_matrix, 2);
699 }
700
701 /// Compute vector of FE interpolated velocity u at local coordinate s
703 Vector<double>& veloc) const
704 {
705 // Find number of nodes
706 unsigned n_node = nnode();
707 // Local shape function
708 Shape psi(n_node);
709 // Find values of shape function
710 shape(s, psi);
711
712 for (unsigned i = 0; i < 3; i++)
713 {
714 // Index at which the nodal value is stored
715 unsigned u_nodal_index = u_index_spherical_nst(i);
716 // Initialise value of u
717 veloc[i] = 0.0;
718 // Loop over the local nodes and sum
719 for (unsigned l = 0; l < n_node; l++)
720 {
721 veloc[i] += nodal_value(l, u_nodal_index) * psi[l];
722 }
723 }
724 }
725
726 /// Return FE interpolated velocity u[i] at local coordinate s
728 const unsigned& i) const
729 {
730 // Find number of nodes
731 unsigned n_node = nnode();
732 // Local shape function
733 Shape psi(n_node);
734 // Find values of shape function
735 shape(s, psi);
736
737 // Get nodal index at which i-th velocity is stored
738 unsigned u_nodal_index = u_index_spherical_nst(i);
739
740 // Initialise value of u
741 double interpolated_u = 0.0;
742 // Loop over the local nodes and sum
743 for (unsigned l = 0; l < n_node; l++)
744 {
745 interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
746 }
747
748 return (interpolated_u);
749 }
750
751 /// Return matrix entry dudx(i,j) of the FE interpolated velocity derivative
752 /// at local coordinate s
754 const unsigned& i,
755 const unsigned& j) const
756 {
757 // Find number of nodes
758 unsigned n_node = nnode();
759 // Local shape function
760 Shape psi(n_node);
761 DShape dpsidx(n_node, 2);
762 // Find values of shape function
763 (void)dshape_eulerian(s, psi, dpsidx);
764
765 // Get nodal index at which i-th velocity is stored
766 unsigned u_nodal_index = u_index_spherical_nst(i);
767
768 // Initialise value of u
769 double interpolated_dudx = 0.0;
770 // Loop over the local nodes and sum
771 for (unsigned l = 0; l < n_node; l++)
772 {
773 interpolated_dudx += nodal_value(l, u_nodal_index) * dpsidx(l, j);
774 }
775
776 return (interpolated_dudx);
777 }
778
779 /// Return FE interpolated pressure at local coordinate s
781 {
782 // Find number of nodes
783 unsigned n_pres = npres_spherical_nst();
784 // Local shape function
785 Shape psi(n_pres);
786 // Find values of shape function
788
789 // Initialise value of p
790 double interpolated_p = 0.0;
791 // Loop over the local nodes and sum
792 for (unsigned l = 0; l < n_pres; l++)
793 {
794 interpolated_p += p_spherical_nst(l) * psi[l];
795 }
796
797 return (interpolated_p);
798 }
799 };
800
801 /// ///////////////////////////////////////////////////////////////////////////
802 /// ///////////////////////////////////////////////////////////////////////////
803 /// ///////////////////////////////////////////////////////////////////////////
804
805
806 //==========================================================================
807 /// Crouzeix_Raviart elements are Navier--Stokes elements with quadratic
808 /// interpolation for velocities and positions, but a discontinuous linear
809 /// pressure interpolation. They can be used within oomph-lib's
810 /// block preconditioning framework.
811 //==========================================================================
813 : public virtual QElement<2, 3>,
814 public virtual SphericalNavierStokesEquations
815 {
816 private:
817 /// Static array of ints to hold required number of variables at nodes
818 static const unsigned Initial_Nvalue[];
819
820 protected:
821 /// Internal index that indicates at which internal data the pressure
822 /// is stored
824
825
826 /// Velocity shape and test functions and their derivs
827 /// w.r.t. to global coords at local coordinate s (taken from geometry)
828 /// Return Jacobian of mapping between local and global coordinates.
830 const Vector<double>& s,
831 Shape& psi,
832 DShape& dpsidx,
833 Shape& test,
834 DShape& dtestdx) const;
835
836 /// Velocity shape and test functions and their derivs
837 /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
838 /// Return Jacobian of mapping between local and global coordinates.
840 const unsigned& ipt,
841 Shape& psi,
842 DShape& dpsidx,
843 Shape& test,
844 DShape& dtestdx) const;
845
846 /// Pressure shape functions at local coordinate s
847 inline void pshape_spherical_nst(const Vector<double>& s, Shape& psi) const;
848
849 /// Pressure shape and test functions at local coordinte s
850 inline void pshape_spherical_nst(const Vector<double>& s,
851 Shape& psi,
852 Shape& test) const;
853
854 /// Return the local equation numbers for the pressure values.
855 inline int p_local_eqn(const unsigned& n) const
856 {
857 return this->internal_local_eqn(P_spherical_nst_internal_index, n);
858 }
859
860 public:
861 /// Constructor, there are 3 internal values (for the pressure)
864 {
865 // Allocate and add one Internal data object that stored 2+1 pressure
866 // values;
868 }
869
870 /// Number of values (pinned or dofs) required at local node n.
871 inline unsigned required_nvalue(const unsigned& n) const
872 {
873 return 3;
874 }
875
876
877 /// Return the i-th pressure value
878 /// (Discontinous pressure interpolation -- no need to cater for hanging
879 /// nodes).
880 double p_spherical_nst(const unsigned& i) const
881 {
882 return this->internal_data_pt(P_spherical_nst_internal_index)->value(i);
883 }
884
885 /// Return number of pressure values
886 unsigned npres_spherical_nst() const
887 {
888 return 3;
889 }
890
891 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
892 void fix_pressure(const unsigned& p_dof, const double& p_value)
893 {
894 this->internal_data_pt(P_spherical_nst_internal_index)->pin(p_dof);
895 this->internal_data_pt(P_spherical_nst_internal_index)
896 ->set_value(p_dof, p_value);
897 }
898
899 /// Add to the set \c paired_load_data pairs containing
900 /// - the pointer to a Data object
901 /// and
902 /// - the index of the value in that Data object
903 /// .
904 /// for all values (pressures, velocities) that affect the
905 /// load computed in the \c get_load(...) function.
907 std::set<std::pair<Data*, unsigned>>& paired_load_data);
908
909 /// Add to the set \c paired_pressure_data pairs containing
910 /// - the pointer to a Data object
911 /// and
912 /// - the index of the value in that Data object
913 /// .
914 /// for pressure values that affect the
915 /// load computed in the \c get_load(...) function.
917 std::set<std::pair<Data*, unsigned>>& paired_pressure_data);
918
919 /// Redirect output to SphericalNavierStokesEquations output
920 void output(std::ostream& outfile)
921 {
923 }
924
925 /// Redirect output to SphericalNavierStokesEquations output
926 void output(std::ostream& outfile, const unsigned& nplot)
927 {
929 }
930
931
932 /// Redirect output to SphericalNavierStokesEquations output
933 void output(FILE* file_pt)
934 {
936 }
937
938 /// Redirect output to SphericalNavierStokesEquations output
939 void output(FILE* file_pt, const unsigned& nplot)
940 {
942 }
943
944
945 /// Full output function:
946 /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
947 /// in tecplot format. Default number of plot points
948 void full_output(std::ostream& outfile)
949 {
951 }
952
953 /// Full output function:
954 /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
955 /// in tecplot format. nplot points in each coordinate direction
956 void full_output(std::ostream& outfile, const unsigned& nplot)
957 {
959 }
960
961
962 /// The number of "DOF types" that degrees of freedom in this element
963 /// are sub-divided into: Velocity (three comp) and pressure.
964 unsigned ndof_types() const
965 {
966 return 4;
967 }
968
969 /// Create a list of pairs for all unknowns in this element,
970 /// so that the first entry in each pair contains the global equation
971 /// number of the unknown, while the second one contains the number
972 /// of the "DOF type" that this unknown is associated with.
973 /// (Function can obviously only be called if the equation numbering
974 /// scheme has been set up.)
976 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const;
977 };
978
979 // Inline functions
980
981 //=======================================================================
982 /// Derivatives of the shape functions and test functions w.r.t. to global
983 /// (Eulerian) coordinates. Return Jacobian of mapping between
984 /// local and global coordinates.
985 //=======================================================================
988 Shape& psi,
989 DShape& dpsidx,
990 Shape& test,
991 DShape& dtestdx) const
992 {
993 // Call the geometrical shape functions and derivatives
994 double J = this->dshape_eulerian(s, psi, dpsidx);
995 // Loop over the test functions and derivatives and set them equal to the
996 // shape functions
997 for (unsigned i = 0; i < 9; i++)
998 {
999 test[i] = psi[i];
1000 dtestdx(i, 0) = dpsidx(i, 0);
1001 dtestdx(i, 1) = dpsidx(i, 1);
1002 }
1003 // Return the jacobian
1004 return J;
1005 }
1006
1007 //=======================================================================
1008 /// Derivatives of the shape functions and test functions w.r.t. to global
1009 /// (Eulerian) coordinates. Return Jacobian of mapping between
1010 /// local and global coordinates.
1011 //=======================================================================
1014 Shape& psi,
1015 DShape& dpsidx,
1016 Shape& test,
1017 DShape& dtestdx) const
1018 {
1019 // Call the geometrical shape functions and derivatives
1020 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
1021 // Loop over the test functions and derivatives and set them equal to the
1022 // shape functions
1023 for (unsigned i = 0; i < 9; i++)
1024 {
1025 test[i] = psi[i];
1026 dtestdx(i, 0) = dpsidx(i, 0);
1027 dtestdx(i, 1) = dpsidx(i, 1);
1028 }
1029 // Return the jacobian
1030 return J;
1031 }
1032
1033
1034 //=======================================================================
1035 /// Pressure shape functions
1036 //=======================================================================
1038 const Vector<double>& s, Shape& psi) const
1039 {
1040 psi[0] = 1.0;
1041 psi[1] = s[0];
1042 psi[2] = s[1];
1043 }
1044
1045 /// Define the pressure shape and test functions
1047 const Vector<double>& s, Shape& psi, Shape& test) const
1048 {
1049 // Call the pressure shape functions
1050 pshape_spherical_nst(s, psi);
1051 // Loop over the test functions and set them equal to the shape functions
1052 for (unsigned i = 0; i < 3; i++) test[i] = psi[i];
1053 }
1054
1055 //=======================================================================
1056 /// Face geometry of the Spherical Crouzeix_Raviart elements
1057 //=======================================================================
1058 template<>
1060 : public virtual QElement<1, 3>
1061 {
1062 public:
1063 FaceGeometry() : QElement<1, 3>() {}
1064 };
1065
1066 //=======================================================================
1067 /// Face geometry of the FaceGeometry of the Spherical
1068 /// Crouzeix_Raviart elements
1069 //=======================================================================
1070 template<>
1072 : public virtual PointElement
1073 {
1074 public:
1076 };
1077
1078
1079 /// /////////////////////////////////////////////////////////////////////////
1080 /// /////////////////////////////////////////////////////////////////////////
1081
1082
1083 //=======================================================================
1084 /// Taylor--Hood elements are Navier--Stokes elements
1085 /// with quadratic interpolation for velocities and positions and
1086 /// continous linear pressure interpolation. They can be used
1087 /// within oomph-lib's block-preconditioning framework.
1088 //=======================================================================
1090 : public virtual QElement<2, 3>,
1091 public virtual SphericalNavierStokesEquations
1092 {
1093 private:
1094 /// Static array of ints to hold number of variables at node
1095 static const unsigned Initial_Nvalue[];
1096
1097 protected:
1098 /// Static array of ints to hold conversion from pressure
1099 /// node numbers to actual node numbers
1100 static const unsigned Pconv[];
1101
1102 /// Velocity shape and test functions and their derivs
1103 /// w.r.t. to global coords at local coordinate s (taken from geometry)
1104 /// Return Jacobian of mapping between local and global coordinates.
1106 const Vector<double>& s,
1107 Shape& psi,
1108 DShape& dpsidx,
1109 Shape& test,
1110 DShape& dtestdx) const;
1111
1112 /// Velocity shape and test functions and their derivs
1113 /// w.r.t. to global coords at local coordinate s (taken from geometry)
1114 /// Return Jacobian of mapping between local and global coordinates.
1116 const unsigned& ipt,
1117 Shape& psi,
1118 DShape& dpsidx,
1119 Shape& test,
1120 DShape& dtestdx) const;
1121
1122 /// Pressure shape functions at local coordinate s
1123 inline void pshape_spherical_nst(const Vector<double>& s, Shape& psi) const;
1124
1125 /// Pressure shape and test functions at local coordinte s
1126 inline void pshape_spherical_nst(const Vector<double>& s,
1127 Shape& psi,
1128 Shape& test) const;
1129
1130 /// Return the local equation numbers for the pressure values.
1131 inline int p_local_eqn(const unsigned& n) const
1132 {
1133 return this->nodal_local_eqn(Pconv[n], p_nodal_index_spherical_nst());
1134 }
1135
1136 public:
1137 /// Constructor, no internal data points
1140 {
1141 }
1142
1143 /// Number of values (pinned or dofs) required at node n. Can
1144 /// be overwritten for hanging node version
1145 inline virtual unsigned required_nvalue(const unsigned& n) const
1146 {
1147 return Initial_Nvalue[n];
1148 }
1149
1150 /// Set the value at which the pressure is stored in the nodes
1151 /// In this case the third index because there are three velocity components
1153 {
1154 return 3;
1155 }
1156
1157 /// Access function for the pressure values at local pressure
1158 /// node n_p (const version)
1159 double p_spherical_nst(const unsigned& n_p) const
1160 {
1161 return this->nodal_value(Pconv[n_p], this->p_nodal_index_spherical_nst());
1162 }
1163
1164 /// Return number of pressure values
1165 unsigned npres_spherical_nst() const
1166 {
1167 return 4;
1168 }
1169
1170 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
1171 void fix_pressure(const unsigned& p_dof, const double& p_value)
1172 {
1173 this->node_pt(Pconv[p_dof])->pin(this->p_nodal_index_spherical_nst());
1174 this->node_pt(Pconv[p_dof])
1175 ->set_value(this->p_nodal_index_spherical_nst(), p_value);
1176 }
1177
1178
1179 /// Add to the set \c paired_load_data pairs containing
1180 /// - the pointer to a Data object
1181 /// and
1182 /// - the index of the value in that Data object
1183 /// .
1184 /// for all values (pressures, velocities) that affect the
1185 /// load computed in the \c get_load(...) function.
1186 void identify_load_data(
1187 std::set<std::pair<Data*, unsigned>>& paired_load_data);
1188
1189
1190 /// Add to the set \c paired_pressure_data pairs containing
1191 /// - the pointer to a Data object
1192 /// and
1193 /// - the index of the value in that Data object
1194 /// .
1195 /// for pressure values that affect the
1196 /// load computed in the \c get_load(...) function.
1198 std::set<std::pair<Data*, unsigned>>& paired_pressure_data);
1199
1200 /// Redirect output to SphericalNavierStokesEquations output
1201 void output(std::ostream& outfile)
1202 {
1204 }
1205
1206 /// Redirect output to SphericalNavierStokesEquations output
1207 void output(std::ostream& outfile, const unsigned& nplot)
1208 {
1210 }
1211
1212 /// Redirect output to SphericalNavierStokesEquations output
1213 void output(FILE* file_pt)
1214 {
1216 }
1217
1218 /// Redirect output to SphericalNavierStokesEquations output
1219 void output(FILE* file_pt, const unsigned& nplot)
1220 {
1222 }
1223
1224
1225 /// The number of "DOF types" that degrees of freedom in this element
1226 /// are sub-divided into: Velocity (3 components) and pressure.
1227 unsigned ndof_types() const
1228 {
1229 return 4;
1230 }
1231
1232 /// Create a list of pairs for all unknowns in this element,
1233 /// so that the first entry in each pair contains the global equation
1234 /// number of the unknown, while the second one contains the number
1235 /// of the "Dof type" that this unknown is associated with.
1236 /// (Function can obviously only be called if the equation numbering
1237 /// scheme has been set up.)
1239 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const;
1240 };
1241
1242 // Inline functions
1243
1244 //==========================================================================
1245 /// 2D :
1246 /// Derivatives of the shape functions and test functions w.r.t to
1247 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1248 /// local and global coordinates.
1249 //==========================================================================
1252 Shape& psi,
1253 DShape& dpsidx,
1254 Shape& test,
1255 DShape& dtestdx) const
1256 {
1257 // Call the geometrical shape functions and derivatives
1258 double J = this->dshape_eulerian(s, psi, dpsidx);
1259 // Loop over the test functions and derivatives and set them equal to the
1260 // shape functions
1261 for (unsigned i = 0; i < 9; i++)
1262 {
1263 test[i] = psi[i];
1264 dtestdx(i, 0) = dpsidx(i, 0);
1265 dtestdx(i, 1) = dpsidx(i, 1);
1266 }
1267 // Return the jacobian
1268 return J;
1269 }
1270
1271
1272 //==========================================================================
1273 /// Derivatives of the shape functions and test functions w.r.t to
1274 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1275 /// local and global coordinates.
1276 //==========================================================================
1279 Shape& psi,
1280 DShape& dpsidx,
1281 Shape& test,
1282 DShape& dtestdx) const
1283 {
1284 // Call the geometrical shape functions and derivatives
1285 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
1286 // Loop over the test functions and derivatives and set them equal to the
1287 // shape functions
1288 for (unsigned i = 0; i < 9; i++)
1289 {
1290 test[i] = psi[i];
1291 dtestdx(i, 0) = dpsidx(i, 0);
1292 dtestdx(i, 1) = dpsidx(i, 1);
1293 }
1294 // Return the jacobian
1295 return J;
1296 }
1297
1298 //==========================================================================
1299 /// Pressure shape functions
1300 //==========================================================================
1302 const Vector<double>& s, Shape& psi) const
1303 {
1304 // Local storage
1305 double psi1[2], psi2[2];
1306 // Call the OneDimensional Shape functions
1307 OneDimLagrange::shape<2>(s[0], psi1);
1308 OneDimLagrange::shape<2>(s[1], psi2);
1309
1310 // Now let's loop over the nodal points in the element
1311 // s1 is the "x" coordinate, s2 the "y"
1312 for (unsigned i = 0; i < 2; i++)
1313 {
1314 for (unsigned j = 0; j < 2; j++)
1315 {
1316 /*Multiply the two 1D functions together to get the 2D function*/
1317 psi[2 * i + j] = psi2[i] * psi1[j];
1318 }
1319 }
1320 }
1321
1322
1323 //==========================================================================
1324 /// Pressure shape and test functions
1325 //==========================================================================
1327 const Vector<double>& s, Shape& psi, Shape& test) const
1328 {
1329 // Call the pressure shape functions
1330 pshape_spherical_nst(s, psi);
1331 // Loop over the test functions and set them equal to the shape functions
1332 for (unsigned i = 0; i < 4; i++) test[i] = psi[i];
1333 }
1334
1335
1336 //=======================================================================
1337 /// Face geometry of the Spherical Taylor_Hood elements
1338 //=======================================================================
1339 template<>
1341 : public virtual QElement<1, 3>
1342 {
1343 public:
1344 FaceGeometry() : QElement<1, 3>() {}
1345 };
1346
1347 //=======================================================================
1348 /// Face geometry of the FaceGeometry of the 2D Taylor Hoodelements
1349 //=======================================================================
1350 template<>
1352 : public virtual PointElement
1353 {
1354 public:
1356 };
1357
1358
1359} // namespace oomph
1360
1361#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:271
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Definition: nodes.h:293
/////////////////////////////////////////////////////////////////////////
Definition: fsi.h:63
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
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 double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3325
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
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 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
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3298
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
Definition: elements.cc:62
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:227
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data.
Definition: elements.h:267
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:192
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:5231
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: elements.h:3439
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void full_output(std::ostream &outfile, const unsigned &nplot)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format....
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
void output(std::ostream &outfile, const unsigned &nplot)
Redirect output to SphericalNavierStokesEquations output.
void output(std::ostream &outfile)
Redirect output to SphericalNavierStokesEquations output.
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
void output(FILE *file_pt, const unsigned &nplot)
Redirect output to SphericalNavierStokesEquations output.
int p_local_eqn(const unsigned &n) const
Return the local equation numbers for the pressure values.
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.
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
double dshape_and_dtest_eulerian_at_knot_spherical_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at ipt-th integation point...
void full_output(std::ostream &outfile)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format....
double p_spherical_nst(const unsigned &i) const
Return the i-th pressure value (Discontinous pressure interpolation – no need to cater for hanging no...
void output(FILE *file_pt)
Redirect output to SphericalNavierStokesEquations output.
QSphericalCrouzeixRaviartElement()
Constructor, there are 3 internal values (for the pressure)
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
unsigned npres_spherical_nst() const
Return number of pressure values.
void pshape_spherical_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocity (thr...
unsigned P_spherical_nst_internal_index
Internal index that indicates at which internal data the pressure is stored.
double dshape_and_dtest_eulerian_spherical_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...
/////////////////////////////////////////////////////////////////////////
QSphericalTaylorHoodElement()
Constructor, no internal data points.
double p_spherical_nst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version)
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocity (3 c...
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 get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
int p_local_eqn(const unsigned &n) const
Return the local equation numbers for the pressure values.
int p_nodal_index_spherical_nst() const
Set the value at which the pressure is stored in the nodes In this case the third index because there...
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.
void output(FILE *file_pt)
Redirect output to SphericalNavierStokesEquations output.
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
void output(std::ostream &outfile, const unsigned &nplot)
Redirect output to SphericalNavierStokesEquations output.
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.
void pshape_spherical_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
double dshape_and_dtest_eulerian_at_knot_spherical_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
double dshape_and_dtest_eulerian_spherical_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 output(FILE *file_pt, const unsigned &nplot)
Redirect output to SphericalNavierStokesEquations output.
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
void output(std::ostream &outfile)
Redirect output to SphericalNavierStokesEquations output.
unsigned npres_spherical_nst() const
Return number of pressure values.
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 elements that solve the Navier–Stokes equations, in axisymmetric spherical polar coordina...
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
double d_kin_energy_dt() const
Get integral of time derivative of kinetic energy over element.
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....
void get_vorticity(const Vector< double > &s, Vector< double > &vorticity) const
Compute the vorticity vector at local coordinate s.
SphericalNavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
double interpolated_dudx_spherical_nst(const Vector< double > &s, const unsigned &i, const unsigned &j) const
Return matrix entry dudx(i,j) of the FE interpolated velocity derivative at local coordinate s.
virtual double dshape_and_dtest_eulerian_at_knot_spherical_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...
double pressure_integral() const
Integral of pressure over element.
virtual void pshape_spherical_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
void compute_error_e(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, FiniteElement::SteadyExactSolutionFctPt exact_soln_dr_pt, FiniteElement::SteadyExactSolutionFctPt exact_soln_dtheta_pt, double &u_error, double &u_norm, double &p_error, double &p_norm)
Validate against exact solution. Solution is provided direct from exact_soln function....
Vector< double > actual_dth(const Vector< double > &x)
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
double *& re_invro_pt()
Pointer to global inverse Froude number.
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 unsigned u_index_spherical_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
SphericalNavierStokesSourceFctPt source_fct_pt() const
Access function for the source-function pointer. Const version.
virtual double dshape_and_dtest_eulerian_spherical_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....
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
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,...
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
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)
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
double u_spherical_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_sp...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Add the elemental contribution to the jacobian matrix. and the residuals vector. Note that this funct...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element's residual Vector.
virtual double p_spherical_nst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes.
double *& re_invfr_pt()
Pointer to global inverse Froude number.
Vector< double > * G_pt
Pointer to global gravity Vector.
Vector< double > actual_dr(const Vector< double > &x)
virtual void get_body_force_spherical_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...
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
double du_dt_spherical_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.
double * ReInvRo_pt
Pointer to global Reynolds number x inverse Rossby number (used when in a rotating frame)
void output(FILE *file_pt)
C-style output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
virtual int p_nodal_index_spherical_nst() const
Return the index at which the pressure is stored if it is stored at the nodes. If not stored at the n...
double * Re_pt
Pointer to global Reynolds number.
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
SphericalNavierStokesSourceFctPt & source_fct_pt()
Access function for the source-function pointer.
virtual void pshape_spherical_nst(const Vector< double > &s, Shape &psi, Shape &test) const =0
Compute the pressure shape and test functions at local coordinate s.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
SphericalNavierStokesBodyForceFctPt body_force_fct_pt() const
Access function for the body-force pointer. Const version.
double cot(const double &th) const
Include a cot function to simplify the NS momentum and jacobian expressions.
const Vector< double > & g() const
Vector of gravitational components.
double interpolated_p_spherical_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
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....
double kin_energy() const
Get integral of kinetic energy over element.
void output_vorticity(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z], [omega_x,omega_y,[and/or omega_z]] in tecplot format. nplot points in each ...
const double & re_invfr() const
Global inverse Froude number.
double dissipation() const
Return integral of dissipation over element.
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
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....
double(* SphericalNavierStokesSourceFctPt)(const double &time, const Vector< double > &x)
Function pointer to source function fct(t,x) x is a Vector!
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
double *& re_pt()
Pointer to Reynolds number.
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number)
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Add the element's contribution to its residuals vector, jacobian matrix and mass matrix.
virtual void fill_in_generic_residual_contribution_spherical_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Compute the residuals for the Navier–Stokes equations; flag=1(or 0): do (or don't) compute the Jacobi...
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
SphericalNavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
SphericalNavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
double u_spherical_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 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...
SphericalNavierStokesEquations()
Constructor: NULL the body force and source function and make sure the ALE terms are included by defa...
virtual unsigned npres_spherical_nst() const =0
Function to return number of pressure degrees of freedom.
double *& density_ratio_pt()
Pointer to Density ratio.
Vector< double > actual(const Vector< double > &x)
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.
double interpolated_u_spherical_nst(const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
void output_veloc(std::ostream &outfile, const unsigned &nplot, const unsigned &t)
Output function: x,y,[z],u,v,[w] in tecplot format. nplot points in each coordinate direction at time...
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...
void(* SphericalNavierStokesBodyForceFctPt)(const double &time, const Vector< double > &x, Vector< double > &body_force)
Function pointer to body force function fct(t,x,f(x)) x is a Vector!
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative....
const double & re() const
Reynolds number.
void interpolated_u_spherical_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
const double & re_invro() const
Global Reynolds number multiplied by inverse Rossby number.
virtual double get_source_spherical_nst(double time, const unsigned &ipt, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
std::string type() const
Return string that indicates the type of the timestepper (e.g. "BDF", "Newmark", etc....
Definition: timesteppers.h:490
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:608
//////////////////////////////////////////////////////////////////// ////////////////////////////////...