From 9ebd7914d1c361acd4dfa8525c203b33ee76631d Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 21 Apr 2016 13:51:27 +0200 Subject: [PATCH 01/16] Added: .gitignore file. Added: Linkage to poro-elasticity module. Changed: Updated README.md to reflect the dependency on PoroElasticity. --- .gitignore | 3 +++ CMakeLists.txt | 46 ++++++++++++++++++++++++++++++++-------------- README.md | 25 ++++++++++++++----------- 3 files changed, 49 insertions(+), 25 deletions(-) create mode 100644 .gitignore 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..d227af8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -21,30 +21,48 @@ 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}) endif() file(GLOB FracEl_SRCS [A-Z]*.C) 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}) @@ -88,7 +106,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/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 From b48370ceab491b30ac3a0934b7b2a2f65a705af0 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 21 Apr 2016 17:01:08 +0200 Subject: [PATCH 02/16] Added: New integrand class PoroFracture which inherits PoroElasticity and have a (pointer to a) FractureElasticity object as member. Changed: Let the parent class of SIMDynElasticity be a template argument. Added: Constructor for mixed problems. --- CMakeLists.txt | 6 +- PoroFracture.C | 166 +++++++++++++++++++++++++++++++++++++++++++++ PoroFracture.h | 98 ++++++++++++++++++++++++++ SIMDynElasticity.h | 76 +++++++++++---------- 4 files changed, 308 insertions(+), 38 deletions(-) create mode 100644 PoroFracture.C create mode 100644 PoroFracture.h diff --git a/CMakeLists.txt b/CMakeLists.txt index d227af8..5975bd5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,9 +46,13 @@ if(EXISTS ${POROELASTIC_DIR}) add_subdirectory(${POROELASTIC_DIR} PoroElasticity) endif() 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 CahnHilliard.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}) diff --git a/PoroFracture.C b/PoroFracture.C new file mode 100644 index 0000000..1e2ee56 --- /dev/null +++ b/PoroFracture.C @@ -0,0 +1,166 @@ +// $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 "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) +{ + m_mode = 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 + return fracEl->evalInt(elMat,fe,X); +} diff --git a/PoroFracture.h b/PoroFracture.h new file mode 100644 index 0000000..9d728fd --- /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 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); + + //! \brief Returns a pointer to the Gauss-point tensile energy array. + virtual const RealArray* getTensileEnergy() 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; + +private: + FractureElasticity* fracEl; //!< Integrand for tangent stiffness evaluation +}; + +#endif diff --git a/SIMDynElasticity.h b/SIMDynElasticity.h index 616d1d5..8ef17f0 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; } @@ -46,22 +54,12 @@ class SIMDynElasticity : public SIMElasticity if (++ncall == 1) // Avoiding infinite recursive calls dSim.printProblem(); else - this->SIMElasticity::printProblem(); + this->Sim::printProblem(); --ncall; } - //! \brief Registers fields for data output. - void registerFields(DataExporter& exporter) - { - int flag = DataExporter::PRIMARY; - if (!Dim::opt.pSolOnly) - flag |= DataExporter::SECONDARY; - exporter.registerField("u","solid displacement",DataExporter::SIM,flag); - exporter.setFieldValue("u",this,&dSim.getSolution()); - } - //! \brief Initializes the problem. - bool init(const TimeStep&) + virtual bool init(const TimeStep&) { dSim.initPrm(); dSim.initSol(dynamic_cast(&dSim) ? 3 : 1); @@ -72,21 +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) + virtual bool saveStep(const TimeStep& tp, int& nBlock) { double old = utl::zero_print_tol; utl::zero_print_tol = 1e-16; @@ -156,7 +143,7 @@ class SIMDynElasticity : public SIMElasticity } //! \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..."; @@ -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,6 +232,8 @@ class SIMDynElasticity : public SIMElasticity } } + //! \brief Returns a const reference to current solution vector. + 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(); } @@ -315,12 +304,25 @@ class SIMDynElasticity : public SIMElasticity static short int ncall = 0; if (++ncall == 1) // Avoiding infinite recursive calls result = dSim.parse(elem); - else if (!strcasecmp(elem->Value(),"elasticity") && !Dim::myProblem) + else if (!strcasecmp(elem->Value(),SIMElasticity::myContext.c_str())) { - std::string form("voigt"); - if (utl::getAttribute(elem,"formulation",form,true) && form != "voigt") - Dim::myProblem = new FractureElasticity(Dim::dimension); - result = this->SIMElasticity::parse(elem); + if (!Dim::myProblem) + { + if (this->getName() == "PoroElasticity") +#ifdef IFEM_HAS_POROELASTIC + Dim::myProblem = new PoroFracture(Dim::dimension); +#else + return false; +#endif + else + { + std::string formulation("voigt"); + utl::getAttribute(elem,"formulation",formulation,true); + if (formulation != "voigt") + Dim::myProblem = new FractureElasticity(Dim::dimension); + } + } + result = this->Sim::parse(elem); const TiXmlElement* child = elem->FirstChildElement(); for (; child; child = child->NextSiblingElement()) @@ -328,7 +330,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; } From 14de7c7e504f516756d8c90ace3af1ccc2480ace Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Fri, 22 Apr 2016 07:11:18 +0200 Subject: [PATCH 03/16] Added: Command-line option -poro to run porelastic fracture solver --- CMakeLists.txt | 3 ++- FractureArgs.C | 7 ++++++ FractureArgs.h | 1 + main_FractureDynamics.C | 48 ++++++++++++++++++++++++++++++----------- 4 files changed, 46 insertions(+), 13 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5975bd5..7a2acb0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -49,7 +49,8 @@ if(EXISTS ${POROELASTIC_DIR}) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DIFEM_HAS_POROELASTIC") endif() -set(FracEl_SRCS CahnHilliard.C FractureElasticity.C FractureElasticityVoigt.C) +set(FracEl_SRCS CahnHilliard.C FractureArgs.C + FractureElasticity.C FractureElasticityVoigt.C) if(EXISTS ${POROELASTIC_DIR}) set(FracEl_SRCS ${FracEl_SRCS} PoroFracture.C) endif() diff --git a/FractureArgs.C b/FractureArgs.C index 2b715f7..f7450cf 100644 --- a/FractureArgs.C +++ b/FractureArgs.C @@ -20,6 +20,7 @@ FractureArgs::FractureArgs () : SIMargsBase("fracturedynamics") { inpfile = nullptr; integrator = coupling = 1; + poroEl = false; } @@ -41,6 +42,10 @@ bool FractureArgs::parseArg (const char* argv) integrator = 4; else if (!strcmp(argv,"-oldHHT")) integrator = 5; + else if (!strncmp(argv,"-poro",5)) + poroEl = true; + else if (!strncmp(argv,"-explcr",7)) + expPhase = true; else if (!strncmp(argv,"-noadap",7)) adap = false; else @@ -78,6 +83,8 @@ bool FractureArgs::parse (const TiXmlElement* elem) integrator = 4; else if (!strcasecmp(child->Value(),"hilberhughestaylor")) integrator = 4; + else if (!strcasecmp(child->Value(),"poroelastic")) + poroEl = true; } return this->SIMargsBase::parse(elem); diff --git a/FractureArgs.h b/FractureArgs.h index a59e649..883b500 100644 --- a/FractureArgs.h +++ b/FractureArgs.h @@ -30,6 +30,7 @@ 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 //! \brief Default constructor. FractureArgs(); diff --git a/main_FractureDynamics.C b/main_FractureDynamics.C index bc8a965..f92bf2d 100644 --- a/main_FractureDynamics.C +++ b/main_FractureDynamics.C @@ -17,6 +17,9 @@ #include "SIMDynElasticity.h" #include "SIMPhaseField.h" #include "SIMFractureQstatic.h" +#ifdef IFEM_HAS_POROELASTIC +#include "SIMPoroElasticity.h" +#endif #include "SIMCoupledSI.h" #include "SIMSolver.h" #include "SIMSolverTS.h" @@ -134,15 +137,11 @@ int runCombined (char* infile, const char* context) \brief Determines whether the quasi-static semi-implicit driver is to be used. */ -template class Cpl, template class Solver=SIMSolver> -int runSimulator3 (const FractureArgs& args) +int runSimulator4 (const FractureArgs& args, const char* contx) { - typedef SIMDynElasticity ElSolver; - typedef SIMPhaseField PhaseSolver; - - const char* contx = Integrator::inputContext; if (args.integrator == 3 && args.coupling == 2) { typedef SIMFractureQstatic Coupler; @@ -211,13 +210,38 @@ int runStandAlone (char* infile, const char* 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; + typedef SIMPhaseField PhaseSolver; + + const char* context = Integrator::inputContext; + if (args.adap) - return runSimulator3(args); + return runSimulator4(args,context); + + return runSimulator4(args,context); +} - return runSimulator3(args); + +/*! + \brief Creates the combined fracture simulator and launches the simulation. +*/ + +template class Cpl> +int runSimulator2 (const FractureArgs& args) +{ + if (args.poroEl) +#ifdef IFEM_HAS_POROELASTIC + return runSimulator3,Cpl>(args); +#else + return 99; // Built without the poroelastic coupling +#endif + + return runSimulator3,Cpl>(args); } @@ -314,8 +338,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] [-[l|q]static|-GA|-HHT] [-adaptive]\n" + <<" [-lag|-spec|-LR] [-2D] [-nGauss ]\n [-nocrack|" + <<"-semiimplicit] [-[l|q]static|-GA|-HHT] [-poro] [-adaptive]\n" <<" [-vtf [-nviz ] [-nu ] [-nv ]]\n [-hdf5] [-principal]\n"; return 0; From 420c3e5c90d62654fc268d1e755ebba627610c14 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 26 Apr 2016 13:56:38 +0200 Subject: [PATCH 04/16] Added: FractureElasticity::setMode to set correct matrix indices --- FractureElasticity.C | 13 +++++++++++++ FractureElasticity.h | 4 ++++ PoroFracture.C | 2 +- 3 files changed, 18 insertions(+), 1 deletion(-) diff --git a/FractureElasticity.C b/FractureElasticity.C index 3573515..b50898f 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 diff --git a/FractureElasticity.h b/FractureElasticity.h index efdb46a..c422ea2 100644 --- a/FractureElasticity.h +++ b/FractureElasticity.h @@ -46,6 +46,10 @@ class FractureElasticity : public Elasticity //! \brief Sets the number of solution variables per node. void setVar(unsigned short int n) { npv = n; } + //! \brief Defines the solution mode before the element assembly is started. + //! \param[in] mode The solution mode to use + virtual void setMode(SIM::SolutionMode mode); + //! \brief Returns whether this norm has explicit boundary contributions. virtual bool hasBoundaryTerms() const { return m_mode != SIM::RECOVERY; } diff --git a/PoroFracture.C b/PoroFracture.C index 1e2ee56..62a6a13 100644 --- a/PoroFracture.C +++ b/PoroFracture.C @@ -51,7 +51,7 @@ void PoroFracture::setMaterial (Material* mat) void PoroFracture::setMode (SIM::SolutionMode mode) { - m_mode = mode; + this->PoroElasticity::setMode(mode); fracEl->setMode(mode); primsol.resize(fracEl->getNoSolutions()); } From c63b935f8ecdd61dd1a3840fd3ef01e9fb9a4084 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 26 Apr 2016 15:22:08 +0200 Subject: [PATCH 05/16] Added: Command-line argument -mixed and the option to run pure poroelasticity without the phase-field coupling, but still with the strain energy split --- main_FractureDynamics.C | 47 +++++++++++++++++++++++++++++++---------- 1 file changed, 36 insertions(+), 11 deletions(-) diff --git a/main_FractureDynamics.C b/main_FractureDynamics.C index f92bf2d..d303197 100644 --- a/main_FractureDynamics.C +++ b/main_FractureDynamics.C @@ -29,6 +29,7 @@ #include "NewmarkNLSIM.h" #include "NonLinSIM.h" #include "ASMstruct.h" +#include "ASMmxBase.h" /*! @@ -140,7 +141,7 @@ int runCombined (char* infile, const char* context) template class Cpl, template class Solver=SIMSolver> -int runSimulator4 (const FractureArgs& args, const char* contx) +int runSimulator5 (const FractureArgs& args, const char* contx) { if (args.integrator == 3 && args.coupling == 2) { @@ -166,10 +167,10 @@ int runSimulator4 (const FractureArgs& args, const char* contx) \param[in] context Input-file context for the time integrator */ -template +template int runStandAlone (char* infile, const char* context) { - typedef SIMDynElasticity SIMElastoDynamics; + typedef SIMDynElasticity SIMElastoDynamics; IFEM::cout <<"\n\n0. Parsing input file(s)." <<"\n========================="<< std::endl; @@ -206,6 +207,27 @@ int runStandAlone (char* infile, const char* context) } +/*! + \brief Determines whether the poroelastic simulation driver is to be used. +*/ + +template +int runSimulator4 (const FractureArgs& args, + const char* context = "newmarksolver") +{ + if (args.poroEl) +#ifdef IFEM_HAS_POROELASTIC + return runStandAlone>(args.inpfile, + context); +#else + return 99; // Built without the poroelastic coupling +#endif + + return runStandAlone>(args.inpfile, + context); +} + + /*! \brief Determines whether the adaptive simulation driver is to be used. */ @@ -220,9 +242,9 @@ int runSimulator3 (const FractureArgs& args) const char* context = Integrator::inputContext; if (args.adap) - return runSimulator4(args,context); + return runSimulator5(args,context); - return runSimulator4(args,context); + return runSimulator5(args,context); } @@ -259,7 +281,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); } } @@ -290,7 +312,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: @@ -318,6 +340,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++) @@ -325,6 +348,8 @@ int main (int argc, char** argv) ; // ignore the input file on the second pass else if (SIMoptions::ignoreOldOptions(argc,argv,i)) ; // ignore the obsolete option + else if (!strcmp(argv[i],"-mixed")) + ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1; else if (!strcmp(argv[i],"-principal")) Elasticity::wantPrincipalStress = true; else if (!strcmp(argv[i],"-dbgElm") && i < argc-1) @@ -337,10 +362,10 @@ int main (int argc, char** argv) if (!args.inpfile) { std::cout <<"usage: "<< argv[0] - <<" [-dense|-spr|-superlu[]|-samg|-petsc]\n" - <<" [-lag|-spec|-LR] [-2D] [-nGauss ]\n [-nocrack|" - <<"-semiimplicit] [-[l|q]static|-GA|-HHT] [-poro] [-adaptive]\n" - <<" [-vtf [-nviz ] [-nu ] [-nv [-dense|-spr|-superlu[]|-samg|-petsc]\n " + <<" [-lag|-spec|-LR] [-2D] [-mixed] [-nGauss ]\n [-nocra" + <<"ck|-semiimplicit] [-[l|q]static|-GA|-HHT] [-poro] [-adaptive]" + <<"\n [-vtf [-nviz ] [-nu ] [-nv ]]\n [-hdf5] [-principal]\n"; return 0; } From 656f8baba17d179fe36eb125b8209104ec574f41 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 28 Apr 2016 17:31:43 +0200 Subject: [PATCH 06/16] Added: Methods for evaluating the phase field and gradient at current point, and the broken permeability tensor based on phase field and deformation state --- FractureElasticity.C | 43 ++++++++++++++++++++++ FractureElasticity.h | 9 +++++ PoroFracture.C | 87 +++++++++++++++++++++++++++++++++++++++++++- PoroFracture.h | 16 ++++++++ 4 files changed, 154 insertions(+), 1 deletion(-) diff --git a/FractureElasticity.C b/FractureElasticity.C index b50898f..218983e 100644 --- a/FractureElasticity.C +++ b/FractureElasticity.C @@ -634,6 +634,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 c422ea2..33b0bda 100644 --- a/FractureElasticity.h +++ b/FractureElasticity.h @@ -100,6 +100,15 @@ class FractureElasticity : public Elasticity virtual bool evalSol(Vector& s, const Vectors& eV, const FiniteElement& fe, const Vec3& X, bool toLocal, Vec3* pdir) 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 index 62a6a13..033d125 100644 --- a/PoroFracture.C +++ b/PoroFracture.C @@ -13,12 +13,22 @@ #include "PoroFracture.h" #include "FractureElasticityVoigt.h" +#include "FiniteElement.h" +#include "Tensor.h" +#include "Vec3Oper.h" #include "Utilities.h" +#include "IFEM.h" +#include "tinyxml.h" PoroFracture::PoroFracture (unsigned short int n) : PoroElasticity(n) { fracEl = new FractureElasticityVoigt(this,n); + + L_per = 0.01; + d_min = 0.1; + Kc = 83.0; + eps = 50.0; } @@ -30,7 +40,21 @@ PoroFracture::~PoroFracture() bool PoroFracture::parse (const TiXmlElement* elem) { - return this->PoroElasticity::parse(elem) & fracEl->parse(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; } @@ -164,3 +188,64 @@ bool PoroFracture::evalElasticityMatrices (ElmMats& elMat, const Matrix&, // Evaluate tangent stiffness matrix and internal forces return fracEl->evalInt(elMat,fe,X); } + + +/*! + This method calculates the anisotropic peremeability for the broken solid + based on a Poiseuielle 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." +*/ + +double PoroFracture::formCrackedPermeabilityTensor (SymmTensor& Kcrack, + const Vectors& eV, + const FiniteElement& fe, + const Vec3& X) const +{ + 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.zero(); + 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 the symmetric tensor C^-1 - (C^-1*n0)otimes(C^-1*n0) + // (the term in the bracket [] of eq. (108) in Miehe2015pfm3) + Vec3 CigD = Kcrack*gradD; // C^-1*gradD + for (unsigned short int j = 1; j <= nsd; j++) + for (unsigned short int i = 1; i <= j; i++) + Kcrack(i,j) -= CigD(i)*CigD(j)/d2; + + // Compute the perpendicular crack stretch + // lambda = gradD*gradD / gradD*C^-1*gradD (see eq. (106)) + double lambda = sqrt(d2 / (gradD*CigD)); + double w = lambda*L_per - L_per; // Crack opening (see eq. (107)) + + // Compute the permeability tensor, scale by d^eps*Kc*w^2*J (see eq. (108)) + Kcrack *= pow(d,eps)*Kc*w*w*F.det(); + return w; +} diff --git a/PoroFracture.h b/PoroFracture.h index 9d728fd..275e525 100644 --- a/PoroFracture.h +++ b/PoroFracture.h @@ -91,8 +91,24 @@ class PoroFracture : public PoroElasticity 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 From 3c2131c39cd70714cbce69376e21eb6624cdc453 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Wed, 4 May 2016 07:23:12 +0200 Subject: [PATCH 07/16] Added: virtual method for computing the total permeability tensor --- PoroFracture.C | 20 ++++++++++++++++++++ PoroFracture.h | 10 ++++++++++ 2 files changed, 30 insertions(+) diff --git a/PoroFracture.C b/PoroFracture.C index 033d125..3a384f4 100644 --- a/PoroFracture.C +++ b/PoroFracture.C @@ -12,6 +12,7 @@ //============================================================================== #include "PoroFracture.h" +#include "PoroMaterial.h" #include "FractureElasticityVoigt.h" #include "FiniteElement.h" #include "Tensor.h" @@ -249,3 +250,22 @@ double PoroFracture::formCrackedPermeabilityTensor (SymmTensor& Kcrack, Kcrack *= pow(d,eps)*Kc*w*w*F.det(); return w; } + + +bool PoroFracture::formPermeabilityTensor (SymmTensor& K, + const Vectors& eV, + const FiniteElement& fe, + const Vec3& X) const +{ + if (this->formCrackedPermeabilityTensor(K,eV,fe,X) < 0.0) + return false; + + const PoroMaterial* pmat = dynamic_cast(material); + if (!pmat) return false; + + Vec3 permeability = pmat->getPermeability(X); + for (size_t i = 1; i <= K.dim(); i++) + K(i,i) += permeability(i); + + return true; +} diff --git a/PoroFracture.h b/PoroFracture.h index 275e525..4a05a1e 100644 --- a/PoroFracture.h +++ b/PoroFracture.h @@ -91,6 +91,16 @@ class PoroFracture : public PoroElasticity 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 From f525bfe36e1d8194c129507e0bd5bd7b96fa1da5 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 10 May 2016 10:58:07 +0200 Subject: [PATCH 08/16] Added: VTF-output of crack opening and permeability coefficients --- FractureElasticity.C | 28 ++++++++++++++++----------- FractureElasticity.h | 5 +++++ PoroFracture.C | 46 ++++++++++++++++++++++++++++++++++++++++++++ PoroFracture.h | 18 +++++++++++++++++ 4 files changed, 86 insertions(+), 11 deletions(-) diff --git a/FractureElasticity.C b/FractureElasticity.C index 218983e..ad80d8e 100644 --- a/FractureElasticity.C +++ b/FractureElasticity.C @@ -521,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()); @@ -535,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); } diff --git a/FractureElasticity.h b/FractureElasticity.h index 33b0bda..e3c96f7 100644 --- a/FractureElasticity.h +++ b/FractureElasticity.h @@ -100,6 +100,11 @@ 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 + 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 diff --git a/PoroFracture.C b/PoroFracture.C index 3a384f4..d7f8e4e 100644 --- a/PoroFracture.C +++ b/PoroFracture.C @@ -269,3 +269,49 @@ bool PoroFracture::formPermeabilityTensor (SymmTensor& K, return true; } + + +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)); + + const PoroMaterial* pmat = dynamic_cast(material); + if (!pmat) return false; + + Vec3 permeability = pmat->getPermeability(X); + for (size_t i = 1; i <= Kc.dim(); i++) + s.push_back(Kc(i,i)+permeability(i)); + + return true; +} diff --git a/PoroFracture.h b/PoroFracture.h index 4a05a1e..393a812 100644 --- a/PoroFracture.h +++ b/PoroFracture.h @@ -79,9 +79,27 @@ class PoroFracture : public PoroElasticity 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 From 5d79cf2369784045d701f43a652749642b3bc85f Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Wed, 25 May 2016 07:14:21 +0200 Subject: [PATCH 09/16] Fixed: Ignore crack stretch slightly less than zero. Fixed: Removed extraction of element velocity and accelereation in PoroFracture::initElement, it is done by the parent class method. --- PoroFracture.C | 64 +++++++++++++++----------------------------------- 1 file changed, 19 insertions(+), 45 deletions(-) diff --git a/PoroFracture.C b/PoroFracture.C index d7f8e4e..12bb5bd 100644 --- a/PoroFracture.C +++ b/PoroFracture.C @@ -113,33 +113,12 @@ bool PoroFracture::initElement (const std::vector& MNPC, // (1 = pressure, 2 = pressure rate, 3 = phase field) elmInt.vec.resize(primsol.size()+3); - // Extract element displacement and pressure vectors + // Extract element displacement, velocity, acceleration 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; + // Extract element phase-field vector + return fracEl->initElement(MNPC,elmInt); } @@ -152,27 +131,13 @@ bool PoroFracture::initElement (const std::vector& MNPC, // (1 = pressure, 2 = pressure rate, 3 = phase field) elmInt.vec.resize(primsol.size()+3); - // Extract element displacement and pressure vectors + // Extract element displacement, velocity, acceleration and pressure vectors if (!this->PoroElasticity::initElement(MNPC,elem_sizes,basis_sizes,elmInt)) return false; - // Extract element phase-field, velocity and acceleration vectors + // Extract element phase-field vector 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; + return fracEl->initElement(std::vector(MNPC.begin(),uEnd),elmInt); } @@ -192,8 +157,8 @@ bool PoroFracture::evalElasticityMatrices (ElmMats& elMat, const Matrix&, /*! - This method calculates the anisotropic peremeability for the broken solid - based on a Poiseuielle flow approximation of the fluid flow in the crack. + 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." */ @@ -244,11 +209,20 @@ double PoroFracture::formCrackedPermeabilityTensor (SymmTensor& Kcrack, // Compute the perpendicular crack stretch // lambda = gradD*gradD / gradD*C^-1*gradD (see eq. (106)) double lambda = sqrt(d2 / (gradD*CigD)); - double w = lambda*L_per - L_per; // Crack opening (see eq. (107)) +#if INT_DEBUG > 3 + std::cout <<"PoroFracture::formCrackedPermeabilityTensor(X = "<< X + <<") lambda = "<< lambda << std::endl; +#endif + if (lambda <= 1.0) + { + Kcrack.zero(); + return 0.0; + } // Compute the permeability tensor, scale by d^eps*Kc*w^2*J (see eq. (108)) + double w = lambda*L_per - L_per; // Crack opening (see eq. (107)) Kcrack *= pow(d,eps)*Kc*w*w*F.det(); - return w; + return w < 0.0 ? 0.0 : w; } From 908cd7a47e16fe27485b61d855f3731c8a2dfd94 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Fri, 27 May 2016 19:06:28 +0200 Subject: [PATCH 10/16] Added: class SIMExplPhaseField representing an explicit phase field, and new command-line options -explcrack to activate it --- CMakeLists.txt | 2 +- FractureArgs.C | 4 +- FractureArgs.h | 1 + SIMExplPhaseField.C | 95 +++++++++++++++++++++++++++++++++++ SIMExplPhaseField.h | 108 ++++++++++++++++++++++++++++++++++++++++ main_FractureDynamics.C | 33 +++++++++--- 6 files changed, 233 insertions(+), 10 deletions(-) create mode 100644 SIMExplPhaseField.C create mode 100644 SIMExplPhaseField.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 7a2acb0..fdd50eb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -49,7 +49,7 @@ if(EXISTS ${POROELASTIC_DIR}) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DIFEM_HAS_POROELASTIC") endif() -set(FracEl_SRCS CahnHilliard.C FractureArgs.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) diff --git a/FractureArgs.C b/FractureArgs.C index f7450cf..7229a65 100644 --- a/FractureArgs.C +++ b/FractureArgs.C @@ -20,7 +20,7 @@ FractureArgs::FractureArgs () : SIMargsBase("fracturedynamics") { inpfile = nullptr; integrator = coupling = 1; - poroEl = false; + poroEl = expPhase = false; } @@ -85,6 +85,8 @@ bool FractureArgs::parse (const TiXmlElement* elem) 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 883b500..ac0ef5c 100644 --- a/FractureArgs.h +++ b/FractureArgs.h @@ -31,6 +31,7 @@ class FractureArgs : public SIMargsBase char integrator; char coupling; //!< Coupling flag (0: none, 1: staggered, 2: semi-implicit) bool poroEl; //!< If \e true, use the poroelastic solver + bool expPhase; //!< If \e true, use an explicit phase field //! \brief Default constructor. FractureArgs(); diff --git a/SIMExplPhaseField.C b/SIMExplPhaseField.C new file mode 100644 index 0000000..5c33717 --- /dev/null +++ b/SIMExplPhaseField.C @@ -0,0 +1,95 @@ +// $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; + } + + size_t numNod = myOwner->getNoNodes(1); + phaseField.resize(numNod); + + for (size_t n = 1; n <= numNod; n++) + { + // Works since we are using basis 1 + Vec4 X = myOwner->getNodeCoord(n); + X.t = tp.time.t; + phaseField[n-1] = (*phaseFunc)(X); + } + + return true; +} + + +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/main_FractureDynamics.C b/main_FractureDynamics.C index d303197..ce1d28b 100644 --- a/main_FractureDynamics.C +++ b/main_FractureDynamics.C @@ -16,6 +16,7 @@ #include "SIM3D.h" #include "SIMDynElasticity.h" #include "SIMPhaseField.h" +#include "SIMExplPhaseField.h" #include "SIMFractureQstatic.h" #ifdef IFEM_HAS_POROELASTIC #include "SIMPoroElasticity.h" @@ -141,7 +142,7 @@ int runCombined (char* infile, const char* context) template class Cpl, template class Solver=SIMSolver> -int runSimulator5 (const FractureArgs& args, const char* contx) +int runSimulator6 (const FractureArgs& args, const char* contx) { if (args.integrator == 3 && args.coupling == 2) { @@ -161,6 +162,22 @@ int runSimulator5 (const FractureArgs& args, const char* contx) } +/*! + \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 @@ -237,14 +254,13 @@ template DynElSolver; - typedef SIMPhaseField PhaseSolver; const char* context = Integrator::inputContext; if (args.adap) - return runSimulator5(args,context); + return runSimulator5(args,context); - return runSimulator5(args,context); + return runSimulator5(args,context); } @@ -362,10 +378,11 @@ int main (int argc, char** argv) if (!args.inpfile) { std::cout <<"usage: "<< argv[0] - <<" [-dense|-spr|-superlu[]|-samg|-petsc]\n " - <<" [-lag|-spec|-LR] [-2D] [-mixed] [-nGauss ]\n [-nocra" - <<"ck|-semiimplicit] [-[l|q]static|-GA|-HHT] [-poro] [-adaptive]" - <<"\n [-vtf [-nviz ] [-nu ] [-nv [-dense|-spr|-superlu[]|-samg|-petsc]\n" + <<" [-lag|-spec|-LR] [-2D] [-mixed] [-nGauss ]\n" + <<" [-nocrack|-explcrack|-semiimplicit]" + <<" [-[l|q]static|-GA|-HHT] [-poro] [-adaptive]\n" + <<" [-vtf [-nviz ] [-nu ] [-nv ]]\n [-hdf5] [-principal]\n"; return 0; } From 7cd00cf0f5474f0a461451d8d730137d94b67131 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 31 May 2016 18:11:22 +0200 Subject: [PATCH 11/16] Added: Command-line option -stopTime overriding the input file setting --- FractureArgs.C | 1 + FractureArgs.h | 3 ++- main_FractureDynamics.C | 43 +++++++++++++++++++++++++++++++---------- 3 files changed, 36 insertions(+), 11 deletions(-) diff --git a/FractureArgs.C b/FractureArgs.C index 7229a65..d1d5606 100644 --- a/FractureArgs.C +++ b/FractureArgs.C @@ -21,6 +21,7 @@ FractureArgs::FractureArgs () : SIMargsBase("fracturedynamics") inpfile = nullptr; integrator = coupling = 1; poroEl = expPhase = false; + stopT = -1.0; } diff --git a/FractureArgs.h b/FractureArgs.h index ac0ef5c..5e25afb 100644 --- a/FractureArgs.h +++ b/FractureArgs.h @@ -32,6 +32,7 @@ class FractureArgs : public SIMargsBase 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(); @@ -42,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/main_FractureDynamics.C b/main_FractureDynamics.C index ce1d28b..3af9399 100644 --- a/main_FractureDynamics.C +++ b/main_FractureDynamics.C @@ -22,7 +22,6 @@ #include "SIMPoroElasticity.h" #endif #include "SIMCoupledSI.h" -#include "SIMSolver.h" #include "SIMSolverTS.h" #include "FractureArgs.h" #include "HHTSIM.h" @@ -49,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. @@ -81,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; @@ -109,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; @@ -142,22 +148,28 @@ int runCombined (char* infile, const char* context) template class Cpl, template class Solver=SIMSolver> -int runSimulator6 (const FractureArgs& args, const char* contx) +int runSimulator6 (const FractureArgs& args, const char* context) { 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); } } @@ -181,11 +193,12 @@ int runSimulator5 (const FractureArgs& args, const char* 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) +int runStandAlone (char* infile, double stopTime, const char* context) { typedef SIMDynElasticity SIMElastoDynamics; @@ -202,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; @@ -235,12 +251,14 @@ int runSimulator4 (const FractureArgs& args, 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); } @@ -368,8 +386,10 @@ int main (int argc, char** argv) 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 @@ -383,7 +403,7 @@ int main (int argc, char** argv) <<" [-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; } @@ -391,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); From 637b1aac9d49f45a22edc61e0a7ac9286f9c2503 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 31 May 2016 18:11:22 +0200 Subject: [PATCH 12/16] Added: Test case from Miehe2015pfm3 --- CMakeLists.txt | 7 + Test/Miehe71-explcrack_FD.reg | 91 ++++++++ Test/Miehe71-explcrack_FDporo.reg | 98 ++++++++ Test/Miehe71-nocrack_FD.reg | 209 +++++++++++++++++ Test/Miehe71-nocrack_FDporo.reg | 147 ++++++++++++ Test/Miehe71.xinp | 66 ++++++ Test/Miehe71_FD.reg | 357 ++++++++++++++++++++++++++++++ Test/Miehe71_FDporo.reg | 128 +++++++++++ 8 files changed, 1103 insertions(+) create mode 100644 Test/Miehe71-explcrack_FD.reg create mode 100644 Test/Miehe71-explcrack_FDporo.reg create mode 100644 Test/Miehe71-nocrack_FD.reg create mode 100644 Test/Miehe71-nocrack_FDporo.reg create mode 100644 Test/Miehe71.xinp create mode 100644 Test/Miehe71_FD.reg create mode 100644 Test/Miehe71_FDporo.reg diff --git a/CMakeLists.txt b/CMakeLists.txt index fdd50eb..cc9aaa2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -89,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() @@ -104,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) 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..9efba89 --- /dev/null +++ b/Test/Miehe71_FDporo.reg @@ -0,0 +1,128 @@ +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 From 18571135cf0bbae0ed4e9220136f464cd92ebb99 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 7 Jul 2016 09:01:00 +0200 Subject: [PATCH 13/16] Added: Projection of explicit phase field function --- SIMExplPhaseField.C | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/SIMExplPhaseField.C b/SIMExplPhaseField.C index 5c33717..0c6ee18 100644 --- a/SIMExplPhaseField.C +++ b/SIMExplPhaseField.C @@ -70,18 +70,9 @@ bool SIMExplPhaseField::solveStep (TimeStep& tp, bool) return false; } - size_t numNod = myOwner->getNoNodes(1); - phaseField.resize(numNod); - - for (size_t n = 1; n <= numNod; n++) - { - // Works since we are using basis 1 - Vec4 X = myOwner->getNodeCoord(n); - X.t = tp.time.t; - phaseField[n-1] = (*phaseFunc)(X); - } - - return true; + phaseField.resize(myOwner->getNoNodes(1)); + return myOwner->project(phaseField,phaseFunc,1,0,1, + SIMoptions::GLOBAL,tp.time.t); } From 7aa4cb2ea576bf492d4f756ede9d99d1de7faa97 Mon Sep 17 00:00:00 2001 From: Eivind Fonn Date: Mon, 23 Jan 2017 12:00:03 +0100 Subject: [PATCH 14/16] Fixed: Computation of cracked permeability tensor. Filter out the pressure variables in the getElementSolution method. --- FractureElasticity.C | 1 + FractureElasticity.h | 3 +- PoroFracture.C | 105 ++++++++++++++++++++++++++++++---------- Test/Miehe71_FDporo.reg | 15 ++++++ 4 files changed, 97 insertions(+), 27 deletions(-) diff --git a/FractureElasticity.C b/FractureElasticity.C index ad80d8e..38c226f 100644 --- a/FractureElasticity.C +++ b/FractureElasticity.C @@ -52,6 +52,7 @@ FractureElasticity::FractureElasticity (IntegrandBase* parent, parent->registerVector("phasefield",&myCVec); // Assuming second vector is pressure, third vector is pressure velocity eC = 3; // and fourth vector is the phase field + npv = nsd+1; // Account for pressure variables in the primary solution vector } diff --git a/FractureElasticity.h b/FractureElasticity.h index e3c96f7..4603d8e 100644 --- a/FractureElasticity.h +++ b/FractureElasticity.h @@ -103,7 +103,8 @@ class FractureElasticity : public Elasticity //! \brief Retrieves the element solution vectors. //! \param[out] eV Element solution vectors //! \param[in] MNPC Nodal point correspondance for the basis function values - bool getElementSolution(Vectors& eV, const std::vector& MNPC) const; + 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 diff --git a/PoroFracture.C b/PoroFracture.C index 12bb5bd..79bc613 100644 --- a/PoroFracture.C +++ b/PoroFracture.C @@ -22,9 +22,55 @@ #include "tinyxml.h" +/*! + \brief Class representing the integrand of elasticity problems with fracture. + + \details This sub-class overrides the getElementSolution method, to account + for that the primary solution vector, as extracted from the patch level, + also contains the pressure variables in addition to the displacements. +*/ + +class PoroFractureElasticity : public FractureElasticityVoigt +{ +public: + //! \brief Constructor for integrands with a parent integrand. + //! \param parent The parent integrand of this one + //! \param[in] n Number of spatial dimensions + PoroFractureElasticity(IntegrandBase* parent, unsigned short int n) + : FractureElasticityVoigt(parent,n) {} + //! \brief Empty destructor. + virtual ~PoroFractureElasticity() {} + + //! \brief Retrieves the element solution vectors. + //! \param[out] eV Element solution vectors + //! \param[in] MNPC Nodal point correspondance for the basis function values + virtual bool getElementSolution(Vectors& eV, + const std::vector& MNPC) const + { + if (!this->FractureElasticityVoigt::getElementSolution(eV,MNPC)) + return false; + + // Filter out the pressure components + // FIXME: Mixed + size_t nen = MNPC.size(); + if (eV.front().size() == nen*(nsd+1)) + { + Vector temp(eV.front()); + Vector& actual = eV.front(); + actual.resize(nen*nsd); + for (size_t n = 0; n < nen; n++) + for (unsigned short int i = 0; i < nsd; i++) + actual[nsd*n+i] = temp[(nsd+1)*n+i]; + } + + return true; + } +}; + + PoroFracture::PoroFracture (unsigned short int n) : PoroElasticity(n) { - fracEl = new FractureElasticityVoigt(this,n); + fracEl = new PoroFractureElasticity(this,n); L_per = 0.01; d_min = 0.1; @@ -161,6 +207,7 @@ bool PoroFracture::evalElasticityMatrices (ElmMats& elMat, const Matrix&, 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, @@ -168,6 +215,13 @@ double PoroFracture::formCrackedPermeabilityTensor (SymmTensor& Kcrack, 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) @@ -179,7 +233,7 @@ double PoroFracture::formCrackedPermeabilityTensor (SymmTensor& Kcrack, else if (d < d_min) { // The crack does not affect the permeability tensor at this point - Kcrack.zero(); + Kcrack = perm; return 0.0; } @@ -199,15 +253,24 @@ double PoroFracture::formCrackedPermeabilityTensor (SymmTensor& Kcrack, if (Kcrack.rightCauchyGreen(F).inverse() == 0.0) return -3.0; - // Compute the symmetric tensor C^-1 - (C^-1*n0)otimes(C^-1*n0) + // 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) -= CigD(i)*CigD(j)/d2; + Kcrack(i,j) -= a2od2*CigD(i)*CigD(j); // Compute the perpendicular crack stretch - // lambda = gradD*gradD / gradD*C^-1*gradD (see eq. (106)) + // 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 @@ -215,14 +278,19 @@ double PoroFracture::formCrackedPermeabilityTensor (SymmTensor& Kcrack, #endif if (lambda <= 1.0) { - Kcrack.zero(); + Kcrack = perm; return 0.0; } // Compute the permeability tensor, scale by d^eps*Kc*w^2*J (see eq. (108)) - double w = lambda*L_per - L_per; // Crack opening (see eq. (107)) - Kcrack *= pow(d,eps)*Kc*w*w*F.det(); - return w < 0.0 ? 0.0 : w; + // 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; } @@ -231,17 +299,7 @@ bool PoroFracture::formPermeabilityTensor (SymmTensor& K, const FiniteElement& fe, const Vec3& X) const { - if (this->formCrackedPermeabilityTensor(K,eV,fe,X) < 0.0) - return false; - - const PoroMaterial* pmat = dynamic_cast(material); - if (!pmat) return false; - - Vec3 permeability = pmat->getPermeability(X); - for (size_t i = 1; i <= K.dim(); i++) - K(i,i) += permeability(i); - - return true; + return this->formCrackedPermeabilityTensor(K,eV,fe,X) >= 0.0; } @@ -279,13 +337,8 @@ bool PoroFracture::evalSol (Vector& s, const FiniteElement& fe, SymmTensor Kc(nsd); s.push_back(this->formCrackedPermeabilityTensor(Kc,eV,fe,X)); - - const PoroMaterial* pmat = dynamic_cast(material); - if (!pmat) return false; - - Vec3 permeability = pmat->getPermeability(X); for (size_t i = 1; i <= Kc.dim(); i++) - s.push_back(Kc(i,i)+permeability(i)); + s.push_back(Kc(i,i)); return true; } diff --git a/Test/Miehe71_FDporo.reg b/Test/Miehe71_FDporo.reg index 9efba89..28443d7 100644 --- a/Test/Miehe71_FDporo.reg +++ b/Test/Miehe71_FDporo.reg @@ -126,3 +126,18 @@ Number of unknowns 441 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 From 1645c9a26070bf66d65fa8649e78497b262a256b Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 5 Apr 2018 14:40:29 +0200 Subject: [PATCH 15/16] Changed: Set npv correctly for mixed problems in the constructor, such that we don't have to override the getLocalIntegral methods. Using override keyword for reimplemented virtual methods. --- FractureElasticity.C | 1 - FractureElasticity.h | 3 --- PoroFracture.C | 30 +++++++----------------------- PoroFracture.h | 14 ++------------ SIMDynElasticity.h | 25 +++++++++++++------------ 5 files changed, 22 insertions(+), 51 deletions(-) diff --git a/FractureElasticity.C b/FractureElasticity.C index 38c226f..ad80d8e 100644 --- a/FractureElasticity.C +++ b/FractureElasticity.C @@ -52,7 +52,6 @@ FractureElasticity::FractureElasticity (IntegrandBase* parent, parent->registerVector("phasefield",&myCVec); // Assuming second vector is pressure, third vector is pressure velocity eC = 3; // and fourth vector is the phase field - npv = nsd+1; // Account for pressure variables in the primary solution vector } diff --git a/FractureElasticity.h b/FractureElasticity.h index 4603d8e..dbcecec 100644 --- a/FractureElasticity.h +++ b/FractureElasticity.h @@ -43,9 +43,6 @@ 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); diff --git a/PoroFracture.C b/PoroFracture.C index 79bc613..b6df31b 100644 --- a/PoroFracture.C +++ b/PoroFracture.C @@ -35,9 +35,11 @@ class PoroFractureElasticity : public FractureElasticityVoigt public: //! \brief Constructor for integrands with a parent integrand. //! \param parent The parent integrand of this one - //! \param[in] n Number of spatial dimensions - PoroFractureElasticity(IntegrandBase* parent, unsigned short int n) - : FractureElasticityVoigt(parent,n) {} + //! \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() {} @@ -68,9 +70,9 @@ public: }; -PoroFracture::PoroFracture (unsigned short int n) : PoroElasticity(n) +PoroFracture::PoroFracture (unsigned short int n, bool m) : PoroElasticity(n,m) { - fracEl = new PoroFractureElasticity(this,n); + fracEl = new PoroFractureElasticity(this, n, m ? n : n+1); L_per = 0.01; d_min = 0.1; @@ -134,24 +136,6 @@ void PoroFracture::initIntegration (size_t nGp, size_t 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) { diff --git a/PoroFracture.h b/PoroFracture.h index 393a812..3d4b20a 100644 --- a/PoroFracture.h +++ b/PoroFracture.h @@ -30,7 +30,8 @@ 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); + //! \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(); @@ -53,17 +54,6 @@ class PoroFracture : public PoroElasticity //! \param[in] nBp Total number of boundary integration points virtual void initIntegration(size_t nGp, size_t nBp); - //! \brief Returns a local integral contribution object for the given element. - //! \param[in] nen Number of nodes on element - //! \param[in] neumann Whether or not we are assembling Neumann BCs - virtual LocalIntegral* getLocalIntegral(size_t nen, - size_t, bool neumann) const; - //! \brief Returns a local integral contribution object for the given element. - //! \param[in] nen Number of nodes on element for each basis - //! \param[in] neumann Whether or not we are assembling Neumann BCs - virtual LocalIntegral* getLocalIntegral(const std::vector& nen, - size_t, bool neumann) const; - using PoroElasticity::initElement; //! \brief Initializes current element for numerical integration. //! \param[in] MNPC Matrix of nodal point correspondance for current element diff --git a/SIMDynElasticity.h b/SIMDynElasticity.h index 8ef17f0..282774d 100644 --- a/SIMDynElasticity.h +++ b/SIMDynElasticity.h @@ -48,7 +48,7 @@ class SIMDynElasticity : public Sim 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 @@ -59,7 +59,7 @@ class SIMDynElasticity : public Sim } //! \brief Initializes the problem. - virtual bool init(const TimeStep&) + bool init(const TimeStep&, bool = false) override { dSim.initPrm(); dSim.initSol(dynamic_cast(&dSim) ? 3 : 1); @@ -73,7 +73,7 @@ class SIMDynElasticity : public Sim //! \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 - virtual 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; @@ -127,12 +127,12 @@ class SIMDynElasticity : public Sim } //! \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(); @@ -143,7 +143,7 @@ class SIMDynElasticity : public Sim } //! \brief Computes the solution for the current time step. - virtual bool solveStep(TimeStep& tp) + bool solveStep(TimeStep& tp) override { if (Dim::msgLevel >= 1) IFEM::cout <<"\n Solving the elasto-dynamics problem..."; @@ -233,9 +233,9 @@ class SIMDynElasticity : public Sim } //! \brief Returns a const reference to current solution vector. - virtual const Vector& getSolution(int i) const { return dSim.getSolution(i); } + 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) @@ -289,16 +289,16 @@ class SIMDynElasticity : public Sim 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; @@ -310,7 +310,8 @@ class SIMDynElasticity : public Sim { if (this->getName() == "PoroElasticity") #ifdef IFEM_HAS_POROELASTIC - Dim::myProblem = new PoroFracture(Dim::dimension); + Dim::myProblem = new PoroFracture(Dim::dimension, + this->mixedProblem()); #else return false; #endif From 0ba13d8eb10f2d7ff6fcd693b45160202831bfe3 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 24 Aug 2018 09:03:28 +0200 Subject: [PATCH 16/16] Added: ifem-poroelasticity to jenkins upstreams --- jenkins/build.sh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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