Collaboration diagram for PlasticProblem:

[legend]

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

MoFEMErrorCode	runProblem ()
	[Run problem] More...

Private Member Functions
MoFEMErrorCode	readMesh ()
	[Run problem] More...

MoFEMErrorCode	setupProblem ()
	[Read mesh] More...

MoFEMErrorCode	boundaryCondition ()
	[Set up problem] More...

MoFEMErrorCode	assembleSystem ()
	[Boundary condition] More...

MoFEMErrorCode	solveSystem ()
	[Solve] More...

MoFEMErrorCode	gettingNorms ()
	[Solve] More...

MoFEMErrorCode	outputResults ()
	[Getting norms] More...

MoFEMErrorCode	checkResults ()
	[Postprocessing results] More...

Private Attributes
MoFEM::Interface &	mField

int	approxOrder = 2

int	geom_order = 2

FieldApproximationBase	approxBase = AINSWORTH_LEGENDRE_BASE

boost::shared_ptr< ClosestPointProjection >	cpPtr

Detailed Description

Examples: adolc_plasticity.cpp.

Definition at line 127 of file adolc_plasticity.cpp.

Constructor & Destructor Documentation

◆ PlasticProblem()

PlasticProblem::PlasticProblem ( MoFEM::Interface & m_field )

inline

Definition at line 129 of file adolc_plasticity.cpp.

129 : mField(m_field) {}

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode PlasticProblem::assembleSystem ( )

private

[Boundary condition]

[Push operators to pipeline]

Examples: adolc_plasticity.cpp.

Definition at line 442 of file adolc_plasticity.cpp.

                                               {
   MoFEMFunctionBegin;
   auto *simple = mField.getInterface<Simple>();
   auto *pipeline_mng = mField.getInterface<PipelineManager>();
  
   // assemble operator to the right hand side
   auto add_domain_ops_lhs = [&](auto &pip) {
     MoFEMFunctionBegin;
     // push forward finite element bases from reference to physical element
     CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
                                                           "GEOMETRY");
     // create local common data
     auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
     // calculate displacement gradients at integration points
     pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
         "U", common_data_ptr->getGardAtGaussPtsPtr()));
     // assemble tangent operator
     CHKERR opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
         mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr);
     MoFEMFunctionReturn(0);
   };
  
   auto add_domain_ops_rhs = [&](auto &pip) {
     MoFEMFunctionBegin;
     // push forward finite element bases from reference to physical element
     CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
                                                           "GEOMETRY");
     // create local common data
     auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
     // calculate displacement gradients at integration points
     pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
         "U", common_data_ptr->getGardAtGaussPtsPtr()));
     // assemble residual
     CHKERR opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
         mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr, Sev::inform);
 #ifdef ADD_CONTACT
     CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
         pip, "SIGMA", "U");
 #endif // ADD_CONTACT
     MoFEMFunctionReturn(0);
   };
  
   // Push operators to the left hand side pipeline. Indicate that domain (i.e.
   // volume/interior) element is used.
   CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
   // Push operators to the right hand side pipeline.
   CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
  
   MoFEMFunctionReturn(0);
 }

◆ boundaryCondition()

MoFEMErrorCode PlasticProblem::boundaryCondition ( )

private

[Set up problem]

[Boundary condition]

Add body froces

Examples: adolc_plasticity.cpp.

Definition at line 359 of file adolc_plasticity.cpp.

                                                  {
   MoFEMFunctionBegin;
   auto *pipeline_mng = mField.getInterface<PipelineManager>();
   auto simple = mField.getInterface<Simple>();
   auto bc_mng = mField.getInterface<BcManager>();
   auto time_scale = boost::make_shared<TimeScale>();
  
   auto rule = [](int, int, int p) { return 2 * p; };
   CHKERR pipeline_mng->setDomainRhsIntegrationRule(cpPtr->integrationRule);
   CHKERR pipeline_mng->setDomainLhsIntegrationRule(cpPtr->integrationRule);
   CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule);
   CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(rule);
  
   CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
       pipeline_mng->getOpBoundaryLhsPipeline(), {HDIV}, "GEOMETRY");
   pipeline_mng->getOpBoundaryLhsPipeline().push_back(
       new OpSetHOWeightsOnSubDim<SPACE_DIM>());
  
   CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
       pipeline_mng->getOpBoundaryRhsPipeline(), {HDIV}, "GEOMETRY");
   pipeline_mng->getOpBoundaryRhsPipeline().push_back(
       new OpSetHOWeightsOnSubDim<SPACE_DIM>());
  
   // Add Natural BCs to RHS
   CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
       pipeline_mng->getOpBoundaryRhsPipeline(), mField, "U", {time_scale},
       Sev::inform);
   // Add Natural BCs to LHS
   CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
       pipeline_mng->getOpBoundaryLhsPipeline(), mField, "U", Sev::inform);
  
 #ifdef ADD_CONTACT
   CHKERR
   ContactOps::opFactoryBoundaryLhs<SPACE_DIM, PETSC, GAUSS, BoundaryEleOp>(
       pipeline_mng->getOpBoundaryLhsPipeline(), "SIGMA", "U");
   CHKERR
   ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEle>(
       mField, pipeline_mng->getOpBoundaryLhsPipeline(),
       simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
       cpPtr->integrationRule);
  
   CHKERR
   ContactOps::opFactoryBoundaryRhs<SPACE_DIM, PETSC, GAUSS, BoundaryEleOp>(
       pipeline_mng->getOpBoundaryRhsPipeline(), "SIGMA", "U");
 #endif // ADD_CONTACT
  
   //! Add body froces
   CHKERR DomainNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
       pipeline_mng->getOpDomainRhsPipeline(), mField, "U", {time_scale},
       "BODY_FORCE", Sev::inform);
  
   // Essential BC
   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);
 #ifdef ADD_CONTACT
   for (auto b : {"FIX_X", "REMOVE_X"})
     CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                              "SIGMA", 0, 0, false, true);
   for (auto b : {"FIX_Y", "REMOVE_Y"})
     CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                              "SIGMA", 1, 1, false, true);
   for (auto b : {"FIX_Z", "REMOVE_Z"})
     CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                              "SIGMA", 2, 2, false, true);
   for (auto b : {"FIX_ALL", "REMOVE_ALL"})
     CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                              "SIGMA", 0, 3, false, true);
   CHKERR bc_mng->removeBlockDOFsOnEntities(
       simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
 #endif
  
   CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
       simple->getProblemName(), "U");
  
   MoFEMFunctionReturn(0);
 }

◆ checkResults()

MoFEMErrorCode PlasticProblem::checkResults ( )

private

[Postprocessing results]

[Check]

Examples: adolc_plasticity.cpp.

Definition at line 1178 of file adolc_plasticity.cpp.

                                             {
   MoFEMFunctionBegin;
 #ifdef ADD_CONTACT
   ContactOps::CommonData::totalTraction.reset();
 #endif
   MoFEMFunctionReturn(0);
 }

◆ gettingNorms()

MoFEMErrorCode PlasticProblem::gettingNorms ( )

private

[Solve]

[Getting norms]

Examples: adolc_plasticity.cpp.

Definition at line 1095 of file adolc_plasticity.cpp.

                                             {
   MoFEMFunctionBegin;
  
   auto simple = mField.getInterface<Simple>();
   auto dm = simple->getDM();
  
   auto T = createDMVector(simple->getDM());
   CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
                                  SCATTER_FORWARD);
   double nrm2;
   CHKERR VecNorm(T, NORM_2, &nrm2);
   MOFEM_LOG("PlasticPrb", Sev::inform) << "Solution norm " << nrm2;
  
   auto post_proc_norm_fe = boost::make_shared<DomainEle>(mField);
  
   post_proc_norm_fe->getRuleHook = cpPtr->integrationRule;
  
   CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
       post_proc_norm_fe->getOpPtrVector(), {H1});
  
   enum NORMS { U_NORM_L2 = 0, PIOLA_NORM, LAST_NORM };
   auto norms_vec =
       createVectorMPI(mField.get_comm(),
                       (mField.get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
   CHKERR VecZeroEntries(norms_vec);
  
   auto u_ptr = boost::make_shared<MatrixDouble>();
   post_proc_norm_fe->getOpPtrVector().push_back(
       new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
  
   post_proc_norm_fe->getOpPtrVector().push_back(
       new OpCalcNormL2Tensor1<SPACE_DIM>(u_ptr, norms_vec, U_NORM_L2));
  
   CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                   post_proc_norm_fe);
  
   CHKERR VecAssemblyBegin(norms_vec);
   CHKERR VecAssemblyEnd(norms_vec);
  
   MOFEM_LOG_CHANNEL("SELF"); // Clear channel from old tags
   if (mField.get_comm_rank() == 0) {
     const double *norms;
     CHKERR VecGetArrayRead(norms_vec, &norms);
     MOFEM_TAG_AND_LOG("SELF", Sev::inform, "example")
         << "norm_u: " << std::scientific << std::sqrt(norms[U_NORM_L2]);
     CHKERR VecRestoreArrayRead(norms_vec, &norms);
   }
  
   MoFEMFunctionReturn(0);
 }

◆ outputResults()

MoFEMErrorCode PlasticProblem::outputResults ( )

private

[Getting norms]

[Postprocessing results]

Examples: adolc_plasticity.cpp.

Definition at line 1148 of file adolc_plasticity.cpp.

                                              {
   MoFEMFunctionBegin;
   PetscInt test_nb = 0;
   PetscBool test_flg = PETSC_FALSE;
   CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-test", &test_nb, &test_flg);
  
   if (test_flg) {
     auto simple = mField.getInterface<Simple>();
     auto T = createDMVector(simple->getDM());
     CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
                                    SCATTER_FORWARD);
     double nrm2;
     CHKERR VecNorm(T, NORM_2, &nrm2);
     MOFEM_LOG("PlasticPrb", Sev::verbose) << "Regression norm " << nrm2;
     double regression_value = 0;
     switch (test_nb) {
     default:
       SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong test number.");
       break;
     }
     if (fabs(nrm2 - regression_value) > 1e-2)
       SETERRQ2(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                "Regression test field; wrong norm value. %6.4e != %6.4e", nrm2,
                regression_value);
   }
   MoFEMFunctionReturn(0);
 }

◆ readMesh()

MoFEMErrorCode PlasticProblem::readMesh ( )

private

[Run problem]

[Read mesh]

Examples: adolc_plasticity.cpp.

Definition at line 173 of file adolc_plasticity.cpp.

                                         {
   MoFEMFunctionBegin;
   auto simple = mField.getInterface<Simple>();
   CHKERR simple->getOptions();
   CHKERR simple->loadFile();
  
   if (simple->getDim() != SPACE_DIM)
     SETERRQ2(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
              "Wrong dimension of mesh %d != %d", simple->getDim(), SPACE_DIM);
  
   MoFEMFunctionReturn(0);
 }

◆ runProblem()

MoFEMErrorCode PlasticProblem::runProblem ( )

[Run problem]

Examples: adolc_plasticity.cpp.

Definition at line 158 of file adolc_plasticity.cpp.

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

◆ setupProblem()

MoFEMErrorCode PlasticProblem::setupProblem ( )

private

[Read mesh]

[Set up problem]

[Material models selection]

FIXME: Here only array of material models is initalized. Each model has unique set of the ADOL-C tags. Pointer is attached based on block name to which entity belongs. That will enable heterogeneity of the model, in addition of heterogeneity of the material properties.

[Material models selection]

Examples: adolc_plasticity.cpp.

Definition at line 188 of file adolc_plasticity.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_NULL, NULL, "-base", list_bases,
                               LASBASETOPT, &choice_base_value, PETSC_NULL);
   CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &approxOrder, PETSC_NULL);
   CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
                             PETSC_NULL);
  
   switch (choice_base_value) {
   case AINSWORTH:
     approxBase = AINSWORTH_LEGENDRE_BASE;
     MOFEM_LOG("WORLD", Sev::inform)
         << "Set AINSWORTH_LEGENDRE_BASE for displacements";
     break;
   case DEMKOWICZ:
     approxBase = DEMKOWICZ_JACOBI_BASE;
     MOFEM_LOG("WORLD", Sev::inform)
         << "Set DEMKOWICZ_JACOBI_BASE for displacements";
     break;
   default:
     approxBase = LASTBASE;
     break;
   }
  
   // Add field
   CHKERR simple->addDomainField("U", H1, approxBase, SPACE_DIM);
   CHKERR simple->addBoundaryField("U", H1, approxBase, SPACE_DIM);
   CHKERR simple->addDataField("GEOMETRY", H1, approxBase, SPACE_DIM);
  
 #ifdef ADD_CONTACT
   CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
                                 SPACE_DIM);
   CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
                                   SPACE_DIM);
   auto get_skin = [&]() {
     Range body_ents;
     CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
     Skinner skin(&mField.get_moab());
     Range skin_ents;
     CHKERR skin.find_skin(0, body_ents, false, skin_ents);
     return skin_ents;
   };
  
   auto filter_blocks = [&](auto skin) {
     bool is_contact_block = false;
     Range contact_range;
     for (auto m :
          mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
  
              (boost::format("%s(.*)") % "CONTACT").str()
  
                  ))
  
     ) {
       is_contact_block =
           true; ///< bloks interation is collectibe, so that is set irrespective
                 ///< if there are enerities in given rank or not in the block
       MOFEM_LOG("CONTACT", Sev::inform)
           << "Find contact block set:  " << m->getName();
       auto meshset = m->getMeshset();
       Range contact_meshset_range;
       CHKERR mField.get_moab().get_entities_by_dimension(
           meshset, SPACE_DIM - 1, contact_meshset_range, true);
  
       CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
           contact_meshset_range);
       contact_range.merge(contact_meshset_range);
     }
     if (is_contact_block) {
       MOFEM_LOG("SYNC", Sev::inform)
           << "Nb entities in contact surface: " << contact_range.size();
       MOFEM_LOG_SYNCHRONISE(mField.get_comm());
       skin = intersect(skin, contact_range);
     }
     return skin;
   };
  
   auto filter_true_skin = [&](auto skin) {
     Range boundary_ents;
     ParallelComm *pcomm =
         ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
     CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                  PSTATUS_NOT, -1, &boundary_ents);
     return boundary_ents;
   };
  
   auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
 #ifdef ADD_CONTACT
   CHKERR simple->setFieldOrder("SIGMA", 0);
   CHKERR simple->setFieldOrder("SIGMA", approxOrder - 1, &boundary_ents);
   CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn_contact",
                                &ContactOps::cn_contact, PETSC_NULL);
   MOFEM_LOG("CONTACT", Sev::inform) << "cn_contact " << ContactOps::cn_contact;
  
 #endif
  
 #ifdef PYTHON_SDF
   sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
   CHKERR sdfPythonPtr->sdfInit("sdf.py");
   ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
 #endif
 #endif
  
   CHKERR simple->setFieldOrder("U", approxOrder);
   CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
   CHKERR simple->setUp();
  
   auto project_ho_geometry = [&]() {
     Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
     return mField.loop_dofs("GEOMETRY", ent_method);
   };
   PetscBool project_geometry = PETSC_TRUE;
   CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-project_geometry",
                              &project_geometry, PETSC_NULL);
   if (project_geometry) {
     CHKERR project_ho_geometry();
   }
  
   //! [Material models selection]
  
   /**
    * FIXME: Here only array of material models is initalized. Each model has
    * unique set of the ADOL-C tags. Pointer is attached based on block name to
    * which entity belongs. That will enable heterogeneity of the model, in
    * addition of heterogeneity of the material properties.
    */
  
   enum MaterialModel {
     VonMisses,
     VonMissesPlaneStress,
     Paraboloidal,
     LastModel
   };
   const char *list_materials[LastModel] = {"VonMisses", "VonMissesPlaneStress",
                                            "Paraboloidal"};
   PetscInt choice_material = VonMisses;
   CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-material", list_materials,
                               LastModel, &choice_material, PETSC_NULL);
  
   switch (choice_material) {
   case VonMisses:
     cpPtr = createMaterial<J2Plasticity<3>>();
     break;
   case VonMissesPlaneStress:
     if (SPACE_DIM != 2)
       SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
               "VonMissesPlaneStrain is only for 2D case");
     cpPtr = createMaterial<J2Plasticity<2>>();
     break;
   case Paraboloidal:
     cpPtr = createMaterial<ParaboloidalPlasticity>();
     break;
   default:
     SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Wrong material model");
   }
  
   if (approxBase == DEMKOWICZ_JACOBI_BASE && SPACE_DIM == 2)
     cpPtr->integrationRule = [](int, int, int p) { return 2 * (p - 2); };
  
   //! [Material models selection]
  
   MoFEMFunctionReturn(0);
 }

◆ solveSystem()

MoFEMErrorCode PlasticProblem::solveSystem ( )

private

[Solve]

[DG projection]

[Set up post-step]

[TS post-step function]

[Set up post-step]

[Create TS]

Examples: adolc_plasticity.cpp.

Definition at line 693 of file adolc_plasticity.cpp.

                                            {
   MoFEMFunctionBegin;
   auto *simple = mField.getInterface<Simple>();
   auto *pipeline_mng = mField.getInterface<PipelineManager>();
  
   auto dm = simple->getDM();
   auto time_scale = boost::make_shared<TimeScale>();
   auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
  
   // Setup postprocessing
   auto create_post_proc_fe = [dm, this, simple]() {
     //! [DG projection]
     // Post-process domain, i.e. volume elements. If 3d case creates stresses
     // are evaluated at points which are trace of the volume element on boundary
     auto post_proc_ele_domain = [this](auto &pip_domain, auto &fe_name) {
       // Push bases on reference element to the phyiscal element
       CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
           pip_domain, {H1, HDIV}, "GEOMETRY");
  
       // calculate displacement gradients at nodes of post processing mesh. For
       // gradient DG projection is obsolete, since displacements can be
       // evaluated at arbitrary points.
       auto grad_ptr = boost::make_shared<MatrixDouble>();
       pip_domain.push_back(
           new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
                                                                    grad_ptr));
  
       // Create operator to run (neted) pipeline in operator. InteriorElem
       // element will have iteration rule and bases as domain element used to
       // solve problem, and remember history variables.
       using InteriorElem = DomainEle;
       auto op_this = new OpLoopThis<InteriorElem>(mField, fe_name, Sev::noisy);
       // add interior element to post-processing pipeline
       pip_domain.push_back(op_this);
  
       // get pointer to interior element, which is in essence pipeline which
       // pipeline.
       auto fe_physics = op_this->getThisFEPtr();
  
       auto evaluate_stress_on_physical_element = [&]() {
         // evaluate stress and plastic strain at Gauss points
         fe_physics->getRuleHook = cpPtr->integrationRule;
         CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
             fe_physics->getOpPtrVector(), {H1});
         auto common_data_ptr =
             boost::make_shared<ADOLCPlasticity::CommonData>();
         // note that gradient of displacements is evaluate again, at
         // physical nodes
         fe_physics->getOpPtrVector().push_back(
             new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
                 "U", common_data_ptr->getGardAtGaussPtsPtr()));
  
         CHKERR cpPtr->addMatBlockOps(mField, fe_physics->getOpPtrVector(),
                                      "ADOLCMAT", Sev::noisy);
         fe_physics->getOpPtrVector().push_back(
             getRawPtrOpCalculateStress<SPACE_DIM>(mField, common_data_ptr,
                                                   cpPtr, false));
         return common_data_ptr;
       };
  
       auto dg_projection_froward_and_back = [&](auto &common_data_ptr) {
         // here we do actual projection of stress and plastic strain to DG space
         int dg_order = approxOrder - 1;
         auto entity_data_l2 =
             boost::make_shared<EntitiesFieldData>(MBENTITYSET);
         auto mass_ptr =
             boost::make_shared<MatrixDouble>(); //< projection operator (mass
                                                 // matrix)
         fe_physics->getOpPtrVector().push_back(new OpDGProjectionMassMatrix(
             dg_order, mass_ptr, entity_data_l2, approxBase, L2));
  
         auto coeffs_ptr_stress = boost::make_shared<MatrixDouble>();
         // project stress on DG space on physical element
         fe_physics->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
             common_data_ptr->getStressMatrixPtr(), coeffs_ptr_stress, mass_ptr,
             entity_data_l2, approxBase, L2));
         // project strains plastic strains DG space on physical element
         auto coeffs_ptr_plastic_strain = boost::make_shared<MatrixDouble>();
         fe_physics->getOpPtrVector().push_back(new OpDGProjectionCoefficients(
             common_data_ptr->getPlasticStrainMatrixPtr(),
             coeffs_ptr_plastic_strain, mass_ptr, entity_data_l2, approxBase,
             L2));
  
         // project stress and plastic strain for DG space on post-process
         // element
         pip_domain.push_back(new OpDGProjectionEvaluation(
             common_data_ptr->getStressMatrixPtr(), coeffs_ptr_stress,
             entity_data_l2, approxBase, L2));
         pip_domain.push_back(new OpDGProjectionEvaluation(
             common_data_ptr->getPlasticStrainMatrixPtr(),
             coeffs_ptr_plastic_strain, entity_data_l2, approxBase, L2));
       };
  
       auto common_data_ptr = evaluate_stress_on_physical_element();
       dg_projection_froward_and_back(common_data_ptr);
  
       return boost::make_tuple(grad_ptr, common_data_ptr->getStressMatrixPtr(),
                                common_data_ptr->getPlasticStrainMatrixPtr());
     };
     //! [DG projection]
  
     // Create tags on post-processing mesh, i.e. those tags are visible in
     // Preview
     auto add_post_proc_map = [&](auto post_proc_fe, auto u_ptr, auto grad_ptr,
                                  auto stress_ptr, auto plastic_strain_ptr,
                                  auto contact_stress_ptr, auto X_ptr) {
       using OpPPMapSPACE_DIM = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
       post_proc_fe->getOpPtrVector().push_back(
  
           new OpPPMapSPACE_DIM(
  
               post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
  
               {},
  
               {{"U", u_ptr}, {"GEOMETRY", X_ptr}},
  
               {{"GRAD", grad_ptr}, {"SIGMA", contact_stress_ptr}},
  
               {}
  
               )
  
       );
  
       using OpPPMap3D = OpPostProcMapInMoab<3, 3>;
       post_proc_fe->getOpPtrVector().push_back(
  
           new OpPPMap3D(
  
               post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
  
               {},
  
               {},
  
               {},
  
               {{"STRESS", stress_ptr}, {"PLASTIC_STRAIN", plastic_strain_ptr}}
  
               )
  
       );
  
       return post_proc_fe;
     };
  
     auto vol_post_proc = [this, simple, post_proc_ele_domain,
                           add_post_proc_map]() {
       PetscBool post_proc_vol = PETSC_FALSE;
       if (SPACE_DIM == 2)
         post_proc_vol = PETSC_TRUE;
  
       CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",
                                  &post_proc_vol, PETSC_NULL);
       if (post_proc_vol == PETSC_FALSE)
         return boost::shared_ptr<PostProcEleDomain>();
       auto post_proc_fe = boost::make_shared<PostProcEleDomain>(mField);
       auto u_ptr = boost::make_shared<MatrixDouble>();
       post_proc_fe->getOpPtrVector().push_back(
           new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
       auto X_ptr = boost::make_shared<MatrixDouble>();
       post_proc_fe->getOpPtrVector().push_back(
           new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
       auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
 #ifdef ADD_CONTACT
       post_proc_fe->getOpPtrVector().push_back(
           new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
               "SIGMA", contact_stress_ptr));
 #else
       contact_stress_ptr = nullptr;
 #endif
       auto [grad_ptr, stress_ptr, plastic_strain_ptr] = post_proc_ele_domain(
           post_proc_fe->getOpPtrVector(), simple->getDomainFEName());
  
       return add_post_proc_map(post_proc_fe, u_ptr, grad_ptr, stress_ptr,
                                plastic_strain_ptr, contact_stress_ptr, X_ptr);
     };
  
     auto skin_post_proc = [this, simple, post_proc_ele_domain,
                            add_post_proc_map]() {
       // create skin of the volume mesh for post-processing,
       // i.e. boundary of the volume mesh
       PetscBool post_proc_skin = PETSC_TRUE;
       if (SPACE_DIM == 2)
         post_proc_skin = PETSC_FALSE;
  
       CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
                                  &post_proc_skin, PETSC_NULL);
  
       if (post_proc_skin == PETSC_FALSE)
         return boost::shared_ptr<PostProcEleBdy>();
  
       auto post_proc_fe = boost::make_shared<PostProcEleBdy>(mField);
       auto u_ptr = boost::make_shared<MatrixDouble>();
       post_proc_fe->getOpPtrVector().push_back(
           new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
       auto X_ptr = boost::make_shared<MatrixDouble>();
       post_proc_fe->getOpPtrVector().push_back(
           new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
       auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
       auto op_loop_side =
           new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
       auto [grad_ptr, stress_ptr, plastic_strain_ptr] = post_proc_ele_domain(
           op_loop_side->getOpPtrVector(), simple->getDomainFEName());
 #ifdef ADD_CONTACT
       op_loop_side->getOpPtrVector().push_back(
           new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
               "SIGMA", contact_stress_ptr));
 #else
       contact_stress_ptr = nullptr;
 #endif
       post_proc_fe->getOpPtrVector().push_back(op_loop_side);
  
       return add_post_proc_map(post_proc_fe, u_ptr, grad_ptr, stress_ptr,
                                plastic_strain_ptr, contact_stress_ptr, X_ptr);
     };
  
     return std::make_pair(vol_post_proc(), skin_post_proc());
   };
  
   auto create_reaction_fe = [&]() {
     auto fe_ptr = boost::make_shared<DomainEle>(mField);
     fe_ptr->getRuleHook = cpPtr->integrationRule;
  
     auto &pip = fe_ptr->getOpPtrVector();
     CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
     auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
     pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
         "U", common_data_ptr->getGardAtGaussPtsPtr()));
     CHKERR opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
         mField, "U", pip, "ADOLCMAT", common_data_ptr, cpPtr, Sev::noisy);
  
     return fe_ptr;
   };
  
   auto add_extra_finite_elements_to_ksp_solver_pipelines = [&]() {
     MoFEMFunctionBegin;
  
     auto pre_proc_ptr = boost::make_shared<FEMethod>();
     auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
     auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
  
     auto get_bc_hook_rhs = [this, pre_proc_ptr, time_scale]() {
       MoFEMFunctionBeginHot;
       CHKERR EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_ptr,
                                                        {time_scale}, false)();
       MoFEMFunctionReturnHot(0);
     };
  
     pre_proc_ptr->preProcessHook = get_bc_hook_rhs;
  
     auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
       MoFEMFunctionBegin;
       CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
           mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
       CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
           mField, post_proc_rhs_ptr, 1.)();
       MoFEMFunctionReturn(0);
     };
     auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
       MoFEMFunctionBegin;
       CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(
           mField, post_proc_lhs_ptr, 1.)();
       MoFEMFunctionReturn(0);
     };
     post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
     post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs;
  
     // This is low level pushing finite elements (pipelines) to solver
     auto ts_ctx_ptr = getDMTsCtx(simple->getDM());
     ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
     ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
     ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
     ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
     MoFEMFunctionReturn(0);
   };
  
   // Add extra finite elements to SNES solver pipelines to resolve essential
   // boundary conditions
   CHKERR add_extra_finite_elements_to_ksp_solver_pipelines();
  
   auto create_monitor_fe = [dm, time_scale](auto &&post_proc_fe,
                                             auto &&reaction_fe) {
     return boost::make_shared<Monitor>(
         dm, post_proc_fe, reaction_fe,
         std::vector<boost::shared_ptr<ScalingMethod>>{time_scale});
   };
  
   //! [Set up post-step]
   auto set_up_post_step = [&](auto ts) {
     MoFEMFunctionBegin;
  
     // create finite element (pipeline) and populate it with operators to
     // update history variables
     auto create_update_ptr = [&]() {
       auto fe_ptr = boost::make_shared<DomainEle>(mField);
       fe_ptr->getRuleHook = cpPtr->integrationRule;
       AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(fe_ptr->getOpPtrVector(),
                                                      {H1});
       auto common_data_ptr = boost::make_shared<ADOLCPlasticity::CommonData>();
       fe_ptr->getOpPtrVector().push_back(
           new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
               "U", common_data_ptr->getGardAtGaussPtsPtr()));
       CHK_THROW_MESSAGE(
           opFactoryDomainUpdate<SPACE_DIM>(mField, fe_ptr->getOpPtrVector(),
                                            "ADOLCMAT", common_data_ptr, cpPtr),
           "push update ops");
       return fe_ptr;
     };
  
     // ts_update_ptr is global static variable
     ts_update_ptr =
         createTSUpdate(simple->getDomainFEName(), create_update_ptr());
  
     //! [TS post-step function]
     // This is pure "C" function which we can to the TS solver
     auto ts_step_post_proc = [](TS ts) {
       MoFEMFunctionBegin;
       if (ts_update_ptr)
         CHKERR ts_update_ptr->postProcess(ts);
       MoFEMFunctionReturn(0);
     };
     //! [TS post-step function]
  
     // finally set up post-step
     CHKERR TSSetPostStep(ts, ts_step_post_proc);
     MoFEMFunctionReturn(0);
   };
   //! [Set up post-step]
  
   // Set monitor which postprocessing results and saves them to the hard drive
   auto set_up_monitor = [&](auto ts) {
     MoFEMFunctionBegin;
     boost::shared_ptr<FEMethod> null_fe;
     auto monitor_ptr =
         create_monitor_fe(create_post_proc_fe(), create_reaction_fe());
     CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
                                null_fe, monitor_ptr);
     MoFEMFunctionReturn(0);
   };
  
   auto set_section_monitor = [&](auto solver) {
     MoFEMFunctionBegin;
     SNES snes;
     CHKERR TSGetSNES(solver, &snes);
     CHKERR SNESMonitorSet(snes,
                           (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
                                              void *))MoFEMSNESMonitorFields,
                           (void *)(snes_ctx_ptr.get()), nullptr);
     MoFEMFunctionReturn(0);
   };
  
   auto set_up_adapt = [&](auto ts) {
     MoFEMFunctionBegin;
     CHKERR TSAdaptRegister(TSADAPTMOFEM, TSAdaptCreateMoFEM);
     TSAdapt adapt;
     CHKERR TSGetAdapt(ts, &adapt);
     MoFEMFunctionReturn(0);
   };
  
   //! [Create TS]
   auto ts = pipeline_mng->createTSIM();
  
   // Set time solver
   double ftime = 1;
   CHKERR TSSetMaxTime(ts, ftime);
   CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
   auto D = createDMVector(simple->getDM());
 #ifdef ADD_CONTACT
   ContactOps::CommonData::createTotalTraction(mField);
 #endif
   CHKERR TSSetSolution(ts, D);
   CHKERR set_section_monitor(ts);
  
   // Set monitor, step adaptivity, and post-step to update history variables
   CHKERR set_up_monitor(ts);
   CHKERR set_up_post_step(ts);
   CHKERR set_up_adapt(ts);
   CHKERR TSSetFromOptions(ts);
  
   CHKERR TSSolve(ts, NULL);
   //! [Create TS]
  
   CHKERR TSGetTime(ts, &ftime);
  
   PetscInt steps, snesfails, rejects, nonlinits, linits;
   CHKERR TSGetStepNumber(ts, &steps);
   CHKERR TSGetSNESFailures(ts, &snesfails);
   CHKERR TSGetStepRejections(ts, &rejects);
   CHKERR TSGetSNESIterations(ts, &nonlinits);
   CHKERR TSGetKSPIterations(ts, &linits);
   MOFEM_LOG_C("PlasticPrb", Sev::inform,
               "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
               "%d, linits %d",
               steps, rejects, snesfails, ftime, nonlinits, linits);
  
   MoFEMFunctionReturn(0);
 }

Member Data Documentation

◆ approxBase

FieldApproximationBase PlasticProblem::approxBase = AINSWORTH_LEGENDRE_BASE

private

Definition at line 147 of file adolc_plasticity.cpp.

◆ approxOrder

int PlasticProblem::approxOrder = 2

private

Definition at line 145 of file adolc_plasticity.cpp.

◆ cpPtr

boost::shared_ptr<ClosestPointProjection> PlasticProblem::cpPtr

private

Closest point projection

Examples: adolc_plasticity.cpp.

Definition at line 148 of file adolc_plasticity.cpp.

◆ geom_order

int PlasticProblem::geom_order = 2

private

Definition at line 146 of file adolc_plasticity.cpp.

◆ mField

MoFEM::Interface& PlasticProblem::mField

private

Examples: adolc_plasticity.cpp.

Definition at line 134 of file adolc_plasticity.cpp.

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

users_modules/adolc-plasticity/adolc_plasticity.cpp

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ PlasticProblem()

Member Function Documentation

◆ assembleSystem()

◆ boundaryCondition()

◆ checkResults()

◆ gettingNorms()

◆ outputResults()

◆ readMesh()

◆ runProblem()

◆ setupProblem()

◆ solveSystem()

Member Data Documentation

◆ approxBase

◆ approxOrder

◆ cpPtr

◆ geom_order

◆ mField