refineable_spherical_navier_stokes_elements.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for refineable axisymmetric quad Navier Stokes elements
27#ifndef OOMPH_REFINEABLE_SPHERICAL_NAVIER_STOKES_ELEMENTS_HEADER
28#define OOMPH_REFINEABLE_SPHERICAL_NAVIER_STOKES_ELEMENTS_HEADER
29
30// Config header generated by autoconfig
31#ifdef HAVE_CONFIG_H
32#include <oomph-lib-config.h>
33#endif
34
35// Oomph lib includes
36#include "../generic/refineable_quad_element.h"
37#include "../generic/error_estimator.h"
39
40namespace oomph
41{
42 //======================================================================
43 /// Refineable version of the Spherical Navier--Stokes equations
44 //======================================================================
46 : public virtual SphericalNavierStokesEquations,
47 public virtual RefineableElement,
48 public virtual ElementWithZ2ErrorEstimator
49 {
50 protected:
51 /// Pointer to n_p-th pressure node (Default: NULL,
52 /// indicating that pressure is not based on nodal interpolation).
53 virtual Node* pressure_node_pt(const unsigned& n_p)
54 {
55 return NULL;
56 }
57
58 /// Unpin all pressure dofs in the element
60
61 /// Pin unused nodal pressure dofs (empty by default, because
62 /// by default pressure dofs are not associated with nodes)
64
65 public:
66 /// Empty Constructor
71 {
72 }
73
74 /// Number of 'flux' terms for Z2 error estimation
76 {
77 // 3 diagonal strain rates, 3 off diagonal
78 return 6;
79 }
80
81 /// Get 'flux' for Z2 error recovery: Upper triangular entries
82 /// in strain rate tensor.
84 {
85 // Specify the number of velocity dimensions
86 unsigned DIM = 3;
87#ifdef PARANOID
88 unsigned num_entries = DIM + ((DIM * DIM) - DIM) / 2;
89 if (flux.size() < num_entries)
90 {
91 std::ostringstream error_message;
92 error_message << "The flux vector is too small, size " << flux.size()
93 << ", whereas it should be " << num_entries << std::endl;
94 throw OomphLibError(error_message.str(),
95 OOMPH_CURRENT_FUNCTION,
96 OOMPH_EXCEPTION_LOCATION);
97 }
98#endif
99 // Get strain rate matrix
100 DenseMatrix<double> strainrate(DIM);
101 this->strain_rate(s, strainrate);
102
103 // Pack into flux Vector
104 unsigned icount = 0;
105
106 // Start with diagonal terms
107 for (unsigned i = 0; i < DIM; i++)
108 {
109 flux[icount] = strainrate(i, i);
110 icount++;
111 }
112
113 // Off diagonals row by row
114 for (unsigned i = 0; i < DIM; i++)
115 {
116 for (unsigned j = i + 1; j < DIM; j++)
117 {
118 flux[icount] = strainrate(i, j);
119 icount++;
120 }
121 }
122 }
123
124 /// Fill in the geometric Jacobian, which in this case is r*r*sin(theta)
126 {
127 return x[0] * x[0] * sin(x[1]);
128 }
129
130
131 /// Further build: pass the pointers down to the sons
133 {
134 // Find the father element
135 RefineableSphericalNavierStokesEquations* cast_father_element_pt =
137 this->father_element_pt());
138
139 // Set the viscosity ratio pointer
140 this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
141 // Set the density ratio pointer
142 this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
143 // Set pointer to global Reynolds number
144 this->Re_pt = cast_father_element_pt->re_pt();
145 // Set pointer to global Reynolds number x Strouhal number (=Womersley)
146 this->ReSt_pt = cast_father_element_pt->re_st_pt();
147 // Set pointer to the global Reynolds number x inverse Rossby number
148 this->ReInvRo_pt = cast_father_element_pt->re_invro_pt();
149 // Set pointer to global Reynolds number x inverse Froude number
150 this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
151 // Set pointer to global gravity Vector
152 this->G_pt = cast_father_element_pt->g_pt();
153
154 // Set pointer to body force function
155 this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
156
157 // Set pointer to volumetric source function
158 this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
159
160 // Set the ALE flag
161 this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
162 }
163
164 /// Loop over all elements in Vector (which typically contains
165 /// all the elements in a fluid mesh) and pin the nodal pressure degrees
166 /// of freedom that are not being used. Function uses
167 /// the member function
168 /// - \c RefineableSphericalNavierStokesEquations::
169 /// pin_all_nodal_pressure_dofs()
170 /// .
171 /// which is empty by default and should be implemented for
172 /// elements with nodal pressure degrees of freedom
173 /// (e.g. for refineable Taylor-Hood.)
175 const Vector<GeneralisedElement*>& element_pt)
176 {
177 // Loop over all elements to brutally pin all nodal pressure degrees of
178 // freedom
179 unsigned n_element = element_pt.size();
180 for (unsigned e = 0; e < n_element; e++)
181 {
182 dynamic_cast<RefineableSphericalNavierStokesEquations*>(element_pt[e])
184 }
185 }
186
187 /// Unpin all pressure dofs in elements listed in vector.
189 const Vector<GeneralisedElement*>& element_pt)
190 {
191 // Loop over all elements to brutally unpin all nodal pressure degrees of
192 // freedom and internal pressure degrees of freedom
193 unsigned n_element = element_pt.size();
194 for (unsigned e = 0; e < n_element; e++)
195 {
196 dynamic_cast<RefineableSphericalNavierStokesEquations*>(element_pt[e])
198 }
199 }
200
201
202 private:
203 /// Add element's contribution to the elemental residual vector
204 /// and/or Jacobian matrix
205 /// flag=1: compute both
206 /// flag=0: compute only residual vector
208 Vector<double>& residuals,
209 DenseMatrix<double>& jacobian,
210 DenseMatrix<double>& mass_matrix,
211 unsigned flag);
212 };
213
214 //======================================================================
215 /// Refineable version of Spherical Quad Taylor Hood elements.
216 /// (note that unlike the cartesian version this is not scale-able
217 /// to higher dimensions!)
218 //======================================================================
222 public virtual RefineableQElement<2>
223 {
224 private:
225 /// Pointer to n_p-th pressure node
226 Node* pressure_node_pt(const unsigned& n_p)
227 {
228 return this->node_pt(this->Pconv[n_p]);
229 }
230
231 /// Unpin all pressure dofs
233 {
234 unsigned n_node = this->nnode();
235 int p_index = this->p_nodal_index_spherical_nst();
236 // loop over nodes
237 for (unsigned n = 0; n < n_node; n++)
238 {
239 this->node_pt(n)->unpin(p_index);
240 }
241 }
242
243 /// Unpin the proper nodal pressure dofs
245 {
246 // Loop over all nodes
247 unsigned n_node = this->nnode();
248 int p_index = this->p_nodal_index_spherical_nst();
249 // loop over all nodes and pin all the nodal pressures
250 for (unsigned n = 0; n < n_node; n++)
251 {
252 this->node_pt(n)->pin(p_index);
253 }
254
255 // Loop over all actual pressure nodes and unpin if they're not hanging
256 unsigned n_pres = this->npres_spherical_nst();
257 for (unsigned l = 0; l < n_pres; l++)
258 {
259 Node* nod_pt = this->node_pt(this->Pconv[l]);
260 if (!nod_pt->is_hanging(p_index))
261 {
262 nod_pt->unpin(p_index);
263 }
264 }
265 }
266
267 public:
268 /// Constructor:
274 {
275 }
276
277 /// Number of values (pinned or dofs) required at node n.
278 // Bumped up to 4 so we don't have to worry if a hanging mid-side node
279 // gets shared by a corner node (which has extra degrees of freedom)
280 unsigned required_nvalue(const unsigned& n) const
281 {
282 return 4;
283 }
284
285 /// Number of continuously interpolated values: 4 (3 velocities + 1
286 /// pressure)
288 {
289 return 4;
290 }
291
292 /// Rebuild from sons: empty
293 void rebuild_from_sons(Mesh*& mesh_pt) {}
294
295 /// Order of recovery shape functions for Z2 error estimation:
296 /// Same order as shape functions.
298 {
299 return 2;
300 }
301
302 /// Number of vertex nodes in the element
303 unsigned nvertex_node() const
304 {
306 }
307
308 /// Pointer to the j-th vertex node in the element
309 Node* vertex_node_pt(const unsigned& j) const
310 {
312 }
313
314 /// Get the function value u in Vector.
315 /// Note: Given the generality of the interface (this function
316 /// is usually called from black-box documentation or interpolation
317 /// routines), the values Vector sets its own size in here.
319 Vector<double>& values)
320 {
321 // Set size of values Vector: u,w,v,p and initialise to zero
322 values.resize(4, 0.0);
323
324 // Calculate three velocities: values[0],...
325 for (unsigned i = 0; i < 3; i++)
326 {
327 values[i] = interpolated_u_spherical_nst(s, i);
328 }
329
330 // Calculate pressure: values[3]
331 values[3] = interpolated_p_spherical_nst(s);
332 }
333
334 /// Get the function value u in Vector.
335 /// Note: Given the generality of the interface (this function
336 /// is usually called from black-box documentation or interpolation
337 /// routines), the values Vector sets its own size in here.
338 void get_interpolated_values(const unsigned& t,
339 const Vector<double>& s,
340 Vector<double>& values)
341 {
342 // Set size of Vector: u,w,v,p
343 values.resize(4);
344
345 // Initialise
346 for (unsigned i = 0; i < 4; i++)
347 {
348 values[i] = 0.0;
349 }
350
351 // Find out how many nodes there are
352 unsigned n_node = nnode();
353
354 // Shape functions
355 Shape psif(n_node);
356 shape(s, psif);
357
358 // Calculate velocities: values[0],...
359 for (unsigned i = 0; i < 3; i++)
360 {
361 // Get the nodal coordinate of the velocity
362 const unsigned u_nodal_index = u_index_spherical_nst(i);
363 for (unsigned l = 0; l < n_node; l++)
364 {
365 values[i] += nodal_value(t, l, u_nodal_index) * psif[l];
366 }
367 }
368
369 // Calculate pressure: values[3]
370 //(no history is carried in the pressure)
371 values[3] = interpolated_p_spherical_nst(s);
372 }
373
374 /// Perform additional hanging node procedures for variables
375 /// that are not interpolated by all nodes. The pressures are stored
376 /// at the 3rd location in each node
378 {
379 this->setup_hang_for_value(3);
380 }
381
382 /// The velocities are isoparametric and so the "nodes" interpolating
383 /// the velocities are the geometric nodes. The pressure "nodes" are a
384 /// subset of the nodes, so when n_value==3, the n-th pressure
385 /// node is returned.
386 Node* interpolating_node_pt(const unsigned& n, const int& n_value)
387
388 {
389 // The only different nodes are the pressure nodes
390 if (n_value == 3)
391 {
392 return this->node_pt(this->Pconv[n]);
393 }
394 // The other variables are interpolated via the usual nodes
395 else
396 {
397 return this->node_pt(n);
398 }
399 }
400
401 /// The pressure nodes are the corner nodes, so when n_value==DIM,
402 /// the fraction is the same as the 1d node number, 0 or 1.
404 const unsigned& i,
405 const int& n_value)
406 {
407 if (n_value == 3)
408 {
409 // The pressure nodes are just located on the boundaries at 0 or 1
410 return double(n1d);
411 }
412 // Otherwise the velocity nodes are the same as the geometric ones
413 else
414 {
415 return this->local_one_d_fraction_of_node(n1d, i);
416 }
417 }
418
419 /// The velocity nodes are the same as the geometric nodes. The
420 /// pressure nodes must be calculated by using the same methods as
421 /// the geometric nodes, but by recalling that there are only two pressure
422 /// nodes per edge.
424 const int& n_value)
425 {
426 // If we are calculating pressure nodes
427 if (n_value == 3)
428 {
429 // Storage for the index of the pressure node
430 unsigned total_index = 0;
431 // The number of nodes along each 1d edge is 2.
432 unsigned NNODE_1D = 2;
433 // Storage for the index along each boundary
434 // Note that it's only a 2D spatial element
435 Vector<int> index(2);
436 // Loop over the coordinates
437 for (unsigned i = 0; i < 2; i++)
438 {
439 // If we are at the lower limit, the index is zero
440 if (s[i] == -1.0)
441 {
442 index[i] = 0;
443 }
444 // If we are at the upper limit, the index is the number of nodes
445 // minus 1
446 else if (s[i] == 1.0)
447 {
448 index[i] = NNODE_1D - 1;
449 }
450 // Otherwise, we have to calculate the index in general
451 else
452 {
453 // For uniformly spaced nodes the 0th node number would be
454 double float_index = 0.5 * (1.0 + s[i]) * (NNODE_1D - 1);
455 index[i] = int(float_index);
456 // What is the excess. This should be safe because the
457 // taking the integer part rounds down
458 double excess = float_index - index[i];
459 // If the excess is bigger than our tolerance there is no node,
460 // return null
463 {
464 return 0;
465 }
466 }
467 /// Construct the general pressure index from the components.
468 total_index +=
469 index[i] * static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
470 static_cast<int>(i)));
471 }
472 // If we've got here we have a node, so let's return a pointer to it
473 return this->node_pt(this->Pconv[total_index]);
474 }
475 // Otherwise velocity nodes are the same as pressure nodes
476 else
477 {
478 return this->get_node_at_local_coordinate(s);
479 }
480 }
481
482
483 /// The number of 1d pressure nodes is 2, the number of 1d velocity
484 /// nodes is the same as the number of 1d geometric nodes.
485 unsigned ninterpolating_node_1d(const int& n_value)
486 {
487 if (n_value == 3)
488 {
489 return 2;
490 }
491 else
492 {
493 return this->nnode_1d();
494 }
495 }
496
497 /// The number of pressure nodes is 4. The number of
498 /// velocity nodes is the same as the number of geometric nodes.
499 unsigned ninterpolating_node(const int& n_value)
500 {
501 if (n_value == 3)
502 {
503 return 4;
504 }
505 else
506 {
507 return this->nnode();
508 }
509 }
510
511 /// The basis interpolating the pressure is given by pshape().
512 /// / The basis interpolating the velocity is shape().
514 Shape& psi,
515 const int& n_value) const
516 {
517 if (n_value == 3)
518 {
519 return this->pshape_spherical_nst(s, psi);
520 }
521 else
522 {
523 return this->shape(s, psi);
524 }
525 }
526 };
527
528
529 /// ////////////////////////////////////////////////////////////////////////
530 /// ////////////////////////////////////////////////////////////////////////
531 /// ////////////////////////////////////////////////////////////////////////
532
533 //=======================================================================
534 /// Face geometry of the RefineableQuadQTaylorHoodElements
535 //=======================================================================
536 template<>
538 : public virtual FaceGeometry<QSphericalTaylorHoodElement>
539 {
540 public:
542 };
543
544
545 //=======================================================================
546 /// Face geometry of the RefineableQuadQTaylorHoodElements
547 //=======================================================================
548 template<>
550 : public virtual FaceGeometry<FaceGeometry<QSphericalTaylorHoodElement>>
551 {
552 public:
554 {
555 }
556 };
557
558
559 //======================================================================
560 /// Refineable version of Spherical Quad Crouzeix Raviart elements
561 /// (note that unlike the cartesian version this is not scale-able
562 /// to higher dimensions!)
563 //======================================================================
567 public virtual RefineableQElement<2>
568 {
569 private:
570 /// Unpin all the internal pressure freedoms
572 {
573 unsigned n_pres = this->npres_spherical_nst();
574 // loop over pressure dofs and unpin
575 for (unsigned l = 0; l < n_pres; l++)
576 {
578 }
579 }
580
581 public:
582 /// Constructor:
588 {
589 }
590
591 /// Number of continuously interpolated values: 3 (velocities)
593 {
594 return 3;
595 }
596
597 /// Rebuild from sons: Reconstruct pressure from the (merged) sons
598 void rebuild_from_sons(Mesh*& mesh_pt)
599 {
600 using namespace QuadTreeNames;
601
602 // Central pressure value:
603 //-----------------------
604
605 // Use average of the sons central pressure values
606 // Other options: Take average of the four (discontinuous)
607 // pressure values at the father's midpoint]
608
609 double av_press = 0.0;
610
611 // Loop over the sons
612 for (unsigned ison = 0; ison < 4; ison++)
613 {
614 // Add the sons midnode pressure
615 av_press += quadtree_pt()
616 ->son_pt(ison)
617 ->object_pt()
619 ->value(0);
620 }
621
622 // Use the average
624 ->set_value(0, 0.25 * av_press);
625
626
627 // Slope in s_0 direction
628 //----------------------
629
630 // Use average of the 2 FD approximations based on the
631 // elements central pressure values
632 // [Other options: Take average of the four
633 // pressure derivatives]
634
635 double slope1 = quadtree_pt()
636 ->son_pt(SE)
637 ->object_pt()
639 ->value(0) -
641 ->son_pt(SW)
642 ->object_pt()
644 ->value(0);
645
646 double slope2 = quadtree_pt()
647 ->son_pt(NE)
648 ->object_pt()
650 ->value(0) -
652 ->son_pt(NW)
653 ->object_pt()
655 ->value(0);
656
657
658 // Use the average
660 ->set_value(1, 0.5 * (slope1 + slope2));
661
662
663 // Slope in s_1 direction
664 //----------------------
665
666 // Use average of the 2 FD approximations based on the
667 // elements central pressure values
668 // [Other options: Take average of the four
669 // pressure derivatives]
670
671 slope1 = quadtree_pt()
672 ->son_pt(NE)
673 ->object_pt()
675 ->value(0) -
677 ->son_pt(SE)
678 ->object_pt()
680 ->value(0);
681
682 slope2 = quadtree_pt()
683 ->son_pt(NW)
684 ->object_pt()
686 ->value(0) -
688 ->son_pt(SW)
689 ->object_pt()
691 ->value(0);
692
693
694 // Use the average
696 ->set_value(2, 0.5 * (slope1 + slope2));
697 }
698
699 /// Order of recovery shape functions for Z2 error estimation:
700 /// Same order as shape functions.
702 {
703 return 2;
704 }
705
706 /// Number of vertex nodes in the element
707 unsigned nvertex_node() const
708 {
710 }
711
712 /// Pointer to the j-th vertex node in the element
713 Node* vertex_node_pt(const unsigned& j) const
714 {
716 }
717
718 /// Get the function value u in Vector.
719 /// Note: Given the generality of the interface (this function
720 /// is usually called from black-box documentation or interpolation
721 /// routines), the values Vector sets its own size in here.
723 Vector<double>& values)
724 {
725 // Set size of Vector: u,w,v and initialise to zero
726 values.resize(3, 0.0);
727
728 // Calculate velocities: values[0],...
729 for (unsigned i = 0; i < 3; i++)
730 {
731 values[i] = interpolated_u_spherical_nst(s, i);
732 }
733 }
734
735 /// Get all function values [u,v..,p] at previous timestep t
736 /// (t=0: present; t>0: previous timestep).
737 /// \n
738 /// Note: Given the generality of the interface (this function is
739 /// usually called from black-box documentation or interpolation routines),
740 /// the values Vector sets its own size in here.
741 /// \n
742 /// Note: No pressure history is kept, so pressure is always
743 /// the current value.
744 void get_interpolated_values(const unsigned& t,
745 const Vector<double>& s,
746 Vector<double>& values)
747 {
748 // Set size of Vector: u,w,v
749 values.resize(3);
750
751 // Initialise
752 for (unsigned i = 0; i < 3; i++)
753 {
754 values[i] = 0.0;
755 }
756
757 // Find out how many nodes there are
758 unsigned n_node = nnode();
759
760 // Shape functions
761 Shape psif(n_node);
762 shape(s, psif);
763
764 // Calculate velocities: values[0],...
765 for (unsigned i = 0; i < 3; i++)
766 {
767 // Get the local index at which the i-th velocity is stored
768 unsigned u_local_index = u_index_spherical_nst(i);
769 for (unsigned l = 0; l < n_node; l++)
770 {
771 values[i] += nodal_value(t, l, u_local_index) * psif[l];
772 }
773 }
774 }
775
776 /// Perform additional hanging node procedures for variables
777 /// that are not interpolated by all nodes. Empty
779
780 /// Further build for Crouzeix_Raviart interpolates the internal
781 /// pressure dofs from father element: Make sure pressure values and
782 /// dp/ds agree between fathers and sons at the midpoints of the son
783 /// elements.
785 {
786 // Call the generic further build
788
789 using namespace QuadTreeNames;
790
791 // What type of son am I? Ask my quadtree representation...
792 int son_type = quadtree_pt()->son_type();
793
794 // Pointer to my father (in element impersonation)
795 RefineableElement* father_el_pt = quadtree_pt()->father_pt()->object_pt();
796
797 Vector<double> s_father(2);
798
799 // Son midpoint is located at the following coordinates in father element:
800
801 // South west son
802 if (son_type == SW)
803 {
804 s_father[0] = -0.5;
805 s_father[1] = -0.5;
806 }
807 // South east son
808 else if (son_type == SE)
809 {
810 s_father[0] = 0.5;
811 s_father[1] = -0.5;
812 }
813 // North east son
814 else if (son_type == NE)
815 {
816 s_father[0] = 0.5;
817 s_father[1] = 0.5;
818 }
819
820 // North west son
821 else if (son_type == NW)
822 {
823 s_father[0] = -0.5;
824 s_father[1] = 0.5;
825 }
826
827 // Pressure value in father element
829 dynamic_cast<RefineableQSphericalCrouzeixRaviartElement*>(father_el_pt);
830 double press = cast_father_el_pt->interpolated_p_spherical_nst(s_father);
831
832 // Pressure value gets copied straight into internal dof:
834
835
836 // The slopes get copied from father
838 ->set_value(1,
839 0.5 * cast_father_el_pt
841 ->value(1));
842
844 ->set_value(2,
845 0.5 * cast_father_el_pt
847 ->value(2));
848 }
849 };
850
851
852 //=======================================================================
853 /// Face geometry of the RefineableQuadQCrouzeixRaviartElements
854 //=======================================================================
855 template<>
857 : public virtual FaceGeometry<QSphericalCrouzeixRaviartElement>
858 {
859 public:
861 };
862
863
864 //=======================================================================
865 /// Face geometry of the RefineableQuadQCrouzeixRaviartElements
866 //=======================================================================
867 template<>
869 : public virtual FaceGeometry<
870 FaceGeometry<QSphericalCrouzeixRaviartElement>>
871 {
872 public:
875 {
876 }
877 };
878
879
880} // namespace oomph
881
882#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
void unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:391
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
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual Node * get_node_at_local_coordinate(const Vector< double > &s) const
If there is a node at this local coordinate, return the pointer to the node.
Definition: elements.cc:3882
virtual unsigned nvertex_node() const
Return the number of vertex nodes in this element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2491
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
Definition: elements.h:1374
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...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2218
virtual Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2500
virtual double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
Get the local fraction of any node in the n-th position in a one dimensional expansion along the i-th...
Definition: elements.h:1858
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
A general mesh class.
Definition: mesh.h:67
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
unsigned npres_spherical_nst() const
Return number of pressure values.
unsigned P_spherical_nst_internal_index
Internal index that indicates at which internal data the pressure is stored.
/////////////////////////////////////////////////////////////////////////
int p_nodal_index_spherical_nst() const
Set the value at which the pressure is stored in the nodes In this case the third index because there...
void pshape_spherical_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
unsigned npres_spherical_nst() const
Return number of pressure values.
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
void setup_hang_for_value(const int &value_id)
Internal helper function that is used to construct the hanging node schemes for the value_id-th inter...
QuadTree * quadtree_pt()
Pointer to quadtree representation of this element.
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
Refineable version of Spherical Quad Crouzeix Raviart elements (note that unlike the cartesian versio...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get all function values [u,v..,p] at previous timestep t (t=0: present; t>0: previous timestep)....
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
void unpin_elemental_pressure_dofs()
Unpin all the internal pressure freedoms.
void further_build()
Further build for Crouzeix_Raviart interpolates the internal pressure dofs from father element: Make ...
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons.
unsigned nvertex_node() const
Number of vertex nodes in the element.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 3 (velocities)
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
Refineable version of Spherical Quad Taylor Hood elements. (note that unlike the cartesian version th...
void pin_elemental_redundant_nodal_pressure_dofs()
Unpin the proper nodal pressure dofs.
Node * interpolating_node_pt(const unsigned &n, const int &n_value)
The velocities are isoparametric and so the "nodes" interpolating the velocities are the geometric no...
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
void interpolating_basis(const Vector< double > &s, Shape &psi, const int &n_value) const
The basis interpolating the pressure is given by pshape(). / The basis interpolating the velocity is ...
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 4 (3 velocities + 1 pressure)
Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node.
Node * get_interpolating_node_at_local_coordinate(const Vector< double > &s, const int &n_value)
The velocity nodes are the same as the geometric nodes. The pressure nodes must be calculated by usin...
unsigned ninterpolating_node(const int &n_value)
The number of pressure nodes is 4. The number of velocity nodes is the same as the number of geometri...
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
unsigned nvertex_node() const
Number of vertex nodes in the element.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned ninterpolating_node_1d(const int &n_value)
The number of 1d pressure nodes is 2, the number of 1d velocity nodes is the same as the number of 1d...
double local_one_d_fraction_of_interpolating_node(const unsigned &n1d, const unsigned &i, const int &n_value)
The pressure nodes are the corner nodes, so when n_value==DIM, the fraction is the same as the 1d nod...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at node n.
Refineable version of the Spherical Navier–Stokes equations.
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Pin unused nodal pressure dofs (empty by default, because by default pressure dofs are not associated...
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...
double geometric_jacobian(const Vector< double > &x)
Fill in the geometric Jacobian, which in this case is r*r*sin(theta)
static void pin_redundant_nodal_pressures(const Vector< GeneralisedElement * > &element_pt)
Loop over all elements in Vector (which typically contains all the elements in a fluid mesh) and pin ...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
static void unpin_all_pressure_dofs(const Vector< GeneralisedElement * > &element_pt)
Unpin all pressure dofs in elements listed in vector.
void fill_in_generic_residual_contribution_spherical_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add element's contribution to the elemental residual vector and/or Jacobian matrix flag=1: compute bo...
void further_build()
Further build: pass the pointers down to the sons.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
A class for elements that solve the Navier–Stokes equations, in axisymmetric spherical polar coordina...
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
SphericalNavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
double *& re_invro_pt()
Pointer to global inverse Froude number.
virtual unsigned u_index_spherical_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
double *& re_invfr_pt()
Pointer to global inverse Froude number.
Vector< double > * G_pt
Pointer to global gravity Vector.
double * ReInvRo_pt
Pointer to global Reynolds number x inverse Rossby number (used when in a rotating frame)
double * Re_pt
Pointer to global Reynolds number.
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
SphericalNavierStokesSourceFctPt & source_fct_pt()
Access function for the source-function pointer.
double interpolated_p_spherical_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
double *& re_pt()
Pointer to Reynolds number.
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number)
SphericalNavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
SphericalNavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
double *& density_ratio_pt()
Pointer to Density ratio.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
void interpolated_u_spherical_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
int son_type() const
Return son type.
Definition: tree.h:214
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Definition: tree.h:88
Tree * son_pt(const int &son_index) const
Return pointer to the son for a given index. Note that to aid code readability specific enums have be...
Definition: tree.h:103
Tree * father_pt() const
Return pointer to father: NULL if it's a root node.
Definition: tree.h:235
//////////////////////////////////////////////////////////////////// ////////////////////////////////...