From 4c26be0d37977cfd9d22a0d7b2eb11e531a9df9b Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Tue, 22 Aug 2017 18:06:20 +0200 Subject: [PATCH 01/13] cleaning function prototypes --- .../materialLaw/mlawNonLocalDamageGurson.cpp | 88 ++++++++++--------- .../materialLaw/mlawNonLocalDamageGurson.h | 6 +- 2 files changed, 51 insertions(+), 43 deletions(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index 74916d1e8..6de39f4ec 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -244,6 +244,7 @@ void mlawNonLocalDamageGurson::constitutive( double fVstar = 0.; double tildefVstar = ipvcur->getNonLocalCorrectedPorosity(); + // avoid unexpected values if(tildefVstar < 0. ) tildefVstar = 0.; if(tildefVstar > _ffstar*(1-1.e-2)) tildefVstar = _ffstar*(1-1.e-2); @@ -255,7 +256,7 @@ void mlawNonLocalDamageGurson::constitutive( static STensor3 Fe1, kcor; static STensor43 Le; static STensor63 dLe; - predictorCorrector(tildefVstar,F0, Fn, ipvprev->getConstRefToFp(), Fp1, ipvprev, ipvcur, P, fVstar, stiff, + predictorCorrector(tildefVstar, Fn, ipvprev->getConstRefToFp(), Fp1, ipvprev, ipvcur, P, fVstar, stiff, dDeltaGammadC, dDeltaHatPdC, dtrNpdC, dDeltafVdC, dNpDevdC, dKCordC, dFpdC, dDeltaGammadtildefVstar, dDeltaHatPdtildefVstar, dtrNpdtildefVstar, dDeltafVdtildefVstar,dNpDevdtildefVstar, dKCordtildefVstar, dFpdtildefVstar, @@ -267,18 +268,18 @@ void mlawNonLocalDamageGurson::constitutive( if (_tangentByPerturbation) { this->tangent_full_perturbation(Tangent, dLocalCorrectedPorosityDStrain, - dStressDNonLocalCorrectedPorosity, dLocalCorrectedPorosityDNonLocalCorrectedPorosity, - P, F0, Fn, ipvprev->getConstRefToFp(),Fp1,tildefVstar, fVstar, ipvprev,ipvcur); + dStressDNonLocalCorrectedPorosity, dLocalCorrectedPorosityDNonLocalCorrectedPorosity, + P, Fn, ipvprev->getConstRefToFp(),Fp1,tildefVstar, fVstar, ipvprev,ipvcur); } else { this->tangent_full_analytic(Tangent, dLocalCorrectedPorosityDStrain, - dStressDNonLocalCorrectedPorosity, dLocalCorrectedPorosityDNonLocalCorrectedPorosity, - P,Fn, ipvprev->getConstRefToFp(),Fp1,tildefVstar, fVstar, ipvprev,ipvcur, - dKCordC, dFpdC,dDeltafVdC,dKCordtildefVstar, dFpdtildefVstar,dDeltafVdtildefVstar, - Fe1, kcor, Le, dLe); + dStressDNonLocalCorrectedPorosity, dLocalCorrectedPorosityDNonLocalCorrectedPorosity, + P,Fn, ipvprev->getConstRefToFp(),Fp1,tildefVstar, fVstar, ipvprev,ipvcur, + dKCordC, dFpdC,dDeltafVdC,dKCordtildefVstar, dFpdtildefVstar,dDeltafVdtildefVstar, + Fe1, kcor, Le, dLe); /* this->tangent_full_perturbation(Tangent, dLocalCorrectedPorosityDStrain, dStressDNonLocalCorrectedPorosity, dLocalCorrectedPorosityDNonLocalCorrectedPorosity, P,Fn, ipvprev->getConstRefToFp(),Fp1,tildefVstar, fVstar, ipvprev,ipvcur);*/ @@ -287,7 +288,7 @@ void mlawNonLocalDamageGurson::constitutive( } -void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STensor3& F0_, const STensor3& F1_, const STensor3& Fp0_, +void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STensor3& F1_, const STensor3& Fp0_, STensor3& Fp1_, const IPNonLocalDamageGurson *q0_, IPNonLocalDamageGurson *q1_, STensor3&P_, double &fVstar, bool stiff, STensor3 &dDeltaGammadC, STensor3 &dDeltaHatPdC, STensor3 &dtrNpdC, @@ -313,6 +314,12 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen double &fV = q1_->getRefToLocalPorosity(); fV = fV0; double DeltafV = 0.; + + // Plastic flow direction + double trNp = 0.; + static STensor3 NpDev,DGNp; + NpDev *= 0.; + DGNp *= 0.; // Derivatives dDeltaGammadC*=0.; @@ -326,8 +333,6 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen - - /* compute elastic predictor */ static STensor3 Fp1inv; static STensor3 Fe1pr; @@ -407,14 +412,11 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen }*/ - static fullMatrix<double> Jinv(4,4); Jinv.setAll(0.); - static fullMatrix<double> J(4,4); J.setAll(0.); - double trNp = 0.; - static STensor3 NpDev,DGNp; - NpDev *= 0.; - DGNp *= 0.; /* Test plasticity */ + static fullMatrix<double> Jinv(4,4); Jinv.setAll(0.); + static fullMatrix<double> J(4,4); J.setAll(0.); + double f = yieldFunction(kprDev, ppr, yield, tildefVstar); if(f > _tol) { @@ -452,10 +454,9 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen if (fV0 + DeltafV < 0.5*fV0){DeltafV = -0.5*fV0;} else if (fV0 + DeltafV > _ff-1.e-10){DeltafV = (_ff-1.e-10-fV0)*0.5;} } - bool CorrectedDeltafV_Flag(false); double DeltaGammapr = DeltaGamma; double Dfvinit = DeltafV; - + bool CorrectedDeltafV_Flag(false); // for porosity growth control // First residual int ite = 0; @@ -643,6 +644,7 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen //Msg::Info("Number of iterations ite = %d , with DeltaGamma = %e ",ite, DeltaGamma); /* final state */ + // Plastic correction DeltaGamma*Np NpDev = kprDev; NpDev*= (3./(yield*yield+6.*_mu*DeltaGamma)); DGNp=NpDev; @@ -650,7 +652,8 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen DGNp(1,1) +=trNp/3.; DGNp(2,2) +=trNp/3.; DGNp*=DeltaGamma; - + + // Plastic increment static STensor3 dFp; STensorOperation::expSTensor3(DGNp,_order,dFp, &ENp); @@ -665,11 +668,18 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen STensorOperation::multSTensor3FirstTranspose(Fe1,Fe1,Ce); //Ce =Fe1.transpose(); //Ce *= Fe1; - STensorOperation::logSTensor3(Ce,_order,logsqrtCe,&Le,&dLe); logsqrtCe*=0.5; // corrotational Kirchhoff stress + // we use the expression used to derive the stiffness as the log is not exact + kcor = NpDev; + kcor *= (-2.*_mu*DeltaGamma); + kcor +=kprDev; + kcor(0,0) +=ppr - _K*DeltaGamma*trNp; + kcor(1,1) +=ppr - _K*DeltaGamma*trNp; + kcor(2,2) +=ppr - _K*DeltaGamma*trNp; + /*static STensor3 kDev; double p; KirchhoffStress(kDev, p, logsqrtCe); @@ -678,13 +688,7 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen kcor(0,0)+=p; kcor(1,1)+=p; kcor(2,2)+=p;*/ - // we use the expression used to derive the stiffness as the log is not exact - kcor = NpDev; - kcor *= (-2.*_mu*DeltaGamma); - kcor +=kprDev; - kcor(0,0) +=ppr - _K*DeltaGamma*trNp; - kcor(1,1) +=ppr - _K*DeltaGamma*trNp; - kcor(2,2) +=ppr - _K*DeltaGamma*trNp; + } //compute here the extended damage growth to be used at next time step // this is an approximation to facilitate the convergence @@ -710,8 +714,8 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen yield, h, yield0, DeltaHatP, DeltaGamma, trNp, fV0, DeltafV, Lpr, ENp); } - else - { + else + { //for debug dump the internal variables to check with perturbation /*dDeltaGammadtildefVstar = DeltaGamma; dDeltaHatPdtildefVstar = DeltaHatP; @@ -720,7 +724,7 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen dNpDevdtildefVstar = kprDev*(3./(yield*yield+6.*_mu*DeltaGamma)); dKCordtildefVstar = kcor; dFpdtildefVstar = Fp1_;*/ - } + } //update internal variables hatP = hatP0 + DeltaHatP; @@ -728,22 +732,26 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen q1_->getRefToFp() = Fp1_; q1_->getRefToElasticEnergy()=deformationEnergy(Ce); + // coalescence fVstar = fV; if(fVstar > _fC) { fVstar = _fC+(_ffstar-_fC)/(_ff-_fC)*(fV-_fC); } - //corotational to PK1 + // corotational to PK1 static STensor3 S; static STensor3 FeS; - - S*=0; - for(int i=0; i<3; i++) - for(int j=0; j<3; j++) - for(int k=0; k<3; k++) - for(int l=0; l<3; l++) + for(int i=0; i<3; i++){ + for(int j=0; j<3; j++){ + S(i,j) = 0.; + for(int k=0; k<3; k++){ + for(int l=0; l<3; l++){ S(i,j)+= Le(i,j,k,l)*kcor(k,l); + } + } + } + } STensorOperation::multSTensor3(Fe1,S,FeS); STensorOperation::multSTensor3SecondTranspose(FeS,Fp1inv,P_); //P_=Fe1; @@ -907,7 +915,7 @@ void mlawNonLocalDamageGurson::tangent_full_perturbation(STensor43 &T_, STensor3 &dLocalCorrectedPorosityDStrain, STensor3 &dStressDNonLocalCorrectedPorosity, double &dLocalCorrectedPorosityDNonLocalCorrectedPorosity, - const STensor3 &P_,const STensor3& F0_,const STensor3 &F_,const STensor3 &Fp0_,const STensor3 &Fp1_, + const STensor3 &P_, const STensor3 &F_,const STensor3 &Fp0_,const STensor3 &Fp1_, double tildefVstar, double fVstar, const IPNonLocalDamageGurson* ipvprev,const IPNonLocalDamageGurson* ipvcur) const { @@ -971,7 +979,7 @@ void mlawNonLocalDamageGurson::tangent_full_perturbation(STensor43 &T_, ipvp.operator=((const IPVariable &)*ipvcur); Fpp(i,j)+=_perturbationfactor; fVstarp =fVstar; - this->predictorCorrector(tildefVstar, F0_, Fpp, ipvprev->getConstRefToFp(), Fp1p, ipvprev, &ipvp, Pp, fVstarp, false, + this->predictorCorrector(tildefVstar, Fpp, ipvprev->getConstRefToFp(), Fp1p, ipvprev, &ipvp, Pp, fVstarp, false, dDeltaGammadC, dDeltaHatPdC, dtrNpdC, dDeltafVdC, dNpDevdC, dKCordC, dFpdC, dDeltaGammadtildefVstar, dDeltaHatPdtildefVstar, dtrNpdtildefVstar, dDeltafVdtildefVstar, dNpDevdtildefVstar, dKCordtildefVstar, dFpdtildefVstar, @@ -1000,7 +1008,7 @@ void mlawNonLocalDamageGurson::tangent_full_perturbation(STensor43 &T_, Fpp = F_; ipvp.operator=((const IPVariable &)*ipvcur); double tildefVstarp = tildefVstar+_perturbationfactor; - this->predictorCorrector(tildefVstarp,F0_, Fpp, ipvprev->getConstRefToFp(), Fp1p, ipvprev, &ipvp, Pp, fVstarp, false, + this->predictorCorrector(tildefVstarp, Fpp, ipvprev->getConstRefToFp(), Fp1p, ipvprev, &ipvp, Pp, fVstarp, false, dDeltaGammadC, dDeltaHatPdC, dtrNpdC, dDeltafVdC, dNpDevdC, dKCordC, dFpdC, dDeltaGammadtildefVstar, dDeltaHatPdtildefVstar, dtrNpdtildefVstar, dDeltafVdtildefVstar, dNpDevdtildefVstar, dKCordtildefVstar, dFpdtildefVstar, diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h index 0155a0e22..28d7a1da8 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h @@ -94,7 +94,7 @@ class mlawNonLocalDamageGurson : public materialLaw Msg::Info("InitialGuessMethod = %d for 1st NR guess with max %d steps and theta = %f",_initialGuessMethod, _maxite, _theta); // Check for meaningfull values if(_maxite < 10){Msg::Warning("Not enough iteration for NR solving in Gurson law");} - if(_theta<0. or _theta > 1.){Msg::Fatal("Wrong theta value for mid-point parameter");} + if(_theta<0. or _theta > 1.){Msg::Fatal("Wrong theta value for mid-point integration parameter");} }; virtual void setPerzynaViscoPlasticLaw(const int num, const double eta, const double m) { @@ -176,7 +176,7 @@ class mlawNonLocalDamageGurson : public materialLaw protected: virtual double deformationEnergy(const STensor3 &C) const ; - virtual void predictorCorrector(double tildefVstar, const STensor3& F0_, const STensor3& F1_, const STensor3& Fp0_, + virtual void predictorCorrector(double tildefVstar, const STensor3& F1_, const STensor3& Fp0_, STensor3& Fp1_, const IPNonLocalDamageGurson *q0_, IPNonLocalDamageGurson *q1, STensor3&P, double &fVstar, bool stiff, STensor3 &dDeltaGammadC, STensor3 &dDeltaHatPdC, STensor3 &dtrNpdC, @@ -210,7 +210,7 @@ class mlawNonLocalDamageGurson : public materialLaw STensor3 &dLocalCorrectedPorosityDStrain, STensor3 &dStressDNonLocalCorrectedPorosity, double &dLocalCorrectedPorosityDNonLocalCorrectedPorosity, - const STensor3 &P_,const STensor3 &F0_, const STensor3 &F_,const STensor3 &Fp0_,const STensor3 &Fp1_, + const STensor3 &P_, const STensor3 &F_,const STensor3 &Fp0_,const STensor3 &Fp1_, double tildefVstar, double fVstar, const IPNonLocalDamageGurson* q0_,const IPNonLocalDamageGurson* q1_) const; virtual void tangent_full_analytic(STensor43 &T_, -- GitLab From 5756c192689f85aba6c55907f45e4f2d76b05995 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Tue, 22 Aug 2017 18:14:38 +0200 Subject: [PATCH 02/13] checkpoint --- .../materialLaw/mlawNonLocalDamageGurson.cpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index 6de39f4ec..7793b6b02 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -1048,11 +1048,11 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, const STensor3 & Fe1, const STensor3 & kCor, const STensor43 &Le, const STensor63 &dLe) const { - // from d /dC -> d /dF + // from d KCor /dC, d Fp /dC -> d KCor /dF, d Fp /dF static STensor43 dKCordF, dFpdF; - for( int i=0; i<3; i++) - for( int j=0; j<3; j++) - for( int k=0; k<3; k++) + for( int i=0; i<3; i++){ + for( int j=0; j<3; j++){ + for( int k=0; k<3; k++){ for( int l=0; l<3; l++) { dKCordF(i,j,k,l) = 0.; @@ -1063,8 +1063,11 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, dFpdF(i,j,k,l) += dFpdC(i,j,l,m)*F_(k,m)+dFpdC(i,j,m,l)*F_(k,m); } } + } + } + } static STensor3 dDeltafVdF; - for( int i=0; i<3; i++) + for( int i=0; i<3; i++){ for( int j=0; j<3; j++) { dDeltafVdF(i,j) = 0.; @@ -1073,6 +1076,7 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, dDeltafVdF(i,j) += dDeltafVdC(i,m)*F_(j,m)+dDeltafVdC(m,i)*F_(j,m); } } + } // d Fp-1 / d static STensor3 Fpinv; -- GitLab From ea856373a6e50bb9659c7bc2f44f8a7534b42869 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Wed, 23 Aug 2017 16:53:30 +0200 Subject: [PATCH 03/13] optimisation --- .../materialLaw/mlawNonLocalDamageGurson.cpp | 69 ++++++++++++++++--- 1 file changed, 60 insertions(+), 9 deletions(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index 7793b6b02..6c556676e 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -1291,9 +1291,31 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre //FpprinvT = Fpprinv.transpose(); static STensor3 kprDevIdev; - kprDevIdev = kprDev*_Idev; + //kprDevIdev = kprDev*_Idev; + for( int i=0; i<3; i++){ + for( int j=0; j<3; j++){ + kprDevIdev(i,j) = 0.; + for( int k=0; k<3; k++){ + for( int l=0; l<3; l++){ + kprDevIdev(i,j) += kprDev(k,l)*_Idev(l,k,i,j); + } + } + } + } + static STensor3 kprDevIdevL; - kprDevIdevL = kprDevIdev*L; + //kprDevIdevL = kprDevIdev*L; + for( int i=0; i<3; i++){ + for( int j=0; j<3; j++){ + kprDevIdevL(i,j) = 0.; + for( int k=0; k<3; k++){ + for( int l=0; l<3; l++){ + kprDevIdevL(i,j) += kprDevIdev(k,l)*L(l,k,i,j); + } + } + } + } + static STensor3 FpinvkprDevIdevL; static STensor3 FpinvkprDevIdevLFpinvT; @@ -1304,10 +1326,24 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre //FpinvkprDevIdevLFpinvT *= FpprinvT; static STensor3 IL; - IL = _I*L; + //IL = _I*L; + for( int i=0; i<3; i++){ + for( int j=0; j<3; j++){ + IL(i,j) = 0.; + for( int k=0; k<3; k++){ + for( int l=0; l<3; l++){ + IL(i,j) += _I(k,l)*L(l,k,i,j); + } + } + } + } + + + static STensor3 FpinvIL; static STensor3 FpinvILFpinvT; + static STensor3 FpinvILFpinvT_rat; STensorOperation::multSTensor3(Fpprinv,IL,FpinvIL); STensorOperation::multSTensor3(FpinvIL,FpprinvT,FpinvILFpinvT); //FpinvILFpinvT = Fpprinv; @@ -1316,21 +1352,36 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre double rat1= 6.*DeltaGamma*yield*_mu/(yield*yield+6*_mu*DeltaGamma)/(yield*yield+6*_mu*DeltaGamma); double rat2= DeltaGamma*trNp*_K/2./yield; - dres0dC = FpinvkprDevIdevLFpinvT*rat1; - dres0dC += FpinvILFpinvT*rat2; + //dres0dC = FpinvkprDevIdevLFpinvT*rat1; + //dres0dC += FpinvILFpinvT*rat2; + dres0dC = FpinvkprDevIdevLFpinvT; + dres0dC *= rat1; + + FpinvILFpinvT_rat = FpinvILFpinvT; + FpinvILFpinvT_rat *= rat2; + dres0dC += FpinvILFpinvT_rat; double rat3= 3.*yield*yield*_mu/(yield*yield+6*_mu*DeltaGamma)/(yield*yield+6*_mu*DeltaGamma); double rat4= 3.*tildefVstar*_q1*_q2*_K/2./yield*sinh(3.*_q2*(ppr-_K*DeltaGamma*trNp)/2./yield); - dres1dC = FpinvkprDevIdevLFpinvT*rat3; - dres1dC += FpinvILFpinvT*rat4; + //dres1dC = FpinvkprDevIdevLFpinvT*rat3; + //dres1dC += FpinvILFpinvT*rat4; + dres1dC = FpinvkprDevIdevLFpinvT; + dres1dC *= rat3; + FpinvILFpinvT_rat = FpinvILFpinvT; + FpinvILFpinvT_rat *= rat4; + dres1dC += FpinvILFpinvT_rat; - double rat6= 9.*tildefVstar*_q1*_q2*_q2*_K/4./yield/yield*cosh(3.*_q2*(ppr-_K*DeltaGamma*trNp)/2./yield); - dres2dC = FpinvILFpinvT*rat6; + + double rat6= 9.*tildefVstar*_q1*_q2*_q2*_K/4./yield/yield*cosh(3.*_q2*(ppr-_K*DeltaGamma*trNp)/2./yield); + //dres2dC = FpinvILFpinvT*rat6; + dres2dC = FpinvILFpinvT; + dres2dC *= rat6; dres2dC *=yield0; //normalize dres3dC *=0.; + dres0dtildefVstar =0; dres1dtildefVstar =2.*_q1*cosh(3.*_q2*(ppr-_K*DeltaGamma*trNp)/2./yield)-2*_q3*_q3*tildefVstar; -- GitLab From 24f715da0740b529dc2e1d0e083279de9d9c002d Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Wed, 23 Aug 2017 17:11:06 +0200 Subject: [PATCH 04/13] optimisation again --- .../materialLaw/mlawNonLocalDamageGurson.cpp | 49 +++++++++++++++---- 1 file changed, 40 insertions(+), 9 deletions(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index 6c556676e..c21a52bd2 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -1156,9 +1156,10 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, } static STensor3 LedKCordtildefVstar; //LedKCordtildefVstar = Le*dKCordtildefVstar; - LedKCordtildefVstar = 0.; + //LedKCordtildefVstar = 0.; for( int i=0; i<3; i++){ for( int j=0; j<3; j++){ + LedKCordtildefVstar(i,j) = 0.; for( int k=0; k<3; k++){ for( int l=0; l<3; l++){ LedKCordtildefVstar(i,j) += Le(i,j,k,l)*dKCordtildefVstar(l,k); @@ -1227,6 +1228,16 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, } dLedtildefVstarKCor = dLedtildefVstar*kCor; + for( int i=0; i<3; i++){ + for( int j=0; j<3; j++){ + dLedtildefVstarKCor(i,j) = 0.; + for( int k=0; k<3; k++){ + for( int l=0; l<3; l++){ + dLedtildefVstarKCor(i,j) += dLedtildefVstar(i,j,k,l)*kCor(l,k); + } + } + } + } } @@ -1350,17 +1361,18 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre //FpinvILFpinvT *= IL; //FpinvILFpinvT *= FpprinvT; + double rat1= 6.*DeltaGamma*yield*_mu/(yield*yield+6*_mu*DeltaGamma)/(yield*yield+6*_mu*DeltaGamma); double rat2= DeltaGamma*trNp*_K/2./yield; //dres0dC = FpinvkprDevIdevLFpinvT*rat1; //dres0dC += FpinvILFpinvT*rat2; dres0dC = FpinvkprDevIdevLFpinvT; dres0dC *= rat1; - FpinvILFpinvT_rat = FpinvILFpinvT; FpinvILFpinvT_rat *= rat2; dres0dC += FpinvILFpinvT_rat; + double rat3= 3.*yield*yield*_mu/(yield*yield+6*_mu*DeltaGamma)/(yield*yield+6*_mu*DeltaGamma); double rat4= 3.*tildefVstar*_q1*_q2*_K/2./yield*sinh(3.*_q2*(ppr-_K*DeltaGamma*trNp)/2./yield); //dres1dC = FpinvkprDevIdevLFpinvT*rat3; @@ -1372,16 +1384,17 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre dres1dC += FpinvILFpinvT_rat; - double rat6= 9.*tildefVstar*_q1*_q2*_q2*_K/4./yield/yield*cosh(3.*_q2*(ppr-_K*DeltaGamma*trNp)/2./yield); //dres2dC = FpinvILFpinvT*rat6; dres2dC = FpinvILFpinvT; dres2dC *= rat6; - dres2dC *=yield0; //normalize + dres2dC *= yield0; //normalize + dres3dC *=0.; + dres0dtildefVstar =0; dres1dtildefVstar =2.*_q1*cosh(3.*_q2*(ppr-_K*DeltaGamma*trNp)/2./yield)-2*_q3*_q3*tildefVstar; @@ -1498,7 +1511,18 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd dNpDevdtildefVstar*= rat2*dDeltaHatPdtildefVstar+rat3*dDeltaGammadtildefVstar; static STensor3 IL; - IL = _I*Lpr; + //IL = _I*Lpr; + for( int i=0; i<3; i++){ + for( int j=0; j<3; j++){ + IL(i,j) = 0.; + for( int k=0; k<3; k++){ + for( int l=0; l<3; l++){ + IL(i,j) += _I(k,l)*Lpr(l,k,i,j); + } + } + } + } + static STensor3 NpDev; NpDev = kprDev; NpDev*= (3./(yield*yield+6.*_mu*DeltaGamma)); @@ -1514,12 +1538,19 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd for(int n=0; n<3; n++) dKCordC(i,j,k,l)+= (_I(i,j)*_K/2.*_I(m,n)*LFPprinvTFPprinv(m,n,k,l))+_mu*_Idev(i,j,m,n)*LFPprinvTFPprinv(m,n,k,l); } + + //dKCordtildefVstar = _I*(-_K*DeltaGamma*dtrNpdtildefVstar-_K*trNp*dDeltaGammadtildefVstar); + //dKCordtildefVstar += NpDev*(-2.*_mu*dDeltaGammadtildefVstar); + //dKCordtildefVstar += dNpDevdtildefVstar*(-2.*_mu*DeltaGamma); + for( int i=0; i<3; i++){ + for( int j=0; j<3; j++){ + dKCordtildefVstar(i,j) = _I(i,j)*(-_K*DeltaGamma*dtrNpdtildefVstar-_K*trNp*dDeltaGammadtildefVstar); + dKCordtildefVstar(i,j) += NpDev(i,j)*(-2.*_mu*dDeltaGammadtildefVstar); + dKCordtildefVstar(i,j) += dNpDevdtildefVstar(i,j)*(-2.*_mu*DeltaGamma); + } + } - dKCordtildefVstar = _I*(-_K*DeltaGamma*dtrNpdtildefVstar-_K*trNp*dDeltaGammadtildefVstar); - dKCordtildefVstar += NpDev*(-2.*_mu*dDeltaGammadtildefVstar); - dKCordtildefVstar += dNpDevdtildefVstar*(-2.*_mu*DeltaGamma); - static STensor43 EdGammaNpdC; for(int i=0; i<3; i++) -- GitLab From c3671c4838444cedd3573c15d20957e55cbad45d Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Wed, 23 Aug 2017 17:15:32 +0200 Subject: [PATCH 05/13] end optimisation --- NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index c21a52bd2..e5c40fdcb 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -1227,7 +1227,7 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, dLedFKCor(i,j,k,l) += dLedF(i,j,m,n,k,l)*kCor(m,n); } - dLedtildefVstarKCor = dLedtildefVstar*kCor; + //dLedtildefVstarKCor = dLedtildefVstar*kCor; for( int i=0; i<3; i++){ for( int j=0; j<3; j++){ dLedtildefVstarKCor(i,j) = 0.; -- GitLab From f9bf4de73aab852194ca8f2a29dda3d84acfc999 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Thu, 24 Aug 2017 11:16:54 +0200 Subject: [PATCH 06/13] correct mistake in stiffness --- NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index e5c40fdcb..55bd1bad6 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -1073,7 +1073,8 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, dDeltafVdF(i,j) = 0.; for( int m=0; m<3; m++) { - dDeltafVdF(i,j) += dDeltafVdC(i,m)*F_(j,m)+dDeltafVdC(m,i)*F_(j,m); + //dDeltafVdF(i,j) += dDeltafVdC(i,m)*F_(j,m)+dDeltafVdC(m,i)*F_(j,m); + dDeltafVdF(i,j) += dDeltafVdC(j,m)*F_(i,m)+dDeltafVdC(m,j)*F_(i,m); } } } -- GitLab From 393918204bba0c9ce35c02992b6e71f928f5a85e Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Thu, 24 Aug 2017 17:34:57 +0200 Subject: [PATCH 07/13] add variables for scatter in intitial Gurson porosity --- .../materialLaw/mlawNonLocalDamageGurson.cpp | 7 +++++-- .../materialLaw/mlawNonLocalDamageGurson.h | 18 ++++++++++++++---- 2 files changed, 19 insertions(+), 6 deletions(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index 55bd1bad6..cad239856 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -27,7 +27,7 @@ mlawNonLocalDamageGurson::mlawNonLocalDamageGurson(const int num,const double E, _perturbationfactor(pert), _tangentByPerturbation(matrixbyPerturbation ), _q1(q1), _q2(q2), _q3(q3), _fC(fC), _ff(ff), _ffstar(ffstar), - _fVinitial(fVinitial), + _fVinitial(fVinitial), _randMaxbound(1.), _randMinbound(1.), _eta(0.), _m(0.), _I4(1.,1.), _I(1.) { @@ -93,7 +93,8 @@ mlawNonLocalDamageGurson::mlawNonLocalDamageGurson(const mlawNonLocalDamageGurso _tangentByPerturbation(source._tangentByPerturbation), _q1(source._q1), _q2(source._q2), _q3(source._q3), _fC(source._fC), _ff(source._ff), _ffstar(source._ffstar), - _fVinitial(source._fVinitial), + _fVinitial(source._fVinitial), _randMaxbound(source._randMaxbound), + _randMinbound(source._randMinbound), _eta(source._eta), _m(source._m), Cel(source.Cel), _I4(source._I4), _I(source._I), _Idev(source._Idev) @@ -174,6 +175,8 @@ mlawNonLocalDamageGurson& mlawNonLocalDamageGurson::operator = (const materialLa _ff=src->_ff; _ffstar=src->_ffstar; _fVinitial=src->_fVinitial; + _randMaxbound=src->_randMaxbound; + _randMinbound=src->_randMinbound; Cel = src->Cel; _I4 = src->_I4; diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h index 28d7a1da8..2363326f6 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h @@ -51,6 +51,8 @@ class mlawNonLocalDamageGurson : public materialLaw // Initial porosity double _fVinitial; + double _randMaxbound; + double _randMinbound; // Numerical parameters int _order; // log and exp approximation order @@ -84,17 +86,17 @@ class mlawNonLocalDamageGurson : public materialLaw virtual void setOrderForLogExp(const int newOrder) { _order = newOrder; - Msg::Info("order = %d is used to approximate log and exp"); + Msg::Info("Gurson law :: order = %d is used to approximate log and exp"); }; virtual void setPredictorCorrectorParameters(const int InitialGuessMethod, const int maxite=20, const double theta=1.) { _initialGuessMethod = InitialGuessMethod; _maxite = maxite; _theta = theta; - Msg::Info("InitialGuessMethod = %d for 1st NR guess with max %d steps and theta = %f",_initialGuessMethod, _maxite, _theta); + Msg::Info("Gurson law :: InitialGuessMethod = %d for 1st NR guess with max %d steps and theta = %f",_initialGuessMethod, _maxite, _theta); // Check for meaningfull values - if(_maxite < 10){Msg::Warning("Not enough iteration for NR solving in Gurson law");} - if(_theta<0. or _theta > 1.){Msg::Fatal("Wrong theta value for mid-point integration parameter");} + if(_maxite < 10){Msg::Warning("Gurson law :: Not enough iteration for NR solving in Gurson law");} + if(_theta<0. or _theta > 1.){Msg::Fatal("Gurson law ::Wrong theta value for mid-point integration parameter");} }; virtual void setPerzynaViscoPlasticLaw(const int num, const double eta, const double m) { @@ -105,6 +107,14 @@ class mlawNonLocalDamageGurson : public materialLaw { _gdnLawContainer.push_back(added_gdnLaw.clone()); }; + virtual void addScatterredInitialPorosity(double f0min, double f0max) + { + Msg::Info("Gurson law :: Add scatter in initial Gurson porosity"); + _randMinbound = f0min; + if(f0min < 0.){Msg::Fatal("Gurson law ::Wrong scatter parameter fV0_min");} + _randMaxbound = f0max; + if(f0max < f0min){Msg::Fatal("Gurson law ::Wrong scatter parameter fV0_max");} + }; #ifndef SWIG // Constructor-destructor mlawNonLocalDamageGurson(const mlawNonLocalDamageGurson &source); -- GitLab From 982d545757f900a076a6d1da502f3bb0afdf7959 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Tue, 5 Sep 2017 09:37:28 +0200 Subject: [PATCH 08/13] clean unnecessary arguments and computations --- .../materialLaw/mlawNonLocalDamageGurson.cpp | 49 ++++++++++--------- .../materialLaw/mlawNonLocalDamageGurson.h | 35 +++++++------ 2 files changed, 44 insertions(+), 40 deletions(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index cad239856..c2d050296 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -706,8 +706,8 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen { computeDInternalVariables(dDeltaGammadC, dDeltaHatPdC, dtrNpdC, dDeltafVdC, dDeltaGammadtildefVstar, dDeltaHatPdtildefVstar, dtrNpdtildefVstar,dDeltafVdtildefVstar, - tildefVstar, Fp0_, kprDev, ppr, yield, h, yield0, DeltaHatP, DeltaGamma, trNp, - fV0, DeltafV, Lpr, Jinv); + tildefVstar, Fp0_, kprDev, ppr, yield, yield0, DeltaHatP, DeltaGamma, trNp, + Lpr, Jinv); computeDNpDevAndDKCorAndDFp(dNpDevdC, dKCordC,dFpdC,dNpDevdtildefVstar,dKCordtildefVstar,dFpdtildefVstar, dDeltaGammadC,dDeltaHatPdC,dtrNpdC, dDeltafVdC, @@ -1293,18 +1293,21 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dres1dC, STensor3 &dres2dC, STensor3 &dres3dC, double &dres0dtildefVstar, double &dres1dtildefVstar, double &dres2dtildefVstar, double &dres3dtildefVstar, - double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, double ppr, - double yield, double h, double yield0, double DeltaHatP, - double DeltaGamma, double trNp, - double fVn, double DeltafV, const STensor43 &L) const + const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, + const double yield, const double yield0, const double DeltaHatP, + const double DeltaGamma, const double trNp, + const STensor43 &L) const { + /* Compute residual derivatives on strains C */ + + // Compute dKirchoff_pr_dev / dC static STensor3 Fpprinv; STensorOperation::inverseSTensor3(Fppr,Fpprinv); //Fpprinv = Fppr.invert(); static STensor3 FpprinvT; STensorOperation::transposeSTensor3(Fpprinv,FpprinvT); //FpprinvT = Fpprinv.transpose(); - + static STensor3 kprDevIdev; //kprDevIdev = kprDev*_Idev; for( int i=0; i<3; i++){ @@ -1330,8 +1333,7 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre } } } - - + static STensor3 FpinvkprDevIdevL; static STensor3 FpinvkprDevIdevLFpinvT; STensorOperation::multSTensor3(Fpprinv,kprDevIdevL,FpinvkprDevIdevL); @@ -1340,6 +1342,9 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre //FpinvkprDevIdevLFpinvT *= kprDevIdevL; //FpinvkprDevIdevLFpinvT *= FpprinvT; + + + // Compute d(Trace of Kirchoff_pr) / dC static STensor3 IL; //IL = _I*L; for( int i=0; i<3; i++){ @@ -1352,9 +1357,6 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre } } } - - - static STensor3 FpinvIL; static STensor3 FpinvILFpinvT; @@ -1366,6 +1368,7 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre //FpinvILFpinvT *= FpprinvT; + // Compute dres / dC double rat1= 6.*DeltaGamma*yield*_mu/(yield*yield+6*_mu*DeltaGamma)/(yield*yield+6*_mu*DeltaGamma); double rat2= DeltaGamma*trNp*_K/2./yield; //dres0dC = FpinvkprDevIdevLFpinvT*rat1; @@ -1398,7 +1401,7 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre dres3dC *=0.; - + // Compute dres / dtilde_fv dres0dtildefVstar =0; dres1dtildefVstar =2.*_q1*cosh(3.*_q2*(ppr-_K*DeltaGamma*trNp)/2./yield)-2*_q3*_q3*tildefVstar; @@ -1416,18 +1419,17 @@ void mlawNonLocalDamageGurson::computeDInternalVariables( double &dDeltaGammadtildefVstar, double &dDeltaHatPdtildefVstar, double &dtrNpdtildefVstar, double &dDeltafVdtildefVstar, - double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, double ppr, - double yield, double h, double yield0, double DeltaHatP, - double DeltaGamma, double trNp, - double fVn, double DeltafV, const STensor43 &L, - const fullMatrix<double> &scaledJinv) const + const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, + const double yield, const double yield0, const double DeltaHatP, + const double DeltaGamma, const double trNp, + const STensor43 &L, const fullMatrix<double> &scaledJinv) const { static STensor3 dres0dC, dres1dC, dres2dC, dres3dC; double dres0dtildefVstar, dres1dtildefVstar, dres2dtildefVstar, dres3dtildefVstar; computeDResidual(dres0dC,dres1dC,dres2dC, dres3dC, dres0dtildefVstar, dres1dtildefVstar, dres2dtildefVstar, - dres3dtildefVstar, tildefVstar, Fppr, kprDev, ppr, yield, h, yield0, DeltaHatP, - DeltaGamma, trNp, fVn, DeltafV, L); + dres3dtildefVstar, tildefVstar, Fppr, kprDev, ppr, yield, yield0, DeltaHatP, + DeltaGamma, trNp, L); static fullVector<double> residual(4); static fullVector<double> sol(4); @@ -1497,6 +1499,7 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd LFPprinvTFPprinv(i,j,k,l)+= Lpr(i,j,m,n)*Fpprinv(k,m)*Fpprinv(l,n); } + // Compute dNp_dev / dC, dtilde_fv double rat1 = 3.*_mu/(yield*yield+6.*_mu*DeltaGamma); double rat2 = -6.*h*yield/(yield*yield+6.*_mu*DeltaGamma)/(yield*yield+6.*_mu*DeltaGamma); double rat3 = -18.*_mu/(yield*yield+6.*_mu*DeltaGamma)/(yield*yield+6.*_mu*DeltaGamma); @@ -1510,10 +1513,11 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd for(int n=0; n<3; n++) dNpDevdC(i,j,k,l)+= _Idev(i,j,m,n)*LFPprinvTFPprinv(m,n,k,l)*rat1 ; } - + dNpDevdtildefVstar = kprDev; dNpDevdtildefVstar*= rat2*dDeltaHatPdtildefVstar+rat3*dDeltaGammadtildefVstar; + /* static STensor3 IL; //IL = _I*Lpr; for( int i=0; i<3; i++){ @@ -1525,8 +1529,9 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd } } } - } + }*/ + // Compute d Kirchoff / dC static STensor3 NpDev; NpDev = kprDev; NpDev*= (3./(yield*yield+6.*_mu*DeltaGamma)); diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h index 2363326f6..9269707de 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h @@ -77,12 +77,14 @@ class mlawNonLocalDamageGurson : public materialLaw STensor43 _Idev; - public: +public: + // Default constructor mlawNonLocalDamageGurson(const int num,const double E,const double nu, const double rho, const double q1, const double q2, const double q3, const double fC, const double ff, const double ffstar, const double fVinitial, const J2IsotropicHardening &j2IH, const CLengthLaw &cLLaw, const double tol=1.e-8, const bool matrixbyPerturbation = false, const double pert = 1e-8); + // Options setting functions virtual void setOrderForLogExp(const int newOrder) { _order = newOrder; @@ -236,24 +238,21 @@ class mlawNonLocalDamageGurson : public materialLaw const STensor3 & Fe1, const STensor3 & kCor, const STensor43 &Le, const STensor63 &dLe) const; virtual void computeDResidual(STensor3 &dres0dC, STensor3 &dres1dC, STensor3 &dres2dC, STensor3 &dres3dC, - double &dres0dtildefVstar, double &dres1dtildefVstar, double &dres2dtildefVstar, - double &dres3dtildefVstar, - double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, double ppr, - double yield, double h, double yield0, double DeltaHatP, - double DeltaGamma, double trNp, - double fVn, double DeltafV, const STensor43 &L) const; + double &dres0dtildefVstar, double &dres1dtildefVstar, double &dres2dtildefVstar, + double &dres3dtildefVstar, + const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, + const double yield, const double yield0, const double DeltaHatP, + const double DeltaGamma, const double trNp, + const STensor43 &L) const; - virtual void computeDInternalVariables( - STensor3 &dDeltaGammadC, STensor3 &dDeltaHatPdC, STensor3 &dtrNpdC, - STensor3 &dDeltafVdC, - double &dDeltaGammadtildefVstar, double &dDeltaHatPdtildefVstar, - double &dtrNpdtildefVstar, - double &dDeltafVdtildefVstar, - double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, double ppr, - double yield, double h, double yield0, double DeltaHatP, - double DeltaGamma, double trNp, - double fVn, double DeltafV, const STensor43 &L, - const fullMatrix<double> &scaledJinv) const; + virtual void computeDInternalVariables(STensor3 &dDeltaGammadC, STensor3 &dDeltaHatPdC, STensor3 &dtrNpdC, + STensor3 &dDeltafVdC, + double &dDeltaGammadtildefVstar, double &dDeltaHatPdtildefVstar, + double &dtrNpdtildefVstar, double &dDeltafVdtildefVstar, + const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, + const double yield, const double yield0, const double DeltaHatP, + const double DeltaGamma, const double trNp, + const STensor43 &L, const fullMatrix<double> &scaledJinv) const; virtual void computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevdC, STensor43 &dKCordC, STensor43 &dFpdC, STensor3 &dNpDevdtildefVstar, STensor3 &dKCordtildefVstar, STensor3 &dFpdtildefVstar, -- GitLab From e19c82be2e0b1905d63e3538f3c06fc841c5a7be Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Tue, 5 Sep 2017 11:29:49 +0200 Subject: [PATCH 09/13] add debug for stiffness and leaning functions --- .../materialLaw/mlawNonLocalDamageGurson.cpp | 39 +++++++++++++------ .../materialLaw/mlawNonLocalDamageGurson.h | 6 +-- 2 files changed, 29 insertions(+), 16 deletions(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index c2d050296..cbd129b77 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -283,9 +283,27 @@ void mlawNonLocalDamageGurson::constitutive( P,Fn, ipvprev->getConstRefToFp(),Fp1,tildefVstar, fVstar, ipvprev,ipvcur, dKCordC, dFpdC,dDeltafVdC,dKCordtildefVstar, dFpdtildefVstar,dDeltafVdtildefVstar, Fe1, kcor, Le, dLe); - /* this->tangent_full_perturbation(Tangent, dLocalCorrectedPorosityDStrain, - dStressDNonLocalCorrectedPorosity, dLocalCorrectedPorosityDNonLocalCorrectedPorosity, - P,Fn, ipvprev->getConstRefToFp(),Fp1,tildefVstar, fVstar, ipvprev,ipvcur);*/ + + /* + // Debug stiffness + double max_value = 0.; + double max_error = 0.; + static STensor43 Tangent_bis; + this->tangent_full_perturbation(Tangent_bis, dLocalCorrectedPorosityDStrain, + dStressDNonLocalCorrectedPorosity, dLocalCorrectedPorosityDNonLocalCorrectedPorosity, + P, Fn, ipvprev->getConstRefToFp(),Fp1,tildefVstar, fVstar, ipvprev,ipvcur); + for( int i=0; i<3; i++){ + for( int j=0; j<3; j++){ + for( int k=0; k<3; k++){ + for( int l=0; l<3; l++){ + if (fabs(Tangent_bis(i,j,k,l)-Tangent(i,j,k,l)) > max_error){max_error = fabs(Tangent_bis(i,j,k,l)-Tangent(i,j,k,l));} + if (fabs(Tangent(i,j,k,l)) > max_value){max_value = fabs(Tangent_bis(i,j,k,l));} + } + } + } + } + if (max_error/max_value > 0.01) {Msg::Info("Relative error on stiffness: %f", max_error/max_value);} + */ } } @@ -706,7 +724,7 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen { computeDInternalVariables(dDeltaGammadC, dDeltaHatPdC, dtrNpdC, dDeltafVdC, dDeltaGammadtildefVstar, dDeltaHatPdtildefVstar, dtrNpdtildefVstar,dDeltafVdtildefVstar, - tildefVstar, Fp0_, kprDev, ppr, yield, yield0, DeltaHatP, DeltaGamma, trNp, + tildefVstar, Fp0_, kprDev, ppr, yield, yield0, DeltaGamma, trNp, Lpr, Jinv); computeDNpDevAndDKCorAndDFp(dNpDevdC, dKCordC,dFpdC,dNpDevdtildefVstar,dKCordtildefVstar,dFpdtildefVstar, @@ -1294,8 +1312,7 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre double &dres0dtildefVstar, double &dres1dtildefVstar, double &dres2dtildefVstar, double &dres3dtildefVstar, const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, - const double yield, const double yield0, const double DeltaHatP, - const double DeltaGamma, const double trNp, + const double yield, const double yield0, const double DeltaGamma, const double trNp, const STensor43 &L) const { /* Compute residual derivatives on strains C */ @@ -1417,18 +1434,16 @@ void mlawNonLocalDamageGurson::computeDInternalVariables( STensor3 &dDeltaGammadC, STensor3 &dDeltaHatPdC, STensor3 &dtrNpdC, STensor3 &dDeltafVdC, double &dDeltaGammadtildefVstar, double &dDeltaHatPdtildefVstar, - double &dtrNpdtildefVstar, - double &dDeltafVdtildefVstar, + double &dtrNpdtildefVstar, double &dDeltafVdtildefVstar, const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, - const double yield, const double yield0, const double DeltaHatP, - const double DeltaGamma, const double trNp, + const double yield, const double yield0, const double DeltaGamma, const double trNp, const STensor43 &L, const fullMatrix<double> &scaledJinv) const { static STensor3 dres0dC, dres1dC, dres2dC, dres3dC; double dres0dtildefVstar, dres1dtildefVstar, dres2dtildefVstar, dres3dtildefVstar; computeDResidual(dres0dC,dres1dC,dres2dC, dres3dC, dres0dtildefVstar, dres1dtildefVstar, dres2dtildefVstar, - dres3dtildefVstar, tildefVstar, Fppr, kprDev, ppr, yield, yield0, DeltaHatP, + dres3dtildefVstar, tildefVstar, Fppr, kprDev, ppr, yield, yield0, DeltaGamma, trNp, L); static fullVector<double> residual(4); @@ -1531,7 +1546,7 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd } }*/ - // Compute d Kirchoff / dC + // Compute d Kirchoff corr / dC, dtilde fv static STensor3 NpDev; NpDev = kprDev; NpDev*= (3./(yield*yield+6.*_mu*DeltaGamma)); diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h index 9269707de..47595e7a7 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h @@ -241,8 +241,7 @@ public: double &dres0dtildefVstar, double &dres1dtildefVstar, double &dres2dtildefVstar, double &dres3dtildefVstar, const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, - const double yield, const double yield0, const double DeltaHatP, - const double DeltaGamma, const double trNp, + const double yield, const double yield0, const double DeltaGamma, const double trNp, const STensor43 &L) const; virtual void computeDInternalVariables(STensor3 &dDeltaGammadC, STensor3 &dDeltaHatPdC, STensor3 &dtrNpdC, @@ -250,8 +249,7 @@ public: double &dDeltaGammadtildefVstar, double &dDeltaHatPdtildefVstar, double &dtrNpdtildefVstar, double &dDeltafVdtildefVstar, const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, - const double yield, const double yield0, const double DeltaHatP, - const double DeltaGamma, const double trNp, + const double yield, const double yield0, const double DeltaGamma, const double trNp, const STensor43 &L, const fullMatrix<double> &scaledJinv) const; virtual void computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevdC, STensor43 &dKCordC, STensor43 &dFpdC, STensor3 &dNpDevdtildefVstar, -- GitLab From 0318bf7deae5bf2498b282825a77bff65949c870 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Tue, 5 Sep 2017 12:00:37 +0200 Subject: [PATCH 10/13] correct mistake in stiffness + cleaning internal function arguments --- .../materialLaw/mlawNonLocalDamageGurson.cpp | 44 +++++++++++-------- .../materialLaw/mlawNonLocalDamageGurson.h | 22 +++++----- 2 files changed, 37 insertions(+), 29 deletions(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index cbd129b77..4958b2365 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -302,7 +302,7 @@ void mlawNonLocalDamageGurson::constitutive( } } } - if (max_error/max_value > 0.01) {Msg::Info("Relative error on stiffness: %f", max_error/max_value);} + if (max_error/max_value > 0.001) {Msg::Info("Relative error on stiffness: %f", max_error/max_value);} */ } } @@ -1439,6 +1439,7 @@ void mlawNonLocalDamageGurson::computeDInternalVariables( const double yield, const double yield0, const double DeltaGamma, const double trNp, const STensor43 &L, const fullMatrix<double> &scaledJinv) const { + /* Compute internal variables derivatives from residual derivatives */ static STensor3 dres0dC, dres1dC, dres2dC, dres3dC; double dres0dtildefVstar, dres1dtildefVstar, dres2dtildefVstar, dres3dtildefVstar; @@ -1446,6 +1447,7 @@ void mlawNonLocalDamageGurson::computeDInternalVariables( dres3dtildefVstar, tildefVstar, Fppr, kprDev, ppr, yield, yield0, DeltaGamma, trNp, L); + // dC static fullVector<double> residual(4); static fullVector<double> sol(4); for(int i=0; i<3; i++) @@ -1468,6 +1470,8 @@ void mlawNonLocalDamageGurson::computeDInternalVariables( dDeltafVdC(i,j) = sol(3); } } + + // dtilde_fv residual(0) = -dres0dtildefVstar; residual(1) = -dres1dtildefVstar; residual(2) = -dres2dtildefVstar; @@ -1485,19 +1489,20 @@ void mlawNonLocalDamageGurson::computeDInternalVariables( } -void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevdC, STensor43 &dKCordC, STensor43 &dFpdC, +void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp(STensor43 &dNpDevdC, STensor43 &dKCordC, STensor43 &dFpdC, STensor3 &dNpDevdtildefVstar, STensor3 &dKCordtildefVstar, STensor3 &dFpdtildefVstar, const STensor3 &dDeltaGammadC,const STensor3 &dDeltaHatPdC, const STensor3 &dtrNpdC, const STensor3 &dDeltafVdC, - double dDeltaGammadtildefVstar, double dDeltaHatPdtildefVstar, - double dtrNpdtildefVstar, double dDeltafVdtildefVstar, - double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, double ppr, - double yield, double h, double yield0, double DeltaHatP, - double DeltaGamma, double trNp, - double fVn, double DeltafV, const STensor43 &Lpr, const STensor43 &ENp) const + const double dDeltaGammadtildefVstar, const double dDeltaHatPdtildefVstar, + const double dtrNpdtildefVstar, const double dDeltafVdtildefVstar, + const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, + const double yield, const double h, const double yield0, const double DeltaHatP, + const double DeltaGamma, const double trNp, + const double fVn, const double DeltafV, const STensor43 &Lpr, const STensor43 &ENp) const { + /* Compute NpDev, Kirchoff corr and Fp derivatives */ static STensor3 Fpprinv; STensorOperation::inverseSTensor3(Fppr,Fpprinv); //Fpprinv = Fppr.invert(); @@ -1515,9 +1520,10 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd } // Compute dNp_dev / dC, dtilde_fv - double rat1 = 3.*_mu/(yield*yield+6.*_mu*DeltaGamma); - double rat2 = -6.*h*yield/(yield*yield+6.*_mu*DeltaGamma)/(yield*yield+6.*_mu*DeltaGamma); - double rat3 = -18.*_mu/(yield*yield+6.*_mu*DeltaGamma)/(yield*yield+6.*_mu*DeltaGamma); + double rat_den1 = (yield*yield+6.*_mu*DeltaGamma); + double rat1 = 3.*_mu/rat_den1; + double rat2 = -6.*h*yield/rat_den1/rat_den1; + double rat3 = -18.*_mu/rat_den1/rat_den1; for(int i=0; i<3; i++) for(int j=0; j<3; j++) for(int k=0; k<3; k++) @@ -1556,7 +1562,7 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd for(int k=0; k<3; k++) for(int l=0; l<3; l++) { - dKCordC(i,j,k,l) = -_K*DeltaGamma*dtrNpdC(k,l)-_K*trNp*dDeltaGammadC(k,l); + dKCordC(i,j,k,l) = -_I(i,j)*_K*DeltaGamma*dtrNpdC(k,l)-_I(i,j)*_K*trNp*dDeltaGammadC(k,l); dKCordC(i,j,k,l) += -2.*_mu* NpDev(i,j)*dDeltaGammadC(k,l)-2.*_mu*DeltaGamma*dNpDevdC(i,j,k,l); for(int m=0; m<3; m++) for(int n=0; n<3; n++) @@ -1575,7 +1581,7 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd } - + // Compute d Fp / dC, dtildefv static STensor43 EdGammaNpdC; for(int i=0; i<3; i++) for(int j=0; j<3; j++) @@ -1588,6 +1594,8 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd EdGammaNpdC(i,j,k,l)+= ENp(i,j,m,n)*(NpDev(m,n)*dDeltaGammadC(k,l)+DeltaGamma*dNpDevdC(m,n,k,l)+ DeltaGamma/3.*_I(m,n)*dtrNpdC(k,l)); } + + static STensor3 EdGammaNpdftildeVstar; for(int i=0; i<3; i++) for(int j=0; j<3; j++) @@ -1597,7 +1605,7 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd for(int n=0; n<3; n++) EdGammaNpdftildeVstar(i,j)+= ENp(i,j,m,n)*(NpDev(m,n)*dDeltaGammadtildefVstar+DeltaGamma*dNpDevdtildefVstar(m,n)+ DeltaGamma/3.*_I(m,n)*dtrNpdtildefVstar); - } + } for(int i=0; i<3; i++) for(int j=0; j<3; j++) @@ -1606,16 +1614,16 @@ void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevd { dFpdC(j,i,k,l) = 0.; for(int m=0; m<3; m++) - dFpdC(j,i,k,l)+= Fppr(m,i)*EdGammaNpdC(m,j,k,l); + dFpdC(j,i,k,l) += Fppr(m,i)*EdGammaNpdC(m,j,k,l); } for(int i=0; i<3; i++) for(int j=0; j<3; j++) { dFpdtildefVstar(j,i) = 0.; - for(int m=0; m<3; m++) - dFpdtildefVstar(j,i)+= Fppr(m,i)*EdGammaNpdftildeVstar(m,j); - } + for(int m=0; m<3; m++) + dFpdtildefVstar(j,i) += Fppr(m,i)*EdGammaNpdftildeVstar(m,j); + } } diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h index 47595e7a7..a38f2e5d9 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h @@ -251,17 +251,17 @@ public: const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, const double yield, const double yield0, const double DeltaGamma, const double trNp, const STensor43 &L, const fullMatrix<double> &scaledJinv) const; - virtual void computeDNpDevAndDKCorAndDFp( STensor43 &dNpDevdC, STensor43 &dKCordC, STensor43 &dFpdC, - STensor3 &dNpDevdtildefVstar, - STensor3 &dKCordtildefVstar, STensor3 &dFpdtildefVstar, - const STensor3 &dDeltaGammadC,const STensor3 &dDeltaHatPdC, - const STensor3 &dtrNpdC, const STensor3 &dDeltafVdC, - double dDeltaGammadtildefVstar, double dDeltaHatPdtildefVstar, - double dtrNpdtildefVstar, double dDeltafVdtildefVstar, - double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, double ppr, - double yield, double h, double yield0, double DeltaHatP, - double DeltaGamma, double trNp, - double fVn, double DeltafV, const STensor43 &Lpr, const STensor43 &ENp) const; + virtual void computeDNpDevAndDKCorAndDFp(STensor43 &dNpDevdC, STensor43 &dKCordC, STensor43 &dFpdC, + STensor3 &dNpDevdtildefVstar, + STensor3 &dKCordtildefVstar, STensor3 &dFpdtildefVstar, + const STensor3 &dDeltaGammadC,const STensor3 &dDeltaHatPdC, + const STensor3 &dtrNpdC, const STensor3 &dDeltafVdC, + const double dDeltaGammadtildefVstar, const double dDeltaHatPdtildefVstar, + const double dtrNpdtildefVstar, const double dDeltafVdtildefVstar, + const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, + const double yield, const double h, const double yield0, const double DeltaHatP, + const double DeltaGamma, const double trNp, + const double fVn, const double DeltafV, const STensor43 &Lpr, const STensor43 &ENp) const; #endif // SWIG }; -- GitLab From 4bd9388c5bd7c327232ac74f23748ed171c38f40 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Tue, 5 Sep 2017 12:11:41 +0200 Subject: [PATCH 11/13] cleaninbg internal function arguments --- .../materialLaw/mlawNonLocalDamageGurson.cpp | 21 ++++++++----------- .../materialLaw/mlawNonLocalDamageGurson.h | 6 ++---- 2 files changed, 11 insertions(+), 16 deletions(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index 4958b2365..b45931992 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -728,12 +728,11 @@ void mlawNonLocalDamageGurson::predictorCorrector(double tildefVstar, const STen Lpr, Jinv); computeDNpDevAndDKCorAndDFp(dNpDevdC, dKCordC,dFpdC,dNpDevdtildefVstar,dKCordtildefVstar,dFpdtildefVstar, - dDeltaGammadC,dDeltaHatPdC,dtrNpdC, dDeltafVdC, - dDeltaGammadtildefVstar, dDeltaHatPdtildefVstar, - dtrNpdtildefVstar, dDeltafVdtildefVstar, - tildefVstar, Fp0_, kprDev, ppr, - yield, h, yield0, DeltaHatP, - DeltaGamma, trNp, fV0, DeltafV, Lpr, ENp); + dDeltaGammadC,dDeltaHatPdC,dtrNpdC, dDeltafVdC, + dDeltaGammadtildefVstar, dDeltaHatPdtildefVstar, + dtrNpdtildefVstar, dDeltafVdtildefVstar, + Fp0_, kprDev, yield, h, + DeltaGamma, trNp, Lpr, ENp); } else { @@ -1490,16 +1489,14 @@ void mlawNonLocalDamageGurson::computeDInternalVariables( } void mlawNonLocalDamageGurson::computeDNpDevAndDKCorAndDFp(STensor43 &dNpDevdC, STensor43 &dKCordC, STensor43 &dFpdC, - STensor3 &dNpDevdtildefVstar, - STensor3 &dKCordtildefVstar, STensor3 &dFpdtildefVstar, + STensor3 &dNpDevdtildefVstar, + STensor3 &dKCordtildefVstar, STensor3 &dFpdtildefVstar, const STensor3 &dDeltaGammadC,const STensor3 &dDeltaHatPdC, const STensor3 &dtrNpdC, const STensor3 &dDeltafVdC, const double dDeltaGammadtildefVstar, const double dDeltaHatPdtildefVstar, const double dtrNpdtildefVstar, const double dDeltafVdtildefVstar, - const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, - const double yield, const double h, const double yield0, const double DeltaHatP, - const double DeltaGamma, const double trNp, - const double fVn, const double DeltafV, const STensor43 &Lpr, const STensor43 &ENp) const + const STensor3 &Fppr, const STensor3 &kprDev, const double yield, const double h, + const double DeltaGamma, const double trNp,const STensor43 &Lpr, const STensor43 &ENp) const { /* Compute NpDev, Kirchoff corr and Fp derivatives */ diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h index a38f2e5d9..1fc1fd190 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h @@ -258,10 +258,8 @@ public: const STensor3 &dtrNpdC, const STensor3 &dDeltafVdC, const double dDeltaGammadtildefVstar, const double dDeltaHatPdtildefVstar, const double dtrNpdtildefVstar, const double dDeltafVdtildefVstar, - const double tildefVstar, const STensor3 &Fppr, const STensor3 &kprDev, const double ppr, - const double yield, const double h, const double yield0, const double DeltaHatP, - const double DeltaGamma, const double trNp, - const double fVn, const double DeltafV, const STensor43 &Lpr, const STensor43 &ENp) const; + const STensor3 &Fppr, const STensor3 &kprDev, const double yield, const double h, + const double DeltaGamma, const double trNp, const STensor43 &Lpr, const STensor43 &ENp) const; #endif // SWIG }; -- GitLab From 639c67cfaa6be40a2ec00d440e3fd128ded2c17b Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Tue, 5 Sep 2017 16:36:14 +0200 Subject: [PATCH 12/13] correct index mistakes in stiffness --- .../materialLaw/mlawNonLocalDamageGurson.cpp | 27 ++++++++++++++----- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index b45931992..176d02556 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -1161,8 +1161,22 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, } static STensor3 LeKCor; - STensorOperation::multSTensor43STensor3SecondTranspose(Le,kCor,LeKCor); + //STensorOperation::multSTensor43STensor3SecondTranspose(Le,kCor,LeKCor); //LeKCor = Le*kCor; + + for (int i=0; i<3; i++){ + for (int j=0; j<3; j++){ + LeKCor(i,j) = 0.; + for (int k=0; k<3; k++){ + for (int l=0; l<3; l++){ + LeKCor(i,j) += Le(i,j,k,l)*kCor(k,l); + } + } + } + } + + + static STensor43 LedKCordF; for( int i=0; i<3; i++) @@ -1183,7 +1197,7 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, LedKCordtildefVstar(i,j) = 0.; for( int k=0; k<3; k++){ for( int l=0; l<3; l++){ - LedKCordtildefVstar(i,j) += Le(i,j,k,l)*dKCordtildefVstar(l,k); + LedKCordtildefVstar(i,j) += Le(i,j,k,l)*dKCordtildefVstar(k,l); } } } @@ -1191,7 +1205,6 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, // dLe / d - // Here we assume Le constant given the Log approximation static STensor43 dLedFKCor; static STensor3 dLedtildefVstarKCor; @@ -1254,7 +1267,7 @@ void mlawNonLocalDamageGurson::tangent_full_analytic(STensor43 &T_, dLedtildefVstarKCor(i,j) = 0.; for( int k=0; k<3; k++){ for( int l=0; l<3; l++){ - dLedtildefVstarKCor(i,j) += dLedtildefVstar(i,j,k,l)*kCor(l,k); + dLedtildefVstarKCor(i,j) += dLedtildefVstar(i,j,k,l)*kCor(k,l); } } } @@ -1331,7 +1344,7 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre kprDevIdev(i,j) = 0.; for( int k=0; k<3; k++){ for( int l=0; l<3; l++){ - kprDevIdev(i,j) += kprDev(k,l)*_Idev(l,k,i,j); + kprDevIdev(i,j) += kprDev(k,l)*_Idev(k,l,i,j); } } } @@ -1344,7 +1357,7 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre kprDevIdevL(i,j) = 0.; for( int k=0; k<3; k++){ for( int l=0; l<3; l++){ - kprDevIdevL(i,j) += kprDevIdev(k,l)*L(l,k,i,j); + kprDevIdevL(i,j) += kprDevIdev(k,l)*L(k,l,i,j); } } } @@ -1368,7 +1381,7 @@ void mlawNonLocalDamageGurson::computeDResidual(STensor3 &dres0dC, STensor3 &dre IL(i,j) = 0.; for( int k=0; k<3; k++){ for( int l=0; l<3; l++){ - IL(i,j) += _I(k,l)*L(l,k,i,j); + IL(i,j) += _I(k,l)*L(k,l,i,j); } } } -- GitLab From 11ec8f7a6be64baa72dfcb057f0b95e2c7f7b308 Mon Sep 17 00:00:00 2001 From: jleclerc <julien.leclerc@ulg.ac.be> Date: Fri, 8 Sep 2017 15:23:54 +0200 Subject: [PATCH 13/13] add scatter in porosity (first version - work only with cg) --- .../internalPoints/ipNonLocalDamageGurson.h | 7 +++++++ .../materialLaw/mlawNonLocalDamageGurson.cpp | 13 +++++++++---- .../materialLaw/mlawNonLocalDamageGurson.h | 15 ++++++++++----- dG3D/src/FractureCohesiveDG3DMaterialLaw.cpp | 2 ++ dG3D/src/dG3DIPVariable.cpp | 8 ++++++++ dG3D/src/dG3DIPVariable.h | 7 ++++--- dG3D/src/dG3DMaterialLaw.cpp | 10 ++++++---- dG3D/src/dG3DMaterialLaw.h | 3 +++ 8 files changed, 49 insertions(+), 16 deletions(-) diff --git a/NonLinearSolver/internalPoints/ipNonLocalDamageGurson.h b/NonLinearSolver/internalPoints/ipNonLocalDamageGurson.h index a803ebbe0..f6a9669ee 100644 --- a/NonLinearSolver/internalPoints/ipNonLocalDamageGurson.h +++ b/NonLinearSolver/internalPoints/ipNonLocalDamageGurson.h @@ -72,6 +72,13 @@ class IPNonLocalDamageGurson : public IPVariableMechanics virtual void restart(); virtual IPVariable* clone() const {return new IPNonLocalDamageGurson(*this);}; + #if defined(HAVE_MPI) + // value to be commmunicated between processes + virtual int numberValuesToCommunicateMPI()const{return 1;} + virtual void fillValuesToCommunicateMPI(double *arrayMPI)const{arrayMPI[0]=_fV;}; + virtual void getValuesFromMPI(const double *arrayMPI){_fV=arrayMPI[0];}; + #endif // HAVE_MPI + // Damage managing virtual void blockDamage(const bool fl){_damageBlocked = fl;}; diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp index 176d02556..c3e4df279 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp @@ -13,7 +13,6 @@ #include <math.h> #include "MInterfaceElement.h" - mlawNonLocalDamageGurson::mlawNonLocalDamageGurson(const int num,const double E,const double nu, const double rho, const double q1, const double q2, const double q3, const double fC, const double ff, const double ffstar, const double fVinitial, @@ -191,9 +190,10 @@ void mlawNonLocalDamageGurson::createIPState(IPStateBase* &ips,const bool* state //bool inter=true; //const MInterfaceElement *iele = static_cast<const MInterfaceElement*>(ele); //if(iele==NULL) inter=false; - IPVariable* ipvi = new IPNonLocalDamageGurson(_fVinitial,_j2IH,_cLLaw,&_gdnLawContainer); - IPVariable* ipv1 = new IPNonLocalDamageGurson(_fVinitial,_j2IH,_cLLaw,&_gdnLawContainer); - IPVariable* ipv2 = new IPNonLocalDamageGurson(_fVinitial,_j2IH,_cLLaw,&_gdnLawContainer); + double fvInit = getInitialPorosity(); + IPVariable* ipvi = new IPNonLocalDamageGurson(fvInit,_j2IH,_cLLaw,&_gdnLawContainer); + IPVariable* ipv1 = new IPNonLocalDamageGurson(fvInit,_j2IH,_cLLaw,&_gdnLawContainer); + IPVariable* ipv2 = new IPNonLocalDamageGurson(fvInit,_j2IH,_cLLaw,&_gdnLawContainer); if(ips != NULL) delete ips; ips = new IP3State(state_,ipvi,ipv1,ipv2); } @@ -218,6 +218,11 @@ double mlawNonLocalDamageGurson::soundSpeed() const return sqrt(_E*factornu/_rho); } + +double mlawNonLocalDamageGurson::frand(const double a,const double b){ + return (rand()/(double)RAND_MAX) * (b-a) + a; +} + void mlawNonLocalDamageGurson::constitutive(const STensor3& F0,const STensor3& Fn,STensor3 &P,const IPNonLocalDamageGurson *q0, IPNonLocalDamageGurson *q1,STensor43 &Tangent, diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h index 1fc1fd190..a8a88b7b0 100644 --- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h +++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.h @@ -113,11 +113,15 @@ public: { Msg::Info("Gurson law :: Add scatter in initial Gurson porosity"); _randMinbound = f0min; - if(f0min < 0.){Msg::Fatal("Gurson law ::Wrong scatter parameter fV0_min");} + if(f0min < 0.){Msg::Fatal("Gurson law :: Wrong scatter parameter fV0_min");} _randMaxbound = f0max; - if(f0max < f0min){Msg::Fatal("Gurson law ::Wrong scatter parameter fV0_max");} + if(f0max < f0min){Msg::Fatal("Gurson law :: Wrong scatter parameter fV0_max");} }; + #ifndef SWIG + static double frand(const double a,const double b); + + // Constructor-destructor mlawNonLocalDamageGurson(const mlawNonLocalDamageGurson &source); mlawNonLocalDamageGurson& operator=(const materialLaw &source); @@ -156,10 +160,11 @@ public: virtual const double bulkModulus() const{return _K;} virtual const double shearModulus() const{return _mu;} virtual const double poissonRatio() const{return _nu;} - virtual const double getInitialPorosity() const {return _fVinitial;} + virtual const double getInitialPorosity() const {return _fVinitial*frand(_randMinbound,_randMaxbound);} // for Ipv creation only! // Specific functions - public: +public: + virtual void constitutive( const STensor3& F0, // initial deformation gradient (input @ time n) const STensor3& Fn, // updated deformation gradient (input @ time n+1) @@ -184,7 +189,7 @@ public: double &dLocalCorrectedPorosityDNonLocalCorrectedPorosity, const bool stiff // if true compute the tangents ) const; - + protected: virtual double deformationEnergy(const STensor3 &C) const ; diff --git a/dG3D/src/FractureCohesiveDG3DMaterialLaw.cpp b/dG3D/src/FractureCohesiveDG3DMaterialLaw.cpp index 227b7be53..2ad785f0e 100644 --- a/dG3D/src/FractureCohesiveDG3DMaterialLaw.cpp +++ b/dG3D/src/FractureCohesiveDG3DMaterialLaw.cpp @@ -57,6 +57,8 @@ void FractureByCohesive3DLaw::createIPState(IPStateBase* &ips,const bool* state_ ipv=NULL; _mbulk->createIPVariable(ipv,ele,nbFF_,GP,gpt); ipvf->setIPvBulk(ipv); + //to have the same random variables at the three states + // Why cohesive IPvariable are created here ? ipvf = static_cast<FractureCohesive3DIPVariable*>(ipvi); diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp index 333514160..3e6037bf5 100644 --- a/dG3D/src/dG3DIPVariable.cpp +++ b/dG3D/src/dG3DIPVariable.cpp @@ -1176,6 +1176,14 @@ void nonLocalDamageJ2HyperDG3DIPVariable::restart() return; } // + +nonLocalDamageGursonDG3DIPVariable::nonLocalDamageGursonDG3DIPVariable(const mlawNonLocalDamageGurson &_gursonlaw,const double fvInitial, const bool oninter): + nonLocalDamageDG3DIPVariableBase(oninter,1) +{ + _nldGursonipv = new IPNonLocalDamageGurson(fvInitial,_gursonlaw.getJ2IsotropicHardening(),_gursonlaw.getCLengthLaw(), + _gursonlaw.getGursonDamageNucleationContainer()); +} + nonLocalDamageGursonDG3DIPVariable::nonLocalDamageGursonDG3DIPVariable(const mlawNonLocalDamageGurson &_law, const bool oninter) : nonLocalDamageDG3DIPVariableBase(oninter,1) { diff --git a/dG3D/src/dG3DIPVariable.h b/dG3D/src/dG3DIPVariable.h index 8ec6efd6f..a853589c7 100644 --- a/dG3D/src/dG3DIPVariable.h +++ b/dG3D/src/dG3DIPVariable.h @@ -848,12 +848,11 @@ class nonLocalDamageJ2HyperDG3DIPVariable : public nonLocalDamageDG3DIPVariableB class nonLocalDamageGursonDG3DIPVariable : public nonLocalDamageDG3DIPVariableBase { - protected: IPNonLocalDamageGurson *_nldGursonipv; public: - nonLocalDamageGursonDG3DIPVariable(const mlawNonLocalDamageGurson &_gursonlaw, const bool oninter=false); + nonLocalDamageGursonDG3DIPVariable(const mlawNonLocalDamageGurson &_gursonlaw,const double fvInitial, const bool oninter=false); nonLocalDamageGursonDG3DIPVariable(const nonLocalDamageGursonDG3DIPVariable &source); virtual nonLocalDamageGursonDG3DIPVariable& operator=(const IPVariable &source); @@ -904,7 +903,9 @@ class nonLocalDamageGursonDG3DIPVariable : public nonLocalDamageDG3DIPVariableBa virtual void restart(); - +private: + nonLocalDamageGursonDG3DIPVariable(const mlawNonLocalDamageGurson &_gursonlaw, const bool oninter=false); + }; class TransverseIsotropicDG3DIPVariable : public dG3DIPVariable // or store data in a different way diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp index 030165bbc..90e586909 100644 --- a/dG3D/src/dG3DMaterialLaw.cpp +++ b/dG3D/src/dG3DMaterialLaw.cpp @@ -1515,9 +1515,10 @@ void NonLocalDamageGursonDG3DMaterialLaw::createIPState(IPStateBase* &ips,const bool inter=true; const MInterfaceElement *iele = dynamic_cast<const MInterfaceElement*>(ele); if(iele==NULL) inter=false; - IPVariable* ipvi = new nonLocalDamageGursonDG3DIPVariable(*_nldGursonlaw,inter); - IPVariable* ipv1 = new nonLocalDamageGursonDG3DIPVariable(*_nldGursonlaw,inter); - IPVariable* ipv2 = new nonLocalDamageGursonDG3DIPVariable(*_nldGursonlaw,inter); + double fvInit = _nldGursonlaw->getInitialPorosity(); + IPVariable* ipvi = new nonLocalDamageGursonDG3DIPVariable(*_nldGursonlaw,fvInit,inter); + IPVariable* ipv1 = new nonLocalDamageGursonDG3DIPVariable(*_nldGursonlaw,fvInit,inter); + IPVariable* ipv2 = new nonLocalDamageGursonDG3DIPVariable(*_nldGursonlaw,fvInit,inter); if(ips != NULL) delete ips; ips = new IP3State(state_,ipvi,ipv1,ipv2); _nldGursonlaw->createIPState((static_cast <nonLocalDamageGursonDG3DIPVariable*> (ipvi))->getIPNonLocalDamageGurson(), @@ -1534,7 +1535,8 @@ void NonLocalDamageGursonDG3DMaterialLaw::createIPVariable(IPVariable* &ipv,cons if(iele == NULL){ inter=false; } - ipv = new nonLocalDamageGursonDG3DIPVariable(*_nldGursonlaw,inter); + double fvInit = _nldGursonlaw->getInitialPorosity(); + ipv = new nonLocalDamageGursonDG3DIPVariable(*_nldGursonlaw,fvInit,inter); IPNonLocalDamageGurson * ipvnl = static_cast <nonLocalDamageGursonDG3DIPVariable*>(ipv)->getIPNonLocalDamageGurson(); _nldGursonlaw->createIPVariable(ipvnl,ele,nbFF_,GP,gpt); } diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h index b105dd2b9..b3b8e9e84 100644 --- a/dG3D/src/dG3DMaterialLaw.h +++ b/dG3D/src/dG3DMaterialLaw.h @@ -948,6 +948,9 @@ class NonLocalDamageGursonDG3DMaterialLaw : public dG3DMaterialLaw {_nldGursonlaw->setPerzynaViscoPlasticLaw(num,eta,m);}; virtual void addNucleationLaw(const GursonDamageNucleation& added_GdnLaw) {_nldGursonlaw->addNucleationLaw(added_GdnLaw);}; + virtual void addScatterredInitialPorosity(const double fmin, const double fmax) + {_nldGursonlaw->addScatterredInitialPorosity(fmin,fmax);}; + #ifndef SWIG NonLocalDamageGursonDG3DMaterialLaw(const NonLocalDamageGursonDG3DMaterialLaw &source); // set the time of _nldlaw -- GitLab