Poroelastic fracture, with broken permeability tensor

.gitignore

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

CMakeLists.txt

@@ -21,30 +21,53 @@ set(EXECUTABLE_OUTPUT_PATH ${CMAKE_BINARY_DIR}/bin)
 set(EXTRA_DOXY_PATHS "${PROJECT_SOURCE_DIR} ${PROJECT_BINARY_DIR}")
 
 if(EXISTS ${PROJECT_SOURCE_DIR}/../IFEM-Elasticity)
-  set(EXTRA_DOXY_PATHS "${EXTRA_DOXY_PATHS} \\
-                        ${PROJECT_SOURCE_DIR}/../IFEM-Elasticity")
-  if(NOT TARGET Elasticity)
-    add_subdirectory(../IFEM-Elasticity Elasticity)
-  endif()
-  include_directories(../IFEM-Elasticity)
+  set(ELASTICITY_DIR ../IFEM-Elasticity)
 elseif(EXISTS ${PROJECT_SOURCE_DIR}/../Elasticity)
-  set(EXTRA_DOXY_PATHS "${EXTRA_DOXY_PATHS} \\
-                        ${PROJECT_SOURCE_DIR}/../Elasticity")
-  if(NOT TARGET Elasticity)
-    add_subdirectory(../Elasticity Elasticity)
+  set(ELASTICITY_DIR ../Elasticity)
+else()
+  message(FATAL_ERROR "Need IFEM-Elasticity in a sibling directory.")
+endif()
+
+set(EXTRA_DOXY_PATHS "${EXTRA_DOXY_PATHS} \\
+                      ${PROJECT_SOURCE_DIR}/${ELASTICITY_DIR}")
+if(NOT TARGET Elasticity)
+  add_subdirectory(${ELASTICITY_DIR} Elasticity)
+endif()
+include_directories(${ELASTICITY_DIR})
+
+set(POROELASTIC_DIR ${PROJECT_SOURCE_DIR}/../IFEM-PoroElasticity/PoroElastic)
+if(NOT EXISTS ${POROELASTIC_DIR})
+  set(POROELASTIC_DIR ${PROJECT_SOURCE_DIR}/../PoroElasticity/PoroElastic)
+endif()
+
+if(EXISTS ${POROELASTIC_DIR})
+  set(EXTRA_DOXY_PATHS "${EXTRA_DOXY_PATHS} ${POROELASTIC_DIR}")
+  if(NOT TARGET PoroElastic)
+    add_subdirectory(${POROELASTIC_DIR} PoroElasticity)
   endif()
-  include_directories(../Elasticity)
+  include_directories(${POROELASTIC_DIR})
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DIFEM_HAS_POROELASTIC")
 endif()
 
-file(GLOB FracEl_SRCS [A-Z]*.C)
+set(FracEl_SRCS FractureArgs.C CahnHilliard.C SIMExplPhaseField.C
+                FractureElasticity.C FractureElasticityVoigt.C)
+if(EXISTS ${POROELASTIC_DIR})
+  set(FracEl_SRCS ${FracEl_SRCS} PoroFracture.C)
+endif()
 
 add_library(CommonFrac STATIC ${FracEl_SRCS})
+set(Common_LIBRARIES CommonFrac Elasticity ${IFEM_LIBRARIES})
+if(EXISTS ${POROELASTIC_DIR})
+  set(Fracture_LIBRARIES CommonFrac PoroElastic Elasticity ${IFEM_LIBRARIES})
+else()
+  set(Fracture_LIBRARIES ${Common_LIBRARIES})
+endif()
 
 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})
 
 list(APPEND CHECK_SOURCES main_CahnHilliard.C main_FractureDynamics.C ${FracEl_SRCS})
 
@@ -66,6 +89,8 @@ else()
                          ${PROJECT_SOURCE_DIR}/Test/*_CH.reg)
   file(GLOB FD_TESTFILES RELATIVE ${PROJECT_SOURCE_DIR}/Test
                          ${PROJECT_SOURCE_DIR}/Test/*_FD.reg)
+  file(GLOB FP_TESTFILES RELATIVE ${PROJECT_SOURCE_DIR}/Test
+                         ${PROJECT_SOURCE_DIR}/Test/*_FDporo.reg)
   foreach(TESTFILE ${CH_TESTFILES})
     ifem_add_test(${TESTFILE} CahnHilliard)
   endforeach()
@@ -81,14 +106,19 @@ else()
   foreach(TESTFILE ${FD_TESTFILES})
     ifem_add_test(${TESTFILE} FractureDynamics)
   endforeach()
+  if(EXISTS ${POROELASTIC_DIR})
+    foreach(TESTFILE ${FP_TESTFILES})
+      ifem_add_test(${TESTFILE} FractureDynamics)
+    endforeach()
+  endif()
 endif()
 list(APPEND TEST_APPS CahnHilliard FractureDynamics)
 
 # Unit tests
 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)

FractureArgs.C

@@ -20,6 +20,8 @@ FractureArgs::FractureArgs () : SIMargsBase("fracturedynamics")
 {
   inpfile = nullptr;
   integrator = coupling = 1;
+  poroEl = expPhase = false;
+  stopT = -1.0;
 }
 
 
@@ -41,6 +43,10 @@ bool FractureArgs::parseArg (const char* argv)
     integrator = 4;
   else if (!strcmp(argv,"-oldHHT"))
     integrator = 5;
+  else if (!strncmp(argv,"-poro",5))
+    poroEl = true;
+  else if (!strncmp(argv,"-explcr",7))
+    expPhase = true;
   else if (!strncmp(argv,"-noadap",7))
     adap = false;
   else
@@ -78,6 +84,10 @@ bool FractureArgs::parse (const TiXmlElement* elem)
         integrator = 4;
       else if (!strcasecmp(child->Value(),"hilberhughestaylor"))
         integrator = 4;
+      else if (!strcasecmp(child->Value(),"poroelastic"))
+        poroEl = true;
+      else if (!strcasecmp(child->Value(),"explicitphase"))
+        expPhase = true;
   }
 
   return this->SIMargsBase::parse(elem);

FractureArgs.h

@@ -30,6 +30,9 @@ class FractureArgs : public SIMargsBase
   //! 4=nonlinear Hilber-Hughes-Taylor)
   char integrator;
   char coupling; //!< Coupling flag (0: none, 1: staggered, 2: semi-implicit)
+  bool poroEl;   //!< If \e true, use the poroelastic solver
+  bool expPhase; //!< If \e true, use an explicit phase field
+  double stopT;  //!< Stop time of the simulation
 
   //! \brief Default constructor.
   FractureArgs();
@@ -40,7 +43,7 @@ class FractureArgs : public SIMargsBase
   void parseFile(const char* argv, int& iarg);
 
 protected:
-  //! \brief Parse an element from the input file
+  //! \brief Parses an element from the input file.
   virtual bool parse(const TiXmlElement* elem);
 };
 

FractureElasticity.C

@@ -121,6 +121,19 @@ void FractureElasticity::printLog () const
 }
 
 
+void FractureElasticity::setMode (SIM::SolutionMode mode)
+{
+  this->Elasticity::setMode(mode);
+  if (eC <= 1) return; // no parent integrand
+
+  eKg = 0; // No geometric stiffness (assuming linear behaviour)
+  eM = eS = 0; // Inertia and external forces are calculated by parent integrand
+  if (eKm) eKm = 2; // Index for stiffness matrix in parent integrand
+  if (iS) iS = 2; // Index for internal force vector in parent integrand
+  eC = mode == SIM::DYNAMIC ? 5 : 3; // include velocity & acceleration vectors
+}
+
+
 void FractureElasticity::initIntegration (size_t nGp, size_t)
 {
   // Initialize internal tensile energy buffer
@@ -508,12 +521,11 @@ bool FractureElasticity::evalInt (LocalIntegral& elmInt,
 }
 
 
-bool FractureElasticity::evalSol (Vector& s,
-                                  const FiniteElement& fe, const Vec3& X,
-                                  const std::vector<int>& MNPC) const
+bool FractureElasticity::getElementSolution (Vectors& eV,
+                                             const std::vector<int>& MNPC) const
 {
   // Extract element displacements
-  Vectors eV(1+eC);
+  eV.resize(1+eC);
   int ierr = 0;
   if (!mySol.empty() && !mySol.front().empty())
     ierr = utl::gather(MNPC,nsd,mySol.front(),eV.front());
@@ -522,14 +534,21 @@ bool FractureElasticity::evalSol (Vector& s,
   if (!myCVec.empty() && ierr == 0)
     ierr = utl::gather(MNPC,1,myCVec,eV[eC]);
 
-  if (ierr > 0)
-  {
-    std::cerr <<" *** FractureElasticity::evalSol: Detected "<< ierr
-              <<" node numbers out of range."<< std::endl;
-    return false;
-  }
+  if (ierr == 0)
+    return true;
 
-  return this->evalSol2(s,eV,fe,X);
+  std::cerr <<" *** FractureElasticity::getElementSolution: Detected "<< ierr
+            <<" node numbers out of range."<< std::endl;
+  return false;
+}
+
+
+bool FractureElasticity::evalSol (Vector& s,
+                                  const FiniteElement& fe, const Vec3& X,
+                                  const std::vector<int>& MNPC) const
+{
+  Vectors eV(1+eC);
+  return this->getElementSolution(eV,MNPC) && this->evalSol2(s,eV,fe,X);
 }
 
 
@@ -621,6 +640,49 @@ bool FractureElasticity::evalSol (Vector& s, const Vectors& eV,
 }
 
 
+double FractureElasticity::evalPhaseField (Vec3& gradD,
+                                           const Vectors& eV,
+                                           const Vector& N,
+                                           const Matrix& dNdX) const
+{
+  if (eV.size() <= eC)
+  {
+    std::cerr <<" *** FractureElasticity::evalPhaseField:"
+              <<" Missing phase field solution vector."<< std::endl;
+    return -1.1;
+
+  }
+  else if (eV[eC].empty()) // No phase field ==> no crack yet
+  {
+    gradD = 0.0;
+    return 0.0;
+  }
+  else if (eV[eC].size() != N.size())
+  {
+    std::cerr <<" *** FractureElasticity::evalPhaseField:"
+              <<" Invalid phase field vector.\n     size(eC) = "
+              << eV[eC].size() <<"   size(N) = "<< N.size() << std::endl;
+    return -1.2;
+  }
+
+  // Invert the nodal phase field values, D = 1 - C,
+  // since that is the convention used in the Miehe papers
+  Vector eD(eV[eC].size()), tmp(nsd);
+  for (size_t i = 0; i < eD.size(); i++)
+    eD[i] = 1.0 - eV[eC][i];
+
+  // Compute the phase field gradient, dD/dX
+  if (dNdX.multiply(eD,tmp,true))
+    gradD = tmp;
+  else
+    return -2.0;
+
+  // Compute the phase field value, filtering out values outside [0.0,1.0]
+  double d = eD.dot(N);
+  return d > 1.0 ? 1.0 : (d < 0.0 ? 0.0 : d);
+}
+
+
 size_t FractureElasticity::getNoFields (int fld) const
 {
   return this->Elasticity::getNoFields(fld) + (fld == 2 ? 4 : 0);

FractureElasticity.h

@@ -43,8 +43,9 @@ class FractureElasticity : public Elasticity
   //! \brief Prints out the problem definition to the log stream.
   virtual void printLog() const;
 
-  //! \brief Sets the number of solution variables per node.
-  void setVar(unsigned short int n) { npv = n; }
+  //! \brief Defines the solution mode before the element assembly is started.
+  //! \param[in] mode The solution mode to use
+  virtual void setMode(SIM::SolutionMode mode);
 
   //! \brief Returns whether this norm has explicit boundary contributions.
   virtual bool hasBoundaryTerms() const { return m_mode != SIM::RECOVERY; }
@@ -96,6 +97,21 @@ class FractureElasticity : public Elasticity
   virtual bool evalSol(Vector& s, const Vectors& eV, const FiniteElement& fe,
                        const Vec3& X, bool toLocal, Vec3* pdir) const;
 
+  //! \brief Retrieves the element solution vectors.
+  //! \param[out] eV Element solution vectors
+  //! \param[in] MNPC Nodal point correspondance for the basis function values
+  virtual bool getElementSolution(Vectors& eV,
+                                  const std::vector<int>& MNPC) const;
+
+  //! \brief Evaluates the phase field and gradient at current point.
+  //! \param[out] gradD Phase field gradient at current point
+  //! \param[in] eV Element solution vectors
+  //! \param[in] N Basis function values at current point
+  //! \param[in] dNdX Basis function gradients at current point
+  //! \return Phase field value at current point
+  double evalPhaseField(Vec3& gradD, const Vectors& eV,
+                        const Vector& N, const Matrix& dNdX) const;
+
   //! \brief Returns a pointer to the Gauss-point tensile energy array.
   virtual const RealArray* getTensileEnergy() const { return &myPhi; }