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