refineable_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 axisymmetric quad Navier Stokes elements
27#ifndef OOMPH_REFINEABLE_AXISYMMETRIC_NAVIER_STOKES_ELEMENTS_HEADER
28#define OOMPH_REFINEABLE_AXISYMMETRIC_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 Axisymmetric Navier--Stokes equations
44 //======================================================================
46 : public virtual AxisymmetricNavierStokesEquations,
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
88#ifdef PARANOID
89 unsigned num_entries = DIM + ((DIM * DIM) - DIM) / 2;
90 if (flux.size() < num_entries)
91 {
92 std::ostringstream error_message;
93 error_message << "The flux vector has the wrong number of entries, "
94 << flux.size() << ", whereas it should be " << num_entries
95 << std::endl;
96 throw OomphLibError(error_message.str(),
97 OOMPH_CURRENT_FUNCTION,
98 OOMPH_EXCEPTION_LOCATION);
99 }
100#endif
101
102
103 // Get strain rate matrix
104 DenseMatrix<double> strainrate(DIM);
105 this->strain_rate(s, strainrate);
106
107 // Pack into flux Vector
108 unsigned icount = 0;
109
110 // Start with diagonal terms
111 for (unsigned i = 0; i < DIM; i++)
112 {
113 flux[icount] = strainrate(i, i);
114 icount++;
115 }
116
117 // Off diagonals row by row
118 for (unsigned i = 0; i < DIM; i++)
119 {
120 for (unsigned j = i + 1; j < DIM; j++)
121 {
122 flux[icount] = strainrate(i, j);
123 icount++;
124 }
125 }
126 }
127
128 /// Fill in the geometric Jacobian, which in this case is r
130 {
131 return x[0];
132 }
133
134
135 /// Further build: pass the pointers down to the sons
137 {
138 // Find the father element
139 RefineableAxisymmetricNavierStokesEquations* cast_father_element_pt =
141 this->father_element_pt());
142
143 // Set the viscosity ratio pointer
144 this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
145 // Set the density ratio pointer
146 this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
147 // Set pointer to global Reynolds number
148 this->Re_pt = cast_father_element_pt->re_pt();
149 // Set pointer to global Reynolds number x Strouhal number (=Womersley)
150 this->ReSt_pt = cast_father_element_pt->re_st_pt();
151 // Set pointer to global Reynolds number x inverse Froude number
152 this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
153 // Set pointer to the global Reynolds number x inverse Rossby number
154 this->ReInvRo_pt = cast_father_element_pt->re_invro_pt();
155 // Set pointer to global gravity Vector
156 this->G_pt = cast_father_element_pt->g_pt();
157
158 // Set pointer to body force function
159 this->Body_force_fct_pt =
160 cast_father_element_pt->axi_nst_body_force_fct_pt();
161
162 // Set pointer to volumetric source function
163 this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
164
165 // Set the ALE_is_disabled flag
166 this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
167 }
168
169 /// Compute the derivatives of the i-th component of
170 /// velocity at point s with respect
171 /// to all data that can affect its value. In addition, return the global
172 /// equation numbers corresponding to the data.
173 /// Overload the non-refineable version to take account of hanging node
174 /// information
176 const unsigned& i,
177 Vector<double>& du_ddata,
178 Vector<unsigned>& global_eqn_number)
179 {
180 // Find number of nodes
181 unsigned n_node = this->nnode();
182 // Local shape function
183 Shape psi(n_node);
184 // Find values of shape function at the given local coordinate
185 this->shape(s, psi);
186
187 // Find the index at which the velocity component is stored
188 const unsigned u_nodal_index = this->u_index_axi_nst(i);
189
190 // Storage for hang info pointer
191 HangInfo* hang_info_pt = 0;
192 // Storage for global equation
193 int global_eqn = 0;
194
195 // Find the number of dofs associated with interpolated u
196 unsigned n_u_dof = 0;
197 for (unsigned l = 0; l < n_node; l++)
198 {
199 unsigned n_master = 1;
200
201 // Local bool (is the node hanging)
202 bool is_node_hanging = this->node_pt(l)->is_hanging();
203
204 // If the node is hanging, get the number of master nodes
205 if (is_node_hanging)
206 {
207 hang_info_pt = this->node_pt(l)->hanging_pt();
208 n_master = hang_info_pt->nmaster();
209 }
210 // Otherwise there is just one master node, the node itself
211 else
212 {
213 n_master = 1;
214 }
215
216 // Loop over the master nodes
217 for (unsigned m = 0; m < n_master; m++)
218 {
219 // Get the equation number
220 if (is_node_hanging)
221 {
222 // Get the equation number from the master node
223 global_eqn =
224 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
225 }
226 else
227 {
228 // Global equation number
229 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
230 }
231
232 // If it's positive add to the count
233 if (global_eqn >= 0)
234 {
235 ++n_u_dof;
236 }
237 }
238 }
239
240 // Now resize the storage schemes
241 du_ddata.resize(n_u_dof, 0.0);
242 global_eqn_number.resize(n_u_dof, 0);
243
244 // Loop over th nodes again and set the derivatives
245 unsigned count = 0;
246 // Loop over the local nodes and sum
247 for (unsigned l = 0; l < n_node; l++)
248 {
249 unsigned n_master = 1;
250 double hang_weight = 1.0;
251
252 // Local bool (is the node hanging)
253 bool is_node_hanging = this->node_pt(l)->is_hanging();
254
255 // If the node is hanging, get the number of master nodes
256 if (is_node_hanging)
257 {
258 hang_info_pt = this->node_pt(l)->hanging_pt();
259 n_master = hang_info_pt->nmaster();
260 }
261 // Otherwise there is just one master node, the node itself
262 else
263 {
264 n_master = 1;
265 }
266
267 // Loop over the master nodes
268 for (unsigned m = 0; m < n_master; m++)
269 {
270 // If the node is hanging get weight from master node
271 if (is_node_hanging)
272 {
273 // Get the hang weight from the master node
274 hang_weight = hang_info_pt->master_weight(m);
275 }
276 else
277 {
278 // Node contributes with full weight
279 hang_weight = 1.0;
280 }
281
282 // Get the equation number
283 if (is_node_hanging)
284 {
285 // Get the equation number from the master node
286 global_eqn =
287 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
288 }
289 else
290 {
291 // Local equation number
292 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
293 }
294
295 if (global_eqn >= 0)
296 {
297 // Set the global equation number
298 global_eqn_number[count] = global_eqn;
299 // Set the derivative with respect to the unknown
300 du_ddata[count] = psi[l] * hang_weight;
301 // Increase the counter
302 ++count;
303 }
304 }
305 }
306 }
307
308
309 /// Loop over all elements in Vector (which typically contains
310 /// all the elements in a fluid mesh) and pin the nodal pressure degrees
311 /// of freedom that are not being used. Function uses
312 /// the member function
313 /// - \c RefineableAxisymmetricNavierStokesEquations::
314 /// pin_all_nodal_pressure_dofs()
315 /// .
316 /// which is empty by default and should be implemented for
317 /// elements with nodal pressure degrees of freedom
318 /// (e.g. for refineable Taylor-Hood.)
320 const Vector<GeneralisedElement*>& element_pt)
321 {
322 // Loop over all elements to brutally pin all nodal pressure degrees of
323 // freedom
324 unsigned n_element = element_pt.size();
325 for (unsigned e = 0; e < n_element; e++)
326 {
328 element_pt[e])
330 }
331 }
332
333 /// Unpin all pressure dofs in elements listed in vector.
335 const Vector<GeneralisedElement*>& element_pt)
336 {
337 // Loop over all elements to brutally unpin all nodal pressure degrees of
338 // freedom and internal pressure degrees of freedom
339 unsigned n_element = element_pt.size();
340 for (unsigned e = 0; e < n_element; e++)
341 {
343 element_pt[e])
345 }
346 }
347
348 /// Compute derivatives of elemental residual vector with respect to
349 /// nodal coordinates. This function computes these terms analytically and
350 /// overwrites the default implementation in the FiniteElement base class.
351 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
353 RankThreeTensor<double>& dresidual_dnodal_coordinates);
354
355 private:
356 /// Add element's contribution to the elemental residual vector
357 /// and/or Jacobian matrix and mass matrix
358 /// flag=2: compute all
359 /// flag=1: compute both residual and Jacobian
360 /// flag=0: compute only residual vector
362 Vector<double>& residuals,
363 DenseMatrix<double>& jacobian,
364 DenseMatrix<double>& mass_matrix,
365 unsigned flag);
366
367 /// Add element's contribution to the derivative of the
368 /// elemental residual vector
369 /// and/or Jacobian matrix and/or mass matrix
370 /// flag=2: compute all
371 /// flag=1: compute both residual and Jacobian
372 /// flag=0: compute only residual vector
374 double* const& parameter_pt,
375 Vector<double>& dres_dparam,
376 DenseMatrix<double>& djac_dparam,
377 DenseMatrix<double>& dmass_matrix_dparam,
378 unsigned flag)
379 {
380 throw OomphLibError("Not yet implemented\n",
381 OOMPH_CURRENT_FUNCTION,
382 OOMPH_EXCEPTION_LOCATION);
383 }
384
385 /// Compute the hessian tensor vector products required to
386 /// perform continuation of bifurcations analytically
388 Vector<double> const& Y,
389 DenseMatrix<double> const& C,
390 DenseMatrix<double>& product)
391 {
392 throw OomphLibError("Not yet implemented\n",
393 OOMPH_CURRENT_FUNCTION,
394 OOMPH_EXCEPTION_LOCATION);
395 }
396 };
397
398 //======================================================================
399 /// Refineable version of Axisymmetric Quad Taylor Hood elements.
400 /// (note that unlike the cartesian version this is not scale-able
401 /// to higher dimensions!)
402 //======================================================================
406 public virtual RefineableQElement<2>
407 {
408 private:
409 /// Pointer to n_p-th pressure node
410 Node* pressure_node_pt(const unsigned& n_p)
411 {
412 return this->node_pt(this->Pconv[n_p]);
413 }
414
415 /// Unpin all pressure dofs
417 {
418 unsigned n_node = this->nnode();
419 int p_index = this->p_nodal_index_axi_nst();
420 // loop over nodes
421 for (unsigned n = 0; n < n_node; n++)
422 {
423 this->node_pt(n)->unpin(p_index);
424 }
425 }
426
427 /// Unpin the proper nodal pressure dofs
429 {
430 // Loop over all nodes
431 unsigned n_node = this->nnode();
432 int p_index = this->p_nodal_index_axi_nst();
433 // loop over all nodes and pin all the nodal pressures
434 for (unsigned n = 0; n < n_node; n++)
435 {
436 this->node_pt(n)->pin(p_index);
437 }
438
439 // Loop over all actual pressure nodes and unpin if they're not hanging
440 unsigned n_pres = this->npres_axi_nst();
441 for (unsigned l = 0; l < n_pres; l++)
442 {
443 Node* nod_pt = this->node_pt(this->Pconv[l]);
444 if (!nod_pt->is_hanging(p_index))
445 {
446 nod_pt->unpin(p_index);
447 }
448 }
449 }
450
451 public:
452 /// Constructor:
458 {
459 }
460
461 /// Number of values (pinned or dofs) required at node n.
462 // Bumped up to 4 so we don't have to worry if a hanging mid-side node
463 // gets shared by a corner node (which has extra degrees of freedom)
464 unsigned required_nvalue(const unsigned& n) const
465 {
466 return 4;
467 }
468
469 /// Number of continuously interpolated values: 4 (3 velocities + 1
470 /// pressure)
472 {
473 return 4;
474 }
475
476 /// Rebuild from sons: empty
477 void rebuild_from_sons(Mesh*& mesh_pt) {}
478
479 /// Order of recovery shape functions for Z2 error estimation:
480 /// Same order as shape functions.
482 {
483 return 2;
484 }
485
486 /// Number of vertex nodes in the element
487 unsigned nvertex_node() const
488 {
490 }
491
492 /// Pointer to the j-th vertex node in the element
493 Node* vertex_node_pt(const unsigned& j) const
494 {
496 }
497
498 /// Get the function value u in Vector.
499 /// Note: Given the generality of the interface (this function
500 /// is usually called from black-box documentation or interpolation
501 /// routines), the values Vector sets its own size in here.
503 Vector<double>& values)
504 {
505 // Set the velocity dimension of the element
506 unsigned DIM = 3;
507
508 // Set size of values Vector: u,w,v,p and initialise to zero
509 values.resize(DIM + 1, 0.0);
510
511 // Calculate velocities: values[0],...
512 for (unsigned i = 0; i < DIM; i++)
513 {
514 values[i] = interpolated_u_axi_nst(s, i);
515 }
516
517 // Calculate pressure: values[DIM]
518 values[DIM] = interpolated_p_axi_nst(s);
519 }
520
521 /// Get the function value u in Vector.
522 /// Note: Given the generality of the interface (this function
523 /// is usually called from black-box documentation or interpolation
524 /// routines), the values Vector sets its own size in here.
525 void get_interpolated_values(const unsigned& t,
526 const Vector<double>& s,
527 Vector<double>& values)
528 {
529 unsigned DIM = 3;
530
531 // Set size of Vector: u,w,v,p
532 values.resize(DIM + 1);
533
534 // Initialise
535 for (unsigned i = 0; i < DIM + 1; i++)
536 {
537 values[i] = 0.0;
538 }
539
540 // Find out how many nodes there are
541 unsigned n_node = nnode();
542
543 // Shape functions
544 Shape psif(n_node);
545 shape(s, psif);
546
547 // Calculate velocities: values[0],...
548 for (unsigned i = 0; i < DIM; i++)
549 {
550 // Get the nodal coordinate of the velocity
551 unsigned u_nodal_index = u_index_axi_nst(i);
552 for (unsigned l = 0; l < n_node; l++)
553 {
554 values[i] += nodal_value(t, l, u_nodal_index) * psif[l];
555 }
556 }
557
558 // Calculate pressure: values[DIM]
559 //(no history is carried in the pressure)
560 values[DIM] = interpolated_p_axi_nst(s);
561 }
562
563 /// Perform additional hanging node procedures for variables
564 /// that are not interpolated by all nodes. The pressures are stored
565 /// at the 3rd location in each node
567 {
568 int DIM = 3;
569 this->setup_hang_for_value(DIM);
570 }
571
572 /// The velocities are isoparametric and so the "nodes" interpolating
573 /// the velocities are the geometric nodes. The pressure "nodes" are a
574 /// subset of the nodes, so when n_value==DIM, the n-th pressure
575 /// node is returned.
576 Node* interpolating_node_pt(const unsigned& n, const int& n_value)
577
578 {
579 int DIM = 3;
580 // The only different nodes are the pressure nodes
581 if (n_value == DIM)
582 {
583 return this->node_pt(this->Pconv[n]);
584 }
585 // The other variables are interpolated via the usual nodes
586 else
587 {
588 return this->node_pt(n);
589 }
590 }
591
592 /// The pressure nodes are the corner nodes, so when n_value==DIM,
593 /// the fraction is the same as the 1d node number, 0 or 1.
595 const unsigned& i,
596 const int& n_value)
597 {
598 int DIM = 3;
599 if (n_value == DIM)
600 {
601 // The pressure nodes are just located on the boundaries at 0 or 1
602 return double(n1d);
603 }
604 // Otherwise the velocity nodes are the same as the geometric ones
605 else
606 {
607 return this->local_one_d_fraction_of_node(n1d, i);
608 }
609 }
610
611 /// The velocity nodes are the same as the geometric nodes. The
612 /// pressure nodes must be calculated by using the same methods as
613 /// the geometric nodes, but by recalling that there are only two pressure
614 /// nodes per edge.
616 const int& n_value)
617 {
618 int DIM = 3;
619 // If we are calculating pressure nodes
620 if (n_value == static_cast<int>(DIM))
621 {
622 // Storage for the index of the pressure node
623 unsigned total_index = 0;
624 // The number of nodes along each 1d edge is 2.
625 unsigned NNODE_1D = 2;
626 // Storage for the index along each boundary
627 // Note that it's only a 2D spatial element
628 Vector<int> index(2);
629 // Loop over the coordinates
630 for (unsigned i = 0; i < 2; i++)
631 {
632 // If we are at the lower limit, the index is zero
633 if (s[i] == -1.0)
634 {
635 index[i] = 0;
636 }
637 // If we are at the upper limit, the index is the number of nodes
638 // minus 1
639 else if (s[i] == 1.0)
640 {
641 index[i] = NNODE_1D - 1;
642 }
643 // Otherwise, we have to calculate the index in general
644 else
645 {
646 // For uniformly spaced nodes the 0th node number would be
647 double float_index = 0.5 * (1.0 + s[i]) * (NNODE_1D - 1);
648 index[i] = int(float_index);
649 // What is the excess. This should be safe because the
650 // taking the integer part rounds down
651 double excess = float_index - index[i];
652 // If the excess is bigger than our tolerance there is no node,
653 // return null
656 {
657 return 0;
658 }
659 }
660 /// Construct the general pressure index from the components.
661 total_index +=
662 index[i] * static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
663 static_cast<int>(i)));
664 }
665 // If we've got here we have a node, so let's return a pointer to it
666 return this->node_pt(this->Pconv[total_index]);
667 }
668 // Otherwise velocity nodes are the same as pressure nodes
669 else
670 {
671 return this->get_node_at_local_coordinate(s);
672 }
673 }
674
675
676 /// The number of 1d pressure nodes is 2, the number of 1d velocity
677 /// nodes is the same as the number of 1d geometric nodes.
678 unsigned ninterpolating_node_1d(const int& n_value)
679 {
680 int DIM = 3;
681 if (n_value == DIM)
682 {
683 return 2;
684 }
685 else
686 {
687 return this->nnode_1d();
688 }
689 }
690
691 /// The number of pressure nodes is 4. The number of
692 /// velocity nodes is the same as the number of geometric nodes.
693 unsigned ninterpolating_node(const int& n_value)
694 {
695 int DIM = 3;
696 if (n_value == DIM)
697 {
698 return 4;
699 }
700 else
701 {
702 return this->nnode();
703 }
704 }
705
706 /// The basis interpolating the pressure is given by pshape().
707 /// / The basis interpolating the velocity is shape().
709 Shape& psi,
710 const int& n_value) const
711 {
712 int DIM = 3;
713 if (n_value == DIM)
714 {
715 return this->pshape_axi_nst(s, psi);
716 }
717 else
718 {
719 return this->shape(s, psi);
720 }
721 }
722 };
723
724
725 /// ////////////////////////////////////////////////////////////////////////
726 /// ////////////////////////////////////////////////////////////////////////
727 /// ////////////////////////////////////////////////////////////////////////
728
729 //=======================================================================
730 /// Face geometry of the RefineableQuadQTaylorHoodElements
731 //=======================================================================
732 template<>
734 : public virtual FaceGeometry<AxisymmetricQTaylorHoodElement>
735 {
736 public:
738 };
739
740
741 //=======================================================================
742 /// Face geometry of the RefineableQuadQTaylorHoodElements
743 //=======================================================================
744 template<>
746 : public virtual FaceGeometry<FaceGeometry<AxisymmetricQTaylorHoodElement>>
747 {
748 public:
751 {
752 }
753 };
754
755
756 //======================================================================
757 /// Refineable version of Axisymmetric Quad Crouzeix Raviart elements
758 /// (note that unlike the cartesian version this is not scale-able
759 /// to higher dimensions!)
760 //======================================================================
764 public virtual RefineableQElement<2>
765 {
766 private:
767 /// Unpin all the internal pressure freedoms
769 {
770 unsigned n_pres = this->npres_axi_nst();
771 // loop over pressure dofs and unpin
772 for (unsigned l = 0; l < n_pres; l++)
773 {
775 }
776 }
777
778 public:
779 /// Constructor:
785 {
786 }
787
788 /// Number of continuously interpolated values: 3 (velocities)
790 {
791 return 3;
792 }
793
794 /// Rebuild from sons: Reconstruct pressure from the (merged) sons
795 void rebuild_from_sons(Mesh*& mesh_pt)
796 {
797 using namespace QuadTreeNames;
798
799 // Central pressure value:
800 //-----------------------
801
802 // Use average of the sons central pressure values
803 // Other options: Take average of the four (discontinuous)
804 // pressure values at the father's midpoint]
805
806 double av_press = 0.0;
807
808 // Loop over the sons
809 for (unsigned ison = 0; ison < 4; ison++)
810 {
811 // Add the sons midnode pressure
812 av_press += quadtree_pt()
813 ->son_pt(ison)
814 ->object_pt()
816 ->value(0);
817 }
818
819 // Use the average
821
822
823 // Slope in s_0 direction
824 //----------------------
825
826 // Use average of the 2 FD approximations based on the
827 // elements central pressure values
828 // [Other options: Take average of the four
829 // pressure derivatives]
830
831 double slope1 = quadtree_pt()
832 ->son_pt(SE)
833 ->object_pt()
835 ->value(0) -
837 ->son_pt(SW)
838 ->object_pt()
840 ->value(0);
841
842 double slope2 = quadtree_pt()
843 ->son_pt(NE)
844 ->object_pt()
846 ->value(0) -
848 ->son_pt(NW)
849 ->object_pt()
851 ->value(0);
852
853
854 // Use the average
856 ->set_value(1, 0.5 * (slope1 + slope2));
857
858
859 // Slope in s_1 direction
860 //----------------------
861
862 // Use average of the 2 FD approximations based on the
863 // elements central pressure values
864 // [Other options: Take average of the four
865 // pressure derivatives]
866
867 slope1 = quadtree_pt()
868 ->son_pt(NE)
869 ->object_pt()
871 ->value(0) -
873 ->son_pt(SE)
874 ->object_pt()
876 ->value(0);
877
878 slope2 = quadtree_pt()
879 ->son_pt(NW)
880 ->object_pt()
882 ->value(0) -
884 ->son_pt(SW)
885 ->object_pt()
887 ->value(0);
888
889
890 // Use the average
892 ->set_value(2, 0.5 * (slope1 + slope2));
893 }
894
895 /// Order of recovery shape functions for Z2 error estimation:
896 /// Same order as shape functions.
898 {
899 return 2;
900 }
901
902 /// Number of vertex nodes in the element
903 unsigned nvertex_node() const
904 {
906 }
907
908 /// Pointer to the j-th vertex node in the element
909 Node* vertex_node_pt(const unsigned& j) const
910 {
912 }
913
914 /// Get the function value u in Vector.
915 /// Note: Given the generality of the interface (this function
916 /// is usually called from black-box documentation or interpolation
917 /// routines), the values Vector sets its own size in here.
919 Vector<double>& values)
920 {
921 unsigned DIM = 3;
922
923 // Set size of Vector: u,w,v and initialise to zero
924 values.resize(DIM, 0.0);
925
926 // Calculate velocities: values[0],...
927 for (unsigned i = 0; i < DIM; i++)
928 {
929 values[i] = interpolated_u_axi_nst(s, i);
930 }
931 }
932
933 /// Get all function values [u,v..,p] at previous timestep t
934 /// (t=0: present; t>0: previous timestep).
935 /// \n
936 /// Note: Given the generality of the interface (this function is
937 /// usually called from black-box documentation or interpolation routines),
938 /// the values Vector sets its own size in here.
939 /// \n
940 /// Note: No pressure history is kept, so pressure is always
941 /// the current value.
942 void get_interpolated_values(const unsigned& t,
943 const Vector<double>& s,
944 Vector<double>& values)
945 {
946 unsigned DIM = 3;
947
948 // Set size of Vector: u,w,v
949 values.resize(DIM);
950
951 // Initialise
952 for (unsigned i = 0; i < DIM; i++)
953 {
954 values[i] = 0.0;
955 }
956
957 // Find out how many nodes there are
958 unsigned n_node = nnode();
959
960 // Shape functions
961 Shape psif(n_node);
962 shape(s, psif);
963
964 // Calculate velocities: values[0],...
965 for (unsigned i = 0; i < DIM; i++)
966 {
967 // Get the local index at which the i-th velocity is stored
968 unsigned u_local_index = u_index_axi_nst(i);
969 for (unsigned l = 0; l < n_node; l++)
970 {
971 values[i] += nodal_value(t, l, u_local_index) * psif[l];
972 }
973 }
974 }
975
976 /// Perform additional hanging node procedures for variables
977 /// that are not interpolated by all nodes. Empty
979
980 /// Further build for Crouzeix_Raviart interpolates the internal
981 /// pressure dofs from father element: Make sure pressure values and
982 /// dp/ds agree between fathers and sons at the midpoints of the son
983 /// elements.
985 {
986 // Call the generic further build
988
989 using namespace QuadTreeNames;
990
991 // What type of son am I? Ask my quadtree representation...
992 int son_type = quadtree_pt()->son_type();
993
994 // Pointer to my father (in element impersonation)
995 RefineableElement* father_el_pt = quadtree_pt()->father_pt()->object_pt();
996
997 Vector<double> s_father(2);
998
999 // Son midpoint is located at the following coordinates in father element:
1000
1001 // South west son
1002 if (son_type == SW)
1003 {
1004 s_father[0] = -0.5;
1005 s_father[1] = -0.5;
1006 }
1007 // South east son
1008 else if (son_type == SE)
1009 {
1010 s_father[0] = 0.5;
1011 s_father[1] = -0.5;
1012 }
1013 // North east son
1014 else if (son_type == NE)
1015 {
1016 s_father[0] = 0.5;
1017 s_father[1] = 0.5;
1018 }
1019
1020 // North west son
1021 else if (son_type == NW)
1022 {
1023 s_father[0] = -0.5;
1024 s_father[1] = 0.5;
1025 }
1026
1027 // Pressure value in father element
1030 father_el_pt);
1031 double press = cast_father_el_pt->interpolated_p_axi_nst(s_father);
1032
1033 // Pressure value gets copied straight into internal dof:
1035
1036
1037 // The slopes get copied from father
1039 ->set_value(
1040 1,
1041 0.5 * cast_father_el_pt->internal_data_pt(P_axi_nst_internal_index)
1042 ->value(1));
1043
1045 ->set_value(
1046 2,
1047 0.5 * cast_father_el_pt->internal_data_pt(P_axi_nst_internal_index)
1048 ->value(2));
1049 }
1050 };
1051
1052
1053 //=======================================================================
1054 /// Face geometry of the RefineableQuadQCrouzeixRaviartElements
1055 //=======================================================================
1056 template<>
1058 : public virtual FaceGeometry<AxisymmetricQCrouzeixRaviartElement>
1059 {
1060 public:
1062 };
1063
1064
1065 //=======================================================================
1066 /// Face geometry of the RefineableQuadQCrouzeixRaviartElements
1067 //=======================================================================
1068 template<>
1071 : public virtual FaceGeometry<
1072 FaceGeometry<AxisymmetricQCrouzeixRaviartElement>>
1073 {
1074 public:
1077 {
1078 }
1079 };
1080
1081
1082} // namespace oomph
1083
1084#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
A class for elements that solve the unsteady axisymmetric Navier–Stokes equations in cylindrical pola...
double(* Source_fct_pt)(const double &time, const Vector< double > &x)
Pointer to volumetric source function.
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
Vector< double > * G_pt
Pointer to global gravity Vector.
void(*&)(const double &time, const Vector< double > &x, Vector< double > &f) axi_nst_body_force_fct_pt()
Access function for the body-force pointer.
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number)
double *& re_pt()
Pointer to Reynolds number.
double *& density_ratio_pt()
Pointer to Density ratio.
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
double(*&)(const double &time, const Vector< double > &x) source_fct_pt()
Access function for the source-function pointer.
double * Re_pt
Pointer to global Reynolds number.
void interpolated_u_axi_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
void(* Body_force_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to body force function.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed....
double *& re_invro_pt()
Pointer to global inverse Froude number.
double *& re_invfr_pt()
Pointer to global inverse Froude number.
double interpolated_p_axi_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
double * ReInvRo_pt
Pointer to global Reynolds number x inverse Rossby number (used when in a rotating frame)
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
virtual unsigned u_index_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
unsigned P_axi_nst_internal_index
Internal index that indicates at which internal data the pressure is stored.
unsigned npres_axi_nst() const
Return number of pressure values.
/////////////////////////////////////////////////////////////////////////
virtual int p_nodal_index_axi_nst() const
Which nodal value represents the pressure?
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
unsigned npres_axi_nst() const
Return number of pressure values.
void pshape_axi_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
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
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
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
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
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
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
Refineable version of the Axisymmetric Navier–Stokes equations.
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
double geometric_jacobian(const Vector< double > &x)
Fill in the geometric Jacobian, which in this case is r.
static void unpin_all_pressure_dofs(const Vector< GeneralisedElement * > &element_pt)
Unpin all pressure dofs in elements listed in vector.
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 ...
void fill_in_generic_dresidual_contribution_axi_nst(double *const &parameter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, unsigned flag)
Add element's contribution to the derivative of the elemental residual vector and/or Jacobian matrix ...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
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_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 and mass matrix fl...
void further_build()
Further build: pass the pointers down to the sons.
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
void dinterpolated_u_axi_nst_ddata(const Vector< double > &s, const unsigned &i, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
Compute the derivatives of the i-th component of velocity at point s with respect to all data that ca...
void fill_in_contribution_to_hessian_vector_products(Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Compute the hessian tensor vector products required to perform continuation of bifurcations analytica...
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...
Refineable version of Axisymmetric Quad Crouzeix Raviart elements (note that unlike the cartesian ver...
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 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 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 nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
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)....
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 3 (velocities)
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void further_build()
Further build for Crouzeix_Raviart interpolates the internal pressure dofs from father element: Make ...
Refineable version of Axisymmetric Quad Taylor Hood elements. (note that unlike the cartesian version...
void pin_elemental_redundant_nodal_pressure_dofs()
Unpin the proper nodal pressure dofs.
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 * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at node n.
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 rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
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...
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 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.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
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 * 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...
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...
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
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...