Taxisym_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 triangular/tetrahedaral Axisymmetric Navier Stokes elements
31 
32 #ifndef OOMPH_TAXISYM_NAVIER_STOKES_ELEMENTS_HEADER
33 #define OOMPH_TAXISYM_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 
41 //OOMPH-LIB headers
42 #include "../generic/Telements.h"
44 #include "../generic/error_estimator.h"
45 
46 namespace oomph
47 {
48 
49 //////////////////////////////////////////////////////////////////////////////
50 //////////////////////////////////////////////////////////////////////////////
51 // NOTE: TRI/TET CROZIER RAVIARTS REQUIRE BUBBLE FUNCTIONS! THEY'RE NOT
52 // STRAIGHTFORWARD GENERALISATIONS OF THE Q-EQUIVALENTS (WHICH ARE
53 // LBB UNSTABLE!)
54 //////////////////////////////////////////////////////////////////////////////
55 //////////////////////////////////////////////////////////////////////////////
56 
57 
58 //==========================================================================
59 ///AxisymmetricTCrouzeix_Raviart elements are Navier--Stokes elements with quadratic
60 ///interpolation for velocities and positions enriched by a single cubic
61 ///bubble function, but a discontinuous linear
62 ///pressure interpolation
63 //==========================================================================
65  public virtual TBubbleEnrichedElement<2,3>,
66  public virtual AxisymmetricNavierStokesEquations,
67  public virtual ElementWithZ2ErrorEstimator
68 {
69  protected:
70 
71 
72  /// Internal index that indicates at which internal datum the pressure is
73  /// stored
75 
76 
77  /// \short Velocity shape and test functions and their derivs
78  /// w.r.t. to global coords at local coordinate s (taken from geometry)
79  ///Return Jacobian of mapping between local and global coordinates.
81  Shape &psi,
82  DShape &dpsidx, Shape &test,
83  DShape &dtestdx) const;
84 
85  /// \short Velocity shape and test functions and their derivs
86  /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
87  ///Return Jacobian of mapping between local and global coordinates.
88  inline double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt,
89  Shape &psi,
90  DShape &dpsidx,
91  Shape &test,
92  DShape &dtestdx)
93  const;
94 
95  /// \short Shape/test functions and derivs w.r.t. to global coords at
96  /// integration point ipt; return Jacobian of mapping (J). Also compute
97  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
99  const unsigned &ipt,
100  Shape &psi,
101  DShape &dpsidx,
102  RankFourTensor<double> &d_dpsidx_dX,
103  Shape &test,
104  DShape &dtestdx,
105  RankFourTensor<double> &d_dtestdx_dX,
106  DenseMatrix<double> &djacobian_dX) const;
107 
108  /// \short Pressure shape and test functions and their derivs
109  /// w.r.t. to global coords at local coordinate s (taken from geometry)
110  /// Return Jacobian of mapping between local and global coordinates.
112  Shape &ppsi,
113  DShape &dppsidx,
114  Shape &ptest,
115  DShape &dptestdx) const;
116 
117 public:
118 
119  /// Pressure shape functions at local coordinate s
120  inline void pshape_axi_nst(const Vector<double> &s, Shape &psi) const;
121 
122  /// Pressure shape and test functions at local coordinte s
123  inline void pshape_axi_nst(const Vector<double> &s, Shape &psi, Shape &test)
124  const;
125 
126  /// Unpin all internal pressure dofs
128 
129  /// Return the local equation numbers for the pressure values.
130  inline int p_local_eqn(const unsigned &n) const
132 
133 public:
134 
135  /// Constructor, there are 3 internal values (for the pressure)
138  {
139  //Allocate and a single internal datum with 3 entries for the
140  //pressure
142  }
143 
144  /// Broken copy constructor
146  const
148  {
149  BrokenCopy::broken_copy("AxisymmetricTCrouzeixRaviartElement");
150  }
151 
152  /// Broken assignment operator
153 //Commented out broken assignment operator because this can lead to a conflict warning
154 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
155 //realise that two separate implementations of the broken function are the same and so,
156 //quite rightly, it shouts.
157  /*void operator=(const AxisymmetricTCrouzeixRaviartElement&)
158  {
159  BrokenCopy::broken_assign("AxisymmetricTCrouzeixRaviartElement");
160  }*/
161 
162 
163  /// \short Number of values (pinned or dofs) required at local node n.
164  inline virtual unsigned required_nvalue(const unsigned &n) const
165  {return 3;}
166 
167 
168  /// \short Return the pressure values at internal dof i_internal
169  /// (Discontinous pressure interpolation -- no need to cater for hanging
170  /// nodes).
171  double p_axi_nst(const unsigned &i) const
172  {return this->internal_data_pt(P_axi_nst_internal_index)->value(i);}
173 
174  /// Return number of pressure values
175  unsigned npres_axi_nst() const {return 3;}
176 
177  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
178  void fix_pressure(const unsigned &p_dof, const double &p_value)
179  {
181  this->internal_data_pt(P_axi_nst_internal_index)->set_value(p_dof, p_value);
182  }
183 
184  /// \short Build FaceElements that apply the Robin boundary condition
185  /// to the pressure advection diffusion problem required by
186  /// Fp preconditioner
187  //void build_fp_press_adv_diff_robin_bc_element(const unsigned&
188  // face_index)
189  // {
190  // this->Pressure_advection_diffusion_robin_element_pt.push_back(
191  // new FpPressureAdvDiffRobinBCElement<AxisymmetricTCrouzeixRaviartElement<DIM> >(
192  // this, face_index));
193  // }
194 
195  /// \short Add to the set paired_load_data
196  /// pairs of pointers to data objects and unsignedegers that
197  /// index the values in the data object that affect the load (traction),
198  /// as specified in the get_load() function.
199  void identify_load_data(std::set<std::pair<Data*,unsigned> >
200  &paired_load_data);
201 
202  /// \short Add to the set \c paired_pressure_data pairs
203  /// containing
204  /// - the pointer to a Data object
205  /// and
206  /// - the index of the value in that Data object
207  /// .
208  /// for all pressure values that affect the
209  /// load computed in the \c get_load(...) function.
211  std::set<std::pair<Data*,unsigned> > &paired_pressure_data);
212 
213  /// Redirect output to NavierStokesEquations output
214  void output(std::ostream &outfile)
216 
217  /// Redirect output to NavierStokesEquations output
218  void output(std::ostream &outfile, const unsigned &nplot)
220 
221  /// Redirect output to NavierStokesEquations output
222  void output(FILE* file_pt)
224 
225  /// Redirect output to NavierStokesEquations output
226  void output(FILE* file_pt, const unsigned &n_plot)
228 
229 
230  /// \short Order of recovery shape functions for Z2 error estimation:
231  /// Same order as unenriched shape functions.
232  unsigned nrecovery_order() {return 2;}
233 
234  /// \short Number of vertex nodes in the element
235  unsigned nvertex_node() const
236  {return 3;}
237 
238  /// \short Pointer to the j-th vertex node in the element
239  Node* vertex_node_pt(const unsigned& j) const
240  {return node_pt(j);}
241 
242  /// Number of 'flux' terms for Z2 error estimation
243  unsigned num_Z2_flux_terms()
244  {
245  // 3 diagonal strain rates, 3 off diagonal
246  return 6;
247  }
248 
249  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
250  /// in strain rate tensor.
252  {
253 #ifdef PARANOID
254  unsigned num_entries=6;
255  if (flux.size() < num_entries)
256  {
257  std::ostringstream error_message;
258  error_message << "The flux vector has the wrong number of entries, "
259  << flux.size() << ", whereas it should be at least "
260  << num_entries << std::endl;
261  throw OomphLibError(error_message.str(),
262  OOMPH_CURRENT_FUNCTION,
263  OOMPH_EXCEPTION_LOCATION);
264  }
265 #endif
266 
267  // Get strain rate matrix
268  DenseMatrix<double> strainrate(3);
269  this->strain_rate(s,strainrate);
270 
271  // Pack into flux Vector
272  unsigned icount=0;
273 
274  // Start with diagonal terms
275  for(unsigned i=0;i<3;i++)
276  {
277  flux[icount]=strainrate(i,i);
278  icount++;
279  }
280 
281  //Off diagonals row by row
282  for(unsigned i=0;i<3;i++)
283  {
284  for(unsigned j=i+1;j<3;j++)
285  {
286  flux[icount]=strainrate(i,j);
287  icount++;
288  }
289  }
290  }
291 
292 
293  /// \short The number of "DOF types" that degrees of freedom in this element
294  /// are sub-divided into: Velocity (3 components) and pressure.
295  unsigned ndof_types() const
296  {
297  return 4;
298  }
299 
300  /// \short Create a list of pairs for all unknowns in this element,
301  /// so that the first entry in each pair contains the global equation
302  /// number of the unknown, while the second one contains the number
303  /// of the "DOF type" that this unknown is associated with.
304  /// (Function can obviously only be called if the equation numbering
305  /// scheme has been set up.)
307  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
308  {
309  // number of nodes
310  unsigned n_node = this->nnode();
311 
312  // number of pressure values
313  unsigned n_press = this->npres_axi_nst();
314 
315  // temporary pair (used to store dof lookup prior to being added to list)
316  std::pair<unsigned,unsigned> dof_lookup;
317 
318  // pressure dof number
319  unsigned pressure_dof_number = 3;
320 
321  // loop over the pressure values
322  for (unsigned n = 0; n < n_press; n++)
323  {
324  // determine local eqn number
325  int local_eqn_number = this->p_local_eqn(n);
326 
327  // ignore pinned values - far away degrees of freedom resulting
328  // from hanging nodes can be ignored since these are be dealt
329  // with by the element containing their master nodes
330  if (local_eqn_number >= 0)
331  {
332  // store dof lookup in temporary pair: First entry in pair
333  // is global equation number; second entry is dof type
334  dof_lookup.first = this->eqn_number(local_eqn_number);
335  dof_lookup.second = pressure_dof_number;
336 
337  // add to list
338  dof_lookup_list.push_front(dof_lookup);
339  }
340  }
341 
342  // loop over the nodes
343  for (unsigned n = 0; n < n_node; n++)
344  {
345  // find the number of values at this node
346  unsigned nv = this->node_pt(n)->nvalue();
347 
348  //loop over these values
349  for (unsigned v = 0; v < nv; v++)
350  {
351  // determine local eqn number
352  int local_eqn_number = this->nodal_local_eqn(n, v);
353 
354  // ignore pinned values
355  if (local_eqn_number >= 0)
356  {
357  // store dof lookup in temporary pair: First entry in pair
358  // is global equation number; second entry is dof type
359  dof_lookup.first = this->eqn_number(local_eqn_number);
360  dof_lookup.second = v;
361 
362  // add to list
363  dof_lookup_list.push_front(dof_lookup);
364 
365  }
366  }
367  }
368  }
369 
370 
371 };
372 
373 //Inline functions
374 
375 //=======================================================================
376 /// Derivatives of the shape functions and test functions w.r.t. to global
377 /// (Eulerian) coordinates. Return Jacobian of mapping between
378 /// local and global coordinates.
379 //=======================================================================
382  const Vector<double> &s, Shape &psi,
383  DShape &dpsidx, Shape &test,
384  DShape &dtestdx) const
385  {
386  //Call the geometrical shape functions and derivatives
387  double J = this->dshape_eulerian(s,psi,dpsidx);
388  //The test functions are equal to the shape functions
389  test = psi;
390  dtestdx = dpsidx;
391  //Return the jacobian
392  return J;
393  }
394 
395 
396 //=======================================================================
397 /// Derivatives of the shape functions and test functions w.r.t. to global
398 /// (Eulerian) coordinates. Return Jacobian of mapping between
399 /// local and global coordinates.
400 //=======================================================================
403  const unsigned &ipt, Shape &psi,
404  DShape &dpsidx, Shape &test,
405  DShape &dtestdx) const
406  {
407  //Call the geometrical shape functions and derivatives
408  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
409  //The test functions are the shape functions
410  test = psi;
411  dtestdx = dpsidx;
412  //Return the jacobian
413  return J;
414  }
415 
416 
417 //=======================================================================
418 /// Define the shape functions (psi) and test functions (test) and
419 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
420 /// and return Jacobian of mapping (J). Additionally compute the
421 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
422 ///
423 /// Galerkin: Test functions = shape functions
424 //=======================================================================
427  const unsigned &ipt, Shape &psi, DShape &dpsidx,
428  RankFourTensor<double> &d_dpsidx_dX,
429  Shape &test, DShape &dtestdx,
430  RankFourTensor<double> &d_dtestdx_dX,
431  DenseMatrix<double> &djacobian_dX) const
432  {
433  // Call the geometrical shape functions and derivatives
434  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
435  djacobian_dX,d_dpsidx_dX);
436 
437  // Set the test functions equal to the shape functions
438  test = psi;
439  dtestdx = dpsidx;
440  d_dtestdx_dX = d_dpsidx_dX;
441 
442  // Return the jacobian
443  return J;
444  }
445 
446 
447 //=======================================================================
448 /// Pressure shape functions
449 //=======================================================================
451  const Vector<double> &s,
452  Shape &psi) const
453  {
454  psi[0] = 1.0;
455  psi[1] = s[0];
456  psi[2] = s[1];
457  }
458 
459 //=======================================================================
460 /// Pressure shape and test functions
461 //=======================================================================
463  const Vector<double> &s,
464  Shape &psi,
465  Shape &test) const
466  {
467  //Call the pressure shape functions
468  this->pshape_axi_nst(s,psi);
469  //The test functions are the shape functions
470  test = psi;
471  }
472 
473 //==========================================================================
474 /// 2D :
475 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
476 /// Return Jacobian of mapping between local and global coordinates.
477 //==========================================================================
480  const Vector<double> &s,
481  Shape &ppsi,
482  DShape &dppsidx,
483  Shape &ptest,
484  DShape &dptestdx) const
485  {
486 
487  // Initalise with shape fcts and derivs. w.r.t. to local coordinates
488  ppsi[0] = 1.0;
489  ppsi[1] = s[0];
490  ppsi[2] = s[1];
491 
492  dppsidx(0,0) = 0.0;
493  dppsidx(1,0) = 1.0;
494  dppsidx(2,0) = 0.0;
495 
496  dppsidx(0,1) = 0.0;
497  dppsidx(1,1) = 0.0;
498  dppsidx(2,1) = 1.0;
499 
500 
501  //Get the values of the shape functions and their local derivatives
502  Shape psi(7);
503  DShape dpsi(7,2);
504  dshape_local(s,psi,dpsi);
505 
506  //Allocate memory for the inverse 2x2 jacobian
507  DenseMatrix<double> inverse_jacobian(2);
508 
509  //Now calculate the inverse jacobian
510  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
511 
512  //Now set the values of the derivatives to be derivs w.r.t. to the
513  // Eulerian coordinates
514  transform_derivatives(inverse_jacobian,dppsidx);
515 
516  //The test functions are equal to the shape functions
517  ptest = ppsi;
518  dptestdx = dppsidx;
519 
520  //Return the determinant of the jacobian
521  return det;
522  }
523 
524 
525 //=======================================================================
526 /// Face geometry of the 2D Crouzeix_Raviart elements
527 //=======================================================================
528 template<>
530  public virtual TElement<1,3>
531  {
532  public:
533  FaceGeometry() : TElement<1,3>() {}
534  };
535 
536 
537 //=======================================================================
538 /// Face geometry of the FaceGeometry of the 2D CrouzeixRaviart elements
539 //=======================================================================
540 template<>
542 public virtual PointElement
543 {
544  public:
546 };
547 
548 
549 ////////////////////////////////////////////////////////////////////////////
550 ////////////////////////////////////////////////////////////////////////////
551 ////////////////////////////////////////////////////////////////////////////
552 
553 
554 
555 //=======================================================================
556 /// Taylor--Hood elements are Navier--Stokes elements
557 /// with quadratic interpolation for velocities and positions and
558 /// continous linear pressure interpolation
559 //=======================================================================
561 public virtual TElement<2,3>,
562 public virtual AxisymmetricNavierStokesEquations,
563 public virtual ElementWithZ2ErrorEstimator
564 
565 {
566  private:
567 
568  /// Static array of ints to hold number of variables at node
569  static const unsigned Initial_Nvalue[];
570 
571  protected:
572 
573  /// \short Static array of ints to hold conversion from pressure
574  /// node numbers to actual node numbers
575  static const unsigned Pconv[];
576 
577  /// \short Velocity shape and test functions and their derivs
578  /// w.r.t. to global coords at local coordinate s (taken from geometry)
579  /// Return Jacobian of mapping between local and global coordinates.
580  inline double dshape_and_dtest_eulerian_axi_nst(const Vector<double> &s,
581  Shape &psi,
582  DShape &dpsidx, Shape &test,
583  DShape &dtestdx) const;
584 
585  /// \short Velocity shape and test functions and their derivs
586  /// w.r.t. to global coords at local coordinate s (taken from geometry)
587  /// Return Jacobian of mapping between local and global coordinates.
588  inline double dshape_and_dtest_eulerian_at_knot_axi_nst(const unsigned &ipt,
589  Shape &psi,
590  DShape &dpsidx,
591  Shape &test,
592  DShape &dtestdx)
593  const;
594 
595  /// \short Shape/test functions and derivs w.r.t. to global coords at
596  /// integration point ipt; return Jacobian of mapping (J). Also compute
597  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
599  const unsigned &ipt,
600  Shape &psi,
601  DShape &dpsidx,
602  RankFourTensor<double> &d_dpsidx_dX,
603  Shape &test,
604  DShape &dtestdx,
605  RankFourTensor<double> &d_dtestdx_dX,
606  DenseMatrix<double> &djacobian_dX) const;
607 
608  /// \short Compute the pressure shape and test functions and derivatives
609  /// w.r.t. global coords at local coordinate s.
610  /// Return Jacobian of mapping between local and global coordinates.
612  Shape &ppsi,
613  DShape &dppsidx,
614  Shape &ptest,
615  DShape &dptestdx) const;
616 
617  /// Unpin all pressure dofs
619 
620  /// Pin all nodal pressure dofs
622 
623  /// Unpin the proper nodal pressure dofs
625 
626 
627 public:
628 
629  /// Constructor, no internal data points
632 
633 
634  /// Broken copy constructor
636  {
637  BrokenCopy::broken_copy("AxisymmetricTTaylorHoodElement");
638  }
639 
640  /// Broken assignment operator
641  /*void operator=(const AxisymmetricTTaylorHoodElement&)
642  {
643  BrokenCopy::broken_assign("AxisymmetricTTaylorHoodElement");
644  }*/
645 
646  /// \short Number of values (pinned or dofs) required at node n. Can
647  /// be overwritten for hanging node version
648  inline virtual unsigned required_nvalue(const unsigned &n) const
649  {return Initial_Nvalue[n];}
650 
651  /// Test whether the pressure dof p_dof hanging or not?
652  //bool pressure_dof_is_hanging(const unsigned& p_dof)
653  // {return this->node_pt(Pconv[p_dof])->is_hanging(DIM);}
654 
655 
656  /// Pressure shape functions at local coordinate s
657  inline void pshape_axi_nst(const Vector<double> &s, Shape &psi) const;
658 
659  /// Pressure shape and test functions at local coordinte s
660  inline void pshape_axi_nst(const Vector<double> &s, Shape &psi,
661  Shape &test) const;
662 
663  /// \short Which nodal value represents the pressure?
664  unsigned p_index_axi_nst() {return 3;}
665 
666  /// \short Pointer to n_p-th pressure node
667  //Node* pressure_node_pt(const unsigned &n_p)
668  //{return this->Node_pt[Pconv[n_p]];}
669 
670  /// Return the local equation numbers for the pressure values.
671  inline int p_local_eqn(const unsigned &n) const
672  {return this->nodal_local_eqn(Pconv[n],3);}
673 
674  /// \short Access function for the pressure values at local pressure
675  /// node n_p (const version)
676  double p_axi_nst(const unsigned &n_p) const
677  {return this->nodal_value(Pconv[n_p],3);}
678 
679  /// \short Set the value at which the pressure is stored in the nodes
680  int p_nodal_index_axi_nst() const {return static_cast<int>(3);}
681 
682  /// Return number of pressure values
683  unsigned npres_axi_nst() const {return 3;}
684 
685  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
686  void fix_pressure(const unsigned &p_dof, const double &p_value)
687  {
688  this->node_pt(Pconv[p_dof])->pin(3);
689  this->node_pt(Pconv[p_dof])->set_value(3,p_value);
690  }
691 
692 
693  /// \short Build FaceElements that apply the Robin boundary condition
694  /// to the pressure advection diffusion problem required by
695  /// Fp preconditioner
696 // void build_fp_press_adv_diff_robin_bc_element(const unsigned&
697 // face_index)
698 // {
699 // this->Pressure_advection_diffusion_robin_element_pt.push_back(
700 // new FpPressureAdvDiffRobinBCElement<AxisymmetricTTaylorHoodElement<DIM> >(
701 // this, face_index));
702 // }
703 
704  /// \short Add to the set \c paired_load_data pairs containing
705  /// - the pointer to a Data object
706  /// and
707  /// - the index of the value in that Data object
708  /// .
709  /// for all values (pressures, velocities) that affect the
710  /// load computed in the \c get_load(...) function.
711  void identify_load_data(
712  std::set<std::pair<Data*,unsigned> > &paired_load_data);
713 
714  /// \short Add to the set \c paired_pressure_data pairs
715  /// containing
716  /// - the pointer to a Data object
717  /// and
718  /// - the index of the value in that Data object
719  /// .
720  /// for all pressure values that affect the
721  /// load computed in the \c get_load(...) function.
723  std::set<std::pair<Data*,unsigned> > &paired_pressure_data);
724 
725  /// Redirect output to NavierStokesEquations output
726  void output(std::ostream &outfile)
728 
729  /// Redirect output to NavierStokesEquations output
730  void output(std::ostream &outfile, const unsigned &nplot)
732 
733  /// Redirect output to NavierStokesEquations output
734  void output(FILE* file_pt)
736 
737  /// Redirect output to NavierStokesEquations output
738  void output(FILE* file_pt, const unsigned &n_plot)
740 
741  /// \short Order of recovery shape functions for Z2 error estimation:
742  /// Same order as shape functions.
743  unsigned nrecovery_order() {return 2;}
744 
745  /// \short Number of vertex nodes in the element
746  unsigned nvertex_node() const
747  {return 3;}
748 
749  /// \short Pointer to the j-th vertex node in the element
750  Node* vertex_node_pt(const unsigned& j) const
751  {return node_pt(j);}
752 
753 
754  /// Number of 'flux' terms for Z2 error estimation
755  unsigned num_Z2_flux_terms()
756  {
757  // 3 diagonal strain rates, 3 off diagonal rates
758  return 6;
759  }
760 
761  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
762  /// in strain rate tensor.
764  {
765 #ifdef PARANOID
766  unsigned num_entries=6;
767  if (flux.size() < num_entries)
768  {
769  std::ostringstream error_message;
770  error_message << "The flux vector has the wrong number of entries, "
771  << flux.size() << ", whereas it should be at least "
772  << num_entries << std::endl;
773  throw OomphLibError(error_message.str(),
774  OOMPH_CURRENT_FUNCTION,
775  OOMPH_EXCEPTION_LOCATION);
776  }
777 #endif
778 
779  // Get strain rate matrix
780  DenseMatrix<double> strainrate(3);
781  this->strain_rate(s,strainrate);
782 
783  // Pack into flux Vector
784  unsigned icount=0;
785 
786  // Start with diagonal terms
787  for(unsigned i=0;i<3;i++)
788  {
789  flux[icount]=strainrate(i,i);
790  icount++;
791  }
792 
793  //Off diagonals row by row
794  for(unsigned i=0;i<3;i++)
795  {
796  for(unsigned j=i+1;j<3;j++)
797  {
798  flux[icount]=strainrate(i,j);
799  icount++;
800  }
801  }
802  }
803 
804  /// \short The number of "DOF types" that degrees of freedom in this element
805  /// are sub-divided into: Velocities (3 components) and pressure.
806  unsigned ndof_types() const
807  {
808  return 4;
809  }
810 
811  /// \short Create a list of pairs for all unknowns in this element,
812  /// so that the first entry in each pair contains the global equation
813  /// number of the unknown, while the second one contains the number
814  /// of the "DOF type" that this unknown is associated with.
815  /// (Function can obviously only be called if the equation numbering
816  /// scheme has been set up.)
818  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
819  {
820  // number of nodes
821  unsigned n_node = this->nnode();
822 
823  // temporary pair (used to store dof lookup prior to being added to list)
824  std::pair<unsigned,unsigned> dof_lookup;
825 
826  // loop over the nodes
827  for (unsigned n = 0; n < n_node; n++)
828  {
829  // find the number of Navier Stokes values at this node
830  unsigned nv = this->required_nvalue(n);
831 
832  //loop over these values
833  for (unsigned v = 0; v < nv; v++)
834  {
835  // determine local eqn number
836  int local_eqn_number = this->nodal_local_eqn(n, v);
837 
838  // ignore pinned values - far away degrees of freedom resulting
839  // from hanging nodes can be ignored since these are be dealt
840  // with by the element containing their master nodes
841  if (local_eqn_number >= 0)
842  {
843  // store dof lookup in temporary pair: Global equation number
844  // is the first entry in pair
845  dof_lookup.first = this->eqn_number(local_eqn_number);
846 
847  // set dof numbers: Dof number is the second entry in pair
848  dof_lookup.second = v;
849 
850  // add to list
851  dof_lookup_list.push_front(dof_lookup);
852  }
853  }
854  }
855  }
856 };
857 
858 
859 
860 
861 //Inline functions
862 
863 //==========================================================================
864 /// Derivatives of the shape functions and test functions w.r.t to
865 /// global (Eulerian) coordinates. Return Jacobian of mapping between
866 /// local and global coordinates.
867 //==========================================================================
870  const Vector<double> &s,
871  Shape &psi,
872  DShape &dpsidx, Shape &test,
873  DShape &dtestdx) const
874 {
875  //Call the geometrical shape functions and derivatives
876  double J = this->dshape_eulerian(s,psi,dpsidx);
877  //Test functions are the shape functions
878  test = psi;
879  dtestdx = dpsidx;
880  //Return the jacobian
881  return J;
882 }
883 
884 
885 //==========================================================================
886 /// Derivatives of the shape functions and test functions w.r.t to
887 /// global (Eulerian) coordinates. Return Jacobian of mapping between
888 /// local and global coordinates.
889 //==========================================================================
892  const unsigned &ipt,Shape &psi, DShape &dpsidx, Shape &test,
893  DShape &dtestdx) const
894 {
895  //Call the geometrical shape functions and derivatives
896  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
897  //Test functions are the shape functions
898  test = psi;
899  dtestdx = dpsidx;
900  //Return the jacobian
901  return J;
902 }
903 
904 //==========================================================================
905 /// Define the shape functions (psi) and test functions (test) and
906 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
907 /// and return Jacobian of mapping (J). Additionally compute the
908 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
909 ///
910 /// Galerkin: Test functions = shape functions
911 //==========================================================================
914  const unsigned &ipt, Shape &psi, DShape &dpsidx,
915  RankFourTensor<double> &d_dpsidx_dX,
916  Shape &test, DShape &dtestdx,
917  RankFourTensor<double> &d_dtestdx_dX,
918  DenseMatrix<double> &djacobian_dX) const
919  {
920  // Call the geometrical shape functions and derivatives
921  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
922  djacobian_dX,d_dpsidx_dX);
923 
924  // Set the test functions equal to the shape functions
925  test = psi;
926  dtestdx = dpsidx;
927  d_dtestdx_dX = d_dpsidx_dX;
928 
929  // Return the jacobian
930  return J;
931  }
932 
933 //==========================================================================
934 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
935 /// Return Jacobian of mapping between local and global coordinates.
936 //==========================================================================
939  const Vector<double> &s,
940  Shape &ppsi,
941  DShape &dppsidx,
942  Shape &ptest,
943  DShape &dptestdx) const
944  {
945 
946  ppsi[0] = s[0];
947  ppsi[1] = s[1];
948  ppsi[2] = 1.0-s[0]-s[1];
949 
950  dppsidx(0,0)=1.0;
951  dppsidx(0,1)=0.0;
952 
953  dppsidx(1,0)=0.0;
954  dppsidx(1,1)=1.0;
955 
956  dppsidx(2,0)=-1.0;
957  dppsidx(2,1)=-1.0;
958 
959  // Allocate memory for the inverse 2x2 jacobian
960  DenseMatrix<double> inverse_jacobian(2);
961 
962 
963  //Get the values of the shape functions and their local derivatives
964  Shape psi(6);
965  DShape dpsi(6,2);
966  dshape_local(s,psi,dpsi);
967 
968  // Now calculate the inverse jacobian
969  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
970 
971  // Now set the values of the derivatives to be derivs w.r.t. to the
972  // Eulerian coordinates
973  transform_derivatives(inverse_jacobian,dppsidx);
974 
975  //Test functions are shape functions
976  ptest = ppsi;
977  dptestdx = dppsidx;
978 
979  // Return the determinant of the jacobian
980  return det;
981 
982  }
983 
984 
985 //==========================================================================
986 /// Pressure shape functions
987 //==========================================================================
990 const
991 {
992  psi[0] = s[0];
993  psi[1] = s[1];
994  psi[2] = 1.0-s[0]-s[1];
995 }
996 
997 //==========================================================================
998 /// Pressure shape and test functions
999 //==========================================================================
1002  Shape &psi,
1003  Shape &test) const
1004 {
1005  //Call the pressure shape functions
1006  this->pshape_axi_nst(s,psi);
1007  //Test functions are shape functions
1008  test = psi;
1009 }
1010 
1011 
1012 //=======================================================================
1013 /// Face geometry of the Axisymmetric Taylor_Hood elements
1014 //=======================================================================
1015 template<>
1017 public virtual TElement<1,3>
1018 {
1019  public:
1020 
1021  /// Constructor: Call constructor of base
1022  FaceGeometry() : TElement<1,3>() {}
1023 };
1024 
1025 
1026 
1027 //=======================================================================
1028 /// Face geometry of the FaceGeometry of the Axisymmetric TaylorHood elements
1029 //=======================================================================
1030 template<>
1032 public virtual PointElement
1033 {
1034  public:
1036 };
1037 
1038 
1039 }
1040 
1041 #endif
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers...
void pshape_axi_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
int internal_local_eqn(const unsigned &i, const unsigned &j) const
Return the local equation number corresponding to the j-th value stored at the i-th internal data...
Definition: elements.h:268
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
virtual unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
double dpshape_and_dptest_eulerian_axi_nst(const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const
Pressure shape and test functions and their derivs w.r.t. to global coords at local coordinate s (tak...
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
Definition: elements.cc:3225
cstr elem_len * i
Definition: cfortran.h:607
AxisymmetricTTaylorHoodElement(const AxisymmetricTTaylorHoodElement &dummy)
Broken copy constructor.
void pin_all_nodal_pressure_dofs()
Pin all nodal pressure dofs.
int p_nodal_index_axi_nst() const
Set the value at which the pressure is stored in the nodes.
int local_eqn_number(const unsigned long &ieqn_global) const
Return the local equation number corresponding to the ieqn_global-th global equation number...
Definition: elements.h:731
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
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
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
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as unenriched shape functions...
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
virtual void transform_derivatives(const DenseMatrix< double > &inverse_jacobian, DShape &dbasis) const
Convert derivative w.r.t.local coordinates to derivatives w.r.t the coordinates used to assemble the ...
Definition: elements.cc:2762
void unpin_all_nodal_pressure_dofs()
Unpin all pressure dofs.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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.
virtual unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
A Rank 4 Tensor class.
Definition: matrices.h:1625
void unpin_proper_nodal_pressure_dofs()
Unpin the proper nodal pressure dofs.
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
AxisymmetricTCrouzeixRaviartElement()
Constructor, there are 3 internal values (for the pressure)
void pshape_axi_nst(const Vector< double > &s, Shape &psi) const
Test whether the pressure dof p_dof hanging or not?
unsigned npres_axi_nst() const
Return number of pressure values.
virtual double local_to_eulerian_mapping(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates, given the derivatives of the shape function...
Definition: elements.h:1461
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
unsigned p_index_axi_nst()
Which nodal value represents the pressure?
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
AxisymmetricTCrouzeixRaviartElement(const AxisymmetricTCrouzeixRaviartElement &dummy)
Broken copy constructor.
int p_local_eqn(const unsigned &n) const
Return the local equation numbers for the pressure values.
double dshape_and_dtest_eulerian_at_knot_axi_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 get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor. ...
static char t char * s
Definition: cfortran.h:572
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
double dshape_and_dtest_eulerian_at_knot_axi_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 local coordinate s (tak...
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
AxisymmetricTTaylorHoodElement()
Constructor, no internal data points.
void unpin_all_internal_pressure_dofs()
Unpin all internal pressure dofs.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
int p_local_eqn(const unsigned &n) const
Pointer to n_p-th pressure node.
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node...
Definition: elements.h:1383
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocities (3...
virtual double dpshape_and_dptest_eulerian_axi_nst(const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const
Compute the pressure shape and test functions and derivatives w.r.t. global coords at local coordinat...
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
unsigned npres_axi_nst() const
Return number of pressure values.
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
double dshape_and_dtest_eulerian_axi_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...
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocity (3 c...
double p_axi_nst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version) ...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
unsigned nvertex_node() const
Number of vertex nodes in the element.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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 add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
Definition: elements.cc:66
double p_axi_nst(const unsigned &i) const
Return the pressure values at internal dof i_internal (Discontinous pressure interpolation – no need ...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
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 output(std::ostream &outfile, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
void output(std::ostream &outfile, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
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.
double dshape_and_dtest_eulerian_axi_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 void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
Definition: elements.h:1914
unsigned nvertex_node() const
Number of vertex nodes in the element.
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3252