From 7547aee6666f6d36325a39ceac56269599ec6754 Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Wed, 13 Jul 2016 15:28:02 +0200
Subject: [PATCH] Added: Input tag <noEnergySplit/> for isotropic strain energy
 degrading. Added: Treat phase field values larger than 1.0 as equal to 1.0.

---
 FractureElasticity.C      | 10 +++++++
 FractureElasticity.h      |  2 ++
 FractureElasticityVoigt.C | 61 +++++++++++++++++++++++++++------------
 3 files changed, 55 insertions(+), 18 deletions(-)

diff --git a/FractureElasticity.C b/FractureElasticity.C
index 7bfc71c..1004aeb 100644
--- a/FractureElasticity.C
+++ b/FractureElasticity.C
@@ -32,6 +32,7 @@
 FractureElasticity::FractureElasticity (unsigned short int n)
   : Elasticity(n), mySol(primsol)
 {
+  noSplit = false;
   alpha = 0.0;
   this->registerVector("phasefield",&myCVec);
   eC = 1; // Assuming second vector is phase field
@@ -42,6 +43,7 @@ FractureElasticity::FractureElasticity (IntegrandBase* parent,
                                         unsigned short int n)
   : Elasticity(n), mySol(parent->getSolutions())
 {
+  noSplit = false;
   alpha = 0.0;
   parent->registerVector("phasefield",&myCVec);
   // Assuming second vector is pressure, third vector is pressure velocity
@@ -62,6 +64,8 @@ bool FractureElasticity::parse (const TiXmlElement* elem)
   const char* value = utl::getValue(elem,"stabilization");
   if (value)
     alpha = atof(value);
+  else if (!strcasecmp(elem->Value(),"noEnergySplit"))
+    noSplit = true;
   else
     return this->Elasticity::parse(elem);
 
@@ -75,6 +79,11 @@ void FractureElasticity::printLog () const
 
   if (alpha != 0.0)
     IFEM::cout <<"\tStabilization parameter: "<< alpha << std::endl;
+
+  if (noSplit)
+    IFEM::cout <<"\tIsotropic degrading of strain energy density."<< std::endl;
+  else
+    IFEM::cout <<"\tDegrading of tensile strain energy density."<< std::endl;
 }
 
 
@@ -269,6 +278,7 @@ double FractureElasticity::getStressDegradation (const Vector& N,
 
   // Evaluate the crack phase field function, c(X)
   double c = eV[eC].empty() ? 1.0 : N.dot(eV[eC]);
+  if (c > 1.0) c = 1.0; // Ignore values larger than 1.0
 
   if (derivative > 0) // Evaluate the 1st derivative of g(c)
     return c > 0.0 ? (1.0-alpha)*c*2.0 : 0.0;
diff --git a/FractureElasticity.h b/FractureElasticity.h
index d85d9ab..d55e63e 100644
--- a/FractureElasticity.h
+++ b/FractureElasticity.h
@@ -121,6 +121,8 @@ class FractureElasticity : public Elasticity
   Vector myCVec; //!< Crack phase field values at control (nodal) points
 
 protected:
+  bool noSplit; //!< If \e true, no strain energy split, just isotropic scaling
+
   unsigned short int eC; //!< Zero-based index to element phase field vector
 
   mutable RealArray myPhi; //!< Tensile energy density at integration points
diff --git a/FractureElasticityVoigt.C b/FractureElasticityVoigt.C
index 39ef3d5..5ea638e 100644
--- a/FractureElasticityVoigt.C
+++ b/FractureElasticityVoigt.C
@@ -262,22 +262,35 @@ bool FractureElasticityVoigt::evalInt (LocalIntegral& elmInt,
     std::cout <<"\nFractureElasticity::evalInt(X = "<< X <<")\nBmat ="<< Bmat;
 #endif
 
-    // Evaluate the material parameters at this point
-    double lambda, mu;
-    if (!material->evaluate(lambda,mu,fe,X))
-      return false;
-
     // Evaluate the stress degradation function
     double Gc = this->getStressDegradation(fe.N,elmInt.vec);
+    if (noSplit)
+    {
+      // Evaluate the constitutive matrix and the stress tensor at this point
+      if (!material->evaluate(dSdE,sigma,myPhi[fe.iGP],fe,X,eps,eps))
+	return false;
+
+      // Degrade the stresses and strain energy isotropically
+      myPhi[fe.iGP] *= Gc;
+      sigma *= Gc;
+    }
+    else
+    {
+      // Evaluate the material parameters at this point
+      double lambda, mu;
+      if (!material->evaluate(lambda,mu,fe,X))
+        return false;
+
 #if INT_DEBUG > 3
-    std::cout <<"lambda = "<< lambda <<" mu = "<< mu <<" G(c) = "<< Gc <<"\n";
-    if (lHaveStrains) std::cout <<"eps =\n"<< eps;
+      std::cout <<"lambda = "<< lambda <<" mu = "<< mu <<" G(c) = "<< Gc <<"\n";
+      if (lHaveStrains) std::cout <<"eps =\n"<< eps;
 #endif
 
-    // Evaluate the stress state at this point
-    if (!this->evalStress(lambda,mu,Gc,eps,&myPhi[fe.iGP],&sigma,
-                          eKm ? &dSdE : nullptr))
-      return false;
+      // Evaluate the stress state at this point, with degraded tensile part
+      if (!this->evalStress(lambda,mu,Gc,eps,&myPhi[fe.iGP],&sigma,
+                            eKm ? &dSdE : nullptr))
+	return false;
+    }
   }
 
   if (eKm)
@@ -396,11 +409,6 @@ bool FractureElasticNorm::evalInt (LocalIntegral& elmInt,
       for (unsigned short int j = i+1; j <= eps.dim(); j++)
         eps(i,j) *= 0.5;
 
-  // Evaluate the material parameters at this point
-  double lambda, mu;
-  if (!p.material->evaluate(lambda,mu,fe,X))
-    return false;
-
   bool printElm = fe.iel == dbgElm;
   if (printElm)
     std::cout <<"\nFractureElasticNorm::evalInt: iel,ip,X = "
@@ -409,8 +417,25 @@ bool FractureElasticNorm::evalInt (LocalIntegral& elmInt,
   // Evaluate the strain energy at this point
   double Phi[4];
   double Gc = p.getStressDegradation(fe.N,elmInt.vec);
-  if (!p.evalStress(lambda,mu,Gc,eps,Phi,nullptr,nullptr,true,printElm))
-    return false;
+  if (p.noSplit)
+  {
+    // Evaluate the strain energy density at this point
+    SymmTensor sigma(eps.dim());
+    if (!p.material->evaluate(Bmat,sigma,Phi[2],fe,X,eps,eps,3))
+      return false;
+    Phi[2] *= Gc; // Isotropic scaling
+    Phi[0] = Phi[1] = Phi[3] = 0.0;
+  }
+  else
+  {
+    // Evaluate the material parameters at this point
+    double lambda, mu;
+    if (!p.material->evaluate(lambda,mu,fe,X))
+      return false;
+    // Evaluate the tensile-degraded strain energy
+    if (!p.evalStress(lambda,mu,Gc,eps,Phi,nullptr,nullptr,true,printElm))
+      return false;
+  }
 
   // Integrate the total elastic energy
   pnorm[0] += Phi[2]*fe.detJxW;