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/CMakeLists.txt b/CMakeLists.txt index 6b18706..1daa584 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -21,30 +21,54 @@ 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(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} \\ - ${PROJECT_SOURCE_DIR}/../Elasticity") - if(NOT TARGET Elasticity) - add_subdirectory(../Elasticity Elasticity) + ${PROJECT_SOURCE_DIR}/${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}) @@ -86,7 +110,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}) + ${Fracture_LIBRARIES}) if(IFEM_COMMON_APP_BUILD) set(TEST_APPS ${TEST_APPS} PARENT_SCOPE) diff --git a/FractureArgs.C b/FractureArgs.C index 2b715f7..7229a65 100644 --- a/FractureArgs.C +++ b/FractureArgs.C @@ -20,6 +20,7 @@ FractureArgs::FractureArgs () : SIMargsBase("fracturedynamics") { inpfile = nullptr; integrator = coupling = 1; + poroEl = expPhase = false; } @@ -41,6 +42,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 +83,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); diff --git a/FractureArgs.h b/FractureArgs.h index a59e649..ac0ef5c 100644 --- a/FractureArgs.h +++ b/FractureArgs.h @@ -30,6 +30,8 @@ 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 //! \brief Default constructor. FractureArgs(); diff --git a/FractureElasticity.C b/FractureElasticity.C index 3573515..38c226f 100644 --- a/FractureElasticity.C +++ b/FractureElasticity.C @@ -52,6 +52,7 @@ FractureElasticity::FractureElasticity (IntegrandBase* parent, parent->registerVector("phasefield",&myCVec); // Assuming second vector is pressure, third vector is pressure velocity eC = 3; // and fourth vector is the phase field + npv = nsd+1; // Account for pressure variables in the primary solution vector } @@ -121,6 +122,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 +522,11 @@ bool FractureElasticity::evalInt (LocalIntegral& elmInt, } -bool FractureElasticity::evalSol (Vector& s, - const FiniteElement& fe, const Vec3& X, - const std::vector& MNPC) const +bool FractureElasticity::getElementSolution (Vectors& eV, + const std::vector& 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 +535,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& MNPC) const +{ + Vectors eV(1+eC); + return this->getElementSolution(eV,MNPC) && this->evalSol2(s,eV,fe,X); } @@ -621,6 +641,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); diff --git a/FractureElasticity.h b/FractureElasticity.h index 66f7092..d1bc950 100644 --- a/FractureElasticity.h +++ b/FractureElasticity.h @@ -46,6 +46,10 @@ class FractureElasticity : public Elasticity //! \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 +100,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& 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; } diff --git a/PoroFracture.C b/PoroFracture.C new file mode 100644 index 0000000..79bc613 --- /dev/null +++ b/PoroFracture.C @@ -0,0 +1,344 @@ +// $Id$ +//============================================================================== +//! +//! \file PoroFracture.C +//! +//! \date Apr 15 2016 +//! +//! \author Knut Morten Okstad / SINTEF +//! +//! \brief Integrand implementations for poroelasticity problems with fracture. +//! +//============================================================================== + +#include "PoroFracture.h" +#include "PoroMaterial.h" +#include "FractureElasticityVoigt.h" +#include "FiniteElement.h" +#include "Tensor.h" +#include "Vec3Oper.h" +#include "Utilities.h" +#include "IFEM.h" +#include "tinyxml.h" + + +/*! + \brief Class representing the integrand of elasticity problems with fracture. + + \details This sub-class overrides the getElementSolution method, to account + for that the primary solution vector, as extracted from the patch level, + also contains the pressure variables in addition to the displacements. +*/ + +class PoroFractureElasticity : public FractureElasticityVoigt +{ +public: + //! \brief Constructor for integrands with a parent integrand. + //! \param parent The parent integrand of this one + //! \param[in] n Number of spatial dimensions + PoroFractureElasticity(IntegrandBase* parent, unsigned short int n) + : FractureElasticityVoigt(parent,n) {} + //! \brief Empty destructor. + virtual ~PoroFractureElasticity() {} + + //! \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& MNPC) const + { + if (!this->FractureElasticityVoigt::getElementSolution(eV,MNPC)) + return false; + + // Filter out the pressure components + // FIXME: Mixed + size_t nen = MNPC.size(); + if (eV.front().size() == nen*(nsd+1)) + { + Vector temp(eV.front()); + Vector& actual = eV.front(); + actual.resize(nen*nsd); + for (size_t n = 0; n < nen; n++) + for (unsigned short int i = 0; i < nsd; i++) + actual[nsd*n+i] = temp[(nsd+1)*n+i]; + } + + return true; + } +}; + + +PoroFracture::PoroFracture (unsigned short int n) : PoroElasticity(n) +{ + fracEl = new PoroFractureElasticity(this,n); + + L_per = 0.01; + d_min = 0.1; + Kc = 83.0; + eps = 50.0; +} + + +PoroFracture::~PoroFracture() +{ + delete fracEl; +} + + +bool PoroFracture::parse (const TiXmlElement* elem) +{ + if (strcasecmp(elem->Value(),"crack")) + return this->PoroElasticity::parse(elem) & fracEl->parse(elem); + + IFEM::cout <<"\tCrack parameters:"; + if (utl::getAttribute(elem,"Kc",Kc)) + IFEM::cout <<" Kc = "<< Kc; + if (utl::getAttribute(elem,"eps",eps)) + IFEM::cout <<" eps = "<< eps; + if (utl::getAttribute(elem,"L_per",L_per)) + IFEM::cout <<" L_per = "<< L_per; + if (utl::getAttribute(elem,"d_min",d_min)) + IFEM::cout <<" d_min = "<< d_min; + IFEM::cout << std::endl; + + return true; +} + + +Material* PoroFracture::parseMatProp (const TiXmlElement* elem, bool) +{ + this->PoroElasticity::parseMatProp(elem,true); + fracEl->setMaterial(material); + return material; +} + + +void PoroFracture::setMaterial (Material* mat) +{ + this->PoroElasticity::setMaterial(mat); + fracEl->setMaterial(mat); +} + + +void PoroFracture::setMode (SIM::SolutionMode mode) +{ + this->PoroElasticity::setMode(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& 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& 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, velocity, acceleration and pressure vectors + if (!this->PoroElasticity::initElement(MNPC,elmInt)) + return false; + + // Extract element phase-field vector + return fracEl->initElement(MNPC,elmInt); +} + + +bool PoroFracture::initElement (const std::vector& MNPC, + const std::vector& elem_sizes, + const std::vector& basis_sizes, + 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, velocity, acceleration and pressure vectors + if (!this->PoroElasticity::initElement(MNPC,elem_sizes,basis_sizes,elmInt)) + return false; + + // Extract element phase-field vector + std::vector::const_iterator uEnd = MNPC.begin() + elem_sizes.front(); + return fracEl->initElement(std::vector(MNPC.begin(),uEnd),elmInt); +} + + +const RealArray* PoroFracture::getTensileEnergy () const +{ + return fracEl->getTensileEnergy(); +} + + +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); +} + + +/*! + This method calculates the anisotropic permeability for the broken solid + based on a Poiseuille flow approximation of the fluid flow in the crack. + See Section 5.5.2 (eqs. (104)-(109)) of Miehe and Maute (2015) + "Phase field modeling of fracture in multi-physics problems. Part III." + as well as clarifying note by Fonn and Kvamsdal (2016). +*/ + +double PoroFracture::formCrackedPermeabilityTensor (SymmTensor& Kcrack, + const Vectors& eV, + const FiniteElement& fe, + const Vec3& X) const +{ + // Base permeability + // FIXME: What to do for non-isotropic permeability? + const PoroMaterial* pmat = dynamic_cast(material); + if (!pmat) + return -4.0; + double perm = pmat->getPermeability(X)[0]; + + Vec3 gradD; // Evaluate the phase field value and gradient + double d = fracEl->evalPhaseField(gradD,eV,fe.N,fe.dNdX); + if (d < 0.0) + { + std::cerr <<" *** PoroFracture::formCrackedPermeabilityTensor(X = "<< X + <<")\n Invalid phase field value: "<< d << std::endl; + return d; + } + else if (d < d_min) + { + // The crack does not affect the permeability tensor at this point + Kcrack = perm; + return 0.0; + } + + double d2 = gradD.length2(); + if (d2 <= 0.0) + { + std::cerr <<" *** PoroFracture::formCrackedPermeabilityTensor(X = "<< X + <<")\n Zero phase field gradient: "<< gradD << std::endl; + return -1.0; + } + + Tensor F(nsd); // Calculate the deformation gradient + if (!this->formDefGradient(eV.front(),fe.N,fe.dNdX,X.x,F)) + return -2.0; + + // Compute the inverse right Cauchy-Green tensor (C^-1) + if (Kcrack.rightCauchyGreen(F).inverse() == 0.0) + return -3.0; + + // Compute alpha = |grad D| / |F^-T grad D|, see F&K eq. (44) + Tensor Finv(F, true); + Finv.inverse(); + Vec3 FtgradD = Finv * gradD; + double alpha = gradD.length() / FtgradD.length(); + SymmTensor kCinv = perm * Kcrack; // k * C^-1 + + // Compute the symmetric tensor C^-1 - alpha * (C^-1*n0)otimes(C^-1*n0) + // (the term in the bracket [] of eq. (108) in Miehe2015pfm3) + // See also F&K eq. (45) + Vec3 CigD = Kcrack*gradD; // C^-1*gradD + double a2od2 = alpha*alpha/d2; + for (unsigned short int j = 1; j <= nsd; j++) + for (unsigned short int i = 1; i <= j; i++) + Kcrack(i,j) -= a2od2*CigD(i)*CigD(j); + + // Compute the perpendicular crack stretch + // lambda = gradD*gradD / gradD*C^-1*gradD (see M&M eq. (106), F&K eq. (36)) + double lambda = sqrt(d2 / (gradD*CigD)); +#if INT_DEBUG > 3 + std::cout <<"PoroFracture::formCrackedPermeabilityTensor(X = "<< X + <<") lambda = "<< lambda << std::endl; +#endif + if (lambda <= 1.0) + { + Kcrack = perm; + return 0.0; + } + + // Compute the permeability tensor, scale by d^eps*Kc*w^2*J (see eq. (108)) + // See also F&K eq. (45) + double w = lambda*L_per - L_per; // Crack opening (see M&M eq. (107)) + if (w < 0.0) w = 0.0; // See F&K eq. (37) + Kcrack *= pow(d,eps) * (w*w/12.0 - perm); + Kcrack += kCinv; + Kcrack *= F.det(); + + return w; +} + + +bool PoroFracture::formPermeabilityTensor (SymmTensor& K, + const Vectors& eV, + const FiniteElement& fe, + const Vec3& X) const +{ + return this->formCrackedPermeabilityTensor(K,eV,fe,X) >= 0.0; +} + + +size_t PoroFracture::getNoFields (int fld) const +{ + size_t nK = fld == 2 ? 1+nsd : 0; + return this->PoroElasticity::getNoFields(fld) + nK; +} + + +std::string PoroFracture::getField2Name (size_t i, const char* prefix) const +{ + if (i < this->PoroElasticity::getNoFields(2)) + return this->PoroElasticity::getField2Name(i,prefix); + + i -= this->PoroElasticity::getNoFields(2); + + static const char* s[4] = { "w", "K_xx", "K_yy", "K_zz" }; + + if (!prefix) return s[i%4]; + + return prefix + std::string(" ") + s[i%4]; +} + + +bool PoroFracture::evalSol (Vector& s, const FiniteElement& fe, + const Vec3& X, const std::vector& MNPC) const +{ + if (!this->PoroElasticity::evalSol(s,fe,X,MNPC)) + return false; + + Vectors eV; + if (!fracEl->getElementSolution(eV,MNPC)) + return false; + + SymmTensor Kc(nsd); + s.push_back(this->formCrackedPermeabilityTensor(Kc,eV,fe,X)); + for (size_t i = 1; i <= Kc.dim(); i++) + s.push_back(Kc(i,i)); + + return true; +} diff --git a/PoroFracture.h b/PoroFracture.h new file mode 100644 index 0000000..2915663 --- /dev/null +++ b/PoroFracture.h @@ -0,0 +1,142 @@ +// $Id$ +//============================================================================== +//! +//! \file PoroFracture.h +//! +//! \date Apr 15 2016 +//! +//! \author Knut Morten Okstad / SINTEF +//! +//! \brief Integrand implementations for poroelasticity 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 in addition 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); + + using PoroElasticity::parseMatProp; + //! \brief Parses material properties from an XML-element. + virtual Material* parseMatProp(const TiXmlElement* elem, bool); + + //! Defines the material properties. + virtual void setMaterial(Material* mat); + + //! \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& nen, + size_t, bool neumann) const; + + using PoroElasticity::initElement; + //! \brief Initializes current element for numerical integration. + //! \param[in] MNPC Matrix of nodal point correspondance for current element + //! \param elmInt The local integral object for current element + virtual bool initElement(const std::vector& 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& MNPC, + const std::vector& elem_sizes, + const std::vector& basis_sizes, + LocalIntegral& elmInt); + + using PoroElasticity::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& MNPC) const; + + //! \brief Returns a pointer to the Gauss-point tensile energy array. + virtual const RealArray* getTensileEnergy() const; + + //! \brief Returns the number of primary/secondary solution field components. + //! \param[in] fld Which field set to consider (1=primary,2=secondary) + virtual size_t getNoFields(int fld) const; + + //! \brief Returns the name of a secondary solution field component. + //! \param[in] i Field component index + //! \param[in] prefix Name prefix for all components + virtual std::string getField2Name(size_t i, const char* prefix) const; + +protected: + //! \brief Computes the elasticity matrices at 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; + + //! \brief Evaluates the permeability tensor at a quadrature point. + //! \param[out] K The permeability tensor + //! \param[in] eV Element solution vectors + //! \param[in] fe Finite element data of current integration point + //! \param[in] X Cartesian coordinates of current integration point + virtual bool formPermeabilityTensor(SymmTensor& K, + const Vectors& eV, + const FiniteElement& fe, + const Vec3& X) const; + + //! \brief Computes the permeability tensor of the broken material. + //! \param[out] Kcrack Permeability tensor of the broken material + //! \param[in] eV Element solution vectors + //! \param[in] fe Finite element data of current integration point + //! \param[in] X Cartesian coordinates of current integration point + //! \return Estimated crack opening + double formCrackedPermeabilityTensor(SymmTensor& Kcrack, + const Vectors& eV, + const FiniteElement& fe, + const Vec3& X) const; + +private: + FractureElasticity* fracEl; //!< Integrand for tangent stiffness evaluation + + double L_per; //!< Characteristic length of crack-perpendicular line + double d_min; //!< Crack phase field value for incorporating Poiseuille flow + double Kc; //!< Spatial permeability in fracture + double eps; //!< Permeability transition exponent +}; + +#endif diff --git a/README.md b/README.md index dbd8c33..960830e 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ ## Introduction -This module contains Fracture Dynamics libraries and applications built +This module contains the Fracture Dynamics library and applications built using the IFEM library. ### Getting all dependencies @@ -12,32 +12,35 @@ using the IFEM library. ### Getting the code -This is done by first navigating to the folder in which you want IFEM installed and typing +This is done by first navigating to a folder `` in which you want +the application and typing git clone https://github.com/OPM/IFEM-Elasticity + git clone https://github.com/OPM/IFEM-PoroElasticity git clone https://github.com/OPM/IFEM-OpenFrac The build system uses sibling directory logic to locate the IFEM-Elasticity -module. +and IFEM-PoroElasticity modules. ### Compiling the code To compile, first navigate to the root catalogue ``. -1. `cd ` +1. `cd IFEM-OpenFrac` 2. `mkdir Debug` 3. `cd Debug` 5. `cmake -DCMAKE_BUILD_TYPE=Debug ..` -6. `make ` +6. `make` -this will compile the library and the fracture dynamics applications. -The binaries can be found in the 'bin' subfolder. -Change all instances of `Debug` with `Release` to drop debug-symbols, -but get faster running code. +This will compile the library and applications. +The executables can be found in the 'bin' sub-folder. +Change all instances of `Debug` with `Release` to drop debug-symbols, +and get a faster running code. ### Testing the code -IFEM is using cmake test system. To compile run all regression- and unit-tests, navigate to your build -folder (i.e. `/Debug`) and type +IFEM uses the cmake test system. +To compile and run all regression- and unit-tests, navigate to your build +folder (i.e. `/IFEM-OpenFrac/Debug`) and type make check diff --git a/SIMDynElasticity.h b/SIMDynElasticity.h index 04b8449..badc32e 100644 --- a/SIMDynElasticity.h +++ b/SIMDynElasticity.h @@ -15,9 +15,11 @@ #define _SIM_DYN_ELASTICITY_H_ #include "NewmarkSIM.h" -#include "SIMElasticity.h" +#include "SIMElasticityWrap.h" #include "FractureElasticityVoigt.h" -#include "DataExporter.h" +#ifdef IFEM_HAS_POROELASTIC +#include "PoroFracture.h" +#endif #include @@ -25,14 +27,19 @@ \brief Driver class for dynamic elasticity problems with fracture. */ -template -class SIMDynElasticity : public SIMElasticity +template< class Dim, class DynSIM=NewmarkSIM, class Sim=SIMElasticityWrap > +class SIMDynElasticity : public Sim { public: //! \brief Default constructor. SIMDynElasticity() : dSim(*this) { - Dim::myHeading = "Elasticity solver"; + vtfStep = 0; + } + + //! \brief Constructor for mixed problems. + SIMDynElasticity(const std::vector& nf) : Sim(nf), dSim(*this) + { vtfStep = 0; } @@ -46,22 +53,12 @@ class SIMDynElasticity : public SIMElasticity if (++ncall == 1) // Avoiding infinite recursive calls dSim.printProblem(); else - this->SIMElasticity::printProblem(); + this->Sim::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&) + virtual bool init(const TimeStep&) { dSim.initPrm(); dSim.initSol(dynamic_cast(&dSim) ? 3 : 1); @@ -72,21 +69,10 @@ class SIMDynElasticity : public SIMElasticity return this->setInitialConditions() && ok; } - //! \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); - } - //! \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; @@ -143,7 +129,7 @@ class SIMDynElasticity : public SIMElasticity virtual bool advanceStep(TimeStep& tp) { return dSim.advanceStep(tp,false); } //! \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..."; @@ -211,7 +197,7 @@ class SIMDynElasticity : public SIMElasticity } //! \brief Returns the tensile energy in gauss points. - const RealArray* getTensileEnergy() const + virtual const RealArray* getTensileEnergy() const { return static_cast(Dim::myProblem)->getTensileEnergy(); } @@ -233,7 +219,7 @@ class SIMDynElasticity : public SIMElasticity } //! \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(); } @@ -287,6 +273,12 @@ class SIMDynElasticity : public SIMElasticity return fel ? fel->getCrackPressure() != nullptr : false; } + //! \brief Add additional MADOF for explicit phase field + void addMADOF() + { + this->Dim::addMADOF(1, 1, false); + } + protected: //! \brief Returns the actual integrand. virtual Elasticity* getIntegrand() @@ -304,12 +296,25 @@ class SIMDynElasticity : public SIMElasticity 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(),SIMElasticity::myContext.c_str())) { - std::string form("voigt"); - if (utl::getAttribute(elem,"formulation",form,true) && form != "voigt") - Dim::myProblem = new FractureElasticity(Dim::dimension); - result = this->SIMElasticity::parse(elem); + if (!Dim::myProblem) + { + if (this->getName() == "PoroElasticity") +#ifdef IFEM_HAS_POROELASTIC + Dim::myProblem = new PoroFracture(Dim::dimension); +#else + return false; +#endif + else + { + std::string formulation("voigt"); + utl::getAttribute(elem,"formulation",formulation,true); + if (formulation != "voigt") + Dim::myProblem = new FractureElasticity(Dim::dimension); + } + } + result = this->Sim::parse(elem); const TiXmlElement* child = elem->FirstChildElement(); for (; child; child = child->NextSiblingElement()) @@ -317,7 +322,7 @@ class SIMDynElasticity : public SIMElasticity this->setIntegrationPrm(3,1); // Disable geometric stiffness } else - result = this->SIMElasticity::parse(elem); + result = this->Dim::parse(elem); --ncall; return result; } diff --git a/SIMExplPhaseField.C b/SIMExplPhaseField.C new file mode 100644 index 0000000..0c6ee18 --- /dev/null +++ b/SIMExplPhaseField.C @@ -0,0 +1,86 @@ +// $Id$ +//============================================================================== +//! +//! \file SIMExplPhaseField.C +//! +//! \date May 27 2016 +//! +//! \author Knut Morten Okstad / SINTEF +//! +//! \brief Solution driver representing an explicit phase-field. +//! +//============================================================================== + +#include "SIMExplPhaseField.h" +#include "SIMoutput.h" +#include "TimeStep.h" +#include "DataExporter.h" +#include "Functions.h" +#include "Utilities.h" +#include "IFEM.h" +#include "tinyxml.h" + + +SIMExplPhaseField::SIMExplPhaseField (SIMoutput* gridOwner) +{ + myHeading = "Explicit phase field"; + myOwner = gridOwner; + phaseFunc = nullptr; + myStep = 0; +} + + +void SIMExplPhaseField::registerFields (DataExporter& exporter) +{ + exporter.registerField("c","phase field",DataExporter::SIM, + DataExporter::PRIMARY); + exporter.setFieldValue("c",this,&phaseField); +} + + +bool SIMExplPhaseField::init (const TimeStep&) +{ + this->registerField("phasefield",phaseField); + return true; +} + + +bool SIMExplPhaseField::parse (const TiXmlElement* elem) +{ + const char* value = utl::getValue(elem,"phasefield"); + if (!value) + return this->SIMbase::parse(elem); + + std::string type; + utl::getAttribute(elem,"type",type); + IFEM::cout <<"\tPhase-field function"; + phaseFunc = utl::parseRealFunc(value,type); + IFEM::cout << std::endl; + + return true; +} + + +bool SIMExplPhaseField::solveStep (TimeStep& tp, bool) +{ + if (!phaseFunc) + { + std::cerr <<" *** SIMExplPhaseField::solveStep: No phase field function." + << std::endl; + return false; + } + + phaseField.resize(myOwner->getNoNodes(1)); + return myOwner->project(phaseField,phaseFunc,1,0,1, + SIMoptions::GLOBAL,tp.time.t); +} + + +bool SIMExplPhaseField::saveStep (const TimeStep& tp, int& nBlock) +{ + if (tp.step%opt.saveInc == 0 && opt.format >= 0) + if (myOwner->writeGlvS1(phaseField,++myStep,nBlock,tp.time.t, + "phase",6,1,true) < 0) return false; + + return true; +} diff --git a/SIMExplPhaseField.h b/SIMExplPhaseField.h new file mode 100644 index 0000000..a52e55d --- /dev/null +++ b/SIMExplPhaseField.h @@ -0,0 +1,106 @@ +// $Id$ +//============================================================================== +//! +//! \file SIMExplPhaseField.h +//! +//! \date May 27 2016 +//! +//! \author Knut Morten Okstad / SINTEF +//! +//! \brief Solution driver representing an explicit phase-field. +//! +//============================================================================== + +#ifndef _SIM_EXPL_PHASE_FIELD_H +#define _SIM_EXPL_PHASE_FIELD_H + +#include "SIMbase.h" +#include "SIMdummy.h" +#include "SIMenums.h" + +class SIMoutput; +class DataExporter; +class TimeStep; +class VTF; +namespace LR { class LRSpline; class RefineData; } + + +/*! + \brief Driver class for an explicit phase-field. +*/ + +class SIMExplPhaseField : public SIMdummy +{ +public: + //! \brief Default constructor. + SIMExplPhaseField(SIMoutput* gridOwner = nullptr); + //! \brief Empty destructor. + virtual ~SIMExplPhaseField() {} + + //! \brief Registers fields for data output. + void registerFields(DataExporter& exporter); + + //! \brief Initializes the problem. + bool init(const TimeStep&); + + //! \brief Saves the converged results of a given time step to VTF file. + bool saveStep(const TimeStep& tp, int& nBlock); + + //! \brief Computes the solution for the current time step. + bool solveStep(TimeStep& tp, bool = true); + + //! \brief Dummy method. + bool postSolve(TimeStep&) { return true; } + //! \brief Dummy method. + bool advanceStep(TimeStep&) { return true; } + //! \brief Dummy method. + bool serialize(std::map&) { return false; } + //! \brief Dummy method. + bool dumpGeometry(std::ostream& os) const { return false; } + //! \brief Dummy method. + bool saveResidual(const TimeStep&, const Vector&, int&) { return true; } + //! \brief Dummy method. + void setTensileEnergy(const RealArray*) {} + //! \brief Dummy method. + void setVTF(VTF*) {} + //! \brief Dummy method. + double getNorm(Vector&, size_t = 0) const { return 0.0; } + //! \brief Dummy method. + const Vector& getGlobalNorms() const { static Vector g; return g; } + //! \brief Returns a const reference to the current solution. + const Vector& getSolution(int = 0) const { return phaseField; } + //! \brief Updates the solution vector. + void setSolution(const Vector& vec) { phaseField = vec; } + //! \brief Returns the maximum number of iterations (unlimited). + int getMaxit() const { return 9999; } + //! \brief Dummy method. + int getInitRefine() const { return 0; } + //! \brief Dummy method. + SIM::ConvStatus solveIteration(TimeStep&) { return SIM::CONVERGED; } + //! \brief Dummy method. + Vector getHistoryField() const { return Vector(); } + //! \brief Dummy method. + void setHistoryField(const RealArray&) {} + //! \brief Dummy method. + bool refine(const LR::RefineData&) { return false; } + //! \brief Dummy method. + bool refine(const LR::RefineData&,Vector&) { return false; } + //! \brief Dummy method. + void getBasis(std::vector&) {} + //! \brief Dummy method. + bool transferHistory(const RealArray&, + std::vector&) { return true; } + +protected: + using SIMbase::parse; + //! \brief Parses a data section from an XML element. + virtual bool parse(const TiXmlElement* elem); + +private: + SIMoutput* myOwner; //!< The FE mesh holder + int myStep; //!< VTF file step counter + RealFunc* phaseFunc; //!< Explicit phase field function + Vector phaseField; //!< Nodal phase field values +}; + +#endif diff --git a/SIMFractureDynamics.h b/SIMFractureDynamics.h index 18f3492..0149588 100644 --- a/SIMFractureDynamics.h +++ b/SIMFractureDynamics.h @@ -50,6 +50,10 @@ class SIMFracture : public Coupling virtual void setupDependencies() { this->S1.registerDependency(&this->S2,"phasefield",1); + + if (this->S2.getHeading().find("Explicit") != std::string::npos) + this->S1.addMADOF(); + // The tensile energy is defined on integration points and not nodal points. // It is a global buffer array across all patches in the model. // Use an explicit call instead of normal couplings for this. diff --git a/Test/Miehe71.xinp b/Test/Miehe71.xinp new file mode 100644 index 0000000..cfc428f --- /dev/null +++ b/Test/Miehe71.xinp @@ -0,0 +1,66 @@ + + + + + + + + + 1 + + + 2 + + + 3 4 + + + + + + 1.0 + 0.1 + abs(x-2.5) + 1 + + + + C=1.0-exp(-0.08)/1.08; + l=0.25; xx=abs(x-2.5); xi=2.0*(xx-l)/l; + if(above(xx,l),1.0-exp(-xi)/(1.0+xi),C*xx/l) + + + + + + + + 3.0 10.0 + + + + + + + + + 3.0 10.0 + + + + + 2 + + + + + constant displacement + + + + 1.0 + 122 + + + diff --git a/main_FractureDynamics.C b/main_FractureDynamics.C index d09804a..6b81cc0 100644 --- a/main_FractureDynamics.C +++ b/main_FractureDynamics.C @@ -16,7 +16,11 @@ #include "SIM3D.h" #include "SIMDynElasticity.h" #include "SIMPhaseField.h" +#include "SIMExplPhaseField.h" #include "SIMFractureQstatic.h" +#ifdef IFEM_HAS_POROELASTIC +#include "SIMPoroElasticity.h" +#endif #include "SIMCoupledSI.h" #include "SIMSolver.h" #include "SIMSolverTS.h" @@ -26,6 +30,7 @@ #include "NewmarkNLSIM.h" #include "NonLinSIM.h" #include "ASMstruct.h" +#include "ASMmxBase.h" /*! @@ -136,15 +141,11 @@ int runCombined (char* infile, const char* context) \brief Determines whether the quasi-static semi-implicit driver is to be used. */ -template class Cpl, template class Solver=SIMSolver> -int runSimulator3 (const FractureArgs& args) +int runSimulator6 (const FractureArgs& args, const char* contx) { - typedef SIMDynElasticity ElSolver; - typedef SIMPhaseField PhaseSolver; - - const char* contx = Integrator::inputContext; if (args.integrator == 3 && args.coupling == 2) { typedef SIMFractureQstatic Coupler; @@ -163,16 +164,32 @@ int runSimulator3 (const FractureArgs& args) } +/*! + \brief Determines whether the explicit phase-field driver is to be used. +*/ + +template class Cpl, + template class Solver=SIMSolver> +int runSimulator5 (const FractureArgs& args, const char* context) +{ + if (args.expPhase) + return runSimulator6(args,context); + else + return runSimulator6,Cpl,Solver>(args,context); +} + + /*! \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 */ -template +template int runStandAlone (char* infile, const char* context) { - typedef SIMDynElasticity SIMElastoDynamics; + typedef SIMDynElasticity SIMElastoDynamics; utl::profiler->start("Model input"); IFEM::cout <<"\n\n0. Parsing input file(s)." @@ -211,17 +228,62 @@ int runStandAlone (char* infile, const char* context) } +/*! + \brief Determines whether the poroelastic simulation driver is to be used. +*/ + +template +int runSimulator4 (const FractureArgs& args, + const char* context = "newmarksolver") +{ + if (args.poroEl) +#ifdef IFEM_HAS_POROELASTIC + return runStandAlone>(args.inpfile, + context); +#else + return 99; // Built without the poroelastic coupling +#endif + + return runStandAlone>(args.inpfile, + context); +} + + /*! \brief Determines whether the adaptive simulation driver is to be used. */ -template class Cpl> -int runSimulator2 (const FractureArgs& args) +template class Cpl> +int runSimulator3 (const FractureArgs& args) { + typedef SIMDynElasticity DynElSolver; + + const char* context = Integrator::inputContext; + if (args.adap) - return runSimulator3(args); + return runSimulator5(args,context); + + return runSimulator5(args,context); +} - return runSimulator3(args); + +/*! + \brief Creates the combined fracture simulator and launches the simulation. +*/ + +template class Cpl> +int runSimulator2 (const FractureArgs& args) +{ + if (args.poroEl) +#ifdef IFEM_HAS_POROELASTIC + return runSimulator3,Cpl>(args); +#else + return 99; // Built without the poroelastic coupling +#endif + + return runSimulator3,Cpl>(args); } @@ -239,7 +301,7 @@ int runSimulator1 (const FractureArgs& args) case 2: return runSimulator2(args); default: // No phase field coupling - return runStandAlone(args.inpfile,Integrator::inputContext); + return runSimulator4(args,Integrator::inputContext); } } @@ -267,7 +329,7 @@ int runSimulator (const FractureArgs& args) { switch (args.integrator) { case 0: - return runStandAlone(args.inpfile,"staticsolver"); + return runSimulator4(args,"staticsolver"); case 1: return runSimulator1(args); case 2: @@ -295,6 +357,7 @@ int main (int argc, char** argv) FractureArgs args; SIMElasticity::planeStrain = true; + ASMmxBase::Type = ASMmxBase::NONE; IFEM::Init(argc,argv,"Fracture dynamics solver"); for (int i = 1; i < argc; i++) @@ -302,6 +365,8 @@ int main (int argc, char** argv) ; // ignore the input file on the second pass else if (SIMoptions::ignoreOldOptions(argc,argv,i)) ; // ignore the obsolete option + else if (!strcmp(argv[i],"-mixed")) + ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1; else if (!strcmp(argv[i],"-principal")) Elasticity::wantPrincipalStress = true; else if (!strcmp(argv[i],"-dbgElm") && i < argc-1) @@ -315,8 +380,9 @@ int main (int argc, char** argv) { std::cout <<"usage: "<< argv[0] <<" [-dense|-spr|-superlu[]|-samg|-petsc]\n" - <<" [-lag|-spec|-LR] [-2D] [-nGauss ]\n " - <<"[-nocrack|-semiimplicit] [-[l|q]static|-GA|-HHT] [-adaptive]\n" + <<" [-lag|-spec|-LR] [-2D] [-mixed] [-nGauss ]\n" + <<" [-nocrack|-explcrack|-semiimplicit]" + <<" [-[l|q]static|-GA|-HHT] [-poro] [-adaptive]\n" <<" [-vtf [-nviz ] [-nu ] [-nv ]]\n [-hdf5] [-principal]\n"; return 0;