Added: Integrand for monolithic fracture elasticity solver

FractureElasticityMonol.C

@@ -0,0 +1,234 @@
+// $Id$
+//==============================================================================
+//!
+//! \file FractureElasticityMonol.C
+//!
+//! \date Jul 10 2016
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Integrand implementation for monolithic fracture elasticity problems.
+//!
+//==============================================================================
+
+#include "FractureElasticityMonol.h"
+#include "FiniteElement.h"
+#include "BlockElmMats.h"
+#include "ElmNorm.h"
+#include "Functions.h"
+#include "Utilities.h"
+#include "IFEM.h"
+#include "tinyxml.h"
+
+
+FractureElasticityMonol::FractureElasticityMonol (unsigned short int n)
+  : FractureElasticityVoigt(n)
+{
+  Gc = smearing = 1.0;
+  scale2nd = 4.0;
+  stabk = 0.0;
+}
+
+
+bool FractureElasticityMonol::parse (const TiXmlElement* elem)
+{
+  const char* value = utl::getValue(elem,"Gc");
+  if (value)
+    Gc = atof(value);
+  else if ((value = utl::getValue(elem,"smearing")))
+    smearing = atof(value);
+  else if ((value = utl::getValue(elem,"stabilization")))
+    stabk = atof(value);
+  else
+    return this->FractureElasticityVoigt::parse(elem);
+
+  return true;
+}
+
+
+void FractureElasticityMonol::printLog () const
+{
+  this->FractureElasticityVoigt::printLog();
+
+  IFEM::cout <<"\tCritical fracture energy density: "<< Gc
+             <<"\n\tSmearing factor: "<< smearing;
+  if (stabk != 0.0)
+    IFEM::cout <<"\n\tStabilization parameter: "<< stabk;
+  if (scale2nd == 2.0)
+    IFEM::cout <<"\n\tUsing fourth-order phase field.";
+  IFEM::cout << std::endl;
+}
+
+
+void FractureElasticityMonol::setMode (SIM::SolutionMode mode)
+{
+  m_mode = mode;
+
+  eM = eKm = eKg = eKc = eAcc = 0;
+  eS = iS = eBc = 0;
+
+  switch (mode)
+    {
+    case SIM::STATIC:
+      eKm  = eKg = 2;
+      eS   = iS  = 2;
+      eAcc = 3;
+      eBc  = 3;
+      eKc  = 4;
+      if (intPrm[3] > 0.0)
+        eKg = 0; // Linear analysis, no geometric stiffness
+      primsol.resize(nSV);
+      break;
+
+    case SIM::RHS_ONLY:
+      eS  = iS = 2;
+      eBc = 3;
+
+    case SIM::RECOVERY:
+      primsol.resize(1);
+      break;
+
+    default:
+      primsol.clear();
+    }
+}
+
+
+LocalIntegral* FractureElasticityMonol::getLocalIntegral (size_t nen, size_t,
+                                                          bool neumann) const
+{
+  BlockElmMats* result = new BlockElmMats(2);
+
+  switch (m_mode)
+  {
+    case SIM::STATIC:
+      result->rhsOnly = neumann;
+      result->withLHS = !neumann;
+      result->resize(neumann ? 0 : 4, 3);
+      break;
+
+    case SIM::RHS_ONLY:
+      result->resize(neumann ? 0 : 4, 3);
+
+    case SIM::RECOVERY:
+      result->rhsOnly = true;
+      result->withLHS = false;
+      break;
+
+    default:
+      ; // 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);
+  */
+  return result;
+}
+
+
+bool FractureElasticityMonol::initElement (const std::vector<int>& MNPC,
+                                           LocalIntegral& elmInt)
+{
+  if (mySol.empty())
+  {
+    std::cerr <<" *** FractureElasticityMonol::initElement:"
+              <<" No primary solution vectors."<< std::endl;
+    return false;
+  }
+  else if (m_mode == SIM::DYNAMIC)
+  {
+    std::cerr <<" *** FractureElasticityMonol::initElement:"
+              <<" Monolithic coupling is not available"
+              <<" for dynamic simulation."<< std::endl;
+    return false;
+  }
+
+  size_t nsol = 1 + eC;
+  if (elmInt.vec.size() < nsol)
+    elmInt.vec.resize(nsol);
+
+  // Extract the element level solution vectors from the monolithic solution
+  Matrix temp(npv,MNPC.size());
+  int ierr = utl::gather(MNPC,npv,primsol.front(),temp);
+  elmInt.vec[eC] = temp.getRow(npv);
+  elmInt.vec.front() = temp.expandRows(-1);
+
+  if (ierr == 0) return true;
+
+  std::cerr <<" *** FractureElasticityMonol::initElement: Detected "
+            << ierr <<" node numbers out of range."<< std::endl;
+  return false;
+}
+
+
+bool FractureElasticityMonol::evalInt (LocalIntegral& elmInt,
+                                       const FiniteElement& fe,
+                                       const Vec3& X) const
+{
+  if (!this->FractureElasticityVoigt::evalInt(elmInt,fe,X))
+    return false;
+
+  if (eAcc)
+  {
+    Matrix& A = static_cast<ElmMats&>(elmInt).A[eAcc-1];
+
+    double PhiPl = myPhi[fe.iGP];
+    double scale = 1.0 + 4.0*smearing*(1.0-stabk)*PhiPl/Gc;
+    double s1JxW = scale*fe.detJxW;
+    double s2JxW = scale2nd*smearing*smearing*fe.detJxW;
+
+    for (size_t i = 1; i <= fe.N.size(); i++)
+      for (size_t j = 1; j <= fe.N.size(); j++)
+      {
+        double grad = 0.0;
+        for (size_t k = 1; k <= nsd; k++)
+          grad += fe.dNdX(i,k)*fe.dNdX(j,k);
+        A(i,j) += fe.N(i)*fe.N(j)*s1JxW + grad*s2JxW;
+      }
+
+    if (scale2nd == 2.0)
+    {
+      // Forth-order phase field
+      double s4JxW = pow(smearing,4.0)*fe.detJxW;
+
+      for (size_t i = 1; i <= fe.N.size(); i++)
+        for (size_t j = 1; j <= fe.N.size(); j++)
+        {
+          double grad = 0.0;
+          for (unsigned short int k = 1; k <= nsd; k++)
+            grad += fe.d2NdX2(i,k,k)*fe.d2NdX2(j,k,k);
+          A(i,j) += grad*s4JxW;
+        }
+    }
+  }
+
+  if (eBc)
+    static_cast<ElmMats&>(elmInt).b[eBc-1].add(fe.N,fe.detJxW);
+
+  return true;
+}
+
+
+bool FractureElasticityMonol::evalSol (Vector& s,
+                                       const FiniteElement& fe, const Vec3& X,
+                                       const std::vector<int>& MNPC) const
+{
+  // Extract element displacements and phase field from the monolithic solution
+  Vectors eV(1+eC);
+  Matrix temp(npv,MNPC.size());
+  int ierr = utl::gather(MNPC,npv,mySol.front(),temp);
+  eV[eC] = temp.getRow(npv);
+  eV.front() = temp.expandRows(-1);
+
+  if (ierr > 0)
+  {
+    std::cerr <<" *** FractureElasticityMonol::evalSol: Detected "<< ierr
+              <<" node numbers out of range."<< std::endl;
+    return false;
+  }
+
+  return this->evalSol2(s,eV,fe,X);
+}

FractureElasticityMonol.h

@@ -0,0 +1,82 @@
+// $Id$
+//==============================================================================
+//!
+//! \file FractureElasticityMonol.h
+//!
+//! \date Jul 10 2016
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Integrand implementation for monolithic fracture elasticity problems.
+//!
+//==============================================================================
+
+#ifndef _FRACTURE_ELASTICITY_MONOL_H
+#define _FRACTURE_ELASTICITY_MONOL_H
+
+#include "FractureElasticityVoigt.h"
+
+
+/*!
+  \brief Class representing the integrand of elasticity problems with fracture.
+  \details This sub-class models the monolithic coupling to the phase field.
+*/
+
+class FractureElasticityMonol : public FractureElasticityVoigt
+{
+public:
+  //! \brief The constructor invokes the parent class constructor only.
+  //! \param[in] n Number of spatial dimensions
+  FractureElasticityMonol(unsigned short int n);
+  //! \brief Empty destructor.
+  virtual ~FractureElasticityMonol() {}
+
+  //! \brief Parses a data section from an XML element.
+  virtual bool parse(const TiXmlElement* elem);
+
+  //! \brief Prints out the problem definition to the log stream.
+  virtual void printLog() const;
+
+  //! \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);
+
+  using FractureElasticityVoigt::getLocalIntegral;
+  //! \brief Returns a local integral container for the given element.
+  //! \param[in] nen Number of nodes on element
+  //! \param[in] neumann Whether or not we are assembling Neumann BC's
+  virtual LocalIntegral* getLocalIntegral(size_t nen, size_t,
+                                          bool neumann) const;
+
+  //! \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 Evaluates the integrand at an interior point.
+  //! \param elmInt The local integral 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 evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
+                       const Vec3& X) const;
+
+  using FractureElasticityVoigt::evalSol;
+  //! \brief Evaluates the secondary solution at a result point.
+  //! \param[out] s Array of solution field values at current point
+  //! \param[in] fe Finite element data at current point
+  //! \param[in] X Cartesian coordinates of current point
+  //! \param[in] MNPC Nodal point correspondance for the basis function values
+  virtual bool evalSol(Vector& s, const FiniteElement& fe,
+                       const Vec3& X, const std::vector<int>& MNPC) const;
+
+protected:
+  unsigned short int eAcc; //!< Zero-based index to element phase-field matrix
+  unsigned short int eBc;  //!< Zero-based index to element phase-field vector
+
+  double Gc;       //!< Fracture energy density
+  double smearing; //!< Smearing factor in crack
+  double stabk;    //!< Stabilization parameter
+  double scale2nd; //!< Scaling factor in front of second order term
+};
+
+#endif