#include "tutorials/vec-0_elasticity/src/ElasticExample.hpp"

Inheritance diagram for ElasticExample:

Collaboration diagram for ElasticExample:

Public Member Functions
	ElasticExample (MoFEM::Interface &m_field)

MoFEMErrorCode	runProblem ()
	[Run problem]

Protected Member Functions
MoFEMErrorCode	readMesh ()
	[Run problem]

MoFEMErrorCode	setupProblem ()
	[Read mesh]

MoFEMErrorCode	boundaryCondition ()
	[Set up problem]

MoFEMErrorCode	assembleSystem ()
	[Boundary condition]

MoFEMErrorCode	solveSystem ()
	[Solve]

MoFEMErrorCode	outputResults ()
	[Solve]

MoFEMErrorCode	checkResults ()
	[Postprocess results]

virtual MoFEMErrorCode	kspSetUpAndSolve (SmartPetscObj< KSP > solver)
	[Push operators to pipeline]

Protected Attributes
MoFEM::Interface &	mField

boost::shared_ptr< MatrixDouble >	vectorFieldPtr = nullptr

Detailed Description

Definition at line 18 of file ElasticExample.hpp.

Constructor & Destructor Documentation

◆ ElasticExample()

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

inline

Definition at line 20 of file ElasticExample.hpp.

20: mField(m_field) {}

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode ElasticExample::assembleSystem ( )

protected

[Boundary condition]

[Push operators to pipeline]

[Integration rule]

[Push domain stiffness matrix to pipeline]

[Push Body forces]

[Push natural boundary conditions]

Examples: mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 138 of file ElasticExample.hpp.

                                              {
  MoFEMFunctionBegin;
  auto pip = mField.getInterface<PipelineManager>();
 
  //! [Integration rule]
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
 
  auto integration_rule_bc = [](int, int, int approx_order) {
    return 2 * approx_order + 1;
  };
 
  CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
  CHKERR pip->setDomainLhsIntegrationRule(integration_rule);
  CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
  CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
  //! [Integration rule]
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
 
  //! [Push domain stiffness matrix to pipeline]
  // Add LHS operator for elasticity (stiffness matrix)
  CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
      mField, pip->getOpDomainLhsPipeline(), "U", "MAT_ELASTIC", Sev::verbose);
  //! [Push domain stiffness matrix to pipeline]
 
  // Add RHS operator for internal forces
  CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
      mField, pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC", Sev::verbose);
 
  //! [Push Body forces]
  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
      pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
  //! [Push Body forces]
 
  //! [Push natural boundary conditions]
  // Add force boundary condition
  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
      pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
  // Add case for mix type of BCs
  CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
      pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::verbose);
  //! [Push natural boundary conditions]
  MoFEMFunctionReturn(0);
}

◆ boundaryCondition()

MoFEMErrorCode ElasticExample::boundaryCondition ( )

protected

[Set up problem]

[Boundary condition]

Examples: mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 117 of file ElasticExample.hpp.

                                                 {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
 
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
                                           "U", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
                                           "U", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
                                           "U", 2, 2);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                           "REMOVE_ALL", "U", 0, 3);
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "U");
 
  MoFEMFunctionReturn(0);
}

◆ checkResults()

MoFEMErrorCode ElasticExample::checkResults ( )

protected

[Postprocess results]

[Check]

Examples: mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 343 of file ElasticExample.hpp.

                                            {
  MOFEM_LOG_CHANNEL("WORLD");
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>();
  MoFEMFunctionBegin;
  pip->getDomainRhsFE().reset();
  pip->getDomainLhsFE().reset();
  pip->getBoundaryRhsFE().reset();
  pip->getBoundaryLhsFE().reset();
 
  auto integration_rule = [](int, int, int p_data) { return 2 * p_data + 1; };
  CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
  CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule);
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pip->getOpBoundaryRhsPipeline(), {}, "GEOMETRY");
 
  // Add RHS operators for internal forces
  CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
      mField, pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC", Sev::verbose);
 
  CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
      pip->getOpDomainRhsPipeline(), mField, "U", Sev::verbose);
 
  CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
      pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::verbose);
 
  auto dm = simple->getDM();
  auto res = createDMVector(dm);
  CHKERR VecSetDM(res, PETSC_NULLPTR);
 
  pip->getDomainRhsFE()->f = res;
  pip->getBoundaryRhsFE()->f = res;
 
  CHKERR VecZeroEntries(res);
 
  CHKERR pip->loopFiniteElements();
  // CHKERR mField.getInterface<FieldBlas>()->fieldScale(-1, "U");
 
  CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
  CHKERR VecAssemblyBegin(res);
  CHKERR VecAssemblyEnd(res);
 
  auto zero_residual_at_constrains = [&]() {
    MoFEMFunctionBegin;
    auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
    auto get_post_proc_hook_rhs = [&]() {
      MoFEMFunctionBegin;
      CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
          mField, fe_post_proc_ptr, res)();
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
          mField, fe_post_proc_ptr, 0, res)();
      MoFEMFunctionReturn(0);
    };
    fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
    CHKERR DMoFEMPostProcessFiniteElements(dm, fe_post_proc_ptr.get());
    MoFEMFunctionReturn(0);
  };
 
  CHKERR zero_residual_at_constrains();
 
  double nrm2;
  CHKERR VecNorm(res, NORM_2, &nrm2);
  MOFEM_LOG_CHANNEL("WORLD");
  MOFEM_LOG_C("WORLD", Sev::inform, "residual = %3.4e\n", nrm2);
 
  int test = 0;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-test", &test, PETSC_NULLPTR);
  if (test > 0) {
    auto post_proc_residual = [&](auto dm, auto f_res, auto out_name) {
      MoFEMFunctionBegin;
      auto post_proc_fe =
          boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
      auto u_vec = boost::make_shared<MatrixDouble>();
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_vec, f_res));
      post_proc_fe->getOpPtrVector().push_back(
 
          new OpPPMap(
 
              post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
              {},
 
              {{"RES", u_vec}},
 
              {}, {})
 
      );
 
      CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                      post_proc_fe);
      post_proc_fe->writeFile(out_name);
      MoFEMFunctionReturn(0);
    };
 
    CHKERR post_proc_residual(simple->getDM(), res, "res.h5m");
 
    constexpr double eps = 1e-8;
    if (nrm2 > eps)
      SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
              "Residual is not zero");
  }
  if (test == 2) {
    if (!vectorFieldPtr || vectorFieldPtr->size1() == 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "atom test %d failed: Field Evaluator did not provide result",
              test);
    }
    auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
    double Ux_ref = 0.46;
    double Uy_ref = -0.03;
    constexpr double eps = 1e-8;
    if (fabs(t_disp(0) - Ux_ref) > eps || fabs(t_disp(1) - Uy_ref) > eps) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "atom test %d failed: Ux_ref = %3.6e, computed = %3.6e, Uy_ref "
              "= %3.6e, computed = %3.6e",
              test, Ux_ref, t_disp(0), Uy_ref, t_disp(1));
    }
  }
  MoFEMFunctionReturn(0);
}

◆ kspSetUpAndSolve()

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

protectedvirtual

[Push operators to pipeline]

Reimplemented in ElasticSchurExample.

Definition at line 194 of file ElasticExample.hpp.

                                                                         {
  MoFEMFunctionBegin;
 
  MOFEM_LOG_CHANNEL("TIMER");
  MOFEM_LOG_TAG("TIMER", "timer");
 
  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
  CHKERR KSPSetUp(solver);
  MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
 
  DM dm;
  CHKERR KSPGetDM(solver, &dm);
  auto D = createDMVector(dm);
  auto F = vectorDuplicate(D);
 
  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 ElasticExample::outputResults ( )

protected

[Solve]

[Postprocess results]

Examples: mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 332 of file ElasticExample.hpp.

                                             {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  CHKERR ElasticPostProc::postProcessElasticResults<SPACE_DIM, DomainEle,
                                                    BoundaryEle, SideEle>(
      mField, simple->getDM(), simple->getDomainFEName(), "out_elastic.h5m");
  MoFEMFunctionReturn(0);
}

◆ readMesh()

MoFEMErrorCode ElasticExample::readMesh ( )

protected

[Run problem]

[Read mesh]

Examples: mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 55 of file ElasticExample.hpp.

                                        {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
  // Add meshsets if config file provided
  CHKERR mField.getInterface<MeshsetsManager>()->setMeshsetFromFile();
  MoFEMFunctionReturn(0);
}

◆ runProblem()

MoFEMErrorCode ElasticExample::runProblem ( )

[Run problem]

Definition at line 41 of file ElasticExample.hpp.

                                          {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR outputResults();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}

◆ setupProblem()

MoFEMErrorCode ElasticExample::setupProblem ( )

protected

[Read mesh]

[Set up problem]

Examples: mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 67 of file ElasticExample.hpp.

                                            {
  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;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set AINSWORTH_LEGENDRE_BASE for displacements";
    break;
  case DEMKOWICZ:
    base = DEMKOWICZ_JACOBI_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set DEMKOWICZ_JACOBI_BASE for displacements";
    break;
  default:
    base = LASTBASE;
    break;
  }
 
  // Add field
  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
  int order = 3;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
 
  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
 
  CHKERR simple->setFieldOrder("U", order);
  CHKERR simple->setFieldOrder("GEOMETRY", 2);
  CHKERR simple->setUp();
 
  auto project_ho_geometry = [&]() {
    Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
    return mField.loop_dofs("GEOMETRY", ent_method);
  };
  CHKERR project_ho_geometry();
 
  MoFEMFunctionReturn(0);
}

◆ solveSystem()

MoFEMErrorCode ElasticExample::solveSystem ( )

protected

[Solve]

Examples: mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 223 of file ElasticExample.hpp.

                                           {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>();
  auto solver = pip->createKSP();
  CHKERR KSPSetFromOptions(solver);
 
  auto set_essential_bc = [this]() {
    MoFEMFunctionBegin;
    // This is low level pushing finite elements (pipelines) to solver
    auto simple = mField.getInterface<Simple>();
    auto dm = simple->getDM();
    auto ksp_ctx_ptr = getDMKspCtx(dm);
 
    auto pre_proc_rhs = boost::make_shared<FEMethod>();
    auto post_proc_rhs = boost::make_shared<FEMethod>();
    auto post_proc_lhs = boost::make_shared<FEMethod>();
 
    auto get_pre_proc_hook = [this, pre_proc_rhs]() {
      return EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_rhs,
                                                       {});
    };
    pre_proc_rhs->preProcessHook = get_pre_proc_hook();
 
    auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
                                                           post_proc_rhs, 1.)();
      MoFEMFunctionReturnHot(0);
    };
 
    auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(mField,
                                                           post_proc_lhs, 1.)();
      MoFEMFunctionReturnHot(0);
    };
 
    post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
    post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
 
    ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
    ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
    ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
    MoFEMFunctionReturn(0);
  };
 
  auto evaluate_field_at_the_point = [&]() {
    MoFEMFunctionBegin;
 
    int coords_dim = 3;
    std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
    PetscBool do_eval_field = PETSC_FALSE;
    CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
                                    field_eval_coords.data(), &coords_dim,
                                    &do_eval_field);
 
    if (do_eval_field) {
 
      vectorFieldPtr = boost::make_shared<MatrixDouble>();
      auto field_eval_data =
          mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
 
      CHKERR mField.getInterface<FieldEvaluatorInterface>()
          ->buildTree<SPACE_DIM>(field_eval_data, simple->getDomainFEName());
 
      field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
      auto no_rule = [](int, int, int) { return -1; };
      auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
      field_eval_fe_ptr->getRuleHook = no_rule;
 
      field_eval_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", vectorFieldPtr));
 
      CHKERR mField.getInterface<FieldEvaluatorInterface>()
          ->evalFEAtThePoint<SPACE_DIM>(
              field_eval_coords.data(), 1e-12, simple->getProblemName(),
              simple->getDomainFEName(), field_eval_data,
              mField.get_comm_rank(), mField.get_comm_rank(), nullptr, MF_EXIST,
              QUIET);
 
      if (vectorFieldPtr->size1()) {
        auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
        if constexpr (SPACE_DIM == 2)
          MOFEM_LOG("FieldEvaluator", Sev::inform)
              << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
        else
          MOFEM_LOG("FieldEvaluator", Sev::inform)
              << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
              << " U_Z: " << t_disp(2);
      }
 
      MOFEM_LOG_SYNCHRONISE(mField.get_comm());
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_essential_bc();
  CHKERR kspSetUpAndSolve(solver);
  CHKERR evaluate_field_at_the_point();
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

◆ mField

MoFEM::Interface& ElasticExample::mField

protected

Examples: mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 25 of file ElasticExample.hpp.

◆ vectorFieldPtr

boost::shared_ptr<MatrixDouble> ElasticExample::vectorFieldPtr = nullptr

protected

Definition at line 26 of file ElasticExample.hpp.

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

tutorials/vec-0_elasticity/src/ElasticExample.hpp

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ ElasticExample()

Member Function Documentation

◆ assembleSystem()

◆ boundaryCondition()

◆ checkResults()

◆ kspSetUpAndSolve()

◆ outputResults()

◆ readMesh()

◆ runProblem()

◆ setupProblem()

◆ solveSystem()

Member Data Documentation

◆ mField

◆ vectorFieldPtr