23#ifndef EXECUTABLE_DIMENSION
24#define EXECUTABLE_DIMENSION 3
67 GAUSS>::OpSource<1, SPACE_DIM>;
72 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
86 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
198 auto get_command_line_parameters = [&]() {
221 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Hardening " <<
H;
222 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Viscous hardening " <<
visH;
223 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Saturation yield stress " <<
Qinf;
224 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Saturation exponent " <<
b_iso;
227 PetscBool is_scale = PETSC_TRUE;
244 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Scaled Hardening " <<
H;
245 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Scaled Viscous hardening " <<
visH;
247 <<
"Scaled Saturation yield stress " <<
Qinf;
253 auto set_matrial_stiffness = [&]() {
269 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
271 auto t_D_axiator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
273 auto t_D_deviator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
276 constexpr double third = boost::math::constants::third<double>();
277 t_D_axiator(
i,
j,
k,
l) =
A *
280 t_D_deviator(
i,
j,
k,
l) =
282 t_D(
i,
j,
k,
l) = t_D_axiator(
i,
j,
k,
l) + t_D_deviator(
i,
j,
k,
l);
300 CHKERR get_command_line_parameters();
301 CHKERR set_matrial_stiffness();
326 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
328 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
330 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
332 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
333 "REMOVE_ALL",
"U", 0, 3);
335 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_X",
"U",
337 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_Y",
"U",
339 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_Z",
"U",
341 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_ALL",
344 auto &bc_map = bc_mng->getBcMapByBlockName();
345 boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{
"FIX_"});
347 CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"REACTION",
350 for (
auto bc : bc_map)
351 MOFEM_LOG(
"EXAMPLE", Sev::verbose) <<
"Marker " << bc.first;
356 std::string reaction_block_set;
357 for (
auto bc : bc_map) {
358 if (bc_mng->checkBlock(bc,
"REACTION")) {
359 reaction_block_set = bc.first;
364 if (
auto bc = bc_mng->popMarkDOFsOnEntities(reaction_block_set)) {
371 ProblemsManager::MarkOP::AND, nodes,
375 MOFEM_LOG(
"EXAMPLE", Sev::warning) <<
"REACTION blockset does not exist";
379 MOFEM_LOG(
"EXAMPLE", Sev::warning) <<
"REACTION blockset does not exist";
401 auto integration_rule_deviator = [](
int o_row,
int o_col,
int approx_order) {
408 auto add_domain_base_ops = [&](
auto &pipeline) {
411 auto det_ptr = boost::make_shared<VectorDouble>();
412 auto jac_ptr = boost::make_shared<MatrixDouble>();
413 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
435 auto add_domain_stress_ops = [&](
auto &pipeline,
auto m_D_ptr) {
442 "Wrong pointer for grad");
470 auto add_domain_ops_lhs_mechanical = [&](
auto &pipeline,
auto m_D_ptr) {
482 pipeline.push_back(
new OpKCauchy(
"U",
"U", m_D_ptr));
491 auto add_domain_ops_rhs_mechanical = [&](
auto &pipeline) {
496 const std::string block_name =
"BODY_FORCE";
497 if (it->getName().compare(0, block_name.size(), block_name) == 0) {
498 std::vector<double> attr;
499 CHKERR it->getAttributes(attr);
500 if (attr.size() == 3) {
505 "Should be three atributes in BODYFORCE blockset, but is %d",
516 auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
517 const auto time = fe_domain_rhs->ts_t;
520 t_source(
i) += (
scale * t_b(
i)) * time;
524 pipeline.push_back(
new OpBodyForce(
"U", get_body_force));
539 auto add_domain_ops_lhs_constrain = [&](
auto &pipeline,
auto m_D_ptr) {
570 auto add_domain_ops_rhs_constrain = [&](
auto &pipeline) {
582 auto add_boundary_ops_lhs_mechanical = [&](
auto &pipeline) {
585 for (
auto bc : bc_map) {
586 if (bc_mng->checkBlock(bc,
"FIX_")){
588 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
590 "U",
"U", [](
double,
double,
double) {
return 1.; },
591 bc.second->getBcEntsPtr()));
598 auto add_boundary_ops_rhs_mechanical = [&](
auto &pipeline) {
604 return fe_domain_rhs->ts_t;
610 return fe_domain_rhs->ts_t *
scale;
616 return -fe_domain_rhs->ts_t;
622 return -fe_domain_rhs->ts_t;
628 if (it->getName().compare(0, 5,
"FORCE") == 0) {
630 std::vector<double> attr_vec;
631 CHKERR it->getMeshsetIdEntitiesByDimension(
633 it->getAttributes(attr_vec);
636 "Wrong size of boundary attributes vector. Set right block "
638 auto force_vec_ptr = boost::make_shared<MatrixDouble>(
SPACE_DIM, 1);
639 std::copy(&attr_vec[0], &attr_vec[
SPACE_DIM],
640 force_vec_ptr->data().begin());
643 boost::make_shared<Range>(force_edges)));
649 auto u_mat_ptr = boost::make_shared<MatrixDouble>();
654 if (bc_mng->checkBlock(bc,
"FIX_")) {
656 new OpSetBc(
"U",
false, bc.second->getBcMarkersPtr()));
657 auto attr_vec = boost::make_shared<MatrixDouble>(
SPACE_DIM, 1);
659 if (bc.second->bcAttributes.size() <
SPACE_DIM)
661 "Wrong size of boundary attributes vector. Set right block "
662 "size attributes. Size of attributes %d",
663 bc.second->bcAttributes.size());
664 std::copy(&bc.second->bcAttributes[0],
666 attr_vec->data().begin());
668 pipeline.push_back(
new OpBoundaryVec(
"U", attr_vec, time_scaled,
669 bc.second->getBcEntsPtr()));
671 "U", u_mat_ptr, [](
double,
double,
double) {
return 1.; },
672 bc.second->getBcEntsPtr()));
682 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
683 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
685 CHKERR add_domain_ops_lhs_mechanical(pipeline_mng->getOpDomainLhsPipeline(),
687 CHKERR add_domain_ops_lhs_constrain(pipeline_mng->getOpDomainLhsPipeline(),
689 CHKERR add_boundary_ops_lhs_mechanical(
690 pipeline_mng->getOpBoundaryLhsPipeline());
698 CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
699 CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
701 CHKERR add_domain_ops_rhs_mechanical(pipeline_mng->getOpDomainRhsPipeline());
702 CHKERR add_domain_ops_rhs_constrain(pipeline_mng->getOpDomainRhsPipeline());
703 CHKERR add_boundary_ops_rhs_mechanical(
704 pipeline_mng->getOpBoundaryRhsPipeline());
711 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule_deviator);
712 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule_deviator);
714 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
715 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
717 auto create_reaction_pipeline = [&](
auto &pipeline) {
722 auto det_ptr = boost::make_shared<VectorDouble>();
723 auto jac_ptr = boost::make_shared<MatrixDouble>();
724 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
742 "Wrong pointer for grad");
779 reactionFe->getRuleHook = integration_rule_deviator;
797 auto set_section_monitor = [&](
auto solver) {
800 CHKERR TSGetSNES(solver, &snes);
801 CHKERR SNESMonitorSet(snes,
804 (
void *)(snes_ctx_ptr.get()),
nullptr);
808 auto create_post_process_element = [&]() {
813 auto det_ptr = boost::make_shared<VectorDouble>();
814 auto jac_ptr = boost::make_shared<MatrixDouble>();
815 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
836 "Wrong pointer for grad");
875 auto scatter_create = [&](
auto D,
auto coeff) {
877 CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
878 ROW,
"U", coeff, coeff, is);
880 CHKERR ISGetLocalSize(is, &loc_size);
884 CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
889 auto set_time_monitor = [&](
auto dm,
auto solver) {
891 boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(
893 boost::shared_ptr<ForcesAndSourcesCore> null;
895 monitor_ptr, null, null);
899 auto set_fieldsplit_preconditioner = [&](
auto solver) {
903 CHKERR TSGetSNES(solver, &snes);
905 CHKERR SNESGetKSP(snes, &ksp);
907 CHKERR KSPGetPC(ksp, &pc);
908 PetscBool is_pcfs = PETSC_FALSE;
909 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
912 if (is_pcfs == PETSC_TRUE) {
915 auto name_prb =
simple->getProblemName();
916 auto is_all_bc = bc_mng->getBlockIS(name_prb,
"FIX_X",
"U", 0, 0);
917 is_all_bc = bc_mng->getBlockIS(name_prb,
"FIX_Y",
"U", 1, 1, is_all_bc);
918 is_all_bc = bc_mng->getBlockIS(name_prb,
"FIX_Z",
"U", 2, 2, is_all_bc);
919 is_all_bc = bc_mng->getBlockIS(name_prb,
"FIX_ALL",
"U", 0, 2, is_all_bc);
922 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
924 <<
"Field split block size " << is_all_bc_size;
926 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
933 auto dm =
simple->getDM();
936 boost::shared_ptr<FEMethod> null;
942 CHKERR create_post_process_element();
948 auto solver = pipeline_mng->createTSIM();
950 CHKERR TSSetSolution(solver,
D);
951 CHKERR set_section_monitor(solver);
952 CHKERR set_time_monitor(dm, solver);
953 CHKERR TSSetSolution(solver,
D);
954 CHKERR TSSetFromOptions(solver);
955 CHKERR set_fieldsplit_preconditioner(solver);
957 CHKERR TSSolve(solver, NULL);
959 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
960 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
969int main(
int argc,
char *argv[]) {
976 auto core_log = logging::core::get();
978 LogManager::createSink(LogManager::getStrmWorld(),
"EXAMPLE"));
979 LogManager::setLog(
"EXAMPLE");
985 DMType dm_name =
"DMMOFEM";
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
void temp(int x, int y=10)
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, 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 TS Jacobian evaluation function
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
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.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
auto smartGetDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
const double D
diffusivity
int main(int argc, char *argv[])
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
EntitiesFieldData::EntData EntData
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
static char help[]
[Solve]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
#define EXECUTABLE_DIMENSION
long double hardening(long double tau, double temp)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
long double hardening_dtau(long double tau, double temp)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
PetscBool is_large_strains
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
[Body force]
constexpr EntityType boundary_ent
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
boost::shared_ptr< DomainEle > feAxiatorRhs
boost::shared_ptr< DomainEle > reactionFe
boost::shared_ptr< std::vector< unsigned char > > reactionMarker
boost::shared_ptr< PostProcEle > postProcFe
boost::shared_ptr< DomainEle > feAxiatorLhs
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
FieldApproximationBase base
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
MoFEMErrorCode createCommonData()
[Set up problem]
std::vector< FTensor::Tensor1< double, 3 > > bodyForces
Example(MoFEM::Interface &m_field)
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode bC()
[Create common data]
boost::shared_ptr< PlasticOps::CommonData > commonPlasticDataPtr
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Simple interface for fast problem set-up.
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
default operator for EDGE element
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
Section manager is used to create indexes and sections.
Get value at integration points for scalar field.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Scale base functions by inverses of measure of element.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.