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