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