Poro fracture WIP

.gitignore

@@ -0,0 +1,3 @@
+*.log
+*.vtf
+*~

CMakeLists.txt

@@ -26,6 +26,17 @@ if(NOT WIN32)
   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall")
 endif(NOT WIN32)
 
+if(EXISTS ${PROJECT_SOURCE_DIR}/../IFEM-PoroElasticity)
+  if(NOT TARGET PoroElastic)
+    add_subdirectory(../IFEM-PoroElasticity/PoroElastic PoroElasticity)
+  endif()
+  include_directories(../IFEM-PoroElasticity/PoroElastic)
+elseif(EXISTS ${PROJECT_SOURCE_DIR}/../PoroElasticity)
+  if(NOT TARGET PoroElastic)
+    add_subdirectory(../PoroElasticity/PoroElastic PoroElasticity)
+  endif()
+  include_directories(../PoroElasticity/PoroElastic)
+endif()
 if(EXISTS ${PROJECT_SOURCE_DIR}/../IFEM-Elasticity)
   if(NOT TARGET Elasticity)
     add_subdirectory(../IFEM-Elasticity Elasticity)
@@ -41,12 +52,14 @@ endif()
 file(GLOB FracEl_SRCS [A-Z]*.C)
 
 add_library(CommonFrac STATIC ${FracEl_SRCS})
+set(Common_LIBRARIES CommonFrac Elasticity ${IFEM_LIBRARIES})
+set(Fracture_LIBRARIES CommonFrac PoroElastic Elasticity ${IFEM_LIBRARIES})
 
 add_executable(CahnHilliard main_CahnHilliard.C)
 add_executable(FractureDynamics main_FractureDynamics.C)
 
-target_link_libraries(CahnHilliard CommonFrac Elasticity ${IFEM_LIBRARIES})
-target_link_libraries(FractureDynamics CommonFrac Elasticity ${IFEM_LIBRARIES})
+target_link_libraries(CahnHilliard ${Common_LIBRARIES})
+target_link_libraries(FractureDynamics ${Fracture_LIBRARIES})
 
 # Installation
 install(TARGETS CahnHilliard FractureDynamics DESTINATION bin)
@@ -79,7 +92,7 @@ list(APPEND TEST_APPS CahnHilliard FractureDynamics)
 IFEM_add_test_app(${PROJECT_SOURCE_DIR}/Test/*.C
                   ${PROJECT_SOURCE_DIR}/Test
                   FractureDynamics
-                  CommonFrac Elasticity ${IFEM_LIBRARIES})
+                  ${Fracture_LIBRARIES})
 
 if(IFEM_COMMON_APP_BUILD)
   set(TEST_APPS ${TEST_APPS} PARENT_SCOPE)

FractureElasticity.C

@@ -27,10 +27,34 @@
 #endif
 
 
-FractureElasticity::FractureElasticity (unsigned short int n) : Elasticity(n)
+FractureElasticity::FractureElasticity (unsigned short int n)
+  : Elasticity(n), mySol(primsol)
 {
   alpha = 0.0;
   this->registerVector("phasefield",&myCVec);
+  eC = 1; // Assuming second vector is phase field 
+}
+
+
+FractureElasticity::FractureElasticity (IntegrandBase* parent,
+                                        unsigned short int n)
+  : Elasticity(n), mySol(parent->getSolutions())
+{
+  alpha = 0.0;
+  parent->registerVector("phasefield",&myCVec);
+  // Assuming second vector is pressure, third vector is pressure velocity
+  eC = 3; // and fourth vector is the phase field
+}
+
+
+void FractureElasticity::setMode (SIM::SolutionMode mode)
+{
+  this->Elasticity::setMode(mode);
+  if (eC != 3) return; // no parent integrand
+
+  eKg = 0; // No geometric stiffness (assuming linear behaviour)
+  eM = eS = 0; // Inertia and external forces are calulated by parent integrand
+  if (eKm) eKm = 2; // Index for stiffness matrix in parent integrand
 }
 
 
@@ -44,28 +68,30 @@ void FractureElasticity::initIntegration (size_t nGp, size_t)
 bool FractureElasticity::initElement (const std::vector<int>& MNPC,
                                       LocalIntegral& elmInt)
 {
-  if (primsol.empty())
+  if (mySol.empty())
   {
     std::cerr <<" *** FractureElasticity::initElement:"
               <<" No primary solution vectors."<< std::endl;
     return false;
   }
 
   int ierr = 0;
-  elmInt.vec.resize(primsol.size()+1);
-
-  // Extract displacement vector for this element
-  if (!primsol.front().empty())
-    ierr = utl::gather(MNPC,npv,primsol.front(),elmInt.vec.front());
+  if (elmInt.vec.empty())
+  {
+    // Extract displacement vector for this element
+    elmInt.vec.resize(mySol.size()+eC);
+    if (!mySol.front().empty())
+      ierr = utl::gather(MNPC,npv,mySol.front(),elmInt.vec.front());
+  }
 
   // Extract crack phase field vector for this element
   if (ierr == 0 && !myCVec.empty())
-    ierr = utl::gather(MNPC,1,myCVec,elmInt.vec[1]);
+    ierr = utl::gather(MNPC,1,myCVec,elmInt.vec[eC]);
 
   // Extract velocity and acceleration vectors for this element
-  for (size_t i = 1; i < primsol.size() && ierr == 0; i++)
-    if (!primsol[i].empty())
-      ierr = utl::gather(MNPC,npv,primsol[i],elmInt.vec[1+i]);
+  for (size_t i = 1; i < mySol.size() && ierr == 0; i++)
+    if (!mySol[i].empty())
+      ierr = utl::gather(MNPC,npv,mySol[i],elmInt.vec[eC+i]);
 
   if (ierr == 0) return true;
 
@@ -210,11 +236,19 @@ bool FractureElasticity::evalStress (double lambda, double mu, double Gc,
 }
 
 
+bool FractureElasticity::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);
+}
+
+
 double FractureElasticity::getStressDegradation (const Vector& N,
-                                                 const Vector& eC) const
+                                                 const Vectors& eV) const
 {
   // Evaluate the crack phase field function, c(X)
-  double c = eC.empty() ? 1.0 : N.dot(eC);
+  double c = eV[eC].empty() ? 1.0 : N.dot(eV[eC]);
   // Evaluate the stress degradation function, g(c), ignoring negative values
   return c > 0.0 ? (1.0-alpha)*c*c + alpha : alpha;
 }
@@ -254,7 +288,7 @@ bool FractureElasticity::evalInt (LocalIntegral& elmInt,
       return false;
 
     // Evaluate the stress degradation function
-    double Gc = this->getStressDegradation(fe.N,elmInt.vec[1]);
+    double Gc = this->getStressDegradation(fe.N,elmInt.vec);
 #if INT_DEBUG > 3
     std::cout <<"lambda = "<< lambda <<" mu = "<< mu <<" G(c) = "<< Gc <<"\n";
     if (lHaveStrains) std::cout <<"eps =\n"<< eps;
@@ -328,32 +362,57 @@ bool FractureElasticity::evalInt (LocalIntegral& elmInt,
 }
 
 
-bool FractureElasticity::evalBou (LocalIntegral& elmInt,
-                                  const FiniteElement& fe, const Vec3& X,
-                                  const Vec3& normal) const
+bool FractureElasticity::checkSolVec (const Vectors& eV, size_t nen,
+                                      const char* name) const
 {
-  if (!tracFld && !fluxFld)
-  {
-    std::cerr <<" *** FractureElasticity::evalBou: No tractions."<< std::endl;
-    return false;
-  }
-  else if (!eS)
-  {
-    if (m_mode == SIM::RECOVERY) return true;
-    std::cerr <<" *** FractureElasticity::evalBou: No load vector."<< std::endl;
-    return false;
-  }
+  if (eV.size() <= eC)
+    std::cerr <<" *** FractureElasticity::"<< name
+              <<": Missing solution vector."<< std::endl;
+  else if (!eV.front().empty() && eV.front().size() != nsd*nen)
+    std::cerr <<" *** FractureElasticity::"<< name
+              <<": Invalid displacement vector.\n     size(eV) = "
+              << eV.front().size() <<"   size(dNdX) = "<< nsd*nen << std::endl;
+  else if (!eV[eC].empty() && eV[eC].size() != nen)
+    std::cerr <<" *** FractureElasticity::"<< name
+              <<": Invalid phase field vector.\n     size(eC) = "
+              << eV[eC].size() <<"   size(N) = "<< nen << std::endl;
+  else
+    return true;
 
-  // Evaluate the surface traction
-  Vec3 T = this->getTraction(X,normal);
+  return false;
+}
 
-  // Integrate the force vector
-  Vector& ES = static_cast<ElmMats&>(elmInt).b[eS-1];
-  for (size_t a = 1; a <= fe.N.size(); a++)
-    for (unsigned short int i = 1; i <= nsd; i++)
-      ES(nsd*(a-1)+i) += T[i-1]*fe.N(a)*fe.detJxW;
 
-  return true;
+double FractureElasticity::crackStretch (const Vectors& eV,
+                                         const FiniteElement& fe,
+                                         const Vec3& X) const
+{
+  if (!this->checkSolVec(eV,fe.N.size(),"crackStretch"))
+    return -1.0;
+  else if (eV[eC].empty())
+    return 1.0; // No phase field ==> unit stretch
+
+  // Compute the deformation gradient (F)
+  Tensor F(nDF), C(nDF);
+  if (!this->formDefGradient(eV.front(),fe.N,fe.dNdX,X.x,F))
+    return -2.0;
+
+  // Compute the inverse right Cauchy-Green tensor (C^-1)
+  SymmTensor Ci(nDF);
+  if (Ci.rightCauchyGreen(F).inverse() == 0.0)
+    return -3.0;
+
+  // Compute the phase field gradient (dD/dX where D = 1 - C)
+  Vector eD(eV[eC]), tmp(nsd);
+  for (size_t i = 0; i < eD.size(); i++)
+    eD[i] = 1.0 - eD[i]; // d = 1 - c
+  if (!fe.dNdX.multiply(eD,tmp,true))
+    return -4.0;
+
+  // Calculate the perpendicular crack stretch
+  // lambda = gradD*gradD / gradD*Ci*gradD (see eq. (106) in Miehe2015pfm3)
+  Vec3 gradD(tmp);
+  return gradD.length2() / (gradD*(Ci*gradD));
 }
 
 
@@ -362,14 +421,14 @@ bool FractureElasticity::evalSol (Vector& s,
                                   const std::vector<int>& MNPC) const
 {
   // Extract element displacements
-  Vectors eV(2);
+  Vectors eV(1+eC);
   int ierr = 0;
-  if (!primsol.empty() && !primsol.front().empty())
-    ierr = utl::gather(MNPC,nsd,primsol.front(),eV.front());
+  if (!mySol.empty() && !mySol.front().empty())
+    ierr = utl::gather(MNPC,nsd,mySol.front(),eV.front());
 
   // Extract crack phase field vector for this element
   if (!myCVec.empty() && ierr == 0)
-    ierr = utl::gather(MNPC,1,myCVec,eV.back());
+    ierr = utl::gather(MNPC,1,myCVec,eV[eC]);
 
   if (ierr > 0)
   {
@@ -378,7 +437,6 @@ bool FractureElasticity::evalSol (Vector& s,
     return false;
   }
 
-
   return this->evalSol2(s,eV,fe,X);
 }
 
@@ -389,26 +447,8 @@ bool FractureElasticity::evalSol (Vector& s, const Vectors& eV,
 {
   PROFILE3("FractureEl::evalSol");
 
-  if (eV.size() < 2)
-  {
-    std::cerr <<" *** FractureElasticity::evalSol: Missing solution vector."
-              << std::endl;
-    return false;
-  }
-  else if (!eV.front().empty() && eV.front().size() != fe.dNdX.rows()*nsd)
-  {
-    std::cerr <<" *** FractureElasticity::evalSol: Invalid displacement vector."
-              <<"\n     size(eV) = "<< eV.front().size() <<"   size(dNdX) = "
-              << fe.dNdX.rows() <<","<< fe.dNdX.cols() << std::endl;
-    return false;
-  }
-  else if (!eV[1].empty() && eV[1].size() != fe.N.size())
-  {
-    std::cerr <<" *** FractureElasticity::evalSol: Invalid phase field vector."
-              <<"\n     size(eC) = "<< eV[1].size() <<"   size(N) = "
-              << fe.N.size() << std::endl;
+  if (!this->checkSolVec(eV,fe.N.size(),"evalSol"))
     return false;
-  }
 
   // Evaluate the symmetric strain tensor, eps
   Matrix Bmat;
@@ -428,7 +468,7 @@ bool FractureElasticity::evalSol (Vector& s, const Vectors& eV,
   // Evaluate the stress state at this point
   SymmTensor sigma(nsd);
   double Phi[3];
-  double Gc = this->getStressDegradation(fe.N,eV[1]);
+  double Gc = this->getStressDegradation(fe.N,eV);
   if (!this->evalStress(lambda,mu,Gc,eps,Phi,sigma))
     return false;
 
@@ -471,14 +511,6 @@ bool FractureElasticity::evalSol (Vector& s, const Vectors& eV,
 }
 
 
-bool FractureElasticity::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);
-}
-
-
 size_t FractureElasticity::getNoFields (int fld) const
 {
   return this->Elasticity::getNoFields(fld) + (fld == 2 ? 4 : 0);