Tnavier_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 Navier Stokes elements
31 
32 #ifndef OOMPH_TNAVIER_STOKES_ELEMENTS_HEADER
33 #define OOMPH_TNAVIER_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"
43 #include "navier_stokes_elements.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 ///TCrouzeix_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 //==========================================================================
64 template <unsigned DIM>
65 class TCrouzeixRaviartElement : public virtual TBubbleEnrichedElement<DIM,3>,
66  public virtual NavierStokesEquations<DIM>,
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.
80  inline double dshape_and_dtest_eulerian_nst(const Vector<double> &s,
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_nst(const unsigned &ipt,
89  Shape &psi,
90  DShape &dpsidx, Shape &test,
91  DShape &dtestdx) const;
92 
93  /// \short Shape/test functions and derivs w.r.t. to global coords at
94  /// integration point ipt; return Jacobian of mapping (J). Also compute
95  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
97  const unsigned &ipt,
98  Shape &psi,
99  DShape &dpsidx,
100  RankFourTensor<double> &d_dpsidx_dX,
101  Shape &test,
102  DShape &dtestdx,
103  RankFourTensor<double> &d_dtestdx_dX,
104  DenseMatrix<double> &djacobian_dX) const;
105 
106  /// \short Pressure shape and test functions and their derivs
107  /// w.r.t. to global coords at local coordinate s (taken from geometry)
108  /// Return Jacobian of mapping between local and global coordinates.
109  inline double dpshape_and_dptest_eulerian_nst(const Vector<double> &s,
110  Shape &ppsi,
111  DShape &dppsidx,
112  Shape &ptest,
113  DShape &dptestdx) const;
114 
115 public:
116 
117  /// Pressure shape functions at local coordinate s
118  inline void pshape_nst(const Vector<double> &s, Shape &psi) const;
119 
120  /// Pressure shape and test functions at local coordinte s
121  inline void pshape_nst(const Vector<double> &s, Shape &psi, Shape &test)
122  const;
123 
124  /// Unpin all internal pressure dofs
126 
127  /// Return the local equation numbers for the pressure values.
128  inline int p_local_eqn(const unsigned &n) const
129  {return this->internal_local_eqn(P_nst_internal_index,n);}
130 
131 public:
132 
133  /// Constructor, there are DIM+1 internal values (for the pressure)
135  NavierStokesEquations<DIM>()
136  {
137  //Allocate and a single internal datum with DIM+1 entries for the
138  //pressure
139  P_nst_internal_index = this->add_internal_data(new Data(DIM+1));
140  }
141 
142  /// Broken copy constructor
144  {
145  BrokenCopy::broken_copy("TCrouzeixRaviartElement");
146  }
147 
148  /// Broken assignment operator
149 //Commented out broken assignment operator because this can lead to a conflict warning
150 //when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
151 //realise that two separate implementations of the broken function are the same and so,
152 //quite rightly, it shouts.
153  /*void operator=(const TCrouzeixRaviartElement<DIM>&)
154  {
155  BrokenCopy::broken_assign("TCrouzeixRaviartElement");
156  }*/
157 
158 
159  /// \short Number of values (pinned or dofs) required at local node n.
160  inline virtual unsigned required_nvalue(const unsigned &n) const
161  {return DIM;}
162 
163 
164  /// \short Return the pressure values at internal dof i_internal
165  /// (Discontinous pressure interpolation -- no need to cater for hanging
166  /// nodes).
167  double p_nst(const unsigned &i) const
168  {return this->internal_data_pt(P_nst_internal_index)->value(i);}
169 
170  /// \short Return the pressure values at internal dof i_internal
171  /// (Discontinous pressure interpolation -- no need to cater for hanging
172  /// nodes).
173  double p_nst(const unsigned &t, const unsigned &i) const
174  {return this->internal_data_pt(P_nst_internal_index)->value(t,i);}
175 
176  /// Return number of pressure values
177  unsigned npres_nst() const {return DIM+1;}
178 
179  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
180  void fix_pressure(const unsigned &p_dof, const double &p_value)
181  {
183  this->internal_data_pt(P_nst_internal_index)->set_value(p_dof, p_value);
184  }
185 
186  /// \short Build FaceElements that apply the Robin boundary condition
187  /// to the pressure advection diffusion problem required by
188  /// Fp preconditioner
190  face_index)
191  {
194  this, face_index));
195  }
196 
197  /// \short Add to the set paired_load_data
198  /// pairs of pointers to data objects and unsignedegers that
199  /// index the values in the data object that affect the load (traction),
200  /// as specified in the get_load() function.
201  void identify_load_data(std::set<std::pair<Data*,unsigned> >
202  &paired_load_data);
203 
204  /// \short Add to the set \c paired_pressure_data pairs
205  /// containing
206  /// - the pointer to a Data object
207  /// and
208  /// - the index of the value in that Data object
209  /// .
210  /// for all pressure values that affect the
211  /// load computed in the \c get_load(...) function.
213  std::set<std::pair<Data*,unsigned> > &paired_pressure_data);
214 
215  /// Redirect output to NavierStokesEquations output
216  void output(std::ostream &outfile)
218 
219  /// Redirect output to NavierStokesEquations output
220  void output(std::ostream &outfile, const unsigned &nplot)
221  {NavierStokesEquations<DIM>::output(outfile,nplot);}
222 
223  /// Redirect output to NavierStokesEquations output
224  void output(FILE* file_pt) {NavierStokesEquations<DIM>::output(file_pt);}
225 
226  /// Redirect output to NavierStokesEquations output
227  void output(FILE* file_pt, const unsigned &n_plot)
228  {NavierStokesEquations<DIM>::output(file_pt,n_plot);}
229 
230 
231  /// \short Order of recovery shape functions for Z2 error estimation:
232  /// Same order as unenriched shape functions.
233  unsigned nrecovery_order() {return 2;}
234 
235  /// \short Number of vertex nodes in the element
236  unsigned nvertex_node() const
237  {return DIM+1;}
238 
239  /// \short Pointer to the j-th vertex node in the element
240  Node* vertex_node_pt(const unsigned& j) const
241  {return node_pt(j);}
242 
243  /// Number of 'flux' terms for Z2 error estimation
244  unsigned num_Z2_flux_terms()
245  {
246  // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
247  return DIM + (DIM*(DIM-1))/2;
248  }
249 
250  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
251  /// in strain rate tensor.
253  {
254 #ifdef PARANOID
255  unsigned num_entries=DIM+(DIM*(DIM-1))/2;
256  if (flux.size() < num_entries)
257  {
258  std::ostringstream error_message;
259  error_message << "The flux vector has the wrong number of entries, "
260  << flux.size() << ", whereas it should be at least "
261  << num_entries << std::endl;
262  throw OomphLibError(error_message.str(),
263  OOMPH_CURRENT_FUNCTION,
264  OOMPH_EXCEPTION_LOCATION);
265  }
266 #endif
267 
268  // Get strain rate matrix
269  DenseMatrix<double> strainrate(DIM);
270  this->strain_rate(s,strainrate);
271 
272  // Pack into flux Vector
273  unsigned icount=0;
274 
275  // Start with diagonal terms
276  for(unsigned i=0;i<DIM;i++)
277  {
278  flux[icount]=strainrate(i,i);
279  icount++;
280  }
281 
282  //Off diagonals row by row
283  for(unsigned i=0;i<DIM;i++)
284  {
285  for(unsigned j=i+1;j<DIM;j++)
286  {
287  flux[icount]=strainrate(i,j);
288  icount++;
289  }
290  }
291  }
292 
293 
294 
295  /// \short Full output function:
296  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
297  /// in tecplot format. Default number of plot points
298  void full_output(std::ostream &outfile)
300 
301  /// \short Full output function:
302  /// x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation
303  /// in tecplot format. nplot points in each coordinate direction
304  void full_output(std::ostream &outfile, const unsigned &nplot)
306 
307 /// \short The number of "DOF types" that degrees of freedom in this element
308  /// are sub-divided into: Velocity and pressure.
309  unsigned ndof_types() const
310  {
311  return DIM+1;
312  }
313 
314  /// \short Create a list of pairs for all unknowns in this element,
315  /// so that the first entry in each pair contains the global equation
316  /// number of the unknown, while the second one contains the number
317  /// of the "DOF types" that this unknown is associated with.
318  /// (Function can obviously only be called if the equation numbering
319  /// scheme has been set up.) Velocity=0; Pressure=1
321  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const;
322 
323 
324 };
325 
326 //Inline functions
327 
328 //=======================================================================
329 /// Derivatives of the shape functions and test functions w.r.t. to global
330 /// (Eulerian) coordinates. Return Jacobian of mapping between
331 /// local and global coordinates.
332 //=======================================================================
333 template<unsigned DIM>
335  const Vector<double> &s, Shape &psi,
336  DShape &dpsidx, Shape &test,
337  DShape &dtestdx) const
338 {
339  //Call the geometrical shape functions and derivatives
340  double J = this->dshape_eulerian(s,psi,dpsidx);
341  //The test functions are equal to the shape functions
342  test = psi;
343  dtestdx = dpsidx;
344  //Return the jacobian
345  return J;
346 }
347 
348 
349 //=======================================================================
350 /// Derivatives of the shape functions and test functions w.r.t. to global
351 /// (Eulerian) coordinates. Return Jacobian of mapping between
352 /// local and global coordinates.
353 //=======================================================================
354 template<unsigned DIM>
355 inline double TCrouzeixRaviartElement<DIM>::
357  const unsigned &ipt, Shape &psi,
358  DShape &dpsidx, Shape &test,
359  DShape &dtestdx) const
360 {
361  //Call the geometrical shape functions and derivatives
362  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
363  //The test functions are the shape functions
364  test = psi;
365  dtestdx = dpsidx;
366  //Return the jacobian
367  return J;
368 }
369 
370 
371 //=======================================================================
372 /// 2D
373 /// Define the shape functions (psi) and test functions (test) and
374 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
375 /// and return Jacobian of mapping (J). Additionally compute the
376 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
377 ///
378 /// Galerkin: Test functions = shape functions
379 //=======================================================================
380 template<unsigned DIM>
381 inline double TCrouzeixRaviartElement<DIM>::
383  const unsigned &ipt, Shape &psi, DShape &dpsidx,
384  RankFourTensor<double> &d_dpsidx_dX,
385  Shape &test, DShape &dtestdx,
386  RankFourTensor<double> &d_dtestdx_dX,
387  DenseMatrix<double> &djacobian_dX) const
388  {
389  // Call the geometrical shape functions and derivatives
390  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
391  djacobian_dX,
392  d_dpsidx_dX);
393 
394  // Loop over the test functions and derivatives and set them equal to the
395  // shape functions
396  for(unsigned i=0;i<9;i++)
397  {
398  test[i] = psi[i];
399 
400  for(unsigned k=0;k<2;k++)
401  {
402  dtestdx(i,k) = dpsidx(i,k);
403 
404  for(unsigned p=0;p<2;p++)
405  {
406  for(unsigned q=0;q<9;q++)
407  {
408  d_dtestdx_dX(p,q,i,k) = d_dpsidx_dX(p,q,i,k);
409  }
410  }
411  }
412  }
413 
414  // Return the jacobian
415  return J;
416  }
417 
418 
419 //=======================================================================
420 /// 2D :
421 /// Pressure shape functions
422 //=======================================================================
423 template<>
425  Shape &psi) const
426 {
427  psi[0] = 1.0;
428  psi[1] = s[0];
429  psi[2] = s[1];
430 }
431 
432 //=======================================================================
433 /// Pressure shape and test functions
434 //=======================================================================
435 template<>
437  Shape &psi,
438  Shape &test) const
439 {
440  //Call the pressure shape functions
441  this->pshape_nst(s,psi);
442  //The test functions are the shape functions
443  test = psi;
444 }
445 
446 
447 //=======================================================================
448 /// 3D :
449 /// Pressure shape functions
450 //=======================================================================
451 template<>
453  Shape &psi)
454 const
455 {
456  psi[0] = 1.0;
457  psi[1] = s[0];
458  psi[2] = s[1];
459  psi[3] = s[2];
460 }
461 
462 
463 //=======================================================================
464 /// Pressure shape and test functions
465 //=======================================================================
466 template<>
468  Shape &psi,
469  Shape &test) const
470 {
471  //Call the pressure shape functions
472  this->pshape_nst(s,psi);
473  //The test functions are the shape functions
474  test = psi;
475 }
476 
477 
478 //==========================================================================
479 /// 2D :
480 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
481 /// Return Jacobian of mapping between local and global coordinates.
482 //==========================================================================
483 template<>
485  const Vector<double> &s,
486  Shape &ppsi,
487  DShape &dppsidx,
488  Shape &ptest,
489  DShape &dptestdx) const
490  {
491 
492  // Initalise with shape fcts and derivs. w.r.t. to local coordinates
493  ppsi[0] = 1.0;
494  ppsi[1] = s[0];
495  ppsi[2] = s[1];
496 
497  dppsidx(0,0) = 0.0;
498  dppsidx(1,0) = 1.0;
499  dppsidx(2,0) = 0.0;
500 
501  dppsidx(0,1) = 0.0;
502  dppsidx(1,1) = 0.0;
503  dppsidx(2,1) = 1.0;
504 
505 
506  //Get the values of the shape functions and their local derivatives
507  Shape psi(7);
508  DShape dpsi(7,2);
509  dshape_local(s,psi,dpsi);
510 
511  //Allocate memory for the inverse 2x2 jacobian
512  DenseMatrix<double> inverse_jacobian(2);
513 
514  //Now calculate the inverse jacobian
515  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
516 
517  //Now set the values of the derivatives to be derivs w.r.t. to the
518  // Eulerian coordinates
519  transform_derivatives(inverse_jacobian,dppsidx);
520 
521  //The test functions are equal to the shape functions
522  ptest = ppsi;
523  dptestdx = dppsidx;
524 
525  //Return the determinant of the jacobian
526  return det;
527  }
528 
529 
530 
531 //==========================================================================
532 /// 3D :
533 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
534 /// Return Jacobian of mapping between local and global coordinates.
535 //==========================================================================
536 template<>
538  const Vector<double> &s,
539  Shape &ppsi,
540  DShape &dppsidx,
541  Shape &ptest,
542  DShape &dptestdx) const
543  {
544 
545  // Initalise with shape fcts and derivs. w.r.t. to local coordinates
546  ppsi[0] = 1.0;
547  ppsi[1] = s[0];
548  ppsi[2] = s[1];
549  ppsi[3] = s[2];
550 
551  dppsidx(0,0) = 0.0;
552  dppsidx(1,0) = 1.0;
553  dppsidx(2,0) = 0.0;
554  dppsidx(3,0) = 0.0;
555 
556  dppsidx(0,1) = 0.0;
557  dppsidx(1,1) = 0.0;
558  dppsidx(2,1) = 1.0;
559  dppsidx(3,1) = 0.0;
560 
561  dppsidx(0,2) = 0.0;
562  dppsidx(1,2) = 0.0;
563  dppsidx(2,2) = 0.0;
564  dppsidx(3,2) = 1.0;
565 
566 
567  //Get the values of the shape functions and their local derivatives
568  Shape psi(11);
569  DShape dpsi(11,3);
570  dshape_local(s,psi,dpsi);
571 
572  // Allocate memory for the inverse 3x3 jacobian
573  DenseMatrix<double> inverse_jacobian(3);
574 
575  // Now calculate the inverse jacobian
576  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
577 
578  // Now set the values of the derivatives to be derivs w.r.t. to the
579  // Eulerian coordinates
580  transform_derivatives(inverse_jacobian,dppsidx);
581 
582  //The test functions are equal to the shape functions
583  ptest = ppsi;
584  dptestdx = dppsidx;
585 
586  // Return the determinant of the jacobian
587  return det;
588 
589  }
590 
591 
592 //=======================================================================
593 /// Face geometry of the 2D Crouzeix_Raviart elements
594 //=======================================================================
595 template<>
596 class FaceGeometry<TCrouzeixRaviartElement<2> >: public virtual TElement<1,3>
597 {
598 public:
599  FaceGeometry() : TElement<1,3>() {}
600 };
601 
602 //=======================================================================
603 /// Face geometry of the 3D Crouzeix_Raviart elements
604 //=======================================================================
605 template<>
607 public virtual TBubbleEnrichedElement<2,3>
608 {
609 
610  public:
612 };
613 
614 
615 //=======================================================================
616 /// Face geometry of the FaceGeometry of the 2D CrouzeixRaviart elements
617 //=======================================================================
618 template<>
620 public virtual PointElement
621 {
622  public:
624 };
625 
626 
627 //=======================================================================
628 /// Face geometry of the FaceGeometry of the 3D Crouzeix_Raviart elements
629 //=======================================================================
630 template<>
632 public virtual TElement<1,3>
633 {
634  public:
635  FaceGeometry() : TElement<1,3>() {}
636 };
637 
638 
639 
640 //=============================================================================
641 /// Create a list of pairs for all unknowns in this element,
642 /// so that the first entry in each pair contains the global equation
643 /// number of the unknown, while the second one contains the number
644 /// of the DOF that this unknown is associated with.
645 /// (Function can obviously only be called if the equation numbering
646 /// scheme has been set up.)
647 //=============================================================================
648 template<unsigned DIM>
650  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
651 {
652  // number of nodes
653  unsigned n_node = this->nnode();
654 
655  // number of pressure values
656  unsigned n_press = this->npres_nst();
657 
658  // temporary pair (used to store dof lookup prior to being added to list)
659  std::pair<unsigned,unsigned> dof_lookup;
660 
661  // pressure dof number
662  unsigned pressure_dof_number = DIM;
663 
664  // loop over the pressure values
665  for (unsigned n = 0; n < n_press; n++)
666  {
667  // determine local eqn number
668  int local_eqn_number = this->p_local_eqn(n);
669 
670  // ignore pinned values - far away degrees of freedom resulting
671  // from hanging nodes can be ignored since these are be dealt
672  // with by the element containing their master nodes
673  if (local_eqn_number >= 0)
674  {
675  // store dof lookup in temporary pair: First entry in pair
676  // is global equation number; second entry is dof type
677  dof_lookup.first = this->eqn_number(local_eqn_number);
678  dof_lookup.second = pressure_dof_number;
679 
680  // add to list
681  dof_lookup_list.push_front(dof_lookup);
682  }
683  }
684 
685  // loop over the nodes
686  for (unsigned n = 0; n < n_node; n++)
687  {
688  // find the number of values at this node
689  unsigned nv = this->node_pt(n)->nvalue();
690 
691  //loop over these values
692  for (unsigned v = 0; v < nv; v++)
693  {
694  // determine local eqn number
695  int local_eqn_number = this->nodal_local_eqn(n, v);
696 
697  // ignore pinned values
698  if (local_eqn_number >= 0)
699  {
700  // store dof lookup in temporary pair: First entry in pair
701  // is global equation number; second entry is dof type
702  dof_lookup.first = this->eqn_number(local_eqn_number);
703  dof_lookup.second = v;
704 
705  // add to list
706  dof_lookup_list.push_front(dof_lookup);
707 
708  }
709  }
710  }
711 }
712 
713 ////////////////////////////////////////////////////////////////////////////
714 ////////////////////////////////////////////////////////////////////////////
715 ////////////////////////////////////////////////////////////////////////////
716 
717 
718 
719 //=======================================================================
720 /// Taylor--Hood elements are Navier--Stokes elements
721 /// with quadratic interpolation for velocities and positions and
722 /// continous linear pressure interpolation
723 //=======================================================================
724 template <unsigned DIM>
725 class TTaylorHoodElement : public virtual TElement<DIM,3>,
726  public virtual NavierStokesEquations<DIM>,
727  public virtual ElementWithZ2ErrorEstimator
728 
729 {
730  private:
731 
732  /// Static array of ints to hold number of variables at node
733  static const unsigned Initial_Nvalue[];
734 
735  protected:
736 
737  /// \short Static array of ints to hold conversion from pressure
738  /// node numbers to actual node numbers
739  static const unsigned Pconv[];
740 
741  /// \short Velocity shape and test functions and their derivs
742  /// w.r.t. to global coords at local coordinate s (taken from geometry)
743  /// Return Jacobian of mapping between local and global coordinates.
744  inline double dshape_and_dtest_eulerian_nst(const Vector<double> &s,
745  Shape &psi,
746  DShape &dpsidx, Shape &test,
747  DShape &dtestdx) const;
748 
749  /// \short Velocity shape and test functions and their derivs
750  /// w.r.t. to global coords at local coordinate s (taken from geometry)
751  /// Return Jacobian of mapping between local and global coordinates.
752  inline double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt,
753  Shape &psi,
754  DShape &dpsidx,
755  Shape &test,
756  DShape &dtestdx) const;
757 
758  /// \short Shape/test functions and derivs w.r.t. to global coords at
759  /// integration point ipt; return Jacobian of mapping (J). Also compute
760  /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
762  const unsigned &ipt,
763  Shape &psi,
764  DShape &dpsidx,
765  RankFourTensor<double> &d_dpsidx_dX,
766  Shape &test,
767  DShape &dtestdx,
768  RankFourTensor<double> &d_dtestdx_dX,
769  DenseMatrix<double> &djacobian_dX) const;
770 
771  /// \short Compute the pressure shape and test functions and derivatives
772  /// w.r.t. global coords at local coordinate s.
773  /// Return Jacobian of mapping between local and global coordinates.
774  virtual double dpshape_and_dptest_eulerian_nst(const Vector<double> &s,
775  Shape &ppsi,
776  DShape &dppsidx,
777  Shape &ptest,
778  DShape &dptestdx) const;
779 
780  /// Unpin all pressure dofs
782 
783  /// Pin all nodal pressure dofs
785 
786  /// Unpin the proper nodal pressure dofs
788 
789 
790 public:
791 
792  /// Constructor, no internal data points
794 
795 
796  /// Broken copy constructor
798  {
799  BrokenCopy::broken_copy("TTaylorHoodElement");
800  }
801 
802  /// Broken assignment operator
803  /*void operator=(const TTaylorHoodElement<DIM>&)
804  {
805  BrokenCopy::broken_assign("TTaylorHoodElement");
806  }*/
807 
808  /// \short Number of values (pinned or dofs) required at node n. Can
809  /// be overwritten for hanging node version
810  inline virtual unsigned required_nvalue(const unsigned &n) const
811  {return Initial_Nvalue[n];}
812 
813  /// Test whether the pressure dof p_dof hanging or not?
814  //bool pressure_dof_is_hanging(const unsigned& p_dof)
815  // {return this->node_pt(Pconv[p_dof])->is_hanging(DIM);}
816 
817 
818  /// Pressure shape functions at local coordinate s
819  inline void pshape_nst(const Vector<double> &s, Shape &psi) const;
820 
821  /// Pressure shape and test functions at local coordinte s
822  inline void pshape_nst(const Vector<double> &s, Shape &psi,
823  Shape &test) const;
824 
825  /// \short Which nodal value represents the pressure?
826  unsigned p_index_nst() {return DIM;}
827 
828  /// \short Pointer to n_p-th pressure node
829  //Node* pressure_node_pt(const unsigned &n_p)
830  //{return this->Node_pt[Pconv[n_p]];}
831 
832  /// Return the local equation numbers for the pressure values.
833  inline int p_local_eqn(const unsigned &n) const
834  {return this->nodal_local_eqn(Pconv[n],DIM);}
835 
836  /// \short Access function for the pressure values at local pressure
837  /// node n_p (const version)
838  double p_nst(const unsigned &n_p) const
839  {return this->nodal_value(Pconv[n_p],DIM);}
840 
841  /// \short Access function for the pressure values at local pressure
842  /// node n_p (const version)
843  double p_nst(const unsigned &t, const unsigned &n_p) const
844  {return this->nodal_value(t,Pconv[n_p],DIM);}
845 
846  /// \short Set the value at which the pressure is stored in the nodes
847  int p_nodal_index_nst() const {return static_cast<int>(DIM);}
848 
849  /// Return number of pressure values
850  unsigned npres_nst() const;
851 
852  /// Pin p_dof-th pressure dof and set it to value specified by p_value.
853  void fix_pressure(const unsigned &p_dof, const double &p_value)
854  {
855  this->node_pt(Pconv[p_dof])->pin(DIM);
856  this->node_pt(Pconv[p_dof])->set_value(DIM,p_value);
857  }
858 
859 
860  /// \short Build FaceElements that apply the Robin boundary condition
861  /// to the pressure advection diffusion problem required by
862  /// Fp preconditioner
864  face_index)
865  {
868  this, face_index));
869  }
870 
871  /// \short Add to the set \c paired_load_data pairs containing
872  /// - the pointer to a Data object
873  /// and
874  /// - the index of the value in that Data object
875  /// .
876  /// for all values (pressures, velocities) that affect the
877  /// load computed in the \c get_load(...) function.
878  void identify_load_data(
879  std::set<std::pair<Data*,unsigned> > &paired_load_data);
880 
881  /// \short Add to the set \c paired_pressure_data pairs
882  /// containing
883  /// - the pointer to a Data object
884  /// and
885  /// - the index of the value in that Data object
886  /// .
887  /// for all pressure values that affect the
888  /// load computed in the \c get_load(...) function.
890  std::set<std::pair<Data*,unsigned> > &paired_pressure_data);
891 
892  /// Redirect output to NavierStokesEquations output
893  void output(std::ostream &outfile)
895 
896  /// Redirect output to NavierStokesEquations output
897  void output(std::ostream &outfile, const unsigned &nplot)
898  {NavierStokesEquations<DIM>::output(outfile,nplot);}
899 
900  /// Redirect output to NavierStokesEquations output
901  void output(FILE* file_pt) {NavierStokesEquations<DIM>::output(file_pt);}
902 
903  /// Redirect output to NavierStokesEquations output
904  void output(FILE* file_pt, const unsigned &n_plot)
905  {NavierStokesEquations<DIM>::output(file_pt,n_plot);}
906 
907  /// \short Order of recovery shape functions for Z2 error estimation:
908  /// Same order as shape functions.
909  unsigned nrecovery_order() {return 2;}
910 
911  /// \short Number of vertex nodes in the element
912  unsigned nvertex_node() const
913  {return DIM+1;}
914 
915  /// \short Pointer to the j-th vertex node in the element
916  Node* vertex_node_pt(const unsigned& j) const
917  {return node_pt(j);}
918 
919 
920  /// Number of 'flux' terms for Z2 error estimation
921  unsigned num_Z2_flux_terms()
922  {
923  // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
924  return DIM + (DIM*(DIM-1))/2;
925  }
926 
927  /// \short Get 'flux' for Z2 error recovery: Upper triangular entries
928  /// in strain rate tensor.
930  {
931 #ifdef PARANOID
932  unsigned num_entries=DIM+(DIM*(DIM-1))/2;
933  if (flux.size() < num_entries)
934  {
935  std::ostringstream error_message;
936  error_message << "The flux vector has the wrong number of entries, "
937  << flux.size() << ", whereas it should be at least "
938  << num_entries << std::endl;
939  throw OomphLibError(error_message.str(),
940  OOMPH_CURRENT_FUNCTION,
941  OOMPH_EXCEPTION_LOCATION);
942  }
943 #endif
944 
945  // Get strain rate matrix
946  DenseMatrix<double> strainrate(DIM);
947  this->strain_rate(s,strainrate);
948 
949  // Pack into flux Vector
950  unsigned icount=0;
951 
952  // Start with diagonal terms
953  for(unsigned i=0;i<DIM;i++)
954  {
955  flux[icount]=strainrate(i,i);
956  icount++;
957  }
958 
959  //Off diagonals row by row
960  for(unsigned i=0;i<DIM;i++)
961  {
962  for(unsigned j=i+1;j<DIM;j++)
963  {
964  flux[icount]=strainrate(i,j);
965  icount++;
966  }
967  }
968  }
969 
970  /// \short The number of "DOF types" that degrees of freedom in this element
971  /// are sub-divided into: Velocity and pressure.
972  unsigned ndof_types() const
973  {
974  return DIM+1;
975  }
976 
977  /// \short Create a list of pairs for all unknowns in this element,
978  /// so that the first entry in each pair contains the global equation
979  /// number of the unknown, while the second one contains the number
980  /// of the "DOF type" that this unknown is associated with.
981  /// (Function can obviously only be called if the equation numbering
982  /// scheme has been set up.) Velocity=0; Pressure=1
984  std::list<std::pair<unsigned long,unsigned> >& dof_lookup_list) const
985  {
986  // number of nodes
987  unsigned n_node = this->nnode();
988 
989  // temporary pair (used to store dof lookup prior to being added to list)
990  std::pair<unsigned,unsigned> dof_lookup;
991 
992  // loop over the nodes
993  for (unsigned n = 0; n < n_node; n++)
994  {
995  // find the number of Navier Stokes values at this node
996  unsigned nv = this->required_nvalue(n);
997 
998  //loop over these values
999  for (unsigned v = 0; v < nv; v++)
1000  {
1001  // determine local eqn number
1002  int local_eqn_number = this->nodal_local_eqn(n, v);
1003 
1004  // ignore pinned values - far away degrees of freedom resulting
1005  // from hanging nodes can be ignored since these are be dealt
1006  // with by the element containing their master nodes
1007  if (local_eqn_number >= 0)
1008  {
1009  // store dof lookup in temporary pair: Global equation number
1010  // is the first entry in pair
1011  dof_lookup.first = this->eqn_number(local_eqn_number);
1012 
1013  // set dof numbers: Dof number is the second entry in pair
1014  dof_lookup.second = v;
1015 
1016  // add to list
1017  dof_lookup_list.push_front(dof_lookup);
1018  }
1019  }
1020  }
1021  }
1022 };
1023 
1024 
1025 
1026 
1027 //Inline functions
1028 
1029 //==========================================================================
1030 /// 2D :
1031 /// Number of pressure values
1032 //==========================================================================
1033 template<>
1034 inline unsigned TTaylorHoodElement<2>::npres_nst() const
1035 {
1036  return 3;
1037 }
1038 
1039 //==========================================================================
1040 /// 3D :
1041 /// Number of pressure values
1042 //==========================================================================
1043 template<>
1044 inline unsigned TTaylorHoodElement<3>::npres_nst() const
1045 {
1046  return 4;
1047 }
1048 
1049 
1050 
1051 //==========================================================================
1052 /// 2D :
1053 /// Derivatives of the shape functions and test functions w.r.t to
1054 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1055 /// local and global coordinates.
1056 //==========================================================================
1057 template<unsigned DIM>
1059  const Vector<double> &s,
1060  Shape &psi,
1061  DShape &dpsidx, Shape &test,
1062  DShape &dtestdx) const
1063 {
1064  //Call the geometrical shape functions and derivatives
1065  double J = this->dshape_eulerian(s,psi,dpsidx);
1066  //Test functions are the shape functions
1067  test = psi;
1068  dtestdx = dpsidx;
1069  //Return the jacobian
1070  return J;
1071 }
1072 
1073 
1074 //==========================================================================
1075 /// Derivatives of the shape functions and test functions w.r.t to
1076 /// global (Eulerian) coordinates. Return Jacobian of mapping between
1077 /// local and global coordinates.
1078 //==========================================================================
1079 template<unsigned DIM>
1081  const unsigned &ipt,Shape &psi, DShape &dpsidx, Shape &test,
1082  DShape &dtestdx) const
1083 {
1084  //Call the geometrical shape functions and derivatives
1085  double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx);
1086  //Test functions are the shape functions
1087  test = psi;
1088  dtestdx = dpsidx;
1089  //Return the jacobian
1090  return J;
1091 }
1092 
1093 //==========================================================================
1094 /// 2D :
1095 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
1096 /// Return Jacobian of mapping between local and global coordinates.
1097 //==========================================================================
1098 template<>
1100  const Vector<double> &s,
1101  Shape &ppsi,
1102  DShape &dppsidx,
1103  Shape &ptest,
1104  DShape &dptestdx) const
1105  {
1106 
1107  ppsi[0] = s[0];
1108  ppsi[1] = s[1];
1109  ppsi[2] = 1.0-s[0]-s[1];
1110 
1111  dppsidx(0,0)=1.0;
1112  dppsidx(0,1)=0.0;
1113 
1114  dppsidx(1,0)=0.0;
1115  dppsidx(1,1)=1.0;
1116 
1117  dppsidx(2,0)=-1.0;
1118  dppsidx(2,1)=-1.0;
1119 
1120  // Allocate memory for the inverse 2x2 jacobian
1121  DenseMatrix<double> inverse_jacobian(2);
1122 
1123 
1124  //Get the values of the shape functions and their local derivatives
1125  Shape psi(6);
1126  DShape dpsi(6,2);
1127  dshape_local(s,psi,dpsi);
1128 
1129  // Now calculate the inverse jacobian
1130  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
1131 
1132  // Now set the values of the derivatives to be derivs w.r.t. to the
1133  // Eulerian coordinates
1134  transform_derivatives(inverse_jacobian,dppsidx);
1135 
1136  //Test functions are shape functions
1137  ptest = ppsi;
1138  dptestdx = dppsidx;
1139 
1140  // Return the determinant of the jacobian
1141  return det;
1142 
1143  }
1144 
1145 
1146 //==========================================================================
1147 /// 3D :
1148 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
1149 /// Return Jacobian of mapping between local and global coordinates.
1150 //==========================================================================
1151 template<>
1153  const Vector<double> &s,
1154  Shape &ppsi,
1155  DShape &dppsidx,
1156  Shape &ptest,
1157  DShape &dptestdx) const
1158  {
1159 
1160  ppsi[0] = s[0];
1161  ppsi[1] = s[1];
1162  ppsi[2] = s[2];
1163  ppsi[3] = 1.0-s[0]-s[1]-s[2];
1164 
1165  dppsidx(0,0)=1.0;
1166  dppsidx(0,1)=0.0;
1167  dppsidx(0,2)=0.0;
1168 
1169  dppsidx(1,0)=0.0;
1170  dppsidx(1,1)=1.0;
1171  dppsidx(1,2)=0.0;
1172 
1173  dppsidx(2,0)=0.0;
1174  dppsidx(2,1)=0.0;
1175  dppsidx(2,2)=1.0;
1176 
1177  dppsidx(3,0)=-1.0;
1178  dppsidx(3,1)=-1.0;
1179  dppsidx(3,2)=-1.0;
1180 
1181 
1182  //Get the values of the shape functions and their local derivatives
1183  Shape psi(10);
1184  DShape dpsi(10,3);
1185  dshape_local(s,psi,dpsi);
1186 
1187  // Allocate memory for the inverse 3x3 jacobian
1188  DenseMatrix<double> inverse_jacobian(3);
1189 
1190  // Now calculate the inverse jacobian
1191  const double det = local_to_eulerian_mapping(dpsi,inverse_jacobian);
1192 
1193  // Now set the values of the derivatives to be derivs w.r.t. to the
1194  // Eulerian coordinates
1195  transform_derivatives(inverse_jacobian,dppsidx);
1196 
1197  //Test functions are shape functions
1198  ptest = ppsi;
1199  dptestdx = dppsidx;
1200 
1201  // Return the determinant of the jacobian
1202  return det;
1203 
1204  }
1205 
1206 
1207 
1208 //==========================================================================
1209 /// 2D :
1210 /// Define the shape functions (psi) and test functions (test) and
1211 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
1212 /// and return Jacobian of mapping (J). Additionally compute the
1213 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
1214 ///
1215 /// Galerkin: Test functions = shape functions
1216 //==========================================================================
1217 template<>
1218  inline double TTaylorHoodElement<2>::
1220  const unsigned &ipt, Shape &psi, DShape &dpsidx,
1221  RankFourTensor<double> &d_dpsidx_dX,
1222  Shape &test, DShape &dtestdx,
1223  RankFourTensor<double> &d_dtestdx_dX,
1224  DenseMatrix<double> &djacobian_dX) const
1225  {
1226  // Call the geometrical shape functions and derivatives
1227  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
1228  djacobian_dX,
1229  d_dpsidx_dX);
1230 
1231  // Loop over the test functions and derivatives and set them equal to the
1232  // shape functions
1233  for(unsigned i=0;i<6;i++)
1234  {
1235  test[i] = psi[i];
1236 
1237  for(unsigned k=0;k<2;k++)
1238  {
1239  dtestdx(i,k) = dpsidx(i,k);
1240 
1241  for(unsigned p=0;p<2;p++)
1242  {
1243  for(unsigned q=0;q<6;q++)
1244  {
1245  d_dtestdx_dX(p,q,i,k) = d_dpsidx_dX(p,q,i,k);
1246  }
1247  }
1248  }
1249  }
1250 
1251  // Return the jacobian
1252  return J;
1253  }
1254 
1255 
1256 //==========================================================================
1257 /// 3D :
1258 /// Define the shape functions (psi) and test functions (test) and
1259 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
1260 /// and return Jacobian of mapping (J). Additionally compute the
1261 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
1262 ///
1263 /// Galerkin: Test functions = shape functions
1264 //==========================================================================
1265 template<>
1266  inline double TTaylorHoodElement<3>::
1268  const unsigned &ipt, Shape &psi, DShape &dpsidx,
1269  RankFourTensor<double> &d_dpsidx_dX,
1270  Shape &test, DShape &dtestdx,
1271  RankFourTensor<double> &d_dtestdx_dX,
1272  DenseMatrix<double> &djacobian_dX) const
1273  {
1274  // Call the geometrical shape functions and derivatives
1275  const double J = this->dshape_eulerian_at_knot(ipt,psi,dpsidx,
1276  djacobian_dX,
1277  d_dpsidx_dX);
1278 
1279  // Loop over the test functions and derivatives and set them equal to the
1280  // shape functions
1281  for(unsigned i=0;i<10;i++)
1282  {
1283  test[i] = psi[i];
1284 
1285  for(unsigned k=0;k<3;k++)
1286  {
1287  dtestdx(i,k) = dpsidx(i,k);
1288 
1289  for(unsigned p=0;p<3;p++)
1290  {
1291  for(unsigned q=0;q<10;q++)
1292  {
1293  d_dtestdx_dX(p,q,i,k) = d_dpsidx_dX(p,q,i,k);
1294  }
1295  }
1296  }
1297  }
1298 
1299  // Return the jacobian
1300  return J;
1301  }
1302 
1303 
1304 
1305 //==========================================================================
1306 /// 2D :
1307 /// Pressure shape functions
1308 //==========================================================================
1309 template<>
1311 const
1312 {
1313  psi[0] = s[0];
1314  psi[1] = s[1];
1315  psi[2] = 1.0-s[0]-s[1];
1316 }
1317 
1318 //==========================================================================
1319 /// 3D :
1320 /// Pressure shape functions
1321 //==========================================================================
1322 template<>
1324 const
1325 {
1326  psi[0] = s[0];
1327  psi[1] = s[1];
1328  psi[2] = s[2];
1329  psi[3] = 1.0-s[0]-s[1]-s[2];
1330 }
1331 
1332 
1333 //==========================================================================
1334 /// Pressure shape and test functions
1335 //==========================================================================
1336 template<unsigned DIM>
1338  Shape &psi,
1339  Shape &test) const
1340 {
1341  //Call the pressure shape functions
1342  this->pshape_nst(s,psi);
1343  //Test functions are shape functions
1344  test = psi;
1345 }
1346 
1347 
1348 //=======================================================================
1349 /// Face geometry of the 2D Taylor_Hood elements
1350 //=======================================================================
1351 template<>
1352 class FaceGeometry<TTaylorHoodElement<2> >: public virtual TElement<1,3>
1353 {
1354  public:
1355 
1356  /// Constructor: Call constructor of base
1357  FaceGeometry() : TElement<1,3>() {}
1358 };
1359 
1360 
1361 
1362 //=======================================================================
1363 /// Face geometry of the 3D Taylor_Hood elements
1364 //=======================================================================
1365 template<>
1366 class FaceGeometry<TTaylorHoodElement<3> >: public virtual TElement<2,3>
1367 {
1368 
1369  public:
1370 
1371  /// Constructor: Call constructor of base
1372  FaceGeometry() : TElement<2,3>() {}
1373 };
1374 
1375 
1376 
1377 //=======================================================================
1378 /// Face geometry of the FaceGeometry of the 2D TaylorHood elements
1379 //=======================================================================
1380 template<>
1382 public virtual PointElement
1383 {
1384  public:
1386 };
1387 
1388 
1389 //=======================================================================
1390 /// Face geometry of the FaceGeometry of the 3D Crouzeix_Raviart elements
1391 //=======================================================================
1392 template<>
1394 public virtual TElement<1,3>
1395 {
1396  public:
1397  FaceGeometry() : TElement<1,3>() {}
1398 };
1399 
1400 
1401 }
1402 
1403 #endif
unsigned npres_nst() const
Return number of pressure values.
TTaylorHoodElement(const TTaylorHoodElement< DIM > &dummy)
Broken copy constructor.
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
unsigned nvertex_node() const
Number of vertex nodes in the element.
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 broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
void pin_all_nodal_pressure_dofs()
Pin all nodal pressure dofs.
void identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
void full_output(std::ostream &outfile)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format...
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
Vector< FpPressureAdvDiffRobinBCElementBase * > Pressure_advection_diffusion_robin_element_pt
Storage for FaceElements that apply Robin BC for pressure adv diff equation used in Fp preconditioner...
FaceGeometry()
Constructor: Call constructor of base.
double dpshape_and_dptest_eulerian_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 output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor. ...
void pshape_nst(const Vector< double > &s, Shape &psi) const
Test whether the pressure dof p_dof hanging or not?
FaceGeometry()
Constructor: Call constructor of base.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
cstr elem_len * i
Definition: cfortran.h:607
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element.
void output(std::ostream &outfile, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
unsigned nvertex_node() const
Number of vertex nodes in the element.
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 full_output(std::ostream &outfile)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format...
unsigned npres_nst() const
Return number of pressure values.
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
int p_nodal_index_nst() const
Set the value at which the pressure is stored in the nodes.
char t
Definition: cfortran.h:572
unsigned p_index_nst()
Which nodal value represents the pressure?
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(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
virtual unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
void unpin_all_nodal_pressure_dofs()
Unpin all pressure dofs.
void full_output(std::ostream &outfile, const unsigned &nplot)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format...
A Rank 4 Tensor class.
Definition: matrices.h:1625
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs of pointers to data objects and unsignedegers that index the va...
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocity and ...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void unpin_proper_nodal_pressure_dofs()
Unpin the proper nodal pressure dofs.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
double p_nst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version) ...
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: 1/2 (du_i/dx_j + du_j/dx_i)
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
void output(FILE *file_pt, const unsigned &n_plot)
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.
void output(std::ostream &outfile, const unsigned &nplot)
Redirect output to NavierStokesEquations output.
static char t char * s
Definition: cfortran.h:572
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 set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:623
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor. ...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as unenriched shape functions...
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double p_nst(const unsigned &t, const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version) ...
virtual unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
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
virtual double dpshape_and_dptest_eulerian_nst(const Vector< double > &s, Shape &ppsi, DShape &dppsidx, Shape &ptest, DShape &dptestdx) const
Compute the pressure shape and test functions and derivatives.
TCrouzeixRaviartElement(const TCrouzeixRaviartElement< DIM > &dummy)
Broken copy constructor.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
double dshape_and_dtest_eulerian_at_knot_nst(const unsigned &ipt, Shape &psi, DShape &dpsidx, Shape &test, DShape &dtestdx) const
Velocity shape and test functions and their derivs w.r.t. to global coords at ipt-th integation point...
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
void unpin_all_internal_pressure_dofs()
Unpin all internal pressure dofs.
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 local coordinate s (tak...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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...
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
TCrouzeixRaviartElement()
Constructor, there are DIM+1 internal values (for the pressure)
int p_local_eqn(const unsigned &n) const
Return the local equation numbers for the pressure values.
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 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...
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...
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 identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
double p_nst(const unsigned &t, const unsigned &i) const
Return the pressure values at internal dof i_internal (Discontinous pressure interpolation – no need ...
double p_nst(const unsigned &i) const
Return the pressure values at internal dof i_internal (Discontinous pressure interpolation – no need ...
int p_local_eqn(const unsigned &n) const
Pointer to n_p-th pressure node.
TTaylorHoodElement()
Constructor, no internal data points.