between_meshes_dg_projection.cpp

Go to the documentation of this file.

/**
 * @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";
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);
  MoFEMErrorCode outputResults(std::string field_name);
  MoFEMErrorCode edgeFlips(BitRefLevel parent_bit, BitRefLevel child_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)
      : DomainEleOp(NOSPACE, OPSPACE), dataPtr(data_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
 
    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;
};
 
//! [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);
 
  CHKERR edgeFlips(parent_bit, child_bit);
  CHKERR reSetupProblem(BitRefLevel().set());
  CHKERR projectResults(parent_bit, child_bit);
 
  CHKERR reSetupProblem(child_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, L2,
                                         AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
 
  CHKERR simpleInterface->setFieldOrder(FIELD_NAME, 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, FIELD_NAME, beta));
 
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpDomainSource(FIELD_NAME, 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) {
  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(), child_bit, child_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(
        child_bit, child_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_filed_eval_op = [&](auto fe_child_ptr) {
    auto field_eval_ptr = mField.getInterface<FieldEvaluatorInterface>();
 
    // Get pointer of FieldEvaluator data. Note finite element and method
    // set integrating 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);
 
    // You can add more fields here
    auto data_ptr = boost::make_shared<MatrixDouble>();
    auto eval_data_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,
                                                      eval_data_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_ptr, data_ptr}}, simpleInterface->getDomainFEName(),
        field_eval_data, 0, mField.get_comm_size(), parent_bit, parent_bit,
        MF_EXIST, QUIET);
 
    fe_child_ptr->getOpPtrVector().push_back(op_ptr);    
    return data_ptr;
  };
 
  // calculate coefficients on child (flipped) elements
  auto dg_projection_base = [&](auto fe_child_ptr, auto data_ptr) {
    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>();
 
    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, there 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();
      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().begin(), coeffs_ptr->data().end(),
                  field_ent_data.begin());
      }
      MoFEMFunctionReturn(0);
    };
    fe_child_ptr->getOpPtrVector().push_back(op_set_coeffs);
 
    MoFEMFunctionReturn(0);
  };
 
  // 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_filed_eval_op, 
  // pass data_ptr values used to project and calculate coefficients
  CHKERR dg_projection_base(fe_child_ptr, get_filed_eval_op(fe_child_ptr));
 
  // That is to test, if projection works, and coefficients are set in correctly
  auto test_data_ptr = boost::make_shared<MatrixDouble>();
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpCalculateVectorFieldValues<FIELD_DIM>(FIELD_NAME, test_data_ptr));
  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpError(test_data_ptr));
 
  CHKERR pipeline_mng->loopFiniteElements();
 
  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>();
  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues(FIELD_NAME, u_ptr));
 
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldGradient<SPACE_DIM>(FIELD_NAME, grad_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(),
 
          OpPPMap::DataMapVec{{FIELD_NAME, u_ptr}},
 
          OpPPMap::DataMapMat{{"GRAD_" + std::string(FIELD_NAME), grad_u_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);
    };
 
    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;
        }
      }
    }
 
    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);
 
#ifndef NDEBUG
    CHKERR save_range(moab, "new_tris.vtk", new_tris);
#endif // NDEBUG
 
    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 adj_skin_tris;
  CHKERR moab.get_adjacencies(skin_edges, SPACE_DIM, false, adj_skin_tris);
 
  Range edges;
  CHKERR moab.get_entities_by_dimension(0, SPACE_DIM - 1, edges);
  CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
      parent_bit, BitRefLevel().set(), edges);
 
  Range new_tris, flipped_tris;
  int flip_count = 0;
  for (auto edge : edges) {
    Range adjacent_tris;
    CHKERR moab.get_adjacencies(&edge, 1, SPACE_DIM, false, adjacent_tris);
    CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
        parent_bit, BitRefLevel().set(), adjacent_tris);
    adjacent_tris = subtract(adjacent_tris, adj_skin_tris);
    adjacent_tris = subtract(adjacent_tris, flipped_tris);
    if (adjacent_tris.size() == 2) {
      CHKERR make_edge_flip(edge, adjacent_tris, new_tris);
      flipped_tris.merge(adjacent_tris);
      ++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]
 
//! [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();
}