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