squash! Added: Tangent matrix for monolithic coupling of elasticity a…

…nd phasefield. Added: Option to return derivatives of the stress degradation function. Added: Parsing and printing of the stabilization parameter alpha, this is the same quantity as stabk in the CahnHilliard integrand.

The monolithic coupling matrix was incorrect

FractureElasticityVoigt.C

@@ -28,8 +28,9 @@
 
 bool FractureElasticityVoigt::evalStress (double lambda, double mu, double Gc,
                                           const SymmTensor& epsil, double* Phi,
-                                          SymmTensor* sigma, Matrix* dSdE,
-                                          bool postProc, bool printElm) const
+                                          SymmTensor* sigma, SymmTensor* sigP,
+                                          Matrix* dSdE, bool postProc,
+                                          bool printElm) const
 {
   PROFILE3("FractureEl::evalStress");
 
@@ -90,6 +91,12 @@ bool FractureElasticityVoigt::evalStress (double lambda, double mu, double Gc,
     *sigma = C0*trEps;
     *sigma += 2.0*mu*(Gc*ePos + eNeg);
   }
+  if (sigP)
+  {
+    // Evaluate the non-degraded tensile stress tensor
+    *sigP = lambda*trEps;
+    *sigP += 2.0*mu*ePos;
+  }
 
   // Evaluate the tensile energy
   Phi[0] = mu*(ePos*ePos).trace();
@@ -113,6 +120,7 @@ bool FractureElasticityVoigt::evalStress (double lambda, double mu, double Gc,
     std::cout <<"M("<< 1+a <<") =\n"<< M[a];
   std::cout <<"ePos =\n"<< ePos <<"eNeg =\n"<< eNeg;
   if (sigma) std::cout <<"sigma =\n"<< *sigma;
+  if (sigP) std::cout <<"sigP =\n"<< *sigP;
   std::cout <<"Phi = "<< Phi[0];
   if (postProc) std::cout <<" "<< Phi[1] <<" "<< Phi[2] <<" "<< Phi[3];
   std::cout << std::endl;
@@ -242,7 +250,7 @@ bool FractureElasticityVoigt::evalInt (LocalIntegral& elmInt,
 
   size_t nstrc = (nsd+1)*nsd/2;
   Matrix Bmat, dSdE(nstrc,nstrc);
-  SymmTensor eps(nsd), sigma(nsd);
+  SymmTensor eps(nsd), sigma(nsd), sigP(nsd);
   bool lHaveStrains = false;
 
   if (eKm || eKg || iS || m_mode == SIM::RECOVERY)
@@ -269,12 +277,12 @@ bool FractureElasticityVoigt::evalInt (LocalIntegral& elmInt,
     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))
+      if (!material->evaluate(dSdE,sigP,myPhi[fe.iGP],fe,X,eps,eps))
         return false;
 
       // Degrade the stresses and strain energy isotropically
       myPhi[fe.iGP] *= Gc;
-      sigma *= Gc;
+      sigma = Gc*sigP;
     }
     else
     {
@@ -290,7 +298,7 @@ bool FractureElasticityVoigt::evalInt (LocalIntegral& elmInt,
 
       // 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))
+                            eKc ? &sigP : nullptr, eKm ? &dSdE : nullptr))
         return false;
     }
   }
@@ -312,6 +320,14 @@ bool FractureElasticityVoigt::evalInt (LocalIntegral& elmInt,
   if (eM) // Integrate the mass matrix
     this->formMassMatrix(elMat.A[eM-1],fe.N,X,fe.detJxW);
 
+  if (iS && lHaveStrains)
+  {
+    // Integrate the internal forces
+    sigma *= -fe.detJxW;
+    if (!Bmat.multiply(sigma,elMat.b[iS-1],true,true)) // ES -= B^T*sigma
+      return false;
+  }
+
   if (eS)
   {
     // Integrate the load vector due to gravitation and other body forces
@@ -320,28 +336,14 @@ bool FractureElasticityVoigt::evalInt (LocalIntegral& elmInt,
     this->formCrackForce(elMat.b[eS-1],elMat.vec,fe,X);
   }
 
-  if (!lHaveStrains)
-    return true;
-
-  Vector BtSig;
-  sigma *= -fe.detJxW;
-  if (eKc && !Bmat.multiply(sigma,BtSig,true)) // BtSig = -B^T*sigma
-    return false;
-
-  if (iS)
-  {
-    // Integrate the internal forces
-    if (eKc)
-      elMat.b[iS-1].add(BtSig); // ES += BtSig
-    else if (!Bmat.multiply(sigma,elMat.b[iS-1],true,true)) // ES -= B^T*sigma
-      return false;
-  }
-
-  if (eKc)
+  if (eKc && lHaveStrains)
   {
     // Integrate the coupling matrix between displacement and phase-field DOFs
-    BtSig *= -this->getStressDegradation(fe.N,elmInt.vec,1); // g'(c)
-    if (!elMat.A[eKc-1].outer_product(BtSig,fe.N,true)) // Kc += B^T*sigma*N*dGc
+    Vector BtSig;
+    sigP *= this->getStressDegradation(fe.N,elmInt.vec,1)*fe.detJxW; // *= g'(c)
+    if (!Bmat.multiply(sigP,BtSig,true)) // BtSig = B^T*g'(c)*sigP
+      return false;
+    if (!elMat.A[eKc-1].outer_product(BtSig,fe.N,true)) // Kc += BtSig*N
       return false;
   }
 
@@ -353,7 +355,7 @@ bool FractureElasticityVoigt::evalStress (double lambda, double mu, double Gc,
                                           const SymmTensor& epsilon,
                                           double* Phi, SymmTensor& sigma) const
 {
-  return this->evalStress(lambda,mu,Gc,epsilon,Phi,&sigma,nullptr,true);
+  return this->evalStress(lambda,mu,Gc,epsilon,Phi,&sigma,nullptr,nullptr,true);
 }
 
 
@@ -445,7 +447,7 @@ bool FractureElasticNorm::evalInt (LocalIntegral& elmInt,
     // Evaluate the tensile-degraded strain energy
     if (!p.evalStress(lambda,mu,Gc,eps,Phi,
                       dirDer ? &sigma : nullptr,
-                      nullptr,true,printElm))
+                      nullptr,nullptr,true,printElm))
       return false;
   }
 

FractureElasticityVoigt.h

@@ -62,7 +62,7 @@ class FractureElasticityVoigt : public FractureElasticity
   //! \brief Evaluates the stress tensor and its derivative w.r.t. the strains.
   bool evalStress(double lambda, double mu, double Gc,
                   const SymmTensor& epsilon, double* Phi,
-                  SymmTensor* sigma, Matrix* dSdE,
+                  SymmTensor* sigma, SymmTensor* sigP, Matrix* dSdE,
                   bool postProc = false, bool printElm = false) const;
 
   unsigned short int eKc; //!< Zero-based index to element coupling matrix

Test/TestFractureElasticity.C

@@ -29,7 +29,7 @@ public:
   virtual ~FracEl() {}
   bool calcStress(double lambda, double mu, double Gc, const SymmTensor& eps,
                   double& Phi, SymmTensor& sigma, Matrix& dSdE) const
-  { return evalStress(lambda,mu,Gc,eps,&Phi,&sigma,&dSdE); }
+  { return evalStress(lambda,mu,Gc,eps,&Phi,&sigma,nullptr,&dSdE); }
   bool calcStress(double lambda, double mu, double Gc, const SymmTensor& eps,
                   double& Phi, SymmTensor& sigma, Tensor4& dSdE) const
   { return FractureElasticity::evalStress(lambda,mu,Gc,eps,&Phi,sigma,&dSdE); }