ADV-6: Map solution between meshes using DG projection

Note

Prerequisites of this tutorial include FUN-0: Hello world

Note

Intended learning outcome:

• Field evaluator
• DG projection
• User defined integration rules

/**
 * @file between_meshes_dg_projection.cpp
 * @example mofem/tutorials/adv-6/between_meshes_dg_projection.cpp
 * 
 * @brief Testing Discontinuous Galerkin (DG) projection operators
 *
 *
 */
 
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
static char help[] = "DG Projection Test - validates discontinuous Galerkin "
                     "projection accuracy\n\n";
 
constexpr char FIELD_NAME_U[] = "U";
constexpr char FIELD_NAME_S[] = "S";
constexpr int BASE_DIM = 1;
constexpr int FIELD_DIM = 1;
constexpr int SPACE_DIM = 2;
constexpr int order = 2;
 
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
using EntData = EntitiesFieldData::EntData;
 
using PostProcFaceEle = PostProcBrokenMeshInMoab<DomainEle>;
 
auto fun = [](const double x, const double y, const double z) {
  return x + y + x * x + y * y;
};
 
using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::BiLinearForm<GAUSS>::OpMass<BASE_DIM, FIELD_DIM>;
 
using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::LinearForm<GAUSS>::OpSource<BASE_DIM, FIELD_DIM>;
 
struct Example {
 
  Example(MoFEM::Interface &m_field) : mField(m_field) {}
 
  MoFEMErrorCode runProblem();
 
private:
  MoFEM::Interface &mField;
  Simple *simpleInterface;
 
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode assembleSystem();
  MoFEMErrorCode solveSystem();
 
  MoFEMErrorCode projectResults(BitRefLevel parent_bit, BitRefLevel child_bit,
                                BitRefLevel refine_bit);
  MoFEMErrorCode outputResults(std::string FIELD_NAME_U);
  MoFEMErrorCode edgeFlips(BitRefLevel parent_bit, BitRefLevel child_bit);
  MoFEMErrorCode refineSkin(BitRefLevel parent_bit, BitRefLevel refine_bit);
  MoFEMErrorCode reSetupProblem(BitRefLevel child_bit);
 
  struct CommonData {
    boost::shared_ptr<MatrixDouble> invJacPtr;
    boost::shared_ptr<VectorDouble> approxVals;
    boost::shared_ptr<MatrixDouble> approxGradVals;
    boost::shared_ptr<MatrixDouble> approxHessianVals;
    SmartPetscObj<Vec> L2Vec;
  };
 
  struct OpError;
};
 
auto save_range = [](moab::Interface &moab, const std::string name,
                     const Range r, std::vector<Tag> tags = {}) {
  MoFEMFunctionBegin;
  auto out_meshset = get_temp_meshset_ptr(moab);
  CHKERR moab.add_entities(*out_meshset, r);
  if (r.size()) {
    CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1,
                           tags.data(), tags.size());
  } else {
    MOFEM_LOG("SELF", Sev::warning) << "Empty range for " << name;
  }
  MoFEMFunctionReturn(0);
};
 
struct Example::OpError : public DomainEleOp {
  boost::shared_ptr<CommonData> commonDataPtr;
 
  OpError(boost::shared_ptr<MatrixDouble> data_ptr,
          BitRefLevel bits = BitRefLevel(), BitRefLevel mask = BitRefLevel())
      : DomainEleOp(NOSPACE, OPSPACE), dataPtr(data_ptr), bitsEle(bits),
        maskEle(mask) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
 
    auto fe_ptr = getNumeredEntFiniteElementPtr();
    auto fe_bit = fe_ptr->getBitRefLevel();
    if ((fe_bit & bitsEle).any() && ((fe_bit & maskEle) == fe_bit)) {
      const int nb_integration_pts = getGaussPts().size2();
 
      auto t_val = getFTensor1FromMat<1>(*(dataPtr));
      auto t_coords = getFTensor1CoordsAtGaussPts();
 
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
        double projected_value = t_val(0);
        double analytical_value = fun(t_coords(0), t_coords(1), t_coords(2));
        double error = projected_value - analytical_value;
 
        constexpr double eps = 1e-8;
        if (std::abs(error) > eps) {
          MOFEM_LOG("SELF", Sev::error)
              << "Projection error too large: " << error << " at point ("
              << t_coords(0) << ", " << t_coords(1) << ")"
              << " projected=" << projected_value
              << " analytical=" << analytical_value;
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "DG projection failed accuracy test");
        }
 
        ++t_val;
        ++t_coords;
      }
 
      MOFEM_LOG("SELF", Sev::noisy)
          << "DG projection accuracy validation passed";
    }
 
    MoFEMFunctionReturn(0);
  }
 
private:
  boost::shared_ptr<MatrixDouble> dataPtr;
  BitRefLevel bitsEle;
  BitRefLevel maskEle;
};
 
//! [Run programme]
MoFEMErrorCode Example::runProblem() {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR outputResults("out_initial.h5m");
 
  auto parent_bit = BitRefLevel().set(0);
  auto child_bit = BitRefLevel().set(1);
  auto refine_bit = BitRefLevel().set(2);
 
  CHKERR edgeFlips(parent_bit, child_bit);
  CHKERR refineSkin(child_bit, refine_bit);
  CHKERR reSetupProblem(BitRefLevel().set());
  CHKERR projectResults(parent_bit, child_bit, refine_bit);
 
  CHKERR reSetupProblem(refine_bit);
  CHKERR outputResults("out_projected.h5m");
 
  MoFEMFunctionReturn(0);
}
//! [Run programme]
 
//! [Read mesh]
MoFEMErrorCode Example::readMesh() {
  MoFEMFunctionBegin;
 
  CHKERR mField.getInterface(simpleInterface);
 
  CHKERR simpleInterface->getOptions();
 
  char mesh_File_Name[255];
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-file_name",
                               mesh_File_Name, 255, PETSC_NULLPTR);
  CHKERR simpleInterface->loadFile("", mesh_File_Name);
 
  MoFEMFunctionReturn(0);
}
//! [Read mesh]
 
//! [Set up problem]
MoFEMErrorCode Example::setupProblem() {
  MoFEMFunctionBegin;
 
  CHKERR simpleInterface->addDomainField(FIELD_NAME_U, L2,
                                         AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
  CHKERR simpleInterface->addDomainField(FIELD_NAME_S, H1,
                                         AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
 
  CHKERR simpleInterface->setFieldOrder(FIELD_NAME_U, order);
  CHKERR simpleInterface->setFieldOrder(FIELD_NAME_S, order);
 
  CHKERR simpleInterface->setUp(PETSC_FALSE);
 
  MoFEMFunctionReturn(0);
}
//! [Set up problem]
 
//! [Push operators to pipeline]
MoFEMErrorCode Example::assembleSystem() {
  MoFEMFunctionBegin;
 
  auto rule = [](int, int, int p) -> int { return 2 * p; };
 
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
 
  CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
 
  auto beta = [](const double, const double, const double) { return 1; };
 
  pipeline_mng->getOpDomainLhsPipeline().push_back(
      new OpDomainMass(FIELD_NAME_U, FIELD_NAME_U, beta));
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpDomainSource(FIELD_NAME_U, fun));
 
  pipeline_mng->getOpDomainLhsPipeline().push_back(
      new OpDomainMass(FIELD_NAME_S, FIELD_NAME_S, beta));
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpDomainSource(FIELD_NAME_S, fun));
 
 
  MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
 
//! [Solve]
MoFEMErrorCode Example::solveSystem() {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
 
  MOFEM_LOG("WORLD", Sev::inform) << "Solving DG projection system";
 
  auto solver = pipeline_mng->createKSP();
  CHKERR KSPSetFromOptions(solver);
  CHKERR KSPSetUp(solver);
 
  auto dm = simpleInterface->getDM();
  auto D = createDMVector(dm);
  auto F = vectorDuplicate(D);
 
  CHKERR KSPSolve(solver, F, D);
 
  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]
 
//! [Project results]
MoFEMErrorCode Example::projectResults(BitRefLevel parent_bit,
                                       BitRefLevel child_bit, 
                                       BitRefLevel refine_bit) {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pipeline_mng = mField.getInterface<PipelineManager>();
 
  pipeline_mng->getDomainLhsFE().reset();
  pipeline_mng->getDomainRhsFE().reset();
  pipeline_mng->getOpDomainRhsPipeline().clear();
 
  auto rule = [](int, int, int p) -> int { return 2 * p; };
  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
 
  // OpLoopThis, is child operator, and is use to execute 
  // fe_child_ptr, only on bit ref level and mask
  // for child elements
  auto get_child_op = [&](auto &pip) {
    auto op_this_child =
        new OpLoopThis<DomainEle>(mField, simple->getDomainFEName(), refine_bit,
                                  child_bit | refine_bit, Sev::noisy);
    auto fe_child_ptr = op_this_child->getThisFEPtr();
    fe_child_ptr->getRuleHook = [] (int, int, int p) { return -1; };
    Range child_edges;
    CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        refine_bit, child_bit | refine_bit, MBEDGE, child_edges);
    // set integration rule, such that integration points are not on flipped edge
    CHKERR setDGSetIntegrationPoints<SPACE_DIM>(
        fe_child_ptr->setRuleHook, [](int, int, int p) { return 2 * p; },
        boost::make_shared<Range>(child_edges));
    pip.push_back(op_this_child);
    return fe_child_ptr;
  };
 
  // Use field evaluator to calculate field values on parent bitref level, 
  // i.e. elements which were flipped.
  auto get_field_eval_op = [&](auto fe_child_ptr) {
    auto field_eval_ptr = mField.getInterface<FieldEvaluatorInterface>();
 
    // Get pointer of FieldEvaluator data. Note finite element and method
    // set integration points is destroyed when this pointer is releases
    auto field_eval_data = field_eval_ptr->getData<DomainEle>();
    // Build tree for particular element
    CHKERR field_eval_ptr->buildTree<SPACE_DIM>(
        field_eval_data, simpleInterface->getDomainFEName(), parent_bit,
        parent_bit | child_bit);
 
    // You can add more fields here
    auto data_U_ptr = boost::make_shared<MatrixDouble>();
    auto eval_data_U_ptr = boost::make_shared<MatrixDouble>();
    auto data_S_ptr = boost::make_shared<MatrixDouble>();
    auto eval_data_S_ptr = boost::make_shared<MatrixDouble>();
 
 
    if (auto fe_eval_ptr = field_eval_data->feMethodPtr.lock()) {
      fe_eval_ptr->getRuleHook = [] (int, int, int p) { return -1; };
      fe_eval_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<FIELD_DIM>(FIELD_NAME_U,
                                                      eval_data_U_ptr));
      fe_eval_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<FIELD_DIM>(FIELD_NAME_S,
                                                      eval_data_S_ptr));
 
      auto op_test = new DomainEleOp(NOSPACE, DomainEleOp::OPSPACE);
      op_test->doWorkRhsHook =
          [](DataOperator *base_op_ptr, int side, EntityType type,
             EntitiesFieldData::EntData &data) -> MoFEMErrorCode {
        MoFEMFunctionBegin;
 
        auto op_ptr = static_cast<DomainEleOp *>(base_op_ptr);
        MOFEM_LOG_CHANNEL("SELF");
        MOFEM_LOG("SELF", Sev::noisy)
            << "Field evaluator method pointer is valid";
        MOFEM_LOG("SELF", Sev::noisy)
            << op_ptr->getGaussPts();
        MOFEM_LOG("SELF", Sev::noisy)
            << "Loop size " << op_ptr->getPtrFE()->getLoopSize();
        MOFEM_LOG("SELF", Sev::noisy)
            << "Coords at gauss pts: " << op_ptr->getCoordsAtGaussPts();
 
        MoFEMFunctionReturn(0);
      };
 
      fe_eval_ptr->getOpPtrVector().push_back(op_test);
 
    } else {
      CHK_THROW_MESSAGE(MOFEM_ATOM_TEST_INVALID,
              "Field evaluator method pointer is expired");
    }
 
    auto op_ptr = field_eval_ptr->getDataOperator<SPACE_DIM>(
        {{eval_data_U_ptr, data_U_ptr}, {eval_data_S_ptr, data_S_ptr}},
        simpleInterface->getDomainFEName(), field_eval_data, 0,
        mField.get_comm_size(), parent_bit, parent_bit | child_bit, MF_EXIST,
        QUIET);
 
    fe_child_ptr->getOpPtrVector().push_back(op_ptr);
    return std::make_pair(
 
        std::vector<std::pair<std::string, boost::shared_ptr<MatrixDouble>>>{
            {FIELD_NAME_U, data_U_ptr}},
 
        std::vector<std::pair<std::string, boost::shared_ptr<MatrixDouble>>>{
            {FIELD_NAME_S, data_S_ptr}}
 
    );
 
  };
 
  // calculate coefficients on child (flipped) elements
  auto dg_projection_base = [&](auto fe_child_ptr, auto vec_data_ptr, auto mat,
                                auto vec) {
    MoFEMFunctionBegin;
    constexpr int projection_order = order;
    auto entity_data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
    auto mass_ptr = boost::make_shared<MatrixDouble>();
    auto coeffs_ptr = boost::make_shared<MatrixDouble>();
 
    // project L2 (directly from coefficients)
    for (auto &p : vec_data_ptr.first) {
      auto field_name = p.first;
      auto data_ptr = p.second;
 
      fe_child_ptr->getOpPtrVector().push_back(new OpDGProjectionMassMatrix(
          projection_order, mass_ptr, entity_data_l2, AINSWORTH_LEGENDRE_BASE,
          L2));
      fe_child_ptr->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
          data_ptr, coeffs_ptr, mass_ptr, entity_data_l2,
          AINSWORTH_LEGENDRE_BASE, L2));
 
      // next two lines are only for testing if projection is correct, they are not
      // essential
      fe_child_ptr->getOpPtrVector().push_back(new OpDGProjectionEvaluation(
          data_ptr, coeffs_ptr, entity_data_l2, AINSWORTH_LEGENDRE_BASE, L2));
      fe_child_ptr->getOpPtrVector().push_back(new OpError(data_ptr));
 
      // set coefficients to flipped element
      auto op_set_coeffs = new DomainEleOp(field_name, DomainEleOp::OPROW);
      op_set_coeffs->doWorkRhsHook =
          [coeffs_ptr](DataOperator *base_op_ptr, int side, EntityType type,
                       EntitiesFieldData::EntData &data) -> MoFEMErrorCode {
        MoFEMFunctionBegin;
        auto field_ents = data.getFieldEntities();
        auto nb_dofs = data.getIndices().size();
        if (!field_ents.size())
          MoFEMFunctionReturnHot(0);
        if (auto e_ptr = field_ents[0]) {
          auto field_ent_data = e_ptr->getEntFieldData();
          std::copy(coeffs_ptr->data().data(),
                    coeffs_ptr->data().data() + nb_dofs,
                    field_ent_data.begin());
        }
        MoFEMFunctionReturn(0);
      };
      fe_child_ptr->getOpPtrVector().push_back(op_set_coeffs);
    }
 
    // project H1 (via coefficients)
    for (auto &p : vec_data_ptr.second) {
      auto field_name = p.first;
      auto data_ptr = p.second;
 
      fe_child_ptr->getOpPtrVector().push_back(new OpDGProjectionMassMatrix(
          projection_order, mass_ptr, entity_data_l2, AINSWORTH_LEGENDRE_BASE,
          L2));
      fe_child_ptr->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
          data_ptr, coeffs_ptr, mass_ptr, entity_data_l2,
          AINSWORTH_LEGENDRE_BASE, L2));
 
      // next two lines are only for testing if projection is correct, they are not
      // essential
      fe_child_ptr->getOpPtrVector().push_back(new OpDGProjectionEvaluation(
          data_ptr, coeffs_ptr, entity_data_l2, AINSWORTH_LEGENDRE_BASE, L2));
      fe_child_ptr->getOpPtrVector().push_back(new OpError(data_ptr));
 
      // assemble to global matrix, since this is H1 (you will do the shame for Hcurl of Hdiv)
      auto beta = [](const double, const double, const double) { return 1; };
      fe_child_ptr->getOpPtrVector().push_back(
          new OpDomainMass(FIELD_NAME_S, FIELD_NAME_S, beta));
      using OpVec = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
          GAUSS>::OpBaseTimesVector<1, FIELD_DIM, 1>;
      fe_child_ptr->getOpPtrVector().push_back(
          new OpVec(FIELD_NAME_S, data_ptr, beta));
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto dm = simple->getDM();
  auto sub_dm = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMCreateSubDM(sub_dm, dm, "SUB");
  CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
  CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
 
  // get only refinement bit DOFs
  auto ref_entities_ptr = boost::make_shared<Range>();
  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByRefLevel(
      refine_bit, child_bit | refine_bit, *ref_entities_ptr);
  Range verts;
  CHKERR mField.get_moab().get_connectivity(*ref_entities_ptr, verts, true);
  ref_entities_ptr->merge(verts);
 
  CHKERR DMMoFEMAddSubFieldRow(sub_dm, FIELD_NAME_S, ref_entities_ptr);
  CHKERR DMMoFEMAddSubFieldCol(sub_dm, FIELD_NAME_S, ref_entities_ptr);
  CHKERR DMSetUp(sub_dm);
 
  auto mat = createDMMatrix(sub_dm);
  auto vec = createDMVector(sub_dm);
 
  // create child operator, and fe_child_ptr element in it
  auto fe_child_ptr = get_child_op(pipeline_mng->getOpDomainRhsPipeline());
  // run dg projection, note that get_field_eval_op, 
  // pass data_ptr values used to project and calculate coefficients
  CHKERR dg_projection_base(fe_child_ptr, get_field_eval_op(fe_child_ptr), mat,
                            vec);
 
  // That is to test, if projection works, and coefficients are set in correctly
  // Note: FIELD_S is not tested, it is in H1, so we have to solve KSP problem first
  auto test_U_data_ptr = boost::make_shared<MatrixDouble>();
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpCalculateVectorFieldValues<FIELD_DIM>(FIELD_NAME_U,
                                                  test_U_data_ptr));
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpError(test_U_data_ptr, refine_bit, BitRefLevel().set()));
 
  auto fe_rhs = pipeline_mng->getCastDomainRhsFE<DomainEle>();
  fe_rhs->ksp_A = mat;
  fe_rhs->ksp_B = mat;
  fe_rhs->ksp_f = vec;
  fe_rhs->data_ctx =
      PetscData::CtxSetF | PetscData::CtxSetA | PetscData::CtxSetB;
  CHKERR pipeline_mng->loopFiniteElements(sub_dm);
  
  CHKERR VecAssemblyBegin(vec);
  CHKERR VecAssemblyEnd(vec);
  CHKERR VecGhostUpdateBegin(vec, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateEnd(vec, ADD_VALUES, SCATTER_REVERSE); 
  CHKERR MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
 
  auto ksp = createKSP(mField.get_comm());
  CHKERR KSPSetOperators(ksp, mat, mat);
  CHKERR KSPSetFromOptions(ksp);
 
  auto sol = createDMVector(sub_dm);
  CHKERR KSPSolve(ksp, vec, sol);
  CHKERR VecGhostUpdateBegin(sol, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(sol, INSERT_VALUES, SCATTER_FORWARD); 
  CHKERR DMoFEMMeshToLocalVector(sub_dm, sol, INSERT_VALUES, SCATTER_REVERSE);
 
  pipeline_mng->getOpDomainRhsPipeline().clear();
  auto test_S_data_ptr = boost::make_shared<MatrixDouble>();
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpCalculateVectorFieldValues<FIELD_DIM>(FIELD_NAME_S,
                                                  test_S_data_ptr));
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpError(test_S_data_ptr, refine_bit, BitRefLevel().set()));
 
  MoFEMFunctionReturn(0);
}
//! [Project results]
 
//! [Output results]
MoFEMErrorCode Example::outputResults(std::string file_name) {
  MoFEMFunctionBegin;
 
  auto pipeline_mng = mField.getInterface<PipelineManager>();
  pipeline_mng->getDomainLhsFE().reset();
 
  auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      post_proc_fe->getOpPtrVector(), {H1});
 
  auto u_ptr = boost::make_shared<VectorDouble>();
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues(FIELD_NAME_U, u_ptr));
  auto s_ptr = boost::make_shared<VectorDouble>();    
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues(FIELD_NAME_S, s_ptr));
 
  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldGradient<SPACE_DIM>(FIELD_NAME_U, grad_u_ptr));
  auto grad_s_ptr = boost::make_shared<MatrixDouble>();
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldGradient<SPACE_DIM>(FIELD_NAME_S, grad_s_ptr));
 
 
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
  post_proc_fe->getOpPtrVector().push_back(
 
      new OpPPMap(
          post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
          OpPPMap::DataMapVec{{FIELD_NAME_U, u_ptr}, {FIELD_NAME_S, s_ptr}},
 
          OpPPMap::DataMapMat{
            
            {"GRAD_" + std::string(FIELD_NAME_U), grad_u_ptr},
            {"GRAD_" + std::string(FIELD_NAME_S), grad_s_ptr}
          
          },
 
          OpPPMap::DataMapMat{},
 
          OpPPMap::DataMapMat{}
 
          )
 
  );
 
  pipeline_mng->getDomainRhsFE() = post_proc_fe;
  CHKERR pipeline_mng->loopFiniteElements();
  CHKERR post_proc_fe->writeFile(file_name);
 
  MoFEMFunctionReturn(0);
}
//! [Output results]
 
//! [Edge flips]
MoFEMErrorCode Example::edgeFlips(BitRefLevel parent_bit,
                                  BitRefLevel child_bit) {
  MoFEMFunctionBegin;
 
  moab::Interface &moab = mField.get_moab();
 
  auto make_edge_flip = [&](auto edge, auto adj_faces, Range &new_tris) {
    MoFEMFunctionBegin;
 
    auto get_conn = [&](EntityHandle e, EntityHandle *conn_cpy) {
      MoFEMFunctionBegin;
      const EntityHandle *conn;
      int num_nodes;
      CHKERR moab.get_connectivity(e, conn, num_nodes, true);
      std::copy(conn, conn + num_nodes, conn_cpy);
      MoFEMFunctionReturn(0);
    };
 
    auto get_tri_normals = [&](auto &conn) {
      std::array<double, 18> coords;
      CHKERR moab.get_coords(conn.data(), 6, coords.data());
      std::array<FTensor::Tensor1<double, 3>, 2> tri_normals;
      for (int t = 0; t != 2; ++t) {
        CHKERR Tools::getTriNormal(&coords[9 * t], &tri_normals[t](0));
      }
      return tri_normals;
    };
 
    auto test_flip = [&](auto &&t_normals) {
      FTENSOR_INDEX(3, i);
      if (t_normals[0](i) * t_normals[1](i) <
          std::numeric_limits<float>::epsilon())
        return false;
      return true;
    };
 
    std::array<EntityHandle, 6> adj_conn;
    CHKERR get_conn(adj_faces[0], &adj_conn[0]);
    CHKERR get_conn(adj_faces[1], &adj_conn[3]);
    std::array<EntityHandle, 2> edge_conn;
    CHKERR get_conn(edge, edge_conn.data());
    std::array<EntityHandle, 2> new_edge_conn;
 
    int j = 1;
    for (int i = 0; i != 6; ++i) {
      if (adj_conn[i] != edge_conn[0] && adj_conn[i] != edge_conn[1]) {
        new_edge_conn[j] = adj_conn[i];
        --j;
      }
    }
 
    auto &new_conn = adj_conn; //< just alias this
    for (int t = 0; t != 2; ++t) {
      for (int i = 0; i != 3; ++i) {
        if (
 
            (adj_conn[3 * t + i % 3] == edge_conn[0] &&
             adj_conn[3 * t + (i + 1) % 3] == edge_conn[1])
 
            ||
 
            (adj_conn[3 * t + i % 3] == edge_conn[1] &&
             adj_conn[3 * t + (i + 1) % 3] == edge_conn[0])
 
        ) {
          new_conn[3 * t + (i + 1) % 3] = new_edge_conn[t];
          break;
        }
      }
    }
 
    if (test_flip(get_tri_normals(new_conn))) {
      for (int t = 0; t != 2; ++t) {
        Range rtri;
        CHKERR moab.get_adjacencies(&new_conn[3 * t], SPACE_DIM + 1, SPACE_DIM,
                                    false, rtri);
        if (!rtri.size()) {
          EntityHandle tri;
          CHKERR moab.create_element(MBTRI, &new_conn[3 * t], SPACE_DIM + 1,
                                     tri);
          new_tris.insert(tri);
        } else {
#ifndef NDEBUG
          if (rtri.size() != 1) {
            MOFEM_LOG("SELF", Sev::error)
                << "Multiple tries created during edge flip for edge " << edge
                << " adjacent faces " << std::endl
                << rtri;
            SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                    "Multiple tries created during edge flip");
          }
#endif // NDEBUG
          new_tris.merge(rtri);
        }
      }
 
      Range new_edges;
      CHKERR moab.get_adjacencies(new_tris, SPACE_DIM - 1, true, new_edges,
                                  moab::Interface::UNION);
    } else {
 
      MOFEM_LOG_CHANNEL("SELF");
      MOFEM_LOG("SELF", Sev::warning)
          << "Edge flip rejected for edge " << edge << " adjacent faces "
          << adj_faces;      
    }
 
    MoFEMFunctionReturn(0);
  };
 
  Range tris;
  CHKERR moab.get_entities_by_dimension(0, SPACE_DIM, tris);
  CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
      parent_bit, BitRefLevel().set(), tris);
  Skinner skin(&moab);
  Range skin_edges;
  CHKERR skin.find_skin(0, tris, false, skin_edges);
 
  Range edges;
  CHKERR moab.get_entities_by_dimension(0, SPACE_DIM - 1, edges);
  edges = subtract(edges, skin_edges);
  CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
      parent_bit, BitRefLevel().set(), edges);
 
  Range new_tris, flipped_tris, forbidden_tris;
  int flip_count = 0;
  for (auto edge : edges) {
    Range adjacent_tris;
    CHKERR moab.get_adjacencies(&edge, 1, SPACE_DIM, true, adjacent_tris);
 
    adjacent_tris = intersect(adjacent_tris, tris);
    adjacent_tris = subtract(adjacent_tris, forbidden_tris);
    if (adjacent_tris.size() == 2) {
 
#ifndef NDEBUG
      int side_number0, sense0, offset0;
      CHKERR mField.get_moab().side_number(adjacent_tris[0], edge, side_number0,
                                           sense0, offset0);
      int side_number1, sense1, offset1;
      CHKERR mField.get_moab().side_number(adjacent_tris[1], edge, side_number1,
                                           sense1, offset1);
      if (sense0 * sense1 > 0)
        SETERRQ(
            PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Cannot flip edge with same orientation in both adjacent faces");
#endif // NDEBUG
 
      Range new_flipped_tris;
      CHKERR make_edge_flip(edge, adjacent_tris, new_flipped_tris);
      if (new_flipped_tris.size()) {
        flipped_tris.merge(adjacent_tris);
        forbidden_tris.merge(adjacent_tris);
        new_tris.merge(new_flipped_tris);
 
#ifndef NDEBUG
        CHKERR save_range(moab,
                          "flipped_tris_" + std::to_string(flip_count) + ".vtk",
                          adjacent_tris);
        CHKERR save_range(
            moab, "new_flipped_tris_" + std::to_string(flip_count) + ".vtk",
            new_flipped_tris);
 
#endif // NDEBUG
 
        ++flip_count;
      }
    }
  }
 
  Range all_tris;
  CHKERR moab.get_entities_by_dimension(0, SPACE_DIM, all_tris);
  Range not_flipped_tris = subtract(all_tris, flipped_tris);
 
  MOFEM_LOG("SELF", Sev::noisy)
      << "Flipped " << flip_count << " edges with two adjacent faces.";
  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevel(not_flipped_tris,
                                                              child_bit);
  CHKERR mField.getInterface<BitRefManager>()->setBitRefLevel(new_tris,
                                                              child_bit);
  CHKERR mField.getInterface<BitRefManager>()->writeBitLevel(
      child_bit, BitRefLevel().set(), "edge_flips_before_refinement.vtk", "VTK",
      "");
 
  MoFEMFunctionReturn(0);
}
//! [Edge flips]
 
//! [Refine skin]
MoFEMErrorCode Example::refineSkin(BitRefLevel parent_bit,
                                   BitRefLevel child_bit) {
  MoFEMFunctionBegin;
 
  moab::Interface &moab = mField.get_moab();
  Range tris;
  CHKERR moab.get_entities_by_dimension(0, SPACE_DIM, tris);
  CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
      parent_bit, BitRefLevel().set(), tris);
 
  Skinner skin(&moab);
  Range skin_edges;
  CHKERR skin.find_skin(0, tris, false, skin_edges);      
 
  auto refine = mField.getInterface<MeshRefinement>();
  CHKERR refine->addVerticesInTheMiddleOfEdges(skin_edges, child_bit);
#ifndef NDEBUG  
  auto debug = true;
#else
  auto debug = false;
#endif
  CHKERR refine->refineTris(tris, child_bit, QUIET, debug);
 
  CHKERR mField.getInterface<BitRefManager>()->writeBitLevel(
      child_bit, BitRefLevel().set(), "edge_flips_after_refinement.vtk", "VTK",
      "");
 
  MoFEMFunctionReturn(0);
}
//! [Refine skin]
 
//! [Re-setup problem after mesh modification
MoFEMErrorCode Example::reSetupProblem(BitRefLevel child_bit) {
  MoFEMFunctionBegin;
  simpleInterface->getBitRefLevel() = child_bit;
  CHKERR simpleInterface->reSetUp(false);
  MoFEMFunctionReturn(0);
}
//! [Re-setup problem after mesh modification]
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc, &argv, NULL, help);
 
  try {
 
    //! [Register MoFEM discrete manager in PETSc]
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    //! [Register MoFEM discrete manager in PETSc]
 
    //! [Create MoAB]
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    //! [Create MoAB]
 
    //! [Create MoFEM]
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
    //! [Create MoFEM]
 
    //! [Execute DG Projection Test]
    Example ex(m_field);
    CHKERR ex.runProblem();
    //! [Execute DG Projection Test]
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
}