v0.14.0
VEC-3: Nonlinear dynamic elastic
Note
Prerequisites of this tutorial include MSH-1: Create a 2D mesh from Gmsh


Note
Intended learning outcome:
  • general structure of a program developed using MoFEM
  • idea of Simple Interface in MoFEM and how to use it
  • second order equation in time
  • idea of domain element in MoFEM and how to use it
  • Use default forms integrals
  • how to push the developed UDOs to the Pipeline
  • utilisation of tools to convert outputs (MOAB) and visualise them (Paraview)
/**
* \file dynamic_elastic.cpp
* \example dynamic_elastic.cpp
*
* Plane stress elastic dynamic problem
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
};
template <> struct ElementsAndOps<3> {
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
PETSC>::LinearForm<GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 1>;
constexpr double rho = 1;
constexpr double omega = 2.4;
constexpr double young_modulus = 1;
constexpr double poisson_ratio = 0.25;
#include <HenckyOps.hpp>
using namespace HenckyOps;
static double *ts_time_ptr;
static double *ts_aa_ptr;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode boundaryCondition();
MoFEMErrorCode assembleSystem();
MoFEMErrorCode solveSystem();
MoFEMErrorCode outputResults();
MoFEMErrorCode checkResults();
};
//! [Run problem]
CHKERR readMesh();
CHKERR setupProblem();
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]
// Add field
CHKERR simple->addDataField("GEOMETRY", H1, AINSWORTH_LEGENDRE_BASE,
int order = 2;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("U", order);
CHKERR simple->setFieldOrder("GEOMETRY", order);
CHKERR simple->setUp();
auto project_ho_geometry = [&]() {
Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
return mField.loop_dofs("GEOMETRY", ent_method);
};
CHKERR project_ho_geometry();
}
//! [Set up problem]
//! [Boundary condition]
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
}
//! [Boundary condition]
//! [Push operators to pipeline]
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto get_body_force = [this](const double, const double, const double) {
t_source(i) = 0.;
t_source(0) = 0.1;
t_source(1) = 1.;
return t_source;
};
auto rho_ptr = boost::make_shared<double>(rho);
auto add_rho_block = [this, rho_ptr](auto &pip, auto block_name, Sev sev) {
struct OpMatRhoBlocks : public DomainEleOp {
OpMatRhoBlocks(boost::shared_ptr<double> rho_ptr,
MoFEM::Interface &m_field, Sev sev,
std::vector<const CubitMeshSets *> meshset_vec_ptr)
: DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), rhoPtr(rho_ptr) {
CHK_THROW_MESSAGE(extractRhoData(m_field, meshset_vec_ptr, sev),
"Can not get data from block");
}
MoFEMErrorCode doWork(int side, EntityType type,
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
*rhoPtr = b.rho;
}
}
*rhoPtr = rho;
}
private:
struct BlockData {
double rho;
Range blockEnts;
};
std::vector<BlockData> blockData;
extractRhoData(MoFEM::Interface &m_field,
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (auto m : meshset_vec_ptr) {
MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Rho Block") << *m;
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() < 1) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Expected that block has two attributes");
}
auto get_block_ents = [&]() {
Range ents;
m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
return ents;
};
MOFEM_TAG_AND_LOG("WORLD", sev, "Mat Rho Block")
<< m->getName() << ": ro = " << block_data[0];
blockData.push_back({block_data[0], get_block_ents()});
}
}
boost::shared_ptr<double> rhoPtr;
};
pip.push_back(new OpMatRhoBlocks(
rho_ptr, mField, sev,
// Get blockset using regular expression
(boost::format("%s(.*)") % block_name).str()
))
));
};
// Get pointer to U_tt shift in domain element
auto get_rho = [rho_ptr](const double, const double, const double) {
return *rho_ptr;
};
// specific time scaling
auto get_time_scale = [this](const double time) {
return sin(time * omega * M_PI);
};
auto apply_lhs = [&](auto &pip) {
"GEOMETRY");
CHKERR HenckyOps::opFactoryDomainLhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
mField, pip, "U", "MAT_ELASTIC", Sev::verbose);
CHKERR add_rho_block(pip, "MAT_RHO", Sev::verbose);
pip.push_back(new OpMass("U", "U", get_rho));
static_cast<OpMass &>(pip.back()).feScalingFun =
[](const FEMethod *fe_ptr) { return fe_ptr->ts_aa; };
};
auto apply_rhs = [&](auto &pip) {
pipeline_mng->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
CHKERR HenckyOps::opFactoryDomainRhs<SPACE_DIM, PETSC, GAUSS, DomainEleOp>(
mField, pip, "U", "MAT_ELASTIC", Sev::inform);
CHKERR add_rho_block(pip, "MAT_RHO", Sev::inform);
auto mat_acceleration_ptr = boost::make_shared<MatrixDouble>();
// Apply inertia
"U", mat_acceleration_ptr));
pip.push_back(new OpInertiaForce("U", mat_acceleration_ptr, get_rho));
pip, mField, "U", {},
{boost::make_shared<TimeScaleVector<SPACE_DIM>>("-time_vector_file",
true)},
"BODY_FORCE", Sev::inform);
};
CHKERR apply_lhs(pipeline_mng->getOpDomainLhsPipeline());
CHKERR apply_rhs(pipeline_mng->getOpDomainRhsPipeline());
auto integration_rule = [](int, int, int approx_order) {
return 2 * approx_order;
};
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
auto get_bc_hook = [&]() {
mField, pipeline_mng->getDomainRhsFE(),
{boost::make_shared<TimeScale>()});
return hook;
};
pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook();
}
//! [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");
}
}
private:
boost::shared_ptr<PostProcEle> postProc;
};
//! [Solve]
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto dm = simple->getDM();
ts = pipeline_mng->createTSIM2();
// Setup postprocessing
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
post_proc_fe->getOpPtrVector(), {H1});
auto common_ptr = commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
mField, post_proc_fe->getOpPtrVector(), "U", "MAT_ELASTIC", Sev::inform);
auto u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
auto X_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{},
{{"U", u_ptr}, {"GEOMETRY", X_ptr}},
{{"GRAD", common_ptr->matGradPtr},
{"FIRST_PIOLA", common_ptr->getMatFirstPiolaStress()}},
{}
)
);
// 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 TSSetMaxTime(ts, ftime);
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
auto T = createDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
SCATTER_FORWARD);
auto TT = vectorDuplicate(T);
CHKERR TS2SetSolution(ts, T, TT);
CHKERR TSSetFromOptions(ts);
CHKERR TSSolve(ts, NULL);
CHKERR TSGetTime(ts, &ftime);
PetscInt steps, snesfails, rejects, nonlinits, linits;
#if PETSC_VERSION_GE(3, 8, 0)
CHKERR TSGetStepNumber(ts, &steps);
#else
CHKERR TSGetTimeStepNumber(ts, &steps);
#endif
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 = createDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
SCATTER_FORWARD);
double nrm2;
CHKERR VecNorm(T, NORM_2, &nrm2);
MOFEM_LOG("EXAMPLE", Sev::inform) << "Regression norm " << nrm2;
constexpr double regression_value = 0.0194561;
if (fabs(nrm2 - regression_value) > 1e-2)
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"Regression test failed; 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";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
// Add logging channel for example
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::setLog("EXAMPLE");
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 interface
//! [Create MoFEM]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}
/**
* \file HenkyOps.hpp
* \example HenkyOps.hpp
*
* @copyright Copyright (c) 2023
*/
#ifndef __HENKY_OPS_HPP__
#define __HENKY_OPS_HPP__
namespace HenckyOps {
constexpr double eps = std::numeric_limits<float>::epsilon();
auto f = [](double v) { return 0.5 * std::log(v); };
auto d_f = [](double v) { return 0.5 / v; };
auto dd_f = [](double v) { return -0.5 / (v * v); };
struct isEq {
static inline auto check(const double &a, const double &b) {
return std::abs(a - b) / absMax < eps;
}
static double absMax;
};
double isEq::absMax = 1;
inline auto is_eq(const double &a, const double &b) {
return isEq::check(a, b);
};
template <int DIM> inline auto get_uniq_nb(double *ptr) {
std::array<double, DIM> tmp;
std::copy(ptr, &ptr[DIM], tmp.begin());
std::sort(tmp.begin(), tmp.end());
isEq::absMax = std::max(std::abs(tmp[0]), std::abs(tmp[DIM - 1]));
return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
};
template <int DIM>
std::max(std::max(std::abs(eig(0)), std::abs(eig(1))), std::abs(eig(2)));
int i = 0, j = 1, k = 2;
if (is_eq(eig(0), eig(1))) {
i = 0;
j = 2;
k = 1;
} else if (is_eq(eig(0), eig(2))) {
i = 0;
j = 1;
k = 2;
} else if (is_eq(eig(1), eig(2))) {
i = 1;
j = 0;
k = 2;
}
eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
FTensor::Tensor1<double, 3> eig_c{eig(i), eig(j), eig(k)};
{
eig(i) = eig_c(i);
eigen_vec(i, j) = eigen_vec_c(i, j);
}
};
struct CommonData : public boost::enable_shared_from_this<CommonData> {
boost::shared_ptr<MatrixDouble> matGradPtr;
boost::shared_ptr<MatrixDouble> matDPtr;
boost::shared_ptr<MatrixDouble> matLogCPlastic;
inline auto getMatFirstPiolaStress() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
inline auto getMatHenckyStress() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
inline auto getMatLogC() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matLogC);
}
inline auto getMatTangent() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &matTangent);
}
};
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpCalculateEigenValsImpl;
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpCalculateLogCImpl;
template <int DIM, IntegrationType I, typename DomainEleOp>
struct OpCalculateLogC_dCImpl;
template <int DIM, IntegrationType I, typename DomainEleOp, int S>
struct OpCalculateHenckyStressImpl;
template <int DIM, IntegrationType I, typename DomainEleOp, int S>
struct OpCalculateHenckyPlasticStressImpl;
template <int DIM, IntegrationType I, typename DomainEleOp, int S>
struct OpCalculatePiolaStressImpl;
template <int DIM, IntegrationType I, typename DomainEleOp, int S>
struct OpHenckyTangentImpl;
template <int DIM, typename DomainEleOp>
struct OpCalculateEigenValsImpl<DIM, GAUSS, DomainEleOp> : public DomainEleOp {
OpCalculateEigenValsImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data)
commonDataPtr(common_data) {
std::fill(&DomainEleOp::doEntities[MBEDGE],
&DomainEleOp::doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
commonDataPtr->matEigVal.resize(DIM, nb_gauss_pts, false);
commonDataPtr->matEigVec.resize(DIM * DIM, nb_gauss_pts, false);
auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
F(i, j) = t_grad(i, j) + t_kd(i, j);
C(i, j) = F(k, i) ^ F(k, j);
for (int ii = 0; ii != DIM; ii++)
for (int jj = 0; jj != DIM; jj++)
eigen_vec(ii, jj) = C(ii, jj);
for (auto ii = 0; ii != DIM; ++ii)
eig(ii) = std::max(std::numeric_limits<double>::epsilon(), eig(ii));
// rare case when two eigen values are equal
auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
if constexpr (DIM == 3) {
if (nb_uniq == 2) {
sort_eigen_vals<DIM>(eig, eigen_vec);
}
}
t_eig_val(i) = eig(i);
t_eig_vec(i, j) = eigen_vec(i, j);
++t_grad;
++t_eig_val;
++t_eig_vec;
}
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename DomainEleOp>
struct OpCalculateLogCImpl<DIM, GAUSS, DomainEleOp> : public DomainEleOp {
OpCalculateLogCImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data)
commonDataPtr(common_data) {
std::fill(&DomainEleOp::doEntities[MBEDGE],
&DomainEleOp::doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
commonDataPtr->matLogC.resize(size_symm, nb_gauss_pts, false);
auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
// rare case when two eigen values are equal
auto nb_uniq = get_uniq_nb<DIM>(&(t_eig_val(0)));
eig(i) = t_eig_val(i);
eigen_vec(i, j) = t_eig_vec(i, j);
auto logC = EigenMatrix::getMat(eig, eigen_vec, f);
t_logC(i, j) = logC(i, j);
++t_eig_val;
++t_eig_vec;
++t_logC;
}
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename DomainEleOp>
struct OpCalculateLogC_dCImpl<DIM, GAUSS, DomainEleOp> : public DomainEleOp {
OpCalculateLogC_dCImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data)
commonDataPtr(common_data) {
std::fill(&DomainEleOp::doEntities[MBEDGE],
&DomainEleOp::doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
commonDataPtr->matLogCdC.resize(size_symm * size_symm, nb_gauss_pts, false);
auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
eig(i) = t_eig_val(i);
eigen_vec(i, j) = t_eig_vec(i, j);
// rare case when two eigen values are equal
auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
auto dlogC_dC = EigenMatrix::getDiffMat(eig, eigen_vec, f, d_f, nb_uniq);
dlogC_dC(i, j, k, l) *= 2;
t_logC_dC(i, j, k, l) = dlogC_dC(i, j, k, l);
++t_logC_dC;
++t_eig_val;
++t_eig_vec;
}
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename DomainEleOp, int S>
struct OpCalculateHenckyStressImpl<DIM, GAUSS, DomainEleOp, S>
: public DomainEleOp {
OpCalculateHenckyStressImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data)
commonDataPtr(common_data) {
std::fill(&DomainEleOp::doEntities[MBEDGE],
&DomainEleOp::doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*commonDataPtr->matDPtr);
auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
commonDataPtr->matHenckyStress.resize(size_symm, nb_gauss_pts, false);
auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
t_T(i, j) = t_D(i, j, k, l) * t_logC(k, l);
++t_logC;
++t_T;
++t_D;
}
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename DomainEleOp, int S>
struct OpCalculateHenckyPlasticStressImpl<DIM, GAUSS, DomainEleOp, S>
: public DomainEleOp {
OpCalculateHenckyPlasticStressImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data,
boost::shared_ptr<MatrixDouble> mat_D_ptr,
const double scale = 1)
: DomainEleOp(field_name, DomainEleOp::OPROW), commonDataPtr(common_data),
scaleStress(scale), matDPtr(mat_D_ptr) {
std::fill(&DomainEleOp::doEntities[MBEDGE],
&DomainEleOp::doEntities[MBMAXTYPE], false);
matLogCPlastic = commonDataPtr->matLogCPlastic;
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*matDPtr);
auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
auto t_logCPlastic = getFTensor2SymmetricFromMat<DIM>(*matLogCPlastic);
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
commonDataPtr->matHenckyStress.resize(size_symm, nb_gauss_pts, false);
auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
t_T(i, j) = t_D(i, j, k, l) * (t_logC(k, l) - t_logCPlastic(k, l));
t_T(i, j) /= scaleStress;
++t_logC;
++t_T;
++t_D;
++t_logCPlastic;
}
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> matDPtr;
boost::shared_ptr<MatrixDouble> matLogCPlastic;
const double scaleStress;
};
template <int DIM, typename DomainEleOp, int S>
struct OpCalculatePiolaStressImpl<DIM, GAUSS, DomainEleOp, S>
: public DomainEleOp {
OpCalculatePiolaStressImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data)
commonDataPtr(common_data) {
std::fill(&DomainEleOp::doEntities[MBEDGE],
&DomainEleOp::doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
#ifdef HENCKY_SMALL_STRAIN
auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*commonDataPtr->matDPtr);
#endif
auto t_logC = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matLogC);
auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
commonDataPtr->matFirstPiolaStress.resize(DIM * DIM, nb_gauss_pts, false);
commonDataPtr->matSecondPiolaStress.resize(size_symm, nb_gauss_pts, false);
auto t_P = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matFirstPiolaStress);
auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
auto t_S =
getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
#ifdef HENCKY_SMALL_STRAIN
t_P(i, j) = t_D(i, j, k, l) * t_grad(k, l);
#else
t_F(i, j) = t_grad(i, j) + t_kd(i, j);
t_S(k, l) = t_T(i, j) * t_logC_dC(i, j, k, l);
t_P(i, l) = t_F(i, k) * t_S(k, l);
#endif
++t_grad;
++t_logC;
++t_logC_dC;
++t_P;
++t_T;
++t_S;
#ifdef HENCKY_SMALL_STRAIN
++t_D;
#endif
}
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename DomainEleOp, int S>
struct OpHenckyTangentImpl<DIM, GAUSS, DomainEleOp, S> : public DomainEleOp {
OpHenckyTangentImpl(const std::string field_name,
boost::shared_ptr<CommonData> common_data,
boost::shared_ptr<MatrixDouble> mat_D_ptr = nullptr)
commonDataPtr(common_data) {
std::fill(&DomainEleOp::doEntities[MBEDGE],
&DomainEleOp::doEntities[MBMAXTYPE], false);
if (mat_D_ptr)
matDPtr = mat_D_ptr;
else
matDPtr = commonDataPtr->matDPtr;
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = DomainEleOp::getGaussPts().size2();
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
commonDataPtr->matTangent.resize(DIM * DIM * DIM * DIM, nb_gauss_pts);
auto dP_dF =
getFTensor4FromMat<DIM, DIM, DIM, DIM, 1>(commonDataPtr->matTangent);
auto t_D = getFTensor4DdgFromMat<DIM, DIM, S>(*matDPtr);
auto t_eig_val = getFTensor1FromMat<DIM>(commonDataPtr->matEigVal);
auto t_eig_vec = getFTensor2FromMat<DIM, DIM>(commonDataPtr->matEigVec);
auto t_T = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matHenckyStress);
auto t_S =
getFTensor2SymmetricFromMat<DIM>(commonDataPtr->matSecondPiolaStress);
auto t_grad = getFTensor2FromMat<DIM, DIM>(*(commonDataPtr->matGradPtr));
auto t_logC_dC = getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->matLogCdC);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
#ifdef HENCKY_SMALL_STRAIN
dP_dF(i, j, k, l) = t_D(i, j, k, l);
#else
t_F(i, j) = t_grad(i, j) + t_kd(i, j);
eig(i) = t_eig_val(i);
eigen_vec(i, j) = t_eig_vec(i, j);
T(i, j) = t_T(i, j);
// rare case when two eigen values are equal
auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
dC_dF(i, j, k, l) = (t_kd(i, l) * t_F(k, j)) + (t_kd(j, l) * t_F(k, i));
auto TL =
EigenMatrix::getDiffDiffMat(eig, eigen_vec, f, d_f, dd_f, T, nb_uniq);
TL(i, j, k, l) *= 4;
P_D_P_plus_TL(i, j, k, l) =
TL(i, j, k, l) +
(t_logC_dC(i, j, o, p) * t_D(o, p, m, n)) * t_logC_dC(m, n, k, l);
P_D_P_plus_TL(i, j, k, l) *= 0.5;
dP_dF(i, j, m, n) = t_kd(i, m) * (t_kd(k, n) * t_S(k, j));
dP_dF(i, j, m, n) +=
t_F(i, k) * (P_D_P_plus_TL(k, j, o, p) * dC_dF(o, p, m, n));
#endif
++dP_dF;
++t_grad;
++t_eig_val;
++t_eig_vec;
++t_logC_dC;
++t_S;
++t_T;
++t_D;
}
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> matDPtr;
};
template <typename DomainEleOp> struct HenkyIntegrators {
template <int DIM, IntegrationType I>
using OpCalculateEigenVals = OpCalculateEigenValsImpl<DIM, I, DomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculateLogC = OpCalculateLogCImpl<DIM, I, DomainEleOp>;
template <int DIM, IntegrationType I>
using OpCalculateLogC_dC = OpCalculateLogC_dCImpl<DIM, I, DomainEleOp>;
template <int DIM, IntegrationType I, int S>
OpCalculateHenckyStressImpl<DIM, I, DomainEleOp, S>;
template <int DIM, IntegrationType I, int S = 0>
OpCalculateHenckyPlasticStressImpl<DIM, I, DomainEleOp, S>;
template <int DIM, IntegrationType I, int S>
OpCalculatePiolaStressImpl<DIM, GAUSS, DomainEleOp, S>;
template <int DIM, IntegrationType I, int S>
using OpHenckyTangent = OpHenckyTangentImpl<DIM, GAUSS, DomainEleOp, S>;
};
template <int DIM>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string block_name,
boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev,
double scale = 1) {
struct OpMatBlocks : public DomainEleOp {
OpMatBlocks(boost::shared_ptr<MatrixDouble> m, double bulk_modulus_K,
double shear_modulus_G, MoFEM::Interface &m_field, Sev sev,
std::vector<const CubitMeshSets *> meshset_vec_ptr,
double scale)
: DomainEleOp(NOSPACE, DomainEleOp::OPSPACE), matDPtr(m),
bulkModulusKDefault(bulk_modulus_K),
shearModulusGDefault(shear_modulus_G), scaleYoungModulus(scale) {
CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
"Can not get data from block");
}
MoFEMErrorCode doWork(int side, EntityType type,
for (auto &b : blockData) {
if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
CHKERR getMatDPtr(matDPtr, b.bulkModulusK * scaleYoungModulus,
b.shearModulusG * scaleYoungModulus);
}
}
CHKERR getMatDPtr(matDPtr, bulkModulusKDefault * scaleYoungModulus,
shearModulusGDefault * scaleYoungModulus);
}
private:
boost::shared_ptr<MatrixDouble> matDPtr;
const double scaleYoungModulus;
struct BlockData {
double bulkModulusK;
double shearModulusG;
Range blockEnts;
};
double bulkModulusKDefault;
double shearModulusGDefault;
std::vector<BlockData> blockData;
extractBlockData(MoFEM::Interface &m_field,
std::vector<const CubitMeshSets *> meshset_vec_ptr,
Sev sev) {
for (auto m : meshset_vec_ptr) {
MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
std::vector<double> block_data;
CHKERR m->getAttributes(block_data);
if (block_data.size() != 2) {
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Expected that block has two attribute");
}
auto get_block_ents = [&]() {
Range ents;
m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
return ents;
};
double young_modulus = block_data[0];
double poisson_ratio = block_data[1];
double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
<< "E = " << young_modulus << " nu = " << poisson_ratio;
blockData.push_back(
{bulk_modulus_K, shear_modulus_G, get_block_ents()});
}
}
MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
double bulk_modulus_K, double shear_modulus_G) {
//! [Calculate elasticity tensor]
auto set_material_stiffness = [&]() {
double A = (DIM == 2)
: 1;
auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
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);
};
//! [Calculate elasticity tensor]
constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
mat_D_ptr->resize(size_symm * size_symm, 1);
set_material_stiffness();
}
};
double E = young_modulus;
double nu = poisson_ratio;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "hencky_", "", "none");
CHKERR PetscOptionsScalar("-young_modulus", "Young modulus", "", E, &E,
PETSC_NULL);
CHKERR PetscOptionsScalar("-poisson_ratio", "poisson ratio", "", nu, &nu,
PETSC_NULL);
ierr = PetscOptionsEnd();
double bulk_modulus_K = E / (3 * (1 - 2 * nu));
double shear_modulus_G = E / (2 * (1 + nu));
pip.push_back(new OpMatBlocks(
mat_D_Ptr, bulk_modulus_K, shear_modulus_G, m_field, sev,
// Get blockset using regular expression
m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
(boost::format("%s(.*)") % block_name).str()
)),
));
}
template <int DIM, IntegrationType I, typename DomainEleOp>
MoFEM::Interface &m_field,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string field_name, std::string block_name, Sev sev, double scale = 1) {
auto common_ptr = boost::make_shared<HenckyOps::CommonData>();
common_ptr->matDPtr = boost::make_shared<MatrixDouble>();
common_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
CHK_THROW_MESSAGE(addMatBlockOps<DIM>(m_field, pip, block_name,
common_ptr->matDPtr, sev, scale),
"addMatBlockOps");
using H = HenkyIntegrators<DomainEleOp>;
field_name, common_ptr->matGradPtr));
pip.push_back(new typename H::template OpCalculateEigenVals<DIM, I>(
field_name, common_ptr));
pip.push_back(
new typename H::template OpCalculateLogC<DIM, I>(field_name, common_ptr));
pip.push_back(new typename H::template OpCalculateLogC_dC<DIM, I>(
field_name, common_ptr));
// Assumes constant D matrix per entity
pip.push_back(new typename H::template OpCalculateHenckyStress<DIM, I, 0>(
field_name, common_ptr));
pip.push_back(new typename H::template OpCalculatePiolaStress<DIM, I, 0>(
field_name, common_ptr));
return common_ptr;
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEM::Interface &m_field,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
Sev sev) {
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template LinearForm<I>;
typename B::template OpGradTimesTensor<1, DIM, DIM>;
pip.push_back(
new OpInternalForcePiola("U", common_ptr->getMatFirstPiolaStress()));
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEM::Interface &m_field,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string field_name, std::string block_name, Sev sev, double scale = 1) {
auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
m_field, pip, field_name, block_name, sev, scale);
CHKERR opFactoryDomainRhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
common_ptr, sev);
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEM::Interface &m_field,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string field_name, boost::shared_ptr<HenckyOps::CommonData> common_ptr,
Sev sev) {
using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
A>::template BiLinearForm<I>;
using OpKPiola = typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
using H = HenkyIntegrators<DomainEleOp>;
// Assumes constant D matrix per entity
pip.push_back(
new typename H::template OpHenckyTangent<DIM, I, 0>(field_name, common_ptr));
pip.push_back(
new OpKPiola(field_name, field_name, common_ptr->getMatTangent()));
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
MoFEM::Interface &m_field,
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string field_name, std::string block_name, Sev sev, double scale = 1) {
auto common_ptr = commonDataFactory<DIM, I, DomainEleOp>(
m_field, pip, field_name, block_name, sev, scale);
CHKERR opFactoryDomainLhs<DIM, A, I, DomainEleOp>(m_field, pip, field_name,
common_ptr, sev);
}
} // namespace HenckyOps
#endif // __HENKY_OPS_HPP__
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::NaturalBC::Assembly::LinearForm
Definition: Natural.hpp:67
EigenMatrix::getMat
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
Definition: MatrixFunction.cpp:53
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
Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: dynamic_first_order_con_law.cpp:1205
HenckyOps::CommonData::matLogCPlastic
boost::shared_ptr< MatrixDouble > matLogCPlastic
Definition: HenckyOps.hpp:83
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
omega
constexpr double omega
Save field DOFS on vertices/tags.
Definition: dynamic_first_order_con_law.cpp:93
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:647
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
HenckyOps::CommonData::getMatFirstPiolaStress
auto getMatFirstPiolaStress()
Definition: HenckyOps.hpp:94
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
rho
double rho
Definition: plastic.cpp:140
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
HenckyOps::HenkyIntegrators::OpCalculateLogC_dC
OpCalculateLogC_dCImpl< DIM, I, DomainEleOp > OpCalculateLogC_dC
Definition: HenckyOps.hpp:581
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
E
HenckyOps::isEq::absMax
static double absMax
Definition: HenckyOps.hpp:23
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
HenckyOps::CommonData::matEigVec
MatrixDouble matEigVec
Definition: HenckyOps.hpp:86
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
HenckyOps::CommonData::getMatLogC
auto getMatLogC()
Definition: HenckyOps.hpp:104
HenckyOps::opFactoryDomainLhs
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
Definition: HenckyOps.hpp:811
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
OpKPiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: seepage.cpp:64
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
bulk_modulus_K
double bulk_modulus_K
Definition: dynamic_first_order_con_law.cpp:96
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
HenckyOps::CommonData::matEigVal
MatrixDouble matEigVal
Definition: HenckyOps.hpp:85
HenckyOps
Definition: HenckyOps.hpp:11
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
scale
double scale
Definition: plastic.cpp:119
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
HenckyOps.hpp
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
HenckyOps::CommonData::matLogC
MatrixDouble matLogC
Definition: HenckyOps.hpp:87
HenckyOps::HenkyIntegrators::OpCalculateLogC
OpCalculateLogCImpl< DIM, I, DomainEleOp > OpCalculateLogC
Definition: HenckyOps.hpp:578
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
HenckyOps::isEq::check
static auto check(const double &a, const double &b)
Definition: HenckyOps.hpp:20
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
HenckyOps::is_eq
auto is_eq(const double &a, const double &b)
Definition: HenckyOps.hpp:28
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: dynamic_first_order_con_law.cpp:463
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: dynamic_first_order_con_law.cpp:63
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
MoFEM::Sev
MoFEM::LogManager::SeverityLevel Sev
Definition: CoreTemplates.hpp:17
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::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: seepage.cpp:66
BlockData
Definition: ElasticityMixedFormulation.hpp:12
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
Example
[Example]
Definition: plastic.cpp:177
double
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
HenckyOps::eps
constexpr double eps
Definition: HenckyOps.hpp:13
PlasticOps::Monitor
Definition: PlasticOpsMonitor.hpp:9
Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: dynamic_first_order_con_law.cpp:893
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
HenckyOps::CommonData::getMatHenckyStress
auto getMatHenckyStress()
Definition: HenckyOps.hpp:99
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MatrixFunction.hpp
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
HenckyOps::CommonData::matFirstPiolaStress
MatrixDouble matFirstPiolaStress
Definition: HenckyOps.hpp:89
HenckyOps::CommonData::matSecondPiolaStress
MatrixDouble matSecondPiolaStress
Definition: HenckyOps.hpp:90
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:17
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
MoFEM::LogManager::SeverityLevel
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::AddFluxToRhsPipelineImpl
Definition: Natural.hpp:46
size_symm
constexpr auto size_symm
Definition: plastic.cpp:42
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
main
int main(int argc, char *argv[])
Definition: activate_deactivate_dofs.cpp:15
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', SPACE_DIM >
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
HenckyOps::HenkyIntegrators::OpCalculateHenckyPlasticStress
OpCalculateHenckyPlasticStressImpl< DIM, I, DomainEleOp, S > OpCalculateHenckyPlasticStress
Definition: HenckyOps.hpp:589
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
convert.n
n
Definition: convert.py:82
HenckyOps::CommonData::matLogCdC
MatrixDouble matLogCdC
Definition: HenckyOps.hpp:88
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
ElementsAndOps
Definition: child_and_parent.cpp:18
Range
DomainEleOp
HenckyOps::sort_eigen_vals
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
Definition: HenckyOps.hpp:41
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
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
HenckyOps::CommonData::matGradPtr
boost::shared_ptr< MatrixDouble > matGradPtr
Definition: HenckyOps.hpp:81
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition: Natural.hpp:65
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
HenckyOps::CommonData::matDPtr
boost::shared_ptr< MatrixDouble > matDPtr
Definition: HenckyOps.hpp:82
shear_modulus_G
double shear_modulus_G
Definition: dynamic_first_order_con_law.cpp:97
HenckyOps::CommonData::getMatTangent
auto getMatTangent()
Definition: HenckyOps.hpp:108
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
HenckyOps::HenkyIntegrators::OpCalculateHenckyStress
OpCalculateHenckyStressImpl< DIM, I, DomainEleOp, S > OpCalculateHenckyStress
Definition: HenckyOps.hpp:585
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
MoFEM::TSMethod::ts_aa
PetscReal ts_aa
shift for U_tt shift for U_tt
Definition: LoopMethods.hpp:161
HenckyOps::HenkyIntegrators::OpCalculateEigenVals
OpCalculateEigenValsImpl< DIM, I, DomainEleOp > OpCalculateEigenVals
Definition: HenckyOps.hpp:575
HenckyOps::opFactoryDomainRhs
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
Definition: HenckyOps.hpp:778
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:225
ts_time_ptr
static double * ts_time_ptr
Definition: nonlinear_dynamic_elastic.cpp:50
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
HenckyOps::HenkyIntegrators::OpHenckyTangent
OpHenckyTangentImpl< DIM, GAUSS, DomainEleOp, S > OpHenckyTangent
Definition: HenckyOps.hpp:596
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: dynamic_first_order_con_law.cpp:536
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FTensor::Ddg
Definition: Ddg_value.hpp:7
approx_order
int approx_order
Definition: test_broken_space.cpp:50
CommonData
Definition: continuity_check_on_skeleton_with_simple_2d_for_h1.cpp:22
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
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
EigenMatrix::getDiffDiffMat
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
Definition: MatrixFunction.cpp:78
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::BlockData
Definition: MeshsetsManager.cpp:755
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:213
HenckyOps::commonDataFactory
auto commonDataFactory(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
Definition: HenckyOps.hpp:745
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
HenckyOps::HenkyIntegrators::OpCalculatePiolaStress
OpCalculatePiolaStressImpl< DIM, GAUSS, DomainEleOp, S > OpCalculatePiolaStress
Definition: HenckyOps.hpp:593
HenckyOps::addMatBlockOps
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev, double scale=1)
Definition: HenckyOps.hpp:601
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::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< DM >
HenckyOps::CommonData::matTangent
MatrixDouble matTangent
Definition: HenckyOps.hpp:92
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
Example::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition: dynamic_first_order_con_law.cpp:1183
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
EigenMatrix::getDiffMat
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
Definition: MatrixFunction.cpp:64
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
HenckyOps::CommonData::matHenckyStress
MatrixDouble matHenckyStress
Definition: HenckyOps.hpp:91
HenckyOps::get_uniq_nb
auto get_uniq_nb(double *ptr)
Definition: HenckyOps.hpp:32
H
double H
Hardening.
Definition: plastic.cpp:124
ts_aa_ptr
static double * ts_aa_ptr
Definition: nonlinear_dynamic_elastic.cpp:51
OpCalculateVectorFieldGradient
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::computeEigenValuesSymmetric
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1452
F
@ F
Definition: free_surface.cpp:394
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698