26#ifndef OOMPH_RAVIART_THOMAS_DARCY_HEADER
27#define OOMPH_RAVIART_THOMAS_DARCY_HEADER
31#include <oomph-lib-config.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;
216 const unsigned&
edge,
const unsigned&
n)
const = 0;
230 virtual double p_value(
const unsigned&
n)
const = 0;
280 for (
unsigned i = 0;
i <
DIM;
i++)
437 const unsigned&
nplot,
444 const unsigned&
nplot,
516 template<
class DARCY_ELEMENT>
549 unsigned nvalue =
dat_pt->nvalue();
550 for (
unsigned i = 0;
i < nvalue;
i++)
559 for (
unsigned j = 0;
j <
n;
j++)
562 for (
unsigned i = 0;
i < nvalue;
i++)
567 if (this->nq_basis_internal() > 0)
571 for (
unsigned i = 0;
i < nvalue;
i++)
596 error_stream <<
"Darcy elements only store two fields so fld = " <<
fld
622 error_stream <<
"Darcy elements only store two fields so fld = " <<
fld
634 const unsigned n_q_basis = this->nq_basis();
635 const unsigned n_p_basis = this->np_basis();
642 double J = this->shape_basis_test_local(
s,
654 for (
unsigned i = 0;
i <
n;
i++)
664 for (
unsigned i = 0;
i <
n;
i++)
666 for (
unsigned j = 0;
j <
m;
j++)
687 error_stream <<
"Darcy elements only store two fields so fld = " <<
fld
719 error_stream <<
"Darcy elements only store two fields so fld = " <<
fld
747 error_stream <<
"Darcy elements only store two fields so fld = " <<
fld
762 unsigned nedge = this->nq_basis_edge();
785 return std::max(this->Initial_Nvalue[
n],
unsigned(1));
795 for (
unsigned j = 0;
j < 3;
j++)
806 const unsigned&
flag)
834 for (
unsigned i = 0;
i <
n_dim;
i++)
860 throw OomphLibError(
"Trying to project Lagrangian coordinates in "
861 "non-SolidElement\n",
985 "positions not in SolidElement\n",
1098 this->interpolated_q(
s,
q_proj);
1122 for (
unsigned i = 0;
i <
n_dim;
i++)
1148 "Wrong field specified. This should never happen\n",
1157 throw OomphLibError(
"Wrong type specificied in Projection_type. "
1158 "This should never happen\n",
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 Vector< Data * > q_edge_data_pt() const =0
Return vector of pointers to the Data objects that store the edge flux values.
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.
virtual Data * q_internal_data_pt() const =0
Return pointer to the Data object that stores the internal flux values.
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.
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.
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 Data * p_data_pt() const =0
Return pointer to the Data object that stores the pressure values.
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 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.
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 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.
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.
SourceFctPt & source_fct_pt()
Access function: Pointer to body force function.
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 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...
FaceGeometry class definition: This policy class is used to allow construction of face elements that ...
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.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational 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...
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 .
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...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
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....
Darcy upgraded to become projectable.
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...
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 ...
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.
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...
SolidFiniteElement class.
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).