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 010bf84..cc9aaa2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -21,30 +21,53 @@ 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(EXTRA_DOXY_PATHS "${EXTRA_DOXY_PATHS} \\ - ${PROJECT_SOURCE_DIR}/../Elasticity") - if(NOT TARGET Elasticity) - add_subdirectory(../Elasticity 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} ${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}) @@ -66,6 +89,8 @@ else() ${PROJECT_SOURCE_DIR}/Test/*_CH.reg) file(GLOB FD_TESTFILES RELATIVE ${PROJECT_SOURCE_DIR}/Test ${PROJECT_SOURCE_DIR}/Test/*_FD.reg) + file(GLOB FP_TESTFILES RELATIVE ${PROJECT_SOURCE_DIR}/Test + ${PROJECT_SOURCE_DIR}/Test/*_FDporo.reg) foreach(TESTFILE ${CH_TESTFILES}) ifem_add_test(${TESTFILE} CahnHilliard) endforeach() @@ -81,6 +106,11 @@ else() foreach(TESTFILE ${FD_TESTFILES}) ifem_add_test(${TESTFILE} FractureDynamics) endforeach() + if(EXISTS ${POROELASTIC_DIR}) + foreach(TESTFILE ${FP_TESTFILES}) + ifem_add_test(${TESTFILE} FractureDynamics) + endforeach() + endif() endif() list(APPEND TEST_APPS CahnHilliard FractureDynamics) @@ -88,7 +118,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..d1d5606 100644 --- a/FractureArgs.C +++ b/FractureArgs.C @@ -20,6 +20,8 @@ FractureArgs::FractureArgs () : SIMargsBase("fracturedynamics") { inpfile = nullptr; integrator = coupling = 1; + poroEl = expPhase = false; + stopT = -1.0; } @@ -41,6 +43,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 +84,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..5e25afb 100644 --- a/FractureArgs.h +++ b/FractureArgs.h @@ -30,6 +30,9 @@ 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 + double stopT; //!< Stop time of the simulation //! \brief Default constructor. FractureArgs(); @@ -40,7 +43,7 @@ class FractureArgs : public SIMargsBase void parseFile(const char* argv, int& iarg); protected: - //! \brief Parse an element from the input file + //! \brief Parses an element from the input file. virtual bool parse(const TiXmlElement* elem); }; diff --git a/FractureElasticity.C b/FractureElasticity.C index 3573515..ad80d8e 100644 --- a/FractureElasticity.C +++ b/FractureElasticity.C @@ -121,6 +121,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 +521,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 +534,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 +640,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 efdb46a..dbcecec 100644 --- a/FractureElasticity.h +++ b/FractureElasticity.h @@ -43,8 +43,9 @@ class FractureElasticity : public Elasticity //! \brief Prints out the problem definition to the log stream. virtual void printLog() const; - //! \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 +97,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..b6df31b --- /dev/null +++ b/PoroFracture.C @@ -0,0 +1,328 @@ +// $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] nd Number of spatial dimensions + //! \param[in] nv Number of primary solution variables per node + PoroFractureElasticity(IntegrandBase* parent, + unsigned short int nd, unsigned short int nv) + : FractureElasticityVoigt(parent,nd) { npv = nv; } + //! \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, bool m) : PoroElasticity(n,m) +{ + fracEl = new PoroFractureElasticity(this, n, m ? n : n+1); + + 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); +} + + +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..3d4b20a --- /dev/null +++ b/PoroFracture.h @@ -0,0 +1,132 @@ +// $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 + //! \param[in] m If \e true, a mixed formulation is used + explicit PoroFracture(unsigned short int n, bool m = false); + //! \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); + + 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 elMat 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 616d1d5..282774d 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,20 @@ \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. + explicit SIMDynElasticity(const std::vector& nf) + : Sim(nf), dSim(*this) + { vtfStep = 0; } @@ -40,28 +48,18 @@ class SIMDynElasticity : public SIMElasticity virtual ~SIMDynElasticity() {} //! \brief Prints out problem-specific data to the log stream. - virtual void printProblem() const + void printProblem() const override { static short int ncall = 0; 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&) + bool init(const TimeStep&, bool = false) override { dSim.initPrm(); dSim.initSol(dynamic_cast(&dSim) ? 3 : 1); @@ -72,21 +70,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) + bool saveStep(const TimeStep& tp, int& nBlock) override { double old = utl::zero_print_tol; utl::zero_print_tol = 1e-16; @@ -140,12 +127,12 @@ class SIMDynElasticity : public SIMElasticity } //! \brief Advances the time step one step forward. - virtual bool advanceStep(TimeStep& tp) { return dSim.advanceStep(tp,false); } + bool advanceStep(TimeStep& tp) override { return dSim.advanceStep(tp,false); } //! \brief Prints out time step identification. //! \param[in] istep Time step counter //! \param[in] time Parameters for time-dependent simulations - virtual void printStep(int istep, const TimeDomain& time) const + void printStep(int istep, const TimeDomain& time) const override { Dim::adm.cout <<"\n step="<< istep <<" time="<< time.t; RealFunc* p = this->haveCrackPressure(); @@ -156,7 +143,7 @@ class SIMDynElasticity : public SIMElasticity } //! \brief Computes the solution for the current time step. - bool solveStep(TimeStep& tp) + bool solveStep(TimeStep& tp) override { if (Dim::msgLevel >= 1) IFEM::cout <<"\n Solving the elasto-dynamics problem..."; @@ -224,7 +211,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(); } @@ -245,8 +232,10 @@ class SIMDynElasticity : public SIMElasticity } } + //! \brief Returns a const reference to current solution vector. + const Vector& getSolution(int i) const override { return dSim.getSolution(i);} //! \brief Returns a const reference to the solution vectors. - const Vectors& getSolutions() const { return dSim.getSolutions(); } + const Vectors& getSolutions() const override { return dSim.getSolutions(); } //! \brief Updates the solution vectors. void setSolutions(const Vectors& dvec) @@ -300,27 +289,41 @@ class SIMDynElasticity : public SIMElasticity protected: //! \brief Returns the actual integrand. - virtual Elasticity* getIntegrand() + Elasticity* getIntegrand() override { if (!Dim::myProblem) // Using the Voigt formulation by default Dim::myProblem = new FractureElasticityVoigt(Dim::dimension); return static_cast(Dim::myProblem); } - using SIMElasticity::parse; + using SIMElasticityWrap::parse; //! \brief Parses a data section from an XML element. - virtual bool parse(const TiXmlElement* elem) + bool parse(const TiXmlElement* elem) override { bool result = true; 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, + this->mixedProblem()); +#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()) @@ -328,7 +331,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..dd1410d --- /dev/null +++ b/SIMExplPhaseField.h @@ -0,0 +1,108 @@ +// $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. + explicit 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. + bool checkStopCriterion () const { return false; } + //! \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/Test/Miehe71-explcrack_FD.reg b/Test/Miehe71-explcrack_FD.reg new file mode 100644 index 0000000..bcde142 --- /dev/null +++ b/Test/Miehe71-explcrack_FD.reg @@ -0,0 +1,91 @@ +Miehe71.xinp -qstatic -explcrack + +Input file: Miehe71.xinp +Equation solver: 2 +Number of Gauss points: 4 +0. Parsing input file(s). +1. Elasticity solver +Parsing input file Miehe71.xinp +Parsing +Parsing + Generating linear geometry on unit parameter domain \[0,1]^2 + Length in X = 5 + Length in Y = 5 + Parsing + Parsing + Topology sets: Left (1,1,1D) + Right (1,2,1D) + TopBottom (1,3,1D) (1,4,1D) + Parsing + Refining P1 19 19 + Parsing +Parsing +Parsing +Parsing +Parsing + Parsing + Material code 0: 254.8 0.3 700 + Parsing + Dirichlet code 2: (fixed) + Parsing + Dirichlet code 1: (fixed) + Parsing + Dirichlet code 1000001 (ramp): 3 \* Ramp(t,10) +Parsing +Parsing +Parsing input file succeeded. +Equation solver: 2 +Number of Gauss points: 2 +2. Explicit phase field +Parsing input file Miehe71.xinp +Parsing +Parsing +Parsing + Phase-field function: 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) +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing input file succeeded. +Equation solver: 2 +Number of Gauss points: 4 +3. Time integration driver +Parsing input file Miehe71.xinp +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing input file succeeded. +10. Preprocessing the finite element model: +11. Elasticity solver +Problem definition: +Elasticity: 2D, gravity = 0 0 +LinIsotropic: E = 254.8, nu = 0.3, rho = 700, alpha = 1.2e-07 + Degrading of tensile strain energy density. +Resolving Dirichlet boundary conditions + Constraining P1 E3 in direction(s) 2 + Constraining P1 E4 in direction(s) 2 + Constraining P1 E1 in direction(s) 1 + Constraining P1 E2 in direction(s) 1 code = 1000001 + >>> SAM model summary <<< +Number of elements 400 +Number of nodes 441 +Number of dofs 882 +Number of constraints 21 +Number of unknowns 798 +100. Starting the simulation + Solving the elasto-dynamics problem... + step=1 time=1 + Primary solution summary: L2-norm : 0.123996 + Max X-component : 0.3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -30.87 + Elastic strain energy: eps_e : 15.435 + Bulk energy: eps_b : 15.435 + Tensile & compressive energies : 15.435 + External energy: ((f,u^h)+(t,u^h))^0.5 : -3.92874 diff --git a/Test/Miehe71-explcrack_FDporo.reg b/Test/Miehe71-explcrack_FDporo.reg new file mode 100644 index 0000000..00e9820 --- /dev/null +++ b/Test/Miehe71-explcrack_FDporo.reg @@ -0,0 +1,98 @@ +Miehe71.xinp -qstatic -explcrack -poro + +Input file: Miehe71.xinp +Equation solver: 2 +Number of Gauss points: 4 +0. Parsing input file(s). +1. Poroelasticity solver +Parsing input file Miehe71.xinp +Parsing +Parsing + Generating linear geometry on unit parameter domain \[0,1]^2 + Length in X = 5 + Length in Y = 5 + Parsing + Parsing + Topology sets: Left (1,1,1D) + Right (1,2,1D) + TopBottom (1,3,1D) (1,4,1D) + Parsing + Refining P1 19 19 + Parsing +Parsing +Parsing +Parsing + Parsing + Material code 0: Poroelastic material, see "Problem definition:" below. + Parsing + Dirichlet code 2: (fixed) + Parsing + Dirichlet code 1: (fixed) + Parsing + Dirichlet code 1000001 (ramp): 3 \* Ramp(t,10) +Parsing +Parsing +Parsing +Parsing input file succeeded. +Equation solver: 2 +Number of Gauss points: 2 +2. Explicit phase field +Parsing input file Miehe71.xinp +Parsing +Parsing +Parsing + Phase-field function: 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) +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing input file succeeded. +Equation solver: 2 +Number of Gauss points: 4 +3. Time integration driver +Parsing input file Miehe71.xinp +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing input file succeeded. +10. Preprocessing the finite element model: +11. Poroelasticity solver +Problem definition: +PoroElasticity: scaling = 0 useDynCoupling = false +Elasticity: 2D, gravity = 0 9.81 + Constitutive Properties: + Young's Modulus, E = 254.8 + Poisson's Ratio, nu = 0.3 + Densities: + Density of Fluid, rhof = 1000 + Density of Solid, rhos = 666.67 + Bulk Moduli: + Biot's coefficient, alpha = 1 + Biot's inverse modulus, M^-1 = 0.01 + Porosity, n = 0.1 + Permeability, K = 2.07e-09 2.07e-09 0 +Resolving Dirichlet boundary conditions + Constraining P1 E3 in direction(s) 2 + Constraining P1 E4 in direction(s) 2 + Constraining P1 E1 in direction(s) 1 + Constraining P1 E2 in direction(s) 1 code = 1000001 + >>> SAM model summary <<< +Number of elements 400 +Number of nodes 441 +Number of dofs 1323 +Number of constraints 21 +Number of unknowns 1239 +100. Starting the simulation + Solving the elasto-dynamics problem... + step=1 time=1 + Primary solution summary: L2-norm : 0.101242 + Max X-displacement : 0.3 + Max pressure : 1.72664e-07 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -30.87 diff --git a/Test/Miehe71-nocrack_FD.reg b/Test/Miehe71-nocrack_FD.reg new file mode 100644 index 0000000..5678031 --- /dev/null +++ b/Test/Miehe71-nocrack_FD.reg @@ -0,0 +1,209 @@ +Miehe71.xinp -qstatic -nocrack -stopTime 50 + +Input file: Miehe71.xinp +Equation solver: 2 +Number of Gauss points: 4 +Simulation stop time: 50 +0. Parsing input file(s). +1. Elasticity solver +Parsing input file Miehe71.xinp +Parsing +Parsing + Generating linear geometry on unit parameter domain \[0,1]^2 + Length in X = 5 + Length in Y = 5 + Parsing + Parsing + Topology sets: Left (1,1,1D) + Right (1,2,1D) + TopBottom (1,3,1D) (1,4,1D) + Parsing + Refining P1 19 19 + Parsing +Parsing +Parsing +Parsing +Parsing + Parsing + Material code 0: 254.8 0.3 700 + Parsing + Dirichlet code 2: (fixed) + Parsing + Dirichlet code 1: (fixed) + Parsing + Dirichlet code 1000001 (ramp): 3 \* Ramp(t,10) +Parsing +Parsing +Parsing input file succeeded. +Equation solver: 2 +Number of Gauss points: 2 +2. Time integration driver +Parsing input file Miehe71.xinp +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing input file succeeded. +10. Preprocessing the finite element model: +11. Elasticity solver +Problem definition: +Elasticity: 2D, gravity = 0 0 +LinIsotropic: E = 254.8, nu = 0.3, rho = 700, alpha = 1.2e-07 + Degrading of tensile strain energy density. +Resolving Dirichlet boundary conditions + Constraining P1 E3 in direction(s) 2 + Constraining P1 E4 in direction(s) 2 + Constraining P1 E1 in direction(s) 1 + Constraining P1 E2 in direction(s) 1 code = 1000001 + >>> SAM model summary <<< +Number of elements 400 +Number of nodes 441 +Number of dofs 882 +Number of constraints 21 +Number of unknowns 798 +100. Starting the simulation + Solving the elasto-dynamics problem... + step=1 time=1 + Primary solution summary: L2-norm : 0.123996 + Max X-component : 0.3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -30.87 + Elastic strain energy: eps_e : 15.435 + Bulk energy: eps_b : 15.435 + Tensile & compressive energies : 15.435 + External energy: ((f,u^h)+(t,u^h))^0.5 : -3.92874 + Solving the elasto-dynamics problem... + step=2 time=2 + Primary solution summary: L2-norm : 0.247992 + Max X-component : 0.6 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -123.48 + Elastic strain energy: eps_e : 61.74 + Bulk energy: eps_b : 61.74 + Tensile & compressive energies : 61.74 + External energy: ((f,u^h)+(t,u^h))^0.5 : -7.85748 + Solving the elasto-dynamics problem... + step=3 time=3 + Primary solution summary: L2-norm : 0.371988 + Max X-component : 0.9 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -277.83 + Elastic strain energy: eps_e : 138.915 + Bulk energy: eps_b : 138.915 + Tensile & compressive energies : 138.915 + External energy: ((f,u^h)+(t,u^h))^0.5 : -11.7862 + Solving the elasto-dynamics problem... + step=4 time=4 + Primary solution summary: L2-norm : 0.495984 + Max X-component : 1.2 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -493.92 + Elastic strain energy: eps_e : 246.96 + Bulk energy: eps_b : 246.96 + Tensile & compressive energies : 246.96 + External energy: ((f,u^h)+(t,u^h))^0.5 : -15.715 + Solving the elasto-dynamics problem... + step=5 time=5 + Primary solution summary: L2-norm : 0.61998 + Max X-component : 1.5 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -771.75 + Elastic strain energy: eps_e : 385.875 + Bulk energy: eps_b : 385.875 + Tensile & compressive energies : 385.875 + External energy: ((f,u^h)+(t,u^h))^0.5 : -19.6437 + Solving the elasto-dynamics problem... + step=6 time=6 + Primary solution summary: L2-norm : 0.743976 + Max X-component : 1.8 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -1111.32 + Elastic strain energy: eps_e : 555.66 + Bulk energy: eps_b : 555.66 + Tensile & compressive energies : 555.66 + External energy: ((f,u^h)+(t,u^h))^0.5 : -23.5724 + Solving the elasto-dynamics problem... + step=7 time=7 + Primary solution summary: L2-norm : 0.867972 + Max X-component : 2.1 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -1512.63 + Elastic strain energy: eps_e : 756.315 + Bulk energy: eps_b : 756.315 + Tensile & compressive energies : 756.315 + External energy: ((f,u^h)+(t,u^h))^0.5 : -27.5012 + Solving the elasto-dynamics problem... + step=8 time=8 + Primary solution summary: L2-norm : 0.991968 + Max X-component : 2.4 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -1975.68 + Elastic strain energy: eps_e : 987.84 + Bulk energy: eps_b : 987.84 + Tensile & compressive energies : 987.84 + External energy: ((f,u^h)+(t,u^h))^0.5 : -31.4299 + Solving the elasto-dynamics problem... + step=9 time=9 + Primary solution summary: L2-norm : 1.11596 + Max X-component : 2.7 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -2500.47 + Elastic strain energy: eps_e : 1250.23 + Bulk energy: eps_b : 1250.23 + Tensile & compressive energies : 1250.23 + External energy: ((f,u^h)+(t,u^h))^0.5 : -35.3587 + Solving the elasto-dynamics problem... + step=10 time=10 + Primary solution summary: L2-norm : 1.23996 + Max X-component : 3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -3087 + Elastic strain energy: eps_e : 1543.5 + Bulk energy: eps_b : 1543.5 + Tensile & compressive energies : 1543.5 + External energy: ((f,u^h)+(t,u^h))^0.5 : -39.2874 + Solving the elasto-dynamics problem... + step=11 time=20 + Primary solution summary: L2-norm : 1.23996 + Max X-component : 3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -3087 + Elastic strain energy: eps_e : 1543.5 + Bulk energy: eps_b : 1543.5 + Tensile & compressive energies : 1543.5 + External energy: ((f,u^h)+(t,u^h))^0.5 : -39.2874 + Solving the elasto-dynamics problem... + step=12 time=30 + Primary solution summary: L2-norm : 1.23996 + Max X-component : 3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -3087 + Elastic strain energy: eps_e : 1543.5 + Bulk energy: eps_b : 1543.5 + Tensile & compressive energies : 1543.5 + External energy: ((f,u^h)+(t,u^h))^0.5 : -39.2874 + Solving the elasto-dynamics problem... + step=13 time=40 + Primary solution summary: L2-norm : 1.23996 + Max X-component : 3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -3087 + Elastic strain energy: eps_e : 1543.5 + Bulk energy: eps_b : 1543.5 + Tensile & compressive energies : 1543.5 + External energy: ((f,u^h)+(t,u^h))^0.5 : -39.2874 + Solving the elasto-dynamics problem... + step=14 time=50 + Primary solution summary: L2-norm : 1.23996 + Max X-component : 3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -3087 + Elastic strain energy: eps_e : 1543.5 + Bulk energy: eps_b : 1543.5 + Tensile & compressive energies : 1543.5 + External energy: ((f,u^h)+(t,u^h))^0.5 : -39.2874 + Time integration completed. diff --git a/Test/Miehe71-nocrack_FDporo.reg b/Test/Miehe71-nocrack_FDporo.reg new file mode 100644 index 0000000..59a75f8 --- /dev/null +++ b/Test/Miehe71-nocrack_FDporo.reg @@ -0,0 +1,147 @@ +Miehe71.xinp -qstatic -nocrack -poro + +Input file: Miehe71.xinp +Equation solver: 2 +Number of Gauss points: 4 +0. Parsing input file(s). +1. Poroelasticity solver +Parsing input file Miehe71.xinp +Parsing +Parsing + Generating linear geometry on unit parameter domain \[0,1]^2 + Length in X = 5 + Length in Y = 5 + Parsing + Parsing + Topology sets: Left (1,1,1D) + Right (1,2,1D) + TopBottom (1,3,1D) (1,4,1D) + Parsing + Refining P1 19 19 + Parsing +Parsing +Parsing +Parsing + Parsing + Material code 0: Poroelastic material, see "Problem definition:" below. + Parsing + Dirichlet code 2: (fixed) + Parsing + Dirichlet code 1: (fixed) + Parsing + Dirichlet code 1000001 (ramp): 3 \* Ramp(t,10) +Parsing +Parsing +Parsing +Parsing input file succeeded. +Equation solver: 2 +Number of Gauss points: 2 +2. Time integration driver +Parsing input file Miehe71.xinp +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing input file succeeded. +10. Preprocessing the finite element model: +11. Poroelasticity solver +Problem definition: +PoroElasticity: scaling = 0 useDynCoupling = false +Elasticity: 2D, gravity = 0 9.81 + Constitutive Properties: + Young's Modulus, E = 254.8 + Poisson's Ratio, nu = 0.3 + Densities: + Density of Fluid, rhof = 1000 + Density of Solid, rhos = 666.67 + Bulk Moduli: + Biot's coefficient, alpha = 1 + Biot's inverse modulus, M^-1 = 0.01 + Porosity, n = 0.1 + Permeability, K = 2.07e-09 2.07e-09 0 +Resolving Dirichlet boundary conditions + Constraining P1 E3 in direction(s) 2 + Constraining P1 E4 in direction(s) 2 + Constraining P1 E1 in direction(s) 1 + Constraining P1 E2 in direction(s) 1 code = 1000001 + >>> SAM model summary <<< +Number of elements 400 +Number of nodes 441 +Number of dofs 1323 +Number of constraints 21 +Number of unknowns 1239 +100. Starting the simulation + Solving the elasto-dynamics problem... + step=1 time=1 + Primary solution summary: L2-norm : 0.101242 + Max X-displacement : 0.3 + Max pressure : 1.7266.e-07 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -30.87 + Solving the elasto-dynamics problem... + step=2 time=2 + Primary solution summary: L2-norm : 0.202485 + Max X-displacement : 0.6 + Max pressure : 3.4532.e-07 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -123.48 + Solving the elasto-dynamics problem... + step=3 time=3 + Primary solution summary: L2-norm : 0.303727 + Max X-displacement : 0.9 + Max pressure : 5.1799.e-07 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -277.83 + Solving the elasto-dynamics problem... + step=4 time=4 + Primary solution summary: L2-norm : 0.404969 + Max X-displacement : 1.2 + Max pressure : 6.9065.e-07 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -493.92 + Solving the elasto-dynamics problem... + step=5 time=5 + Primary solution summary: L2-norm : 0.506211 + Max X-displacement : 1.5 + Max pressure : 8.6332.e-07 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -771.75 + Solving the elasto-dynamics problem... + step=6 time=6 + Primary solution summary: L2-norm : 0.607454 + Max X-displacement : 1.8 + Max pressure : 1.03598e-06 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -1111.32 + Solving the elasto-dynamics problem... + step=7 time=7 + Primary solution summary: L2-norm : 0.708696 + Max X-displacement : 2.1 + Max pressure : 1.20865e-06 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -1512.63 + Solving the elasto-dynamics problem... + step=8 time=8 + Primary solution summary: L2-norm : 0.809938 + Max X-displacement : 2.4 + Max pressure : 1.38131e-06 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -1975.68 + Solving the elasto-dynamics problem... + step=9 time=9 + Primary solution summary: L2-norm : 0.911181 + Max X-displacement : 2.7 + Max pressure : 1.55398e-06 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -2500.47 + Solving the elasto-dynamics problem... + step=10 time=10 + Primary solution summary: L2-norm : 1.01242 + Max X-displacement : 3 + Max pressure : 1.72664e-06 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -3087 diff --git a/Test/Miehe71.xinp b/Test/Miehe71.xinp new file mode 100644 index 0000000..52cc9b8 --- /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 + 121 + + + diff --git a/Test/Miehe71_FD.reg b/Test/Miehe71_FD.reg new file mode 100644 index 0000000..1b26739 --- /dev/null +++ b/Test/Miehe71_FD.reg @@ -0,0 +1,357 @@ +Miehe71.xinp -qstatic -stopTime 50 + +Input file: Miehe71.xinp +Equation solver: 2 +Number of Gauss points: 4 +Simulation stop time: 50 +0. Parsing input file(s). +1. Elasticity solver +Parsing input file Miehe71.xinp +Parsing +Parsing + Generating linear geometry on unit parameter domain \[0,1]^2 + Length in X = 5 + Length in Y = 5 + Parsing + Parsing + Topology sets: Left (1,1,1D) + Right (1,2,1D) + TopBottom (1,3,1D) (1,4,1D) + Parsing + Refining P1 19 19 + Parsing +Parsing +Parsing +Parsing +Parsing + Parsing + Material code 0: 254.8 0.3 700 + Parsing + Dirichlet code 2: (fixed) + Parsing + Dirichlet code 1: (fixed) + Parsing + Dirichlet code 1000001 (ramp): 3 \* Ramp(t,10) +Parsing +Parsing +Parsing input file succeeded. +Equation solver: 2 +Number of Gauss points: 2 +2. Cahn-Hilliard solver +Parsing input file Miehe71.xinp +Parsing +Parsing + Parsing + Parsing + Topology sets: Left (1,1,1D) + Right (1,2,1D) + TopBottom (1,3,1D) (1,4,1D) + Parsing + Refining P1 19 19 + Parsing +Parsing + Initial crack function: abs(x-2.5) +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing input file succeeded. +Equation solver: 2 +Number of Gauss points: 2 +3. Time integration driver +Parsing input file Miehe71.xinp +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing input file succeeded. +10. Preprocessing the finite element model: +11. Elasticity solver +Problem definition: +Elasticity: 2D, gravity = 0 0 +LinIsotropic: E = 254.8, nu = 0.3, rho = 700, alpha = 1.2e-07 + Degrading of tensile strain energy density. +Resolving Dirichlet boundary conditions + Constraining P1 E3 in direction(s) 2 + Constraining P1 E4 in direction(s) 2 + Constraining P1 E1 in direction(s) 1 + Constraining P1 E2 in direction(s) 1 code = 1000001 + >>> SAM model summary <<< +Number of elements 400 +Number of nodes 441 +Number of dofs 882 +Number of constraints 21 +Number of unknowns 798 +12. Cahn-Hilliard solver +Problem definition: +Cahn-Hilliard: 2D + Critical fracture energy density: 1 + Smearing factor: 0.1 + Max value in crack: 0.001 + Initial crack specified as a function. + Enforcing crack irreversibility using history buffer. +Resolving Dirichlet boundary conditions + >>> SAM model summary <<< +Number of elements 400 +Number of nodes 441 +Number of dofs 441 +Number of unknowns 441 +100. Starting the simulation + Solving the elasto-dynamics problem... + step=1 time=1 + Primary solution summary: L2-norm : 0.123996 + Max X-component : 0.3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -30.87 + Elastic strain energy: eps_e : 15.435 + Bulk energy: eps_b : 15.435 + Tensile & compressive energies : 15.435 + External energy: ((f,u^h)+(t,u^h))^0.5 : -3.92874 + Solving crack phase field at step=1 time=1 + Primary solution summary: L2-norm : 0.749993 + Max value : 0.80195 + Min value : -0.125837 + Range : 1.12584 + Dissipated energy: eps_d : 9.21268 + L1-norm: |c^h| = (|c^h|) : 17.9004 + Normalized L1-norm: |c^h|/V : 0.716017 + Solving the elasto-dynamics problem... + step=2 time=2 + Primary solution summary: L2-norm : 0.266817 + Max X-component : 0.6 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -40.3582 + Elastic strain energy: eps_e : 20.1512 + Bulk energy: eps_b : 20.1512 + Tensile & compressive energies : 172.847 + External energy: ((f,u^h)+(t,u^h))^0.5 : -4.49211 + Solving crack phase field at step=2 time=2 + Primary solution summary: L2-norm : 0.712248 + Max value : 0.79706 + Min value : -0.028883 + Range : 1.02888 + Dissipated energy: eps_d : 12.1754 + L1-norm: |c^h| = (|c^h|) : 16.4857 + Normalized L1-norm: |c^h|/V : 0.659428 + Solving the elasto-dynamics problem... + step=3 time=3 + Primary solution summary: L2-norm : 0.438842 + Max X-component : 0.9 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -9.64493 + Elastic strain energy: eps_e : 4.78831 + Bulk energy: eps_b : 4.78831 + Tensile & compressive energies : 1201.4 + External energy: ((f,u^h)+(t,u^h))^0.5 : -2.19601 + Solving crack phase field at step=3 time=3 + Primary solution summary: L2-norm : 0.709708 + Max value : 0.79706 + Min value : -0.003712 + Range : 1.00371 + Dissipated energy: eps_d : 13.4453 + L1-norm: |c^h| = (|c^h|) : 16.2243 + Normalized L1-norm: |c^h|/V : 0.648973 + Solving the elasto-dynamics problem... + step=4 time=4 + Primary solution summary: L2-norm : 0.592569 + Max X-component : 1.2 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -0.526115 + Elastic strain energy: eps_e : 0.261561 + Bulk energy: eps_b : 0.261561 + Tensile & compressive energies : 2457.94 + External energy: ((f,u^h)+(t,u^h))^0.5 : -0.512891 + Solving crack phase field at step=4 time=4 + Primary solution summary: L2-norm : 0.709526 + Max value : 0.79706 + Min value : -0.001158 + Range : 1.00116 + Dissipated energy: eps_d : 13.5885 + L1-norm: |c^h| = (|c^h|) : 16.1978 + Normalized L1-norm: |c^h|/V : 0.647912 + Solving the elasto-dynamics problem... + step=5 time=5 + Primary solution summary: L2-norm : 0.740939 + Max X-component : 1.5 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -0.308136 + Elastic strain energy: eps_e : 0.103966 + Bulk energy: eps_b : 0.103966 + Tensile & compressive energies : 3854.06 + External energy: ((f,u^h)+(t,u^h))^0.5 : -0.392515 + Solving crack phase field at step=5 time=5 + Primary solution summary: L2-norm : 0.709465 + Max value : 0.79706 + Min value : -0.000588 + Range : 1.00059 + Dissipated energy: eps_d : 13.6409 + L1-norm: |c^h| = (|c^h|) : 16.1882 + Normalized L1-norm: |c^h|/V : 0.647527 + Solving the elasto-dynamics problem... + step=6 time=6 + Primary solution summary: L2-norm : 0.889182 + Max X-component : 1.8 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -0.18709 + Elastic strain energy: eps_e : 0.0623575 + Bulk energy: eps_b : 0.0623575 + Tensile & compressive energies : 5553.78 + External energy: ((f,u^h)+(t,u^h))^0.5 : -0.305851 + Solving crack phase field at step=6 time=6 + Primary solution summary: L2-norm : 0.709434 + Max value : 0.79706 + Min value : -0.000421 + Range : 1.00042 + Dissipated energy: eps_d : 13.6699 + L1-norm: |c^h| = (|c^h|) : 16.1829 + Normalized L1-norm: |c^h|/V : 0.647315 + Solving the elasto-dynamics problem... + step=7 time=7 + Primary solution summary: L2-norm : 1.0374 + Max X-component : 2.1 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -0.124697 + Elastic strain energy: eps_e : 0.0411738 + Bulk energy: eps_b : 0.0411738 + Tensile & compressive energies : 7561.28 + External energy: ((f,u^h)+(t,u^h))^0.5 : -0.249697 + Solving crack phase field at step=7 time=7 + Primary solution summary: L2-norm : 0.709415 + Max value : 0.79706 + Min value : -0.000315 + Range : 1.00032 + Dissipated energy: eps_d : 13.6875 + L1-norm: |c^h| = (|c^h|) : 16.1797 + Normalized L1-norm: |c^h|/V : 0.647186 + Solving the elasto-dynamics problem... + step=8 time=8 + Primary solution summary: L2-norm : 1.18561 + Max X-component : 2.4 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -0.0888788 + Elastic strain energy: eps_e : 0.0291413 + Bulk energy: eps_b : 0.0291413 + Tensile & compressive energies : 9877.08 + External energy: ((f,u^h)+(t,u^h))^0.5 : -0.210807 + Solving crack phase field at step=8 time=8 + Primary solution summary: L2-norm : 0.709403 + Max value : 0.79706 + Min value : -0.000244 + Range : 1.00024 + Dissipated energy: eps_d : 13.699 + L1-norm: |c^h| = (|c^h|) : 16.1776 + Normalized L1-norm: |c^h|/V : 0.647102 + Solving the elasto-dynamics problem... + step=9 time=9 + Primary solution summary: L2-norm : 1.33382 + Max X-component : 2.7 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -0.0664736 + Elastic strain energy: eps_e : 0.0216767 + Bulk energy: eps_b : 0.0216767 + Tensile & compressive energies : 12501.4 + External energy: ((f,u^h)+(t,u^h))^0.5 : -0.18231 + Solving crack phase field at step=9 time=9 + Primary solution summary: L2-norm : 0.709394 + Max value : 0.79706 + Min value : -0.000195 + Range : 1.0002 + Dissipated energy: eps_d : 13.7069 + L1-norm: |c^h| = (|c^h|) : 16.1761 + Normalized L1-norm: |c^h|/V : 0.647045 + Solving the elasto-dynamics problem... + step=10 time=10 + Primary solution summary: L2-norm : 1.48203 + Max X-component : 3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -0.0515521 + Elastic strain energy: eps_e : 0.0167381 + Bulk energy: eps_b : 0.0167381 + Tensile & compressive energies : 15434.2 + External energy: ((f,u^h)+(t,u^h))^0.5 : -0.160549 + Solving crack phase field at step=10 time=10 + Primary solution summary: L2-norm : 0.709388 + Max value : 0.79706 + Min value : -0.000159 + Range : 1.00016 + Dissipated energy: eps_d : 13.7125 + L1-norm: |c^h| = (|c^h|) : 16.1751 + Normalized L1-norm: |c^h|/V : 0.647003 + Solving the elasto-dynamics problem... + step=11 time=20 + Primary solution summary: L2-norm : 1.48203 + Max X-component : 3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -0.0219946 + Elastic strain energy: eps_e : 0.0109967 + Bulk energy: eps_b : 0.0109967 + Tensile & compressive energies : 15434.5 + External energy: ((f,u^h)+(t,u^h))^0.5 : -0.104868 + Solving crack phase field at step=11 time=20 + Primary solution summary: L2-norm : 0.709388 + Max value : 0.79706 + Min value : -0.000159 + Range : 1.00016 + Dissipated energy: eps_d : 13.7125 + L1-norm: |c^h| = (|c^h|) : 16.1751 + Normalized L1-norm: |c^h|/V : 0.647003 + Solving the elasto-dynamics problem... + step=12 time=30 + Primary solution summary: L2-norm : 1.48203 + Max X-component : 3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -0.0219926 + Elastic strain energy: eps_e : 0.0109963 + Bulk energy: eps_b : 0.0109963 + Tensile & compressive energies : 15434.5 + External energy: ((f,u^h)+(t,u^h))^0.5 : -0.104863 + Solving crack phase field at step=12 time=30 + Primary solution summary: L2-norm : 0.709388 + Max value : 0.79706 + Min value : -0.000159 + Range : 1.00016 + Dissipated energy: eps_d : 13.7125 + L1-norm: |c^h| = (|c^h|) : 16.1751 + Normalized L1-norm: |c^h|/V : 0.647003 + Solving the elasto-dynamics problem... + step=13 time=40 + Primary solution summary: L2-norm : 1.48203 + Max X-component : 3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -0.0219926 + Elastic strain energy: eps_e : 0.0109963 + Bulk energy: eps_b : 0.0109963 + Tensile & compressive energies : 15434.5 + External energy: ((f,u^h)+(t,u^h))^0.5 : -0.104863 + Solving crack phase field at step=13 time=40 + Primary solution summary: L2-norm : 0.709388 + Max value : 0.79706 + Min value : -0.000159 + Range : 1.00016 + Dissipated energy: eps_d : 13.7125 + L1-norm: |c^h| = (|c^h|) : 16.1751 + Normalized L1-norm: |c^h|/V : 0.647003 + Solving the elasto-dynamics problem... + step=14 time=50 + Primary solution summary: L2-norm : 1.48203 + Max X-component : 3 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -0.0219926 + Elastic strain energy: eps_e : 0.0109963 + Bulk energy: eps_b : 0.0109963 + Tensile & compressive energies : 15434.5 + External energy: ((f,u^h)+(t,u^h))^0.5 : -0.104863 + Solving crack phase field at step=14 time=50 + Primary solution summary: L2-norm : 0.709388 + Max value : 0.79706 + Min value : -0.000159 + Range : 1.00016 + Dissipated energy: eps_d : 13.7125 + L1-norm: |c^h| = (|c^h|) : 16.1751 + Normalized L1-norm: |c^h|/V : 0.647003 + Time integration completed. diff --git a/Test/Miehe71_FDporo.reg b/Test/Miehe71_FDporo.reg new file mode 100644 index 0000000..28443d7 --- /dev/null +++ b/Test/Miehe71_FDporo.reg @@ -0,0 +1,143 @@ +Miehe71.xinp -qstatic -poro + +Input file: Miehe71.xinp +Equation solver: 2 +Number of Gauss points: 4 +0. Parsing input file(s). +1. Poroelasticity solver +Parsing input file Miehe71.xinp +Parsing +Parsing + Generating linear geometry on unit parameter domain \[0,1]^2 + Length in X = 5 + Length in Y = 5 + Parsing + Parsing + Topology sets: Left (1,1,1D) + Right (1,2,1D) + TopBottom (1,3,1D) (1,4,1D) + Parsing + Refining P1 19 19 + Parsing +Parsing +Parsing +Parsing + Parsing + Material code 0: Poroelastic material, see "Problem definition:" below. + Parsing + Dirichlet code 2: (fixed) + Parsing + Dirichlet code 1: (fixed) + Parsing + Dirichlet code 1000001 (ramp): 3 \* Ramp(t,10) +Parsing +Parsing +Parsing +Parsing input file succeeded. +Equation solver: 2 +Number of Gauss points: 2 +2. Cahn-Hilliard solver +Parsing input file Miehe71.xinp +Parsing +Parsing + Parsing + Parsing + Topology sets: Left (1,1,1D) + Right (1,2,1D) + TopBottom (1,3,1D) (1,4,1D) + Parsing + Refining P1 19 19 + Parsing +Parsing + Initial crack function: abs(x-2.5) +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing input file succeeded. +Equation solver: 2 +Number of Gauss points: 2 +3. Time integration driver +Parsing input file Miehe71.xinp +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing +Parsing input file succeeded. +10. Preprocessing the finite element model: +11. Poroelasticity solver +Problem definition: +PoroElasticity: scaling = 0 useDynCoupling = false +Elasticity: 2D, gravity = 0 9.81 + Constitutive Properties: + Young's Modulus, E = 254.8 + Poisson's Ratio, nu = 0.3 + Densities: + Density of Fluid, rhof = 1000 + Density of Solid, rhos = 666.67 + Bulk Moduli: + Biot's coefficient, alpha = 1 + Biot's inverse modulus, M^-1 = 0.01 + Porosity, n = 0.1 + Permeability, K = 2.07e-09 2.07e-09 0 +Resolving Dirichlet boundary conditions + Constraining P1 E3 in direction(s) 2 + Constraining P1 E4 in direction(s) 2 + Constraining P1 E1 in direction(s) 1 + Constraining P1 E2 in direction(s) 1 code = 1000001 + >>> SAM model summary <<< +Number of elements 400 +Number of nodes 441 +Number of dofs 1323 +Number of constraints 21 +Number of unknowns 1239 +12. Cahn-Hilliard solver +Problem definition: +Cahn-Hilliard: 2D + Critical fracture energy density: 1 + Smearing factor: 0.1 + Max value in crack: 0.001 + Initial crack specified as a function. + Enforcing crack irreversibility using history buffer. +Resolving Dirichlet boundary conditions + >>> SAM model summary <<< +Number of elements 400 +Number of nodes 441 +Number of dofs 441 +Number of unknowns 441 +100. Starting the simulation + Solving the elasto-dynamics problem... + step=1 time=1 + Primary solution summary: L2-norm : 0.101242 + Max X-displacement : 0.3 + Max pressure : 1.72664e-07 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -30.87 + Solving crack phase field at step=1 time=1 + Primary solution summary: L2-norm : 0.749993 + Max value : 0.80195 + Min value : -0.125837 + Range : 1.12584 + Dissipated energy: eps_d : 9.21268 + L1-norm: |c^h| = (|c^h|) : 17.9004 + Normalized L1-norm: |c^h|/V : 0.716017 + Solving the elasto-dynamics problem... + step=2 time=2 + Primary solution summary: L2-norm : 0.217789 + Max X-displacement : 0.6 + Max pressure : 2.10171e-06 + Total reaction forces: Sum(R) : 0 0 + displacement\*reactions: (R,u) : -40.5157 + Solving crack phase field at step=2 time=2 + Primary solution summary: L2-norm : 0.711321 + Max value : 0.796012 + Min value : -0.0289993 + Range : 1.029 + Dissipated energy: eps_d : 12.1926 + L1-norm: |c^h| = (|c^h|) : 16.4654 + Normalized L1-norm: |c^h|/V : 0.658615 diff --git a/jenkins/build.sh b/jenkins/build.sh index 83175d3..9afe4f1 100755 --- a/jenkins/build.sh +++ b/jenkins/build.sh @@ -18,10 +18,11 @@ function clone_ifem { # Upstreams and revisions declare -a upstream -upstreams=(IFEM-Elasticity) +upstreams=(IFEM-Elasticity IFEM-PoroElasticity) declare -A upstreamRev upstreamRev[IFEM-Elasticity]=master +upstreamRev[IFEM-PoroElasticity]=master IFEM_REVISION=master if grep -qi "ifem=" <<< $ghprbCommentBody diff --git a/main_FractureDynamics.C b/main_FractureDynamics.C index bc8a965..3af9399 100644 --- a/main_FractureDynamics.C +++ b/main_FractureDynamics.C @@ -16,9 +16,12 @@ #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" #include "FractureArgs.h" #include "HHTSIM.h" @@ -26,6 +29,7 @@ #include "NewmarkNLSIM.h" #include "NonLinSIM.h" #include "ASMstruct.h" +#include "ASMmxBase.h" /*! @@ -44,6 +48,9 @@ public: //! \brief Empty destructor. virtual ~SIMDriver() {} + //! \brief Overrides the stop time that was read from the input file. + void setStopTime(double t) { Solver::tp.stopTime = t; } + protected: using Solver::parse; //! \brief Parses a data section from an XML element. @@ -76,12 +83,13 @@ private: /*! \brief Creates the combined fracture simulator and launches the simulation. \param[in] infile The input file to parse + \param[in] stopTime Stop time of the simulation (if non-negative) \param[in] context Input-file context for the time integrator */ template class Solver=SIMSolver> -int runCombined (char* infile, const char* context) +int runCombined (char* infile, double stopTime, const char* context) { IFEM::cout <<"\n\n0. Parsing input file(s)." <<"\n========================="<< std::endl; @@ -104,6 +112,9 @@ int runCombined (char* infile, const char* context) if (!readModel(solver,infile)) return 1; + if (stopTime >= 0.0) + solver.setStopTime(stopTime); + IFEM::cout <<"\n\n10. Preprocessing the finite element model:" <<"\n==========================================="<< std::endl; @@ -134,43 +145,62 @@ 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* context) { - typedef SIMDynElasticity ElSolver; - typedef SIMPhaseField PhaseSolver; - - const char* contx = Integrator::inputContext; if (args.integrator == 3 && args.coupling == 2) { typedef SIMFractureQstatic Coupler; - return runCombined(args.inpfile,contx); + return runCombined(args.inpfile, + args.stopT, + context); } else if (args.integrator == 3 && args.coupling == 3) { typedef SIMFractureMiehe Coupler; - return runCombined(args.inpfile,contx); + return runCombined(args.inpfile, + args.stopT, + context); } else { typedef SIMFracture Coupler; - return runCombined(args.inpfile,contx); + return runCombined(args.inpfile, + args.stopT, + context); } } +/*! + \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] stopTime Stop time of the simulation (if non-negative) \param[in] context Input-file context for the time integrator */ -template -int runStandAlone (char* infile, const char* context) +template +int runStandAlone (char* infile, double stopTime, const char* context) { - typedef SIMDynElasticity SIMElastoDynamics; + typedef SIMDynElasticity SIMElastoDynamics; IFEM::cout <<"\n\n0. Parsing input file(s)." <<"\n========================="<< std::endl; @@ -185,6 +215,9 @@ int runStandAlone (char* infile, const char* context) if (!readModel(solver,infile)) return 1; + if (stopTime >= 0.0) + solver.setStopTime(stopTime); + IFEM::cout <<"\n\n10. Preprocessing the finite element model:" <<"\n==========================================="<< std::endl; @@ -207,17 +240,64 @@ 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, + args.stopT, + context); +#else + return 99; // Built without the poroelastic coupling +#endif + + return runStandAlone>(args.inpfile, + args.stopT, + 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 runSimulator3(args); + return runSimulator5(args,context); +} + + +/*! + \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); } @@ -235,7 +315,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); } } @@ -266,7 +346,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: @@ -294,6 +374,7 @@ int main (int argc, char** argv) FractureArgs args; Elastic::planeStrain = true; + ASMmxBase::Type = ASMmxBase::NONE; IFEM::Init(argc,argv,"Fracture dynamics solver"); for (int i = 1; i < argc; i++) @@ -301,10 +382,14 @@ 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) + else if (!strncmp(argv[i],"-dbgEl",6) && i < argc-1) FractureElasticNorm::dbgElm = atoi(argv[++i]); + else if (!strncmp(argv[i],"-stopT",6) && i < argc-1) + args.stopT = atof(argv[++i]); else if (!args.inpfile) args.parseFile(argv[i],i); else @@ -314,10 +399,11 @@ 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"; + <<" [-nw ]]\n [-hdf5] [-principal] [-stopTime ]\n"; return 0; } @@ -325,7 +411,10 @@ int main (int argc, char** argv) IFEM::getOptions().discretization = ASM::LRSpline; IFEM::cout <<"\nInput file: "<< args.inpfile; - IFEM::getOptions().print(IFEM::cout) << std::endl; + IFEM::getOptions().print(IFEM::cout); + if (args.stopT >= 0.0) + IFEM::cout <<"\nSimulation stop time: "<< args.stopT; + IFEM::cout << std::endl; if (args.dim == 2) return runSimulator(args);