SolutionMapping.hpp

 
 
/** \file SolutionMapping.hpp
 * \example SolutionMapping.hpp
 */
 
struct SetEnrichedIntegration {
 
  SetEnrichedIntegration(boost::shared_ptr<Range> new_els);
 
  MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
                            int order_col, int order_data);
 
private:
  struct Fe : public ForcesAndSourcesCore {
    using ForcesAndSourcesCore::dataOnElement;
 
  private:
    using ForcesAndSourcesCore::ForcesAndSourcesCore;
  };
 
  boost::shared_ptr<Range> newEls;
  static inline std::map<long int, MatrixDouble> mapRefCoords;
};
 
struct SolutionMapping {
 
  SolutionMapping(MoFEM::Interface &m_field);
 
  MoFEMErrorCode mapFields(SmartPetscObj<DM> &intermediateDM,
                           Range elsToRemove = Range(), Range elsToAdd = Range(),
                           Range elsToFlip = Range(), PetscInt numVecs = 1,
                           Vec vecsToMapFrom[] = PETSC_NULLPTR,
                           Vec vecsToMapTo[] = PETSC_NULLPTR);
  MoFEMErrorCode reSetUp(std::vector<std::string> fIelds,
                         std::vector<int> oRders, BitRefLevel bIt);
 
private:
  MoFEM::Interface &m_field;
  Range elsToRemove;
  Range elsToAdd;
  Range elsToFlip;
  SmartPetscObj<DM> subDM, intermediateDM;
 
  boost::shared_ptr<DomainEle> feRhs;
  boost::shared_ptr<DomainEle> feLhs;
 
  MoFEMErrorCode setUp(SmartPetscObj<DM> &subDM,
                       SmartPetscObj<DM> &intermediateDM, Range elsToAdd,
                       Range elsToFlip);
  MoFEMErrorCode removeDofs(Range elsToRemove, Range elsToAdd);
  MoFEMErrorCode pushOperators(Range elsToFlip, Range elsToAdd,
                               Vec vecIn = PETSC_NULLPTR);
  MoFEMErrorCode solveMap(SmartPetscObj<DM> &subDM,
                          SmartPetscObj<DM> &intermediateDM, PetscInt numVecs,
                          Vec vecsToMapFrom[], Vec vecsToMapTo[]);
  MoFEMErrorCode postProcess();
};
 
SolutionMapping::SolutionMapping(MoFEM::Interface &m_field) : m_field(m_field) {
  feRhs = boost::make_shared<DomainEle>(m_field);
  feLhs = boost::make_shared<DomainEle>(m_field);
}
 
MoFEMErrorCode SolutionMapping::mapFields(SmartPetscObj<DM> &intermediate_dm,
                                          Range els_to_remove, Range els_to_add,
                                          Range els_to_flip, PetscInt num_vecs,
                                          Vec vecs_to_map_from[],
                                          Vec vecs_to_map_to[]) {
  MoFEMFunctionBegin;
 
  CHKERR pushOperators(els_to_flip, els_to_add);
  CHKERR setUp(subDM, intermediate_dm, els_to_add, els_to_flip);
  CHKERR solveMap(subDM, intermediate_dm, num_vecs, vecs_to_map_from,
                  vecs_to_map_to);
 
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode SolutionMapping::setUp(SmartPetscObj<DM> &sub_dm,
                                      SmartPetscObj<DM> &intermediate_dm,
                                      Range els_to_add, Range els_to_flip) {
  MoFEMFunctionBegin;
 
  MOFEM_LOG("REMESHING", Sev::inform) << "Set up sub DM for mapping";
 
  Simple *simple = m_field.getInterface<Simple>();
 
  Range sub_dm_elems;
 
  auto create_sub_dm = [&]() {
    MoFEMFunctionBegin;
    sub_dm = createDM(m_field.get_comm(), "DMMOFEM");
    auto create_impl = [&]() {
      MoFEMFunctionBegin;
      CHKERR DMSetType(sub_dm, "DMMOFEM");
      CHKERR DMMoFEMCreateSubDM(sub_dm, intermediate_dm, "SUB_MAPPING");
      CHKERR DMMoFEMAddElement(sub_dm, simple->getDomainFEName());
      CHKERR DMMoFEMSetSquareProblem(sub_dm, PETSC_TRUE);
      CHKERR DMMoFEMSetDestroyProblem(sub_dm, PETSC_TRUE);
 
      MOFEM_LOG("REMESHING", Sev::inform) << "Created sub DM";
 
      for (auto f : {"T", "TAU", "EP"}) {
        CHKERR DMMoFEMAddSubFieldRow(sub_dm, f,
                                     boost::make_shared<Range>(els_to_add));
        CHKERR DMMoFEMAddSubFieldCol(sub_dm, f,
                                     boost::make_shared<Range>(els_to_add));
      }
 
      auto bit = BitRefLevel().set();
      CHKERR m_field.modify_problem_ref_level_set_bit("SUB_MAPPING", bit);
      CHKERR m_field.modify_problem_mask_ref_level_set_bit("SUB_MAPPING", bit);
 
      CHKERR DMMoFEMSetIsPartitioned(sub_dm, is_distributed_mesh);
      CHKERR DMSubDMSetUp_MoFEM(sub_dm);
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR create_impl();
    MoFEMFunctionReturn(0);
  };
 
  CHK_THROW_MESSAGE(create_sub_dm(), "failed to create sub dm");
 
  MOFEM_LOG("REMESHING", Sev::inform) << "Set up sub DM for mapping <= done";
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode SolutionMapping::reSetUp(std::vector<std::string> fields,
                                        std::vector<int> orders,
                                        BitRefLevel bit) {
  MoFEMFunctionBegin;
 
  Simple *simple = m_field.getInterface<Simple>();
 
  auto add_ents_order_to_field =
      [&](MoFEM::Interface &m_field,
          const std::initializer_list<std::string> &f,
          const std::initializer_list<EntityType> &t, int _order) {
        MoFEMFunctionBegin;
        for (auto _f : f) {
          for (auto _t : t) {
            CHKERR m_field.add_ents_to_field_by_type(0, _t, _f);
            CHKERR m_field.set_field_order(0, _t, _f, _order);
          }
        }
        MoFEMFunctionReturn(0);
      };
 
  for (size_t i = 0; i < fields.size(); ++i) {
    CHKERR add_ents_order_to_field(m_field, {fields[i]}, {MBTRI}, orders[i]);
  }
 
  CHKERR m_field.build_fields();
  CHKERR m_field.add_ents_to_finite_element_by_type(0, MBTRI,
                                                    simple->getDomainFEName());
  CHKERR m_field.build_finite_elements();
  CHKERR m_field.build_adjacencies(bit);
 
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode SolutionMapping::removeDofs(Range els_to_remove,
                                           Range els_to_add) {
  MoFEMFunctionBegin;
 
  const Problem *prb_ptr;
  CHKERR m_field.get_problem("SUB_MAPPING", &prb_ptr);
  int nb_init_row_dofs = prb_ptr->getNbDofsRow();
 
  MOFEM_LOG("REMESHING", Sev::verbose)
      << "Number of initial row dofs: " << nb_init_row_dofs;
 
  if (is_distributed_mesh == PETSC_FALSE)
    CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
        els_to_remove);
 
  auto prb_mng = m_field.getInterface<ProblemsManager>();
 
  auto remove_dofs_on_ents = [&](std::string field, const Range &ents) {
    if (is_distributed_mesh == PETSC_TRUE)
      return prb_mng->removeDofsOnEntities("SUB_MAPPING", field, ents);
    else
      return prb_mng->removeDofsOnEntitiesNotDistributed("SUB_MAPPING", field,
                                                         ents);
  };
 
  for (std::string field : {"EP", "TAU", "T"}) {
    remove_dofs_on_ents(field, els_to_remove);
 
    nb_init_row_dofs = prb_ptr->getNbDofsRow();
 
    MOFEM_LOG("REMESHING", Sev::verbose)
        << "Number of row dofs after removing " << field
        << " field: " << nb_init_row_dofs;
  }
 
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode SolutionMapping::pushOperators(Range els_to_flip,
                                              Range els_to_add, Vec vec_in) {
  MoFEMFunctionBegin;
 
  SmartPetscObj<Vec> smart_vec_in;
  if (vec_in)
    smart_vec_in = SmartPetscObj<Vec>(vec_in, true);
 
  MOFEM_LOG("REMESHING", Sev::inform) << "Got past smart vector";
 
  feRhs->getOpPtrVector().clear();
 
  MOFEM_LOG("REMESHING", Sev::inform) << "Cleared RHS FE";
 
  Simple *simple = m_field.getInterface<Simple>();
 
  auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order; };
 
  // feRhs->getRuleHook = vol_rule;
  // feLhs->getRuleHook = vol_rule;
  feRhs->getRuleHook = [](int, int, int) { return -1; };
  feRhs->setRuleHook =
      SetEnrichedIntegration(boost::make_shared<Range>(els_to_add));
  feLhs->getRuleHook = [](int, int, int) { return -1; };
  feLhs->setRuleHook =
      SetEnrichedIntegration(boost::make_shared<Range>(els_to_add));
 
  auto T_ptr = boost::make_shared<VectorDouble>();
  auto TAU_ptr = boost::make_shared<VectorDouble>();
  auto EP_ptr = boost::make_shared<MatrixDouble>();
 
  auto set_domain_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {L2});
 
    auto old_T_ptr = boost::make_shared<VectorDouble>();
    auto old_TAU_ptr = boost::make_shared<VectorDouble>();
    auto old_EP_ptr = boost::make_shared<MatrixDouble>();
 
    auto old_el_ops = [&](auto &pip) {
      MoFEMFunctionBegin;
 
      // get old solution
      if (vec_in) {
        pip.push_back(
            new OpCalculateScalarFieldValues("T", old_T_ptr, smart_vec_in));
        pip.push_back(
            new OpCalculateScalarFieldValues("TAU", old_TAU_ptr, smart_vec_in));
        pip.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
            "EP", old_EP_ptr, smart_vec_in));
      } else {
        pip.push_back(new OpCalculateScalarFieldValues("T", old_T_ptr));
        pip.push_back(new OpCalculateScalarFieldValues("TAU", old_TAU_ptr));
        pip.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
            "EP", old_EP_ptr));
      }
      // pip.push_back(new EdgeFlippingOps::filterScalarSolution<AT, IT>(
      //     "T", old_T_ptr, nullptr, nullptr, nullptr));
      // pip.push_back(new EdgeFlippingOps::filterScalarSolution<AT, IT>(
      //     "TAU", old_TAU_ptr, nullptr, nullptr, nullptr));
      // pip.push_back(new EdgeFlippingOps::filterSymTensorSolution<AT, IT,
      // SPACE_DIM>("EP", old_EP_ptr, nullptr, nullptr, nullptr));
 
      MoFEMFunctionReturn(0);
    };
 
    auto new_el_ops = [&](auto &pp_fe, auto &&p) {
      MoFEMFunctionBegin;
 
      auto new_T_ptr = boost::make_shared<VectorDouble>();
      auto new_TAU_ptr = boost::make_shared<VectorDouble>();
      auto new_EP_ptr = boost::make_shared<MatrixDouble>();
      if (vec_in) {
        pip.push_back(
            new OpCalculateScalarFieldValues("T", new_T_ptr, smart_vec_in));
        pip.push_back(
            new OpCalculateScalarFieldValues("TAU", new_TAU_ptr, smart_vec_in));
        pip.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
            "EP", new_EP_ptr, smart_vec_in));
      } else {
        MOFEM_LOG("REMESHING", Sev::inform) << "In other branch of child ops";
        pip.push_back(new OpCalculateScalarFieldValues("T", new_T_ptr));
        pip.push_back(new OpCalculateScalarFieldValues("TAU", new_TAU_ptr));
        pip.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
            "EP", new_EP_ptr));
      }
 
      pip.push_back(
          new EdgeFlippingOps::OpRhsScalarLeastSquaresProj<AT, IT, SPACE_DIM>(
              "T", old_T_ptr, nullptr, new_T_ptr, nullptr, m_field));
      pip.push_back(
          new EdgeFlippingOps::OpRhsScalarLeastSquaresProj<AT, IT, SPACE_DIM>(
              "TAU", old_TAU_ptr, nullptr, new_TAU_ptr, nullptr, m_field));
      pip.push_back(
          new EdgeFlippingOps::OpRhsSymTensorLeastSquaresProj<AT, IT,
                                                              SPACE_DIM>(
              "EP", old_EP_ptr, nullptr, new_EP_ptr, nullptr, m_field));
 
      MoFEMFunctionReturn(0);
    };
 
    // create an OpLoopRange object
    auto op_range = new OpLoopRange<DomainEleOnRange>(
        m_field, simple->getDomainFEName(),
        boost::make_shared<Range>(els_to_flip), Sev::noisy);
 
    op_range->getRangeFEPtr()->getRuleHook = [](int, int, int) { return -1; };
 
    // push operator
    pip.push_back(op_range);
 
    // call the parent ops from each child pipeline
    CHK_MOAB_THROW(new_el_ops(pip, old_el_ops(op_range->getOpPtrVector())),
                   "failed to loop over child ops");
 
    MoFEMFunctionReturn(0);
  };
 
  auto set_domain_lhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {L2});
 
    using OpLhsScalarLeastSquaresProj = FormsIntegrators<DomainEleOp>::Assembly<
        AT>::BiLinearForm<IT>::OpMass<1, 1>;
    pip.push_back(new OpLhsScalarLeastSquaresProj("T", "T"));
    pip.push_back(new OpLhsScalarLeastSquaresProj("TAU", "TAU"));
    pip.push_back(
        new EdgeFlippingOps::OpLhsSymTensorLeastSquaresProj<SPACE_DIM, IT>(
            "EP", "EP"));
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_domain_rhs(feRhs->getOpPtrVector());
  CHKERR set_domain_lhs(feLhs->getOpPtrVector());
 
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode SolutionMapping::solveMap(SmartPetscObj<DM> &sub_dm,
                                         SmartPetscObj<DM> &intermediate_dm,
                                         PetscInt num_vecs,
                                         Vec vecs_to_map_from[],
                                         Vec vecs_to_map_to[]) {
  MoFEMFunctionBegin;
 
  auto simple = m_field.getInterface<Simple>();
  auto is_mgr = m_field.getInterface<ISManager>();
 
  auto post_proc = [&](auto dm, std::string file_name) {
    MoFEMFunctionBegin;
 
    auto T_ptr = boost::make_shared<VectorDouble>();
    auto TAU_ptr = boost::make_shared<VectorDouble>();
    auto EP_ptr = boost::make_shared<MatrixDouble>();
 
    auto post_proc_fe = boost::make_shared<PostProcEle>(m_field);
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("T", T_ptr));
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateScalarFieldValues("TAU", TAU_ptr));
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>("EP", EP_ptr));
    // post_proc_fe->getOpPtrVector().push_back(
    //     new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(post_proc_fe->getPostProcMesh(),
                    post_proc_fe->getMapGaussPts(),
 
                    {{"T", T_ptr}, {"TAU", TAU_ptr}},
 
                    {},
 
                    {},
 
                    {{"EP", EP_ptr}}
 
                    )
 
    );
 
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                    post_proc_fe);
    CHKERR post_proc_fe->writeFile(file_name);
 
    MoFEMFunctionReturn(0);
  };
 
  auto null_fe = boost::shared_ptr<FEMethod>();
 
  CHKERR DMMoFEMKSPSetComputeOperators(sub_dm, simple->getDomainFEName(), feLhs,
                                       null_fe, null_fe);
  // CHKERR DMMoFEMKSPSetComputeRHS(sub_dm, simple->getDomainFEName(), feRhs,
  //                                null_fe, null_fe);
 
  auto ksp = MoFEM::createKSP(m_field.get_comm());
  CHKERR KSPSetDM(ksp, sub_dm);
  CHKERR KSPSetFromOptions(ksp);
 
  // We create a temp_D to store current state of the mesh. We doing to preserve
  // whatever is currently stored on entities, we will restore it to keep it as
  // it was before.
 
  auto tmp_D = createDMVector(intermediate_dm);
  CHKERR DMoFEMMeshToLocalVector(intermediate_dm, tmp_D, INSERT_VALUES,
                                 SCATTER_FORWARD);
  CHKERR VecGhostUpdateBegin(tmp_D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(tmp_D, INSERT_VALUES, SCATTER_FORWARD);
 
  auto D = createDMVector(sub_dm);
  auto F = vectorDuplicate(D);
 
  for (int i = 0; i < num_vecs; ++i) {
    CHKERR VecGhostUpdateBegin(vecs_to_map_from[i], INSERT_VALUES,
                               SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(vecs_to_map_from[i], INSERT_VALUES,
                             SCATTER_FORWARD);
    CHKERR DMoFEMMeshToLocalVector(intermediate_dm, vecs_to_map_from[i],
                                   INSERT_VALUES, SCATTER_REVERSE);
    // map
    CHKERR VecZeroEntries(F);
    CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
    feRhs->ksp_f = F;
    CHKERR DMoFEMLoopFiniteElements(sub_dm, simple->getDomainFEName(), feRhs);
    CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecAssemblyBegin(F);
    CHKERR VecAssemblyEnd(F);
    CHKERR KSPSolve(ksp, F, D);
 
    // update data
    CHKERR DMoFEMMeshToLocalVector(sub_dm, D, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR DMoFEMMeshToLocalVector(intermediate_dm, vecs_to_map_to[i],
                                   INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateBegin(vecs_to_map_to[i], INSERT_VALUES,
                               SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(vecs_to_map_to[i], INSERT_VALUES, SCATTER_FORWARD);
 
    CHKERR post_proc(intermediate_dm, "out_map_" + std::to_string(i) + ".h5m");
  }
 
  CHKERR DMoFEMMeshToLocalVector(intermediate_dm, tmp_D, INSERT_VALUES,
                                 SCATTER_REVERSE);
 
  MOFEM_LOG("REMESHING", Sev::inform) << "Mapped thermoplastic state variables";
 
  MoFEMFunctionReturn(0);
};
 
// Use example
 
// CHKERR SolutionMapping(m_field).mapFields();