ElasticTieMeshExample Struct Reference

Inheritance diagram for ElasticTieMeshExample:

Collaboration diagram for ElasticTieMeshExample:

Public Member Functions
MoFEMErrorCode	readMesh () override
	[Run problem]
MoFEMErrorCode	setupProblem () override
	[Read mesh]
MoFEMErrorCode	setUpSolver (SmartPetscObj< KSP > solver) override
	[Push operators to pipeline]
MoFEMErrorCode	kspSetUpAndSolve (SmartPetscObj< KSP > solver) override
MoFEMErrorCode	checkResults () override
	[Postprocess results]
	ElasticExample (MoFEM::Interface &m_field)
Public Member Functions inherited from ElasticExample
	ElasticExample (MoFEM::Interface &m_field)
virtual	~ElasticExample ()=default
MoFEMErrorCode	runProblem ()
	[Run problem]

Public Attributes
boost::shared_ptr< ElasticTie::CommonData >	tieData = nullptr
Range	physicalBoundaryFaces

Private Member Functions
MoFEMErrorCode	checkStressError (const int test)

Additional Inherited Members
Protected Member Functions inherited from ElasticExample
virtual MoFEMErrorCode	boundaryCondition ()
	[Set up problem]
virtual MoFEMErrorCode	assembleSystem ()
	[Boundary condition]
virtual MoFEMErrorCode	solveSystem ()
	[Solve]
virtual MoFEMErrorCode	outputResults ()
	[Solve]
Protected Attributes inherited from ElasticExample
MoFEM::Interface &	mField
boost::shared_ptr< MatrixDouble >	vectorFieldPtr = nullptr

Detailed Description

Examples

mofem/tutorials/vec-11_elastic_tie_mesh/elastic_tie_mesh.cpp.

Definition at line 76 of file elastic_tie_mesh.cpp.

Member Function Documentation

checkResults()

MoFEMErrorCode ElasticTieMeshExample::checkResults ( )

overridevirtual

[Postprocess results]

[Check]

Reimplemented from ElasticExample.

Examples

mofem/tutorials/vec-11_elastic_tie_mesh/elastic_tie_mesh.cpp.

Definition at line 366 of file elastic_tie_mesh.cpp.

                                                   {
  MoFEMFunctionBegin;
 
  int test = 0;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-test", &test, PETSC_NULLPTR);
  if (test == 10 || test == 11)
    CHKERR checkStressError(test);
  MoFEMFunctionReturn(0);
}

checkStressError()

MoFEMErrorCode ElasticTieMeshExample::checkStressError ( const int  test )

private

Examples

mofem/tutorials/vec-11_elastic_tie_mesh/elastic_tie_mesh.cpp.

Definition at line 280 of file elastic_tie_mesh.cpp.

                                                                     {
  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);
}

ElasticExample()

ElasticExample::ElasticExample ( MoFEM::Interface &  m_field )

inline

Definition at line 20 of file ElasticExample.hpp.

: mField(m_field) {}

kspSetUpAndSolve()

MoFEMErrorCode ElasticTieMeshExample::kspSetUpAndSolve ( SmartPetscObj< KSP >  solver )

overridevirtual

Reimplemented from ElasticExample.

Examples

mofem/tutorials/vec-11_elastic_tie_mesh/elastic_tie_mesh.cpp.

Definition at line 272 of file elastic_tie_mesh.cpp.

                                                                 {
  MoFEMFunctionBegin;
  CHKERR ElasticExample::kspSetUpAndSolve(solver);
  CHKERR ElasticTie::logTieDisplacementNorms(mField, tieData, "U");
 
  MoFEMFunctionReturn(0);
}

readMesh()

MoFEMErrorCode ElasticTieMeshExample::readMesh ( )

overridevirtual

[Run problem]

[Read mesh]

Reimplemented from ElasticExample.

Examples

mofem/tutorials/vec-11_elastic_tie_mesh/elastic_tie_mesh.cpp.

Definition at line 93 of file elastic_tie_mesh.cpp.

                                               {
  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);
}

setupProblem()

MoFEMErrorCode ElasticTieMeshExample::setupProblem ( )

overridevirtual

[Read mesh]

[Set up problem]

Reimplemented from ElasticExample.

Examples

mofem/tutorials/vec-11_elastic_tie_mesh/elastic_tie_mesh.cpp.

Definition at line 201 of file elastic_tie_mesh.cpp.

                                                   {
  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);
}

setUpSolver()

MoFEMErrorCode ElasticTieMeshExample::setUpSolver ( SmartPetscObj< KSP >  solver )

overridevirtual

[Push operators to pipeline]

Reimplemented from ElasticExample.

Examples

mofem/tutorials/vec-11_elastic_tie_mesh/elastic_tie_mesh.cpp.

Definition at line 265 of file elastic_tie_mesh.cpp.

                                                                           {
  MoFEMFunctionBegin;
  CHKERR ElasticTie::setupTieSolver(mField, tieData, solver);
  MoFEMFunctionReturn(0);
}

Member Data Documentation

physicalBoundaryFaces

Range ElasticTieMeshExample::physicalBoundaryFaces

Examples

mofem/tutorials/vec-11_elastic_tie_mesh/elastic_tie_mesh.cpp.

Definition at line 81 of file elastic_tie_mesh.cpp.

tieData

boost::shared_ptr<ElasticTie::CommonData> ElasticTieMeshExample::tieData = nullptr

Examples

mofem/tutorials/vec-11_elastic_tie_mesh/elastic_tie_mesh.cpp.

Definition at line 80 of file elastic_tie_mesh.cpp.

---

The documentation for this struct was generated from the following file:

• tutorials/vec-11_elastic_tie_mesh/elastic_tie_mesh.cpp