diff --git a/docs/examples/DynamicsDiagonalPeriodic/laminate/eta=00.10__alpha=0.20/CMakeLists.txt b/docs/examples/DynamicsDiagonalPeriodic/laminate/eta=00.10__alpha=0.20/CMakeLists.txt index e69c369..e1dbb6a 100644 --- a/docs/examples/DynamicsDiagonalPeriodic/laminate/eta=00.10__alpha=0.20/CMakeLists.txt +++ b/docs/examples/DynamicsDiagonalPeriodic/laminate/eta=00.10__alpha=0.20/CMakeLists.txt @@ -1,43 +1,53 @@ cmake_minimum_required(VERSION 3.1) # define a project name project(example) # set optimization level set(CMAKE_BUILD_TYPE Release) add_definitions(-DNDEBUG) # set C++ standard set(CMAKE_CXX_STANDARD 14) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) +# option to enable OpenMP +option(OPENMP "Use OpenMP" ON) +if(OPENMP) + find_package(OpenMP) + if(OPENMP_FOUND) + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + endif() +endif() + # load pkg-config find_package(PkgConfig) # load HDF5 find_package(HDF5 COMPONENTS CXX REQUIRED) include_directories(${HDF5_INCLUDE_DIRS}) set(LINK_LIBS ${HDF5_LIBS} ${HDF5_LIBRARIES}) # load header-only modules using pkg-config # - cppmat pkg_check_modules(CPPMAT REQUIRED cppmat) include_directories(${CPPMAT_INCLUDE_DIRS}) # - GooseFEM pkg_check_modules(GOOSEFEM REQUIRED GooseFEM) include_directories(${GOOSEFEM_INCLUDE_DIRS}) # - GooseMaterial pkg_check_modules(GOOSEMATERIAL REQUIRED GooseMaterial) include_directories(${GOOSEMATERIAL_INCLUDE_DIRS}) # - eigen3 pkg_check_modules(EIGEN3 REQUIRED eigen3) include_directories(${EIGEN3_INCLUDE_DIRS}) # - HDF5pp pkg_check_modules(HDF5PP REQUIRED HDF5pp) include_directories(${HDF5PP_INCLUDE_DIRS}) # create executable add_executable(${PROJECT_NAME} main.cpp) # link libraries target_link_libraries(${PROJECT_NAME} ${LINK_LIBS}) diff --git a/docs/examples/DynamicsDiagonalPeriodic/laminate/eta=00.10__alpha=0.20/main.cpp b/docs/examples/DynamicsDiagonalPeriodic/laminate/eta=00.10__alpha=0.20/main.cpp index d976896..6f2f0ce 100644 --- a/docs/examples/DynamicsDiagonalPeriodic/laminate/eta=00.10__alpha=0.20/main.cpp +++ b/docs/examples/DynamicsDiagonalPeriodic/laminate/eta=00.10__alpha=0.20/main.cpp @@ -1,184 +1,217 @@ -#define _USE_MATH_DEFINES // to use "M_PI" from "math.h" - -#include -#include #include #include #include -#include +#include +#include // ------------------------------------------------------------------------------------------------- -using MatS = GooseFEM::MatS; -using MatD = GooseFEM::MatD; -using ColD = GooseFEM::ColD; - -using vec = cppmat::cartesian2d::vector ; -using T2 = cppmat::cartesian2d::tensor2 ; -using T2s = cppmat::cartesian2d::tensor2s; -using T2d = cppmat::cartesian2d::tensor2d; - -namespace GM = GooseMaterial::AmorphousSolid::LinearStrain::Elastic::Cartesian2d; +using ColD = GooseFEM::ColD; +using T2 = cppmat::cartesian2d::tensor2; +using Mat = GooseMaterial::AmorphousSolid::LinearStrain::Elastic::Cartesian2d::Material; // ================================================================================================= -class Quadrature +class Material { public: - T2s eps, epsdot, sig; - size_t nhard; - GM::Material hard, soft, rayleigh; + // class variables + // --------------- + + // strain/stress, parameters, output variables + cppmat::matrix eps, epsdot, sig, rho, alpha, Epot; + + // dimensions + size_t nelem, nne=4, ndim=2, nip=4; + + // constitutive response + std::vector material; + + // damping + Mat rayleigh; + + // class functions + // --------------- - double Ebar, Vbar; + // constructor + Material(size_t nelem, size_t nhard); - Quadrature(size_t nhard, double eta); + // compute stress for one integration point + void updated_eps (size_t e, size_t k); + void updated_epsdot(size_t e, size_t k); - double density (size_t elem, size_t k, double V); - void stressStrain (size_t elem, size_t k, double V); - void stressStrainRate (size_t elem, size_t k, double V); - void stressStrainPost (size_t elem, size_t k, double V); - void stressStrainRatePost(size_t elem, size_t k, double V); + // compute post variables + void post(); }; // ------------------------------------------------------------------------------------------------- -Quadrature::Quadrature(size_t _nhard, double eta) +Material::Material(size_t _nelem, size_t _nhard) { - nhard = _nhard; - hard = GM::Material(100.,10.); - soft = GM::Material(100., 1.); - rayleigh = GM::Material(eta ,eta); -} + // copy from input + nelem = _nelem; -// ------------------------------------------------------------------------------------------------- + // allocate symmetric tensors and scalars of each integration point + epsdot.resize({nelem,nip,3}); + eps .resize({nelem,nip,3}); + sig .resize({nelem,nip,3}); + Epot .resize({nelem,nip }); -double Quadrature::density(size_t elem, size_t k, double V) -{ - return 1.0; -} -// ------------------------------------------------------------------------------------------------- + // allocate and set density and non-Galilean damping coefficient of each nodal integration point + rho .resize({nelem,nne}); rho .setConstant(1.0); + alpha.resize({nelem,nne}); alpha.setConstant(0.2); -void Quadrature::stressStrain(size_t elem, size_t k, double V) -{ - if ( elem < nhard ) sig = hard.stress(eps); - else sig = soft.stress(eps); + // constitutive response per element + for ( size_t e = 0 ; e < nelem ; ++e ) + { + if ( e < _nhard ) + { + Mat mat = Mat(100.,10.); + material.push_back(mat); + } + else + { + Mat mat = Mat(100.,1.); + material.push_back(mat); + } + } + + // damping + rayleigh = Mat(10.,.1); } + // ------------------------------------------------------------------------------------------------- -void Quadrature::stressStrainRate(size_t elem, size_t k, double V) +void Material::updated_eps(size_t e, size_t k) { - sig = rayleigh.stress(epsdot); + // local views of the global arrays (speeds up indexing, and increases readability) + cppmat::cartesian2d::tensor2s Eps, Epsdot, Sig; + + // pointer to stress/strain + Epsdot.map(&epsdot(e,k)); + Eps .map(&eps (e,k)); + Sig .map(&sig (e,k)); + + // compute stress + Sig.copy( material[e].stress(Eps) + rayleigh.stress(Epsdot) ); } + // ------------------------------------------------------------------------------------------------- -void Quadrature::stressStrainPost(size_t elem, size_t k, double V) +void Material::updated_epsdot(size_t e, size_t k) { - Vbar += V; + // local views of the global arrays (speeds up indexing, and increases readability) + cppmat::cartesian2d::tensor2s Eps, Epsdot, Sig; + + // pointer to stress/strain + Epsdot.map(&epsdot(e,k)); + Eps .map(&eps (e,k)); + Sig .map(&sig (e,k)); - if ( elem < nhard ) Ebar += hard.energy(eps) * V; - else Ebar += soft.energy(eps) * V; + // compute stress + Sig.copy( material[e].stress(Eps) + rayleigh.stress(Epsdot) ); } // ------------------------------------------------------------------------------------------------- -void Quadrature::stressStrainRatePost(size_t elem, size_t k, double V) +void Material::post() { +#pragma omp parallel +{ + // local views of the global arrays (speeds up indexing, and increases readability) + cppmat::cartesian2d::tensor2s Eps; + + // loop over all elements / integration points + #pragma omp for + for ( size_t e = 0 ; e < nelem ; ++e ) + { + for ( size_t k = 0 ; k < nip ; ++k ) + { + // pointer to stress/strain + Eps.map(&eps(e,k)); + + // compute energy + Epot(e,k) = material[e].energy(Eps); + } + } +} } // ================================================================================================= +using Mesh = GooseFEM::Mesh::Quad4::Regular; +using Element = GooseFEM::Element::Diagonal::SmallStrain::Quad4; +using Simulation = GooseFEM::Dynamics::Diagonal::Periodic; + +// ================================================================================================= + int main() { // set simulation parameters - double T = 5000.; // total time - double dt = 1.e-1; // time increment + double T = 20. ; // total time + double dt = 1.e-2; // time increment size_t nx = 40 ; // number of elements in both directions double gamma = .05 ; // displacement step - double eta = 0.1 ; // damping coefficient - double Gbar = 1.0 ; // equivalent shear modulus - double rho = 1.0 ; // density - // background damping - double alpha = std::pow(2.,.5)*2.*M_PI/static_cast(nx)*std::pow(Gbar/rho,.5)*rho; - std::cout << alpha << std::endl; // class which provides the mesh - GooseFEM::Mesh::Quad4::Regular mesh(nx,nx,1.); - // reference node + Mesh mesh(nx,nx,1.); + // extract information size_t nodeOrigin = mesh.nodeOrigin(); + size_t nelem = mesh.nelem(); + size_t nhard = nx*nx/4; - // class which provides the constitutive response at each quadrature point - auto quadrature = std::make_shared(nx*nx/4,eta); - - // class which provides the response of each element - using Elem = GooseFEM::Dynamics::Diagonal::LinearStrain::Quad4; - auto elem = std::make_shared(quadrature); - - // class which provides the system and an increment - GooseFEM::Dynamics::Diagonal::Periodic sim( - elem, + // simulation class + Simulation sim = Simulation( + std::make_unique(std::make_unique(nelem,nhard),nelem), mesh.coor(), mesh.conn(), mesh.dofsPeriodic(), - dt, - alpha + dt ); // define update in macroscopic deformation gradient T2 dFbar(0.); dFbar(0,1) = gamma; // update the displacement according to the macroscopic deformation gradient update for ( size_t i = 0 ; i < sim.nnode ; ++i ) for ( size_t j = 0 ; j < sim.ndim ; ++j ) for ( size_t k = 0 ; k < sim.ndim ; ++k ) - sim.u(i,j) += dFbar(j,k) * ( sim.x0(i,k) - sim.x0(nodeOrigin,k) ); - - // compute the externally applied strain energy - double Eext = std::pow(40.,2.) * (.25*100.+.75*1.) * std::pow(.5*gamma,2.); - size_t inc; + sim.u(i,j) += dFbar(j,k) * ( sim.x(i,k) - sim.x(nodeOrigin,k) ); // output variables - ColD Epot(static_cast(T/dt)); Epot.setZero(); // potential energy + ColD Epot(static_cast(T/dt)); Epot.setZero(); ColD Ekin(static_cast(T/dt)); Ekin.setZero(); ColD t (static_cast(T/dt)); t .setZero(); // loop over increments - for ( inc = 0 ; inc < static_cast(Epot.size()) ; ++inc ) + for ( size_t inc = 0 ; inc < static_cast(Epot.size()) ; ++inc ) { // - compute increment sim.velocityVerlet(); - // - post: energy based on nodes - // -- store time + // - store time t(inc) = sim.t; - // -- kinetic energy + + // - store total kinetic energy for ( size_t i = 0 ; i < sim.ndof ; ++i ) - Ekin(inc) += .5 / sim.Minv(i) * std::pow( sim.V(i) , 2. ); - // -- potential energy - quadrature->Ebar = 0.0; - quadrature->Vbar = 0.0; - sim.post(); - Epot(inc) = quadrature->Ebar; - - // check stopping criterion - if ( inc >= 10 ) if ( Ekin(inc)/Eext < 1.e-6 and Ekin(inc-10)/Eext < 1.e-6 ) break; - } + Ekin(inc) += .5 * sim.M(i) * std::pow(sim.V(i),2.); - // truncate to simulation size - Epot.conservativeResize(inc); - Ekin.conservativeResize(inc); - t .conservativeResize(inc); + // - store total potential energy + sim.elem->mat->post(); + Epot(inc) = sim.elem->mat->Epot.average(sim.elem->V) * sim.elem->V.sum(); + } // write to output file - H5p::File f = H5p::File("example.hdf5"); + H5p::File f = H5p::File("example.hdf5","w"); f.write("/global/Epot",Epot ); f.write("/global/Ekin",Ekin ); f.write("/global/t" ,t ); f.write("/mesh/coor" ,mesh.coor()); f.write("/mesh/conn" ,mesh.conn()); f.write("/mesh/disp" ,sim.u ); return 0; } diff --git a/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/Verlet/CMakeLists.txt b/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/Verlet/CMakeLists.txt index e69c369..e1dbb6a 100644 --- a/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/Verlet/CMakeLists.txt +++ b/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/Verlet/CMakeLists.txt @@ -1,43 +1,53 @@ cmake_minimum_required(VERSION 3.1) # define a project name project(example) # set optimization level set(CMAKE_BUILD_TYPE Release) add_definitions(-DNDEBUG) # set C++ standard set(CMAKE_CXX_STANDARD 14) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) +# option to enable OpenMP +option(OPENMP "Use OpenMP" ON) +if(OPENMP) + find_package(OpenMP) + if(OPENMP_FOUND) + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") + endif() +endif() + # load pkg-config find_package(PkgConfig) # load HDF5 find_package(HDF5 COMPONENTS CXX REQUIRED) include_directories(${HDF5_INCLUDE_DIRS}) set(LINK_LIBS ${HDF5_LIBS} ${HDF5_LIBRARIES}) # load header-only modules using pkg-config # - cppmat pkg_check_modules(CPPMAT REQUIRED cppmat) include_directories(${CPPMAT_INCLUDE_DIRS}) # - GooseFEM pkg_check_modules(GOOSEFEM REQUIRED GooseFEM) include_directories(${GOOSEFEM_INCLUDE_DIRS}) # - GooseMaterial pkg_check_modules(GOOSEMATERIAL REQUIRED GooseMaterial) include_directories(${GOOSEMATERIAL_INCLUDE_DIRS}) # - eigen3 pkg_check_modules(EIGEN3 REQUIRED eigen3) include_directories(${EIGEN3_INCLUDE_DIRS}) # - HDF5pp pkg_check_modules(HDF5PP REQUIRED HDF5pp) include_directories(${HDF5PP_INCLUDE_DIRS}) # create executable add_executable(${PROJECT_NAME} main.cpp) # link libraries target_link_libraries(${PROJECT_NAME} ${LINK_LIBS}) diff --git a/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/Verlet/main.cpp b/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/Verlet/main.cpp index cee2242..5f729e9 100644 --- a/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/Verlet/main.cpp +++ b/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/Verlet/main.cpp @@ -1,193 +1,193 @@ #include #include #include #include #include // ------------------------------------------------------------------------------------------------- using ColD = GooseFEM::ColD; using T2 = cppmat::cartesian2d::tensor2; using Mat = GooseMaterial::AmorphousSolid::LinearStrain::Elastic::Cartesian2d::Material; // ================================================================================================= class Material { public: // class variables // --------------- // strain/stress, parameters, output variables cppmat::matrix eps, epsdot, sig, rho, alpha, Epot; // dimensions size_t nelem, nne=4, ndim=2, nip=4; // constitutive response std::vector material; // class functions // --------------- // constructor Material(size_t nelem, size_t nhard); // compute stress for one integration point void updated_eps (size_t e, size_t k); void updated_epsdot(size_t e, size_t k){}; // compute post variables void post(); }; // ------------------------------------------------------------------------------------------------- Material::Material(size_t _nelem, size_t _nhard) { // copy from input nelem = _nelem; // allocate symmetric tensors and scalars of each integration point eps .resize({nelem,nip,3}); sig .resize({nelem,nip,3}); Epot.resize({nelem,nip }); // allocate and set density and non-Galilean damping coefficient of each nodal integration point rho .resize({nelem,nne}); rho .setConstant(1.); alpha.resize({nelem,nne}); alpha.setConstant(0.); // constitutive response per element for ( size_t e = 0 ; e < nelem ; ++e ) { if ( e < _nhard ) { Mat mat = Mat(100.,10.); material.push_back(mat); } else { Mat mat = Mat(100.,1.); material.push_back(mat); } } } // ------------------------------------------------------------------------------------------------- void Material::updated_eps(size_t e, size_t k) { // local views of the global arrays (speeds up indexing, and increases readability) cppmat::cartesian2d::tensor2s Eps, Sig; // pointer to stress/strain Eps.map(&eps(e,k)); Sig.map(&sig(e,k)); // compute stress Sig.copy( material[e].stress(Eps) ); } // ------------------------------------------------------------------------------------------------- void Material::post() { #pragma omp parallel { // local views of the global arrays (speeds up indexing, and increases readability) cppmat::cartesian2d::tensor2s Eps; // loop over all elements / integration points #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) { for ( size_t k = 0 ; k < nip ; ++k ) { // pointer to stress/strain Eps.map(&eps(e,k)); // compute energy Epot(e,k) = material[e].energy(Eps); } } } } // ================================================================================================= using Mesh = GooseFEM::Mesh::Quad4::Regular; -using Element = GooseFEM::Element::Diagonal::LinearStrain::Quad4; +using Element = GooseFEM::Element::Diagonal::SmallStrain::Quad4; using Simulation = GooseFEM::Dynamics::Diagonal::Periodic; // ================================================================================================= int main() { // set simulation parameters double T = 20. ; // total time double dt = 1.e-2; // time increment size_t nx = 40 ; // number of elements in both directions double gamma = .05 ; // displacement step // class which provides the mesh Mesh mesh(nx,nx,1.); // extract information size_t nodeOrigin = mesh.nodeOrigin(); size_t nelem = mesh.nelem(); size_t nhard = nx*nx/4; // simulation class Simulation sim = Simulation( std::make_unique(std::make_unique(nelem,nhard),nelem), mesh.coor(), mesh.conn(), mesh.dofsPeriodic(), dt ); // define update in macroscopic deformation gradient T2 dFbar(0.); dFbar(0,1) = gamma; // update the displacement according to the macroscopic deformation gradient update for ( size_t i = 0 ; i < sim.nnode ; ++i ) for ( size_t j = 0 ; j < sim.ndim ; ++j ) for ( size_t k = 0 ; k < sim.ndim ; ++k ) sim.u(i,j) += dFbar(j,k) * ( sim.x(i,k) - sim.x(nodeOrigin,k) ); // output variables ColD Epot(static_cast(T/dt)); Epot.setZero(); ColD Ekin(static_cast(T/dt)); Ekin.setZero(); ColD t (static_cast(T/dt)); t .setZero(); // loop over increments for ( size_t inc = 0 ; inc < static_cast(Epot.size()) ; ++inc ) { // - compute increment sim.Verlet(); // - store time t(inc) = sim.t; // - store total kinetic energy for ( size_t i = 0 ; i < sim.ndof ; ++i ) Ekin(inc) += .5 * sim.M(i) * std::pow(sim.V(i),2.); // - store total potential energy sim.elem->mat->post(); Epot(inc) = sim.elem->mat->Epot.average(sim.elem->V) * sim.elem->V.sum(); } // write to output file H5p::File f = H5p::File("example.hdf5","w"); f.write("/global/Epot",Epot ); f.write("/global/Ekin",Ekin ); f.write("/global/t" ,t ); f.write("/mesh/coor" ,mesh.coor()); f.write("/mesh/conn" ,mesh.conn()); f.write("/mesh/disp" ,sim.u ); return 0; } diff --git a/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/velocityVerlet/main.cpp b/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/velocityVerlet/main.cpp index 461d450..8d16df0 100644 --- a/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/velocityVerlet/main.cpp +++ b/docs/examples/DynamicsDiagonalPeriodic/laminate/no_damping/velocityVerlet/main.cpp @@ -1,193 +1,193 @@ #include #include #include #include #include // ------------------------------------------------------------------------------------------------- using ColD = GooseFEM::ColD; using T2 = cppmat::cartesian2d::tensor2; using Mat = GooseMaterial::AmorphousSolid::LinearStrain::Elastic::Cartesian2d::Material; // ================================================================================================= class Material { public: // class variables // --------------- // strain/stress, parameters, output variables cppmat::matrix eps, epsdot, sig, rho, alpha, Epot; // dimensions size_t nelem, nne=4, ndim=2, nip=4; // constitutive response std::vector material; // class functions // --------------- // constructor Material(size_t nelem, size_t nhard); // compute stress for one integration point void updated_eps (size_t e, size_t k); void updated_epsdot(size_t e, size_t k){}; // compute post variables void post(); }; // ------------------------------------------------------------------------------------------------- Material::Material(size_t _nelem, size_t _nhard) { // copy from input nelem = _nelem; // allocate symmetric tensors and scalars of each integration point eps .resize({nelem,nip,3}); sig .resize({nelem,nip,3}); Epot.resize({nelem,nip }); // allocate and set density and non-Galilean damping coefficient of each nodal integration point rho .resize({nelem,nne}); rho .setConstant(1.); alpha.resize({nelem,nne}); alpha.setConstant(0.); // constitutive response per element for ( size_t e = 0 ; e < nelem ; ++e ) { if ( e < _nhard ) { Mat mat = Mat(100.,10.); material.push_back(mat); } else { Mat mat = Mat(100.,1.); material.push_back(mat); } } } // ------------------------------------------------------------------------------------------------- void Material::updated_eps(size_t e, size_t k) { // local views of the global arrays (speeds up indexing, and increases readability) cppmat::cartesian2d::tensor2s Eps, Sig; // pointer to stress/strain Eps.map(&eps(e,k)); Sig.map(&sig(e,k)); // compute stress Sig.copy( material[e].stress(Eps) ); } // ------------------------------------------------------------------------------------------------- void Material::post() { #pragma omp parallel { // local views of the global arrays (speeds up indexing, and increases readability) cppmat::cartesian2d::tensor2s Eps; // loop over all elements / integration points #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) { for ( size_t k = 0 ; k < nip ; ++k ) { // pointer to stress/strain Eps.map(&eps(e,k)); // compute energy Epot(e,k) = material[e].energy(Eps); } } } } // ================================================================================================= using Mesh = GooseFEM::Mesh::Quad4::Regular; -using Element = GooseFEM::Element::Diagonal::LinearStrain::Quad4; +using Element = GooseFEM::Element::Diagonal::SmallStrain::Quad4; using Simulation = GooseFEM::Dynamics::Diagonal::Periodic; // ================================================================================================= int main() { // set simulation parameters double T = 20. ; // total time double dt = 1.e-2; // time increment size_t nx = 40 ; // number of elements in both directions double gamma = .05 ; // displacement step // class which provides the mesh Mesh mesh(nx,nx,1.); // extract information size_t nodeOrigin = mesh.nodeOrigin(); size_t nelem = mesh.nelem(); size_t nhard = nx*nx/4; // simulation class Simulation sim = Simulation( std::make_unique(std::make_unique(nelem,nhard),nelem), mesh.coor(), mesh.conn(), mesh.dofsPeriodic(), dt ); // define update in macroscopic deformation gradient T2 dFbar(0.); dFbar(0,1) = gamma; // update the displacement according to the macroscopic deformation gradient update for ( size_t i = 0 ; i < sim.nnode ; ++i ) for ( size_t j = 0 ; j < sim.ndim ; ++j ) for ( size_t k = 0 ; k < sim.ndim ; ++k ) sim.u(i,j) += dFbar(j,k) * ( sim.x(i,k) - sim.x(nodeOrigin,k) ); // output variables ColD Epot(static_cast(T/dt)); Epot.setZero(); ColD Ekin(static_cast(T/dt)); Ekin.setZero(); ColD t (static_cast(T/dt)); t .setZero(); // loop over increments for ( size_t inc = 0 ; inc < static_cast(Epot.size()) ; ++inc ) { // - compute increment sim.velocityVerlet(); // - store time t(inc) = sim.t; // - store total kinetic energy for ( size_t i = 0 ; i < sim.ndof ; ++i ) Ekin(inc) += .5 * sim.M(i) * std::pow(sim.V(i),2.); // - store total potential energy sim.elem->mat->post(); Epot(inc) = sim.elem->mat->Epot.average(sim.elem->V) * sim.elem->V.sum(); } // write to output file H5p::File f = H5p::File("example.hdf5","w"); f.write("/global/Epot",Epot ); f.write("/global/Ekin",Ekin ); f.write("/global/t" ,t ); f.write("/mesh/coor" ,mesh.coor()); f.write("/mesh/conn" ,mesh.conn()); f.write("/mesh/disp" ,sim.u ); return 0; } diff --git a/src/GooseFEM/DynamicsDiagonalLinearStrainQuad4.h b/src/GooseFEM/DynamicsDiagonalLinearStrainQuad4.h index cbbdf7e..8bd80c3 100644 --- a/src/GooseFEM/DynamicsDiagonalLinearStrainQuad4.h +++ b/src/GooseFEM/DynamicsDiagonalLinearStrainQuad4.h @@ -1,416 +1,416 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ -#ifndef GOOSEFEM_DYNAMICS_DIAGONAL_QUAD4_H -#define GOOSEFEM_DYNAMICS_DIAGONAL_QUAD4_H +#ifndef GOOSEFEM_DYNAMICS_DIAGONAL_SMALLSTRAIN_QUAD4_H +#define GOOSEFEM_DYNAMICS_DIAGONAL_SMALLSTRAIN_QUAD4_H // ------------------------------------------------------------------------------------------------- #include "Macros.h" #include // ------------------------------------------------------------------------------------------------- namespace GooseFEM { namespace Element { namespace Diagonal { -namespace LinearStrain { +namespace SmallStrain { // ========================== N.B. most loops are unrolled for efficiency ========================== template class Quad4 { public: // class variables // --------------- // constitutive response per integration point, per element std::unique_ptr mat; // matrices to store the element data cppmat::matrix x, u, v, f, M, D, dNx, dNxi, w, V, dNxi_n, w_n, V_n; // dimensions size_t nelem, nip=4, nne=4, ndim=2; // signals for solvers bool changed_f=false, changed_M=true, changed_D=true; // class functions // --------------- // constructor Quad4(std::unique_ptr mat, size_t nelem); // recompute relevant quantities after "x", "u", or "v" have been externally updated void updated_x(); void updated_u(); void updated_v(); }; // ================================================================================================= template Quad4::Quad4(std::unique_ptr _mat, size_t _nelem) { // copy from input nelem = _nelem; mat = std::move(_mat); // allocate matrices // - x .resize({nelem,nne,ndim}); u .resize({nelem,nne,ndim}); v .resize({nelem,nne,ndim}); f .resize({nelem,nne,ndim}); // - M .resize({nelem,nne*ndim,nne*ndim}); D .resize({nelem,nne*ndim,nne*ndim}); // - dNx .resize({nelem,nip,nne,ndim}); // - V .resize({nelem,nip}); V_n .resize({nelem,nne}); // - dNxi .resize({nip,nne,ndim}); dNxi_n.resize({nne,nne,ndim}); // - w .resize({nip}); w_n .resize({nne}); // zero initialize matrices (only diagonal is written, no new zero-initialization necessary) M.setZero(); D.setZero(); // shape function gradient at all Gauss points, in local coordinates // - k == 0 dNxi(0,0,0) = -.25*(1.+1./std::sqrt(3.)); dNxi(0,0,1) = -.25*(1.+1./std::sqrt(3.)); dNxi(0,1,0) = +.25*(1.+1./std::sqrt(3.)); dNxi(0,1,1) = -.25*(1.-1./std::sqrt(3.)); dNxi(0,2,0) = +.25*(1.-1./std::sqrt(3.)); dNxi(0,2,1) = +.25*(1.-1./std::sqrt(3.)); dNxi(0,3,0) = -.25*(1.-1./std::sqrt(3.)); dNxi(0,3,1) = +.25*(1.+1./std::sqrt(3.)); // - k == 1 dNxi(1,0,0) = -.25*(1.+1./std::sqrt(3.)); dNxi(1,0,1) = -.25*(1.-1./std::sqrt(3.)); dNxi(1,1,0) = +.25*(1.+1./std::sqrt(3.)); dNxi(1,1,1) = -.25*(1.+1./std::sqrt(3.)); dNxi(1,2,0) = +.25*(1.-1./std::sqrt(3.)); dNxi(1,2,1) = +.25*(1.+1./std::sqrt(3.)); dNxi(1,3,0) = -.25*(1.-1./std::sqrt(3.)); dNxi(1,3,1) = +.25*(1.-1./std::sqrt(3.)); // - k == 2 dNxi(2,0,0) = -.25*(1.-1./std::sqrt(3.)); dNxi(2,0,1) = -.25*(1.-1./std::sqrt(3.)); dNxi(2,1,0) = +.25*(1.-1./std::sqrt(3.)); dNxi(2,1,1) = -.25*(1.+1./std::sqrt(3.)); dNxi(2,2,0) = +.25*(1.+1./std::sqrt(3.)); dNxi(2,2,1) = +.25*(1.+1./std::sqrt(3.)); dNxi(2,3,0) = -.25*(1.+1./std::sqrt(3.)); dNxi(2,3,1) = +.25*(1.-1./std::sqrt(3.)); // - k == 3 dNxi(3,0,0) = -.25*(1.-1./std::sqrt(3.)); dNxi(3,0,1) = -.25*(1.+1./std::sqrt(3.)); dNxi(3,1,0) = +.25*(1.-1./std::sqrt(3.)); dNxi(3,1,1) = -.25*(1.-1./std::sqrt(3.)); dNxi(3,2,0) = +.25*(1.+1./std::sqrt(3.)); dNxi(3,2,1) = +.25*(1.-1./std::sqrt(3.)); dNxi(3,3,0) = -.25*(1.+1./std::sqrt(3.)); dNxi(3,3,1) = +.25*(1.+1./std::sqrt(3.)); // integration point weight at all Gauss points w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; // shape function gradient at all nodes, in local coordinates // - k == 0 dNxi_n(0,0,0) = -0.5; dNxi_n(0,0,1) = -0.5; dNxi_n(0,1,0) = +0.5; dNxi_n(0,1,1) = 0.0; dNxi_n(0,2,0) = 0.0; dNxi_n(0,2,1) = 0.0; dNxi_n(0,3,0) = 0.0; dNxi_n(0,3,1) = +0.5; // - k == 1 dNxi_n(1,0,0) = -0.5; dNxi_n(1,0,1) = 0.0; dNxi_n(1,1,0) = +0.5; dNxi_n(1,1,1) = -0.5; dNxi_n(1,2,0) = 0.0; dNxi_n(1,2,1) = +0.5; dNxi_n(1,3,0) = 0.0; dNxi_n(1,3,1) = 0.0; // - k == 2 dNxi_n(2,0,0) = 0.0; dNxi_n(2,0,1) = 0.0; dNxi_n(2,1,0) = 0.0; dNxi_n(2,1,1) = -0.5; dNxi_n(2,2,0) = +0.5; dNxi_n(2,2,1) = +0.5; dNxi_n(2,3,0) = -0.5; dNxi_n(2,3,1) = 0.0; // - k == 3 dNxi_n(3,0,0) = 0.0; dNxi_n(3,0,1) = -0.5; dNxi_n(3,1,0) = 0.0; dNxi_n(3,1,1) = 0.0; dNxi_n(3,2,0) = +0.5; dNxi_n(3,2,1) = 0.0; dNxi_n(3,3,0) = -0.5; dNxi_n(3,3,1) = +0.5; // integration point weight at all nodes w_n(0) = 1.; w_n(1) = 1.; w_n(2) = 1.; w_n(3) = 1.; // Note, the above is a specialization of the following: // // - Shape function gradients // // dNxi(0,0) = -.25*(1.-xi(k,1)); dNxi(0,1) = -.25*(1.-xi(k,0)); // dNxi(1,0) = +.25*(1.-xi(k,1)); dNxi(1,1) = -.25*(1.+xi(k,0)); // dNxi(2,0) = +.25*(1.+xi(k,1)); dNxi(2,1) = +.25*(1.+xi(k,0)); // dNxi(3,0) = -.25*(1.+xi(k,1)); dNxi(3,1) = +.25*(1.-xi(k,0)); // // - Gauss point coordinates and weights // // xi(0,0) = -1./std::sqrt(3.); xi(0,1) = -1./std::sqrt(3.); w(0) = 1.; // xi(1,0) = +1./std::sqrt(3.); xi(1,1) = -1./std::sqrt(3.); w(1) = 1.; // xi(2,0) = +1./std::sqrt(3.); xi(2,1) = +1./std::sqrt(3.); w(2) = 1.; // xi(3,0) = -1./std::sqrt(3.); xi(3,1) = +1./std::sqrt(3.); w(3) = 1.; // // - Nodal coordinates and weights // // xi(0,0) = -1.; xi(0,1) = -1.; w(0) = 1.; // xi(1,0) = +1.; xi(1,1) = -1.; w(1) = 1.; // xi(2,0) = +1.; xi(2,1) = +1.; w(2) = 1.; // xi(3,0) = -1.; xi(3,1) = +1.; w(3) = 1.; } // ================================================================================================= template void Quad4::updated_x() { #pragma omp parallel { // intermediate quantities cppmat::cartesian2d::tensor2 J_, Jinv_; double Jdet_; // local views of the global arrays (speeds up indexing, and increases readability) cppmat::tiny::matrix2 M_, D_; cppmat::tiny::matrix2 dNxi_, dNx_, x_; cppmat::tiny::vector w_, V_; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) { // pointer to element positions x_.map(&x(e)); // nodal quadrature // ---------------- // pointer to element mass/damping matrix, nodal volume, and integration weight M_.map(&M (e)); D_.map(&D (e)); V_.map(&V_n(e)); w_.map(&w_n(0)); // loop over nodes, i.e. the integration points for ( size_t k = 0 ; k < nne ; ++k ) { // - pointer to the shape function gradients (local coordinates) dNxi_.map(&dNxi_n(k)); // - Jacobian // J(i,j) += dNxi(m,i) * xe(m,j) J_(0,0) = dNxi_(0,0)*x_(0,0) + dNxi_(1,0)*x_(1,0) + dNxi_(2,0)*x_(2,0) + dNxi_(3,0)*x_(3,0); J_(0,1) = dNxi_(0,0)*x_(0,1) + dNxi_(1,0)*x_(1,1) + dNxi_(2,0)*x_(2,1) + dNxi_(3,0)*x_(3,1); J_(1,0) = dNxi_(0,1)*x_(0,0) + dNxi_(1,1)*x_(1,0) + dNxi_(2,1)*x_(2,0) + dNxi_(3,1)*x_(3,0); J_(1,1) = dNxi_(0,1)*x_(0,1) + dNxi_(1,1)*x_(1,1) + dNxi_(2,1)*x_(2,1) + dNxi_(3,1)*x_(3,1); // - determinant of the Jacobian Jdet_ = J_.det(); // - integration point volume V_(k) = w_(k) * Jdet_; // - assemble element mass matrix // M(m+i,n+i) += N(m) * rho * V * N(n); M_(k*2 ,k*2 ) = mat->rho(e,k) * V_(k); M_(k*2+1,k*2+1) = mat->rho(e,k) * V_(k); // - assemble element non-Galilean damping matrix // D(m+i,n+i) += N(m) * alpha * V * N(n); D_(k*2 ,k*2 ) = mat->alpha(e,k) * V_(k); D_(k*2+1,k*2+1) = mat->alpha(e,k) * V_(k); } // Gaussian quadrature // ------------------- // pointer to element integration volume and weight V_.map(&V(e)); w_.map(&w(0)); // loop over Gauss points for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to the shape function gradients (local/global coordinates) dNxi_.map(&dNxi( k)); dNx_ .map(&dNx (e,k)); // - Jacobian // J(i,j) += dNxi(m,i) * xe(m,j) J_(0,0) = dNxi_(0,0)*x_(0,0) + dNxi_(1,0)*x_(1,0) + dNxi_(2,0)*x_(2,0) + dNxi_(3,0)*x_(3,0); J_(0,1) = dNxi_(0,0)*x_(0,1) + dNxi_(1,0)*x_(1,1) + dNxi_(2,0)*x_(2,1) + dNxi_(3,0)*x_(3,1); J_(1,0) = dNxi_(0,1)*x_(0,0) + dNxi_(1,1)*x_(1,0) + dNxi_(2,1)*x_(2,0) + dNxi_(3,1)*x_(3,0); J_(1,1) = dNxi_(0,1)*x_(0,1) + dNxi_(1,1)*x_(1,1) + dNxi_(2,1)*x_(2,1) + dNxi_(3,1)*x_(3,1); // - determinant and inverse of the Jacobian Jdet_ = J_.det(); Jinv_ = J_.inv(); // - integration point volume V_(k) = w_(k) * Jdet_; // - shape function gradients (global coordinates) // dNx(m,i) += Jinv(i,j) * dNxi(m,j) for ( size_t m = 0 ; m < nne ; ++m ) { dNx_(m,0) = Jinv_(0,0) * dNxi_(m,0) + Jinv_(0,1) * dNxi_(m,1); dNx_(m,1) = Jinv_(1,0) * dNxi_(m,0) + Jinv_(1,1) * dNxi_(m,1); } } } // set signals changed_f = false; changed_M = true; changed_D = true; } // #pragma omp parallel } // ================================================================================================= template void Quad4::updated_u() { #pragma omp parallel { // intermediate quantities cppmat::cartesian2d::tensor2 gradu_; cppmat::cartesian2d::tensor2s eps_, sig_; // local views of the global arrays (speeds up indexing, and increases readability) cppmat::tiny::matrix2 dNx_, u_, f_; cppmat::tiny::vector V_; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) { // pointer to element forces, displacements, and integration volume f_.map(&f(e)); u_.map(&u(e)); V_.map(&V(e)); // zero initialize forces f_.setZero(); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to the shape function gradients, strain and stress tensor (stored symmetric) dNx_.map(&dNx (e,k)); eps_.map(&mat->eps(e,k)); sig_.map(&mat->sig(e,k)); // - displacement gradient // gradu_(i,j) += dNx(m,i) * ue(m,j) gradu_(0,0) = dNx_(0,0)*u_(0,0) + dNx_(1,0)*u_(1,0) + dNx_(2,0)*u_(2,0) + dNx_(3,0)*u_(3,0); gradu_(0,1) = dNx_(0,0)*u_(0,1) + dNx_(1,0)*u_(1,1) + dNx_(2,0)*u_(2,1) + dNx_(3,0)*u_(3,1); gradu_(1,0) = dNx_(0,1)*u_(0,0) + dNx_(1,1)*u_(1,0) + dNx_(2,1)*u_(2,0) + dNx_(3,1)*u_(3,0); gradu_(1,1) = dNx_(0,1)*u_(0,1) + dNx_(1,1)*u_(1,1) + dNx_(2,1)*u_(2,1) + dNx_(3,1)*u_(3,1); // - strain (stored symmetric) // eps(i,j) = .5 * ( gradu_(i,j) + gradu_(j,i) ) eps_(0,0) = gradu_(0,0); eps_(0,1) = .5 * ( gradu_(0,1) + gradu_(1,0) ); eps_(1,1) = gradu_(1,1); // - constitutive response mat->updated_eps(e,k); // - assemble to element force // f(m,j) += dNx(m,i) * sig(i,j) * V; for ( size_t m = 0 ; m < nne ; ++m ) { f_(m,0) += dNx_(m,0) * sig_(0,0) * V_(k) + dNx_(m,1) * sig_(1,0) * V_(k); f_(m,1) += dNx_(m,0) * sig_(0,1) * V_(k) + dNx_(m,1) * sig_(1,1) * V_(k); } } } // set signals changed_f = true; changed_M = false; changed_D = false; } } // ================================================================================================= template void Quad4::updated_v() { #pragma omp parallel { // intermediate quantities cppmat::cartesian2d::tensor2 gradv_; cppmat::cartesian2d::tensor2s epsdot_, sig_; // local views of the global arrays (speeds up indexing, and increases readability) cppmat::tiny::matrix2 dNx_, v_, f_; cppmat::tiny::vector V_; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) { // pointer to element forces, displacements, and integration volume f_.map(&f(e)); v_.map(&u(e)); V_.map(&V(e)); // zero initialize forces f_.setZero(); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to the shape function gradients, strain-rate and stress tensor (stored symmetric) dNx_ .map(&dNx (e,k)); epsdot_.map(&mat->eps(e,k)); sig_ .map(&mat->sig(e,k)); // - displacement gradient // gradv_(i,j) += dNx(m,i) * ue(m,j) gradv_(0,0) = dNx_(0,0)*v_(0,0) + dNx_(1,0)*v_(1,0) + dNx_(2,0)*v_(2,0) + dNx_(3,0)*v_(3,0); gradv_(0,1) = dNx_(0,0)*v_(0,1) + dNx_(1,0)*v_(1,1) + dNx_(2,0)*v_(2,1) + dNx_(3,0)*v_(3,1); gradv_(1,0) = dNx_(0,1)*v_(0,0) + dNx_(1,1)*v_(1,0) + dNx_(2,1)*v_(2,0) + dNx_(3,1)*v_(3,0); gradv_(1,1) = dNx_(0,1)*v_(0,1) + dNx_(1,1)*v_(1,1) + dNx_(2,1)*v_(2,1) + dNx_(3,1)*v_(3,1); // - strain (stored symmetric) // epsdot(i,j) = .5 * ( gradv_(i,j) + gradv_(j,i) ) epsdot_(0,0) = gradv_(0,0); epsdot_(0,1) = .5 * ( gradv_(0,1) + gradv_(1,0) ); epsdot_(1,1) = gradv_(1,1); // - constitutive response mat->updated_epsdot(e,k); // - assemble to element force // f(m,j) += dNdx(m,i) * sig(i,j) * V; for ( size_t m = 0 ; m < nne ; ++m ) { f_(m,0) += dNx_(m,0) * sig_(0,0) * V_(k) + dNx_(m,1) * sig_(1,0) * V_(k); f_(m,1) += dNx_(m,0) * sig_(0,1) * V_(k) + dNx_(m,1) * sig_(1,1) * V_(k); } } } // set signals changed_f = true; changed_M = false; changed_D = false; } } // ================================================================================================= }}}} // namespace ... #endif