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 1f24de8..aef2af6 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,14 @@ endif() file(GLOB FracEl_SRCS [A-Z]*.C) add_library(CommonFrac STATIC ${FracEl_SRCS}) +set(Common_LIBRARIES CommonFrac Elasticity ${IFEM_LIBRARIES}) +set(Fracture_LIBRARIES CommonFrac PoroElastic 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 ${Fracture_LIBRARIES}) # Installation install(TARGETS CahnHilliard FractureDynamics DESTINATION bin) @@ -79,7 +92,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/FractureElasticity.C b/FractureElasticity.C index 397fa58..7e5ef08 100644 --- a/FractureElasticity.C +++ b/FractureElasticity.C @@ -27,10 +27,34 @@ #endif -FractureElasticity::FractureElasticity (unsigned short int n) : Elasticity(n) +FractureElasticity::FractureElasticity (unsigned short int n) + : Elasticity(n), mySol(primsol) { alpha = 0.0; this->registerVector("phasefield",&myCVec); + eC = 1; // Assuming second vector is phase field +} + + +FractureElasticity::FractureElasticity (IntegrandBase* parent, + unsigned short int n) + : Elasticity(n), mySol(parent->getSolutions()) +{ + alpha = 0.0; + parent->registerVector("phasefield",&myCVec); + // Assuming second vector is pressure, third vector is pressure velocity + eC = 3; // and fourth vector is the phase field +} + + +void FractureElasticity::setMode (SIM::SolutionMode mode) +{ + this->Elasticity::setMode(mode); + if (eC != 3) return; // no parent integrand + + eKg = 0; // No geometric stiffness (assuming linear behaviour) + eM = eS = 0; // Inertia and external forces are calulated by parent integrand + if (eKm) eKm = 2; // Index for stiffness matrix in parent integrand } @@ -44,7 +68,7 @@ void FractureElasticity::initIntegration (size_t nGp, size_t) bool FractureElasticity::initElement (const std::vector& MNPC, LocalIntegral& elmInt) { - if (primsol.empty()) + if (mySol.empty()) { std::cerr <<" *** FractureElasticity::initElement:" <<" No primary solution vectors."<< std::endl; @@ -52,20 +76,22 @@ bool FractureElasticity::initElement (const std::vector& 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; @@ -210,11 +236,19 @@ bool FractureElasticity::evalStress (double lambda, double mu, double Gc, } +bool FractureElasticity::evalStress (double lambda, double mu, double Gc, + const SymmTensor& epsilon, + double* Phi, SymmTensor& sigma) const +{ + return this->evalStress(lambda,mu,Gc,epsilon,Phi,sigma,nullptr,true); +} + + double FractureElasticity::getStressDegradation (const Vector& N, - const Vector& eC) const + const Vectors& eV) const { // Evaluate the crack phase field function, c(X) - double c = eC.empty() ? 1.0 : N.dot(eC); + double c = eV[eC].empty() ? 1.0 : N.dot(eV[eC]); // Evaluate the stress degradation function, g(c), ignoring negative values return c > 0.0 ? (1.0-alpha)*c*c + alpha : alpha; } @@ -254,7 +288,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; @@ -328,32 +362,57 @@ bool FractureElasticity::evalInt (LocalIntegral& elmInt, } -bool FractureElasticity::evalBou (LocalIntegral& elmInt, - const FiniteElement& fe, const Vec3& X, - const Vec3& normal) const +bool FractureElasticity::checkSolVec (const Vectors& eV, size_t nen, + const char* name) const { - if (!tracFld && !fluxFld) - { - std::cerr <<" *** FractureElasticity::evalBou: No tractions."<< std::endl; - return false; - } - else if (!eS) - { - if (m_mode == SIM::RECOVERY) return true; - std::cerr <<" *** FractureElasticity::evalBou: No load vector."<< std::endl; - return false; - } + if (eV.size() <= eC) + std::cerr <<" *** FractureElasticity::"<< name + <<": Missing solution vector."<< std::endl; + else if (!eV.front().empty() && eV.front().size() != nsd*nen) + std::cerr <<" *** FractureElasticity::"<< name + <<": Invalid displacement vector.\n size(eV) = " + << eV.front().size() <<" size(dNdX) = "<< nsd*nen << std::endl; + else if (!eV[eC].empty() && eV[eC].size() != nen) + std::cerr <<" *** FractureElasticity::"<< name + <<": Invalid phase field vector.\n size(eC) = " + << eV[eC].size() <<" size(N) = "<< nen << std::endl; + else + return true; - // Evaluate the surface traction - Vec3 T = this->getTraction(X,normal); + return false; +} - // Integrate the force vector - Vector& ES = static_cast(elmInt).b[eS-1]; - for (size_t a = 1; a <= fe.N.size(); a++) - for (unsigned short int i = 1; i <= nsd; i++) - ES(nsd*(a-1)+i) += T[i-1]*fe.N(a)*fe.detJxW; - return true; +double FractureElasticity::crackStretch (const Vectors& eV, + const FiniteElement& fe, + const Vec3& X) const +{ + if (!this->checkSolVec(eV,fe.N.size(),"crackStretch")) + return -1.0; + else if (eV[eC].empty()) + return 1.0; // No phase field ==> unit stretch + + // Compute the deformation gradient (F) + Tensor F(nDF), C(nDF); + if (!this->formDefGradient(eV.front(),fe.N,fe.dNdX,X.x,F)) + return -2.0; + + // Compute the inverse right Cauchy-Green tensor (C^-1) + SymmTensor Ci(nDF); + if (Ci.rightCauchyGreen(F).inverse() == 0.0) + return -3.0; + + // Compute the phase field gradient (dD/dX where D = 1 - C) + Vector eD(eV[eC]), tmp(nsd); + for (size_t i = 0; i < eD.size(); i++) + eD[i] = 1.0 - eD[i]; // d = 1 - c + if (!fe.dNdX.multiply(eD,tmp,true)) + return -4.0; + + // Calculate the perpendicular crack stretch + // lambda = gradD*gradD / gradD*Ci*gradD (see eq. (106) in Miehe2015pfm3) + Vec3 gradD(tmp); + return gradD.length2() / (gradD*(Ci*gradD)); } @@ -362,14 +421,14 @@ bool FractureElasticity::evalSol (Vector& s, const std::vector& MNPC) const { // Extract element displacements - Vectors eV(2); + Vectors eV(1+eC); 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) - ierr = utl::gather(MNPC,1,myCVec,eV.back()); + ierr = utl::gather(MNPC,1,myCVec,eV[eC]); if (ierr > 0) { @@ -378,7 +437,6 @@ bool FractureElasticity::evalSol (Vector& s, return false; } - return this->evalSol2(s,eV,fe,X); } @@ -389,26 +447,8 @@ bool FractureElasticity::evalSol (Vector& s, const Vectors& eV, { PROFILE3("FractureEl::evalSol"); - if (eV.size() < 2) - { - std::cerr <<" *** FractureElasticity::evalSol: Missing solution vector." - << std::endl; - return false; - } - else if (!eV.front().empty() && eV.front().size() != fe.dNdX.rows()*nsd) - { - std::cerr <<" *** FractureElasticity::evalSol: Invalid displacement vector." - <<"\n size(eV) = "<< eV.front().size() <<" size(dNdX) = " - << fe.dNdX.rows() <<","<< fe.dNdX.cols() << std::endl; - return false; - } - else if (!eV[1].empty() && eV[1].size() != fe.N.size()) - { - std::cerr <<" *** FractureElasticity::evalSol: Invalid phase field vector." - <<"\n size(eC) = "<< eV[1].size() <<" size(N) = " - << fe.N.size() << std::endl; + if (!this->checkSolVec(eV,fe.N.size(),"evalSol")) return false; - } // Evaluate the symmetric strain tensor, eps Matrix Bmat; @@ -428,7 +468,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; @@ -471,14 +511,6 @@ bool FractureElasticity::evalSol (Vector& s, const Vectors& eV, } -bool FractureElasticity::evalStress (double lambda, double mu, double Gc, - const SymmTensor& epsilon, - double* Phi, SymmTensor& sigma) const -{ - return this->evalStress(lambda,mu,Gc,epsilon,Phi,sigma,nullptr,true); -} - - 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 7e7dede..1b75439 100644 --- a/FractureElasticity.h +++ b/FractureElasticity.h @@ -29,9 +29,23 @@ class FractureElasticity : public Elasticity //! \brief The constructor invokes the parent class constructor only. //! \param[in] n Number of spatial dimensions FractureElasticity(unsigned short int n); + //! \brief Constructor for integrands with a parent integrand. + //! \param parent The parent integrand of this one + //! \param[in] n Number of spatial dimensions + FractureElasticity(IntegrandBase* parent, unsigned short int n); //! \brief Empty destructor. virtual ~FractureElasticity() {} + //! \brief Sets the number of solution variables per node. + void setVar(unsigned short int n) { npv = n; } + + //! \brief 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; } + //! \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); @@ -48,14 +62,6 @@ class FractureElasticity : public Elasticity virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe, const Vec3& X) const; - //! \brief Evaluates the integrand at a boundary 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 - //! \param[in] normal Boundary normal vector at current integration point - virtual bool evalBou(LocalIntegral& elmInt, const FiniteElement& fe, - const Vec3& X, const Vec3& normal) const; - using Elasticity::evalSol; //! \brief Evaluates the secondary solution at a result point. //! \param[out] s Array of solution field values at current point @@ -63,7 +69,7 @@ class FractureElasticity : public Elasticity //! \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; + const Vec3& X, const std::vector& MNPC) const; //! \brief Evaluates the finite element (FE) solution at an integration point. //! \param[out] s The FE stress values at current point @@ -76,7 +82,15 @@ class FractureElasticity : public Elasticity const Vec3& X, bool toLocal, Vec3* pdir) const; //! \brief Returns a pointer to the Gauss-point tensile energy array. - const RealArray* getTensileEnergy() const { return &myPhi; } + virtual const RealArray* getTensileEnergy() const { return &myPhi; } + + //! \brief Evaluates the crack stretch at an integration point. + //! \param[in] eV Element solution vectors + //! \param[in] fe Finite element data at current point + //! \param[in] X Cartesian coordinates of current point + //! \return Perpendicular relative crack opening meassure + double crackStretch(const Vectors& eV, const FiniteElement& fe, + const Vec3& X) const; //! \brief Returns the number of primary/secondary solution field components. //! \param[in] fld which field set to consider (1=primary, 2=secondary) @@ -103,13 +117,20 @@ 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; + + //! \brief Internal helper used to check solution vector dimensions. + bool checkSolVec(const Vectors& eV, size_t nen, const char* name) 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/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..f8709f4 --- /dev/null +++ b/PoroFracture.C @@ -0,0 +1,174 @@ +// $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" +#include "Utilities.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); +} + + +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 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.expandRows(-1); // Assign velocity vector + eMat.resize(nsd+1,MNPC.size()); + eMat = elmInt.vec[uA]; // Interleaved acceleration vector + elmInt.vec[uA] = eMat.expandRows(-1); // Assign acceleration vector + // We don't need the pressure acceleration + + return true; +} + + +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 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::const_iterator uEnd = MNPC.begin() + elem_sizes.front(); + if (!fracEl->initElement(std::vector(MNPC.begin(),uEnd),elmInt)) + return false; + + if (primsol.size() < 2) + return true; // Quasi-static simulation + + // Extract the element level pressure rate vector + std::vector 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; +} + + +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 + if (!fracEl->evalInt(elMat,fe,X)) + return false; + + PoroElasticity::Mats* pelm = dynamic_cast(&elMat); + if (!pelm) return false; + + // Evaluate the perpendicular crack stretch + pelm->setCrackStretch(fracEl->crackStretch(elMat.vec,fe,X)); + return true; +} diff --git a/PoroFracture.h b/PoroFracture.h new file mode 100644 index 0000000..de25c14 --- /dev/null +++ b/PoroFracture.h @@ -0,0 +1,98 @@ +// $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 Elasticity::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); + + //! \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/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 5fac336..9a293c7 100644 --- a/SIMDynElasticity.h +++ b/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. //! //============================================================================== @@ -16,24 +15,25 @@ #define _SIM_DYN_ELASTICITY_H_ #include "NewmarkSIM.h" -#include "SIMElasticity.h" +#include "SIMElasticityWrap.h" #include "FractureElasticityVoigt.h" -#include "DataExporter.h" +#include "PoroFracture.h" /*! - \brief Driver wrapping an elasticity solver with an ISolver interface. + \brief Driver class for dynamic elasticity problems with fracture. */ -template -class SIMDynElasticity : public SIMElasticity +template> +class SIMDynElasticity : public Sim { public: //! \brief Default constructor. - SIMDynElasticity() : SIMElasticity(false), dSim(*this), vtfStep(0) - { - Dim::myHeading = "Elasticity solver"; - } + SIMDynElasticity() : dSim(*this), vtfStep(0) {} + + //! \brief Constructor for mixed problems. + SIMDynElasticity(const std::vector& flds) + : Sim(flds), dSim(*this), vtfStep(0) {} //! \brief Empty destructor. virtual ~SIMDynElasticity() {} @@ -45,46 +45,25 @@ 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(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; @@ -113,7 +92,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..."; @@ -162,9 +141,9 @@ 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(); + return static_cast(Dim::myProblem)->getTensileEnergy(); } //! \brief Returns a const reference to the global norms. @@ -174,7 +153,7 @@ class SIMDynElasticity : public SIMElasticity 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 +179,7 @@ class SIMDynElasticity : public SIMElasticity //! \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(Dim::myProblem); } @@ -212,15 +191,21 @@ 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) + { + std::string form("voigt"); + utl::getAttribute(elem,"formulation",form,true); + if (this->getName() == "PoroElasticity") + Dim::myProblem = new PoroFracture(Dim::dimension); + else if (form != "voigt") + Dim::myProblem = new FractureElasticity(Dim::dimension); + } + result = this->Sim::parse(elem); } else - result = this->SIMElasticity::parse(elem); + result = this->Dim::parse(elem); --ncall; return result; } diff --git a/main_FractureDynamics.C b/main_FractureDynamics.C index 24d68e6..e6af8f0 100644 --- a/main_FractureDynamics.C +++ b/main_FractureDynamics.C @@ -17,6 +17,7 @@ #include "SIMDynElasticity.h" #include "SIMPhaseField.h" #include "SIMFractureDynamics.h" +#include "SIMPoroElasticity.h" #include "SIMCoupledSI.h" #include "SIMSolver.h" #include "SIMSolverTS.h" @@ -72,13 +73,13 @@ private: \param[in] infile The input file to parse */ -template class Cpl, template class Solver=SIMSolver> int runSimulator2 (char* infile) { - typedef SIMDynElasticity SIMElastoDynamics; - typedef SIMPhaseField SIMCrackField; + typedef SIMDynElasticity SIMElastoDynamics; + typedef SIMPhaseField SIMCrackField; typedef SIMFracture SIMFractureDynamics; @@ -140,10 +141,10 @@ int runSimulator2 (char* infile) \param[in] context Input-file context for the time integrator */ -template -int runSimulator3 (char* infile, const char* context = "newmarksolver") +template +int runSimulator3 (char* infile, const char* context) { - typedef SIMDynElasticity SIMElastoDynamics; + typedef SIMDynElasticity SIMElastoDynamics; utl::profiler->start("Model input"); IFEM::cout <<"\n\n0. Parsing input file(s)." @@ -186,19 +187,56 @@ int runSimulator3 (char* infile, const char* context = "newmarksolver") } +/*! + \brief Creates and launches a stand-alone elasticity simulator (no coupling). + \param[in] infile The input file to parse + \param[in] poroel Use poroelastic solver + \param[in] context Input-file context for the time integrator +*/ + +template +int runSimulator4 (char* infile, bool poroel, + const char* context = "newmarksolver") +{ + if (poroel) + return runSimulator3>(infile,context); + + return runSimulator3>(infile,context); +} + + /*! \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 Cpl> +template class Cpl> int runSolver (char* infile, bool timeslabs) { if (timeslabs) - return runSimulator2(infile); + return runSimulator2(infile); + + return runSimulator2(infile); +} + + +/*! + \brief Creates the combined fracture simulator and launches the simulation. + \param[in] infile The input file to parse + \param[in] poroel Use poroelastic solver + \param[in] adaptiv Use time-slab adaptive solver +*/ + +template class Cpl> +int runSolver (char* infile, bool poroel, bool adaptiv) +{ + if (poroel) + return runSolver,Cpl>(infile,adaptiv); - return runSimulator2(infile); + return runSolver,Cpl>(infile,adaptiv); } @@ -206,18 +244,19 @@ int runSolver (char* infile, bool timeslabs) \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] poroel Use poroelastic solver \param[in] timeslabs Use time-slab adaptive solver */ template -int runSimulator1 (char* infile, char coupling, bool timeslabs) +int runSimulator1 (char* infile, char coupling, bool poroel, bool timeslabs) { if (coupling == 1) - return runSolver(infile,timeslabs); + return runSolver(infile,poroel,timeslabs); else if (coupling == 2) - return runSolver(infile,timeslabs); + return runSolver(infile,poroel,timeslabs); else - return runSimulator3(infile); // No phase field coupling + return runSimulator4(infile,poroel); // No phase field } @@ -241,18 +280,20 @@ public: \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] poroel Use poroelastic solver \param[in] timeslabs Use time-slab adaptive solver */ template -int runSimulator (char* infile, char integrator, char coupling, bool timeslabs) +int runSimulator (char* infile, char integrator, char coupling, + bool poroel, bool timeslabs) { if (integrator == 2) - return runSimulator1(infile,coupling,timeslabs); + return runSimulator1(infile,coupling,poroel,timeslabs); else if (integrator > 0) - return runSimulator1(infile,coupling,timeslabs); + return runSimulator1(infile,coupling,poroel,timeslabs); else - return runSimulator3(infile,"staticsolver"); + return runSimulator4(infile,poroel,"staticsolver"); } @@ -269,7 +310,9 @@ int main (int argc, char** argv) char coupling = 1; char integrator = 1; bool twoD = false; + bool poroel = false; bool adaptive = false; + ASMmxBase::Type = ASMmxBase::NONE; IFEM::Init(argc,argv); @@ -278,6 +321,8 @@ int main (int argc, char** argv) ; // ignore the obsolete option else if (!strcmp(argv[i],"-2D")) twoD = SIMElasticity::planeStrain = true; + else if (!strcmp(argv[i],"-mixed")) + ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1; else if (!strcmp(argv[i],"-nocrack")) coupling = 0; else if (!strcmp(argv[i],"-semiimplicit")) @@ -290,6 +335,8 @@ int main (int argc, char** argv) Elasticity::wantPrincipalStress = true; else if (!strcmp(argv[i],"-dbgElm") && i < argc-1) FractureElasticNorm::dbgElm = atoi(argv[++i]); + else if (!strncmp(argv[i],"-poro",5)) + poroel = true; else if (!strncmp(argv[i],"-adap",5)) adaptive = true; else if (!infile) @@ -301,8 +348,8 @@ 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] [-static|-GA] [-adaptive]\n" + <<" [-lag|-spec|-LR] [-2D] [-mixed] [-nGauss ]\n " + <<" [-nocrack|-semiimplicit] [-static|-GA] [-poro] [-adaptive]\n" <<" [-vtf [-nviz ] [-nu ] [-nv ]] [-hdf5] [-principal]\n"<< std::endl; return 0; @@ -320,7 +367,7 @@ int main (int argc, char** argv) IFEM::cout << std::endl; if (twoD) - return runSimulator(infile,integrator,coupling,adaptive); + return runSimulator(infile,integrator,coupling,poroel,adaptive); else - return runSimulator(infile,integrator,coupling,adaptive); + return runSimulator(infile,integrator,coupling,poroel,adaptive); }