FractureElasticityMonol.C

@@ -14,8 +14,6 @@
 #include "FractureElasticityMonol.h"
 #include "FiniteElement.h"
 #include "BlockElmMats.h"
-#include "ElmNorm.h"
-#include "Functions.h"
 #include "Utilities.h"
 #include "IFEM.h"
 #include "tinyxml.h"
@@ -119,12 +117,11 @@ LocalIntegral* FractureElasticityMonol::getLocalIntegral (size_t nen, size_t,
       ; // Other solution modes not supported
   }
 
-  /*TODO
-  result->redim(npv*nen,0);
-  result->redim(nsd*nen,1);
-  result->redim(nen,nen,2);
-  result->redim(nsd*nen,nen,eKc-1);
-  */
+  result->redim(eKm-1,nen,nsd);
+  result->redim(eAcc-1,nen,1);
+  result->redimOffDiag(eKc-1);
+  result->redimNewtonMat();
+
   return result;
 }
 

SIMDynElasticity.h

@@ -3,12 +3,11 @@
 //!
 //! \file SIMDynElasticity.h
 //!
-//! \date Jul 13 2015
+//! \date Dec 04 2015
 //!
-//! \author Arne Morten Kvarving / SINTEF
+//! \author Knut Morten Okstad / SINTEF
 //!
-//! \brief Wrapper equipping the elasticity solver with
-//! time-stepping support and phase field coupling.
+//! \brief Dynamic simulation driver for elasticity problems with fracture.
 //!
 //==============================================================================
 
@@ -17,22 +16,24 @@
 
 #include "NewmarkSIM.h"
 #include "SIMElasticity.h"
-#include "FractureElasticityVoigt.h"
+#include "FractureElasticityMonol.h"
 #include "DataExporter.h"
 
 
 /*!
-  \brief Driver wrapping an elasticity solver with an ISolver interface.
+  \brief Driver class for dynamic elasticity problems with fracture.
 */
 
 template<class Dim, class DynSIM=NewmarkSIM>
 class SIMDynElasticity : public SIMElasticity<Dim>
 {
 public:
   //! \brief Default constructor.
-  SIMDynElasticity() : SIMElasticity<Dim>(false), dSim(*this), vtfStep(0)
+  SIMDynElasticity(bool mon = false) : SIMElasticity<Dim>(false), dSim(*this)
   {
     Dim::myHeading = "Elasticity solver";
+    monolithic = mon;
+    vtfStep = 0;
   }
 
   //! \brief Empty destructor.
@@ -174,9 +175,9 @@ class SIMDynElasticity : public SIMElasticity<Dim>
   }
 
   //! \brief Returns the tensile energy in gauss points.
-  const RealArray* getTensileEnergy() const
+  virtual 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.
@@ -212,8 +213,13 @@ class SIMDynElasticity : public SIMElasticity<Dim>
   //! \brief Returns the actual integrand.
   virtual Elasticity* getIntegrand()
   {
-    if (!Dim::myProblem) // Using the Voigt formulation by default
-      Dim::myProblem = new FractureElasticityVoigt(Dim::dimension);
+    if (!Dim::myProblem)
+    {
+      if (monolithic)
+        Dim::myProblem = new FractureElasticityMonol(Dim::dimension);
+      else // Using the Voigt elasticity formulation by default
+        Dim::myProblem = new FractureElasticityVoigt(Dim::dimension);
+    }
     return static_cast<Elasticity*>(Dim::myProblem);
   }
 
@@ -224,11 +230,20 @@ 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(),"cahnhilliard") && monolithic)
+    {
+      const TiXmlElement* child = elem->FirstChildElement();
+      for (; child; child = child->NextSiblingElement())
+        this->getIntegrand()->parse(child);
+    }
     else if (!strcasecmp(elem->Value(),"elasticity") && !Dim::myProblem)
     {
-      std::string form("voigt");
-      if (utl::getAttribute(elem,"formulation",form,true) && form != "voigt")
-        Dim::myProblem = new FractureElasticity(Dim::dimension);
+      if (!monolithic)
+      {
+        std::string form("voigt");
+        if (utl::getAttribute(elem,"formulation",form,true) && form != "voigt")
+          Dim::myProblem = new FractureElasticity(Dim::dimension);
+      }
       result = this->SIMElasticity<Dim>::parse(elem);
     }
     else
@@ -242,7 +257,9 @@ class SIMDynElasticity : public SIMElasticity<Dim>
   Matrix projSol; //!< Projected secondary solution fields
   Matrix eNorm;   //!< Element norm values
   Vector gNorm;   //!< Global norm values
-  int    vtfStep; //!< VTF file step counter
+
+  int  vtfStep;    //!< VTF file step counter
+  bool monolithic; //!< If \e true, use the monolithic integrand
 };
 
 #endif

main_FractureDynamics.C

@@ -157,18 +157,19 @@ int runSimulator3 (char* infile)
   \brief Creates and launches a stand-alone elasticity simulator (no coupling).
   \param[in] infile The input file to parse
   \param[in] context Input-file context for the time integrator
+  \param[in] monolithic If \e true, use monolithic coupling to phase-field
 */
 
 template<class Dim, class Integrator=NewmarkSIM>
-int runStandAlone (char* infile, const char* context)
+int runStandAlone (char* infile, const char* context, bool monolithic = false)
 {
   typedef SIMDynElasticity<Dim,Integrator> SIMElastoDynamics;
 
   utl::profiler->start("Model input");
   IFEM::cout <<"\n\n0. Parsing input file(s)."
              <<"\n========================="<< std::endl;
 
-  SIMElastoDynamics elastoSim;
+  SIMElastoDynamics elastoSim(monolithic);
   if (!elastoSim.read(infile))
     return 1;
 
@@ -231,6 +232,9 @@ int runSimulator1 (const FDargs& args)
     return runSimulator2<Dim,Integrator,SIMCoupled>(args);
   case 2:
     return runSimulator2<Dim,Integrator,SIMCoupledSI>(args);
+  case 3: // Monolithic coupling, single integrand
+    return runStandAlone<Dim,Integrator>(args.inpfile,
+                                         Integrator::inputContext,true);
   default: // No phase field coupling
     return runStandAlone<Dim,Integrator>(args.inpfile,Integrator::inputContext);
   }
@@ -300,6 +304,8 @@ int main (int argc, char** argv)
       args.coupling = 0;
     else if (!strcmp(argv[i],"-semiimplicit"))
       args.coupling = 2;
+    else if (!strcmp(argv[i],"-monolithic"))
+      args.coupling = 3;
     else if (!strcmp(argv[i],"-lstatic"))
       args.integrator = 0;
     else if (!strcmp(argv[i],"-GA"))
@@ -325,10 +331,11 @@ int main (int argc, char** argv)
   {
     std::cout <<"usage: "<< argv[0]
               <<" <inputfile> [-dense|-spr|-superlu[<nt>]|-samg|-petsc]\n"
-              <<"       [-lag|-spec|-LR] [-2D] [-nGauss <n>]\n       "
-              <<"[-nocrack|-semiimplicit] [-[l|q]static|-GA|-HHT] [-adaptive]\n"
+              <<"       [-lag|-spec|-LR] [-2D] [-nGauss <n>]\n"
+              <<"       [-nocrack|-semiimplicit|-monolithic]"
+              <<" [-[l|q]static|-GA|-HHT] [-adaptive]\n"
               <<"       [-vtf <format> [-nviz <nviz>] [-nu <nu>] [-nv <nv]"
-              <<" [-nw <nw>]] [-hdf5] [-principal]\n"<< std::endl;
+              <<" [-nw <nw>]]\n       [-hdf5] [-principal]\n"<< std::endl;
     return 0;
   }