Plane stress elastic dynamic problem

* \file nonlinear_elastic.cpp
* \example nonlinear_elastic.cpp
* Plane stress elastic dynamic problem
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* License for more details.
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <MoFEM.hpp>
using namespace MoFEM;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
template <> struct ElementsAndOps<3> {
constexpr int SPACE_DIM = 3; //< Space dimension of problem, mesh
using OpK = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
using OpInternalForce = FormsIntegrators<DomainEleOp>::Assembly<
using OpBodyForce = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
GAUSS>::OpSource<1, SPACE_DIM>;
using OpBoundaryVec = FormsIntegrators<BoundaryEleOp>::Assembly<
constexpr double young_modulus = 100;
constexpr double poisson_ratio = 0.3;
constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
#include <HenckyOps.hpp>
using namespace HenckyOps;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode createCommonData();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
MoFEMErrorCode checkResults();
boost::shared_ptr<MatrixDouble> matGradPtr;
boost::shared_ptr<MatrixDouble> matStrainPtr;
boost::shared_ptr<MatrixDouble> matStressPtr;
boost::shared_ptr<MatrixDouble> matDPtr;
boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
//! [Create common data]
auto set_material_stiffness = [&]() {
constexpr double A =
(SPACE_DIM == 2) ? 2 * shear_modulus_G /
: 1;
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
t_kd(i, j) * t_kd(k, l);
matGradPtr = boost::make_shared<MatrixDouble>();
matStrainPtr = boost::make_shared<MatrixDouble>();
matStressPtr = boost::make_shared<MatrixDouble>();
matDPtr = boost::make_shared<MatrixDouble>();
commonHenckyDataPtr = boost::make_shared<HenckyOps::CommonData>();
commonHenckyDataPtr->matGradPtr = matGradPtr;
commonHenckyDataPtr->matDPtr = matDPtr;
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
matDPtr->resize(size_symm * size_symm, 1);
CHKERR set_material_stiffness();
//! [Create common data]
//! [Run problem]
CHKERR readMesh();
CHKERR setupProblem();
CHKERR createCommonData();
CHKERR boundaryCondition();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR outputResults();
CHKERR checkResults();
//! [Run problem]
//! [Read mesh]
auto simple = mField.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
//! [Read mesh]
//! [Set up problem]
Simple *simple = mField.getInterface<Simple>();
// Add field
int order = 2;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("U", order);
CHKERR simple->setUp();
//! [Set up problem]
//! [Boundary condition]
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
"U", 0, 3);
auto get_time = [&](double, double, double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_rhs = pipeline_mng->getDomainRhsFE();
return fe_domain_rhs->ts_t;
if (it->getName().compare(0, 5, "FORCE") == 0) {
Range force_edges;
std::vector<double> attr_vec;
CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
force_edges, true);
if (attr_vec.size() < SPACE_DIM)
"Wrong size of boundary attributes vector. Set right block "
"size attributes.");
auto force_vec_ptr = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
std::copy(&attr_vec[0], &attr_vec[SPACE_DIM],
new OpBoundaryVec("U", force_vec_ptr, get_time,
//! [Define gravity vector]
auto get_body_force = [this](const double, const double, const double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
const auto time = fe_domain_rhs->ts_t;
// hardcoded gravity load in y direction
t_source(i) = 0;
t_source(1) = 1 * time;
return t_source;
//! [Define gravity vector]
new OpBodyForce("U", get_body_force));
//! [Boundary condition]
//! [Push operators to pipeline]
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto add_domain_base_ops = [&](auto &pipeline) {
if (SPACE_DIM == 2) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
"U", matGradPtr));
new OpCalculateEigenVals<SPACE_DIM>("U", commonHenckyDataPtr));
new OpCalculateLogC<SPACE_DIM>("U", commonHenckyDataPtr));
new OpCalculateLogC_dC<SPACE_DIM>("U", commonHenckyDataPtr));
new OpCalculateHenckyStress<SPACE_DIM>("U", commonHenckyDataPtr));
new OpCalculatePiolaStress<SPACE_DIM>("U", commonHenckyDataPtr));
auto add_domain_ops_lhs = [&](auto &pipeline) {
new OpHenckyTangent<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(new OpK("U", "U", commonHenckyDataPtr->getMatTangent()));
auto add_domain_ops_rhs = [&](auto &pipeline) {
// Calculate internal forece
pipeline.push_back(new OpInternalForce(
"U", commonHenckyDataPtr->getMatFirstPiolaStress()));
CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
CHKERR add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
auto integration_rule = [](int, int, int approx_order) {
return 2 * (approx_order - 1);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
//! [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, boost::shared_ptr<PostProcEle> post_proc)
: dM(dm), postProc(post_proc){};
MoFEMErrorCode postProcess() {
constexpr int save_every_nth_step = 1;
if (ts_step % save_every_nth_step == 0) {
CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
CHKERR postProc->writeFile(
"out_step_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
SmartPetscObj<DM> dM;
boost::shared_ptr<PostProcEle> postProc;
//! [Solve]
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto dm = simple->getDM();
auto ts = pipeline_mng->createTSIM();
// Setup postprocessing
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
if (SPACE_DIM == 2) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
new OpCalculateHOJacForFace(jac_ptr));
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
new OpSetInvJacH1ForFace(inv_jac_ptr));
"U", commonHenckyDataPtr->matGradPtr));
new OpCalculateEigenVals<SPACE_DIM>("U", commonHenckyDataPtr));
new OpCalculateLogC<SPACE_DIM>("U", commonHenckyDataPtr));
new OpCalculateLogC_dC<SPACE_DIM>("U", commonHenckyDataPtr));
new OpCalculateHenckyStress<SPACE_DIM>("U", commonHenckyDataPtr));
new OpCalculatePiolaStress<SPACE_DIM>("U", commonHenckyDataPtr));
post_proc_fe->getOpPtrVector().push_back(new OpPostProcHencky<SPACE_DIM>(
"U", post_proc_fe->postProcMesh, post_proc_fe->mapGaussPts,
// Add monitor to time solver
boost::shared_ptr<FEMethod> null_fe;
auto monitor_ptr = boost::make_shared<Monitor>(dm, post_proc_fe);
CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
null_fe, monitor_ptr);
double ftime = 1;
CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
auto D = smartCreateDMVector(simple->getDM());
CHKERR TSSetSolution(ts, D);
CHKERR TSSetFromOptions(ts);
CHKERR TSGetTime(ts, &ftime);
PetscInt steps, snesfails, rejects, nonlinits, linits;
CHKERR TSGetTimeStepNumber(ts, &steps);
CHKERR TSGetSNESFailures(ts, &snesfails);
CHKERR TSGetStepRejections(ts, &rejects);
CHKERR TSGetSNESIterations(ts, &nonlinits);
CHKERR TSGetKSPIterations(ts, &linits);
MOFEM_LOG_C("EXAMPLE", Sev::inform,
"steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
"%D, linits %D\n",
steps, rejects, snesfails, ftime, nonlinits, linits);
//! [Solve]
//! [Postprocess results]
PetscBool test_flg = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
if (test_flg) {
auto *simple = mField.getInterface<Simple>();
auto T = smartCreateDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
double nrm2;
CHKERR VecNorm(T, NORM_2, &nrm2);
MOFEM_LOG("EXAMPLE", Sev::inform) << "Regression norm " << nrm2;
constexpr double regression_value = 1.09572;
if (fabs(nrm2 - regression_value) > 1e-2)
"Regression test faileed; wrong norm value.");
//! [Postprocess results]
//! [Check]
//! [Check]
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
// Add logging channel for example
auto core_log = logging::core::get();
MOFEM_LOG_TAG("EXAMPLE", "example");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [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 insterface
//! [Create MoFEM]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main(int argc, char *argv[])
static char help[]
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
Catch errors.
Definition: definitions.h:385
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
Definition: definitions.h:161
Definition: definitions.h:53
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto integration_rule
constexpr auto t_kd
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:544
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const double T
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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: DMMMoFEM.cpp:994
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
double A
int save_every_nth_step
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:95
static constexpr int approx_order
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode outputResults()
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
Postprocess on face.
Post processing.
EntitiesFieldData::EntData EntData
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpK
constexpr int SPACE_DIM
constexpr double poisson_ratio
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
constexpr double young_modulus