v0.14.0
PoissonOperators.hpp
/**
* \file PoissonOperators.hpp
* \example PoissonOperators.hpp
*
*/
#ifndef __POISSONOPERATORS_HPP__
#define __POISSONOPERATORS_HPP__
#include <MoFEM.hpp>
using namespace MoFEM;
namespace PoissonExample {
/**
* \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.
*
*/
OpK(bool symm = true)
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,
// 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,
// 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,
// 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,
// 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)>
FSource;
OpF(FSource f_source)
: OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator>("U"),
fSource(f_source) {}
protected:
FSource fSource;
VectorDouble locVec;
/**
* \brief Integrate local entity vector
* @param data entity data on element row
* @return error code
*/
// 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
*/
// 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
*
*/
OpC(const bool assemble_transpose)
false),
assembleTranspose(assemble_transpose) {}
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
EntityType col_type,
// 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,
// 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,
// 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;
if (assembleTranspose) {
// 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)>
FVal;
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
*/
// 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
*/
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)>
UVal;
typedef boost::function<FTensor::Tensor1<double, 3>(
const double, const double, const double)>
GVal;
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,
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
*/
// 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
*/
// 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,
// 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
*/
// 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
*/
// 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
*/
// 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
*/
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 elements element instances
domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
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>(
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
*/
MoFEMErrorCode creatFEToPostProcessResults(
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));
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,
bool trans = true) const {
// Create elements element instances
domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
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__
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
OpK
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradSymTensorGrad< BASE_DIM, SPACE_DIM, SPACE_DIM, 0 > OpK
[Define entities]
Definition: elastic.cpp:39
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
OpError
Definition: initial_diffusion.cpp:61
FTensor::Tensor1< double, 3 >
OpK
Definition: simple_elasticity.cpp:16
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
PoissonExample
Definition: AuxPoissonFunctions.hpp:12
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1254
MoFEM.hpp
VolRule
Set integration rule to volume elements.
Definition: simple_interface.cpp:88
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
OPBASE
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PoissonExample::PostProcFE
PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore > PostProcFE
Definition: PoissonOperators.hpp:15
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MoFEM::OpSetHOInvJacToScalarBases
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:73
double
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:105
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
MoFEM::PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:670
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MoFEM::OpCalculateHOJac< 3 >
Definition: HODataOperators.hpp:269
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1269
FaceRule
Set integration rule to boundary elements.
Definition: simple_interface.cpp:91
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
FTensor::Tensor0
Definition: Tensor0.hpp:16
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3249
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::ForcesAndSourcesCore::RuleHookFun
boost::function< int(int order_row, int order_col, int order_data)> RuleHookFun
Definition: ForcesAndSourcesCore.hpp:28
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698