ElasticPostProc.hpp

Go to the documentation of this file.

#ifndef __ELASTIC_POST_PROC_HPP__
#define __ELASTIC_POST_PROC_HPP__
 
#include <HookeOps.hpp>
 
namespace ElasticPostProc {
 
template <int SPACE_DIM, typename DomainEle, typename BoundaryEle,
          typename SideEle>
MoFEMErrorCode postProcessElasticResults(
    MoFEM::Interface &mField, SmartPetscObj<DM> dm,
    const std::string &domain_fe_name, const std::string &out_file_name,
    std::vector<std::pair<std::string, SmartPetscObj<Vec>>> extra_vectors = {},
    const std::vector<std::string> &tags_to_transfer = {},
    const Sev hooke_ops_sev = Sev::verbose) {
  using namespace MoFEM;
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>();
 
  using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
  using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
  pip->getDomainRhsFE().reset();
  pip->getDomainLhsFE().reset();
  pip->getBoundaryRhsFE().reset();
  pip->getBoundaryLhsFE().reset();
 
  auto post_proc_mesh = boost::make_shared<moab::Core>();
  auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
      mField, post_proc_mesh);
  auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
      mField, post_proc_mesh);
 
  using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
 
  auto calculate_stress_ops = [&](auto &ops) {
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(ops, {H1}, "GEOMETRY");
 
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEle>(
        mField, ops, "U", "MAT_ELASTIC", hooke_ops_sev);
 
    auto u_ptr = boost::make_shared<MatrixDouble>();
    ops.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
    auto x_ptr = boost::make_shared<MatrixDouble>();
    ops.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
 
    DataMapMat vec_map = {{"U", u_ptr}, {"GEOMETRY", x_ptr}};
    if (!extra_vectors.empty()) {
      for (const auto &vec_pair : extra_vectors) {
        auto extra_ptr = boost::make_shared<MatrixDouble>();
        ops.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
            "U", extra_ptr, vec_pair.second));
        vec_map[vec_pair.first] = extra_ptr;
      }
    }
 
    DataMapMat sym_map = {{"STRAIN", common_ptr->getMatStrain()},
                          {"STRESS", common_ptr->getMatCauchyStress()}};
 
    return boost::make_tuple(common_ptr->getMatCauchyStress(), vec_map, sym_map);
  };
 
  auto get_tags_on_pmesh = [&](bool post_proc_skin) {
    int def_val_int = 0;
    Tag tag_mat;
    CHKERR mField.get_moab().tag_get_handle(
        "MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
        MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
 
    auto meshset_vec_ptr =
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
            std::regex((boost::format("%s(.*)") % "MAT_ELASTIC").str()));
 
    Range skin_ents;
    std::unique_ptr<Skinner> skin_ptr;
    if (post_proc_skin) {
      skin_ptr = std::make_unique<Skinner>(&mField.get_moab());
      auto boundary_meshset = simple->getBoundaryMeshSet();
      CHKERR mField.get_moab().get_entities_by_handle(boundary_meshset,
                                                      skin_ents, true);
    }
 
    for (auto m : meshset_vec_ptr) {
      Range ents_3d;
      CHKERR mField.get_moab().get_entities_by_handle(m->getMeshset(), ents_3d,
                                                      true);
      int const id = m->getMeshsetId();
      ents_3d = ents_3d.subset_by_dimension(SPACE_DIM);
      if (post_proc_skin) {
        Range skin_faces;
        CHKERR skin_ptr->find_skin(0, ents_3d, false, skin_faces);
        ents_3d = intersect(skin_ents, skin_faces);
      }
      CHKERR mField.get_moab().tag_clear_data(tag_mat, ents_3d, &id);
    }
 
    std::vector<Tag> tags{tag_mat};
    for (const auto &tag_name : tags_to_transfer) {
      Tag tag = 0;
      if (mField.get_moab().tag_get_handle(tag_name.c_str(), tag) ==
          MB_SUCCESS) {
        tags.push_back(tag);
      }
    }
    return tags;
  };
 
  auto post_proc_domain = [&](auto pmesh) {
    auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField, pmesh);
    auto [mat_stress_ptr, vec_map, sym_map] =
        calculate_stress_ops(post_proc_fe->getOpPtrVector());
    static_cast<void>(mat_stress_ptr);
 
    post_proc_fe->getOpPtrVector().push_back(new OpPPMap(
        post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(), {},
        vec_map, {}, sym_map));
 
    post_proc_fe->setTagsToTransfer(get_tags_on_pmesh(false));
    return post_proc_fe;
  };
 
  auto post_proc_boundary = [&](auto pmesh) {
    auto post_proc_fe = boost::make_shared<PostProcEleBdy>(mField, pmesh);
    AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
 
    auto op_loop_side =
        new OpLoopSide<SideEle>(mField, domain_fe_name, SPACE_DIM);
    auto [mat_stress_ptr, vec_map, sym_map] =
        calculate_stress_ops(op_loop_side->getOpPtrVector());
    post_proc_fe->getOpPtrVector().push_back(op_loop_side);
 
    auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
    post_proc_fe->getOpPtrVector().push_back(
        new ElasticOps::OpCalculateTraction(mat_stress_ptr, mat_traction_ptr));
    vec_map["T"] = mat_traction_ptr;
 
    post_proc_fe->getOpPtrVector().push_back(new OpPPMap(
        post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(), {},
        vec_map, {}, sym_map));
 
    post_proc_fe->setTagsToTransfer(get_tags_on_pmesh(true));
    return post_proc_fe;
  };
 
  PetscBool post_proc_skin_only = PETSC_FALSE;
  if (SPACE_DIM == 3) {
    post_proc_skin_only = PETSC_TRUE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin_only",
                               &post_proc_skin_only, PETSC_NULLPTR);
  }
 
  if (post_proc_skin_only == PETSC_FALSE) {
    pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
  } else {
    pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
  }
 
  CHKERR DMoFEMPreProcessFiniteElements(dm, post_proc_begin->getFEMethod());
  CHKERR pip->loopFiniteElements();
  CHKERR DMoFEMPostProcessFiniteElements(dm, post_proc_end->getFEMethod());
 
  CHKERR post_proc_end->writeFile(out_file_name);
  MoFEMFunctionReturn(0);
}
 
} // namespace ElasticPostProc
 
#endif // __ELASTIC_POST_PROC_HPP__