refineable_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 refineable 2D quad Navier Stokes elements
27
28#ifndef OOMPH_REFINEABLE_NAVIER_STOKES_ELEMENTS_HEADER
29#define OOMPH_REFINEABLE_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// Oomph-lib headers
37#include "../generic/refineable_quad_element.h"
38#include "../generic/refineable_brick_element.h"
39#include "../generic/hp_refineable_elements.h"
40#include "../generic/error_estimator.h"
42
43namespace oomph
44{
45 /// ////////////////////////////////////////////////////////////////////
46 /// ////////////////////////////////////////////////////////////////////
47 /// ////////////////////////////////////////////////////////////////////
48
49
50 //======================================================================
51 /// A class for elements that allow the imposition of Robin boundary
52 /// conditions for the pressure advection diffusion problem in the
53 /// Fp preconditioner.
54 /// The geometrical information can be read from the FaceGeometery<ELEMENT>
55 /// class and and thus, we can be generic enough without the need to have
56 /// a separate equations class
57 //======================================================================
58 template<class ELEMENT>
60 : public virtual FpPressureAdvDiffRobinBCElement<ELEMENT>
61 {
62 public:
63 /// Constructor, which takes a "bulk" element and the value of the index
64 /// and its limit
66 const int& face_index)
67 : FaceGeometry<ELEMENT>(),
69 FpPressureAdvDiffRobinBCElement<ELEMENT>(element_pt, face_index, true)
70 {
71 }
72
73 /// This function returns the residuals for the
74 /// traction function. flag=1 (or 0): do (or don't) compute the
75 /// Jacobian as well.
77 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
78 };
79
80
81 //============================================================================
82 /// Get residuals and Jacobian of Robin boundary conditions in pressure
83 /// advection diffusion problem in Fp preconditoner. Refineable version.
84 //============================================================================
85 template<class ELEMENT>
88 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
89 {
90 // Pointers to hang info objects
91 HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
92
93 // Storage for local coordinates in FaceElement and associted bulk element
94 unsigned my_dim = this->dim();
95 Vector<double> s(my_dim);
96 Vector<double> s_bulk(my_dim + 1);
97
98 // Storage for outer unit normal
99 Vector<double> unit_normal(my_dim + 1);
100
101 // Storage for veloc in bulk element
102 Vector<double> veloc(my_dim + 1);
103
104 // Set the value of n_intpt
105 unsigned n_intpt = this->integral_pt()->nweight();
106
107 // Integers to store local equation numbers
108 int local_eqn = 0;
109 int local_unknown = 0;
110
111 // Get cast bulk element
112 ELEMENT* bulk_el_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
113
114 // Find out how many pressure dofs there are in the bulk element
115 unsigned n_pres = bulk_el_pt->npres_nst();
116
117
118 // Which nodal value represents the pressure? (Negative if pressure
119 // is not based on nodal interpolation).
120 int p_index = bulk_el_pt->p_nodal_index_nst();
121
122 // Local array of booleans that are true if the l-th pressure value is
123 // hanging (avoid repeated virtual function calls)
124 bool pressure_dof_is_hanging[n_pres];
125 // If the pressure is stored at a node
126 if (p_index >= 0)
127 {
128 // Read out whether the pressure is hanging
129 for (unsigned l = 0; l < n_pres; ++l)
130 {
131 pressure_dof_is_hanging[l] =
132 bulk_el_pt->pressure_node_pt(l)->is_hanging(p_index);
133 }
134 }
135 // Otherwise the pressure is not stored at a node and so cannot hang
136 else
137 {
138 // pressure advection diffusion doesn't work for this one!
139 throw OomphLibError(
140 "Pressure advection diffusion does not work in this case\n",
141 OOMPH_CURRENT_FUNCTION,
142 OOMPH_EXCEPTION_LOCATION);
143
144 for (unsigned l = 0; l < n_pres; ++l)
145 {
146 pressure_dof_is_hanging[l] = false;
147 }
148 }
149
150 // Get the Reynolds number from the bulk element
151 double re = bulk_el_pt->re();
152
153 // Set up memory for pressure shape and test functions
154 Shape psip(n_pres), testp(n_pres);
155
156 // Loop over the integration points
157 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
158 {
159 // Get the integral weight
160 double w = this->integral_pt()->weight(ipt);
161
162 // Assign values of local coordinate in FaceElement
163 for (unsigned i = 0; i < my_dim; i++)
164 s[i] = this->integral_pt()->knot(ipt, i);
165
166 // Find corresponding coordinate in the the bulk element
167 s_bulk = this->local_coordinate_in_bulk(s);
168
169 /// Get outer unit normal
170 this->outer_unit_normal(ipt, unit_normal);
171
172 // Get velocity in bulk element
173 bulk_el_pt->interpolated_u_nst(s_bulk, veloc);
174
175 // Get normal component of veloc
176 double flux = 0.0;
177 for (unsigned i = 0; i < my_dim + 1; i++)
178 {
179 flux += veloc[i] * unit_normal[i];
180 }
181
182 // Modify bc: If outflow (flux>0) apply Neumann condition instead
183 if (flux > 0.0) flux = 0.0;
184
185 // Get pressure
186 double interpolated_press = bulk_el_pt->interpolated_p_nst(s_bulk);
187
188 // Call the pressure shape and test functions in bulk element
189 bulk_el_pt->pshape_nst(s_bulk, psip, testp);
190
191 // Find the Jacobian of the mapping within the FaceElement
192 double J = this->J_eulerian(s);
193
194 // Premultiply the weights and the Jacobian
195 double W = w * J;
196
197
198 // Number of master nodes and storage for the weight of the shape function
199 unsigned n_master = 1;
200 double hang_weight = 1.0;
201
202 // Loop over the pressure shape functions
203 for (unsigned l = 0; l < n_pres; l++)
204 {
205 // If the pressure dof is hanging
206 if (pressure_dof_is_hanging[l])
207 {
208 // Pressure dof is hanging so it must be nodal-based
209 // Get the hang info object
210 hang_info_pt = bulk_el_pt->pressure_node_pt(l)->hanging_pt(p_index);
211
212 // Get the number of master nodes from the pressure node
213 n_master = hang_info_pt->nmaster();
214 }
215 // Otherwise the node is its own master
216 else
217 {
218 n_master = 1;
219 }
220
221 // Loop over the master nodes
222 for (unsigned m = 0; m < n_master; m++)
223 {
224 // Get the number of the unknown
225 // If the pressure dof is hanging
226 if (pressure_dof_is_hanging[l])
227 {
228 // Get the local equation from the master node
229 local_eqn = bulk_el_pt->local_hang_eqn(
230 hang_info_pt->master_node_pt(m), p_index);
231 // Get the weight for the node
232 hang_weight = hang_info_pt->master_weight(m);
233 }
234 else
235 {
236 local_eqn = bulk_el_pt->p_local_eqn(l);
237 hang_weight = 1.0;
238 }
239
240 // If the equation is not pinned
241 if (local_eqn >= 0)
242 {
243 residuals[local_eqn] -=
244 re * flux * interpolated_press * testp[l] * W * hang_weight;
245
246 // Jacobian too?
247 if (flag)
248 {
249 // Number of master nodes and weights
250 unsigned n_master2 = 1;
251 double hang_weight2 = 1.0;
252
253 // Loop over the pressure shape functions
254 for (unsigned l2 = 0; l2 < n_pres; l2++)
255 {
256 // If the pressure dof is hanging
257 if (pressure_dof_is_hanging[l2])
258 {
259 hang_info2_pt =
260 bulk_el_pt->pressure_node_pt(l2)->hanging_pt(p_index);
261 // Pressure dof is hanging so it must be nodal-based
262 // Get the number of master nodes from the pressure node
263 n_master2 = hang_info2_pt->nmaster();
264 }
265 // Otherwise the node is its own master
266 else
267 {
268 n_master2 = 1;
269 }
270
271 // Loop over the master nodes
272 for (unsigned m2 = 0; m2 < n_master2; m2++)
273 {
274 // Get the number of the unknown
275 // If the pressure dof is hanging
276 if (pressure_dof_is_hanging[l2])
277 {
278 // Get the unknown from the master node
279 local_unknown = bulk_el_pt->local_hang_eqn(
280 hang_info2_pt->master_node_pt(m2), p_index);
281 // Get the weight from the hanging object
282 hang_weight2 = hang_info2_pt->master_weight(m2);
283 }
284 else
285 {
286 local_unknown = bulk_el_pt->p_local_eqn(l2);
287 hang_weight2 = 1.0;
288 }
289
290 // If the unknown is not pinned
291 if (local_unknown >= 0)
292 {
293 jacobian(local_eqn, local_unknown) -=
294 re * flux * psip[l2] * testp[l] * W * hang_weight *
295 hang_weight2;
296 }
297 }
298 }
299 }
300 }
301 } /*End of Jacobian calculation*/
302 } // End of if not boundary condition
303 } // End of loop over l
304 }
305
306
307 /// ////////////////////////////////////////////////////////////////////
308 /// ////////////////////////////////////////////////////////////////////
309 /// ////////////////////////////////////////////////////////////////////
310
311
312 //======================================================================
313 /// Refineable version of the Navier--Stokes equations
314 ///
315 ///
316 //======================================================================
317 template<unsigned DIM>
319 : public virtual NavierStokesEquations<DIM>,
320 public virtual RefineableElement,
321 public virtual ElementWithZ2ErrorEstimator
322 {
323 protected:
324 /// Unpin all pressure dofs in the element
326
327 /// Pin unused nodal pressure dofs (empty by default, because
328 /// by default pressure dofs are not associated with nodes)
330
331 public:
332 /// Constructor
334 : NavierStokesEquations<DIM>(),
337 {
338 }
339
340
341 /// Loop over all elements in Vector (which typically contains
342 /// all the elements in a fluid mesh) and pin the nodal pressure degrees
343 /// of freedom that are not being used. Function uses
344 /// the member function
345 /// - \c RefineableNavierStokesEquations::
346 /// pin_elemental_redundant_nodal_pressure_dofs()
347 /// .
348 /// which is empty by default and should be implemented for
349 /// elements with nodal pressure degrees of freedom
350 /// (e.g. for refineable Taylor-Hood.)
352 const Vector<GeneralisedElement*>& element_pt)
353 {
354 // Loop over all elements and call the function that pins their
355 // unused nodal pressure data
356 unsigned n_element = element_pt.size();
357 for (unsigned e = 0; e < n_element; e++)
358 {
359 dynamic_cast<RefineableNavierStokesEquations<DIM>*>(element_pt[e])
361 }
362 }
363
364 /// Unpin all pressure dofs in elements listed in vector.
366 const Vector<GeneralisedElement*>& element_pt)
367 {
368 // Loop over all elements
369 unsigned n_element = element_pt.size();
370 for (unsigned e = 0; e < n_element; e++)
371 {
372 dynamic_cast<RefineableNavierStokesEquations<DIM>*>(element_pt[e])
374 }
375 }
376
377 /// Pointer to n_p-th pressure node (Default: NULL,
378 /// indicating that pressure is not based on nodal interpolation).
379 virtual Node* pressure_node_pt(const unsigned& n_p)
380 {
381 return NULL;
382 }
383
384 /// Compute the diagonal of the velocity/pressure mass matrices.
385 /// If which one=0, both are computed, otherwise only the pressure
386 /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
387 /// LSC version of the preconditioner only needs that one)
389 Vector<double>& press_mass_diag,
390 Vector<double>& veloc_mass_diag,
391 const unsigned& which_one = 0);
392
393
394 /// Number of 'flux' terms for Z2 error estimation
396 {
397 // DIM diagonal strain rates, DIM(DIM -1) /2 off diagonal rates
398 return DIM + (DIM * (DIM - 1)) / 2;
399 }
400
401 /// Get 'flux' for Z2 error recovery: Upper triangular entries
402 /// in strain rate tensor.
404 {
405#ifdef PARANOID
406 unsigned num_entries = DIM + (DIM * (DIM - 1)) / 2;
407 if (flux.size() < num_entries)
408 {
409 std::ostringstream error_message;
410 error_message << "The flux vector has the wrong number of entries, "
411 << flux.size() << ", whereas it should be at least "
412 << num_entries << std::endl;
413 throw OomphLibError(error_message.str(),
414 OOMPH_CURRENT_FUNCTION,
415 OOMPH_EXCEPTION_LOCATION);
416 }
417#endif
418
419 // Get strain rate matrix
420 DenseMatrix<double> strainrate(DIM);
421 this->strain_rate(s, strainrate);
422
423 // Pack into flux Vector
424 unsigned icount = 0;
425
426 // Start with diagonal terms
427 for (unsigned i = 0; i < DIM; i++)
428 {
429 flux[icount] = strainrate(i, i);
430 icount++;
431 }
432
433 // Off diagonals row by row
434 for (unsigned i = 0; i < DIM; i++)
435 {
436 for (unsigned j = i + 1; j < DIM; j++)
437 {
438 flux[icount] = strainrate(i, j);
439 icount++;
440 }
441 }
442 }
443
444 /// Further build, pass the pointers down to the sons
446 {
447 // Find the father element
448 RefineableNavierStokesEquations<DIM>* cast_father_element_pt =
450 this->father_element_pt());
451
452 // Set the viscosity ratio pointer
453 this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
454 // Set the density ratio pointer
455 this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
456 // Set pointer to global Reynolds number
457 this->Re_pt = cast_father_element_pt->re_pt();
458 // Set pointer to global Reynolds number x Strouhal number (=Womersley)
459 this->ReSt_pt = cast_father_element_pt->re_st_pt();
460 // Set pointer to global Reynolds number x inverse Froude number
461 this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
462 // Set pointer to global gravity Vector
463 this->G_pt = cast_father_element_pt->g_pt();
464
465 // Set pointer to body force function
466 this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
467
468 // Set pointer to volumetric source function
469 this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
470
471 // Set the ALE flag
472 this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
473 }
474
475
476 /// Compute the derivatives of the i-th component of
477 /// velocity at point s with respect
478 /// to all data that can affect its value. In addition, return the global
479 /// equation numbers corresponding to the data.
480 /// Overload the non-refineable version to take account of hanging node
481 /// information
483 const unsigned& i,
484 Vector<double>& du_ddata,
485 Vector<unsigned>& global_eqn_number)
486 {
487 // Find number of nodes
488 unsigned n_node = this->nnode();
489 // Local shape function
490 Shape psi(n_node);
491 // Find values of shape function at the given local coordinate
492 this->shape(s, psi);
493
494 // Find the index at which the velocity component is stored
495 const unsigned u_nodal_index = this->u_index_nst(i);
496
497 // Storage for hang info pointer
498 HangInfo* hang_info_pt = 0;
499 // Storage for global equation
500 int global_eqn = 0;
501
502 // Find the number of dofs associated with interpolated u
503 unsigned n_u_dof = 0;
504 for (unsigned l = 0; l < n_node; l++)
505 {
506 unsigned n_master = 1;
507
508 // Local bool (is the node hanging)
509 bool is_node_hanging = this->node_pt(l)->is_hanging();
510
511 // If the node is hanging, get the number of master nodes
512 if (is_node_hanging)
513 {
514 hang_info_pt = this->node_pt(l)->hanging_pt();
515 n_master = hang_info_pt->nmaster();
516 }
517 // Otherwise there is just one master node, the node itself
518 else
519 {
520 n_master = 1;
521 }
522
523 // Loop over the master nodes
524 for (unsigned m = 0; m < n_master; m++)
525 {
526 // Get the equation number
527 if (is_node_hanging)
528 {
529 // Get the equation number from the master node
530 global_eqn =
531 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
532 }
533 else
534 {
535 // Global equation number
536 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
537 }
538
539 // If it's positive add to the count
540 if (global_eqn >= 0)
541 {
542 ++n_u_dof;
543 }
544 }
545 }
546
547 // Now resize the storage schemes
548 du_ddata.resize(n_u_dof, 0.0);
549 global_eqn_number.resize(n_u_dof, 0);
550
551 // Loop over th nodes again and set the derivatives
552 unsigned count = 0;
553 // Loop over the local nodes and sum
554 for (unsigned l = 0; l < n_node; l++)
555 {
556 unsigned n_master = 1;
557 double hang_weight = 1.0;
558
559 // Local bool (is the node hanging)
560 bool is_node_hanging = this->node_pt(l)->is_hanging();
561
562 // If the node is hanging, get the number of master nodes
563 if (is_node_hanging)
564 {
565 hang_info_pt = this->node_pt(l)->hanging_pt();
566 n_master = hang_info_pt->nmaster();
567 }
568 // Otherwise there is just one master node, the node itself
569 else
570 {
571 n_master = 1;
572 }
573
574 // Loop over the master nodes
575 for (unsigned m = 0; m < n_master; m++)
576 {
577 // If the node is hanging get weight from master node
578 if (is_node_hanging)
579 {
580 // Get the hang weight from the master node
581 hang_weight = hang_info_pt->master_weight(m);
582 }
583 else
584 {
585 // Node contributes with full weight
586 hang_weight = 1.0;
587 }
588
589 // Get the equation number
590 if (is_node_hanging)
591 {
592 // Get the equation number from the master node
593 global_eqn =
594 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
595 }
596 else
597 {
598 // Local equation number
599 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
600 }
601
602 if (global_eqn >= 0)
603 {
604 // Set the global equation number
605 global_eqn_number[count] = global_eqn;
606 // Set the derivative with respect to the unknown
607 du_ddata[count] = psi[l] * hang_weight;
608 // Increase the counter
609 ++count;
610 }
611 }
612 }
613 }
614
615
616 protected:
617 /// Add element's contribution to elemental residual vector and/or
618 /// Jacobian matrix
619 /// flag=1: compute both
620 /// flag=0: compute only residual vector
622 Vector<double>& residuals,
623 DenseMatrix<double>& jacobian,
624 DenseMatrix<double>& mass_matrix,
625 unsigned flag);
626
627 /// Compute the residuals for the associated pressure advection
628 /// diffusion problem. Used by the Fp preconditioner.
629 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
631 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag);
632
633
634 /// Compute derivatives of elemental residual vector with respect
635 /// to nodal coordinates. Overwrites default implementation in
636 /// FiniteElement base class.
637 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
639 RankThreeTensor<double>& dresidual_dnodal_coordinates);
640 };
641
642
643 //======================================================================
644 /// Refineable version of Taylor Hood elements. These classes
645 /// can be written in total generality.
646 //======================================================================
647 template<unsigned DIM>
649 : public QTaylorHoodElement<DIM>,
650 public virtual RefineableNavierStokesEquations<DIM>,
651 public virtual RefineableQElement<DIM>
652 {
653 private:
654 /// Unpin all pressure dofs
656 {
657 // find the index at which the pressure is stored
658 int p_index = this->p_nodal_index_nst();
659 unsigned n_node = this->nnode();
660 // loop over nodes
661 for (unsigned n = 0; n < n_node; n++)
662 {
663 this->node_pt(n)->unpin(p_index);
664 }
665 }
666
667 /// Pin all nodal pressure dofs that are not required
669 {
670 // Find the pressure index
671 int p_index = this->p_nodal_index_nst();
672 // Loop over all nodes
673 unsigned n_node = this->nnode();
674 // loop over all nodes and pin all the nodal pressures
675 for (unsigned n = 0; n < n_node; n++)
676 {
677 this->node_pt(n)->pin(p_index);
678 }
679
680 // Loop over all actual pressure nodes and unpin if they're not hanging
681 unsigned n_pres = this->npres_nst();
682 for (unsigned l = 0; l < n_pres; l++)
683 {
684 Node* nod_pt = this->pressure_node_pt(l);
685 if (!nod_pt->is_hanging(p_index))
686 {
687 nod_pt->unpin(p_index);
688 }
689 }
690 }
691
692 public:
693 /// Constructor
697 RefineableQElement<DIM>(),
698 QTaylorHoodElement<DIM>()
699 {
700 }
701
702 /// Number of values required at local node n. In order to simplify
703 /// matters, we allocate storage for pressure variables at all the nodes
704 /// and then pin those that are not used.
705 unsigned required_nvalue(const unsigned& n) const
706 {
707 return DIM + 1;
708 }
709
710 /// Number of continuously interpolated values: (DIM velocities + 1
711 /// pressure)
713 {
714 return DIM + 1;
715 }
716
717 /// Rebuild from sons: empty
718 void rebuild_from_sons(Mesh*& mesh_pt) {}
719
720 /// Order of recovery shape functions for Z2 error estimation:
721 /// Same order as shape functions.
723 {
724 return 2;
725 }
726
727 /// Number of vertex nodes in the element
728 unsigned nvertex_node() const
729 {
731 }
732
733 /// Pointer to the j-th vertex node in the element
734 Node* vertex_node_pt(const unsigned& j) const
735 {
737 }
738
739 /// Get the function value u in Vector.
740 /// Note: Given the generality of the interface (this function
741 /// is usually called from black-box documentation or interpolation
742 /// routines), the values Vector sets its own size in here.
744 Vector<double>& values)
745 {
746 // Set size of Vector: u,v,p and initialise to zero
747 values.resize(DIM + 1, 0.0);
748
749 // Calculate velocities: values[0],...
750 for (unsigned i = 0; i < DIM; i++)
751 {
752 values[i] = this->interpolated_u_nst(s, i);
753 }
754
755 // Calculate pressure: values[DIM]
756 values[DIM] = this->interpolated_p_nst(s);
757 }
758
759 /// Get the function value u in Vector.
760 /// Note: Given the generality of the interface (this function
761 /// is usually called from black-box documentation or interpolation
762 /// routines), the values Vector sets its own size in here.
763 void get_interpolated_values(const unsigned& t,
764 const Vector<double>& s,
765 Vector<double>& values)
766 {
767 // Set size of Vector: u,v,p
768 values.resize(DIM + 1);
769
770 // Initialise
771 for (unsigned i = 0; i < DIM + 1; i++)
772 {
773 values[i] = 0.0;
774 }
775
776 // Find out how many nodes there are
777 unsigned n_node = this->nnode();
778
779 // Shape functions
780 Shape psif(n_node);
781 this->shape(s, psif);
782
783 // Calculate velocities: values[0],...
784 for (unsigned i = 0; i < DIM; i++)
785 {
786 // Get the index at which the i-th velocity is stored
787 unsigned u_nodal_index = this->u_index_nst(i);
788 for (unsigned l = 0; l < n_node; l++)
789 {
790 values[i] += this->nodal_value(t, l, u_nodal_index) * psif[l];
791 }
792 }
793
794 // Calculate pressure: values[DIM]
795 //(no history is carried in the pressure)
796 values[DIM] = this->interpolated_p_nst(s);
797 }
798
799
800 /// Perform additional hanging node procedures for variables
801 /// that are not interpolated by all nodes. The pressures are stored
802 /// at the p_nodal_index_nst-th location in each node
804 {
805 this->setup_hang_for_value(this->p_nodal_index_nst());
806 }
807
808 /// Pointer to n_p-th pressure node
809 Node* pressure_node_pt(const unsigned& n_p)
810 {
811 return this->node_pt(this->Pconv[n_p]);
812 }
813
814 /// The velocities are isoparametric and so the "nodes" interpolating
815 /// the velocities are the geometric nodes. The pressure "nodes" are a
816 /// subset of the nodes, so when value_id==DIM, the n-th pressure
817 /// node is returned.
818 Node* interpolating_node_pt(const unsigned& n, const int& value_id)
819
820 {
821 // The only different nodes are the pressure nodes
822 if (value_id == DIM)
823 {
824 return this->pressure_node_pt(n);
825 }
826 // The other variables are interpolated via the usual nodes
827 else
828 {
829 return this->node_pt(n);
830 }
831 }
832
833 /// The pressure nodes are the corner nodes, so when n_value==DIM,
834 /// the fraction is the same as the 1d node number, 0 or 1.
836 const unsigned& i,
837 const int& value_id)
838 {
839 if (value_id == DIM)
840 {
841 // The pressure nodes are just located on the boundaries at 0 or 1
842 return double(n1d);
843 }
844 // Otherwise the velocity nodes are the same as the geometric ones
845 else
846 {
847 return this->local_one_d_fraction_of_node(n1d, i);
848 }
849 }
850
851 /// The velocity nodes are the same as the geometric nodes. The
852 /// pressure nodes must be calculated by using the same methods as
853 /// the geometric nodes, but by recalling that there are only two pressure
854 /// nodes per edge.
856 const int& value_id)
857 {
858 // If we are calculating pressure nodes
859 if (value_id == DIM)
860 {
861 // Storage for the index of the pressure node
862 unsigned total_index = 0;
863 // The number of nodes along each 1d edge is 2.
864 unsigned NNODE_1D = 2;
865 // Storage for the index along each boundary
866 Vector<int> index(DIM);
867 // Loop over the coordinates
868 for (unsigned i = 0; i < DIM; i++)
869 {
870 // If we are at the lower limit, the index is zero
871 if (s[i] == -1.0)
872 {
873 index[i] = 0;
874 }
875 // If we are at the upper limit, the index is the number of nodes
876 // minus 1
877 else if (s[i] == 1.0)
878 {
879 index[i] = NNODE_1D - 1;
880 }
881 // Otherwise, we have to calculate the index in general
882 else
883 {
884 // For uniformly spaced nodes the 0th node number would be
885 double float_index = 0.5 * (1.0 + s[i]) * (NNODE_1D - 1);
886 index[i] = int(float_index);
887 // What is the excess. This should be safe because the
888 // taking the integer part rounds down
889 double excess = float_index - index[i];
890 // If the excess is bigger than our tolerance there is no node,
891 // return null
894 {
895 return 0;
896 }
897 }
898 /// Construct the general pressure index from the components.
899 total_index +=
900 index[i] * static_cast<unsigned>(pow(static_cast<float>(NNODE_1D),
901 static_cast<int>(i)));
902 }
903 // If we've got here we have a node, so let's return a pointer to it
904 return this->pressure_node_pt(total_index);
905 }
906 // Otherwise velocity nodes are the same as pressure nodes
907 else
908 {
909 return this->get_node_at_local_coordinate(s);
910 }
911 }
912
913
914 /// The number of 1d pressure nodes is 2, the number of 1d velocity
915 /// nodes is the same as the number of 1d geometric nodes.
916 unsigned ninterpolating_node_1d(const int& value_id)
917 {
918 if (value_id == DIM)
919 {
920 return 2;
921 }
922 else
923 {
924 return this->nnode_1d();
925 }
926 }
927
928 /// The number of pressure nodes is 2^DIM. The number of
929 /// velocity nodes is the same as the number of geometric nodes.
930 unsigned ninterpolating_node(const int& value_id)
931 {
932 if (value_id == DIM)
933 {
934 return static_cast<unsigned>(pow(2.0, static_cast<int>(DIM)));
935 }
936 else
937 {
938 return this->nnode();
939 }
940 }
941
942 /// The basis interpolating the pressure is given by pshape().
943 /// / The basis interpolating the velocity is shape().
945 Shape& psi,
946 const int& value_id) const
947 {
948 if (value_id == DIM)
949 {
950 return this->pshape_nst(s, psi);
951 }
952 else
953 {
954 return this->shape(s, psi);
955 }
956 }
957
958
959 /// Build FaceElements that apply the Robin boundary condition
960 /// to the pressure advection diffusion problem required by
961 /// Fp preconditioner
962 void build_fp_press_adv_diff_robin_bc_element(const unsigned& face_index)
963 {
966 RefineableQTaylorHoodElement<DIM>>(this, face_index));
967 }
968
969
970 /// Add to the set \c paired_load_data pairs containing
971 /// - the pointer to a Data object
972 /// and
973 /// - the index of the value in that Data object
974 /// .
975 /// for all values (pressures, velocities) that affect the
976 /// load computed in the \c get_load(...) function.
977 /// (Overloads non-refineable version and takes hanging nodes
978 /// into account)
980 std::set<std::pair<Data*, unsigned>>& paired_load_data)
981 {
982 // Get the nodal indices at which the velocities are stored
983 unsigned u_index[DIM];
984 for (unsigned i = 0; i < DIM; i++)
985 {
986 u_index[i] = this->u_index_nst(i);
987 }
988
989 // Loop over the nodes
990 unsigned n_node = this->nnode();
991 for (unsigned n = 0; n < n_node; n++)
992 {
993 // Pointer to current node
994 Node* nod_pt = this->node_pt(n);
995
996 // Check if it's hanging:
997 if (nod_pt->is_hanging())
998 {
999 // It's hanging -- get number of master nodes
1000 unsigned nmaster = nod_pt->hanging_pt()->nmaster();
1001
1002 // Loop over masters
1003 for (unsigned j = 0; j < nmaster; j++)
1004 {
1005 Node* master_nod_pt = nod_pt->hanging_pt()->master_node_pt(j);
1006
1007 // Loop over the velocity components and add pointer to their data
1008 // and indices to the vectors
1009 for (unsigned i = 0; i < DIM; i++)
1010 {
1011 paired_load_data.insert(
1012 std::make_pair(master_nod_pt, u_index[i]));
1013 }
1014 }
1015 }
1016 // Not hanging
1017 else
1018 {
1019 // Loop over the velocity components and add pointer to their data
1020 // and indices to the vectors
1021 for (unsigned i = 0; i < DIM; i++)
1022 {
1023 paired_load_data.insert(
1024 std::make_pair(this->node_pt(n), u_index[i]));
1025 }
1026 }
1027 }
1028
1029 // Get the nodal index at which the pressure is stored
1030 int p_index = this->p_nodal_index_nst();
1031
1032 // Loop over the pressure data
1033 unsigned n_pres = this->npres_nst();
1034 for (unsigned l = 0; l < n_pres; l++)
1035 {
1036 // Get the pointer to the nodal pressure
1037 Node* pres_node_pt = this->pressure_node_pt(l);
1038 // Check if the pressure dof is hanging
1039 if (pres_node_pt->is_hanging(p_index))
1040 {
1041 // Get the pointer to the hang info object
1042 // (pressure is stored as p_index--th nodal dof).
1043 HangInfo* hang_info_pt = pres_node_pt->hanging_pt(p_index);
1044
1045 // Get number of pressure master nodes (pressure is stored
1046 unsigned nmaster = hang_info_pt->nmaster();
1047
1048 // Loop over pressure master nodes
1049 for (unsigned m = 0; m < nmaster; m++)
1050 {
1051 // The p_index-th entry in each nodal data is the pressure, which
1052 // affects the traction
1053 paired_load_data.insert(
1054 std::make_pair(hang_info_pt->master_node_pt(m), p_index));
1055 }
1056 }
1057 // It's not hanging
1058 else
1059 {
1060 // The p_index-th entry in each nodal data is the pressure, which
1061 // affects the traction
1062 paired_load_data.insert(std::make_pair(pres_node_pt, p_index));
1063 }
1064 }
1065 }
1066 };
1067
1068
1069 //=======================================================================
1070 /// Face geometry of the RefineableQTaylorHoodElements is the
1071 /// same as the Face geometry of the QTaylorHoodElements.
1072 //=======================================================================
1073 template<unsigned DIM>
1075 : public virtual FaceGeometry<QTaylorHoodElement<DIM>>
1076 {
1077 public:
1079 };
1080
1081
1082 //=======================================================================
1083 /// Face geometry of the face geometry of
1084 /// the RefineableQTaylorHoodElements is the
1085 /// same as the Face geometry of the Face geometry of QTaylorHoodElements.
1086 //=======================================================================
1087 template<unsigned DIM>
1089 : public virtual FaceGeometry<FaceGeometry<QTaylorHoodElement<DIM>>>
1090 {
1091 public:
1093 };
1094
1095
1096 /// ////////////////////////////////////////////////////////////////////////
1097 /// ////////////////////////////////////////////////////////////////////////
1098 /// ////////////////////////////////////////////////////////////////////////
1099
1100
1101 //======================================================================
1102 /// Refineable version of Crouzeix Raviart elements. Generic class definitions
1103 //======================================================================
1104 template<unsigned DIM>
1106 : public QCrouzeixRaviartElement<DIM>,
1107 public virtual RefineableNavierStokesEquations<DIM>,
1108 public virtual RefineableQElement<DIM>
1109 {
1110 private:
1111 /// Unpin all internal pressure dofs
1113 {
1114 unsigned n_pres = this->npres_nst();
1115 // loop over pressure dofs and unpin them
1116 for (unsigned l = 0; l < n_pres; l++)
1117 {
1119 }
1120 }
1121
1122 public:
1123 /// Constructor
1127 RefineableQElement<DIM>(),
1129 {
1130 }
1131
1132
1133 /// Broken copy constructor
1135 const RefineableQCrouzeixRaviartElement<DIM>& dummy) = delete;
1136
1137 /// Broken assignment operator
1138 // Commented out broken assignment operator because this can lead to a
1139 // conflict warning when used in the virtual inheritence hierarchy.
1140 // Essentially the compiler doesn't realise that two separate
1141 // implementations of the broken function are the same and so, quite
1142 // rightly, it shouts.
1143 /*void operator=(const RefineableQCrouzeixRaviartElement<DIM>&) = delete;*/
1144
1145 /// Number of continuously interpolated values: DIM (velocities)
1147 {
1148 return DIM;
1149 }
1150
1151 /// Rebuild from sons: Reconstruct pressure from the (merged) sons
1152 /// This must be specialised for each dimension.
1153 inline void rebuild_from_sons(Mesh*& mesh_pt);
1154
1155 /// Order of recovery shape functions for Z2 error estimation:
1156 /// Same order as shape functions.
1158 {
1159 return 2;
1160 }
1161
1162 /// Number of vertex nodes in the element
1163 unsigned nvertex_node() const
1164 {
1166 }
1167
1168 /// Pointer to the j-th vertex node in the element
1169 Node* vertex_node_pt(const unsigned& j) const
1170 {
1172 }
1173
1174 /// Get the function value u in Vector.
1175 /// Note: Given the generality of the interface (this function
1176 /// is usually called from black-box documentation or interpolation
1177 /// routines), the values Vector sets its own size in here.
1179 Vector<double>& values)
1180 {
1181 // Set size of Vector: u,v,p and initialise to zero
1182 values.resize(DIM, 0.0);
1183
1184 // Calculate velocities: values[0],...
1185 for (unsigned i = 0; i < DIM; i++)
1186 {
1187 values[i] = this->interpolated_u_nst(s, i);
1188 }
1189 }
1190
1191 /// Get all function values [u,v..,p] at previous timestep t
1192 /// (t=0: present; t>0: previous timestep).
1193 ///
1194 /// Note: Given the generality of the interface (this function
1195 /// is usually called from black-box documentation or interpolation
1196 /// routines), the values Vector sets its own size in here.
1197 ///
1198 /// Note: No pressure history is kept, so pressure is always
1199 /// the current value.
1200 void get_interpolated_values(const unsigned& t,
1201 const Vector<double>& s,
1202 Vector<double>& values)
1203 {
1204 // Set size of Vector: u,v,p
1205 values.resize(DIM);
1206
1207 // Initialise
1208 for (unsigned i = 0; i < DIM; i++)
1209 {
1210 values[i] = 0.0;
1211 }
1212
1213 // Find out how many nodes there are
1214 unsigned n_node = this->nnode();
1215
1216 // Shape functions
1217 Shape psif(n_node);
1218 this->shape(s, psif);
1219
1220 // Calculate velocities: values[0],...
1221 for (unsigned i = 0; i < DIM; i++)
1222 {
1223 // Get the nodal index at which the i-th velocity component is stored
1224 unsigned u_nodal_index = this->u_index_nst(i);
1225 for (unsigned l = 0; l < n_node; l++)
1226 {
1227 values[i] += this->nodal_value(t, l, u_nodal_index) * psif[l];
1228 }
1229 }
1230 }
1231
1232 /// Perform additional hanging node procedures for variables
1233 /// that are not interpolated by all nodes. Empty
1235
1236 /// Further build for Crouzeix_Raviart interpolates the internal
1237 /// pressure dofs from father element: Make sure pressure values and
1238 /// dp/ds agree between fathers and sons at the midpoints of the son
1239 /// elements. This must be specialised for each dimension.
1240 inline void further_build();
1241
1242
1243 /// Build FaceElements that apply the Robin boundary condition
1244 /// to the pressure advection diffusion problem required by
1245 /// Fp preconditioner
1246 void build_fp_press_adv_diff_robin_bc_element(const unsigned& face_index)
1247 {
1250 RefineableQCrouzeixRaviartElement<DIM>>(this, face_index));
1251 }
1252
1253
1254 /// Add to the set \c paired_load_data pairs containing
1255 /// - the pointer to a Data object
1256 /// and
1257 /// - the index of the value in that Data object
1258 /// .
1259 /// for all values (pressures, velocities) that affect the
1260 /// load computed in the \c get_load(...) function.
1261 /// (Overloads non-refineable version and takes hanging nodes
1262 /// into account)
1264 std::set<std::pair<Data*, unsigned>>& paired_load_data)
1265 {
1266 // Get the nodal indices at which the velocities are stored
1267 unsigned u_index[DIM];
1268 for (unsigned i = 0; i < DIM; i++)
1269 {
1270 u_index[i] = this->u_index_nst(i);
1271 }
1272
1273 // Loop over the nodes
1274 unsigned n_node = this->nnode();
1275 for (unsigned n = 0; n < n_node; n++)
1276 {
1277 // Pointer to current node
1278 Node* nod_pt = this->node_pt(n);
1279
1280 // Check if it's hanging:
1281 if (nod_pt->is_hanging())
1282 {
1283 // It's hanging -- get number of master nodes
1284 unsigned nmaster = nod_pt->hanging_pt()->nmaster();
1285
1286 // Loop over masters
1287 for (unsigned j = 0; j < nmaster; j++)
1288 {
1289 Node* master_nod_pt = nod_pt->hanging_pt()->master_node_pt(j);
1290
1291 // Loop over the velocity components and add pointer to their data
1292 // and indices to the vectors
1293 for (unsigned i = 0; i < DIM; i++)
1294 {
1295 paired_load_data.insert(
1296 std::make_pair(master_nod_pt, u_index[i]));
1297 }
1298 }
1299 }
1300 // Not hanging
1301 else
1302 {
1303 // Loop over the velocity components and add pointer to their data
1304 // and indices to the vectors
1305 for (unsigned i = 0; i < DIM; i++)
1306 {
1307 paired_load_data.insert(
1308 std::make_pair(this->node_pt(n), u_index[i]));
1309 }
1310 }
1311 }
1312
1313
1314 // Loop over the pressure data (can't be hanging!)
1315 unsigned n_pres = this->npres_nst();
1316 for (unsigned l = 0; l < n_pres; l++)
1317 {
1318 // The entries in the internal data at P_nst_internal_index
1319 // are the pressures, which affect the traction
1320 paired_load_data.insert(std::make_pair(
1321 this->internal_data_pt(this->P_nst_internal_index), l));
1322 }
1323 }
1324 };
1325
1326
1327 //======================================================================
1328 /// p-refineable version of Crouzeix Raviart elements. Generic class
1329 /// definitions
1330 //======================================================================
1331 template<unsigned DIM>
1333 : public QCrouzeixRaviartElement<DIM>,
1334 public virtual RefineableNavierStokesEquations<DIM>,
1335 public virtual PRefineableQElement<DIM, 3>
1336 {
1337 private:
1338 /// Unpin all internal pressure dofs
1340 {
1341 unsigned n_pres = this->npres_nst();
1342 n_pres = this->internal_data_pt(this->P_nst_internal_index)->nvalue();
1343 // loop over pressure dofs and unpin them
1344 for (unsigned l = 0; l < n_pres; l++)
1345 {
1347 }
1348 }
1349
1350 public:
1351 /// Constructor
1355 PRefineableQElement<DIM, 3>(),
1357 {
1358 // Set the p-order
1359 this->p_order() = 3;
1360
1361 // Set integration scheme
1362 // (To avoid memory leaks in pre-build and p-refine where new
1363 // integration schemes are created)
1365
1366 // Resize pressure storage
1367 // (Constructor for QCrouzeixRaviartElement sets up DIM+1 pressure values)
1368 if (this->internal_data_pt(this->P_nst_internal_index)->nvalue() <=
1369 this->npres_nst())
1370 {
1372 ->resize(this->npres_nst());
1373 }
1374 else
1375 {
1376 Data* new_data_pt = new Data(this->npres_nst());
1377 delete this->internal_data_pt(this->P_nst_internal_index);
1378 this->internal_data_pt(this->P_nst_internal_index) = new_data_pt;
1379 }
1380 }
1381
1382 /// Destructor
1384 {
1385 delete this->integral_pt();
1386 }
1387
1388
1389 /// Broken copy constructor
1391 const PRefineableQCrouzeixRaviartElement<DIM>& dummy) = delete;
1392
1393 /// Broken assignment operator
1394 /*void operator=(const PRefineableQCrouzeixRaviartElement<DIM>&) =
1395 * delete;*/
1396
1397 /// Return the i-th pressure value
1398 /// (Discontinous pressure interpolation -- no need to cater for hanging
1399 /// nodes).
1400 double p_nst(const unsigned& i) const
1401 {
1402 return this->internal_data_pt(this->P_nst_internal_index)->value(i);
1403 }
1404
1405 double p_nst(const unsigned& t, const unsigned& i) const
1406 {
1407 return this->internal_data_pt(this->P_nst_internal_index)->value(t, i);
1408 }
1409
1410 /// // Return number of pressure values
1411 unsigned npres_nst() const
1412 {
1413 return (this->p_order() - 2) * (this->p_order() - 2);
1414 }
1415
1416 /// Pin p_dof-th pressure dof and set it to value specified by p_value.
1417 void fix_pressure(const unsigned& p_dof, const double& p_value)
1418 {
1419 this->internal_data_pt(this->P_nst_internal_index)->pin(p_dof);
1421 ->set_value(p_dof, p_value);
1422 }
1423
1424 unsigned required_nvalue(const unsigned& n) const
1425 {
1426 return DIM;
1427 }
1428
1429 /// Number of continuously interpolated values: DIM (velocities)
1431 {
1432 return DIM;
1433 }
1434
1435 /// Rebuild from sons: Reconstruct pressure from the (merged) sons
1436 /// This must be specialised for each dimension.
1437 void rebuild_from_sons(Mesh*& mesh_pt)
1438 {
1439 // Do p-refineable version
1441 // Do Crouzeix-Raviart version
1442 // Need to reconstruct pressure manually!
1443 for (unsigned p = 0; p < npres_nst(); p++)
1444 {
1445 // BENFLAG: Set to zero for now -- don't do projection problem yet
1446 this->internal_data_pt(this->P_nst_internal_index)->set_value(p, 0.0);
1447 }
1448 }
1449
1450 /// Order of recovery shape functions for Z2 error estimation:
1451 /// - Same order as shape functions.
1452 // unsigned nrecovery_order()
1453 // {
1454 // if(this->nnode_1d() < 4) {return (this->nnode_1d()-1);}
1455 // else {return 3;}
1456 // }
1457 /// - Constant recovery order, since recovery order of the first element
1458 /// is used for the whole mesh.
1460 {
1461 return 3;
1462 }
1463
1464 /// Number of vertex nodes in the element
1465 unsigned nvertex_node() const
1466 {
1468 }
1469
1470 /// Pointer to the j-th vertex node in the element
1471 Node* vertex_node_pt(const unsigned& j) const
1472 {
1474 }
1475
1476 /// Velocity shape and test functions and their derivs
1477 /// w.r.t. to global coords at local coordinate s (taken from geometry)
1478 /// Return Jacobian of mapping between local and global coordinates.
1480 Shape& psi,
1481 DShape& dpsidx,
1482 Shape& test,
1483 DShape& dtestdx) const;
1484
1485 /// Velocity shape and test functions and their derivs
1486 /// w.r.t. to global coords at ipt-th integation point (taken from geometry)
1487 /// Return Jacobian of mapping between local and global coordinates.
1488 inline double dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
1489 Shape& psi,
1490 DShape& dpsidx,
1491 Shape& test,
1492 DShape& dtestdx) const;
1493
1494 /// Pressure shape functions at local coordinate s
1495 inline void pshape_nst(const Vector<double>& s, Shape& psi) const;
1496
1497 /// Pressure shape and test functions at local coordinte s
1498 inline void pshape_nst(const Vector<double>& s,
1499 Shape& psi,
1500 Shape& test) const;
1501
1502 /// Get the function value u in Vector.
1503 /// Note: Given the generality of the interface (this function
1504 /// is usually called from black-box documentation or interpolation
1505 /// routines), the values Vector sets its own size in here.
1507 Vector<double>& values)
1508 {
1509 // Set size of Vector: u,v,p and initialise to zero
1510 values.resize(DIM, 0.0);
1511
1512 // Calculate velocities: values[0],...
1513 for (unsigned i = 0; i < DIM; i++)
1514 {
1515 values[i] = this->interpolated_u_nst(s, i);
1516 }
1517 }
1518
1519 /// Get all function values [u,v..,p] at previous timestep t
1520 /// (t=0: present; t>0: previous timestep).
1521 ///
1522 /// Note: Given the generality of the interface (this function
1523 /// is usually called from black-box documentation or interpolation
1524 /// routines), the values Vector sets its own size in here.
1525 ///
1526 /// Note: No pressure history is kept, so pressure is always
1527 /// the current value.
1528 void get_interpolated_values(const unsigned& t,
1529 const Vector<double>& s,
1530 Vector<double>& values)
1531 {
1532 // Set size of Vector: u,v,p
1533 values.resize(DIM);
1534
1535 // Initialise
1536 for (unsigned i = 0; i < DIM; i++)
1537 {
1538 values[i] = 0.0;
1539 }
1540
1541 // Find out how many nodes there are
1542 unsigned n_node = this->nnode();
1543
1544 // Shape functions
1545 Shape psif(n_node);
1546 this->shape(s, psif);
1547
1548 // Calculate velocities: values[0],...
1549 for (unsigned i = 0; i < DIM; i++)
1550 {
1551 // Get the nodal index at which the i-th velocity component is stored
1552 unsigned u_nodal_index = this->u_index_nst(i);
1553 for (unsigned l = 0; l < n_node; l++)
1554 {
1555 values[i] += this->nodal_value(t, l, u_nodal_index) * psif[l];
1556 }
1557 }
1558 }
1559
1560 /// Perform additional hanging node procedures for variables
1561 /// that are not interpolated by all nodes. Empty
1563
1564 /// Further build for Crouzeix_Raviart interpolates the internal
1565 /// pressure dofs from father element: Make sure pressure values and
1566 /// dp/ds agree between fathers and sons at the midpoints of the son
1567 /// elements. This must be specialised for each dimension.
1569 };
1570
1571
1572 //=======================================================================
1573 /// Face geometry of the RefineableQuadQCrouzeixRaviartElements
1574 //=======================================================================
1575 template<unsigned DIM>
1577 : public virtual FaceGeometry<QCrouzeixRaviartElement<DIM>>
1578 {
1579 public:
1581 };
1582
1583 //======================================================================
1584 /// Face geometry of the face geometry of
1585 /// the RefineableQCrouzeixRaviartElements is the
1586 /// same as the Face geometry of the Face geometry of
1587 /// QCrouzeixRaviartElements.
1588 //=======================================================================
1589 template<unsigned DIM>
1591 : public virtual FaceGeometry<FaceGeometry<QCrouzeixRaviartElement<DIM>>>
1592 {
1593 public:
1595 {
1596 }
1597 };
1598
1599
1600 // Inline functions
1601
1602 //=====================================================================
1603 /// 2D Rebuild from sons: Reconstruct pressure from the (merged) sons
1604 //=====================================================================
1605 template<>
1607 Mesh*& mesh_pt)
1608 {
1609 using namespace QuadTreeNames;
1610
1611 // Central pressure value:
1612 //-----------------------
1613
1614 // Use average of the sons central pressure values
1615 // Other options: Take average of the four (discontinuous)
1616 // pressure values at the father's midpoint]
1617
1618 double av_press = 0.0;
1619
1620 // Loop over the sons
1621 for (unsigned ison = 0; ison < 4; ison++)
1622 {
1623 // Add the sons midnode pressure
1624 // Note that we can assume that the pressure is stored at the same
1625 // location because these are EXACTLY the same type of elements
1626 av_press += quadtree_pt()
1627 ->son_pt(ison)
1628 ->object_pt()
1629 ->internal_data_pt(this->P_nst_internal_index)
1630 ->value(0);
1631 }
1632
1633 // Use the average
1634 internal_data_pt(this->P_nst_internal_index)->set_value(0, 0.25 * av_press);
1635
1636
1637 // Slope in s_0 direction
1638 //----------------------
1639
1640 // Use average of the 2 FD approximations based on the
1641 // elements central pressure values
1642 // [Other options: Take average of the four
1643 // pressure derivatives]
1644
1645 double slope1 = quadtree_pt()
1646 ->son_pt(SE)
1647 ->object_pt()
1648 ->internal_data_pt(this->P_nst_internal_index)
1649 ->value(0) -
1650 quadtree_pt()
1651 ->son_pt(SW)
1652 ->object_pt()
1653 ->internal_data_pt(this->P_nst_internal_index)
1654 ->value(0);
1655
1656 double slope2 = quadtree_pt()
1657 ->son_pt(NE)
1658 ->object_pt()
1659 ->internal_data_pt(this->P_nst_internal_index)
1660 ->value(0) -
1661 quadtree_pt()
1662 ->son_pt(NW)
1663 ->object_pt()
1664 ->internal_data_pt(this->P_nst_internal_index)
1665 ->value(0);
1666
1667
1668 // Use the average
1669 internal_data_pt(this->P_nst_internal_index)
1670 ->set_value(1, 0.5 * (slope1 + slope2));
1671
1672
1673 // Slope in s_1 direction
1674 //----------------------
1675
1676 // Use average of the 2 FD approximations based on the
1677 // elements central pressure values
1678 // [Other options: Take average of the four
1679 // pressure derivatives]
1680
1681 slope1 = quadtree_pt()
1682 ->son_pt(NE)
1683 ->object_pt()
1684 ->internal_data_pt(this->P_nst_internal_index)
1685 ->value(0) -
1686 quadtree_pt()
1687 ->son_pt(SE)
1688 ->object_pt()
1689 ->internal_data_pt(this->P_nst_internal_index)
1690 ->value(0);
1691
1692 slope2 = quadtree_pt()
1693 ->son_pt(NW)
1694 ->object_pt()
1695 ->internal_data_pt(this->P_nst_internal_index)
1696 ->value(0) -
1697 quadtree_pt()
1698 ->son_pt(SW)
1699 ->object_pt()
1700 ->internal_data_pt(this->P_nst_internal_index)
1701 ->value(0);
1702
1703
1704 // Use the average
1705 internal_data_pt(this->P_nst_internal_index)
1706 ->set_value(2, 0.5 * (slope1 + slope2));
1707 }
1708
1709
1710 //=================================================================
1711 /// 3D Rebuild from sons: Reconstruct pressure from the (merged) sons
1712 //=================================================================
1713 template<>
1715 Mesh*& mesh_pt)
1716 {
1717 using namespace OcTreeNames;
1718
1719 // Central pressure value:
1720 //-----------------------
1721
1722 // Use average of the sons central pressure values
1723 // Other options: Take average of the four (discontinuous)
1724 // pressure values at the father's midpoint]
1725
1726 double av_press = 0.0;
1727
1728 // Loop over the sons
1729 for (unsigned ison = 0; ison < 8; ison++)
1730 {
1731 // Add the sons midnode pressure
1732 av_press += octree_pt()
1733 ->son_pt(ison)
1734 ->object_pt()
1735 ->internal_data_pt(this->P_nst_internal_index)
1736 ->value(0);
1737 }
1738
1739 // Use the average
1740 internal_data_pt(this->P_nst_internal_index)
1741 ->set_value(0, 0.125 * av_press);
1742
1743
1744 // Slope in s_0 direction
1745 //----------------------
1746
1747 // Use average of the 4 FD approximations based on the
1748 // elements central pressure values
1749 // [Other options: Take average of the four
1750 // pressure derivatives]
1751
1752 double slope1 = octree_pt()
1753 ->son_pt(RDF)
1754 ->object_pt()
1755 ->internal_data_pt(this->P_nst_internal_index)
1756 ->value(0) -
1757 octree_pt()
1758 ->son_pt(LDF)
1759 ->object_pt()
1760 ->internal_data_pt(this->P_nst_internal_index)
1761 ->value(0);
1762
1763 double slope2 = octree_pt()
1764 ->son_pt(RUF)
1765 ->object_pt()
1766 ->internal_data_pt(this->P_nst_internal_index)
1767 ->value(0) -
1768 octree_pt()
1769 ->son_pt(LUF)
1770 ->object_pt()
1771 ->internal_data_pt(this->P_nst_internal_index)
1772 ->value(0);
1773
1774 double slope3 = octree_pt()
1775 ->son_pt(RDB)
1776 ->object_pt()
1777 ->internal_data_pt(this->P_nst_internal_index)
1778 ->value(0) -
1779 octree_pt()
1780 ->son_pt(LDB)
1781 ->object_pt()
1782 ->internal_data_pt(this->P_nst_internal_index)
1783 ->value(0);
1784
1785 double slope4 = octree_pt()
1786 ->son_pt(RUB)
1787 ->object_pt()
1788 ->internal_data_pt(this->P_nst_internal_index)
1789 ->value(0) -
1790 octree_pt()
1791 ->son_pt(LUB)
1792 ->object_pt()
1793 ->internal_data_pt(this->P_nst_internal_index)
1794 ->value(0);
1795
1796
1797 // Use the average
1798 internal_data_pt(this->P_nst_internal_index)
1799 ->set_value(1, 0.25 * (slope1 + slope2 + slope3 + slope4));
1800
1801
1802 // Slope in s_1 direction
1803 //----------------------
1804
1805 // Use average of the 4 FD approximations based on the
1806 // elements central pressure values
1807 // [Other options: Take average of the four
1808 // pressure derivatives]
1809
1810 slope1 = octree_pt()
1811 ->son_pt(LUB)
1812 ->object_pt()
1813 ->internal_data_pt(this->P_nst_internal_index)
1814 ->value(0) -
1815 octree_pt()
1816 ->son_pt(LDB)
1817 ->object_pt()
1818 ->internal_data_pt(this->P_nst_internal_index)
1819 ->value(0);
1820
1821 slope2 = octree_pt()
1822 ->son_pt(RUB)
1823 ->object_pt()
1824 ->internal_data_pt(this->P_nst_internal_index)
1825 ->value(0) -
1826 octree_pt()
1827 ->son_pt(RDB)
1828 ->object_pt()
1829 ->internal_data_pt(this->P_nst_internal_index)
1830 ->value(0);
1831
1832 slope3 = octree_pt()
1833 ->son_pt(LUF)
1834 ->object_pt()
1835 ->internal_data_pt(this->P_nst_internal_index)
1836 ->value(0) -
1837 octree_pt()
1838 ->son_pt(LDF)
1839 ->object_pt()
1840 ->internal_data_pt(this->P_nst_internal_index)
1841 ->value(0);
1842
1843 slope4 = octree_pt()
1844 ->son_pt(RUF)
1845 ->object_pt()
1846 ->internal_data_pt(this->P_nst_internal_index)
1847 ->value(0) -
1848 octree_pt()
1849 ->son_pt(RDF)
1850 ->object_pt()
1851 ->internal_data_pt(this->P_nst_internal_index)
1852 ->value(0);
1853
1854
1855 // Use the average
1856 internal_data_pt(this->P_nst_internal_index)
1857 ->set_value(2, 0.25 * (slope1 + slope2 + slope3 + slope4));
1858
1859
1860 // Slope in s_2 direction
1861 //----------------------
1862
1863 // Use average of the 4 FD approximations based on the
1864 // elements central pressure values
1865 // [Other options: Take average of the four
1866 // pressure derivatives]
1867
1868 slope1 = octree_pt()
1869 ->son_pt(LUF)
1870 ->object_pt()
1871 ->internal_data_pt(this->P_nst_internal_index)
1872 ->value(0) -
1873 octree_pt()
1874 ->son_pt(LUB)
1875 ->object_pt()
1876 ->internal_data_pt(this->P_nst_internal_index)
1877 ->value(0);
1878
1879 slope2 = octree_pt()
1880 ->son_pt(RUF)
1881 ->object_pt()
1882 ->internal_data_pt(this->P_nst_internal_index)
1883 ->value(0) -
1884 octree_pt()
1885 ->son_pt(RUB)
1886 ->object_pt()
1887 ->internal_data_pt(this->P_nst_internal_index)
1888 ->value(0);
1889
1890 slope3 = octree_pt()
1891 ->son_pt(LDF)
1892 ->object_pt()
1893 ->internal_data_pt(this->P_nst_internal_index)
1894 ->value(0) -
1895 octree_pt()
1896 ->son_pt(LDB)
1897 ->object_pt()
1898 ->internal_data_pt(this->P_nst_internal_index)
1899 ->value(0);
1900
1901 slope4 = octree_pt()
1902 ->son_pt(RDF)
1903 ->object_pt()
1904 ->internal_data_pt(this->P_nst_internal_index)
1905 ->value(0) -
1906 octree_pt()
1907 ->son_pt(RDB)
1908 ->object_pt()
1909 ->internal_data_pt(this->P_nst_internal_index)
1910 ->value(0);
1911
1912 // Use the average
1913 internal_data_pt(this->P_nst_internal_index)
1914 ->set_value(3, 0.25 * (slope1 + slope2 + slope3 + slope4));
1915 }
1916
1917
1918 //======================================================================
1919 /// 2D Further build for Crouzeix_Raviart interpolates the internal
1920 /// pressure dofs from father element: Make sure pressure values and
1921 /// dp/ds agree between fathers and sons at the midpoints of the son
1922 /// elements.
1923 //======================================================================
1924 template<>
1926 {
1927 // Call the generic further build
1929
1930 using namespace QuadTreeNames;
1931
1932 // What type of son am I? Ask my quadtree representation...
1933 int son_type = quadtree_pt()->son_type();
1934
1935 // Pointer to my father (in element impersonation)
1936 RefineableElement* father_el_pt = quadtree_pt()->father_pt()->object_pt();
1937
1938 Vector<double> s_father(2);
1939
1940 // Son midpoint is located at the following coordinates in father element:
1941
1942 // South west son
1943 if (son_type == SW)
1944 {
1945 s_father[0] = -0.5;
1946 s_father[1] = -0.5;
1947 }
1948 // South east son
1949 else if (son_type == SE)
1950 {
1951 s_father[0] = 0.5;
1952 s_father[1] = -0.5;
1953 }
1954 // North east son
1955 else if (son_type == NE)
1956 {
1957 s_father[0] = 0.5;
1958 s_father[1] = 0.5;
1959 }
1960
1961 // North west son
1962 else if (son_type == NW)
1963 {
1964 s_father[0] = -0.5;
1965 s_father[1] = 0.5;
1966 }
1967
1968 // Pressure value in father element
1969 RefineableQCrouzeixRaviartElement<2>* cast_father_element_pt =
1970 dynamic_cast<RefineableQCrouzeixRaviartElement<2>*>(father_el_pt);
1971
1972 double press = cast_father_element_pt->interpolated_p_nst(s_father);
1973
1974 // Pressure value gets copied straight into internal dof:
1975 internal_data_pt(this->P_nst_internal_index)->set_value(0, press);
1976
1977 // The slopes get copied from father
1978 for (unsigned i = 1; i < 3; i++)
1979 {
1980 double half_father_slope =
1981 0.5 *
1982 cast_father_element_pt->internal_data_pt(this->P_nst_internal_index)
1983 ->value(i);
1984 // Set the value in the son
1985 internal_data_pt(this->P_nst_internal_index)
1986 ->set_value(i, half_father_slope);
1987 }
1988 }
1989
1990
1991 //=======================================================================
1992 /// 3D Further build for Crouzeix_Raviart interpolates the internal
1993 /// pressure dofs from father element: Make sure pressure values and
1994 /// dp/ds agree between fathers and sons at the midpoints of the son
1995 /// elements.
1996 //=======================================================================
1997 template<>
1999 {
2001
2002 using namespace OcTreeNames;
2003
2004 // What type of son am I? Ask my octree representation...
2005 int son_type = octree_pt()->son_type();
2006
2007 // Pointer to my father (in element impersonation)
2008 RefineableQElement<3>* father_el_pt = dynamic_cast<RefineableQElement<3>*>(
2009 octree_pt()->father_pt()->object_pt());
2010
2011 Vector<double> s_father(3);
2012
2013 // Son midpoint is located at the following coordinates in father element:
2014 for (unsigned i = 0; i < 3; i++)
2015 {
2016 s_father[i] = 0.5 * OcTree::Direction_to_vector[son_type][i];
2017 }
2018
2019 // Pressure value in father element
2020 RefineableQCrouzeixRaviartElement<3>* cast_father_element_pt =
2021 dynamic_cast<RefineableQCrouzeixRaviartElement<3>*>(father_el_pt);
2022
2023 double press = cast_father_element_pt->interpolated_p_nst(s_father);
2024
2025 // Pressure value gets copied straight into internal dof:
2026 internal_data_pt(this->P_nst_internal_index)->set_value(0, press);
2027
2028 // The slopes get copied from father
2029 for (unsigned i = 1; i < 4; i++)
2030 {
2031 double half_father_slope =
2032 0.5 *
2033 cast_father_element_pt->internal_data_pt(this->P_nst_internal_index)
2034 ->value(i);
2035 // Set the value
2036 internal_data_pt(this->P_nst_internal_index)
2037 ->set_value(i, half_father_slope);
2038 }
2039 }
2040
2041 //=======================================================================
2042 /// 2D
2043 /// Derivatives of the shape functions and test functions w.r.t. to global
2044 /// (Eulerian) coordinates. Return Jacobian of mapping between
2045 /// local and global coordinates.
2046 //=======================================================================
2047 template<>
2049 2>::dshape_and_dtest_eulerian_nst(const Vector<double>& s,
2050 Shape& psi,
2051 DShape& dpsidx,
2052 Shape& test,
2053 DShape& dtestdx) const
2054 {
2055 // Call the geometrical shape functions and derivatives
2056 double J = this->dshape_eulerian(s, psi, dpsidx);
2057
2058 // Loop over the test functions and derivatives and set them equal to the
2059 // shape functions
2060 for (unsigned i = 0; i < nnode_1d() * nnode_1d(); i++)
2061 {
2062 test[i] = psi[i];
2063 dtestdx(i, 0) = dpsidx(i, 0);
2064 dtestdx(i, 1) = dpsidx(i, 1);
2065 }
2066
2067 // Return the jacobian
2068 return J;
2069 }
2070
2071 //=======================================================================
2072 /// 2D
2073 /// Derivatives of the shape functions and test functions w.r.t. to global
2074 /// (Eulerian) coordinates. Return Jacobian of mapping between
2075 /// local and global coordinates.
2076 //=======================================================================
2077 template<>
2079 2>::dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
2080 Shape& psi,
2081 DShape& dpsidx,
2082 Shape& test,
2083 DShape& dtestdx) const
2084 {
2085 // Call the geometrical shape functions and derivatives
2086 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
2087
2088 // Loop over the test functions and derivatives and set them equal to the
2089 // shape functions
2090 for (unsigned i = 0; i < nnode_1d() * nnode_1d(); i++)
2091 {
2092 test[i] = psi[i];
2093 dtestdx(i, 0) = dpsidx(i, 0);
2094 dtestdx(i, 1) = dpsidx(i, 1);
2095 }
2096
2097 // Return the jacobian
2098 return J;
2099 }
2100
2101 //=======================================================================
2102 /// 3D
2103 /// Derivatives of the shape functions and test functions w.r.t. to global
2104 /// (Eulerian) coordinates. Return Jacobian of mapping between
2105 /// local and global coordinates.
2106 //=======================================================================
2107 template<>
2109 3>::dshape_and_dtest_eulerian_nst(const Vector<double>& s,
2110 Shape& psi,
2111 DShape& dpsidx,
2112 Shape& test,
2113 DShape& dtestdx) const
2114 {
2115 // Call the geometrical shape functions and derivatives
2116 double J = this->dshape_eulerian(s, psi, dpsidx);
2117
2118 // Loop over the test functions and derivatives and set them equal to the
2119 // shape functions
2120 for (unsigned i = 0; i < nnode_1d() * nnode_1d() * nnode_1d(); i++)
2121 {
2122 test[i] = psi[i];
2123 dtestdx(i, 0) = dpsidx(i, 0);
2124 dtestdx(i, 1) = dpsidx(i, 1);
2125 dtestdx(i, 2) = dpsidx(i, 2);
2126 }
2127
2128 // Return the jacobian
2129 return J;
2130 }
2131
2132 //=======================================================================
2133 /// 3D
2134 /// Derivatives of the shape functions and test functions w.r.t. to global
2135 /// (Eulerian) coordinates. Return Jacobian of mapping between
2136 /// local and global coordinates.
2137 //=======================================================================
2138 template<>
2140 3>::dshape_and_dtest_eulerian_at_knot_nst(const unsigned& ipt,
2141 Shape& psi,
2142 DShape& dpsidx,
2143 Shape& test,
2144 DShape& dtestdx) const
2145 {
2146 // Call the geometrical shape functions and derivatives
2147 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
2148
2149 // Loop over the test functions and derivatives and set them equal to the
2150 // shape functions
2151 for (unsigned i = 0; i < nnode_1d() * nnode_1d() * nnode_1d(); i++)
2152 {
2153 test[i] = psi[i];
2154 dtestdx(i, 0) = dpsidx(i, 0);
2155 dtestdx(i, 1) = dpsidx(i, 1);
2156 dtestdx(i, 2) = dpsidx(i, 2);
2157 }
2158
2159 // Return the jacobian
2160 return J;
2161 }
2162
2163 //=======================================================================
2164 /// 2D :
2165 /// Pressure shape functions
2166 //=======================================================================
2167 template<>
2169 const Vector<double>& s, Shape& psi) const
2170 {
2171 unsigned npres = this->npres_nst();
2172 if (npres == 1)
2173 {
2174 psi[0] = 1.0;
2175 }
2176 else
2177 {
2178 // Get number of pressure modes
2179 unsigned npres_1d = (int)std::sqrt((double)npres);
2180
2181 // Local storage
2182 // Call the one-dimensional modal shape functions
2183 OneDimensionalModalShape psi1(npres_1d, s[0]);
2184 OneDimensionalModalShape psi2(npres_1d, s[1]);
2185
2186 // Now let's loop over the nodal points in the element
2187 // s1 is the "x" coordinate, s2 the "y"
2188 for (unsigned i = 0; i < npres_1d; i++)
2189 {
2190 for (unsigned j = 0; j < npres_1d; j++)
2191 {
2192 // Multiply the two 1D functions together to get the 2D function
2193 psi[i * npres_1d + j] = psi2[i] * psi1[j];
2194 }
2195 }
2196 }
2197 }
2198
2199 /// Define the pressure shape and test functions
2200 template<>
2202 const Vector<double>& s, Shape& psi, Shape& test) const
2203 {
2204 // Call the pressure shape functions
2205 pshape_nst(s, psi);
2206
2207 // Loop over the test functions and set them equal to the shape functions
2208 if (this->npres_nst() == 1)
2209 {
2210 test[0] = psi[0];
2211 }
2212 else
2213 {
2214 for (unsigned i = 0; i < this->npres_nst(); i++) test[i] = psi[i];
2215 }
2216 }
2217
2218 //=======================================================================
2219 /// 3D :
2220 /// Pressure shape functions
2221 //=======================================================================
2222 template<>
2224 const Vector<double>& s, Shape& psi) const
2225 {
2226 unsigned npres = this->npres_nst();
2227 if (npres == 1)
2228 {
2229 psi[0] = 1.0;
2230 }
2231 else
2232 {
2233 // Get number of pressure modes
2234 unsigned npres_1d = (int)std::sqrt((double)npres);
2235
2236 // Local storage
2237 // Call the one-dimensional modal shape functions
2238 OneDimensionalModalShape psi1(npres_1d, s[0]);
2239 OneDimensionalModalShape psi2(npres_1d, s[1]);
2240 OneDimensionalModalShape psi3(npres_1d, s[2]);
2241
2242 // Now let's loop over the nodal points in the element
2243 // s1 is the "x" coordinate, s2 the "y"
2244 for (unsigned i = 0; i < npres_1d; i++)
2245 {
2246 for (unsigned j = 0; j < npres_1d; j++)
2247 {
2248 for (unsigned k = 0; k < npres_1d; k++)
2249 {
2250 // Multiply the two 1D functions together to get the 2D function
2251 psi[i * npres_1d * npres_1d + j * npres_1d + k] =
2252 psi3[i] * psi2[j] * psi1[k];
2253 }
2254 }
2255 }
2256 }
2257 }
2258
2259 /// Define the pressure shape and test functions
2260 template<>
2262 const Vector<double>& s, Shape& psi, Shape& test) const
2263 {
2264 // Call the pressure shape functions
2265 pshape_nst(s, psi);
2266
2267 // Loop over the test functions and set them equal to the shape functions
2268 if (this->npres_nst() == 1)
2269 {
2270 test[0] = psi[0];
2271 }
2272 else
2273 {
2274 for (unsigned i = 0; i < this->npres_nst(); i++) test[i] = psi[i];
2275 }
2276 }
2277
2278} // namespace oomph
2279
2280#endif
e
Definition: cfortran.h:571
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 unpin(const unsigned &i)
Unpin the i-th stored variable.
Definition: nodes.h:391
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
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
virtual void resize(const unsigned &n_value)
Change (increase) the number of values that may be stored.
Definition: nodes.cc:1002
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
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 Node * get_node_at_local_coordinate(const Vector< double > &s) const
If there is a node at this local coordinate, return the pointer to the node.
Definition: elements.cc:3882
virtual unsigned nvertex_node() const
Return the number of vertex nodes in this element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2491
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
Definition: elements.h:1374
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
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
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2218
virtual void set_integration_scheme(Integral *const &integral_pt)
Set the spatial integration scheme.
Definition: elements.cc:3210
virtual Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element. Broken virtual function in "pure" finite elements.
Definition: elements.h:2500
virtual double local_one_d_fraction_of_node(const unsigned &n1d, const unsigned &i)
Get the local fraction of any node in the n-th position in a one dimensional expansion along the i-th...
Definition: elements.h:1858
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Class for multidimensional Gauss Lobatto Legendre integration rules empty - just establishes template...
Definition: integral.h:1281
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
Definition: elements.h:622
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
A general mesh class.
Definition: mesh.h:67
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
virtual unsigned u_index_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
double * Re_pt
Pointer to global Reynolds number.
void interpolated_u_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
double *& density_ratio_pt()
Pointer to Density ratio.
double *& re_invfr_pt()
Pointer to global inverse Froude number.
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
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)
double *& re_pt()
Pointer to Reynolds number.
Vector< double > * G_pt
Pointer to global gravity Vector.
NavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number)
NavierStokesSourceFctPt & source_fct_pt()
Access function for the source-function pointer.
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
Vector< FpPressureAdvDiffRobinBCElementBase * > Pressure_advection_diffusion_robin_element_pt
Storage for FaceElements that apply Robin BC for pressure adv diff equation used in Fp preconditioner...
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
static Vector< Vector< int > > Direction_to_vector
For each direction, i.e. a son_type (vertex), a face or an edge, this defines a vector that indicates...
Definition: octree.h:353
Non-templated class that returns modal hierachical shape functions based on Legendre polynomials.
Definition: shape.h:1349
An OomphLibError object which should be thrown when an run-time error is encountered....
p-refineable version of Crouzeix Raviart elements. Generic class definitions
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons This must be specialised for each dime...
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
double p_nst(const unsigned &t, const unsigned &i) const
Return the i-th pressure value (Discontinous pressure interpolation – no need to cater for hanging no...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
double p_nst(const unsigned &i) const
Broken assignment operator.
void further_build()
Further build for Crouzeix_Raviart interpolates the internal pressure dofs from father element: Make ...
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 nvertex_node() const
Number of vertex nodes in the element.
PRefineableQCrouzeixRaviartElement(const PRefineableQCrouzeixRaviartElement< DIM > &dummy)=delete
Broken copy constructor.
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at local node n.
void pshape_nst(const Vector< double > &s, Shape &psi, Shape &test) const
Pressure shape and test functions at local coordinte s.
unsigned npres_nst() const
// Return number of pressure values
void unpin_elemental_pressure_dofs()
Unpin all internal pressure dofs.
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 further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation:
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: DIM (velocities)
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 get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get all function values [u,v..,p] at previous timestep t (t=0: present; t>0: previous timestep).
A class that is used to template the p-refineable Q elements by dimension. It's really nothing more t...
Definition: Qelements.h:2274
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
unsigned npres_nst() const
Return number of pressure values.
unsigned P_nst_internal_index
Internal index that indicates at which internal data the pressure is stored.
/////////////////////////////////////////////////////////////////////////
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
virtual int p_nodal_index_nst() const
Set the value at which the pressure is stored in the nodes.
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
unsigned npres_nst() const
Return number of pressure values.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1370
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
virtual RefineableElement * father_element_pt() const
Return a pointer to the father element.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
RefineableFpPressureAdvDiffRobinBCElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
virtual void fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
This function returns the residuals for the traction function. flag=1 (or 0): do (or don't) compute t...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
static void unpin_all_pressure_dofs(const Vector< GeneralisedElement * > &element_pt)
Unpin all pressure dofs in elements listed in vector.
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Pin unused nodal pressure dofs (empty by default, because by default pressure dofs are not associated...
void fill_in_generic_residual_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
void further_build()
Further build, pass the pointers down to the sons.
void dinterpolated_u_nst_ddata(const Vector< double > &s, const unsigned &i, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
Compute the derivatives of the i-th component of velocity at point s with respect to all data that ca...
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
static void pin_redundant_nodal_pressures(const Vector< GeneralisedElement * > &element_pt)
Loop over all elements in Vector (which typically contains all the elements in a fluid mesh) and pin ...
void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed,...
void fill_in_generic_pressure_advection_diffusion_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
virtual Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interp...
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 get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get all function values [u,v..,p] at previous timestep t (t=0: present; t>0: previous timestep).
unsigned ncont_interpolated_values() const
Broken assignment operator.
void unpin_elemental_pressure_dofs()
Unpin all internal pressure dofs.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void further_build()
Further build for Crouzeix_Raviart interpolates the internal pressure dofs from father element: Make ...
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons This must be specialised for each dime...
unsigned nvertex_node() const
Number of vertex nodes in the element.
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
RefineableQCrouzeixRaviartElement(const RefineableQCrouzeixRaviartElement< DIM > &dummy)=delete
Broken copy constructor.
Refineable version of QElement<3,NNODE_1D>.
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
Definition: Qelements.h:2259
Refineable version of Taylor Hood elements. These classes can be written in total generality.
unsigned ninterpolating_node_1d(const int &value_id)
The number of 1d pressure nodes is 2, the number of 1d velocity nodes is the same as the number of 1d...
unsigned ninterpolating_node(const int &value_id)
The number of pressure nodes is 2^DIM. The number of velocity nodes is the same as the number of geom...
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
void unpin_elemental_pressure_dofs()
Unpin all pressure dofs.
Node * interpolating_node_pt(const unsigned &n, const int &value_id)
The velocities are isoparametric and so the "nodes" interpolating the velocities are the geometric no...
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
unsigned nvertex_node() const
Number of vertex nodes in the element.
void pin_elemental_redundant_nodal_pressure_dofs()
Pin all nodal pressure dofs that are not required.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
unsigned required_nvalue(const unsigned &n) const
Number of values required at local node n. In order to simplify matters, we allocate storage for pres...
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
double local_one_d_fraction_of_interpolating_node(const unsigned &n1d, const unsigned &i, const int &value_id)
The pressure nodes are the corner nodes, so when n_value==DIM, the fraction is the same as the 1d nod...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: (DIM velocities + 1 pressure)
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
Node * get_interpolating_node_at_local_coordinate(const Vector< double > &s, const int &value_id)
The velocity nodes are the same as the geometric nodes. The pressure nodes must be calculated by usin...
void build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node.
void interpolating_basis(const Vector< double > &s, Shape &psi, const int &value_id) const
The basis interpolating the pressure is given by pshape(). / The basis interpolating the velocity is ...
void identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...