46 return cos(2 * x * M_PI) * cos(2 * y * M_PI);
55 -2 * M_PI * cos(2 * M_PI * y) * sin(2 * M_PI * x),
56 -2 * M_PI * cos(2 * M_PI * x) * sin(2 * M_PI * y)
63 return -(2 * x * x + 2 * y * y);
65 return 8 * M_PI * M_PI * cos(2 * x * M_PI) * cos(2 * y * M_PI);
71static char help[] =
"...\n\n";
101 : domainField(
"U"), mField(m_field), oRder(4) {}
147 auto save_shared = [&](
auto meshset, std::string prefix) {
178 auto add_base_ops = [&](
auto &pipeline) {
179 auto det_ptr = boost::make_shared<VectorDouble>();
180 auto jac_ptr = boost::make_shared<MatrixDouble>();
181 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
187 add_base_ops(pipeline_mng->getOpDomainLhsPipeline());
190 [](
const double,
const double,
const double) {
return 1; }));
191 pipeline_mng->getOpDomainRhsPipeline().push_back(
195 auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
196 add_base_ops(side_fe_ptr->getOpPtrVector());
197 side_fe_ptr->getOpPtrVector().push_back(
201 pipeline_mng->getOpSkeletonLhsPipeline().push_back(
205 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
207 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
218 auto rule_lhs = [](int, int,
int p) ->
int {
return 2 *
p; };
219 auto rule_rhs = [](int, int,
int p) ->
int {
return 2 *
p; };
220 auto rule_2 = [
this](int, int, int) {
return 2 *
oRder; };
223 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
224 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
226 CHKERR pipeline_mng->setSkeletonLhsIntegrationRule(rule_2);
227 CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
228 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(rule_2);
229 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule_2);
241 auto ksp_solver = pipeline_mng->
createKSP();
242 CHKERR KSPSetFromOptions(ksp_solver);
243 CHKERR KSPSetUp(ksp_solver);
253 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
254 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
266 pipeline_mng->getDomainLhsFE().reset();
267 pipeline_mng->getSkeletonRhsFE().reset();
268 pipeline_mng->getSkeletonLhsFE().reset();
269 pipeline_mng->getBoundaryRhsFE().reset();
270 pipeline_mng->getBoundaryLhsFE().reset();
272 auto rule = [](int, int,
int p) ->
int {
return 2 *
p; };
273 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
274 auto rule_2 = [
this](int, int, int) {
return 2 *
oRder; };
275 CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(rule_2);
277 auto add_base_ops = [&](
auto &pipeline) {
278 auto det_ptr = boost::make_shared<VectorDouble>();
279 auto jac_ptr = boost::make_shared<MatrixDouble>();
280 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
286 auto u_vals_ptr = boost::make_shared<VectorDouble>();
287 auto u_grad_ptr = boost::make_shared<MatrixDouble>();
289 add_base_ops(pipeline_mng->getOpDomainRhsPipeline());
290 pipeline_mng->getOpDomainRhsPipeline().push_back(
292 pipeline_mng->getOpDomainRhsPipeline().push_back(
295 enum NORMS {
L2 = 0, SEMI_NORM,
H1, LAST_NORM };
296 std::array<int, LAST_NORM> error_indices;
297 for (
int i = 0;
i != LAST_NORM; ++
i)
298 error_indices[
i] =
i;
310 if (
const size_t nb_dofs = data.getIndices().size()) {
312 const int nb_integration_pts =
o->getGaussPts().size2();
313 auto t_w =
o->getFTensor0IntegrationWeight();
315 auto t_grad = getFTensor1FromMat<2>(*u_grad_ptr);
316 auto t_coords =
o->getFTensor1CoordsAtGaussPts();
318 std::array<double, LAST_NORM> error;
319 std::fill(error.begin(), error.end(), 0);
321 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
323 const double alpha = t_w *
o->getMeasure();
325 t_val -
u_exact(t_coords(0), t_coords(1), t_coords(2));
327 auto t_exact_grad =
u_grad_exact(t_coords(0), t_coords(1), t_coords(2));
329 const double diff_grad =
330 (t_grad(
i) - t_exact_grad(
i)) * (t_grad(
i) - t_exact_grad(
i));
332 error[
L2] += alpha * pow(diff, 2);
333 error[SEMI_NORM] += alpha * diff_grad;
341 error[
H1] = error[
L2] + error[SEMI_NORM];
350 auto side_fe_ptr = boost::make_shared<FaceSideEle>(
mField);
351 add_base_ops(side_fe_ptr->getOpPtrVector());
352 side_fe_ptr->getOpPtrVector().push_back(
354 std::array<VectorDouble, 2> side_vals;
355 std::array<double, 2> area_map;
358 side_vals_op->doWorkRhsHook = [&](
DataOperator *op_ptr,
int side,
364 const auto nb_in_loop =
o->getFEMethod()->nInTheLoop;
365 area_map[nb_in_loop] =
o->getMeasure();
366 side_vals[nb_in_loop] = *u_vals_ptr;
369 side_vals[1].clear();
374 side_fe_ptr->getOpPtrVector().push_back(side_vals_op);
381 CHKERR o->loopSideFaces(
"dFE", *side_fe_ptr);
382 const auto in_the_loop = side_fe_ptr->nInTheLoop;
385 const std::array<std::string, 2> ele_type_name = {
"BOUNDARY",
"SKELETON"};
387 <<
"do_work_rhs_error in_the_loop " << ele_type_name[in_the_loop];
390 const double s =
o->getMeasure() / (area_map[0] + area_map[1]);
393 constexpr std::array<int, 2> sign_array{1, -1};
395 std::array<double, LAST_NORM> error;
396 std::fill(error.begin(), error.end(), 0);
398 const int nb_integration_pts =
o->getGaussPts().size2();
401 side_vals[1].resize(nb_integration_pts,
false);
402 auto t_coords =
o->getFTensor1CoordsAtGaussPts();
404 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
405 t_val_m =
u_exact(t_coords(0), t_coords(1), t_coords(2));
413 auto t_w =
o->getFTensor0IntegrationWeight();
415 for (
size_t gg = 0; gg != nb_integration_pts; ++gg) {
417 const double alpha =
o->getMeasure() * t_w;
418 const auto diff = t_val_p - t_val_m;
419 error[SEMI_NORM] += alpha *
p * diff * diff;
427 error[
H1] = error[SEMI_NORM];
435 skeleton_error_op->doWorkRhsHook = do_work_rhs_error;
437 boundary_error_op->doWorkRhsHook = do_work_rhs_error;
439 pipeline_mng->getOpDomainRhsPipeline().push_back(error_op);
440 pipeline_mng->getOpSkeletonLhsPipeline().push_back(skeleton_error_op);
441 pipeline_mng->getOpBoundaryRhsPipeline().push_back(boundary_error_op);
443 CHKERR pipeline_mng->loopFiniteElements();
445 CHKERR VecAssemblyBegin(l2_vec);
446 CHKERR VecAssemblyEnd(l2_vec);
450 CHKERR VecGetArrayRead(l2_vec, &array);
451 MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm L2 %6.4e",
452 std::sqrt(array[
L2]));
453 MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm Energetic %6.4e",
454 std::sqrt(array[SEMI_NORM]));
455 MOFEM_LOG_C(
"SELF", Sev::inform,
"Error Norm H1 %6.4e",
456 std::sqrt(array[
H1]));
459 constexpr double eps = 1e-12;
460 if (std::sqrt(array[
H1]) >
eps)
464 CHKERR VecRestoreArrayRead(l2_vec, &array);
482 pipeline_mng->getSkeletonRhsFE().reset();
483 pipeline_mng->getSkeletonLhsFE().reset();
484 pipeline_mng->getBoundaryRhsFE().reset();
485 pipeline_mng->getBoundaryLhsFE().reset();
487 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
489 auto u_ptr = boost::make_shared<VectorDouble>();
490 post_proc_fe->getOpPtrVector().push_back(
495 post_proc_fe->getOpPtrVector().push_back(
499 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
511 pipeline_mng->getDomainRhsFE() = post_proc_fe;
512 CHKERR pipeline_mng->loopFiniteElements();
513 CHKERR post_proc_fe->writeFile(
"out_result.h5m");
537int main(
int argc,
char *argv[]) {
546 DMType dm_name =
"DMMOFEM";
550 moab::Core mb_instance;
551 moab::Interface &moab = mb_instance;
#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)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
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.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data 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.
bool & getAddBoundaryFE()
Get the addSkeletonFE.
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
bool & getAddSkeletonFE()
Get the addSkeletonFE.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
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 tp collect data from elements on the side of Edge/Face.
Operator to evaluate Dirichlet boundary conditions using DG.
Operator the left hand side matrix.