v0.14.0
adolc_plasticity.cpp

\brife Implementation of plasticity problem using ADOLC-C

/**
* \file adolc_plasticity.cpp
* \ingroup adoc_plasticity
* \example adolc_plasticity.cpp
*
* \brife Implementation of plasticity problem using ADOLC-C
*
*/
#include <MoFEM.hpp>
#include <MatrixFunction.hpp>
using namespace MoFEM;
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
using EntData = EntitiesFieldData::EntData;
/** Select finite element type for integration on domain based on problem
* dimension
*/
using DomainEle = PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::DomainEle;
/** Select finite element type for integrate on boundary based on problem
* dimension
*/
using BoundaryEle =
PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
/** Operators used to assemble domain integrals
*/
using DomainEleOp = DomainEle::UserDataOperator;
/** Operators used to assemble boundary integrals
*/
using BoundaryEleOp = BoundaryEle::UserDataOperator;
/** Linear forms used to integrate body forces
*/
using DomainNaturalBC =
NaturalBC<DomainEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
/** Select linear froms reading data from blockest (e.g. "BODY_FORCE") and
* applying body force.
*/
using OpBodyForce =
DomainNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, SPACE_DIM>;
// Select natural BC boundary condition
// Use natural boundary conditions implemented with adv-1 example which includes
// spring BC
#include <ContactNaturalBC.hpp>
/** Use Gauss quadrature rule and PETSc assembly to integrate neural boundary
* conditions. Select linear forms operators.
*/
using BoundaryRhsBCs =
NaturalBC<BoundaryEleOp>::Assembly<PETSC>::LinearForm<GAUSS>;
/* Use specialization from adv-1 integrating boundary conditions on forces and
* with springs. OP is used to integrate the right hand side
*/
using OpBoundaryRhsBCs =
BoundaryRhsBCs::OpFlux<ContactOps::BoundaryBCs, 1, SPACE_DIM>;
/** Use Gauss quadrature rule and PETSc assembly to integrate neural boundary
* conditions. Select bi-linear forms operators.
*/
using BoundaryLhsBCs =
NaturalBC<BoundaryEleOp>::Assembly<PETSC>::BiLinearForm<GAUSS>;
/** Use specialization from adv-1 integrating boundary conditions on forces and
* with springs
*/
using OpBoundaryLhsBCs =
BoundaryLhsBCs::OpFlux<ContactOps::BoundaryBCs, 1, SPACE_DIM>;
// Note: We create only skin post-processing mesh.
template <int DIM>
struct ElementsAndOps; //< Select finite elements and operators used for
// postprocessing and contact space
template <> struct ElementsAndOps<2> {
using PostProcEleDomain = PostProcBrokenMeshInMoabBase<DomainEle>;
using PostProcEleBdy = PostProcBrokenMeshInMoabBase<BoundaryEle>;
using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
static constexpr FieldSpace CONTACT_SPACE = HCURL;
};
template <> struct ElementsAndOps<3> {
using PostProcEleDomain = PostProcBrokenMeshInMoabBase<DomainEle>;
using PostProcEleBdy = PostProcBrokenMeshInMoabBase<BoundaryEle>;
using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
static constexpr FieldSpace CONTACT_SPACE = HDIV;
};
// Define contact space by dimension
constexpr FieldSpace ElementsAndOps<2>::CONTACT_SPACE;
constexpr FieldSpace ElementsAndOps<3>::CONTACT_SPACE;
constexpr FieldSpace CONTACT_SPACE = ElementsAndOps<SPACE_DIM>::CONTACT_SPACE;
using PostProcEleDomain = ElementsAndOps<SPACE_DIM>::PostProcEleDomain;
using SideEle = ElementsAndOps<SPACE_DIM>::SideEle;
using PostProcEleBdy = ElementsAndOps<SPACE_DIM>::PostProcEleBdy;
#include <ADOLCPlasticity.hpp>
#include <ADOLCPlasticityMaterialModels.hpp>
using namespace ADOLCPlasticity;
double scale = 1.;
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#endif
namespace ContactOps {
double cn_contact = 0.1;
}; // namespace ContactOps
#include <ContactOps.hpp>
#endif // ADD_CONTACT
struct PlasticProblem {
PlasticProblem(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEM::Interface &mField;
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode gettingNorms();
MoFEMErrorCode outputResults();
MoFEMErrorCode checkResults();
int approxOrder = 2;
int geom_order = 2;
FieldApproximationBase approxBase = AINSWORTH_LEGENDRE_BASE;
boost::shared_ptr<ClosestPointProjection> cpPtr; ///< Closest point
///< projection
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
boost::shared_ptr<ContactOps::SDFPython> sdfPythonPtr;
#endif
#endif // ADD_CONTACT
};
//! [Run problem]
MoFEMErrorCode PlasticProblem::runProblem() {
MoFEMFunctionBegin;
CHKERR readMesh();
CHKERR setupProblem();
CHKERR boundaryCondition();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR gettingNorms();
CHKERR outputResults();
CHKERR checkResults();
MoFEMFunctionReturn(0);
}
//! [Run problem]
//! [Read mesh]
MoFEMErrorCode PlasticProblem::readMesh() {
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);
}
//! [Read mesh]
//! [Set up problem]
MoFEMErrorCode PlasticProblem::setupProblem() {
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);
}
//! [Set up problem]
//! [Boundary condition]
MoFEMErrorCode PlasticProblem::boundaryCondition() {
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);
}
//! [Boundary condition]
//! [Push operators to pipeline]
MoFEMErrorCode PlasticProblem::assembleSystem() {
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);
}
//! [Push operators to pipeline]
/**
* @brief Monitor solution
*
* This functions is called by TS solver at the end of each step. It is used
* to output results to the hard drive.
*/
struct Monitor : public FEMethod {
Monitor(SmartPetscObj<DM> dm,
std::pair<boost::shared_ptr<PostProcEleDomain>,
boost::shared_ptr<PostProcEleBdy>>
pair_post_proc_fe,
boost::shared_ptr<DomainEle> reaction_fe,
std::vector<boost::shared_ptr<ScalingMethod>> smv)
: dM(dm), reactionFE(reaction_fe), vecOfTimeScalingMethods(smv) {
fRes = createDMVector(dM);
domainPostProcFe = pair_post_proc_fe.first;
skinPostProcFe = pair_post_proc_fe.second;
MoFEM::Interface *m_field_ptr;
CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
#ifdef ADD_CONTACT
auto get_integrate_traction = [&]() {
auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
CHK_THROW_MESSAGE(
(AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
integrate_traction->getOpPtrVector(), {HDIV}, "GEOMETRY")),
"Apply transform");
// We have to integrate on curved face geometry, thus integration weight
// have to adjusted.
integrate_traction->getOpPtrVector().push_back(
new OpSetHOWeightsOnSubDim<SPACE_DIM>());
integrate_traction->getRuleHook = [](int, int, int approx_order) {
return 2 * approx_order + 2 - 1;
};
CHK_THROW_MESSAGE(
(ContactOps::opFactoryCalculateTraction<SPACE_DIM, GAUSS,
BoundaryEleOp>(
integrate_traction->getOpPtrVector(), "SIGMA", 0)),
"push operators to calculate traction");
return integrate_traction;
};
integrateTraction = get_integrate_traction();
#endif
}
MoFEMErrorCode postProcess() {
MoFEMFunctionBegin;
int save_every_nth_step = 1;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every_nth_step",
&save_every_nth_step, PETSC_NULL);
#ifdef ADD_CONTACT
auto print_traction = [&](const std::string msg) {
MoFEMFunctionBegin;
MoFEM::Interface *m_field_ptr;
CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
if (!m_field_ptr->get_comm_rank()) {
const double *t_ptr;
CHKERR VecGetArrayRead(ContactOps::CommonData::totalTraction, &t_ptr);
MOFEM_LOG_C("CONTACT", Sev::inform, "%s time %3.4e %3.4e %3.4e %3.4e",
msg.c_str(), ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
CHKERR VecRestoreArrayRead(ContactOps::CommonData::totalTraction,
&t_ptr);
}
MoFEMFunctionReturn(0);
};
#endif
auto make_vtk = [&]() {
MoFEMFunctionBegin;
if (domainPostProcFe) {
CHKERR DMoFEMLoopFiniteElements(dM, "dFE", domainPostProcFe,
getCacheWeakPtr());
CHKERR domainPostProcFe->writeFile(
"out_plastic_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
}
if (skinPostProcFe) {
CHKERR DMoFEMLoopFiniteElements(dM, "bFE", skinPostProcFe,
getCacheWeakPtr());
CHKERR skinPostProcFe->writeFile(
"out_skin_plastic_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
}
MoFEMFunctionReturn(0);
};
if (!(ts_step % save_every_nth_step)) {
CHKERR make_vtk();
}
if (reactionFE) {
CHKERR VecZeroEntries(fRes);
reactionFE->f = fRes;
CHKERR DMoFEMLoopFiniteElements(dM, "dFE", reactionFE, getCacheWeakPtr());
CHKERR VecAssemblyBegin(fRes);
CHKERR VecAssemblyEnd(fRes);
CHKERR VecGhostUpdateBegin(fRes, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(fRes, ADD_VALUES, SCATTER_REVERSE);
MoFEM::Interface *m_field_ptr;
CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
*m_field_ptr, reactionFE, fRes)();
double nrm;
CHKERR VecNorm(fRes, NORM_2, &nrm);
MOFEM_LOG("PlasticPrb", Sev::verbose)
<< "Residual norm " << nrm << " at time step " << ts_step;
}
#ifdef ADD_CONTACT
auto calculate_traction = [&] {
MoFEMFunctionBegin;
CHKERR VecZeroEntries(ContactOps::CommonData::totalTraction);
CHKERR DMoFEMLoopFiniteElements(dM, "bFE", integrateTraction);
CHKERR VecAssemblyBegin(ContactOps::CommonData::totalTraction);
CHKERR VecAssemblyEnd(ContactOps::CommonData::totalTraction);
MoFEMFunctionReturn(0);
};
#endif
auto get_min_max_displacement = [&]() {
MoFEMFunctionBegin;
MoFEM::Interface *m_field_ptr;
CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
std::array<double, 4> a_min = {DBL_MAX, DBL_MAX, DBL_MAX, 0};
std::array<double, 4> a_max = {-DBL_MAX, -DBL_MAX, -DBL_MAX, 0};
auto get_min_max = [&](boost::shared_ptr<FieldEntity> field_entity_ptr) {
MoFEMFunctionBegin;
int d = 0;
for (auto v : field_entity_ptr->getEntFieldData()) {
a_min[d] = std::min(a_min[d], v);
a_max[d] = std::max(a_max[d], v);
++d;
}
MoFEMFunctionReturn(0);
};
a_min[SPACE_DIM] = 0;
a_max[SPACE_DIM] = 0;
Range verts;
CHKERR m_field_ptr->get_field_entities_by_type("U", MBVERTEX, verts);
CHKERR m_field_ptr->getInterface<FieldBlas>()->fieldLambdaOnEntities(
get_min_max, "U", &verts);
auto mpi_reduce = [&](auto &a, auto op) {
std::array<double, 3> a_mpi = {0, 0, 0};
MPI_Allreduce(a.data(), a_mpi.data(), 3, MPI_DOUBLE, op,
m_field_ptr->get_comm());
return a_mpi;
};
auto a_min_mpi = mpi_reduce(a_min, MPI_MIN);
auto a_max_mpi = mpi_reduce(a_max, MPI_MAX);
MOFEM_LOG("PlasticPrb", Sev::inform)
<< "Min displacement " << a_min_mpi[0] << " " << a_min_mpi[1] << " "
<< a_min_mpi[2];
MOFEM_LOG("PlasticPrb", Sev::inform)
<< "Max displacement " << a_max_mpi[0] << " " << a_max_mpi[1] << " "
<< a_max_mpi[2];
MoFEMFunctionReturn(0);
};
CHKERR get_min_max_displacement();
#ifdef ADD_CONTACT
CHKERR calculate_traction();
CHKERR print_traction("Contact force");
#endif
double scale = 1;
for (auto s : vecOfTimeScalingMethods) {
scale *= s->getScale(this->ts_t);
}
MOFEM_LOG("PlasticPrb", Sev::inform)
<< "Time: " << this->ts_t << " scale: " << scale;
MoFEMFunctionReturn(0);
}
private:
SmartPetscObj<DM> dM;
boost::shared_ptr<PostProcEleDomain> domainPostProcFe;
boost::shared_ptr<PostProcEleBdy> skinPostProcFe;
boost::shared_ptr<FEMethod> reactionFE;
boost::shared_ptr<BoundaryEle> integrateTraction;
SmartPetscObj<Vec> fRes;
VecOfTimeScalingMethods vecOfTimeScalingMethods;
};
static boost::shared_ptr<TSUpdate> ts_update_ptr = nullptr;
//! [Solve]
MoFEMErrorCode PlasticProblem::solveSystem() {
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);
}
//! [Solve]
//! [Getting norms]
MoFEMErrorCode PlasticProblem::gettingNorms() {
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);
}
//! [Getting norms]
//! [Postprocessing results]
MoFEMErrorCode PlasticProblem::outputResults() {
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);
}
//! [Postprocessing results]
//! [Check]
MoFEMErrorCode PlasticProblem::checkResults() {
MoFEMFunctionBegin;
#ifdef ADD_CONTACT
ContactOps::CommonData::totalTraction.reset();
#endif
MoFEMFunctionReturn(0);
}
//! [Check]
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
Py_Initialize();
np::initialize();
#endif
#endif // ADD_CONTACT
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "PlasticPrb"));
LogManager::setLog("PlasticPrb");
MOFEM_LOG_TAG("PlasticPrb", "PlasticPrb");
MOFEM_LOG("PlasticPrb", Sev::inform) << "SPACE_DIM " << SPACE_DIM;
#ifdef ADD_CONTACT
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
LogManager::setLog("CONTACT");
MOFEM_LOG_TAG("CONTACT", "Contact");
#endif // ADD_CONTACT
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
CHKERR DMRegister_MoFEM(dm_name);
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database interface
//! [Create MoFEM]
//! [PlasticProblem]
PlasticProblem ex(m_field);
CHKERR ex.runProblem();
//! [PlasticProblem]
}
CATCH_ERRORS;
CHKERR MoFEM::Core::Finalize();
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
#endif // ADD_CONTACT
return 0;
}
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
SPACE_DIM
constexpr int SPACE_DIM
Definition: adolc_plasticity.cpp:14
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
PlasticProblem::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: adolc_plasticity.cpp:693
MoFEM::EssentialPreProcReaction< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:157
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
scale
double scale
Definition: adolc_plasticity.cpp:109
LASTBASE
@ LASTBASE
Definition: definitions.h:69
MoFEM::OpLoopThis
Execute "this" element in the operator.
Definition: DGProjection.hpp:68
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
ADOLCPlasticity::createTSUpdate
boost::shared_ptr< TSUpdate > createTSUpdate(std::string fe_name, boost::shared_ptr< FEMethod > fe_ptr)
ContactOps
Definition: contact.cpp:99
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::EssentialPostProcLhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:134
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::NaturalBC::Assembly::BiLinearForm
Definition: Natural.hpp:74
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition: dynamic_first_order_con_law.cpp:53
ts_update_ptr
static boost::shared_ptr< TSUpdate > ts_update_ptr
Definition: adolc_plasticity.cpp:690
MoFEM::EssentialPostProcRhs< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:113
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
MoFEM::OpDGProjectionCoefficients
Definition: DGProjection.hpp:119
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
help
static char help[]
[Check]
Definition: adolc_plasticity.cpp:1187
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
geom_order
int geom_order
Order if fixed.
Definition: plastic.cpp:137
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpFluxLhsImpl
Definition: Natural.hpp:43
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
PlasticProblem
Definition: adolc_plasticity.cpp:127
CONTACT_SPACE
constexpr FieldSpace CONTACT_SPACE
Definition: adolc_plasticity.cpp:100
PlasticProblem::outputResults
MoFEMErrorCode outputResults()
[Getting norms]
Definition: adolc_plasticity.cpp:1148
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
PlasticProblem::gettingNorms
MoFEMErrorCode gettingNorms()
[Solve]
Definition: adolc_plasticity.cpp:1095
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::EssentialPreProc< DisplacementCubitBcData >
Specialization for DisplacementCubitBcData.
Definition: EssentialDisplacementCubitBcData.hpp:91
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
BoundaryEleOp
ContactOps::CommonData::totalTraction
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:29
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
PlasticProblem::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: adolc_plasticity.cpp:188
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2572
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MatrixFunction.hpp
MoFEM::AddFluxToLhsPipelineImpl
Definition: Natural.hpp:49
MoFEM::OpCalcNormL2Tensor1
Get norm of input MatrixDouble for Tensor1.
Definition: NormsOperators.hpp:44
TSADAPTMOFEM
#define TSADAPTMOFEM
Definition: TsCtx.hpp:10
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
ADOLCPlasticity
Definition: ADOLCPlasticity.hpp:24
PostProcEleBdy
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
Definition: dynamic_first_order_con_law.cpp:55
PlasticProblem::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: adolc_plasticity.cpp:173
BiLinearForm
MoFEM::OpSetHOWeightsOnSubDim
Definition: HODataOperators.hpp:145
MoFEM::PostProcBrokenMeshInMoabBase
Definition: PostProcBrokenMeshInMoabBase.hpp:94
ContactOps.hpp
PlasticProblem::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: adolc_plasticity.cpp:359
PlasticProblem::assembleSystem
MoFEMErrorCode assembleSystem()
[Boundary condition]
Definition: adolc_plasticity.cpp:442
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
ContactOps::CommonData::createTotalTraction
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:31
ContactOps::cn_contact
double cn_contact
Definition: contact.cpp:100
PlasticProblem::checkResults
MoFEMErrorCode checkResults()
[Postprocessing results]
Definition: adolc_plasticity.cpp:1178
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
ElementsAndOps
Definition: child_and_parent.cpp:18
DomainEle
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomainEle DomainEle
Definition: adolc_plasticity.cpp:22
Range
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
EntData
EntitiesFieldData::EntData EntData
Definition: adolc_plasticity.cpp:17
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
main
int main(int argc, char *argv[])
Definition: adolc_plasticity.cpp:1189
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
MoFEM::OpDGProjectionEvaluation
Definition: DGProjection.hpp:136
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
MoFEM::DMoFEMGetInterfacePtr
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:414
ADOLCPlasticityMaterialModels.hpp
Matetial models for plasticity.
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
ADOLCPlasticity.hpp
Monitor
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:783
approx_order
int approx_order
Definition: test_broken_space.cpp:50
MoFEM::VecOfTimeScalingMethods
std::vector< boost::shared_ptr< ScalingMethod > > VecOfTimeScalingMethods
Vector of time scaling methods.
Definition: Natural.hpp:20
MoFEM::DMMoFEMTSSetMonitor
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::MoFEMSNESMonitorFields
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:232
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
PlasticProblem::cpPtr
boost::shared_ptr< ClosestPointProjection > cpPtr
Definition: adolc_plasticity.cpp:148
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
ContactNaturalBC.hpp
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::getDMSnesCtx
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1127
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:578
MoFEM::SmartPetscObj< DM >
PlasticProblem::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: adolc_plasticity.cpp:158
MoFEM::CoreInterface::get_field_entities_by_type
virtual MoFEMErrorCode get_field_entities_by_type(const std::string name, EntityType type, Range &ents) const =0
get entities in the field by type
ContactOps::opFactoryCalculateTraction
MoFEMErrorCode opFactoryCalculateTraction(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1364
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1141
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1290
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
PlasticProblem::mField
MoFEM::Interface & mField
Definition: adolc_plasticity.cpp:134
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::TSAdaptCreateMoFEM
PetscErrorCode TSAdaptCreateMoFEM(TSAdapt adapt)
Craete MOFEM adapt.
Definition: TsCtx.cpp:797
MoFEM::OpDGProjectionMassMatrix
Definition: DGProjection.hpp:106
Generated by Doxygen 1.8.17 and hosted at University of Glasgow