Added: Calculation of f'(alpha) by direct integration together with f…

…(alpha)

FractureElasticityMonol.C

@@ -27,6 +27,7 @@ FractureElasticityMonol::FractureElasticityMonol (unsigned short int n, int ord)
   crtol = 0.0;
   use4th = ord == 4;
   npv = nsd + 1;
+  eC = 2; // Assuming third vector is phase field
 }
 
 
@@ -87,9 +88,11 @@ void FractureElasticityMonol::setMode (SIM::SolutionMode mode)
     case SIM::INT_FORCES:
       iS  = 2;
       eBc = 3;
+      primsol.resize(1);
+      break;
 
     case SIM::RECOVERY:
-      primsol.resize(1);
+      primsol.resize(FractureElasticNorm::dirDer ? 2 : 1);
       break;
 
     default:
@@ -151,6 +154,8 @@ bool FractureElasticityMonol::initElement (const std::vector<int>& MNPC,
   }
 
   size_t nsol = 1 + eC;
+  if (m_mode == SIM::RECOVERY && FractureElasticNorm::dirDer)
+    nsol++;
   if (elmInt.vec.size() < nsol)
     elmInt.vec.resize(nsol);
 
@@ -255,6 +260,26 @@ bool FractureElasticityMonol::getSolution (Vectors& eV,
   std::cout <<"Element displacement vector:"<< eV.front()
             <<"Element phase field vector:"<< eV[eC];
 #endif
+
+  if (m_mode == SIM::RECOVERY && FractureElasticNorm::dirDer)
+  {
+    // Extract element solution correction vectors
+    ierr = utl::gather(MNPC,npv,mySol.back(),temp);
+    if (ierr > 0)
+    {
+      std::cerr <<" *** FractureElasticityMonol::getSolution: Detected "
+                << ierr <<" node numbers out of range."<< std::endl;
+      return false;
+    }
+
+    eV[eC+1] = temp.getRow(npv);
+    eV[1]    = temp.expandRows(-1);
+#if INT_DEBUG > 2
+    std::cout <<"Element displacement correction vector:"<< eV[1]
+              <<"Element phase field correction vector:"<< eV[eC+1];
+#endif
+  }
+
   return true;
 }
 
@@ -321,5 +346,19 @@ bool FractureElasticMoNorm::evalInt (LocalIntegral& elmInt,
   // Integrate the dissipated energy
   pnorm[5] += scale*fe.detJxW;
 
+  if (dirDer)
+  {
+    Vector gradCc; // Compute the gradient of the phase field correction
+    if (!fe.dNdX.multiply(elmInt.vec[p.eC+1],gradCc,true))
+      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);
+    if (C < p.crtol) scale += p.gammaInv*C*Cc;
+
+    // Integrate directional derivative of the dissipated energy
+    pnorm[4] += scale*fe.detJxW;
+  }
+
   return true;
 }

FractureElasticityVoigt.C

@@ -363,7 +363,9 @@ NormBase* FractureElasticityVoigt::getNormIntegrand (AnaSol*) const
 }
 
 
-int FractureElasticNorm::dbgElm = 0;
+bool FractureElasticNorm::extEnr = true;
+bool FractureElasticNorm::dirDer = false;
+int  FractureElasticNorm::dbgElm = 0;
 
 
 FractureElasticNorm::FractureElasticNorm (FractureElasticityVoigt& p)
@@ -424,11 +426,12 @@ bool FractureElasticNorm::evalInt (LocalIntegral& elmInt,
   // Evaluate the strain energy at this point
   double Phi[4];
   double Gc = p.getStressDegradation(fe.N,elmInt.vec);
+  SymmTensor sigma(eps.dim());
   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))
+    Matrix Cmat;
+    if (!p.material->evaluate(Cmat,sigma,Phi[2],fe,X,eps,eps,3))
       return false;
     Phi[2] *= Gc; // Isotropic scaling
     Phi[0] = Phi[1] = Phi[3] = 0.0;
@@ -440,14 +443,16 @@ bool FractureElasticNorm::evalInt (LocalIntegral& elmInt,
     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))
+    if (!p.evalStress(lambda,mu,Gc,eps,Phi,
+                      dirDer ? &sigma : nullptr,
+                      nullptr,true,printElm))
       return false;
   }
 
   // Integrate the total elastic energy
   pnorm[0] += Phi[2]*fe.detJxW;
 
-  if (p.haveLoads())
+  if (extEnr && p.haveLoads())
   {
     // Evaluate the body load
     Vec3 f = p.getBodyforce(X);
@@ -460,8 +465,16 @@ bool FractureElasticNorm::evalInt (LocalIntegral& elmInt,
   // Integrate the tensile and compressive energies
   pnorm[2] += Phi[0]*fe.detJxW;
   pnorm[3] += Phi[1]*fe.detJxW;
-  // Integrate the bulk energy
-  pnorm[4] += Phi[3]*fe.detJxW;
+  if (dirDer)
+  {
+    // Integrate directional derivative of the elastic energy
+    Bmat.multiply(elmInt.vec[1],eps);
+    double Cc = fe.N.dot(elmInt.vec[p.eC+1]);
+    pnorm[4] += (sigma.innerProd(eps) + 2.0*Phi[0]*Cc)*fe.detJxW;
+  }
+  else
+    // Integrate the bulk energy
+    pnorm[4] += Phi[3]*fe.detJxW;
 
   return true;
 }

FractureElasticityVoigt.h

@@ -97,7 +97,9 @@ class FractureElasticNorm : public ElasticityNorm
   //! \brief Returns the name of a norm quantity.
   virtual std::string getName(size_t, size_t j, const char*) const;
 
-  static int dbgElm; //!< Element index for debug output
+  static bool extEnr; //!< Flag for integration of external energy
+  static bool dirDer; //!< Flag for calculation of directional derivatives
+  static int  dbgElm; //!< Element index for debug output
 };
 
 #endif

QuasiStaticSIM.C

@@ -13,6 +13,7 @@
 
 #include "QuasiStaticSIM.h"
 #include "HermiteInterpolator.h"
+#include "FractureElasticityVoigt.h"
 #include "SIMoutput.h"
 #include "TimeStep.h"
 #include "Utilities.h"
@@ -30,6 +31,7 @@
 QuasiStaticSIM::QuasiStaticSIM (SIMbase& sim) : NonLinSIM(sim)
 {
   numPt = numPtPos = nDump = 0;
+  version = 1;
 }
 
 
@@ -66,6 +68,7 @@ bool QuasiStaticSIM::parse (const TiXmlElement* elem)
         numPt = params.size();
       }
       utl::getAttribute(child,"nDump",nDump);
+      utl::getAttribute(child,"version",version);
     }
 
   return this->NonLinSIM::parse(elem);
@@ -84,6 +87,10 @@ void QuasiStaticSIM::printProblem () const
       IFEM::cout << (i%10 ? " " : "\n\t            ") << params[i];
     if (nDump > 0)
       IFEM::cout <<"\n\tDumping f(alpha) at "<< nDump <<" sampling points.";
+    if (version == 1)
+      IFEM::cout <<"\n\tEvaluating f'(alpha) as {F}_int * {u}_corr.";
+    else if (version == 2)
+      IFEM::cout <<"\n\tEvaluating f'(alpha) by direct integration.";
   }
   IFEM::cout << std::endl;
 }
@@ -93,8 +100,10 @@ bool QuasiStaticSIM::evalEnergyFunctional (const TimeDomain& time,
                                            const Vectors& pSol,
                                            double* fVal, double* fDer)
 {
-  if (fDer)
+  if (fDer && version == 1)
   {
+    // Integrate f'(alpha) as the dot-product between internal force vector
+    // and the solution correction vector
     if (!model.setMode(SIM::INT_FORCES))
       return false;
 
@@ -107,16 +116,23 @@ bool QuasiStaticSIM::evalEnergyFunctional (const TimeDomain& time,
     *fDer = -residual.dot(linsol);
   }
 
-  if (fVal)
+  if (fVal || version == 2)
   {
+    FractureElasticNorm::dirDer = version == 2;
+
     if (!model.setMode(SIM::RECOVERY))
       return false;
 
     Vectors gNorm;
     if (!model.solutionNorms(time,pSol,gNorm))
       return false;
 
-    *fVal = gNorm.front()(1) + gNorm.front()(6);
+    if (fVal)
+      *fVal = gNorm.front()(1) + gNorm.front()(6);
+    if (fDer && version == 2)
+      *fDer = gNorm.front()(5); // f'(alpha) is integrated by the norm class
+
+    FractureElasticNorm::dirDer = false;
   }
 
   return true;
@@ -129,8 +145,11 @@ bool QuasiStaticSIM::lineSearch (TimeStep& param)
   if (numPtPos < 2)
     return true; // No line search
 
+  FractureElasticNorm::extEnr = false; // Disable external energy calculation,
+  // since it uses a path integral valid at the converged configuration only
+
   // Make a temporary copy of the primary solution
-  Vectors tmpSol(1,solution.front());
+  Vectors tmpSol(version == 2 ? 2 : 1, solution.front());
   Vector& solVec = tmpSol.front();
   double curr, prev;
 
@@ -148,7 +167,13 @@ bool QuasiStaticSIM::lineSearch (TimeStep& param)
   std::cout <<"\tLine search? curr: "<< curr <<" prev: "<< prev << std::endl;
 #endif
   if (curr < prev)
+  {
+    FractureElasticNorm::extEnr = true; // Enable external energy calculation
     return true; // No line search needed in this iteration
+  }
+
+  if (version == 2) // Store the solution correction in tmpSol
+    tmpSol.back() = linsol;
 
   // Evaluate f'(alpha) at alpha=0.0
   if (!this->evalEnergyFunctional(param.time,tmpSol,nullptr,&curr))
@@ -171,6 +196,7 @@ bool QuasiStaticSIM::lineSearch (TimeStep& param)
     if (!this->evalEnergyFunctional(param.time,tmpSol,&values[i],&derivs[i]))
       return false;
   }
+  FractureElasticNorm::extEnr = true; // Enable external energy calculation
 
 #if SP_DEBUG > 1
   std::cout <<"\nParameters:\n";

QuasiStaticSIM.h

@@ -52,6 +52,7 @@ class QuasiStaticSIM : public NonLinSIM
   size_t    numPt;    // Total number of alpha-values
   size_t    numPtPos; // Number of non-negative alpha-values
   size_t    nDump;    // Number of grid-points to dump f(alpha) to file for
+  int       version;  // Algorithm version, for internal testing
 };
 
 #endif

SIMDynElasticity.h

@@ -283,6 +283,15 @@ class SIMDynElasticity : public SIMElasticity<Dim>
     return fel ? fel->getCrackPressure() != 0.0 : false;
   }
 
+  //! \brief Computes (possibly problem-dependent) external energy contribution.
+  virtual double externalEnergy(const Vectors& psol) const
+  {
+    if (FractureElasticNorm::extEnr)
+      return this->Dim::externalEnergy(psol);
+
+    return 0.0; // No external energy integration
+  }
+
 protected:
   //! \brief Returns the actual integrand.
   virtual Elasticity* getIntegrand()