ElasticAdaptiveExample Struct Reference

#include "tutorials/vec-8_elasticity_adaptive/src/ElasticAdaptiveExample.hpp"

Inheritance diagram for ElasticAdaptiveExample:

Collaboration diagram for ElasticAdaptiveExample:

Public Member Functions
	ElasticAdaptiveExample (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	[Run problem]
Public Member Functions inherited from ElasticExample
	ElasticExample (MoFEM::Interface &m_field)
virtual	~ElasticExample ()=default
MoFEMErrorCode	runProblem ()
	[Run problem]

Private Member Functions
MoFEMErrorCode	readMesh ()
	[Run problem]
MoFEMErrorCode	setupAdaptivity ()
	[Read mesh]
MoFEMErrorCode	computeErrorNorms ()
	[Postprocess results]
MoFEMErrorCode	outputResults (int ref_iter=-1)
	[Postprocess results]
MoFEMErrorCode	refineOrder (int ref_level)
MoFEMErrorCode	computeErrorIndicators ()
MoFEMErrorCode	refineGeometry ()
MoFEMErrorCode	calculateGeometryError ()
	Calculate geometry error and save on a tag.
MoFEMErrorCode	kspSetUpAndSolve (SmartPetscObj< KSP > solver) override
	[Set up problem]

Private Attributes
int	oRder = 2
PetscBool	refGeom = PETSC_FALSE
int	refAdaptNum = 0
double	refAdaptThreshold = 1
double	refGeomThreshold = 1
std::array< double, 2 >	meanError = {0.0, 0.0}
double	errorIndicMean = 0.0
int	atomTest = 0

Additional Inherited Members
Protected Types inherited from ElasticExample
using	Pipeline = boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator >
Protected Member Functions inherited from ElasticExample
virtual MoFEMErrorCode	setupProblem ()
	[Read mesh]
virtual MoFEMErrorCode	boundaryCondition ()
	[Set up problem]
virtual MoFEMErrorCode	assembleSystem ()
virtual MoFEMErrorCode	solveSystem ()
	[Solve]
virtual MoFEMErrorCode	outputResults ()
	[Solve]
virtual MoFEMErrorCode	checkResults ()
	[Postprocess results]
virtual MoFEMErrorCode	addDomainMaterialOps (Pipeline &lhs_pipeline, Pipeline &rhs_pipeline)
	Add constitutive operators to the domain assembly pipelines.
virtual MoFEMErrorCode	setUpSolver (SmartPetscObj< KSP > solver)
	[Push operators to pipeline]
Protected Attributes inherited from ElasticExample
MoFEM::Interface &	mField
boost::shared_ptr< MatrixDouble >	vectorFieldPtr = nullptr

Detailed Description

Definition at line 12 of file ElasticAdaptiveExample.hpp.

Constructor & Destructor Documentation

ElasticAdaptiveExample()

ElasticAdaptiveExample::ElasticAdaptiveExample ( MoFEM::Interface &  m_field )

inline

Definition at line 14 of file ElasticAdaptiveExample.hpp.

: ElasticExample(m_field) {}

Member Function Documentation

calculateGeometryError()

MoFEMErrorCode ElasticAdaptiveExample::calculateGeometryError ( )

private

Calculate geometry error and save on a tag.

Definition at line 14 of file HighOrderGeometryRef.hpp.

                                                              {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  double def_val_tag = 0.0;
  Tag tag_err;
  CHKERR mField.get_moab().tag_get_handle("GEOMETRY_ERROR", 1, MB_TYPE_DOUBLE,
                                          tag_err, MB_TAG_CREAT | MB_TAG_SPARSE,
                                          &def_val_tag);
 
  auto LoCoordsPtr = boost::make_shared<MatrixDouble>();
  auto *op_ptr = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
  op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
                              EntityType type,
                              EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
    auto op_ptr = static_cast<BoundaryEleOp *>(base_op_ptr);
    *LoCoordsPtr = op_ptr->getCoordsAtGaussPts();
    MoFEMFunctionReturn(0);
  };
 
  auto *op_ptr_ho = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
  op_ptr_ho->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
                                 EntityType type,
                                 EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
    auto op_ptr_ho = static_cast<BoundaryEleOp *>(base_op_ptr);
    auto fe_ent = op_ptr_ho->getFEEntityHandle();
    FTensor::Index<'i', SPACE_DIM> i;
    MatrixDouble &lo_mat = *LoCoordsPtr;
    MatrixDouble &ho_mat = op_ptr_ho->getCoordsAtGaussPts();
    auto t_lo_coords = getFTensor1FromMat<SPACE_DIM>(lo_mat);
    auto t_ho_coords = getFTensor1FromMat<SPACE_DIM>(ho_mat);
    auto t_w = op_ptr_ho->getFTensor0IntegrationWeight();
    double err_int = 0;
    const auto nb_gauss_pts = op_ptr_ho->getGaussPts().size2();
    // E = sqrt( int (x-xh)*(x-xh) da )  / face_area
    for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
      err_int += t_w * ((t_lo_coords(i) - t_ho_coords(i)) *
                        (t_lo_coords(i) - t_ho_coords(i)));
      ++t_lo_coords;
      ++t_ho_coords;
      ++t_w;
    }
    const double err_measure = sqrt(err_int) / op_ptr_ho->getMeasure();
    if (err_measure > std::numeric_limits<float>::epsilon()) {
      meanError[IndicInds::ERROR_TOT] += err_measure;
      meanError[IndicInds::COUNT_ENTS] += 1.0;
    }
 
    CHKERR mField.get_moab().tag_set_data(tag_err, &fe_ent, 1, &err_measure);
    Range adj_tet;
    CHKERR mField.get_moab().get_adjacencies(&fe_ent, 1, SPACE_DIM, false,
                                             adj_tet);
    if (adj_tet.size() != 1)
      CHK_THROW_MESSAGE(
          MOFEM_DATA_INCONSISTENCY,
          "There should be only one ent adjacent to the boundary ent");
    auto tet_ent = *adj_tet.begin();
    double tet_err_measure = 0;
    CHKERR mField.get_moab().tag_get_data(tag_err, &tet_ent, 1,
                                          &tet_err_measure);
    if (err_measure > tet_err_measure) {
      CHKERR mField.get_moab().tag_set_data(tag_err, &tet_ent, 1, &err_measure);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto fe = boost::make_shared<BoundaryEle>(mField);
  fe->getOpPtrVector().push_back(op_ptr);
 
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      fe->getOpPtrVector(), {NOSPACE}, "GEOMETRY");
  fe->getOpPtrVector().push_back(op_ptr_ho);
  auto dm = simple->getDM();
  CHK_THROW_MESSAGE(
      DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(), fe),
      "calc metric");
 
  MoFEMFunctionReturn(0);
}

computeErrorIndicators()

MoFEMErrorCode ElasticAdaptiveExample::computeErrorIndicators ( )

private

Definition at line 230 of file AdaptiveOrderRef.hpp.

                                                              {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
  auto simple = mField.getInterface<Simple>();
 
  Tag th_error_ind;
  CHKERR getTagHandle(mField, "SOLUTION_ERROR", MB_TYPE_DOUBLE, th_error_ind);
 
  auto integration_rule = [](int, int, int p_data) { return 2 * p_data + 1; };
  CHKERR pipeline_mng->setSkeletonRhsIntegrationRule(integration_rule);
  auto &skeleton_pipeline = pipeline_mng->getOpSkeletonRhsPipeline();
 
  if (skeleton_pipeline.empty()) {
 
    auto op_loop_side =
        new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        op_loop_side->getOpPtrVector(), {H1}, "GEOMETRY");
 
    // Use HookeOps commonDataFactory to get matStrainPtr and matCauchyStressPtr
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEle>(
        mField, op_loop_side->getOpPtrVector(), "U", "MAT_ELASTIC",
        Sev::verbose);
 
    auto mat_jump_ptr = boost::make_shared<MatrixDouble>();
    op_loop_side->getOpPtrVector().push_back(new OpCalculateTractionJump(
        "U", common_ptr->getMatCauchyStress(), mat_jump_ptr));
 
    auto *op_ptr = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
    op_ptr->doWorkRhsHook = [this, th_error_ind,
                             mat_jump_ptr](DataOperator *base_op_ptr, int side,
                                           EntityType type,
                                           EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
      FTensor::Index<'i', SPACE_DIM> i;
      auto op_ptr = static_cast<BoundaryEleOp *>(base_op_ptr);
      auto fe_ent = op_ptr->getFEEntityHandle();
 
      auto t_jump = getFTensor1FromMat<SPACE_DIM>(*mat_jump_ptr);
      auto t_w = op_ptr->getFTensor0IntegrationWeight();
      double err_int = 0;
      const auto nb_gauss_pts = op_ptr->getGaussPts().size2();
      for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
        err_int += t_w * (t_jump(i) * t_jump(i));
        ++t_jump;
        ++t_w;
      }
      double err_measure = sqrt(err_int);
      CHKERR mField.get_moab().tag_set_data(th_error_ind, &fe_ent, 1,
                                            &err_measure);
 
      MoFEMFunctionReturn(0);
    };
 
    skeleton_pipeline.push_back(op_loop_side);
    skeleton_pipeline.push_back(op_ptr);
  }
 
  Range fe_ents;
  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM - 1, fe_ents);
  double zero = 0;
  CHKERR mField.get_moab().tag_clear_data(th_error_ind, fe_ents, &zero);
 
  CHKERR DMoFEMLoopFiniteElements(simple->getDM(), simple->getSkeletonFEName(),
                                  pipeline_mng->getSkeletonRhsFE());
 
  Range domain_ents_part;
  CHKERR getDomainEntsPart(mField, domain_ents_part);
 
  ParallelComm *pcomm =
      ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
  Range empty_range;
  CHKERR pcomm->reduce_tags(th_error_ind, MPI_SUM, empty_range);
 
  std::array<double, 3> error_indic_loc = {0.0, 0.0, 0.0};
 
  for (auto domain_ent : domain_ents_part) {
    Range face_ents;
    CHKERR mField.get_moab().get_adjacencies(
        &domain_ent, 1, SPACE_DIM - 1, true, face_ents, moab::Interface::UNION);
 
    double err_indic_sum = 0;
    for (auto face_ent : face_ents) {
      double err_indic = 0;
      CHKERR mField.get_moab().tag_get_data(th_error_ind, &face_ent, 1,
                                            &err_indic);
      err_indic_sum += err_indic;
    }
    CHKERR mField.get_moab().tag_set_data(th_error_ind, &domain_ent, 1,
                                          &err_indic_sum);
 
    error_indic_loc[IndicInds::ERROR_TOT] += err_indic_sum;
 
    const EntityHandle *vert_conn;
    int vert_num;
    CHKERR mField.get_moab().get_connectivity(domain_ent, vert_conn, vert_num,
                                              true);
    std::vector<double> vpos(3 * vert_num);
    CHKERR mField.get_moab().get_coords(vert_conn, vert_num, &vpos[0]);
    double vol = Tools::tetVolume(vpos.data());
 
    error_indic_loc[IndicInds::ERROR_NORM] += err_indic_sum * vol;
  }
 
  CHKERR copyTagOnSkin(mField, "SOLUTION_ERROR", MB_TYPE_DOUBLE);
 
  error_indic_loc[IndicInds::COUNT_ENTS] =
      static_cast<double>(domain_ents_part.size());
 
  std::array<double, 3> error_indic_glob = {0.0, 0.0, 0.0};
 
  MPI_Allreduce(&error_indic_loc[0], &error_indic_glob[0], 3, MPI_DOUBLE,
                MPI_SUM, mField.get_comm());
 
  errorIndicMean = error_indic_glob[IndicInds::ERROR_TOT] /
                   error_indic_glob[IndicInds::COUNT_ENTS];
 
  MOFEM_LOG("WORLD", Sev::inform) << "Global error indicator (norm): "
                                  << error_indic_glob[IndicInds::ERROR_NORM];
  MOFEM_LOG("WORLD", Sev::inform)
      << "Global error indicator (mean): " << errorIndicMean;
 
  MoFEMFunctionReturn(0);
}

computeErrorNorms()

MoFEMErrorCode ElasticAdaptiveExample::computeErrorNorms ( )

private

[Postprocess results]

[Error Norms]

Definition at line 216 of file ElasticAdaptiveExample.hpp.

                                                         {
  MOFEM_LOG_CHANNEL("WORLD");
  auto simple = mField.getInterface<Simple>();
  MoFEMFunctionBegin;
 
  auto integration_rule = [](int, int, int p_data) { return p_data + 1; };
  auto evaluation_fe = boost::make_shared<DomainEle>(mField);
  evaluation_fe->getRuleHook = integration_rule;
 
  auto &evaluation_pipeline = evaluation_fe->getOpPtrVector();
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(evaluation_pipeline,
                                                        {H1}, "GEOMETRY");
 
  auto u_ptr = boost::make_shared<MatrixDouble>();
  evaluation_pipeline.push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
 
  auto analytical_disp_ptr = boost::make_shared<MatrixDouble>();
 
  evaluation_pipeline.push_back(new OpGetTensor1fromFunc<SPACE_DIM, SPACE_DIM>(
      analytical_disp_ptr, analyticalDisplacement));
 
  auto normsVec = createVectorMPI(mField.get_comm(),
                                  (mField.get_comm_rank() == 0) ? 1 : 0, 1);
 
  evaluation_pipeline.push_back(new OpCalcNormL2Tensor1<SPACE_DIM>(
      u_ptr, normsVec, 0, analytical_disp_ptr));
 
  CHKERR VecZeroEntries(normsVec);
  CHKERR DMoFEMLoopFiniteElements(simple->getDM(), simple->getDomainFEName(),
                                  evaluation_fe);
  CHKERR VecAssemblyBegin(normsVec);
  CHKERR VecAssemblyEnd(normsVec);
 
  if (mField.get_comm_rank() == 0) {
    const double *norms;
    CHKERR VecGetArrayRead(normsVec, &norms);
    MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
        << "Displacement error L2 norm: " << std::scientific
        << std::sqrt(norms[0]);
 
    if (atomTest == 3 && std::sqrt(norms[0]) > 2.195e-06) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "atom test failed: displacement error L2 norm is too high: "
              "%3.6e",
              std::sqrt(norms[0]));
    }
 
    if (atomTest == 4 && std::sqrt(norms[0]) > 1.995e-07) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "atom test failed: displacement error L2 norm is too high: "
              "%3.6e",
              std::sqrt(norms[0]));
    }
 
    CHKERR VecRestoreArrayRead(normsVec, &norms);
  }
 
  MoFEMFunctionReturn(0);
}

kspSetUpAndSolve()

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

overrideprivatevirtual

[Set up problem]

Reimplemented from ElasticExample.

Definition at line 155 of file ElasticAdaptiveExample.hpp.

                                                                  {
  MoFEMFunctionBegin;
 
  MOFEM_LOG_CHANNEL("TIMER");
  MOFEM_LOG_TAG("TIMER", "timer");
 
  DM dm;
  CHKERR KSPGetDM(solver, &dm);
  auto D = createDMVector(dm);
  auto F = vectorDuplicate(D);
 
  CHKERR VecZeroEntries(D);
  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
  CHKERR KSPSetUp(solver);
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
  CHKERR KSPSolve(solver, F, D);
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
 
  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
  MoFEMFunctionReturn(0);
}

outputResults()

MoFEMErrorCode ElasticAdaptiveExample::outputResults ( int  ref_iter = -1 )

private

[Postprocess results]

Definition at line 196 of file ElasticAdaptiveExample.hpp.

                                                                 {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  std::ostringstream strm;
  strm << "out_elastic";
  if (ref_iter >= 0) {
    strm << "_" << ref_iter;
  }
  strm << ".h5m";
 
  CHKERR ElasticPostProc::postProcessElasticResults<SPACE_DIM, DomainEle,
                                                    BoundaryEle, SideEle>(
      mField, simple->getDM(), simple->getDomainFEName(), strm.str(), {},
      {"GEOMETRY_ERROR", "SOLUTION_ERROR", "ORDER"}, Sev::verbose);
  MoFEMFunctionReturn(0);
}

readMesh()

MoFEMErrorCode ElasticAdaptiveExample::readMesh ( )

privatevirtual

[Run problem]

[Read mesh]

Reimplemented from ElasticExample.

Definition at line 79 of file ElasticAdaptiveExample.hpp.

                                                {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  CHKERR simple->getOptions();
 
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-ref_geom", &refGeom,
                             PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-ref_adapt_num", &refAdaptNum,
                            PETSC_NULLPTR);
 
  if ((refAdaptNum > 0) && refGeom) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Adaptive order and geometry refinement cannot be used together");
  }
 
  simple->getAddSkeletonFE() = true;
  CHKERR simple->loadFile();
 
  // Add meshsets if config file provided
  CHKERR mField.getInterface<MeshsetsManager>()->setMeshsetFromFile();
 
  auto update_ghost_ents = [&]() {
    MoFEMFunctionBegin;
 
    for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
             std::regex((boost::format("%s(.*)") % "MAT_ELASTIC").str()))
 
    ) {
      Range ents;
      CHKERR mField.get_moab().get_entities_by_handle(m->getMeshset(), ents,
                                                      true);
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ents);
      CHKERR mField.get_moab().add_entities(m->getMeshset(), ents);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR update_ghost_ents();
 
  MoFEMFunctionReturn(0);
}

refineGeometry()

MoFEMErrorCode ElasticAdaptiveExample::refineGeometry ( )

private

Definition at line 97 of file HighOrderGeometryRef.hpp.

                                                      {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-ref_geom_threshold",
                             &refGeomThreshold, PETSC_NULLPTR);
 
  if (Tag tag_err = 0; mField.get_moab().tag_get_handle(
                           "GEOMETRY_ERROR", tag_err) == MB_SUCCESS) {
 
    std::array<double, 2> mean_error_glob = {0.0, 0.0};
    MPI_Allreduce(&meanError[0], &mean_error_glob[0], 2, MPI_DOUBLE, MPI_SUM,
                  mField.get_comm());
    if (mean_error_glob[IndicInds::COUNT_ENTS] > 0) {
      refGeomThreshold *= (mean_error_glob[IndicInds::ERROR_TOT] /
                           mean_error_glob[IndicInds::COUNT_ENTS]);
    }
 
    MOFEM_LOG("WORLD", Sev::inform)
        << "Geometry refinement threshold = " << refGeomThreshold;
 
    auto get_ents_to_refine = [&]() {
      Range outer_ref_ents;
 
      auto *op_ptr_ref = new BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE);
      op_ptr_ref->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
                                      EntityType type,
                                      EntitiesFieldData::EntData &data) {
        MoFEMFunctionBegin;
        auto op_ptr_ref = static_cast<BoundaryEleOp *>(base_op_ptr);
        auto fe_ent = op_ptr_ref->getFEEntityHandle();
        double err_measure = 0;
        CHKERR mField.get_moab().tag_get_data(tag_err, &fe_ent, 1,
                                              &err_measure);
        if (err_measure > refGeomThreshold)
          outer_ref_ents.insert(fe_ent);
 
        MoFEMFunctionReturn(0);
      };
 
      auto fe = boost::make_shared<BoundaryEle>(mField);
      fe->getOpPtrVector().push_back(op_ptr_ref);
      auto dm = simple->getDM();
      CHK_THROW_MESSAGE(
          DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(), fe),
          "ref ents");
      return outer_ref_ents;
    };
 
    auto ents_to_ref = get_ents_to_refine();
 
    Range all_domain_ents;
    CHKERR mField.get_moab().get_entities_by_handle(simple->getMeshset(),
                                                    all_domain_ents);
    Range all_boundary_ents;
    CHKERR mField.get_moab().get_entities_by_handle(
        simple->getBoundaryMeshSet(), all_boundary_ents);
 
    Range tets_to_refine[2];
    CHKERR mField.get_moab().get_adjacencies(ents_to_ref, SPACE_DIM, false,
                                             tets_to_refine[1],
                                             moab::Interface::UNION);
    Range ent_conn;
    CHKERR mField.get_moab().get_adjacencies(tets_to_refine[1], SPACE_DIM - 1,
                                             true, ent_conn,
                                             moab::Interface::UNION);
    CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ent_conn);
 
    CHKERR mField.get_moab().get_adjacencies(
        ent_conn, SPACE_DIM, true, tets_to_refine[0], moab::Interface::UNION);
    tets_to_refine[0] = subtract(tets_to_refine[0], tets_to_refine[1]);
 
#ifndef NDEBUG
    CHKERR CutMeshInterface::SaveData(mField.get_moab())(
        "ents_to_refine_ho_2.vtk", tets_to_refine[0]);
    CHKERR CutMeshInterface::SaveData(mField.get_moab())(
        "ents_to_refine_ho_3.vtk", tets_to_refine[1]);
#endif
 
    for (int i = 0; i < 2; ++i) {
      Range ents_to_refine;
      ents_to_refine.merge(tets_to_refine[i]);
      for (int d = 1; d < SPACE_DIM; ++d) {
        Range lower_dim_ents;
        CHKERR mField.get_moab().get_adjacencies(
            tets_to_refine[i], d, true, lower_dim_ents, moab::Interface::UNION);
        ents_to_refine.merge(lower_dim_ents);
      }
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          ents_to_refine);
 
      CHKERR simple->setFieldOrder("U", oRder + i + 1, &ents_to_refine);
    }
 
    CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities("U");
    CHKERR simple->reSetUp();
 
    // save order on tags (volumes and boundary)
    int orders[] = {oRder, oRder + 1, oRder + 2};
    Tag tag_order;
 
    CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, tag_order);
 
    for (int i = 0; i < 2; ++i) {
      Range adj_faces;
      CHKERR mField.get_moab().get_adjacencies(tets_to_refine[i], SPACE_DIM - 1,
                                               true, adj_faces,
                                               moab::Interface::UNION);
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          adj_faces);
 
      adj_faces = intersect(adj_faces, all_boundary_ents);
      CHKERR mField.get_moab().tag_clear_data(tag_order, tets_to_refine[i],
                                              &orders[i + 1]);
      CHKERR mField.get_moab().tag_clear_data(tag_order, adj_faces,
                                              &orders[i + 1]);
    }
  }
  MoFEMFunctionReturn(0);
}

refineOrder()

MoFEMErrorCode ElasticAdaptiveExample::refineOrder ( int  ref_level )

private

Definition at line 147 of file AdaptiveOrderRef.hpp.

                                                                {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto comm = mField.getInterface<CommInterface>();
 
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-ref_adapt_threshold",
                             &refAdaptThreshold, PETSC_NULLPTR);
 
  Tag th_error_ind, th_order;
  CHKERR getTagHandle(mField, "SOLUTION_ERROR", MB_TYPE_DOUBLE, th_error_ind);
  CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, th_order);
 
  std::vector<Range> refinement_levels(ref_level + 2);
 
  Range domain_ents_part;
  CHKERR getDomainEntsPart(mField, domain_ents_part);
 
  for (auto ent : domain_ents_part) {
    double err_indic = 0;
    CHKERR mField.get_moab().tag_get_data(th_error_ind, &ent, 1, &err_indic);
    int order;
    CHKERR mField.get_moab().tag_get_data(th_order, &ent, 1, &order);
 
    if (err_indic > refAdaptThreshold * errorIndicMean) {
      order += 1;
      refinement_levels[order - oRder].insert(ent);
    }
  }
 
  for (int ll = refinement_levels.size() - 1; ll > 1; ll--) {
    CHKERR comm->synchroniseEntities(refinement_levels[ll]);
 
    Range face_ents;
    CHKERR mField.get_moab().get_adjacencies(refinement_levels[ll],
                                             SPACE_DIM - 1, true, face_ents,
                                             moab::Interface::UNION);
    Range tets;
    CHKERR mField.get_moab().get_adjacencies(face_ents, SPACE_DIM, true, tets,
                                             moab::Interface::UNION);
    tets = subtract(tets, refinement_levels[ll]);
    refinement_levels[ll - 1].merge(tets);
  }
 
  for (int ll = 1; ll < refinement_levels.size(); ll++) {
 
    if (oRder + ll > 8) {
      MOFEM_LOG("WORLD", Sev::warning)
          << "setting approximation order higher than 8 is not currently "
             "supported"
          << endl;
    } else {
 
      for (auto ent : refinement_levels[ll]) {
        int order = oRder + ll;
        CHKERR mField.get_moab().tag_set_data(th_order, &ent, 1, &order);
      }
 
      Range ents_to_refine;
      ents_to_refine.merge(refinement_levels[ll]);
      for (int d = 1; d < SPACE_DIM; ++d) {
        Range lower_dim_ents;
        CHKERR mField.get_moab().get_adjacencies(refinement_levels[ll], d, true,
                                                 lower_dim_ents,
                                                 moab::Interface::UNION);
        ents_to_refine.merge(lower_dim_ents);
      }
 
      CHKERR comm->synchroniseEntities(ents_to_refine);
 
      CHKERR simple->setFieldOrder("U", oRder + ll, &ents_to_refine);
    }
  }
 
  CHKERR copyTagOnSkin(mField, "ORDER", MB_TYPE_INTEGER);
 
  CHKERR comm->synchroniseFieldEntities("U");
  CHKERR simple->reSetUp();
  
  CHKERR boundaryCondition();
 
  MoFEMFunctionReturn(0);
}

runProblem()

MoFEMErrorCode ElasticAdaptiveExample::runProblem

(

)

[Run problem]

Definition at line 46 of file ElasticAdaptiveExample.hpp.

                                                  {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR setupAdaptivity();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  if (atomTest == 3)
    CHKERR computeErrorNorms();
  CHKERR checkResults();
 
  if (refAdaptNum > 0) {
    CHKERR computeErrorIndicators();
    CHKERR outputResults(0);
  } else {
    CHKERR outputResults();
  }
 
  for (int it = 0; it < refAdaptNum; ++it) {
    CHKERR refineOrder(it + 1);
    CHKERR solveSystem();
    CHKERR checkResults();
    if (it == refAdaptNum - 1 && atomTest == 4)
      CHKERR computeErrorNorms();
    CHKERR computeErrorIndicators();
    CHKERR outputResults(it + 1);
  }
  MoFEMFunctionReturn(0);
}

setupAdaptivity()

MoFEMErrorCode ElasticAdaptiveExample::setupAdaptivity ( )

private

[Read mesh]

[Set up problem]

Definition at line 125 of file ElasticAdaptiveExample.hpp.

                                                       {
  MoFEMFunctionBegin;
 
  Range domainEntities;
  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
                                                     domainEntities);
  Tag th_order;
  CHKERR getTagHandle(mField, "ORDER", MB_TYPE_INTEGER, th_order);
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &oRder, PETSC_NULLPTR);
  for (auto ent : domainEntities) {
    CHKERR mField.get_moab().tag_set_data(th_order, &ent, 1, &oRder);
  }
 
  CHKERR copyTagOnSkin(mField, "ORDER", MB_TYPE_INTEGER);
 
  if (refGeom) {
    CHKERR calculateGeometryError();
    CHKERR refineGeometry(); // optional ho refinement
  }
 
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-test", &atomTest,
                            PETSC_NULLPTR);
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

atomTest

int ElasticAdaptiveExample::atomTest = 0

private

Definition at line 37 of file ElasticAdaptiveExample.hpp.

errorIndicMean

double ElasticAdaptiveExample::errorIndicMean = 0.0

private

Definition at line 35 of file ElasticAdaptiveExample.hpp.

meanError

std::array<double, 2> ElasticAdaptiveExample::meanError = {0.0, 0.0}

private

Definition at line 34 of file ElasticAdaptiveExample.hpp.

{0.0, 0.0};

oRder

int ElasticAdaptiveExample::oRder = 2

private

Definition at line 29 of file ElasticAdaptiveExample.hpp.

refAdaptNum

int ElasticAdaptiveExample::refAdaptNum = 0

private

Definition at line 31 of file ElasticAdaptiveExample.hpp.

refAdaptThreshold

double ElasticAdaptiveExample::refAdaptThreshold = 1

private

Definition at line 32 of file ElasticAdaptiveExample.hpp.

refGeom

PetscBool ElasticAdaptiveExample::refGeom = PETSC_FALSE

private

Definition at line 30 of file ElasticAdaptiveExample.hpp.

refGeomThreshold

double ElasticAdaptiveExample::refGeomThreshold = 1

private

Definition at line 33 of file ElasticAdaptiveExample.hpp.

---

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

• tutorials/vec-8_elasticity_adaptive/src/ElasticAdaptiveExample.hpp
• tutorials/vec-8_elasticity_adaptive/src/AdaptiveOrderRef.hpp
• tutorials/vec-8_elasticity_adaptive/src/HighOrderGeometryRef.hpp