From 0e3c9fb2f38c3c30540f5da89d23fa62c3ff422c Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Tue, 12 Apr 2016 16:29:50 +0200
Subject: [PATCH 1/4] Added: .gitignore Changed: Moved template functions to
 separate file such that they can be reused.

---
 .gitignore              |   3 +
 SIMDriver.h             |  59 +++++++++++++++++++
 main_FractureDynamics.C | 123 ++--------------------------------------
 runTimeInt.h            |  92 ++++++++++++++++++++++++++++++
 4 files changed, 160 insertions(+), 117 deletions(-)
 create mode 100644 .gitignore
 create mode 100644 SIMDriver.h
 create mode 100644 runTimeInt.h

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..9eeb9cf
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+*.log
+*.vtf
+*~
diff --git a/SIMDriver.h b/SIMDriver.h
new file mode 100644
index 0000000..58bbff3
--- /dev/null
+++ b/SIMDriver.h
@@ -0,0 +1,59 @@
+// $Id$
+//==============================================================================
+//!
+//! \file SIMDriver.h
+//!
+//! \date Nov 15 2015
+//!
+//! \author Knut Morten Okstad
+//!
+//! \brief Simulation driver class template.
+//!
+//==============================================================================
+
+#ifndef _SIM_DRIVER_H
+#define _SIM_DRIVER_H
+
+
+/*!
+  \brief Simulation driver template class.
+
+  \details Only the parse method is reimplemented here to handle that the
+  time stepping parameters may be located within the specified context, and
+  for obtaining the file name for global energy output from child simulator.
+*/
+
+template<class T, template<class S> class Solver>
+class SIMDriver : public Solver<T>
+{
+public:
+  //! \brief The constructor initializes the reference to the actual solver.
+  SIMDriver(T& s, const char* c = nullptr) : Solver<T>(s), context(c) {}
+  //! \brief Empty destructor.
+  virtual ~SIMDriver() {}
+
+protected:
+  //! \brief Parses a data section from an XML element.
+  virtual bool parse(const TiXmlElement* elem)
+  {
+    if (!strcasecmp(elem->Value(),context))
+    {
+      const TiXmlElement* child = elem->FirstChildElement();
+      for (; child; child = child->NextSiblingElement())
+        this->SIMSolver<T>::parse(child);
+    }
+    else if (!strcasecmp(elem->Value(),"postprocessing"))
+    {
+      const TiXmlElement* child = elem->FirstChildElement("energyfile");
+      if (child && child->FirstChild())
+        this->S1.setEnergyFile(child->FirstChild()->Value());
+    }
+
+    return this->Solver<T>::parse(elem);
+  }
+
+private:
+  const char* context; //!< XML-tag to search for time-stepping input within
+};
+
+#endif
diff --git a/main_FractureDynamics.C b/main_FractureDynamics.C
index 24d68e6..2271a9e 100644
--- a/main_FractureDynamics.C
+++ b/main_FractureDynamics.C
@@ -17,56 +17,12 @@
 #include "SIMDynElasticity.h"
 #include "SIMPhaseField.h"
 #include "SIMFractureDynamics.h"
-#include "SIMCoupledSI.h"
 #include "SIMSolver.h"
-#include "SIMSolverTS.h"
-#include "GenAlphaSIM.h"
-#include "NonLinSIM.h"
+#include "SIMDriver.h"
 #include "ASMstruct.h"
 #include "AppCommon.h"
 
 
-/*!
-  \brief Dynamic simulation driver.
-
-  \details Only the parse method is reimplemented here to handle that the
-  time stepping parameters may be located within the specified context.
-*/
-
-template<class T, template<class S1> class Solver>
-class SIMDriver : public Solver<T>
-{
-public:
-  //! \brief The constructor initializes the reference to the actual solver.
-  SIMDriver(T& s, const char* c = nullptr) : Solver<T>(s), context(c) {}
-  //! \brief Empty destructor.
-  virtual ~SIMDriver() {}
-
-protected:
-  //! \brief Parses a data section from an XML element.
-  virtual bool parse(const TiXmlElement* elem)
-  {
-    if (!strcasecmp(elem->Value(),context))
-    {
-      const TiXmlElement* child = elem->FirstChildElement();
-      for (; child; child = child->NextSiblingElement())
-        this->SIMSolver<T>::parse(child);
-    }
-    else if (!strcasecmp(elem->Value(),"postprocessing"))
-    {
-      const TiXmlElement* child = elem->FirstChildElement("energyfile");
-      if (child && child->FirstChild())
-        this->S1.setEnergyFile(child->FirstChild()->Value());
-    }
-
-    return this->Solver<T>::parse(elem);
-  }
-
-private:
-  const char* context; //!< XML-tag to search for time-stepping input within
-};
-
-
 /*!
   \brief Creates the combined fracture simulator and launches the simulation.
   \param[in] infile The input file to parse
@@ -75,7 +31,7 @@ private:
 template<class Dim, class Integrator,
          template<class T1, class T2> class Cpl,
          template<class T1> class Solver=SIMSolver>
-int runSimulator2 (char* infile)
+int runCplSimulator (char* infile)
 {
   typedef SIMDynElasticity<Dim,Integrator> SIMElastoDynamics;
   typedef SIMPhaseField<Dim>               SIMCrackField;
@@ -141,7 +97,7 @@ int runSimulator2 (char* infile)
 */
 
 template<class Dim, class Integrator=NewmarkSIM>
-int runSimulator3 (char* infile, const char* context = "newmarksolver")
+int runSimulator (char* infile, const char* context = "newmarksolver")
 {
   typedef SIMDynElasticity<Dim,Integrator> SIMElastoDynamics;
 
@@ -186,74 +142,7 @@ int runSimulator3 (char* infile, const char* context = "newmarksolver")
 }
 
 
-/*!
-  \brief Creates the combined fracture simulator and launches the simulation.
-  \param[in] infile The input file to parse
-  \param[in] timeslabs Use time-slab adaptive solver
-*/
-
-template<class Dim, class Integrator, template<class T1, class T2> class Cpl>
-int runSolver (char* infile, bool timeslabs)
-{
-  if (timeslabs)
-    return runSimulator2<Dim,Integrator,Cpl,SIMSolverTS>(infile);
-
-  return runSimulator2<Dim,Integrator,Cpl>(infile);
-}
-
-
-/*!
-  \brief Creates the combined fracture simulator and launches the simulation.
-  \param[in] infile The input file to parse
-  \param[in] coupling Coupling flag (0: none, 1: staggered, 2: semi-implicit)
-  \param[in] timeslabs Use time-slab adaptive solver
-*/
-
-template<class Dim, class Integrator=NewmarkSIM>
-int runSimulator1 (char* infile, char coupling, bool timeslabs)
-{
-  if (coupling == 1)
-    return runSolver<Dim,Integrator,SIMCoupled>(infile,timeslabs);
-  else if (coupling == 2)
-    return runSolver<Dim,Integrator,SIMCoupledSI>(infile,timeslabs);
-  else
-    return runSimulator3<Dim,Integrator>(infile); // No phase field coupling
-}
-
-
-/*!
-  \brief Linear quasi-static solution driver.
-*/
-
-class LinSIM : public NonLinSIM
-{
-public:
-  //! \brief The constructor forwards to the parent class constructor.
-  LinSIM(SIMbase& sim) : NonLinSIM(sim,NonLinSIM::NONE) { fromIni = false; }
-  //! \brief Empty destructor.
-  virtual ~LinSIM() {}
-};
-
-
-/*!
-  \brief Creates the combined fracture simulator and launches the simulation.
-  \param[in] infile The input file to parse
-  \param[in] integrator The time integrator to use (0=linear quasi-static,
-             no phase-field coupling, 1=linear Newmark, 2=Generalized alpha)
-  \param[in] coupling Coupling flag (0: none, 1: staggered, 2: semi-implicit)
-  \param[in] timeslabs Use time-slab adaptive solver
-*/
-
-template<class Dim>
-int runSimulator (char* infile, char integrator, char coupling, bool timeslabs)
-{
-  if (integrator == 2)
-    return runSimulator1<Dim,GenAlphaSIM>(infile,coupling,timeslabs);
-  else if (integrator > 0)
-    return runSimulator1<Dim>(infile,coupling,timeslabs);
-  else
-    return runSimulator3<Dim,LinSIM>(infile,"staticsolver");
-}
+#include "runTimeInt.h"
 
 
 /*!
@@ -320,7 +209,7 @@ int main (int argc, char** argv)
   IFEM::cout << std::endl;
 
   if (twoD)
-    return runSimulator<SIM2D>(infile,integrator,coupling,adaptive);
+    return runTimeInt<SIM2D>(infile,integrator,coupling,adaptive);
   else
-    return runSimulator<SIM3D>(infile,integrator,coupling,adaptive);
+    return runTimeInt<SIM3D>(infile,integrator,coupling,adaptive);
 }
diff --git a/runTimeInt.h b/runTimeInt.h
new file mode 100644
index 0000000..4143680
--- /dev/null
+++ b/runTimeInt.h
@@ -0,0 +1,92 @@
+// $Id$
+//==============================================================================
+//!
+//! \file runTimeInt.h
+//!
+//! \date Dec 1 2015
+//!
+//! \author Knut Morten Okstad
+//!
+//! \brief Coupled simulation drivers with time evolution.
+//!
+//==============================================================================
+
+#ifndef _RUN_TIME_INT_H
+#define _RUN_TIME_INT_H
+
+#include "SIMCoupledSI.h"
+#include "SIMSolverTS.h"
+#include "GenAlphaSIM.h"
+#include "NonLinSIM.h"
+
+
+/*!
+  \brief Linear quasi-static solution driver.
+*/
+
+class LinSIM : public NonLinSIM
+{
+public:
+  //! \brief The constructor forwards to the parent class constructor.
+  LinSIM(SIMbase& sim) : NonLinSIM(sim,NonLinSIM::NONE) { fromIni = false; }
+  //! \brief Empty destructor.
+  virtual ~LinSIM() {}
+};
+
+
+/*!
+  \brief Creates an adaptive simulator and launches the simulation.
+  \param[in] infile The input file to parse
+  \param[in] adaptive If \e true, use time-slab adaptive solver
+*/
+
+template<class Dim, class Integrator, template<class T1, class T2> class Cpl>
+int runAdapSolver (char* infile, bool adaptive)
+{
+  if (adaptive)
+    return runCplSimulator<Dim,Integrator,Cpl,SIMSolverTS>(infile);
+
+  return runCplSimulator<Dim,Integrator,Cpl>(infile);
+}
+
+
+/*!
+  \brief Creates a coupled fracture simulator and launches the simulation.
+  \param[in] infile The input file to parse
+  \param[in] coupling Coupling flag (0: none, 1: staggered, 2: semi-implicit)
+  \param[in] adaptive If \e true, use time-slab adaptive solver
+*/
+
+template<class Dim, class Integrator=NewmarkSIM>
+int runCplSolver (char* infile, char coupling, bool adaptive)
+{
+  if (coupling == 1)
+    return runAdapSolver<Dim,Integrator,SIMCoupled>(infile,adaptive);
+  else if (coupling == 2)
+    return runAdapSolver<Dim,Integrator,SIMCoupledSI>(infile,adaptive);
+  else
+    return runSimulator<Dim,Integrator>(infile); // Single field, no coupling
+}
+
+
+/*!
+  \brief Creates a simulator with time evolution and launches the simulation.
+  \param[in] infile The input file to parse
+  \param[in] integrator The time integrator to use (0=linear quasi-static,
+             no phase-field coupling, 1=linear Newmark, 2=Generalized alpha)
+  \param[in] coupling Coupling flag (0: none, 1: staggered, 2: semi-implicit)
+  \param[in] adaptive If \e true, use time-slab adaptive solver
+*/
+
+template<class Dim>
+int runTimeInt (char* infile, char integrator, char coupling, bool adaptive)
+{
+  if (integrator == 2)
+    return runCplSolver<Dim,GenAlphaSIM>(infile,coupling,adaptive);
+  else if (integrator > 0)
+    return runCplSolver<Dim>(infile,coupling,adaptive);
+  else
+    return runSimulator<Dim,LinSIM>(infile,"staticsolver"); // Single field
+}
+
+#endif

From 33a52dfac769e6cb28be0ca8105803b2dd9b0ff9 Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Fri, 15 Apr 2016 10:47:32 +0200
Subject: [PATCH 2/4] Changed: Let SIMDynElasticity inherit SIMPoroElasticity

---
 CMakeLists.txt     | 18 ++++++++--
 SIMDynElasticity.h | 88 ++++++++++++++++++++++++----------------------
 2 files changed, 60 insertions(+), 46 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1f24de8..c520594 100644
--- a/CMakeLists.txt
+++ b/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)
diff --git a/SIMDynElasticity.h b/SIMDynElasticity.h
index 5fac336..ec5f9f7 100644
--- a/SIMDynElasticity.h
+++ b/SIMDynElasticity.h
@@ -16,9 +16,8 @@
 #define _SIM_DYN_ELASTICITY_H_
 
 #include "NewmarkSIM.h"
-#include "SIMElasticity.h"
+#include "SIMPoroElasticity.h"
 #include "FractureElasticityVoigt.h"
-#include "DataExporter.h"
 
 
 /*!
@@ -26,18 +25,28 @@
 */
 
 template<class Dim, class DynSIM=NewmarkSIM>
-class SIMDynElasticity : public SIMElasticity<Dim>
+class SIMDynElasticity : public SIMPoroElasticity<Dim>
 {
 public:
   //! \brief Default constructor.
-  SIMDynElasticity() : SIMElasticity<Dim>(false), dSim(*this), vtfStep(0)
-  {
-    Dim::myHeading = "Elasticity solver";
-  }
+  SIMDynElasticity() : dSim(*this), vtfStep(0) {}
+
+  //! \brief Constructor for mixed poroelastic problems.
+  SIMDynElasticity(const std::vector<unsigned char>& flds)
+    : SIMPoroElasticity<Dim>(flds), dSim(*this), vtfStep(0) {}
 
   //! \brief Empty destructor.
   virtual ~SIMDynElasticity() {}
 
+  //! \brief Returns the name of this simulator (for use in the HDF5 export).
+  virtual std::string getName() const
+  {
+    if (SIMElasticity<Dim>::myContext == "poroelasticity")
+      return "PoroElasticity";
+    else
+      return "Elasticity";
+  }
+
   //! \brief Prints out problem-specific data to the log stream.
   virtual void printProblem() const
   {
@@ -45,46 +54,25 @@ class SIMDynElasticity : public SIMElasticity<Dim>
     if (++ncall == 1) // Avoiding infinite recursive calls
       dSim.printProblem();
     else
-      this->SIMElasticity<Dim>::printProblem();
+      this->Dim::printProblem();
     --ncall;
   }
 
-  //! \brief Registers fields for data output.
-  void registerFields(DataExporter& exporter)
-  {
-    int flag = DataExporter::PRIMARY;
-    if (!Dim::opt.pSolOnly)
-      flag |= DataExporter::SECONDARY;
-    exporter.registerField("u","solid displacement",DataExporter::SIM,flag);
-    exporter.setFieldValue("u",this,&dSim.getSolution());
-  }
-
-  //! \brief Initializes the problem.
-  bool init(const TimeStep&)
+  //! \brief Initializes the solution vectors.
+  virtual bool init(const TimeStep&)
   {
     dSim.initPrm();
     dSim.initSol(3);
 
-    this->setMode(SIM::INIT);
+    bool ok = this->setMode(SIM::INIT);
     this->setQuadratureRule(Dim::opt.nGauss[0],true);
-    return true;
-  }
-
-  //! \brief Opens a new VTF-file and writes the model geometry to it.
-  //! \param[in] fileName File name used to construct the VTF-file name from
-  //! \param geoBlk Running geometry block counter
-  //! \param nBlock Running result block counter
-  bool saveModel(char* fileName, int& geoBlk, int& nBlock)
-  {
-    if (Dim::opt.format < 0) return true;
-
-    return dSim.saveModel(geoBlk,nBlock,fileName);
+    return ok;
   }
 
   //! \brief Saves the converged results of a given time step to VTF file.
   //! \param[in] tp Time stepping parameters
   //! \param nBlock Running result block counter
-  bool saveStep(const TimeStep& tp, int& nBlock)
+  virtual bool saveStep(const TimeStep& tp, int& nBlock)
   {
     double old = utl::zero_print_tol;
     utl::zero_print_tol = 1e-16;
@@ -110,10 +98,14 @@ class SIMDynElasticity : public SIMElasticity<Dim>
   }
 
   //! \brief Advances the time step one step forward.
-  virtual bool advanceStep(TimeStep& tp) { return dSim.advanceStep(tp,false); }
+  virtual bool advanceStep(TimeStep& tp)
+  {
+    dSim.advanceStep(tp,false);
+    return this->SIMElasticity<Dim>::advanceStep(tp);
+  }
 
   //! \brief Computes the solution for the current time step.
-  bool solveStep(TimeStep& tp)
+  virtual bool solveStep(TimeStep& tp)
   {
     if (Dim::msgLevel >= 1)
       IFEM::cout <<"\n  Solving the elasto-dynamics problem...";
@@ -164,7 +156,7 @@ class SIMDynElasticity : public SIMElasticity<Dim>
   //! \brief Returns the tensile energy in gauss points.
   const RealArray* getTensileEnergy() const
   {
-    return static_cast<FractureElasticity*>(Dim::myProblem)->getTensileEnergy();
+    return static_cast<Elasticity*>(Dim::myProblem)->getTensileEnergy();
   }
 
   //! \brief Returns a const reference to the global norms.
@@ -174,7 +166,7 @@ class SIMDynElasticity : public SIMElasticity<Dim>
   void setEnergyFile(const std::string&) {}
 
   //! \brief Returns a const reference to current solution vector.
-  const Vector& getSolution(int idx = 0) const { return dSim.getSolution(idx); }
+  virtual const Vector& getSolution(int i) const { return dSim.getSolution(i); }
   //! \brief Returns a const reference to the solution vectors.
   const Vectors& getSolutions() const { return dSim.getSolutions(); }
 
@@ -200,7 +192,7 @@ class SIMDynElasticity : public SIMElasticity<Dim>
   //! \brief Returns the actual integrand.
   virtual Elasticity* getIntegrand()
   {
-    if (!Dim::myProblem) // Using the Voigt formulation by default
+    if (!Dim::myProblem) // Using the Voigt elasticity formulation by default
       Dim::myProblem = new FractureElasticityVoigt(Dim::dimension);
     return static_cast<Elasticity*>(Dim::myProblem);
   }
@@ -212,15 +204,25 @@ class SIMDynElasticity : public SIMElasticity<Dim>
     static short int ncall = 0;
     if (++ncall == 1) // Avoiding infinite recursive calls
       result = dSim.parse(elem);
-    else if (!strcasecmp(elem->Value(),"elasticity") && !Dim::myProblem)
+    else if (!strcasecmp(elem->Value(),"elasticity"))
     {
       std::string form("voigt");
-      if (utl::getAttribute(elem,"formulation",form,true) && form != "voigt")
-        Dim::myProblem = new FractureElasticity(Dim::dimension);
+      if (!Dim::myProblem)
+        if (utl::getAttribute(elem,"formulation",form,true) && form != "voigt")
+          Dim::myProblem = new FractureElasticity(Dim::dimension);
+      SIMElasticity<Dim>::myContext = "elasticity";
       result = this->SIMElasticity<Dim>::parse(elem);
     }
-    else
+    else if (!strcasecmp(elem->Value(),"poroelasticity"))
+    {
+      if (!Dim::myProblem)
+        Dim::myProblem = new PoroElasticity(Dim::dimension);
+      SIMElasticity<Dim>::myContext = "poroelasticity";
       result = this->SIMElasticity<Dim>::parse(elem);
+    }
+    else
+      result = this->Dim::parse(elem);
+
     --ncall;
     return result;
   }

From 820360b0455d198707c652e913ec6804ecc0f9b3 Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Fri, 15 Apr 2016 17:29:39 +0200
Subject: [PATCH 3/4] Added: New integrand class PoroFracture

---
 FractureElasticity.C      |  8 +++++
 FractureElasticity.h      |  4 +++
 FractureElasticityVoigt.h |  7 +++-
 PoroFracture.C            | 61 +++++++++++++++++++++++++++++++++++
 PoroFracture.h            | 67 +++++++++++++++++++++++++++++++++++++++
 SIMDynElasticity.h        |  3 +-
 6 files changed, 148 insertions(+), 2 deletions(-)
 create mode 100644 PoroFracture.C
 create mode 100644 PoroFracture.h

diff --git a/FractureElasticity.C b/FractureElasticity.C
index 397fa58..3591563 100644
--- a/FractureElasticity.C
+++ b/FractureElasticity.C
@@ -34,6 +34,14 @@ FractureElasticity::FractureElasticity (unsigned short int n) : Elasticity(n)
 }
 
 
+FractureElasticity::FractureElasticity (IntegrandBase* parent,
+                                        unsigned short int n) : Elasticity(n)
+{
+  alpha = 0.0;
+  parent->registerVector("phasefield",&myCVec);
+}
+
+
 void FractureElasticity::initIntegration (size_t nGp, size_t)
 {
   // Initialize internal tensile energy buffer
diff --git a/FractureElasticity.h b/FractureElasticity.h
index 7e7dede..395443f 100644
--- a/FractureElasticity.h
+++ b/FractureElasticity.h
@@ -29,6 +29,10 @@ 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() {}
 
diff --git a/FractureElasticityVoigt.h b/FractureElasticityVoigt.h
index 7bfe2a7..f7cf838 100644
--- a/FractureElasticityVoigt.h
+++ b/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() {}
 
diff --git a/PoroFracture.C b/PoroFracture.C
new file mode 100644
index 0000000..f001727
--- /dev/null
+++ b/PoroFracture.C
@@ -0,0 +1,61 @@
+// $Id$
+//==============================================================================
+//!
+//! \file PoroFracture.C
+//!
+//! \date Apr 15 2016
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Integrand implementations for elasticity problems with fracture.
+//!
+//==============================================================================
+
+#include "PoroFracture.h"
+#include "FractureElasticityVoigt.h"
+
+
+PoroFracture::PoroFracture (unsigned short int n) : PoroElasticity(n)
+{
+  fracEl = new FractureElasticityVoigt(this,n);
+}
+
+
+PoroFracture::~PoroFracture()
+{
+  delete fracEl;
+}
+
+
+bool PoroFracture::parse (const TiXmlElement* elem)
+{
+  return this->PoroElasticity::parse(elem) & fracEl->parse(elem);
+}
+
+
+void PoroFracture::initIntegration (size_t nGp, size_t nBp)
+{
+  fracEl->initIntegration(nGp,nBp);
+}
+
+
+bool PoroFracture::initElement (const std::vector<int>& MNPC,
+                                LocalIntegral& elmInt)
+{
+  return this->PoroElasticity::initElement(MNPC,elmInt) &&
+         fracEl->initElement(MNPC,elmInt);
+}
+
+
+const RealArray* PoroFracture::getTensileEnergy () const
+{
+  return fracEl->getTensileEnergy();
+}
+
+
+bool PoroFracture::evalElasticityMatrices (ElmMats& elMat, const Matrix&,
+                                           const FiniteElement& fe,
+                                           const Vec3& X) const
+{
+  return fracEl->evalInt(elMat,fe,X);
+}
diff --git a/PoroFracture.h b/PoroFracture.h
new file mode 100644
index 0000000..c13c8e9
--- /dev/null
+++ b/PoroFracture.h
@@ -0,0 +1,67 @@
+// $Id$
+//==============================================================================
+//!
+//! \file PoroFracture.h
+//!
+//! \date Apr 15 2016
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Integrand implementations for elasticity problems with fracture.
+//!
+//==============================================================================
+
+#ifndef _PORO_FRACTURE_H
+#define _PORO_FRACTURE_H
+
+#include "PoroElasticity.h"
+
+class FractureElasticity;
+
+
+/*!
+  \brief Class representing the integrand of poroelasticity with fracture.
+
+  \details This class inherits PoroElasticity and uses elements from
+  FractureElasticity through a private member.
+*/
+
+class PoroFracture : public PoroElasticity
+{
+public:
+  //! \brief The constructor allocates the internal FractureElasticy object.
+  //! \param[in] n Number of spatial dimensions
+  PoroFracture(unsigned short int n);
+  //! \brief The destructor deletes the internal FractureElasticy object.
+  virtual ~PoroFracture();
+
+  //! \brief Parses a data section from an XML-element.
+  virtual bool parse(const TiXmlElement* elem);
+
+  //! \brief Initializes the integrand with the number of integration points.
+  //! \param[in] nGp Total number of interior integration points
+  //! \param[in] nBp Total number of boundary integration points
+  virtual void initIntegration(size_t nGp, size_t nBp);
+
+  //! \brief Initializes current element for numerical integration.
+  //! \param[in] MNPC Matrix of nodal point correspondance for current element
+  //! \param elmInt Local integral for element
+  virtual bool initElement(const std::vector<int>& MNPC, LocalIntegral& elmInt);
+
+  //! \brief Returns a pointer to the Gauss-point tensile energy array.
+  virtual const RealArray* getTensileEnergy() const;
+
+protected:
+  //! \brief Computes the elasticity matrices for a quadrature point.
+  //! \param elmInt The element matrix object to receive the contributions
+  //! \param[in] fe Finite element data of current integration point
+  //! \param[in] X Cartesian coordinates of current integration point
+  virtual bool evalElasticityMatrices(ElmMats& elMat, const Matrix&,
+                                      const FiniteElement& fe,
+                                      const Vec3& X) const;
+
+private:
+  FractureElasticity* fracEl; //!< Integrand for tangent stiffness evaluation
+};
+
+#endif
diff --git a/SIMDynElasticity.h b/SIMDynElasticity.h
index ec5f9f7..fac272a 100644
--- a/SIMDynElasticity.h
+++ b/SIMDynElasticity.h
@@ -18,6 +18,7 @@
 #include "NewmarkSIM.h"
 #include "SIMPoroElasticity.h"
 #include "FractureElasticityVoigt.h"
+#include "PoroFracture.h"
 
 
 /*!
@@ -216,7 +217,7 @@ class SIMDynElasticity : public SIMPoroElasticity<Dim>
     else if (!strcasecmp(elem->Value(),"poroelasticity"))
     {
       if (!Dim::myProblem)
-        Dim::myProblem = new PoroElasticity(Dim::dimension);
+        Dim::myProblem = new PoroFracture(Dim::dimension);
       SIMElasticity<Dim>::myContext = "poroelasticity";
       result = this->SIMElasticity<Dim>::parse(elem);
     }

From 183b73f59dc9cf37d9da7564d87b0492dfd663e4 Mon Sep 17 00:00:00 2001
From: Knut Morten Okstad <knut.morten.okstad@sintef.no>
Date: Tue, 19 Apr 2016 19:15:11 +0200
Subject: [PATCH 4/4] fixup! Added: New integrand class PoroFracture

---
 FractureElasticity.C      | 49 +++++++++++---------
 FractureElasticity.h      |  9 +++-
 FractureElasticityVoigt.C |  4 +-
 PoroFracture.C            | 94 ++++++++++++++++++++++++++++++++++++++-
 PoroFracture.h            | 32 +++++++++++--
 5 files changed, 158 insertions(+), 30 deletions(-)

diff --git a/FractureElasticity.C b/FractureElasticity.C
index 3591563..73c20b1 100644
--- a/FractureElasticity.C
+++ b/FractureElasticity.C
@@ -27,18 +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)
+                                        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
 }
 
 
@@ -52,7 +57,7 @@ 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;
@@ -60,20 +65,22 @@ bool FractureElasticity::initElement (const std::vector<int>& MNPC,
   }
 
   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;
 
@@ -219,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;
 }
@@ -262,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;
@@ -372,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)
@@ -397,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;
@@ -410,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;
   }
@@ -436,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;
 
diff --git a/FractureElasticity.h b/FractureElasticity.h
index 395443f..a09939f 100644
--- a/FractureElasticity.h
+++ b/FractureElasticity.h
@@ -36,6 +36,9 @@ class FractureElasticity : public Elasticity
   //! \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);
@@ -107,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
diff --git a/FractureElasticityVoigt.C b/FractureElasticityVoigt.C
index a93129e..220b11f 100644
--- a/FractureElasticityVoigt.C
+++ b/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;
 
diff --git a/PoroFracture.C b/PoroFracture.C
index f001727..c9d7d21 100644
--- a/PoroFracture.C
+++ b/PoroFracture.C
@@ -13,6 +13,7 @@
 
 #include "PoroFracture.h"
 #include "FractureElasticityVoigt.h"
+#include "Utilities.h"
 
 
 PoroFracture::PoroFracture (unsigned short int n) : PoroElasticity(n)
@@ -33,17 +34,105 @@ bool PoroFracture::parse (const TiXmlElement* elem)
 }
 
 
+void PoroFracture::setMode (SIM::SolutionMode mode)
+{
+  m_mode = mode;
+  fracEl->setMode(mode);
+  primsol.resize(fracEl->getNoSolutions());
+}
+
+
 void PoroFracture::initIntegration (size_t nGp, size_t nBp)
 {
   fracEl->initIntegration(nGp,nBp);
 }
 
 
+LocalIntegral* PoroFracture::getLocalIntegral (size_t nen,
+                                               size_t, bool neumann) const
+{
+  LocalIntegral* elmInt = this->PoroElasticity::getLocalIntegral(nen,0,neumann);
+  fracEl->setVar(nsd+1);
+  return elmInt;
+}
+
+
+LocalIntegral* PoroFracture::getLocalIntegral (const std::vector<size_t>& nen,
+                                               size_t, bool neumann) const
+{
+  LocalIntegral* elmInt = this->PoroElasticity::getLocalIntegral(nen,0,neumann);
+  fracEl->setVar(nsd);
+  return elmInt;
+}
+
+
+bool PoroFracture::initElement (const std::vector<int>& MNPC,
+                                LocalIntegral& elmInt)
+{
+  // Allocating three more vectors on the element level compared to global level
+  // (1 = pressure, 2 = pressure rate, 3 = phase field)
+  elmInt.vec.resize(primsol.size()+3);
+
+  // Extract element displacement and pressure vectors
+  if (!this->PoroElasticity::initElement(MNPC,elmInt))
+    return false;
+
+  // Extract element phase-field, velocity and acceleration vectors
+  if (!fracEl->initElement(MNPC,elmInt))
+    return false;
+
+  if (primsol.size() < 2)
+    return true; // Quasi-static simulation
+
+  // For the standard (non-mixed) formulation, the displacement and pressure
+  // variables are assumed stored interleaved in the element solution vector,
+  // now we need to separate them (for velocities and accelerations)
+  size_t uA = elmInt.vec.size() - 2; // Index to element acceleration vector
+  size_t uV = uA - 1;                // Index to element velocity vector
+  size_t pV = 2;                     // Index to element pressure rate vector
+  Matrix eMat(nsd+1,MNPC.size());
+  eMat = elmInt.vec[uV]; // Interleaved velocity vector
+  elmInt.vec[pV] = eMat.getRow(nsd+1);  // Assign pressure rate vector
+  elmInt.vec[uV] = eMat.expandCols(-1); // Assign velocity vector
+  eMat.resize(nsd+1,MNPC.size());
+  eMat = elmInt.vec[uA]; // Interleaved acceleration vector
+  elmInt.vec[uA] = eMat.expandCols(-1); // Assign acceleration vector
+  // We don't need the pressure acceleration
+
+  return true;
+}
+
+
 bool PoroFracture::initElement (const std::vector<int>& MNPC,
+                                const std::vector<size_t>& elem_sizes,
+                                const std::vector<size_t>& basis_sizes,
                                 LocalIntegral& elmInt)
 {
-  return this->PoroElasticity::initElement(MNPC,elmInt) &&
-         fracEl->initElement(MNPC,elmInt);
+  // Allocating three more vectors on the element level compared to global level
+  // (1 = pressure, 2 = pressure rate, 3 = phase field)
+  elmInt.vec.resize(primsol.size()+3);
+
+  // Extract element displacement and pressure vectors
+  if (!this->PoroElasticity::initElement(MNPC,elem_sizes,basis_sizes,elmInt))
+    return false;
+
+  // Extract element phase-field, velocity and acceleration vectors
+  std::vector<int>::const_iterator uEnd = MNPC.begin() + elem_sizes.front();
+  if (!fracEl->initElement(std::vector<int>(MNPC.begin(),uEnd),elmInt))
+    return false;
+
+  if (primsol.size() < 2)
+    return true; // Quasi-static simulation
+
+  // Extract the element level pressure rate vector
+  std::vector<int> MNPCp(uEnd,MNPC.end());
+  int ierr = utl::gather(MNPCp, 0,1, primsol[primsol.size()-2], elmInt.vec[2],
+                         nsd*basis_sizes.front(), basis_sizes.front());
+  if (ierr == 0) return true;
+
+  std::cerr <<" *** PoroFracture::initElement: Detected "<< ierr
+            <<" node numbers out of range."<< std::endl;
+  return false;
 }
 
 
@@ -57,5 +146,6 @@ bool PoroFracture::evalElasticityMatrices (ElmMats& elMat, const Matrix&,
                                            const FiniteElement& fe,
                                            const Vec3& X) const
 {
+  // Evaluate tangent stiffness matrix and internal forces
   return fracEl->evalInt(elMat,fe,X);
 }
diff --git a/PoroFracture.h b/PoroFracture.h
index c13c8e9..119206d 100644
--- a/PoroFracture.h
+++ b/PoroFracture.h
@@ -7,7 +7,7 @@
 //!
 //! \author Knut Morten Okstad / SINTEF
 //!
-//! \brief Integrand implementations for elasticity problems with fracture.
+//! \brief Integrand implementations for poroelasticity problems with fracture.
 //!
 //==============================================================================
 
@@ -21,9 +21,8 @@ class FractureElasticity;
 
 /*!
   \brief Class representing the integrand of poroelasticity with fracture.
-
   \details This class inherits PoroElasticity and uses elements from
-  FractureElasticity through a private member.
+  FractureElasticity in addition through a private member.
 */
 
 class PoroFracture : public PoroElasticity
@@ -38,15 +37,40 @@ class PoroFracture : public PoroElasticity
   //! \brief Parses a data section from an XML-element.
   virtual bool parse(const TiXmlElement* elem);
 
+  //! \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 Initializes the integrand with the number of integration points.
   //! \param[in] nGp Total number of interior integration points
   //! \param[in] nBp Total number of boundary integration points
   virtual void initIntegration(size_t nGp, size_t nBp);
 
+  //! \brief Returns a local integral contribution object for the given element.
+  //! \param[in] nen Number of nodes on element
+  //! \param[in] neumann Whether or not we are assembling Neumann BCs
+  virtual LocalIntegral* getLocalIntegral(size_t nen,
+                                          size_t, bool neumann) const;
+  //! \brief Returns a local integral contribution object for the given element.
+  //! \param[in] nen Number of nodes on element for each basis
+  //! \param[in] neumann Whether or not we are assembling Neumann BCs
+  virtual LocalIntegral* getLocalIntegral(const std::vector<size_t>& nen,
+                                          size_t, bool neumann) const;
+
+  using Elasticity::initElement;
   //! \brief Initializes current element for numerical integration.
   //! \param[in] MNPC Matrix of nodal point correspondance for current element
-  //! \param elmInt Local integral for element
+  //! \param elmInt The local integral object for current element
   virtual bool initElement(const std::vector<int>& MNPC, LocalIntegral& elmInt);
+  //! \brief Initializes current element for numerical integration.
+  //! \param[in] MNPC Matrix of nodal point correspondance for current element
+  //! \param[in] elem_sizes Size of each basis on the element
+  //! \param[in] basis_sizes Size of each basis on the patch level
+  //! \param elmInt The local integral object for current element
+  virtual bool initElement(const std::vector<int>& MNPC,
+                           const std::vector<size_t>& elem_sizes,
+                           const std::vector<size_t>& basis_sizes,
+                           LocalIntegral& elmInt);
 
   //! \brief Returns a pointer to the Gauss-point tensile energy array.
   virtual const RealArray* getTensileEnergy() const;