v0.13.1
Loading...
Searching...
No Matches
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>
#include <MatrixFunction.hpp>
using namespace MoFEM;
template <int DIM> struct ElementsAndOps {};
template <> struct ElementsAndOps<2> {
using DomainEle = FaceElementForcesAndSourcesCore;
};
template <> struct ElementsAndOps<3> {
using DomainEle = VolumeElementForcesAndSourcesCore;
};
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
using EntData = EntitiesFieldData::EntData;
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
using OpK = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
using OpMass = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
GAUSS>::OpMass<1, SPACE_DIM>;
using OpInternalForce = FormsIntegrators<DomainEleOp>::Assembly<
PETSC>::LinearForm<GAUSS>::OpGradTimesTensor<1, SPACE_DIM, SPACE_DIM>;
using OpBodyForce = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::LinearForm<
GAUSS>::OpSource<1, SPACE_DIM>;
using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
PETSC>::LinearForm<GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 1>;
constexpr bool is_quasi_static = false;
constexpr double rho = 1;
constexpr double omega = 2.4;
constexpr double young_modulus = 1;
constexpr double poisson_ratio = 0.25;
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;
static double *ts_time_ptr;
static double *ts_aa_ptr;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
MoFEMErrorCode runProblem();
private:
MoFEM::Interface &mField;
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> matAccelerationPtr;
boost::shared_ptr<MatrixDouble> matInertiaPtr;
boost::shared_ptr<MatrixDouble> matDPtr;
boost::shared_ptr<MatrixDouble> matTangentPtr;
};
//! [Create common data]
MoFEMErrorCode Example::createCommonData() {
MoFEMFunctionBegin;
auto set_matrial_stiffness = [&]() {
FTensor::Index<'i', SPACE_DIM> i;
FTensor::Index<'j', SPACE_DIM> j;
FTensor::Index<'k', SPACE_DIM> k;
FTensor::Index<'l', SPACE_DIM> l;
constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
MoFEMFunctionBegin;
constexpr double A =
(SPACE_DIM == 2) ? 2 * shear_modulus_G /
(bulk_modulus_K + (4. / 3.) * 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);
MoFEMFunctionReturn(0);
};
matGradPtr = boost::make_shared<MatrixDouble>();
matStrainPtr = boost::make_shared<MatrixDouble>();
matStressPtr = boost::make_shared<MatrixDouble>();
matAccelerationPtr = boost::make_shared<MatrixDouble>();
matInertiaPtr = boost::make_shared<MatrixDouble>();
matDPtr = boost::make_shared<MatrixDouble>();
matTangentPtr = boost::make_shared<MatrixDouble>();
constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
matDPtr->resize(size_symm * size_symm, 1);
CHKERR set_matrial_stiffness();
MoFEMFunctionReturn(0);
}
//! [Create common data]
//! [Run problem]
MoFEMErrorCode Example::runProblem() {
MoFEMFunctionBegin;
CHKERR readMesh();
CHKERR setupProblem();
CHKERR createCommonData();
CHKERR boundaryCondition();
CHKERR assembleSystem();
CHKERR solveSystem();
CHKERR outputResults();
CHKERR checkResults();
MoFEMFunctionReturn(0);
}
//! [Run problem]
//! [Read mesh]
MoFEMErrorCode Example::readMesh() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
CHKERR simple->getOptions();
CHKERR simple->loadFile();
MoFEMFunctionReturn(0);
}
//! [Read mesh]
//! [Set up problem]
MoFEMErrorCode Example::setupProblem() {
MoFEMFunctionBegin;
Simple *simple = mField.getInterface<Simple>();
// Add field
CHKERR simple->addDomainField("U", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
SPACE_DIM);
int order = 2;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
CHKERR simple->setFieldOrder("U", order);
CHKERR simple->setUp();
MoFEMFunctionReturn(0);
}
//! [Set up problem]
//! [Boundary condition]
MoFEMErrorCode Example::boundaryCondition() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
MoFEMFunctionReturn(0);
}
//! [Boundary condition]
//! [Push operators to pipeline]
MoFEMErrorCode Example::assembleSystem() {
MoFEMFunctionBegin;
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto add_base_transform = [&](auto &pipeline) {
pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
pipeline.push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(
new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
pipeline.push_back(new OpSetHOWeights(det_ptr));
};
add_base_transform(pipeline_mng->getOpDomainRhsPipeline());
add_base_transform(pipeline_mng->getOpDomainLhsPipeline());
// Get pointer to U_tt shift in domain element
auto get_rho = [this](const double, const double, const double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
return rho * fe_domain_lhs->ts_aa;
};
auto get_body_force = [this](const double, const double, const double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
FTensor::Index<'i', SPACE_DIM> i;
FTensor::Tensor1<double, SPACE_DIM> t_source;
t_source(i) = 0.;
t_source(0) = 0.1;
t_source(1) = 1.;
const auto time = fe_domain_rhs->ts_t;
t_source(i) *= sin(time * omega * M_PI);
return t_source;
};
auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
henky_common_data_ptr->matGradPtr = matGradPtr;
henky_common_data_ptr->matDPtr = matDPtr;
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
matGradPtr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpHenckyTangent<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpK("U", "U", henky_common_data_ptr->getMatTangent()));
if (!is_quasi_static)
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpMass("U", "U", get_rho));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpBodyForce("U", get_body_force));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
matGradPtr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpInternalForce(
"U", henky_common_data_ptr->getMatFirstPiolaStress()));
if (!is_quasi_static) {
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>("U",
matAccelerationPtr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpScaleMatrix("U", rho, matAccelerationPtr, matInertiaPtr));
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpInertiaForce(
"U", matInertiaPtr, [](double, double, double) { return 1.; }));
}
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 = [&]() {
EssentialPreProc<DisplacementCubitBcData> hook(
mField, pipeline_mng->getDomainRhsFE(),
{boost::make_shared<TimeScale>()});
return hook;
};
pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook();
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, boost::shared_ptr<PostProcEle> post_proc)
: dM(dm), postProc(post_proc){};
MoFEMErrorCode postProcess() {
MoFEMFunctionBegin;
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");
}
MoFEMFunctionReturn(0);
}
private:
SmartPetscObj<DM> dM;
boost::shared_ptr<PostProcEle> postProc;
};
//! [Solve]
MoFEMErrorCode Example::solveSystem() {
MoFEMFunctionBegin;
auto *simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto dm = simple->getDM();
MoFEM::SmartPetscObj<TS> ts;
if (is_quasi_static)
ts = pipeline_mng->createTSIM();
else
ts = pipeline_mng->createTSIM2();
// Setup postprocessing
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>("U",
matGradPtr));
post_proc_fe->getOpPtrVector().push_back(
new OpSymmetrizeTensor<SPACE_DIM>("U", matGradPtr, matStrainPtr));
post_proc_fe->getOpPtrVector().push_back(
new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
"U", matStrainPtr, matStressPtr, matDPtr));
auto u_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{},
{{"U", u_ptr}},
{},
{{"STRAIN", matStrainPtr}, {"STRESS", matStressPtr}}
)
);
// 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 = smartCreateDMVector(simple->getDM());
CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
SCATTER_FORWARD);
if (is_quasi_static) {
CHKERR TSSetSolution(ts, T);
CHKERR TSSetFromOptions(ts);
} else {
auto TT = smartVectorDuplicate(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);
MoFEMFunctionReturn(0);
}
//! [Solve]
//! [Postprocess results]
MoFEMErrorCode Example::outputResults() {
MoFEMFunctionBegin;
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,
SCATTER_FORWARD);
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)
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
"Regression test faileed; wrong norm value.");
}
MoFEMFunctionReturn(0);
}
//! [Postprocess results]
//! [Check]
MoFEMErrorCode Example::checkResults() {
MoFEMFunctionBegin;
MoFEMFunctionReturn(0);
}
//! [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::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
MOFEM_LOG_TAG("EXAMPLE", "example");
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 insterface
//! [Create MoFEM]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
CATCH_ERRORS;
CHKERR MoFEM::Core::Finalize();
}
param_file
std::string param_file
Definition: DefaultParams.hpp:113
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
help
static char help[]
Definition: activate_deactivate_dofs.cpp:13
main
int main()
Definition: adol-c_atom.cpp:46
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
is_quasi_static
constexpr bool is_quasi_static
Definition: contact.cpp:117
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
OpInternalForce
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
Definition: elastic.cpp:44
shear_modulus_G
constexpr double shear_modulus_G
Definition: elastic.cpp:63
bulk_modulus_K
constexpr double bulk_modulus_K
Definition: elastic.cpp:62
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
Definition: fluid_structure_eigenproblem.cpp:59
poisson_ratio
auto poisson_ratio
Definition: fluid_structure_eigenproblem.cpp:89
PostProcEle
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
Definition: fluid_structure_eigenproblem.cpp:48
order
int order
Definition: fluid_structure_eigenproblem.cpp:84
integration_rule
auto integration_rule
Definition: free_surface.cpp:111
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:73
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: DMMMoFEM.cpp:470
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
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: DMMMoFEM.cpp:533
MoFEM::smartCreateDMVector
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:97
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:90
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
T
const double T
Definition: mixed_reac_diff.cpp:158
OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
Definition: multifield_plasticity.cpp:73
OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
Definition: multifield_plasticity.cpp:50
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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: DMMMoFEM.cpp:1003
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::smartVectorDuplicate
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:200
ts_aa_ptr
static double * ts_aa_ptr
Definition: nonlinear_dynamic_elastic.cpp:54
omega
constexpr double omega
Definition: nonlinear_dynamic_elastic.cpp:44
ts_time_ptr
static double * ts_time_ptr
Definition: nonlinear_dynamic_elastic.cpp:53
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
young_modulus
double young_modulus
Definition: plastic.cpp:96
rho
double rho
Definition: plastic.cpp:98
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:10
size_symm
constexpr auto size_symm
Definition: plastic.cpp:33
A
constexpr AssemblyType A
Definition: plastic.cpp:35
approx_order
static constexpr int approx_order
Definition: prism_polynomial_approximation.cpp:14
l
FTensor::Index< 'l', 3 > l
Definition: shallow_wave.cpp:79
j
FTensor::Index< 'j', 3 > j
Definition: shallow_wave.cpp:78
omega
constexpr double omega
Definition: shallow_wave.cpp:62
i
FTensor::Index< 'i', 3 > i
Definition: shallow_wave.cpp:77
Example
[Example]
Definition: plastic.cpp:137
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: helmholtz.cpp:106
Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Applying essential BC]
Definition: helmholtz.cpp:154
Example::readMesh
MoFEMErrorCode readMesh()
[run problem]
Definition: helmholtz.cpp:73
Example::matAccelerationPtr
boost::shared_ptr< MatrixDouble > matAccelerationPtr
Definition: nonlinear_dynamic_elastic.cpp:77
Example::matInertiaPtr
boost::shared_ptr< MatrixDouble > matInertiaPtr
Definition: nonlinear_dynamic_elastic.cpp:78
Example::matGradPtr
boost::shared_ptr< MatrixDouble > matGradPtr
Definition: nonlinear_elastic.cpp:76
Example::matStressPtr
boost::shared_ptr< MatrixDouble > matStressPtr
Definition: nonlinear_elastic.cpp:78
Example::matStrainPtr
boost::shared_ptr< MatrixDouble > matStrainPtr
Definition: nonlinear_elastic.cpp:77
Example::checkResults
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: contact.cpp:659
Example::solveSystem
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Definition: helmholtz.cpp:235
Example::createCommonData
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:257
Example::matTangentPtr
boost::shared_ptr< MatrixDouble > matTangentPtr
Definition: nonlinear_dynamic_elastic.cpp:81
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:172
Example::mField
MoFEM::Interface & mField
Definition: plastic.cpp:144
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:184
Example::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition: helmholtz.cpp:255
Example::matDPtr
boost::shared_ptr< MatrixDouble > matDPtr
Definition: eigen_elastic.cpp:65
HenckyOps::OpCalculateLogC_dC
Definition: HenckyOps.hpp:208
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:23
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
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::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:72
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
MoFEM::EssentialPreProc
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:39
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field valuse for given petsc vector.
Definition: UserDataOperators.hpp:547
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:419
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:583
MoFEM::OpSetHOInvJacToScalarBases
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:73
MoFEM::OpSetHOWeights
Set inverse jacobian to base functions.
Definition: HODataOperators.hpp:133
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:552
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition: PetscSmartObj.hpp:81
MoFEM::TSMethod::ts_step
PetscInt ts_step
time step number
Definition: LoopMethods.hpp:159
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
Monitor
Monitor solution.
Definition: heat_equation.cpp:68
Monitor::dM
SmartPetscObj< DM > dM
Definition: heat_equation.cpp:89
Monitor::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: heat_equation.cpp:78
Monitor::postProc
boost::shared_ptr< PostProcEle > postProc
Definition: heat_equation.cpp:90
namespace HenckyOps {
constexpr double eps = std::numeric_limits<float>::epsilon();
auto f = [](double v) { return 0.5 * std::log(static_cast<long double>(v)); };
auto d_f = [](double v) { return 0.5 / static_cast<long double>(v); };
auto dd_f = [](double v) {
return -0.5 / (static_cast<long double>(v) * static_cast<long double>(v));
};
inline auto is_eq(const double &a, const double &b) {
return std::abs(a - b) < eps;
};
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());
return std::distance(tmp.begin(), std::unique(tmp.begin(), tmp.end(), is_eq));
};
template <int DIM>
inline auto sort_eigen_vals(FTensor::Tensor1<double, DIM> &eig,
FTensor::Tensor2<double, DIM, DIM> &eigen_vec) {
std::array<std::pair<double, size_t>, DIM> tmp;
auto is_eq_pair = [](const auto &a, const auto &b) {
return a.first < b.first;
};
for (size_t i = 0; i != DIM; ++i)
tmp[i] = std::make_pair(eig(i), i);
std::sort(tmp.begin(), tmp.end(), is_eq_pair);
int i, j, k;
if (is_eq_pair(tmp[0], tmp[1])) {
i = tmp[0].second;
j = tmp[2].second;
k = tmp[1].second;
} else {
i = tmp[0].second;
j = tmp[1].second;
k = tmp[2].second;
}
FTensor::Tensor2<double, 3, 3> eigen_vec_c{
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)};
{
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
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;
MatrixDouble matEigVal;
MatrixDouble matEigVec;
MatrixDouble matLogC;
MatrixDouble matLogCdC;
MatrixDouble matFirstPiolaStress;
MatrixDouble matSecondPiolaStress;
MatrixDouble matHenckyStress;
MatrixDouble matTangent;
inline auto getMatFirstPiolaStress() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
&matFirstPiolaStress);
}
inline auto getMatHenckyStress() {
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
&matHenckyStress);
}
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> struct OpCalculateEigenVals : public DomainEleOp {
OpCalculateEigenVals(const std::string field_name,
boost::shared_ptr<CommonData> common_data)
: DomainEleOp(field_name, DomainEleOp::OPROW),
commonDataPtr(common_data) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = 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) {
FTensor::Tensor2<double, DIM, DIM> F;
FTensor::Tensor2_symmetric<double, DIM> C;
FTensor::Tensor1<double, DIM> eig;
FTensor::Tensor2<double, DIM, DIM> eigen_vec;
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);
CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
// rare case when two eigen values are equal
auto nb_uniq = get_uniq_nb<DIM>(&eig(0));
if (DIM == 3 && 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;
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM> struct OpCalculateLogC : public DomainEleOp {
OpCalculateLogC(const std::string field_name,
boost::shared_ptr<CommonData> common_data)
: DomainEleOp(field_name, DomainEleOp::OPROW),
commonDataPtr(common_data) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = 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)));
FTensor::Tensor1<double, DIM> eig;
FTensor::Tensor2<double, DIM, DIM> eigen_vec;
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;
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM> struct OpCalculateLogC_dC : public DomainEleOp {
OpCalculateLogC_dC(const std::string field_name,
boost::shared_ptr<CommonData> common_data)
: DomainEleOp(field_name, DomainEleOp::OPROW),
commonDataPtr(common_data) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
FTensor::Index<'l', DIM> l;
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = 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) {
FTensor::Tensor1<double, DIM> eig;
FTensor::Tensor2<double, DIM, DIM> eigen_vec;
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;
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM> struct OpCalculateHenckyStress : public DomainEleOp {
OpCalculateHenckyStress(const std::string field_name,
boost::shared_ptr<CommonData> common_data)
: DomainEleOp(field_name, DomainEleOp::OPROW),
commonDataPtr(common_data) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
FTensor::Index<'l', DIM> l;
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = getGaussPts().size2();
auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*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;
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM> struct OpCalculateHenckyPlasticStress : public DomainEleOp {
OpCalculateHenckyPlasticStress(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(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
matLogCPlastic = commonDataPtr->matLogCPlastic;
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
FTensor::Index<'l', DIM> l;
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = getGaussPts().size2();
auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*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;
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> matDPtr;
boost::shared_ptr<MatrixDouble> matLogCPlastic;
const double scaleStress;
};
template <int DIM> struct OpCalculatePiolaStress : public DomainEleOp {
OpCalculatePiolaStress(const std::string field_name,
boost::shared_ptr<CommonData> common_data)
: DomainEleOp(field_name, DomainEleOp::OPROW),
commonDataPtr(common_data) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
FTensor::Index<'l', DIM> l;
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = getGaussPts().size2();
auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->matDPtr);
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) {
FTensor::Tensor2<double, DIM, DIM> t_F;
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);
++t_grad;
++t_logC;
++t_logC_dC;
++t_P;
++t_T;
++t_S;
++t_D;
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM> struct OpHenckyTangent : public DomainEleOp {
OpHenckyTangent(const std::string field_name,
boost::shared_ptr<CommonData> common_data,
boost::shared_ptr<MatrixDouble> mat_D_ptr = nullptr)
: DomainEleOp(field_name, DomainEleOp::OPROW),
commonDataPtr(common_data) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
if (mat_D_ptr)
matDPtr = mat_D_ptr;
else
matDPtr = commonDataPtr->matDPtr;
}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
FTensor::Index<'l', DIM> l;
FTensor::Index<'m', DIM> m;
FTensor::Index<'n', DIM> n;
FTensor::Index<'o', DIM> o;
FTensor::Index<'p', DIM> p;
constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
// const size_t nb_gauss_pts = matGradPtr->size2();
const size_t nb_gauss_pts = 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, 0>(*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) {
FTensor::Tensor2<double, DIM, DIM> t_F;
t_F(i, j) = t_grad(i, j) + t_kd(i, j);
FTensor::Tensor1<double, DIM> eig;
FTensor::Tensor2<double, DIM, DIM> eigen_vec;
FTensor::Tensor2_symmetric<double, DIM> T;
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));
FTensor::Tensor4<double, DIM, DIM, DIM, DIM> dC_dF;
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;
FTensor::Ddg<double, DIM, DIM> P_D_P_plus_TL;
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));
++dP_dF;
++t_grad;
++t_eig_val;
++t_eig_vec;
++t_logC_dC;
++t_S;
++t_T;
++t_D;
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<CommonData> commonDataPtr;
boost::shared_ptr<MatrixDouble> matDPtr;
};
template <int DIM> struct OpPostProcHencky : public DomainEleOp {
OpPostProcHencky(const std::string field_name,
moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts,
boost::shared_ptr<CommonData> common_data_ptr);
MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
private:
moab::Interface &postProcMesh;
std::vector<EntityHandle> &mapGaussPts;
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM>
OpPostProcHencky<DIM>::OpPostProcHencky(
const std::string field_name, moab::Interface &post_proc_mesh,
std::vector<EntityHandle> &map_gauss_pts,
boost::shared_ptr<CommonData> common_data_ptr)
: DomainEleOp(field_name, DomainEleOp::OPROW), postProcMesh(post_proc_mesh),
mapGaussPts(map_gauss_pts), commonDataPtr(common_data_ptr) {
// Opetor is only executed for vertices
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
//! [Postprocessing]
template <int DIM>
MoFEMErrorCode OpPostProcHencky<DIM>::doWork(int side, EntityType type,
EntData &data) {
MoFEMFunctionBegin;
FTensor::Index<'i', DIM> i;
FTensor::Index<'j', DIM> j;
FTensor::Index<'k', DIM> k;
std::array<double, 9> def;
std::fill(def.begin(), def.end(), 0);
auto get_tag = [&](const std::string name, size_t size) {
Tag th;
CHKERR postProcMesh.tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
MB_TAG_CREAT | MB_TAG_SPARSE,
def.data());
return th;
};
MatrixDouble3by3 mat(3, 3);
auto set_matrix = [&](auto &t) -> MatrixDouble3by3 & {
mat.clear();
for (size_t r = 0; r != SPACE_DIM; ++r)
for (size_t c = 0; c != SPACE_DIM; ++c)
mat(r, c) = t(r, c);
return mat;
};
auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
&*mat.data().begin());
};
auto th_grad = get_tag("GRAD", 9);
auto th_stress = get_tag("STRESS", 9);
auto t_piola =
getFTensor2FromMat<DIM, DIM>(commonDataPtr->matFirstPiolaStress);
auto t_grad = getFTensor2FromMat<DIM, DIM>(*commonDataPtr->matGradPtr);
FTensor::Tensor2<double, DIM, DIM> F;
FTensor::Tensor2_symmetric<double, DIM> cauchy_stress;
for (int gg = 0; gg != commonDataPtr->matGradPtr->size2(); ++gg) {
F(i, j) = t_grad(i, j) + kronecker_delta(i, j);
double inv_t_detF;
if (DIM == 2)
determinantTensor2by2(F, inv_t_detF);
else
determinantTensor3by3(F, inv_t_detF);
inv_t_detF = 1. / inv_t_detF;
cauchy_stress(i, j) = inv_t_detF * (t_piola(i, k) ^ F(j, k));
CHKERR set_tag(th_grad, gg, set_matrix(t_grad));
CHKERR set_tag(th_stress, gg, set_matrix(cauchy_stress));
++t_piola;
++t_grad;
}
MoFEMFunctionReturn(0);
}
} // namespace HenckyOps
o
static Index< 'o', 3 > o
Definition: BasicFeTools.hpp:89
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
a
constexpr double a
Definition: approx_sphere.cpp:30
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:14
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
m
FTensor::Index< 'm', SPACE_DIM > m
Definition: fluid_structure_eigenproblem.cpp:65
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
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
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
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
FTensor::kronecker_delta
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
Definition: Kronecker_Delta.hpp:81
HenckyOps::get_uniq_nb
auto get_uniq_nb(double *ptr)
Definition: HenckyOps.hpp:15
HenckyOps::sort_eigen_vals
auto sort_eigen_vals(FTensor::Tensor1< double, DIM > &eig, FTensor::Tensor2< double, DIM, DIM > &eigen_vec)
Definition: HenckyOps.hpp:23
HenckyOps::f
auto f
Definition: HenckyOps.hpp:5
HenckyOps::eps
constexpr double eps
Definition: HenckyOps.hpp:3
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:7
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:6
HenckyOps::is_eq
auto is_eq(const double &a, const double &b)
Definition: HenckyOps.hpp:11
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::Types::MatrixDouble3by3
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
MoFEM::determinantTensor2by2
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:1353
MoFEM::determinantTensor3by3
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:1342
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:21
scale
double scale
Definition: plastic.cpp:94
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:17
HenckyOps::CommonData::getMatFirstPiolaStress
auto getMatFirstPiolaStress()
Definition: HenckyOps.hpp:76
HenckyOps::CommonData::matEigVal
MatrixDouble matEigVal
Definition: HenckyOps.hpp:67
HenckyOps::CommonData::matLogCdC
MatrixDouble matLogCdC
Definition: HenckyOps.hpp:70
HenckyOps::CommonData::matEigVec
MatrixDouble matEigVec
Definition: HenckyOps.hpp:68
HenckyOps::CommonData::matTangent
MatrixDouble matTangent
Definition: HenckyOps.hpp:74
HenckyOps::CommonData::matLogCPlastic
boost::shared_ptr< MatrixDouble > matLogCPlastic
Definition: HenckyOps.hpp:65
HenckyOps::CommonData::matFirstPiolaStress
MatrixDouble matFirstPiolaStress
Definition: HenckyOps.hpp:71
HenckyOps::CommonData::matLogC
MatrixDouble matLogC
Definition: HenckyOps.hpp:69
HenckyOps::CommonData::matDPtr
boost::shared_ptr< MatrixDouble > matDPtr
Definition: HenckyOps.hpp:64
HenckyOps::CommonData::matSecondPiolaStress
MatrixDouble matSecondPiolaStress
Definition: HenckyOps.hpp:72
HenckyOps::CommonData::matHenckyStress
MatrixDouble matHenckyStress
Definition: HenckyOps.hpp:73
HenckyOps::CommonData::matGradPtr
boost::shared_ptr< MatrixDouble > matGradPtr
Definition: HenckyOps.hpp:63
HenckyOps::OpCalculateEigenVals::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:155
HenckyOps::OpCalculateEigenVals::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:104
HenckyOps::OpCalculateEigenVals::OpCalculateEigenVals
OpCalculateEigenVals(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:97
HenckyOps::OpCalculateHenckyPlasticStress::OpCalculateHenckyPlasticStress
OpCalculateHenckyPlasticStress(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr, const double scale=1)
Definition: HenckyOps.hpp:302
HenckyOps::OpCalculateHenckyPlasticStress::matDPtr
boost::shared_ptr< MatrixDouble > matDPtr
Definition: HenckyOps.hpp:346
HenckyOps::OpCalculateHenckyPlasticStress::matLogCPlastic
boost::shared_ptr< MatrixDouble > matLogCPlastic
Definition: HenckyOps.hpp:347
HenckyOps::OpCalculateHenckyPlasticStress::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:345
HenckyOps::OpCalculateHenckyPlasticStress::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:313
HenckyOps::OpCalculateHenckyStress::OpCalculateHenckyStress
OpCalculateHenckyStress(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:261
HenckyOps::OpCalculateHenckyStress::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:297
HenckyOps::OpCalculateHenckyStress::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:268
HenckyOps::OpCalculateLogC_dC::OpCalculateLogC_dC
OpCalculateLogC_dC(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:210
HenckyOps::OpCalculateLogC_dC::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:217
HenckyOps::OpCalculateLogC_dC::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:256
HenckyOps::OpCalculateLogC::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:205
HenckyOps::OpCalculateLogC::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:167
HenckyOps::OpCalculateLogC::OpCalculateLogC
OpCalculateLogC(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:160
HenckyOps::OpCalculatePiolaStress::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:360
HenckyOps::OpCalculatePiolaStress::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:403
HenckyOps::OpCalculatePiolaStress::OpCalculatePiolaStress
OpCalculatePiolaStress(const std::string field_name, boost::shared_ptr< CommonData > common_data)
Definition: HenckyOps.hpp:353
HenckyOps::OpHenckyTangent::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:494
HenckyOps::OpHenckyTangent::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: HenckyOps.hpp:419
HenckyOps::OpHenckyTangent::matDPtr
boost::shared_ptr< MatrixDouble > matDPtr
Definition: HenckyOps.hpp:495
HenckyOps::OpHenckyTangent::OpHenckyTangent
OpHenckyTangent(const std::string field_name, boost::shared_ptr< CommonData > common_data, boost::shared_ptr< MatrixDouble > mat_D_ptr=nullptr)
Definition: HenckyOps.hpp:407
HenckyOps::OpPostProcHencky::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
Definition: HenckyOps.hpp:524
HenckyOps::OpPostProcHencky::OpPostProcHencky
OpPostProcHencky(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
Definition: HenckyOps.hpp:512
HenckyOps::OpPostProcHencky::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: HenckyOps.hpp:507
HenckyOps::OpPostProcHencky::postProcMesh
moab::Interface & postProcMesh
Definition: HenckyOps.hpp:506
HenckyOps::OpPostProcHencky::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: HenckyOps.hpp:508
MoFEM::DataOperator::doEntities
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Definition: DataOperators.hpp:85
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:564
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1247
Generated by Doxygen 1.9.5 and hosted at University of Glasgow