43struct CommonData :
public boost::enable_shared_from_this<CommonData> {
61 return boost::shared_ptr<BlockParams>(shared_from_this(), &
blockParams);
65 boost::shared_ptr<MatrixDouble>
mDPtr;
88 return boost::shared_ptr<VectorDouble>(shared_from_this(), &
plasticSurface);
91 return boost::shared_ptr<VectorDouble>(shared_from_this(), &
plasticTau);
94 return boost::shared_ptr<VectorDouble>(shared_from_this(), &
plasticTauDot);
97 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
plasticStrain);
100 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
104 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
plasticFlow);
120template <
int DIM, IntegrationType I,
typename DomainEleOp>
123template <
int DIM, IntegrationType I,
typename DomainEleOp>
126template <
int DIM, IntegrationType I,
typename DomainEleOp>
129template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
132template <IntegrationType I,
typename AssemblyDomainEleOp>
135template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
138template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
141template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
144template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
147template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
150template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
153template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
156template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
159template <
int DIM, IntegrationType I,
typename AssemblyDomainEleOp>
162template <IntegrationType I,
typename AssemblyDomainEleOp>
166 template <
int DIM, IntegrationType I>
170 template <
int DIM, IntegrationType I>
173 template <
int DIM, IntegrationType I>
181 template <
int DIM, IntegrationType I>
185 template <IntegrationType I>
189 template <
int DIM, IntegrationType I>
193 template <
int DIM, IntegrationType I>
198 template <
int DIM, IntegrationType I>
202 template <
int DIM, IntegrationType I>
206 template <
int DIM, IntegrationType I>
210 template <
int DIM, IntegrationType I>
214 template <
int DIM, IntegrationType I>
218 template <
int DIM, IntegrationType I>
222 template <
int DIM, IntegrationType I>
226 template <IntegrationType I>
242using Pip = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>;
249 boost::shared_ptr<MatrixDouble> mat_D_Ptr,
250 boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
251 double scale_value, Sev sev) {
255 OpMatBlocks(boost::shared_ptr<MatrixDouble> m_D_ptr,
256 boost::shared_ptr<CommonData::BlockParams> mat_params_ptr,
258 std::vector<const CubitMeshSets *> meshset_vec_ptr)
260 matParamsPtr(mat_params_ptr), scaleVal(scale_value) {
262 "Can not get data from block");
265 MoFEMErrorCode doWork(
int side,
EntityType type,
266 EntitiesFieldData::EntData &data) {
269 auto getK = [](
auto &
p) {
275 auto getG = [](
auto &
p) {
281 auto scale_fun = [
this](
auto &
p) {
282 p[CommonData::YOUNG_MODULUS] *= scaleVal;
283 p[CommonData::SIGMA_Y] *= scaleVal;
284 p[CommonData::H] *= scaleVal;
285 p[CommonData::VIS_H] *= scaleVal;
286 p[CommonData::QINF] *= scaleVal;
287 p[CommonData::C1_k] *= scaleVal;
290 for (
auto &b : blockData) {
291 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
292 *matParamsPtr = b.bParams;
293 scale_fun(*matParamsPtr);
294 CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
301 scale_fun(*matParamsPtr);
302 CHKERR getMatDPtr(matDPtr, getK(*matParamsPtr), getG(*matParamsPtr));
308 boost::shared_ptr<MatrixDouble> matDPtr;
309 boost::shared_ptr<CommonData::BlockParams> matParamsPtr;
310 const double scaleVal;
313 std::array<double, CommonData::LAST_PARAM> bParams;
316 std::vector<BlockData> blockData;
328 std::vector<const CubitMeshSets *> meshset_vec_ptr,
332 for (
auto m : meshset_vec_ptr) {
334 std::vector<double> block_data;
335 CHKERR m->getAttributes(block_data);
336 if (block_data.size() != CommonData::LAST_PARAM) {
338 "Wrong number of block attribute");
340 auto get_block_ents = [&]() {
348 for (
auto i = 0;
i != CommonData::LAST_PARAM; ++
i) {
349 block_params[
i] = block_data[
i];
353 <<
"E = " << block_params[CommonData::YOUNG_MODULUS]
354 <<
" nu = " << block_params[CommonData::POISSON_RATIO];
357 <<
"sigma_y = " << block_params[CommonData::SIGMA_Y] << std::endl
358 <<
"h = " << block_params[CommonData::H] << std::endl
359 <<
"vis_h = " << block_params[CommonData::VIS_H] << std::endl
360 <<
"qinf = " << block_params[CommonData::QINF] << std::endl
361 <<
"biso = " << block_params[CommonData::BISO] << std::endl
362 <<
"C1_k = " << block_params[CommonData::C1_k] << std::endl;
364 blockData.push_back({block_params, get_block_ents()});
381 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
385 auto set_material_stiffness = [&]() {
391 double A = (DIM == 2)
395 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*mat_D_ptr);
402 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
404 set_material_stiffness();
410 pip.push_back(
new OpMatBlocks(
411 mat_D_Ptr, mat_params_ptr, scale_value, m_field, sev,
414 m_field.
getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
416 (boost::format(
"%s(.*)") % block_name).str()
425template <
int DIM, IntegrationType I,
typename DomainEleOp>
428 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
429 std::string u, std::string ep, std::string tau,
double scale, Sev sev) {
433 auto common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
434 common_plastic_ptr = boost::make_shared<PlasticOps::CommonData>();
436 constexpr auto size_symm = (DIM * (DIM + 1)) / 2;
437 auto make_d_mat = []() {
441 common_plastic_ptr->mDPtr = make_d_mat();
442 common_plastic_ptr->mGradPtr = boost::make_shared<MatrixDouble>();
443 common_plastic_ptr->mStrainPtr = boost::make_shared<MatrixDouble>();
444 common_plastic_ptr->mStressPtr = boost::make_shared<MatrixDouble>();
446 auto m_D_ptr = common_plastic_ptr->mDPtr;
449 common_plastic_ptr->getParamsPtr(),
451 "add mat block ops");
453 pip.push_back(
new OpCalculateScalarFieldValues(
454 tau, common_plastic_ptr->getPlasticTauPtr()));
455 pip.push_back(
new OpCalculateTensor2SymmetricFieldValues<DIM>(
456 ep, common_plastic_ptr->getPlasticStrainPtr()));
458 u, common_plastic_ptr->mGradPtr));
463 common_henky_ptr = boost::make_shared<HenckyOps::CommonData>();
464 common_henky_ptr->matGradPtr = common_plastic_ptr->mGradPtr;
465 common_henky_ptr->matDPtr = common_plastic_ptr->mDPtr;
466 common_henky_ptr->matLogCPlastic =
467 common_plastic_ptr->getPlasticStrainPtr();
468 common_plastic_ptr->mStrainPtr = common_henky_ptr->getMatLogC();
469 common_plastic_ptr->mStressPtr = common_henky_ptr->getMatHenckyStress();
473 pip.push_back(
new typename H::template OpCalculateEigenVals<DIM, I>(
474 u, common_henky_ptr));
476 new typename H::template OpCalculateLogC<DIM, I>(u, common_henky_ptr));
477 pip.push_back(
new typename H::template OpCalculateLogC_dC<DIM, I>(
478 u, common_henky_ptr));
480 typename H::template OpCalculateHenckyPlasticStress<DIM, I>(
481 u, common_henky_ptr, m_D_ptr));
482 pip.push_back(
new typename H::template OpCalculatePiolaStress<DIM, I>(
483 u, common_henky_ptr));
486 pip.push_back(
new OpSymmetrizeTensor<SPACE_DIM>(
487 u, common_plastic_ptr->mGradPtr, common_plastic_ptr->mStrainPtr));
488 pip.push_back(
new typename P::template OpPlasticStress<DIM, I>(
489 u, common_plastic_ptr, m_D_ptr));
492 pip.push_back(
new typename P::template OpCalculatePlasticSurface<DIM, I>(
493 u, common_plastic_ptr));
495 return std::make_tuple(common_plastic_ptr, common_henky_ptr);
498template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
501 std::string u, std::string ep, std::string tau) {
504 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
507 typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
508 using OpInternalForcePiola =
513 auto [common_plastic_ptr, common_henky_ptr] =
514 createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
515 ep, tau,
scale, Sev::inform);
517 auto m_D_ptr = common_plastic_ptr->mDPtr;
519 pip.push_back(
new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
520 ep, common_plastic_ptr->getPlasticStrainDotPtr()));
521 pip.push_back(
new OpCalculateScalarFieldValuesDot(
522 tau, common_plastic_ptr->getPlasticTauDotPtr()));
523 pip.push_back(
new typename P::template OpCalculatePlasticity<DIM, I>(
524 u, common_plastic_ptr, m_D_ptr));
527 if (common_henky_ptr) {
528 pip.push_back(
new OpInternalForcePiola(
529 u, common_henky_ptr->getMatFirstPiolaStress()));
536 typename P::template Assembly<A>::template OpCalculateConstraintsRhs<I>(
537 tau, common_plastic_ptr, m_D_ptr));
539 new typename P::template Assembly<A>::template OpCalculatePlasticFlowRhs<
540 DIM,
I>(ep, common_plastic_ptr, m_D_ptr));
545template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
548 std::string u, std::string ep, std::string tau) {
553 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
555 using OpKPiola =
typename B::template OpGradTensorGrad<1, DIM, DIM, 1>;
556 using OpKCauchy =
typename B::template OpGradSymTensorGrad<1, DIM, DIM, 0>;
560 auto [common_plastic_ptr, common_henky_ptr] =
561 createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
562 ep, tau,
scale, Sev::verbose);
564 auto m_D_ptr = common_plastic_ptr->mDPtr;
566 pip.push_back(
new OpCalculateTensor2SymmetricFieldValuesDot<DIM>(
567 ep, common_plastic_ptr->getPlasticStrainDotPtr()));
568 pip.push_back(
new OpCalculateScalarFieldValuesDot(
569 tau, common_plastic_ptr->getPlasticTauDotPtr()));
570 pip.push_back(
new typename P::template OpCalculatePlasticity<DIM, I>(
571 u, common_plastic_ptr, m_D_ptr));
573 if (common_henky_ptr) {
575 pip.push_back(
new typename H::template OpHenckyTangent<DIM, I>(
576 u, common_henky_ptr, m_D_ptr));
577 pip.push_back(
new OpKPiola(u, u, common_henky_ptr->getMatTangent()));
579 new typename P::template Assembly<A>::
580 template OpCalculatePlasticInternalForceLhs_LogStrain_dEP<DIM, I>(
581 u, ep, common_plastic_ptr, common_henky_ptr, m_D_ptr));
583 pip.push_back(
new OpKCauchy(u, u, m_D_ptr));
584 pip.push_back(
new typename P::template Assembly<A>::
585 template OpCalculatePlasticInternalForceLhs_dEP<DIM, I>(
586 u, ep, common_plastic_ptr, m_D_ptr));
589 if (common_henky_ptr) {
591 new typename P::template Assembly<A>::
592 template OpCalculateConstraintsLhs_LogStrain_dU<DIM, I>(
593 tau, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
595 new typename P::template Assembly<A>::
596 template OpCalculatePlasticFlowLhs_LogStrain_dU<DIM, I>(
597 ep, u, common_plastic_ptr, common_henky_ptr, m_D_ptr));
601 typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dU<
602 DIM,
I>(tau, u, common_plastic_ptr, m_D_ptr));
605 typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dU<
606 DIM,
I>(ep, u, common_plastic_ptr, m_D_ptr));
611 typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dEP<
612 DIM,
I>(ep, ep, common_plastic_ptr, m_D_ptr));
615 typename P::template Assembly<A>::template OpCalculatePlasticFlowLhs_dTAU<
616 DIM,
I>(ep, tau, common_plastic_ptr, m_D_ptr));
619 typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dEP<
620 DIM,
I>(tau, ep, common_plastic_ptr, m_D_ptr));
623 typename P::template Assembly<A>::template OpCalculateConstraintsLhs_dTAU<
624 I>(tau, tau, common_plastic_ptr));
629template <
int DIM, AssemblyType A, IntegrationType I,
typename DomainEleOp>
631 std::string block_name,
Pip &pip,
632 std::string u, std::string ep,
636 using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
639 typename B::template OpGradTimesSymTensor<1, DIM, DIM>;
640 using OpInternalForcePiola =
643 auto [common_plastic_ptr, common_henky_ptr] =
644 createCommonPlasticOps<DIM, I, DomainEleOp>(m_field, block_name, pip, u,
645 ep, tau, 1., Sev::inform);
648 if (common_henky_ptr) {
649 pip.push_back(
new OpInternalForcePiola(
650 u, common_henky_ptr->getMatFirstPiolaStress()));
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#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.
FTensor::Index< 'm', SPACE_DIM > m
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto createCommonPlasticOps(MoFEM::Interface &m_field, std::string block_name, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string u, std::string ep, std::string tau, double scale, Sev sev)
MoFEMErrorCode opFactoryDomainReactions(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > Pip
FTensor::Index< 'M', 3 > M
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, std::string block_name, Pip &pip, boost::shared_ptr< MatrixDouble > mat_D_Ptr, boost::shared_ptr< CommonData::BlockParams > mat_params_ptr, double scale_value, Sev sev)
boost::shared_ptr< PlasticOps::CommonData > CommonPlasticPtr
FTensor::Index< 'N', 3 > N
FTensor::Index< 'I', 3 > I
[Common data]
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
FTensor::Index< 'J', 3 > J
boost::shared_ptr< HenckyOps::CommonData > CommonHenkyPtr
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
double young_modulus
Young modulus.
double C1_k
Kinematic hardening.
double Qinf
Saturation yield stress.
double visH
Viscous hardening.
double b_iso
Saturation exponent.
PetscBool is_large_strains
Large strains.
double sigmaY
Yield stress.
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
@ OPSPACE
operator do Work is execute on space data
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
boost::shared_ptr< MatrixDouble > mStrainPtr
MatrixDouble resFlowDstrainDot
auto getPlasticStrainPtr()
auto getPlasticStrainDotPtr()
boost::shared_ptr< MatrixDouble > mGradPtr
static std::array< int, 5 > activityData
MatrixDouble resCdPlasticStrain
VectorDouble plasticTauDot
boost::shared_ptr< MatrixDouble > mDPtr
[Common data set externally]
VectorDouble plasticSurface
[Common data set externally]
auto getPlasticTauDotPtr()
MatrixDouble plasticStrain
boost::shared_ptr< MatrixDouble > mStressPtr
MatrixDouble plasticStrainDot
MatrixDouble resFlowDstrain
auto getPlasticSurfacePtr()
std::array< double, LAST_PARAM > BlockParams
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Linear elastic problem]