#ifndef __POISSONOPERATORS_HPP__
#define __POISSONOPERATORS_HPP__
symm) {}
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
if (!nbRows)
if (!nbCols)
nbIntegrationPts = getGaussPts().size2();
if (row_side == col_side && row_type == col_type) {
isDiag = true;
} else {
isDiag = false;
}
CHKERR iNtegrate(row_data, col_data);
CHKERR aSsemble(row_data, col_data);
}
protected:
int nbRows;
int nbCols;
int nbIntegrationPts;
bool isDiag;
locMat.resize(nbRows, nbCols, false);
locMat.clear();
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
const double alpha = t_w * vol;
for (int rr = 0; rr != nbRows; rr++) {
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);
if (!isDiag && sYmm) {
locMat = trans(locMat);
&*locMat.data().begin(), ADD_VALUES);
}
}
};
template <
typename OPBASE>
struct OpBaseRhs :
public OPBASE {
if (!nbRows)
nbIntegrationPts = OPBASE::getGaussPts().size2();
}
protected:
int nbRows;
int nbIntegrationPts;
};
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;
locVec.resize(nbRows, false);
locVec.clear();
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
const double alpha =
vol * t_w * fSource(t_coords(NX), t_coords(NY), t_coords(NZ));
&*locVec.data().begin());
for (int rr = 0; rr != nbRows; rr++) {
t_a -= alpha * t_v;
++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),
assembleTranspose(assemble_transpose) {}
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
if (!nbRows)
if (!nbCols)
nbIntegrationPts = getGaussPts().size2();
CHKERR iNtegrate(row_data, col_data);
CHKERR aSsemble(row_data, col_data);
}
private:
int nbRows;
int nbCols;
int nbIntegrationPts;
const bool assembleTranspose;
locMat.resize(nbRows, nbCols, false);
locMat.clear();
const double area = getArea();
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
const double alpha = area * t_w;
&*locMat.data().begin());
for (int rr = 0; rr != nbRows; rr++) {
for (int cc = 0; cc != nbCols; cc++) {
c += alpha * t_row * t_col;
++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);
if (assembleTranspose) {
locMat = trans(locMat);
&*locMat.data().begin(), ADD_VALUES);
}
}
};
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:
FVal fValue;
const double bEta;
locVec.resize(nbRows, false);
locVec.clear();
const double area = getArea() * bEta;
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
double alpha =
area * t_w * fValue(t_coords(NX), t_coords(NY), t_coords(NZ));
&*locVec.data().begin());
for (int rr = 0; rr != nbRows; rr++) {
t_a += alpha * t_l;
++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)>
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) {}
if (!nbRows)
nbIntegrationPts = getGaussPts().size2();
}
private:
UVal uValue;
GVal gValue;
boost::shared_ptr<VectorDouble> fieldVals;
boost::shared_ptr<MatrixDouble> gradVals;
const double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_grad = getFTensor1FromMat<3>(*gradVals);
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
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);
++t_w;
++t_u;
++t_grad;
++t_coords;
}
}
}
};
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:
locMat.resize(nbRows, nbCols, false);
locMat.clear();
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_grad = getFTensor1FromMat<3>(*gradVals);
for (int gg = 0; gg != nbIntegrationPts; gg++) {
const double alpha = t_w * vol;
const double beta = alpha *
A(t_u);
t_gamma(
i) = (alpha * diffA(t_u)) * t_grad(
i);
&*locMat.data().begin());
for (int rr = 0; rr != nbRows; rr++) {
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<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:
locVec.resize(nbRows, false);
locVec.clear();
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
auto t_grad = getFTensor1FromMat<3>(*gradVals);
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
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);
&*locVec.data().begin());
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<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;
locVec.resize(nbRows, false);
locVec.clear();
const double area = getArea() * bEta;
auto t_w = getFTensor0IntegrationWeight();
auto t_coords = getFTensor1CoordsAtGaussPts();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
double alpha = area * t_w;
&*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;
locVec.resize(nbRows, false);
locVec.clear();
const double area = getArea() * bEta;
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nbIntegrationPts; gg++) {
double alpha = area * t_w;
&*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;
}
}
};
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(
domain_error->getOpPtrVector().push_back(
domain_error->getOpPtrVector().push_back(
new OpError(f_u, g_u, values_at_integration_ptr,
grad_at_integration_ptr, global_error));
}
boost::shared_ptr<PostProcFE> &post_proc_volume) const {
post_proc_volume =
boost::shared_ptr<PostProcFE>(
new PostProcFE(mField));
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(
post_proc_volume->getOpPtrVector().push_back(
post_proc_volume->getOpPtrVector().push_back(
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(
post_proc_volume->getOpPtrVector().push_back(
post_proc_volume->getOpPtrVector().push_back(
post_proc_volume->getOpPtrVector().push_back(
post_proc_volume->getPostProcMesh(),
post_proc_volume->getMapGaussPts(),
{{"U", u_ptr}, {"ERROR", e_ptr}},
{{"GRAD", grad_ptr}},
{},
{}
)
);
}
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 {
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(
domain_lhs_fe->getOpPtrVector().push_back(
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(
domain_rhs_fe->getOpPtrVector().push_back(
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(
boundary_rhs_fe->getOpPtrVector().push_back(
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 //__POISSONOPERATORS_HPP__