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// Version 1.0; svn revision $LastChangedRevision: 1194 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-27 08:44:43 +0100 (Fri, 27 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for refineable axisymmetric quad Navier Stokes elements
31 #ifndef OOMPH_GENERALISED_NEWTONIAN_REFINEABLE_AXISYMMETRIC_NAVIER_STOKES_ELEMENTS_HEADER
32 #define OOMPH_GENERALISED_NEWTONIAN_REFINEABLE_AXISYMMETRIC_NAVIER_STOKES_ELEMENTS_HEADER
33 
34 // Config header generated by autoconfig
35 #ifdef HAVE_CONFIG_H
36 #include <oomph-lib-config.h>
37 #endif
38 
39 //Oomph lib includes
40 #include "../generic/refineable_quad_element.h"
41 #include "../generic/error_estimator.h"
43 
44 namespace oomph
45 {
46 
47 
48 //======================================================================
49 ///Refineable version of the Axisymmetric Navier--Stokes equations
50 //======================================================================
53 public virtual RefineableElement,
54 public virtual ElementWithZ2ErrorEstimator
55 {
56  protected:
57 
58  /// \short Pointer to n_p-th pressure node (Default: NULL,
59  /// indicating that pressure is not based on nodal interpolation).
60  virtual Node* pressure_node_pt(const unsigned& n_p) {return NULL;}
61 
62  /// \short Unpin all pressure dofs in the element
63  virtual void unpin_elemental_pressure_dofs()=0;
64 
65 /// \short Pin unused nodal pressure dofs (empty by default, because
66  /// by default pressure dofs are not associated with nodes)
68 
69  public:
70 
71  /// \short Empty Constructor
76 
77  /// Number of 'flux' terms for Z2 error estimation
78  unsigned num_Z2_flux_terms()
79  {
80  // 3 diagonal strain rates, 3 off diagonal
81  return 6;
82  }
83 
84  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
85  /// in strain rate tensor.
87  {
88  //Specify the number of velocity dimensions
89  unsigned DIM=3;
90 
91 #ifdef PARANOID
92  unsigned num_entries=DIM+((DIM*DIM)-DIM)/2;
93  if (flux.size() < num_entries)
94  {
95  std::ostringstream error_message;
96  error_message << "The flux vector has the wrong number of entries, "
97  << flux.size() << ", whereas it should be "
98  << num_entries << std::endl;
99  throw OomphLibError(
100  error_message.str(),
101  OOMPH_CURRENT_FUNCTION,
102  OOMPH_EXCEPTION_LOCATION);
103  }
104 #endif
105 
106 
107  // Get strain rate matrix
108  DenseMatrix<double> strainrate(DIM);
109  this->strain_rate(s,strainrate);
110 
111  // Pack into flux Vector
112  unsigned icount=0;
113 
114  // Start with diagonal terms
115  for(unsigned i=0;i<DIM;i++)
116  {
117  flux[icount]=strainrate(i,i);
118  icount++;
119  }
120 
121  //Off diagonals row by row
122  for(unsigned i=0;i<DIM;i++)
123  {
124  for(unsigned j=i+1;j<DIM;j++)
125  {
126  flux[icount]=strainrate(i,j);
127  icount++;
128  }
129  }
130  }
131 
132  /// Fill in the geometric Jacobian, which in this case is r
133  double geometric_jacobian(const Vector<double> &x) {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
142  = dynamic_cast<
144  (this->father_element_pt());
145 
146  //Set the viscosity ratio pointer
147  this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
148  //Set the density ratio pointer
149  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
150  //Set pointer to global Reynolds number
151  this->Re_pt = cast_father_element_pt->re_pt();
152  //Set pointer to global Reynolds number x Strouhal number (=Womersley)
153  this->ReSt_pt = cast_father_element_pt->re_st_pt();
154  //Set pointer to global Reynolds number x inverse Froude number
155  this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
156  //Set pointer to the global Reynolds number x inverse Rossby number
157  this->ReInvRo_pt = cast_father_element_pt->re_invro_pt();
158  //Set pointer to global gravity Vector
159  this->G_pt = cast_father_element_pt->g_pt();
160 
161  //Set pointer to body force function
162  this->Body_force_fct_pt =
163  cast_father_element_pt->axi_nst_body_force_fct_pt();
164 
165  //Set pointer to volumetric source function
166  this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
167 
168  //Set the ALE_is_disabled flag
169  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
170  }
171 
172  /// \short Compute the derivatives of the i-th component of
173  /// velocity at point s with respect
174  /// to all data that can affect its value. In addition, return the global
175  /// equation numbers corresponding to the data.
176  /// Overload the non-refineable version to take account of hanging node
177  /// information
179  const unsigned &i,
180  Vector<double> &du_ddata,
181  Vector<unsigned> &global_eqn_number)
182  {
183  //Find number of nodes
184  unsigned n_node = this->nnode();
185  //Local shape function
186  Shape psi(n_node);
187  //Find values of shape function at the given local coordinate
188  this->shape(s,psi);
189 
190  //Find the index at which the velocity component is stored
191  const unsigned u_nodal_index = this->u_index_axi_nst(i);
192 
193  //Storage for hang info pointer
194  HangInfo* hang_info_pt=0;
195  //Storage for global equation
196  int global_eqn = 0;
197 
198  //Find the number of dofs associated with interpolated u
199  unsigned n_u_dof=0;
200  for(unsigned l=0;l<n_node;l++)
201  {
202  unsigned n_master = 1;
203 
204  //Local bool (is the node hanging)
205  bool is_node_hanging = this->node_pt(l)->is_hanging();
206 
207  //If the node is hanging, get the number of master nodes
208  if(is_node_hanging)
209  {
210  hang_info_pt = this->node_pt(l)->hanging_pt();
211  n_master = hang_info_pt->nmaster();
212  }
213  //Otherwise there is just one master node, the node itself
214  else
215  {
216  n_master = 1;
217  }
218 
219  //Loop over the master nodes
220  for(unsigned m=0;m<n_master;m++)
221  {
222  //Get the equation number
223  if(is_node_hanging)
224  {
225  //Get the equation number from the master node
226  global_eqn = hang_info_pt->master_node_pt(m)->
227  eqn_number(u_nodal_index);
228  }
229  else
230  {
231  // Global equation number
232  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
233  }
234 
235  //If it's positive add to the count
236  if(global_eqn >= 0) {++n_u_dof;}
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 = hang_info_pt->master_node_pt(m)->
287  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 
310  /// \short 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  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<
331  }
332  }
333 
334  /// \short Unpin all pressure dofs in elements listed in vector.
336  element_pt)
337  {
338  // Loop over all elements to brutally unpin all nodal pressure degrees of
339  // freedom and internal pressure degrees of freedom
340  unsigned n_element = element_pt.size();
341  for(unsigned e=0;e<n_element;e++)
342  {
343  dynamic_cast<
345  (element_pt[e])->unpin_elemental_pressure_dofs();
346  }
347  }
348 
349  /// \short Compute derivatives of elemental residual vector with respect to
350  /// nodal coordinates. This function computes these terms analytically and
351  /// overwrites the default implementation in the FiniteElement base class.
352  /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
354  dresidual_dnodal_coordinates);
355 
356  private:
357 
358  /// \short Add element's contribution to the elemental residual vector
359  /// and/or Jacobian matrix and mass matrix
360  /// flag=2: compute all
361  /// flag=1: compute both residual and Jacobian
362  /// flag=0: compute only residual vector
364  Vector<double> &residuals,
365  DenseMatrix<double> &jacobian,
366  DenseMatrix<double> &mass_matrix,
367  unsigned flag);
368 
369  /// \short Add element's contribution to the derivative of the
370  /// elemental residual vector
371  /// and/or Jacobian matrix and/or mass matrix
372  /// flag=2: compute all
373  /// flag=1: compute both residual and Jacobian
374  /// flag=0: compute only residual vector
376  double* const &parameter_pt,
377  Vector<double> &dres_dparam,
378  DenseMatrix<double> &djac_dparam,
379  DenseMatrix<double> &dmass_matrix_dparam,
380  unsigned flag)
381  {
382  throw OomphLibError("Not yet implemented\n",
383  OOMPH_CURRENT_FUNCTION,
384  OOMPH_EXCEPTION_LOCATION);
385  }
386 
387  /// \short Compute the hessian tensor vector products required to
388  /// perform continuation of bifurcations analytically
390  Vector<double> const &Y,
391  DenseMatrix<double> const &C,
392  DenseMatrix<double> &product)
393  {
394  throw OomphLibError(
395  "Not yet implemented\n",
396  OOMPH_CURRENT_FUNCTION,
397  OOMPH_EXCEPTION_LOCATION);
398  }
399 
400 
401 
402 };
403 
404 //======================================================================
405 /// Refineable version of Axisymmetric Quad Taylor Hood elements.
406 /// (note that unlike the cartesian version this is not scale-able
407 /// to higher dimensions!)
408 //======================================================================
412 public virtual RefineableQElement<2>
413 {
414  private:
415 
416  /// \short Pointer to n_p-th pressure node
417  Node* pressure_node_pt(const unsigned &n_p)
418  {return this->node_pt(this->Pconv[n_p]);}
419 
420  /// Unpin all pressure dofs
422  {
423  unsigned n_node = this->nnode();
424  int p_index = this->p_nodal_index_axi_nst();
425  // loop over nodes
426  for(unsigned n=0;n<n_node;n++) {this->node_pt(n)->unpin(p_index);}
427  }
428 
429  /// Unpin the proper nodal pressure dofs
431  {
432  //Loop over all nodes
433  unsigned n_node = this->nnode();
434  int p_index = this->p_nodal_index_axi_nst();
435  // loop over all nodes and pin all the nodal pressures
436  for(unsigned n=0;n<n_node;n++) {this->node_pt(n)->pin(p_index);}
437 
438  // Loop over all actual pressure nodes and unpin if they're not hanging
439  unsigned n_pres = this->npres_axi_nst();
440  for(unsigned l=0;l<n_pres;l++)
441  {
442  Node* nod_pt = this->node_pt(this->Pconv[l]);
443  if (!nod_pt->is_hanging(p_index)) {nod_pt->unpin(p_index);}
444  }
445  }
446 
447  public:
448 
449  /// \short Constructor:
453  RefineableQElement<2>(),
455 
456  /// Number of values (pinned or dofs) required at node n.
457  // Bumped up to 4 so we don't have to worry if a hanging mid-side node
458  // gets shared by a corner node (which has extra degrees of freedom)
459  unsigned required_nvalue(const unsigned &n) const {return 4;}
460 
461  /// Number of continuously interpolated values: 4 (3 velocities + 1 pressure)
462  unsigned ncont_interpolated_values() const {return 4;}
463 
464  /// Rebuild from sons: empty
465  void rebuild_from_sons(Mesh* &mesh_pt) {}
466 
467  /// \short Order of recovery shape functions for Z2 error estimation:
468  /// Same order as shape functions.
469  unsigned nrecovery_order() {return 2;}
470 
471  /// \short Number of vertex nodes in the element
472  unsigned nvertex_node() const
474 
475  /// \short Pointer to the j-th vertex node in the element
476  Node* vertex_node_pt(const unsigned& j) const
477  {
479  }
480 
481 /// \short Get the function value u in Vector.
482 /// Note: Given the generality of the interface (this function
483 /// is usually called from black-box documentation or interpolation routines),
484 /// the values Vector sets its own size in here.
486  {
487  //Set the velocity dimension of the element
488  unsigned DIM=3;
489 
490  // Set size of values Vector: u,w,v,p and initialise to zero
491  values.resize(DIM+1,0.0);
492 
493  //Calculate velocities: values[0],...
494  for(unsigned i=0;i<DIM;i++) {values[i] = interpolated_u_axi_nst(s,i);}
495 
496  //Calculate pressure: values[DIM]
497  values[DIM] = interpolated_p_axi_nst(s);
498  }
499 
500 /// \short Get the function value u in Vector.
501 /// Note: Given the generality of the interface (this function
502 /// is usually called from black-box documentation or interpolation routines),
503 /// the values Vector sets its own size in here.
504  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
505  Vector<double>& values)
506  {
507  unsigned DIM=3;
508 
509  // Set size of Vector: u,w,v,p
510  values.resize(DIM+1);
511 
512  // Initialise
513  for(unsigned i=0;i<DIM+1;i++) {values[i]=0.0;}
514 
515  //Find out how many nodes there are
516  unsigned n_node = nnode();
517 
518  // Shape functions
519  Shape psif(n_node);
520  shape(s,psif);
521 
522  //Calculate velocities: values[0],...
523  for(unsigned i=0;i<DIM;i++)
524  {
525  //Get the nodal coordinate of the velocity
526  unsigned u_nodal_index = u_index_axi_nst(i);
527  for(unsigned l=0;l<n_node;l++)
528  {
529  values[i] += nodal_value(t,l,u_nodal_index)*psif[l];
530  }
531  }
532 
533  //Calculate pressure: values[DIM]
534  //(no history is carried in the pressure)
535  values[DIM] = interpolated_p_axi_nst(s);
536  }
537 
538  /// \short Perform additional hanging node procedures for variables
539  /// that are not interpolated by all nodes. The pressures are stored
540  /// at the 3rd location in each node
542  {
543  int DIM=3;
544  this->setup_hang_for_value(DIM);
545  }
546 
547  /// \short The velocities are isoparametric and so the "nodes" interpolating
548  /// the velocities are the geometric nodes. The pressure "nodes" are a
549  /// subset of the nodes, so when n_value==DIM, the n-th pressure
550  /// node is returned.
551  Node* interpolating_node_pt(const unsigned &n,
552  const int &n_value)
553 
554  {
555  int DIM=3;
556  //The only different nodes are the pressure nodes
557  if(n_value==DIM) {return this->node_pt(this->Pconv[n]);}
558  //The other variables are interpolated via the usual nodes
559  else {return this->node_pt(n);}
560  }
561 
562  /// \short The pressure nodes are the corner nodes, so when n_value==DIM,
563  /// the fraction is the same as the 1d node number, 0 or 1.
564  double local_one_d_fraction_of_interpolating_node(const unsigned &n1d,
565  const unsigned &i,
566  const int &n_value)
567  {
568  int DIM=3;
569  if(n_value==DIM)
570  {
571  //The pressure nodes are just located on the boundaries at 0 or 1
572  return double(n1d);
573  }
574  //Otherwise the velocity nodes are the same as the geometric ones
575  else
576  {
577  return this->local_one_d_fraction_of_node(n1d,i);
578  }
579  }
580 
581  /// \short The velocity nodes are the same as the geometric nodes. The
582  /// pressure nodes must be calculated by using the same methods as
583  /// the geometric nodes, but by recalling that there are only two pressure
584  /// nodes per edge.
586  const int &n_value)
587  {
588  int DIM=3;
589  //If we are calculating pressure nodes
590  if(n_value==static_cast<int>(DIM))
591  {
592  //Storage for the index of the pressure node
593  unsigned total_index=0;
594  //The number of nodes along each 1d edge is 2.
595  unsigned NNODE_1D = 2;
596  //Storage for the index along each boundary
597  //Note that it's only a 2D spatial element
598  Vector<int> index(2);
599  //Loop over the coordinates
600  for(unsigned i=0;i<2;i++)
601  {
602  //If we are at the lower limit, the index is zero
603  if(s[i]==-1.0)
604  {
605  index[i] = 0;
606  }
607  //If we are at the upper limit, the index is the number of nodes minus 1
608  else if(s[i] == 1.0)
609  {
610  index[i] = NNODE_1D-1;
611  }
612  //Otherwise, we have to calculate the index in general
613  else
614  {
615  //For uniformly spaced nodes the 0th node number would be
616  double float_index = 0.5*(1.0 + s[i])*(NNODE_1D-1);
617  index[i] = int(float_index);
618  //What is the excess. This should be safe because the
619  //taking the integer part rounds down
620  double excess = float_index - index[i];
621  //If the excess is bigger than our tolerance there is no node,
622  //return null
623  if((excess > FiniteElement::Node_location_tolerance) && (
624  (1.0 - excess) > FiniteElement::Node_location_tolerance))
625  {
626  return 0;
627  }
628  }
629  ///Construct the general pressure index from the components.
630  total_index += index[i]*
631  static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
632  static_cast<int>(i)));
633  }
634  //If we've got here we have a node, so let's return a pointer to it
635  return this->node_pt(this->Pconv[total_index]);
636  }
637  //Otherwise velocity nodes are the same as pressure nodes
638  else
639  {
640  return this->get_node_at_local_coordinate(s);
641  }
642  }
643 
644 
645  /// \short The number of 1d pressure nodes is 2, the number of 1d velocity
646  /// nodes is the same as the number of 1d geometric nodes.
647  unsigned ninterpolating_node_1d(const int &n_value)
648  {
649  int DIM=3;
650  if(n_value==DIM) {return 2;}
651  else {return this->nnode_1d();}
652  }
653 
654  /// \short The number of pressure nodes is 4. The number of
655  /// velocity nodes is the same as the number of geometric nodes.
656  unsigned ninterpolating_node(const int &n_value)
657  {
658  int DIM=3;
659  if(n_value==DIM) {return 4;}
660  else {return this->nnode();}
661  }
662 
663  /// \short The basis interpolating the pressure is given by pshape().
664  //// The basis interpolating the velocity is shape().
666  Shape &psi,
667  const int &n_value) const
668  {
669  int DIM=3;
670  if(n_value==DIM) {return this->pshape_axi_nst(s,psi);}
671  else {return this->shape(s,psi);}
672  }
673 
674 };
675 
676 
677 ///////////////////////////////////////////////////////////////////////////
678 ///////////////////////////////////////////////////////////////////////////
679 ///////////////////////////////////////////////////////////////////////////
680 
681 //=======================================================================
682 /// Face geometry of the RefineableQuadQTaylorHoodElements
683 //=======================================================================
684 template<>
687 public virtual FaceGeometry<GeneralisedNewtonianAxisymmetricQTaylorHoodElement>
688 {
689  public:
692 };
693 
694 
695 //=======================================================================
696 /// Face geometry of the RefineableQuadQTaylorHoodElements
697 //=======================================================================
698 template<>
701 public virtual FaceGeometry<FaceGeometry<
702 GeneralisedNewtonianAxisymmetricQTaylorHoodElement> >
703 {
704  public:
708 };
709 
710 
711 //======================================================================
712 /// Refineable version of Axisymmetric Quad Crouzeix Raviart elements
713 /// (note that unlike the cartesian version this is not scale-able
714 /// to higher dimensions!)
715 //======================================================================
719 public virtual RefineableQElement<2>
720 {
721  private:
722 
723  ///Unpin all the internal pressure freedoms
725  {
726  unsigned n_pres = this->npres_axi_nst();
727  // loop over pressure dofs and unpin
728  for(unsigned l=0;l<n_pres;l++)
730  }
731 
732  public:
733 
734 /// \short Constructor:
738  RefineableQElement<2>(),
740 
741  /// Number of continuously interpolated values: 3 (velocities)
742  unsigned ncont_interpolated_values() const {return 3;}
743 
744  /// Rebuild from sons: Reconstruct pressure from the (merged) sons
745  void rebuild_from_sons(Mesh* &mesh_pt)
746  {
747  using namespace QuadTreeNames;
748 
749  //Central pressure value:
750  //-----------------------
751 
752  // Use average of the sons central pressure values
753  // Other options: Take average of the four (discontinuous)
754  // pressure values at the father's midpoint]
755 
756  double av_press=0.0;
757 
758  // Loop over the sons
759  for (unsigned ison=0;ison<4;ison++)
760  {
761  // Add the sons midnode pressure
762  av_press += quadtree_pt()->son_pt(ison)->object_pt()->
764  }
765 
766  // Use the average
768 
769 
770  //Slope in s_0 direction
771  //----------------------
772 
773  // Use average of the 2 FD approximations based on the
774  // elements central pressure values
775  // [Other options: Take average of the four
776  // pressure derivatives]
777 
778  double slope1= quadtree_pt()->son_pt(SE)->object_pt()->
780  quadtree_pt()->son_pt(SW)->object_pt()->
782 
783  double slope2 = quadtree_pt()->son_pt(NE)->object_pt()->
785  quadtree_pt()->son_pt(NW)->object_pt()->
787 
788 
789  // Use the average
791  set_value(1,0.5*(slope1+slope2));
792 
793 
794  //Slope in s_1 direction
795  //----------------------
796 
797  // Use average of the 2 FD approximations based on the
798  // elements central pressure values
799  // [Other options: Take average of the four
800  // pressure derivatives]
801 
802  slope1 = quadtree_pt()->son_pt(NE)->object_pt()->
804  quadtree_pt()->son_pt(SE)->object_pt()->
806 
807  slope2 = quadtree_pt()->son_pt(NW)->object_pt()->
809  quadtree_pt()->son_pt(SW)->object_pt()->
811 
812 
813  // Use the average
815  set_value(2,0.5*(slope1+slope2));
816 
817 
818  }
819 
820 /// \short Order of recovery shape functions for Z2 error estimation:
821 /// Same order as shape functions.
822  unsigned nrecovery_order()
823  {
824  return 2;
825  }
826 
827  /// \short Number of vertex nodes in the element
828  unsigned nvertex_node() const
830  nvertex_node();}
831 
832  /// \short Pointer to the j-th vertex node in the element
833  Node* vertex_node_pt(const unsigned& j) const
835  vertex_node_pt(j);}
836 
837 /// \short Get the function value u in Vector.
838 /// Note: Given the generality of the interface (this function
839 /// is usually called from black-box documentation or interpolation routines),
840 /// the values Vector sets its own size in here.
842  {
843  unsigned DIM=3;
844 
845  // Set size of Vector: u,w,v and initialise to zero
846  values.resize(DIM,0.0);
847 
848  //Calculate velocities: values[0],...
849  for(unsigned i=0;i<DIM;i++) {values[i] = interpolated_u_axi_nst(s,i);}
850  }
851 
852  /// \short Get all function values [u,v..,p] at previous timestep t
853  /// (t=0: present; t>0: previous timestep).
854  /// \n
855  /// Note: Given the generality of the interface (this function is
856  /// usually called from black-box documentation or interpolation routines),
857  /// the values Vector sets its own size in here.
858  /// \n
859  /// Note: No pressure history is kept, so pressure is always
860  /// the current value.
861  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
862  Vector<double>& values)
863  {
864  unsigned DIM=3;
865 
866  // Set size of Vector: u,w,v
867  values.resize(DIM);
868 
869  // Initialise
870  for(unsigned i=0;i<DIM;i++) {values[i]=0.0;}
871 
872  //Find out how many nodes there are
873  unsigned n_node = nnode();
874 
875  // Shape functions
876  Shape psif(n_node);
877  shape(s,psif);
878 
879  //Calculate velocities: values[0],...
880  for(unsigned i=0;i<DIM;i++)
881  {
882  //Get the local index at which the i-th velocity is stored
883  unsigned u_local_index = u_index_axi_nst(i);
884  for(unsigned l=0;l<n_node;l++)
885  {
886  values[i] += nodal_value(t,l,u_local_index)*psif[l];
887  }
888  }
889  }
890 
891  /// \short Perform additional hanging node procedures for variables
892  /// that are not interpolated by all nodes. Empty
894 
895  /// Further build for Crouzeix_Raviart interpolates the internal
896  /// pressure dofs from father element: Make sure pressure values and
897  /// dp/ds agree between fathers and sons at the midpoints of the son
898  /// elements.
900  {
901  //Call the generic further build
904 
905  using namespace QuadTreeNames;
906 
907  // What type of son am I? Ask my quadtree representation...
908  int son_type=quadtree_pt()->son_type();
909 
910  // Pointer to my father (in element impersonation)
911  RefineableElement* father_el_pt=
913 
914  Vector<double> s_father(2);
915 
916  // Son midpoint is located at the following coordinates in father element:
917 
918  // South west son
919  if (son_type==SW)
920  {
921  s_father[0]=-0.5;
922  s_father[1]=-0.5;
923  }
924  // South east son
925  else if (son_type==SE)
926  {
927  s_father[0]= 0.5;
928  s_father[1]=-0.5;
929  }
930  // North east son
931  else if (son_type==NE)
932  {
933  s_father[0]=0.5;
934  s_father[1]=0.5;
935  }
936 
937  // North west son
938  else if (son_type==NW)
939  {
940  s_father[0]=-0.5;
941  s_father[1]= 0.5;
942  }
943 
944  // Pressure value in father element
946  cast_father_el_pt=
947  dynamic_cast<
949  (father_el_pt);
950  double press=cast_father_el_pt->interpolated_p_axi_nst(s_father);
951 
952  // Pressure value gets copied straight into internal dof:
954 
955 
956  // The slopes get copied from father
958  set_value(1,0.5*cast_father_el_pt
960  ->value(1));
961 
963  set_value(2,0.5*cast_father_el_pt
965  ->value(2));
966 }
967 
968 };
969 
970 
971 
972 //=======================================================================
973 /// Face geometry of the RefineableQuadQCrouzeixRaviartElements
974 //=======================================================================
975 template<>
978 public virtual FaceGeometry<
979 GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement>
980 {
981  public:
984 };
985 
986 
987 //=======================================================================
988 /// Face geometry of the RefineableQuadQCrouzeixRaviartElements
989 //=======================================================================
990 template<>
992  FaceGeometry<
994 public virtual FaceGeometry<FaceGeometry<
995 GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement> >
996 {
997  public:
1001 };
1002 
1003 
1004 }
1005 
1006 #endif
void unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:386
static void unpin_all_pressure_dofs(const Vector< GeneralisedElement * > &element_pt)
Unpin all pressure dofs in elements listed in vector.
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number) ...
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 4 (3 velocities + 1 pressure)
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Definition: tree.h:102
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 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...
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. This function computes these terms analytically and overwrites the default implementation in the FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 3 (velocities)
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
cstr elem_len * i
Definition: cfortran.h:607
virtual unsigned u_index_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number...
Definition: elements.h:709
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 ...
double * ReInvRo_pt
Pointer to global Reynolds number x inverse Rossby number (used when in a rotating frame) ...
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:2458
double(* Source_fct_pt)(const double &time, const Vector< double > &x)
Pointer to volumetric source function.
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...
char t
Definition: cfortran.h:572
void(* Body_force_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to body force function.
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...
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
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...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
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...
e
Definition: cfortran.h:575
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...
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:2139
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
double geometric_jacobian(const Vector< double > &x)
Fill in the geometric Jacobian, which in this case is r.
double(*&)(const double &time, const Vector< double > &x) source_fct_pt()
Access function for the source-function pointer.
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number) ...
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at node n.
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 pshape_axi_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
QuadTree * quadtree_pt()
Pointer to quadtree representation of this element.
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
A Rank 3 Tensor class.
Definition: matrices.h:1337
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Pin unused nodal pressure dofs (empty by default, because.
virtual int p_nodal_index_axi_nst() const
Which nodal value represents the pressure?
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
static char t char * s
Definition: cfortran.h:572
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:267
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
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...
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:3804
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 interpolated_u_axi_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
double interpolated_p_axi_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
Definition: elements.h:1334
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor. ...
Class that contains data for hanging nodes.
Definition: nodes.h:684
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
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:2371
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
int son_type() const
Return son type.
Definition: tree.h:209
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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 ...
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 *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
virtual unsigned nvertex_node() const
Definition: elements.h:2362
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...
Tree * father_pt() const
Return pointer to father: NULL if it's a root node.
Definition: tree.h:221
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:1810
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:291
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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...
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:114
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
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...
A general mesh class.
Definition: mesh.h:74
void(*&)(const double &time, const Vector< double > &x, Vector< double > &f) axi_nst_body_force_fct_pt()
Access function for the body-force pointer.
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...