using EntData = EntitiesFieldData::EntData;
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
-1;
auto u_exact = [](
const double x,
const double y,
const double) {
return x * x * y * y;
else
return cos(2 * x * M_PI) * cos(2 * y * M_PI);
};
else
-2 * M_PI * cos(2 * M_PI * y) * sin(2 * M_PI * x),
-2 * M_PI * cos(2 * M_PI * x) * sin(2 * M_PI * y)
};
};
auto source = [](
const double x,
const double y,
const double) {
return -(2 * x * x + 2 * y * y);
else
return 8 * M_PI * M_PI * cos(2 * x * M_PI) * cos(2 * y * M_PI);
};
static char help[] =
"...\n\n";
public:
private:
};
: domainField("U"), mField(m_field), oRder(4) {}
}
PETSC_NULL);
PETSC_NULL);
auto save_shared = [&](auto meshset, std::string prefix) {
auto file_name =
prefix + "_" +
1);
};
}
}
auto add_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
};
add_base_ops(pipeline_mng->getOpDomainLhsPipeline());
[](const double, const double, const double) { return 1; }));
pipeline_mng->getOpDomainRhsPipeline().push_back(
auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
add_base_ops(side_fe_ptr->getOpPtrVector());
side_fe_ptr->getOpPtrVector().push_back(
pipeline_mng->getOpSkeletonLhsPipeline().push_back(
pipeline_mng->getOpBoundaryLhsPipeline().push_back(
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
}
auto rule_lhs = [](
int,
int,
int p) ->
int {
return 2 *
p; };
auto rule_rhs = [](
int,
int,
int p) ->
int {
return 2 *
p; };
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
CHKERR pipeline_mng->setSkeletonLhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(rule_2);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule_2);
}
CHKERR KSPSetFromOptions(ksp_solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getSkeletonRhsFE().reset();
pipeline_mng->getSkeletonLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
auto add_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
};
auto u_vals_ptr = boost::make_shared<VectorDouble>();
auto u_grad_ptr = boost::make_shared<MatrixDouble>();
add_base_ops(pipeline_mng->getOpDomainRhsPipeline());
pipeline_mng->getOpDomainRhsPipeline().push_back(
pipeline_mng->getOpDomainRhsPipeline().push_back(
enum NORMS {
L2 = 0, SEMI_NORM,
H1, LAST_NORM };
std::array<int, LAST_NORM> error_indices;
for (
int i = 0;
i != LAST_NORM; ++
i)
if (const size_t nb_dofs = data.getIndices().size()) {
const int nb_integration_pts =
o->getGaussPts().size2();
auto t_w =
o->getFTensor0IntegrationWeight();
auto t_val = getFTensor0FromVec(*u_vals_ptr);
auto t_grad = getFTensor1FromMat<2>(*u_grad_ptr);
auto t_coords =
o->getFTensor1CoordsAtGaussPts();
std::array<double, LAST_NORM> error;
std::fill(error.begin(), error.end(), 0);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w *
o->getMeasure();
const double diff =
t_val -
u_exact(t_coords(0), t_coords(1), t_coords(2));
auto t_exact_grad =
u_grad_exact(t_coords(0), t_coords(1), t_coords(2));
const double diff_grad =
(t_grad(
i) - t_exact_grad(
i)) * (t_grad(
i) - t_exact_grad(
i));
error[SEMI_NORM] +=
alpha * diff_grad;
++t_w;
++t_val;
++t_grad;
++t_coords;
}
error[
H1] = error[
L2] + error[SEMI_NORM];
ADD_VALUES);
}
};
auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
add_base_ops(side_fe_ptr->getOpPtrVector());
side_fe_ptr->getOpPtrVector().push_back(
std::array<VectorDouble, 2> side_vals;
std::array<double, 2> area_map;
side_vals_op->doWorkRhsHook = [&](
DataOperator *op_ptr,
int side,
const auto nb_in_loop =
o->getFEMethod()->nInTheLoop;
area_map[nb_in_loop] =
o->getMeasure();
side_vals[nb_in_loop] = *u_vals_ptr;
if (!nb_in_loop) {
area_map[1] = 0;
side_vals[1].clear();
}
};
side_fe_ptr->getOpPtrVector().push_back(side_vals_op);
CHKERR o->loopSideFaces(
"dFE", *side_fe_ptr);
const auto in_the_loop = side_fe_ptr->nInTheLoop;
#ifndef NDEBUG
const std::array<std::string, 2> ele_type_name = {"BOUNDARY", "SKELETON"};
<< "do_work_rhs_error in_the_loop " << ele_type_name[in_the_loop];
#endif
const double s =
o->getMeasure() / (area_map[0] + area_map[1]);
constexpr std::array<int, 2> sign_array{1, -1};
std::array<double, LAST_NORM> error;
std::fill(error.begin(), error.end(), 0);
const int nb_integration_pts =
o->getGaussPts().size2();
if (!in_the_loop) {
side_vals[1].resize(nb_integration_pts, false);
auto t_coords =
o->getFTensor1CoordsAtGaussPts();
auto t_val_m = getFTensor0FromVec(side_vals[1]);
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
t_val_m =
u_exact(t_coords(0), t_coords(1), t_coords(2));
++t_coords;
++t_val_m;
}
}
auto t_val_p = getFTensor0FromVec(side_vals[0]);
auto t_val_m = getFTensor0FromVec(side_vals[1]);
auto t_w =
o->getFTensor0IntegrationWeight();
for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha =
o->getMeasure() * t_w;
const auto diff = t_val_p - t_val_m;
error[SEMI_NORM] +=
alpha *
p * diff * diff;
++t_w;
++t_val_p;
++t_val_m;
}
error[
H1] = error[SEMI_NORM];
ADD_VALUES);
};
skeleton_error_op->doWorkRhsHook = do_work_rhs_error;
boundary_error_op->doWorkRhsHook = do_work_rhs_error;
pipeline_mng->getOpDomainRhsPipeline().push_back(error_op);
pipeline_mng->getOpSkeletonLhsPipeline().push_back(skeleton_error_op);
pipeline_mng->getOpBoundaryRhsPipeline().push_back(boundary_error_op);
CHKERR pipeline_mng->loopFiniteElements();
CHKERR VecAssemblyBegin(l2_vec);
CHKERR VecAssemblyEnd(l2_vec);
const double *array;
CHKERR VecGetArrayRead(l2_vec, &array);
MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm L2 %6.4e",
MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm Energetic %6.4e",
std::sqrt(array[SEMI_NORM]));
MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm H1 %6.4e",
constexpr double eps = 1e-12;
if (std::sqrt(array[
H1]) >
eps)
}
CHKERR VecRestoreArrayRead(l2_vec, &array);
}
}
pipeline_mng->getSkeletonRhsFE().reset();
pipeline_mng->getSkeletonLhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
pipeline_mng->getBoundaryLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{{"U", u_ptr}},
{},
{},
{})
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile(
"out_result.h5m");
}
}
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
CHKERR poisson_problem.runProgram();
}
return 0;
}
#define MOFEM_LOG_C(channel, severity, format,...)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
ElementsAndOps< SPACE_DIM >::FaceSideEle FaceSideEle
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
auto createDMVector(DM dm)
Get smart vector from DM.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
DomainEle::UserDataOperator DomainEleOp
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< BASE_DIM, FIELD_DIM, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
BoundaryEle::UserDataOperator BoundaryEleOp
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
keeps basic data about problem
DofIdx getNbDofsRow() const
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator tp collect data from elements on the side of Edge/Face.
MoFEM::Interface & mField
MoFEMErrorCode boundaryCondition()
[Setup problem]
MoFEMErrorCode checkResults()
[Solve system]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode setIntegrationRules()
[Assemble system]
MoFEMErrorCode assembleSystem()
[Boundary condition]
Poisson2DiscontGalerkin(MoFEM::Interface &m_field)
MoFEMErrorCode solveSystem()
[Set integration rules]
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode outputResults()
[Output results]
MoFEMErrorCode runProgram()
[Output results]
Operator to evaluate Dirichlet boundary conditions using DG.
Operator the left hand side matrix.