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