refineable_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: 1182 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-20 16:50:20 +0100 (Fri, 20 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 2D quad Navier Stokes elements
31 
32 #ifndef OOMPH_REFINEABLE_NAVIER_STOKES_ELEMENTS_HEADER
33 #define OOMPH_REFINEABLE_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 headers
41 #include "../generic/refineable_quad_element.h"
42 #include "../generic/refineable_brick_element.h"
43 #include "../generic/hp_refineable_elements.h"
44 #include "../generic/error_estimator.h"
45 #include "navier_stokes_elements.h"
46 
47 namespace oomph
48 {
49 
50 ///////////////////////////////////////////////////////////////////////
51 ///////////////////////////////////////////////////////////////////////
52 ///////////////////////////////////////////////////////////////////////
53 
54 
55 //======================================================================
56 /// A class for elements that allow the imposition of Robin boundary
57 /// conditions for the pressure advection diffusion problem in the
58 /// Fp preconditioner.
59 /// The geometrical information can be read from the FaceGeometery<ELEMENT>
60 /// class and and thus, we can be generic enough without the need to have
61 /// a separate equations class
62 //======================================================================
63 template <class ELEMENT>
66 {
67 
68  public:
69 
70  ///Constructor, which takes a "bulk" element and the value of the index
71  ///and its limit
73  const int &face_index) :
74  FaceGeometry<ELEMENT>(), FaceElement(),
75  FpPressureAdvDiffRobinBCElement<ELEMENT>(element_pt,face_index,true)
76  {}
77 
78  /// \short This function returns the residuals for the
79  /// traction function. flag=1 (or 0): do (or don't) compute the
80  /// Jacobian as well.
82  Vector<double> &residuals,
83  DenseMatrix<double> &jacobian,
84  unsigned flag);
85 
86 
87 };
88 
89 
90 
91 
92 
93 
94 
95 
96 
97 
98 //============================================================================
99 /// Get residuals and Jacobian of Robin boundary conditions in pressure
100 /// advection diffusion problem in Fp preconditoner. Refineable version.
101 //============================================================================
102 template<class ELEMENT>
105  Vector<double> &residuals,
106  DenseMatrix<double> &jacobian,
107  unsigned flag)
108 {
109  //Pointers to hang info objects
110  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
111 
112  //Storage for local coordinates in FaceElement and associted bulk element
113  unsigned my_dim=this->dim();
114  Vector<double> s(my_dim);
115  Vector<double> s_bulk(my_dim+1);
116 
117  // Storage for outer unit normal
118  Vector<double> unit_normal(my_dim+1);
119 
120  //Storage for veloc in bulk element
121  Vector<double> veloc(my_dim+1);
122 
123  //Set the value of n_intpt
124  unsigned n_intpt = this->integral_pt()->nweight();
125 
126  //Integers to store local equation numbers
127  int local_eqn=0;
128  int local_unknown=0;
129 
130  // Get cast bulk element
131  ELEMENT* bulk_el_pt=dynamic_cast<ELEMENT*>(this->bulk_element_pt());
132 
133  //Find out how many pressure dofs there are in the bulk element
134  unsigned n_pres = bulk_el_pt->npres_nst();
135 
136 
137  // Which nodal value represents the pressure? (Negative if pressure
138  // is not based on nodal interpolation).
139  int p_index = bulk_el_pt->p_nodal_index_nst();
140 
141  // Local array of booleans that are true if the l-th pressure value is
142  // hanging (avoid repeated virtual function calls)
143  bool pressure_dof_is_hanging[n_pres];
144  //If the pressure is stored at a node
145  if(p_index >= 0)
146  {
147  //Read out whether the pressure is hanging
148  for(unsigned l=0;l<n_pres;++l)
149  {
150  pressure_dof_is_hanging[l] =
151  bulk_el_pt->pressure_node_pt(l)->is_hanging(p_index);
152  }
153  }
154  //Otherwise the pressure is not stored at a node and so cannot hang
155  else
156  {
157  // pressure advection diffusion doesn't work for this one!
158  throw OomphLibError(
159  "Pressure advection diffusion does not work in this case\n",
160  OOMPH_CURRENT_FUNCTION,
161  OOMPH_EXCEPTION_LOCATION);
162 
163  for(unsigned l=0;l<n_pres;++l)
164  {pressure_dof_is_hanging[l] = false;}
165  }
166 
167  // Get the Reynolds number from the bulk element
168  double re = bulk_el_pt->re();
169 
170  //Set up memory for pressure shape and test functions
171  Shape psip(n_pres), testp(n_pres);
172 
173  //Loop over the integration points
174  for(unsigned ipt=0;ipt<n_intpt;ipt++)
175  {
176  //Get the integral weight
177  double w = this->integral_pt()->weight(ipt);
178 
179  //Assign values of local coordinate in FaceElement
180  for(unsigned i=0;i<my_dim;i++) s[i] = this->integral_pt()->knot(ipt,i);
181 
182  // Find corresponding coordinate in the the bulk element
183  s_bulk=this->local_coordinate_in_bulk(s);
184 
185  /// Get outer unit normal
186  this->outer_unit_normal(ipt,unit_normal);
187 
188  // Get velocity in bulk element
189  bulk_el_pt->interpolated_u_nst(s_bulk,veloc);
190 
191  // Get normal component of veloc
192  double flux=0.0;
193  for (unsigned i=0;i<my_dim+1;i++)
194  {
195  flux+=veloc[i]*unit_normal[i];
196  }
197 
198  // Modify bc: If outflow (flux>0) apply Neumann condition instead
199  if (flux>0.0) flux=0.0;
200 
201  // Get pressure
202  double interpolated_press=bulk_el_pt->interpolated_p_nst(s_bulk);
203 
204  //Call the pressure shape and test functions in bulk element
205  bulk_el_pt->pshape_nst(s_bulk,psip,testp);
206 
207  //Find the Jacobian of the mapping within the FaceElement
208  double J = this->J_eulerian(s);
209 
210  //Premultiply the weights and the Jacobian
211  double W = w*J;
212 
213 
214  //Number of master nodes and storage for the weight of the shape function
215  unsigned n_master=1; double hang_weight=1.0;
216 
217  //Loop over the pressure shape functions
218  for(unsigned l=0;l<n_pres;l++)
219  {
220  //If the pressure dof is hanging
221  if(pressure_dof_is_hanging[l])
222  {
223  // Pressure dof is hanging so it must be nodal-based
224  // Get the hang info object
225  hang_info_pt = bulk_el_pt->pressure_node_pt(l)->hanging_pt(p_index);
226 
227  //Get the number of master nodes from the pressure node
228  n_master = hang_info_pt->nmaster();
229  }
230  //Otherwise the node is its own master
231  else
232  {
233  n_master = 1;
234  }
235 
236  //Loop over the master nodes
237  for(unsigned m=0;m<n_master;m++)
238  {
239  //Get the number of the unknown
240  //If the pressure dof is hanging
241  if(pressure_dof_is_hanging[l])
242  {
243  //Get the local equation from the master node
244  local_eqn =
245  bulk_el_pt->local_hang_eqn(hang_info_pt->master_node_pt(m),
246  p_index);
247  //Get the weight for the node
248  hang_weight = hang_info_pt->master_weight(m);
249  }
250  else
251  {
252  local_eqn = bulk_el_pt->p_local_eqn(l);
253  hang_weight = 1.0;
254  }
255 
256  //If the equation is not pinned
257  if(local_eqn >= 0)
258  {
259  residuals[local_eqn] -=
260  re*flux*interpolated_press*testp[l]*W*hang_weight;
261 
262  // Jacobian too?
263  if(flag)
264  {
265  //Number of master nodes and weights
266  unsigned n_master2=1; double hang_weight2=1.0;
267 
268  //Loop over the pressure shape functions
269  for(unsigned l2=0;l2<n_pres;l2++)
270  {
271  //If the pressure dof is hanging
272  if(pressure_dof_is_hanging[l2])
273  {
274  hang_info2_pt =
275  bulk_el_pt->pressure_node_pt(l2)->hanging_pt(p_index);
276  // Pressure dof is hanging so it must be nodal-based
277  //Get the number of master nodes from the pressure node
278  n_master2 = hang_info2_pt->nmaster();
279  }
280  //Otherwise the node is its own master
281  else
282  {
283  n_master2 = 1;
284  }
285 
286  //Loop over the master nodes
287  for(unsigned m2=0;m2<n_master2;m2++)
288  {
289  //Get the number of the unknown
290  //If the pressure dof is hanging
291  if(pressure_dof_is_hanging[l2])
292  {
293  //Get the unknown from the master node
294  local_unknown =
295  bulk_el_pt->local_hang_eqn(
296  hang_info2_pt->master_node_pt(m2),
297  p_index);
298  //Get the weight from the hanging object
299  hang_weight2 = hang_info2_pt->master_weight(m2);
300  }
301  else
302  {
303  local_unknown = bulk_el_pt->p_local_eqn(l2);
304  hang_weight2 = 1.0;
305  }
306 
307  //If the unknown is not pinned
308  if(local_unknown >= 0)
309  {
310  jacobian(local_eqn,local_unknown)-=
311  re*flux*psip[l2]*testp[l]*W*hang_weight*hang_weight2;
312  }
313  }
314  }
315  }
316  }
317  } /*End of Jacobian calculation*/
318  } //End of if not boundary condition
319  }//End of loop over l
320 }
321 
322 
323 
324 
325 
326 
327 
328 ///////////////////////////////////////////////////////////////////////
329 ///////////////////////////////////////////////////////////////////////
330 ///////////////////////////////////////////////////////////////////////
331 
332 
333 
334 //======================================================================
335 /// Refineable version of the Navier--Stokes equations
336 ///
337 ///
338 //======================================================================
339 template<unsigned DIM>
341 public virtual NavierStokesEquations<DIM>,
342 public virtual RefineableElement,
343 public virtual ElementWithZ2ErrorEstimator
344 {
345  protected:
346 
347  /// \short Unpin all pressure dofs in the element
348  virtual void unpin_elemental_pressure_dofs()=0;
349 
350  /// \short Pin unused nodal pressure dofs (empty by default, because
351  /// by default pressure dofs are not associated with nodes)
353 
354  public:
355 
356  /// \short Constructor
358  NavierStokesEquations<DIM>(),
361 
362 
363  /// \short Loop over all elements in Vector (which typically contains
364  /// all the elements in a fluid mesh) and pin the nodal pressure degrees
365  /// of freedom that are not being used. Function uses
366  /// the member function
367  /// - \c RefineableNavierStokesEquations::
368  /// pin_elemental_redundant_nodal_pressure_dofs()
369  /// .
370  /// which is empty by default and should be implemented for
371  /// elements with nodal pressure degrees of freedom
372  /// (e.g. for refineable Taylor-Hood.)
374  element_pt)
375  {
376  // Loop over all elements and call the function that pins their
377  // unused nodal pressure data
378  unsigned n_element = element_pt.size();
379  for(unsigned e=0;e<n_element;e++)
380  {
381  dynamic_cast<RefineableNavierStokesEquations<DIM>*>(element_pt[e])->
383  }
384  }
385 
386  /// \short Unpin all pressure dofs in elements listed in vector.
388  element_pt)
389  {
390  // Loop over all elements
391  unsigned n_element = element_pt.size();
392  for(unsigned e=0;e<n_element;e++)
393  {
394  dynamic_cast<RefineableNavierStokesEquations<DIM>*>(element_pt[e])->
396  }
397  }
398 
399  /// \short Pointer to n_p-th pressure node (Default: NULL,
400  /// indicating that pressure is not based on nodal interpolation).
401  virtual Node* pressure_node_pt(const unsigned& n_p) {return NULL;}
402 
403  /// \short Compute the diagonal of the velocity/pressure mass matrices.
404  /// If which one=0, both are computed, otherwise only the pressure
405  /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
406  /// LSC version of the preconditioner only needs that one)
408  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
409  const unsigned& which_one=0);
410 
411 
412  /// Number of 'flux' terms for Z2 error estimation
413  unsigned num_Z2_flux_terms()
414  {
415  // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
416  return DIM + (DIM*(DIM-1))/2;
417  }
418 
419  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
420  /// in strain rate tensor.
422  {
423 #ifdef PARANOID
424  unsigned num_entries=DIM+(DIM*(DIM-1))/2;
425  if (flux.size() < num_entries)
426  {
427  std::ostringstream error_message;
428  error_message << "The flux vector has the wrong number of entries, "
429  << flux.size() << ", whereas it should be at least "
430  << num_entries << std::endl;
431  throw OomphLibError(error_message.str(),
432  OOMPH_CURRENT_FUNCTION,
433  OOMPH_EXCEPTION_LOCATION);
434  }
435 #endif
436 
437  // Get strain rate matrix
438  DenseMatrix<double> strainrate(DIM);
439  this->strain_rate(s,strainrate);
440 
441  // Pack into flux Vector
442  unsigned icount=0;
443 
444  // Start with diagonal terms
445  for(unsigned i=0;i<DIM;i++)
446  {
447  flux[icount]=strainrate(i,i);
448  icount++;
449  }
450 
451  //Off diagonals row by row
452  for(unsigned i=0;i<DIM;i++)
453  {
454  for(unsigned j=i+1;j<DIM;j++)
455  {
456  flux[icount]=strainrate(i,j);
457  icount++;
458  }
459  }
460  }
461 
462  /// Further build, pass the pointers down to the sons
464  {
465  //Find the father element
466  RefineableNavierStokesEquations<DIM>* cast_father_element_pt
467  = dynamic_cast<RefineableNavierStokesEquations<DIM>*>
468  (this->father_element_pt());
469 
470  //Set the viscosity ratio pointer
471  this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
472  //Set the density ratio pointer
473  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
474  //Set pointer to global Reynolds number
475  this->Re_pt = cast_father_element_pt->re_pt();
476  //Set pointer to global Reynolds number x Strouhal number (=Womersley)
477  this->ReSt_pt = cast_father_element_pt->re_st_pt();
478  //Set pointer to global Reynolds number x inverse Froude number
479  this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
480  //Set pointer to global gravity Vector
481  this->G_pt = cast_father_element_pt->g_pt();
482 
483  //Set pointer to body force function
484  this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
485 
486  //Set pointer to volumetric source function
487  this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
488 
489  //Set the ALE flag
490  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
491  }
492 
493 
494  /// \short Compute the derivatives of the i-th component of
495  /// velocity at point s with respect
496  /// to all data that can affect its value. In addition, return the global
497  /// equation numbers corresponding to the data.
498  /// Overload the non-refineable version to take account of hanging node
499  /// information
501  const unsigned &i,
502  Vector<double> &du_ddata,
503  Vector<unsigned> &global_eqn_number)
504  {
505  //Find number of nodes
506  unsigned n_node = this->nnode();
507  //Local shape function
508  Shape psi(n_node);
509  //Find values of shape function at the given local coordinate
510  this->shape(s,psi);
511 
512  //Find the index at which the velocity component is stored
513  const unsigned u_nodal_index = this->u_index_nst(i);
514 
515  //Storage for hang info pointer
516  HangInfo* hang_info_pt=0;
517  //Storage for global equation
518  int global_eqn = 0;
519 
520  //Find the number of dofs associated with interpolated u
521  unsigned n_u_dof=0;
522  for(unsigned l=0;l<n_node;l++)
523  {
524  unsigned n_master = 1;
525 
526  //Local bool (is the node hanging)
527  bool is_node_hanging = this->node_pt(l)->is_hanging();
528 
529  //If the node is hanging, get the number of master nodes
530  if(is_node_hanging)
531  {
532  hang_info_pt = this->node_pt(l)->hanging_pt();
533  n_master = hang_info_pt->nmaster();
534  }
535  //Otherwise there is just one master node, the node itself
536  else
537  {
538  n_master = 1;
539  }
540 
541  //Loop over the master nodes
542  for(unsigned m=0;m<n_master;m++)
543  {
544  //Get the equation number
545  if(is_node_hanging)
546  {
547  //Get the equation number from the master node
548  global_eqn = hang_info_pt->master_node_pt(m)->
549  eqn_number(u_nodal_index);
550  }
551  else
552  {
553  // Global equation number
554  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
555  }
556 
557  //If it's positive add to the count
558  if(global_eqn >= 0) {++n_u_dof;}
559  }
560  }
561 
562  //Now resize the storage schemes
563  du_ddata.resize(n_u_dof,0.0);
564  global_eqn_number.resize(n_u_dof,0);
565 
566  //Loop over th nodes again and set the derivatives
567  unsigned count=0;
568  //Loop over the local nodes and sum
569  for(unsigned l=0;l<n_node;l++)
570  {
571  unsigned n_master = 1;
572  double hang_weight = 1.0;
573 
574  //Local bool (is the node hanging)
575  bool is_node_hanging = this->node_pt(l)->is_hanging();
576 
577  //If the node is hanging, get the number of master nodes
578  if(is_node_hanging)
579  {
580  hang_info_pt = this->node_pt(l)->hanging_pt();
581  n_master = hang_info_pt->nmaster();
582  }
583  //Otherwise there is just one master node, the node itself
584  else
585  {
586  n_master = 1;
587  }
588 
589  //Loop over the master nodes
590  for(unsigned m=0;m<n_master;m++)
591  {
592  //If the node is hanging get weight from master node
593  if(is_node_hanging)
594  {
595  //Get the hang weight from the master node
596  hang_weight = hang_info_pt->master_weight(m);
597  }
598  else
599  {
600  // Node contributes with full weight
601  hang_weight = 1.0;
602  }
603 
604  //Get the equation number
605  if(is_node_hanging)
606  {
607  //Get the equation number from the master node
608  global_eqn = hang_info_pt->master_node_pt(m)->
609  eqn_number(u_nodal_index);
610  }
611  else
612  {
613  // Local equation number
614  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
615  }
616 
617  if (global_eqn >= 0)
618  {
619  //Set the global equation number
620  global_eqn_number[count] = global_eqn;
621  //Set the derivative with respect to the unknown
622  du_ddata[count] = psi[l]*hang_weight;
623  //Increase the counter
624  ++count;
625  }
626  }
627  }
628  }
629 
630 
631 
632  protected:
633 
634 
635 /// \short Add element's contribution to elemental residual vector and/or
636 /// Jacobian matrix
637 /// flag=1: compute both
638 /// flag=0: compute only residual vector
640  Vector<double> &residuals,
641  DenseMatrix<double> &jacobian,
642  DenseMatrix<double> &mass_matrix,
643  unsigned flag);
644 
645  /// \short Compute the residuals for the associated pressure advection
646  /// diffusion problem. Used by the Fp preconditioner.
647  /// flag=1(or 0): do (or don't) compute the Jacobian as well.
649  Vector<double> &residuals, DenseMatrix<double> &jacobian, unsigned flag);
650 
651 
652  /// \short Compute derivatives of elemental residual vector with respect
653  /// to nodal coordinates. Overwrites default implementation in
654  /// FiniteElement base class.
655  /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
657  dresidual_dnodal_coordinates);
658 
659 };
660 
661 
662 //======================================================================
663 /// Refineable version of Taylor Hood elements. These classes
664 /// can be written in total generality.
665 //======================================================================
666 template<unsigned DIM>
668 public QTaylorHoodElement<DIM>,
669 public virtual RefineableNavierStokesEquations<DIM>,
670 public virtual RefineableQElement<DIM>
671 {
672  private:
673 
674  /// Unpin all pressure dofs
676  {
677  //find the index at which the pressure is stored
678  int p_index = this->p_nodal_index_nst();
679  unsigned n_node = this->nnode();
680  // loop over nodes
681  for(unsigned n=0;n<n_node;n++)
682  {this->node_pt(n)->unpin(p_index);}
683  }
684 
685  /// Pin all nodal pressure dofs that are not required
687  {
688  //Find the pressure index
689  int p_index = this->p_nodal_index_nst();
690  //Loop over all nodes
691  unsigned n_node = this->nnode();
692  // loop over all nodes and pin all the nodal pressures
693  for(unsigned n=0;n<n_node;n++) {this->node_pt(n)->pin(p_index);}
694 
695  // Loop over all actual pressure nodes and unpin if they're not hanging
696  unsigned n_pres = this->npres_nst();
697  for(unsigned l=0;l<n_pres;l++)
698  {
699  Node* nod_pt = this->pressure_node_pt(l);
700  if (!nod_pt->is_hanging(p_index)) {nod_pt->unpin(p_index);}
701  }
702  }
703 
704  public:
705 
706  /// \short Constructor
710  RefineableQElement<DIM>(),
711  QTaylorHoodElement<DIM>() {}
712 
713  /// \short Number of values required at local node n. In order to simplify
714  /// matters, we allocate storage for pressure variables at all the nodes
715  /// and then pin those that are not used.
716  unsigned required_nvalue(const unsigned &n) const {return DIM+1;}
717 
718  /// Number of continuously interpolated values: (DIM velocities + 1 pressure)
719  unsigned ncont_interpolated_values() const {return DIM+1;}
720 
721  /// Rebuild from sons: empty
722  void rebuild_from_sons(Mesh* &mesh_pt) {}
723 
724  /// \short Order of recovery shape functions for Z2 error estimation:
725  /// Same order as shape functions.
726  unsigned nrecovery_order() {return 2;}
727 
728  /// \short Number of vertex nodes in the element
729  unsigned nvertex_node() const
731 
732  /// \short Pointer to the j-th vertex node in the element
733  Node* vertex_node_pt(const unsigned& j) const
735 
736 /// \short Get the function value u in Vector.
737 /// Note: Given the generality of the interface (this function
738 /// is usually called from black-box documentation or interpolation routines),
739 /// the values Vector sets its own size in here.
741  {
742  // Set size of Vector: u,v,p and initialise to zero
743  values.resize(DIM+1,0.0);
744 
745  //Calculate velocities: values[0],...
746  for(unsigned i=0;i<DIM;i++) {values[i] = this->interpolated_u_nst(s,i);}
747 
748  //Calculate pressure: values[DIM]
749  values[DIM] = this->interpolated_p_nst(s);
750  }
751 
752 /// \short Get the function value u in Vector.
753 /// Note: Given the generality of the interface (this function
754 /// is usually called from black-box documentation or interpolation routines),
755 /// the values Vector sets its own size in here.
756  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
757  Vector<double>& values)
758  {
759  // Set size of Vector: u,v,p
760  values.resize(DIM+1);
761 
762  // Initialise
763  for(unsigned i=0;i<DIM+1;i++) {values[i]=0.0;}
764 
765  //Find out how many nodes there are
766  unsigned n_node = this->nnode();
767 
768  // Shape functions
769  Shape psif(n_node);
770  this->shape(s,psif);
771 
772  //Calculate velocities: values[0],...
773  for(unsigned i=0;i<DIM;i++)
774  {
775  //Get the index at which the i-th velocity is stored
776  unsigned u_nodal_index = this->u_index_nst(i);
777  for(unsigned l=0;l<n_node;l++)
778  {
779  values[i] += this->nodal_value(t,l,u_nodal_index)*psif[l];
780  }
781  }
782 
783  //Calculate pressure: values[DIM]
784  //(no history is carried in the pressure)
785  values[DIM] = this->interpolated_p_nst(s);
786  }
787 
788 
789  /// \short Perform additional hanging node procedures for variables
790  /// that are not interpolated by all nodes. The pressures are stored
791  /// at the p_nodal_index_nst-th location in each node
793  {
794  this->setup_hang_for_value(this->p_nodal_index_nst());
795  }
796 
797  /// \short Pointer to n_p-th pressure node
798  Node* pressure_node_pt(const unsigned &n_p)
799  {return this->node_pt(this->Pconv[n_p]);}
800 
801  /// \short The velocities are isoparametric and so the "nodes" interpolating
802  /// the velocities are the geometric nodes. The pressure "nodes" are a
803  /// subset of the nodes, so when value_id==DIM, the n-th pressure
804  /// node is returned.
805  Node* interpolating_node_pt(const unsigned &n,
806  const int &value_id)
807 
808  {
809  //The only different nodes are the pressure nodes
810  if(value_id==DIM) {return this->pressure_node_pt(n);}
811  //The other variables are interpolated via the usual nodes
812  else {return this->node_pt(n);}
813  }
814 
815  /// \short The pressure nodes are the corner nodes, so when n_value==DIM,
816  /// the fraction is the same as the 1d node number, 0 or 1.
817  double local_one_d_fraction_of_interpolating_node(const unsigned &n1d,
818  const unsigned &i,
819  const int &value_id)
820  {
821  if(value_id==DIM)
822  {
823  //The pressure nodes are just located on the boundaries at 0 or 1
824  return double(n1d);
825  }
826  //Otherwise the velocity nodes are the same as the geometric ones
827  else
828  {
829  return this->local_one_d_fraction_of_node(n1d,i);
830  }
831  }
832 
833  /// \short The velocity nodes are the same as the geometric nodes. The
834  /// pressure nodes must be calculated by using the same methods as
835  /// the geometric nodes, but by recalling that there are only two pressure
836  /// nodes per edge.
838  const int &value_id)
839  {
840  //If we are calculating pressure nodes
841  if(value_id==DIM)
842  {
843  //Storage for the index of the pressure node
844  unsigned total_index=0;
845  //The number of nodes along each 1d edge is 2.
846  unsigned NNODE_1D = 2;
847  //Storage for the index along each boundary
848  Vector<int> index(DIM);
849  //Loop over the coordinates
850  for(unsigned i=0;i<DIM;i++)
851  {
852  //If we are at the lower limit, the index is zero
853  if(s[i]==-1.0)
854  {
855  index[i] = 0;
856  }
857  //If we are at the upper limit, the index is the number of nodes minus 1
858  else if(s[i] == 1.0)
859  {
860  index[i] = NNODE_1D-1;
861  }
862  //Otherwise, we have to calculate the index in general
863  else
864  {
865  //For uniformly spaced nodes the 0th node number would be
866  double float_index = 0.5*(1.0 + s[i])*(NNODE_1D-1);
867  index[i] = int(float_index);
868  //What is the excess. This should be safe because the
869  //taking the integer part rounds down
870  double excess = float_index - index[i];
871  //If the excess is bigger than our tolerance there is no node,
872  //return null
873  if((excess > FiniteElement::Node_location_tolerance) && (
874  (1.0 - excess) > FiniteElement::Node_location_tolerance))
875  {
876  return 0;
877  }
878  }
879  ///Construct the general pressure index from the components.
880  total_index += index[i]*
881  static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
882  static_cast<int>(i)));
883  }
884  //If we've got here we have a node, so let's return a pointer to it
885  return this->pressure_node_pt(total_index);
886  }
887  //Otherwise velocity nodes are the same as pressure nodes
888  else
889  {
890  return this->get_node_at_local_coordinate(s);
891  }
892  }
893 
894 
895  /// \short The number of 1d pressure nodes is 2, the number of 1d velocity
896  /// nodes is the same as the number of 1d geometric nodes.
897  unsigned ninterpolating_node_1d(const int &value_id)
898  {
899  if(value_id==DIM) {return 2;}
900  else {return this->nnode_1d();}
901  }
902 
903  /// \short The number of pressure nodes is 2^DIM. The number of
904  /// velocity nodes is the same as the number of geometric nodes.
905  unsigned ninterpolating_node(const int &value_id)
906  {
907  if(value_id==DIM)
908  {return static_cast<unsigned>(pow(2.0,static_cast<int>(DIM)));}
909  else {return this->nnode();}
910  }
911 
912  /// \short The basis interpolating the pressure is given by pshape().
913  //// The basis interpolating the velocity is shape().
915  Shape &psi,
916  const int &value_id) const
917  {
918  if(value_id==DIM) {return this->pshape_nst(s,psi);}
919  else {return this->shape(s,psi);}
920  }
921 
922 
923 
924  /// \short Build FaceElements that apply the Robin boundary condition
925  /// to the pressure advection diffusion problem required by
926  /// Fp preconditioner
928  face_index)
929  {
932  this, face_index));
933  }
934 
935 
936  /// \short Add to the set \c paired_load_data pairs containing
937  /// - the pointer to a Data object
938  /// and
939  /// - the index of the value in that Data object
940  /// .
941  /// for all values (pressures, velocities) that affect the
942  /// load computed in the \c get_load(...) function.
943  /// (Overloads non-refineable version and takes hanging nodes
944  /// into account)
946  std::set<std::pair<Data*,unsigned> > &paired_load_data)
947  {
948  //Get the nodal indices at which the velocities are stored
949  unsigned u_index[DIM];
950  for(unsigned i=0;i<DIM;i++) {u_index[i] = this->u_index_nst(i);}
951 
952  //Loop over the nodes
953  unsigned n_node = this->nnode();
954  for(unsigned n=0;n<n_node;n++)
955  {
956  // Pointer to current node
957  Node* nod_pt=this->node_pt(n);
958 
959  // Check if it's hanging:
960  if (nod_pt->is_hanging())
961  {
962  // It's hanging -- get number of master nodes
963  unsigned nmaster=nod_pt->hanging_pt()->nmaster();
964 
965  // Loop over masters
966  for (unsigned j=0;j<nmaster;j++)
967  {
968  Node* master_nod_pt=nod_pt->hanging_pt()->master_node_pt(j);
969 
970  //Loop over the velocity components and add pointer to their data
971  //and indices to the vectors
972  for(unsigned i=0;i<DIM;i++)
973  {
974  paired_load_data.insert(std::make_pair(master_nod_pt,u_index[i]));
975  }
976  }
977  }
978  //Not hanging
979  else
980  {
981  //Loop over the velocity components and add pointer to their data
982  //and indices to the vectors
983  for(unsigned i=0;i<DIM;i++)
984  {
985  paired_load_data.insert(std::make_pair(this->node_pt(n),u_index[i]));
986  }
987  }
988  }
989 
990  //Get the nodal index at which the pressure is stored
991  int p_index = this->p_nodal_index_nst();
992 
993  //Loop over the pressure data
994  unsigned n_pres= this->npres_nst();
995  for(unsigned l=0;l<n_pres;l++)
996  {
997  //Get the pointer to the nodal pressure
998  Node* pres_node_pt = this->pressure_node_pt(l);
999  // Check if the pressure dof is hanging
1000  if(pres_node_pt->is_hanging(p_index))
1001  {
1002  //Get the pointer to the hang info object
1003  // (pressure is stored as p_index--th nodal dof).
1004  HangInfo* hang_info_pt = pres_node_pt->hanging_pt(p_index);
1005 
1006  // Get number of pressure master nodes (pressure is stored
1007  unsigned nmaster = hang_info_pt->nmaster();
1008 
1009  // Loop over pressure master nodes
1010  for(unsigned m=0;m<nmaster;m++)
1011  {
1012  //The p_index-th entry in each nodal data is the pressure, which
1013  //affects the traction
1014  paired_load_data.insert(
1015  std::make_pair(hang_info_pt->master_node_pt(m),p_index));
1016  }
1017  }
1018  // It's not hanging
1019  else
1020  {
1021  //The p_index-th entry in each nodal data is the pressure, which
1022  //affects the traction
1023  paired_load_data.insert(std::make_pair(pres_node_pt,p_index));
1024  }
1025  }
1026 
1027  }
1028 
1029 
1030 };
1031 
1032 
1033 //=======================================================================
1034 /// \short Face geometry of the RefineableQTaylorHoodElements is the
1035 /// same as the Face geometry of the QTaylorHoodElements.
1036 //=======================================================================
1037 template<unsigned DIM>
1039 public virtual FaceGeometry<QTaylorHoodElement<DIM> >
1040 {
1041  public:
1043 };
1044 
1045 
1046 //=======================================================================
1047 /// \short Face geometry of the face geometry of
1048 /// the RefineableQTaylorHoodElements is the
1049 /// same as the Face geometry of the Face geometry of QTaylorHoodElements.
1050 //=======================================================================
1051 template<unsigned DIM>
1053 public virtual FaceGeometry<FaceGeometry<QTaylorHoodElement<DIM> > >
1054 {
1055  public:
1057  {}
1058 };
1059 
1060 
1061 ///////////////////////////////////////////////////////////////////////////
1062 ///////////////////////////////////////////////////////////////////////////
1063 ///////////////////////////////////////////////////////////////////////////
1064 
1065 
1066 //======================================================================
1067 /// Refineable version of Crouzeix Raviart elements. Generic class definitions
1068 //======================================================================
1069 template<unsigned DIM>
1071 public QCrouzeixRaviartElement<DIM>,
1072 public virtual RefineableNavierStokesEquations<DIM>,
1073 public virtual RefineableQElement<DIM>
1074 {
1075  private:
1076 
1077  /// Unpin all internal pressure dofs
1079  {
1080  unsigned n_pres = this->npres_nst();
1081  // loop over pressure dofs and unpin them
1082  for(unsigned l=0;l<n_pres;l++)
1083  {this->internal_data_pt(this->P_nst_internal_index)->unpin(l);}
1084  }
1085 
1086  public:
1087 
1088  /// \short Constructor
1092  RefineableQElement<DIM>(),
1093  QCrouzeixRaviartElement<DIM>() {}
1094 
1095 
1096  /// Broken copy constructor
1099  {
1100  BrokenCopy::broken_copy("RefineableQCrouzeixRaviartElement");
1101  }
1102 
1103  /// Broken assignment operator
1104 //Commented out broken assignment operator because this can lead to a conflict warning
1105 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
1106 //realise that two separate implementations of the broken function are the same and so,
1107 //quite rightly, it shouts.
1108  /*void operator=(const RefineableQCrouzeixRaviartElement<DIM>&)
1109  {
1110  BrokenCopy::broken_assign("RefineableQCrouzeixRaviartElement");
1111  }*/
1112 
1113  /// Number of continuously interpolated values: DIM (velocities)
1114  unsigned ncont_interpolated_values() const {return DIM;}
1115 
1116  /// \short Rebuild from sons: Reconstruct pressure from the (merged) sons
1117  /// This must be specialised for each dimension.
1118  inline void rebuild_from_sons(Mesh* &mesh_pt);
1119 
1120  /// \short Order of recovery shape functions for Z2 error estimation:
1121  /// Same order as shape functions.
1122  unsigned nrecovery_order() {return 2;}
1123 
1124  /// \short Number of vertex nodes in the element
1125  unsigned nvertex_node() const
1127 
1128  /// \short Pointer to the j-th vertex node in the element
1129  Node* vertex_node_pt(const unsigned& j) const
1130  {
1132  }
1133 
1134 /// \short Get the function value u in Vector.
1135 /// Note: Given the generality of the interface (this function
1136 /// is usually called from black-box documentation or interpolation routines),
1137 /// the values Vector sets its own size in here.
1139  {
1140  // Set size of Vector: u,v,p and initialise to zero
1141  values.resize(DIM,0.0);
1142 
1143  //Calculate velocities: values[0],...
1144  for(unsigned i=0;i<DIM;i++) {values[i] = this->interpolated_u_nst(s,i);}
1145  }
1146 
1147  /// \short Get all function values [u,v..,p] at previous timestep t
1148  /// (t=0: present; t>0: previous timestep).
1149  ///
1150  /// Note: Given the generality of the interface (this function
1151  /// is usually called from black-box documentation or interpolation routines),
1152  /// the values Vector sets its own size in here.
1153  ///
1154  /// Note: No pressure history is kept, so pressure is always
1155  /// the current value.
1156  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
1157  Vector<double>& values)
1158  {
1159  // Set size of Vector: u,v,p
1160  values.resize(DIM);
1161 
1162  // Initialise
1163  for(unsigned i=0;i<DIM;i++) {values[i]=0.0;}
1164 
1165  //Find out how many nodes there are
1166  unsigned n_node = this->nnode();
1167 
1168  // Shape functions
1169  Shape psif(n_node);
1170  this->shape(s,psif);
1171 
1172  //Calculate velocities: values[0],...
1173  for(unsigned i=0;i<DIM;i++)
1174  {
1175  //Get the nodal index at which the i-th velocity component is stored
1176  unsigned u_nodal_index = this->u_index_nst(i);
1177  for(unsigned l=0;l<n_node;l++)
1178  {
1179  values[i] += this->nodal_value(t,l,u_nodal_index)*psif[l];
1180  }
1181  }
1182  }
1183 
1184  /// \short Perform additional hanging node procedures for variables
1185  /// that are not interpolated by all nodes. Empty
1187 
1188  /// Further build for Crouzeix_Raviart interpolates the internal
1189  /// pressure dofs from father element: Make sure pressure values and
1190  /// dp/ds agree between fathers and sons at the midpoints of the son
1191  /// elements. This must be specialised for each dimension.
1192  inline void further_build();
1193 
1194 
1195 
1196  /// \short Build FaceElements that apply the Robin boundary condition
1197  /// to the pressure advection diffusion problem required by
1198  /// Fp preconditioner
1200  face_index)
1201  {
1204  this, face_index));
1205  }
1206 
1207 
1208  /// \short Add to the set \c paired_load_data pairs containing
1209  /// - the pointer to a Data object
1210  /// and
1211  /// - the index of the value in that Data object
1212  /// .
1213  /// for all values (pressures, velocities) that affect the
1214  /// load computed in the \c get_load(...) function.
1215  /// (Overloads non-refineable version and takes hanging nodes
1216  /// into account)
1218  std::set<std::pair<Data*,unsigned> > &paired_load_data)
1219  {
1220  //Get the nodal indices at which the velocities are stored
1221  unsigned u_index[DIM];
1222  for(unsigned i=0;i<DIM;i++) {u_index[i] = this->u_index_nst(i);}
1223 
1224  //Loop over the nodes
1225  unsigned n_node = this->nnode();
1226  for(unsigned n=0;n<n_node;n++)
1227  {
1228  // Pointer to current node
1229  Node* nod_pt=this->node_pt(n);
1230 
1231  // Check if it's hanging:
1232  if (nod_pt->is_hanging())
1233  {
1234  // It's hanging -- get number of master nodes
1235  unsigned nmaster=nod_pt->hanging_pt()->nmaster();
1236 
1237  // Loop over masters
1238  for (unsigned j=0;j<nmaster;j++)
1239  {
1240  Node* master_nod_pt=nod_pt->hanging_pt()->master_node_pt(j);
1241 
1242  //Loop over the velocity components and add pointer to their data
1243  //and indices to the vectors
1244  for(unsigned i=0;i<DIM;i++)
1245  {
1246  paired_load_data.insert(std::make_pair(master_nod_pt,u_index[i]));
1247  }
1248  }
1249  }
1250  //Not hanging
1251  else
1252  {
1253  //Loop over the velocity components and add pointer to their data
1254  //and indices to the vectors
1255  for(unsigned i=0;i<DIM;i++)
1256  {
1257  paired_load_data.insert(std::make_pair(this->node_pt(n),u_index[i]));
1258  }
1259  }
1260  }
1261 
1262 
1263  //Loop over the pressure data (can't be hanging!)
1264  unsigned n_pres = this->npres_nst();
1265  for(unsigned l=0;l<n_pres;l++)
1266  {
1267  //The entries in the internal data at P_nst_internal_index
1268  //are the pressures, which affect the traction
1269  paired_load_data.insert(
1270  std::make_pair(this->internal_data_pt(this->P_nst_internal_index),l));
1271  }
1272  }
1273 
1274 
1275 };
1276 
1277 
1278 //======================================================================
1279 /// p-refineable version of Crouzeix Raviart elements. Generic class
1280 /// definitions
1281 //======================================================================
1282 template<unsigned DIM>
1284 public QCrouzeixRaviartElement<DIM>,
1285 public virtual RefineableNavierStokesEquations<DIM>,
1286 public virtual PRefineableQElement<DIM,3>
1287 {
1288  private:
1289 
1290  /// Unpin all internal pressure dofs
1292  {
1293  unsigned n_pres = this->npres_nst();
1294  n_pres = this->internal_data_pt(this->P_nst_internal_index)->nvalue();
1295  // loop over pressure dofs and unpin them
1296  for(unsigned l=0;l<n_pres;l++)
1297  {this->internal_data_pt(this->P_nst_internal_index)->unpin(l);}
1298  }
1299 
1300  public:
1301 
1302  /// \short Constructor
1306  PRefineableQElement<DIM,3>(),
1308  {
1309  // Set the p-order
1310  this->p_order()=3;
1311 
1312  // Set integration scheme
1313  // (To avoid memory leaks in pre-build and p-refine where new
1314  // integration schemes are created)
1316 
1317  // Resize pressure storage
1318  // (Constructor for QCrouzeixRaviartElement sets up DIM+1 pressure values)
1319  if (this->internal_data_pt(this->P_nst_internal_index)->nvalue()
1320  <= this->npres_nst())
1321  {
1323  ->resize(this->npres_nst());
1324  }
1325  else
1326  {
1327  Data* new_data_pt = new Data(this->npres_nst());
1328  delete this->internal_data_pt(this->P_nst_internal_index);
1329  this->internal_data_pt(this->P_nst_internal_index) = new_data_pt;
1330  }
1331  }
1332 
1333  /// \short Destructor
1335  {
1336  delete this->integral_pt();
1337  }
1338 
1339 
1340  /// Broken copy constructor
1343  {
1344  BrokenCopy::broken_copy("PRefineableQCrouzeixRaviartElement");
1345  }
1346 
1347  /// Broken assignment operator
1348  /*void operator=(const PRefineableQCrouzeixRaviartElement<DIM>&)
1349  {
1350  BrokenCopy::broken_assign("PRefineableQCrouzeixRaviartElement");
1351  }*/
1352 
1353  /// \short Return the i-th pressure value
1354  /// (Discontinous pressure interpolation -- no need to cater for hanging
1355  /// nodes).
1356  double p_nst(const unsigned &i) const
1357  {return this->internal_data_pt(this->P_nst_internal_index)->value(i);}
1358 
1359  double p_nst(const unsigned &t, const unsigned &i) const
1360  {return this->internal_data_pt(this->P_nst_internal_index)->value(t,i);}
1361 
1362  ///// Return number of pressure values
1363  unsigned npres_nst() const {return (this->p_order()-2)*(this->p_order()-2);}
1364 
1365  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
1366  void fix_pressure(const unsigned &p_dof, const double &p_value)
1367  {
1368  this->internal_data_pt(this->P_nst_internal_index)->pin(p_dof);
1369  this->internal_data_pt(this->P_nst_internal_index)->set_value(p_dof,p_value);
1370  }
1371 
1372  unsigned required_nvalue(const unsigned& n) const
1373  {return DIM;}
1374 
1375  /// Number of continuously interpolated values: DIM (velocities)
1376  unsigned ncont_interpolated_values() const {return DIM;}
1377 
1378  /// \short Rebuild from sons: Reconstruct pressure from the (merged) sons
1379  /// This must be specialised for each dimension.
1380  void rebuild_from_sons(Mesh* &mesh_pt)
1381  {
1382  //Do p-refineable version
1384  //Do Crouzeix-Raviart version
1385  //Need to reconstruct pressure manually!
1386  for(unsigned p=0; p<npres_nst(); p++)
1387  {
1388  //BENFLAG: Set to zero for now -- don't do projection problem yet
1389  this->internal_data_pt(this->P_nst_internal_index)->set_value(p,0.0);
1390  }
1391 
1392  }
1393 
1394  /// \short Order of recovery shape functions for Z2 error estimation:
1395  /// - Same order as shape functions.
1396  //unsigned nrecovery_order()
1397  // {
1398  // if(this->nnode_1d() < 4) {return (this->nnode_1d()-1);}
1399  // else {return 3;}
1400  // }
1401  /// - Constant recovery order, since recovery order of the first element
1402  /// is used for the whole mesh.
1403  unsigned nrecovery_order() {return 3;}
1404 
1405  /// \short Number of vertex nodes in the element
1406  unsigned nvertex_node() const
1408 
1409  /// \short Pointer to the j-th vertex node in the element
1410  Node* vertex_node_pt(const unsigned& j) const
1411  {
1413  }
1414 
1415  /// \short Velocity shape and test functions and their derivs
1416  /// w.r.t. to global coords at local coordinate s (taken from geometry)
1417  ///Return Jacobian of mapping between local and global coordinates.
1418  inline double dshape_and_dtest_eulerian_nst(const Vector<double> &s,
1419  Shape &psi,
1420  DShape &dpsidx, Shape &test,
1421  DShape &dtestdx) const;
1422 
1423  /// \short Velocity shape and test functions and their derivs
1424  /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
1425  ///Return Jacobian of mapping between local and global coordinates.
1426  inline double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt,
1427  Shape &psi,
1428  DShape &dpsidx,
1429  Shape &test,
1430  DShape &dtestdx) const;
1431 
1432  /// Pressure shape functions at local coordinate s
1433  inline void pshape_nst(const Vector<double> &s, Shape &psi) const;
1434 
1435  /// Pressure shape and test functions at local coordinte s
1436  inline void pshape_nst(const Vector<double> &s, Shape &psi, Shape &test) const;
1437 
1438  /// \short Get the function value u in Vector.
1439  /// Note: Given the generality of the interface (this function
1440  /// is usually called from black-box documentation or interpolation routines),
1441  /// the values Vector sets its own size in here.
1443  {
1444  // Set size of Vector: u,v,p and initialise to zero
1445  values.resize(DIM,0.0);
1446 
1447  //Calculate velocities: values[0],...
1448  for(unsigned i=0;i<DIM;i++) {values[i] = this->interpolated_u_nst(s,i);}
1449  }
1450 
1451  /// \short Get all function values [u,v..,p] at previous timestep t
1452  /// (t=0: present; t>0: previous timestep).
1453  ///
1454  /// Note: Given the generality of the interface (this function
1455  /// is usually called from black-box documentation or interpolation routines),
1456  /// the values Vector sets its own size in here.
1457  ///
1458  /// Note: No pressure history is kept, so pressure is always
1459  /// the current value.
1460  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
1461  Vector<double>& values)
1462  {
1463  // Set size of Vector: u,v,p
1464  values.resize(DIM);
1465 
1466  // Initialise
1467  for(unsigned i=0;i<DIM;i++) {values[i]=0.0;}
1468 
1469  //Find out how many nodes there are
1470  unsigned n_node = this->nnode();
1471 
1472  // Shape functions
1473  Shape psif(n_node);
1474  this->shape(s,psif);
1475 
1476  //Calculate velocities: values[0],...
1477  for(unsigned i=0;i<DIM;i++)
1478  {
1479  //Get the nodal index at which the i-th velocity component is stored
1480  unsigned u_nodal_index = this->u_index_nst(i);
1481  for(unsigned l=0;l<n_node;l++)
1482  {
1483  values[i] += this->nodal_value(t,l,u_nodal_index)*psif[l];
1484  }
1485  }
1486  }
1487 
1488  /// \short Perform additional hanging node procedures for variables
1489  /// that are not interpolated by all nodes. Empty
1491 
1492  /// Further build for Crouzeix_Raviart interpolates the internal
1493  /// pressure dofs from father element: Make sure pressure values and
1494  /// dp/ds agree between fathers and sons at the midpoints of the son
1495  /// elements. This must be specialised for each dimension.
1496  void further_build();
1497 
1498 
1499 };
1500 
1501 
1502 //=======================================================================
1503 /// Face geometry of the RefineableQuadQCrouzeixRaviartElements
1504 //=======================================================================
1505 template<unsigned DIM>
1507 public virtual FaceGeometry<QCrouzeixRaviartElement<DIM> >
1508 {
1509  public:
1511 };
1512 
1513 //======================================================================
1514 /// \short Face geometry of the face geometry of
1515 /// the RefineableQCrouzeixRaviartElements is the
1516 /// same as the Face geometry of the Face geometry of
1517 /// QCrouzeixRaviartElements.
1518 //=======================================================================
1519 template<unsigned DIM>
1521 public virtual FaceGeometry<FaceGeometry<QCrouzeixRaviartElement<DIM> > >
1522 {
1523  public:
1525  {}
1526 };
1527 
1528 
1529 //Inline functions
1530 
1531 //=====================================================================
1532 /// 2D Rebuild from sons: Reconstruct pressure from the (merged) sons
1533 //=====================================================================
1534 template<>
1537 {
1538  using namespace QuadTreeNames;
1539 
1540  //Central pressure value:
1541  //-----------------------
1542 
1543  // Use average of the sons central pressure values
1544  // Other options: Take average of the four (discontinuous)
1545  // pressure values at the father's midpoint]
1546 
1547  double av_press=0.0;
1548 
1549  // Loop over the sons
1550  for (unsigned ison=0;ison<4;ison++)
1551  {
1552  // Add the sons midnode pressure
1553  // Note that we can assume that the pressure is stored at the same
1554  // location because these are EXACTLY the same type of elements
1555  av_press += quadtree_pt()->son_pt(ison)->object_pt()->
1556  internal_data_pt(this->P_nst_internal_index)->value(0);
1557  }
1558 
1559  // Use the average
1560  internal_data_pt(this->P_nst_internal_index)->set_value(0,0.25*av_press);
1561 
1562 
1563  //Slope in s_0 direction
1564  //----------------------
1565 
1566  // Use average of the 2 FD approximations based on the
1567  // elements central pressure values
1568  // [Other options: Take average of the four
1569  // pressure derivatives]
1570 
1571  double slope1=
1572  quadtree_pt()->son_pt(SE)->object_pt()->
1573  internal_data_pt(this->P_nst_internal_index)->value(0) -
1574  quadtree_pt()->son_pt(SW)->object_pt()->
1575  internal_data_pt(this->P_nst_internal_index)->value(0);
1576 
1577  double slope2 =
1578  quadtree_pt()->son_pt(NE)->object_pt()->
1579  internal_data_pt(this->P_nst_internal_index)->value(0) -
1580  quadtree_pt()->son_pt(NW)->object_pt()->
1581  internal_data_pt(this->P_nst_internal_index)->value(0);
1582 
1583 
1584  // Use the average
1585  internal_data_pt(this->P_nst_internal_index)->
1586  set_value(1,0.5*(slope1+slope2));
1587 
1588 
1589 
1590  //Slope in s_1 direction
1591  //----------------------
1592 
1593  // Use average of the 2 FD approximations based on the
1594  // elements central pressure values
1595  // [Other options: Take average of the four
1596  // pressure derivatives]
1597 
1598  slope1 =
1599  quadtree_pt()->son_pt(NE)->object_pt()
1600  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1601  quadtree_pt()->son_pt(SE)->object_pt()
1602  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1603 
1604  slope2=
1605  quadtree_pt()->son_pt(NW)->object_pt()
1606  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1607  quadtree_pt()->son_pt(SW)->object_pt()
1608  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1609 
1610 
1611  // Use the average
1612  internal_data_pt(this->P_nst_internal_index)->
1613  set_value(2,0.5*(slope1+slope2));
1614 }
1615 
1616 
1617 
1618 //=================================================================
1619 /// 3D Rebuild from sons: Reconstruct pressure from the (merged) sons
1620 //=================================================================
1621 template<>
1624 {
1625  using namespace OcTreeNames;
1626 
1627  //Central pressure value:
1628  //-----------------------
1629 
1630  // Use average of the sons central pressure values
1631  // Other options: Take average of the four (discontinuous)
1632  // pressure values at the father's midpoint]
1633 
1634  double av_press=0.0;
1635 
1636  // Loop over the sons
1637  for (unsigned ison=0;ison<8;ison++)
1638  {
1639  // Add the sons midnode pressure
1640  av_press += octree_pt()->son_pt(ison)->object_pt()->
1641  internal_data_pt(this->P_nst_internal_index)->value(0);
1642  }
1643 
1644  // Use the average
1645  internal_data_pt(this->P_nst_internal_index)->
1646  set_value(0,0.125*av_press);
1647 
1648 
1649  //Slope in s_0 direction
1650  //----------------------
1651 
1652  // Use average of the 4 FD approximations based on the
1653  // elements central pressure values
1654  // [Other options: Take average of the four
1655  // pressure derivatives]
1656 
1657  double slope1 =
1658  octree_pt()->son_pt(RDF)->object_pt()->
1659  internal_data_pt(this->P_nst_internal_index)->value(0) -
1660  octree_pt()->son_pt(LDF)->object_pt()->
1661  internal_data_pt(this->P_nst_internal_index)->value(0);
1662 
1663  double slope2 =
1664  octree_pt()->son_pt(RUF)->object_pt()->
1665  internal_data_pt(this->P_nst_internal_index)->value(0) -
1666  octree_pt()->son_pt(LUF)->object_pt()->
1667  internal_data_pt(this->P_nst_internal_index)->value(0);
1668 
1669  double slope3 =
1670  octree_pt()->son_pt(RDB)->object_pt()->
1671  internal_data_pt(this->P_nst_internal_index)->value(0) -
1672  octree_pt()->son_pt(LDB)->object_pt()->
1673  internal_data_pt(this->P_nst_internal_index)->value(0);
1674 
1675  double slope4 =
1676  octree_pt()->son_pt(RUB)->object_pt()->
1677  internal_data_pt(this->P_nst_internal_index)->value(0) -
1678  octree_pt()->son_pt(LUB)->object_pt()->
1679  internal_data_pt(this->P_nst_internal_index)->value(0);
1680 
1681 
1682  // Use the average
1683  internal_data_pt(this->P_nst_internal_index)->
1684  set_value(1,0.25*(slope1+slope2+slope3+slope4));
1685 
1686 
1687  //Slope in s_1 direction
1688  //----------------------
1689 
1690  // Use average of the 4 FD approximations based on the
1691  // elements central pressure values
1692  // [Other options: Take average of the four
1693  // pressure derivatives]
1694 
1695  slope1 =
1696  octree_pt()->son_pt(LUB)->object_pt()
1697  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1698  octree_pt()->son_pt(LDB)->object_pt()
1699  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1700 
1701  slope2 =
1702  octree_pt()->son_pt(RUB)->object_pt()
1703  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1704  octree_pt()->son_pt(RDB)->object_pt()
1705  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1706 
1707  slope3 =
1708  octree_pt()->son_pt(LUF)->object_pt()
1709  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1710  octree_pt()->son_pt(LDF)->object_pt()
1711  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1712 
1713  slope4 =
1714  octree_pt()->son_pt(RUF)->object_pt()
1715  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1716  octree_pt()->son_pt(RDF)->object_pt()
1717  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1718 
1719 
1720  // Use the average
1721  internal_data_pt(this->P_nst_internal_index)->
1722  set_value(2,0.25*(slope1+slope2+slope3+slope4));
1723 
1724 
1725  //Slope in s_2 direction
1726  //----------------------
1727 
1728  // Use average of the 4 FD approximations based on the
1729  // elements central pressure values
1730  // [Other options: Take average of the four
1731  // pressure derivatives]
1732 
1733  slope1 =
1734  octree_pt()->son_pt(LUF)->object_pt()
1735  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1736  octree_pt()->son_pt(LUB)->object_pt()
1737  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1738 
1739  slope2 =
1740  octree_pt()->son_pt(RUF)->object_pt()
1741  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1742  octree_pt()->son_pt(RUB)->object_pt()
1743  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1744 
1745  slope3 =
1746  octree_pt()->son_pt(LDF)->object_pt()
1747  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1748  octree_pt()->son_pt(LDB)->object_pt()
1749  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1750 
1751  slope4 =
1752  octree_pt()->son_pt(RDF)->object_pt()
1753  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1754  octree_pt()->son_pt(RDB)->object_pt()
1755  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1756 
1757  // Use the average
1758  internal_data_pt(this->P_nst_internal_index)->
1759  set_value(3,0.25*(slope1+slope2+slope3+slope4));
1760 
1761 }
1762 
1763 
1764 //======================================================================
1765 /// 2D Further build for Crouzeix_Raviart interpolates the internal
1766 /// pressure dofs from father element: Make sure pressure values and
1767 /// dp/ds agree between fathers and sons at the midpoints of the son
1768 /// elements.
1769 //======================================================================
1770 template<>
1772 {
1773  //Call the generic further build
1775 
1776  using namespace QuadTreeNames;
1777 
1778  // What type of son am I? Ask my quadtree representation...
1779  int son_type=quadtree_pt()->son_type();
1780 
1781  // Pointer to my father (in element impersonation)
1782  RefineableElement* father_el_pt= quadtree_pt()->father_pt()->object_pt();
1783 
1784  Vector<double> s_father(2);
1785 
1786  // Son midpoint is located at the following coordinates in father element:
1787 
1788  // South west son
1789  if (son_type==SW)
1790  {
1791  s_father[0]=-0.5;
1792  s_father[1]=-0.5;
1793  }
1794  // South east son
1795  else if (son_type==SE)
1796  {
1797  s_father[0]= 0.5;
1798  s_father[1]=-0.5;
1799  }
1800  // North east son
1801  else if (son_type==NE)
1802  {
1803  s_father[0]=0.5;
1804  s_father[1]=0.5;
1805  }
1806 
1807  // North west son
1808  else if (son_type==NW)
1809  {
1810  s_father[0]=-0.5;
1811  s_father[1]= 0.5;
1812  }
1813 
1814  // Pressure value in father element
1815  RefineableQCrouzeixRaviartElement<2>* cast_father_element_pt=
1816  dynamic_cast<RefineableQCrouzeixRaviartElement<2>*>(father_el_pt);
1817 
1818  double press=cast_father_element_pt->interpolated_p_nst(s_father);
1819 
1820  // Pressure value gets copied straight into internal dof:
1821  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1822 
1823  // The slopes get copied from father
1824  for(unsigned i=1;i<3;i++)
1825  {
1826  double half_father_slope = 0.5*cast_father_element_pt->
1827  internal_data_pt(this->P_nst_internal_index)->value(i);
1828  //Set the value in the son
1829  internal_data_pt(this->P_nst_internal_index)->
1830  set_value(i,half_father_slope);
1831  }
1832 }
1833 
1834 
1835 //=======================================================================
1836 /// 3D Further build for Crouzeix_Raviart interpolates the internal
1837 /// pressure dofs from father element: Make sure pressure values and
1838 /// dp/ds agree between fathers and sons at the midpoints of the son
1839 /// elements.
1840 //=======================================================================
1841 template<>
1843 {
1845 
1846  using namespace OcTreeNames;
1847 
1848  // What type of son am I? Ask my octree representation...
1849  int son_type=octree_pt()->son_type();
1850 
1851  // Pointer to my father (in element impersonation)
1852  RefineableQElement<3>* father_el_pt=
1853  dynamic_cast<RefineableQElement<3>*>(
1854  octree_pt()->father_pt()->object_pt());
1855 
1856  Vector<double> s_father(3);
1857 
1858  // Son midpoint is located at the following coordinates in father element:
1859  for(unsigned i=0;i<3;i++)
1860  {
1861  s_father[i]=0.5*OcTree::Direction_to_vector[son_type][i];
1862  }
1863 
1864  // Pressure value in father element
1865  RefineableQCrouzeixRaviartElement<3>* cast_father_element_pt=
1866  dynamic_cast<RefineableQCrouzeixRaviartElement<3>*>(father_el_pt);
1867 
1868  double press=cast_father_element_pt->interpolated_p_nst(s_father);
1869 
1870  // Pressure value gets copied straight into internal dof:
1871  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1872 
1873  // The slopes get copied from father
1874  for(unsigned i=1;i<4;i++)
1875  {
1876  double half_father_slope = 0.5*cast_father_element_pt
1877  ->internal_data_pt(this->P_nst_internal_index)->value(i);
1878  //Set the value
1879  internal_data_pt(this->P_nst_internal_index)
1880  ->set_value(i,half_father_slope);
1881  }
1882 }
1883 
1884 //=======================================================================
1885 /// 2D
1886 /// Derivatives of the shape functions and test functions w.r.t. to global
1887 /// (Eulerian) coordinates. Return Jacobian of mapping between
1888 /// local and global coordinates.
1889 //=======================================================================
1890 template<>
1893  DShape &dpsidx, Shape &test,
1894  DShape &dtestdx) const
1895 {
1896  //Call the geometrical shape functions and derivatives
1897  double J = this->dshape_eulerian(s,psi,dpsidx);
1898 
1899  //Loop over the test functions and derivatives and set them equal to the
1900  //shape functions
1901  for(unsigned i=0;i<nnode_1d()*nnode_1d();i++)
1902  {
1903  test[i] = psi[i];
1904  dtestdx(i,0) = dpsidx(i,0);
1905  dtestdx(i,1) = dpsidx(i,1);
1906  }
1907 
1908  //Return the jacobian
1909  return J;
1910 }
1911 
1912 //=======================================================================
1913 /// 2D
1914 /// Derivatives of the shape functions and test functions w.r.t. to global
1915 /// (Eulerian) coordinates. Return Jacobian of mapping between
1916 /// local and global coordinates.
1917 //=======================================================================
1918 template<>
1921  DShape &dpsidx, Shape &test,
1922  DShape &dtestdx) const
1923 {
1924  //Call the geometrical shape functions and derivatives
1925  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
1926 
1927  //Loop over the test functions and derivatives and set them equal to the
1928  //shape functions
1929  for(unsigned i=0;i<nnode_1d()*nnode_1d();i++)
1930  {
1931  test[i] = psi[i];
1932  dtestdx(i,0) = dpsidx(i,0);
1933  dtestdx(i,1) = dpsidx(i,1);
1934  }
1935 
1936  //Return the jacobian
1937  return J;
1938 }
1939 
1940 //=======================================================================
1941 /// 3D
1942 /// Derivatives of the shape functions and test functions w.r.t. to global
1943 /// (Eulerian) coordinates. Return Jacobian of mapping between
1944 /// local and global coordinates.
1945 //=======================================================================
1946 template<>
1949  DShape &dpsidx, Shape &test,
1950  DShape &dtestdx) const
1951 {
1952  //Call the geometrical shape functions and derivatives
1953  double J = this->dshape_eulerian(s,psi,dpsidx);
1954 
1955  //Loop over the test functions and derivatives and set them equal to the
1956  //shape functions
1957  for(unsigned i=0;i<nnode_1d()*nnode_1d()*nnode_1d();i++)
1958  {
1959  test[i] = psi[i];
1960  dtestdx(i,0) = dpsidx(i,0);
1961  dtestdx(i,1) = dpsidx(i,1);
1962  dtestdx(i,2) = dpsidx(i,2);
1963  }
1964 
1965  //Return the jacobian
1966  return J;
1967 }
1968 
1969 //=======================================================================
1970 /// 3D
1971 /// Derivatives of the shape functions and test functions w.r.t. to global
1972 /// (Eulerian) coordinates. Return Jacobian of mapping between
1973 /// local and global coordinates.
1974 //=======================================================================
1975 template<>
1978  DShape &dpsidx, Shape &test,
1979  DShape &dtestdx) const
1980 {
1981  //Call the geometrical shape functions and derivatives
1982  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
1983 
1984  //Loop over the test functions and derivatives and set them equal to the
1985  //shape functions
1986  for(unsigned i=0;i<nnode_1d()*nnode_1d()*nnode_1d();i++)
1987  {
1988  test[i] = psi[i];
1989  dtestdx(i,0) = dpsidx(i,0);
1990  dtestdx(i,1) = dpsidx(i,1);
1991  dtestdx(i,2) = dpsidx(i,2);
1992  }
1993 
1994  //Return the jacobian
1995  return J;
1996 }
1997 
1998 //=======================================================================
1999 /// 2D :
2000 /// Pressure shape functions
2001 //=======================================================================
2002 template<>
2004 pshape_nst(const Vector<double> &s, Shape &psi) const
2005 {
2006  unsigned npres = this->npres_nst();
2007  if (npres==1)
2008  {
2009  psi[0] = 1.0;
2010  }
2011  else
2012  {
2013  // Get number of pressure modes
2014  unsigned npres_1d = (int) std::sqrt((double)npres);
2015 
2016  //Local storage
2017  //Call the one-dimensional modal shape functions
2018  OneDimensionalModalShape psi1(npres_1d,s[0]);
2019  OneDimensionalModalShape psi2(npres_1d,s[1]);
2020 
2021  //Now let's loop over the nodal points in the element
2022  //s1 is the "x" coordinate, s2 the "y"
2023  for(unsigned i=0;i<npres_1d;i++)
2024  {
2025  for(unsigned j=0;j<npres_1d;j++)
2026  {
2027  //Multiply the two 1D functions together to get the 2D function
2028  psi[i*npres_1d + j] = psi2[i]*psi1[j];
2029  }
2030  }
2031  }
2032 }
2033 
2034 ///Define the pressure shape and test functions
2035 template<>
2037 pshape_nst(const Vector<double> &s, Shape &psi, Shape &test) const
2038 {
2039  //Call the pressure shape functions
2040  pshape_nst(s,psi);
2041 
2042  //Loop over the test functions and set them equal to the shape functions
2043  if (this->npres_nst()==1)
2044  {
2045  test[0] = psi[0];
2046  }
2047  else
2048  {
2049  for(unsigned i=0;i<this->npres_nst();i++) test[i] = psi[i];
2050  }
2051 }
2052 
2053 //=======================================================================
2054 /// 3D :
2055 /// Pressure shape functions
2056 //=======================================================================
2057 template<>
2059 pshape_nst(const Vector<double> &s, Shape &psi) const
2060 {
2061  unsigned npres = this->npres_nst();
2062  if (npres==1)
2063  {
2064  psi[0] = 1.0;
2065  }
2066  else
2067  {
2068  // Get number of pressure modes
2069  unsigned npres_1d = (int) std::sqrt((double)npres);
2070 
2071  //Local storage
2072  //Call the one-dimensional modal shape functions
2073  OneDimensionalModalShape psi1(npres_1d,s[0]);
2074  OneDimensionalModalShape psi2(npres_1d,s[1]);
2075  OneDimensionalModalShape psi3(npres_1d,s[2]);
2076 
2077  //Now let's loop over the nodal points in the element
2078  //s1 is the "x" coordinate, s2 the "y"
2079  for(unsigned i=0;i<npres_1d;i++)
2080  {
2081  for(unsigned j=0;j<npres_1d;j++)
2082  {
2083  for(unsigned k=0;k<npres_1d;k++)
2084  {
2085  //Multiply the two 1D functions together to get the 2D function
2086  psi[i*npres_1d*npres_1d + j*npres_1d + k] = psi3[i]*psi2[j]*psi1[k];
2087  }
2088  }
2089  }
2090  }
2091 }
2092 
2093 ///Define the pressure shape and test functions
2094 template<>
2096 pshape_nst(const Vector<double> &s, Shape &psi, Shape &test) const
2097 {
2098  //Call the pressure shape functions
2099  pshape_nst(s,psi);
2100 
2101  //Loop over the test functions and set them equal to the shape functions
2102  if (this->npres_nst()==1)
2103  {
2104  test[0] = psi[0];
2105  }
2106  else
2107  {
2108  for(unsigned i=0;i<this->npres_nst();i++) test[i] = psi[i];
2109  }
2110 }
2111 
2112 }
2113 
2114 #endif
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
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 unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:386
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor. ...
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. Overwrites default implementation in FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}.
Vector< FpPressureAdvDiffRobinBCElementBase * > Pressure_advection_diffusion_robin_element_pt
Storage for FaceElements that apply Robin BC for pressure adv diff equation used in Fp preconditioner...
void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
unsigned nvertex_node() const
Number of vertex nodes in the element.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
cstr elem_len * i
Definition: cfortran.h:607
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4367
void dinterpolated_u_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 rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons This must be specialised for each dime...
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3148
void unpin_elemental_pressure_dofs()
Unpin all internal pressure dofs.
double * Re_pt
Pointer to global Reynolds number.
NavierStokesSourceFctPt & source_fct_pt()
Access function for the source-function pointer.
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
unsigned nvertex_node() const
Number of vertex nodes in the element.
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
A general Finite Element class.
Definition: elements.h:1271
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...
double p_nst(const unsigned &i) const
Broken assignment operator.
char t
Definition: cfortran.h:572
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number) ...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
double *& re_pt()
Pointer to Reynolds number.
void unpin_elemental_pressure_dofs()
Unpin all pressure dofs.
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
virtual void fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
This function returns the residuals for the traction function. flag=1 (or 0): do (or don't) compute t...
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
e
Definition: cfortran.h:575
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
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 ncont_interpolated_values() const
Broken assignment operator.
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
virtual int p_nodal_index_nst() const
Set the value at which the pressure is stored in the nodes.
void unpin_elemental_pressure_dofs()
Unpin all internal pressure dofs.
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
Node * interpolating_node_pt(const unsigned &n, const int &value_id)
The velocities are isoparametric and so the "nodes" interpolating the velocities are the geometric no...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
unsigned required_nvalue(const unsigned &n) const
Number of values required at local node n. In order to simplify matters, we allocate storage for pres...
double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at ipt-th integation point...
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
void pin_elemental_redundant_nodal_pressure_dofs()
Pin all nodal pressure dofs that are not required.
A Rank 3 Tensor class.
Definition: matrices.h:1337
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 further_build()
Further build, pass the pointers down to the sons.
Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node.
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
double p_nst(const unsigned &t, const unsigned &i) const
Return the i-th pressure value (Discontinous pressure interpolation – no need to cater for hanging no...
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)...
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
static Vector< Vector< int > > Direction_to_vector
For each direction, i.e. a son_type (vertex), a face or an edge, this defines a vector that indicates...
Definition: octree.h:342
static char t char * s
Definition: cfortran.h:572
RefineableFpPressureAdvDiffRobinBCElement(FiniteElement *const &element_pt, const int &face_index)
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 *& re_invfr_pt()
Pointer to global inverse Froude number.
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
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)...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
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 further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
Refineable version of Crouzeix Raviart elements. Generic class definitions.
Vector< double > * G_pt
Pointer to global gravity Vector.
unsigned npres_nst() const
Return number of pressure values.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: DIM (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...
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 const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
double *& density_ratio_pt()
Pointer to Density ratio.
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
void pshape_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:684
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
unsigned ninterpolating_node(const int &value_id)
The number of pressure nodes is 2^DIM. The number of velocity nodes is the same as the number of geom...
unsigned nvertex_node() const
Number of vertex nodes in the element.
double dshape_and_dtest_eulerian_nst(const Vector< double > &s, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
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
void interpolating_basis(const Vector< double > &s, Shape &psi, const int &value_id) const
The basis interpolating the pressure is given by pshape(). / The basis interpolating the velocity is ...
double local_one_d_fraction_of_interpolating_node(const unsigned &n1d, const unsigned &i, const int &value_id)
The pressure nodes are the corner nodes, so when n_value==DIM, the fraction is the same as the 1d nod...
void fix_pressure(const unsigned &p_dof, const double &p_value)
Pin p_dof-th pressure dof and set it to value specified by p_value.
void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed...
NavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
Node * get_interpolating_node_at_local_coordinate(const Vector< double > &s, const int &value_id)
The velocity nodes are the same as the geometric nodes. The pressure nodes must be calculated by usin...
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
unsigned npres_nst() const
Return number of pressure values.
void fill_in_generic_residual_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
virtual void resize(const unsigned &n_value)
Change (increase) the number of values that may be stored.
Definition: nodes.cc:961
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
void fill_in_generic_pressure_advection_diffusion_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number) ...
void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
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
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons This must be specialised for each dime...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned npres_nst() const
Return number of pressure values.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation:
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
static void unpin_all_pressure_dofs(const Vector< GeneralisedElement * > &element_pt)
Unpin all pressure dofs in elements listed in vector.
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
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 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 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...
unsigned ninterpolating_node_1d(const int &value_id)
The number of 1d pressure nodes is 2, the number of 1d velocity nodes is the same as the number of 1d...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: (DIM velocities + 1 pressure)
void interpolated_u_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
virtual unsigned u_index_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.