v0.14.0
Loading...
Searching...
No Matches
PoissonOperators.hpp
/**
* \file PoissonOperators.hpp
* \example PoissonOperators.hpp
*
*/
#ifndef __POISSONOPERATORS_HPP__
#define __POISSONOPERATORS_HPP__
namespace PoissonExample {
using PostProcFE = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
/**
* \brief Calculate the grad-grad operator and assemble matrix
*
* Calculate
* \f[
* \mathbf{K}=\int_\Omega \nabla \boldsymbol\phi \cdot \nabla \boldsymbol\phi
* \textrm{d}\Omega \f] and assemble to global matrix.
*
* This operator is executed on element for each unique combination of entities.
*
*/
struct OpK : public VolumeElementForcesAndSourcesCore::UserDataOperator {
OpK(bool symm = true)
: VolumeElementForcesAndSourcesCore::UserDataOperator("U", "U", OPROWCOL,
symm) {}
/**
* \brief Do calculations for give operator
* @param row_side row side number (local number) of entity on element
* @param col_side column side number (local number) of entity on element
* @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
* @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
* @param row_data data for row
* @param col_data data for column
* @return error code
*/
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
// get number of dofs on row
nbRows = row_data.getIndices().size();
// if no dofs on row, exit that work, nothing to do here
if (!nbRows)
// get number of dofs on column
nbCols = col_data.getIndices().size();
// if no dofs on Columbia, exit nothing to do here
if (!nbCols)
// get number of integration points
nbIntegrationPts = getGaussPts().size2();
// check if entity block is on matrix diagonal
if (row_side == col_side && row_type == col_type) {
isDiag = true; // yes, it is
} else {
isDiag = false;
}
// integrate local matrix for entity block
CHKERR iNtegrate(row_data, col_data);
// assemble local matrix
CHKERR aSsemble(row_data, col_data);
}
protected:
///< error code
int nbRows; ///< number of dofs on rows
int nbCols; ///< number if dof on column
int nbIntegrationPts; ///< number of integration points
bool isDiag; ///< true if this block is on diagonal
FTensor::Index<'i', 3> i; ///< summit Index
MatrixDouble locMat; ///< local entity block matrix
/**
* \brief Integrate grad-grad operator
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
// set size of local entity bock
locMat.resize(nbRows, nbCols, false);
// clear matrix
locMat.clear();
// get element volume
double vol = getVolume();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get base function gradient on rows
auto t_row_grad = row_data.getFTensor1DiffN<3>();
// loop over integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
// take into account Jacobean
const double alpha = t_w * vol;
// noalias(locMat) +=
// alpha*prod(row_data.getDiffN(gg),trans(col_data.getDiffN(gg))); take
// fist element to local matrix
FTensor::Tensor0<double *> a(&*locMat.data().begin());
// loop over rows base functions
for (int rr = 0; rr != nbRows; rr++) {
// get column base functions gradient at gauss point gg
auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
// loop over columns
for (int cc = 0; cc != nbCols; cc++) {
// calculate element of local matrix
a += alpha * (t_row_grad(i) * t_col_grad(i));
++t_col_grad; // move to another gradient of base function on column
++a; // move to another element of local matrix in column
}
++t_row_grad; // move to another element of gradient of base function on
// row
}
++t_w; // move to another integration weight
}
}
/**
* \brief Assemble local entity block matrix
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
// get pointer to first global index on row
const int *row_indices = &*row_data.getIndices().data().begin();
// get pointer to first global index on column
const int *col_indices = &*col_data.getIndices().data().begin();
Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
: getFEMethod()->snes_B;
// assemble local matrix
CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
&*locMat.data().begin(), ADD_VALUES);
if (!isDiag && sYmm) {
// if not diagonal term and since global matrix is symmetric assemble
// transpose term.
locMat = trans(locMat);
CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
&*locMat.data().begin(), ADD_VALUES);
}
}
};
/**
* \brief template class for integration oh the right hand side
*/
template <typename OPBASE> struct OpBaseRhs : public OPBASE {
OpBaseRhs(const std::string field_name) : OPBASE(field_name, OPBASE::OPROW) {}
/**
* \brief This function is called by finite element
*
* Do work is composed from two operations, integrate and assembly. Also,
* it set values nbRows, and nbIntegrationPts.
*
*/
MoFEMErrorCode doWork(int row_side, EntityType row_type,
EntitiesFieldData::EntData &row_data) {
// get number of dofs on row
nbRows = row_data.getIndices().size();
if (!nbRows)
// get number of integration points
nbIntegrationPts = OPBASE::getGaussPts().size2();
// integrate local vector
CHKERR iNtegrate(row_data);
// assemble local vector
CHKERR aSsemble(row_data);
}
/**
* \brief Class dedicated to integrate operator
* @param data entity data on element row
* @return error code
*/
virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) = 0;
/**
* \brief Class dedicated to assemble operator to global system vector
* @param data entity data (indices, base functions, etc. ) on element row
* @return error code
*/
virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data) = 0;
protected:
///< error code
int nbRows; ///< number of dofs on row
int nbIntegrationPts; ///< number of integration points
};
/**
* \brief Operator calculate source term,
*
* \f[
* \mathbf{F} = \int_\Omega \boldsymbol\phi f \textrm{d}\Omega
* \f]
*
*/
struct OpF
: public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
typedef boost::function<double(const double, const double, const double)>
OpF(FSource f_source)
: OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator>("U"),
fSource(f_source) {}
protected:
VectorDouble locVec;
/**
* \brief Integrate local entity vector
* @param data entity data on element row
* @return error code
*/
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
// set size of local vector
locVec.resize(nbRows, false);
// clear local entity vector
locVec.clear();
// get finite element volume
double vol = getVolume();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get base functions on entity
auto t_v = data.getFTensor0N();
// get coordinates at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// loop over all integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
// evaluate constant term
const double alpha =
vol * t_w * fSource(t_coords(NX), t_coords(NY), t_coords(NZ));
// get element of local vector
&*locVec.data().begin());
// loop over base functions
for (int rr = 0; rr != nbRows; rr++) {
// add to local vector source term
t_a -= alpha * t_v;
++t_a; // move to next element of local vector
++t_v; // move to next base function
}
++t_w; // move to next integration weight
++t_coords; // move to next physical coordinates at integration point
}
}
/**
* \brief assemble local entity vector to the global right hand side
* @param data entity data, i.e. global indices of local vector
* @return error code
*/
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data) {
// get global indices of local vector
const int *indices = &*data.getIndices().data().begin();
// get values from local vector
const double *vals = &*locVec.data().begin();
Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
: getFEMethod()->snes_f;
// assemble vector
CHKERR VecSetValues(f, nbRows, indices, vals, ADD_VALUES);
}
};
/**
* \brief Calculate constrains matrix
*
* \f[
* \mathbf{C} = \int_{\partial\Omega} \boldsymbol\psi \boldsymbol\phi
* \textrm{d}\partial\Omega \f] where \f$\lambda \f$ is base function on
* boundary
*
*/
struct OpC : public FaceElementForcesAndSourcesCore::UserDataOperator {
OpC(const bool assemble_transpose)
: FaceElementForcesAndSourcesCore::UserDataOperator("L", "U", OPROWCOL,
false),
assembleTranspose(assemble_transpose) {}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
// get number of dofs on row
nbRows = row_data.getIndices().size();
// exit here if no dofs on row, nothing to do
if (!nbRows)
// get number of dofs on column,
nbCols = col_data.getIndices().size();
// exit here if no dofs on roe, nothing to do
if (!nbCols)
// get number of integration points
nbIntegrationPts = getGaussPts().size2();
// integrate local constrains matrix
CHKERR iNtegrate(row_data, col_data);
// assemble local constrains matrix
CHKERR aSsemble(row_data, col_data);
}
private:
///< error code
int nbRows; ///< number of dofs on row
int nbCols; ///< number of dofs on column
int nbIntegrationPts; ///< number of integration points
const bool assembleTranspose; ///< assemble transpose, i.e. CT if set to true
MatrixDouble locMat; ///< local constrains matrix
/** \brief Integrate local constrains matrix
*/
inline MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
// set size of local constrains matrix
locMat.resize(nbRows, nbCols, false);
// clear matrix
locMat.clear();
// get area of element
const double area = getArea();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get base functions on entity
auto t_row = row_data.getFTensor0N();
// run over integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
const double alpha = area * t_w;
// get element of local matrix
&*locMat.data().begin());
// run over base functions on rows
for (int rr = 0; rr != nbRows; rr++) {
// get first base functions on column for integration point gg
auto t_col = col_data.getFTensor0N(gg, 0);
// run over base function on column
for (int cc = 0; cc != nbCols; cc++) {
// integrate element of constrains matrix
c += alpha * t_row * t_col;
++t_col; // move to next base function on column
++c; // move to next element of local matrix
}
++t_row; // move to next base function on row
}
++t_w; // move to next integrate weight
}
}
/**
* \brief integrate local constrains matrix
*/
inline MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
// get indices on row
const int *row_indices = &*row_data.getIndices().data().begin();
// get indices on column
const int *col_indices = &*col_data.getIndices().data().begin();
Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
: getFEMethod()->snes_B;
// assemble local matrix
CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
&*locMat.data().begin(), ADD_VALUES);
// cerr << locMat << endl;
// assemble transpose of local matrix
locMat = trans(locMat);
CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
&*locMat.data().begin(), ADD_VALUES);
}
}
};
/**
* \brief Assemble constrains vector
*
* \f[
* \mathbf{g} = \int_{\partial\Omega} \boldsymbol\psi \overline{u}
* \textrm{d}\partial\Omega \f]
*
*/
struct Op_g
: public OpBaseRhs<FaceElementForcesAndSourcesCore::UserDataOperator> {
typedef boost::function<double(const double, const double, const double)>
Op_g(FVal f_value, const string field_name = "L", const double beta = 1)
: OpBaseRhs<FaceElementForcesAndSourcesCore::UserDataOperator>(
fValue(f_value), bEta(beta) {}
protected:
FTensor::Number<0> NX; ///< x-direction index
FTensor::Number<1> NY; ///< y-direction index
FTensor::Number<2> NZ; ///< z-direction index
FVal fValue; ///< Function pointer evaluating values of "U" at the boundary
VectorDouble locVec;
const double bEta;
/**
* \brief Integrate local constrains vector
*/
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
// set size to local vector
locVec.resize(nbRows, false);
// clear local vector
locVec.clear();
// get face area
const double area = getArea() * bEta;
// get integration weight
auto t_w = getFTensor0IntegrationWeight();
// get base function
auto t_l = data.getFTensor0N();
// get coordinates at integration point
auto t_coords = getFTensor1CoordsAtGaussPts();
// make loop over integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
// evaluate function on boundary and scale it by area and integration
// weight
double alpha =
area * t_w * fValue(t_coords(NX), t_coords(NY), t_coords(NZ));
// get element of vector
&*locVec.data().begin());
//
for (int rr = 0; rr != nbRows; rr++) {
t_a += alpha * t_l;
++t_a;
++t_l;
}
++t_w;
++t_coords;
}
}
/**
* \brief assemble constrains vectors
*/
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data) {
const int *indices = &*data.getIndices().data().begin();
const double *vals = &*locVec.data().begin();
Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
: getFEMethod()->snes_f;
CHKERR VecSetValues(f, nbRows, indices, &*vals, ADD_VALUES);
}
};
/**
* \brief Evaluate error
*/
struct OpError
: public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
typedef boost::function<double(const double, const double, const double)>
typedef boost::function<FTensor::Tensor1<double, 3>(
const double, const double, const double)>
OpError(UVal u_value, GVal g_value,
boost::shared_ptr<VectorDouble> &field_vals,
boost::shared_ptr<MatrixDouble> &grad_vals, Vec global_error)
: OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator>("ERROR"),
globalError(global_error), uValue(u_value), gValue(g_value),
fieldVals(field_vals), gradVals(grad_vals) {}
MoFEMErrorCode doWork(int row_side, EntityType row_type,
EntitiesFieldData::EntData &row_data) {
nbRows = row_data.getFieldData().size();
if (!nbRows)
nbIntegrationPts = getGaussPts().size2();
CHKERR iNtegrate(row_data);
CHKERR aSsemble(row_data);
}
private:
Vec globalError; ///< ghost vector with global (integrated over volume) error
UVal uValue; ///< function with exact solution
GVal gValue; ///< function with exact solution for gradient
boost::shared_ptr<VectorDouble> fieldVals;
boost::shared_ptr<MatrixDouble> gradVals;
/**
* \brief Integrate error
*/
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
// clear field dofs
data.getFieldData().clear();
// get volume of element
const double vol = getVolume();
// get integration weight
auto t_w = getFTensor0IntegrationWeight();
// get solution at integration point
auto t_u = getFTensor0FromVec(*fieldVals);
// get solution at integration point
auto t_grad = getFTensor1FromMat<3>(*gradVals);
// get coordinates at integration point
auto t_coords = getFTensor1CoordsAtGaussPts();
// keep exact gradient and error or gradient
FTensor::Tensor1<double, 3> t_exact_grad, t_error_grad;
// integrate over
for (int gg = 0; gg != nbIntegrationPts; gg++) {
double alpha = vol * t_w;
// evaluate exact value
double exact_u = uValue(t_coords(NX), t_coords(NY), t_coords(NZ));
// evaluate exact gradient
t_exact_grad = gValue(t_coords(NX), t_coords(NY), t_coords(NZ));
// calculate gradient errro
t_error_grad(i) = t_grad(i) - t_exact_grad(i);
// error
double error = pow(t_u - exact_u, 2) + t_error_grad(i) * t_error_grad(i);
// iterate over base functions
data.getFieldData()[0] += alpha * error;
++t_w; // move to next integration point
++t_u; // next value of function at integration point
++t_grad; // next gradient at integration point
++t_coords; // next coordinate at integration point
}
}
/**
* \brief Assemble error
*/
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data) {
// set error on mesh
data.getFieldDofs()[0]->getFieldData() = sqrt(data.getFieldData()[0]);
// assemble vector to global error
CHKERR VecSetValue(globalError, 0, data.getFieldData()[0], ADD_VALUES);
}
};
struct OpKt : public OpK {
OpKt(boost::function<double(const double)> a,
boost::function<double(const double)> diff_a,
boost::shared_ptr<VectorDouble> &field_vals,
boost::shared_ptr<MatrixDouble> &grad_vals)
: OpK(false), A(a), diffA(diff_a), fieldVals(field_vals),
gradVals(grad_vals) {}
protected:
/**
* \brief Integrate grad-grad operator
* @param row_data row data (consist base functions on row entity)
* @param col_data column data (consist base functions on column entity)
* @return error code
*/
inline MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
// set size of local entity bock
locMat.resize(nbRows, nbCols, false);
// clear matrix
locMat.clear();
// get element volume
double vol = getVolume();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get solution at integration point
auto t_u = getFTensor0FromVec(*fieldVals);
// get solution at integration point
auto t_grad = getFTensor1FromMat<3>(*gradVals);
// get base function gradient on rows
auto t_row_grad = row_data.getFTensor1DiffN<3>();
// loop over integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
// take into account Jacobian
const double alpha = t_w * vol;
const double beta = alpha * A(t_u);
t_gamma(i) = (alpha * diffA(t_u)) * t_grad(i);
// take fist element to local matrix
&*locMat.data().begin());
// loop over rows base functions
for (int rr = 0; rr != nbRows; rr++) {
// get column base function
auto t_col = col_data.getFTensor0N(gg, 0);
// get column base functions gradient at gauss point gg
auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
// loop over columns
for (int cc = 0; cc != nbCols; cc++) {
// calculate element of local matrix
a += (t_row_grad(i) * beta) * t_col_grad(i) +
t_row_grad(i) * (t_gamma(i) * t_col);
++t_col; // move to next base function
++t_col_grad; // move to another gradient of base function on column
++a; // move to another element of local matrix in column
}
++t_row_grad; // move to another element of gradient of base function on
// row
}
++t_w; // move to another integration weight
++t_u; // move to next value at integration point
++t_grad; // move to next gradient value
}
}
boost::function<double(const double)> A;
boost::function<double(const double)> diffA;
boost::shared_ptr<VectorDouble> fieldVals;
boost::shared_ptr<MatrixDouble> gradVals;
};
struct OpResF_Domain : public OpF {
OpResF_Domain(FSource f_source, boost::function<double(const double)> a,
boost::shared_ptr<VectorDouble> &field_vals,
boost::shared_ptr<MatrixDouble> &grad_vals)
: OpF(f_source), A(a), fieldVals(field_vals), gradVals(grad_vals) {}
protected:
/**
* \brief Integrate local entity vector
* @param data entity data on element row
* @return error code
*/
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
// set size of local vector
locVec.resize(nbRows, false);
// clear local entity vector
locVec.clear();
// get finite element volume
double vol = getVolume();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get solution at integration point
auto t_u = getFTensor0FromVec(*fieldVals);
// get solution at integration point
auto t_grad = getFTensor1FromMat<3>(*gradVals);
// get base functions on entity
auto t_v = data.getFTensor0N();
// get base function gradient on rows
auto t_v_grad = data.getFTensor1DiffN<3>();
// get coordinates at integration points
auto t_coords = getFTensor1CoordsAtGaussPts();
// loop over all integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
// evaluate constant term
const double alpha = vol * t_w;
const double source_term =
alpha * fSource(t_coords(NX), t_coords(NY), t_coords(NZ));
grad_term(i) = (alpha * A(t_u)) * t_grad(i);
// get element of local vector
&*locVec.data().begin());
// loop over base functions
for (int rr = 0; rr != nbRows; rr++) {
// add to local vector source term
t_a += t_v_grad(i) * grad_term(i) + t_v * source_term;
++t_a; // move to next element of local vector
++t_v; // move to next base function
++t_v_grad; // move to next gradient of base function
}
++t_w; // move to next integration weights
++t_u; // move to next value
++t_grad; // move to next gradient value
++t_coords; // move to next physical coordinates at integration point
}
}
boost::function<double(const double)> A;
boost::shared_ptr<VectorDouble> fieldVals;
boost::shared_ptr<MatrixDouble> gradVals;
};
struct OpRes_g : public Op_g {
OpRes_g(FVal f_value, boost::shared_ptr<VectorDouble> &field_vals)
: Op_g(f_value, "L", 1), fieldVals(field_vals) {}
protected:
boost::shared_ptr<VectorDouble> fieldVals;
/**
* \brief Integrate local constrains vector
*/
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
// set size to local vector
locVec.resize(nbRows, false);
// clear local vector
locVec.clear();
// get face area
const double area = getArea() * bEta;
// get integration weight
auto t_w = getFTensor0IntegrationWeight();
// get base function
auto t_l = data.getFTensor0N();
// get solution at integration point
auto t_u = getFTensor0FromVec(*fieldVals);
// get coordinates at integration point
auto t_coords = getFTensor1CoordsAtGaussPts();
// make loop over integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
// evaluate function on boundary and scale it by area and integration
// weight
double alpha = area * t_w;
// get element of vector
&*locVec.data().begin());
for (int rr = 0; rr != nbRows; rr++) {
t_a += alpha * t_l *
(t_u - fValue(t_coords(NX), t_coords(NY), t_coords(NZ)));
++t_a;
++t_l;
}
++t_w;
++t_u;
++t_coords;
}
}
};
struct OpResF_Boundary : public Op_g {
OpResF_Boundary(boost::shared_ptr<VectorDouble> &lambda_vals)
: Op_g(FVal(), "U", 1), lambdaVals(lambda_vals) {}
protected:
boost::shared_ptr<VectorDouble> lambdaVals;
/**
* \brief Integrate local constrains vector
*/
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
// set size to local vector
locVec.resize(nbRows, false);
// clear local vector
locVec.clear();
// get face area
const double area = getArea() * bEta;
// get integration weight
auto t_w = getFTensor0IntegrationWeight();
// get base function
auto t_u = data.getFTensor0N();
// get solution at integration point
auto t_lambda = getFTensor0FromVec(*lambdaVals);
// make loop over integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
// evaluate function on boundary and scale it by area and integration
// weight
double alpha = area * t_w;
// get element of vector
&*locVec.data().begin());
for (int rr = 0; rr != nbRows; rr++) {
t_a += alpha * t_u * t_lambda;
++t_a;
++t_u;
}
++t_w;
++t_lambda;
}
}
};
/**
* \brief Set integration rule to volume elements
*
* This rule is used to integrate \f$\nabla v \cdot \nabla u\f$, thus
* if approximation field and testing field are polynomial order "p",
* then rule for exact integration is 2*(p-1).
*
* Integration rule is order of polynomial which is calculated exactly. Finite
* element selects integration method based on return of this function.
*
*/
struct VolRule {
int operator()(int, int, int p) const { return 2 * (p - 1); }
};
/**
* \brief Set integration rule to boundary elements
*
* This is uses to integrate values on the face. Is used to integrate
* \f$(\mathbf{n} \cdot \lambda) u\f$, where Lagrange multiplayer
* is order "p_row" and approximate function is order "p_col".
*
* Integration rule is order of polynomial which is calculated exactly. Finite
* element selects integration method based on return of this function.
*
*/
struct FaceRule {
int operator()(int p_row, int p_col, int p_data) const {
return 2 * p_data + 1;
}
};
/**
* \brief Create finite elements instances
*
* Create finite element instances and add operators to finite elements.
*
*/
struct CreateFiniteElements {
CreateFiniteElements(MoFEM::Interface &m_field) : mField(m_field) {}
/**
* \brief Create finite element to calculate matrix and vectors
*/
boost::function<double(const double, const double, const double)> f_u,
boost::function<double(const double, const double, const double)>
f_source,
boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
bool trans = true) const {
// Create elements element instances
domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
// Set integration rule to elements instances
domain_lhs_fe->getRuleHook = VolRule();
domain_rhs_fe->getRuleHook = VolRule();
boundary_lhs_fe->getRuleHook = FaceRule();
boundary_rhs_fe->getRuleHook = FaceRule();
// Add operators to element instances
// Add operator grad-grad for calculate matrix
domain_lhs_fe->getOpPtrVector().push_back(new OpK());
// Add operator to calculate source terms
domain_rhs_fe->getOpPtrVector().push_back(new OpF(f_source));
// Add operator calculating constrains matrix
boundary_lhs_fe->getOpPtrVector().push_back(new OpC(trans));
// Add operator calculating constrains vector
boundary_rhs_fe->getOpPtrVector().push_back(new Op_g(f_u));
}
/**
* \brief Create finite element to calculate error
*/
MoFEMErrorCode createFEToEvaluateError(
boost::function<double(const double, const double, const double)> f_u,
boost::function<FTensor::Tensor1<double, 3>(const double, const double,
const double)>
g_u,
Vec global_error,
boost::shared_ptr<ForcesAndSourcesCore> &domain_error) const {
// Create finite element instance to calculate error
domain_error = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
domain_error->getRuleHook = VolRule();
// Set integration rule
// Crate shared vector storing values of field "u" on integration points on
// element. element is local and is used to exchange data between operators.
boost::shared_ptr<VectorDouble> values_at_integration_ptr =
boost::make_shared<VectorDouble>();
// Storing gradients of field
boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
boost::make_shared<MatrixDouble>();
// Add default operator to calculate field values at integration points
domain_error->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
// Add default operator to calculate field gradient at integration points
domain_error->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
// Add operator to integrate error element by element.
domain_error->getOpPtrVector().push_back(
new OpError(f_u, g_u, values_at_integration_ptr,
grad_at_integration_ptr, global_error));
}
/**
* \brief Create finite element to post-process results
*/
boost::shared_ptr<PostProcFE> &post_proc_volume) const {
// Note that user can stack together arbitrary number of operators to
// compose complex PDEs.
// Post-process results. This is standard element, with functionality
// enabling refining mesh for post-processing. In addition in
// PostProcOnRefMesh.hpp are implanted set of users operators to
// post-processing fields. Here using simplified mechanism for
// post-processing finite element, we add operators to save data from field
// on mesh tags for ParaView visualization.
post_proc_volume =
boost::shared_ptr<PostProcFE>(new PostProcFE(mField));
// Add operators to the elements, starting with some generic
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_volume->getOpPtrVector().push_back(
new OpCalculateHOJac<3>(jac_ptr));
post_proc_volume->getOpPtrVector().push_back(
new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
post_proc_volume->getOpPtrVector().push_back(
new OpSetHOInvJacToScalarBases<3>(H1, inv_jac_ptr));
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
auto e_ptr = boost::make_shared<VectorDouble>();
post_proc_volume->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", u_ptr));
post_proc_volume->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<3>("U", grad_ptr));
post_proc_volume->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("ERROR", e_ptr));
using OpPPMap = OpPostProcMapInMoab<3, 3>;
post_proc_volume->getOpPtrVector().push_back(
new OpPPMap(
post_proc_volume->getPostProcMesh(),
post_proc_volume->getMapGaussPts(),
{{"U", u_ptr}, {"ERROR", e_ptr}},
{{"GRAD", grad_ptr}},
{},
{}
)
);
}
/**
* \brief Create finite element to calculate matrix and vectors
*/
MoFEMErrorCode createFEToAssembleMatrixAndVectorForNonlinearProblem(
boost::function<double(const double, const double, const double)> f_u,
boost::function<double(const double, const double, const double)>
f_source,
boost::function<double(const double)> a,
boost::function<double(const double)> diff_a,
boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
ForcesAndSourcesCore::RuleHookFun vol_rule,
ForcesAndSourcesCore::RuleHookFun face_rule = FaceRule(),
bool trans = true) const {
// Create elements element instances
domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
new VolumeElementForcesAndSourcesCore(mField));
boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
// Set integration rule to elements instances
domain_lhs_fe->getRuleHook = vol_rule;
domain_rhs_fe->getRuleHook = vol_rule;
boundary_lhs_fe->getRuleHook = face_rule;
boundary_rhs_fe->getRuleHook = face_rule;
// Set integration rule
// Crate shared vector storing values of field "u" on integration points on
// element. element is local and is used to exchange data between operators.
boost::shared_ptr<VectorDouble> values_at_integration_ptr =
boost::make_shared<VectorDouble>();
// Storing gradients of field
boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
boost::make_shared<MatrixDouble>();
// multipliers values
boost::shared_ptr<VectorDouble> multiplier_at_integration_ptr =
boost::make_shared<VectorDouble>();
// Add operators to element instances
// Add default operator to calculate field values at integration points
domain_lhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
// Add default operator to calculate field gradient at integration points
domain_lhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
// Add operator grad-(1+u^2)grad for calculate matrix
domain_lhs_fe->getOpPtrVector().push_back(new OpKt(
a, diff_a, values_at_integration_ptr, grad_at_integration_ptr));
// Add default operator to calculate field values at integration points
domain_rhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
// Add default operator to calculate field gradient at integration points
domain_rhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
// Add operator to calculate source terms
domain_rhs_fe->getOpPtrVector().push_back(new OpResF_Domain(
f_source, a, values_at_integration_ptr, grad_at_integration_ptr));
// Add operator calculating constrains matrix
boundary_lhs_fe->getOpPtrVector().push_back(new OpC(trans));
// Add default operator to calculate field values at integration points
boundary_rhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
// Add default operator to calculate values of Lagrange multipliers
boundary_rhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("L", multiplier_at_integration_ptr));
// Add operator calculating constrains vector
boundary_rhs_fe->getOpPtrVector().push_back(
new OpRes_g(f_u, values_at_integration_ptr));
boundary_rhs_fe->getOpPtrVector().push_back(
new OpResF_Boundary(multiplier_at_integration_ptr));
}
private:
};
} // namespace PoissonExample
#endif //__POISSONOPERATORS_HPP__
static Index< 'p', 3 > p
constexpr double a
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
const double c
speed of light (cm/ns)
PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore > PostProcFE
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
Set integration rule to boundary elements.
Deprecated interface functions.
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< PostProcFE > &post_proc_volume) const
Create finite element to post-process results.
MoFEMErrorCode createFEToEvaluateError(boost::function< double(const double, const double, const double)> f_u, boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> g_u, Vec global_error, boost::shared_ptr< ForcesAndSourcesCore > &domain_error) const
Create finite element to calculate error.
MoFEMErrorCode createFEToAssembleMatrixAndVector(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::shared_ptr< ForcesAndSourcesCore > &domain_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &domain_rhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_rhs_fe, bool trans=true) const
Create finite element to calculate matrix and vectors.
int operator()(int p_row, int p_col, int p_data) const
FTensor::Number< 2 > NZ
z-direction index
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
assemble constrains vectors
boost::function< double(const double, const double, const double)> FVal
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
FTensor::Number< 0 > NX
x-direction index
FTensor::Number< 1 > NY
y-direction index
FVal fValue
Function pointer evaluating values of "U" at the boundary.
int nbIntegrationPts
number of integration points
virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)=0
Class dedicated to assemble operator to global system vector.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
This function is called by finite element.
virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)=0
Class dedicated to integrate operator.
MatrixDouble locMat
local constrains matrix
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate local constrains matrix
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate local constrains matrix.
int nbIntegrationPts
number of integration points
int nbRows
< error code
const bool assembleTranspose
assemble transpose, i.e. CT if set to true
int nbCols
number of dofs on column
Vec globalError
ghost vector with global (integrated over volume) error
UVal uValue
function with exact solution
boost::shared_ptr< MatrixDouble > gradVals
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
boost::function< double(const double, const double, const double)> UVal
FTensor::Number< 2 > NZ
FTensor::Index< 'i', 3 > i
GVal gValue
function with exact solution for gradient
boost::shared_ptr< VectorDouble > fieldVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate error.
FTensor::Number< 1 > NY
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> GVal
FTensor::Number< 0 > NX
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
Assemble error.
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
assemble local entity vector to the global right hand side
FTensor::Number< 2 > NZ
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local entity vector.
boost::function< double(const double, const double, const double)> FSource
FTensor::Number< 1 > NY
FTensor::Number< 0 > NX
virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
int nbRows
< error code
int nbCols
number if dof on column
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations for give operator.
bool isDiag
true if this block is on diagonal
MatrixDouble locMat
local entity block matrix
int nbIntegrationPts
number of integration points
FTensor::Index< 'i', 3 > i
summit Index
virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble local entity block matrix.
boost::function< double(const double)> diffA
boost::shared_ptr< MatrixDouble > gradVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< VectorDouble > fieldVals
boost::function< double(const double)> A
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
boost::shared_ptr< VectorDouble > fieldVals
boost::shared_ptr< VectorDouble > lambdaVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
boost::shared_ptr< MatrixDouble > gradVals
boost::shared_ptr< VectorDouble > fieldVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local entity vector.
boost::function< double(const double)> A
FTensor::Index< 'i', 3 > i
int operator()(int, int, int p) const
Set integration rule.