#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
};
};
GAUSS>::OpSource<1, SPACE_DIM>;
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
}
}
private:
boost::shared_ptr<PlasticOps::CommonData> commonPlasticDataPtr;
boost::shared_ptr<HenckyOps::CommonData> commonHenckyDataPtr;
boost::shared_ptr<PostProcEle> postProcFe;
boost::shared_ptr<DomainEle> reactionFe;
boost::shared_ptr<DomainEle> feAxiatorRhs;
boost::shared_ptr<DomainEle> feAxiatorLhs;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
boost::shared_ptr<std::vector<unsigned char>> reactionMarker;
std::vector<FTensor::Tensor1<double, 3>> bodyForces;
};
}
Simple *
simple = mField.getInterface<Simple>();
}
auto get_command_line_parameters = [&]() {
PETSC_NULL);
PETSC_NULL);
MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Hardening " <<
H;
MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Viscous hardening " <<
visH;
MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Saturation yield stress " <<
Qinf;
PetscBool is_scale = PETSC_TRUE;
PETSC_NULL);
if (is_scale) {
MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Scaled Hardening " <<
H;
MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Scaled Viscous hardening " <<
visH;
<<
"Scaled Saturation yield stress " <<
Qinf;
}
};
auto set_matrial_stiffness = [&]() {
: 1;
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
*commonPlasticDataPtr->mDPtr);
auto t_D_axiator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
*commonPlasticDataPtr->mDPtr_Axiator);
auto t_D_deviator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
*commonPlasticDataPtr->mDPtr_Deviator);
constexpr
double third = boost::math::constants::third<double>();
t_D_axiator(
i,
j,
k,
l) =
A *
t_D_deviator(
i,
j,
k,
l) =
t_D(
i,
j,
k,
l) = t_D_axiator(
i,
j,
k,
l) + t_D_deviator(
i,
j,
k,
l);
};
commonPlasticDataPtr = boost::make_shared<PlasticOps::CommonData>();
commonPlasticDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mDPtr->resize(size_symm * size_symm, 1);
commonPlasticDataPtr->mDPtr_Axiator = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mDPtr_Axiator->resize(size_symm * size_symm, 1);
commonPlasticDataPtr->mDPtr_Deviator = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mDPtr_Deviator->resize(size_symm * size_symm, 1);
commonPlasticDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
commonPlasticDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
CHKERR get_command_line_parameters();
CHKERR set_matrial_stiffness();
commonHenckyDataPtr = boost::make_shared<HenckyOps::CommonData>();
commonHenckyDataPtr->matGradPtr = commonPlasticDataPtr->mGradPtr;
commonHenckyDataPtr->matDPtr = commonPlasticDataPtr->mDPtr;
commonHenckyDataPtr->matLogCPlastic =
commonPlasticDataPtr->getPlasticStrainPtr();
commonPlasticDataPtr->mStrainPtr = commonHenckyDataPtr->getMatLogC();
commonPlasticDataPtr->mStressPtr =
commonHenckyDataPtr->getMatHenckyStress();
}
}
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
auto prb_mng = mField.getInterface<ProblemsManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_X",
"U",
0, 0);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_Y",
"U",
1, 1);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_Z",
"U",
2, 2);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_ALL",
"U", 0, 3);
auto &bc_map = bc_mng->getBcMapByBlockName();
boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"REACTION",
"U", 0, 3);
for (auto bc : bc_map)
MOFEM_LOG(
"EXAMPLE", Sev::verbose) <<
"Marker " << bc.first;
std::string reaction_block_set;
for (auto bc : bc_map) {
if (bc_mng->checkBlock(bc, "REACTION")) {
reaction_block_set = bc.first;
break;
}
}
if (auto bc = bc_mng->popMarkDOFsOnEntities(reaction_block_set)) {
reactionMarker = bc->getBcMarkersPtr();
Range nodes;
CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes,
true);
ProblemsManager::MarkOP::AND, nodes,
*reactionMarker);
} else {
MOFEM_LOG(
"EXAMPLE", Sev::warning) <<
"REACTION blockset does not exist";
}
if (!reactionMarker) {
MOFEM_LOG(
"EXAMPLE", Sev::warning) <<
"REACTION blockset does not exist";
}
}
auto pipeline_mng = mField.getInterface<PipelineManager>();
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
feAxiatorLhs = boost::make_shared<DomainEle>(mField);
feAxiatorRhs = boost::make_shared<DomainEle>(mField);
};
feAxiatorLhs->getRuleHook = integration_rule_axiator;
feAxiatorRhs->getRuleHook = integration_rule_axiator;
auto integration_rule_deviator = [](
int o_row,
int o_col,
int approx_order) {
};
};
auto add_domain_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
}
"TAU", commonPlasticDataPtr->getPlasticTauDotPtr()));
pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
pipeline.push_back(new OpCalculateTensor2SymmetricFieldValuesDot<SPACE_DIM>(
"EP", commonPlasticDataPtr->getPlasticStrainDotPtr()));
"U", commonPlasticDataPtr->mGradPtr));
pipeline.push_back(new OpCalculateScalarFieldValues(
"TAU", commonPlasticDataPtr->getPlasticTauPtr()));
};
auto add_domain_stress_ops = [&](auto &pipeline, auto m_D_ptr) {
if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
"Wrong pointer for grad");
pipeline.push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(
new OpCalculateLogC<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
"U", commonHenckyDataPtr, m_D_ptr));
pipeline.push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", commonHenckyDataPtr));
} else {
pipeline.push_back(
new OpSymmetrizeTensor<SPACE_DIM>("U", commonPlasticDataPtr->mGradPtr,
commonPlasticDataPtr->mStrainPtr));
pipeline.push_back(
new OpPlasticStress("U", commonPlasticDataPtr, m_D_ptr, 1));
}
if (m_D_ptr != commonPlasticDataPtr->mDPtr_Axiator)
pipeline.push_back(
new OpCalculatePlasticSurface("U", commonPlasticDataPtr));
};
auto add_domain_ops_lhs_mechanical = [&](auto &pipeline, auto m_D_ptr) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
new OpHenckyTangent<SPACE_DIM>("U", commonHenckyDataPtr, m_D_ptr));
pipeline.push_back(
new OpKPiola(
"U",
"U", commonHenckyDataPtr->getMatTangent()));
pipeline.push_back(new OpCalculatePlasticInternalForceLhs_LogStrain_dEP(
"U", "EP", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
} else {
pipeline.push_back(
new OpKCauchy(
"U",
"U", m_D_ptr));
pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP(
"U", "EP", commonPlasticDataPtr, m_D_ptr));
}
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_ops_rhs_mechanical = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
const std::string block_name = "BODY_FORCE";
if (it->getName().compare(0, block_name.size(), block_name) == 0) {
std::vector<double> attr;
CHKERR it->getAttributes(attr);
if (attr.size() == 3) {
bodyForces.push_back(
} else {
"Should be three atributes in BODYFORCE blockset, but is %d",
attr.size());
}
}
}
auto get_body_force = [this](const double, const double, const double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
const auto time = fe_domain_rhs->ts_t;
for (auto &t_b : bodyForces)
t_source(
i) += (
scale * t_b(
i)) * time;
return t_source;
};
pipeline.push_back(
new OpBodyForce(
"U", get_body_force));
"U", commonHenckyDataPtr->getMatFirstPiolaStress()));
} else {
pipeline.push_back(
}
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_ops_lhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
new OpHenckyTangent<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(new OpCalculateContrainsLhs_LogStrain_dU(
"TAU", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
pipeline.push_back(new OpCalculatePlasticFlowLhs_LogStrain_dU(
"EP", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
} else {
pipeline.push_back(new OpCalculatePlasticFlowLhs_dU(
"EP", "U", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(new OpCalculateContrainsLhs_dU(
"TAU", "U", commonPlasticDataPtr, m_D_ptr));
}
pipeline.push_back(new OpCalculatePlasticFlowLhs_dEP(
"EP", "EP", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(
new OpCalculatePlasticFlowLhs_dTAU("EP", "TAU", commonPlasticDataPtr));
pipeline.push_back(new OpCalculateContrainsLhs_dEP(
"TAU", "EP", commonPlasticDataPtr, m_D_ptr));
pipeline.push_back(
new OpCalculateContrainsLhs_dTAU("TAU", "TAU", commonPlasticDataPtr));
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_ops_rhs_constrain = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
new OpCalculatePlasticFlowRhs("EP", commonPlasticDataPtr));
pipeline.push_back(
new OpCalculateContrainsRhs("TAU", commonPlasticDataPtr));
};
auto add_boundary_ops_lhs_mechanical = [&](auto &pipeline) {
auto &bc_map = mField.getInterface<BcManager>()->getBcMapByBlockName();
for (auto bc : bc_map) {
if (bc_mng->checkBlock(bc, "FIX_")){
pipeline.push_back(
new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
"U", "U", [](double, double, double) { return 1.; },
bc.second->getBcEntsPtr()));
pipeline.push_back(new OpUnSetBc("U"));
}
}
};
auto add_boundary_ops_rhs_mechanical = [&](auto &pipeline) {
auto get_time = [&](double, double, double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
return fe_domain_rhs->ts_t;
};
auto get_time_scaled = [&](double, double, double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
return fe_domain_rhs->ts_t *
scale;
};
auto get_minus_time = [&](double, double, double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
return -fe_domain_rhs->ts_t;
};
auto time_scaled = [&](double, double, double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_rhs = pipeline_mng->getBoundaryRhsFE();
return -fe_domain_rhs->ts_t;
};
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
if (it->getName().compare(0, 5, "FORCE") == 0) {
Range force_edges;
std::vector<double> attr_vec;
CHKERR it->getMeshsetIdEntitiesByDimension(
mField.get_moab(),
SPACE_DIM - 1, force_edges,
true);
it->getAttributes(attr_vec);
"Wrong size of boundary attributes vector. Set right block "
"size attributes.");
auto force_vec_ptr = boost::make_shared<MatrixDouble>(
SPACE_DIM, 1);
std::copy(&attr_vec[0], &attr_vec[
SPACE_DIM],
force_vec_ptr->data().begin());
pipeline.push_back(
boost::make_shared<Range>(force_edges)));
}
}
pipeline.push_back(new OpUnSetBc("U"));
auto u_mat_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_mat_ptr));
for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
if (bc_mng->checkBlock(bc, "FIX_")) {
pipeline.push_back(
new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
auto attr_vec = boost::make_shared<MatrixDouble>(
SPACE_DIM, 1);
attr_vec->clear();
if (bc.second->bcAttributes.size() <
SPACE_DIM)
"Wrong size of boundary attributes vector. Set right block "
"size attributes. Size of attributes %d",
bc.second->bcAttributes.size());
std::copy(&bc.second->bcAttributes[0],
attr_vec->data().begin());
bc.second->getBcEntsPtr()));
"U", u_mat_ptr, [](double, double, double) { return 1.; },
bc.second->getBcEntsPtr()));
pipeline.push_back(new OpUnSetBc("U"));
}
}
};
CHKERR add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainLhsPipeline(),
commonPlasticDataPtr->mDPtr_Deviator);
CHKERR add_domain_ops_lhs_mechanical(pipeline_mng->getOpDomainLhsPipeline(),
commonPlasticDataPtr->mDPtr_Deviator);
CHKERR add_domain_ops_lhs_constrain(pipeline_mng->getOpDomainLhsPipeline(),
commonPlasticDataPtr->mDPtr_Deviator);
CHKERR add_boundary_ops_lhs_mechanical(
pipeline_mng->getOpBoundaryLhsPipeline());
CHKERR add_domain_base_ops(feAxiatorLhs->getOpPtrVector());
CHKERR add_domain_stress_ops(feAxiatorLhs->getOpPtrVector(),
commonPlasticDataPtr->mDPtr_Axiator);
CHKERR add_domain_ops_lhs_mechanical(feAxiatorLhs->getOpPtrVector(),
commonPlasticDataPtr->mDPtr_Axiator);
CHKERR add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_domain_stress_ops(pipeline_mng->getOpDomainRhsPipeline(),
commonPlasticDataPtr->mDPtr_Deviator);
CHKERR add_domain_ops_rhs_mechanical(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_domain_ops_rhs_constrain(pipeline_mng->getOpDomainRhsPipeline());
CHKERR add_boundary_ops_rhs_mechanical(
pipeline_mng->getOpBoundaryRhsPipeline());
CHKERR add_domain_base_ops(feAxiatorRhs->getOpPtrVector());
CHKERR add_domain_stress_ops(feAxiatorRhs->getOpPtrVector(),
commonPlasticDataPtr->mDPtr_Axiator);
CHKERR add_domain_ops_rhs_mechanical(feAxiatorRhs->getOpPtrVector());
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule_deviator);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule_deviator);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
auto create_reaction_pipeline = [&](auto &pipeline) {
if (reactionMarker) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
}
pipeline.push_back(
"U", commonPlasticDataPtr->mGradPtr));
pipeline.push_back(new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
"Wrong pointer for grad");
pipeline.push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(
new OpCalculateLogC<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", commonHenckyDataPtr));
pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
"U", commonHenckyDataPtr, commonPlasticDataPtr->mDPtr));
pipeline.push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", commonHenckyDataPtr));
} else {
pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
"U", commonPlasticDataPtr->mGradPtr,
commonPlasticDataPtr->mStrainPtr));
pipeline.push_back(new OpPlasticStress("U", commonPlasticDataPtr,
commonPlasticDataPtr->mDPtr, 1));
}
pipeline.push_back(new OpSetBc("U", false, reactionMarker));
"U", commonHenckyDataPtr->getMatFirstPiolaStress()));
} else {
pipeline.push_back(
}
pipeline.push_back(new OpUnSetBc("U"));
}
};
reactionFe = boost::make_shared<DomainEle>(mField);
reactionFe->getRuleHook = integration_rule_deviator;
CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
}
Simple *
simple = mField.getInterface<Simple>();
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
ISManager *is_manager = mField.getInterface<ISManager>();
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
snes,
(
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
};
auto create_post_process_element = [&]() {
postProcFe = boost::make_shared<PostProcEle>(mField);
postProcFe->generateReferenceElementMesh();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
postProcFe->getOpPtrVector().push_back(
postProcFe->getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
postProcFe->getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
}
postProcFe->getOpPtrVector().push_back(
"U", commonPlasticDataPtr->mGradPtr));
postProcFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
"TAU", commonPlasticDataPtr->getPlasticTauPtr()));
postProcFe->getOpPtrVector().push_back(
new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
"EP", commonPlasticDataPtr->getPlasticStrainPtr()));
if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
"Wrong pointer for grad");
postProcFe->getOpPtrVector().push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", commonHenckyDataPtr));
postProcFe->getOpPtrVector().push_back(
new OpCalculateLogC<SPACE_DIM>("U", commonHenckyDataPtr));
postProcFe->getOpPtrVector().push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", commonHenckyDataPtr));
postProcFe->getOpPtrVector().push_back(
new OpCalculateHenckyPlasticStress<SPACE_DIM>(
"U", commonHenckyDataPtr, commonPlasticDataPtr->mDPtr,
scale));
postProcFe->getOpPtrVector().push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", commonHenckyDataPtr));
postProcFe->getOpPtrVector().push_back(new OpPostProcHencky<SPACE_DIM>(
"U", postProcFe->postProcMesh, postProcFe->mapGaussPts,
commonHenckyDataPtr));
} else {
postProcFe->getOpPtrVector().push_back(
new OpSymmetrizeTensor<SPACE_DIM>("U", commonPlasticDataPtr->mGradPtr,
commonPlasticDataPtr->mStrainPtr));
postProcFe->getOpPtrVector().push_back(new OpPlasticStress(
"U", commonPlasticDataPtr, commonPlasticDataPtr->mDPtr,
scale));
postProcFe->getOpPtrVector().push_back(
"U", postProcFe->postProcMesh, postProcFe->mapGaussPts,
commonPlasticDataPtr->mStrainPtr,
commonPlasticDataPtr->mStressPtr));
}
postProcFe->getOpPtrVector().push_back(
new OpCalculatePlasticSurface("U", commonPlasticDataPtr));
postProcFe->getOpPtrVector().push_back(
new OpPostProcPlastic("U", postProcFe->postProcMesh,
postProcFe->mapGaussPts, commonPlasticDataPtr));
postProcFe->addFieldValuesPostProc("U");
};
auto scatter_create = [&](
auto D,
auto coeff) {
SmartPetscObj<IS> is;
CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
ROW,
"U", coeff, coeff, is);
int loc_size;
CHKERR ISGetLocalSize(is, &loc_size);
CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &
v);
VecScatter scatter;
CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
return std::make_tuple(SmartPetscObj<Vec>(
v),
SmartPetscObj<VecScatter>(scatter));
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(
dm, postProcFe, reactionFe, uXScatter, uYScatter, uZScatter));
boost::shared_ptr<ForcesAndSourcesCore> null;
monitor_ptr, null, null);
};
auto set_fieldsplit_preconditioner = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto bc_mng = mField.getInterface<BcManager>();
auto name_prb =
simple->getProblemName();
auto is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_X", "U", 0, 0);
is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Y", "U", 1, 1, is_all_bc);
is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Z", "U", 2, 2, is_all_bc);
is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_ALL", "U", 0, 2, is_all_bc);
int is_all_bc_size;
CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
<< "Field split block size " << is_all_bc_size;
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc);
}
};
boost::shared_ptr<FEMethod> null;
null, null);
null, null);
CHKERR create_post_process_element();
uXScatter = scatter_create(
D, 0);
uYScatter = scatter_create(
D, 1);
uZScatter = scatter_create(
D, 2);
auto solver = pipeline_mng->createTSIM();
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetFromOptions(solver);
CHKERR set_fieldsplit_preconditioner(solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
try {
DMType dm_name = "DMMOFEM";
}
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
#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 ...
void temp(int x, int y=10)
auto smartCreateDMVector
Get smart vector from DM.
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
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#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
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)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
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
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode bC()
[Create common data]
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Scale base functions by inverses of measure of element.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]