refineable_linearised_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 linearised axisymmetric Navier-Stokes elements
27
28#ifndef OOMPH_REFINEABLE_LINEARISED_NAVIER_STOKES_ELEMENTS_HEADER
29#define OOMPH_REFINEABLE_LINEARISED_NAVIER_STOKES_ELEMENTS_HEADER
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36// oomph-lib includes
37#include "../generic/refineable_quad_element.h"
38#include "../generic/error_estimator.h"
40
41namespace oomph
42{
43 //=======================================================================
44 /// Refineable version of the linearised axisymmetric
45 /// Navier--Stokes equations
46 //=======================================================================
48 : public virtual LinearisedNavierStokesEquations,
49 public virtual RefineableElement,
50 public virtual ElementWithZ2ErrorEstimator
51 {
52 protected:
53 /// Pointer to n_p-th pressure node (Default: NULL,
54 /// indicating that pressure is not based on nodal interpolation).
55 virtual Node* pressure_node_pt(const unsigned& n_p)
56 {
57 return NULL;
58 }
59
60 /// Unpin all pressure dofs in the element
62
63 /// Pin unused nodal pressure dofs (empty by default, because
64 /// by default pressure dofs are not associated with nodes)
66
67 public:
68 /// Empty Constructor
73 {
74 }
75
76 /// Number of 'flux' terms for Z2 error estimation
78 {
79 // 3 diagonal strain rates, 3 off diagonal real and imaginary parts
80 return 2 * (DIM + ((DIM * DIM) - DIM) / 2);
81 }
82
83 /// Get 'flux' for Z2 error recovery: Upper triangular entries
84 /// in strain rate tensor.
86 {
87#ifdef PARANOID
88 unsigned num_entries = 2 * (DIM + ((DIM * DIM) - DIM) / 2);
89 if (flux.size() != num_entries)
90 {
91 std::ostringstream error_message;
92 error_message << "The flux vector has the wrong number of entries, "
93 << flux.size() << ", whereas it should be " << num_entries
94 << std::endl;
95 throw OomphLibError(
96 error_message.str(),
97 "RefineableLinearisedNavierStokesEquations::get_Z2_flux()",
98 OOMPH_EXCEPTION_LOCATION);
99 }
100#endif
101
102 // Get strain rate matrix
103 DenseMatrix<double> real_strainrate(DIM);
104 this->strain_rate(s, real_strainrate, 0);
105 DenseMatrix<double> imag_strainrate(DIM);
106 this->strain_rate(s, imag_strainrate, 1);
107
108
109 // Pack into flux Vector
110 unsigned icount = 0;
111
112 // Start with diagonal terms
113 for (unsigned i = 0; i < DIM; i++)
114 {
115 flux[icount] = real_strainrate(i, i);
116 icount++;
117 }
118
119 // Off diagonals row by row
120 for (unsigned i = 0; i < DIM; i++)
121 {
122 for (unsigned j = i + 1; j < DIM; j++)
123 {
124 flux[icount] = real_strainrate(i, j);
125 icount++;
126 }
127 }
128
129 // Start with diagonal terms
130 for (unsigned i = 0; i < DIM; i++)
131 {
132 flux[icount] = imag_strainrate(i, i);
133 icount++;
134 }
135
136 // Off diagonals row by row
137 for (unsigned i = 0; i < DIM; i++)
138 {
139 for (unsigned j = i + 1; j < DIM; j++)
140 {
141 flux[icount] = imag_strainrate(i, j);
142 icount++;
143 }
144 }
145 }
146
147 /// Further build: pass the pointers down to the sons
149 {
150 // Find the father element
151 RefineableLinearisedNavierStokesEquations* cast_father_element_pt =
153 this->father_element_pt());
154
155 // Set the viscosity ratio pointer
156 this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
157
158 // Set the density ratio pointer
159 this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
160
161 // Set pointer to global Reynolds number
162 this->Re_pt = cast_father_element_pt->re_pt();
163
164 // Set pointer to global Reynolds number x Strouhal number (=Womersley)
165 this->ReSt_pt = cast_father_element_pt->re_st_pt();
166
167 if (cast_father_element_pt->normalisation_element_pt() != 0)
168 {
169 // Pass down the normalisation element
171 cast_father_element_pt->normalisation_element_pt());
172 }
173
174 // Set the ALE_is_disabled flag
175 this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
176 }
177
178 /// Loop over all elements in Vector (which typically contains
179 /// all the elements in a fluid mesh) and pin the nodal pressure degrees
180 /// of freedom that are not being used. Function uses
181 /// the member function
182 /// - \c RefineableLinearisedNavierStokesEquations::
183 /// pin_all_nodal_pressure_dofs()
184 /// .
185 /// which is empty by default and should be implemented for
186 /// elements with nodal pressure degrees of freedom
187 /// (e.g. for refineable Taylor-Hood.)
189 const Vector<GeneralisedElement*>& element_pt)
190 {
191 // Loop over all elements to brutally pin all nodal pressure degrees of
192 // freedom
193 const unsigned n_element = element_pt.size();
194 for (unsigned e = 0; e < n_element; e++)
195 {
196 dynamic_cast<RefineableLinearisedNavierStokesEquations*>(element_pt[e])
198 }
199 }
200
201 /// Unpin all pressure dofs in elements listed in Vector
203 const Vector<GeneralisedElement*>& element_pt)
204 {
205 // Loop over all elements to brutally unpin all nodal pressure degrees of
206 // freedom and internal pressure degrees of freedom
207 const unsigned n_element = element_pt.size();
208 for (unsigned e = 0; e < n_element; e++)
209 {
210 dynamic_cast<RefineableLinearisedNavierStokesEquations*>(element_pt[e])
212 }
213 }
214
215
216 private:
217 /// Add element's contribution to the elemental residual vector
218 /// and/or Jacobian matrix
219 /// flag=1: compute both
220 /// flag=0: compute only residual vector
222 Vector<double>& residuals,
223 DenseMatrix<double>& jacobian,
224 DenseMatrix<double>& mass_matrix,
225 unsigned flag);
226
227 }; // End of RefineableLinearisedNavierStokesEquations class defn
228
229
230 /// ///////////////////////////////////////////////////////////////////////////
231 /// ///////////////////////////////////////////////////////////////////////////
232 /// ///////////////////////////////////////////////////////////////////////////
233
234
235 //=======================================================================
236 /// Refineable version of linearised axisymmetric quadratic
237 /// Crouzeix-Raviart elements
238 //=======================================================================
242 public virtual RefineableQElement<2>
243 {
244 private:
245 /// Unpin all the internal pressure freedoms
247 {
248 const unsigned n_pres = this->npres_linearised_nst();
249 // Loop over pressure dofs and unpin
250 for (unsigned l = 0; l < n_pres; l++)
251 {
252 // There are two pressure components
253 for (unsigned i = 0; i < 2; i++)
254 {
256 }
257 }
258 }
259
260 public:
261 /// Constructor
267 {
268 }
269
270 /// Number of continuously interpolated values: 4*DIM (velocities)
272 {
273 return 4 * DIM;
274 }
275
276 /// Rebuild from sons: Reconstruct pressure from the (merged) sons
277 void rebuild_from_sons(Mesh*& mesh_pt)
278 {
279 using namespace QuadTreeNames;
280
281 // Central pressure value:
282 // -----------------------
283
284 // Use average of the sons central pressure values
285 // Other options: Take average of the four (discontinuous)
286 // pressure values at the father's midpoint]
287
288 Vector<double> av_press(2, 0.0);
289
290 // Loop over the sons
291 for (unsigned ison = 0; ison < 4; ison++)
292 {
293 // Loop over the two pressure components
294 for (unsigned i = 0; i < 2; i++)
295 {
296 // Add the sons midnode pressure
297 av_press[i] +=
299 ->son_pt(ison)
300 ->object_pt()
302 ->value(0);
303 }
304 }
305
306 // Use the average
307 for (unsigned i = 0; i < 2; i++)
308 {
310 ->set_value(0, 0.25 * av_press[i]);
311 }
312
313 Vector<double> slope1(2, 0.0), slope2(2, 0.0);
314
315 // Loop over pressure components
316 for (unsigned i = 0; i < 2; i++)
317 {
318 // Slope in s_0 direction
319 // ----------------------
320
321 // Use average of the 2 FD approximations based on the
322 // elements central pressure values
323 // [Other options: Take average of the four
324 // pressure derivatives]
325
326 slope1[i] = quadtree_pt()
327 ->son_pt(SE)
328 ->object_pt()
330 ->value(0) -
332 ->son_pt(SW)
333 ->object_pt()
335 ->value(0);
336
337 slope2[i] = quadtree_pt()
338 ->son_pt(NE)
339 ->object_pt()
341 ->value(0) -
343 ->son_pt(NW)
344 ->object_pt()
346 ->value(0);
347
348 // Use the average
350 ->set_value(1, 0.5 * (slope1[i] + slope2[i]));
351
352 // Slope in s_1 direction
353 // ----------------------
354
355 // Use average of the 2 FD approximations based on the
356 // elements central pressure values
357 // [Other options: Take average of the four
358 // pressure derivatives]
359
360 slope1[i] = quadtree_pt()
361 ->son_pt(NE)
362 ->object_pt()
364 ->value(0) -
366 ->son_pt(SE)
367 ->object_pt()
369 ->value(0);
370
371 slope2[i] = quadtree_pt()
372 ->son_pt(NW)
373 ->object_pt()
375 ->value(0) -
377 ->son_pt(SW)
378 ->object_pt()
380 ->value(0);
381
382 // Use the average
384 ->set_value(2, 0.5 * (slope1[i] + slope2[i]));
385 } // End of loop over pressure components
386 }
387
388 /// Order of recovery shape functions for Z2 error estimation:
389 /// Same order as shape functions.
391 {
392 return 2;
393 }
394
395 /// Number of vertex nodes in the element
396 unsigned nvertex_node() const
397 {
399 }
400
401 /// Pointer to the j-th vertex node in the element
402 Node* vertex_node_pt(const unsigned& j) const
403 {
405 }
406
407 /// Get the function value u in Vector.
408 /// Note: Given the generality of the interface (this function
409 /// is usually called from black-box documentation or interpolation
410 /// routines), the values Vector sets its own size in here.
412 Vector<double>& values)
413 {
414 // Determine size of values Vector: U^C, U^S, W^C, W^S, V^C, V^S
415 const unsigned n_values = 4 * DIM;
416
417 // Set size of values Vector and initialise to zero
418 values.resize(n_values, 0.0);
419
420 // Calculate velocities: values[0],...
421 for (unsigned i = 0; i < n_values; i++)
422 {
424 }
425 }
426
427 /// Get all function values [U^C,U^S,...,P^S] at previous timestep t
428 /// (t=0: present; t>0: previous timestep).
429 /// \n
430 /// Note: Given the generality of the interface (this function is
431 /// usually called from black-box documentation or interpolation
432 /// routines), the values Vector sets its own size in here.
433 /// \n
434 /// Note: No pressure history is kept, so pressure is always
435 /// the current value.
436 void get_interpolated_values(const unsigned& t,
437 const Vector<double>& s,
438 Vector<double>& values)
439 {
440 // Set size of Vector: U^C, U^S, W^C, W^S, V^C, V^S
441 values.resize(4 * DIM);
442
443 // Initialise to zero
444 for (unsigned i = 0; i < 4 * DIM; i++)
445 {
446 values[i] = 0.0;
447 }
448
449 // Determine number of nodes in element
450 const unsigned n_node = nnode();
451
452 // Shape functions
453 Shape psif(n_node);
454 shape(s, psif);
455
456 // Calculate velocities: values[0],...
457 for (unsigned i = 0; i < (4 * DIM); i++)
458 {
459 // Get the local index at which the i-th velocity is stored
460 const unsigned u_local_index = u_index_linearised_nst(i);
461 for (unsigned l = 0; l < n_node; l++)
462 {
463 values[i] += nodal_value(t, l, u_local_index) * psif[l];
464 }
465 }
466 }
467
468 /// Perform additional hanging node procedures for variables
469 /// that are not interpolated by all nodes. Empty
471
472 /// Further build for Crouzeix_Raviart interpolates the internal
473 /// pressure dofs from father element: Make sure pressure values and
474 /// dp/ds agree between fathers and sons at the midpoints of the son
475 /// elements.
477 {
478 // Call the generic further build
480
481 using namespace QuadTreeNames;
482
483 // What type of son am I? Ask my quadtree representation...
484 int son_type = quadtree_pt()->son_type();
485
486 // Pointer to my father (in element impersonation)
487 RefineableElement* father_el_pt = quadtree_pt()->father_pt()->object_pt();
488
489 Vector<double> s_father(2);
490
491 // Son midpoint is located at the following coordinates in father element:
492
493 // South west son
494 if (son_type == SW)
495 {
496 s_father[0] = -0.5;
497 s_father[1] = -0.5;
498 }
499 // South east son
500 else if (son_type == SE)
501 {
502 s_father[0] = 0.5;
503 s_father[1] = -0.5;
504 }
505 // North east son
506 else if (son_type == NE)
507 {
508 s_father[0] = 0.5;
509 s_father[1] = 0.5;
510 }
511
512 // North west son
513 else if (son_type == NW)
514 {
515 s_father[0] = -0.5;
516 s_father[1] = 0.5;
517 }
518
519 // Pressure values in father element
520 // ---------------------------------
521
522 // Find pointer to father element
525 father_el_pt);
526
527 // Set up storage for pressure in father element
528 Vector<double> press(2, 0.0);
529
530 // Loop over pressure components
531 for (unsigned i = 0; i < 2; i++)
532 {
533 // Get pressure from father element
534 press[i] =
535 cast_father_el_pt->interpolated_p_linearised_nst(s_father, i);
536
537 // Pressure value gets copied straight into internal dof:
539 ->set_value(0, press[i]);
540
541 // The slopes get copied from father
543 ->set_value(1,
544 0.5 *
545 cast_father_el_pt
547 ->value(1));
548
550 ->set_value(2,
551 0.5 *
552 cast_father_el_pt
554 ->value(2));
555 } // End of loop over pressure components
556 }
557
558 }; // End of RefineableLinearisedQCrouzeixRaviartElement defn
559
560
561 /// ////////////////////////////////////////////////////////////////////////
562 /// ////////////////////////////////////////////////////////////////////////
563 /// ////////////////////////////////////////////////////////////////////////
564
565
566 //=======================================================================
567 /// Face geometry of the refineable linearised axisym
568 /// Crouzeix-Raviart elements
569 //=======================================================================
570 template<>
572 : public virtual FaceGeometry<LinearisedQCrouzeixRaviartElement>
573 {
574 public:
576 };
577
578
579 //=======================================================================
580 /// Face geometry of face geometric of the refineable linearised
581 /// axisym Crouzeix-Raviart elements
582 //=======================================================================
583 template<>
585 : public virtual FaceGeometry<
586 FaceGeometry<LinearisedQCrouzeixRaviartElement>>
587 {
588 public:
591 {
592 }
593 };
594
595
596 /// ///////////////////////////////////////////////////////////////////////////
597 /// ///////////////////////////////////////////////////////////////////////////
598 /// ///////////////////////////////////////////////////////////////////////////
599
600
601 //=======================================================================
602 /// Refineable version of linearised axisymmetric quadratic
603 /// Taylor-Hood elements
604 //=======================================================================
608 public virtual RefineableQElement<2>
609 {
610 private:
611 /// Pointer to n_p-th pressure node
612 Node* pressure_node_pt(const unsigned& n_p)
613 {
614 return this->node_pt(this->Pconv[n_p]);
615 }
616
617 /// Unpin all pressure dofs
619 {
620 // Determine number of nodes in element
621 const unsigned n_node = this->nnode();
622
623 // Get nodal indeces of the two pressure components
624 Vector<int> p_index(2);
625 for (unsigned i = 0; i < 2; i++)
626 {
627 p_index[i] = this->p_index_linearised_nst(i);
628 }
629
630 // Loop over nodes and unpin both pressure components
631 for (unsigned n = 0; n < n_node; n++)
632 {
633 for (unsigned i = 0; i < 2; i++)
634 {
635 this->node_pt(n)->unpin(p_index[i]);
636 }
637 }
638 }
639
640 /// Unpin the proper nodal pressure dofs
642 {
643 // Determine number of nodes in element
644 const unsigned n_node = this->nnode();
645
646 // Get nodal indeces of the two pressure components
647 Vector<int> p_index(2);
648 for (unsigned i = 0; i < 2; i++)
649 {
650 p_index[i] = this->p_index_linearised_nst(i);
651 }
652
653 // Loop over all nodes and pin all the nodal pressures
654 for (unsigned n = 0; n < n_node; n++)
655 {
656 for (unsigned i = 0; i < 2; i++)
657 {
658 this->node_pt(n)->pin(p_index[i]);
659 }
660 }
661
662 // Loop over all actual pressure nodes and unpin if they're not hanging
663 const unsigned n_pres = this->npres_linearised_nst();
664 for (unsigned l = 0; l < n_pres; l++)
665 {
666 Node* nod_pt = this->node_pt(this->Pconv[l]);
667 for (unsigned i = 0; i < 2; i++)
668 {
669 if (!nod_pt->is_hanging(p_index[i]))
670 {
671 nod_pt->unpin(p_index[i]);
672 }
673 }
674 }
675 }
676
677 public:
678 /// Constructor:
684 {
685 }
686
687 /// Number of values (pinned or dofs) required at node n.
688 /// Bumped up to 8 so we don't have to worry if a hanging mid-side node
689 /// gets shared by a corner node (which has extra degrees of freedom)
690 unsigned required_nvalue(const unsigned& n) const
691 {
692 return 8;
693 }
694
695 /// Number of continuously interpolated values: 8
696 /// (6 velocities + 2 pressures)
698 {
699 return 8;
700 }
701
702 /// Rebuild from sons: empty
703 void rebuild_from_sons(Mesh*& mesh_pt) {}
704
705 /// Order of recovery shape functions for Z2 error estimation:
706 /// Same order as shape functions.
708 {
709 return 2;
710 }
711
712 /// Number of vertex nodes in the element
713 unsigned nvertex_node() const
714 {
716 }
717
718 /// Pointer to the j-th vertex node in the element
719 Node* vertex_node_pt(const unsigned& j) const
720 {
722 }
723
724 /// Get the function value u in Vector.
725 /// Note: Given the generality of the interface (this function
726 /// is usually called from black-box documentation or interpolation
727 /// routines), the values Vector sets its own size in here.
729 Vector<double>& values)
730 {
731 // Determine size of values Vector:
732 // U^C, U^S, W^C, W^S, V^C, V^S, P^C, P^S
733 const unsigned n_values = 2 * (DIM + 1);
734
735 // Set size of values Vector and initialise to zero
736 values.resize(n_values, 0.0);
737
738 // Calculate velocities: values[0],...
739 for (unsigned i = 0; i < (2 * DIM); i++)
740 {
742 }
743
744 // Calculate pressure: values[DIM], values[DIM+1]
745 for (unsigned i = 0; i < 2; i++)
746 {
747 values[DIM + i] = interpolated_p_linearised_nst(s, i);
748 }
749 }
750
751 /// Get the function value u in Vector.
752 /// Note: Given the generality of the interface (this function
753 /// is usually called from black-box documentation or interpolation
754 /// routines), the values Vector sets its own size in here.
755 void get_interpolated_values(const unsigned& t,
756 const Vector<double>& s,
757 Vector<double>& values)
758 {
759 // Set size of values Vector: U^C, U^S, W^C, W^S, V^C, V^S, P^C, P^S
760 values.resize(2 * (DIM + 1));
761
762 // Initialise all entries to zero
763 for (unsigned i = 0; i < 2 * (DIM + 1); i++)
764 {
765 values[i] = 0.0;
766 }
767
768 // Determine number of nodes in element
769 const unsigned n_node = nnode();
770
771 // Shape functions
772 Shape psif(n_node);
773 shape(s, psif);
774
775 // Calculate velocities: values[0],...
776 for (unsigned i = 0; i < (2 * DIM); i++)
777 {
778 // Get the nodal coordinate of the velocity
779 const unsigned u_nodal_index = u_index_linearised_nst(i);
780 for (unsigned l = 0; l < n_node; l++)
781 {
782 values[i] += nodal_value(t, l, u_nodal_index) * psif[l];
783 }
784 }
785
786 // Calculate pressure: values[DIM], values[DIM+1]
787 // (no history is carried in the pressure)
788 for (unsigned i = 0; i < 2; i++)
789 {
790 values[DIM + i] = interpolated_p_linearised_nst(s, i);
791 }
792 }
793
794 /// Perform additional hanging node procedures for variables
795 /// that are not interpolated by all nodes. The two pressure components
796 /// are stored at the 6th and 7th location in each node
798 {
799 for (unsigned i = 0; i < 2; i++)
800 {
801 this->setup_hang_for_value((2 * DIM) + i);
802 }
803 }
804
805 /// The velocities are isoparametric and so the "nodes"
806 /// interpolating the velocities are the geometric nodes. The
807 /// pressure "nodes" are a subset of the nodes, so when n_value==6
808 /// or 7, the n-th pressure node is returned.
809 Node* interpolating_node_pt(const unsigned& n, const int& n_value)
810
811 {
812 // The only different nodes are the pressure nodes
813 if (n_value == (2 * DIM) || n_value == ((2 * DIM) + 1))
814 {
815 return this->node_pt(this->Pconv[n]);
816 }
817
818 // The other variables are interpolated via the usual nodes
819 else
820 {
821 return this->node_pt(n);
822 }
823 }
824
825 /// The pressure nodes are the corner nodes, so when n_value==6
826 /// or 7, the fraction is the same as the 1d node number, 0 or 1.
828 const unsigned& i,
829 const int& n_value)
830 {
831 if (n_value == (2 * DIM) || n_value == ((2 * DIM) + 1))
832 {
833 // The pressure nodes are just located on the boundaries at 0 or 1
834 return double(n1d);
835 }
836 // Otherwise the velocity nodes are the same as the geometric ones
837 else
838 {
839 return this->local_one_d_fraction_of_node(n1d, i);
840 }
841 }
842
843 /// The velocity nodes are the same as the geometric nodes.
844 /// The pressure nodes must be calculated by using the same methods
845 /// as the geometric nodes, but by recalling that there are only two
846 /// pressure nodes per edge.
848 const int& n_value)
849 {
850 // If we are calculating pressure nodes
851 if (n_value == static_cast<int>(2 * DIM) ||
852 n_value == static_cast<int>((2 * DIM) + 1))
853 {
854 // Storage for the index of the pressure node
855 unsigned total_index = 0;
856 // The number of nodes along each 1d edge is 2.
857 const unsigned NNODE_1D = 2;
858 // Storage for the index along each boundary
859 // Note that it's only a 2D spatial element
860 Vector<int> index(2);
861 // Loop over the coordinates
862 for (unsigned i = 0; i < 2; i++)
863 {
864 // If we are at the lower limit, the index is zero
865 if (s[i] == -1.0)
866 {
867 index[i] = 0;
868 }
869
870 // If we are at the upper limit, the index is the number of nodes
871 // minus 1
872 else if (s[i] == 1.0)
873 {
874 index[i] = NNODE_1D - 1;
875 }
876
877 // Otherwise, we have to calculate the index in general
878 else
879 {
880 // For uniformly spaced nodes the 0th node number would be
881 double float_index = 0.5 * (1.0 + s[i]) * (NNODE_1D - 1);
882 index[i] = int(float_index);
883 // What is the excess. This should be safe because the
884 // taking the integer part rounds down
885 double excess = float_index - index[i];
886 // If the excess is bigger than our tolerance there is no node,
887 // return null
890 {
891 return 0;
892 }
893 }
894 /// Construct the general pressure index from the components.
895 total_index +=
896 index[i] * static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
897 static_cast<int>(i)));
898 }
899 // If we've got here we have a node, so let's return a pointer to it
900 return this->node_pt(this->Pconv[total_index]);
901 }
902 // Otherwise velocity nodes are the same as pressure nodes
903 else
904 {
905 return this->get_node_at_local_coordinate(s);
906 }
907 }
908
909
910 /// The number of 1d pressure nodes is 2, the number of 1d
911 /// velocity nodes is the same as the number of 1d geometric nodes.
912 unsigned ninterpolating_node_1d(const int& n_value)
913 {
914 if (n_value == (2 * DIM) || n_value == ((2 * DIM) + 1))
915 {
916 return 2;
917 }
918 else
919 {
920 return this->nnode_1d();
921 }
922 }
923
924 /// The number of pressure nodes is 4. The number of
925 /// velocity nodes is the same as the number of geometric nodes.
926 unsigned ninterpolating_node(const int& n_value)
927 {
928 if (n_value == (2 * DIM) || n_value == ((2 * DIM) + 1))
929 {
930 return 4;
931 }
932 else
933 {
934 return this->nnode();
935 }
936 }
937
938 /// The basis interpolating the pressure is given by pshape().
939 /// / The basis interpolating the velocity is shape().
941 Shape& psi,
942 const int& n_value) const
943 {
944 if (n_value == (2 * DIM) || n_value == ((2 * DIM) + 1))
945 {
946 return this->pshape_linearised_nst(s, psi);
947 }
948 else
949 {
950 return this->shape(s, psi);
951 }
952 }
953
954 }; // End of RefineableLinearisedQTaylorHoodElement class defn
955
956
957 /// ////////////////////////////////////////////////////////////////////////
958 /// ////////////////////////////////////////////////////////////////////////
959 /// ////////////////////////////////////////////////////////////////////////
960
961
962 //=======================================================================
963 /// Face geometry of the refineable linearised axisym
964 /// Taylor-Hood elements
965 //=======================================================================
966 template<>
968 : public virtual FaceGeometry<LinearisedQTaylorHoodElement>
969 {
970 public:
972 };
973
974
975 //=======================================================================
976 /// Face geometry of face geometric of the refineable linearised
977 /// axisym Taylor-Hood elements
978 //=======================================================================
979 template<>
981 : public virtual FaceGeometry<FaceGeometry<LinearisedQTaylorHoodElement>>
982 {
983 public:
985 {
986 }
987 };
988
989
990} // End of namespace oomph
991
992#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 class for elements that solve the linearised version of the unsteady Navier–Stokes equations in cyl...
virtual unsigned u_index_linearised_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate, const unsigned &real) const
Strain-rate tensor: where (in that order)
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed....
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
void set_eigenfunction_normalisation_element(LinearisedNavierStokesEigenfunctionNormalisationElement *const &normalisation_el_pt)
the boolean flag check_nodal_data is set to false.
double *& re_pt()
Pointer to Reynolds number.
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
double interpolated_p_linearised_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated pressure p[i] at local coordinate s.
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
double * Re_pt
Pointer to global Reynolds number.
double *& viscosity_ratio_pt()
Pointer to the viscosity ratio.
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
double *& density_ratio_pt()
Pointer to the density ratio.
double interpolated_u_linearised_nst(const Vector< double > &s, const unsigned &i) const
Compute the element's residual Vector and the jacobian matrix. Virtual function can be overloaded by ...
LinearisedNavierStokesEigenfunctionNormalisationElement * normalisation_element_pt()
Pointer to normalisation element.
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
unsigned npres_linearised_nst() const
Return number of pressure values corresponding to a single pressure component.
Vector< unsigned > P_linearised_nst_internal_index
Internal indices that indicate at which internal data the pressure values are stored....
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
unsigned npres_linearised_nst() const
Return number of pressure values corresponding to a single pressure component.
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
void pshape_linearised_nst(const Vector< double > &s, Shape &psi) const
Compute the pressure shape functions at local coordinate s.
virtual int p_index_linearised_nst(const unsigned &i) const
Which nodal value represents the pressure? Overload version in base class which returns static int "P...
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....
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.
Refineable version of the linearised axisymmetric Navier–Stokes equations.
void further_build()
Further build: pass the pointers down to the sons.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
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 ...
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...
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 void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor.
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_linearised_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 for Crouzeix_Raviart interpolates the internal pressure dofs from father element: Make ...
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get all function values [U^C,U^S,...,P^S] at previous timestep t (t=0: present; t>0: previous timeste...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
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.
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons.
void unpin_elemental_pressure_dofs()
Unpin all the internal pressure freedoms.
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 ncont_interpolated_values() const
Number of continuously interpolated values: 4*DIM (velocities)
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
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==6 or 7, the fraction is the same as the 1d ...
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at node n. Bumped up to 8 so we don't have to worry if a h...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 8 (6 velocities + 2 pressures)
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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...
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 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 pin_elemental_redundant_nodal_pressure_dofs()
Unpin the proper nodal pressure dofs.
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...
unsigned nvertex_node() const
Number of vertex nodes 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...
Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node.
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 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 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
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...