30 #ifndef OOMPH_SPACETIME_FLUID_TRACTION_ELEMENTS_HEADER
31 #define OOMPH_SPACETIME_FLUID_TRACTION_ELEMENTS_HEADER
35 #include <oomph-lib-config.h>
51 template<
class ELEMENT>
52 class NavierStokesSpaceTimeTractionElement
53 :
public virtual FaceGeometry<ELEMENT>,
54 public virtual FaceElement
59 const Vector<double>& x,
60 const Vector<double>& n,
61 Vector<double>& result);
69 const bool& called_from_refineable_constructor =
false)
79 if (!called_from_refineable_constructor)
82 if (element_pt->
dim() == 3)
95 std::ostringstream error_message_stream;
99 <<
"This flux element will not work correctly "
100 <<
"if nodes are hanging" << std::endl;
104 OOMPH_CURRENT_FUNCTION,
105 OOMPH_EXCEPTION_LOCATION);
152 residuals, jacobian, 1);
164 void output(std::ostream& outfile,
const unsigned& nplot)
175 const unsigned& flag);
183 const unsigned&
i)
const
226 for (
unsigned i = 0;
i <
Dim - 1;
i++)
236 (*Traction_fct_pt)(time, x, n, result);
253 template<
class ELEMENT>
256 Vector<double>& residuals,
257 DenseMatrix<double>& jacobian,
258 const unsigned& flag)
261 unsigned n_node = nnode();
270 unsigned n_intpt = integral_pt()->nweight();
276 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
279 double w = integral_pt()->weight(ipt);
283 double J = shape_and_test_at_knot(ipt, psif, testf);
289 double interpolated_t = 0.0;
292 Vector<double> interpolated_x(
Dim - 1, 0.0);
295 Vector<double> interpolated_spacetime_n(
Dim, 0.0);
298 Vector<double> interpolated_n(
Dim - 1, 0.0);
301 Vector<double> traction(
Dim - 1, 0.0);
304 for (
unsigned l = 0; l < n_node; l++)
307 interpolated_t += raw_nodal_position(l,
Dim - 1) * psif(l);
310 for (
unsigned i = 0;
i <
Dim - 1;
i++)
313 interpolated_x[
i] += nodal_position(l,
i) * psif(l);
318 outer_unit_normal(ipt, interpolated_spacetime_n);
323 if (std::fabs(interpolated_spacetime_n[
Dim - 1]) > 1.0e-15)
326 std::ostringstream error_message_stream;
330 <<
"The component of the outer unit normal in the "
331 <<
"time-direction\nshould be zero but the actual "
332 <<
"value is: " << interpolated_spacetime_n[
Dim - 1] << std::endl;
335 throw OomphLibError(error_message_stream.str(),
336 OOMPH_CURRENT_FUNCTION,
337 OOMPH_EXCEPTION_LOCATION);
342 for (
unsigned i = 0;
i <
Dim - 1;
i++)
345 interpolated_n[
i] = interpolated_spacetime_n[
i];
349 get_traction(interpolated_t, interpolated_x, interpolated_n, traction);
354 for (
unsigned l = 0; l < n_node; l++)
357 for (
unsigned i = 0;
i <
Dim - 1;
i++)
360 local_eqn = u_local_eqn(l,
i);
366 residuals[local_eqn] += traction[
i] * testf[l] *
W;
390 template<
class ELEMENT>
391 class RefineableNavierStokesSpaceTimeTractionElement
392 :
public virtual NavierStokesSpaceTimeTractionElement<ELEMENT>,
393 public virtual NonRefineableElementWithHangingNodes
432 residuals, jacobian, 1);
442 const unsigned& flag);
453 template<
class ELEMENT>
458 const unsigned& flag)
461 unsigned u_nodal_index[this->Dim - 1];
464 for (
unsigned i = 0;
i < this->Dim - 1;
i++)
468 dynamic_cast<ELEMENT*
>(this->bulk_element_pt())->u_index_nst(
i);
472 unsigned n_node = nnode();
481 unsigned n_intpt = integral_pt()->nweight();
487 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
490 double w = integral_pt()->weight(ipt);
494 double J = this->shape_and_test_at_knot(ipt, psif, testf);
500 double interpolated_t = 0.0;
503 Vector<double> interpolated_x(this->Dim - 1, 0.0);
506 Vector<double> interpolated_spacetime_n(this->Dim, 0.0);
509 Vector<double> interpolated_n(this->Dim - 1, 0.0);
512 Vector<double> traction(this->Dim - 1, 0.0);
515 for (
unsigned l = 0; l < n_node; l++)
518 interpolated_t += raw_nodal_position(l, this->Dim - 1) * psif(l);
521 for (
unsigned i = 0;
i < this->Dim - 1;
i++)
524 interpolated_x[
i] += this->nodal_position(l,
i) * psif(l);
529 this->outer_unit_normal(ipt, interpolated_spacetime_n);
534 if (std::fabs(interpolated_n[this->Dim - 1]) > 1.0e-15)
537 std::ostringstream error_message_stream;
540 error_message_stream <<
"The component of the outer unit normal in the "
541 <<
"time-direction\nshould be zero but the actual "
542 <<
"value is: " << interpolated_n[this->Dim - 1]
546 throw OomphLibError(error_message_stream.str(),
547 OOMPH_CURRENT_FUNCTION,
548 OOMPH_EXCEPTION_LOCATION);
553 for (
unsigned i = 0;
i < this->Dim - 1;
i++)
556 interpolated_n[
i] = interpolated_spacetime_n[
i];
561 interpolated_t, interpolated_x, interpolated_n, traction);
566 unsigned n_master = 1;
569 double hang_weight = 1.0;
572 HangInfo* hang_info_pt = 0;
576 for (
unsigned l = 0; l < n_node; l++)
579 bool is_node_hanging = this->node_pt(l)->is_hanging();
585 hang_info_pt = this->node_pt(l)->hanging_pt();
588 n_master = hang_info_pt->nmaster();
598 for (
unsigned m = 0; m < n_master; m++)
601 for (
unsigned i = 0;
i < this->Dim - 1;
i++)
609 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
613 hang_weight = hang_info_pt->master_weight(m);
619 local_eqn = this->nodal_local_eqn(l, u_nodal_index[
i]);
629 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...
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 void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
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...
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,...
unsigned Dim
The highest dimension of the problem.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns just the residuals.
NavierStokesSpaceTimeTractionElement(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.
void output(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,[z],u,v,[w],p in tecplot format.
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_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the jacobian.
void output(std::ostream &outfile)
Overload the output function.
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 ...
~NavierStokesSpaceTimeTractionElement()
Destructor should not delete anything.
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...
void(* Traction_fct_pt)(const double &time, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Pointer to an imposed traction function.
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(*&)(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 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_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the Jacobian.
RefineableNavierStokesSpaceTimeTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the face index.
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...
~RefineableNavierStokesSpaceTimeTractionElement()
Destructor should not delete anything.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns just the residuals.
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 ...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...