PoroElasticity integration

.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,13 @@ endif()
 file(GLOB FracEl_SRCS [A-Z]*.C)
 
 add_library(CommonFrac STATIC ${FracEl_SRCS})
+set(Common_LIBRARIES CommonFrac 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 PoroElastic ${Common_LIBRARIES})
 
 # Installation
 install(TARGETS CahnHilliard FractureDynamics DESTINATION bin)
@@ -79,7 +91,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})
+                  PoroElastic ${Common_LIBRARIES})
 
 if(IFEM_COMMON_APP_BUILD)
   set(TEST_APPS ${TEST_APPS} PARENT_SCOPE)

FractureElasticity.C

@@ -27,10 +27,23 @@
 #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
 }
 
 
@@ -44,28 +57,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;
 
@@ -211,10 +226,10 @@ bool FractureElasticity::evalStress (double lambda, double mu, double Gc,
 
 
 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 +269,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;
@@ -364,8 +379,8 @@ bool FractureElasticity::evalSol (Vector& s,
   // Extract element displacements
   Vectors eV(2);
   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)
@@ -389,7 +404,7 @@ bool FractureElasticity::evalSol (Vector& s, const Vectors& eV,
 {
   PROFILE3("FractureEl::evalSol");
 
-  if (eV.size() < 2)
+  if (eV.size() <= eC)
   {
     std::cerr <<" *** FractureElasticity::evalSol: Missing solution vector."
               << std::endl;
@@ -402,10 +417,10 @@ bool FractureElasticity::evalSol (Vector& s, const Vectors& eV,
               << fe.dNdX.rows() <<","<< fe.dNdX.cols() << std::endl;
     return false;
   }
-  else if (!eV[1].empty() && eV[1].size() != fe.N.size())
+  else if (!eV[eC].empty() && eV[eC].size() != fe.N.size())
   {
     std::cerr <<" *** FractureElasticity::evalSol: Invalid phase field vector."
-              <<"\n     size(eC) = "<< eV[1].size() <<"   size(N) = "
+              <<"\n     size(eC) = "<< eV[eC].size() <<"   size(N) = "
               << fe.N.size() << std::endl;
     return false;
   }
@@ -428,7 +443,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;
 

FractureElasticity.h

@@ -29,9 +29,16 @@ class FractureElasticity : public Elasticity
   //! \brief The constructor invokes the parent class constructor only.
   //! \param[in] n Number of spatial dimensions
   FractureElasticity(unsigned short int n);
+  //! \brief Constructor for integrands with a parent integrand.
+  //! \param parent The parent integrand of this one
+  //! \param[in] n Number of spatial dimensions
+  FractureElasticity(IntegrandBase* parent, unsigned short int n);
   //! \brief Empty destructor.
   virtual ~FractureElasticity() {}
 
+  //! \brief Sets the number of solution variables per node.
+  void setVar(unsigned short int n) { npv = n; }
+
   //! \brief Initializes the integrand with the number of integration points.
   //! \param[in] nGp Total number of interior integration points
   virtual void initIntegration(size_t nGp, size_t);
@@ -103,13 +110,17 @@ class FractureElasticity : public Elasticity
                   bool postProc = false) const;
 
   //! \brief Evaluates the stress degradation function \a g(c) at current point.
-  double getStressDegradation(const Vector& N, const Vector& eC) const;
+  double getStressDegradation(const Vector& N, const Vectors& eV) const;
+
+private:
+  unsigned short int eC; //!< Zero-based index to element phase field vector
 
 protected:
   double alpha;  //!< Relaxation factor for the crack phase field
   Vector myCVec; //!< Crack phase field values at nodal points
 
   mutable RealArray myPhi; //!< Tensile energy density at integration points
+  Vectors&          mySol; //!< Primary solution vectors for current patch
 };
 
 #endif

FractureElasticityVoigt.C

@@ -268,7 +268,7 @@ bool FractureElasticityVoigt::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;
@@ -391,7 +391,7 @@ bool FractureElasticNorm::evalInt (LocalIntegral& elmInt,
 
   // Evaluate the strain energy at this point
   double Phi[4];
-  double Gc = p.getStressDegradation(fe.N,elmInt.vec[1]);
+  double Gc = p.getStressDegradation(fe.N,elmInt.vec);
   if (!p.evalStress(lambda,mu,Gc,eps,Phi,nullptr,nullptr,true,printElm))
     return false;
 

FractureElasticityVoigt.h

@@ -12,7 +12,7 @@
 //==============================================================================
 
 #ifndef _FRACTURE_ELASTICITY_VOIGT_H
-#define _FRACTURE_ELASTICITY_VOIGH_H
+#define _FRACTURE_ELASTICITY_VOIGT_H
 
 #include "FractureElasticity.h"
 
@@ -30,6 +30,11 @@ class FractureElasticityVoigt : public FractureElasticity
   //! \brief The constructor invokes the parent class constructor only.
   //! \param[in] n Number of spatial dimensions
   FractureElasticityVoigt(unsigned short int n) : FractureElasticity(n) {}
+  //! \brief Constructor for integrands with a parent integrand.
+  //! \param parent The parent integrand of this one
+  //! \param[in] n Number of spatial dimensions
+  FractureElasticityVoigt(IntegrandBase* parent, unsigned short int n)
+    : FractureElasticity(parent,n) {}
   //! \brief Empty destructor.
   virtual ~FractureElasticityVoigt() {}