27 #ifndef OOMPH_GENERALISED_NEWTONIAN_REFINEABLE_AXISYMMETRIC_NAVIER_STOKES_ELEMENTS_HEADER
28 #define OOMPH_GENERALISED_NEWTONIAN_REFINEABLE_AXISYMMETRIC_NAVIER_STOKES_ELEMENTS_HEADER
32 #include <oomph-lib-config.h>
36 #include "../generic/refineable_quad_element.h"
37 #include "../generic/error_estimator.h"
89 unsigned num_entries = DIM + ((DIM * DIM) - DIM) / 2;
90 if (flux.size() < num_entries)
92 std::ostringstream error_message;
93 error_message <<
"The flux vector has the wrong number of entries, "
94 << flux.size() <<
", whereas it should be " << num_entries
97 OOMPH_CURRENT_FUNCTION,
98 OOMPH_EXCEPTION_LOCATION);
111 for (
unsigned i = 0;
i < DIM;
i++)
113 flux[icount] = strainrate(
i,
i);
118 for (
unsigned i = 0;
i < DIM;
i++)
120 for (
unsigned j =
i + 1; j < DIM; j++)
122 flux[icount] = strainrate(
i, j);
140 cast_father_element_pt =
dynamic_cast<
149 this->
Re_pt = cast_father_element_pt->
re_pt();
157 this->
G_pt = cast_father_element_pt->
g_pt();
182 unsigned n_node = this->
nnode();
197 unsigned n_u_dof = 0;
198 for (
unsigned l = 0; l < n_node; l++)
200 unsigned n_master = 1;
209 n_master = hang_info_pt->
nmaster();
218 for (
unsigned m = 0; m < n_master; m++)
242 du_ddata.resize(n_u_dof, 0.0);
243 global_eqn_number.resize(n_u_dof, 0);
248 for (
unsigned l = 0; l < n_node; l++)
250 unsigned n_master = 1;
251 double hang_weight = 1.0;
260 n_master = hang_info_pt->
nmaster();
269 for (
unsigned m = 0; m < n_master; m++)
299 global_eqn_number[count] = global_eqn;
301 du_ddata[count] = psi[l] * hang_weight;
325 unsigned n_element = element_pt.size();
326 for (
unsigned e = 0;
e < n_element;
e++)
341 unsigned n_element = element_pt.size();
342 for (
unsigned e = 0;
e < n_element;
e++)
377 double*
const& parameter_pt,
384 OOMPH_CURRENT_FUNCTION,
385 OOMPH_EXCEPTION_LOCATION);
396 OOMPH_CURRENT_FUNCTION,
397 OOMPH_EXCEPTION_LOCATION);
421 unsigned n_node = this->
nnode();
424 for (
unsigned n = 0; n < n_node; n++)
434 unsigned n_node = this->
nnode();
437 for (
unsigned n = 0; n < n_node; n++)
444 for (
unsigned l = 0; l < n_pres; l++)
449 nod_pt->
unpin(p_index);
513 values.resize(DIM + 1, 0.0);
516 for (
unsigned i = 0;
i < DIM;
i++)
536 values.resize(DIM + 1);
539 for (
unsigned i = 0;
i < DIM + 1;
i++)
545 unsigned n_node =
nnode();
552 for (
unsigned i = 0;
i < DIM;
i++)
556 for (
unsigned l = 0; l < n_node; l++)
624 if (n_value ==
static_cast<int>(DIM))
627 unsigned total_index = 0;
629 unsigned NNODE_1D = 2;
634 for (
unsigned i = 0;
i < 2;
i++)
643 else if (
s[
i] == 1.0)
645 index[
i] = NNODE_1D - 1;
651 double float_index = 0.5 * (1.0 +
s[
i]) * (NNODE_1D - 1);
652 index[
i] = int(float_index);
655 double excess = float_index - index[
i];
666 index[
i] *
static_cast<unsigned>(pow(
static_cast<float>(NNODE_1D),
667 static_cast<int>(
i)));
706 return this->
nnode();
714 const int& n_value)
const
723 return this->
shape(s, psi);
740 GeneralisedNewtonianAxisymmetricQTaylorHoodElement>
757 FaceGeometry<GeneralisedNewtonianAxisymmetricQTaylorHoodElement>>
784 for (
unsigned l = 0; l < n_pres; l++)
809 using namespace QuadTreeNames;
818 double av_press = 0.0;
821 for (
unsigned ison = 0; ison < 4; ison++)
938 values.resize(DIM, 0.0);
941 for (
unsigned i = 0;
i < DIM;
i++)
966 for (
unsigned i = 0;
i < DIM;
i++)
972 unsigned n_node =
nnode();
979 for (
unsigned i = 0;
i < DIM;
i++)
983 for (
unsigned l = 0; l < n_node; l++)
1004 using namespace QuadTreeNames;
1023 else if (son_type ==
SE)
1029 else if (son_type ==
NE)
1036 else if (son_type ==
NW)
1044 cast_father_el_pt =
dynamic_cast<
1076 GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement>
1093 FaceGeometry<GeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement>>
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.
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
double value(const unsigned &i) const
Return i-th stored value. This function is not virtual so that it can be inlined. This means that if ...
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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.
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...
Data *& internal_data_pt(const unsigned &i)
Return a pointer to i-th internal data object.
A class for elements that solve the unsteady axisymmetric Navier–Stokes equations in cylindrical pola...
double *& re_st_pt()
Pointer to product of Reynolds and Strouhal number (=Womersley number)
virtual unsigned u_index_axi_nst(const unsigned &i) const
Return the index at which the i-th unknown velocity component is stored. The default value,...
void(*&)(const double &time, const Vector< double > &x, Vector< double > &f) axi_nst_body_force_fct_pt()
Access function for the body-force pointer.
double * ReInvRo_pt
Pointer to global Reynolds number x inverse Rossby number (used when in a rotating frame)
double *& viscosity_ratio_pt()
Pointer to Viscosity Ratio.
Vector< double > *& g_pt()
Pointer to Vector of gravitational components.
double *& density_ratio_pt()
Pointer to Density ratio.
double * Viscosity_Ratio_pt
Pointer to the viscosity ratio (relative to the viscosity used in the definition of the Reynolds numb...
double * ReInvFr_pt
Pointer to global Reynolds number x inverse Froude number (= Bond number / Capillary number)
double(*&)(const double &time, const Vector< double > &x) source_fct_pt()
Access function for the source-function pointer.
double * Re_pt
Pointer to global Reynolds number.
double * Density_Ratio_pt
Pointer to the density ratio (relative to the density used in the definition of the Reynolds number)
void(* Body_force_fct_pt)(const double &time, const Vector< double > &x, Vector< double > &result)
Pointer to body force function.
double *& re_invro_pt()
Pointer to global inverse Froude number.
double * ReSt_pt
Pointer to global Reynolds number x Strouhal number (=Womersley)
double(* Source_fct_pt)(const double &time, const Vector< double > &x)
Pointer to volumetric source function.
double *& re_invfr_pt()
Pointer to global inverse Froude number.
double interpolated_p_axi_nst(const Vector< double > &s) const
Return FE interpolated pressure at local coordinate s.
Vector< double > * G_pt
Pointer to global gravity Vector.
double *& re_pt()
Pointer to Reynolds number.
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when the time-derivatives are computed....
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: where (in that order)
void interpolated_u_axi_nst(const Vector< double > &s, Vector< double > &veloc) const
Compute vector of FE interpolated velocity u at local coordinate s.
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
unsigned npres_axi_nst() const
Return number of pressure values.
unsigned P_axi_nst_internal_index
Internal index that indicates at which internal data the pressure is stored.
/////////////////////////////////////////////////////////////////////////
static const unsigned Pconv[]
Static array of ints to hold conversion from pressure node numbers to actual node numbers.
virtual int p_nodal_index_axi_nst() const
Which nodal value represents the pressure?
unsigned npres_axi_nst() const
Return number of pressure values.
void pshape_axi_nst(const Vector< double > &s, Shape &psi) const
Pressure shape functions at local coordinate s.
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....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
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.
Refineable version of the Axisymmetric Navier–Stokes equations.
double geometric_jacobian(const Vector< double > &x)
Fill in the geometric Jacobian, which in this case is r.
virtual void pin_elemental_redundant_nodal_pressure_dofs()
Pin unused nodal pressure dofs (empty by default, because by default pressure dofs are not associated...
void fill_in_generic_residual_contribution_axi_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, unsigned flag)
Add element's contribution to the elemental residual vector and/or Jacobian matrix and mass matrix fl...
void dinterpolated_u_axi_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_dresidual_contribution_axi_nst(double *const ¶meter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, unsigned flag)
Add element's contribution to the derivative of the elemental residual vector and/or Jacobian matrix ...
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 ...
virtual Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node (Default: NULL, indicating that pressure is not based on nodal interp...
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Get 'flux' for Z2 error recovery: Upper triangular entries in strain rate tensor.
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_contribution_to_hessian_vector_products(Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Compute the hessian tensor vector products required to perform continuation of bifurcations analytica...
static void unpin_all_pressure_dofs(const Vector< GeneralisedElement * > &element_pt)
Unpin all pressure dofs in elements listed in vector.
unsigned num_Z2_flux_terms()
Number of 'flux' terms for Z2 error estimation.
virtual void unpin_elemental_pressure_dofs()=0
Unpin all pressure dofs in the element.
RefineableGeneralisedNewtonianAxisymmetricNavierStokesEquations()
Empty Constructor.
void further_build()
Further build: pass the pointers down to the sons.
Refineable version of Axisymmetric Quad Crouzeix Raviart elements (note that unlike the cartesian ver...
void further_build()
Further build for Crouzeix_Raviart interpolates the internal pressure dofs from father element: Make ...
void get_interpolated_values(const unsigned &t, const Vector< double > &s, Vector< double > &values)
Get all function values [u,v..,p] at previous timestep t (t=0: present; t>0: previous timestep)....
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: Reconstruct pressure from the (merged) sons.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 3 (velocities)
RefineableGeneralisedNewtonianAxisymmetricQCrouzeixRaviartElement()
Constructor:
unsigned nvertex_node() const
Number of vertex nodes in the element.
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
void get_interpolated_values(const Vector< double > &s, Vector< double > &values)
Get the function value u in Vector. Note: Given the generality of the interface (this function is usu...
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void unpin_elemental_pressure_dofs()
Unpin all the internal pressure freedoms.
Refineable version of Axisymmetric Quad Taylor Hood elements. (note that unlike the cartesian version...
unsigned required_nvalue(const unsigned &n) const
Number of values (pinned or dofs) required at node n.
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...
unsigned nvertex_node() const
Number of vertex nodes in the element.
void pin_elemental_redundant_nodal_pressure_dofs()
Unpin the proper nodal pressure dofs.
Node * interpolating_node_pt(const unsigned &n, const int &n_value)
The velocities are isoparametric and so the "nodes" interpolating the velocities are the geometric no...
Node * pressure_node_pt(const unsigned &n_p)
Pointer to n_p-th pressure node.
RefineableGeneralisedNewtonianAxisymmetricQTaylorHoodElement()
Constructor:
unsigned ninterpolating_node(const int &n_value)
The number of pressure nodes is 4. The number of velocity nodes is the same as the number of geometri...
Node * vertex_node_pt(const unsigned &j) const
Pointer to the j-th vertex node in the element.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: 4 (3 velocities + 1 pressure)
Node * get_interpolating_node_at_local_coordinate(const Vector< double > &s, const int &n_value)
The velocity nodes are the same as the geometric nodes. The pressure nodes must be calculated by usin...
void rebuild_from_sons(Mesh *&mesh_pt)
Rebuild from sons: empty.
void unpin_elemental_pressure_dofs()
Unpin all pressure dofs.
double local_one_d_fraction_of_interpolating_node(const unsigned &n1d, const unsigned &i, const int &n_value)
The pressure nodes are the corner nodes, so when n_value==DIM, the fraction is the same as the 1d nod...
unsigned ninterpolating_node_1d(const int &n_value)
The number of 1d pressure nodes is 2, the number of 1d velocity nodes is the same as the number of 1d...
void further_setup_hanging_nodes()
Perform additional hanging node procedures for variables that are not interpolated by all nodes....
unsigned nrecovery_order()
Order of recovery shape functions for Z2 error estimation: Same order as shape functions.
void interpolating_basis(const Vector< double > &s, Shape &psi, const int &n_value) const
The basis interpolating the pressure is given by pshape(). / The basis interpolating the velocity is ...
QuadTree * quadtree_pt()
Pointer to quadtree representation of this element.
void setup_hang_for_value(const int &value_id)
Internal helper function that is used to construct the hanging node schemes for the value_id-th inter...
A class that is used to template the refineable Q elements by dimension. It's really nothing more tha...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
RefineableElement * object_pt() const
Return the pointer to the object (RefineableElement) represented by the tree.
Tree * father_pt() const
Return pointer to father: NULL if it's a root node.
int son_type() const
Return son type.
Tree * son_pt(const int &son_index) const
Return pointer to the son for a given index. Note that to aid code readability specific enums have be...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...