generalised_newtonian_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: 1194 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-27 08:44:43 +0100 (Fri, 27 May 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header file for refineable 2D quad Navier Stokes elements
31 
32 #ifndef OOMPH_GENERALISED_NEWTONIAN_REFINEABLE_NAVIER_STOKES_ELEMENTS_HEADER
33 #define OOMPH_GENERALISED_NEWTONIAN_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"
46 
47 namespace oomph
48 {
49 
50 ///////////////////////////////////////////////////////////////////////
51 ///////////////////////////////////////////////////////////////////////
52 ///////////////////////////////////////////////////////////////////////
53 
54 
55 
56 
57 
58 //======================================================================
59 /// Refineable version of the Navier--Stokes equations
60 ///
61 ///
62 //======================================================================
63 template<unsigned DIM>
66 public virtual RefineableElement,
67 public virtual ElementWithZ2ErrorEstimator
68 {
69  protected:
70 
71  /// \short Unpin all pressure dofs in the element
72  virtual void unpin_elemental_pressure_dofs()=0;
73 
74  /// \short Pin unused nodal pressure dofs (empty by default, because
75  /// by default pressure dofs are not associated with nodes)
77 
78  public:
79 
80  /// \short Constructor
85 
86 
87  /// \short Loop over all elements in Vector (which typically contains
88  /// all the elements in a fluid mesh) and pin the nodal pressure degrees
89  /// of freedom that are not being used. Function uses
90  /// the member function
91  /// - \c RefineableNavierStokesEquations::
92  /// pin_elemental_redundant_nodal_pressure_dofs()
93  /// .
94  /// which is empty by default and should be implemented for
95  /// elements with nodal pressure degrees of freedom
96  /// (e.g. for refineable Taylor-Hood.)
98  element_pt)
99  {
100  // Loop over all elements and call the function that pins their
101  // unused nodal pressure data
102  unsigned n_element = element_pt.size();
103  for(unsigned e=0;e<n_element;e++)
104  {
106  (element_pt[e])->
108  }
109  }
110 
111  /// \short Unpin all pressure dofs in elements listed in vector.
113  element_pt)
114  {
115  // Loop over all elements
116  unsigned n_element = element_pt.size();
117  for(unsigned e=0;e<n_element;e++)
118  {
120  (element_pt[e])->
122  }
123  }
124 
125  /// \short Pointer to n_p-th pressure node (Default: NULL,
126  /// indicating that pressure is not based on nodal interpolation).
127  virtual Node* pressure_node_pt(const unsigned& n_p) {return NULL;}
128 
129  /// \short Compute the diagonal of the velocity/pressure mass matrices.
130  /// If which one=0, both are computed, otherwise only the pressure
131  /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
132  /// LSC version of the preconditioner only needs that one)
134  Vector<double> &press_mass_diag, Vector<double> &veloc_mass_diag,
135  const unsigned& which_one=0);
136 
137 
138  /// Number of 'flux' terms for Z2 error estimation
139  unsigned num_Z2_flux_terms()
140  {
141  // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
142  return DIM + (DIM*(DIM-1))/2;
143  }
144 
145  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
146  /// in strain rate tensor.
148  {
149 #ifdef PARANOID
150  unsigned num_entries=DIM+(DIM*(DIM-1))/2;
151  if (flux.size() < num_entries)
152  {
153  std::ostringstream error_message;
154  error_message << "The flux vector has the wrong number of entries, "
155  << flux.size() << ", whereas it should be at least "
156  << num_entries << std::endl;
157  throw OomphLibError(error_message.str(),
158  OOMPH_CURRENT_FUNCTION,
159  OOMPH_EXCEPTION_LOCATION);
160  }
161 #endif
162 
163  // Get strain rate matrix
164  DenseMatrix<double> strainrate(DIM);
165  this->strain_rate(s,strainrate);
166 
167  // Pack into flux Vector
168  unsigned icount=0;
169 
170  // Start with diagonal terms
171  for(unsigned i=0;i<DIM;i++)
172  {
173  flux[icount]=strainrate(i,i);
174  icount++;
175  }
176 
177  //Off diagonals row by row
178  for(unsigned i=0;i<DIM;i++)
179  {
180  for(unsigned j=i+1;j<DIM;j++)
181  {
182  flux[icount]=strainrate(i,j);
183  icount++;
184  }
185  }
186  }
187 
188  /// Further build, pass the pointers down to the sons
190  {
191  //Find the father element
193  cast_father_element_pt
195  (this->father_element_pt());
196 
197  //Set the viscosity ratio pointer
198  this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
199  //Set the density ratio pointer
200  this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
201  //Set pointer to global Reynolds number
202  this->Re_pt = cast_father_element_pt->re_pt();
203  //Set pointer to global Reynolds number x Strouhal number (=Womersley)
204  this->ReSt_pt = cast_father_element_pt->re_st_pt();
205  //Set pointer to global Reynolds number x inverse Froude number
206  this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
207  //Set pointer to global gravity Vector
208  this->G_pt = cast_father_element_pt->g_pt();
209 
210  //Set pointer to body force function
211  this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
212 
213  //Set pointer to volumetric source function
214  this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
215 
216  //Set the ALE flag
217  this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
218  }
219 
220 
221  /// \short Compute the derivatives of the i-th component of
222  /// velocity at point s with respect
223  /// to all data that can affect its value. In addition, return the global
224  /// equation numbers corresponding to the data.
225  /// Overload the non-refineable version to take account of hanging node
226  /// information
228  const unsigned &i,
229  Vector<double> &du_ddata,
230  Vector<unsigned> &global_eqn_number)
231  {
232  //Find number of nodes
233  unsigned n_node = this->nnode();
234  //Local shape function
235  Shape psi(n_node);
236  //Find values of shape function at the given local coordinate
237  this->shape(s,psi);
238 
239  //Find the index at which the velocity component is stored
240  const unsigned u_nodal_index = this->u_index_nst(i);
241 
242  //Storage for hang info pointer
243  HangInfo* hang_info_pt=0;
244  //Storage for global equation
245  int global_eqn = 0;
246 
247  //Find the number of dofs associated with interpolated u
248  unsigned n_u_dof=0;
249  for(unsigned l=0;l<n_node;l++)
250  {
251  unsigned n_master = 1;
252 
253  //Local bool (is the node hanging)
254  bool is_node_hanging = this->node_pt(l)->is_hanging();
255 
256  //If the node is hanging, get the number of master nodes
257  if(is_node_hanging)
258  {
259  hang_info_pt = this->node_pt(l)->hanging_pt();
260  n_master = hang_info_pt->nmaster();
261  }
262  //Otherwise there is just one master node, the node itself
263  else
264  {
265  n_master = 1;
266  }
267 
268  //Loop over the master nodes
269  for(unsigned m=0;m<n_master;m++)
270  {
271  //Get the equation number
272  if(is_node_hanging)
273  {
274  //Get the equation number from the master node
275  global_eqn = hang_info_pt->master_node_pt(m)->
276  eqn_number(u_nodal_index);
277  }
278  else
279  {
280  // Global equation number
281  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
282  }
283 
284  //If it's positive add to the count
285  if(global_eqn >= 0) {++n_u_dof;}
286  }
287  }
288 
289  //Now resize the storage schemes
290  du_ddata.resize(n_u_dof,0.0);
291  global_eqn_number.resize(n_u_dof,0);
292 
293  //Loop over th nodes again and set the derivatives
294  unsigned count=0;
295  //Loop over the local nodes and sum
296  for(unsigned l=0;l<n_node;l++)
297  {
298  unsigned n_master = 1;
299  double hang_weight = 1.0;
300 
301  //Local bool (is the node hanging)
302  bool is_node_hanging = this->node_pt(l)->is_hanging();
303 
304  //If the node is hanging, get the number of master nodes
305  if(is_node_hanging)
306  {
307  hang_info_pt = this->node_pt(l)->hanging_pt();
308  n_master = hang_info_pt->nmaster();
309  }
310  //Otherwise there is just one master node, the node itself
311  else
312  {
313  n_master = 1;
314  }
315 
316  //Loop over the master nodes
317  for(unsigned m=0;m<n_master;m++)
318  {
319  //If the node is hanging get weight from master node
320  if(is_node_hanging)
321  {
322  //Get the hang weight from the master node
323  hang_weight = hang_info_pt->master_weight(m);
324  }
325  else
326  {
327  // Node contributes with full weight
328  hang_weight = 1.0;
329  }
330 
331  //Get the equation number
332  if(is_node_hanging)
333  {
334  //Get the equation number from the master node
335  global_eqn = hang_info_pt->master_node_pt(m)->
336  eqn_number(u_nodal_index);
337  }
338  else
339  {
340  // Local equation number
341  global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
342  }
343 
344  if (global_eqn >= 0)
345  {
346  //Set the global equation number
347  global_eqn_number[count] = global_eqn;
348  //Set the derivative with respect to the unknown
349  du_ddata[count] = psi[l]*hang_weight;
350  //Increase the counter
351  ++count;
352  }
353  }
354  }
355  }
356 
357 
358 
359  protected:
360 
361 
362 /// \short Add element's contribution to elemental residual vector and/or
363 /// Jacobian matrix
364 /// flag=1: compute both
365 /// flag=0: compute only residual vector
367  Vector<double> &residuals,
368  DenseMatrix<double> &jacobian,
369  DenseMatrix<double> &mass_matrix,
370  unsigned flag);
371 
372  /// \short Compute derivatives of elemental residual vector with respect
373  /// to nodal coordinates. Overwrites default implementation in
374  /// FiniteElement base class.
375  /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
377  dresidual_dnodal_coordinates);
378 
379 };
380 
381 
382 //======================================================================
383 /// Refineable version of Taylor Hood elements. These classes
384 /// can be written in total generality.
385 //======================================================================
386 template<unsigned DIM>
390 public virtual RefineableQElement<DIM>
391 {
392  private:
393 
394  /// Unpin all pressure dofs
396  {
397  //find the index at which the pressure is stored
398  int p_index = this->p_nodal_index_nst();
399  unsigned n_node = this->nnode();
400  // loop over nodes
401  for(unsigned n=0;n<n_node;n++)
402  {this->node_pt(n)->unpin(p_index);}
403  }
404 
405  /// Pin all nodal pressure dofs that are not required
407  {
408  //Find the pressure index
409  int p_index = this->p_nodal_index_nst();
410  //Loop over all nodes
411  unsigned n_node = this->nnode();
412  // loop over all nodes and pin all the nodal pressures
413  for(unsigned n=0;n<n_node;n++) {this->node_pt(n)->pin(p_index);}
414 
415  // Loop over all actual pressure nodes and unpin if they're not hanging
416  unsigned n_pres = this->npres_nst();
417  for(unsigned l=0;l<n_pres;l++)
418  {
419  Node* nod_pt = this->pressure_node_pt(l);
420  if (!nod_pt->is_hanging(p_index)) {nod_pt->unpin(p_index);}
421  }
422  }
423 
424  public:
425 
426  /// \short Constructor
430  RefineableQElement<DIM>(),
432 
433  /// \short Number of values required at local node n. In order to simplify
434  /// matters, we allocate storage for pressure variables at all the nodes
435  /// and then pin those that are not used.
436  unsigned required_nvalue(const unsigned &n) const {return DIM+1;}
437 
438  /// Number of continuously interpolated values: (DIM velocities + 1 pressure)
439  unsigned ncont_interpolated_values() const {return DIM+1;}
440 
441  /// Rebuild from sons: empty
442  void rebuild_from_sons(Mesh* &mesh_pt) {}
443 
444  /// \short Order of recovery shape functions for Z2 error estimation:
445  /// Same order as shape functions.
446  unsigned nrecovery_order() {return 2;}
447 
448  /// \short Number of vertex nodes in the element
449  unsigned nvertex_node() const
451 
452  /// \short Pointer to the j-th vertex node in the element
453  Node* vertex_node_pt(const unsigned& j) const
455 
456 /// \short Get the function value u in Vector.
457 /// Note: Given the generality of the interface (this function
458 /// is usually called from black-box documentation or interpolation routines),
459 /// the values Vector sets its own size in here.
461  {
462  // Set size of Vector: u,v,p and initialise to zero
463  values.resize(DIM+1,0.0);
464 
465  //Calculate velocities: values[0],...
466  for(unsigned i=0;i<DIM;i++) {values[i] = this->interpolated_u_nst(s,i);}
467 
468  //Calculate pressure: values[DIM]
469  values[DIM] = this->interpolated_p_nst(s);
470  }
471 
472 /// \short Get the function value u in Vector.
473 /// Note: Given the generality of the interface (this function
474 /// is usually called from black-box documentation or interpolation routines),
475 /// the values Vector sets its own size in here.
476  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
477  Vector<double>& values)
478  {
479  // Set size of Vector: u,v,p
480  values.resize(DIM+1);
481 
482  // Initialise
483  for(unsigned i=0;i<DIM+1;i++) {values[i]=0.0;}
484 
485  //Find out how many nodes there are
486  unsigned n_node = this->nnode();
487 
488  // Shape functions
489  Shape psif(n_node);
490  this->shape(s,psif);
491 
492  //Calculate velocities: values[0],...
493  for(unsigned i=0;i<DIM;i++)
494  {
495  //Get the index at which the i-th velocity is stored
496  unsigned u_nodal_index = this->u_index_nst(i);
497  for(unsigned l=0;l<n_node;l++)
498  {
499  values[i] += this->nodal_value(t,l,u_nodal_index)*psif[l];
500  }
501  }
502 
503  //Calculate pressure: values[DIM]
504  //(no history is carried in the pressure)
505  values[DIM] = this->interpolated_p_nst(s);
506  }
507 
508 
509  /// \short Perform additional hanging node procedures for variables
510  /// that are not interpolated by all nodes. The pressures are stored
511  /// at the p_nodal_index_nst-th location in each node
513  {
514  this->setup_hang_for_value(this->p_nodal_index_nst());
515  }
516 
517  /// \short Pointer to n_p-th pressure node
518  Node* pressure_node_pt(const unsigned &n_p)
519  {return this->node_pt(this->Pconv[n_p]);}
520 
521  /// \short The velocities are isoparametric and so the "nodes" interpolating
522  /// the velocities are the geometric nodes. The pressure "nodes" are a
523  /// subset of the nodes, so when value_id==DIM, the n-th pressure
524  /// node is returned.
525  Node* interpolating_node_pt(const unsigned &n,
526  const int &value_id)
527 
528  {
529  //The only different nodes are the pressure nodes
530  if(value_id==DIM) {return this->pressure_node_pt(n);}
531  //The other variables are interpolated via the usual nodes
532  else {return this->node_pt(n);}
533  }
534 
535  /// \short The pressure nodes are the corner nodes, so when n_value==DIM,
536  /// the fraction is the same as the 1d node number, 0 or 1.
537  double local_one_d_fraction_of_interpolating_node(const unsigned &n1d,
538  const unsigned &i,
539  const int &value_id)
540  {
541  if(value_id==DIM)
542  {
543  //The pressure nodes are just located on the boundaries at 0 or 1
544  return double(n1d);
545  }
546  //Otherwise the velocity nodes are the same as the geometric ones
547  else
548  {
549  return this->local_one_d_fraction_of_node(n1d,i);
550  }
551  }
552 
553  /// \short The velocity nodes are the same as the geometric nodes. The
554  /// pressure nodes must be calculated by using the same methods as
555  /// the geometric nodes, but by recalling that there are only two pressure
556  /// nodes per edge.
558  const int &value_id)
559  {
560  //If we are calculating pressure nodes
561  if(value_id==DIM)
562  {
563  //Storage for the index of the pressure node
564  unsigned total_index=0;
565  //The number of nodes along each 1d edge is 2.
566  unsigned NNODE_1D = 2;
567  //Storage for the index along each boundary
568  Vector<int> index(DIM);
569  //Loop over the coordinates
570  for(unsigned i=0;i<DIM;i++)
571  {
572  //If we are at the lower limit, the index is zero
573  if(s[i]==-1.0)
574  {
575  index[i] = 0;
576  }
577  //If we are at the upper limit, the index is the number of nodes minus 1
578  else if(s[i] == 1.0)
579  {
580  index[i] = NNODE_1D-1;
581  }
582  //Otherwise, we have to calculate the index in general
583  else
584  {
585  //For uniformly spaced nodes the 0th node number would be
586  double float_index = 0.5*(1.0 + s[i])*(NNODE_1D-1);
587  index[i] = int(float_index);
588  //What is the excess. This should be safe because the
589  //taking the integer part rounds down
590  double excess = float_index - index[i];
591  //If the excess is bigger than our tolerance there is no node,
592  //return null
593  if((excess > FiniteElement::Node_location_tolerance) && (
594  (1.0 - excess) > FiniteElement::Node_location_tolerance))
595  {
596  return 0;
597  }
598  }
599  ///Construct the general pressure index from the components.
600  total_index += index[i]*
601  static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
602  static_cast<int>(i)));
603  }
604  //If we've got here we have a node, so let's return a pointer to it
605  return this->pressure_node_pt(total_index);
606  }
607  //Otherwise velocity nodes are the same as pressure nodes
608  else
609  {
610  return this->get_node_at_local_coordinate(s);
611  }
612  }
613 
614 
615  /// \short The number of 1d pressure nodes is 2, the number of 1d velocity
616  /// nodes is the same as the number of 1d geometric nodes.
617  unsigned ninterpolating_node_1d(const int &value_id)
618  {
619  if(value_id==DIM) {return 2;}
620  else {return this->nnode_1d();}
621  }
622 
623  /// \short The number of pressure nodes is 2^DIM. The number of
624  /// velocity nodes is the same as the number of geometric nodes.
625  unsigned ninterpolating_node(const int &value_id)
626  {
627  if(value_id==DIM)
628  {return static_cast<unsigned>(pow(2.0,static_cast<int>(DIM)));}
629  else {return this->nnode();}
630  }
631 
632  /// \short The basis interpolating the pressure is given by pshape().
633  //// The basis interpolating the velocity is shape().
635  Shape &psi,
636  const int &value_id) const
637  {
638  if(value_id==DIM) {return this->pshape_nst(s,psi);}
639  else {return this->shape(s,psi);}
640  }
641 
642 
643  /// \short Add to the set \c paired_load_data pairs containing
644  /// - the pointer to a Data object
645  /// and
646  /// - the index of the value in that Data object
647  /// .
648  /// for all values (pressures, velocities) that affect the
649  /// load computed in the \c get_load(...) function.
650  /// (Overloads non-refineable version and takes hanging nodes
651  /// into account)
653  std::set<std::pair<Data*,unsigned> > &paired_load_data)
654  {
655  //Get the nodal indices at which the velocities are stored
656  unsigned u_index[DIM];
657  for(unsigned i=0;i<DIM;i++) {u_index[i] = this->u_index_nst(i);}
658 
659  //Loop over the nodes
660  unsigned n_node = this->nnode();
661  for(unsigned n=0;n<n_node;n++)
662  {
663  // Pointer to current node
664  Node* nod_pt=this->node_pt(n);
665 
666  // Check if it's hanging:
667  if (nod_pt->is_hanging())
668  {
669  // It's hanging -- get number of master nodes
670  unsigned nmaster=nod_pt->hanging_pt()->nmaster();
671 
672  // Loop over masters
673  for (unsigned j=0;j<nmaster;j++)
674  {
675  Node* master_nod_pt=nod_pt->hanging_pt()->master_node_pt(j);
676 
677  //Loop over the velocity components and add pointer to their data
678  //and indices to the vectors
679  for(unsigned i=0;i<DIM;i++)
680  {
681  paired_load_data.insert(std::make_pair(master_nod_pt,u_index[i]));
682  }
683  }
684  }
685  //Not hanging
686  else
687  {
688  //Loop over the velocity components and add pointer to their data
689  //and indices to the vectors
690  for(unsigned i=0;i<DIM;i++)
691  {
692  paired_load_data.insert(std::make_pair(this->node_pt(n),u_index[i]));
693  }
694  }
695  }
696 
697  //Get the nodal index at which the pressure is stored
698  int p_index = this->p_nodal_index_nst();
699 
700  //Loop over the pressure data
701  unsigned n_pres= this->npres_nst();
702  for(unsigned l=0;l<n_pres;l++)
703  {
704  //Get the pointer to the nodal pressure
705  Node* pres_node_pt = this->pressure_node_pt(l);
706  // Check if the pressure dof is hanging
707  if(pres_node_pt->is_hanging(p_index))
708  {
709  //Get the pointer to the hang info object
710  // (pressure is stored as p_index--th nodal dof).
711  HangInfo* hang_info_pt = pres_node_pt->hanging_pt(p_index);
712 
713  // Get number of pressure master nodes (pressure is stored
714  unsigned nmaster = hang_info_pt->nmaster();
715 
716  // Loop over pressure master nodes
717  for(unsigned m=0;m<nmaster;m++)
718  {
719  //The p_index-th entry in each nodal data is the pressure, which
720  //affects the traction
721  paired_load_data.insert(
722  std::make_pair(hang_info_pt->master_node_pt(m),p_index));
723  }
724  }
725  // It's not hanging
726  else
727  {
728  //The p_index-th entry in each nodal data is the pressure, which
729  //affects the traction
730  paired_load_data.insert(std::make_pair(pres_node_pt,p_index));
731  }
732  }
733 
734  }
735 
736 
737 };
738 
739 
740 //=======================================================================
741 /// \short Face geometry of the RefineableQTaylorHoodElements is the
742 /// same as the Face geometry of the QTaylorHoodElements.
743 //=======================================================================
744 template<unsigned DIM>
746 public virtual FaceGeometry<GeneralisedNewtonianQTaylorHoodElement<DIM> >
747 {
748  public:
750 };
751 
752 
753 //=======================================================================
754 /// \short Face geometry of the face geometry of
755 /// the RefineableQTaylorHoodElements is the
756 /// same as the Face geometry of the Face geometry of QTaylorHoodElements.
757 //=======================================================================
758 template<unsigned DIM>
761 public virtual FaceGeometry<FaceGeometry<
762 GeneralisedNewtonianQTaylorHoodElement<DIM> > >
763 {
764  public:
767  {}
768 };
769 
770 
771 ///////////////////////////////////////////////////////////////////////////
772 ///////////////////////////////////////////////////////////////////////////
773 ///////////////////////////////////////////////////////////////////////////
774 
775 
776 //======================================================================
777 /// Refineable version of Crouzeix Raviart elements. Generic class definitions
778 //======================================================================
779 template<unsigned DIM>
783 public virtual RefineableQElement<DIM>
784 {
785  private:
786 
787  /// Unpin all internal pressure dofs
789  {
790  unsigned n_pres = this->npres_nst();
791  // loop over pressure dofs and unpin them
792  for(unsigned l=0;l<n_pres;l++)
793  {this->internal_data_pt(this->P_nst_internal_index)->unpin(l);}
794  }
795 
796  public:
797 
798  /// \short Constructor
802  RefineableQElement<DIM>(),
804 
805 
806  /// Broken copy constructor
809  {
811  "RefineableGeneralisedNewtonianQCrouzeixRaviartElement");
812  }
813 
814  /// Broken assignment operator
815 //Commented out broken assignment operator because this can lead to a conflict warning
816 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
817 //realise that two separate implementations of the broken function are the same and so,
818 //quite rightly, it shouts.
819  /*void operator=(const
820  RefineableGeneralisedNewtonianQCrouzeixRaviartElement<DIM>&)
821  {
822  BrokenCopy::broken_assign(
823  "RefineableGeneralisedNewtonianQCrouzeixRaviartElement");
824  }*/
825 
826  /// Number of continuously interpolated values: DIM (velocities)
827  unsigned ncont_interpolated_values() const {return DIM;}
828 
829  /// \short Rebuild from sons: Reconstruct pressure from the (merged) sons
830  /// This must be specialised for each dimension.
831  inline void rebuild_from_sons(Mesh* &mesh_pt);
832 
833  /// \short Order of recovery shape functions for Z2 error estimation:
834  /// Same order as shape functions.
835  unsigned nrecovery_order() {return 2;}
836 
837  /// \short Number of vertex nodes in the element
838  unsigned nvertex_node() const
840 
841  /// \short Pointer to the j-th vertex node in the element
842  Node* vertex_node_pt(const unsigned& j) const
843  {
845  }
846 
847 /// \short Get the function value u in Vector.
848 /// Note: Given the generality of the interface (this function
849 /// is usually called from black-box documentation or interpolation routines),
850 /// the values Vector sets its own size in here.
852  {
853  // Set size of Vector: u,v,p and initialise to zero
854  values.resize(DIM,0.0);
855 
856  //Calculate velocities: values[0],...
857  for(unsigned i=0;i<DIM;i++) {values[i] = this->interpolated_u_nst(s,i);}
858  }
859 
860  /// \short Get all function values [u,v..,p] at previous timestep t
861  /// (t=0: present; t>0: previous timestep).
862  ///
863  /// Note: Given the generality of the interface (this function
864  /// is usually called from black-box documentation or interpolation routines),
865  /// the values Vector sets its own size in here.
866  ///
867  /// Note: No pressure history is kept, so pressure is always
868  /// the current value.
869  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
870  Vector<double>& values)
871  {
872  // Set size of Vector: u,v,p
873  values.resize(DIM);
874 
875  // Initialise
876  for(unsigned i=0;i<DIM;i++) {values[i]=0.0;}
877 
878  //Find out how many nodes there are
879  unsigned n_node = this->nnode();
880 
881  // Shape functions
882  Shape psif(n_node);
883  this->shape(s,psif);
884 
885  //Calculate velocities: values[0],...
886  for(unsigned i=0;i<DIM;i++)
887  {
888  //Get the nodal index at which the i-th velocity component is stored
889  unsigned u_nodal_index = this->u_index_nst(i);
890  for(unsigned l=0;l<n_node;l++)
891  {
892  values[i] += this->nodal_value(t,l,u_nodal_index)*psif[l];
893  }
894  }
895  }
896 
897  /// \short Perform additional hanging node procedures for variables
898  /// that are not interpolated by all nodes. Empty
900 
901  /// Further build for Crouzeix_Raviart interpolates the internal
902  /// pressure dofs from father element: Make sure pressure values and
903  /// dp/ds agree between fathers and sons at the midpoints of the son
904  /// elements. This must be specialised for each dimension.
905  inline void further_build();
906 
907 
908  /// \short Add to the set \c paired_load_data pairs containing
909  /// - the pointer to a Data object
910  /// and
911  /// - the index of the value in that Data object
912  /// .
913  /// for all values (pressures, velocities) that affect the
914  /// load computed in the \c get_load(...) function.
915  /// (Overloads non-refineable version and takes hanging nodes
916  /// into account)
918  std::set<std::pair<Data*,unsigned> > &paired_load_data)
919  {
920  //Get the nodal indices at which the velocities are stored
921  unsigned u_index[DIM];
922  for(unsigned i=0;i<DIM;i++) {u_index[i] = this->u_index_nst(i);}
923 
924  //Loop over the nodes
925  unsigned n_node = this->nnode();
926  for(unsigned n=0;n<n_node;n++)
927  {
928  // Pointer to current node
929  Node* nod_pt=this->node_pt(n);
930 
931  // Check if it's hanging:
932  if (nod_pt->is_hanging())
933  {
934  // It's hanging -- get number of master nodes
935  unsigned nmaster=nod_pt->hanging_pt()->nmaster();
936 
937  // Loop over masters
938  for (unsigned j=0;j<nmaster;j++)
939  {
940  Node* master_nod_pt=nod_pt->hanging_pt()->master_node_pt(j);
941 
942  //Loop over the velocity components and add pointer to their data
943  //and indices to the vectors
944  for(unsigned i=0;i<DIM;i++)
945  {
946  paired_load_data.insert(std::make_pair(master_nod_pt,u_index[i]));
947  }
948  }
949  }
950  //Not hanging
951  else
952  {
953  //Loop over the velocity components and add pointer to their data
954  //and indices to the vectors
955  for(unsigned i=0;i<DIM;i++)
956  {
957  paired_load_data.insert(std::make_pair(this->node_pt(n),u_index[i]));
958  }
959  }
960  }
961 
962 
963  //Loop over the pressure data (can't be hanging!)
964  unsigned n_pres = this->npres_nst();
965  for(unsigned l=0;l<n_pres;l++)
966  {
967  //The entries in the internal data at P_nst_internal_index
968  //are the pressures, which affect the traction
969  paired_load_data.insert(
970  std::make_pair(this->internal_data_pt(this->P_nst_internal_index),l));
971  }
972  }
973 
974 
975 };
976 
977 
978 //======================================================================
979 /// p-refineable version of Crouzeix Raviart elements. Generic class
980 /// definitions
981 //======================================================================
982 template<unsigned DIM>
986 public virtual PRefineableQElement<DIM,3>
987 {
988  private:
989 
990  /// Unpin all internal pressure dofs
992  {
993  unsigned n_pres = this->npres_nst();
994  n_pres = this->internal_data_pt(this->P_nst_internal_index)->nvalue();
995  // loop over pressure dofs and unpin them
996  for(unsigned l=0;l<n_pres;l++)
997  {this->internal_data_pt(this->P_nst_internal_index)->unpin(l);}
998  }
999 
1000  public:
1001 
1002  /// \short Constructor
1006  PRefineableQElement<DIM,3>(),
1008  {
1009  // Set the p-order
1010  this->p_order()=3;
1011 
1012  // Set integration scheme
1013  // (To avoid memory leaks in pre-build and p-refine where new
1014  // integration schemes are created)
1016 
1017  // Resize pressure storage
1018  // (Constructor for QCrouzeixRaviartElement sets up DIM+1 pressure values)
1019  if (this->internal_data_pt(this->P_nst_internal_index)->nvalue()
1020  <= this->npres_nst())
1021  {
1023  ->resize(this->npres_nst());
1024  }
1025  else
1026  {
1027  Data* new_data_pt = new Data(this->npres_nst());
1028  delete this->internal_data_pt(this->P_nst_internal_index);
1029  this->internal_data_pt(this->P_nst_internal_index) = new_data_pt;
1030  }
1031  }
1032 
1033  /// \short Destructor
1035  {
1036  delete this->integral_pt();
1037  }
1038 
1039 
1040  /// Broken copy constructor
1043  {
1045  "PRefineableGeneralisedNewtonianQCrouzeixRaviartElement");
1046  }
1047 
1048  /// Broken assignment operator
1049  /*void operator=(const
1050  PRefineableGeneralisedNewtonianQCrouzeixRaviartElement<DIM>&)
1051  {
1052  BrokenCopy::broken_assign(
1053  "PRefineableGeneralisedNewtonianQCrouzeixRaviartElement");
1054  }*/
1055 
1056  /// \short Return the i-th pressure value
1057  /// (Discontinous pressure interpolation -- no need to cater for hanging
1058  /// nodes).
1059  double p_nst(const unsigned &i) const
1060  {return this->internal_data_pt(this->P_nst_internal_index)->value(i);}
1061 
1062  double p_nst(const unsigned &t, const unsigned &i) const
1063  {return this->internal_data_pt(this->P_nst_internal_index)->value(t,i);}
1064 
1065  ///// Return number of pressure values
1066  unsigned npres_nst() const {return (this->p_order()-2)*(this->p_order()-2);}
1067 
1068  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
1069  void fix_pressure(const unsigned &p_dof, const double &p_value)
1070  {
1071  this->internal_data_pt(this->P_nst_internal_index)->pin(p_dof);
1072  this->internal_data_pt(this->P_nst_internal_index)->set_value(p_dof,p_value);
1073  }
1074 
1075  unsigned required_nvalue(const unsigned& n) const
1076  {return DIM;}
1077 
1078  /// Number of continuously interpolated values: DIM (velocities)
1079  unsigned ncont_interpolated_values() const {return DIM;}
1080 
1081  /// \short Rebuild from sons: Reconstruct pressure from the (merged) sons
1082  /// This must be specialised for each dimension.
1083  void rebuild_from_sons(Mesh* &mesh_pt)
1084  {
1085  //Do p-refineable version
1087  //Do Crouzeix-Raviart version
1088  //Need to reconstruct pressure manually!
1089  for(unsigned p=0; p<npres_nst(); p++)
1090  {
1091  //BENFLAG: Set to zero for now -- don't do projection problem yet
1092  this->internal_data_pt(this->P_nst_internal_index)->set_value(p,0.0);
1093  }
1094 
1095  }
1096 
1097  /// \short Order of recovery shape functions for Z2 error estimation:
1098  /// - Same order as shape functions.
1099  //unsigned nrecovery_order()
1100  // {
1101  // if(this->nnode_1d() < 4) {return (this->nnode_1d()-1);}
1102  // else {return 3;}
1103  // }
1104  /// - Constant recovery order, since recovery order of the first element
1105  /// is used for the whole mesh.
1106  unsigned nrecovery_order() {return 3;}
1107 
1108  /// \short Number of vertex nodes in the element
1109  unsigned nvertex_node() const
1111 
1112  /// \short Pointer to the j-th vertex node in the element
1113  Node* vertex_node_pt(const unsigned& j) const
1114  {
1116  }
1117 
1118  /// \short Velocity shape and test functions and their derivs
1119  /// w.r.t. to global coords at local coordinate s (taken from geometry)
1120  ///Return Jacobian of mapping between local and global coordinates.
1121  inline double dshape_and_dtest_eulerian_nst(const Vector<double> &s,
1122  Shape &psi,
1123  DShape &dpsidx, Shape &test,
1124  DShape &dtestdx) const;
1125 
1126  /// \short Velocity shape and test functions and their derivs
1127  /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
1128  ///Return Jacobian of mapping between local and global coordinates.
1129  inline double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt,
1130  Shape &psi,
1131  DShape &dpsidx,
1132  Shape &test,
1133  DShape &dtestdx) const;
1134 
1135  /// Pressure shape functions at local coordinate s
1136  inline void pshape_nst(const Vector<double> &s, Shape &psi) const;
1137 
1138  /// Pressure shape and test functions at local coordinte s
1139  inline void pshape_nst(const Vector<double> &s, Shape &psi, Shape &test) const;
1140 
1141  /// \short Get the function value u in Vector.
1142  /// Note: Given the generality of the interface (this function
1143  /// is usually called from black-box documentation or interpolation routines),
1144  /// the values Vector sets its own size in here.
1146  {
1147  // Set size of Vector: u,v,p and initialise to zero
1148  values.resize(DIM,0.0);
1149 
1150  //Calculate velocities: values[0],...
1151  for(unsigned i=0;i<DIM;i++) {values[i] = this->interpolated_u_nst(s,i);}
1152  }
1153 
1154  /// \short Get all function values [u,v..,p] at previous timestep t
1155  /// (t=0: present; t>0: previous timestep).
1156  ///
1157  /// Note: Given the generality of the interface (this function
1158  /// is usually called from black-box documentation or interpolation routines),
1159  /// the values Vector sets its own size in here.
1160  ///
1161  /// Note: No pressure history is kept, so pressure is always
1162  /// the current value.
1163  void get_interpolated_values(const unsigned& t, const Vector<double>&s,
1164  Vector<double>& values)
1165  {
1166  // Set size of Vector: u,v,p
1167  values.resize(DIM);
1168 
1169  // Initialise
1170  for(unsigned i=0;i<DIM;i++) {values[i]=0.0;}
1171 
1172  //Find out how many nodes there are
1173  unsigned n_node = this->nnode();
1174 
1175  // Shape functions
1176  Shape psif(n_node);
1177  this->shape(s,psif);
1178 
1179  //Calculate velocities: values[0],...
1180  for(unsigned i=0;i<DIM;i++)
1181  {
1182  //Get the nodal index at which the i-th velocity component is stored
1183  unsigned u_nodal_index = this->u_index_nst(i);
1184  for(unsigned l=0;l<n_node;l++)
1185  {
1186  values[i] += this->nodal_value(t,l,u_nodal_index)*psif[l];
1187  }
1188  }
1189  }
1190 
1191  /// \short Perform additional hanging node procedures for variables
1192  /// that are not interpolated by all nodes. Empty
1194 
1195  /// Further build for Crouzeix_Raviart interpolates the internal
1196  /// pressure dofs from father element: Make sure pressure values and
1197  /// dp/ds agree between fathers and sons at the midpoints of the son
1198  /// elements. This must be specialised for each dimension.
1199  void further_build();
1200 
1201 
1202 };
1203 
1204 
1205 //=======================================================================
1206 /// Face geometry of the RefineableQuadQCrouzeixRaviartElements
1207 //=======================================================================
1208 template<unsigned DIM>
1211 public virtual FaceGeometry<GeneralisedNewtonianQCrouzeixRaviartElement<DIM> >
1212 {
1213  public:
1216 };
1217 
1218 //======================================================================
1219 /// \short Face geometry of the face geometry of
1220 /// the RefineableQCrouzeixRaviartElements is the
1221 /// same as the Face geometry of the Face geometry of
1222 /// QCrouzeixRaviartElements.
1223 //=======================================================================
1224 template<unsigned DIM>
1227 public virtual FaceGeometry<
1228 FaceGeometry<GeneralisedNewtonianQCrouzeixRaviartElement<DIM> > >
1229 {
1230  public:
1233  {}
1234 };
1235 
1236 
1237 //Inline functions
1238 
1239 //=====================================================================
1240 /// 2D Rebuild from sons: Reconstruct pressure from the (merged) sons
1241 //=====================================================================
1242 template<>
1245 {
1246  using namespace QuadTreeNames;
1247 
1248  //Central pressure value:
1249  //-----------------------
1250 
1251  // Use average of the sons central pressure values
1252  // Other options: Take average of the four (discontinuous)
1253  // pressure values at the father's midpoint]
1254 
1255  double av_press=0.0;
1256 
1257  // Loop over the sons
1258  for (unsigned ison=0;ison<4;ison++)
1259  {
1260  // Add the sons midnode pressure
1261  // Note that we can assume that the pressure is stored at the same
1262  // location because these are EXACTLY the same type of elements
1263  av_press += quadtree_pt()->son_pt(ison)->object_pt()->
1264  internal_data_pt(this->P_nst_internal_index)->value(0);
1265  }
1266 
1267  // Use the average
1268  internal_data_pt(this->P_nst_internal_index)->set_value(0,0.25*av_press);
1269 
1270 
1271  //Slope in s_0 direction
1272  //----------------------
1273 
1274  // Use average of the 2 FD approximations based on the
1275  // elements central pressure values
1276  // [Other options: Take average of the four
1277  // pressure derivatives]
1278 
1279  double slope1=
1280  quadtree_pt()->son_pt(SE)->object_pt()->
1281  internal_data_pt(this->P_nst_internal_index)->value(0) -
1282  quadtree_pt()->son_pt(SW)->object_pt()->
1283  internal_data_pt(this->P_nst_internal_index)->value(0);
1284 
1285  double slope2 =
1286  quadtree_pt()->son_pt(NE)->object_pt()->
1287  internal_data_pt(this->P_nst_internal_index)->value(0) -
1288  quadtree_pt()->son_pt(NW)->object_pt()->
1289  internal_data_pt(this->P_nst_internal_index)->value(0);
1290 
1291 
1292  // Use the average
1293  internal_data_pt(this->P_nst_internal_index)->
1294  set_value(1,0.5*(slope1+slope2));
1295 
1296 
1297 
1298  //Slope in s_1 direction
1299  //----------------------
1300 
1301  // Use average of the 2 FD approximations based on the
1302  // elements central pressure values
1303  // [Other options: Take average of the four
1304  // pressure derivatives]
1305 
1306  slope1 =
1307  quadtree_pt()->son_pt(NE)->object_pt()
1308  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1309  quadtree_pt()->son_pt(SE)->object_pt()
1310  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1311 
1312  slope2=
1313  quadtree_pt()->son_pt(NW)->object_pt()
1314  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1315  quadtree_pt()->son_pt(SW)->object_pt()
1316  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1317 
1318 
1319  // Use the average
1320  internal_data_pt(this->P_nst_internal_index)->
1321  set_value(2,0.5*(slope1+slope2));
1322 }
1323 
1324 
1325 
1326 //=================================================================
1327 /// 3D Rebuild from sons: Reconstruct pressure from the (merged) sons
1328 //=================================================================
1329 template<>
1332 {
1333  using namespace OcTreeNames;
1334 
1335  //Central pressure value:
1336  //-----------------------
1337 
1338  // Use average of the sons central pressure values
1339  // Other options: Take average of the four (discontinuous)
1340  // pressure values at the father's midpoint]
1341 
1342  double av_press=0.0;
1343 
1344  // Loop over the sons
1345  for (unsigned ison=0;ison<8;ison++)
1346  {
1347  // Add the sons midnode pressure
1348  av_press += octree_pt()->son_pt(ison)->object_pt()->
1349  internal_data_pt(this->P_nst_internal_index)->value(0);
1350  }
1351 
1352  // Use the average
1353  internal_data_pt(this->P_nst_internal_index)->
1354  set_value(0,0.125*av_press);
1355 
1356 
1357  //Slope in s_0 direction
1358  //----------------------
1359 
1360  // Use average of the 4 FD approximations based on the
1361  // elements central pressure values
1362  // [Other options: Take average of the four
1363  // pressure derivatives]
1364 
1365  double slope1 =
1366  octree_pt()->son_pt(RDF)->object_pt()->
1367  internal_data_pt(this->P_nst_internal_index)->value(0) -
1368  octree_pt()->son_pt(LDF)->object_pt()->
1369  internal_data_pt(this->P_nst_internal_index)->value(0);
1370 
1371  double slope2 =
1372  octree_pt()->son_pt(RUF)->object_pt()->
1373  internal_data_pt(this->P_nst_internal_index)->value(0) -
1374  octree_pt()->son_pt(LUF)->object_pt()->
1375  internal_data_pt(this->P_nst_internal_index)->value(0);
1376 
1377  double slope3 =
1378  octree_pt()->son_pt(RDB)->object_pt()->
1379  internal_data_pt(this->P_nst_internal_index)->value(0) -
1380  octree_pt()->son_pt(LDB)->object_pt()->
1381  internal_data_pt(this->P_nst_internal_index)->value(0);
1382 
1383  double slope4 =
1384  octree_pt()->son_pt(RUB)->object_pt()->
1385  internal_data_pt(this->P_nst_internal_index)->value(0) -
1386  octree_pt()->son_pt(LUB)->object_pt()->
1387  internal_data_pt(this->P_nst_internal_index)->value(0);
1388 
1389 
1390  // Use the average
1391  internal_data_pt(this->P_nst_internal_index)->
1392  set_value(1,0.25*(slope1+slope2+slope3+slope4));
1393 
1394 
1395  //Slope in s_1 direction
1396  //----------------------
1397 
1398  // Use average of the 4 FD approximations based on the
1399  // elements central pressure values
1400  // [Other options: Take average of the four
1401  // pressure derivatives]
1402 
1403  slope1 =
1404  octree_pt()->son_pt(LUB)->object_pt()
1405  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1406  octree_pt()->son_pt(LDB)->object_pt()
1407  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1408 
1409  slope2 =
1410  octree_pt()->son_pt(RUB)->object_pt()
1411  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1412  octree_pt()->son_pt(RDB)->object_pt()
1413  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1414 
1415  slope3 =
1416  octree_pt()->son_pt(LUF)->object_pt()
1417  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1418  octree_pt()->son_pt(LDF)->object_pt()
1419  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1420 
1421  slope4 =
1422  octree_pt()->son_pt(RUF)->object_pt()
1423  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1424  octree_pt()->son_pt(RDF)->object_pt()
1425  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1426 
1427 
1428  // Use the average
1429  internal_data_pt(this->P_nst_internal_index)->
1430  set_value(2,0.25*(slope1+slope2+slope3+slope4));
1431 
1432 
1433  //Slope in s_2 direction
1434  //----------------------
1435 
1436  // Use average of the 4 FD approximations based on the
1437  // elements central pressure values
1438  // [Other options: Take average of the four
1439  // pressure derivatives]
1440 
1441  slope1 =
1442  octree_pt()->son_pt(LUF)->object_pt()
1443  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1444  octree_pt()->son_pt(LUB)->object_pt()
1445  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1446 
1447  slope2 =
1448  octree_pt()->son_pt(RUF)->object_pt()
1449  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1450  octree_pt()->son_pt(RUB)->object_pt()
1451  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1452 
1453  slope3 =
1454  octree_pt()->son_pt(LDF)->object_pt()
1455  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1456  octree_pt()->son_pt(LDB)->object_pt()
1457  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1458 
1459  slope4 =
1460  octree_pt()->son_pt(RDF)->object_pt()
1461  ->internal_data_pt(this->P_nst_internal_index)->value(0) -
1462  octree_pt()->son_pt(RDB)->object_pt()
1463  ->internal_data_pt(this->P_nst_internal_index)->value(0);
1464 
1465  // Use the average
1466  internal_data_pt(this->P_nst_internal_index)->
1467  set_value(3,0.25*(slope1+slope2+slope3+slope4));
1468 
1469 }
1470 
1471 
1472 //======================================================================
1473 /// 2D Further build for Crouzeix_Raviart interpolates the internal
1474 /// pressure dofs from father element: Make sure pressure values and
1475 /// dp/ds agree between fathers and sons at the midpoints of the son
1476 /// elements.
1477 //======================================================================
1478 template<>
1481 {
1482  //Call the generic further build
1484 
1485  using namespace QuadTreeNames;
1486 
1487  // What type of son am I? Ask my quadtree representation...
1488  int son_type=quadtree_pt()->son_type();
1489 
1490  // Pointer to my father (in element impersonation)
1491  RefineableElement* father_el_pt= quadtree_pt()->father_pt()->object_pt();
1492 
1493  Vector<double> s_father(2);
1494 
1495  // Son midpoint is located at the following coordinates in father element:
1496 
1497  // South west son
1498  if (son_type==SW)
1499  {
1500  s_father[0]=-0.5;
1501  s_father[1]=-0.5;
1502  }
1503  // South east son
1504  else if (son_type==SE)
1505  {
1506  s_father[0]= 0.5;
1507  s_father[1]=-0.5;
1508  }
1509  // North east son
1510  else if (son_type==NE)
1511  {
1512  s_father[0]=0.5;
1513  s_father[1]=0.5;
1514  }
1515 
1516  // North west son
1517  else if (son_type==NW)
1518  {
1519  s_father[0]=-0.5;
1520  s_father[1]= 0.5;
1521  }
1522 
1523  // Pressure value in father element
1525  cast_father_element_pt=
1527  (father_el_pt);
1528 
1529  double press=cast_father_element_pt->interpolated_p_nst(s_father);
1530 
1531  // Pressure value gets copied straight into internal dof:
1532  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1533 
1534  // The slopes get copied from father
1535  for(unsigned i=1;i<3;i++)
1536  {
1537  double half_father_slope = 0.5*cast_father_element_pt->
1538  internal_data_pt(this->P_nst_internal_index)->value(i);
1539  //Set the value in the son
1540  internal_data_pt(this->P_nst_internal_index)->
1541  set_value(i,half_father_slope);
1542  }
1543 }
1544 
1545 
1546 //=======================================================================
1547 /// 3D Further build for Crouzeix_Raviart interpolates the internal
1548 /// pressure dofs from father element: Make sure pressure values and
1549 /// dp/ds agree between fathers and sons at the midpoints of the son
1550 /// elements.
1551 //=======================================================================
1552 template<>
1555 {
1557 
1558  using namespace OcTreeNames;
1559 
1560  // What type of son am I? Ask my octree representation...
1561  int son_type=octree_pt()->son_type();
1562 
1563  // Pointer to my father (in element impersonation)
1564  RefineableQElement<3>* father_el_pt=
1565  dynamic_cast<RefineableQElement<3>*>(
1566  octree_pt()->father_pt()->object_pt());
1567 
1568  Vector<double> s_father(3);
1569 
1570  // Son midpoint is located at the following coordinates in father element:
1571  for(unsigned i=0;i<3;i++)
1572  {
1573  s_father[i]=0.5*OcTree::Direction_to_vector[son_type][i];
1574  }
1575 
1576  // Pressure value in father element
1578  cast_father_element_pt=
1580  (father_el_pt);
1581 
1582  double press=cast_father_element_pt->interpolated_p_nst(s_father);
1583 
1584  // Pressure value gets copied straight into internal dof:
1585  internal_data_pt(this->P_nst_internal_index)->set_value(0,press);
1586 
1587  // The slopes get copied from father
1588  for(unsigned i=1;i<4;i++)
1589  {
1590  double half_father_slope = 0.5*cast_father_element_pt
1591  ->internal_data_pt(this->P_nst_internal_index)->value(i);
1592  //Set the value
1593  internal_data_pt(this->P_nst_internal_index)
1594  ->set_value(i,half_father_slope);
1595  }
1596 }
1597 
1598 //=======================================================================
1599 /// 2D
1600 /// Derivatives of the shape functions and test functions w.r.t. to global
1601 /// (Eulerian) coordinates. Return Jacobian of mapping between
1602 /// local and global coordinates.
1603 //=======================================================================
1604 template<>
1607  DShape &dpsidx, Shape &test,
1608  DShape &dtestdx) const
1609 {
1610  //Call the geometrical shape functions and derivatives
1611  double J = this->dshape_eulerian(s,psi,dpsidx);
1612 
1613  //Loop over the test functions and derivatives and set them equal to the
1614  //shape functions
1615  for(unsigned i=0;i<nnode_1d()*nnode_1d();i++)
1616  {
1617  test[i] = psi[i];
1618  dtestdx(i,0) = dpsidx(i,0);
1619  dtestdx(i,1) = dpsidx(i,1);
1620  }
1621 
1622  //Return the jacobian
1623  return J;
1624 }
1625 
1626 //=======================================================================
1627 /// 2D
1628 /// Derivatives of the shape functions and test functions w.r.t. to global
1629 /// (Eulerian) coordinates. Return Jacobian of mapping between
1630 /// local and global coordinates.
1631 //=======================================================================
1632 template<>
1635  DShape &dpsidx, Shape &test,
1636  DShape &dtestdx) const
1637 {
1638  //Call the geometrical shape functions and derivatives
1639  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
1640 
1641  //Loop over the test functions and derivatives and set them equal to the
1642  //shape functions
1643  for(unsigned i=0;i<nnode_1d()*nnode_1d();i++)
1644  {
1645  test[i] = psi[i];
1646  dtestdx(i,0) = dpsidx(i,0);
1647  dtestdx(i,1) = dpsidx(i,1);
1648  }
1649 
1650  //Return the jacobian
1651  return J;
1652 }
1653 
1654 //=======================================================================
1655 /// 3D
1656 /// Derivatives of the shape functions and test functions w.r.t. to global
1657 /// (Eulerian) coordinates. Return Jacobian of mapping between
1658 /// local and global coordinates.
1659 //=======================================================================
1660 template<>
1663  DShape &dpsidx, Shape &test,
1664  DShape &dtestdx) const
1665 {
1666  //Call the geometrical shape functions and derivatives
1667  double J = this->dshape_eulerian(s,psi,dpsidx);
1668 
1669  //Loop over the test functions and derivatives and set them equal to the
1670  //shape functions
1671  for(unsigned i=0;i<nnode_1d()*nnode_1d()*nnode_1d();i++)
1672  {
1673  test[i] = psi[i];
1674  dtestdx(i,0) = dpsidx(i,0);
1675  dtestdx(i,1) = dpsidx(i,1);
1676  dtestdx(i,2) = dpsidx(i,2);
1677  }
1678 
1679  //Return the jacobian
1680  return J;
1681 }
1682 
1683 //=======================================================================
1684 /// 3D
1685 /// Derivatives of the shape functions and test functions w.r.t. to global
1686 /// (Eulerian) coordinates. Return Jacobian of mapping between
1687 /// local and global coordinates.
1688 //=======================================================================
1689 template<>
1692  DShape &dpsidx, Shape &test,
1693  DShape &dtestdx) const
1694 {
1695  //Call the geometrical shape functions and derivatives
1696  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
1697 
1698  //Loop over the test functions and derivatives and set them equal to the
1699  //shape functions
1700  for(unsigned i=0;i<nnode_1d()*nnode_1d()*nnode_1d();i++)
1701  {
1702  test[i] = psi[i];
1703  dtestdx(i,0) = dpsidx(i,0);
1704  dtestdx(i,1) = dpsidx(i,1);
1705  dtestdx(i,2) = dpsidx(i,2);
1706  }
1707 
1708  //Return the jacobian
1709  return J;
1710 }
1711 
1712 //=======================================================================
1713 /// 2D :
1714 /// Pressure shape functions
1715 //=======================================================================
1716 template<>
1718 pshape_nst(const Vector<double> &s, Shape &psi) const
1719 {
1720  unsigned npres = this->npres_nst();
1721  if (npres==1)
1722  {
1723  psi[0] = 1.0;
1724  }
1725  else
1726  {
1727  // Get number of pressure modes
1728  unsigned npres_1d = (int) std::sqrt((double)npres);
1729 
1730  //Local storage
1731  //Call the one-dimensional modal shape functions
1732  OneDimensionalModalShape psi1(npres_1d,s[0]);
1733  OneDimensionalModalShape psi2(npres_1d,s[1]);
1734 
1735  //Now let's loop over the nodal points in the element
1736  //s1 is the "x" coordinate, s2 the "y"
1737  for(unsigned i=0;i<npres_1d;i++)
1738  {
1739  for(unsigned j=0;j<npres_1d;j++)
1740  {
1741  //Multiply the two 1D functions together to get the 2D function
1742  psi[i*npres_1d + j] = psi2[i]*psi1[j];
1743  }
1744  }
1745  }
1746 }
1747 
1748 ///Define the pressure shape and test functions
1749 template<>
1751 pshape_nst(const Vector<double> &s, Shape &psi, Shape &test) const
1752 {
1753  //Call the pressure shape functions
1754  pshape_nst(s,psi);
1755 
1756  //Loop over the test functions and set them equal to the shape functions
1757  if (this->npres_nst()==1)
1758  {
1759  test[0] = psi[0];
1760  }
1761  else
1762  {
1763  for(unsigned i=0;i<this->npres_nst();i++) test[i] = psi[i];
1764  }
1765 }
1766 
1767 //=======================================================================
1768 /// 3D :
1769 /// Pressure shape functions
1770 //=======================================================================
1771 template<>
1773 pshape_nst(const Vector<double> &s, Shape &psi) const
1774 {
1775  unsigned npres = this->npres_nst();
1776  if (npres==1)
1777  {
1778  psi[0] = 1.0;
1779  }
1780  else
1781  {
1782  // Get number of pressure modes
1783  unsigned npres_1d = (int) std::sqrt((double)npres);
1784 
1785  //Local storage
1786  //Call the one-dimensional modal shape functions
1787  OneDimensionalModalShape psi1(npres_1d,s[0]);
1788  OneDimensionalModalShape psi2(npres_1d,s[1]);
1789  OneDimensionalModalShape psi3(npres_1d,s[2]);
1790 
1791  //Now let's loop over the nodal points in the element
1792  //s1 is the "x" coordinate, s2 the "y"
1793  for(unsigned i=0;i<npres_1d;i++)
1794  {
1795  for(unsigned j=0;j<npres_1d;j++)
1796  {
1797  for(unsigned k=0;k<npres_1d;k++)
1798  {
1799  //Multiply the two 1D functions together to get the 2D function
1800  psi[i*npres_1d*npres_1d + j*npres_1d + k] = psi3[i]*psi2[j]*psi1[k];
1801  }
1802  }
1803  }
1804  }
1805 }
1806 
1807 ///Define the pressure shape and test functions
1808 template<>
1810 pshape_nst(const Vector<double> &s, Shape &psi, Shape &test) const
1811 {
1812  //Call the pressure shape functions
1813  pshape_nst(s,psi);
1814 
1815  //Loop over the test functions and set them equal to the shape functions
1816  if (this->npres_nst()==1)
1817  {
1818  test[0] = psi[0];
1819  }
1820  else
1821  {
1822  for(unsigned i=0;i<this->npres_nst();i++) test[i] = psi[i];
1823  }
1824 }
1825 
1826 }
1827 
1828 #endif
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...
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.
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
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 pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
virtual unsigned u_index_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: (DIM velocities + 1 pressure)
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)
virtual int p_nodal_index_nst() const
Set the value at which the pressure is stored in the nodes.
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
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...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
static void unpin_all_pressure_dofs(const Vector< GeneralisedElement * > &element_pt)
Unpin all pressure dofs in elements listed in vector.
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons This must be specialised for each dime...
cstr elem_len * i
Definition: cfortran.h:607
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Pin unused nodal pressure dofs (empty by default, because by default pressure dofs are not associated...
NavierStokesSourceFctPt & source_fct_pt()
Access function for the source-function pointer.
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3148
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
void pin_elemental_redundant_nodal_pressure_dofs()
Pin all nodal pressure dofs that are not required.
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
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
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 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
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number) ...
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
e
Definition: cfortran.h:575
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 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
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...
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 ...
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes...
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}.
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...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
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...
A Rank 3 Tensor class.
Definition: matrices.h:1337
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...
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
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
void interpolated_u_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
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...
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
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
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...
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
NavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
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
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
Class that contains data for hanging nodes.
Definition: nodes.h:684
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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 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...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation:
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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...
virtual void resize(const unsigned &n_value)
Change (increase) the number of values that may be stored.
Definition: nodes.cc:961
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor. ...
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 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 This must be specialised for each dime...
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 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...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: DIM (velocities)
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 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 ...
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 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)...
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...
Refineable version of Crouzeix Raviart elements. Generic class definitions.
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
double *& re_invfr_pt()
Pointer to global inverse Froude number.
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number) ...
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
unsigned npres_nst() const
Return number of pressure values.