SolutionMapping.hpp

 
 
/** \file SolutionMapping.hpp
 * \example SolutionMapping.hpp
 */
 
#include <DGProjection.hpp>
struct SetEnrichedIntegration {
 
  SetEnrichedIntegration();
 
  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;
  };
 
  static inline std::map<long int, MatrixDouble> mapRefCoords;
};
 
struct SolutionMapping {
 
  SolutionMapping(MoFEM::Interface &m_field);
 
  MoFEMErrorCode mapFields(SmartPetscObj<DM> &intermediateDM,
                           BitRefLevel parentBit, BitRefLevel childBit,
                           PetscInt numVecs = 1,
                           Vec vecsToMapFrom[] = PETSC_NULLPTR,
                           Vec vecsToMapTo[] = PETSC_NULLPTR);
  MoFEMErrorCode reSetUp(std::vector<std::string> L2Fields,
                         std::vector<int> L2Oders,
                         std::vector<std::string> H1Fields,
                         std::vector<int> H1Orders,
                         std::vector<std::string> HdivFields,
                         std::vector<int> HdivOrders, BitRefLevel bIt);
 
private:
  MoFEM::Interface &m_field;
  BitRefLevel parentBit, childBit;
  SmartPetscObj<DM> subDM, intermediateDM;
 
  boost::shared_ptr<DomainEle> feRhs;
  boost::shared_ptr<DomainEle> feLhs;
 
  MoFEMErrorCode setUp(SmartPetscObj<DM> &subDM,
                       SmartPetscObj<DM> &intermediateDM,
                       BitRefLevel parent_bit, BitRefLevel child_bit);
  MoFEMErrorCode projectResults(BitRefLevel parent_bit, BitRefLevel child_bit);
  MoFEMErrorCode solveMap(SmartPetscObj<DM> &subDM,
                          SmartPetscObj<DM> &intermediateDM,
                          BitRefLevel parent_bit, BitRefLevel child_bit,
                          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,
                                          BitRefLevel parent_bit,
                                          BitRefLevel child_bit,
                                          PetscInt num_vecs,
                                          Vec vecs_to_map_from[],
                                          Vec vecs_to_map_to[]) {
  MoFEMFunctionBegin;
 
  CHKERR setUp(subDM, intermediate_dm, parent_bit, child_bit);
  CHKERR solveMap(subDM, intermediate_dm, parent_bit, child_bit, num_vecs,
                  vecs_to_map_from, vecs_to_map_to);
 
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode SolutionMapping::setUp(SmartPetscObj<DM> &sub_dm,
                                      SmartPetscObj<DM> &intermediate_dm,
                                      BitRefLevel parent_bit,
                                      BitRefLevel child_bit) {
  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";
 
      // get only refinement bit DOFs
      auto ref_entities_ptr = boost::make_shared<Range>();
      auto child_mask = parent_bit;
      CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
          child_bit, child_mask.flip(), SPACE_DIM, *ref_entities_ptr);
 
      Range edges;
      CHKERR m_field.get_moab().get_adjacencies(*ref_entities_ptr,
                                                SPACE_DIM - 1, true, edges,
                                                moab::Interface::UNION);
      ref_entities_ptr->merge(edges);
 
      Range verts;
      CHKERR m_field.get_moab().get_adjacencies(*ref_entities_ptr, 0, true,
                                                verts, moab::Interface::UNION);
      ref_entities_ptr->merge(verts);
 
      for (auto f : {"U", "FLUX"}) {
        CHKERR DMMoFEMAddSubFieldRow(sub_dm, f, ref_entities_ptr);
        CHKERR DMMoFEMAddSubFieldCol(sub_dm, f, ref_entities_ptr);
      }
 
      CHKERR m_field.modify_problem_ref_level_set_bit("SUB_MAPPING",
                                                      parent_bit | child_bit);
      CHKERR m_field.modify_problem_mask_ref_level_set_bit("SUB_MAPPING",
                                                           BitRefLevel().set());
 
      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> L2_fields,
                                        std::vector<int> L2_orders,
                                        std::vector<std::string> H1_fields,
                                        std::vector<int> H1_orders,
                                        std::vector<std::string> Hdiv_fields,
                                        std::vector<int> Hdiv_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 < L2_fields.size(); ++i) {
    CHKERR add_ents_order_to_field(m_field, {L2_fields[i]}, {MBTRI},
                                   L2_orders[i]);
  }
  for (size_t i = 0; i < H1_fields.size(); ++i) {
    CHKERR add_ents_order_to_field(m_field, {H1_fields[i]}, {MBTRI, MBEDGE},
                                   H1_orders[i]);
    CHKERR m_field.set_field_order(0, MBVERTEX, {H1_fields[i]},
                                   1); // H1 only supports 1st order on vertices
  }
  for (size_t i = 0; i < Hdiv_fields.size(); ++i) {
    CHKERR add_ents_order_to_field(m_field, {Hdiv_fields[i]}, {MBTRI, MBEDGE},
                                   Hdiv_orders[i]);
  }
 
  Skinner skin(&m_field.get_moab());
  Range faces;
  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
      BitRefLevel().set(FLIPPED_BIT + num_refinement_levels),
      BitRefLevel().set(), MBTRI, faces);
  Range skin_edges;
  CHKERR skin.find_skin(0, faces, false, skin_edges);
 
#ifdef ADD_CONTACT
 
  auto filter_blocks = [&](auto skin) {
    bool is_contact_block = false;
    Range contact_range;
    for (auto m :
         m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % "CONTACT").str()
 
                 ))
 
    ) {
      is_contact_block =
          true; ///< blocs interation is collective, so that is set irrespective
                ///< if there are entities in given rank or not in the block
      MOFEM_LOG("CONTACT", Sev::inform)
          << "Find contact block set:  " << m->getName();
      auto meshset = m->getMeshset();
      Range contact_meshset_range;
      CHKERR m_field.get_moab().get_entities_by_dimension(
          meshset, SPACE_DIM - 1, contact_meshset_range, true);
 
      CHKERR m_field.getInterface<CommInterface>()->synchroniseEntities(
          contact_meshset_range);
      contact_range.merge(contact_meshset_range);
    }
    if (is_contact_block) {
      MOFEM_LOG("SYNC", Sev::inform)
          << "Nb entities in contact surface: " << contact_range.size();
      MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
      skin = intersect(skin, contact_range);
    }
    return skin;
  };
 
  auto sigma_trace_ents = filter_blocks(skin_edges);
  CHKERR m_field.add_ents_to_field_by_type(0, MBTRI, "SIGMA");
  CHKERR m_field.set_field_order(0, MBTRI, "SIGMA", 0);
  CHKERR m_field.add_ents_to_field_by_type(sigma_trace_ents, MBEDGE, "SIGMA");
  CHKERR m_field.set_field_order(sigma_trace_ents, "SIGMA", order - 1);
#endif // ADD_CONTACT
 
  CHKERR m_field.add_ents_to_finite_element_by_type(0, MBTRI,
                                                    simple->getDomainFEName());
  CHKERR m_field.add_ents_to_finite_element_by_type(
      skin_edges, MBEDGE, simple->getBoundaryFEName());
 
  CHKERR m_field.build_fields();
  CHKERR m_field.build_finite_elements();
  CHKERR m_field.build_adjacencies(bit);
 
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode SolutionMapping::projectResults(BitRefLevel parent_bit,
                                               BitRefLevel child_bit) {
  MoFEMFunctionBegin;
 
  feRhs->getOpPtrVector().clear();
  feLhs->getOpPtrVector().clear();
 
  MOFEM_LOG("REMESHING", Sev::inform) << "Cleared FEs";
 
  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;
 
  // 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>(m_field, simple->getDomainFEName(), child_bit,
                                  parent_bit.flip(), Sev::noisy);
    auto fe_child_ptr = op_this_child->getThisFEPtr();
    fe_child_ptr->getRuleHook = [](int, int, int p) { return -1; };
    Range child_edges;
    CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
        child_bit, parent_bit.flip(), MBEDGE, child_edges);
    // set integration rule, such that integration points are not on flipped edge
    CHKERR setDGSetIntegrationPoints<2>( // Currently only implemented for 2D
        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 = m_field.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, simple->getDomainFEName(), parent_bit,
        BitRefLevel().set());
 
    // You can add more fields here
    auto new_T_ptr = boost::make_shared<MatrixDouble>();
    auto old_T_ptr = boost::make_shared<MatrixDouble>();
    auto new_TAU_ptr = boost::make_shared<MatrixDouble>();
    auto old_TAU_ptr = boost::make_shared<MatrixDouble>();
    auto new_EP_ptr = boost::make_shared<MatrixDouble>();
    auto old_EP_ptr = boost::make_shared<MatrixDouble>();
 
    auto new_U_ptr = boost::make_shared<MatrixDouble>();
    auto old_U_ptr = boost::make_shared<MatrixDouble>();
    auto new_flux_ptr = boost::make_shared<MatrixDouble>();
    auto old_flux_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<1>("T", old_T_ptr));
      fe_eval_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<1>("TAU", old_TAU_ptr));
      constexpr int size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
      fe_eval_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<size_symm>("EP", old_EP_ptr));
      fe_eval_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", old_U_ptr));
      fe_eval_ptr->getOpPtrVector().push_back(
          new OpCalculateHVecVectorField<3>("FLUX", old_flux_ptr));
 
    } else {
      CHK_THROW_MESSAGE(MOFEM_ATOM_TEST_INVALID,
                        "Field evaluator method pointer is expired");
    }
 
    auto op_ptr = field_eval_ptr->getDataOperator<SPACE_DIM>(
        {
 
            {old_T_ptr, new_T_ptr},
            {old_TAU_ptr, new_TAU_ptr},
            {old_EP_ptr, new_EP_ptr},
            {old_U_ptr, new_U_ptr},
            {old_flux_ptr, new_flux_ptr}
 
        },
        simple->getDomainFEName(), field_eval_data, 0, m_field.get_comm_size(),
        parent_bit, child_bit.flip(), MF_EXIST, QUIET);
 
    fe_child_ptr->getOpPtrVector().push_back(op_ptr);
 
    using FieldTupleVec = std::vector<
        std::tuple<std::string, boost::shared_ptr<MatrixDouble>, int>>;
 
    return std::make_pair(
 
        std::vector<std::pair<std::string, boost::shared_ptr<MatrixDouble>>>{
            {"T", new_T_ptr}, {"TAU", new_TAU_ptr}, {"EP", new_EP_ptr}},
 
        std::vector<std::pair<std::string, boost::shared_ptr<MatrixDouble>>>{
            {"U", new_U_ptr}, {"FLUX", new_flux_ptr}}
 
    );
  };
 
  // calculate coefficients on child (flipped) elements
  auto dg_projection_base = [&](auto fe_child_ptr, auto vec_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>();
 
    constexpr int projection_order = 2; // use order 2 for DG projection
 
    for (auto &[field_name, data_ptr] : vec_data_ptr.first) {
 
      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));
 
      // 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);
    }
 
    for (auto &p : vec_data_ptr.second) {
      auto field_name = p.first;
      auto data_ptr = p.second;
 
      if (field_name == "U") {
 
        using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
            PETSC>::BiLinearForm<GAUSS>::OpMass<1, SPACE_DIM>;
 
        // 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, field_name, beta));
        using OpVec = FormsIntegrators<DomainEleOp>::Assembly<
            PETSC>::LinearForm<GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 1>;
        fe_child_ptr->getOpPtrVector().push_back(
            new OpVec(field_name, data_ptr, beta));
      }
 
      if (field_name == "FLUX") {
 
        using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
            PETSC>::BiLinearForm<GAUSS>::OpMass<3, 3>;
 
        // 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, field_name, beta));
        using OpVec = FormsIntegrators<DomainEleOp>::Assembly<
            PETSC>::LinearForm<GAUSS>::OpBaseTimesVector<3, 3, 1>;
        fe_child_ptr->getOpPtrVector().push_back(
            new OpVec(field_name, data_ptr, beta));
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  // create child operator, and fe_child_ptr element in it
  auto fe_child_ptr = get_child_op(feRhs->getOpPtrVector());
  // 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));
 
  MOFEM_LOG("REMESHING", Sev::inform) << "Before looping FEs";
  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
      subDM, simple->getDomainFEName(), feRhs, 0, m_field.get_comm_size());
  // CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(subDM, simple->getDomainFEName(),
  //                                             feRhs, m_field.get_comm_rank(),
  //                                             m_field.get_comm_rank());
  MOFEM_LOG("REMESHING", Sev::inform) << "After looping RHS";
  // CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
  //     subDM, simple->getDomainFEName(), feLhs, 0, m_field.get_comm_size());
  // MOFEM_LOG("REMESHING", Sev::inform) << "After looping LHS";
 
  MoFEMFunctionReturn(0);
};
 
MoFEMErrorCode SolutionMapping::solveMap(
    SmartPetscObj<DM> &sub_dm, SmartPetscObj<DM> &intermediate_dm,
    BitRefLevel parent_bit, BitRefLevel child_bit, PetscInt num_vecs,
    Vec vecs_to_map_from[], Vec vecs_to_map_to[]) {
  MoFEMFunctionBegin;
 
  auto simple = m_field.getInterface<Simple>();
 
  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 U_ptr = boost::make_shared<MatrixDouble>();
    auto FLUX_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 OpCalculateHVecVectorField<3, SPACE_DIM>("FLUX", FLUX_ptr));
 
    post_proc_fe->getOpPtrVector().push_back(new OpPPMap(
        post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
        {{"T", T_ptr}, {"TAU", TAU_ptr}},
 
        {{"U", U_ptr}, {"FLUX", FLUX_ptr}},
 
        {},
 
        {{"EP", EP_ptr}}
 
        )
 
    );
 
    MOFEM_LOG("WORLD", Sev::inform) << "Just before looping post_proc_fe";
    {
      CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dm, simple->getDomainFEName(),
                                                  post_proc_fe, 0,
                                                  m_field.get_comm_size());
      auto &post_proc_moab = post_proc_fe->getPostProcMesh();
      auto output_name = "out_all_rank_" +
                         std::to_string(m_field.get_comm_rank()) + "_" +
                         file_name + ".h5m";
      CHKERR post_proc_moab.write_file(output_name.c_str());
    }
    {
      CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
          dm, simple->getDomainFEName(), post_proc_fe, m_field.get_comm_rank(),
          m_field.get_comm_rank());
      auto &post_proc_moab = post_proc_fe->getPostProcMesh();
      auto output_name = "out_only_on_rank" +
                         std::to_string(m_field.get_comm_rank()) + "_" +
                         file_name + ".h5m";
      CHKERR post_proc_moab.write_file(output_name.c_str());
    }
 
    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 vec = createDMVector(sub_dm);
  int size;
  CHKERR VecGetSize(vec, &size);
  SmartPetscObj<Mat> mat;
  if (size)
    mat = createDMMatrix(sub_dm);
  auto sol = createDMVector(sub_dm);
 
  auto ksp = createKSP(m_field.get_comm());
  CHKERR KSPSetOperators(ksp, mat, mat);
  CHKERR KSPSetFromOptions(ksp);
 
  feRhs->ksp_A = mat;
  feRhs->ksp_B = mat;
  feRhs->ksp_f = vec;
  feRhs->data_ctx =
      PetscData::CtxSetF | PetscData::CtxSetA | PetscData::CtxSetB;
 
  // 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);
 
  for (int i = 0; i < num_vecs; ++i) {
    // load data to intermediate dm from vecs_to_map_from[i]
    CHKERR VecGhostUpdateBegin(vecs_to_map_from[i], INSERT_VALUES,
                               SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(vecs_to_map_from[i], INSERT_VALUES,
                             SCATTER_FORWARD);
    CHKERR DMoFEMMeshToGlobalVector(intermediate_dm, vecs_to_map_from[i],
                                    INSERT_VALUES, SCATTER_REVERSE);
 
    if (size)
      CHKERR MatZeroEntries(mat);
    CHKERR VecZeroEntries(vec);
    CHKERR VecGhostUpdateBegin(vec, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(vec, INSERT_VALUES, SCATTER_FORWARD);
    // map
    CHKERR projectResults(parent_bit, child_bit);
 
    if (size) {
      CHKERR VecAssemblyBegin(vec);
      CHKERR VecAssemblyEnd(vec);
      CHKERR MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
      CHKERR MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
 
      CHKERR KSPSolve(ksp, vec, sol);
      CHKERR VecGhostUpdateBegin(sol, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(sol, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR DMoFEMMeshToGlobalVector(sub_dm, sol, INSERT_VALUES,
                                      SCATTER_REVERSE);
    }
 
    // update data
    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, "map_" + std::to_string(i));
  }
 
  CHKERR DMoFEMMeshToLocalVector(intermediate_dm, tmp_D, INSERT_VALUES,
                                 SCATTER_REVERSE);
 
  MOFEM_LOG("REMESHING", Sev::inform) << "Mapped thermoplastic state variables";
 
  MoFEMFunctionReturn(0);
};