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