29 #ifndef OOMPH_MIXED_ORDER_PETROV_GALERKIN_SPACE_TIME_FLUID_TRACTION_ELEMENTS_HEADER
30 #define OOMPH_MIXED_ORDER_PETROV_GALERKIN_SPACE_TIME_FLUID_TRACTION_ELEMENTS_HEADER
34 #include <oomph-lib-config.h>
50 template<
class ELEMENT>
68 const bool& called_from_refineable_constructor =
false)
78 if (!called_from_refineable_constructor)
81 if (element_pt->
dim() == 3)
94 std::ostringstream error_message_stream;
98 <<
"This flux element will not work correctly "
99 <<
"if nodes are hanging" << std::endl;
103 OOMPH_CURRENT_FUNCTION,
104 OOMPH_EXCEPTION_LOCATION);
151 residuals, jacobian, 1);
163 void output(std::ostream& outfile,
const unsigned& nplot)
174 const unsigned& flag);
182 const unsigned&
i)
const
208 const unsigned el_dim =
dim();
215 "Can only deal with 2D space-time traction elements!",
216 OOMPH_CURRENT_FUNCTION,
217 OOMPH_EXCEPTION_LOCATION);
224 for (
unsigned i = 0;
i < el_dim;
i++)
249 for (
unsigned j = 0; j <
nnode_1d; j++)
255 psi[index] = psi_values[0][
i] * psi_values[1][j];
282 for (
unsigned j = 0; j <
nnode_1d; j++)
288 test[index] = test_values[0][
i] * test_values[1][j];
311 for (
unsigned i = 0;
i <
Dim - 1;
i++)
321 (*Traction_fct_pt)(time, x, n, result);
338 template<
class ELEMENT>
343 const unsigned& flag)
346 unsigned n_node = nnode();
355 unsigned n_intpt = integral_pt()->nweight();
361 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
364 double w = integral_pt()->weight(ipt);
368 double J = shape_and_test_at_knot(ipt, psif, testf);
374 double interpolated_t = 0.0;
389 for (
unsigned l = 0; l < n_node; l++)
392 interpolated_t += raw_nodal_position(l,
Dim - 1) * psif(l);
395 for (
unsigned i = 0;
i <
Dim - 1;
i++)
398 interpolated_x[
i] += nodal_position(l,
i) * psif(l);
403 outer_unit_normal(ipt, interpolated_spacetime_n);
408 if (std::fabs(interpolated_spacetime_n[
Dim - 1]) > 1.0e-15)
411 std::ostringstream error_message_stream;
415 <<
"The component of the outer unit normal in the "
416 <<
"time-direction\nshould be zero but the actual "
417 <<
"value is: " << interpolated_spacetime_n[
Dim - 1] << std::endl;
421 OOMPH_CURRENT_FUNCTION,
422 OOMPH_EXCEPTION_LOCATION);
427 for (
unsigned i = 0;
i <
Dim - 1;
i++)
430 interpolated_n[
i] = interpolated_spacetime_n[
i];
434 get_traction(interpolated_t, interpolated_x, interpolated_n, traction);
439 for (
unsigned l = 0; l < n_node; l++)
442 for (
unsigned i = 0;
i <
Dim - 1;
i++)
445 local_eqn = u_local_eqn(l,
i);
451 residuals[local_eqn] += traction[
i] * testf[l] *
W;
475 template<
class ELEMENT>
517 residuals, jacobian, 1);
527 const unsigned& flag);
538 template<
class ELEMENT>
543 const unsigned& flag)
546 unsigned u_nodal_index[this->
Dim - 1];
549 for (
unsigned i = 0;
i < this->
Dim - 1;
i++)
553 dynamic_cast<ELEMENT*
>(this->bulk_element_pt())->u_index_nst(
i);
557 unsigned n_node = nnode();
566 unsigned n_intpt = integral_pt()->nweight();
572 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
575 double w = integral_pt()->weight(ipt);
579 double J = this->shape_and_test_at_knot(ipt, psif, testf);
585 double interpolated_t = 0.0;
600 for (
unsigned l = 0; l < n_node; l++)
603 interpolated_t += raw_nodal_position(l, this->Dim - 1) * psif(l);
606 for (
unsigned i = 0;
i < this->Dim - 1;
i++)
609 interpolated_x[
i] += this->nodal_position(l,
i) * psif(l);
614 this->outer_unit_normal(ipt, interpolated_spacetime_n);
619 if (std::fabs(interpolated_n[this->Dim - 1]) > 1.0e-15)
622 std::ostringstream error_message_stream;
625 error_message_stream <<
"The component of the outer unit normal in the "
626 <<
"time-direction\nshould be zero but the actual "
627 <<
"value is: " << interpolated_n[this->Dim - 1]
632 OOMPH_CURRENT_FUNCTION,
633 OOMPH_EXCEPTION_LOCATION);
638 for (
unsigned i = 0;
i < this->Dim - 1;
i++)
641 interpolated_n[
i] = interpolated_spacetime_n[
i];
646 interpolated_t, interpolated_x, interpolated_n, traction);
651 unsigned n_master = 1;
654 double hang_weight = 1.0;
661 for (
unsigned l = 0; l < n_node; l++)
664 bool is_node_hanging = this->node_pt(l)->is_hanging();
670 hang_info_pt = this->node_pt(l)->hanging_pt();
673 n_master = hang_info_pt->
nmaster();
683 for (
unsigned m = 0; m < n_master; m++)
686 for (
unsigned i = 0;
i < this->Dim - 1;
i++)
704 local_eqn = this->nodal_local_eqn(l, u_nodal_index[
i]);
714 residuals[local_eqn] += traction[
i] * testf[l] *
W * hang_weight;
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)
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
In a FaceElement, the "global" intrinsic coordinate of the element along the boundary,...
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
double J_eulerian_at_knot(const unsigned &ipt) const
Return the Jacobian of the mapping from local to global coordinates at the ipt-th integration point O...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
A general Finite Element class.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement....
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
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.
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
A class for elements that allow the imposition of an applied traction to the Navier–Stokes equations ...
virtual int u_local_eqn(const unsigned &n, const unsigned &i)
Access function that returns the local equation numbers for velocity components. u_local_eqn(n,...
NavierStokesMixedOrderSpaceTimeTractionElement(FiniteElement *const &element_pt, const int &face_index, const bool &called_from_refineable_constructor=false)
Constructor, which takes a "bulk" element and the value of the index and its limit.
double shape_and_test_at_knot(const unsigned &ipt, Shape &psi, Shape &test) const
Function to compute the shape and test functions and to return the Jacobian of mapping.
void fill_in_generic_residual_contribution_fluid_traction(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
This function returns the residuals for the traction function. flag=1(or 0): do (or don't) compute th...
~NavierStokesMixedOrderSpaceTimeTractionElement()
Destructor should not delete anything.
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
void(*&)(const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &result) traction_fct_pt()
Access function for the imposed traction pointer.
unsigned Dim
The highest dimension of the problem.
void output(std::ostream &outfile)
Overload the output function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns just the residuals.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the jacobian.
void get_traction(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Function to calculate the traction applied to the fluid.
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Pointer to an imposed traction function.
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
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 ...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns just the residuals.
RefineableNavierStokesMixedOrderSpaceTimeTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the face index.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the Jacobian.
void refineable_fill_in_generic_residual_contribution_fluid_traction(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
This function returns the residuals for the traction function. flag=1(or 0): do (or don't) compute th...
~RefineableNavierStokesMixedOrderSpaceTimeTractionElement()
Destructor should not delete anything.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values are the same as those in the bulk element.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
unsigned Dim
Dimension of zeta tuples (set by get_dim_helper) – needed because we store the scalar coordinates in ...
void shape< 3 >(const double &s, double *Psi)
1D shape functions specialised to quadratic order (3 Nodes)
void shape< 3 >(const double &s, double *Psi)
1D shape functions specialised to quadratic order (3 Nodes)
void shape< 3 >(const double &s, double *Psi)
1D shape functions specialised to quadratic order (3 Nodes)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...