refineable_linearised_axisym_navier_stokes_elements.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Header file for refineable linearised axisymmetric Navier-Stokes elements
27
28#ifndef OOMPH_REFINEABLE_LINEARISED_AXISYM_NAVIER_STOKES_ELEMENTS_HEADER
29#define OOMPH_REFINEABLE_LINEARISED_AXISYM_NAVIER_STOKES_ELEMENTS_HEADER
30
31// Config header generated by autoconfig
32#ifdef HAVE_CONFIG_H
33#include <oomph-lib-config.h>
34#endif
35
36// oomph-lib includes
37#include "../generic/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 //=======================================================================
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
80 return 6;
81 }
82
83 /// Get 'flux' for Z2 error recovery: Upper triangular entries
84 /// in strain rate tensor.
86 {
87 // Specify the number of velocity dimensions
88 unsigned DIM = 3;
89
90#ifdef PARANOID
91 unsigned num_entries = DIM + ((DIM * DIM) - DIM) / 2;
92 if (flux.size() != num_entries)
93 {
94 std::ostringstream error_message;
95 error_message << "The flux vector has the wrong number of entries, "
96 << flux.size() << ", whereas it should be " << num_entries
97 << std::endl;
98 throw OomphLibError(error_message.str(),
99 "RefineableLinearisedAxisymmetricNavierStokesEquati"
100 "ons::get_Z2_flux()",
101 OOMPH_EXCEPTION_LOCATION);
102 }
103#endif
104
105 // Get strain rate matrix
106 DenseMatrix<double> strainrate(DIM);
107 this->strain_rate(s, strainrate);
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] = 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] = strainrate(i, j);
125 icount++;
126 }
127 }
128 }
129
130 /// Fill in the geometric Jacobian, which in this case is r
132 {
133 return x[0];
134 }
135
136 /// Further build: pass the pointers down to the sons
138 {
139 // Find the father element
141 cast_father_element_pt =
143 this->father_element_pt());
144
145 // Set the viscosity ratio pointer
146 this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
147
148 // Set the density ratio pointer
149 this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
150
151 // Set pointer to global Reynolds number
152 this->Re_pt = cast_father_element_pt->re_pt();
153
154 // Set pointer to global Reynolds number x Strouhal number (=Womersley)
155 this->ReSt_pt = cast_father_element_pt->re_st_pt();
156
157 // Set pointer to azimuthal mode number
159 cast_father_element_pt->azimuthal_mode_number_pt();
160
161 // Set the ALE_is_disabled flag
162 this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
163 }
164
165 /// Loop over all elements in Vector (which typically contains
166 /// all the elements in a fluid mesh) and pin the nodal pressure degrees
167 /// of freedom that are not being used. Function uses
168 /// the member function
169 /// - \c RefineableLinearisedAxisymmetricNavierStokesEquations::
170 /// pin_all_nodal_pressure_dofs()
171 /// .
172 /// which is empty by default and should be implemented for
173 /// elements with nodal pressure degrees of freedom
174 /// (e.g. for refineable Taylor-Hood.)
176 const Vector<GeneralisedElement*>& element_pt)
177 {
178 // Loop over all elements to brutally pin all nodal pressure degrees of
179 // freedom
180 const unsigned n_element = element_pt.size();
181 for (unsigned e = 0; e < n_element; e++)
182 {
184 element_pt[e])
186 }
187 }
188
189 /// Unpin all pressure dofs in elements listed in Vector
191 const Vector<GeneralisedElement*>& element_pt)
192 {
193 // Loop over all elements to brutally unpin all nodal pressure degrees of
194 // freedom and internal pressure degrees of freedom
195 const unsigned n_element = element_pt.size();
196 for (unsigned e = 0; e < n_element; e++)
197 {
199 element_pt[e])
201 }
202 }
203
204
205 private:
206 /// Add element's contribution to the elemental residual vector
207 /// and/or Jacobian matrix
208 /// flag=1: compute both
209 /// flag=0: compute only residual vector
211 Vector<double>& residuals,
212 DenseMatrix<double>& jacobian,
213 DenseMatrix<double>& mass_matrix,
214 unsigned flag);
215
216 }; // End of RefineableLinearisedAxisymmetricNavierStokesEquations class defn
217
218
219 /// ///////////////////////////////////////////////////////////////////////////
220 /// ///////////////////////////////////////////////////////////////////////////
221 /// ///////////////////////////////////////////////////////////////////////////
222
223
224 //=======================================================================
225 /// Refineable version of linearised axisymmetric quadratic
226 /// Crouzeix-Raviart elements
227 //=======================================================================
231 public virtual RefineableQElement<2>
232 {
233 private:
234 /// Unpin all the internal pressure freedoms
236 {
237 const unsigned n_pres = this->npres_linearised_axi_nst();
238 // Loop over pressure dofs and unpin
239 for (unsigned l = 0; l < n_pres; l++)
240 {
241 // There are two pressure components
242 for (unsigned i = 0; i < 2; i++)
243 {
245 ->unpin(l);
246 }
247 }
248 }
249
250 public:
251 /// Constructor
257 {
258 }
259
260 /// Number of continuously interpolated values: 6 (velocities)
262 {
263 return 6;
264 }
265
266 /// Rebuild from sons: Reconstruct pressure from the (merged) sons
267 void rebuild_from_sons(Mesh*& mesh_pt)
268 {
269 using namespace QuadTreeNames;
270
271 // Central pressure value:
272 // -----------------------
273
274 // Use average of the sons central pressure values
275 // Other options: Take average of the four (discontinuous)
276 // pressure values at the father's midpoint]
277
278 Vector<double> av_press(2, 0.0);
279
280 // Loop over the sons
281 for (unsigned ison = 0; ison < 4; ison++)
282 {
283 // Loop over the two pressure components
284 for (unsigned i = 0; i < 2; i++)
285 {
286 // Add the sons midnode pressure
287 av_press[i] +=
289 ->son_pt(ison)
290 ->object_pt()
292 ->value(0);
293 }
294 }
295
296 // Use the average
297 for (unsigned i = 0; i < 2; i++)
298 {
300 ->set_value(0, 0.25 * av_press[i]);
301 }
302
303 Vector<double> slope1(2, 0.0), slope2(2, 0.0);
304
305 // Loop over pressure components
306 for (unsigned i = 0; i < 2; i++)
307 {
308 // Slope in s_0 direction
309 // ----------------------
310
311 // Use average of the 2 FD approximations based on the
312 // elements central pressure values
313 // [Other options: Take average of the four
314 // pressure derivatives]
315
316 slope1[i] = quadtree_pt()
317 ->son_pt(SE)
318 ->object_pt()
320 ->value(0) -
322 ->son_pt(SW)
323 ->object_pt()
325 ->value(0);
326
327 slope2[i] = quadtree_pt()
328 ->son_pt(NE)
329 ->object_pt()
331 ->value(0) -
333 ->son_pt(NW)
334 ->object_pt()
336 ->value(0);
337
338 // Use the average
340 ->set_value(1, 0.5 * (slope1[i] + slope2[i]));
341
342 // Slope in s_1 direction
343 // ----------------------
344
345 // Use average of the 2 FD approximations based on the
346 // elements central pressure values
347 // [Other options: Take average of the four
348 // pressure derivatives]
349
350 slope1[i] = quadtree_pt()
351 ->son_pt(NE)
352 ->object_pt()
354 ->value(0) -
356 ->son_pt(SE)
357 ->object_pt()
359 ->value(0);
360
361 slope2[i] = quadtree_pt()
362 ->son_pt(NW)
363 ->object_pt()
365 ->value(0) -
367 ->son_pt(SW)
368 ->object_pt()
370 ->value(0);
371
372 // Use the average
374 ->set_value(2, 0.5 * (slope1[i] + slope2[i]));
375 } // End of loop over pressure components
376 }
377
378 /// Order of recovery shape functions for Z2 error estimation:
379 /// Same order as shape functions.
381 {
382 return 2;
383 }
384
385 /// Number of vertex nodes in the element
386 unsigned nvertex_node() const
387 {
389 }
390
391 /// Pointer to the j-th vertex node in the element
392 Node* vertex_node_pt(const unsigned& j) const
393 {
395 }
396
397 /// Get the function value u in Vector.
398 /// Note: Given the generality of the interface (this function
399 /// is usually called from black-box documentation or interpolation
400 /// routines), the values Vector sets its own size in here.
402 Vector<double>& values)
403 {
404 // Set the velocity dimension of the element
405 const unsigned DIM = 3;
406
407 // Determine size of values Vector: U^C, U^S, W^C, W^S, V^C, V^S
408 const unsigned n_values = 2 * DIM;
409
410 // Set size of values Vector and initialise to zero
411 values.resize(n_values, 0.0);
412
413 // Calculate velocities: values[0],...
414 for (unsigned i = 0; i < n_values; i++)
415 {
417 }
418 }
419
420 /// Get all function values [U^C,U^S,...,P^S] at previous timestep t
421 /// (t=0: present; t>0: previous timestep).
422 /// \n
423 /// Note: Given the generality of the interface (this function is
424 /// usually called from black-box documentation or interpolation
425 /// routines), the values Vector sets its own size in here.
426 /// \n
427 /// Note: No pressure history is kept, so pressure is always
428 /// the current value.
429 void get_interpolated_values(const unsigned& t,
430 const Vector<double>& s,
431 Vector<double>& values)
432 {
433 const unsigned DIM = 3;
434
435 // Set size of Vector: U^C, U^S, W^C, W^S, V^C, V^S
436 values.resize(2 * DIM);
437
438 // Initialise to zero
439 for (unsigned i = 0; i < 2 * DIM; i++)
440 {
441 values[i] = 0.0;
442 }
443
444 // Determine number of nodes in element
445 const unsigned n_node = nnode();
446
447 // Shape functions
448 Shape psif(n_node);
449 shape(s, psif);
450
451 // Calculate velocities: values[0],...
452 for (unsigned i = 0; i < (2 * DIM); i++)
453 {
454 // Get the local index at which the i-th velocity is stored
455 const unsigned u_local_index = u_index_linearised_axi_nst(i);
456 for (unsigned l = 0; l < n_node; l++)
457 {
458 values[i] += nodal_value(t, l, u_local_index) * psif[l];
459 }
460 }
461 }
462
463 /// Perform additional hanging node procedures for variables
464 /// that are not interpolated by all nodes. Empty
466
467 /// Further build for Crouzeix_Raviart interpolates the internal
468 /// pressure dofs from father element: Make sure pressure values and
469 /// dp/ds agree between fathers and sons at the midpoints of the son
470 /// elements.
472 {
473 // Call the generic further build
475
476 using namespace QuadTreeNames;
477
478 // What type of son am I? Ask my quadtree representation...
479 int son_type = quadtree_pt()->son_type();
480
481 // Pointer to my father (in element impersonation)
482 RefineableElement* father_el_pt = quadtree_pt()->father_pt()->object_pt();
483
484 Vector<double> s_father(2);
485
486 // Son midpoint is located at the following coordinates in father element:
487
488 // South west son
489 if (son_type == SW)
490 {
491 s_father[0] = -0.5;
492 s_father[1] = -0.5;
493 }
494 // South east son
495 else if (son_type == SE)
496 {
497 s_father[0] = 0.5;
498 s_father[1] = -0.5;
499 }
500 // North east son
501 else if (son_type == NE)
502 {
503 s_father[0] = 0.5;
504 s_father[1] = 0.5;
505 }
506
507 // North west son
508 else if (son_type == NW)
509 {
510 s_father[0] = -0.5;
511 s_father[1] = 0.5;
512 }
513
514 // Pressure values in father element
515 // ---------------------------------
516
517 // Find pointer to father element
519 cast_father_el_pt = dynamic_cast<
521 father_el_pt);
522
523 // Set up storage for pressure in father element
524 Vector<double> press(2, 0.0);
525
526 // Loop over pressure components
527 for (unsigned i = 0; i < 2; i++)
528 {
529 // Get pressure from father element
530 press[i] =
531 cast_father_el_pt->interpolated_p_linearised_axi_nst(s_father, i);
532
533 // Pressure value gets copied straight into internal dof:
535 ->set_value(0, press[i]);
536
537 // The slopes get copied from father
539 ->set_value(
540 1,
541 0.5 * cast_father_el_pt
543 ->value(1));
544
546 ->set_value(
547 2,
548 0.5 * cast_father_el_pt
550 ->value(2));
551 } // End of loop over pressure components
552 }
553
554 }; // End of RefineableLinearisedAxisymmetricQCrouzeixRaviartElement defn
555
556
557 /// ////////////////////////////////////////////////////////////////////////
558 /// ////////////////////////////////////////////////////////////////////////
559 /// ////////////////////////////////////////////////////////////////////////
560
561
562 //=======================================================================
563 /// Face geometry of the refineable linearised axisym
564 /// Crouzeix-Raviart elements
565 //=======================================================================
566 template<>
568 : public virtual FaceGeometry<LinearisedAxisymmetricQCrouzeixRaviartElement>
569 {
570 public:
573 {
574 }
575 };
576
577
578 //=======================================================================
579 /// Face geometry of face geometric of the refineable linearised
580 /// axisym Crouzeix-Raviart elements
581 //=======================================================================
582 template<>
585 : public virtual FaceGeometry<
586 FaceGeometry<LinearisedAxisymmetricQCrouzeixRaviartElement>>
587 {
588 public:
590 : FaceGeometry<
592 {
593 }
594 };
595
596
597 /// ///////////////////////////////////////////////////////////////////////////
598 /// ///////////////////////////////////////////////////////////////////////////
599 /// ///////////////////////////////////////////////////////////////////////////
600
601
602 //=======================================================================
603 /// Refineable version of linearised axisymmetric quadratic
604 /// Taylor-Hood elements
605 //=======================================================================
609 public virtual RefineableQElement<2>
610 {
611 private:
612 /// Pointer to n_p-th pressure node
613 Node* pressure_node_pt(const unsigned& n_p)
614 {
615 return this->node_pt(this->Pconv[n_p]);
616 }
617
618 /// Unpin all pressure dofs
620 {
621 // Determine number of nodes in element
622 const unsigned n_node = this->nnode();
623
624 // Get nodal indeces of the two pressure components
625 Vector<int> p_index(2);
626 for (unsigned i = 0; i < 2; i++)
627 {
628 p_index[i] = this->p_index_linearised_axi_nst(i);
629 }
630
631 // Loop over nodes and unpin both pressure components
632 for (unsigned n = 0; n < n_node; n++)
633 {
634 for (unsigned i = 0; i < 2; i++)
635 {
636 this->node_pt(n)->unpin(p_index[i]);
637 }
638 }
639 }
640
641 /// Unpin the proper nodal pressure dofs
643 {
644 // Determine number of nodes in element
645 const unsigned n_node = this->nnode();
646
647 // Get nodal indeces of the two pressure components
648 Vector<int> p_index(2);
649 for (unsigned i = 0; i < 2; i++)
650 {
651 p_index[i] = this->p_index_linearised_axi_nst(i);
652 }
653
654 // Loop over all nodes and pin all the nodal pressures
655 for (unsigned n = 0; n < n_node; n++)
656 {
657 for (unsigned i = 0; i < 2; i++)
658 {
659 this->node_pt(n)->pin(p_index[i]);
660 }
661 }
662
663 // Loop over all actual pressure nodes and unpin if they're not hanging
664 const unsigned n_pres = this->npres_linearised_axi_nst();
665 for (unsigned l = 0; l < n_pres; l++)
666 {
667 Node* nod_pt = this->node_pt(this->Pconv[l]);
668 for (unsigned i = 0; i < 2; i++)
669 {
670 if (!nod_pt->is_hanging(p_index[i]))
671 {
672 nod_pt->unpin(p_index[i]);
673 }
674 }
675 }
676 }
677
678 public:
679 /// Constructor:
685 {
686 }
687
688 /// Number of values (pinned or dofs) required at node n.
689 /// Bumped up to 8 so we don't have to worry if a hanging mid-side node
690 /// gets shared by a corner node (which has extra degrees of freedom)
691 unsigned required_nvalue(const unsigned& n) const
692 {
693 return 8;
694 }
695
696 /// Number of continuously interpolated values: 8
697 /// (6 velocities + 2 pressures)
699 {
700 return 8;
701 }
702
703 /// Rebuild from sons: empty
704 void rebuild_from_sons(Mesh*& mesh_pt) {}
705
706 /// Order of recovery shape functions for Z2 error estimation:
707 /// Same order as shape functions.
709 {
710 return 2;
711 }
712
713 /// Number of vertex nodes in the element
714 unsigned nvertex_node() const
715 {
717 }
718
719 /// Pointer to the j-th vertex node in the element
720 Node* vertex_node_pt(const unsigned& j) const
721 {
723 }
724
725 /// Get the function value u in Vector.
726 /// Note: Given the generality of the interface (this function
727 /// is usually called from black-box documentation or interpolation
728 /// routines), the values Vector sets its own size in here.
730 Vector<double>& values)
731 {
732 // Set the velocity dimension of the element
733 const unsigned DIM = 3;
734
735 // Determine size of values Vector:
736 // U^C, U^S, W^C, W^S, V^C, V^S, P^C, P^S
737 const unsigned n_values = 2 * (DIM + 1);
738
739 // Set size of values Vector and initialise to zero
740 values.resize(n_values, 0.0);
741
742 // Calculate velocities: values[0],...
743 for (unsigned i = 0; i < (2 * DIM); i++)
744 {
746 }
747
748 // Calculate pressure: values[DIM], values[DIM+1]
749 for (unsigned i = 0; i < 2; i++)
750 {
751 values[DIM + i] = interpolated_p_linearised_axi_nst(s, i);
752 }
753 }
754
755 /// Get the function value u in Vector.
756 /// Note: Given the generality of the interface (this function
757 /// is usually called from black-box documentation or interpolation
758 /// routines), the values Vector sets its own size in here.
759 void get_interpolated_values(const unsigned& t,
760 const Vector<double>& s,
761 Vector<double>& values)
762 {
763 // Set the velocity dimension of the element
764 const unsigned DIM = 3;
765
766 // Set size of values Vector: U^C, U^S, W^C, W^S, V^C, V^S, P^C, P^S
767 values.resize(2 * (DIM + 1));
768
769 // Initialise all entries to zero
770 for (unsigned i = 0; i < 2 * (DIM + 1); i++)
771 {
772 values[i] = 0.0;
773 }
774
775 // Determine number of nodes in element
776 const unsigned n_node = nnode();
777
778 // Shape functions
779 Shape psif(n_node);
780 shape(s, psif);
781
782 // Calculate velocities: values[0],...
783 for (unsigned i = 0; i < (2 * DIM); i++)
784 {
785 // Get the nodal coordinate of the velocity
786 const unsigned u_nodal_index = u_index_linearised_axi_nst(i);
787 for (unsigned l = 0; l < n_node; l++)
788 {
789 values[i] += nodal_value(t, l, u_nodal_index) * psif[l];
790 }
791 }
792
793 // Calculate pressure: values[DIM], values[DIM+1]
794 // (no history is carried in the pressure)
795 for (unsigned i = 0; i < 2; i++)
796 {
797 values[DIM + i] = interpolated_p_linearised_axi_nst(s, i);
798 }
799 }
800
801 /// Perform additional hanging node procedures for variables
802 /// that are not interpolated by all nodes. The two pressure components
803 /// are stored at the 6th and 7th location in each node
805 {
806 const unsigned DIM = 3;
807 for (unsigned i = 0; i < 2; i++)
808 {
809 this->setup_hang_for_value((2 * DIM) + i);
810 }
811 }
812
813 /// The velocities are isoparametric and so the "nodes"
814 /// interpolating the velocities are the geometric nodes. The
815 /// pressure "nodes" are a subset of the nodes, so when n_value==6
816 /// or 7, the n-th pressure node is returned.
817 Node* interpolating_node_pt(const unsigned& n, const int& n_value)
818
819 {
820 const int DIM = 3;
821
822 // The only different nodes are the pressure nodes
823 if (n_value == (2 * DIM) || n_value == ((2 * DIM) + 1))
824 {
825 return this->node_pt(this->Pconv[n]);
826 }
827
828 // The other variables are interpolated via the usual nodes
829 else
830 {
831 return this->node_pt(n);
832 }
833 }
834
835 /// The pressure nodes are the corner nodes, so when n_value==6
836 /// or 7, the fraction is the same as the 1d node number, 0 or 1.
838 const unsigned& i,
839 const int& n_value)
840 {
841 const int DIM = 3;
842 if (n_value == (2 * DIM) || n_value == ((2 * DIM) + 1))
843 {
844 // The pressure nodes are just located on the boundaries at 0 or 1
845 return double(n1d);
846 }
847 // Otherwise the velocity nodes are the same as the geometric ones
848 else
849 {
850 return this->local_one_d_fraction_of_node(n1d, i);
851 }
852 }
853
854 /// The velocity nodes are the same as the geometric nodes.
855 /// The pressure nodes must be calculated by using the same methods
856 /// as the geometric nodes, but by recalling that there are only two
857 /// pressure nodes per edge.
859 const int& n_value)
860 {
861 const int DIM = 3;
862
863 // If we are calculating pressure nodes
864 if (n_value == static_cast<int>(2 * DIM) ||
865 n_value == static_cast<int>((2 * DIM) + 1))
866 {
867 // Storage for the index of the pressure node
868 unsigned total_index = 0;
869 // The number of nodes along each 1d edge is 2.
870 const unsigned NNODE_1D = 2;
871 // Storage for the index along each boundary
872 // Note that it's only a 2D spatial element
873 Vector<int> index(2);
874 // Loop over the coordinates
875 for (unsigned i = 0; i < 2; i++)
876 {
877 // If we are at the lower limit, the index is zero
878 if (s[i] == -1.0)
879 {
880 index[i] = 0;
881 }
882
883 // If we are at the upper limit, the index is the number of nodes
884 // minus 1
885 else if (s[i] == 1.0)
886 {
887 index[i] = NNODE_1D - 1;
888 }
889
890 // Otherwise, we have to calculate the index in general
891 else
892 {
893 // For uniformly spaced nodes the 0th node number would be
894 double float_index = 0.5 * (1.0 + s[i]) * (NNODE_1D - 1);
895 index[i] = int(float_index);
896 // What is the excess. This should be safe because the
897 // taking the integer part rounds down
898 double excess = float_index - index[i];
899 // If the excess is bigger than our tolerance there is no node,
900 // return null
903 {
904 return 0;
905 }
906 }
907 /// Construct the general pressure index from the components.
908 total_index +=
909 index[i] * static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
910 static_cast<int>(i)));
911 }
912 // If we've got here we have a node, so let's return a pointer to it
913 return this->node_pt(this->Pconv[total_index]);
914 }
915 // Otherwise velocity nodes are the same as pressure nodes
916 else
917 {
918 return this->get_node_at_local_coordinate(s);
919 }
920 }
921
922
923 /// The number of 1d pressure nodes is 2, the number of 1d
924 /// velocity nodes is the same as the number of 1d geometric nodes.
925 unsigned ninterpolating_node_1d(const int& n_value)
926 {
927 const int DIM = 3;
928 if (n_value == (2 * DIM) || n_value == ((2 * DIM) + 1))
929 {
930 return 2;
931 }
932 else
933 {
934 return this->nnode_1d();
935 }
936 }
937
938 /// The number of pressure nodes is 4. The number of
939 /// velocity nodes is the same as the number of geometric nodes.
940 unsigned ninterpolating_node(const int& n_value)
941 {
942 const int DIM = 3;
943 if (n_value == (2 * DIM) || n_value == ((2 * DIM) + 1))
944 {
945 return 4;
946 }
947 else
948 {
949 return this->nnode();
950 }
951 }
952
953 /// The basis interpolating the pressure is given by pshape().
954 /// / The basis interpolating the velocity is shape().
956 Shape& psi,
957 const int& n_value) const
958 {
959 const int DIM = 3;
960 if (n_value == (2 * DIM) || n_value == ((2 * DIM) + 1))
961 {
962 return this->pshape_linearised_axi_nst(s, psi);
963 }
964 else
965 {
966 return this->shape(s, psi);
967 }
968 }
969
970 }; // End of RefineableLinearisedAxisymmetricQTaylorHoodElement class defn
971
972
973 /// ////////////////////////////////////////////////////////////////////////
974 /// ////////////////////////////////////////////////////////////////////////
975 /// ////////////////////////////////////////////////////////////////////////
976
977
978 //=======================================================================
979 /// Face geometry of the refineable linearised axisym
980 /// Taylor-Hood elements
981 //=======================================================================
982 template<>
984 : public virtual FaceGeometry<LinearisedAxisymmetricQTaylorHoodElement>
985 {
986 public:
988 };
989
990
991 //=======================================================================
992 /// Face geometry of face geometric of the refineable linearised
993 /// axisym Taylor-Hood elements
994 //=======================================================================
995 template<>
998 : public virtual FaceGeometry<
999 FaceGeometry<LinearisedAxisymmetricQTaylorHoodElement>>
1000 {
1001 public:
1004 {
1005 }
1006 };
1007
1008
1009} // End of namespace oomph
1010
1011#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...
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
int * Azimuthal_Mode_Number_pt
Pointer to azimuthal mode number k in e^ik(theta) decomposition.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed....
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
double interpolated_u_linearised_axi_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated velocity u[i] at local coordinate s.
double interpolated_p_linearised_axi_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated pressure p[i] at local coordinate s.
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
int *& azimuthal_mode_number_pt()
Pointer to azimuthal mode number k in e^ik(theta) decomposition.
virtual unsigned u_index_linearised_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
Vector< unsigned > P_linearised_axi_nst_internal_index
Internal indices that indicate at which internal data the pressure values are stored....
unsigned npres_linearised_axi_nst() const
Return number of pressure values corresponding to a single pressure component.
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi) const
Compute the pressure shape functions at local coordinate s.
unsigned npres_linearised_axi_nst() const
Return number of pressure values corresponding to a single pressure component.
virtual int p_index_linearised_axi_nst(const unsigned &i) const
Which nodal value represents the pressure? Overload version in base class which returns static int "P...
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
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.
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 fill_in_generic_residual_contribution_linearised_axi_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...
double geometric_jacobian(const Vector< double > &x)
Fill in the geometric Jacobian, which in this case is r.
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 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 get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor.
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Pin unused nodal pressure dofs (empty by default, because by default pressure dofs are not associated...
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
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...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 6 (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...
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
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 * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
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 ...
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 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 further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
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 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...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 8 (6 velocities + 2 pressures)
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 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...
Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node.
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 ...
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...