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 <ElasticTie.hpp>
#include <AnalyticalSolutions.hpp>
 
struct ElasticTieMeshTutorial : public ElasticTieMeshExample {
 
  using ElasticTieMeshExample::ElasticTieMeshExample;
 
  MoFEMErrorCode checkResults() override;
 
private:
  MoFEMErrorCode checkStressError(const int test);
};
 
MoFEMErrorCode ElasticTieMeshTutorial::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 ElasticTieMeshTutorial::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;
 
    ElasticTieMeshTutorial ex(m_field);
    CHKERR ex.runProblem();
  }
  CATCH_ERRORS;
 
  CHKERR MoFEM::Core::Finalize();
}