Poisson2DiscontGalerkin Struct Reference

Collaboration diagram for Poisson2DiscontGalerkin:

Public Member Functions
	Poisson2DiscontGalerkin (MoFEM::Interface &m_field)
MoFEMErrorCode	runProgram ()
	[Output results]

Private Member Functions
MoFEMErrorCode	readMesh ()
	[Read mesh]
MoFEMErrorCode	setupProblem ()
	[Read mesh]
MoFEMErrorCode	boundaryCondition ()
	[Setup problem]
MoFEMErrorCode	assembleSystem ()
	[Boundary condition]
MoFEMErrorCode	setIntegrationRules ()
	[Assemble system]
MoFEMErrorCode	solveSystem ()
	[Set integration rules]
MoFEMErrorCode	checkResults ()
	[Solve system]
MoFEMErrorCode	outputResults ()
	[Output results]

Private Attributes
MoFEM::Interface &	mField
Simple *	simpleInterface
std::string	domainField
int	oRder

Detailed Description

Definition at line 74 of file poisson_2d_dis_galerkin.cpp.

Constructor & Destructor Documentation

Poisson2DiscontGalerkin()

Poisson2DiscontGalerkin::Poisson2DiscontGalerkin

(

MoFEM::Interface &

m_field

)

Definition at line 101 of file poisson_2d_dis_galerkin.cpp.

    : domainField("U"), mField(m_field), oRder(4) {}

Member Function Documentation

assembleSystem()

MoFEMErrorCode Poisson2DiscontGalerkin::assembleSystem ( )

private

[Boundary condition]

[Assemble system]

Definition at line 174 of file poisson_2d_dis_galerkin.cpp.

                                                       {
  MoFEMFunctionBegin;
 
  auto pipeline_mng = mField.getInterface<PipelineManager>();
 
  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>();
    pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
    pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
    pipeline.push_back(new OpSetInvJacL2ForFace(inv_jac_ptr));
  };
 
  add_base_ops(pipeline_mng->getOpDomainLhsPipeline());
  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpDomainGradGrad(
      domainField, domainField,
      [](const double, const double, const double) { return 1; }));
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpDomainSource(domainField, source));
 
  // Push operators to the Pipeline for Skeleton
  auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
  add_base_ops(side_fe_ptr->getOpPtrVector());
  side_fe_ptr->getOpPtrVector().push_back(
      new OpCalculateSideData(domainField, domainField));
 
  // Push operators to the Pipeline for Skeleton
  pipeline_mng->getOpSkeletonLhsPipeline().push_back(
      new OpL2LhsPenalty(side_fe_ptr));
 
  // Push operators to the Pipeline for Boundary
  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
      new OpL2LhsPenalty(side_fe_ptr));
  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
      new OpL2BoundaryRhs(side_fe_ptr, u_exact));
 
  MoFEMFunctionReturn(0);
}

boundaryCondition()

MoFEMErrorCode Poisson2DiscontGalerkin::boundaryCondition ( )

private

[Setup problem]

[Boundary condition]

Definition at line 166 of file poisson_2d_dis_galerkin.cpp.

                                                          {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}

checkResults()

MoFEMErrorCode Poisson2DiscontGalerkin::checkResults ( )

private

[Solve system]

Definition at line 264 of file poisson_2d_dis_galerkin.cpp.

                                                     {
  MoFEMFunctionBegin;
 
  auto pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto rule = [](int, int, int p) -> int { return 2 * p; };
  CHKERR pipeline_mng->setEvaluationIntegrationRule(rule);
  auto rule_2 = [this](int, int, int) { return 2 * oRder; };
  auto &evaluation_pipeline = pipeline_mng->getOpEvaluationPipeline();
 
  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>();
    pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
    pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
    pipeline.push_back(new OpSetInvJacL2ForFace(inv_jac_ptr));
  };
 
  auto u_vals_ptr = boost::make_shared<VectorDouble>();
  auto u_grad_ptr = boost::make_shared<MatrixDouble>();
 
  add_base_ops(evaluation_pipeline);
  evaluation_pipeline.push_back(
      new OpCalculateScalarFieldValues(domainField, u_vals_ptr));
  evaluation_pipeline.push_back(
      new OpCalculateScalarFieldGradient<2>(domainField, u_grad_ptr));
 
  enum NORMS { L2 = 0, SEMI_NORM, H1, LAST_NORM };
  std::array<int, LAST_NORM> error_indices;
  for (int i = 0; i != LAST_NORM; ++i)
    error_indices[i] = i;
  auto l2_vec = createVectorMPI(
      mField.get_comm(), (!mField.get_comm_rank()) ? LAST_NORM : 0, LAST_NORM);
 
  auto error_op = new DomainEleOp(domainField, DomainEleOp::OPROW);
  error_op->doWorkRhsHook = [&](DataOperator *op_ptr, int side, EntityType type,
                                EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
    auto o = static_cast<DomainEleOp *>(op_ptr);
 
    FTensor::Index<'i', 2> 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[L2] += alpha * pow(diff, 2);
        error[SEMI_NORM] += alpha * diff_grad;
 
        ++t_w;
        ++t_val;
        ++t_grad;
        ++t_coords;
      }
 
      error[H1] = error[L2] + error[SEMI_NORM];
 
      CHKERR VecSetValues(l2_vec, LAST_NORM, error_indices.data(), error.data(),
                          ADD_VALUES);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
  add_base_ops(side_fe_ptr->getOpPtrVector());
  side_fe_ptr->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues(domainField, u_vals_ptr));
  std::array<VectorDouble, 2> side_vals;
  std::array<double, 2> area_map;
  std::array<EntityHandle, 2> side_ents = {0, 0};
 
  auto side_vals_op = new DomainEleOp(domainField, DomainEleOp::OPROW);
  side_vals_op->doWorkRhsHook = [&](DataOperator *op_ptr, int side,
                                    EntityType type,
                                    EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
    auto o = static_cast<FaceSideOp *>(op_ptr);
 
    const auto nb_in_loop = o->getFEMethod()->nInTheLoop;
    area_map[nb_in_loop] = o->getMeasure();
    side_vals[nb_in_loop] = *u_vals_ptr;
    side_ents[nb_in_loop] = o->getFEEntityHandle();
    if (!nb_in_loop) {
      area_map[1] = 0;
      side_vals[1].clear();
      side_ents[1] = 0;
    }
 
    MoFEMFunctionReturn(0);
  };
  side_fe_ptr->getOpPtrVector().push_back(side_vals_op);
 
  auto do_work_rhs_error = [&](DataOperator *op_ptr, int side, EntityType type,
                               EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
    auto o = static_cast<BoundaryEleOp *>(op_ptr);
 
    CHKERR o->loopSideFaces("dFE", *side_fe_ptr);
    const auto in_the_loop = side_fe_ptr->nInTheLoop;
    if (in_the_loop) {
      if (auto side_ptr_fe = o->getSidePtrFE()) {
        if (side_ptr_fe->numeredEntFiniteElementPtr->getEnt() != side_ents[0]) {
          MoFEMFunctionReturnHot(0);
        }
      }
    }
 
#ifndef NDEBUG
    const std::array<std::string, 2> ele_type_name = {"BOUNDARY", "SKELETON"};
    MOFEM_LOG("SELF", Sev::noisy)
        << "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]);
    const double p = penalty * s;
 
    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];
    CHKERR VecSetValues(l2_vec, LAST_NORM, error_indices.data(), error.data(),
                        ADD_VALUES);
 
    MoFEMFunctionReturn(0);
  };
 
  auto skeleton_error_op = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
  skeleton_error_op->doWorkRhsHook = do_work_rhs_error;
  auto boundary_error_op = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
  boundary_error_op->doWorkRhsHook = do_work_rhs_error;
 
  evaluation_pipeline.push_back(error_op);
 
  auto add_side_error_op = [&](const std::string &fe_name, auto error_op) {
    auto op_loop_side = new OpLoopSide<BoundaryEle>(
        mField, fe_name, SPACE_DIM - 1, Sev::verbose);
    op_loop_side->getSideFEPtr()->getRuleHook = rule_2;
    op_loop_side->getOpPtrVector().push_back(error_op);
    evaluation_pipeline.push_back(op_loop_side);
  };
 
  add_side_error_op(simpleInterface->getSkeletonFEName(), skeleton_error_op);
  add_side_error_op(simpleInterface->getBoundaryFEName(), boundary_error_op);
 
  CHKERR pipeline_mng->loopFiniteElementsEvaluation();
 
  CHKERR VecAssemblyBegin(l2_vec);
  CHKERR VecAssemblyEnd(l2_vec);
 
  if (mField.get_comm_rank() == 0) {
    const double *array;
    CHKERR VecGetArrayRead(l2_vec, &array);
    MOFEM_LOG_C("SELF", Sev::inform, "Error Norm L2 %6.4e",
                std::sqrt(array[L2]));
    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",
                std::sqrt(array[H1]));
 
    if(is_test) {
      constexpr double eps = 1e-12;
      if (std::sqrt(array[H1]) > eps)
        SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Error is too big");
    }
                
    CHKERR VecRestoreArrayRead(l2_vec, &array);
    const MoFEM::Problem *problem_ptr;
    CHKERR DMMoFEMGetProblemPtr(simpleInterface->getDM(), &problem_ptr);
    MOFEM_LOG_C("SELF", Sev::inform, "Nb. DOFs %d",
                problem_ptr->getNbDofsRow());
  }
 
  MoFEMFunctionReturn(0);
}

outputResults()

MoFEMErrorCode Poisson2DiscontGalerkin::outputResults ( )

private

[Output results]

Definition at line 489 of file poisson_2d_dis_galerkin.cpp.

                                                      {
  MoFEMFunctionBegin;
 
  auto pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
 
  auto u_ptr = boost::make_shared<VectorDouble>();
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues(domainField, u_ptr));
 
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
  post_proc_fe->getOpPtrVector().push_back(
 
      new OpPPMap(
 
          post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
          {{"U", u_ptr}},
 
          {},
 
          {},
 
          {})
 
  );
 
  pipeline_mng->getDomainPostProcFE() = post_proc_fe;
  CHKERR pipeline_mng->loopFiniteElementsPostProc();
  CHKERR post_proc_fe->writeFile("out_result.h5m");
 
  MoFEMFunctionReturn(0);
}

readMesh()

MoFEMErrorCode Poisson2DiscontGalerkin::readMesh ( )

private

[Read mesh]

Definition at line 105 of file poisson_2d_dis_galerkin.cpp.

                                                 {
  MoFEMFunctionBegin;
 
  CHKERR mField.getInterface(simpleInterface);
  CHKERR simpleInterface->getOptions();
 
  // Only L2 field is set in this example. Two lines bellow forces simple
  // interface to creat lower dimension (edge) elements, despite that fact that
  // there is no field spanning on such elements. We need them for DG method.
  simpleInterface->getAddSkeletonFE() = true;
  simpleInterface->getAddBoundaryFE() = true;
 
  CHKERR simpleInterface->loadFile();
 
  MoFEMFunctionReturn(0);
}

runProgram()

MoFEMErrorCode Poisson2DiscontGalerkin::runProgram

(

)

[Output results]

[Run program]

Definition at line 527 of file poisson_2d_dis_galerkin.cpp.

                                                   {
  MoFEMFunctionBegin;
 
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR setIntegrationRules();
  CHKERR solveSystem();
  CHKERR outputResults();
  CHKERR checkResults();
 
  MoFEMFunctionReturn(0);
}

setIntegrationRules()

MoFEMErrorCode Poisson2DiscontGalerkin::setIntegrationRules ( )

private

[Assemble system]

[Set integration rules]

Definition at line 216 of file poisson_2d_dis_galerkin.cpp.

                                                            {
  MoFEMFunctionBegin;
 
  auto rule_lhs = [](int, int, int p) -> int { return 2 * p; };
  auto rule_rhs = [](int, int, int p) -> int { return 2 * p; };
  auto rule_2 = [this](int, int, int) { return 2 * oRder; };
 
  auto pipeline_mng = mField.getInterface<PipelineManager>();
  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);
 
  MoFEMFunctionReturn(0);
}

setupProblem()

MoFEMErrorCode Poisson2DiscontGalerkin::setupProblem ( )

private

[Read mesh]

[Setup problem]

Definition at line 124 of file poisson_2d_dis_galerkin.cpp.

                                                     {
  MoFEMFunctionBegin;
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &oRder, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-penalty", &penalty,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-phi", &phi, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-nitsche", &nitsche,
                               PETSC_NULLPTR);
  PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_test", &is_test, PETSC_NULLPTR);
 
  MOFEM_LOG("WORLD", Sev::inform) << "Set order: " << oRder;
  MOFEM_LOG("WORLD", Sev::inform) << "Set penalty: " << penalty;
  MOFEM_LOG("WORLD", Sev::inform) << "Set phi: " << phi;
  MOFEM_LOG("WORLD", Sev::inform) << "Set nitche: " << nitsche;
  MOFEM_LOG("WORLD", Sev::inform) << "Set test: " << (is_test == PETSC_TRUE);
 
  CHKERR simpleInterface->addDomainField(domainField, L2,
                                         AINSWORTH_LOBATTO_BASE, 1);
  CHKERR simpleInterface->setFieldOrder(domainField, oRder);
  CHKERR simpleInterface->setUp();
 
  // This is only for debigging and experimentation, to see boundary and edge
  // elements.
  [[maybe_unused]] auto save_shared = [&](auto meshset, std::string prefix) {
    MoFEMFunctionBegin;
    auto file_name =
        prefix + "_" +
        boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk";
    CHKERR mField.get_moab().write_file(file_name.c_str(), "VTK", "", &meshset,
                                        1);
    MoFEMFunctionReturn(0);
  };
 
  // CHKERR save_shared(simpleInterface->getBoundaryMeshSet(), "bdy");
  // CHKERR save_shared(simpleInterface->getSkeletonMeshSet(), "skel");
 
  MoFEMFunctionReturn(0);
}

solveSystem()

MoFEMErrorCode Poisson2DiscontGalerkin::solveSystem ( )

private

[Set integration rules]

[Solve system]

Definition at line 237 of file poisson_2d_dis_galerkin.cpp.

                                                    {
  MoFEMFunctionBegin;
 
  auto pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto ksp_solver = pipeline_mng->createKSP();
  CHKERR KSPSetFromOptions(ksp_solver);
  CHKERR KSPSetUp(ksp_solver);
 
  // Create RHS and solution vectors
  auto dm = simpleInterface->getDM();
  auto F = createDMVector(dm);
  auto D = vectorDuplicate(F);
 
  CHKERR KSPSolve(ksp_solver, F, D);
  CHKERR mField.getInterface<TestingInterface>()->checkSolver(
      ksp_solver, "poisson 2d discontinuous Galerkin");
 
  // Scatter result data on the mesh
  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

domainField

std::string Poisson2DiscontGalerkin::domainField

private

Definition at line 97 of file poisson_2d_dis_galerkin.cpp.

mField

MoFEM::Interface& Poisson2DiscontGalerkin::mField

private

Definition at line 93 of file poisson_2d_dis_galerkin.cpp.

oRder

int Poisson2DiscontGalerkin::oRder

private

Definition at line 98 of file poisson_2d_dis_galerkin.cpp.

simpleInterface

Simple* Poisson2DiscontGalerkin::simpleInterface

private

Definition at line 94 of file poisson_2d_dis_galerkin.cpp.

---

The documentation for this struct was generated from the following file:

• tutorials/scl-11_poisson_2d_dg/poisson_2d_dis_galerkin.cpp