27 #ifndef OOMPH_REFINEABLE_SPACE_TIME_NAVIER_STOKES_ELEMENTS_HEADER
28 #define OOMPH_REFINEABLE_SPACE_TIME_NAVIER_STOKES_ELEMENTS_HEADER
32 #include <oomph-lib-config.h>
52 template<
class ELEMENT>
53 class RefineableFpPressureAdvDiffRobinBCSpaceTimeElement
54 :
public virtual FpPressureAdvDiffRobinBCSpaceTimeElement<ELEMENT>
74 const unsigned& construct_jacobian_flag);
82 template<
class ELEMENT>
91 OOMPH_CURRENT_FUNCTION,
92 OOMPH_EXCEPTION_LOCATION);
102 unsigned my_dim = this->dim();
119 unsigned n_intpt = this->integral_pt()->nweight();
125 int local_unknown = 0;
128 ELEMENT* bulk_el_pt =
dynamic_cast<ELEMENT*
>(this->bulk_element_pt());
131 unsigned n_pres = bulk_el_pt->npres_nst();
135 int p_index = bulk_el_pt->p_nodal_index_nst();
139 bool pressure_dof_is_hanging[n_pres];
145 for (
unsigned l = 0; l < n_pres; ++l)
148 pressure_dof_is_hanging[l] =
149 bulk_el_pt->pressure_node_pt(l)->is_hanging(p_index);
156 std::ostringstream error_message_stream;
159 error_message_stream <<
"Pressure advection diffusion does not work "
160 <<
"in this case!" << std::endl;
163 throw OomphLibError(error_message_stream.str(),
164 OOMPH_CURRENT_FUNCTION,
165 OOMPH_EXCEPTION_LOCATION);
168 for (
unsigned l = 0; l < n_pres; ++l)
171 pressure_dof_is_hanging[l] =
false;
176 double re = bulk_el_pt->re();
179 Shape psip(n_pres), testp(n_pres);
182 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
185 double w = this->integral_pt()->weight(ipt);
188 for (
unsigned i = 0;
i < my_dim;
i++)
189 s[
i] = this->integral_pt()->knot(ipt,
i);
192 s_bulk = this->local_coordinate_in_bulk(
s);
195 this->outer_unit_normal(ipt, unit_normal);
198 bulk_el_pt->interpolated_u_nst(s_bulk, velocity);
202 for (
unsigned i = 0;
i < my_dim + 1;
i++)
204 flux += velocity[
i] * unit_normal[
i];
208 if (flux > 0.0) flux = 0.0;
211 double interpolated_press = bulk_el_pt->interpolated_p_nst(s_bulk);
214 bulk_el_pt->pshape_nst(s_bulk, psip, testp);
217 double J = this->J_eulerian(
s);
224 unsigned n_master = 1;
225 double hang_weight = 1.0;
228 for (
unsigned l = 0; l < n_pres; l++)
231 if (pressure_dof_is_hanging[l])
235 hang_info_pt = bulk_el_pt->pressure_node_pt(l)->hanging_pt(p_index);
238 n_master = hang_info_pt->
nmaster();
247 for (
unsigned m = 0; m < n_master; m++)
251 if (pressure_dof_is_hanging[l])
254 local_eqn = bulk_el_pt->local_hang_eqn(
261 local_eqn = bulk_el_pt->p_local_eqn(l);
268 residuals[local_eqn] -=
269 re * flux * interpolated_press * testp[l] *
W * hang_weight;
275 unsigned n_master2 = 1;
276 double hang_weight2 = 1.0;
279 for (
unsigned l2 = 0; l2 < n_pres; l2++)
282 if (pressure_dof_is_hanging[l2])
285 bulk_el_pt->pressure_node_pt(l2)->hanging_pt(p_index);
288 n_master2 = hang_info2_pt->
nmaster();
297 for (
unsigned m2 = 0; m2 < n_master2; m2++)
301 if (pressure_dof_is_hanging[l2])
304 local_unknown = bulk_el_pt->local_hang_eqn(
311 local_unknown = bulk_el_pt->p_local_eqn(l2);
316 if (local_unknown >= 0)
318 jacobian(local_eqn, local_unknown) -=
319 (re * flux * psip[l2] * testp[l] *
W * hang_weight *
338 template<
unsigned DIM>
339 class RefineableSpaceTimeNavierStokesEquations
340 :
public virtual SpaceTimeNavierStokesEquations<DIM>,
341 public virtual RefineableElement,
342 public virtual ElementWithZ2ErrorEstimator
376 unsigned n_element = element_pt.size();
379 for (
unsigned e = 0;
e < n_element;
e++)
394 unsigned n_element = element_pt.size();
397 for (
unsigned e = 0;
e < n_element;
e++)
423 const unsigned& which_one = 0);
431 return (DIM + (DIM * (DIM - 1)) / 2) + DIM;
441 unsigned num_entries = (DIM + (DIM * (DIM - 1)) / 2) + DIM;
444 if (flux.size() < num_entries)
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;
451 OOMPH_CURRENT_FUNCTION,
452 OOMPH_EXCEPTION_LOCATION);
466 for (
unsigned i = 0;
i < DIM;
i++)
469 flux[icount] = strainrate(
i,
i);
476 for (
unsigned i = 0;
i < DIM;
i++)
479 for (
unsigned j =
i + 1; j < DIM; j++)
482 flux[icount] = strainrate(
i, j);
490 unsigned n_node = this->
nnode();
496 DShape dpsifdx(n_node, DIM + 1);
502 for (
unsigned j = 0; j < DIM; j++)
508 for (
unsigned l = 0; l < n_node; l++)
511 flux[icount] += this->
nodal_value(l, u_nodal_index) * dpsifdx(l, DIM);
535 this->
Re_pt = cast_father_element_pt->re_pt();
538 this->
St_pt = cast_father_element_pt->st_pt();
543 cast_father_element_pt->is_strouhal_stored_as_external_data();
549 unsigned data_index = 0;
553 cast_father_element_pt->external_data_pt(data_index));
557 this->
ReInvFr_pt = cast_father_element_pt->re_invfr_pt();
560 this->
G_pt = cast_father_element_pt->g_pt();
566 this->
Source_fct_pt = cast_father_element_pt->source_fct_pt();
585 unsigned n_node = this->
nnode();
592 const unsigned u_nodal_index = this->
u_index_nst(i);
600 unsigned n_u_dof = 0;
601 for (
unsigned l = 0; l < n_node; l++)
603 unsigned n_master = 1;
612 n_master = hang_info_pt->
nmaster();
621 for (
unsigned m = 0; m < n_master; m++)
645 du_ddata.resize(n_u_dof, 0.0);
646 global_eqn_number.resize(n_u_dof, 0);
652 for (
unsigned l = 0; l < n_node; l++)
655 unsigned n_master = 1;
658 double hang_weight = 1.0;
670 n_master = hang_info_pt->
nmaster();
680 for (
unsigned m = 0; m < n_master; m++)
711 global_eqn_number[count] = global_eqn;
714 du_ddata[count] = psi[l] * hang_weight;
732 const unsigned& compute_jacobian_flag);
741 const unsigned& compute_jacobian_flag);
757 template<
unsigned DIM>
829 values.resize(DIM + 1, 0.0);
832 for (
unsigned i = 0;
i < DIM;
i++)
860 std::ostringstream error_message_stream;
863 error_message_stream <<
"History values don't make sense in "
864 <<
"space-time elements!" << std::endl;
868 OOMPH_CURRENT_FUNCTION,
869 OOMPH_EXCEPTION_LOCATION);
946 unsigned total_index = 0;
955 for (
unsigned i = 0;
i < DIM + 1;
i++)
965 else if (
s[
i] == 1.0)
975 double float_index = 0.5 * (1.0 +
s[
i]) * (
nnode_1d - 1);
978 index[
i] = int(float_index);
982 double excess = float_index - index[
i];
995 index[
i] *
static_cast<unsigned>(pow(
static_cast<float>(
nnode_1d),
996 static_cast<int>(
i)));
1016 if (value_id == DIM)
1035 if (value_id == DIM)
1038 return static_cast<unsigned>(pow(2.0,
static_cast<int>(DIM + 1)));
1044 return this->
nnode();
1053 const int& value_id)
const
1056 if (value_id == DIM)
1065 return this->
shape(s, psi);
1091 std::set<std::pair<Data*, unsigned>>& paired_load_data)
1094 unsigned u_index[DIM];
1097 for (
unsigned i = 0;
i < DIM;
i++)
1104 unsigned n_node = this->
nnode();
1107 for (
unsigned n = 0; n < n_node; n++)
1119 for (
unsigned j = 0; j < nmaster; j++)
1126 for (
unsigned i = 0;
i < DIM;
i++)
1129 paired_load_data.insert(
1130 std::make_pair(master_nod_pt, u_index[
i]));
1139 for (
unsigned i = 0;
i < DIM;
i++)
1142 paired_load_data.insert(
1143 std::make_pair(this->
node_pt(n), u_index[
i]));
1155 for (
unsigned l = 0; l < n_pres; l++)
1168 unsigned nmaster = hang_info_pt->
nmaster();
1171 for (
unsigned m = 0; m < nmaster; m++)
1175 paired_load_data.insert(
1184 paired_load_data.insert(std::make_pair(pres_node_pt, p_index));
1197 unsigned n_node = this->
nnode();
1200 for (
unsigned n = 0; n < n_node; n++)
1215 unsigned n_node = this->
nnode();
1218 for (
unsigned n = 0; n < n_node; n++)
1228 for (
unsigned l = 0; l < n_pres; l++)
1237 nod_pt->
unpin(p_index);
1248 template<
unsigned DIM>
1249 class FaceGeometry<RefineableQTaylorHoodSpaceTimeElement<DIM>>
1250 :
public virtual FaceGeometry<QTaylorHoodSpaceTimeElement<DIM>>
1263 template<
unsigned DIM>
1264 class FaceGeometry<FaceGeometry<RefineableQTaylorHoodSpaceTimeElement<DIM>>>
1265 :
public virtual FaceGeometry<
1266 FaceGeometry<QTaylorHoodSpaceTimeElement<DIM>>>
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
void pin(const unsigned &i)
Pin the i-th stored variable.
void unpin(const unsigned &i)
Unpin the i-th stored variable.
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....
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
FaceGeometry()
Constructor; empty.
FaceGeometry()
Constructor; empty.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
A general Finite Element class.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
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...
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.
virtual unsigned nvertex_node() const
Return the number of vertex nodes in this element. Broken virtual function in "pure" finite elements.
static const double Node_location_tolerance
Default value for the tolerance to be used when locating nodes via local coordinates.
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.
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.
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...
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
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...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Class that contains data for hanging nodes.
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
unsigned nmaster() const
Return the number of master nodes.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
bool is_hanging() const
Test whether the node is geometrically hanging.
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.
unsigned npres_nst() const
Return number of pressure values.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
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...
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...
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...
unsigned ncont_interpolated_values() const
Number of continuously interpolated values.
Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node.
RefineableQTaylorHoodSpaceTimeElement()
Constructor.
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...
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
unsigned nvertex_node() const
Number of vertex nodes in the element.
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...
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 * 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 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...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void identify_load_data(std::set< std::pair< Data *, unsigned >> &paired_load_data)
Add to the set paired_load_data pairs containing.
void unpin_elemental_pressure_dofs()
Unpin all pressure dofs.
void pin_elemental_redundant_nodal_pressure_dofs()
Pin all nodal pressure dofs that are not required.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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 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.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
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,...
RefineableSpaceTimeNavierStokesEquations()
Constructor.
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
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...
void further_build()
Further build, pass the pointers down to the sons.
void dinterpolated_u_nst_ddata(const Vector< double > &s, const unsigned &i, Vector< double > &du_ddata, Vector< unsigned > &global_eqn_number)
Compute the derivatives of the i-th component of velocity at point s with respect to all data that ca...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number)
NavierStokesBodyForceFctPt Body_force_fct_pt
Pointer to body force function.
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)
Vector< double > * G_pt
Pointer to global gravity Vector.
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...
double * Re_pt
Pointer to global Reynolds number.
NavierStokesSourceFctPt Source_fct_pt
Pointer to volumetric source function.
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....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...