#ifndef __POISSONOPERATORS_HPP__
#define __POISSONOPERATORS_HPP__
symm) {}
EntityType col_type,
nbRows = row_data.getIndices().size();
nbCols = col_data.getIndices().size();
if (row_side == col_side && row_type == col_type) {
} else {
}
}
protected:
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_grad = row_data.getFTensor1DiffN<3>();
const double alpha = t_w * vol;
for (
int rr = 0; rr !=
nbRows; rr++) {
auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
for (
int cc = 0; cc !=
nbCols; cc++) {
a +=
alpha * (t_row_grad(
i) * t_col_grad(
i));
++t_col_grad;
}
++t_row_grad;
}
++t_w;
}
}
const int *row_indices = &*row_data.getIndices().data().begin();
const int *col_indices = &*col_data.getIndices().data().begin();
Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
: getFEMethod()->snes_B;
&*
locMat.data().begin(), ADD_VALUES);
&*
locMat.data().begin(), ADD_VALUES);
}
}
};
template <
typename OPBASE>
struct OpBaseRhs :
public OPBASE {
OpBaseRhs(
const std::string field_name) :
OPBASE(field_name, OPBASE::OPROW) {}
nbRows = row_data.getIndices().size();
}
protected:
};
struct OpF
: public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
typedef boost::function<double(const double, const double, const double)>
: OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator>("U"),
protected:
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_v = data.getFTensor0N();
auto t_coords = getFTensor1CoordsAtGaussPts();
for (
int rr = 0; rr !=
nbRows; rr++) {
++t_a;
++t_v;
}
++t_w;
++t_coords;
}
}
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;
}
};
OpC(
const bool assemble_transpose)
false),
EntityType col_type,
nbRows = row_data.getIndices().size();
nbCols = col_data.getIndices().size();
}
private:
const double area = getArea();
auto t_w = getFTensor0IntegrationWeight();
auto t_row = row_data.getFTensor0N();
const double alpha = area * t_w;
for (
int rr = 0; rr !=
nbRows; rr++) {
auto t_col = col_data.getFTensor0N(gg, 0);
for (
int cc = 0; cc !=
nbCols; cc++) {
++t_col;
}
++t_row;
}
++t_w;
}
}
const int *row_indices = &*row_data.getIndices().data().begin();
const int *col_indices = &*col_data.getIndices().data().begin();
Mat
B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
: getFEMethod()->snes_B;
&*
locMat.data().begin(), ADD_VALUES);
&*
locMat.data().begin(), ADD_VALUES);
}
}
};
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>(
field_name),
protected:
const double area = getArea() *
bEta;
auto t_w = getFTensor0IntegrationWeight();
auto t_l = data.getFTensor0N();
auto t_coords = getFTensor1CoordsAtGaussPts();
area * t_w *
fValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
for (
int rr = 0; rr !=
nbRows; rr++) {
++t_a;
++t_l;
}
++t_w;
++t_coords;
}
}
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;
}
};
: 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)>
boost::shared_ptr<VectorDouble> &field_vals,
boost::shared_ptr<MatrixDouble> &grad_vals,
Vec global_error)
: OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator>("ERROR"),
nbRows = row_data.getFieldData().size();
}
private:
boost::shared_ptr<MatrixDouble>
gradVals;
data.getFieldData().clear();
const double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_grad = getFTensor1FromMat<3>(*
gradVals);
auto t_coords = getFTensor1CoordsAtGaussPts();
double alpha = vol * t_w;
double exact_u =
uValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
t_exact_grad =
gValue(t_coords(
NX), t_coords(
NY), t_coords(
NZ));
t_error_grad(
i) = t_grad(
i) - t_exact_grad(
i);
double error = pow(t_u - exact_u, 2) + t_error_grad(
i) * t_error_grad(
i);
data.getFieldData()[0] +=
alpha * error;
++t_w;
++t_u;
++t_grad;
++t_coords;
}
}
data.getFieldDofs()[0]->getFieldData() = sqrt(data.getFieldData()[0]);
}
};
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)
protected:
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_grad = getFTensor1FromMat<3>(*
gradVals);
auto t_row_grad = row_data.getFTensor1DiffN<3>();
const double alpha = t_w * vol;
const double beta =
alpha *
A(t_u);
for (
int rr = 0; rr !=
nbRows; rr++) {
auto t_col = col_data.getFTensor0N(gg, 0);
auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
for (
int cc = 0; cc !=
nbCols; cc++) {
a += (t_row_grad(
i) * beta) * t_col_grad(
i) +
t_row_grad(
i) * (t_gamma(
i) * t_col);
++t_col;
++t_col_grad;
}
++t_row_grad;
}
++t_w;
++t_u;
++t_grad;
}
}
boost::function<double(
const double)>
A;
boost::function<double(
const double)>
diffA;
boost::shared_ptr<MatrixDouble>
gradVals;
};
struct OpResF_Domain : public OpF {
boost::shared_ptr<VectorDouble> &field_vals,
boost::shared_ptr<MatrixDouble> &grad_vals)
protected:
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_grad = getFTensor1FromMat<3>(*
gradVals);
auto t_v = data.getFTensor0N();
auto t_v_grad = data.getFTensor1DiffN<3>();
auto t_coords = getFTensor1CoordsAtGaussPts();
const double alpha = vol * t_w;
const double source_term =
grad_term(
i) = (
alpha *
A(t_u)) * t_grad(
i);
for (
int rr = 0; rr !=
nbRows; rr++) {
t_a += t_v_grad(
i) * grad_term(
i) + t_v * source_term;
++t_a;
++t_v;
++t_v_grad;
}
++t_w;
++t_u;
++t_grad;
++t_coords;
}
}
boost::function<double(
const double)>
A;
boost::shared_ptr<MatrixDouble>
gradVals;
};
struct OpRes_g : public Op_g {
OpRes_g(
FVal f_value, boost::shared_ptr<VectorDouble> &field_vals)
protected:
const double area = getArea() *
bEta;
auto t_w = getFTensor0IntegrationWeight();
auto t_l = data.getFTensor0N();
auto t_coords = getFTensor1CoordsAtGaussPts();
double alpha = area * t_w;
for (
int rr = 0; rr !=
nbRows; rr++) {
++t_a;
++t_l;
}
++t_w;
++t_u;
++t_coords;
}
}
};
struct OpResF_Boundary : public Op_g {
protected:
const double area = getArea() *
bEta;
auto t_w = getFTensor0IntegrationWeight();
auto t_u = data.getFTensor0N();
double alpha = area * t_w;
for (
int rr = 0; rr !=
nbRows; rr++) {
t_a +=
alpha * t_u * t_lambda;
++t_a;
++t_u;
}
++t_w;
++t_lambda;
}
}
};
int operator()(
int,
int,
int p)
const {
return 2 * (
p - 1); }
};
int operator()(
int p_row,
int p_col,
int p_data)
const {
return 2 * p_data + 1;
}
};
struct CreateFiniteElements {
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 {
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>(
domain_lhs_fe->getRuleHook =
VolRule();
domain_rhs_fe->getRuleHook =
VolRule();
boundary_lhs_fe->getRuleHook =
FaceRule();
boundary_rhs_fe->getRuleHook =
FaceRule();
domain_lhs_fe->getOpPtrVector().push_back(
new OpK());
domain_rhs_fe->getOpPtrVector().push_back(new OpF(f_source));
boundary_lhs_fe->getOpPtrVector().push_back(new OpC(trans));
boundary_rhs_fe->getOpPtrVector().push_back(new Op_g(f_u));
}
boost::function<double(const double, const double, const double)> f_u,
const double)>
g_u,
boost::shared_ptr<ForcesAndSourcesCore> &domain_error) const {
domain_error = boost::shared_ptr<ForcesAndSourcesCore>(
domain_error->getRuleHook =
VolRule();
boost::shared_ptr<VectorDouble> values_at_integration_ptr =
boost::make_shared<VectorDouble>();
boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
boost::make_shared<MatrixDouble>();
domain_error->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
domain_error->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
domain_error->getOpPtrVector().push_back(
new OpError(f_u, g_u, values_at_integration_ptr,
grad_at_integration_ptr, global_error));
}
boost::shared_ptr<ForcesAndSourcesCore> &post_proc_volume) const {
post_proc_volume = boost::shared_ptr<ForcesAndSourcesCore>(
CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
post_proc_volume)
->generateReferenceElementMesh();
CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
post_proc_volume)
->addFieldValuesPostProc("U");
CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
post_proc_volume)
->addFieldValuesPostProc("ERROR");
CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
post_proc_volume)
->addFieldValuesGradientPostProc("U");
}
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 {
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>(
domain_lhs_fe->getRuleHook = vol_rule;
domain_rhs_fe->getRuleHook = vol_rule;
boundary_lhs_fe->getRuleHook = face_rule;
boundary_rhs_fe->getRuleHook = face_rule;
boost::shared_ptr<VectorDouble> values_at_integration_ptr =
boost::make_shared<VectorDouble>();
boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
boost::make_shared<MatrixDouble>();
boost::shared_ptr<VectorDouble> multiplier_at_integration_ptr =
boost::make_shared<VectorDouble>();
domain_lhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
domain_lhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
domain_lhs_fe->getOpPtrVector().push_back(new OpKt(
a, diff_a, values_at_integration_ptr, grad_at_integration_ptr));
domain_rhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
domain_rhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
domain_rhs_fe->getOpPtrVector().push_back(new OpResF_Domain(
f_source,
a, values_at_integration_ptr, grad_at_integration_ptr));
boundary_lhs_fe->getOpPtrVector().push_back(new OpC(trans));
boundary_rhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
boundary_rhs_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("L", multiplier_at_integration_ptr));
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:
};
}
#endif
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpK
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
constexpr auto VecSetValues
constexpr auto MatSetValues
DataForcesAndSourcesCore::EntData EntData
Set integration rule to boundary elements.
Deprecated interface functions.
MoFEM::Interface & mField
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< ForcesAndSourcesCore > &post_proc_volume) const
Create finite element to post-process results.
CreateFiniteElements(MoFEM::Interface &m_field)
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 finite element to calculate matrix and vectors.
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
Op_g(FVal f_value, const string field_name="L", const double beta=1)
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.
OpBaseRhs(const std::string field_name)
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
OpC(const bool assemble_transpose)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate local constrains matrix.
int nbIntegrationPts
number of integration points
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
OpError(UVal u_value, GVal g_value, boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &grad_vals, Vec global_error)
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::Index< 'i', 3 > i
GVal gValue
function with exact solution for gradient
boost::shared_ptr< VectorDouble > fieldVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate error.
boost::function< FTensor::Tensor1< double, 3 > const double, const double, const double)> GVal
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
Assemble error.
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
assemble local entity vector to the global right hand side
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local entity vector.
boost::function< double(const double, const double, const double)> FSource
virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
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
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)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
boost::shared_ptr< VectorDouble > fieldVals
OpRes_g(FVal f_value, boost::shared_ptr< VectorDouble > &field_vals)
boost::shared_ptr< VectorDouble > lambdaVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
OpResF_Boundary(boost::shared_ptr< VectorDouble > &lambda_vals)
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
OpResF_Domain(FSource f_source, boost::function< double(const double)> a, boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &grad_vals)
int operator()(int, int, int p) const