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 Axisymmetric Navier Stokes elements
27
28#ifndef OOMPH_TAXISYM_NAVIER_STOKES_ELEMENTS_HEADER
29#define OOMPH_TAXISYM_NAVIER_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 /// AxisymmetricTCrouzeix_Raviart elements are Navier--Stokes elements with
55 /// quadratic interpolation for velocities and positions enriched by a single
56 /// cubic bubble function, but a discontinuous linear pressure interpolation
57 //==========================================================================
59 : public virtual TBubbleEnrichedElement<2, 3>,
61 public virtual ElementWithZ2ErrorEstimator
62 {
63 protected:
64 /// Internal index that indicates at which internal datum the pressure is
65 /// stored
67
68
69 /// Velocity shape and test functions and their derivs
70 /// w.r.t. to global coords at local coordinate s (taken from geometry)
71 /// Return Jacobian of mapping between local and global coordinates.
73 Shape& psi,
74 DShape& dpsidx,
75 Shape& test,
76 DShape& dtestdx) const;
77
78 /// Velocity shape and test functions and their derivs
79 /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
80 /// Return Jacobian of mapping between local and global coordinates.
82 const unsigned& ipt,
83 Shape& psi,
84 DShape& dpsidx,
85 Shape& test,
86 DShape& dtestdx) const;
87
88 /// Shape/test functions and derivs w.r.t. to global coords at
89 /// integration point ipt; return Jacobian of mapping (J). Also compute
90 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
92 const unsigned& ipt,
93 Shape& psi,
94 DShape& dpsidx,
95 RankFourTensor<double>& d_dpsidx_dX,
96 Shape& test,
97 DShape& dtestdx,
98 RankFourTensor<double>& d_dtestdx_dX,
99 DenseMatrix<double>& djacobian_dX) const;
100
101 /// Pressure shape and test functions and their derivs
102 /// w.r.t. to global coords at local coordinate s (taken from geometry)
103 /// Return Jacobian of mapping between local and global coordinates.
105 Shape& ppsi,
106 DShape& dppsidx,
107 Shape& ptest,
108 DShape& dptestdx) const;
109
110 public:
111 /// Pressure shape functions at local coordinate s
112 inline void pshape_axi_nst(const Vector<double>& s, Shape& psi) const;
113
114 /// Pressure shape and test functions at local coordinte s
115 inline void pshape_axi_nst(const Vector<double>& s,
116 Shape& psi,
117 Shape& test) const;
118
119 /// Unpin all internal pressure dofs
121
122 /// Return the local equation numbers for the pressure values.
123 inline int p_local_eqn(const unsigned& n) const
124 {
125 return this->internal_local_eqn(P_axi_nst_internal_index, n);
126 }
127
128 public:
129 /// Constructor, there are 3 internal values (for the pressure)
132 {
133 // Allocate and a single internal datum with 3 entries for the
134 // pressure
136 }
137
138 /// Broken copy constructor
140 const AxisymmetricTCrouzeixRaviartElement& dummy) = delete;
141
142 /// Broken assignment operator
143 // Commented out broken assignment operator because this can lead to a
144 // conflict warning when used in the virtual inheritence hierarchy.
145 // Essentially the compiler doesn't realise that two separate
146 // implementations of the broken function are the same and so, quite
147 // rightly, it shouts.
148 /*void operator=(const AxisymmetricTCrouzeixRaviartElement&) =
149 delete;*/
150
151
152 /// Number of values (pinned or dofs) required at local node n.
153 inline virtual unsigned required_nvalue(const unsigned& n) const
154 {
155 return 3;
156 }
157
158
159 /// Return the pressure values at internal dof i_internal
160 /// (Discontinous pressure interpolation -- no need to cater for hanging
161 /// nodes).
162 double p_axi_nst(const unsigned& i) const
163 {
164 return this->internal_data_pt(P_axi_nst_internal_index)->value(i);
165 }
166
167 /// Return number of pressure values
168 unsigned npres_axi_nst() const
169 {
170 return 3;
171 }
172
173 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
174 void fix_pressure(const unsigned& p_dof, const double& p_value)
175 {
176 this->internal_data_pt(P_axi_nst_internal_index)->pin(p_dof);
177 this->internal_data_pt(P_axi_nst_internal_index)
178 ->set_value(p_dof, p_value);
179 }
180
181 /// Build FaceElements that apply the Robin boundary condition
182 /// to the pressure advection diffusion problem required by
183 /// Fp preconditioner
184 // void build_fp_press_adv_diff_robin_bc_element(const unsigned&
185 // face_index)
186 // {
187 // this->Pressure_advection_diffusion_robin_element_pt.push_back(
188 // new
189 // FpPressureAdvDiffRobinBCElement<AxisymmetricTCrouzeixRaviartElement<DIM>
190 // >(
191 // this, face_index));
192 // }
193
194 /// Add to the set paired_load_data
195 /// pairs of pointers to data objects and unsignedegers that
196 /// index the values in the data object that affect the load (traction),
197 /// as specified in the get_load() function.
199 std::set<std::pair<Data*, unsigned>>& paired_load_data);
200
201 /// Add to the set \c paired_pressure_data pairs
202 /// containing
203 /// - the pointer to a Data object
204 /// and
205 /// - the index of the value in that Data object
206 /// .
207 /// for all pressure values that affect the
208 /// load computed in the \c get_load(...) function.
210 std::set<std::pair<Data*, unsigned>>& paired_pressure_data);
211
212 /// Redirect output to NavierStokesEquations output
213 void output(std::ostream& outfile)
214 {
216 }
217
218 /// Redirect output to NavierStokesEquations output
219 void output(std::ostream& outfile, const unsigned& nplot)
220 {
222 }
223
224 /// Redirect output to NavierStokesEquations output
225 void output(FILE* file_pt)
226 {
228 }
229
230 /// Redirect output to NavierStokesEquations output
231 void output(FILE* file_pt, const unsigned& n_plot)
232 {
234 }
235
236
237 /// Order of recovery shape functions for Z2 error estimation:
238 /// Same order as unenriched shape functions.
240 {
241 return 2;
242 }
243
244 /// Number of vertex nodes in the element
245 unsigned nvertex_node() const
246 {
247 return 3;
248 }
249
250 /// Pointer to the j-th vertex node in the element
251 Node* vertex_node_pt(const unsigned& j) const
252 {
253 return node_pt(j);
254 }
255
256 /// Number of 'flux' terms for Z2 error estimation
258 {
259 // 3 diagonal strain rates, 3 off diagonal
260 return 6;
261 }
262
263 /// Get 'flux' for Z2 error recovery: Upper triangular entries
264 /// in strain rate tensor.
266 {
267#ifdef PARANOID
268 unsigned num_entries = 6;
269 if (flux.size() < num_entries)
270 {
271 std::ostringstream error_message;
272 error_message << "The flux vector has the wrong number of entries, "
273 << flux.size() << ", whereas it should be at least "
274 << num_entries << std::endl;
275 throw OomphLibError(error_message.str(),
276 OOMPH_CURRENT_FUNCTION,
277 OOMPH_EXCEPTION_LOCATION);
278 }
279#endif
280
281 // Get strain rate matrix
282 DenseMatrix<double> strainrate(3);
283 this->strain_rate(s, strainrate);
284
285 // Pack into flux Vector
286 unsigned icount = 0;
287
288 // Start with diagonal terms
289 for (unsigned i = 0; i < 3; i++)
290 {
291 flux[icount] = strainrate(i, i);
292 icount++;
293 }
294
295 // Off diagonals row by row
296 for (unsigned i = 0; i < 3; i++)
297 {
298 for (unsigned j = i + 1; j < 3; j++)
299 {
300 flux[icount] = strainrate(i, j);
301 icount++;
302 }
303 }
304 }
305
306
307 /// The number of "DOF types" that degrees of freedom in this element
308 /// are sub-divided into: Velocity (3 components) and pressure.
309 unsigned ndof_types() const
310 {
311 return 4;
312 }
313
314 /// Create a list of pairs for all unknowns in this element,
315 /// so that the first entry in each pair contains the global equation
316 /// number of the unknown, while the second one contains the number
317 /// of the "DOF type" that this unknown is associated with.
318 /// (Function can obviously only be called if the equation numbering
319 /// scheme has been set up.)
321 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
322 {
323 // number of nodes
324 unsigned n_node = this->nnode();
325
326 // number of pressure values
327 unsigned n_press = this->npres_axi_nst();
328
329 // temporary pair (used to store dof lookup prior to being added to list)
330 std::pair<unsigned, unsigned> dof_lookup;
331
332 // pressure dof number
333 unsigned pressure_dof_number = 3;
334
335 // loop over the pressure values
336 for (unsigned n = 0; n < n_press; n++)
337 {
338 // determine local eqn number
339 int local_eqn_number = this->p_local_eqn(n);
340
341 // ignore pinned values - far away degrees of freedom resulting
342 // from hanging nodes can be ignored since these are be dealt
343 // with by the element containing their master nodes
344 if (local_eqn_number >= 0)
345 {
346 // store dof lookup in temporary pair: First entry in pair
347 // is global equation number; second entry is dof type
348 dof_lookup.first = this->eqn_number(local_eqn_number);
349 dof_lookup.second = pressure_dof_number;
350
351 // add to list
352 dof_lookup_list.push_front(dof_lookup);
353 }
354 }
355
356 // loop over the nodes
357 for (unsigned n = 0; n < n_node; n++)
358 {
359 // find the number of values at this node
360 unsigned nv = this->node_pt(n)->nvalue();
361
362 // loop over these values
363 for (unsigned v = 0; v < nv; v++)
364 {
365 // determine local eqn number
366 int local_eqn_number = this->nodal_local_eqn(n, v);
367
368 // ignore pinned values
369 if (local_eqn_number >= 0)
370 {
371 // store dof lookup in temporary pair: First entry in pair
372 // is global equation number; second entry is dof type
373 dof_lookup.first = this->eqn_number(local_eqn_number);
374 dof_lookup.second = v;
375
376 // add to list
377 dof_lookup_list.push_front(dof_lookup);
378 }
379 }
380 }
381 }
382 };
383
384 // Inline functions
385
386 //=======================================================================
387 /// Derivatives of the shape functions and test functions w.r.t. to global
388 /// (Eulerian) coordinates. Return Jacobian of mapping between
389 /// local and global coordinates.
390 //=======================================================================
393 Shape& psi,
394 DShape& dpsidx,
395 Shape& test,
396 DShape& dtestdx) const
397 {
398 // Call the geometrical shape functions and derivatives
399 double J = this->dshape_eulerian(s, psi, dpsidx);
400 // The test functions are equal to the shape functions
401 test = psi;
402 dtestdx = dpsidx;
403 // Return the jacobian
404 return J;
405 }
406
407
408 //=======================================================================
409 /// Derivatives of the shape functions and test functions w.r.t. to global
410 /// (Eulerian) coordinates. Return Jacobian of mapping between
411 /// local and global coordinates.
412 //=======================================================================
415 Shape& psi,
416 DShape& dpsidx,
417 Shape& test,
418 DShape& dtestdx) const
419 {
420 // Call the geometrical shape functions and derivatives
421 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
422 // The test functions are the shape functions
423 test = psi;
424 dtestdx = dpsidx;
425 // Return the jacobian
426 return J;
427 }
428
429
430 //=======================================================================
431 /// Define the shape functions (psi) and test functions (test) and
432 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
433 /// and return Jacobian of mapping (J). Additionally compute the
434 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
435 ///
436 /// Galerkin: Test functions = shape functions
437 //=======================================================================
440 const unsigned& ipt,
441 Shape& psi,
442 DShape& dpsidx,
443 RankFourTensor<double>& d_dpsidx_dX,
444 Shape& test,
445 DShape& dtestdx,
446 RankFourTensor<double>& d_dtestdx_dX,
447 DenseMatrix<double>& djacobian_dX) const
448 {
449 // Call the geometrical shape functions and derivatives
450 const double J = this->dshape_eulerian_at_knot(
451 ipt, psi, dpsidx, djacobian_dX, d_dpsidx_dX);
452
453 // Set the test functions equal to the shape functions
454 test = psi;
455 dtestdx = dpsidx;
456 d_dtestdx_dX = d_dpsidx_dX;
457
458 // Return the jacobian
459 return J;
460 }
461
462
463 //=======================================================================
464 /// Pressure shape functions
465 //=======================================================================
467 const Vector<double>& s, Shape& psi) const
468 {
469 psi[0] = 1.0;
470 psi[1] = s[0];
471 psi[2] = s[1];
472 }
473
474 //=======================================================================
475 /// Pressure shape and test functions
476 //=======================================================================
478 const Vector<double>& s, Shape& psi, Shape& test) const
479 {
480 // Call the pressure shape functions
481 this->pshape_axi_nst(s, psi);
482 // The test functions are the shape functions
483 test = psi;
484 }
485
486 //==========================================================================
487 /// 2D :
488 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
489 /// Return Jacobian of mapping between local and global coordinates.
490 //==========================================================================
493 Shape& ppsi,
494 DShape& dppsidx,
495 Shape& ptest,
496 DShape& dptestdx) const
497 {
498 // Initalise with shape fcts and derivs. w.r.t. to local coordinates
499 ppsi[0] = 1.0;
500 ppsi[1] = s[0];
501 ppsi[2] = s[1];
502
503 dppsidx(0, 0) = 0.0;
504 dppsidx(1, 0) = 1.0;
505 dppsidx(2, 0) = 0.0;
506
507 dppsidx(0, 1) = 0.0;
508 dppsidx(1, 1) = 0.0;
509 dppsidx(2, 1) = 1.0;
510
511
512 // Get the values of the shape functions and their local derivatives
513 Shape psi(7);
514 DShape dpsi(7, 2);
515 dshape_local(s, psi, dpsi);
516
517 // Allocate memory for the inverse 2x2 jacobian
518 DenseMatrix<double> inverse_jacobian(2);
519
520 // Now calculate the inverse jacobian
521 const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
522
523 // Now set the values of the derivatives to be derivs w.r.t. to the
524 // Eulerian coordinates
525 transform_derivatives(inverse_jacobian, dppsidx);
526
527 // The test functions are equal to the shape functions
528 ptest = ppsi;
529 dptestdx = dppsidx;
530
531 // Return the determinant of the jacobian
532 return det;
533 }
534
535
536 //=======================================================================
537 /// Face geometry of the 2D Crouzeix_Raviart elements
538 //=======================================================================
539 template<>
541 : public virtual TElement<1, 3>
542 {
543 public:
544 FaceGeometry() : TElement<1, 3>() {}
545 };
546
547
548 //=======================================================================
549 /// Face geometry of the FaceGeometry of the 2D CrouzeixRaviart elements
550 //=======================================================================
551 template<>
553 : public virtual PointElement
554 {
555 public:
557 };
558
559
560 /// /////////////////////////////////////////////////////////////////////////
561 /// /////////////////////////////////////////////////////////////////////////
562 /// /////////////////////////////////////////////////////////////////////////
563
564
565 //=======================================================================
566 /// Taylor--Hood elements are Navier--Stokes elements
567 /// with quadratic interpolation for velocities and positions and
568 /// continous linear pressure interpolation
569 //=======================================================================
571 : public virtual TElement<2, 3>,
573 public virtual ElementWithZ2ErrorEstimator
574
575 {
576 private:
577 /// Static array of ints to hold number of variables at node
578 static const unsigned Initial_Nvalue[];
579
580 protected:
581 /// Static array of ints to hold conversion from pressure
582 /// node numbers to actual node numbers
583 static const unsigned Pconv[];
584
585 /// Velocity shape and test functions and their derivs
586 /// w.r.t. to global coords at local coordinate s (taken from geometry)
587 /// Return Jacobian of mapping between local and global coordinates.
589 Shape& psi,
590 DShape& dpsidx,
591 Shape& test,
592 DShape& dtestdx) const;
593
594 /// Velocity shape and test functions and their derivs
595 /// w.r.t. to global coords at local coordinate s (taken from geometry)
596 /// Return Jacobian of mapping between local and global coordinates.
598 const unsigned& ipt,
599 Shape& psi,
600 DShape& dpsidx,
601 Shape& test,
602 DShape& dtestdx) const;
603
604 /// Shape/test functions and derivs w.r.t. to global coords at
605 /// integration point ipt; return Jacobian of mapping (J). Also compute
606 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
608 const unsigned& ipt,
609 Shape& psi,
610 DShape& dpsidx,
611 RankFourTensor<double>& d_dpsidx_dX,
612 Shape& test,
613 DShape& dtestdx,
614 RankFourTensor<double>& d_dtestdx_dX,
615 DenseMatrix<double>& djacobian_dX) const;
616
617 /// Compute the pressure shape and test functions and derivatives
618 /// w.r.t. global coords at local coordinate s.
619 /// Return Jacobian of mapping between local and global coordinates.
621 Shape& ppsi,
622 DShape& dppsidx,
623 Shape& ptest,
624 DShape& dptestdx) const;
625
626 /// Unpin all pressure dofs
628
629 /// Pin all nodal pressure dofs
631
632 /// Unpin the proper nodal pressure dofs
634
635
636 public:
637 /// Constructor, no internal data points
640 {
641 }
642
643
644 /// Broken copy constructor
646 const AxisymmetricTTaylorHoodElement& dummy) = delete;
647
648 /// Broken assignment operator
649 /*void operator=(const AxisymmetricTTaylorHoodElement&) =
650 delete;*/
651
652 /// Number of values (pinned or dofs) required at node n. Can
653 /// be overwritten for hanging node version
654 inline virtual unsigned required_nvalue(const unsigned& n) const
655 {
656 return Initial_Nvalue[n];
657 }
658
659 /// Test whether the pressure dof p_dof hanging or not?
660 // bool pressure_dof_is_hanging(const unsigned& p_dof)
661 // {return this->node_pt(Pconv[p_dof])->is_hanging(DIM);}
662
663
664 /// Pressure shape functions at local coordinate s
665 inline void pshape_axi_nst(const Vector<double>& s, Shape& psi) const;
666
667 /// Pressure shape and test functions at local coordinte s
668 inline void pshape_axi_nst(const Vector<double>& s,
669 Shape& psi,
670 Shape& test) const;
671
672 /// Which nodal value represents the pressure?
674 {
675 return 3;
676 }
677
678 /// Pointer to n_p-th pressure node
679 // Node* pressure_node_pt(const unsigned &n_p)
680 //{return this->Node_pt[Pconv[n_p]];}
681
682 /// Return the local equation numbers for the pressure values.
683 inline int p_local_eqn(const unsigned& n) const
684 {
685 return this->nodal_local_eqn(Pconv[n], 3);
686 }
687
688 /// Access function for the pressure values at local pressure
689 /// node n_p (const version)
690 double p_axi_nst(const unsigned& n_p) const
691 {
692 return this->nodal_value(Pconv[n_p], 3);
693 }
694
695 /// Set the value at which the pressure is stored in the nodes
697 {
698 return static_cast<int>(3);
699 }
700
701 /// Return number of pressure values
702 unsigned npres_axi_nst() const
703 {
704 return 3;
705 }
706
707 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
708 void fix_pressure(const unsigned& p_dof, const double& p_value)
709 {
710 this->node_pt(Pconv[p_dof])->pin(3);
711 this->node_pt(Pconv[p_dof])->set_value(3, p_value);
712 }
713
714
715 /// Build FaceElements that apply the Robin boundary condition
716 /// to the pressure advection diffusion problem required by
717 /// Fp preconditioner
718 // void build_fp_press_adv_diff_robin_bc_element(const unsigned&
719 // face_index)
720 // {
721 // this->Pressure_advection_diffusion_robin_element_pt.push_back(
722 // new FpPressureAdvDiffRobinBCElement<AxisymmetricTTaylorHoodElement<DIM>
723 // >(
724 // this, face_index));
725 // }
726
727 /// Add to the set \c paired_load_data pairs containing
728 /// - the pointer to a Data object
729 /// and
730 /// - the index of the value in that Data object
731 /// .
732 /// for all values (pressures, velocities) that affect the
733 /// load computed in the \c get_load(...) function.
735 std::set<std::pair<Data*, unsigned>>& paired_load_data);
736
737 /// Add to the set \c paired_pressure_data pairs
738 /// containing
739 /// - the pointer to a Data object
740 /// and
741 /// - the index of the value in that Data object
742 /// .
743 /// for all pressure values that affect the
744 /// load computed in the \c get_load(...) function.
746 std::set<std::pair<Data*, unsigned>>& paired_pressure_data);
747
748 /// Redirect output to NavierStokesEquations output
749 void output(std::ostream& outfile)
750 {
752 }
753
754 /// Redirect output to NavierStokesEquations output
755 void output(std::ostream& outfile, const unsigned& nplot)
756 {
758 }
759
760 /// Redirect output to NavierStokesEquations output
761 void output(FILE* file_pt)
762 {
764 }
765
766 /// Redirect output to NavierStokesEquations output
767 void output(FILE* file_pt, const unsigned& n_plot)
768 {
770 }
771
772 /// Order of recovery shape functions for Z2 error estimation:
773 /// Same order as shape functions.
775 {
776 return 2;
777 }
778
779 /// Number of vertex nodes in the element
780 unsigned nvertex_node() const
781 {
782 return 3;
783 }
784
785 /// Pointer to the j-th vertex node in the element
786 Node* vertex_node_pt(const unsigned& j) const
787 {
788 return node_pt(j);
789 }
790
791
792 /// Number of 'flux' terms for Z2 error estimation
794 {
795 // 3 diagonal strain rates, 3 off diagonal rates
796 return 6;
797 }
798
799 /// Get 'flux' for Z2 error recovery: Upper triangular entries
800 /// in strain rate tensor.
802 {
803#ifdef PARANOID
804 unsigned num_entries = 6;
805 if (flux.size() < num_entries)
806 {
807 std::ostringstream error_message;
808 error_message << "The flux vector has the wrong number of entries, "
809 << flux.size() << ", whereas it should be at least "
810 << num_entries << std::endl;
811 throw OomphLibError(error_message.str(),
812 OOMPH_CURRENT_FUNCTION,
813 OOMPH_EXCEPTION_LOCATION);
814 }
815#endif
816
817 // Get strain rate matrix
818 DenseMatrix<double> strainrate(3);
819 this->strain_rate(s, strainrate);
820
821 // Pack into flux Vector
822 unsigned icount = 0;
823
824 // Start with diagonal terms
825 for (unsigned i = 0; i < 3; i++)
826 {
827 flux[icount] = strainrate(i, i);
828 icount++;
829 }
830
831 // Off diagonals row by row
832 for (unsigned i = 0; i < 3; i++)
833 {
834 for (unsigned j = i + 1; j < 3; j++)
835 {
836 flux[icount] = strainrate(i, j);
837 icount++;
838 }
839 }
840 }
841
842 /// The number of "DOF types" that degrees of freedom in this element
843 /// are sub-divided into: Velocities (3 components) and pressure.
844 unsigned ndof_types() const
845 {
846 return 4;
847 }
848
849 /// Create a list of pairs for all unknowns in this element,
850 /// so that the first entry in each pair contains the global equation
851 /// number of the unknown, while the second one contains the number
852 /// of the "DOF type" that this unknown is associated with.
853 /// (Function can obviously only be called if the equation numbering
854 /// scheme has been set up.)
856 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
857 {
858 // number of nodes
859 unsigned n_node = this->nnode();
860
861 // temporary pair (used to store dof lookup prior to being added to list)
862 std::pair<unsigned, unsigned> dof_lookup;
863
864 // loop over the nodes
865 for (unsigned n = 0; n < n_node; n++)
866 {
867 // find the number of Navier Stokes values at this node
868 unsigned nv = this->required_nvalue(n);
869
870 // loop over these values
871 for (unsigned v = 0; v < nv; v++)
872 {
873 // determine local eqn number
874 int local_eqn_number = this->nodal_local_eqn(n, v);
875
876 // ignore pinned values - far away degrees of freedom resulting
877 // from hanging nodes can be ignored since these are be dealt
878 // with by the element containing their master nodes
879 if (local_eqn_number >= 0)
880 {
881 // store dof lookup in temporary pair: Global equation number
882 // is the first entry in pair
883 dof_lookup.first = this->eqn_number(local_eqn_number);
884
885 // set dof numbers: Dof number is the second entry in pair
886 dof_lookup.second = v;
887
888 // add to list
889 dof_lookup_list.push_front(dof_lookup);
890 }
891 }
892 }
893 }
894 };
895
896
897 // Inline functions
898
899 //==========================================================================
900 /// Derivatives of the shape functions and test functions w.r.t to
901 /// global (Eulerian) coordinates. Return Jacobian of mapping between
902 /// local and global coordinates.
903 //==========================================================================
906 Shape& psi,
907 DShape& dpsidx,
908 Shape& test,
909 DShape& dtestdx) const
910 {
911 // Call the geometrical shape functions and derivatives
912 double J = this->dshape_eulerian(s, psi, dpsidx);
913 // Test functions are the shape functions
914 test = psi;
915 dtestdx = dpsidx;
916 // Return the jacobian
917 return J;
918 }
919
920
921 //==========================================================================
922 /// Derivatives of the shape functions and test functions w.r.t to
923 /// global (Eulerian) coordinates. Return Jacobian of mapping between
924 /// local and global coordinates.
925 //==========================================================================
928 Shape& psi,
929 DShape& dpsidx,
930 Shape& test,
931 DShape& dtestdx) const
932 {
933 // Call the geometrical shape functions and derivatives
934 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
935 // Test functions are the shape functions
936 test = psi;
937 dtestdx = dpsidx;
938 // Return the jacobian
939 return J;
940 }
941
942 //==========================================================================
943 /// Define the shape functions (psi) and test functions (test) and
944 /// their derivatives w.r.t. global coordinates (dpsidx and dtestdx)
945 /// and return Jacobian of mapping (J). Additionally compute the
946 /// derivatives of dpsidx, dtestdx and J w.r.t. nodal coordinates.
947 ///
948 /// Galerkin: Test functions = shape functions
949 //==========================================================================
952 const unsigned& ipt,
953 Shape& psi,
954 DShape& dpsidx,
955 RankFourTensor<double>& d_dpsidx_dX,
956 Shape& test,
957 DShape& dtestdx,
958 RankFourTensor<double>& d_dtestdx_dX,
959 DenseMatrix<double>& djacobian_dX) const
960 {
961 // Call the geometrical shape functions and derivatives
962 const double J = this->dshape_eulerian_at_knot(
963 ipt, psi, dpsidx, djacobian_dX, d_dpsidx_dX);
964
965 // Set the test functions equal to the shape functions
966 test = psi;
967 dtestdx = dpsidx;
968 d_dtestdx_dX = d_dpsidx_dX;
969
970 // Return the jacobian
971 return J;
972 }
973
974 //==========================================================================
975 /// Pressure shape and test functions and derivs w.r.t. to Eulerian coords.
976 /// Return Jacobian of mapping between local and global coordinates.
977 //==========================================================================
980 Shape& ppsi,
981 DShape& dppsidx,
982 Shape& ptest,
983 DShape& dptestdx) const
984 {
985 ppsi[0] = s[0];
986 ppsi[1] = s[1];
987 ppsi[2] = 1.0 - s[0] - s[1];
988
989 dppsidx(0, 0) = 1.0;
990 dppsidx(0, 1) = 0.0;
991
992 dppsidx(1, 0) = 0.0;
993 dppsidx(1, 1) = 1.0;
994
995 dppsidx(2, 0) = -1.0;
996 dppsidx(2, 1) = -1.0;
997
998 // Allocate memory for the inverse 2x2 jacobian
999 DenseMatrix<double> inverse_jacobian(2);
1000
1001
1002 // Get the values of the shape functions and their local derivatives
1003 Shape psi(6);
1004 DShape dpsi(6, 2);
1005 dshape_local(s, psi, dpsi);
1006
1007 // Now calculate the inverse jacobian
1008 const double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
1009
1010 // Now set the values of the derivatives to be derivs w.r.t. to the
1011 // Eulerian coordinates
1012 transform_derivatives(inverse_jacobian, dppsidx);
1013
1014 // Test functions are shape functions
1015 ptest = ppsi;
1016 dptestdx = dppsidx;
1017
1018 // Return the determinant of the jacobian
1019 return det;
1020 }
1021
1022
1023 //==========================================================================
1024 /// Pressure shape functions
1025 //==========================================================================
1027 const Vector<double>& s, Shape& psi) const
1028 {
1029 psi[0] = s[0];
1030 psi[1] = s[1];
1031 psi[2] = 1.0 - s[0] - s[1];
1032 }
1033
1034 //==========================================================================
1035 /// Pressure shape and test functions
1036 //==========================================================================
1038 const Vector<double>& s, Shape& psi, Shape& test) const
1039 {
1040 // Call the pressure shape functions
1041 this->pshape_axi_nst(s, psi);
1042 // Test functions are shape functions
1043 test = psi;
1044 }
1045
1046
1047 //=======================================================================
1048 /// Face geometry of the Axisymmetric Taylor_Hood elements
1049 //=======================================================================
1050 template<>
1052 : public virtual TElement<1, 3>
1053 {
1054 public:
1055 /// Constructor: Call constructor of base
1056 FaceGeometry() : TElement<1, 3>() {}
1057 };
1058
1059
1060 //=======================================================================
1061 /// Face geometry of the FaceGeometry of the Axisymmetric TaylorHood elements
1062 //=======================================================================
1063 template<>
1065 : public virtual PointElement
1066 {
1067 public:
1069 };
1070
1071
1072} // namespace oomph
1073
1074#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A class for elements that solve the unsteady axisymmetric Navier–Stokes equations in cylindrical pola...
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
void output(std::ostream &outfile)
Output function: x,y,[z],u,v,[w],p in tecplot format. Default number of plot points.
///////////////////////////////////////////////////////////////////////////
unsigned P_axi_nst_internal_index
Internal index that indicates at which internal datum the pressure is stored.
AxisymmetricTCrouzeixRaviartElement()
Constructor, there are 3 internal values (for the pressure)
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...
unsigned nvertex_node() const
Number of vertex nodes in the element.
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
unsigned npres_axi_nst() const
Return number of pressure values.
virtual unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
double p_axi_nst(const unsigned &i) const
Return the pressure values at internal dof i_internal (Discontinous pressure interpolation – no need ...
void fix_pressure(const unsigned &p_dof, const double &p_value)
Pin p_dof-th pressure dof and set it to value specified by p_value.
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
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...
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor.
AxisymmetricTCrouzeixRaviartElement(const AxisymmetricTCrouzeixRaviartElement &dummy)=delete
Broken copy constructor.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as unenriched shape functions.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocity (3 c...
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
void unpin_all_internal_pressure_dofs()
Unpin all internal pressure dofs.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void output(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.
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 identify_pressure_data(std::set< std::pair< Data *, unsigned > > &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
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.
///////////////////////////////////////////////////////////////////////// ///////////////////////////...
void unpin_all_nodal_pressure_dofs()
Unpin all pressure dofs.
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
void output(FILE *file_pt)
Redirect output to NavierStokesEquations output.
AxisymmetricTTaylorHoodElement()
Constructor, no internal data points.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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 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.
int p_nodal_index_axi_nst() const
Set the value at which the pressure is stored in the nodes.
unsigned nvertex_node() const
Number of vertex nodes in the element.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
unsigned p_index_axi_nst()
Which nodal value represents the pressure?
unsigned npres_axi_nst() const
Return number of pressure values.
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 identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
void output(std::ostream &outfile)
Redirect output to NavierStokesEquations output.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Velocities (3...
int p_local_eqn(const unsigned &n) const
Pointer to n_p-th pressure node.
virtual unsigned required_nvalue(const unsigned &n) const
Broken assignment operator.
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.
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...
void pshape_axi_nst(const Vector< double > &s, Shape &psi) const
Test whether the pressure dof p_dof hanging or not?
static const unsigned Initial_Nvalue[]
Static array of ints to hold number of variables at node.
double p_axi_nst(const unsigned &n_p) const
Access function for the pressure values at local pressure node n_p (const version)
void pin_all_nodal_pressure_dofs()
Pin all nodal pressure dofs.
void output(FILE *file_pt, const unsigned &n_plot)
Redirect output to NavierStokesEquations output.
void unpin_proper_nodal_pressure_dofs()
Unpin the proper nodal pressure dofs.
AxisymmetricTTaylorHoodElement(const AxisymmetricTTaylorHoodElement &dummy)=delete
Broken copy constructor.
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_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...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
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
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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...