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