polar_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// Created (18/03/08) by recopying from my_oomph-codes/polar_navier_stokes.h
27// Now I know what I'm doing it needs updating to current oomph conventions.
28
29// 17/06/08 - Weakform adjusted to "correct" traction version
30
31// "_pnst" added to: - npres()
32// - pshape()
33// - u()
34// - p()
35// - du_dt()
36// - interpolated_u()
37// - interpolated_p()
38// - interpolated_dudx()
39// - dshape_and_dtest_eulerian()
40// - dshape_and_dtest_eulerian_at_knot()
41
42// The following generality is necessary for elements with multiple physics.
43// That is where we can't just assume that values one and two at the nodes
44// are velocities and the third a pressure.
45
46// Then renamed "index_of_nodal_pressure_value to p_nodal_index_pnst.
47
48// Added u_index_pnst (by default just returns the value you give it)
49
50// The internal pressure data in Crouzeix Raviart elements is now stored
51// as one data object with three values:
52// Added P_pnst_internal_index which stores the index of that internal
53// data, specified in the constructor.
54// Adjusted: - p_pnst
55// - fix_pressure
56// - get_load_data
57// To reflect this change in internal data storage.
58
59// Plus added Pressure_not_stored_at_node, default = -100.
60
61// assign_additional_local_eqn_numbers removed in favour of
62// - p_local_eqn
63// And we have:
64// - local_eqn = nodal_local_eqn(l,u_nodal_index[i]);
65// - local_eqn = p_local_eqn(l);
66// in fill_in_generic_residual_contribution
67
68// Also removed the following functions for pinning/unpinning pressures
69// as they're only needed in the refineable case:
70//----------------------------------------------------------------------
71// - unpin_all_nodal_pressure_dofs()
72// - unpin_all_internal_pressure_dofs()
73// - pin_all_nodal_pressure_dofs()
74// - unpin_proper_nodal_pressure_dofs()
75// - pin_redundant_nodal_pressures()
76// - unpin_all_pressure_dofs()
77// - pressure_dof_is_hanging()
78//----------------------------------------------------------------------
79
80// strain_rate adapted to return the polar strain components.
81
82// interpolated positions / velocities now assembled using nodal_position
83// and nodal_value - both of which should cope with hanging nodes.
84
85// Also stokes_output() added to compute convergence rates when exact
86// (Stokes) solution is known.
87
88#ifndef OOMPH_POLAR_NAVIER_STOKES
89#define OOMPH_POLAR_NAVIER_STOKES
90
91// Config header generated by autoconfig
92#ifdef HAVE_CONFIG_H
93#include <oomph-lib-config.h>
94#endif
95
96// OOMPH-LIB headers
97#include "../generic/Qelements.h"
98
99
100namespace oomph
101{
102 //=====================================================================
103 /// A class for elements that solve the polar Navier--Stokes equations,
104 /// This contains the generic maths -- any concrete implementation must
105 /// be derived from this.
106 //======================================================================
108 {
109 public:
110 /// Function pointer to body force function fct(t,x,f(x))
111 /// x is a Vector!
112 typedef void (*NavierStokesBodyForceFctPt)(const double& time,
113 const Vector<double>& x,
114 Vector<double>& body_force);
115
116 /// Function pointer to source function fct(t,x)
117 /// x is a Vector!
118 typedef double (*NavierStokesSourceFctPt)(const double& time,
119 const Vector<double>& x);
120
121 private:
122 /// Static "magic" number that indicates that the pressure is
123 /// not stored at a node
125
126 /// Static default value for the physical constants (all initialised to
127 /// zero)
129
130 /// Static default value for the physical ratios (all are initialised to
131 /// one)
133
134 /// Static default value for the gravity vector
136
137 protected:
138 // Physical constants
139
140 /// Pointer to the viscosity ratio (relative to the
141 /// viscosity used in the definition of the Reynolds number)
143
144 /// Pointer to the density ratio (relative to the density used in the
145 /// definition of the Reynolds number)
147
148 // Pointers to global physical constants
149
150 /// Pointer to the angle alpha
151 double* Alpha_pt;
152
153 /// Pointer to global Reynolds number
154 double* Re_pt;
155
156 /// Pointer to global Reynolds number x Strouhal number (=Womersley)
157 double* ReSt_pt;
158
159 /// Pointer to global Reynolds number x inverse Froude number
160 /// (= Bond number / Capillary number)
161 double* ReInvFr_pt;
162
163 /// Pointer to global gravity Vector
165
166 /// Pointer to body force function
168
169 /// Pointer to volumetric source function
171
172 /// Access function for the local equation number information for
173 /// the pressure.
174 /// p_local_eqn[n] = local equation number or < 0 if pinned
175 virtual int p_local_eqn(const unsigned& n) = 0;
176
177 /// Compute the shape functions and derivatives
178 /// w.r.t. global coords at local coordinate s.
179 /// Return Jacobian of mapping between local and global coordinates.
181 Shape& psi,
182 DShape& dpsidx,
183 Shape& test,
184 DShape& dtestdx) const = 0;
185
186 /// Compute the shape functions and derivatives
187 /// w.r.t. global coords at ipt-th integration point
188 /// Return Jacobian of mapping between local and global coordinates.
190 const unsigned& ipt,
191 Shape& psi,
192 DShape& dpsidx,
193 Shape& test,
194 DShape& dtestdx) const = 0;
195
196 /// Compute the pressure shape functions at local coordinate s
197 virtual void pshape_pnst(const Vector<double>& s, Shape& psi) const = 0;
198
199 /// Compute the pressure shape and test functions
200 /// at local coordinate s
201 virtual void pshape_pnst(const Vector<double>& s,
202 Shape& psi,
203 Shape& test) const = 0;
204
205
206 /// Calculate the body force at a given time and Eulerian position
207 void get_body_force(double time,
208 const Vector<double>& x,
209 Vector<double>& result)
210 {
211 // If the function pointer is zero return zero
212 if (Body_force_fct_pt == 0)
213 {
214 // Loop over dimensions and set body forces to zero
215 for (unsigned i = 0; i < 2; i++)
216 {
217 result[i] = 0.0;
218 }
219 }
220 // Otherwise call the function
221 else
222 {
223 (*Body_force_fct_pt)(time, x, result);
224 }
225 }
226
227 /// Calculate the source fct at given time and
228 /// Eulerian position
229 double get_source_fct(double time, const Vector<double>& x)
230 {
231 // If the function pointer is zero return zero
232 if (Source_fct_pt == 0)
233 {
234 return 0;
235 }
236 // Otherwise call the function
237 else
238 {
239 return (*Source_fct_pt)(time, x);
240 }
241 }
242
243 public:
244 /// Constructor: NULL the body force and source function
246 {
247 // Set all the Physical parameter pointers to the default value zero
252 // Set the Physical ratios to the default value of 1
255 }
256
257 /// Vector to decide whether the stress-divergence form is used or not
258 // N.B. This needs to be public so that the intel compiler gets things
259 // correct somehow the access function messes things up when going to
260 // refineable navier--stokes
262
263 // Access functions for the physical constants
264
265 /// Reynolds number
266 const double& re() const
267 {
268 return *Re_pt;
269 }
270
271 /// Alpha
272 const double& alpha() const
273 {
274 return *Alpha_pt;
275 }
276
277 /// Product of Reynolds and Strouhal number (=Womersley number)
278 const double& re_st() const
279 {
280 return *ReSt_pt;
281 }
282
283 /// Pointer to Reynolds number
284 double*& re_pt()
285 {
286 return Re_pt;
287 }
288
289 /// Pointer to Alpha
290 double*& alpha_pt()
291 {
292 return Alpha_pt;
293 }
294
295 /// Pointer to product of Reynolds and Strouhal number (=Womersley number)
296 double*& re_st_pt()
297 {
298 return ReSt_pt;
299 }
300
301 /// Viscosity ratio for element: Element's viscosity relative to the
302 /// viscosity used in the definition of the Reynolds number
303 const double& viscosity_ratio() const
304 {
305 return *Viscosity_Ratio_pt;
306 }
307
308 /// Pointer to Viscosity Ratio
310 {
311 return Viscosity_Ratio_pt;
312 }
313
314 /// Density ratio for element: Element's density relative to the
315 /// viscosity used in the definition of the Reynolds number
316 const double& density_ratio() const
317 {
318 return *Density_Ratio_pt;
319 }
320
321 /// Pointer to Density ratio
323 {
324 return Density_Ratio_pt;
325 }
326
327 /// Global inverse Froude number
328 const double& re_invfr() const
329 {
330 return *ReInvFr_pt;
331 }
332
333 /// Pointer to global inverse Froude number
334 double*& re_invfr_pt()
335 {
336 return ReInvFr_pt;
337 }
338
339 /// Vector of gravitational components
340 const Vector<double>& g() const
341 {
342 return *G_pt;
343 }
344
345 /// Pointer to Vector of gravitational components
347 {
348 return G_pt;
349 }
350
351 /// Access function for the body-force pointer
353 {
354 return Body_force_fct_pt;
355 }
356
357 /// Access function for the body-force pointer. Const version
359 {
360 return Body_force_fct_pt;
361 }
362
363 /// Access function for the source-function pointer
365 {
366 return Source_fct_pt;
367 }
368
369 /// Access function for the source-function pointer. Const version
371 {
372 return Source_fct_pt;
373 }
374
375 /// Function to return number of pressure degrees of freedom
376 virtual unsigned npres_pnst() const = 0;
377
378 /// Velocity i at local node n. Uses suitably interpolated value
379 /// for hanging nodes.
380 virtual double u_pnst(const unsigned& n, const unsigned& i) const = 0;
381
382 /// Velocity i at local node n at timestep t (t=0: present;
383 /// t>0: previous). Uses suitably interpolated value for hanging nodes.
384 virtual double u_pnst(const unsigned& t,
385 const unsigned& n,
386 const unsigned& i) const = 0;
387
388 /// Return the index at which the i-th unknown velocity component
389 /// is stored. The default value, i, is appropriate for single-physics
390 /// problems.
391 /// In derived multi-physics elements, this function should be overloaded
392 /// to reflect the chosen storage scheme. Note that these equations require
393 /// that the unknowns are always stored at the same indices at each node.
394 virtual inline unsigned u_index_pnst(const unsigned& i) const
395 {
396 return i;
397 }
398
399 /// i-th component of du/dt at local node n.
400 /// Uses suitably interpolated value for hanging nodes.
401 double du_dt_pnst(const unsigned& n, const unsigned& i) const
402 {
403 // Get the data's timestepper
405
406 // Number of timsteps (past & present)
407 unsigned n_time = time_stepper_pt->ntstorage();
408
409 // Initialise dudt
410 double dudt = 0.0;
411 // Loop over the timesteps, if there is a non Steady timestepper
412 if (time_stepper_pt->type() != "Steady")
413 {
414 for (unsigned t = 0; t < n_time; t++)
415 {
416 dudt += time_stepper_pt->weight(1, t) * u_pnst(t, n, i);
417 }
418 }
419
420 return dudt;
421 }
422
423 /// Pressure at local pressure "node" n_p
424 /// Uses suitably interpolated value for hanging nodes.
425 virtual double p_pnst(const unsigned& n_p) const = 0;
426
427 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
428 virtual void fix_pressure(const unsigned& p_dof, const double& p_value) = 0;
429
430 /// Which nodal value represents the pressure? (Default: negative,
431 /// indicating that pressure is not based on nodal interpolation).
432 virtual int p_nodal_index_pnst()
433 {
435 }
436
437 /// Pointer to n_p-th pressure node (Default: NULL,
438 /// indicating that pressure is not based on nodal interpolation).
439 virtual Node* pressure_node_pt(const unsigned& n_p)
440 {
441 return NULL;
442 }
443
444 /// Integral of pressure over element
445 double pressure_integral() const;
446
447 /// Return integral of dissipation over element
448 double dissipation() const;
449
450 /// Return dissipation at local coordinate s
451 double dissipation(const Vector<double>& s) const;
452
453 /// Get integral of kinetic energy over element
454 double kin_energy() const;
455
456 /// Strain-rate tensor: Now returns polar strain
457 void strain_rate(const Vector<double>& s,
459
460 /// Function to return polar strain multiplied by r
463
464 /// Compute traction (on the viscous scale) at local coordinate s
465 /// for outer unit normal N
466 void get_traction(const Vector<double>& s,
467 const Vector<double>& N,
468 Vector<double>& traction);
469
470 /// The potential loading on an external SolidElement is always
471 /// provided by the traction function
473 const Vector<double>& xi,
474 const Vector<double>& x,
475 const Vector<double>& N,
476 Vector<double>& load)
477 {
478 get_traction(s, N, load);
479 }
480
481
482 /// Output functionget_vels(const Vector<double>& x_to_get,
483 /// Vector<double>& vels): x,y,[z],u,v,[w],p in tecplot format. Default
484 /// number of plot points
485 void output(std::ostream& outfile)
486 {
487 unsigned nplot = 5;
488 output(outfile, nplot);
489 }
490
491 /// Output function: x,y,[z],u,v,[w],p
492 /// in tecplot format. nplot points in each coordinate direction
493 void output(std::ostream& outfile, const unsigned& nplot);
494
495 /// C-style output function: x,y,[z],u,v,[w],p
496 /// in tecplot format. Default number of plot points
497 void output(FILE* file_pt)
498 {
499 unsigned nplot = 5;
500 output(file_pt, nplot);
501 }
502
503 /// C-style output function: x,y,[z],u,v,[w],p
504 /// in tecplot format. nplot points in each coordinate direction
505 void output(FILE* file_pt, const unsigned& nplot);
506
507 /// Full output function:
508 /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
509 /// in tecplot format. Default number of plot points
510 void full_output(std::ostream& outfile)
511 {
512 unsigned nplot = 5;
513 full_output(outfile, nplot);
514 }
515
516 /// Full output function:
517 /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
518 /// in tecplot format. nplot points in each coordinate direction
519 void full_output(std::ostream& outfile, const unsigned& nplot);
520
521
522 /// Output function: x,y,[z],u,v,[w] in tecplot format.
523 /// nplot points in each coordinate direction at timestep t
524 /// (t=0: present; t>0: previous timestep)
525 void output_veloc(std::ostream& outfile,
526 const unsigned& nplot,
527 const unsigned& t);
528
529 /// Output exact solution specified via function pointer
530 /// at a given number of plot points. Function prints as
531 /// many components as are returned in solution Vector
532 void output_fct(std::ostream& outfile,
533 const unsigned& nplot,
535
536 /// Output exact solution specified via function pointer
537 /// at a given time and at a given number of plot points.
538 /// Function prints as many components as are returned in solution Vector.
539 void output_fct(std::ostream& outfile,
540 const unsigned& nplot,
541 const double& time,
543
544 /// Validate against exact solution at given time
545 /// Solution is provided via function pointer.
546 /// Plot at a given number of plot points and compute L2 error
547 /// and L2 norm of velocity solution over element
548 void compute_error(std::ostream& outfile,
550 const double& time,
551 double& error,
552 double& norm);
553
554 /// Validate against exact solution.
555 /// Solution is provided via function pointer.
556 /// Plot at a given number of plot points and compute L2 error
557 /// and L2 norm of velocity solution over element
558 void compute_error(std::ostream& outfile,
560 double& error,
561 double& norm);
562
563 /// Compute the element's residual Vector
565 {
566 // Call the generic residuals function with flag set to 0
570 0);
571 }
572
573 /// Compute the element's residual Vector and the jacobian matrix
574 /// Virtual function can be overloaded by hanging-node version
576 DenseMatrix<double>& jacobian)
577 {
578 // Call the generic routine with the flag set to 1
580 residuals, jacobian, GeneralisedElement::Dummy_matrix, 1);
581
582 /*
583 // Only needed for test purposes
584 //Create a new matrix
585 unsigned n_dof = ndof();
586 DenseMatrix<double> jacobian_fd(n_dof,n_dof,0.0);
587 fill_in_jacobian_from_internal_by_fd(residuals,jacobian_fd);
588 fill_in_jacobian_from_nodal_by_fd(residuals,jacobian_fd);
589
590 if(n_dof==21)
591 {
592 for(unsigned i=0;i<n_dof;i++)
593 {
594 for(unsigned j=0;j<n_dof;j++)
595 {
596 if(std::fabs(jacobian(i,j) - jacobian_fd(i,j)) > 1.0e-3)
597 {
598 std::cout << i << " " << j << " " << std::fabs(jacobian(i,j) -
599 jacobian_fd(i,j)) << std::endl;
600 }
601 }
602 }
603 }
604 // Only needed for test purposes
605
606 // Only needed for test purposes
607 //Create a new matrix
608 unsigned n_dof = ndof();
609 DenseMatrix<double>
610 jacobian_new(n_dof,n_dof,0.0),jacobian_old(n_dof,n_dof,0.0);
611 Vector<double> residuals_new(n_dof,0.0),residuals_old(n_dof,0.0);
612
613 //Call the generic routine with the flag set to 1
614 PolarNavierStokesEquations::fill_in_generic_residual_contribution(residuals_old,jacobian_old,
615 GeneralisedElement::Dummy_matrix,1);
616 //Call the generic routine with the flag set to 1
617 fill_in_generic_residual_contribution(residuals_new,jacobian_new,
618 GeneralisedElement::Dummy_matrix,1);
619
620 if(n_dof==21)
621 {
622 for(unsigned i=0;i<n_dof;i++)
623 {
624 if(std::fabs(residuals_new[i] - residuals_old[i])>1.0e-8) std::cout << i
625 << " " << std::fabs(residuals_new[i] - residuals_old[i]) << std::endl;
626 //std::cout << residuals_new[i] << std::endl;
627 for(unsigned j=0;j<n_dof;j++)
628 {
629 if(std::fabs(jacobian_new(i,j) - jacobian_old(i,j)) > 1.0e-8)
630 {
631 std::cout << i << " " << j << " " << std::fabs(jacobian_new(i,j) -
632 jacobian_old(i,j)) << std::endl;
633 }
634 }
635 }
636 }
637 // Only needed for test purposes
638 */
639 } // End of fill_in_contribution_to_jacobian
640
641 /// Compute the element's residual Vector and the jacobian matrix
642 /// Plus the mass matrix especially for eigenvalue problems
644 Vector<double>& residuals,
645 DenseMatrix<double>& jacobian,
646 DenseMatrix<double>& mass_matrix)
647 {
648 // Call the generic routine with the flag set to 2
650 residuals, jacobian, mass_matrix, 2);
651 }
652
653 /// Compute the residuals for the Navier--Stokes equations;
654 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
656 Vector<double>& residuals,
657 DenseMatrix<double>& jacobian,
658 DenseMatrix<double>& mass_matrix,
659 unsigned flag);
660
661 /// Compute vector of FE interpolated velocity u at local coordinate s
663 Vector<double>& veloc) const
664 {
665 // Find number of nodes
666 unsigned n_node = nnode();
667 // Local shape function
668 Shape psi(n_node);
669 // Find values of shape function
670 shape(s, psi);
671
672 for (unsigned i = 0; i < 2; i++)
673 {
674 // Initialise value of u
675 veloc[i] = 0.0;
676 // Loop over the local nodes and sum
677 for (unsigned l = 0; l < n_node; l++)
678 {
679 veloc[i] += u_pnst(l, i) * psi[l];
680 }
681 }
682 }
683
684 /// Return FE interpolated velocity u[i] at local coordinate s
685 double interpolated_u_pnst(const Vector<double>& s, const unsigned& i) const
686 {
687 // Find number of nodes
688 unsigned n_node = nnode();
689 // Local shape function
690 Shape psi(n_node);
691 // Find values of shape function
692 shape(s, psi);
693
694 // Initialise value of u
695 double interpolated_u = 0.0;
696 // Loop over the local nodes and sum
697 for (unsigned l = 0; l < n_node; l++)
698 {
699 interpolated_u += u_pnst(l, i) * psi[l];
700 }
701
702 return (interpolated_u);
703 }
704
705 /// Return FE interpolated pressure at local coordinate s
707 {
708 // Find number of nodes
709 unsigned n_pres = npres_pnst();
710 // Local shape function
711 Shape psi(n_pres);
712 // Find values of shape function
713 pshape_pnst(s, psi);
714
715 // Initialise value of p
716 double interpolated_p = 0.0;
717 // Loop over the local nodes and sum
718 for (unsigned l = 0; l < n_pres; l++)
719 {
720 interpolated_p += p_pnst(l) * psi[l];
721 }
722
723 return (interpolated_p);
724 }
725
726 /// ////////////////////////////////////////////////////////////////
727 /// My own work: ///
728 /// Return FE interpolated velocity derivative du[i]/dx[j] ///
729 /// at local coordinate s ///
730 /// ////////////////////////////////////////////////////////////////
732 const unsigned& i,
733 const unsigned& j) const
734 {
735 // Find number of nodes
736 unsigned n_node = nnode();
737
738 // Set up memory for the shape and test functions
739 Shape psif(n_node);
740 DShape dpsifdx(n_node, 2);
741
742 // double J =
743 dshape_eulerian(s, psif, dpsifdx);
744
745 // Initialise to zero
746 double interpolated_dudx = 0.0;
747
748 // Calculate velocity derivative:
749
750 // Loop over nodes
751 for (unsigned l = 0; l < n_node; l++)
752 {
753 interpolated_dudx += u_pnst(l, i) * dpsifdx(l, j);
754 }
755
756 return (interpolated_dudx);
757 }
758 };
759
760 /// ///////////////////////////////////////////////////////////////////////////
761 /// ///////////////////////////////////////////////////////////////////////////
762 /// ///////////////////////////////////////////////////////////////////////////
763
764
765 //==========================================================================
766 /// Crouzeix_Raviart elements are Navier--Stokes elements with quadratic
767 /// interpolation for velocities and positions, but a discontinuous linear
768 /// pressure interpolation
769 //==========================================================================
770 class PolarCrouzeixRaviartElement : public virtual QElement<2, 3>,
771 public virtual PolarNavierStokesEquations
772 {
773 private:
774 /// Static array of ints to hold required number of variables at nodes
775 static const unsigned Initial_Nvalue[];
776
777 protected:
778 /// Internal index that indicates at which internal data the pressure
779 /// is stored
781
782 /// Velocity shape and test functions and their derivs
783 /// w.r.t. to global coords at local coordinate s (taken from geometry)
784 /// Return Jacobian of mapping between local and global coordinates.
786 Shape& psi,
787 DShape& dpsidx,
788 Shape& test,
789 DShape& dtestdx) const;
790
791 /// Velocity shape and test functions and their derivs
792 /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
793 /// Return Jacobian of mapping between local and global coordinates.
794 inline double dshape_and_dtest_eulerian_at_knot_pnst(const unsigned& ipt,
795 Shape& psi,
796 DShape& dpsidx,
797 Shape& test,
798 DShape& dtestdx) const;
799
800 /// Pressure shape functions at local coordinate s
801 inline void pshape_pnst(const Vector<double>& s, Shape& psi) const;
802
803 /// Pressure shape and test functions at local coordinte s
804 inline void pshape_pnst(const Vector<double>& s,
805 Shape& psi,
806 Shape& test) const;
807
808 /// Return the local equation numbers for the pressure values.
809 inline int p_local_eqn(const unsigned& n)
810 {
811 return this->internal_local_eqn(P_pnst_internal_index, n);
812 }
813
814 public:
815 /// Constructor, there are DIM+1 internal values (for the pressure)
818 {
819 // Allocate and add one Internal data object that stores 3 pressure
820 // values;
822 }
823
824 /// Number of values (pinned or dofs) required at local node n.
825 virtual unsigned required_nvalue(const unsigned& n) const;
826
827 /// Return the velocity component u[i] at local node
828 /// n. Uses suitably interpolated value for hanging nodes.
829 double u_pnst(const unsigned& n, const unsigned& i) const
830 {
831 return this->nodal_value(n, i);
832 }
833
834 /// Return the velocity component u[i] at local node
835 /// n at timestep t (t=0: present; t>0: previous timestep).
836 /// Uses suitably interpolated value for hanging nodes.
837 double u_pnst(const unsigned& t, const unsigned& n, const unsigned& i) const
838 {
839 return this->nodal_value(t, n, i);
840 }
841
842 /// Return the pressure values at internal dof i_internal
843 /// (Discontinous pressure interpolation -- no need to cater for hanging
844 /// nodes).
845 double p_pnst(const unsigned& i_internal) const
846 {
847 return *this->internal_data_pt(P_pnst_internal_index)
848 ->value_pt(i_internal);
849 }
850
851 /// Return number of pressure values
852 unsigned npres_pnst() const
853 {
854 return 3;
855 }
856
857 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
858 void fix_pressure(const unsigned& p_dof, const double& p_value)
859 {
860 this->internal_data_pt(P_pnst_internal_index)->pin(p_dof);
861 *this->internal_data_pt(P_pnst_internal_index)->value_pt(p_dof) = p_value;
862 }
863
864 /// Add to the set paired_load_data
865 /// pairs of pointers to data objects and unsigned integers that
866 /// index the values in the data object that affect the load (traction),
867 /// as specified in the get_load() function.
868 void get_load_data(std::set<std::pair<Data*, unsigned>>& paired_load_data);
869
870 /// Redirect output to NavierStokesEquations output
871 void output(std::ostream& outfile)
872 {
874 }
875
876 /// Redirect output to NavierStokesEquations output
877 void output(std::ostream& outfile, const unsigned& Nplot)
878 {
880 }
881
882
883 /// Redirect output to NavierStokesEquations output
884 void output(FILE* file_pt)
885 {
887 }
888
889 /// Redirect output to NavierStokesEquations output
890 void output(FILE* file_pt, const unsigned& Nplot)
891 {
893 }
894
895
896 /// Full output function:
897 /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
898 /// in tecplot format. Default number of plot points
899 void full_output(std::ostream& outfile)
900 {
902 }
903
904 /// Full output function:
905 /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
906 /// in tecplot format. nplot points in each coordinate direction
907 void full_output(std::ostream& outfile, const unsigned& nplot)
908 {
910 }
911 };
912
913 // Inline functions
914
915 //=======================================================================
916 /// 2D
917 /// Derivatives of the shape functions and test functions w.r.t. to global
918 /// (Eulerian) coordinates. Return Jacobian of mapping between
919 /// local and global coordinates.
920 //=======================================================================
922 const Vector<double>& s,
923 Shape& psi,
924 DShape& dpsidx,
925 Shape& test,
926 DShape& dtestdx) const
927 {
928 // Call the geometrical shape functions and derivatives
929 double J = QElement<2, 3>::dshape_eulerian(s, psi, dpsidx);
930 // Loop over the test functions and derivatives and set them equal to the
931 // shape functions
932 for (unsigned i = 0; i < 9; i++)
933 {
934 test[i] = psi[i];
935 dtestdx(i, 0) = dpsidx(i, 0);
936 dtestdx(i, 1) = dpsidx(i, 1);
937 }
938 // Return the jacobian
939 return J;
940 }
941
942
943 //=======================================================================
944 /// 2D
945 /// Derivatives of the shape functions and test functions w.r.t. to global
946 /// (Eulerian) coordinates. Return Jacobian of mapping between
947 /// local and global coordinates.
948 //=======================================================================
951 Shape& psi,
952 DShape& dpsidx,
953 Shape& test,
954 DShape& dtestdx) const
955 {
956 // Call the geometrical shape functions and derivatives
957 double J = QElement<2, 3>::dshape_eulerian_at_knot(ipt, psi, dpsidx);
958 // Loop over the test functions and derivatives and set them equal to the
959 // shape functions
960 for (unsigned i = 0; i < 9; i++)
961 {
962 test[i] = psi[i];
963 dtestdx(i, 0) = dpsidx(i, 0);
964 dtestdx(i, 1) = dpsidx(i, 1);
965 }
966 // Return the jacobian
967 return J;
968 }
969
970 //=======================================================================
971 /// 2D :
972 /// Pressure shape functions
973 //=======================================================================
975 Shape& psi) const
976 {
977 psi[0] = 1.0;
978 psi[1] = s[0];
979 psi[2] = s[1];
980 }
981
982 /// Define the pressure shape and test functions
984 Shape& psi,
985 Shape& test) const
986 {
987 // Call the pressure shape functions
988 pshape_pnst(s, psi);
989 // Loop over the test functions and set them equal to the shape functions
990 for (unsigned i = 0; i < 3; i++) test[i] = psi[i];
991 }
992
993 //=======================================================================
994 /// Face geometry of the 2D Crouzeix_Raviart elements
995 //=======================================================================
996 template<>
998 : public virtual QElement<1, 3>
999 {
1000 public:
1001 FaceGeometry() : QElement<1, 3>() {}
1002 };
1003
1004 /// /////////////////////////////////////////////////////////////////////////
1005 /// /////////////////////////////////////////////////////////////////////////
1006
1007 //=======================================================================
1008 /// Taylor--Hood elements are Navier--Stokes elements
1009 /// with quadratic interpolation for velocities and positions and
1010 /// continous linear pressure interpolation
1011 //=======================================================================
1012 class PolarTaylorHoodElement : public virtual QElement<2, 3>,
1013 public virtual PolarNavierStokesEquations
1014 {
1015 private:
1016 /// Static array of ints to hold number of variables at node
1017 static const unsigned Initial_Nvalue[];
1018
1019 protected:
1020 /// Static array of ints to hold conversion from pressure
1021 /// node numbers to actual node numbers
1022 static const unsigned Pconv[];
1023
1024 /// Velocity shape and test functions and their derivs
1025 /// w.r.t. to global coords at local coordinate s (taken from geometry)
1026 /// Return Jacobian of mapping between local and global coordinates.
1027 inline double dshape_and_dtest_eulerian_pnst(const Vector<double>& s,
1028 Shape& psi,
1029 DShape& dpsidx,
1030 Shape& test,
1031 DShape& dtestdx) const;
1032
1033 /// Velocity shape and test functions and their derivs
1034 /// w.r.t. to global coords at local coordinate s (taken from geometry)
1035 /// Return Jacobian of mapping between local and global coordinates.
1036 inline double dshape_and_dtest_eulerian_at_knot_pnst(const unsigned& ipt,
1037 Shape& psi,
1038 DShape& dpsidx,
1039 Shape& test,
1040 DShape& dtestdx) const;
1041
1042 /// Pressure shape functions at local coordinate s
1043 inline void pshape_pnst(const Vector<double>& s, Shape& psi) const;
1044
1045 /// Pressure shape and test functions at local coordinte s
1046 inline void pshape_pnst(const Vector<double>& s,
1047 Shape& psi,
1048 Shape& test) const;
1049
1050 /// Return the local equation numbers for the pressure values.
1051 inline int p_local_eqn(const unsigned& n)
1052 {
1053 return this->nodal_local_eqn(Pconv[n], p_nodal_index_pnst());
1054 }
1055
1056 public:
1057 /// Constructor, no internal data points
1059
1060 /// Number of values (pinned or dofs) required at node n. Can
1061 /// be overwritten for hanging node version
1062 inline virtual unsigned required_nvalue(const unsigned& n) const
1063 {
1064 return Initial_Nvalue[n];
1065 }
1066
1067 /// Which nodal value represents the pressure?
1069 {
1070 return 2;
1071 }
1072
1073 /// Pointer to n_p-th pressure node
1074 Node* pressure_node_pt(const unsigned& n_p)
1075 {
1076 return this->node_pt(Pconv[n_p]);
1077 }
1078
1079 /// Return the velocity component u[i] at local node
1080 /// n. Uses suitably interpolated value for hanging nodes.
1081 double u_pnst(const unsigned& n, const unsigned& i) const
1082 {
1083 return this->nodal_value(n, i);
1084 }
1085
1086 /// Return the velocity component u[i] at local node
1087 /// n at timestep t (t=0: present; t>0: previous timestep).
1088 /// Uses suitably interpolated value for hanging nodes.
1089 double u_pnst(const unsigned& t, const unsigned& n, const unsigned& i) const
1090 {
1091 return this->nodal_value(t, n, i);
1092 }
1093
1094 /// Access function for the pressure values at local pressure
1095 /// node n_p (const version)
1096 double p_pnst(const unsigned& n_p) const
1097 {
1098 return this->nodal_value(Pconv[n_p], 2);
1099 }
1100 // {return this->nodal_value(Pconv[n_p],this->p_nodal_index_pnst());}
1101
1102 /// Return number of pressure values
1103 unsigned npres_pnst() const
1104 {
1105 return 4;
1106 }
1107
1108 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
1109 void fix_pressure(const unsigned& p_dof, const double& p_value)
1110 {
1111 this->node_pt(Pconv[p_dof])->pin(this->p_nodal_index_pnst());
1112 *this->node_pt(Pconv[p_dof])->value_pt(this->p_nodal_index_pnst()) =
1113 p_value;
1114 }
1115
1116
1117 /// Add to the set paired_load_data
1118 /// pairs of pointers to data objects and unsigned integers that
1119 /// index the values in the data object that affect the load (traction),
1120 /// as specified in the get_load() function.
1121 void get_load_data(std::set<std::pair<Data*, unsigned>>& paired_load_data);
1122
1123 /// Redirect output to NavierStokesEquations output
1124 void output(std::ostream& outfile)
1125 {
1127 }
1128
1129 /// Redirect output to NavierStokesEquations output
1130 void output(std::ostream& outfile, const unsigned& Nplot)
1131 {
1132 PolarNavierStokesEquations::output(outfile, Nplot);
1133 }
1134
1135 /// Redirect output to NavierStokesEquations output
1136 void output(FILE* file_pt)
1137 {
1139 }
1140
1141 /// Redirect output to NavierStokesEquations output
1142 void output(FILE* file_pt, const unsigned& Nplot)
1143 {
1144 PolarNavierStokesEquations::output(file_pt, Nplot);
1145 }
1146 };
1147
1148 // Inline functions
1149
1150 //==========================================================================
1151 /// 2D :
1152 /// Derivatives of the shape functions and test functions w.r.t to
1153 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1154 /// local and global coordinates.
1155 //==========================================================================
1157 const Vector<double>& s,
1158 Shape& psi,
1159 DShape& dpsidx,
1160 Shape& test,
1161 DShape& dtestdx) const
1162 {
1163 // Call the geometrical shape functions and derivatives
1164 double J = QElement<2, 3>::dshape_eulerian(s, psi, dpsidx);
1165 // Loop over the test functions and derivatives and set them equal to the
1166 // shape functions
1167 for (unsigned i = 0; i < 9; i++)
1168 {
1169 test[i] = psi[i];
1170 dtestdx(i, 0) = dpsidx(i, 0);
1171 dtestdx(i, 1) = dpsidx(i, 1);
1172 }
1173 // Return the jacobian
1174 return J;
1175 }
1176
1177
1178 //==========================================================================
1179 /// 2D :
1180 /// Derivatives of the shape functions and test functions w.r.t to
1181 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1182 /// local and global coordinates.
1183 //==========================================================================
1185 const unsigned& ipt,
1186 Shape& psi,
1187 DShape& dpsidx,
1188 Shape& test,
1189 DShape& dtestdx) const
1190 {
1191 // Call the geometrical shape functions and derivatives
1192 double J = QElement<2, 3>::dshape_eulerian_at_knot(ipt, psi, dpsidx);
1193 // Loop over the test functions and derivatives and set them equal to the
1194 // shape functions
1195 for (unsigned i = 0; i < 9; i++)
1196 {
1197 test[i] = psi[i];
1198 dtestdx(i, 0) = dpsidx(i, 0);
1199 dtestdx(i, 1) = dpsidx(i, 1);
1200 }
1201 // Return the jacobian
1202 return J;
1203 }
1204
1205
1206 //==========================================================================
1207 /// 2D :
1208 /// Pressure shape functions
1209 //==========================================================================
1211 Shape& psi) const
1212 {
1213 // Call the OneDimensional Shape functions
1214 // Local storage
1215 double psi1[2], psi2[2];
1216 // Call the OneDimensional Shape functions
1217 OneDimLagrange::shape<2>(s[0], psi1);
1218 OneDimLagrange::shape<2>(s[1], psi2);
1219
1220 // Now let's loop over the nodal points in the element
1221 // s1 is the "x" coordinate, s2 the "y"
1222 for (unsigned i = 0; i < 2; i++)
1223 {
1224 for (unsigned j = 0; j < 2; j++)
1225 {
1226 /*Multiply the two 1D functions together to get the 2D function*/
1227 psi[2 * i + j] = psi2[i] * psi1[j];
1228 }
1229 }
1230 }
1231
1232
1233 //==========================================================================
1234 /// 2D :
1235 /// Pressure shape and test functions
1236 //==========================================================================
1238 Shape& psi,
1239 Shape& test) const
1240 {
1241 // Call the pressure shape functions
1242 pshape_pnst(s, psi);
1243 // Loop over the test functions and set them equal to the shape functions
1244 for (unsigned i = 0; i < 4; i++) test[i] = psi[i];
1245 }
1246
1247 //=======================================================================
1248 /// Face geometry of the 2D Taylor_Hood elements
1249 //=======================================================================
1250 template<>
1251 class FaceGeometry<PolarTaylorHoodElement> : public virtual QElement<1, 3>
1252 {
1253 public:
1254 FaceGeometry() : QElement<1, 3>() {}
1255 };
1256
1257} // namespace oomph
1258#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
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
Definition: nodes.h:324
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: 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 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
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
double u_pnst(const unsigned &t, const unsigned &n, const unsigned &i) const
Return the velocity component u[i] at local node n at timestep t (t=0: present; t>0: previous timeste...
void output(FILE *file_pt, const unsigned &Nplot)
Redirect output to NavierStokesEquations 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.
int p_local_eqn(const unsigned &n)
Return the local equation numbers for the pressure values.
PolarCrouzeixRaviartElement()
Constructor, there are DIM+1 internal values (for the pressure)
void pshape_pnst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
unsigned npres_pnst() const
Return number of pressure values.
virtual unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
void get_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs of pointers to data objects and unsigned integers that index th...
unsigned P_pnst_internal_index
Internal index that indicates at which internal data the pressure is stored.
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....
double u_pnst(const unsigned &n, const unsigned &i) const
Return the velocity component u[i] at local node n. Uses suitably interpolated value for hanging node...
double dshape_and_dtest_eulerian_at_knot_pnst(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...
double p_pnst(const unsigned &i_internal) const
Return the pressure values at internal dof i_internal (Discontinous pressure interpolation – no need ...
void output(std::ostream &outfile, const unsigned &Nplot)
Redirect output to NavierStokesEquations output.
double dshape_and_dtest_eulerian_pnst(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(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
static const unsigned Initial_Nvalue[]
Static array of ints to hold required number of variables at nodes.
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....
A class for elements that solve the polar Navier–Stokes equations, This contains the generic maths – ...
PolarNavierStokesEquations()
Constructor: NULL the body force and source function.
double dissipation() const
Return integral of dissipation over element.
double *& re_invfr_pt()
Pointer to global inverse Froude number.
double interpolated_u_pnst(const Vector< double > &s, const unsigned &i) const
Return FE interpolated velocity u[i] at local coordinate s.
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute the element's residual Vector.
virtual void pshape_pnst(const Vector< double > &s, Shape &psi) const =0
Compute the pressure shape functions at local coordinate s.
const double & re() const
Reynolds number.
double kin_energy() const
Get integral of kinetic energy over element.
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Compute the element's residual Vector and the jacobian matrix Virtual function can be overloaded by h...
double *& re_pt()
Pointer to Reynolds number.
virtual void pshape_pnst(const Vector< double > &s, Shape &psi, Shape &test) const =0
Compute the pressure shape and test functions at local coordinate s.
const double & re_st() const
Product of Reynolds and Strouhal number (=Womersley number)
Vector< double > * G_pt
Pointer to global gravity Vector.
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not.
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
const double & density_ratio() const
Density ratio for element: Element's density relative to the viscosity used in the definition of the ...
virtual int p_local_eqn(const unsigned &n)=0
Access function for the local equation number information for the pressure. p_local_eqn[n] = local eq...
void get_load(const Vector< double > &s, const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
The potential loading on an external SolidElement is always provided by the traction function.
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
static double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
virtual double p_pnst(const unsigned &n_p) const =0
Pressure at local pressure "node" n_p Uses suitably interpolated value for hanging nodes.
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)
Compute traction (on the viscous scale) at local coordinate s for outer unit normal N.
const double & viscosity_ratio() const
Viscosity ratio for element: Element's viscosity relative to the viscosity used in the definition of ...
virtual double dshape_and_dtest_eulerian_pnst(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....
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: Now returns polar strain.
double *& density_ratio_pt()
Pointer to Density ratio.
double get_source_fct(double time, const Vector< double > &x)
Calculate the source fct at given time and Eulerian position.
NavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
void output(std::ostream &outfile)
Output functionget_vels(const Vector<double>& x_to_get, Vector<double>& vels): x,y,...
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....
const double & re_invfr() const
Global inverse Froude number.
const Vector< double > & g() const
Vector of gravitational components.
virtual int p_nodal_index_pnst()
Which nodal value represents the pressure? (Default: negative, indicating that pressure is not based ...
void get_body_force(double time, const Vector< double > &x, Vector< double > &result)
Calculate the body force at a given time and Eulerian position.
virtual double dshape_and_dtest_eulerian_at_knot_pnst(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...
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 * Re_pt
Pointer to global Reynolds number.
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
double(* NavierStokesSourceFctPt)(const double &time, const Vector< double > &x)
Function pointer to source function fct(t,x) x is a Vector!
NavierStokesSourceFctPt & source_fct_pt()
Access function for the source-function pointer.
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 strain_rate_by_r(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Function to return polar strain multiplied by r.
virtual unsigned u_index_pnst(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)
void output(FILE *file_pt)
C-style output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary 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....
void interpolated_u_pnst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
double du_dt_pnst(const unsigned &n, const unsigned &i) const
i-th component of du/dt at local node n. Uses suitably interpolated value for hanging nodes.
virtual double u_pnst(const unsigned &t, const unsigned &n, const unsigned &i) const =0
Velocity i at local node n at timestep t (t=0: present; t>0: previous). Uses suitably interpolated va...
virtual void fill_in_generic_residual_contribution(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...
double pressure_integral() const
Integral of pressure over element.
virtual Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interp...
void(* NavierStokesBodyForceFctPt)(const double &time, const Vector< double > &x, Vector< double > &body_force)
Function pointer to body force function fct(t,x,f(x)) x is a Vector!
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Compute the element's residual Vector and the jacobian matrix Plus the mass matrix especially for eig...
virtual double u_pnst(const unsigned &n, const unsigned &i) const =0
Velocity i at local node n. Uses suitably interpolated value for hanging nodes.
NavierStokesSourceFctPt source_fct_pt() const
Access function for the source-function pointer. Const version.
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....
double * Alpha_pt
Pointer to the angle alpha.
virtual unsigned npres_pnst() const =0
Function to return number of pressure degrees of freedom.
NavierStokesBodyForceFctPt body_force_fct_pt() const
Access function for the body-force pointer. Const version.
double interpolated_dudx_pnst(const Vector< double > &s, const unsigned &i, const unsigned &j) const
//////////////////////////////////////////////////////////////// My own work: /// Return FE interpola...
double interpolated_p_pnst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
/////////////////////////////////////////////////////////////////////////
void pshape_pnst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
void get_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs of pointers to data objects and unsigned integers that index th...
double p_pnst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version)
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.
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 &Nplot)
Redirect output to NavierStokesEquations output.
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
unsigned npres_pnst() const
Return number of pressure values.
double dshape_and_dtest_eulerian_pnst(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...
double u_pnst(const unsigned &n, const unsigned &i) const
Return the velocity component u[i] at local node n. Uses suitably interpolated value for hanging node...
Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node.
double dshape_and_dtest_eulerian_at_knot_pnst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
int p_local_eqn(const unsigned &n)
Return the local equation numbers for the pressure values.
void output(std::ostream &outfile, const unsigned &Nplot)
Redirect output to NavierStokesEquations output.
double u_pnst(const unsigned &t, const unsigned &n, const unsigned &i) const
Return the velocity component u[i] at local node n at timestep t (t=0: present; t>0: previous timeste...
PolarTaylorHoodElement()
Constructor, no internal data points.
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations 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.
virtual int p_nodal_index_pnst()
Which nodal value represents the pressure?
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
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
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...