v0.15.0
Loading...
Searching...
No Matches
PlasticOps::OpCalculatePlasticityImpl< DIM, GAUSS, DomainEleOp > Struct Template Reference

#include "tutorials/adv-0/src/PlasticOpsGeneric.hpp"

Inheritance diagram for PlasticOps::OpCalculatePlasticityImpl< DIM, GAUSS, DomainEleOp >:
[legend]
Collaboration diagram for PlasticOps::OpCalculatePlasticityImpl< DIM, GAUSS, DomainEleOp >:
[legend]

Public Member Functions

 OpCalculatePlasticityImpl (const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > m_D_ptr)
 
MoFEMErrorCode doWork (int side, EntityType type, EntData &data)
 

Protected Attributes

boost::shared_ptr< CommonDatacommonDataPtr
 
boost::shared_ptr< MatrixDouble > mDPtr
 

Detailed Description

template<int DIM, typename DomainEleOp>
struct PlasticOps::OpCalculatePlasticityImpl< DIM, GAUSS, DomainEleOp >

Definition at line 469 of file PlasticOpsGeneric.hpp.

Constructor & Destructor Documentation

◆ OpCalculatePlasticityImpl()

template<int DIM, typename DomainEleOp >
PlasticOps::OpCalculatePlasticityImpl< DIM, GAUSS, DomainEleOp >::OpCalculatePlasticityImpl ( const std::string field_name,
boost::shared_ptr< CommonData > common_data_ptr,
boost::shared_ptr< MatrixDouble > m_D_ptr )
Examples
PlasticOpsGeneric.hpp.

Definition at line 481 of file PlasticOpsGeneric.hpp.

484 : DomainEleOp(field_name, DomainEleOp::OPROW),
485 commonDataPtr(common_data_ptr), mDPtr(m_D_ptr) {
486 // Opetor is only executed for vertices
487 std::fill(&DomainEleOp::doEntities[MBEDGE],
488 &DomainEleOp::doEntities[MBMAXTYPE], false);
489}
DomainEle::UserDataOperator DomainEleOp
Finire element operator type.
constexpr auto field_name

Member Function Documentation

◆ doWork()

template<int DIM, typename DomainEleOp >
MoFEMErrorCode PlasticOps::OpCalculatePlasticityImpl< DIM, GAUSS, DomainEleOp >::doWork ( int side,
EntityType type,
EntData & data )

< material parameters

Examples
PlasticOpsGeneric.hpp.

Definition at line 492 of file PlasticOpsGeneric.hpp.

493 {
495
496 FTensor::Index<'i', DIM> i;
497 FTensor::Index<'j', DIM> j;
498 FTensor::Index<'k', DIM> k;
499 FTensor::Index<'l', DIM> l;
500 FTensor::Index<'m', DIM> m;
501 FTensor::Index<'n', DIM> n;
502
503 auto &params = commonDataPtr->blockParams; ///< material parameters
504
505 auto nb_gauss_pts = DomainEleOp::getGaussPts().size2();
506 auto t_w = DomainEleOp::getFTensor0IntegrationWeight();
507 auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
508 auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
509 auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
510 auto t_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticFlow);
511 auto t_plastic_strain =
512 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrain);
513 auto t_plastic_strain_dot =
514 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->plasticStrainDot);
515 auto t_stress =
516 getFTensor2SymmetricFromMat<DIM>(*(commonDataPtr->mStressPtr));
517
518 auto t_D = getFTensor4DdgFromMat<DIM, DIM, 0>(*commonDataPtr->mDPtr);
519 auto t_D_Op = getFTensor4DdgFromMat<DIM, DIM, 0>(*mDPtr);
520
521 auto t_diff_plastic_strain = diff_tensor(FTensor::Number<DIM>());
522 auto t_diff_deviator = diff_deviator(diff_tensor(FTensor::Number<DIM>()));
523
524 FTensor::Ddg<double, DIM, DIM> t_flow_dir_dstress;
525 FTensor::Ddg<double, DIM, DIM> t_flow_dir_dstrain;
526 t_flow_dir_dstress(i, j, k, l) =
527 1.5 * (t_diff_deviator(M, N, i, j) * t_diff_deviator(M, N, k, l));
528 t_flow_dir_dstrain(i, j, k, l) =
529 t_flow_dir_dstress(i, j, m, n) * t_D_Op(m, n, k, l);
530
531
532 auto t_alpha_dir =
534
535 commonDataPtr->resC.resize(nb_gauss_pts, false);
536 commonDataPtr->resCdTau.resize(nb_gauss_pts, false);
537 commonDataPtr->resCdStrain.resize(size_symm, nb_gauss_pts, false);
538 commonDataPtr->resCdPlasticStrain.resize(size_symm, nb_gauss_pts, false);
539 commonDataPtr->resFlow.resize(size_symm, nb_gauss_pts, false);
540 commonDataPtr->resFlowDtau.resize(size_symm, nb_gauss_pts, false);
541 commonDataPtr->resFlowDstrain.resize(size_symm * size_symm, nb_gauss_pts,
542 false);
543 commonDataPtr->resFlowDstrainDot.resize(size_symm * size_symm, nb_gauss_pts,
544 false);
545
546 commonDataPtr->resC.clear();
547 commonDataPtr->resCdTau.clear();
548 commonDataPtr->resCdStrain.clear();
549 commonDataPtr->resCdPlasticStrain.clear();
550 commonDataPtr->resFlow.clear();
551 commonDataPtr->resFlowDtau.clear();
552 commonDataPtr->resFlowDstrain.clear();
553 commonDataPtr->resFlowDstrainDot.clear();
554
555 auto t_res_c = getFTensor0FromVec(commonDataPtr->resC);
556 auto t_res_c_dtau = getFTensor0FromVec(commonDataPtr->resCdTau);
557 auto t_res_c_dstrain =
558 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdStrain);
559 auto t_res_c_plastic_strain =
560 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resCdPlasticStrain);
561 auto t_res_flow = getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlow);
562 auto t_res_flow_dtau =
563 getFTensor2SymmetricFromMat<DIM>(commonDataPtr->resFlowDtau);
564 auto t_res_flow_dstrain =
565 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrain);
566 auto t_res_flow_dplastic_strain =
567 getFTensor4DdgFromMat<DIM, DIM>(commonDataPtr->resFlowDstrainDot);
568
569 auto next = [&]() {
570 ++t_tau;
571 ++t_tau_dot;
572 ++t_f;
573 ++t_flow;
574 ++t_plastic_strain;
575 ++t_plastic_strain_dot;
576 ++t_stress;
577 ++t_res_c;
578 ++t_res_c_dtau;
579 ++t_res_c_dstrain;
580 ++t_res_c_plastic_strain;
581 ++t_res_flow;
582 ++t_res_flow_dtau;
583 ++t_res_flow_dstrain;
584 ++t_res_flow_dplastic_strain;
585 ++t_w;
586 };
587
588 auto get_avtive_pts = [&]() {
589 int nb_points_avtive_on_elem = 0;
590 int nb_points_on_elem = 0;
591
592 auto t_tau = getFTensor0FromVec(commonDataPtr->plasticTau);
593 auto t_tau_dot = getFTensor0FromVec(commonDataPtr->plasticTauDot);
594 auto t_f = getFTensor0FromVec(commonDataPtr->plasticSurface);
595 auto t_plastic_strain_dot =
596 getFTensor2SymmetricFromMat<SPACE_DIM>(commonDataPtr->plasticStrainDot);
597
598 auto dt = this->getTStimeStep();
599
600 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
601 auto eqiv = equivalent_strain_dot(t_plastic_strain_dot);
602 const auto ww = w(
603 eqiv, t_tau_dot, t_f,
604 iso_hardening(t_tau, params[CommonData::H], params[CommonData::QINF],
605 params[CommonData::BISO], params[CommonData::SIGMA_Y]),
606 params[CommonData::SIGMA_Y]);
607 const auto sign_ww = constrian_sign(ww, dt);
608
609 ++nb_points_on_elem;
610 if (sign_ww > 0) {
611 ++nb_points_avtive_on_elem;
612 }
613
614 ++t_tau;
615 ++t_tau_dot;
616 ++t_f;
617 ++t_plastic_strain_dot;
618 }
619
620 int &active_points = PlasticOps::CommonData::activityData[0];
621 int &avtive_full_elems = PlasticOps::CommonData::activityData[1];
622 int &avtive_elems = PlasticOps::CommonData::activityData[2];
623 int &nb_points = PlasticOps::CommonData::activityData[3];
624 int &nb_elements = PlasticOps::CommonData::activityData[4];
625
626 ++nb_elements;
627 nb_points += nb_points_on_elem;
628 if (nb_points_avtive_on_elem > 0) {
629 ++avtive_elems;
630 active_points += nb_points_avtive_on_elem;
631 if (nb_points_avtive_on_elem == nb_points_on_elem) {
632 ++avtive_full_elems;
633 }
634 }
635
636 if (nb_points_avtive_on_elem != nb_points_on_elem)
637 return 1;
638 else
639 return 0;
640 };
641
642 if (DomainEleOp::getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
643 get_avtive_pts();
644 }
645
646 auto dt = this->getTStimeStep();
647 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
648
649 auto eqiv = equivalent_strain_dot(t_plastic_strain_dot);
650 auto t_diff_eqiv = diff_equivalent_strain_dot(eqiv, t_plastic_strain_dot,
651 t_diff_plastic_strain,
653
654 const auto sigma_y =
655 iso_hardening(t_tau, params[CommonData::H], params[CommonData::QINF],
656 params[CommonData::BISO], params[CommonData::SIGMA_Y]);
657 const auto d_sigma_y =
658 iso_hardening_dtau(t_tau, params[CommonData::H],
659 params[CommonData::QINF], params[CommonData::BISO]);
660
661 auto ww = w(eqiv, t_tau_dot, t_f, sigma_y, params[CommonData::SIGMA_Y]);
662 auto abs_ww = constrain_abs(ww, dt);
663 auto sign_ww = constrian_sign(ww, dt);
664
665 auto c = constraint(eqiv, t_tau_dot, t_f, sigma_y, abs_ww,
666 params[CommonData::VIS_H], params[CommonData::SIGMA_Y]);
667 auto c_dot_tau = diff_constrain_ddot_tau(sign_ww, eqiv, t_tau_dot,
668 params[CommonData::VIS_H],
669 params[CommonData::SIGMA_Y]);
670 auto c_equiv = diff_constrain_deqiv(sign_ww, eqiv, t_tau_dot,
671 params[CommonData::SIGMA_Y]);
672 auto c_sigma_y = diff_constrain_dsigma_y(sign_ww);
673 auto c_f = diff_constrain_df(sign_ww);
674
675 auto t_dev_stress = deviator(
676
677 t_stress, trace(t_stress),
678
679 kinematic_hardening(t_plastic_strain, params[CommonData::C1_k])
680
681 );
682
684 t_flow_dir(k, l) = 1.5 * (t_dev_stress(I, J) * t_diff_deviator(I, J, k, l));
686 t_flow_dstrain(i, j) = t_flow(k, l) * t_D_Op(k, l, i, j);
687
688 auto get_res_c = [&]() { return c; };
689
690 auto get_res_c_dstrain = [&](auto &t_diff_res) {
691 t_diff_res(i, j) = c_f * t_flow_dstrain(i, j);
692 };
693
694 auto get_res_c_dplastic_strain = [&](auto &t_diff_res) {
695 t_diff_res(i, j) = (this->getTSa() * c_equiv) * t_diff_eqiv(i, j);
696 t_diff_res(k, l) -= c_f * t_flow(i, j) * t_alpha_dir(i, j, k, l);
697 };
698
699 auto get_res_c_dtau = [&]() {
700 return this->getTSa() * c_dot_tau + c_sigma_y * d_sigma_y;
701 };
702
703 auto get_res_c_plastic_strain = [&](auto &t_diff_res) {
704 t_diff_res(k, l) = -c_f * t_flow(i, j) * t_alpha_dir(i, j, k, l);
705 };
706
707 auto get_res_flow = [&](auto &t_res_flow) {
708 const auto a = sigma_y;
709 const auto b = t_tau_dot;
710 t_res_flow(k, l) = a * t_plastic_strain_dot(k, l) - b * t_flow_dir(k, l);
711 };
712
713 auto get_res_flow_dtau = [&](auto &t_res_flow_dtau) {
714 const auto da = d_sigma_y;
715 const auto db = this->getTSa();
716 t_res_flow_dtau(k, l) =
717 da * t_plastic_strain_dot(k, l) - db * t_flow_dir(k, l);
718 };
719
720 auto get_res_flow_dstrain = [&](auto &t_res_flow_dstrain) {
721 const auto b = t_tau_dot;
722 t_res_flow_dstrain(m, n, k, l) = -t_flow_dir_dstrain(m, n, k, l) * b;
723 };
724
725 auto get_res_flow_dplastic_strain = [&](auto &t_res_flow_dplastic_strain) {
726 const auto a = sigma_y;
727 t_res_flow_dplastic_strain(m, n, k, l) =
728 (a * this->getTSa()) * t_diff_plastic_strain(m, n, k, l);
729 const auto b = t_tau_dot;
730 t_res_flow_dplastic_strain(m, n, i, j) +=
731 (t_flow_dir_dstrain(m, n, k, l) * t_alpha_dir(k, l, i, j)) * b;
732 };
733
734 t_res_c = get_res_c();
735 get_res_flow(t_res_flow);
736
737 if (this->getTSCtx() == TSMethod::TSContext::CTX_TSSETIJACOBIAN) {
738 t_res_c_dtau = get_res_c_dtau();
739 get_res_c_dstrain(t_res_c_dstrain);
740 get_res_c_dplastic_strain(t_res_c_plastic_strain);
741 get_res_flow_dtau(t_res_flow_dtau);
742 get_res_flow_dstrain(t_res_flow_dstrain);
743 get_res_flow_dplastic_strain(t_res_flow_dplastic_strain);
744 }
745
746 next();
747 }
748
750}
constexpr double a
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FTensor::Index< 'i', SPACE_DIM > i
double dt
const double c
speed of light (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double diff_constrain_ddot_tau(double sign, double eqiv, double dot_tau, double vis_H, double sigma_Y)
double constrian_sign(double x, double dt)
FTensor::Index< 'M', 3 > M
double diff_constrain_deqiv(double sign, double eqiv, double dot_tau, double sigma_Y)
auto diff_equivalent_strain_dot(const T1 eqiv, T2 &t_plastic_strain_dot, T3 &t_diff_plastic_strain, FTensor::Number< DIM >)
double constraint(double eqiv, double dot_tau, double f, double sigma_y, double abs_w, double vis_H, double sigma_Y)
auto diff_constrain_df(double sign)
FTensor::Index< 'N', 3 > N
auto diff_constrain_dsigma_y(double sign)
FTensor::Index< 'I', 3 > I
[Common data]
FTensor::Index< 'J', 3 > J
auto diff_tensor(FTensor::Number< DIM >)
[Lambda functions]
auto diff_deviator(FTensor::Ddg< double, DIM, DIM > &&t_diff_stress, FTensor::Number< DIM >)
double trace(FTensor::Tensor2_symmetric< T, 2 > &t_stress)
auto deviator(FTensor::Tensor2_symmetric< T, DIM > &t_stress, double trace, FTensor::Tensor2_symmetric< double, DIM > &t_alpha, FTensor::Number< DIM >)
double w(double eqiv, double dot_tau, double f, double sigma_y, double sigma_Y)
auto equivalent_strain_dot(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain_dot)
double constrain_abs(double x, double dt)
auto kinematic_hardening(FTensor::Tensor2_symmetric< T, DIM > &t_plastic_strain, double C1_k)
Definition plastic.cpp:92
double iso_hardening_dtau(double tau, double H, double Qinf, double b_iso)
Definition plastic.cpp:78
constexpr auto size_symm
Definition plastic.cpp:42
double iso_hardening(double tau, double H, double Qinf, double b_iso, double sigmaY)
Definition plastic.cpp:73
auto kinematic_hardening_dplastic_strain(double C1_k)
Definition plastic.cpp:106
FTensor::Index< 'm', 3 > m
static std::array< int, 5 > activityData

Member Data Documentation

◆ commonDataPtr

template<int DIM, typename DomainEleOp >
boost::shared_ptr<CommonData> PlasticOps::OpCalculatePlasticityImpl< DIM, GAUSS, DomainEleOp >::commonDataPtr
protected

Definition at line 476 of file PlasticOpsGeneric.hpp.

◆ mDPtr

template<int DIM, typename DomainEleOp >
boost::shared_ptr<MatrixDouble> PlasticOps::OpCalculatePlasticityImpl< DIM, GAUSS, DomainEleOp >::mDPtr
protected

Definition at line 477 of file PlasticOpsGeneric.hpp.


The documentation for this struct was generated from the following file: