elastic_tie_mesh.cpp

Go to the documentation of this file.

/**
 * @file elastic_tie_mesh.cpp
 * @example mofem/tutorials/vec-11_elastic_tie_mesh/elastic_tie_mesh.cpp
 * @brief Elastic mesh tie example with a tie-aware partitioned mesh loader
 *
 * @copyright Anonymous authors (c) 2025 under the MIT license
 */
 
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
constexpr int SPACE_DIM = EXECUTABLE_DIMENSION;
constexpr AssemblyType A =
    (SCHUR_ASSEMBLE) ? AssemblyType::BLOCK_SCHUR : AssemblyType::PETSC;
constexpr IntegrationType I = IntegrationType::GAUSS;
 
using EntData = EntitiesFieldData::EntData;
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
using BoundaryEle =
    PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
using DomainEleOp = DomainEle::UserDataOperator;
using BoundaryEleOp = BoundaryEle::UserDataOperator;
 
struct DomainBCs {};
struct BoundaryBCs {};
 
using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<A>::LinearForm<I>;
using OpDomainRhsBCs = DomainRhsBCs::OpFlux<DomainBCs, 1, SPACE_DIM>;
using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::LinearForm<I>;
using OpBoundaryRhsBCs = BoundaryRhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::BiLinearForm<I>;
using OpBoundaryLhsBCs = BoundaryLhsBCs::OpFlux<BoundaryBCs, 1, SPACE_DIM>;
using SideEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
 
#include <ElasticExample.hpp>
#include <ElasticTie.hpp>
#include <AnalyticalSolutions.hpp>
 
namespace ElasticTie {
 
MoFEMErrorCode makeClosureMeshset(moab::Interface &moab, const Range &ents,
                                  const int max_dim, EntityHandle &meshset) {
  MoFEMFunctionBegin;
  Range closure(ents);
  for (int d = 0; d < max_dim; ++d) {
    Range adj;
    CHKERR moab.get_adjacencies(ents, d, false, adj, moab::Interface::UNION);
    closure.merge(adj);
  }
  CHKERR moab.create_meshset(MESHSET_SET, meshset);
  CHKERR moab.add_entities(meshset, closure);
  MoFEMFunctionReturn(0);
}
 
Range getPhysicalVolumeEntities(const Range &ents) {
  Range physical_ents = ents.subset_by_type(MBTET);
  physical_ents.merge(ents.subset_by_type(MBHEX));
  return physical_ents;
}
 
MoFEMErrorCode getPrismTriangleFaces(moab::Interface &moab, const Range &ents,
                                     Range &triangle_faces) {
  MoFEMFunctionBegin;
  triangle_faces.clear();
  const auto prisms = ents.subset_by_type(MBPRISM);
  if (!prisms.empty()) {
    CHKERR moab.get_adjacencies(prisms, 2, false, triangle_faces,
                                moab::Interface::UNION);
    triangle_faces = triangle_faces.subset_by_type(MBTRI);
  }
  MoFEMFunctionReturn(0);
}
} // namespace ElasticTie
 
struct ElasticTieMeshExample : public ElasticExample {
 
  using ElasticExample::ElasticExample;
 
  boost::shared_ptr<ElasticTie::CommonData> tieData = nullptr;
  Range physicalBoundaryFaces;
 
  MoFEMErrorCode readMesh() override;
  MoFEMErrorCode setupProblem() override;
  MoFEMErrorCode setUpSolver(SmartPetscObj<KSP> solver) override;
  MoFEMErrorCode kspSetUpAndSolve(SmartPetscObj<KSP> solver) override;
  MoFEMErrorCode checkResults() override;
 
private:
  MoFEMErrorCode checkStressError(const int test);
};
 
MoFEMErrorCode ElasticTieMeshExample::readMesh() {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto meshset_mng = mField.getInterface<MeshsetsManager>();
  auto comm_interface = mField.getInterface<CommInterface>();
  CHKERR simple->getOptions();
  physicalBoundaryFaces.clear();
 
  if (SPACE_DIM == 3 && mField.get_comm_size() > 1) {
    auto load_fun = [&](MoFEM::Interface &mf, const char *file_name,
                        const char *) -> MoFEMErrorCode {
      MoFEMFunctionBegin;
      auto &moab = mf.get_moab();
      CHKERR moab.load_file(file_name, 0,
                            "PARALLEL=BCAST;PARTITION=PARALLEL_PARTITION;");
 
      auto *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
      if (pcomm == nullptr)
        pcomm = new ParallelComm(&moab, mf.get_comm());
 
      const auto rank = pcomm->rank();
      Tag part_tag = pcomm->part_tag();
      Range all_ents, tagged_sets, proc_ents, off_proc_ents;
      CHKERR moab.get_entities_by_handle(0, all_ents, false);
      CHKERR moab.get_entities_by_type_and_tag(
          0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
          moab::Interface::UNION);
 
      for (auto meshset : tagged_sets) {
        int part = -1;
        CHKERR moab.tag_get_data(part_tag, &meshset, 1, &part);
        Range meshset_ents;
        CHKERR moab.get_entities_by_handle(meshset, meshset_ents, true);
        if (part == rank)
          proc_ents.merge(meshset_ents);
        else
          off_proc_ents.merge(meshset_ents);
        CHKERR moab.tag_clear_data(part_tag, meshset_ents, &part);
      }
 
      const auto all_volume_ents = all_ents.subset_by_dimension(3);
      const auto proc_volume_ents = proc_ents.subset_by_dimension(3);
      const auto all_phys_ents = ElasticTie::getPhysicalVolumeEntities(all_ents);
      const auto proc_phys_ents = ElasticTie::getPhysicalVolumeEntities(proc_ents);
 
      std::array<Range, 4> proc_ents_skin;
      Skinner skin(&moab);
      Range all_skin, proc_skin;
      CHKERR skin.find_skin(0, all_volume_ents, false, all_skin);
      CHKERR skin.find_skin(0, proc_volume_ents, false, proc_skin);
      proc_ents_skin[2] = subtract(proc_skin, all_skin);
 
      Range all_phys_skin, local_phys_faces, physical_boundary_faces;
      CHKERR skin.find_skin(0, all_phys_ents.subset_by_dimension(3), false,
                            all_phys_skin);
      CHKERR moab.get_adjacencies(proc_phys_ents, 2, false, local_phys_faces,
                                  moab::Interface::UNION);
      physical_boundary_faces = intersect(all_phys_skin, local_phys_faces);
      if (!physical_boundary_faces.empty()) {
        CHKERR moab.tag_clear_data(pcomm->partition_tag(), physical_boundary_faces,
                                   &rank);
      }
      proc_ents_skin[2].merge(physical_boundary_faces);
 
      Range prism_triangle_faces;
      CHKERR ElasticTie::getPrismTriangleFaces(moab, proc_ents, prism_triangle_faces);
      proc_ents_skin[2].merge(prism_triangle_faces);
      CHKERR moab.get_adjacencies(proc_ents_skin[2], 1, false,
                                  proc_ents_skin[1],
                                  moab::Interface::UNION);
      CHKERR moab.get_connectivity(proc_ents_skin[2], proc_ents_skin[0],
                                   false);
 
      auto to_remove = off_proc_ents;
      for (int d = 2; d >= 0; --d)
        to_remove = subtract(to_remove, proc_ents_skin[d]);
 
      Range all_meshsets;
      CHKERR moab.get_entities_by_type(0, MBENTITYSET, all_meshsets, true);
      for (auto meshset : all_meshsets)
        CHKERR moab.remove_entities(meshset, to_remove);
 
      for (int d = 3; d > 0; --d) {
        Range ents = to_remove.subset_by_dimension(d);
        if (!ents.empty())
          CHKERR moab.delete_entities(ents);
      }
 
      physicalBoundaryFaces = physical_boundary_faces;
      CHKERR pcomm->resolve_shared_ents(0, proc_ents, 3, -1,
                                        proc_ents_skin.data());
 
      MoFEMFunctionReturn(0);
    };
 
    CHKERR simple->loadFile("", "", load_fun);
  } else {
    CHKERR simple->loadFile();
  }
 
  if (SPACE_DIM == 3 && mField.get_comm_size() > 1) {
    CHKERR comm_interface->synchroniseEntities(physicalBoundaryFaces);
  }
 
  CHKERR meshset_mng->setMeshsetFromFile();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ElasticTieMeshExample::setupProblem() {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
  PetscInt choice_base_value = AINSWORTH;
  CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
                              LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
 
  FieldApproximationBase base;
  switch (choice_base_value) {
  case AINSWORTH:
    base = AINSWORTH_LEGENDRE_BASE;
    break;
  case DEMKOWICZ:
    base = DEMKOWICZ_JACOBI_BASE;
    break;
  default:
    base = LASTBASE;
    break;
  }
 
  int order = 3;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
 
  if (SPACE_DIM == 3) {
    Range domain_ents;
    CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
                                                    domain_ents, true);
    auto tet_hex_ents = ElasticTie::getPhysicalVolumeEntities(domain_ents);
 
    if (!tet_hex_ents.empty()) {
      EntityHandle domain_meshset;
      CHKERR ElasticTie::makeClosureMeshset(mField.get_moab(), tet_hex_ents, SPACE_DIM,
                                domain_meshset);
      simple->getMeshset() = domain_meshset;
 
      Range boundary_faces;
      if (mField.get_comm_size() > 1) {
        Skinner skin(&mField.get_moab());
        Range local_skin;
        CHKERR skin.find_skin(0, tet_hex_ents, false, local_skin);
        boundary_faces = intersect(physicalBoundaryFaces, local_skin);
      } else if (!physicalBoundaryFaces.empty()) {
        boundary_faces = physicalBoundaryFaces;
      } else {
        Skinner skin(&mField.get_moab());
        CHKERR skin.find_skin(0, tet_hex_ents, false, boundary_faces);
      }
      EntityHandle boundary_meshset;
      CHKERR ElasticTie::makeClosureMeshset(mField.get_moab(), boundary_faces,
                                SPACE_DIM - 1, boundary_meshset);
      simple->getBoundaryMeshSet() = boundary_meshset;
    }
  }
 
  CHKERR ElasticExample::setupProblem();
  CHKERR ElasticTie::initializeTieConstraints(mField, tieData, "U", base,
                                                 order);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ElasticTieMeshExample::setUpSolver(SmartPetscObj<KSP> solver) {
  MoFEMFunctionBegin;
  CHKERR ElasticTie::setupTieSolver(mField, tieData, solver);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
ElasticTieMeshExample::kspSetUpAndSolve(SmartPetscObj<KSP> solver) {
  MoFEMFunctionBegin;
  CHKERR ElasticExample::kspSetUpAndSolve(solver);
  CHKERR ElasticTie::logTieDisplacementNorms(mField, tieData, "U");
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ElasticTieMeshExample::checkStressError(const int test) {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
 
  auto norm_fe = boost::make_shared<DomainEle>(mField);
  norm_fe->getRuleHook = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      norm_fe->getOpPtrVector(), {H1}, "GEOMETRY");
 
  auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEle>(
      mField, norm_fe->getOpPtrVector(), "U", "MAT_ELASTIC", Sev::verbose);
 
  auto exact_stress_ptr = boost::make_shared<MatrixDouble>();
 
  constexpr double uniaxial_stress_xx = 1e3 * 0.3 / 3.0;
  auto uniaxial_stress = [](const double, const double, const double) {
    MatrixDouble stress((SPACE_DIM * (SPACE_DIM + 1)) / 2, 1);
    stress.clear();
    auto t_stress = getFTensor2SymmetricFromMat<
        SPACE_DIM, -1, DataLayoutTraits<DataLayout::CoeffsByGauss>>(stress);
    t_stress(0, 0) = uniaxial_stress_xx;
    return stress;
  };
 
  const AnalyticalSolutions::HollowCylinderUnderRadialPressure<SPACE_DIM>
      lame_solution(0.5, 1.0, 2.0, 1.0, 1e3, 0.3);
  auto lame_stress = [&lame_solution](const double x, const double y,
                                      const double z) {
    return lame_solution.stress(x, y, z);
  };
 
  MatrixFunc stress_func = (test == 11) ? MatrixFunc(lame_stress)
                                        : MatrixFunc(uniaxial_stress);
  norm_fe->getOpPtrVector().push_back(
      new OpGetTensor2SymmetricfromFunc<SPACE_DIM, SPACE_DIM>(
          exact_stress_ptr, stress_func));
 
  enum Norms { STRESS_ERROR_L2 = 0, STRESS_EXACT_L2, LAST_NORM };
  auto norms_vec =
      createVectorMPI(mField.get_comm(),
                      (mField.get_comm_rank() == 0) ? LAST_NORM : 0,
                      LAST_NORM);
  CHKERR VecZeroEntries(norms_vec);
 
  norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor2Symmetric<SPACE_DIM>(
          exact_stress_ptr, norms_vec, STRESS_EXACT_L2));
  norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor2Symmetric<SPACE_DIM>(
          common_ptr->getMatCauchyStress(), norms_vec, STRESS_ERROR_L2,
          exact_stress_ptr));
 
  CHKERR DMoFEMLoopFiniteElements(simple->getDM(), simple->getDomainFEName(),
                                  norm_fe);
  CHKERR VecAssemblyBegin(norms_vec);
  CHKERR VecAssemblyEnd(norms_vec);
 
  if (mField.get_comm_rank() == 0) {
    const double *norms;
    CHKERR VecGetArrayRead(norms_vec, &norms);
 
    const double error_l2 = std::sqrt(norms[STRESS_ERROR_L2]);
    const double exact_l2 = std::sqrt(norms[STRESS_EXACT_L2]);
 
    MOFEM_LOG_C("WORLD", Sev::inform, "STRESS_ERROR_L2 = %.16e\n", error_l2);
    MOFEM_LOG_C("WORLD", Sev::inform, "STRESS_EXACT_L2 = %.16e\n", exact_l2);
 
    CHKERR VecRestoreArrayRead(norms_vec, &norms);
 
    double max_error_l2 = 1e-2;
    CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-stress_error_l2_tol",
                               &max_error_l2, PETSC_NULLPTR);
    if (error_l2 >= max_error_l2) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Stress L2 error %.16e exceeds tolerance %.16e", error_l2,
              max_error_l2);
    }
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ElasticTieMeshExample::checkResults() {
  MoFEMFunctionBegin;
 
  int test = 0;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-test", &test, PETSC_NULLPTR);
  if (test == 10 || test == 11)
    CHKERR checkStressError(test);
  MoFEMFunctionReturn(0);
}
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  MeshsetsManager::broadcastMeshsetsOn = false;
 
  auto core_log = logging::core::get();
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmSync(), "FieldEvaluator"));
  LogManager::setLog("FieldEvaluator");
  MOFEM_LOG_TAG("FieldEvaluator", "field_eval");
 
  try {
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    DMType dm_name_mg = "DMMOFEM_MG";
    CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
 
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    ElasticTieMeshExample ex(m_field);
    CHKERR ex.runProblem();
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
}