refineable_mixed_order_petrov_galerkin_space_time_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-2023 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_MIXED_ORDER_PETROV_GALERKIN_SPACE_TIME_NAVIER_STOKES_ELEMENTS_HEADER
28 #define OOMPH_REFINEABLE_MIXED_ORDER_PETROV_GALERKIN_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 
42 namespace 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>
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>(),
62  FaceElement(),
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  // Number of master nodes and storage for the weight of the shape function
223  unsigned n_master = 1;
224  double hang_weight = 1.0;
225 
226  // Loop over the pressure shape functions
227  for (unsigned l = 0; l < n_pres; l++)
228  {
229  // If the pressure dof is hanging
230  if (pressure_dof_is_hanging[l])
231  {
232  // Pressure dof is hanging so it must be nodal-based
233  // Get the hang info object
234  hang_info_pt = bulk_el_pt->pressure_node_pt(l)->hanging_pt(p_index);
235 
236  // Get the number of master nodes from the pressure node
237  n_master = hang_info_pt->nmaster();
238  }
239  // Otherwise the node is its own master
240  else
241  {
242  n_master = 1;
243  }
244 
245  // Loop over the master nodes
246  for (unsigned m = 0; m < n_master; m++)
247  {
248  // Get the number of the unknown
249  // If the pressure dof is hanging
250  if (pressure_dof_is_hanging[l])
251  {
252  // Get the local equation from the master node
253  local_eqn = bulk_el_pt->local_hang_eqn(
254  hang_info_pt->master_node_pt(m), p_index);
255  // Get the weight for the node
256  hang_weight = hang_info_pt->master_weight(m);
257  }
258  else
259  {
260  local_eqn = bulk_el_pt->p_local_eqn(l);
261  hang_weight = 1.0;
262  }
263 
264  // If the equation is not pinned
265  if (local_eqn >= 0)
266  {
267  residuals[local_eqn] -=
268  re * flux * interpolated_press * testp[l] * W * hang_weight;
269 
270  // Jacobian too?
271  if (flag)
272  {
273  // Number of master nodes and weights
274  unsigned n_master2 = 1;
275  double hang_weight2 = 1.0;
276 
277  // Loop over the pressure shape functions
278  for (unsigned l2 = 0; l2 < n_pres; l2++)
279  {
280  // If the pressure dof is hanging
281  if (pressure_dof_is_hanging[l2])
282  {
283  hang_info2_pt =
284  bulk_el_pt->pressure_node_pt(l2)->hanging_pt(p_index);
285  // Pressure dof is hanging so it must be nodal-based
286  // Get the number of master nodes from the pressure node
287  n_master2 = hang_info2_pt->nmaster();
288  }
289  // Otherwise the node is its own master
290  else
291  {
292  n_master2 = 1;
293  }
294 
295  // Loop over the master nodes
296  for (unsigned m2 = 0; m2 < n_master2; m2++)
297  {
298  // Get the number of the unknown
299  // If the pressure dof is hanging
300  if (pressure_dof_is_hanging[l2])
301  {
302  // Get the unknown from the master node
303  local_unknown = bulk_el_pt->local_hang_eqn(
304  hang_info2_pt->master_node_pt(m2), p_index);
305  // Get the weight from the hanging object
306  hang_weight2 = hang_info2_pt->master_weight(m2);
307  }
308  else
309  {
310  local_unknown = bulk_el_pt->p_local_eqn(l2);
311  hang_weight2 = 1.0;
312  }
313 
314  // If the unknown is not pinned
315  if (local_unknown >= 0)
316  {
317  jacobian(local_eqn, local_unknown) -=
318  (re * flux * psip[l2] * testp[l] * W * hang_weight *
319  hang_weight2);
320  }
321  } // for (unsigned m2=0;m2<n_master2;m2++)
322  } // for (unsigned l2=0;l2<n_pres;l2++)
323  } // if (flag)
324  } // if (local_eqn>=0)
325  } // for (unsigned m=0;m<n_master;m++)
326  } // for (unsigned l=0;l<n_pres;l++)
327  } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
328  } // End of fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc
329 
330  /// ////////////////////////////////////////////////////////////////////
331  /// ////////////////////////////////////////////////////////////////////
332  /// ////////////////////////////////////////////////////////////////////
333 
334  //======================================================================
335  /// Refineable version of the Navier-Stokes equations
336  //======================================================================
337  template<unsigned DIM>
339  : public virtual SpaceTimeNavierStokesMixedOrderEquations<DIM>,
340  public virtual RefineableElement,
341  public virtual ElementWithZ2ErrorEstimator
342  {
343  protected:
344  /// Unpin all pressure dofs in the element
345  virtual void unpin_elemental_pressure_dofs() = 0;
346 
347  /// Pin unused nodal pressure dofs (empty by default, because
348  /// by default pressure dofs are not associated with nodes)
350 
351  public:
352  /// Constructor
357  {
358  }
359 
360 
361  /// Loop over all elements in Vector (which typically contains
362  /// all the elements in a fluid mesh) and pin the nodal pressure degrees
363  /// of freedom that are not being used. Function uses
364  /// the member function
365  /// - \c RefineableSpaceTimeNavierStokesMixedOrderEquations::
366  /// pin_elemental_redundant_nodal_pressure_dofs()
367  /// .
368  /// which is empty by default and should be implemented for
369  /// elements with nodal pressure degrees of freedom
370  /// (e.g. for refineable Taylor-Hood.)
372  const Vector<GeneralisedElement*>& element_pt)
373  {
374  // How many elements are there?
375  unsigned n_element = element_pt.size();
376 
377  // Loop over all elements
378  for (unsigned e = 0; e < n_element; e++)
379  {
380  // Call the function that pins the unused nodal pressure data
382  element_pt[e])
384  }
385  } // End of pin_redundant_nodal_pressures
386 
387 
388  /// Unpin all pressure dofs in elements listed in vector.
390  const Vector<GeneralisedElement*>& element_pt)
391  {
392  // How many elements are there?
393  unsigned n_element = element_pt.size();
394 
395  // Loop over all elements
396  for (unsigned e = 0; e < n_element; e++)
397  {
398  // Unpin the pressure dofs in the e-th element
400  element_pt[e])
402  }
403  } // End of unpin_all_pressure_dofs
404 
405 
406  /// Pointer to n_p-th pressure node (Default: NULL,
407  /// indicating that pressure is not based on nodal interpolation).
408  virtual Node* pressure_node_pt(const unsigned& n_p)
409  {
410  // Return a null pointer
411  return NULL;
412  } // End of pressure_node_pt
413 
414 
415  /// Compute the diagonal of the velocity/pressure mass matrices.
416  /// If which one=0, both are computed, otherwise only the pressure
417  /// (which_one=1) or the velocity mass matrix (which_one=2 -- the
418  /// LSC version of the preconditioner only needs that one)
420  Vector<double>& press_mass_diag,
421  Vector<double>& veloc_mass_diag,
422  const unsigned& which_one = 0);
423 
424 
425  /// Number of 'flux' terms for Z2 error estimation
426  unsigned num_Z2_flux_terms()
427  {
428  // DIM diagonal strain rates, DIM(DIM-1)/2 off diagonal rates (plus
429  // DIM velocity time-derivatives)
430  return (DIM + (DIM * (DIM - 1)) / 2) + DIM;
431  } // End of num_Z2_flux_terms
432 
433 
434  /// Get 'flux' for Z2 error recovery: Upper triangular entries
435  /// in strain rate tensor.
437  {
438 #ifdef PARANOID
439  // Calculate the number of entries there should be
440  unsigned num_entries = (DIM + (DIM * (DIM - 1)) / 2) + DIM;
441 
442  // Check if the flux vector has the correct size
443  if (flux.size() < num_entries)
444  {
445  std::ostringstream error_message;
446  error_message << "The flux vector has the wrong number of entries, "
447  << flux.size() << ", whereas it should be at least "
448  << num_entries << std::endl;
449  throw OomphLibError(error_message.str(),
450  OOMPH_CURRENT_FUNCTION,
451  OOMPH_EXCEPTION_LOCATION);
452  }
453 #endif
454 
455  // Allocate space for the strain-rate
456  DenseMatrix<double> strainrate(DIM, DIM, 0.0);
457 
458  // Get strain rate matrix
459  this->strain_rate(s, strainrate);
460 
461  // Pack into flux Vector
462  unsigned icount = 0;
463 
464  // Loop over the diagonal entries
465  for (unsigned i = 0; i < DIM; i++)
466  {
467  // Add the next strain rate entry to the flux vector
468  flux[icount] = strainrate(i, i);
469 
470  // Increment the counter
471  icount++;
472  }
473 
474  // Loop over the rows of the matrix
475  for (unsigned i = 0; i < DIM; i++)
476  {
477  // Loop over the lower triangular columns
478  for (unsigned j = i + 1; j < DIM; j++)
479  {
480  // Add the next strain rate entry to the flux vector
481  flux[icount] = strainrate(i, j);
482 
483  // Increment the counter
484  icount++;
485  }
486  } // for (unsigned i=0;i<DIM;i++)
487 
488  // The number of nodes in the element
489  unsigned n_node = this->nnode();
490 
491  // Set up memory for the shape and test functions
492  Shape psif(n_node);
493 
494  // Set up memory for the derivatives of the shape and test functions
495  DShape dpsifdx(n_node, DIM + 1);
496 
497  // Call the derivatives of the shape and test functions
498  dshape_eulerian(s, psif, dpsifdx);
499 
500  // Loop over velocity components
501  for (unsigned j = 0; j < DIM; j++)
502  {
503  // Find the index at which the variable is stored
504  unsigned u_nodal_index = this->u_index_nst(j);
505 
506  // Loop over nodes
507  for (unsigned l = 0; l < n_node; l++)
508  {
509  // Add the time-derivative contribution from the l-th node
510  flux[icount] += this->nodal_value(l, u_nodal_index) * dpsifdx(l, DIM);
511  }
512 
513  // Increment the counter
514  icount++;
515  } // for (unsigned j=0;j<DIM;j++)
516  } // End of get_Z2_flux
517 
518 
519  /// Further build, pass the pointers down to the sons
521  {
522  // Find the father element
524  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->ReSt_pt = cast_father_element_pt->re_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_reynolds_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>
761  public virtual RefineableQElement<DIM + 1>
762  {
763  public:
764  /// Constructor
766  : RefineableElement(),
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
785  unsigned ncont_interpolated_values() const
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.
798  unsigned nrecovery_order()
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.
916  double local_one_d_fraction_of_interpolating_node(const unsigned& n_1d,
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
986  ((1.0 - excess) > FiniteElement::Node_location_tolerance))
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  face_index));
1079  } // End of build_fp_press_adv_diff_robin_bc_element
1080 
1081 
1082  /// Add to the set \c paired_load_data pairs containing
1083  /// - the pointer to a Data object
1084  /// and
1085  /// - the index of the value in that Data object
1086  /// .
1087  /// for all values (pressures, velocities) that affect the
1088  /// load computed in the \c get_load(...) function.
1089  /// (Overloads non-refineable version and takes hanging nodes
1090  /// into account)
1092  std::set<std::pair<Data*, unsigned>>& paired_load_data)
1093  {
1094  // Allocate space for the velocity component indices
1095  unsigned u_index[DIM];
1096 
1097  // Loop over the velocity components
1098  for (unsigned i = 0; i < DIM; i++)
1099  {
1100  // Get the nodal indices at which the velocities are stored
1101  u_index[i] = this->u_index_nst(i);
1102  }
1103 
1104  // Get the number of nodes in the element
1105  unsigned n_node = this->nnode();
1106 
1107  // Loop over the nodes
1108  for (unsigned n = 0; n < n_node; n++)
1109  {
1110  // Pointer to current node
1111  Node* nod_pt = this->node_pt(n);
1112 
1113  // Check if it's hanging:
1114  if (nod_pt->is_hanging())
1115  {
1116  // It's hanging -- get number of master nodes
1117  unsigned nmaster = nod_pt->hanging_pt()->nmaster();
1118 
1119  // Loop over master nodes
1120  for (unsigned j = 0; j < nmaster; j++)
1121  {
1122  // Create a pointer to the j-th master node
1123  Node* master_nod_pt = nod_pt->hanging_pt()->master_node_pt(j);
1124 
1125  // Loop over the velocity components and add a pointer to their data
1126  // and indices to the vectors
1127  for (unsigned i = 0; i < DIM; i++)
1128  {
1129  // Create a pair and add it to the storage
1130  paired_load_data.insert(
1131  std::make_pair(master_nod_pt, u_index[i]));
1132  }
1133  } // for (unsigned j=0;j<nmaster;j++)
1134  }
1135  // If the node is not hanging
1136  else
1137  {
1138  // Loop over the velocity components and add pointer to their data
1139  // and indices to the vectors
1140  for (unsigned i = 0; i < DIM; i++)
1141  {
1142  // Create a pair and add it to the storage
1143  paired_load_data.insert(
1144  std::make_pair(this->node_pt(n), u_index[i]));
1145  }
1146  } // if (nod_pt->is_hanging())
1147  } // for (unsigned n=0;n<n_node;n++)
1148 
1149  // Get the nodal index at which the pressure is stored
1150  int p_index = this->p_nodal_index_nst();
1151 
1152  // Get the number of pressure degrees of freedom
1153  unsigned n_pres = this->npres_nst();
1154 
1155  // Loop over the pressure data
1156  for (unsigned l = 0; l < n_pres; l++)
1157  {
1158  // Get the pointer to the l-th pressure node
1159  Node* pres_node_pt = this->pressure_node_pt(l);
1160 
1161  // Check if the pressure dof is hanging
1162  if (pres_node_pt->is_hanging(p_index))
1163  {
1164  // Get the pointer to the hang info object (pressure is stored
1165  // as p_index-th nodal dof).
1166  HangInfo* hang_info_pt = pres_node_pt->hanging_pt(p_index);
1167 
1168  // Get number of pressure master nodes (pressure is stored
1169  unsigned nmaster = hang_info_pt->nmaster();
1170 
1171  // Loop over pressure master nodes
1172  for (unsigned m = 0; m < nmaster; m++)
1173  {
1174  // The p_index-th entry in each nodal data is the pressure, which
1175  // affects the traction
1176  paired_load_data.insert(
1177  std::make_pair(hang_info_pt->master_node_pt(m), p_index));
1178  }
1179  }
1180  // If the pressure dof is not hanging
1181  else
1182  {
1183  // The p_index-th entry in each nodal data is the pressure, which
1184  // affects the traction
1185  paired_load_data.insert(std::make_pair(pres_node_pt, p_index));
1186  } // if (pres_node_pt->is_hanging(p_index))
1187  } // for (unsigned l=0;l<n_pres;l++)
1188  } // End of identify_load_data
1189 
1190  private:
1191  /// Unpin all pressure dofs
1193  {
1194  // Find the index at which the pressure is stored
1195  int p_index = this->p_nodal_index_nst();
1196 
1197  // Get the number of nodes in the element
1198  unsigned n_node = this->nnode();
1199 
1200  // Loop over nodes
1201  for (unsigned n = 0; n < n_node; n++)
1202  {
1203  // Unpin the p_index-th dof at node n
1204  this->node_pt(n)->unpin(p_index);
1205  }
1206  } // End of unpin_elemental_pressure_dofs
1207 
1208 
1209  /// Pin all nodal pressure dofs that are not required
1211  {
1212  // Find the pressure index
1213  int p_index = this->p_nodal_index_nst();
1214 
1215  // Loop over all nodes
1216  unsigned n_node = this->nnode();
1217 
1218  // Loop over all nodes and pin all the nodal pressures
1219  for (unsigned n = 0; n < n_node; n++)
1220  {
1221  // Pin the p_index-th dof at node n
1222  this->node_pt(n)->pin(p_index);
1223  }
1224 
1225  // Loop over all actual pressure nodes and unpin if they're not hanging
1226  unsigned n_pres = this->npres_nst();
1227 
1228  // Loop over all pressure nodes
1229  for (unsigned l = 0; l < n_pres; l++)
1230  {
1231  // Get a pointer to the l-th pressure node
1232  Node* nod_pt = this->pressure_node_pt(l);
1233 
1234  // If the node isn't hanging
1235  if (!nod_pt->is_hanging(p_index))
1236  {
1237  // Unpin the p_index-th dof at this node
1238  nod_pt->unpin(p_index);
1239  }
1240  } // for (unsigned l=0;l<n_pres;l++)
1241  } // End of pin_elemental_redundant_nodal_pressure_dofs
1242  };
1243 
1244 
1245  //=======================================================================
1246  /// Face geometry of the class:
1247  /// RefineableQTaylorHoodMixedOrderSpaceTimeElements
1248  /// is the same as the Face geometry of:
1249  /// QTaylorHoodMixedOrderSpaceTimeElements.
1250  //=======================================================================
1251  template<unsigned DIM>
1253  : public virtual FaceGeometry<QTaylorHoodMixedOrderSpaceTimeElement<DIM>>
1254  {
1255  public:
1256  /// Constructor; empty
1258  {
1259  }
1260  };
1261 
1262 
1263  //=======================================================================
1264  /// Face geometry of the face geometry of the
1265  /// RefineableQTaylorHoodMixedOrderSpaceTimeElements is the same as the Face
1266  /// geometry of the Face geometry of QTaylorHoodMixedOrderSpaceTimeElements.
1267  //=======================================================================
1268  template<unsigned DIM>
1271  : public virtual FaceGeometry<
1272  FaceGeometry<QTaylorHoodMixedOrderSpaceTimeElement<DIM>>>
1273  {
1274  public:
1275  /// Constructor; empty
1278  {
1279  }
1280  };
1281 
1282 } // End of namespace oomph
1283 
1284 #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
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:367
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
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
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
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
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
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
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 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 that contains data for hanging nodes.
Definition: nodes.h:742
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
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
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
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
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
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.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
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...
RefineableFpPressureAdvDiffRobinBCMixedOrderSpaceTimeElement(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, const unsigned &construct_jacobian_flag)
This function returns the residuals for the traction function. If construct_jacobian_flag=1 (or 0): d...
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 identify_load_data(std::set< std::pair< Data *, unsigned >> &paired_load_data)
Add to the set paired_load_data pairs containing.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
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...
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...
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...
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...
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...
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...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
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...
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 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 further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
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_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....
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
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...
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 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...
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 get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
NavierStokesMixedOrderBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
void interpolated_u_nst(const Vector< double > &s, Vector< double > &velocity) const
Compute vector of FE interpolated velocity u at local coordinate s.
bool ReynoldsStrouhal_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....
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number)
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)
virtual double interpolated_p_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
NavierStokesMixedOrderSourceFctPt Source_fct_pt
Pointer to volumetric source function.
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
void store_reynolds_strouhal_as_external_data(Data *reynolds_strouhal_data_pt)
Function that tells us whether the period is stored as external data.
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,...
Vector< FpPressureAdvDiffRobinBCMixedOrderSpaceTimeElementBase * > Pressure_advection_diffusion_robin_element_pt
Storage for FaceElements that apply Robin BC for pressure adv diff equation used in Fp preconditioner...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...