v0.15.0
Loading...
Searching...
No Matches
Public Member Functions | List of all members
ADOLCPlasticity::OpCalculateStress< DIM, SMALL_STRAIN > Struct Template Reference
Inheritance diagram for ADOLCPlasticity::OpCalculateStress< DIM, SMALL_STRAIN >:
[legend]
Collaboration diagram for ADOLCPlasticity::OpCalculateStress< DIM, SMALL_STRAIN >:
[legend]

Public Member Functions

MoFEMErrorCode doWork (int side, EntityType type, DataForcesAndSourcesCore::EntData &data) override
 
- Public Member Functions inherited from ADOLCPlasticity::OpCalculateStressImpl< DIM, SMALL_STRAIN >
 OpCalculateStressImpl (MoFEM::Interface &m_field, boost::shared_ptr< CommonData > common_data_ptr, boost::shared_ptr< ClosestPointProjection > cp_ptr, bool calc_lhs)
 

Additional Inherited Members

- Protected Member Functions inherited from ADOLCPlasticity::OpCalculateStressImpl< DIM, SMALL_STRAIN >
MoFEMErrorCode getTags ()
 
MoFEMErrorCode setTagsData (const EntityHandle tet, const int nb_gauss_pts, const int nb_internal_variables)
 
- Protected Attributes inherited from ADOLCPlasticity::OpCalculateStressImpl< DIM, SMALL_STRAIN >
MoFEM::Interface & mField
 
boost::shared_ptr< CommonData > commonDataPtr
 
boost::shared_ptr< ClosestPointProjection > cpPtr
 
bool calcLhs
 
Tag thPlasticStrain
 
Tag thInternalVariables
 

Detailed Description

template<int DIM>
struct ADOLCPlasticity::OpCalculateStress< DIM, SMALL_STRAIN >

Definition at line 292 of file ADOLCPlasticity.cpp.

Member Function Documentation

◆ doWork()

template<int DIM>
MoFEMErrorCode ADOLCPlasticity::OpCalculateStress< DIM, SMALL_STRAIN >::doWork ( int side,
EntityType type,
DataForcesAndSourcesCore::EntData & data )
override

Definition at line 600 of file ADOLCPlasticity.cpp.

601 {
602 MoFEMFunctionBegin;
603
604 auto do_work_impl = [this](auto commonDataPtr, auto cpPtr) {
605 MoFEMFunctionBegin;
606
607 int nb_gauss_pts = this->getGaussPts().size2();
608
609 int nb_internal_variables = cpPtr->nbInternalVariables;
610 auto tet = this->getNumeredEntFiniteElementPtr()->getEnt();
611 CHKERR this->setTagsData(tet, nb_gauss_pts, nb_internal_variables);
612
613 commonDataPtr->activeVariablesW.resize(
614 nb_gauss_pts, 12 + cpPtr->nbInternalVariables, false);
615 commonDataPtr->gradientW.resize(nb_gauss_pts,
616 12 + cpPtr->nbInternalVariables, false);
617 commonDataPtr->materialTangentOperator.resize(36, nb_gauss_pts, false);
618 commonDataPtr->deltaGamma.resize(nb_gauss_pts);
619
620 FTensor::Index<'i', 3> i;
621 FTensor::Index<'j', 3> j;
622 FTensor::Index<'k', 3> k;
623 FTensor::Index<'l', 3> l;
624 FTensor::Index<'Z', 6> Z;
625 FTensor::Index<'Y', 6> Y;
626
627 // Code uses trial stress, until first solve of linear system of equations
628 int iter = 1;
629 if (this->getFEMethod()->snes_ctx != SnesMethod::CTX_SNESNONE) {
630 CHKERR SNESGetLinearSolveIterations(this->getFEMethod()->snes, &iter);
631 }
632
633 // b-bar
634 auto ave_tr_strain = calcAveStrain(
635 this->commonDataPtr->bBar, nb_gauss_pts,
636 getFTensor2FromMat<DIM, DIM>(commonDataPtr->gradAtGaussPts),
637 this->getFTensor0IntegrationWeight());
638
639 auto t_voight_stress_op = voight_to_stress_op();
640
641 auto t_grad = getFTensor2FromMat<DIM, DIM>(commonDataPtr->gradAtGaussPts);
642
643 auto t_Cep = getFTensor4DdgFromMat<3, 3, 1>(
644 this->commonDataPtr->materialTangentOperator);
645 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
646
647 auto tmp_active =
648 getVectorAdaptor(&(commonDataPtr->activeVariablesW(gg, 0)),
649 12 + cpPtr->nbInternalVariables);
650 cpPtr->activeVariablesW.swap(tmp_active);
651 auto tmp_gradient = getVectorAdaptor(&(commonDataPtr->gradientW(gg, 0)),
652 12 + cpPtr->nbInternalVariables);
653 cpPtr->gradientW.swap(tmp_gradient);
654
655 auto t_voigt_strain =
656 calcStrain(commonDataPtr->bBar, ave_tr_strain, t_grad,
657 getFTensor1FromPtr<6>(&(cpPtr->activeVariablesW[6])));
658 ++t_grad;
659
660 // Copy plastic strains and internal variables from tags
661 auto copy_plastic_strain_and_internal_variables = [&]() {
662 MoFEMFunctionBegin;
663
664 // Get data stored on Tags
665 VectorAdaptor plastic_strain_tags =
666 VectorAdaptor(6, ublas::shallow_array_adaptor<double>(
667 6, &commonDataPtr->plasticStrainPtr[gg * 6]));
668 VectorAdaptor internal_variables_tags = VectorAdaptor(
669 nb_internal_variables,
670 ublas::shallow_array_adaptor<double>(
671 nb_internal_variables,
672 &commonDataPtr
673 ->internalVariablesPtr[gg * nb_internal_variables]));
674
675 // Set values to closest point projection
676 cpPtr->plasticStrain0.resize(6, false);
677 noalias(cpPtr->plasticStrain0) = plastic_strain_tags;
678 cpPtr->internalVariables0.resize(nb_internal_variables, false);
679 noalias(cpPtr->internalVariables0) = internal_variables_tags;
680
681 auto cp_plastic_strain = cpPtr->getPlasticStrain();
682 auto cp_internal_variables = cpPtr->getInternalVariables();
683 noalias(cp_plastic_strain) = plastic_strain_tags;
684 noalias(cp_internal_variables) = internal_variables_tags;
685 MoFEMFunctionReturn(0);
686 };
687
688 CHKERR copy_plastic_strain_and_internal_variables();
689
690 auto closest_point_projection = [&](auto &recalculate_elastic_tangent) {
691 MoFEMFunctionBegin;
692 for (auto t : {ClosestPointProjection::H, ClosestPointProjection::Y,
693 ClosestPointProjection::W}) {
694 CHKERR cpPtr->setParams(t, recalculate_elastic_tangent);
695 }
696 CHKERR cpPtr->playW_NoHessian();
697 CHKERR cpPtr->playY_NoGradient();
698 double y = cpPtr->y;
699 cpPtr->deltaGamma = 0; // zero lagrange multiplier
700 if (iter > 0 && y > std::numeric_limits<double>::epsilon()) {
701 CHKERR cpPtr->solveClosestProjection();
702 }
703 commonDataPtr->deltaGamma[gg] = cpPtr->deltaGamma;
704 MoFEMFunctionReturn(0);
705 };
706
707 auto evaluate_consistent_tangent_matrix =
708 [&](auto &recalculate_elastic_tangent) {
709 MoFEMFunctionBegin;
710 if (recalculate_elastic_tangent)
711 CHKERR cpPtr->playW();
712
713 if (iter > 0 &&
714 cpPtr->deltaGamma > std::numeric_limits<double>::epsilon()) {
715
716 CHKERR cpPtr->consistentTangent();
717
718 auto &m = cpPtr->Cep;
719 // for generic case of non-associated plasticity consistent
720 // tangent matrix is non-symmetric
721 auto t_voight_cep =
722 FTensor::Tensor2<FTensor::PackPtr<double *, 0>, 6, 6>(
723 &m(0, 0), &m(0, 1), &m(0, 2), &m(0, 3), &m(0, 4),
724 &m(0, 5),
725
726 &m(1, 0), &m(1, 1), &m(1, 2), &m(1, 3), &m(1, 4),
727 &m(1, 5),
728
729 &m(2, 0), &m(2, 1), &m(2, 2), &m(2, 3), &m(2, 4),
730 &m(2, 5),
731
732 &m(3, 0), &m(3, 1), &m(3, 2), &m(3, 3), &m(3, 4),
733 &m(3, 5),
734
735 &m(4, 0), &m(4, 1), &m(4, 2), &m(4, 3), &m(4, 4),
736 &m(4, 5),
737
738 &m(5, 0), &m(5, 1), &m(5, 2), &m(5, 3), &m(5, 4), &m(5, 5)
739
740 );
741
742 t_Cep(i, j, k, l) =
743 t_voight_stress_op(i, j, Z) *
744 (t_voight_cep(Z, Y) * t_voight_stress_op(k, l, Y));
745
746 } else {
747
748 auto &m = cpPtr->C;
749 auto t_voight_cep =
750 FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 0>, 6>(
751 &m(0, 0), &m(0, 1), &m(0, 2), &m(0, 3), &m(0, 4),
752 &m(0, 5),
753
754 &m(1, 1), &m(1, 2), &m(1, 3), &m(1, 4), &m(1, 5),
755
756 &m(2, 2), &m(2, 3), &m(2, 4), &m(2, 5),
757
758 &m(3, 3), &m(3, 4), &m(3, 5),
759
760 &m(4, 4), &m(4, 5),
761
762 &m(5, 5));
763
764 t_Cep(i, j, k, l) =
765 t_voight_stress_op(i, j, Z) *
766 (t_voight_cep(Z, Y) * t_voight_stress_op(k, l, Y));
767 }
768
769 ++t_Cep;
770 MoFEMFunctionReturn(0);
771 };
772
773 // Always recalculate elastic tangent if first element
774 bool recalculate_elastic_tangent =
775 (this->getNinTheLoop() == 0) ? true : false;
776 CHKERR closest_point_projection(recalculate_elastic_tangent);
777 if (this->calcLhs)
778 CHKERR evaluate_consistent_tangent_matrix(recalculate_elastic_tangent);
779 }
780
781 auto calcs_stress_matrix = [&]() {
782 MoFEMFunctionBegin;
783 auto t_voight_stress_op = voight_to_stress_op();
784 commonDataPtr->stressMatrix.resize(6, commonDataPtr->gradientW.size1(),
785 false);
786 auto t_stess =
787 getFTensor2SymmetricFromMat<3>(commonDataPtr->stressMatrix);
788 auto t_voight_stress = commonDataPtr->getFTensor1StressAtGaussPts();
789 for (auto gg = 0; gg != commonDataPtr->gradientW.size1(); ++gg) {
790 t_stess(i, j) = t_voight_stress_op(i, j, Z) * t_voight_stress(Z);
791 ++t_voight_stress;
792 ++t_stess;
793 }
794 MoFEMFunctionReturn(0);
795 };
796
797 auto calcs_plastic_strain_matrix = [&]() {
798 MoFEMFunctionBegin;
799 auto t_voight_strain_op = voight_to_stress_op();
800 commonDataPtr->plasticStrainMatrix.resize(
801 6, commonDataPtr->activeVariablesW.size1(), false);
802 auto t_plastic_strain =
803 getFTensor2SymmetricFromMat<3>(commonDataPtr->plasticStrainMatrix);
804 auto t_voight_plastic_strain =
805 commonDataPtr->getFTensor1PlasticStrainAtGaussPts();
806 for (auto gg = 0; gg != commonDataPtr->activeVariablesW.size1(); ++gg) {
807 t_plastic_strain(i, j) =
808 t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
809 ++t_voight_plastic_strain;
810 ++t_plastic_strain;
811 }
812 MoFEMFunctionReturn(0);
813 };
814
815 CHKERR calcs_stress_matrix();
816 CHKERR calcs_plastic_strain_matrix();
817 MoFEMFunctionReturn(0);
818 };
819
820 CHKERR do_work_impl(this->commonDataPtr, this->cpPtr);
821
822 MoFEMFunctionReturn(0);
823}
FTensor::Tensor2_symmetric
Definition Tensor2_symmetric_value.hpp:14
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
ADOLCPlasticity::voight_to_stress_op
FTensor::Dg< double, 3, 6 > voight_to_stress_op()
Op convert Vight stress vector to stress tensor.
Definition ADOLCPlasticity.hpp:51
ADOLCPlasticity::calcAveStrain
double calcAveStrain(bool b_bar, const int nb_gauss_pts, FTensor::Tensor2< T1, DIM1, DIM2 > &&t_grad, FTensor::Tensor0< T2 > &&t_w)
Definition ADOLCPlasticity.cpp:18
ADOLCPlasticity::calcStrain
FTensor::Tensor1< T2, DIM3 > calcStrain(bool b_bar, double ave_tr_strain, FTensor::Tensor2< T1, DIM1, DIM2 > &t_grad, FTensor::Tensor1< T2, DIM3 > &&t_voight_strain)
Definition ADOLCPlasticity.cpp:43
MoFEM::Types::VectorAdaptor
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition Types.hpp:115
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition Templates.hpp:31
MoFEM::getFTensor4DdgFromMat
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
Definition Templates.hpp:418
MoFEM::getFTensor2FromMat
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
Definition Templates.hpp:278
MoFEM::getFTensor2SymmetricFromMat
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
Definition Templates.hpp:334
MoFEM::getFTensor1FromPtr
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition Templates.hpp:838
t
constexpr double t
plate stiffness
Definition plate.cpp:58
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
ADOLCPlasticity::OpCalculateStressImpl< DIM, SMALL_STRAIN >::setTagsData
MoFEMErrorCode setTagsData(const EntityHandle tet, const int nb_gauss_pts, const int nb_internal_variables)
Definition ADOLCPlasticity.cpp:325
ADOLCPlasticity::OpCalculateStressImpl< DIM, SMALL_STRAIN >::cpPtr
boost::shared_ptr< ClosestPointProjection > cpPtr
Definition ADOLCPlasticity.cpp:262
The documentation for this struct was generated from the following file:
Generated by Doxygen 1.12.0 and hosted at University of Glasgow