SCL-11: Discontinuous Galerkin for Poisson problem

Table of Contents

• Problem formulation
• Source code main code

Note

Prerequisites of this tutorial include SCL-1: Poisson's equation (homogeneous BC)

Note

Intended learning outcome:

• general structure of a program developed using MoFEM
• idea of Simple Interface in MoFEM and how to use it
• integartion skeleton

Problem formulation

\[ -(v_{,j}, \sigma_{j}(\mathbf{u}))_V+(v, \overline{t})_{S^\sigma} = 0 \]

\[ -(v_{j}, \sigma_{j}(\mathbf{u}))_\Omega -\sum_{i} (\left\{[v n_j] - \theta\gamma\sigma_{j}(\mathbf{v})\right\} + \theta\gamma\sigma_{j}(\mathbf{v}), \hat{\sigma}_{j}(\mathbf{u}))_{F_i}+(v_i, \overline{t}_i)_{S^\sigma}= 0 \]

\[ -(v_{j}, \sigma_{j}(\mathbf{u}))_\Omega+(v_i, \overline{t})_{S^\sigma} -\sum_{i} (\left\{[v_i n_j] - \theta\gamma\sigma_{ij}(\mathbf{v})\right\}, \hat{\sigma}_{j}(\mathbf{u}))_{F_i}-\sum_{i} (\theta\gamma\sigma_{j}(\mathbf{v}),\sigma_{j}(\mathbf{u}))_{F_i} = 0 \]

\[ \hat{\sigma} = \gamma^{-1} \left\{ [u n_j] - \gamma \sigma_{j} (\mathbf{u}) \right\} \]

\[ -(v_{j}, \sigma_{j}(\mathbf{u}))_\Omega+(v, \overline{t})_{S^\sigma} -\sum_{i} (\left\{[v n_j] - \theta\gamma\sigma_{j}(\mathbf{v})\right\}, \gamma^{-1} \left\{ [u n_j] - \gamma \sigma_{j}(\mathbf{u})\right\})_{F_i} -\sum_{i} (\theta\gamma\sigma_{j}(\mathbf{v}), \gamma^{-1} \left\{\gamma \sigma_{j}(\mathbf{u})\right\})_{F_i} = 0 \]

Source code main code

/**
 * \file poisson_2d_dis_galerkin.cpp
 * \example poisson_2d_dis_galerkin.cpp
 *
 * Example of implementation for discontinuous Galerkin.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
constexpr int BASE_DIM = 1;
constexpr int FIELD_DIM = 1;
constexpr int SPACE_DIM = 2;
 
using EntData = EntitiesFieldData::EntData;
 
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
 
using BoundaryEle =
    PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
using FaceSideEle =
    PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
using FaceSideOp = FaceSideEle::UserDataOperator;
 
using PostProcEle =  PostProcBrokenMeshInMoab<DomainEle>;
 
static double penalty = 1e6;
static double phi =
    -1; // 1 - symmetric Nitsche, 0 - nonsymmetric, -1 antisymmetrica
static double nitsche = 1;
 
#include <PoissonDiscontinousGalerkin.hpp>
 
using OpDomainGradGrad = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::BiLinearForm<GAUSS>::OpGradGrad<BASE_DIM, FIELD_DIM, SPACE_DIM>;
using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::LinearForm<GAUSS>::OpSource<BASE_DIM, FIELD_DIM>;
 
PetscBool is_test = PETSC_FALSE;
 
auto u_exact = [](const double x, const double y, const double) {
  if (is_test)
    return x * x * y * y;
  else
    return cos(2 * x * M_PI) * cos(2 * y * M_PI);
};
 
auto u_grad_exact = [](const double x, const double y, const double) {
  if (is_test)
    return FTensor::Tensor1<double, 2>{2 * x * y * y, 2 * x * x * y};
  else
    return FTensor::Tensor1<double, 2>{
 
        -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) {
  if (is_test)
    return -(2 * x * x + 2 * y * y);
  else
    return 8 * M_PI * M_PI * cos(2 * x * M_PI) * cos(2 * y * M_PI);
};
 
using namespace MoFEM;
using namespace Poisson2DiscontGalerkinOperators;
 
static char help[] = "...\n\n";
 
struct Poisson2DiscontGalerkin {
public:
  Poisson2DiscontGalerkin(MoFEM::Interface &m_field);
 
  // Declaration of the main function to run analysis
  MoFEMErrorCode runProgram();
 
private:
  // Declaration of other main functions called in runProgram()
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode boundaryCondition();
  MoFEMErrorCode assembleSystem();
  MoFEMErrorCode setIntegrationRules();
  MoFEMErrorCode solveSystem();
  MoFEMErrorCode checkResults();
  MoFEMErrorCode outputResults();
 
  // MoFEM interfaces
  MoFEM::Interface &mField;
  Simple *simpleInterface;
 
  // Field name and approximation order
  std::string domainField;
  int oRder;
};
 
Poisson2DiscontGalerkin::Poisson2DiscontGalerkin(MoFEM::Interface &m_field)
    : domainField("U"), mField(m_field), oRder(4) {}
 
//! [Read mesh]
MoFEMErrorCode Poisson2DiscontGalerkin::readMesh() {
  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);
}
//! [Read mesh]
 
//! [Setup problem]
MoFEMErrorCode Poisson2DiscontGalerkin::setupProblem() {
  MoFEMFunctionBegin;
 
  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-penalty", &penalty,
                               PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-phi", &phi, PETSC_NULL);
  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-nitsche", &nitsche,
                               PETSC_NULL);
  PetscOptionsGetBool(PETSC_NULL, "", "-is_test", &is_test, PETSC_NULL);
 
  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.
  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);
}
//! [Setup problem]
 
//! [Boundary condition]
MoFEMErrorCode Poisson2DiscontGalerkin::boundaryCondition() {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}
 
//! [Boundary condition]
 
//! [Assemble system]
MoFEMErrorCode Poisson2DiscontGalerkin::assembleSystem() {
  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);
}
//! [Assemble system]
 
//! [Set integration rules]
MoFEMErrorCode Poisson2DiscontGalerkin::setIntegrationRules() {
  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);
}
//! [Set integration rules]
 
//! [Solve system]
MoFEMErrorCode Poisson2DiscontGalerkin::solveSystem() {
  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);
 
  // 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);
}
//! [Solve system]
 
MoFEMErrorCode Poisson2DiscontGalerkin::checkResults() {
  MoFEMFunctionBegin;
 
  auto pipeline_mng = mField.getInterface<PipelineManager>();
  pipeline_mng->getDomainRhsFE().reset();
  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);
  auto rule_2 = [this](int, int, int) { return 2 * oRder; };
  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>();
    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(pipeline_mng->getOpDomainRhsPipeline());
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpCalculateScalarFieldValues(domainField, u_vals_ptr));
  pipeline_mng->getOpDomainRhsPipeline().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;
 
  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;
    if (!nb_in_loop) {
      area_map[1] = 0;
      side_vals[1].clear();
    }
 
    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;
 
#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;
 
    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];
    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;
 
  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);
 
  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);
}
 
//! [Output results]
MoFEMErrorCode Poisson2DiscontGalerkin::outputResults() {
  MoFEMFunctionBegin;
 
  auto pipeline_mng = mField.getInterface<PipelineManager>();
  pipeline_mng->getDomainLhsFE().reset();
  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(
      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->getDomainRhsFE() = post_proc_fe;
  CHKERR pipeline_mng->loopFiniteElements();
  CHKERR post_proc_fe->writeFile("out_result.h5m");
 
  MoFEMFunctionReturn(0);
}
//! [Output results]
 
//! [Run program]
MoFEMErrorCode Poisson2DiscontGalerkin::runProgram() {
  MoFEMFunctionBegin;
 
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR setIntegrationRules();
  CHKERR solveSystem();
  CHKERR outputResults();
  CHKERR checkResults();
 
  MoFEMFunctionReturn(0);
}
//! [Run program]
 
//! [Main]
int main(int argc, char *argv[]) {
 
  // Initialisation of MoFEM/PETSc and MOAB data structures
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  // Error handling
  try {
    // Register MoFEM discrete manager in PETSc
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    // Create MOAB instance
    moab::Core mb_instance;              // mesh database
    moab::Interface &moab = mb_instance; // mesh database interface
 
    // Create MoFEM instance
    MoFEM::Core core(moab);           // finite element database
    MoFEM::Interface &m_field = core; // finite element interface
 
    // Run the main analysis
    Poisson2DiscontGalerkin poisson_problem(m_field);
    CHKERR poisson_problem.runProgram();
  }
  CATCH_ERRORS;
 
  // Finish work: cleaning memory, getting statistics, etc.
  MoFEM::Core::Finalize();
 
  return 0;
}
//! [Main]

Source code main code

/**
 * \file PoissonDiscontinousGalerkin.hpp
 * \example PoissonDiscontinousGalerkin.hpp
 *
 * Example of implementation for discontinuous Galerkin.
 */
 
 
 
// Define name if it has not been defined yet
#ifndef __POISSON2DISGALERKIN_HPP__
#define __POISSON2DISGALERKIN_HPP__
 
// Namespace that contains necessary UDOs, will be included in the main program
namespace Poisson2DiscontGalerkinOperators {
 
// Declare FTensor index for 2D problem
FTensor::Index<'i', SPACE_DIM> i;
 
enum ElementSide { LEFT_SIDE = 0, RIGHT_SIDE };
// data for skeleton computation
std::array<VectorInt, 2>
    indicesRowSideMap; ///< indices on rows for left hand-side
std::array<VectorInt, 2>
    indicesColSideMap; ///< indices on columns for left hand-side
std::array<MatrixDouble, 2> rowBaseSideMap; // base functions on rows
std::array<MatrixDouble, 2> colBaseSideMap; // base function  on columns
std::array<MatrixDouble, 2> rowDiffBaseSideMap; // derivative of base functions
std::array<MatrixDouble, 2> colDiffBaseSideMap; // derivative of base functions
std::array<double, 2> areaMap; // area/volume of elements on the side
std::array<int, 2> senseMap; // orientaton of local element edge/face in respect
                             // to global orientation of edge/face
 
/**
 * @brief Operator tp collect data from elements on the side of Edge/Face
 * 
 */
struct OpCalculateSideData : public FaceSideOp {
 
  OpCalculateSideData(std::string field_name, std::string col_field_name)
      : FaceSideOp(field_name, col_field_name, FaceSideOp::OPROWCOL) {
 
    std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
 
    for (auto t = moab::CN::TypeDimensionMap[SPACE_DIM].first;
         t <= moab::CN::TypeDimensionMap[SPACE_DIM].second; ++t)
      doEntities[t] = true;
  }
 
  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
                        EntityType col_type, EntData &row_data,
                        EntData &col_data) {
    MoFEMFunctionBeginHot;
 
    // Note: THat for L2 base data rows, and columns are the same, so operator
    // above can be simpler operator for the right hand side, and data can be
    // stored only for rows, since for columns data are the same. However for
    // complex multi-physics problems that not necessary would be a case. For
    // generality, we keep generic case.
 
    if ((CN::Dimension(row_type) == SPACE_DIM) &&
        (CN::Dimension(col_type) == SPACE_DIM)) {
      const auto nb_in_loop = getFEMethod()->nInTheLoop;
      indicesRowSideMap[nb_in_loop] = row_data.getIndices();
      indicesColSideMap[nb_in_loop] = col_data.getIndices();
      rowBaseSideMap[nb_in_loop] = row_data.getN();
      colBaseSideMap[nb_in_loop] = col_data.getN();
      rowDiffBaseSideMap[nb_in_loop] = row_data.getDiffN();
      colDiffBaseSideMap[nb_in_loop] = col_data.getDiffN();
      areaMap[nb_in_loop] = getMeasure();
      senseMap[nb_in_loop] = getSkeletonSense();
      if (!nb_in_loop) {
        indicesRowSideMap[1].clear();
        indicesColSideMap[1].clear();
        rowBaseSideMap[1].clear();
        colBaseSideMap[1].clear();
        rowDiffBaseSideMap[1].clear();
        colDiffBaseSideMap[1].clear();
        areaMap[1] = 0;
        senseMap[1] = 0;
      }
    } else {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Should not happen");
    }
 
    MoFEMFunctionReturnHot(0);
  }
};
 
template <typename T> inline auto get_ntensor(T &base_mat) {
  return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(
      &*base_mat.data().begin());
};
 
template <typename T> inline auto get_ntensor(T &base_mat, int gg, int bb) {
  double *ptr = &base_mat(gg, bb);
  return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(ptr);
};
 
template <typename T> inline auto get_diff_ntensor(T &base_mat) {
  double *ptr = &*base_mat.data().begin();
  return getFTensor1FromPtr<2>(ptr);
};
 
template <typename T>
inline auto get_diff_ntensor(T &base_mat, int gg, int bb) {
  double *ptr = &base_mat(gg, 2 * bb);
  return getFTensor1FromPtr<2>(ptr);
};
 
/**
 * @brief Operator the left hand side matrix 
 * 
 */
struct OpL2LhsPenalty : public BoundaryEleOp {
public:
 
  /**
   * @brief Construct a new OpL2LhsPenalty 
   * 
   * @param side_fe_ptr pointer to FE to evaluate side elements
   */
  OpL2LhsPenalty(boost::shared_ptr<FaceSideEle> side_fe_ptr)
      : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), sideFEPtr(side_fe_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
 
    // Collect data from side domian elements
    CHKERR loopSideFaces("dFE", *sideFEPtr);
    const auto in_the_loop =
        sideFEPtr->nInTheLoop; // return number of elements on the side
 
#ifndef NDEBUG
    const std::array<std::string, 2> ele_type_name = {"BOUNDARY", "SKELETON"};
    MOFEM_LOG("SELF", Sev::noisy)
        << "OpL2LhsPenalty inTheLoop " << ele_type_name[in_the_loop];
    MOFEM_LOG("SELF", Sev::noisy)
        << "OpL2LhsPenalty sense " << senseMap[0] << " " << senseMap[1];
#endif
 
    // calculate  penalty
    const double s = getMeasure() / (areaMap[0] + areaMap[1]);
    const double p = penalty * s;
 
    // get normal of the face or edge
    auto t_normal = getFTensor1Normal();
    t_normal.normalize();
 
    // get number of integration points on face
    const size_t nb_integration_pts = getGaussPts().size2();
 
    // beta paramter is zero, when penalty method is used, also, takes value 1,
    // when boundary edge/face is evaluated, and 2 when is skeleton edge/face.
    const double beta = static_cast<double>(nitsche) / (in_the_loop + 1);
 
    // iterate the sides rows
    for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
 
      // gent number of DOFs on the right side. 
      const auto nb_rows = indicesRowSideMap[s0].size();
 
      if (nb_rows) {
 
        // get orientation of the local element edge
        const auto sense_row = senseMap[s0];
 
        // iterate the side cols
        const auto nb_row_base_functions = rowBaseSideMap[s0].size2();
        for (auto s1 : {LEFT_SIDE, RIGHT_SIDE}) {
 
          // get orientation of the edge
          const auto sense_col = senseMap[s1];
 
          // number of dofs, for homogenous approximation this value is
          // constant.
          const auto nb_cols = indicesColSideMap[s1].size();
 
          // resize local element matrix
          locMat.resize(nb_rows, nb_cols, false);
          locMat.clear();
 
          // get base functions, and integration weights
          auto t_row_base = get_ntensor(rowBaseSideMap[s0]);
          auto t_diff_row_base = get_diff_ntensor(rowDiffBaseSideMap[s0]);
          auto t_w = getFTensor0IntegrationWeight();
 
          // iterate integration points on face/edge
          for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
 
            // t_w is integration weight, and measure is element area, or
            // volume, depending if problem is in 2d/3d.
            const double alpha = getMeasure() * t_w;
            auto t_mat = locMat.data().begin();
            
            // iterate rows
            size_t rr = 0;
            for (; rr != nb_rows; ++rr) {
 
              // calculate tetting function 
              FTensor::Tensor1<double, SPACE_DIM> t_vn_plus;
              t_vn_plus(i) = beta * (phi * t_diff_row_base(i) / p);
              FTensor::Tensor1<double, SPACE_DIM> t_vn;
              t_vn(i) = t_row_base * t_normal(i) * sense_row - t_vn_plus(i);
 
              // get base functions on columns
              auto t_col_base = get_ntensor(colBaseSideMap[s1], gg, 0);
              auto t_diff_col_base =
                  get_diff_ntensor(colDiffBaseSideMap[s1], gg, 0);
 
              // iterate columns
              for (size_t cc = 0; cc != nb_cols; ++cc) {
 
                // calculate variance of tested function 
                FTensor::Tensor1<double, SPACE_DIM> t_un;
                t_un(i) = -p * (t_col_base * t_normal(i) * sense_col -
                                beta * t_diff_col_base(i) / p);
 
                // assemble matrix
                *t_mat -= alpha * (t_vn(i) * t_un(i));
                *t_mat -= alpha * (t_vn_plus(i) * (beta * t_diff_col_base(i)));
 
                // move to next column base and element of matrix
                ++t_col_base;
                ++t_diff_col_base;
                ++t_mat;
              }
 
              // move to next row base
              ++t_row_base;
              ++t_diff_row_base;
            }
 
            // this is obsolete for this particular example, we keep it for
            // generality. in case of multi-physics problems different fields can
            // chare hierarchical base but use different approx. order, so is
            // possible to have more base functions than DOFs on element.
            for (; rr < nb_row_base_functions; ++rr) {
              ++t_row_base;
              ++t_diff_row_base;
            }
 
            ++t_w;
          }
 
          // assemble system
          CHKERR ::MatSetValues(getKSPB(), indicesRowSideMap[s0].size(),
                                &*indicesRowSideMap[s0].begin(),
                                indicesColSideMap[s1].size(),
                                &*indicesColSideMap[s1].begin(),
                                &*locMat.data().begin(), ADD_VALUES);
 
          // stop of boundary element
          if (!in_the_loop)
            MoFEMFunctionReturnHot(0);
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<FaceSideEle>
      sideFEPtr; ///< pointer to element to get data on edge/face sides
  MatrixDouble locMat; ///< local operator matrix
};
 
/**
 * @brief Operator to evaluate Dirichlet boundary conditions using DG
 *
 */
struct OpL2BoundaryRhs : public BoundaryEleOp {
public:
  OpL2BoundaryRhs(boost::shared_ptr<FaceSideEle> side_fe_ptr, ScalarFun fun)
      : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE), sideFEPtr(side_fe_ptr),
        sourceFun(fun) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
 
    // get normal of the face or edge
    CHKERR loopSideFaces("dFE", *sideFEPtr);
    const auto in_the_loop =
        sideFEPtr->nInTheLoop; // return number of elements on the side
 
    // calculate  penalty
    const double s = getMeasure() / (areaMap[0]);
    const double p = penalty * s;
 
    // get normal of the face or edge
    auto t_normal = getFTensor1Normal();
    t_normal.normalize();
 
    auto t_w = getFTensor0IntegrationWeight();
 
    const size_t nb_rows = indicesRowSideMap[0].size();
 
    if (!nb_rows)
      MoFEMFunctionReturnHot(0);
 
    // resize local operator vector
    rhsF.resize(nb_rows, false);
    rhsF.clear();
 
    const size_t nb_integration_pts = getGaussPts().size2();
    const size_t nb_row_base_functions = rowBaseSideMap[0].size2();
 
    // shape funcs
    auto t_row_base = get_ntensor(rowBaseSideMap[0]);
    auto t_diff_row_base = get_diff_ntensor(rowDiffBaseSideMap[0]);
    auto t_coords = getFTensor1CoordsAtGaussPts();
 
    const auto sense_row = senseMap[0];
    const double beta = static_cast<double>(nitsche) / (in_the_loop + 1);
 
    for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
      const double alpha = getMeasure() * t_w;
 
      const double source_val =
          -p * sourceFun(t_coords(0), t_coords(1), t_coords(2));
 
      auto t_f = rhsF.data().begin();
 
      size_t rr = 0;
      for (; rr != nb_rows; ++rr) {
 
        FTensor::Tensor1<double, SPACE_DIM> t_vn_plus;
        t_vn_plus(i) = beta * (phi * t_diff_row_base(i) / p);
        FTensor::Tensor1<double, SPACE_DIM> t_vn;
        t_vn(i) = t_row_base * t_normal(i) * sense_row - t_vn_plus(i);
 
        // assemble operator local vactor
        *t_f -= alpha * t_vn(i) * (source_val * t_normal(i));
 
        ++t_row_base;
        ++t_diff_row_base;
        ++t_f;
      }
 
      for (; rr < nb_row_base_functions; ++rr) {
        ++t_row_base;
        ++t_diff_row_base;
      }
 
      ++t_w;
      ++t_coords;
    }
 
    // assemble local operator vector to global vector
    CHKERR ::VecSetValues(getKSPf(), indicesRowSideMap[0].size(),
                          &*indicesRowSideMap[0].begin(), &*rhsF.data().begin(),
                          ADD_VALUES);
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<FaceSideEle>
      sideFEPtr;       ///< pointer to element to get data on edge/face sides
  ScalarFun sourceFun; ///< pointer to function to evaluate value of function on boundary
  VectorDouble rhsF;   ///< vector to store local operator right hand side
};
 
}; // namespace Poisson2DiscontGalerkinOperators
 
#endif //__POISSON2DISGALERKIN_HPP__