refineable_linearised_axisym_navier_stokes_elements.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 linearised axisymmetric Navier-Stokes elements
31 
32 #ifndef OOMPH_REFINEABLE_LINEARISED_AXISYM_NAVIER_STOKES_ELEMENTS_HEADER
33 #define OOMPH_REFINEABLE_LINEARISED_AXISYM_NAVIER_STOKES_ELEMENTS_HEADER
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37 #include <oomph-lib-config.h>
38 #endif
39 
40 // oomph-lib includes
41 #include "../generic/refineable_quad_element.h"
42 #include "../generic/error_estimator.h"
44 
45 namespace oomph
46 {
47 
48  //=======================================================================
49  /// \short Refineable version of the linearised axisymmetric
50  /// Navier--Stokes equations
51  //=======================================================================
54  public virtual RefineableElement,
55  public virtual ElementWithZ2ErrorEstimator
56  {
57  protected:
58 
59  /// \short Pointer to n_p-th pressure node (Default: NULL,
60  /// indicating that pressure is not based on nodal interpolation).
61  virtual Node* pressure_node_pt(const unsigned& n_p) {return NULL;}
62 
63  /// \short Unpin all pressure dofs in the element
64  virtual void unpin_elemental_pressure_dofs()=0;
65 
66  /// \short Pin unused nodal pressure dofs (empty by default, because
67  /// by default pressure dofs are not associated with nodes)
69 
70  public:
71 
72  /// \short Empty Constructor
77 
78  /// Number of 'flux' terms for Z2 error estimation
79  unsigned num_Z2_flux_terms()
80  {
81  // 3 diagonal strain rates, 3 off diagonal
82  return 6;
83  }
84 
85  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
86  /// in strain rate tensor.
88  {
89  // Specify the number of velocity dimensions
90  unsigned DIM = 3;
91 
92 #ifdef PARANOID
93  unsigned num_entries=DIM+((DIM*DIM)-DIM)/2;
94  if (flux.size()!=num_entries)
95  {
96  std::ostringstream error_message;
97  error_message << "The flux vector has the wrong number of entries, "
98  << flux.size() << ", whereas it should be "
99  << num_entries << std::endl;
100  throw OomphLibError(
101  error_message.str(),
102  "RefineableLinearisedAxisymmetricNavierStokesEquations::get_Z2_flux()",
103  OOMPH_EXCEPTION_LOCATION);
104  }
105 #endif
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  /// Further build: pass the pointers down to the sons
137  {
138  // Find the father element
140  cast_father_element_pt
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 
147  // Set the density ratio pointer
148  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
149 
150  // Set pointer to global Reynolds number
151  this->Re_pt = cast_father_element_pt->re_pt();
152 
153  // Set pointer to global Reynolds number x Strouhal number (=Womersley)
154  this->ReSt_pt = cast_father_element_pt->re_st_pt();
155 
156  // Set pointer to azimuthal mode number
158  cast_father_element_pt->azimuthal_mode_number_pt();
159 
160  // Set the ALE_is_disabled flag
161  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
162  }
163 
164  /// \short Loop over all elements in Vector (which typically contains
165  /// all the elements in a fluid mesh) and pin the nodal pressure degrees
166  /// of freedom that are not being used. Function uses
167  /// the member function
168  /// - \c RefineableLinearisedAxisymmetricNavierStokesEquations::
169  /// pin_all_nodal_pressure_dofs()
170  /// .
171  /// which is empty by default and should be implemented for
172  /// elements with nodal pressure degrees of freedom
173  /// (e.g. for refineable Taylor-Hood.)
175  const Vector<GeneralisedElement*>& element_pt)
176  {
177  // Loop over all elements to brutally pin all nodal pressure degrees of
178  // freedom
179  const unsigned n_element=element_pt.size();
180  for(unsigned e=0;e<n_element;e++)
181  {
184  }
185  }
186 
187  /// Unpin all pressure dofs in elements listed in Vector
189  element_pt)
190  {
191  // Loop over all elements to brutally unpin all nodal pressure degrees of
192  // freedom and internal pressure degrees of freedom
193  const unsigned n_element = element_pt.size();
194  for(unsigned e=0;e<n_element;e++)
195  {
197  (element_pt[e])->unpin_elemental_pressure_dofs();
198  }
199  }
200 
201 
202  private:
203 
204  /// \short Add element's contribution to the elemental residual vector
205  /// and/or Jacobian matrix
206  /// flag=1: compute both
207  /// flag=0: compute only residual vector
209  Vector<double> &residuals,
210  DenseMatrix<double> &jacobian,
211  DenseMatrix<double> &mass_matrix,
212  unsigned flag);
213 
214  }; // End of RefineableLinearisedAxisymmetricNavierStokesEquations class defn
215 
216 
217 
218 //////////////////////////////////////////////////////////////////////////////
219 //////////////////////////////////////////////////////////////////////////////
220 //////////////////////////////////////////////////////////////////////////////
221 
222 
223  //=======================================================================
224  /// \short Refineable version of linearised axisymmetric quadratic
225  /// Crouzeix-Raviart elements
226  //=======================================================================
230  public virtual RefineableQElement<2>
231  {
232  private:
233 
234  /// Unpin all the internal pressure freedoms
236  {
237  const unsigned n_pres = this->npres_linearised_axi_nst();
238  // Loop over pressure dofs and unpin
239  for(unsigned l=0;l<n_pres;l++)
240  {
241  // There are two pressure components
242  for(unsigned i=0;i<2;i++)
243  {
245  ->unpin(l);
246  }
247  }
248  }
249 
250  public:
251 
252  /// Constructor
256  RefineableQElement<2>(),
258 
259  /// Number of continuously interpolated values: 6 (velocities)
260  unsigned ncont_interpolated_values() const { return 6; }
261 
262  /// Rebuild from sons: Reconstruct pressure from the (merged) sons
263  void rebuild_from_sons(Mesh* &mesh_pt)
264  {
265  using namespace QuadTreeNames;
266 
267  // Central pressure value:
268  // -----------------------
269 
270  // Use average of the sons central pressure values
271  // Other options: Take average of the four (discontinuous)
272  // pressure values at the father's midpoint]
273 
274  Vector<double> av_press(2,0.0);
275 
276  // Loop over the sons
277  for(unsigned ison=0;ison<4;ison++)
278  {
279  // Loop over the two pressure components
280  for(unsigned i=0;i<2;i++)
281  {
282  // Add the sons midnode pressure
283  av_press[i] += quadtree_pt()->son_pt(ison)->object_pt()->
285  }
286  }
287 
288  // Use the average
289  for(unsigned i=0;i<2;i++)
290  {
292  ->set_value(0,0.25*av_press[i]);
293  }
294 
295  Vector<double> slope1(2,0.0), slope2(2,0.0);
296 
297  // Loop over pressure components
298  for(unsigned i=0;i<2;i++)
299  {
300  // Slope in s_0 direction
301  // ----------------------
302 
303  // Use average of the 2 FD approximations based on the
304  // elements central pressure values
305  // [Other options: Take average of the four
306  // pressure derivatives]
307 
308  slope1[i] = quadtree_pt()->son_pt(SE)->object_pt()->
310  quadtree_pt()->son_pt(SW)->object_pt()->
312 
313  slope2[i] = quadtree_pt()->son_pt(NE)->object_pt()->
315  quadtree_pt()->son_pt(NW)->object_pt()->
317 
318  // Use the average
320  set_value(1,0.5*(slope1[i]+slope2[i]));
321 
322  // Slope in s_1 direction
323  // ----------------------
324 
325  // Use average of the 2 FD approximations based on the
326  // elements central pressure values
327  // [Other options: Take average of the four
328  // pressure derivatives]
329 
330  slope1[i] = quadtree_pt()->son_pt(NE)->object_pt()->
332  quadtree_pt()->son_pt(SE)->object_pt()->
334 
335  slope2[i] = quadtree_pt()->son_pt(NW)->object_pt()->
337  quadtree_pt()->son_pt(SW)->object_pt()->
339 
340  // Use the average
342  set_value(2,0.5*(slope1[i]+slope2[i]));
343  } // End of loop over pressure components
344  }
345 
346  /// \short Order of recovery shape functions for Z2 error estimation:
347  /// Same order as shape functions.
348  unsigned nrecovery_order() { return 2; }
349 
350  /// \short Number of vertex nodes in the element
351  unsigned nvertex_node() const
353 
354  /// \short Pointer to the j-th vertex node in the element
355  Node* vertex_node_pt(const unsigned& j) const
356  {
358  }
359 
360  /// \short Get the function value u in Vector.
361  /// Note: Given the generality of the interface (this function
362  /// is usually called from black-box documentation or interpolation
363  /// routines), the values Vector sets its own size in here.
365  Vector<double>& values)
366  {
367  // Set the velocity dimension of the element
368  const unsigned DIM=3;
369 
370  // Determine size of values Vector: U^C, U^S, W^C, W^S, V^C, V^S
371  const unsigned n_values = 2*DIM;
372 
373  // Set size of values Vector and initialise to zero
374  values.resize(n_values,0.0);
375 
376  // Calculate velocities: values[0],...
377  for(unsigned i=0;i<n_values;i++)
378  {
380  }
381  }
382 
383  /// \short Get all function values [U^C,U^S,...,P^S] at previous timestep t
384  /// (t=0: present; t>0: previous timestep).
385  /// \n
386  /// Note: Given the generality of the interface (this function is
387  /// usually called from black-box documentation or interpolation
388  /// routines), the values Vector sets its own size in here.
389  /// \n
390  /// Note: No pressure history is kept, so pressure is always
391  /// the current value.
392  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
393  Vector<double>& values)
394  {
395  const unsigned DIM=3;
396 
397  // Set size of Vector: U^C, U^S, W^C, W^S, V^C, V^S
398  values.resize(2*DIM);
399 
400  // Initialise to zero
401  for(unsigned i=0;i<2*DIM;i++) { values[i]=0.0; }
402 
403  // Determine number of nodes in element
404  const unsigned n_node = nnode();
405 
406  // Shape functions
407  Shape psif(n_node);
408  shape(s,psif);
409 
410  // Calculate velocities: values[0],...
411  for(unsigned i=0;i<(2*DIM);i++)
412  {
413  // Get the local index at which the i-th velocity is stored
414  const unsigned u_local_index = u_index_linearised_axi_nst(i);
415  for(unsigned l=0;l<n_node;l++)
416  {
417  values[i] += nodal_value(t,l,u_local_index)*psif[l];
418  }
419  }
420  }
421 
422  /// \short Perform additional hanging node procedures for variables
423  /// that are not interpolated by all nodes. Empty
425 
426  /// Further build for Crouzeix_Raviart interpolates the internal
427  /// pressure dofs from father element: Make sure pressure values and
428  /// dp/ds agree between fathers and sons at the midpoints of the son
429  /// elements.
431  {
432  // Call the generic further build
434 
435  using namespace QuadTreeNames;
436 
437  // What type of son am I? Ask my quadtree representation...
438  int son_type=quadtree_pt()->son_type();
439 
440  // Pointer to my father (in element impersonation)
441  RefineableElement* father_el_pt=
443 
444  Vector<double> s_father(2);
445 
446  // Son midpoint is located at the following coordinates in father element:
447 
448  // South west son
449  if(son_type==SW)
450  {
451  s_father[0]=-0.5;
452  s_father[1]=-0.5;
453  }
454  // South east son
455  else if(son_type==SE)
456  {
457  s_father[0]= 0.5;
458  s_father[1]=-0.5;
459  }
460  // North east son
461  else if(son_type==NE)
462  {
463  s_father[0]=0.5;
464  s_father[1]=0.5;
465  }
466 
467  // North west son
468  else if(son_type==NW)
469  {
470  s_father[0]=-0.5;
471  s_father[1]= 0.5;
472  }
473 
474  // Pressure values in father element
475  // ---------------------------------
476 
477  // Find pointer to father element
479  cast_father_el_pt=
481  (father_el_pt);
482 
483  // Set up storage for pressure in father element
484  Vector<double> press(2,0.0);
485 
486  // Loop over pressure components
487  for(unsigned i=0;i<2;i++)
488  {
489  // Get pressure from father element
490  press[i]
491  = cast_father_el_pt->interpolated_p_linearised_axi_nst(s_father,i);
492 
493  // Pressure value gets copied straight into internal dof:
495  set_value(0,press[i]);
496 
497  // The slopes get copied from father
499  set_value(1,0.5*cast_father_el_pt
501  ->value(1));
502 
504  set_value(2,0.5*cast_father_el_pt
506  ->value(2));
507  } // End of loop over pressure components
508  }
509 
510  }; // End of RefineableLinearisedAxisymmetricQCrouzeixRaviartElement defn
511 
512 
513 
514 ///////////////////////////////////////////////////////////////////////////
515 ///////////////////////////////////////////////////////////////////////////
516 ///////////////////////////////////////////////////////////////////////////
517 
518 
519 
520  //=======================================================================
521  /// \short Face geometry of the refineable linearised axisym
522  /// Crouzeix-Raviart elements
523  //=======================================================================
524  template<>
526  public virtual FaceGeometry<LinearisedAxisymmetricQCrouzeixRaviartElement>
527  {
528  public:
531  };
532 
533 
534 
535  //=======================================================================
536  /// \short Face geometry of face geometric of the refineable linearised
537  /// axisym Crouzeix-Raviart elements
538  //=======================================================================
539  template<>
542  public virtual FaceGeometry
543  <FaceGeometry<LinearisedAxisymmetricQCrouzeixRaviartElement> >
544  {
545  public:
549  };
550 
551 
552 
553 //////////////////////////////////////////////////////////////////////////////
554 //////////////////////////////////////////////////////////////////////////////
555 //////////////////////////////////////////////////////////////////////////////
556 
557 
558 
559  //=======================================================================
560  /// \short Refineable version of linearised axisymmetric quadratic
561  /// Taylor-Hood elements
562  //=======================================================================
566  public virtual RefineableQElement<2>
567  {
568  private:
569 
570  /// Pointer to n_p-th pressure node
571  Node* pressure_node_pt(const unsigned &n_p)
572  { return this->node_pt(this->Pconv[n_p]); }
573 
574  /// Unpin all pressure dofs
576  {
577  // Determine number of nodes in element
578  const unsigned n_node = this->nnode();
579 
580  // Get nodal indeces of the two pressure components
581  Vector<int> p_index(2);
582  for(unsigned i=0;i<2;i++)
583  {
584  p_index[i] = this->p_index_linearised_axi_nst(i);
585  }
586 
587  // Loop over nodes and unpin both pressure components
588  for(unsigned n=0;n<n_node;n++)
589  {
590  for(unsigned i=0;i<2;i++) { this->node_pt(n)->unpin(p_index[i]); }
591  }
592  }
593 
594  /// Unpin the proper nodal pressure dofs
596  {
597  // Determine number of nodes in element
598  const unsigned n_node = this->nnode();
599 
600  // Get nodal indeces of the two pressure components
601  Vector<int> p_index(2);
602  for(unsigned i=0;i<2;i++)
603  {
604  p_index[i] = this->p_index_linearised_axi_nst(i);
605  }
606 
607  // Loop over all nodes and pin all the nodal pressures
608  for(unsigned n=0;n<n_node;n++)
609  {
610  for(unsigned i=0;i<2;i++) { this->node_pt(n)->pin(p_index[i]); }
611  }
612 
613  // Loop over all actual pressure nodes and unpin if they're not hanging
614  const unsigned n_pres = this->npres_linearised_axi_nst();
615  for(unsigned l=0;l<n_pres;l++)
616  {
617  Node* nod_pt = this->node_pt(this->Pconv[l]);
618  for(unsigned i=0;i<2;i++)
619  {
620  if(!nod_pt->is_hanging(p_index[i])) { nod_pt->unpin(p_index[i]); }
621  }
622  }
623  }
624 
625  public:
626 
627  /// \short Constructor:
631  RefineableQElement<2>(),
633 
634  /// \short Number of values (pinned or dofs) required at node n.
635  /// Bumped up to 8 so we don't have to worry if a hanging mid-side node
636  /// gets shared by a corner node (which has extra degrees of freedom)
637  unsigned required_nvalue(const unsigned &n) const { return 8; }
638 
639  /// \short Number of continuously interpolated values: 8
640  /// (6 velocities + 2 pressures)
641  unsigned ncont_interpolated_values() const { return 8; }
642 
643  /// Rebuild from sons: empty
644  void rebuild_from_sons(Mesh* &mesh_pt) {}
645 
646  /// \short Order of recovery shape functions for Z2 error estimation:
647  /// Same order as shape functions.
648  unsigned nrecovery_order() { return 2; }
649 
650  /// Number of vertex nodes in the element
651  unsigned nvertex_node() const
653 
654  /// Pointer to the j-th vertex node in the element
655  Node* vertex_node_pt(const unsigned& j) const
656  {
658  }
659 
660  /// \short Get the function value u in Vector.
661  /// Note: Given the generality of the interface (this function
662  /// is usually called from black-box documentation or interpolation
663  /// routines), the values Vector sets its own size in here.
665  Vector<double> &values)
666  {
667  // Set the velocity dimension of the element
668  const unsigned DIM = 3;
669 
670  // Determine size of values Vector:
671  // U^C, U^S, W^C, W^S, V^C, V^S, P^C, P^S
672  const unsigned n_values = 2*(DIM+1);
673 
674  // Set size of values Vector and initialise to zero
675  values.resize(n_values,0.0);
676 
677  // Calculate velocities: values[0],...
678  for(unsigned i=0;i<(2*DIM);i++)
679  {
681  }
682 
683  // Calculate pressure: values[DIM], values[DIM+1]
684  for(unsigned i=0;i<2;i++)
685  {
686  values[DIM+i] = interpolated_p_linearised_axi_nst(s,i);
687  }
688  }
689 
690  /// \short Get the function value u in Vector.
691  /// Note: Given the generality of the interface (this function
692  /// is usually called from black-box documentation or interpolation
693  /// routines), the values Vector sets its own size in here.
694  void get_interpolated_values(const unsigned &t,
695  const Vector<double> &s,
696  Vector<double> &values)
697  {
698  // Set the velocity dimension of the element
699  const unsigned DIM=3;
700 
701  // Set size of values Vector: U^C, U^S, W^C, W^S, V^C, V^S, P^C, P^S
702  values.resize(2*(DIM+1));
703 
704  // Initialise all entries to zero
705  for(unsigned i=0;i<2*(DIM+1);i++) { values[i]=0.0; }
706 
707  // Determine number of nodes in element
708  const unsigned n_node = nnode();
709 
710  // Shape functions
711  Shape psif(n_node);
712  shape(s,psif);
713 
714  // Calculate velocities: values[0],...
715  for(unsigned i=0;i<(2*DIM);i++)
716  {
717  // Get the nodal coordinate of the velocity
718  const unsigned u_nodal_index = u_index_linearised_axi_nst(i);
719  for(unsigned l=0;l<n_node;l++)
720  {
721  values[i] += nodal_value(t,l,u_nodal_index)*psif[l];
722  }
723  }
724 
725  // Calculate pressure: values[DIM], values[DIM+1]
726  // (no history is carried in the pressure)
727  for(unsigned i=0;i<2;i++)
728  {
729  values[DIM+i] = interpolated_p_linearised_axi_nst(s,i);
730  }
731  }
732 
733  /// \short Perform additional hanging node procedures for variables
734  /// that are not interpolated by all nodes. The two pressure components
735  /// are stored at the 6th and 7th location in each node
737  {
738  const unsigned DIM = 3;
739  for(unsigned i=0;i<2;i++)
740  {
741  this->setup_hang_for_value((2*DIM)+i);
742  }
743  }
744 
745  /// \short The velocities are isoparametric and so the "nodes"
746  /// interpolating the velocities are the geometric nodes. The
747  /// pressure "nodes" are a subset of the nodes, so when n_value==6
748  /// or 7, the n-th pressure node is returned.
749  Node* interpolating_node_pt(const unsigned &n,
750  const int &n_value)
751 
752  {
753  const int DIM = 3;
754 
755  // The only different nodes are the pressure nodes
756  if(n_value==(2*DIM) || n_value==((2*DIM)+1))
757  {
758  return this->node_pt(this->Pconv[n]);
759  }
760 
761  // The other variables are interpolated via the usual nodes
762  else { return this->node_pt(n); }
763  }
764 
765  /// \short The pressure nodes are the corner nodes, so when n_value==6
766  /// or 7, the fraction is the same as the 1d node number, 0 or 1.
767  double local_one_d_fraction_of_interpolating_node(const unsigned &n1d,
768  const unsigned &i,
769  const int &n_value)
770  {
771  const int DIM=3;
772  if(n_value==(2*DIM) || n_value==((2*DIM)+1))
773  {
774  // The pressure nodes are just located on the boundaries at 0 or 1
775  return double(n1d);
776  }
777  // Otherwise the velocity nodes are the same as the geometric ones
778  else
779  {
780  return this->local_one_d_fraction_of_node(n1d,i);
781  }
782  }
783 
784  /// \short The velocity nodes are the same as the geometric nodes.
785  /// The pressure nodes must be calculated by using the same methods
786  /// as the geometric nodes, but by recalling that there are only two
787  /// pressure nodes per edge.
789  const int &n_value)
790  {
791  const int DIM = 3;
792 
793  // If we are calculating pressure nodes
794  if(n_value==static_cast<int>(2*DIM)
795  || n_value==static_cast<int>((2*DIM)+1))
796  {
797  // Storage for the index of the pressure node
798  unsigned total_index = 0;
799  // The number of nodes along each 1d edge is 2.
800  const unsigned NNODE_1D = 2;
801  // Storage for the index along each boundary
802  // Note that it's only a 2D spatial element
803  Vector<int> index(2);
804  // Loop over the coordinates
805  for(unsigned i=0;i<2;i++)
806  {
807  // If we are at the lower limit, the index is zero
808  if(s[i]==-1.0) { index[i] = 0; }
809 
810  // If we are at the upper limit, the index is the number of nodes
811  // minus 1
812  else if(s[i] == 1.0) { index[i] = NNODE_1D-1; }
813 
814  // Otherwise, we have to calculate the index in general
815  else
816  {
817  // For uniformly spaced nodes the 0th node number would be
818  double float_index = 0.5*(1.0 + s[i])*(NNODE_1D-1);
819  index[i] = int(float_index);
820  //What is the excess. This should be safe because the
821  //taking the integer part rounds down
822  double excess = float_index - index[i];
823  //If the excess is bigger than our tolerance there is no node,
824  //return null
825  if((excess > FiniteElement::Node_location_tolerance) && (
826  (1.0 - excess) > FiniteElement::Node_location_tolerance))
827  {
828  return 0;
829  }
830  }
831  ///Construct the general pressure index from the components.
832  total_index += index[i]*
833  static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
834  static_cast<int>(i)));
835  }
836  // If we've got here we have a node, so let's return a pointer to it
837  return this->node_pt(this->Pconv[total_index]);
838  }
839  // Otherwise velocity nodes are the same as pressure nodes
840  else
841  {
842  return this->get_node_at_local_coordinate(s);
843  }
844  }
845 
846 
847  /// \short The number of 1d pressure nodes is 2, the number of 1d
848  /// velocity nodes is the same as the number of 1d geometric nodes.
849  unsigned ninterpolating_node_1d(const int &n_value)
850  {
851  const int DIM=3;
852  if(n_value==(2*DIM) || n_value==((2*DIM)+1)) { return 2; }
853  else { return this->nnode_1d(); }
854  }
855 
856  /// \short The number of pressure nodes is 4. The number of
857  /// velocity nodes is the same as the number of geometric nodes.
858  unsigned ninterpolating_node(const int &n_value)
859  {
860  const int DIM=3;
861  if(n_value==(2*DIM) || n_value==((2*DIM)+1)) { return 4; }
862  else { return this->nnode(); }
863  }
864 
865  /// \short The basis interpolating the pressure is given by pshape().
866  //// The basis interpolating the velocity is shape().
868  Shape &psi,
869  const int &n_value) const
870  {
871  const int DIM=3;
872  if(n_value==(2*DIM) || n_value==((2*DIM)+1))
873  {
874  return this->pshape_linearised_axi_nst(s,psi);
875  }
876  else { return this->shape(s,psi); }
877  }
878 
879  }; // End of RefineableLinearisedAxisymmetricQTaylorHoodElement class defn
880 
881 
882 
883 ///////////////////////////////////////////////////////////////////////////
884 ///////////////////////////////////////////////////////////////////////////
885 ///////////////////////////////////////////////////////////////////////////
886 
887 
888 
889  //=======================================================================
890  /// \short Face geometry of the refineable linearised axisym
891  /// Taylor-Hood elements
892  //=======================================================================
893  template<>
895  public virtual FaceGeometry<LinearisedAxisymmetricQTaylorHoodElement>
896  {
897  public:
899  };
900 
901 
902 
903  //=======================================================================
904  /// \short Face geometry of face geometric of the refineable linearised
905  /// axisym Taylor-Hood elements
906  //=======================================================================
907  template<>
910  public virtual FaceGeometry<FaceGeometry
911  <LinearisedAxisymmetricQTaylorHoodElement> >
912  {
913  public:
916  };
917 
918 
919 
920 } // End of namespace oomph
921 
922 #endif
void unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:386
Refineable version of linearised axisymmetric quadratic Taylor-Hood elements.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Definition: tree.h:102
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons.
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 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 pshape_linearised_axi_nst(const Vector< double > &s, Shape &psi) const
Compute the pressure shape functions at local coordinate s.
virtual unsigned u_index_linearised_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value...
cstr elem_len * i
Definition: cfortran.h:607
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...
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
Vector< unsigned > P_linearised_axi_nst_internal_index
Internal indices that indicate at which internal data the pressure values are stored. We note that there are two pressure values, corresponding to the functions P^C(r,z,t) and P^S(r,z,t) which multiply the cosine and sine terms respectively.
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
char t
Definition: cfortran.h:572
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...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 8 (6 velocities + 2 pressures)
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor. ...
double geometric_jacobian(const Vector< double > &x)
Fill in the geometric Jacobian, which in this case is r.
double local_one_d_fraction_of_interpolating_node(const unsigned &n1d, const unsigned &i, const int &n_value)
The pressure nodes are the corner nodes, so when n_value==6 or 7, the fraction is the same as the 1d ...
e
Definition: cfortran.h:575
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...
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
unsigned npres_linearised_axi_nst() const
Return number of pressure values corresponding to a single pressure component.
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 all function values [U^C,U^S,...,P^S] at previous timestep t (t=0: present; t>0: previous timeste...
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 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 ...
QuadTree * quadtree_pt()
Pointer to quadtree representation of this element.
double interpolated_p_linearised_axi_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated pressure p[i] at local coordinate s...
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
double interpolated_u_linearised_axi_nst(const Vector< double > &s, const unsigned &i) const
Return the i-th component of the FE interpolated velocity u[i] at local coordinate s...
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
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
void fill_in_generic_residual_contribution_linearised_axi_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add element's contribution to the elemental residual vector and/or Jacobian matrix flag=1: compute bo...
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
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
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)
static void unpin_all_pressure_dofs(const Vector< GeneralisedElement * > &element_pt)
Unpin all pressure dofs in elements listed in Vector.
int *& azimuthal_mode_number_pt()
Pointer to azimuthal mode number k in e^ik(theta) decomposition.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
virtual int p_index_linearised_axi_nst(const unsigned &i) const
Which nodal value represents the pressure?
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
Definition: elements.h:1334
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...
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node.
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
int son_type() const
Return son type.
Definition: tree.h:209
Refineable version of the linearised axisymmetric Navier–Stokes equations.
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at node n. Bumped up to 8 so we don't have to worry if a h...
int * Azimuthal_Mode_Number_pt
Pointer to azimuthal mode number k in e^ik(theta) decomposition.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
virtual unsigned nvertex_node() const
Definition: elements.h:2362
Refineable version of linearised axisymmetric quadratic Crouzeix-Raviart elements.
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Pin unused nodal pressure dofs (empty by default, because by default pressure dofs are not associated...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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...
Tree * father_pt() const
Return pointer to father: NULL if it's a root node.
Definition: tree.h:221
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
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
A class for elements that solve the linearised version of the unsteady Navier–Stokes equations in cyl...
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...
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
static void pin_redundant_nodal_pressures(const Vector< GeneralisedElement * > &element_pt)
Loop over all elements in Vector (which typically contains all the elements in a fluid mesh) and pin ...
virtual void 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
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 6 (velocities)
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
Node * get_interpolating_node_at_local_coordinate(const Vector< double > &s, const int &n_value)
The velocity nodes are the same as the geometric nodes. The pressure nodes must be calculated by usin...