refineable_discontinuous_galerkin_equal_order_pressure_spacetime_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 space-time Navier Stokes elements
27#ifndef OOMPH_REFINEABLE_DISCONTINUOUS_GALERKIN_EQUAL_ORDER_PRESSURE_SPACE_TIME_NAVIER_STOKES_ELEMENTS_HEADER
28#define OOMPH_REFINEABLE_DISCONTINUOUS_GALERKIN_EQUAL_ORDER_PRESSURE_SPACE_TIME_NAVIER_STOKES_ELEMENTS_HEADER
29
30// Config header generated by autoconfig
31#ifdef HAVE_CONFIG_H
32#include <oomph-lib-config.h>
33#endif
34
35// Oomph-lib headers
41
42namespace oomph
43{
44 //======================================================================
45 /// A class for elements that allow the imposition of Robin boundary
46 /// conditions for the pressure advection diffusion problem in the
47 /// Fp preconditioner.
48 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
49 /// class and and thus, we can be generic enough without the need to have
50 /// a separate equations class
51 //======================================================================
52 template<class ELEMENT>
54 : public virtual FpPressureAdvDiffRobinBCSpaceTimeElement<ELEMENT>
55 {
56 public:
57 /// Constructor, which takes a "bulk" element and the value of the
58 /// index and its limit
60 FiniteElement* const& element_pt, const int& face_index)
61 : FaceGeometry<ELEMENT>(),
64 element_pt, face_index, true)
65 {
66 }
67
68 /// This function returns the residuals for the traction
69 /// function. If construct_jacobian_flag=1 (or 0): do (or don't)
70 /// compute the Jacobian as well.
72 Vector<double>& residuals,
73 DenseMatrix<double>& jacobian,
74 const unsigned& construct_jacobian_flag);
75 };
76
77
78 //======================================================================
79 /// Get residuals and Jacobian of Robin boundary conditions in pressure
80 /// advection diffusion problem in Fp preconditoner. Refineable version.
81 //======================================================================
82 template<class ELEMENT>
85 Vector<double>& residuals,
86 DenseMatrix<double>& jacobian,
87 const unsigned& flag)
88 {
89 // Throw a warning
90 throw OomphLibError("You shouldn't be using this just yet!",
91 OOMPH_CURRENT_FUNCTION,
92 OOMPH_EXCEPTION_LOCATION);
93
94 // Pointers to first hang info object
95 HangInfo* hang_info_pt = 0;
96
97 // Pointers to second hang info object
98 HangInfo* hang_info2_pt = 0;
99
100 // Get the dimension of the element
101 // DRAIG: Should be 2 as bulk (space-time) element is 3D...
102 unsigned my_dim = this->dim();
103
104 // Storage for local coordinates in FaceElement
105 Vector<double> s(my_dim, 0.0);
106
107 // Storage for local coordinates in the associated bulk element
108 Vector<double> s_bulk(my_dim + 1, 0.0);
109
110 // Storage for outer unit normal
111 // DRAIG: Need to be careful here...
112 Vector<double> unit_normal(my_dim + 1, 0.0);
113
114 // Storage for velocity in bulk element
115 // DRAIG: Need to be careful here...
116 Vector<double> velocity(my_dim, 0.0);
117
118 // Set the value of n_intpt
119 unsigned n_intpt = this->integral_pt()->nweight();
120
121 // Integer to store local equation number
122 int local_eqn = 0;
123
124 // Integer to store local unknown number
125 int local_unknown = 0;
126
127 // Get upcast bulk element
128 ELEMENT* bulk_el_pt = dynamic_cast<ELEMENT*>(this->bulk_element_pt());
129
130 // Find out how many pressure dofs there are in the bulk element
131 unsigned n_pres = bulk_el_pt->npres_nst();
132
133 // Which nodal value represents the pressure? (Negative if pressure
134 // is not based on nodal interpolation).
135 int p_index = bulk_el_pt->p_nodal_index_nst();
136
137 // Local array of booleans that are true if the l-th pressure value is
138 // hanging (avoid repeated virtual function calls)
139 bool pressure_dof_is_hanging[n_pres];
140
141 // If the pressure is stored at a node
142 if (p_index >= 0)
143 {
144 // Read out whether the pressure is hanging
145 for (unsigned l = 0; l < n_pres; ++l)
146 {
147 // Check the hang status of the pressure node
148 pressure_dof_is_hanging[l] =
149 bulk_el_pt->pressure_node_pt(l)->is_hanging(p_index);
150 }
151 }
152 // Otherwise the pressure is not stored at a node and so cannot hang
153 else
154 {
155 // Create an output stream
156 std::ostringstream error_message_stream;
157
158 // Create an error message
159 error_message_stream << "Pressure advection diffusion does not work "
160 << "in this case!" << std::endl;
161
162 // Throw an error
163 throw OomphLibError(error_message_stream.str(),
164 OOMPH_CURRENT_FUNCTION,
165 OOMPH_EXCEPTION_LOCATION);
166
167 // Loop over the pressure dofs
168 for (unsigned l = 0; l < n_pres; ++l)
169 {
170 // DRAIG: This makes no sense...
171 pressure_dof_is_hanging[l] = false;
172 }
173 } // if (p_index>=0)
174
175 // Get the Reynolds number from the bulk element
176 double re = bulk_el_pt->re();
177
178 // Set up memory for pressure shape and test functions
179 Shape psip(n_pres), testp(n_pres);
180
181 // Loop over the integration points
182 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
183 {
184 // Get the integral weight
185 double w = this->integral_pt()->weight(ipt);
186
187 // Assign values of local coordinate in FaceElement
188 for (unsigned i = 0; i < my_dim; i++)
189 s[i] = this->integral_pt()->knot(ipt, i);
190
191 // Find corresponding coordinate in the the bulk element
192 s_bulk = this->local_coordinate_in_bulk(s);
193
194 /// Get outer unit normal
195 this->outer_unit_normal(ipt, unit_normal);
196
197 // Get velocity in bulk element
198 bulk_el_pt->interpolated_u_nst(s_bulk, velocity);
199
200 // Get normal component of veloc
201 double flux = 0.0;
202 for (unsigned i = 0; i < my_dim + 1; i++)
203 {
204 flux += velocity[i] * unit_normal[i];
205 }
206
207 // Modify bc: If outflow (flux>0) apply Neumann condition instead
208 if (flux > 0.0) flux = 0.0;
209
210 // Get pressure
211 double interpolated_press = bulk_el_pt->interpolated_p_nst(s_bulk);
212
213 // Call the pressure shape and test functions in bulk element
214 bulk_el_pt->pshape_nst(s_bulk, psip, testp);
215
216 // Find the Jacobian of the mapping within the FaceElement
217 double J = this->J_eulerian(s);
218
219 // Premultiply the weights and the Jacobian
220 double W = w * J;
221
222
223 // Number of master nodes and storage for the weight of the shape function
224 unsigned n_master = 1;
225 double hang_weight = 1.0;
226
227 // Loop over the pressure shape functions
228 for (unsigned l = 0; l < n_pres; l++)
229 {
230 // If the pressure dof is hanging
231 if (pressure_dof_is_hanging[l])
232 {
233 // Pressure dof is hanging so it must be nodal-based
234 // Get the hang info object
235 hang_info_pt = bulk_el_pt->pressure_node_pt(l)->hanging_pt(p_index);
236
237 // Get the number of master nodes from the pressure node
238 n_master = hang_info_pt->nmaster();
239 }
240 // Otherwise the node is its own master
241 else
242 {
243 n_master = 1;
244 }
245
246 // Loop over the master nodes
247 for (unsigned m = 0; m < n_master; m++)
248 {
249 // Get the number of the unknown
250 // If the pressure dof is hanging
251 if (pressure_dof_is_hanging[l])
252 {
253 // Get the local equation from the master node
254 local_eqn = bulk_el_pt->local_hang_eqn(
255 hang_info_pt->master_node_pt(m), p_index);
256 // Get the weight for the node
257 hang_weight = hang_info_pt->master_weight(m);
258 }
259 else
260 {
261 local_eqn = bulk_el_pt->p_local_eqn(l);
262 hang_weight = 1.0;
263 }
264
265 // If the equation is not pinned
266 if (local_eqn >= 0)
267 {
268 residuals[local_eqn] -=
269 re * flux * interpolated_press * testp[l] * W * hang_weight;
270
271 // Jacobian too?
272 if (flag)
273 {
274 // Number of master nodes and weights
275 unsigned n_master2 = 1;
276 double hang_weight2 = 1.0;
277
278 // Loop over the pressure shape functions
279 for (unsigned l2 = 0; l2 < n_pres; l2++)
280 {
281 // If the pressure dof is hanging
282 if (pressure_dof_is_hanging[l2])
283 {
284 hang_info2_pt =
285 bulk_el_pt->pressure_node_pt(l2)->hanging_pt(p_index);
286 // Pressure dof is hanging so it must be nodal-based
287 // Get the number of master nodes from the pressure node
288 n_master2 = hang_info2_pt->nmaster();
289 }
290 // Otherwise the node is its own master
291 else
292 {
293 n_master2 = 1;
294 }
295
296 // Loop over the master nodes
297 for (unsigned m2 = 0; m2 < n_master2; m2++)
298 {
299 // Get the number of the unknown
300 // If the pressure dof is hanging
301 if (pressure_dof_is_hanging[l2])
302 {
303 // Get the unknown from the master node
304 local_unknown = bulk_el_pt->local_hang_eqn(
305 hang_info2_pt->master_node_pt(m2), p_index);
306 // Get the weight from the hanging object
307 hang_weight2 = hang_info2_pt->master_weight(m2);
308 }
309 else
310 {
311 local_unknown = bulk_el_pt->p_local_eqn(l2);
312 hang_weight2 = 1.0;
313 }
314
315 // If the unknown is not pinned
316 if (local_unknown >= 0)
317 {
318 jacobian(local_eqn, local_unknown) -=
319 (re * flux * psip[l2] * testp[l] * W * hang_weight *
320 hang_weight2);
321 }
322 } // for (unsigned m2=0;m2<n_master2;m2++)
323 } // for (unsigned l2=0;l2<n_pres;l2++)
324 } // if (flag)
325 } // if (local_eqn>=0)
326 } // for (unsigned m=0;m<n_master;m++)
327 } // for (unsigned l=0;l<n_pres;l++)
328 } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
329 } // End of fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc
330
331 /// ////////////////////////////////////////////////////////////////////
332 /// ////////////////////////////////////////////////////////////////////
333 /// ////////////////////////////////////////////////////////////////////
334
335 //======================================================================
336 /// Refineable version of the Navier-Stokes equations
337 //======================================================================
338 template<unsigned DIM>
340 : public virtual SpaceTimeNavierStokesEquations<DIM>,
341 public virtual RefineableElement,
342 public virtual ElementWithZ2ErrorEstimator
343 {
344 protected:
345 /// Unpin all pressure dofs in the element
347
348 /// Pin unused nodal pressure dofs (empty by default, because
349 /// by default pressure dofs are not associated with nodes)
351
352 public:
353 /// Constructor
358 {
359 }
360
361
362 /// Loop over all elements in Vector (which typically contains
363 /// all the elements in a fluid mesh) and pin the nodal pressure degrees
364 /// of freedom that are not being used. Function uses
365 /// the member function
366 /// - \c RefineableSpaceTimeNavierStokesEquations::
367 /// pin_elemental_redundant_nodal_pressure_dofs()
368 /// .
369 /// which is empty by default and should be implemented for
370 /// elements with nodal pressure degrees of freedom
371 /// (e.g. for refineable Taylor-Hood.)
373 const Vector<GeneralisedElement*>& element_pt)
374 {
375 // How many elements are there?
376 unsigned n_element = element_pt.size();
377
378 // Loop over all elements
379 for (unsigned e = 0; e < n_element; e++)
380 {
381 // Call the function that pins the unused nodal pressure data
383 element_pt[e])
385 }
386 } // End of pin_redundant_nodal_pressures
387
388
389 /// Unpin all pressure dofs in elements listed in vector.
391 const Vector<GeneralisedElement*>& element_pt)
392 {
393 // How many elements are there?
394 unsigned n_element = element_pt.size();
395
396 // Loop over all elements
397 for (unsigned e = 0; e < n_element; e++)
398 {
399 // Unpin the pressure dofs in the e-th element
401 element_pt[e])
403 }
404 } // End of unpin_all_pressure_dofs
405
406
407 /// Pointer to n_p-th pressure node (Default: NULL,
408 /// indicating that pressure is not based on nodal interpolation).
409 virtual Node* pressure_node_pt(const unsigned& n_p)
410 {
411 // Return a null pointer
412 return NULL;
413 } // End of pressure_node_pt
414
415
416 /// Compute the diagonal of the velocity/pressure mass matrices.
417 /// If which one=0, both are computed, otherwise only the pressure
418 /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
419 /// LSC version of the preconditioner only needs that one)
421 Vector<double>& press_mass_diag,
422 Vector<double>& veloc_mass_diag,
423 const unsigned& which_one = 0);
424
425
426 /// Number of 'flux' terms for Z2 error estimation
428 {
429 // DIM diagonal strain rates, DIM(DIM-1)/2 off diagonal rates (plus
430 // DIM velocity time-derivatives)
431 return (DIM + (DIM * (DIM - 1)) / 2) + DIM;
432 } // End of num_Z2_flux_terms
433
434
435 /// Get 'flux' for Z2 error recovery: Upper triangular entries
436 /// in strain rate tensor.
438 {
439#ifdef PARANOID
440 // Calculate the number of entries there should be
441 unsigned num_entries = (DIM + (DIM * (DIM - 1)) / 2) + DIM;
442
443 // Check if the flux vector has the correct size
444 if (flux.size() < num_entries)
445 {
446 std::ostringstream error_message;
447 error_message << "The flux vector has the wrong number of entries, "
448 << flux.size() << ", whereas it should be at least "
449 << num_entries << std::endl;
450 throw OomphLibError(error_message.str(),
451 OOMPH_CURRENT_FUNCTION,
452 OOMPH_EXCEPTION_LOCATION);
453 }
454#endif
455
456 // Allocate space for the strain-rate
457 DenseMatrix<double> strainrate(DIM, DIM, 0.0);
458
459 // Get strain rate matrix
460 this->strain_rate(s, strainrate);
461
462 // Pack into flux Vector
463 unsigned icount = 0;
464
465 // Loop over the diagonal entries
466 for (unsigned i = 0; i < DIM; i++)
467 {
468 // Add the next strain rate entry to the flux vector
469 flux[icount] = strainrate(i, i);
470
471 // Increment the counter
472 icount++;
473 }
474
475 // Loop over the rows of the matrix
476 for (unsigned i = 0; i < DIM; i++)
477 {
478 // Loop over the lower triangular columns
479 for (unsigned j = i + 1; j < DIM; j++)
480 {
481 // Add the next strain rate entry to the flux vector
482 flux[icount] = strainrate(i, j);
483
484 // Increment the counter
485 icount++;
486 }
487 } // for (unsigned i=0;i<DIM;i++)
488
489 // The number of nodes in the element
490 unsigned n_node = this->nnode();
491
492 // Set up memory for the shape and test functions
493 Shape psif(n_node);
494
495 // Set up memory for the derivatives of the shape and test functions
496 DShape dpsifdx(n_node, DIM + 1);
497
498 // Call the derivatives of the shape and test functions
499 dshape_eulerian(s, psif, dpsifdx);
500
501 // Loop over velocity components
502 for (unsigned j = 0; j < DIM; j++)
503 {
504 // Find the index at which the variable is stored
505 unsigned u_nodal_index = this->u_index_nst(j);
506
507 // Loop over nodes
508 for (unsigned l = 0; l < n_node; l++)
509 {
510 // Add the time-derivative contribution from the l-th node
511 flux[icount] += this->nodal_value(l, u_nodal_index) * dpsifdx(l, DIM);
512 }
513
514 // Increment the counter
515 icount++;
516 } // for (unsigned j=0;j<DIM;j++)
517 } // End of get_Z2_flux
518
519
520 /// Further build, pass the pointers down to the sons
522 {
523 // Find the father element
524 RefineableSpaceTimeNavierStokesEquations<DIM>* cast_father_element_pt =
526 this->father_element_pt());
527
528 // Set the viscosity ratio pointer
529 this->Viscosity_Ratio_pt = cast_father_element_pt->viscosity_ratio_pt();
530
531 // Set the density ratio pointer
532 this->Density_Ratio_pt = cast_father_element_pt->density_ratio_pt();
533
534 // Set pointer to global Reynolds number
535 this->Re_pt = cast_father_element_pt->re_pt();
536
537 // Set pointer to global Reynolds number x Strouhal number (=Womersley)
538 this->St_pt = cast_father_element_pt->st_pt();
539
540 // The Strouhal number (which is a proxy for the period here) is not
541 // stored as external data
543 cast_father_element_pt->is_strouhal_stored_as_external_data();
544
545 // If we're storing the Strouhal number as external data
547 {
548 // The index of the external data (which contains the st)
549 unsigned data_index = 0;
550
551 // Get the external data pointer from the father and store it
553 cast_father_element_pt->external_data_pt(data_index));
554 }
555
556 // Set pointer to global Reynolds number x inverse Froude number
557 this->ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
558
559 // Set pointer to global gravity Vector
560 this->G_pt = cast_father_element_pt->g_pt();
561
562 // Set pointer to body force function
563 this->Body_force_fct_pt = cast_father_element_pt->body_force_fct_pt();
564
565 // Set pointer to volumetric source function
566 this->Source_fct_pt = cast_father_element_pt->source_fct_pt();
567
568 // Set the ALE flag
569 this->ALE_is_disabled = cast_father_element_pt->ALE_is_disabled;
570 } // End of further_build
571
572
573 /// Compute the derivatives of the i-th component of
574 /// velocity at point s with respect
575 /// to all data that can affect its value. In addition, return the global
576 /// equation numbers corresponding to the data.
577 /// Overload the non-refineable version to take account of hanging node
578 /// information
580 const unsigned& i,
581 Vector<double>& du_ddata,
582 Vector<unsigned>& global_eqn_number)
583 {
584 // Find the number of nodes in the element
585 unsigned n_node = this->nnode();
586 // Local shape function
587 Shape psi(n_node);
588 // Find values of shape function at the given local coordinate
589 this->shape(s, psi);
590
591 // Find the index at which the velocity component is stored
592 const unsigned u_nodal_index = this->u_index_nst(i);
593
594 // Storage for hang info pointer
595 HangInfo* hang_info_pt = 0;
596 // Storage for global equation
597 int global_eqn = 0;
598
599 // Find the number of dofs associated with interpolated u
600 unsigned n_u_dof = 0;
601 for (unsigned l = 0; l < n_node; l++)
602 {
603 unsigned n_master = 1;
604
605 // Local bool (is the node hanging)
606 bool is_node_hanging = this->node_pt(l)->is_hanging();
607
608 // If the node is hanging, get the number of master nodes
609 if (is_node_hanging)
610 {
611 hang_info_pt = this->node_pt(l)->hanging_pt();
612 n_master = hang_info_pt->nmaster();
613 }
614 // Otherwise there is just one master node, the node itself
615 else
616 {
617 n_master = 1;
618 }
619
620 // Loop over the master nodes
621 for (unsigned m = 0; m < n_master; m++)
622 {
623 // Get the equation number
624 if (is_node_hanging)
625 {
626 // Get the equation number from the master node
627 global_eqn =
628 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
629 }
630 else
631 {
632 // Global equation number
633 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
634 }
635
636 // If it's positive add to the count
637 if (global_eqn >= 0)
638 {
639 ++n_u_dof;
640 }
641 }
642 }
643
644 // Now resize the storage schemes
645 du_ddata.resize(n_u_dof, 0.0);
646 global_eqn_number.resize(n_u_dof, 0);
647
648 // Loop over the nodes again and set the derivatives
649 unsigned count = 0;
650
651 // Loop over the nodes in the element
652 for (unsigned l = 0; l < n_node; l++)
653 {
654 // Initialise the number of master nodes to one
655 unsigned n_master = 1;
656
657 // Initialise the hang weight to one
658 double hang_weight = 1.0;
659
660 // Local bool (is the node hanging)
661 bool is_node_hanging = this->node_pt(l)->is_hanging();
662
663 // If the node is hanging, get the number of master nodes
664 if (is_node_hanging)
665 {
666 // Get the HangInfo pointer associated with the l-th node
667 hang_info_pt = this->node_pt(l)->hanging_pt();
668
669 // How many master nodes does this node have?
670 n_master = hang_info_pt->nmaster();
671 }
672 // Otherwise there is just one master node, the node itself
673 else
674 {
675 // Set n_master to one
676 n_master = 1;
677 }
678
679 // Loop over the master nodes
680 for (unsigned m = 0; m < n_master; m++)
681 {
682 // If the node is hanging get weight from master node
683 if (is_node_hanging)
684 {
685 // Get the hang weight from the master node
686 hang_weight = hang_info_pt->master_weight(m);
687 }
688 else
689 {
690 // Node contributes with full weight
691 hang_weight = 1.0;
692 }
693
694 // Get the equation number
695 if (is_node_hanging)
696 {
697 // Get the equation number from the master node
698 global_eqn =
699 hang_info_pt->master_node_pt(m)->eqn_number(u_nodal_index);
700 }
701 else
702 {
703 // Get the global equation number
704 global_eqn = this->node_pt(l)->eqn_number(u_nodal_index);
705 }
706
707 // If it's a proper degree of freedom
708 if (global_eqn >= 0)
709 {
710 // Set the global equation number
711 global_eqn_number[count] = global_eqn;
712
713 // Set the derivative with respect to the unknown
714 du_ddata[count] = psi[l] * hang_weight;
715
716 // Increase the counter
717 ++count;
718 }
719 } // for (unsigned m=0;m<n_master;m++)
720 } // for (unsigned l=0;l<n_node;l++)
721 } // End of dinterpolated_u_nst_ddata
722
723 protected:
724 /// Add the elements contribution to elemental residual vector
725 /// and/or Jacobian matrix.
726 /// compute_jacobian_flag=0: compute residual vector only
727 /// compute_jacobian_flag=1: compute both
729 Vector<double>& residuals,
730 DenseMatrix<double>& jacobian,
731 DenseMatrix<double>& mass_matrix,
732 const unsigned& compute_jacobian_flag);
733
734
735 /// Compute the residuals for the associated pressure advection
736 /// diffusion problem. Used by the Fp preconditioner.
737 /// flag=1(or 0): do (or don't) compute the Jacobian as well.
739 Vector<double>& residuals,
740 DenseMatrix<double>& jacobian,
741 const unsigned& compute_jacobian_flag);
742
743
744 /// Compute derivatives of elemental residual vector with respect
745 /// to nodal coordinates. Overwrites default implementation in
746 /// FiniteElement base class.
747 /// dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}
749 RankThreeTensor<double>& dresidual_dnodal_coordinates);
750 };
751
752
753 //======================================================================
754 /// Refineable version of Taylor Hood elements. These classes
755 /// can be written in total generality.
756 //======================================================================
757 template<unsigned DIM>
759 : public QTaylorHoodSpaceTimeElement<DIM>,
761 public virtual RefineableQElement<DIM + 1>
762 {
763 public:
764 /// Constructor
768 RefineableQElement<DIM + 1>(),
770 {
771 }
772
773
774 /// Number of values required at local node n. In order to simplify
775 /// matters, we allocate storage for pressure variables at all the nodes
776 /// and then pin those that are not used.
777 unsigned required_nvalue(const unsigned& n) const
778 {
779 // There are DIM velocity components and 1 pressure component
780 return DIM + 1;
781 } // End of required_nvalue
782
783
784 /// Number of continuously interpolated values
786 {
787 // There are DIM velocities and 1 pressure component
788 return DIM + 1;
789 } // End of ncont_interpolated_values
790
791
792 /// Rebuild from sons: empty
793 void rebuild_from_sons(Mesh*& mesh_pt) {}
794
795
796 /// Order of recovery shape functions for Z2 error estimation:
797 /// Same order as shape functions.
799 {
800 // Using quadratic interpolation
801 return 2;
802 } // End of nrecovery_order
803
804
805 /// Number of vertex nodes in the element
806 unsigned nvertex_node() const
807 {
808 // Call the base class implementation of the function
810 } // End of nvertex_node
811
812
813 /// Pointer to the j-th vertex node in the element
814 Node* vertex_node_pt(const unsigned& j) const
815 {
816 // Call the base class implementation of the function
818 } // End of vertex_node_pt
819
820
821 /// Get the function value u in Vector.
822 /// Note: Given the generality of the interface (this function is usually
823 /// called from black-box documentation or interpolation routines), the
824 /// values Vector sets its own size in here.
826 Vector<double>& values)
827 {
828 // Set size of Vector: u,v,p and initialise to zero
829 values.resize(DIM + 1, 0.0);
830
831 // Calculate velocities: values[0],...
832 for (unsigned i = 0; i < DIM; i++)
833 {
834 // Calculate the i-th velocity component at local coordinates s
835 values[i] = this->interpolated_u_nst(s, i);
836 }
837
838 // Calculate the pressure field at local coordinates s
839 values[DIM] = this->interpolated_p_nst(s);
840 } // End of get_interpolated_values
841
842
843 /// Get the function value u in Vector.
844 /// Note: Given the generality of the interface (this function
845 /// is usually called from black-box documentation or interpolation
846 /// routines), the values Vector sets its own size in here.
847 void get_interpolated_values(const unsigned& t,
848 const Vector<double>& s,
849 Vector<double>& values)
850 {
851 // The only time this makes sense (if you can even say that...)
852 if (t == 0)
853 {
854 // Call the other function
855 get_interpolated_values(s, values);
856 }
857 else
858 {
859 // Create an output stream
860 std::ostringstream error_message_stream;
861
862 // Create an error message
863 error_message_stream << "History values don't make sense in "
864 << "space-time elements!" << std::endl;
865
866 // Throw an error message
867 throw OomphLibError(error_message_stream.str(),
868 OOMPH_CURRENT_FUNCTION,
869 OOMPH_EXCEPTION_LOCATION);
870 }
871 } // End of get_interpolated_values
872
873
874 /// Perform additional hanging node procedures for variables
875 /// that are not interpolated by all nodes. The pressures are stored
876 /// at the p_nodal_index_nst-th location in each node
878 {
879 // Set up the hang info for the pressure nodes
880 this->setup_hang_for_value(this->p_nodal_index_nst());
881 } // End of further_setup_hanging_nodes
882
883
884 /// Pointer to n_p-th pressure node
885 Node* pressure_node_pt(const unsigned& n_p)
886 {
887 // Return a pointer to the n_p-th pressure node
888 return this->node_pt(this->Pconv[n_p]);
889 } // End of pressure node_pt
890
891
892 /// The velocities are isoparametric and so the "nodes" interpolating
893 /// the velocities are the geometric nodes. The pressure "nodes" are a
894 /// subset of the nodes, so when value_id==DIM, the n-th pressure
895 /// node is returned.
896 Node* interpolating_node_pt(const unsigned& n, const int& value_id)
897
898 {
899 // The only different nodes are the pressure nodes
900 if (value_id == DIM)
901 {
902 // Return a pointer to the n-th pressure node
903 return this->pressure_node_pt(n);
904 }
905 // The other variables are interpolated via the usual nodes
906 else
907 {
908 // Return a pointer to the n-th regular node
909 return this->node_pt(n);
910 }
911 } // End of interpolating_node_pt
912
913
914 /// The pressure nodes are the corner nodes, so when n_value==DIM,
915 /// the fraction is the same as the 1D node number, 0 or 1.
917 const unsigned& i,
918 const int& value_id)
919 {
920 // If we're dealing with a pressure node
921 if (value_id == DIM)
922 {
923 // The pressure nodes are just located on the boundaries at 0 or 1
924 return double(n_1d);
925 }
926 // Otherwise we're dealing with a velocity node
927 else
928 {
929 // The velocity nodes are the same as the geometric ones
930 return this->local_one_d_fraction_of_node(n_1d, i);
931 }
932 } // End of local_one_d_fraction_of_interpolating_node
933
934
935 /// The velocity nodes are the same as the geometric nodes. The
936 /// pressure nodes must be calculated by using the same methods as
937 /// the geometric nodes, but by recalling that there are only two pressure
938 /// nodes per edge.
940 const int& value_id)
941 {
942 // If we are calculating pressure nodes
943 if (value_id == DIM)
944 {
945 // Storage for the index of the pressure node
946 unsigned total_index = 0;
947
948 // The number of nodes along each 1D edge is 2.
949 unsigned nnode_1d = 2;
950
951 // Storage for the index along each boundary
952 Vector<int> index(DIM + 1);
953
954 // Loop over the coordinate directions
955 for (unsigned i = 0; i < DIM + 1; i++)
956 {
957 // If we are at the lower limit, the index is zero
958 if (s[i] == -1.0)
959 {
960 // We're at the first node
961 index[i] = 0;
962 }
963 // If we are at the upper limit, the index is the number of nodes
964 // minus 1
965 else if (s[i] == 1.0)
966 {
967 // We're on the last node
968 index[i] = nnode_1d - 1;
969 }
970 // Otherwise, we have to calculate the index in general
971 else
972 {
973 // For uniformly spaced nodes this should produce an integer when
974 // s[i] is associated with a nodal location
975 double float_index = 0.5 * (1.0 + s[i]) * (nnode_1d - 1);
976
977 // Take the integer part of the float_index value
978 index[i] = int(float_index);
979
980 // What is the excess. This should be safe because the taking the
981 // integer part rounds down
982 double excess = float_index - index[i];
983
984 // If the excess is bigger than our tolerance there is no node
987 {
988 // As there is no node, return null
989 return 0;
990 }
991 } // if (s[i]==-1.0)
992
993 // Construct the general pressure index from the components.
994 total_index +=
995 index[i] * static_cast<unsigned>(pow(static_cast<float>(nnode_1d),
996 static_cast<int>(i)));
997 } // for (unsigned i=0;i<DIM;i++)
998
999 // If we've got here we have a node, so let's return a pointer to it
1000 return this->pressure_node_pt(total_index);
1001 }
1002 // Otherwise velocity nodes are the same as pressure nodes
1003 else
1004 {
1005 // Call the regular helper function to find the right node
1006 return this->get_node_at_local_coordinate(s);
1007 } // if (value_id==DIM)
1008 } // End of get_interpolating_node_at_local_coordinate
1009
1010
1011 /// The number of 1D pressure nodes is 2, the number of 1D velocity
1012 /// nodes is the same as the number of 1D geometric nodes.
1013 unsigned ninterpolating_node_1d(const int& value_id)
1014 {
1015 // If we're dealing with the pressure nodes
1016 if (value_id == DIM)
1017 {
1018 // Using linear interpolation for the pressure
1019 return 2;
1020 }
1021 // If we're dealing with the velocity nodes
1022 else
1023 {
1024 // Every node is a velocity interpolating node
1025 return this->nnode_1d();
1026 }
1027 } // End of ninterpolating_node_1d
1028
1029
1030 /// The number of pressure nodes is 2^DIM. The number of
1031 /// velocity nodes is the same as the number of geometric nodes.
1032 unsigned ninterpolating_node(const int& value_id)
1033 {
1034 // If we want the pressure basis functions
1035 if (value_id == DIM)
1036 {
1037 // There are 2^{DIM+1} pressure dofs (also interpolating in time)
1038 return static_cast<unsigned>(pow(2.0, static_cast<int>(DIM + 1)));
1039 }
1040 // If we want the velocity basis functions
1041 else
1042 {
1043 // Every node is a velocity interpolating node
1044 return this->nnode();
1045 }
1046 } // End if ninterpolating_node
1047
1048
1049 /// The basis interpolating the pressure is given by pshape().
1050 /// / The basis interpolating the velocity is shape().
1052 Shape& psi,
1053 const int& value_id) const
1054 {
1055 // If we want the pressure basis functions
1056 if (value_id == DIM)
1057 {
1058 // Call the pressure shape function
1059 return this->pshape_nst(s, psi);
1060 }
1061 // If we want the velocity basis functions
1062 else
1063 {
1064 // Call the velocity shape function
1065 return this->shape(s, psi);
1066 }
1067 } // End of interpolating_basis
1068
1069
1070 /// Build FaceElements that apply the Robin boundary condition
1071 /// to the pressure advection diffusion problem required by
1072 /// Fp preconditioner
1073 void build_fp_press_adv_diff_robin_bc_element(const unsigned& face_index)
1074 {
1078 } // End of build_fp_press_adv_diff_robin_bc_element
1079
1080
1081 /// Add to the set \c paired_load_data pairs containing
1082 /// - the pointer to a Data object
1083 /// and
1084 /// - the index of the value in that Data object
1085 /// .
1086 /// for all values (pressures, velocities) that affect the
1087 /// load computed in the \c get_load(...) function.
1088 /// (Overloads non-refineable version and takes hanging nodes
1089 /// into account)
1091 std::set<std::pair<Data*, unsigned>>& paired_load_data)
1092 {
1093 // Allocate space for the velocity component indices
1094 unsigned u_index[DIM];
1095
1096 // Loop over the velocity components
1097 for (unsigned i = 0; i < DIM; i++)
1098 {
1099 // Get the nodal indices at which the velocities are stored
1100 u_index[i] = this->u_index_nst(i);
1101 }
1102
1103 // Get the number of nodes in the element
1104 unsigned n_node = this->nnode();
1105
1106 // Loop over the nodes
1107 for (unsigned n = 0; n < n_node; n++)
1108 {
1109 // Pointer to current node
1110 Node* nod_pt = this->node_pt(n);
1111
1112 // Check if it's hanging:
1113 if (nod_pt->is_hanging())
1114 {
1115 // It's hanging -- get number of master nodes
1116 unsigned nmaster = nod_pt->hanging_pt()->nmaster();
1117
1118 // Loop over master nodes
1119 for (unsigned j = 0; j < nmaster; j++)
1120 {
1121 // Create a pointer to the j-th master node
1122 Node* master_nod_pt = nod_pt->hanging_pt()->master_node_pt(j);
1123
1124 // Loop over the velocity components and add a pointer to their data
1125 // and indices to the vectors
1126 for (unsigned i = 0; i < DIM; i++)
1127 {
1128 // Create a pair and add it to the storage
1129 paired_load_data.insert(
1130 std::make_pair(master_nod_pt, u_index[i]));
1131 }
1132 } // for (unsigned j=0;j<nmaster;j++)
1133 }
1134 // If the node is not hanging
1135 else
1136 {
1137 // Loop over the velocity components and add pointer to their data
1138 // and indices to the vectors
1139 for (unsigned i = 0; i < DIM; i++)
1140 {
1141 // Create a pair and add it to the storage
1142 paired_load_data.insert(
1143 std::make_pair(this->node_pt(n), u_index[i]));
1144 }
1145 } // if (nod_pt->is_hanging())
1146 } // for (unsigned n=0;n<n_node;n++)
1147
1148 // Get the nodal index at which the pressure is stored
1149 int p_index = this->p_nodal_index_nst();
1150
1151 // Get the number of pressure degrees of freedom
1152 unsigned n_pres = this->npres_nst();
1153
1154 // Loop over the pressure data
1155 for (unsigned l = 0; l < n_pres; l++)
1156 {
1157 // Get the pointer to the l-th pressure node
1158 Node* pres_node_pt = this->pressure_node_pt(l);
1159
1160 // Check if the pressure dof is hanging
1161 if (pres_node_pt->is_hanging(p_index))
1162 {
1163 // Get the pointer to the hang info object (pressure is stored
1164 // as p_index-th nodal dof).
1165 HangInfo* hang_info_pt = pres_node_pt->hanging_pt(p_index);
1166
1167 // Get number of pressure master nodes (pressure is stored
1168 unsigned nmaster = hang_info_pt->nmaster();
1169
1170 // Loop over pressure master nodes
1171 for (unsigned m = 0; m < nmaster; m++)
1172 {
1173 // The p_index-th entry in each nodal data is the pressure, which
1174 // affects the traction
1175 paired_load_data.insert(
1176 std::make_pair(hang_info_pt->master_node_pt(m), p_index));
1177 }
1178 }
1179 // If the pressure dof is not hanging
1180 else
1181 {
1182 // The p_index-th entry in each nodal data is the pressure, which
1183 // affects the traction
1184 paired_load_data.insert(std::make_pair(pres_node_pt, p_index));
1185 } // if (pres_node_pt->is_hanging(p_index))
1186 } // for (unsigned l=0;l<n_pres;l++)
1187 } // End of identify_load_data
1188
1189 private:
1190 /// Unpin all pressure dofs
1192 {
1193 // Find the index at which the pressure is stored
1194 int p_index = this->p_nodal_index_nst();
1195
1196 // Get the number of nodes in the element
1197 unsigned n_node = this->nnode();
1198
1199 // Loop over nodes
1200 for (unsigned n = 0; n < n_node; n++)
1201 {
1202 // Unpin the p_index-th dof at node n
1203 this->node_pt(n)->unpin(p_index);
1204 }
1205 } // End of unpin_elemental_pressure_dofs
1206
1207
1208 /// Pin all nodal pressure dofs that are not required
1210 {
1211 // Find the pressure index
1212 int p_index = this->p_nodal_index_nst();
1213
1214 // Loop over all nodes
1215 unsigned n_node = this->nnode();
1216
1217 // Loop over all nodes and pin all the nodal pressures
1218 for (unsigned n = 0; n < n_node; n++)
1219 {
1220 // Pin the p_index-th dof at node n
1221 this->node_pt(n)->pin(p_index);
1222 }
1223
1224 // Loop over all actual pressure nodes and unpin if they're not hanging
1225 unsigned n_pres = this->npres_nst();
1226
1227 // Loop over all pressure nodes
1228 for (unsigned l = 0; l < n_pres; l++)
1229 {
1230 // Get a pointer to the l-th pressure node
1231 Node* nod_pt = this->pressure_node_pt(l);
1232
1233 // If the node isn't hanging
1234 if (!nod_pt->is_hanging(p_index))
1235 {
1236 // Unpin the p_index-th dof at this node
1237 nod_pt->unpin(p_index);
1238 }
1239 } // for (unsigned l=0;l<n_pres;l++)
1240 } // End of pin_elemental_redundant_nodal_pressure_dofs
1241 };
1242
1243
1244 //=======================================================================
1245 /// Face geometry of the RefineableQTaylorHoodSpaceTimeElements is the
1246 /// same as the Face geometry of the QTaylorHoodSpaceTimeElements.
1247 //=======================================================================
1248 template<unsigned DIM>
1250 : public virtual FaceGeometry<QTaylorHoodSpaceTimeElement<DIM>>
1251 {
1252 public:
1253 /// Constructor; empty
1255 };
1256
1257
1258 //=======================================================================
1259 /// Face geometry of the face geometry of the
1260 /// RefineableQTaylorHoodSpaceTimeElements is the same as the Face geometry
1261 /// of the Face geometry of QTaylorHoodSpaceTimeElements.
1262 //=======================================================================
1263 template<unsigned DIM>
1265 : public virtual FaceGeometry<
1266 FaceGeometry<QTaylorHoodSpaceTimeElement<DIM>>>
1267 {
1268 public:
1269 /// Constructor; empty
1272 {
1273 }
1274 };
1275
1276} // End of namespace oomph
1277
1278#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
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
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
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
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
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 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 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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Data *& external_data_pt(const unsigned &i)
Return a pointer to i-th external data object.
Definition: elements.h:659
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
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
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
virtual int p_nodal_index_nst() const
Set the value at which the pressure is stored in the nodes.
void pshape_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
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.
A class for elements that allow the imposition of Robin boundary conditions for the pressure advectio...
virtual void fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &construct_jacobian_flag)
This function returns the residuals for the traction function. If construct_jacobian_flag=1 (or 0): d...
RefineableFpPressureAdvDiffRobinBCSpaceTimeElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the value of the index and its limit.
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.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
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 build_fp_press_adv_diff_robin_bc_element(const unsigned &face_index)
Build FaceElements that apply the Robin boundary condition to the pressure advection diffusion proble...
double local_one_d_fraction_of_interpolating_node(const unsigned &n_1d, 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 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...
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...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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 the function value u in Vector. Note: Given the generality of the interface (this function is usu...
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 identify_load_data(std::set< std::pair< Data *, unsigned > > &paired_load_data)
Add to the set paired_load_data pairs containing.
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...
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...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor.
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, const unsigned &compute_jacobian_flag)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
void fill_in_generic_residual_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, const unsigned &compute_jacobian_flag)
Add the elements contribution to elemental residual vector and/or Jacobian matrix....
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...
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
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 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...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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 * St_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
void store_strouhal_as_external_data(Data *strouhal_data_pt)
Function that tells us whether the period is stored as external data.
Vector< FpPressureAdvDiffRobinBCSpaceTimeElementBase * > Pressure_advection_diffusion_robin_element_pt
Storage for FaceElements that apply Robin BC for pressure adv diff equation used in Fp preconditioner...
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
bool is_strouhal_stored_as_external_data() const
Are we storing the Strouhal number as external data?
NavierStokesBodyForceFctPt & body_force_fct_pt()
Access function for the body-force pointer.
void interpolated_u_nst(const Vector< double > &s, Vector< double > &velocity) const
Compute vector of FE interpolated velocity u at local coordinate s.
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,...
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)
bool Strouhal_is_stored_as_external_data
Boolean to indicate whether or not the Strouhal value is stored as external data (if it's also an unk...
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...