26 #ifndef OOMPH_RAVIART_THOMAS_DARCY_HEADER
27 #define OOMPH_RAVIART_THOMAS_DARCY_HEADER
31 #include <oomph-lib-config.h>
34 #include "../generic/elements.h"
35 #include "../generic/shape.h"
36 #include "../generic/error_estimator.h"
37 #include "../generic/projection.h"
46 template<
unsigned DIM>
93 for (
unsigned i = 0;
i < n;
i++)
101 (*Source_fct_pt)(x, b);
117 (*Mass_source_fct_pt)(x, b);
148 virtual double q_edge(
const unsigned& n)
const = 0;
152 const unsigned& j)
const = 0;
166 virtual void set_q_edge(
const unsigned& n,
const double& value) = 0;
185 Shape& q_basis)
const = 0;
190 Shape& div_q_basis_ds)
const = 0;
195 const unsigned n_node = this->
nnode();
196 Shape psi(n_node, DIM);
197 const unsigned n_q_basis = this->
nq_basis();
198 Shape q_basis_local(n_q_basis, DIM);
216 const unsigned& edge,
const unsigned& n)
const = 0;
221 const unsigned& edge,
const unsigned& n,
Vector<double>& x)
const = 0;
230 virtual double p_value(
const unsigned& n)
const = 0;
242 virtual void set_p_value(
const unsigned& n,
const double& value) = 0;
253 const Shape& q_basis_local,
255 Shape& q_basis)
const;
277 Shape q_basis(n_q_basis, DIM);
280 for (
unsigned i = 0;
i < DIM;
i++)
283 for (
unsigned l = 0; l < n_q_basis_edge; l++)
287 for (
unsigned l = n_q_basis_edge; l < n_q_basis; l++)
300 Shape q_basis(n_q_basis, DIM);
304 for (
unsigned l = 0; l < n_q_basis_edge; l++)
306 q_i +=
q_edge(l) * q_basis(l,
i);
308 for (
unsigned l = n_q_basis_edge; l < n_q_basis; l++)
310 q_i +=
q_internal(l - n_q_basis_edge) * q_basis(l,
i);
323 unsigned n_node =
nnode();
324 const unsigned n_q_basis =
nq_basis();
328 Shape div_q_basis_ds(n_q_basis);
348 for (
unsigned l = 0; l < n_q_basis_edge; l++)
350 div_q += 1.0 / det * div_q_basis_ds(l) *
q_edge(l);
355 for (
unsigned l = n_q_basis_edge; l < n_q_basis; l++)
357 div_q += 1.0 / det * div_q_basis_ds(l) *
q_internal(l - n_q_basis_edge);
381 Shape p_basis(n_p_basis);
390 for (
unsigned l = 0; l < n_p_basis; l++)
431 void output(std::ostream& outfile,
const unsigned& nplot);
437 const unsigned& nplot,
444 const unsigned& nplot,
480 Shape& div_q_basis_ds,
481 Shape& div_q_test_ds)
const = 0;
492 Shape& div_q_basis_ds,
493 Shape& div_q_test_ds)
const = 0;
516 template<
class DARCY_ELEMENT>
534 std::stringstream error_stream;
535 error_stream <<
"Darcy elements only store 2 fields so fld = " << fld
538 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
548 Data* dat_pt = this->p_data_pt();
549 unsigned nvalue = dat_pt->
nvalue();
550 for (
unsigned i = 0;
i < nvalue;
i++)
552 data_values.push_back(std::make_pair(dat_pt,
i));
558 unsigned n = edge_dat_pt.size();
559 for (
unsigned j = 0; j < n; j++)
561 unsigned nvalue = edge_dat_pt[j]->nvalue();
562 for (
unsigned i = 0;
i < nvalue;
i++)
564 data_values.push_back(std::make_pair(edge_dat_pt[j],
i));
567 if (this->nq_basis_internal() > 0)
569 Data* int_dat_pt = this->q_internal_data_pt();
570 unsigned nvalue = int_dat_pt->
nvalue();
571 for (
unsigned i = 0;
i < nvalue;
i++)
573 data_values.push_back(std::make_pair(int_dat_pt,
i));
595 std::stringstream error_stream;
596 error_stream <<
"Darcy elements only store two fields so fld = " << fld
599 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
621 std::stringstream error_stream;
622 error_stream <<
"Darcy elements only store two fields so fld = " << fld
625 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
632 const unsigned n_dim = this->
dim();
633 const unsigned n_node = this->
nnode();
634 const unsigned n_q_basis = this->nq_basis();
635 const unsigned n_p_basis = this->np_basis();
639 Shape psi_geom(n_node), q_basis(n_q_basis, n_dim),
640 q_test(n_q_basis, n_dim), p_basis(n_p_basis), p_test(n_p_basis),
641 div_q_basis_ds(n_q_basis), div_q_test_ds(n_q_basis);
642 double J = this->shape_basis_test_local(
s,
653 unsigned n = p_basis.nindex1();
654 for (
unsigned i = 0;
i < n;
i++)
662 unsigned n = q_basis.
nindex1();
663 unsigned m = q_basis.nindex2();
664 for (
unsigned i = 0;
i < n;
i++)
666 for (
unsigned j = 0; j < m; j++)
668 psi(
i, j) = q_basis(
i, j);
686 std::stringstream error_stream;
687 error_stream <<
"Darcy elements only store two fields so fld = " << fld
690 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
694 double return_value = 0.0;
699 this->interpolated_p(
s, return_value);
704 OOMPH_CURRENT_FUNCTION,
705 OOMPH_EXCEPTION_LOCATION);
718 std::stringstream error_stream;
719 error_stream <<
"Darcy elements only store two fields so fld = " << fld
722 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
726 unsigned return_value = 0;
729 return_value = this->np_basis();
733 return_value = this->nq_basis();
746 std::stringstream error_stream;
747 error_stream <<
"Darcy elements only store two fields so fld = " << fld
750 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
754 int return_value = 0;
758 return_value = this->p_local_eqn(j);
762 unsigned nedge = this->nq_basis_edge();
765 return_value = this->q_edge_local_eqn(j);
769 return_value = this->q_internal_local_eqn(j - nedge);
777 void output(std::ostream& outfile,
const unsigned& nplot)
785 return std::max(this->Initial_Nvalue[n],
unsigned(1));
795 for (
unsigned j = 0; j < 3; j++)
806 const unsigned& flag)
811 unsigned n_dim = this->
dim();
820 const unsigned n_node = this->
nnode();
831 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
834 for (
unsigned i = 0;
i < n_dim;
i++)
850 int local_eqn = 0, local_unknown = 0;
858 if (solid_el_pt == 0)
860 throw OomphLibError(
"Trying to project Lagrangian coordinates in "
861 "non-SolidElement\n",
862 OOMPH_CURRENT_FUNCTION,
863 OOMPH_EXCEPTION_LOCATION);
867 Shape psi(n_node, n_position_type);
882 double interpolated_xi_bar =
887 for (
unsigned l = 0; l < n_node; ++l)
890 for (
unsigned k = 0; k < n_position_type; ++k)
900 residuals[local_eqn] +=
901 (interpolated_xi_proj - interpolated_xi_bar) * psi(l, k) *
907 for (
unsigned l2 = 0; l2 < n_node; l2++)
910 for (
unsigned k2 = 0; k2 < n_position_type; k2++)
914 if (local_unknown >= 0)
916 jacobian(local_eqn, local_unknown) +=
917 psi(l2, k2) * psi(l, k) *
W;
933 Shape psi(n_node, n_position_type);
944 double interpolated_x_proj = 0.0;
946 if (solid_el_pt != 0)
948 interpolated_x_proj =
954 interpolated_x_proj = this->
get_field(0, fld,
s);
962 for (
unsigned l = 0; l < n_node; ++l)
965 for (
unsigned k = 0; k < n_position_type; ++k)
968 if (solid_el_pt != 0)
982 if (n_position_type != 1)
985 "positions not in SolidElement\n",
986 OOMPH_CURRENT_FUNCTION,
987 OOMPH_EXCEPTION_LOCATION);
996 residuals[local_eqn] +=
997 (interpolated_x_proj - interpolated_x_bar) * psi(l, k) *
W;
1002 for (
unsigned l2 = 0; l2 < n_node; l2++)
1005 for (
unsigned k2 = 0; k2 < n_position_type; k2++)
1008 if (solid_el_pt != 0)
1018 if (local_unknown >= 0)
1020 jacobian(local_eqn, local_unknown) +=
1021 psi(l2, k2) * psi(l, k) *
W;
1050 double interpolated_value_proj = this->
get_field(now, fld,
s);
1053 double interpolated_value_bar =
1057 for (
unsigned l = 0; l < n_value; l++)
1063 residuals[local_eqn] +=
1064 (interpolated_value_proj - interpolated_value_bar) *
1070 for (
unsigned l2 = 0; l2 < n_value; l2++)
1073 if (local_unknown >= 0)
1075 jacobian(local_eqn, local_unknown) +=
1076 psi[l2] * psi[l] *
W;
1087 Shape psi(n_value, n_dim);
1098 this->interpolated_q(
s, q_proj);
1102 cast_el_pt->interpolated_q(other_s, q_bar);
1107 std::stringstream error_stream;
1108 error_stream <<
"Darcy elements are steady!\n";
1110 OOMPH_CURRENT_FUNCTION,
1111 OOMPH_EXCEPTION_LOCATION);
1116 for (
unsigned l = 0; l < n_value; l++)
1122 for (
unsigned i = 0;
i < n_dim;
i++)
1125 residuals[local_eqn] +=
1126 (q_proj[
i] - q_bar[
i]) * psi(l,
i) *
W;
1131 for (
unsigned l2 = 0; l2 < n_value; l2++)
1134 if (local_unknown >= 0)
1136 jacobian(local_eqn, local_unknown) +=
1137 psi(l2,
i) * psi(l,
i) *
W;
1148 "Wrong field specified. This should never happen\n",
1149 OOMPH_CURRENT_FUNCTION,
1150 OOMPH_EXCEPTION_LOCATION);
1157 throw OomphLibError(
"Wrong type specificied in Projection_type. "
1158 "This should never happen\n",
1159 OOMPH_CURRENT_FUNCTION,
1160 OOMPH_EXCEPTION_LOCATION);
1174 template<
class ELEMENT>
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Class implementing the generic maths of the Darcy equations using Raviart-Thomas elements with both e...
virtual unsigned required_nvalue(const unsigned &n) const =0
Number of values required at node n.
MassSourceFctPt mass_source_fct_pt() const
Access function: Pointer to mass source function (const version)
MassSourceFctPt Mass_source_fct_pt
Pointer to the mass source function.
virtual unsigned nq_basis_edge() const =0
Return the number of edge basis functions for q.
virtual int q_internal_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th internal degree of freedom.
double interpolated_div_q(const Vector< double > &s)
Calculate the FE representation of div q and return it.
virtual void edge_flux_interpolation_point_global(const unsigned &j, Vector< double > &x) const =0
Compute the global coordinates of the flux_interpolation point associated with the j-th edge-based q ...
virtual unsigned np_basis() const =0
Return the total number of pressure basis functions.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output FE representation of exact soln: x,y,q1,q2,div_q,p at Nplot^DIM plot points.
virtual void scale_basis(Shape &basis) const =0
Scale the edge basis to allow arbitrary edge mappings.
virtual unsigned q_edge_index(const unsigned &n) const =0
Return the nodal index at which the nth edge unknown is stored.
void output_with_projected_flux(std::ostream &outfile, const unsigned &nplot, const Vector< double > unit_normal)
Output incl. projection of fluxes into direction of the specified unit vector.
double interpolated_q(const Vector< double > &s, const unsigned i) const
Calculate the FE representation of the i-th component of q.
virtual void fill_in_generic_residual_contribution(Vector< double > &residuals, DenseMatrix< double > &jacobian, bool flag)
Fill in residuals and, if flag==true, jacobian.
virtual int p_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th pressure degree of freedom.
virtual double shape_basis_test_local(const Vector< double > &s, Shape &psi, Shape &q_basis, Shape &q_test, Shape &p_basis, Shape &p_test, Shape &div_q_basis_ds, Shape &div_q_test_ds) const =0
Returns the geometric basis, and the q, p and divergence basis functions and test functions at local ...
SourceFctPt Source_fct_pt
Pointer to body force function.
virtual double shape_basis_test_local_at_knot(const unsigned &ipt, Shape &psi, Shape &q_basis, Shape &q_test, Shape &p_basis, Shape &p_test, Shape &div_q_basis_ds, Shape &div_q_test_ds) const =0
Returns the geometric basis, and the q, p and divergence basis functions and test functions at integr...
virtual void edge_flux_interpolation_point_global(const unsigned &edge, const unsigned &n, Vector< double > &x) const =0
Returns the global coordinates of the nth flux interpolation point along an edge.
void get_q_basis(const Vector< double > &s, Shape &q_basis) const
Returns the transformed basis at local coordinate s.
virtual unsigned nedge_flux_interpolation_point() const =0
Returns the number of flux interpolation points along each edge of the element.
virtual int q_edge_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th edge (flux) degree of freedom.
double interpolated_p(const Vector< double > &s) const
Calculate the FE representation of p and return it.
DarcyEquations()
Constructor.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Z2 flux (use actual flux)
void interpolated_q(const Vector< double > &s, Vector< double > &q) const
Calculate the FE representation of q.
SourceFctPt source_fct_pt() const
Access function: Pointer to body force function (const version)
virtual void set_q_edge(const unsigned &n, const double &value)=0
Set the values of the edge (flux) degree of freedom.
virtual void face_local_coordinate_of_flux_interpolation_point(const unsigned &edge, const unsigned &n, Vector< double > &s) const =0
Compute the face element coordinates of the nth flux interpolation point along an edge.
virtual void pin_q_internal_value(const unsigned &n)=0
Pin the nth internal q value.
virtual unsigned face_index_of_q_edge_basis_fct(const unsigned &j) const =0
Return the face index associated with edge flux degree of freedom.
virtual double p_value(const unsigned &n) const =0
Return the nth pressure value.
virtual void pin_superfluous_darcy_dofs()
Helper function to pin superfluous dofs (empty; can be overloaded in projectable elements where we in...
double transform_basis(const Vector< double > &s, const Shape &q_basis_local, Shape &psi, Shape &q_basis) const
Performs a div-conserving transformation of the vector basis functions from the reference element to ...
virtual void set_p_value(const unsigned &n, const double &value)=0
Set the nth pressure value.
void interpolated_p(const Vector< double > &s, double &p) const
Calculate the FE representation of p.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in contribution to residuals for the Darcy equations.
unsigned num_Z2_flux_terms()
Number off flux terms for Z2 error estimator (use actual flux)
virtual unsigned q_internal_index() const =0
Return the index of the internal data where the q_internal degrees of freedom are stored.
virtual unsigned q_edge_node_number(const unsigned &n) const =0
Return the number of the node where the nth edge unknown is stored.
virtual void pin_p_value(const unsigned &n)=0
Pin the nth pressure value.
virtual Data * q_internal_data_pt() const =0
Return pointer to the Data object that stores the internal flux values.
virtual void set_q_internal(const unsigned &n, const double &value)=0
Set the values of the internal degree of freedom.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
Compute the error between the FE solution and the exact solution using the H(div) norm for q and L^2 ...
virtual void get_p_basis(const Vector< double > &s, Shape &p_basis) const =0
Return the pressure basis.
virtual double q_edge(const unsigned &n) const =0
Return the values of the n-th edge (flux) degree of freedom.
void interpolated_div_q(const Vector< double > &s, double &div_q) const
Calculate the FE representation of div q.
void(* MassSourceFctPt)(const Vector< double > &x, double &f)
Mass source function pointer typedef.
unsigned nq_basis() const
Return the total number of computational basis functions for q.
SourceFctPt & source_fct_pt()
Access function: Pointer to body force function.
MassSourceFctPt & mass_source_fct_pt()
Access function: Pointer to mass source function.
void output(std::ostream &outfile)
Output with default number of plot points.
void source(const Vector< double > &x, Vector< double > &b) const
Indirect access to the source function - returns 0 if no source function has been set.
virtual Vector< Data * > q_edge_data_pt() const =0
Return vector of pointers to the Data objects that store the edge flux values.
virtual unsigned nq_basis_internal() const =0
Return the number of internal basis functions for q.
virtual void get_div_q_basis_local(const Vector< double > &s, Shape &div_q_basis_ds) const =0
Returns the local form of the q basis and dbasis/ds at local coordinate s.
void(* SourceFctPt)(const Vector< double > &x, Vector< double > &f)
Source function pointer typedef.
unsigned self_test()
Self test – empty for now.
virtual Data * p_data_pt() const =0
Return pointer to the Data object that stores the pressure values.
void mass_source(const Vector< double > &x, double &b) const
Indirect access to the mass source function - returns 0 if no mass source function has been set.
virtual unsigned face_index_of_edge(const unsigned &j) const =0
Return the face index associated with specified edge.
virtual void get_q_basis_local(const Vector< double > &s, Shape &q_basis) const =0
Returns the local form of the q basis at local coordinate s.
virtual Vector< double > edge_flux_interpolation_point(const unsigned &edge, const unsigned &n) const =0
Returns the local coordinate of nth flux interpolation point along an edge.
virtual double q_internal(const unsigned &n) const =0
Return the values of the internal degree of freedom.
A class that represents a collection of data; each Data object may contain many different individual ...
void pin(const unsigned &i)
Pin the i-th stored variable.
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
A general Finite Element class.
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
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...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
unsigned nnode() const
Return the number of nodes.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
virtual double local_to_eulerian_mapping(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates, given the derivatives of the shape function...
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
An OomphLibError object which should be thrown when an run-time error is encountered....
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
void output(std::ostream &outfile, const unsigned &nplot)
Output FE representation of soln as in underlying element.
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
void pin_superfluous_darcy_dofs()
Helper function to pin superfluous dofs; required because we introduce at least one dof per node to a...
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!...
void residual_for_projection(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Residual for the projection step. Flag indicates if we want the Jacobian (1) or not (0)....
ProjectableDarcyElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
unsigned required_nvalue(const unsigned &n) const
Overload required_nvalue to return at least one value.
unsigned nfields_for_projection()
Number of fields to be projected: 2 (pressure and flux)
Template-free Base class for projectable elements.
unsigned Projected_lagrangian
When projecting the Lagrangain coordinates indicate which coordiante is to be projected.
Projection_Type Projection_type
Variable to indicate if we're projecting the history values of the nodal coordinates (Coordinate) the...
unsigned Time_level_for_projection
Time level we are projecting (0=current values; >0: history values)
unsigned Projected_field
Field that is currently being projected.
unsigned Projected_coordinate
When projecting the history values of the nodal coordinates, this is the coordinate we're projecting.
virtual double get_field(const unsigned &time, const unsigned &fld, const Vector< double > &s)=0
Return the fld-th field at local coordinates s at time-level time (time=0: current value; time>0: his...
Wrapper class for projectable elements. Adds "projectability" to the underlying ELEMENT.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
unsigned nindex1() const
Return the range of index 1 of the shape function object.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
int position_local_eqn(const unsigned &n, const unsigned &k, const unsigned &j) const
Access function that returns the local equation number that corresponds to the j-th coordinate of the...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
void output()
Doc the command line arguments.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...