diff --git a/NonLinearSolver/materialLaw/STensorOperations.cpp b/NonLinearSolver/materialLaw/STensorOperations.cpp index 059751a9b5250da0c0079e4e231c1f7ca333be28..2add110c83e40c330c706c6403b121b9fe5822c0 100644 --- a/NonLinearSolver/materialLaw/STensorOperations.cpp +++ b/NonLinearSolver/materialLaw/STensorOperations.cpp @@ -632,7 +632,7 @@ void STensorOperation::alignEigenDecomposition_NormBased(const STensor3& A, cons a3_align = eigenValReal1_permuted(2); } -void STensorOperation::alignEigenDecomposition_NormBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, bool rotate){ +void STensorOperation::alignEigenDecomposition_NormBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, int& rotate){ // Align A to B, where B is fixed // Initialise 1 @@ -684,9 +684,9 @@ void STensorOperation::alignEigenDecomposition_NormBased_indexOnly(const STensor // Assign the aligned indices normamised_norm = norm_min; // /pow(A.dotprod(),0.5); - Msg::Info("normamised_norm of A = %f",normamised_norm); + // Msg::Info("normamised_norm of A = %f",normamised_norm); if(normamised_norm > 1.e-5){ - rotate = true; + rotate = 1; alignedIndex = eigenValReal1_permuted_Index; } } @@ -765,7 +765,7 @@ void STensorOperation::alignEigenDecomposition_TangentBased_indexOnly(const STen alignedIndex = eigenValReal1_permuted_Index; } -void STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, bool rotate){ +void STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, int& rotate){ // Align A to B, where B is fixed // Initialise 1 @@ -825,24 +825,26 @@ void STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly( // Generate permutated eigVal vector static fullVector<double> eigenValReal1_permuted(3); static fullVector<int> eigenValReal1_permuted_Index(3); eigenValReal1_permuted_Index(0) = 0; eigenValReal1_permuted_Index(1) = 1; eigenValReal1_permuted_Index(2) = 2; - for (int p=0; p<3; p++) - for (int q=0; q<3; q++){ - if(q!=p){ - for (int r=0; r<3; r++){ - if(r!=p && r!=q){ - - STensorOperation::zero(A_reformed); - for (int k=0; k<3; k++) - for (int l=0; l<3; l++) - A_reformed(k,l) += eigenValReal1(p)*B1(k,l) + eigenValReal1(q)*B2(k,l) + eigenValReal1(r)*B3(k,l); - A_reformed -= A; - STensorOperation::doubleContractionSTensor3(A_reformed,A_reformed,norm1); - norm1 = pow(norm1,0.5); - - norm2 = theta(p,0) + theta(q,1) + theta(r,2); - - if(norm1 < norm_min1){ // if( ( - norm + norm_min)>1.e-6 ) - if(norm2 < norm_min2){ + double tol(1.e-5); + if ( ( abs(eigenValReal1(0)-eigenValReal1(1)) > tol ) || ( abs(eigenValReal1(0)-eigenValReal1(2)) > tol ) || ( abs(eigenValReal1(1)-eigenValReal1(2)) > tol ) ){ + for (int p=0; p<3; p++) + for (int q=0; q<3; q++){ + if(q!=p){ + for (int r=0; r<3; r++){ + if(r!=p && r!=q){ + + STensorOperation::zero(A_reformed); + for (int k=0; k<3; k++) + for (int l=0; l<3; l++) + A_reformed(k,l) += eigenValReal1(p)*B1(k,l) + eigenValReal1(q)*B2(k,l) + eigenValReal1(r)*B3(k,l); + A_reformed -= A; + STensorOperation::doubleContractionSTensor3(A_reformed,A_reformed,norm1); + norm1 = pow(norm1,0.5); + + norm2 = theta(p,0) + theta(q,1) + theta(r,2); + + if(norm1 < norm_min1){ // if( ( - norm + norm_min)>1.e-6 ) + // if(norm2 < norm_min2){ norm_min1 = norm1; norm_min2 = norm2; eigenValReal1_permuted(0) = eigenValReal1(p); @@ -851,19 +853,69 @@ void STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly( eigenValReal1_permuted_Index(1) = q; eigenValReal1_permuted(2) = eigenValReal1(r); eigenValReal1_permuted_Index(2) = r; + // } } } } + } + } + alignedIndex = eigenValReal1_permuted_Index; + } + else{ + // Compare and pick the smallest theta1 for the 1st aligned Eigenvector wrt Eigenvectors2 + double smol1 = theta(0,0); + double smol2 = theta(1,1); + double smol3 = theta(2,2); + int index1 = 0, index2 = 1, index3 = 2; + + // pick 1 + int temp1 = 0; + for (int l=0; l<3; l++){ + if(temp1==0){ + smol1 = theta(l,0); + index1 = l; + temp1 = 1; + } + else if(theta(l,0)<smol1 && temp1!=0){ + smol1 = theta(l,0); + index1 = l; + temp1 = 1; + } + } + // pick 2 + int temp2 = 0; + for (int l=0; l<3; l++){ + if(l!=index1){ + if(temp1==0){ + smol2 = theta(l,1); + index2 = l; + temp2 = 1; + } + if(theta(l,1)<smol2 && temp2!=0){ + smol2 = theta(l,1); + index2 = l; + temp2 = 1; + } } } - // Assign the aligned indices - if (eigenValReal1_permuted_Index(0) != 0 || eigenValReal1_permuted_Index(1) != 1 || eigenValReal1_permuted_Index(2) != 2){ - rotate = true; - alignedIndex = eigenValReal1_permuted_Index; + // pick 3 + for (int l=0; l<3; l++){ + if(l!=index1 && l!=index2) index3 = l; + } + + // Assign the aligned indices + alignedIndex(0) = index1; + alignedIndex(1) = index2; + alignedIndex(2) = index3; } + /* + if (norm_min2 > pi/2){ + rotate = 1; + alignedIndex = eigenValReal1_permuted_Index; + }*/ } void STensorOperation::VRDecomposition(const STensor3& a, STensor3& V, STensor3&R) diff --git a/NonLinearSolver/materialLaw/STensorOperations.h b/NonLinearSolver/materialLaw/STensorOperations.h index a60a00393f0a49f5c8d3206c4a1a2102c3d64b60..744a5fa167a682dc9accb21112287d63d4ecb293 100644 --- a/NonLinearSolver/materialLaw/STensorOperations.h +++ b/NonLinearSolver/materialLaw/STensorOperations.h @@ -1933,9 +1933,9 @@ namespace STensorOperation{ } - void alignEigenDecomposition_NormBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, bool rotate = false); + void alignEigenDecomposition_NormBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, int& rotate); - void alignEigenDecomposition_NormAndDotProductBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, bool rotate = false); + void alignEigenDecomposition_NormAndDotProductBased_indexOnly(const STensor3& A, const STensor3& B, fullVector<int>& alignedIndex, int& rotate); void alignEigenDecomposition_NormBased(const STensor3& A, const STensor3& B, const double& b1, const double& b2, const double& b3, const STensor3& B1, const STensor3& B2, const STensor3& B3, diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.cpp index 6890f81418092db948ab297712f4ff8d19dfce2a..dc5e62105a8972667d213512ae72bd84fdda98b7 100644 --- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.cpp @@ -482,7 +482,7 @@ void mlawNonLinearTVENonLinearTVP2::predictorCorrector_TVP_nonAssociatedFlow_non // NEW EIGEN static fullVector<int> alignedIndex(3); alignedIndex(0) = 0; alignedIndex(1) = 1; alignedIndex(2) = 2; - bool rotate = false; + int rotate = 0; STensor3& R = q1->_R; STensor43& dRdEe = q1->_dRdEe; STensor43& dRtdEe = q1->_dRtdEe; @@ -502,10 +502,10 @@ void mlawNonLinearTVENonLinearTVP2::predictorCorrector_TVP_nonAssociatedFlow_non } else if(_rotationCorrectionScheme == 2){ Ee0_rot = q0->_Ee; - // STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); + STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); // STensorOperation::alignEigenDecomposition_EigenVectorDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex); // STensorOperation::alignEigenDecomposition_TangentBased_indexOnly(q0->_Ee,Ee,alignedIndex); - STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); + // STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); } // NEW EIGEN @@ -515,7 +515,6 @@ void mlawNonLinearTVENonLinearTVP2::predictorCorrector_TVP_nonAssociatedFlow_non mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm,&Ge_TrEeTensor); } else{ - // mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm,&Ge_TrEeTensor); mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,stiff,alignedIndex,rotate,Bd_stiffnessTerm,&Ge_TrEeTensor,&rotationStiffness); // NEW EIGEN } Ge_TrEeTensor_pr = Ge_TrEeTensor; @@ -526,8 +525,7 @@ void mlawNonLinearTVENonLinearTVP2::predictorCorrector_TVP_nonAssociatedFlow_non mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm); } else{ - // mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm); - mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,stiff,alignedIndex,rotate,Bd_stiffnessTerm,NULL,&rotationStiffness); // NEW EIGEN + mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ce,Ee,Ee0_rot,q0,q1,Ke,Ge,T0,T,stiff,alignedIndex,rotate,Bd_stiffnessTerm,NULL,&rotationStiffness); // NEW EIGEN } rotationStiffness_pr = rotationStiffness; } @@ -808,11 +806,11 @@ void mlawNonLinearTVENonLinearTVP2::predictorCorrector_TVP_nonAssociatedFlow_non if (_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){ // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3 getIterated_DPhi(T0,T,q0,q1,Gamma,Cxtr,Cxdev,Cepr,Eepr,trXn,devXn,Ke,Ge,Ge_Tensor,ptilde,devPhi,Phi,N,expGN,dexpAdA, - Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,alignedIndex,&Ge_TrEeTensor,&rotationStiffness); + Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,alignedIndex,rotate,&Ge_TrEeTensor,&rotationStiffness); } else{ getIterated_DPhi(T0,T,q0,q1,Gamma,Cxtr,Cxdev,Cepr,Eepr,trXn,devXn,Ke,Ge,Ge_Tensor,ptilde,devPhi,Phi,N,expGN,dexpAdA, - Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,alignedIndex,NULL,&rotationStiffness); + Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,alignedIndex,rotate,NULL,&rotationStiffness); } PhiEq = sqrt(1.5*devPhi.dotprod()); @@ -2197,7 +2195,7 @@ void mlawNonLinearTVENonLinearTVP2::getIterated_DPhi(const double& T0, const dou double& Ke, double& Ge, STensor43& Ge_Tensor, double& ptilde, STensor3& devPhi, STensor3& Phi, STensor3& N, STensor3& expGN, STensor43& dexpAdA, - STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv, fullVector<int>& alignedIndex, bool rotate, + STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv, fullVector<int>& alignedIndex, int& rotate, STensor43* Ge_TrEeTensor, STensor43* rotationStiffness) const{ // This function calculates Phi iteratively. @@ -2248,10 +2246,10 @@ void mlawNonLinearTVENonLinearTVP2::getIterated_DPhi(const double& T0, const dou } else if(_rotationCorrectionScheme == 2){ Ee0_rot = q0->_Ee; - // STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); + STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); // STensorOperation::alignEigenDecomposition_EigenVectorDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex); // STensorOperation::alignEigenDecomposition_TangentBased_indexOnly(q0->_Ee,Ee,alignedIndex); - STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); + // STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); } // NEW EIGEN @@ -2494,10 +2492,10 @@ void mlawNonLinearTVENonLinearTVP2::getIterated_DPhi(const double& T0, const dou // EIGEN NEW if(_rotationCorrectionScheme == 2){ Ee0_rot = q0->_Ee; - // STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); + STensorOperation::alignEigenDecomposition_NormBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); // STensorOperation::alignEigenDecomposition_EigenVectorDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex); // STensorOperation::alignEigenDecomposition_TangentBased_indexOnly(q0->_Ee,Ee,alignedIndex); - STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); + // STensorOperation::alignEigenDecomposition_NormAndDotProductBased_indexOnly(q0->_Ee,Ee,alignedIndex,rotate); } // STensorOperation::getEigenDecomposition(Ee,e1,e2,e3,E1,E2,E3); diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.h b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.h index a7f40ccf6503aa404d6ec01a72510c29268c1ea6..4461d31a7757a7fa86bdc1e67520283039e8f252 100644 --- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.h +++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP2.h @@ -31,7 +31,7 @@ class mlawNonLinearTVENonLinearTVP2 : public mlawNonLinearTVP{ double& Ke, double& Ge, STensor43& Ge_Tensor, double& ptilde, STensor3& devPhi, STensor3& Phi, STensor3& N, STensor3& expGN, STensor43& dexpAdA, - STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv, fullVector<int>& alignedIndex, bool rotate = false, + STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv, fullVector<int>& alignedIndex, int& rotate, STensor43* Ge_TrEeTensor = NULL, STensor43* rotationStiffness = NULL) const; virtual void getDho(const double& Gamma, const STensor3& Cepr, const STensor3& Ceinvpr, const STensor3& KS, diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp index 049251b86d3dcf789f4bdb987748ee8870b85cb6..419f6158a2c2ed82782c099b45ddfc96ae904bf6 100644 --- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp @@ -3112,7 +3112,7 @@ void mlawNonLinearTVM::get_Ee0s_Ai0s_ders(const STensor3& Ce, const STensor3& Ee void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ce, const STensor3& Ee, const STensor3& Ee0, const IPNonLinearTVM *q0, IPNonLinearTVM *q1, double& Ke, double& Ge, const double T0, const double T1, - const bool stiff, const fullVector<int>& alignedIndex, const bool rotate, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2, STensor43* rotationStiffness) const{ + const bool stiff, const fullVector<int>& alignedIndex, const int& rotate, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2, STensor43* rotationStiffness) const{ // Caution!! -> For _rotationCorrectionScheme == 0 and 1, Ee0 is already rotated. _R and its derivatives are in the IP @@ -3157,7 +3157,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ce, co // mlawNonLinearTVM::getEe0s2(Ce,Ee,devEe0,devEe0s,dDevEe0sdEe); q1->_A_rot = q0->_A; // mlawNonLinearTVM::getAi0s2(Ce,Ee,q0->_A,q1->_A_rot,dAi0sdEe); - if(rotate){mlawNonLinearTVM::get_Ee0s_Ai0s_ders(Ce,Ee,devEe0,q0->_A,alignedIndex,devEe0s,dDevEe0sdEe,q1->_A_rot,dAi0sdEe);} + mlawNonLinearTVM::get_Ee0s_Ai0s_ders(Ce,Ee,devEe0,q0->_A,alignedIndex,devEe0s,dDevEe0sdEe,q1->_A_rot,dAi0sdEe); // devEe0s.print("after"); // Msg::Info("alignedIndex = %d, %d, %d",alignedIndex(0),alignedIndex(1),alignedIndex(2)); // FLE } @@ -3233,7 +3233,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ce, co q1->_elasticShearPropertyScaleFactor = B; q1->_dAv_dTrEe = dA_dTrEe; q1->_dBd_dDevEe = dB_dDevEe; - + // Add to stress static STensor3 sigExtra; @@ -3764,6 +3764,8 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ce, co q1->_kirchhoff(1,1) += p; q1->_kirchhoff(2,2) += p; q1->_trCorKirinf_TVE = p; + + if(rotate == 0){STensorOperation::zero(*rotationStiffness);} } // ======================================================================================================================================== diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h index d32f1ffef97a3ffcfc676561e4a6c6e0a30ba301..62926d4eb1d478283e71d14d338fb2d60f6ef3ae 100644 --- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h +++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h @@ -137,7 +137,7 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{ void ThermoViscoElasticPredictor_forTVP(const STensor3& Ce, const STensor3& Ee, const STensor3& Ee0, const IPNonLinearTVM *q0, IPNonLinearTVM *q1, double& Ke, double& Ge, const double T0, const double T1, - const bool stiff, const fullVector<int>& alignedIndex, const bool rotate, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2 = NULL, STensor43* rotationStiffness = NULL) const; + const bool stiff, const fullVector<int>& alignedIndex, const int& rotate, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2 = NULL, STensor43* rotationStiffness = NULL) const; void isotropicHookTensor(const double K, const double G, STensor43& L) const; diff --git a/dG3D/benchmarks/nonLinearTVP_cube_Rotation_biaxialShear/cubeTVP.py b/dG3D/benchmarks/nonLinearTVP_cube_Rotation_biaxialShear/cubeTVP.py index efa138a92d38decba8e009ef9c8713a422052b0a..05bfd11d2c7567db8af36e8e47096948c364fddd 100644 --- a/dG3D/benchmarks/nonLinearTVP_cube_Rotation_biaxialShear/cubeTVP.py +++ b/dG3D/benchmarks/nonLinearTVP_cube_Rotation_biaxialShear/cubeTVP.py @@ -1,8 +1,8 @@ #coding-Utf-8-*- from gmshpy import * -# from dG3Dpy import* -from dG3DpyDebug import* +from dG3Dpy import* +# from dG3DpyDebug import* from math import* import csv import numpy as np @@ -93,7 +93,7 @@ law1.setViscoelasticMethod(0) law1.setShiftFactorConstantsWLF(C1,C2) law1.setReferenceTemperature(Tref) -law1.useRotationCorrectionBool(False,2) # bool, int_Scheme = 0 -> R0, 1 -> R1 +law1.useRotationCorrectionBool(True,2) # bool, int_Scheme = 0 -> R0, 1 -> R1 """ law1.setExtraBranchNLType(5)