linearised_axisym_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 linearised axisymmetric Navier-Stokes elements
27
28#ifndef OOMPH_LINEARISED_AXISYM_NAVIER_STOKES_ELEMENTS_HEADER
29#define OOMPH_LINEARISED_AXISYM_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// oomph-lib includes
37#include "../generic/Qelements.h"
38#include "../generic/fsi.h"
39
40namespace oomph
41{
42 //=======================================================================
43 /// A class for elements that solve the linearised version of the
44 /// unsteady Navier--Stokes equations in cylindrical polar coordinates,
45 /// where we have Fourier-decomposed in the azimuthal direction so that
46 /// the theta-dependance is replaced by an azimuthal mode number.
47 //=======================================================================
49 : public virtual FiniteElement
50 {
51 private:
52 /// Static "magic" number that indicates that the pressure is not
53 /// stored at a node
55
56 /// Static default value for the physical constants
57 /// (all initialised to zero)
59
60 /// Static default value for the azimuthal mode number (zero)
62
63 /// Static default value for the physical ratios (all initialised to one)
65
66 protected:
67 // Physical constants
68 // ------------------
69
70 /// Pointer to the viscosity ratio (relative to the
71 /// viscosity used in the definition of the Reynolds number)
73
74 /// Pointer to the density ratio (relative to the
75 /// density used in the definition of the Reynolds number)
77
78 /// Pointer to global Reynolds number
79 double* Re_pt;
80
81 /// Pointer to global Reynolds number x Strouhal number (=Womersley)
82 double* ReSt_pt;
83
84 /// Pointer to azimuthal mode number k in e^ik(theta) decomposition
86
87 /// Pointer to base flow solution (velocity components) function
88 void (*Base_flow_u_fct_pt)(const double& time,
89 const Vector<double>& x,
90 Vector<double>& result);
91
92 /// Pointer to derivatives of base flow solution velocity
93 /// components w.r.t. global coordinates (r and z) function
94 void (*Base_flow_dudx_fct_pt)(const double& time,
95 const Vector<double>& x,
96 DenseMatrix<double>& result);
97
98 /// Boolean flag to indicate if ALE formulation is disabled when
99 /// the time-derivatives are computed. Only set to true if you're sure
100 /// that the mesh is stationary.
102
103 /// Access function for the local equation number
104 /// information for the i-th component of the pressure.
105 /// p_local_eqn[n,i] = local equation number or < 0 if pinned.
106 virtual int p_local_eqn(const unsigned& n, const unsigned& i) = 0;
107
108 /// Compute the shape functions and their derivatives
109 /// w.r.t. global coordinates at local coordinate s.
110 /// Return Jacobian of mapping between local and global coordinates.
112 const Vector<double>& s,
113 Shape& psi,
114 DShape& dpsidx,
115 Shape& test,
116 DShape& dtestdx) const = 0;
117
118 /// Compute the shape functions and their derivatives
119 /// w.r.t. global coordinates at the ipt-th integration point.
120 /// Return Jacobian of mapping between local and global coordinates.
122 const unsigned& ipt,
123 Shape& psi,
124 DShape& dpsidx,
125 Shape& test,
126 DShape& dtestdx) const = 0;
127
128 /// Compute the pressure shape functions at local coordinate s
130 Shape& psi) const = 0;
131
132 /// Compute the pressure shape and test functions at local coordinate s
134 Shape& psi,
135 Shape& test) const = 0;
136
137 /// Calculate the velocity components of the base flow solution
138 /// at a given time and Eulerian position
139 virtual void get_base_flow_u(const double& time,
140 const unsigned& ipt,
141 const Vector<double>& x,
142 Vector<double>& result) const
143 {
144 // If the function pointer is zero return zero
145 if (Base_flow_u_fct_pt == 0)
146 {
147 // Loop over velocity components and set base flow solution to zero
148 for (unsigned i = 0; i < 3; i++)
149 {
150 result[i] = 0.0;
151 }
152 }
153 // Otherwise call the function
154 else
155 {
156 (*Base_flow_u_fct_pt)(time, x, result);
157 }
158 }
159
160 /// Calculate the derivatives of the velocity components of the
161 /// base flow solution w.r.t. global coordinates (r and z) at a given
162 /// time and Eulerian position
163 virtual void get_base_flow_dudx(const double& time,
164 const unsigned& ipt,
165 const Vector<double>& x,
166 DenseMatrix<double>& result) const
167 {
168 // If the function pointer is zero return zero
169 if (Base_flow_dudx_fct_pt == 0)
170 {
171 // Loop over velocity components
172 for (unsigned i = 0; i < 3; i++)
173 {
174 // Loop over coordinate directions and set to zero
175 for (unsigned j = 0; j < 2; j++)
176 {
177 result(i, j) = 0.0;
178 }
179 }
180 }
181 // Otherwise call the function
182 else
183 {
184 (*Base_flow_dudx_fct_pt)(time, x, result);
185 }
186 }
187
188 /// Compute the residuals for the Navier-Stokes equations;
189 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
191 Vector<double>& residuals,
192 DenseMatrix<double>& jacobian,
193 DenseMatrix<double>& mass_matrix,
194 unsigned flag);
195
196 public:
197 /// Constructor: NULL the base flow solution and the
198 /// derivatives of the base flow function
201 {
202 // Set all the physical parameter pointers to the default value of zero
205
206 // Set the azimuthal mode number to default (zero)
208
209 // Set the physical ratios to the default value of one
212 }
213
214 /// Vector to decide whether the stress-divergence form is used or not.
215 // N.B. This needs to be public so that the intel compiler gets things
216 // correct. Somehow the access function messes things up when going to
217 // refineable navier--stokes
219
220 // Access functions for the physical constants
221 // -------------------------------------------
222
223 /// Reynolds number
224 const double& re() const
225 {
226 return *Re_pt;
227 }
228
229 /// Product of Reynolds and Strouhal number (=Womersley number)
230 const double& re_st() const
231 {
232 return *ReSt_pt;
233 }
234
235 /// Azimuthal mode number k in e^ik(theta) decomposition
236 const int& azimuthal_mode_number() const
237 {
239 }
240
241 /// Pointer to Reynolds number
242 double*& re_pt()
243 {
244 return Re_pt;
245 }
246
247 /// Pointer to product of Reynolds and Strouhal number (=Womersley number)
248 double*& re_st_pt()
249 {
250 return ReSt_pt;
251 }
252
253 /// Pointer to azimuthal mode number k in e^ik(theta) decomposition
255 {
257 }
258
259 /// Viscosity ratio for element: element's viscosity relative
260 /// to the viscosity used in the definition of the Reynolds number
261 const double& viscosity_ratio() const
262 {
263 return *Viscosity_Ratio_pt;
264 }
265
266 /// Pointer to the viscosity ratio
268 {
269 return Viscosity_Ratio_pt;
270 }
271
272 /// Density ratio for element: element's density relative
273 /// to the viscosity used in the definition of the Reynolds number
274 const double& density_ratio() const
275 {
276 return *Density_Ratio_pt;
277 }
278
279 /// Pointer to the density ratio
281 {
282 return Density_Ratio_pt;
283 }
284
285 /// Access function for the base flow solution pointer
286 void (*&base_flow_u_fct_pt())(const double& time,
287 const Vector<double>& x,
289 {
290 return Base_flow_u_fct_pt;
291 }
292
293 /// Access function for the derivatives of the base flow
294 /// w.r.t. global coordinates solution pointer
295 void (*&base_flow_dudx_fct_pt())(const double& time,
296 const Vector<double>& x,
298 {
300 }
301
302 /// Return the number of pressure degrees of freedom
303 /// associated with a single pressure component in the element
304 virtual unsigned npres_linearised_axi_nst() const = 0;
305
306 /// Return the index at which the i-th unknown velocity
307 /// component is stored. The default value, i, is appropriate for
308 /// single-physics problems. In derived multi-physics elements, this
309 /// function should be overloaded to reflect the chosen storage scheme.
310 /// Note that these equations require that the unknowns are always
311 /// stored at the same indices at each node.
312 virtual inline unsigned u_index_linearised_axi_nst(const unsigned& i) const
313 {
314 return i;
315 }
316
317 /// Return the i-th component of du/dt at local node n.
318 /// Uses suitably interpolated value for hanging nodes.
319 double du_dt_linearised_axi_nst(const unsigned& n, const unsigned& i) const
320 {
321 // Get the data's timestepper
323
324 // Initialise dudt
325 double dudt = 0.0;
326
327 // Loop over the timesteps, if there is a non-steady timestepper
329 {
330 // Get the index at which the velocity is stored
331 const unsigned u_nodal_index = u_index_linearised_axi_nst(i);
332
333 // Determine number of timsteps (past & present)
334 const unsigned n_time = time_stepper_pt->ntstorage();
335
336 // Add the contributions to the time derivative
337 for (unsigned t = 0; t < n_time; t++)
338 {
339 dudt +=
340 time_stepper_pt->weight(1, t) * nodal_value(t, n, u_nodal_index);
341 }
342 }
343
344 return dudt;
345 }
346
347 /// Disable ALE, i.e. assert the mesh is not moving -- you do this
348 /// at your own risk!
350 {
351 ALE_is_disabled = true;
352 }
353
354 /// (Re-)enable ALE, i.e. take possible mesh motion into account
355 /// when evaluating the time-derivative. Note: By default, ALE is
356 /// enabled, at the expense of possibly creating unnecessary work
357 /// in problems where the mesh is, in fact, stationary.
359 {
360 ALE_is_disabled = false;
361 }
362
363 /// Return the i-th pressure value at local pressure "node" n_p.
364 /// Uses suitably interpolated value for hanging nodes.
365 virtual double p_linearised_axi_nst(const unsigned& n_p,
366 const unsigned& i) const = 0;
367
368 /// Which nodal value represents the pressure?
369 // N.B. This function has return type "int" (rather than "unsigned"
370 // as in the u_index case) so that we can return the "magic" number
371 // "Pressure_not_stored_at_node" ( = -100 )
372 virtual inline int p_index_linearised_axi_nst(const unsigned& i) const
373 {
375 }
376
377 /// Strain-rate tensor: \f$ e_{ij} \f$
378 /// where \f$ i,j = r,z,\theta \f$ (in that order)
379 void strain_rate(const Vector<double>& s,
381
382 /// Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
383 /// in tecplot format. Default number of plot points
384 void output(std::ostream& outfile)
385 {
386 const unsigned nplot = 5;
387 output(outfile, nplot);
388 }
389
390 /// Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
391 /// in tecplot format. nplot points in each coordinate direction
392 void output(std::ostream& outfile, const unsigned& nplot);
393
394 /// Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
395 /// in tecplot format. Default number of plot points
396 void output(FILE* file_pt)
397 {
398 const unsigned nplot = 5;
399 output(file_pt, nplot);
400 }
401
402 /// Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S
403 /// in tecplot format. nplot points in each coordinate direction
404 void output(FILE* file_pt, const unsigned& nplot);
405
406 /// Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S,
407 /// in tecplot format. nplot points in each coordinate direction
408 /// at timestep t (t=0: present; t>0: previous timestep)
409 void output_veloc(std::ostream& outfile,
410 const unsigned& nplot,
411 const unsigned& t);
412
413 /// Compute the element's residual Vector
415 {
416 // Call the generic residuals function with flag set to 0
417 // and using a dummy matrix argument
419 residuals,
422 0);
423 }
424
425 /// Compute the element's residual Vector and the jacobian matrix.
426 /// Virtual function can be overloaded by hanging-node version.
428 DenseMatrix<double>& jacobian)
429 {
430 // Call the generic routine with the flag set to 1
432 residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
433 }
434
435 /// Add the element's contribution to its residuals vector,
436 /// jacobian matrix and mass matrix
438 Vector<double>& residuals,
439 DenseMatrix<double>& jacobian,
440 DenseMatrix<double>& mass_matrix)
441 {
442 // Call the generic routine with the flag set to 2
444 residuals, jacobian, mass_matrix, 2);
445 }
446
447 /// Return the i-th component of the FE interpolated velocity
448 /// u[i] at local coordinate s
450 const unsigned& i) const
451 {
452 // Determine number of nodes in the element
453 const unsigned n_node = nnode();
454
455 // Provide storage for local shape functions
456 Shape psi(n_node);
457
458 // Find values of shape functions
459 shape(s, psi);
460
461 // Get the index at which the velocity is stored
462 const unsigned u_nodal_index = u_index_linearised_axi_nst(i);
463
464 // Initialise value of u
465 double interpolated_u = 0.0;
466
467 // Loop over the local nodes and sum
468 for (unsigned l = 0; l < n_node; l++)
469 {
470 interpolated_u += nodal_value(l, u_nodal_index) * psi[l];
471 }
472
473 return (interpolated_u);
474 }
475
476 /// Return the i-th component of the FE interpolated pressure
477 /// p[i] at local coordinate s
479 const unsigned& i) const
480 {
481 // Determine number of pressure nodes in the element
482 const unsigned n_pressure_nodes = npres_linearised_axi_nst();
483
484 // Provide storage for local shape functions
485 Shape psi(n_pressure_nodes);
486
487 // Find values of shape functions
489
490 // Initialise value of p
491 double interpolated_p = 0.0;
492
493 // Loop over the local nodes and sum
494 for (unsigned l = 0; l < n_pressure_nodes; l++)
495 {
496 // N.B. The pure virtual function p_linearised_axi_nst(...)
497 // automatically calculates the index at which the pressure value
498 // is stored, so we don't need to worry about this here
499 interpolated_p += p_linearised_axi_nst(l, i) * psi[l];
500 }
501
502 return (interpolated_p);
503 }
504
505 }; // End of LinearisedAxisymmetricNavierStokesEquations class definition
506
507
508 /// ///////////////////////////////////////////////////////////////////////////
509 /// ///////////////////////////////////////////////////////////////////////////
510 /// ///////////////////////////////////////////////////////////////////////////
511
512
513 //=======================================================================
514 /// Crouzeix-Raviart elements are Navier-Stokes elements with quadratic
515 /// interpolation for velocities and positions, but a discontinuous
516 /// linear pressure interpolation
517 //=======================================================================
519 : public virtual QElement<2, 3>,
521 {
522 private:
523 /// Static array of ints to hold required number of variables at nodes
524 static const unsigned Initial_Nvalue[];
525
526 protected:
527 /// Internal indices that indicate at which internal data the
528 /// pressure values are stored. We note that there are two pressure
529 /// values, corresponding to the functions P^C(r,z,t) and P^S(r,z,t)
530 /// which multiply the cosine and sine terms respectively.
532
533 /// Velocity shape and test functions and their derivatives
534 /// w.r.t. global coordinates at local coordinate s (taken from geometry).
535 /// Return Jacobian of mapping between local and global coordinates.
537 const Vector<double>& s,
538 Shape& psi,
539 DShape& dpsidx,
540 Shape& test,
541 DShape& dtestdx) const;
542
543 /// Velocity shape and test functions and their derivatives
544 /// w.r.t. global coordinates at the ipt-th integation point
545 /// (taken from geometry).
546 /// Return Jacobian of mapping between local and global coordinates.
548 const unsigned& ipt,
549 Shape& psi,
550 DShape& dpsidx,
551 Shape& test,
552 DShape& dtestdx) const;
553
554 /// Compute the pressure shape functions at local coordinate s
555 inline void pshape_linearised_axi_nst(const Vector<double>& s,
556 Shape& psi) const;
557
558 /// Compute the pressure shape and test functions at local coordinate s
559 inline void pshape_linearised_axi_nst(const Vector<double>& s,
560 Shape& psi,
561 Shape& test) const;
562
563 public:
564 /// Constructor: there are three internal values for each
565 /// of the two pressure components
567 : QElement<2, 3>(),
570 {
571 // Loop over the two pressure components
572 for (unsigned i = 0; i < 2; i++)
573 {
574 // Allocate and add one internal data object for each of the two
575 // pressure components that store the three pressure values
577 this->add_internal_data(new Data(3));
578 }
579 }
580
581 /// Return number of values (pinned or dofs) required at local node n
582 virtual unsigned required_nvalue(const unsigned& n) const;
583
584 /// Return the pressure value i at internal dof i_internal
585 /// (Discontinous pressure interpolation -- no need to cater for
586 /// hanging nodes)
587 double p_linearised_axi_nst(const unsigned& i_internal,
588 const unsigned& i) const
589 {
591 ->value(i_internal);
592 }
593
594 /// Return number of pressure values corresponding to a
595 /// single pressure component
597 {
598 return 3;
599 }
600
601 /// Fix both components of the internal pressure degrees
602 /// of freedom p_dof to pvalue
603 void fix_pressure(const unsigned& p_dof, const double& pvalue)
604 {
605 // Loop over the two pressure components
606 for (unsigned i = 0; i < 2; i++)
607 {
608 this->internal_data_pt(P_linearised_axi_nst_internal_index[i])
609 ->pin(p_dof);
611 ->set_value(p_dof, pvalue);
612 }
613 }
614
615 /// Overload the access function for the i-th component of the
616 /// pressure's local equation numbers
617 inline int p_local_eqn(const unsigned& n, const unsigned& i)
618 {
620 }
621
622 /// Redirect output to NavierStokesEquations output
623 void output(std::ostream& outfile)
624 {
626 }
627
628 /// Redirect output to NavierStokesEquations output
629 void output(std::ostream& outfile, const unsigned& n_plot)
630 {
632 }
633
634 /// Redirect output to NavierStokesEquations output
635 void output(FILE* file_pt)
636 {
638 }
639
640 /// Redirect output to NavierStokesEquations output
641 void output(FILE* file_pt, const unsigned& n_plot)
642 {
644 }
645
646 /// The number of "dof-blocks" that degrees of freedom in this
647 /// element are sub-divided into: Velocity and pressure.
648 unsigned ndof_types() const
649 {
650 return 8;
651 }
652
653 }; // End of LinearisedAxisymmetricQCrouzeixRaviartElement class definition
654
655
656 // Inline functions
657 // ----------------
658
659 //=======================================================================
660 /// Derivatives of the shape functions and test functions w.r.t.
661 /// global (Eulerian) coordinates at local coordinate s.
662 /// Return Jacobian of mapping between local and global coordinates.
663 //=======================================================================
666 Shape& psi,
667 DShape& dpsidx,
668 Shape& test,
669 DShape& dtestdx) const
670 {
671 // Call the geometrical shape functions and derivatives
672 const double J = this->dshape_eulerian(s, psi, dpsidx);
673
674 // Loop over the test functions and derivatives and set them
675 // equal to the shape functions
676 for (unsigned i = 0; i < 9; i++)
677 {
678 test[i] = psi[i];
679 dtestdx(i, 0) = dpsidx(i, 0);
680 dtestdx(i, 1) = dpsidx(i, 1);
681 }
682
683 // Return the Jacobian
684 return J;
685 }
686
687 //=======================================================================
688 /// Derivatives of the shape functions and test functions w.r.t.
689 /// global (Eulerian) coordinates at the ipt-th integration point.
690 /// Return Jacobian of mapping between local and global coordinates.
691 //=======================================================================
694 Shape& psi,
695 DShape& dpsidx,
696 Shape& test,
697 DShape& dtestdx) const
698 {
699 // Call the geometrical shape functions and derivatives
700 const double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
701
702 // Loop over the test functions and derivatives and set them
703 // equal to the shape functions
704 for (unsigned i = 0; i < 9; i++)
705 {
706 test[i] = psi[i];
707 dtestdx(i, 0) = dpsidx(i, 0);
708 dtestdx(i, 1) = dpsidx(i, 1);
709 }
710
711 // Return the Jacobian
712 return J;
713 }
714
715 //=======================================================================
716 /// Pressure shape functions
717 //=======================================================================
720 {
721 psi[0] = 1.0;
722 psi[1] = s[0];
723 psi[2] = s[1];
724 }
725
726 //=======================================================================
727 /// Define the pressure shape and test functions
728 //=======================================================================
731 Shape& psi,
732 Shape& test) const
733 {
734 // Call the pressure shape functions
736
737 // Loop over the test functions and set them equal to the shape functions
738 for (unsigned i = 0; i < 3; i++)
739 {
740 test[i] = psi[i];
741 }
742 }
743
744 //=======================================================================
745 /// Face geometry of the linearised axisym Crouzeix-Raviart elements
746 //=======================================================================
747 template<>
749 : public virtual QElement<1, 3>
750 {
751 public:
752 FaceGeometry() : QElement<1, 3>() {}
753 };
754
755 //=======================================================================
756 /// Face geometry of face geometry of the linearised axisymmetric
757 /// Crouzeix Raviart elements
758 //=======================================================================
759 template<>
762 : public virtual PointElement
763 {
764 public:
766 };
767
768
769 /// ///////////////////////////////////////////////////////////////////////////
770 /// ///////////////////////////////////////////////////////////////////////////
771 /// ///////////////////////////////////////////////////////////////////////////
772
773
774 //=======================================================================
775 /// Taylor--Hood elements are Navier--Stokes elements with quadratic
776 /// interpolation for velocities and positions and continuous linear
777 /// pressure interpolation
778 //=======================================================================
780 : public virtual QElement<2, 3>,
782 {
783 private:
784 /// Static array of ints to hold number of variables at node
785 static const unsigned Initial_Nvalue[];
786
787 protected:
788 /// Static array of ints to hold conversion from pressure
789 /// node numbers to actual node numbers
790 static const unsigned Pconv[];
791
792 /// Velocity shape and test functions and their derivatives
793 /// w.r.t. global coordinates at local coordinate s (taken from geometry).
794 /// Return Jacobian of mapping between local and global coordinates.
796 const Vector<double>& s,
797 Shape& psi,
798 DShape& dpsidx,
799 Shape& test,
800 DShape& dtestdx) const;
801
802 /// Velocity shape and test functions and their derivatives
803 /// w.r.t. global coordinates the ipt-th integation point
804 /// (taken from geometry).
805 /// Return Jacobian of mapping between local and global coordinates.
807 const unsigned& ipt,
808 Shape& psi,
809 DShape& dpsidx,
810 Shape& test,
811 DShape& dtestdx) const;
812
813 /// Compute the pressure shape functions at local coordinate s
814 inline void pshape_linearised_axi_nst(const Vector<double>& s,
815 Shape& psi) const;
816
817 /// Compute the pressure shape and test functions at local coordinte s
818 inline void pshape_linearised_axi_nst(const Vector<double>& s,
819 Shape& psi,
820 Shape& test) const;
821
822 public:
823 /// Constructor, no internal data points
826 {
827 }
828
829 /// Number of values (pinned or dofs) required at node n. Can
830 /// be overwritten for hanging node version
831 inline virtual unsigned required_nvalue(const unsigned& n) const
832 {
833 return Initial_Nvalue[n];
834 }
835
836 /// Which nodal value represents the pressure? Overload version
837 /// in base class which returns static int "Pressure_not_stored_at_node"
838 virtual int p_index_linearised_axi_nst(const unsigned& i) const
839 {
840 return (6 + i);
841 }
842
843 /// Access function for the i-th component of pressure
844 /// at local pressure node n_p (const version).
845 double p_linearised_axi_nst(const unsigned& n_p, const unsigned& i) const
846 {
848 }
849
850 /// Return number of pressure values corresponding to a
851 /// single pressure component
853 {
854 return 4;
855 }
856
857 /// Fix both components of the pressure at local pressure
858 /// node n_p to pvalue
859 void fix_pressure(const unsigned& n_p, const double& pvalue)
860 {
861 // Loop over the two pressure components
862 for (unsigned i = 0; i < 2; i++)
863 {
864 this->node_pt(Pconv[n_p])->pin(p_index_linearised_axi_nst(i));
865 this->node_pt(Pconv[n_p])
867 }
868 }
869
870 /// Overload the access function for the i-th component of the
871 /// pressure's local equation numbers
872 inline int p_local_eqn(const unsigned& n, const unsigned& i)
873 {
875 }
876
877 /// Redirect output to NavierStokesEquations output
878 void output(std::ostream& outfile)
879 {
881 }
882
883 /// Redirect output to NavierStokesEquations output
884 void output(std::ostream& outfile, const unsigned& n_plot)
885 {
887 }
888
889 /// Redirect output to NavierStokesEquations output
890 void output(FILE* file_pt)
891 {
893 }
894
895 /// Redirect output to NavierStokesEquations output
896 void output(FILE* file_pt, const unsigned& n_plot)
897 {
899 }
900
901 /// Returns the number of "dof-blocks" that degrees of freedom
902 /// in this element are sub-divided into: Velocity and pressure.
903 unsigned ndof_types() const
904 {
905 return 8;
906 }
907
908 }; // End of LinearisedAxisymmetricQTaylorHoodElement class definition
909
910
911 // Inline functions
912 // ----------------
913
914 //=======================================================================
915 /// Derivatives of the shape functions and test functions w.r.t
916 /// global (Eulerian) coordinates at local coordinate s.
917 /// Return Jacobian of mapping between local and global coordinates.
918 //=======================================================================
921 Shape& psi,
922 DShape& dpsidx,
923 Shape& test,
924 DShape& dtestdx) const
925 {
926 // Call the geometrical shape functions and derivatives
927 const double J = this->dshape_eulerian(s, psi, dpsidx);
928
929 // Loop over the test functions and derivatives and set them
930 // equal to the shape functions
931 for (unsigned i = 0; i < 9; i++)
932 {
933 test[i] = psi[i];
934 dtestdx(i, 0) = dpsidx(i, 0);
935 dtestdx(i, 1) = dpsidx(i, 1);
936 }
937
938 // Return the Jacobian
939 return J;
940 }
941
942 //=======================================================================
943 /// Derivatives of the shape functions and test functions w.r.t
944 /// global (Eulerian) coordinates at the ipt-th integration point.
945 /// Return Jacobian of mapping between local and global coordinates.
946 //=======================================================================
949 Shape& psi,
950 DShape& dpsidx,
951 Shape& test,
952 DShape& dtestdx) const
953 {
954 // Call the geometrical shape functions and derivatives
955 const double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
956
957 // Loop over the test functions and derivatives and set them
958 // equal to the shape functions
959 for (unsigned i = 0; i < 9; i++)
960 {
961 test[i] = psi[i];
962 dtestdx(i, 0) = dpsidx(i, 0);
963 dtestdx(i, 1) = dpsidx(i, 1);
964 }
965
966 // Return the Jacobian
967 return J;
968 }
969
970 //=======================================================================
971 /// Pressure shape functions
972 //=======================================================================
975 {
976 // Allocate local storage for the pressure shape functions
977 double psi1[2], psi2[2];
978
979 // Call the one-dimensional shape functions
980 OneDimLagrange::shape<2>(s[0], psi1);
981 OneDimLagrange::shape<2>(s[1], psi2);
982
983 // Now let's loop over the nodal points in the element
984 // s1 is the "r" coordinate, s2 the "z"
985 for (unsigned i = 0; i < 2; i++)
986 {
987 for (unsigned j = 0; j < 2; j++)
988 {
989 // Multiply the two 1D functions together to get the 2D function
990 psi[2 * i + j] = psi2[i] * psi1[j];
991 }
992 }
993 }
994
995 //=======================================================================
996 /// Pressure shape and test functions
997 //=======================================================================
1000 Shape& psi,
1001 Shape& test) const
1002 {
1003 // Call the pressure shape functions
1005
1006 // Loop over the test functions and set them equal to the shape functions
1007 for (unsigned i = 0; i < 4; i++)
1008 {
1009 test[i] = psi[i];
1010 }
1011 }
1012
1013 //=======================================================================
1014 /// Face geometry of the linearised axisymmetric Taylor Hood elements
1015 //=======================================================================
1016 template<>
1018 : public virtual QElement<1, 3>
1019 {
1020 public:
1021 FaceGeometry() : QElement<1, 3>() {}
1022 };
1023
1024 //=======================================================================
1025 /// Face geometry of the face geometry of the linearised
1026 /// axisymmetric Taylor Hood elements
1027 //=======================================================================
1028 template<>
1030 : public virtual PointElement
1031 {
1032 public:
1034 };
1035
1036
1037} // namespace oomph
1038
1039#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: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
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
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
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
A class for elements that solve the linearised version of the unsteady Navier–Stokes equations in cyl...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element's residual Vector.
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
virtual double dshape_and_dtest_eulerian_linearised_axi_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and their derivatives w.r.t. global coordinates at local coordinate s....
void output(FILE *file_pt)
Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S in tecplot format. Default number of pl...
void(* Base_flow_dudx_fct_pt)(const double &time, const Vector< double > &x, DenseMatrix< double > &result)
Pointer to derivatives of base flow solution velocity components w.r.t. global coordinates (r and z) ...
const int & azimuthal_mode_number() const
Azimuthal mode number k in e^ik(theta) decomposition.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Add the element's contribution to its residuals vector, jacobian matrix and mass matrix.
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
static int Default_Azimuthal_Mode_Number_Value
Static default value for the azimuthal mode number (zero)
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
int * Azimuthal_Mode_Number_pt
Pointer to azimuthal mode number k in e^ik(theta) decomposition.
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
void disable_ALE()
Disable ALE, i.e. assert the mesh is not moving – you do this at your own risk!
double du_dt_linearised_axi_nst(const unsigned &n, const unsigned &i) const
Return the i-th component of du/dt at local node n. Uses suitably interpolated value for hanging node...
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 ...
void(*&)(const double &time, const Vector< double > &x, DenseMatrix< double > &f) base_flow_dudx_fct_pt()
Access function for the derivatives of the base flow w.r.t. global coordinates solution pointer.
void output(std::ostream &outfile)
Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, P^C, P^S in tecplot format. Default number of pl...
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all initialised to one)
virtual void get_base_flow_u(const double &time, const unsigned &ipt, const Vector< double > &x, Vector< double > &result) const
Calculate the velocity components of the base flow solution at a given time and Eulerian position.
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
virtual int p_local_eqn(const unsigned &n, const unsigned &i)=0
Access function for the local equation number information for the i-th component of the pressure....
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed....
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
void enable_ALE()
(Re-)enable ALE, i.e. take possible mesh motion into account when evaluating the time-derivative....
virtual void pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
double interpolated_u_linearised_axi_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated velocity u[i] at local coordinate s.
void output_veloc(std::ostream &outfile, const unsigned &nplot, const unsigned &t)
Output function: r, z, U^C, U^S, V^C, V^S, W^C, W^S, in tecplot format. nplot points in each coordina...
virtual void get_base_flow_dudx(const double &time, const unsigned &ipt, const Vector< double > &x, DenseMatrix< double > &result) const
Calculate the derivatives of the velocity components of the base flow solution w.r....
virtual double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const =0
Compute the shape functions and their derivatives w.r.t. global coordinates at the ipt-th integration...
virtual unsigned npres_linearised_axi_nst() const =0
Return the number of pressure degrees of freedom associated with a single pressure component in the e...
double interpolated_p_linearised_axi_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated pressure p[i] at local coordinate s.
const double & viscosity_ratio() const
Viscosity ratio for element: element's viscosity relative to the viscosity used in the definition of ...
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
int *& azimuthal_mode_number_pt()
Pointer to azimuthal mode number k in e^ik(theta) decomposition.
virtual unsigned u_index_linearised_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
void(* Base_flow_u_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to base flow solution (velocity components) function.
virtual void pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi, Shape &test) const =0
Compute the pressure shape and test functions at local coordinate s.
virtual void fill_in_generic_residual_contribution_linearised_axi_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...
virtual double p_linearised_axi_nst(const unsigned &n_p, const unsigned &i) const =0
Return the i-th pressure value at local pressure "node" n_p. Uses suitably interpolated value for han...
virtual int p_index_linearised_axi_nst(const unsigned &i) const
Which nodal value represents the pressure?
LinearisedAxisymmetricNavierStokesEquations()
Constructor: NULL the base flow solution and the derivatives of the base flow function.
const double & density_ratio() const
Density ratio for element: element's density relative to the viscosity used in the definition of the ...
void(*&)(const double &time, const Vector< double > &x, Vector< double > &f) base_flow_u_fct_pt()
Access function for the base flow solution pointer.
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void output(std::ostream &outfile, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
unsigned ndof_types() const
The number of "dof-blocks" that degrees of freedom in this element are sub-divided into: Velocity and...
double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivatives w.r.t. global coordinates at the ipt-th integ...
double p_linearised_axi_nst(const unsigned &i_internal, const unsigned &i) const
Return the pressure value i at internal dof i_internal (Discontinous pressure interpolation – no need...
Vector< unsigned > P_linearised_axi_nst_internal_index
Internal indices that indicate at which internal data the pressure values are stored....
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
LinearisedAxisymmetricQCrouzeixRaviartElement()
Constructor: there are three internal values for each of the two pressure components.
double dshape_and_dtest_eulerian_linearised_axi_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivatives w.r.t. global coordinates at local coordinate...
virtual unsigned required_nvalue(const unsigned &n) const
Return number of values (pinned or dofs) required at local node n.
void pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi) const
Compute the pressure shape functions at local coordinate s.
void fix_pressure(const unsigned &p_dof, const double &pvalue)
Fix both components of the internal pressure degrees of freedom p_dof to pvalue.
unsigned npres_linearised_axi_nst() const
Return number of pressure values corresponding to a single pressure component.
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
int p_local_eqn(const unsigned &n, const unsigned &i)
Overload the access function for the i-th component of the pressure's local equation numbers.
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
unsigned ndof_types() const
Returns the number of "dof-blocks" that degrees of freedom in this element are sub-divided into: Velo...
double dshape_and_dtest_eulerian_at_knot_linearised_axi_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivatives w.r.t. global coordinates the ipt-th integati...
void fix_pressure(const unsigned &n_p, const double &pvalue)
Fix both components of the pressure at local pressure node n_p to pvalue.
void pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi) const
Compute the pressure shape functions at local coordinate s.
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.
double p_linearised_axi_nst(const unsigned &n_p, const unsigned &i) const
Access function for the i-th component of pressure at local pressure node n_p (const version).
void output(std::ostream &outfile, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
unsigned npres_linearised_axi_nst() const
Return number of pressure values corresponding to a single pressure component.
virtual int p_index_linearised_axi_nst(const unsigned &i) const
Which nodal value represents the pressure? Overload version in base class which returns static int "P...
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
double dshape_and_dtest_eulerian_linearised_axi_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivatives w.r.t. global coordinates at local coordinate...
LinearisedAxisymmetricQTaylorHoodElement()
Constructor, no internal data points.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
int p_local_eqn(const unsigned &n, const unsigned &i)
Overload the access function for the i-th component of the pressure's local equation numbers.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: elements.h:3439
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Definition: Qelements.h:459
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
////////////////////////////////////////////////////////////////////// //////////////////////////////...
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
bool is_steady() const
Flag to indicate if a timestepper has been made steady (possibly temporarily to switch off time-depen...
Definition: timesteppers.h:389
void shape< 2 >(const double &s, double *Psi)
1D shape functions specialised to linear order (2 Nodes)
Definition: shape.h:608
//////////////////////////////////////////////////////////////////// ////////////////////////////////...