squash! Added: Calculation of f'(alpha) by direct integration togethe…

…r with f(alpha)

There was also an arror in the calculation of f(alpha)
(missing Gc-factor in front of the gradient term)

FractureElasticityMonol.C

@@ -340,7 +340,8 @@ bool FractureElasticMoNorm::evalInt (LocalIntegral& elmInt,
 
   double C = fe.N.dot(elmInt.vec[p.eC]);
   double d = 1.0 - C;
-  double scale = 0.25*p.Gc*d*d/p.smearing + p.smearing*gradC.dot(gradC);
+  double ell = 2.0*p.smearing;
+  double scale = 0.5*p.Gc*(gradC.dot(gradC)*ell + d*d/ell);
   if (C < p.crtol) scale += 0.5*p.gammaInv*C*C;
 
   // Integrate the dissipated energy
@@ -353,7 +354,7 @@ bool FractureElasticMoNorm::evalInt (LocalIntegral& elmInt,
       return false;
 
     double Cc = fe.N.dot(elmInt.vec[p.eC+1]);
-    scale = -0.5*p.Gc*d*Cc/p.smearing + p.smearing*gradC.dot(gradCc);
+    scale = p.Gc*(gradC.dot(gradCc)*ell - d*Cc/ell);
     if (C < p.crtol) scale += p.gammaInv*C*Cc;
 
     // Integrate directional derivative of the dissipated energy

FractureElasticityVoigt.C

@@ -471,8 +471,9 @@ bool FractureElasticNorm::evalInt (LocalIntegral& elmInt,
   {
     // Integrate directional derivative of the elastic energy
     Bmat.multiply(elmInt.vec[1],eps);
+    double dGc = p.getStressDegradation(fe.N,elmInt.vec,1);
     double Cc = fe.N.dot(elmInt.vec[p.eC+1]);
-    pnorm[4] += (sigma.innerProd(eps) + 2.0*Phi[0]*Cc)*fe.detJxW;
+    pnorm[4] += (sigma.innerProd(eps) + dGc*Phi[0]*Cc)*fe.detJxW;
   }
   else
     // Integrate the bulk energy