diff --git a/docs/DynamicsDiagonal.rst b/docs/DynamicsDiagonal.rst index 8982e56..9d859db 100644 --- a/docs/DynamicsDiagonal.rst +++ b/docs/DynamicsDiagonal.rst @@ -1,277 +1,277 @@ **************************** GooseFEM::Dynamics::Diagonal **************************** [:download:`source: DynamicsDiagonal.h <../src/GooseFEM/DynamicsDiagonal.h>`, :download:`source: DynamicsDiagonal.cpp <../src/GooseFEM/DynamicsDiagonal.cpp>`] Overview ======== Principle --------- The philosophy is to provide some structure to efficiently run a finite element simulation which remains customizable. Even more customization can be obtained by copy/pasting the source and modifying it to your need. The idea that is followed involves a hierarchy of three classes, whereby a class that is higher in hierarchy writes to some field of the class that is lower in hierarchy, runs a function, and reads from some field. In general: * **Discretized system** (``GooseFEM::Dynamics::Diagonal::Periodic``, ``GooseFEM::Dynamics::Diagonal::SemiPeriodic``). * Defines the discretized system. * Writes element positions, displacements, and velocity of all elements to the *element definition*. * Assembles the diagonal (inverse) mass matrix, the displacement dependent forces, the velocity dependent forces, and the diagonal damping matrix from the element arrays computed in *element definition*. * Provides time integrators. * **Element definition** (``GooseFEM::Dynamics::Diagonal::SmallStrain::Qaud4``) Provides the element arrays by performing numerical quadrature. At the integration point the strain and strain-rate are computed and constitutive response is probed from the quadrature point definition. * **Quadrature point definition** This class is not provided, and should be provided by the user. Example ------- A simple example is: -[:download:`source: examples/DynamicsDiagonalPeriodic/laminate/no_damping/Verlet/main.cpp `] +[:download:`source: examples/DynamicsDiagonal/Periodic/laminate/no_damping/Verlet/main.cpp `] .. code-block:: cpp #include #include #include #include // ------------------------------------------------------------------- using MatS = GooseFEM::MatS; using MatD = GooseFEM::MatD; using ColD = GooseFEM::ColD; using T2 = cppmat::cartesian2d::tensor2 ; using T2s = cppmat::cartesian2d::tensor2s; namespace GM = GooseMaterial::AmorphousSolid::SmallStrain::Elastic::Cartesian2d; // =================================================================== class Quadrature { public: T2s eps, epsdot, sig; size_t nhard; GM::Material hard, soft; Quadrature(size_t nhard); double density (size_t e, size_t k); void stressStrain (size_t e, size_t k); void stressStrainRate (size_t e, size_t k); void stressStrainPost (size_t e, size_t k); void stressStrainRatePost(size_t e, size_t k); }; // ------------------------------------------------------------------- Quadrature::Quadrature(size_t _nhard) { nhard = _nhard; hard = GM::Material(100.,10.); soft = GM::Material(100., 1.); } // ------------------------------------------------------------------- double Quadrature::density(size_t elem, size_t k, double V) { return 1.0; } // ------------------------------------------------------------------- void Quadrature::stressStrain(size_t elem, size_t k, double V) { if ( elem < nhard ) sig = hard.stress(eps); else sig = soft.stress(eps); } // ------------------------------------------------------------------- void Quadrature::stressStrainRate(size_t elem, size_t k, double V) { } // ------------------------------------------------------------------- void Quadrature::stressStrainPost(size_t elem, size_t k, double V) { Vbar += V; if ( elem < nhard ) Ebar += hard.energy(eps) * V; else Ebar += soft.energy(eps) * V; } // ------------------------------------------------------------------- void Quadrature::stressStrainRatePost(size_t elem, size_t k, double V) { } // =================================================================== int main() { // class which provides the mesh GooseFEM::Mesh::Quad4::Regular mesh(40,40,1.); // class which provides the constitutive response at each quadrature point auto quadrature = std::make_shared(40*40/4); // class which provides the response of each element using Elem = GooseFEM::Dynamics::Diagonal::SmallStrain::Quad4; auto elem = std::make_shared(quadrature); // class which provides the system and an increment GooseFEM::Dynamics::Diagonal::Periodic sim( elem, mesh.coor(), mesh.conn(), mesh.dofsPeriodic(), 1.e-2, 0.0 ); // loop over increments for ( ... ) { // - set displacement of fixed DOFs ... // - compute time increment sim.Verlet(); // - post-process quadrature->Ebar = 0.0; quadrature->Vbar = 0.0; sim.post(); ... } return 0; } Pseudo-code ----------- What is happening inside ``Verlet`` is evaluating the forces (and the mass matrix), and updating the displacements by solving the system. In pseudo-code: * Mass matrix: .. code-block:: python sim.computeMinv(): { for e in elements: sim->elem->xe(i,j) = ... sim->elem->ue(i,j) = ... sim->elem->computeM(e): { for k in integration-points: sim->elem->M(...,...) += ... * sim->elem->quad->density(e,k,V) } M(...) += sim->elem->M(i,i) } * Displacement dependent force: .. code-block:: python sim.computeFu(): { for e in elements: sim->elem->xe(i,j) = ... sim->elem->ue(i,j) = ... sim->elem->computeFu(e): { for k in integration-points: sim->elem->quad->eps(i,j) = ... sim->elem->quad->stressStrain(e,k,V) sim->elem->fu(...) += ... * sim->elem->quad->sig(i,j) } Fu(...) += sim->elem->fu(i) } * Velocity dependent force: .. code-block:: python sim.computeFv(): { for e in elements: sim->elem->xe(i,j) = ... sim->elem->ue(i,j) = ... sim->elem->ve(i,j) = ... sim->elem->computeFv(e): { for k in integration-points: sim->elem->quad->epsdot(i,j) = ... sim->elem->quad->stressStrainRate(e,k,V) sim->elem->fv(...) += ... * sim->elem->quad->sig(i,j) } Fv(...) += sim->elem->fu(i) } Signature --------- From this it is clear that: * ``GooseFEM::Dynamics::Diagonal::Periodic`` requires the following minimal signature from ``GooseFEM::Dynamics::Diagonal::SmallStrain::Qaud4``: .. code-block:: cpp class Element { public: matrix M; // should have operator(i,j) column fu, fv; // should have operator(i) matrix xe, ue, ve; // should have operator(i,j) void computeM (size_t elem); // mass matrix <- quad->density void computeFu(size_t elem); // displacement dependent forces <- quad->stressStrain void computeFv(size_t elem); // displacement dependent forces <- quad->stressStrainRate void post (size_t elem); // post-process <- quad->stressStrain(Rate) } * ``GooseFEM::Dynamics::Diagonal::SmallStrain::Qaud4`` requires the minimal signature from ``Quadrature`` .. code-block:: cpp class Quadrature { public: tensor eps, epsdot, sig; // should have operator(i,j) 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); } diff --git a/docs/examples/DynamicsDiagonal/Periodic/laminate/eta=00.10__alpha=0.00/main.cpp b/docs/examples/DynamicsDiagonal/Periodic/laminate/eta=00.10__alpha=0.00/main.cpp index 86f29f8..2c0861c 100644 --- a/docs/examples/DynamicsDiagonal/Periodic/laminate/eta=00.10__alpha=0.00/main.cpp +++ b/docs/examples/DynamicsDiagonal/Periodic/laminate/eta=00.10__alpha=0.00/main.cpp @@ -1,230 +1,233 @@ #include #include #include #include #include // ------------------------------------------------------------------------------------------------- // introduce aliases using ColD = GooseFEM::ColD; using T2 = cppmat::cartesian2d::tensor2; using Mat = GooseMaterial::AmorphousSolid::LinearStrain::Elastic::Cartesian2d::Material; // ================================================================================================= class Material { public: // class variables // --------------- // strain(-rate), stress, mass density, background damping [mandatory] cppmat::matrix eps, epsdot, sig, rho, alpha; // mesh dimensions size_t nelem, nne=4, ndim=2, nip=4; // constitutive response [customize] std::vector material; // damping [customize] Mat rayleigh; // class functions // --------------- // constructor [customize] Material(size_t nelem, size_t nhard); // compute stress for a certain element "e" and integration point "k" [mandatory] void updated_eps (size_t e, size_t k); void updated_epsdot(size_t e, size_t k); // post process variables / functions [customize] // - potential energy cppmat::matrix Epot; // - compute potential energy void post(); }; // ------------------------------------------------------------------------------------------------- Material::Material(size_t _nelem, size_t _nhard) { // copy from input [customize] nelem = _nelem; // allocate symmetric tensors and scalars per element, per integration point [mandatory] eps .resize({nelem,nip,3}); epsdot.resize({nelem,nip,3}); sig .resize({nelem,nip,3}); rho .resize({nelem,nne }); alpha .resize({nelem,nne }); // set mass-density and background damping [customize] rho .setConstant(1.0); alpha.setConstant(0.0); // post variables [customize] Epot.resize({nelem,nip}); // constitutive response per element (per integration post) [customize] 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); } } // homogeneous damping [customize] rayleigh = Mat(10.,.1); } // ------------------------------------------------------------------------------------------------- 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 Epsdot, Eps, Sig; // point to correct position in matrix of tensors Epsdot.map(&epsdot(e,k)); Eps .map(&eps (e,k)); // compute stress Sig = material[e].stress(Eps) + rayleigh.stress(Epsdot); // copy to matrix of tensors std::copy(Sig.begin(), Sig.end(), &sig(e,k)); } // ------------------------------------------------------------------------------------------------- void Material::updated_epsdot(size_t e, size_t k) { // local views of the global arrays (speeds up indexing, and increases readability) cppmat::cartesian2d::tensor2s Eps, Epsdot, Sig; // point to correct position in matrix of tensors Epsdot.map(&epsdot(e,k)); Eps .map(&eps (e,k)); // compute stress Sig = material[e].stress(Eps) + rayleigh.stress(Epsdot); - + // copy to matrix of tensors std::copy(Sig.begin(), Sig.end(), &sig(e,k)); } // ------------------------------------------------------------------------------------------------- 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 ) { // point to correct position in matrix of tensors Eps.map(&eps(e,k)); // compute energy Epot(e,k) = material[e].energy(Eps); } } } } // ================================================================================================= // introduce aliases using Mesh = GooseFEM::Mesh::Quad4::Regular; using Element = GooseFEM::Dynamics::Diagonal::SmallStrain::Quad4; using Simulation = GooseFEM::Dynamics::Diagonal::Periodic; // ================================================================================================= int main() { // set simulation parameters double T = 60. ; // 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.nodesOrigin(); 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) ); + // process updates for displacement dependent variables + sim.updated_u(); + // 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->vol) * sim.elem->vol.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/DynamicsDiagonal/Periodic/laminate/eta=00.10__alpha=0.20/main.cpp b/docs/examples/DynamicsDiagonal/Periodic/laminate/eta=00.10__alpha=0.20/main.cpp index 5233e56..2b2e93e 100644 --- a/docs/examples/DynamicsDiagonal/Periodic/laminate/eta=00.10__alpha=0.20/main.cpp +++ b/docs/examples/DynamicsDiagonal/Periodic/laminate/eta=00.10__alpha=0.20/main.cpp @@ -1,230 +1,233 @@ #include #include #include #include #include // ------------------------------------------------------------------------------------------------- // introduce aliases using ColD = GooseFEM::ColD; using T2 = cppmat::cartesian2d::tensor2; using Mat = GooseMaterial::AmorphousSolid::LinearStrain::Elastic::Cartesian2d::Material; // ================================================================================================= class Material { public: // class variables // --------------- // strain(-rate), stress, mass density, background damping [mandatory] cppmat::matrix eps, epsdot, sig, rho, alpha; // mesh dimensions size_t nelem, nne=4, ndim=2, nip=4; // constitutive response [customize] std::vector material; // damping [customize] Mat rayleigh; // class functions // --------------- // constructor [customize] Material(size_t nelem, size_t nhard); // compute stress for a certain element "e" and integration point "k" [mandatory] void updated_eps (size_t e, size_t k); void updated_epsdot(size_t e, size_t k); // post process variables / functions [customize] // - potential energy cppmat::matrix Epot; // - compute potential energy void post(); }; // ------------------------------------------------------------------------------------------------- Material::Material(size_t _nelem, size_t _nhard) { // copy from input [customize] nelem = _nelem; // allocate symmetric tensors and scalars per element, per integration point [mandatory] eps .resize({nelem,nip,3}); epsdot.resize({nelem,nip,3}); sig .resize({nelem,nip,3}); rho .resize({nelem,nne }); alpha .resize({nelem,nne }); // set mass-density and background damping [customize] rho .setConstant(1.0); alpha.setConstant(0.2); // post variables [customize] Epot.resize({nelem,nip}); // constitutive response per element (per integration post) [customize] 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); } } // homogeneous damping [customize] rayleigh = Mat(10.,.1); } // ------------------------------------------------------------------------------------------------- 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 Epsdot, Eps, Sig; // point to correct position in matrix of tensors Epsdot.map(&epsdot(e,k)); Eps .map(&eps (e,k)); // compute stress Sig = material[e].stress(Eps) + rayleigh.stress(Epsdot); // copy to matrix of tensors std::copy(Sig.begin(), Sig.end(), &sig(e,k)); } // ------------------------------------------------------------------------------------------------- void Material::updated_epsdot(size_t e, size_t k) { // local views of the global arrays (speeds up indexing, and increases readability) cppmat::cartesian2d::tensor2s Eps, Epsdot, Sig; // point to correct position in matrix of tensors Epsdot.map(&epsdot(e,k)); Eps .map(&eps (e,k)); // compute stress Sig = material[e].stress(Eps) + rayleigh.stress(Epsdot); - + // copy to matrix of tensors std::copy(Sig.begin(), Sig.end(), &sig(e,k)); } // ------------------------------------------------------------------------------------------------- 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 ) { // point to correct position in matrix of tensors Eps.map(&eps(e,k)); // compute energy Epot(e,k) = material[e].energy(Eps); } } } } // ================================================================================================= // introduce aliases using Mesh = GooseFEM::Mesh::Quad4::Regular; using Element = GooseFEM::Dynamics::Diagonal::SmallStrain::Quad4; using Simulation = GooseFEM::Dynamics::Diagonal::Periodic; // ================================================================================================= int main() { // set simulation parameters double T = 60. ; // 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.nodesOrigin(); 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) ); + // process updates for displacement dependent variables + sim.updated_u(); + // 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->vol) * sim.elem->vol.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/DynamicsDiagonal/Periodic/laminate/no_damping/Verlet/main.cpp b/docs/examples/DynamicsDiagonal/Periodic/laminate/no_damping/Verlet/main.cpp index d8c3e2e..50ae9d1 100644 --- a/docs/examples/DynamicsDiagonal/Periodic/laminate/no_damping/Verlet/main.cpp +++ b/docs/examples/DynamicsDiagonal/Periodic/laminate/no_damping/Verlet/main.cpp @@ -1,211 +1,214 @@ #include #include #include #include #include // ------------------------------------------------------------------------------------------------- // introduce aliases using ColD = GooseFEM::ColD; using T2 = cppmat::cartesian2d::tensor2; using Mat = GooseMaterial::AmorphousSolid::LinearStrain::Elastic::Cartesian2d::Material; // ================================================================================================= class Material { public: // class variables // --------------- // strain(-rate), stress, mass density, background damping [mandatory] cppmat::matrix eps, epsdot, sig, rho, alpha; // mesh dimensions size_t nelem, nne=4, ndim=2, nip=4; // constitutive response [customize] std::vector material; // class functions // --------------- // constructor [customize] Material(size_t nelem, size_t nhard); // compute stress for a certain element "e" and integration point "k" [mandatory] void updated_eps (size_t e, size_t k); void updated_epsdot(size_t e, size_t k); // post process variables / functions [customize] // - potential energy cppmat::matrix Epot; // - compute potential energy void post(); }; // ------------------------------------------------------------------------------------------------- Material::Material(size_t _nelem, size_t _nhard) { // copy from input [customize] nelem = _nelem; // allocate symmetric tensors and scalars per element, per integration point [mandatory] eps .resize({nelem,nip,3}); epsdot.resize({nelem,nip,3}); sig .resize({nelem,nip,3}); rho .resize({nelem,nne }); alpha .resize({nelem,nne }); // set mass-density and background damping [customize] rho .setConstant(1.0); alpha.setConstant(0.0); // post variables [customize] Epot.resize({nelem,nip}); // constitutive response per element (per integration post) [customize] 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; // point to correct position in matrix of tensors Eps.map(&eps(e,k)); // compute stress Sig = material[e].stress(Eps); // copy to matrix of tensors std::copy(Sig.begin(), Sig.end(), &sig(e,k)); } // ------------------------------------------------------------------------------------------------- void Material::updated_epsdot(size_t e, size_t k) { } // ------------------------------------------------------------------------------------------------- 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 ) { // point to correct position in matrix of tensors Eps.map(&eps(e,k)); // compute energy Epot(e,k) = material[e].energy(Eps); } } } } // ================================================================================================= // introduce aliases using Mesh = GooseFEM::Mesh::Quad4::Regular; using Element = GooseFEM::Dynamics::Diagonal::SmallStrain::Quad4; using Simulation = GooseFEM::Dynamics::Diagonal::Periodic; // ================================================================================================= int main() { // set simulation parameters double T = 60. ; // 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.nodesOrigin(); 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) ); + // process updates for displacement dependent variables + sim.updated_u(); + // 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->vol) * sim.elem->vol.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/DynamicsDiagonal/Periodic/laminate/no_damping/velocityVerlet/main.cpp b/docs/examples/DynamicsDiagonal/Periodic/laminate/no_damping/velocityVerlet/main.cpp index 874d529..83a4daf 100644 --- a/docs/examples/DynamicsDiagonal/Periodic/laminate/no_damping/velocityVerlet/main.cpp +++ b/docs/examples/DynamicsDiagonal/Periodic/laminate/no_damping/velocityVerlet/main.cpp @@ -1,211 +1,214 @@ #include #include #include #include #include // ------------------------------------------------------------------------------------------------- // introduce aliases using ColD = GooseFEM::ColD; using T2 = cppmat::cartesian2d::tensor2; using Mat = GooseMaterial::AmorphousSolid::LinearStrain::Elastic::Cartesian2d::Material; // ================================================================================================= class Material { public: // class variables // --------------- // strain(-rate), stress, mass density, background damping [mandatory] cppmat::matrix eps, epsdot, sig, rho, alpha; // mesh dimensions size_t nelem, nne=4, ndim=2, nip=4; // constitutive response [customize] std::vector material; // class functions // --------------- // constructor [customize] Material(size_t nelem, size_t nhard); // compute stress for a certain element "e" and integration point "k" [mandatory] void updated_eps (size_t e, size_t k); void updated_epsdot(size_t e, size_t k); // post process variables / functions [customize] // - potential energy cppmat::matrix Epot; // - compute potential energy void post(); }; // ------------------------------------------------------------------------------------------------- Material::Material(size_t _nelem, size_t _nhard) { // copy from input [customize] nelem = _nelem; // allocate symmetric tensors and scalars per element, per integration point [mandatory] eps .resize({nelem,nip,3}); epsdot.resize({nelem,nip,3}); sig .resize({nelem,nip,3}); rho .resize({nelem,nne }); alpha .resize({nelem,nne }); // set mass-density and background damping [customize] rho .setConstant(1.0); alpha.setConstant(0.0); // post variables [customize] Epot.resize({nelem,nip}); // constitutive response per element (per integration post) [customize] 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; // point to correct position in matrix of tensors Eps.map(&eps(e,k)); // compute stress Sig = material[e].stress(Eps); // copy to matrix of tensors std::copy(Sig.begin(), Sig.end(), &sig(e,k)); } // ------------------------------------------------------------------------------------------------- void Material::updated_epsdot(size_t e, size_t k) { } // ------------------------------------------------------------------------------------------------- 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 ) { // point to correct position in matrix of tensors Eps.map(&eps(e,k)); // compute energy Epot(e,k) = material[e].energy(Eps); } } } } // ================================================================================================= // introduce aliases using Mesh = GooseFEM::Mesh::Quad4::Regular; using Element = GooseFEM::Dynamics::Diagonal::SmallStrain::Quad4; using Simulation = GooseFEM::Dynamics::Diagonal::Periodic; // ================================================================================================= int main() { // set simulation parameters double T = 60. ; // 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.nodesOrigin(); 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) ); + // process updates for displacement dependent variables + sim.updated_u(); + // 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->vol) * sim.elem->vol.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/DynamicsDiagonal/Periodic/laminate/overdamped/CMakeLists.txt b/docs/examples/DynamicsDiagonal/Periodic/laminate/overdamped/CMakeLists.txt new file mode 100644 index 0000000..0aafb9a --- /dev/null +++ b/docs/examples/DynamicsDiagonal/Periodic/laminate/overdamped/CMakeLists.txt @@ -0,0 +1,68 @@ +cmake_minimum_required(VERSION 3.1) + +# define a project name +project(example) + +# define empty list of libraries to link +set(PROJECT_LIBS "") + +# set optimization level +set(CMAKE_BUILD_TYPE Release) + +# set C++ standard +set(CMAKE_CXX_STANDARD 14) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) + +# add optimization options: switch on/off assertions +option(ASSERT "Switch on assertions" OFF) +if(NOT ASSERT) + add_definitions(-DNDEBUG) +endif() + +# set warnings on +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -pedantic") + +# OpenMP settings +option(OPENMP "Switch on parallelization" ON) +if(OPENMP) + if (NOT "${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang") + 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() +else() + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-unknown-pragmas") +endif() + +# load pkg-config +find_package(PkgConfig) + +# load HDF5 +find_package(HDF5 COMPONENTS CXX REQUIRED) +include_directories(${HDF5_INCLUDE_DIRS}) +set(PROJECT_LIBS ${PROJECT_LIBS} ${HDF5_LIBS} ${HDF5_LIBRARIES}) + +# load header-only modules using pkg-config +# - eigen3 +pkg_check_modules(EIGEN3 REQUIRED eigen3) +include_directories(SYSTEM ${EIGEN3_INCLUDE_DIRS}) +# - 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}) +# - 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} ${PROJECT_LIBS}) diff --git a/docs/examples/DynamicsDiagonal/Periodic/laminate/overdamped/example.svg b/docs/examples/DynamicsDiagonal/Periodic/laminate/overdamped/example.svg new file mode 100644 index 0000000..f273a98 --- /dev/null +++ b/docs/examples/DynamicsDiagonal/Periodic/laminate/overdamped/example.svg @@ -0,0 +1,11053 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/examples/DynamicsDiagonal/Periodic/laminate/no_damping/Verlet/main.cpp b/docs/examples/DynamicsDiagonal/Periodic/laminate/overdamped/main.cpp similarity index 87% copy from docs/examples/DynamicsDiagonal/Periodic/laminate/no_damping/Verlet/main.cpp copy to docs/examples/DynamicsDiagonal/Periodic/laminate/overdamped/main.cpp index d8c3e2e..9c13a77 100644 --- a/docs/examples/DynamicsDiagonal/Periodic/laminate/no_damping/Verlet/main.cpp +++ b/docs/examples/DynamicsDiagonal/Periodic/laminate/overdamped/main.cpp @@ -1,211 +1,211 @@ #include #include #include #include #include // ------------------------------------------------------------------------------------------------- // introduce aliases using ColD = GooseFEM::ColD; using T2 = cppmat::cartesian2d::tensor2; using Mat = GooseMaterial::AmorphousSolid::LinearStrain::Elastic::Cartesian2d::Material; // ================================================================================================= class Material { public: // class variables // --------------- - // strain(-rate), stress, mass density, background damping [mandatory] - cppmat::matrix eps, epsdot, sig, rho, alpha; + // strain(-rate), stress, background damping [mandatory] + cppmat::matrix eps, sig, alpha; // mesh dimensions size_t nelem, nne=4, ndim=2, nip=4; // constitutive response [customize] std::vector material; // class functions // --------------- // constructor [customize] Material(size_t nelem, size_t nhard); // compute stress for a certain element "e" and integration point "k" [mandatory] - void updated_eps (size_t e, size_t k); - void updated_epsdot(size_t e, size_t k); + void updated_eps(size_t e, size_t k); // post process variables / functions [customize] // - potential energy cppmat::matrix Epot; // - compute potential energy void post(); }; // ------------------------------------------------------------------------------------------------- Material::Material(size_t _nelem, size_t _nhard) { // copy from input [customize] nelem = _nelem; // allocate symmetric tensors and scalars per element, per integration point [mandatory] eps .resize({nelem,nip,3}); - epsdot.resize({nelem,nip,3}); sig .resize({nelem,nip,3}); - rho .resize({nelem,nne }); alpha .resize({nelem,nne }); // set mass-density and background damping [customize] - rho .setConstant(1.0); - alpha.setConstant(0.0); + alpha.setConstant(0.2); // post variables [customize] Epot.resize({nelem,nip}); // constitutive response per element (per integration post) [customize] 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; // point to correct position in matrix of tensors Eps.map(&eps(e,k)); // compute stress Sig = material[e].stress(Eps); // copy to matrix of tensors std::copy(Sig.begin(), Sig.end(), &sig(e,k)); } // ------------------------------------------------------------------------------------------------- -void Material::updated_epsdot(size_t e, size_t k) -{ -} - -// ------------------------------------------------------------------------------------------------- - 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 ) { // point to correct position in matrix of tensors Eps.map(&eps(e,k)); // compute energy Epot(e,k) = material[e].energy(Eps); } } } } // ================================================================================================= // introduce aliases using Mesh = GooseFEM::Mesh::Quad4::Regular; -using Element = GooseFEM::Dynamics::Diagonal::SmallStrain::Quad4; -using Simulation = GooseFEM::Dynamics::Diagonal::Periodic; +using Element = GooseFEM::OverdampedDynamics::Diagonal::SmallStrain::Quad4; +using Simulation = GooseFEM::OverdampedDynamics::Diagonal::Periodic; // ================================================================================================= int main() { // set simulation parameters double T = 60. ; // total time - double dt = 1.e-2; // time increment + double dt = 1.e-3; // time increment size_t nx = 40 ; // number of elements in both directions double gamma = .05 ; // displacement step + // damping and mass-density (to allow conversion) + double alpha = 0.2; + double rho = 1.0; + // class which provides the mesh Mesh mesh(nx,nx,1.); // extract information size_t nodeOrigin = mesh.nodesOrigin(); 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) ); + // process updates for displacement dependent variables + sim.updated_u(); + // 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(); + // mass matrix + ColD M = rho/alpha * sim.D; + // loop over increments for ( size_t inc = 0 ; inc < static_cast(Epot.size()) ; ++inc ) { // - compute increment - sim.Verlet(); + sim.forwardEuler(); // - 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.); + Ekin(inc) += .5 * 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->vol) * sim.elem->vol.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/DynamicsDiagonal/Periodic/laminate/overdamped/plot.py b/docs/examples/DynamicsDiagonal/Periodic/laminate/overdamped/plot.py new file mode 100644 index 0000000..c31473d --- /dev/null +++ b/docs/examples/DynamicsDiagonal/Periodic/laminate/overdamped/plot.py @@ -0,0 +1,57 @@ + +import sys,os,re +import h5py +import numpy as np +import matplotlib.pyplot as plt +import goosempl as gplt +import matplotlib as mpl + +plt.style.use(['goose','goose-latex']) + +# -------------------------------------------------------------------------------------------------- + +name = 'example.hdf5' +f = h5py.File(name,'r') + +fig = plt.figure(figsize=(20,10)) +fig.set_tight_layout(True) + +# -------------------------------------------------------------------------------------------------- + +ax = fig.add_subplot(1,2,1) + +t = f['/global/t' ][:] +Epot = f['/global/Epot'][:] +Ekin = f['/global/Ekin'][:] + +ax.plot(t,Epot+Ekin,color='k',label=r'$E_\mathrm{pot} + E_\mathrm{kin}$') +ax.plot(t,Epot ,color='r',label=r'$E_\mathrm{pot}$') +ax.plot(t,Ekin ,color='b',label=r'$E_\mathrm{kin}$') + +plt.legend(loc='lower right') + +plt.xlabel(r'$t$') +plt.ylabel(r'$E$') + +# -------------------------------------------------------------------------------------------------- + +ax = fig.add_subplot(1,2,2) + +coor = f['/mesh/coor'][:] +conn = f['/mesh/conn'][:] +u = f['/mesh/disp'][:] + +im = gplt.patch(coor=coor+u,conn=conn) + +Lx = np.max(coor[:,0])-np.min(coor[:,0]) +Ly = np.max(coor[:,1])-np.min(coor[:,1]) + +plt.xlim([-0.1*Lx,1.3*Lx]) +plt.ylim([-0.1*Ly,1.3*Ly]) + +ax.set_aspect('equal') + +plt.savefig(re.sub(r'(.*)(\.hdf5)',r'\1.svg',name)) +plt.show() + + diff --git a/docs/theory/.gitignore b/docs/theory/.gitignore new file mode 100644 index 0000000..b6418e5 --- /dev/null +++ b/docs/theory/.gitignore @@ -0,0 +1,221 @@ +## Core latex/pdflatex auxiliary files: +*.aux +*.lof +*.log +*.lot +*.fls +*.out +*.toc +*.fmt +*.fot +*.cb +*.cb2 + +## Intermediate documents: +*.dvi +*.xdv +*-converted-to.* +# these rules might exclude image files for figures etc. +# *.ps +# *.eps +# *.pdf + +## Generated if empty string is given at "Please type another file name for output:" +.pdf + +## Bibliography auxiliary files (bibtex/biblatex/biber): +*.bbl +*.bcf +*.blg +*-blx.aux +*-blx.bib +*.run.xml + +## Build tool auxiliary files: +*.fdb_latexmk +*.synctex +*.synctex(busy) +*.synctex.gz +*.synctex.gz(busy) +*.pdfsync + +## Auxiliary and intermediate files from other packages: +# algorithms +*.alg +*.loa + +# achemso +acs-*.bib + +# amsthm +*.thm + +# beamer +*.nav +*.pre +*.snm +*.vrb + +# changes +*.soc + +# cprotect +*.cpt + +# elsarticle (documentclass of Elsevier journals) +*.spl + +# endnotes +*.ent + +# fixme +*.lox + +# feynmf/feynmp +*.mf +*.mp +*.t[1-9] +*.t[1-9][0-9] +*.tfm + +#(r)(e)ledmac/(r)(e)ledpar +*.end +*.?end +*.[1-9] +*.[1-9][0-9] +*.[1-9][0-9][0-9] +*.[1-9]R +*.[1-9][0-9]R +*.[1-9][0-9][0-9]R +*.eledsec[1-9] +*.eledsec[1-9]R +*.eledsec[1-9][0-9] +*.eledsec[1-9][0-9]R +*.eledsec[1-9][0-9][0-9] +*.eledsec[1-9][0-9][0-9]R + +# glossaries +*.acn +*.acr +*.glg +*.glo +*.gls +*.glsdefs + +# gnuplottex +*-gnuplottex-* + +# gregoriotex +*.gaux +*.gtex + +# hyperref +*.brf + +# knitr +*-concordance.tex +# TODO Comment the next line if you want to keep your tikz graphics files +*.tikz +*-tikzDictionary + +# listings +*.lol + +# makeidx +*.idx +*.ilg +*.ind +*.ist + +# minitoc +*.maf +*.mlf +*.mlt +*.mtc[0-9]* +*.slf[0-9]* +*.slt[0-9]* +*.stc[0-9]* + +# minted +_minted* +*.pyg + +# morewrites +*.mw + +# nomencl +*.nlo + +# pax +*.pax + +# pdfpcnotes +*.pdfpc + +# sagetex +*.sagetex.sage +*.sagetex.py +*.sagetex.scmd + +# scrwfile +*.wrt + +# sympy +*.sout +*.sympy +sympy-plots-for-*.tex/ + +# pdfcomment +*.upa +*.upb + +# pythontex +*.pytxcode +pythontex-files-*/ + +# thmtools +*.loe + +# TikZ & PGF +*.dpth +*.md5 +*.auxlock + +# todonotes +*.tdo + +# easy-todo +*.lod + +# xindy +*.xdy + +# xypic precompiled matrices +*.xyc + +# endfloat +*.ttt +*.fff + +# Latexian +TSWLatexianTemp* + +## Editors: +# WinEdt +*.bak +*.sav + +# Texpad +.texpadtmp + +# Kile +*.backup + +# KBibTeX +*~[0-9]* + +# auto folder when using emacs and auctex +/auto/* + +# expex forward references with \gathertags +*-tags.tex diff --git a/docs/theory/figures/isoparametric-transform.pdf b/docs/theory/figures/isoparametric-transform.pdf new file mode 100644 index 0000000..1ebc1ad Binary files /dev/null and b/docs/theory/figures/isoparametric-transform.pdf differ diff --git a/docs/theory/figures/isoparametric-transform.svg b/docs/theory/figures/isoparametric-transform.svg new file mode 100644 index 0000000..a66a1d2 --- /dev/null +++ b/docs/theory/figures/isoparametric-transform.svg @@ -0,0 +1,684 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/theory/figures/newton-raphson.pdf b/docs/theory/figures/newton-raphson.pdf new file mode 100644 index 0000000..afa1908 Binary files /dev/null and b/docs/theory/figures/newton-raphson.pdf differ diff --git a/docs/theory/figures/newton-raphson.svg b/docs/theory/figures/newton-raphson.svg new file mode 100644 index 0000000..1f141da --- /dev/null +++ b/docs/theory/figures/newton-raphson.svg @@ -0,0 +1,549 @@ + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/theory/figures/problem-discretized.pdf b/docs/theory/figures/problem-discretized.pdf new file mode 100644 index 0000000..e72faac Binary files /dev/null and b/docs/theory/figures/problem-discretized.pdf differ diff --git a/docs/theory/figures/problem-discretized.svg b/docs/theory/figures/problem-discretized.svg new file mode 100644 index 0000000..622e6e3 --- /dev/null +++ b/docs/theory/figures/problem-discretized.svg @@ -0,0 +1,1513 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/theory/figures/problem.pdf b/docs/theory/figures/problem.pdf new file mode 100644 index 0000000..4238181 Binary files /dev/null and b/docs/theory/figures/problem.pdf differ diff --git a/docs/theory/figures/problem.svg b/docs/theory/figures/problem.svg new file mode 100644 index 0000000..50bd0e1 --- /dev/null +++ b/docs/theory/figures/problem.svg @@ -0,0 +1,657 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/theory/figures/shape-functions-1d-element-0.pdf b/docs/theory/figures/shape-functions-1d-element-0.pdf new file mode 100644 index 0000000..627f2a2 Binary files /dev/null and b/docs/theory/figures/shape-functions-1d-element-0.pdf differ diff --git a/docs/theory/figures/shape-functions-1d-element-0.svg b/docs/theory/figures/shape-functions-1d-element-0.svg new file mode 100644 index 0000000..0a0a4fe --- /dev/null +++ b/docs/theory/figures/shape-functions-1d-element-0.svg @@ -0,0 +1,520 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + 0 + 1 + 2 + 3 + 4 + 0 + 1 + 2 + 3 + + + + + + + + diff --git a/docs/theory/figures/shape-functions-1d-element-1.pdf b/docs/theory/figures/shape-functions-1d-element-1.pdf new file mode 100644 index 0000000..10bd96a Binary files /dev/null and b/docs/theory/figures/shape-functions-1d-element-1.pdf differ diff --git a/docs/theory/figures/shape-functions-1d-element-1.svg b/docs/theory/figures/shape-functions-1d-element-1.svg new file mode 100644 index 0000000..a4f07d1 --- /dev/null +++ b/docs/theory/figures/shape-functions-1d-element-1.svg @@ -0,0 +1,520 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + 0 + 1 + 2 + 3 + 4 + 0 + 1 + 2 + 3 + + + + + + + + + + diff --git a/docs/theory/figures/shape-functions-1d-element-2.pdf b/docs/theory/figures/shape-functions-1d-element-2.pdf new file mode 100644 index 0000000..7eba2bf Binary files /dev/null and b/docs/theory/figures/shape-functions-1d-element-2.pdf differ diff --git a/docs/theory/figures/shape-functions-1d-element-2.svg b/docs/theory/figures/shape-functions-1d-element-2.svg new file mode 100644 index 0000000..3185ea7 --- /dev/null +++ b/docs/theory/figures/shape-functions-1d-element-2.svg @@ -0,0 +1,536 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + 0 + 1 + 2 + 3 + 4 + 0 + 1 + 2 + 3 + + + + + + + + + + diff --git a/docs/theory/figures/shape-functions-1d-element-3.pdf b/docs/theory/figures/shape-functions-1d-element-3.pdf new file mode 100644 index 0000000..bf9a0b9 Binary files /dev/null and b/docs/theory/figures/shape-functions-1d-element-3.pdf differ diff --git a/docs/theory/figures/shape-functions-1d-element-3.svg b/docs/theory/figures/shape-functions-1d-element-3.svg new file mode 100644 index 0000000..42fa1fd --- /dev/null +++ b/docs/theory/figures/shape-functions-1d-element-3.svg @@ -0,0 +1,520 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + 0 + 1 + 2 + 3 + 4 + 0 + 1 + 2 + 3 + + + + + + + + + + diff --git a/docs/theory/figures/shape-functions-1d-node-2.pdf b/docs/theory/figures/shape-functions-1d-node-2.pdf new file mode 100644 index 0000000..ff31870 Binary files /dev/null and b/docs/theory/figures/shape-functions-1d-node-2.pdf differ diff --git a/docs/theory/figures/shape-functions-1d-node-2.svg b/docs/theory/figures/shape-functions-1d-node-2.svg new file mode 100644 index 0000000..459ce95 --- /dev/null +++ b/docs/theory/figures/shape-functions-1d-node-2.svg @@ -0,0 +1,473 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + 0 + 1 + 2 + 3 + 4 + 0 + 1 + 2 + 3 + + + + + + + + + diff --git a/docs/theory/figures/shape-functions-1d.pdf b/docs/theory/figures/shape-functions-1d.pdf new file mode 100644 index 0000000..c24a1fb Binary files /dev/null and b/docs/theory/figures/shape-functions-1d.pdf differ diff --git a/docs/theory/figures/shape-functions-1d.svg b/docs/theory/figures/shape-functions-1d.svg new file mode 100644 index 0000000..6f0766d --- /dev/null +++ b/docs/theory/figures/shape-functions-1d.svg @@ -0,0 +1,629 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + 0 + 1 + 2 + 3 + 4 + 0 + 1 + 2 + 3 + + + + + + + + + diff --git a/docs/theory/goose-article.cls b/docs/theory/goose-article.cls new file mode 100644 index 0000000..c3484e4 --- /dev/null +++ b/docs/theory/goose-article.cls @@ -0,0 +1,290 @@ +% ------------------------------------------------------------------------------ +% @LaTeX-class-file{ +% filename = "goose-article.cls", +% version = "0.3", +% date = "16 November 2016", +% codetable = "ISO/ASCII", +% keywords = "LaTeX, goose-article", +% supported = "http://tdegeus.github.io/GooseLaTeX", +% docstring = "A customized article-class." +% ------------------------------------------------------------------------------ + +% ============================================================================== +% Basic settings +% ============================================================================== + +\NeedsTeXFormat{LaTeX2e} +\ProvidesClass{goose-article}[2016/11/16 v0.3 customized article class] + +% load default article class +\LoadClass[a4paper,fleqn,twoside]{article} + +% ============================================================================== +% Options +% ============================================================================== + +% bibliography [namecite] +% ----------------------- + +\newif\if@namecite +\let\if@namecite\iffalse +\DeclareOption{namecite}{\let\if@namecite\iftrue} + +% Review: [narrow,doublespacing] +% ------------------------------ + +\newif\if@narrow +\let\if@narrow\iffalse +\DeclareOption{narrow}{\let\if@narrow\iftrue} + +\newif\if@doublespacing +\let\if@doublespacing\iffalse +\DeclareOption{doublespacing}{\let\if@doublespacing\iftrue} + +% Font: [times,garamond,verdana] +% ------------------------------ + +\newif\if@font +\let\if@font\iffalse + +\newif\if@times +\let\if@times\iffalse +\DeclareOption{times}{\let\if@times\iftrue\let\if@font\iftrue} + +\newif\if@garamond +\let\if@garamond\iffalse +\DeclareOption{garamond}{\let\if@garamond\iftrue\let\if@font\iftrue} + +\newif\if@verdana +\let\if@verdana\iffalse +\DeclareOption{verdana}{\let\if@verdana\iftrue\let\if@font\iftrue} + +% Process options +% --------------- + +\ProcessOptions\relax + +% ============================================================================== +% Packages +% ============================================================================== + +\RequirePackage[left=25mm,right=25mm,top=30mm,bottom=30mm]{geometry} +\if@narrow + \newgeometry{left=40mm,right=40mm,top=30mm,bottom=30mm} +\fi + +\RequirePackage{fix-cm} +\RequirePackage{graphicx} +\RequirePackage[labelformat=simple]{subcaption} +\RequirePackage{amssymb,amsfonts,bm,cancel} +\RequirePackage{enumerate} +\RequirePackage[fleqn]{amsmath} +\RequirePackage[font={small},labelfont=bf,labelsep=period]{caption} +\RequirePackage[bf]{titlesec} +\RequirePackage[dvipsnames]{xcolor} +\RequirePackage{authblk} +\RequirePackage{hyperref} + +\renewcommand\thesubfigure{(\alph{subfigure})} + +% ============================================================================== +% Citation style +% ============================================================================== + +\if@namecite + \RequirePackage{natbib} + \let\oldbibliography\bibliography + \renewcommand{\bibliography}[1]{% + \bibliographystyle{apalike} + \setlength{\bibsep}{3pt plus 0.3ex} + \def\bibfont{\scriptsize} + \oldbibliography{#1} + } +\else + \RequirePackage[square,sort&compress,numbers,comma]{natbib} + \let\oldbibliography\bibliography + \renewcommand{\bibliography}[1]{% + \bibliographystyle{unsrtnat} + \setlength{\bibsep}{3pt plus 0.3ex} + \def\bibfont{\scriptsize} + \oldbibliography{#1} + } +\fi + +% define DOI command +\newcommand*{\doi}[1]{\sloppy\href{http://dx.doi.org/#1}{doi:~#1}} +\newcommand*{\eprint}[1]{\sloppy\href{http://arxiv.org/abs/#1}{arXiv:~#1}} + +% ============================================================================== +% Font +% ============================================================================== + +\if@font + \RequirePackage{mathspec} + \RequirePackage{fontspec} +\fi + +\if@garamond + \AtBeginDocument{% + \setmathfont(Digits)[Numbers=OldStyle]{Garamond}% + \setromanfont[% + BoldFont={Garamond Bold},% + ItalicFont={Garamond Italic},% + Mapping=tex-text% + ]{Garamond}% + }% +\fi + +\if@verdana + \AtBeginDocument{% + \defaultfontfeatures{Mapping=tex-text} + \setmainfont[Mapping=tex-text]{Verdana} + \setsansfont{Verdana} + \renewcommand*{\familydefault}{\sfdefault} + \setmathfont(Greek,Digits,Latin){Verdana} + \setmathrm{Verdana} + }% +\fi + +\if@times + \AtBeginDocument{% + \setmathfont(Digits)[Numbers=OldStyle]{Times New Roman}% + \setromanfont[% + BoldFont={Times New Roman Bold},% + ItalicFont={Times New Roman Italic},% + Mapping=tex-text% + ]{Times New Roman}% + }% +\fi + +% ============================================================================== +% Spacing +% ============================================================================== + +\if@doublespacing + \AtBeginDocument{% + \usepackage{setspace} + \doublespacing + }% +\fi + +% ============================================================================== +% Hyperref +% ============================================================================== + +\AtBeginDocument{% + \hypersetup{% + pdftitle=\@title, + citecolor=Blue, + filecolor=black, + linkcolor=Blue, + urlcolor=Blue, + breaklinks, + hidelinks, + bookmarksopen=true, + }% +} + +% ============================================================================== +% Headers +% ============================================================================== + +% {command}{format}{label}{sep}{before-code}{after-code} +% - section +\titleformat + {\section} + {\normalfont\large\bfseries} + {\thesection} + {1em} + {} + {} +% - subsection +\titleformat + {\subsection} + {\normalfont\normalsize\bfseries} + {\thesubsection} + {1em} + {} + {} +% - paragraph +\titleformat + {\paragraph} + {\normalfont\normalsize\bfseries} + {} + {} + {} + +% spacing {left}{top}{bottom}, of: +% - section +\titlespacing{\section} + {0pt}{1eM}{0.8em} +% - subsection +\titlespacing{\subsection} + {0pt}{0.8em}{0.4em} + +% no paragraph indent after header +\AtBeginDocument{% + \makeatletter + \let\orig@afterheading\@afterheading + \def\@afterheading{% + \@afterindentfalse + \orig@afterheading} + \makeatother +} + +% ============================================================================== +% Title +% ============================================================================== + +% title, authors, affiliation layout +% - contact +\newcommand{\contact}[1]{\date{{#1}}} +% - basic settings +\newcommand{\nl}{\vspace*{-.3em}\protect\\ \fontsize{9}{10.8}\itshape} +\renewcommand\Authfont{\fontsize{11}{14.4}\selectfont} +\renewcommand\Affilfont{\fontsize{9}{10.8}\itshape} +\renewcommand*{\Authsep}{, } +\renewcommand*{\Authand}{, } +\renewcommand*{\Authands}{, } +% - actual layout +\renewcommand{\maketitle}{% + \newpage + \null + \vskip 2em% + \begin{center}% + \let \footnote \thanks + {\Large\bfseries \@title \par}% + \vskip 1.5em% + {\normalsize + \lineskip .5em% + \begin{tabular}[t]{c}% + \@author + \end{tabular}\par}% + \vskip 0.5em% + {\small \@date}% + \end{center}% + \par + \vskip 1.5em + \thispagestyle{fancy} +} + +% ============================================================================== +% Header and footer +% ============================================================================== + +\RequirePackage{fancyhdr} +\setlength{\headheight}{16pt} +\pagestyle{fancy} +\fancyhead{} +\fancyfoot{} +\renewcommand{\headrulewidth}{0pt} +\renewcommand{\footrulewidth}{0pt} +\fancyhead[RO,LE]{\scriptsize \thepage} +\newcommand{\header}[1]{\fancyhead[LO,RE]{\scriptsize\itshape#1}} + +% ============================================================================== +% Abstract and keywords +% ============================================================================== + +\renewenvironment{abstract}{\section*{Abstract}\noindent\ignorespaces}{\ignorespacesafterend} +\newcommand{\keywords}[1]{\vspace*{0.5eM}\noindent\textbf{Keywords: }#1} diff --git a/docs/theory/readme.pdf b/docs/theory/readme.pdf new file mode 100644 index 0000000..8a32390 Binary files /dev/null and b/docs/theory/readme.pdf differ diff --git a/docs/theory/readme.rst b/docs/theory/readme.rst index 37907cd..0532cb5 100644 --- a/docs/theory/readme.rst +++ b/docs/theory/readme.rst @@ -1,938 +1,1144 @@ ********************* Finite Element Method ********************* .. contents:: **Contents** :local: :depth: 2 :backlinks: top In the sequel the theory of the Finite Element is discussed in a compact way. This discussion is by no means comprehensive. Therefore one is invited to dive in more complete textbooks. The key contribution of this reader is that it is supported by many examples that can be easily extended and customized into efficient, production-ready code. To this end, the examples are in C++14, and are specifically written such that they are mostly 'what you see is what you get'. The entire structure is in the main-file, and not hidden somewhere in a library. To simplify your life we do use several libraries, each of which, however, only with a dedicated task, which can be understood, used, and checked independently of the Finite Element Method or any specific application. More specifically we use: * `GooseMaterial `_ Provides the constitutive response (and optionally the constitutive tangent) of several materials. * `cppmat `_ Provides tensor classes and operations. (The number of tensor operations are limited in the main program, and even non-standard, but this library is crucial to compute the material response implemented in `GooseMaterial `_.) * `Eigen3 `_ A linear algebra library. As you will notice, Eigen plays an important role in GooseFEM, and glues everything together since in the end the Finite Element Method is just a way to cast a problem into a set linear or linearized equations. Most of the efficiency of the final program will depend on the efficiency of the implementation of the linear algebra. In several examples we will simplify the structure by using dense matrices together with a simple solver which solves the resulting linear system. In reality one should always use sparse matrices combined with a more efficient solver. As you will notice, many examples need only few changes to be transformed in a production code. .. note:: **Compilation** Unless otherwise mentioned, the examples can be compiled as follows. Provided that ``pkg-config`` is set-up correctly one can use .. code-block:: bash clang++ `pkg-config --cflags Eigen3 cppmat GooseMaterial GooseFEM` -std=c++14 -o example example_name.cpp (whereby ``clang++`` can be replaced by for example ``g++``). If one does not want to use ``pkg-config``, one has to specify ``-I/path/to/library`` for each of the libraries. For further development it is strongly advised to include the options ``-Wpedantic -Wall`` to get on top of mistakes. Once the code is ready, one should compile with optimization (``-O3``) and without assertions (``-DNDEBUG``). The `Eigen3 documentation `_ further recommends the option ``-march=native`` to enable vectorization optimized for you architecture. Statics ======= The conceptual idea ------------------- We begin our discussion by considering a static, solid mechanics, problem. Loosely speaking the goal is to find a deformation map, :math:`\vec{x} = \varphi(\vec{X},t)`, that maps a body :math:`\Omega_0` to a deformed state :math:`\Omega` that satisfies equilibrium and the boundary conditions applied on :math:`\Gamma`. .. image:: problem.svg :width: 550px :align: center As is the case in the illustration, the body can be subjected to two kinds of boundary conditions: * *Essential* or *Dirichlet* boundary conditions on :math:`\Gamma_p`, whereby the displacements are prescribed. * *Natural* or *Neumann* boundary conditions on :math:`\Gamma_u`, whereby the tractions are prescribed. (Whereby traction-free is also perfectly acceptable.) In practice, we are not explicitly looking for the map :math:`\vec{x} = \varphi(\vec{X},t)`, but for the deformation gradient :math:`\bm{F}`, or in fact for a displacement field :math:`\vec{u} = \vec{x} - \vec{X}`. To make things a bit more explicit, the deformation gradient is defined as follows: .. math:: \vec{x} = \bm{F} \cdot \vec{X} hence .. math:: \bm{F} = \frac{\partial \varphi}{\partial \vec{X}} = \big( \vec{\nabla}_0 \, \vec{x} \big)^T = \bm{I} + \big( \vec{\nabla}_0 \, \vec{u} \big)^T Momentum balance ---------------- We start from the linear momentum balance: .. math:: \vec{\nabla} \cdot \bm{\sigma}(\vec{x}) = \vec{0} \qquad \vec{x} \in \Omega where :math:`\bm{\sigma}` is the Cauchy stress which depends on the new position :math:`\vec{x}` and thus on the displacement :math:`\vec{u}`. It has been assumed that all actions are instantaneous (no inertia) and, for simplicity, that there are no body forces. Loosely speaking the interpretation of this equation is that *the sum of all forces vanishes* everywhere in the domain :math:`\Omega` .. note:: The following nomenclature has been used .. math:: \vec{\nabla} \cdot \bm{\sigma} = \frac{ \partial \sigma_{ij} }{ \partial x_i } The crux of the Finite Element Method is that this non-linear differential equation is solved in a weak sense. I.e. .. math:: \int\limits_\Omega \vec{\phi}(\vec{X}) \cdot \big[\, \vec{\nabla} \cdot \bm{\sigma}(\vec{x}) \,\big] \; \mathrm{d}\Omega = 0 \qquad \forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d where :math:`\vec{\phi}` are test functions. For reasons that become obvious below, we apply integration by parts, which results in .. math:: \int\limits_\Omega \big[\, \vec{\nabla} \vec{\phi}(\vec{X}) \,\big] : \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega = \int\limits_\Omega \vec{\nabla} \cdot \big[\, \vec{\phi}(\vec{X}) \cdot \bm{\sigma}(\vec{x}) \,\big] \; \mathrm{d}\Omega \qquad \forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d .. note:: Use has been made of the following chain rule .. math:: \vec{\nabla} \cdot \big[\, \vec{\phi} \cdot \bm{\sigma}^T \,\big] = \big[\, \vec{\nabla} \vec{\phi} \,\big] : \bm{\sigma}^T + \vec{\phi} \cdot \big[\, \vec{\nabla} \cdot \bm{\sigma} \,\big] together with the symmetry of the Cauchy stress .. math:: \bm{\sigma} = \bm{\sigma}^T and the following nomenclature: .. math:: C = \bm{A} : \bm{B} = A_{ij} B_{ji} The right-hand side of this equation can be reduced to an area integral by employing Gauss's divergence theorem. The result reads .. math:: \int\limits_\Omega \big[\, \vec{\nabla} \vec{\phi}(\vec{X}) \,\big] : \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega = \int\limits_\Gamma \vec{\phi}(\vec{X}) \cdot \underbrace{ \vec{n}(\vec{x}) \cdot \bm{\sigma}(\vec{x}) }_{ \vec{t}(\vec{x}) } \; \mathrm{d}\Gamma \qquad \forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d .. note:: Gauss's divergence theorem states that .. math:: \int\limits_\Omega \vec{\nabla} \cdot \vec{a}(\vec{x}) \; \mathrm{d}\Omega = \int\limits_\Gamma \vec{n}(\vec{x}) \cdot \vec{a}(\vec{x}) \; \mathrm{d}\Gamma where :math:`\vec{n}` is the normal along the surface :math:`\Gamma`. Discretization -------------- The problem is now discretized using :math:`n` *nodes* that are connected through *elements*, which define the discretized domain :math:`\Omega^h_0`. `Shape functions`_ :math:`N_i(\vec{X})` are used to extrapolate the nodal quantities throughout the domain :math:`\Omega^h_0` (and :math:`\Omega^h`), as follows: .. math:: \vec{x}(\vec{X},t) \approx \vec{x}^h(\vec{X},t) = \sum_{i=1}^{n} N_i (\vec{X}) \; \vec{x}_i (t) = \underline{N}^\mathsf{T} (\vec{X}) \; \underline{\vec{x}} (t) Following standard Galerkin .. math:: \vec{\phi}(\vec{X}) \approx \vec{\phi}^h(\vec{X}) = \underline{N}^\mathsf{T} (\vec{X}) \; \underline{\vec{\phi}} .. note:: Applied to our problem sketch, a discretization might look like this. The nodes are clearly marked as circles. The lines connecting the nodes clearly mark the elements which are in this case three-node triangles (Tri3 in GooseFEM) .. image:: problem-discretized.svg :width: 300px :align: center Applied to the balance equation we obtain .. math:: \underline{\vec{\phi}}^\mathsf{T} \cdot \int\limits_{\Omega^h} \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big] \cdot \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega = \underline{\vec{\phi}}^\mathsf{T} \cdot \int\limits_{\Gamma^h} \underline{N}(\vec{X}) \cdot \vec{t}(\vec{x}) \; \mathrm{d}\Gamma \qquad \forall \; \underline{\vec{\phi}} \in \mathbb{R}^d_n from which the dependency on :math:`\underline{\vec{\phi}}` can be dropped: .. math:: \int\limits_{\Omega^h} \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big] \cdot \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega = \int\limits_{\Gamma^h} \underline{N}(\vec{X}) \cdot \vec{t}(\vec{x}) \; \mathrm{d}\Gamma This corresponds to (non-linear) set of nodal balance equations: .. math:: \underline{\vec{f}}(\vec{x}) = \underline{\vec{t}}(\vec{x}) with: * *Internal forces* .. math:: \underline{\vec{f}}(\vec{x}) = \int\limits_{\Omega^h} \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big] \cdot \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega * *Boundary tractions* .. math:: \underline{\vec{t}}(\vec{x}) = \int\limits_{\Gamma^h} \underline{N}(\vec{X}) \cdot \vec{t}(\vec{x}) \; \mathrm{d}\Gamma which is zero in the interior of the domain, i.e. in :math:`\Omega^h \bigcap \Gamma^h`, while they can be zero or non-zero in :math:`\Gamma^h` depending on the problem details. Iterative solution -- small strain ---------------------------------- A commonly used strategy to solve the non-linear system is the iterative Newton-Raphson scheme (see inset below). The idea is thereby to formulate an initial guess for the solution, determine possible residual forces, and use these forces to come to a better guess for the solution. This is continued until the solution has been found, i.e. when the residual vanishes. This solution technique is discussed here in the context of small deformations, while it is later generalized. Assuming the deformations to be small allows us to assume that :math:`\Omega = \Omega_0`, and thus that :math:`\nabla = \nabla_0`. Also we define a strain tensor .. math:: \bm{\varepsilon} = \tfrac{1}{2} \left[ \nabla_0 \vec{u} + \big[\, \nabla_0 \vec{u} \,\big]^T \right] = \mathbb{I}_s : \big[\, \nabla_0 \vec{u} \,\big] and use some non-linear relationship between it and the stress .. math:: \bm{\sigma} = \bm{\sigma} \big( \bm{\varepsilon} \big) To simplify our discussion we assume the boundary tractions to be some known constant. Our nodal equilibrium equations now read .. math:: \underline{\vec{r}}(\vec{x}) = \underline{\vec{t}} - \underline{\vec{f}}(\vec{x}) = \underline{\vec{0}} with .. math:: \underline{\vec{f}}(\vec{x}) = \int\limits_{\Omega^h_0} \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big] \cdot \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega To come to an iterative solution, we linearize as this point. This results in .. math:: \int\limits_{\Omega^h_0} \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big] \cdot \mathbb{K}\big(\vec{x}_{(i)}\big) \cdot \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big]^\mathsf{T} \; \mathrm{d}\Omega \cdot \delta \underline{\vec{x}} = \underline{\vec{t}} - \int\limits_{\Omega^h_0} \big[\, \vec{\nabla}_0 \underline{N}(\vec{X}) \,\big] \cdot \bm{\sigma}\big(\vec{x}_{(i)}\big) \; \mathrm{d}\Omega where .. math:: \mathbb{K}\big(\vec{x}_{(i)}\big) = \left. \frac{\partial \bm{\sigma}}{\partial \bm{\varepsilon}} \right|_{\vec{x}_{(i)}} : \mathbb{I}_s where the left part is the *constitutive tangent operator* and the right part comes from the strain definition. Note that this right part, the symmetrization using :math:`\mathbb{I}_s`, can often be omitted as many *constitutive tangent operators* already symmetrize. In a shorter notation, this is our iterative update: .. math:: \underline{\underline{\mathbb{K}}}_{(i)} \cdot \delta \underline{\vec{x}} = \underline{\vec{t}} - \underline{\vec{f}}_{(i)} with .. math:: \underline{\underline{\mathbb{K}}}_{(i)} = \int\limits_{\Omega^h_0} \big[\, \vec{\nabla}_0 \underline{N} \,\big] \cdot \mathbb{K}\big(\vec{x}_{(i)}\big) \cdot \big[\, \vec{\nabla}_0 \underline{N} \,\big]^\mathsf{T} \; \mathrm{d}\Omega and .. math:: \underline{\vec{f}}_{(i)} = \int\limits_{\Omega^h_0} \big[\, \vec{\nabla}_0 \underline{N} \,\big] \cdot \bm{\sigma}\big(\vec{x}_{(i)}\big) \; \mathrm{d}\Omega .. note:: This is a good point to study some examples: * :ref:`fem_examples_small-strain_linear_dense` We slowly work up to an iterative scheme starting from a linear problem, written, however, in such a way that the step towards a non-linear problem is small. * :ref:`fem_examples_small-strain_nonlinear_dense` Here we employ Newton-Raphson to solve the non-linear equilibrium equation. It is easy to see that once the above examples have been understood this step is indeed trivial. .. note:: **Newton-Raphson in one dimension** We try to find :math:`x` such that .. math:: r(x) = 0 We will make a guess for :math:`x` and (hopefully) iteratively improve this guess. This iterative value is denoted using :math:`x_{(i)}`. Therefore we will make use of the following Taylor expansion .. math:: r \big( x_{(i+1)} \big) = r \big( x_{(i)} \big) + \left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \delta x + \mathcal{O} \big( \delta x^2 \big) \approx 0 where .. math:: \delta x = x_{(i+1)} - x_{(i)} We now determine :math:`\delta x` by neglecting higher order terms, which results in .. math:: r \big( x_{(i)} \big) + \left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \delta x = 0 From which we obtain :math:`\delta x` as .. math:: \delta x = - \left[ \left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \right]^{-1} r \big( x_{(i)} \big) Thereafter we set .. math:: x_{(i+1)} = x_{(i)} + \delta x And check if we have reached our solution within a certain accuracy :math:`\epsilon`: .. math:: \left| r \big( x_{(i+1)} \big) \right| < \epsilon If not, we repeat the above steps until we do. The iterative scheme is well understood from the following illustration: .. image:: newton-raphson.svg :width: 300px :align: center Dynamics ======== Momentum balance ---------------- -We continue with our balance equation and add inertia and damping to it: +Weak form +^^^^^^^^^ + +We continue with our balance equation and add inertia to it: + +.. math:: + + \rho\, \ddot{\vec{x}} + = + \vec{\nabla} \cdot + \bm{\sigma}(\vec{x}) + \qquad + \vec{x} \in \Omega + +where :math:`\rho` is the mass density, and second time derivative of the position :math:`\vec{x}` is the acceleration :math:`\vec{a} = \ddot{\vec{x}}`. Note that this function is the continuum equivalent of :math:`\vec{f} = m \vec{a}`. + +Like before, we will solve this equation in a weak sense + +.. math:: + + \int\limits_\Omega + \rho(\vec{x})\; \vec{\phi}(\vec{X}) \cdot \ddot{\vec{x}} \; + \mathrm{d}\Omega + = + \int\limits_\Omega + \vec{\phi}(\vec{X}) + \cdot + \Big[\, + \vec{\nabla} + \cdot + \bm{\sigma}(\vec{x}) + \,\Big] \; + \mathrm{d}\Omega + \qquad + \forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d + +Integration by parts results in + +.. math:: + + \int\limits_\Omega + \rho(\vec{x})\; \vec{\phi}(\vec{X}) \cdot \ddot{\vec{x}} \; + \mathrm{d}\Omega + = + \int\limits_\Gamma + \vec{\phi}(\vec{X}) \cdot \vec{t}(\vec{x}) \; + \mathrm{d}\Gamma + - + \int\limits_\Omega + \big[\, \vec{\nabla} \vec{\phi}(\vec{X}) \,\big] + : + \bm{\sigma}(\vec{x}) \; + \mathrm{d}\Omega + \qquad + \forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d + +Which we will discretize as before: + +.. math:: + + \underline{\vec{\phi}}^\mathsf{T} \cdot + \int\limits_\Omega + \rho(\vec{x})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \; + \mathrm{d}\Omega \; + \underline{\ddot{\vec{x}}} + = + \underline{\vec{\phi}}^\mathsf{T} \cdot + \int\limits_\Gamma + \underline{N}(\vec{X})\; \vec{t}(\vec{x}) \; + \mathrm{d}\Gamma + - + \underline{\vec{\phi}}^\mathsf{T} \cdot + \int\limits_\Omega + \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big] + : + \bm{\sigma}(\vec{x}) \; + \mathrm{d}\Omega + \qquad + \forall \; \underline{\vec{\phi}} \in \mathbb{R}^d_n + +Which is independent of the test functions, hence: + +.. math:: + + \underbrace{ + \int\limits_\Omega + \rho(\vec{x})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \; + \mathrm{d}\Omega + }_{\underline{\underline{M}}(\vec{x})} \; + \underline{\ddot{\vec{x}}} + = + \underbrace{ + \int\limits_\Gamma + \underline{N}(\vec{X})\; \vec{t}(\vec{x}) \; + \mathrm{d}\Gamma + }_{\underline{\vec{t}}(\vec{x})} + - + \underbrace{ + \int\limits_\Omega + \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big] + : + \bm{\sigma}(\vec{x}) \; + \mathrm{d}\Omega + }_{\underline{\vec{f}}(\vec{x})} + +Which we can denote as follows + +.. math:: + + \underline{\underline{M}}(\vec{x})\; \underline{\ddot{\vec{x}}} + = + \underline{\vec{t}}(\vec{x}) + - + \underline{\vec{f}}(\vec{x}) + +Where :math:`\underline{\underline{M}}(\vec{x})` is the *mass matrix*, :math:`\underline{\vec{t}}(\vec{x})` are the *boundary tractions*, and :math:`\underline{\vec{f}}(\vec{x})` are the *internal forces*. + +.. note:: + + For problems where the local volume is conversed (either weakly to incompressible elasticity, or strongly by adding an incompressibility constraint) it make makes sense to assume the mass matrix constant, as any change of volume results in an equivalent change of the density. In that case + + .. math:: + + \int\limits_{\Omega} + \rho(\vec{x}) + \;\mathrm{d}\Omega + = + \int\limits_{\Omega_0} + \rho(\vec{X}) + \;\mathrm{d}\Omega_0 + + Which results in: + + .. math:: + + \underline{\underline{M}}(\vec{X}) + = + \int\limits_{\Omega_0} + \rho(\vec{X})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \; + \mathrm{d}\Omega_0 + = + \mathrm{constant} + +.. note:: + + To enhance computational efficiency, it may be a good option concentrate the mass to point masses on the nodes. This has to strong advantage that the mass matrix becomes diagonal. Consequently, instead of solving a linear system one just has to just do :math:`n_\mathrm{dof}` decoupled inversions (where :math:`n_\mathrm{dof}` is the number of degrees-of-freedom). + + See: :ref:`fem_examples_dynamic_diagonal-mass`. + +Time discretization +^^^^^^^^^^^^^^^^^^^ + +To solve the second order differential equation in time, one typically also discretizes time, however with some finite difference based scheme. To simplify notation below to following notation is used: for the velocity :math:`\vec{v} = \dot{\vec{x}}`, and for the acceleration :math:`\vec{a} = \ddot{\vec{x}}`. + +Verlet +"""""" + +1. Compute the velocity on a dummy time grid: + + .. math:: + + \vec{v}_{n+1/2} = \vec{v}_{n-1/2} + \Delta_t \; a_n + + Note that this involves solving for :math:`a_n`. + +2. Update the positions + + .. math:: + + x_{n+1} = x_n + \Delta_t v_{n + 1/2} + +Note that an important feature is that is time reversible. + +Damping +------- + +The equations such as presented in the previous section can lead to indefinite oscillations. For example when the constitutive response, hidden in :math:`\bm{\sigma}(\vec{x})`, is elasticity, there is no form of dissipation. In the absence of numerical inaccuracies this corresponds to elastic waves that travel through the system indefinitely. + +To avoid this one needs to add damping to the system. Two types of damping are considered: + +* Viscosity: + + .. math:: + + \rho\, \ddot{\vec{x}} + = + \vec{\nabla} \cdot + \bm{\sigma}(\vec{x}) + \textcolor{red}{+ \vec{\nabla} \cdot \bm{\sigma}_\eta (\vec{x})} + \qquad + \vec{x} \in \Omega + + + + +We continue with our balance equation and add inertia to it: + +.. math:: + + \rho\, \ddot{\vec{x}} + = + \vec{\nabla} \cdot + \bm{\sigma}(\vec{x}) + \qquad + \vec{x} \in \Omega + +where :math:`\rho` is the mass density. This function is the continuum equivalent of :math:`\vec{f} = m \vec{a}`. .. math:: \rho\, \ddot{\vec{x}} = \vec{\nabla} \cdot \bm{\sigma}(\vec{x}) + \eta\, \nabla^2\dot{\vec{x}} \qquad \vec{x} \in \Omega where :math:`\rho` is the density and :math:`\eta` the viscosity (a.k.a. the damping coefficient). The first and second time derivative of the position :math:`\vec{x}` are respectively the velocity :math:`\vec{v} = \dot{\vec{x}}` and the acceleration :math:`\vec{a} = \ddot{\vec{x}}`. We can generalize this as follows (which will also simplify our proceedings below) .. math:: \rho(\vec{x})\, \ddot{\vec{x}} = \vec{\nabla} \cdot \big[\, \bm{\sigma}(\vec{x}) + \bm{\sigma}_{\eta}(\vec{\dot{x}} ) \,\big] \qquad \vec{x} \in \Omega .. note:: To retrieve the original form .. math:: \bm{\sigma}_{\eta} = \eta\; \vec{\nabla} \dot{\vec{x}} But, we can now also use other expressions. For example, the damping equivalent of linear elasticity: .. math:: \bm{\sigma}_{\eta} (\vec{x}) = \mathbb{C}_{\eta} (\vec{x}) : \dot{\bm{\varepsilon}} (\vec{x}) with .. math:: \mathbb{C}_{\eta} (\vec{x}) = \kappa (\vec{x}) \bm{I} \otimes \bm{I} + 2 \gamma (\vec{x}) \mathbb{I}_d where :math:`\kappa` is the bulk viscosity while :math:`\gamma` is the shear viscosity. Furthermore .. math:: \dot{\bm{\varepsilon}} (\vec{x}) = \tfrac{1}{2} \big[\, \vec{\nabla} \dot{\vec{x}} + [\, \vec{\nabla} \dot{\vec{x}} \,]^T \,\big] Our original form is retrieved when :math:`\kappa = \tfrac{2}{3} \gamma`, both are independent of :math:`\vec{x}`, and :math:`\dot{\vec{x}}` possesses the necessary symmetries. Like before, we will solve this equation in a weak sense .. math:: \int\limits_\Omega \rho(\vec{x})\; \vec{\phi}(\vec{X}) \cdot \ddot{\vec{x}} \; \mathrm{d}\Omega = \int\limits_\Omega \vec{\phi}(\vec{X}) \cdot \Big[\, \vec{\nabla} \cdot \big[\, \bm{\sigma}(\vec{x}) + \bm{\sigma}_{\eta}(\vec{\dot{x}} ) \,\big] \,\Big] \; \mathrm{d}\Omega \qquad \forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d Integration by parts results in .. math:: \int\limits_\Omega \rho(\vec{x})\; \vec{\phi}(\vec{X}) \cdot \ddot{\vec{x}} \; \mathrm{d}\Omega = \int\limits_\Gamma \vec{\phi}(\vec{X}) \cdot \big[\, \vec{t}(\vec{x}) + \vec{t}_{\eta}(\vec{x}) \,\big] \; \mathrm{d}\Gamma - \int\limits_\Omega \big[\, \vec{\nabla} \vec{\phi}(\vec{X}) \,\big] : \big[\, \bm{\sigma}(\vec{x}) + \bm{\sigma}_{\eta}(\dot{\vec{x}}) \,\big] \; \mathrm{d}\Omega \qquad \forall \; \vec{\phi}(\vec{X}) \in \mathbb{R}^d Which we will discretize as before: .. math:: \underline{\vec{\phi}}^\mathsf{T} \cdot \int\limits_\Omega \rho(\vec{x})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \; \mathrm{d}\Omega \; \underline{\ddot{\vec{x}}} = \underline{\vec{\phi}}^\mathsf{T} \cdot \int\limits_\Gamma \underline{N}(\vec{X})\; \big[\, \vec{t}(\vec{x}) + \vec{t}_{\eta}(\vec{x}) \,\big] \; \mathrm{d}\Gamma - \underline{\vec{\phi}}^\mathsf{T} \cdot \int\limits_\Omega \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big] : \big[\, \bm{\sigma}(\vec{x}) + \bm{\sigma}_{\eta}(\dot{\vec{x}}) \,\big] \; \mathrm{d}\Omega \qquad \forall \; \underline{\vec{\phi}} \in \mathbb{R}^d_n Which is independent of the test functions, hence: .. math:: \int\limits_\Omega \rho(\vec{x})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \; \mathrm{d}\Omega \; \underline{\ddot{\vec{x}}} = \int\limits_\Gamma \underline{N}(\vec{X})\; \big[\, \vec{t}(\vec{x}) + \vec{t}_{\eta}(\vec{x}) \,\big] \; \mathrm{d}\Gamma - \int\limits_\Omega \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big] : \big[\, \bm{\sigma}(\vec{x}) + \bm{\sigma}_{\eta}(\dot{\vec{x}}) \,\big] \; \mathrm{d}\Omega Which we can denote as follows .. math:: \underline{\underline{M}}(\vec{x})\; \underline{\ddot{\vec{x}}} = \underline{\vec{t}}(\vec{x}) + \underline{\vec{t}}_{\eta}(\vec{x}) - \underline{\vec{f}}(\vec{x}) - \underline{\vec{f}}_{\eta}(\vec{x}) whereby we have introduced: * *Mass matrix* .. math:: \underline{\underline{M}}(\vec{x}) = \int\limits_\Omega \rho(\vec{x})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \; \mathrm{d}\Omega * *Boundary tractions* .. math:: \underline{\vec{t}}(\vec{x}) = \int\limits_\Gamma \underline{N}(\vec{X})\; \vec{t}(\vec{x}) \; \mathrm{d}\Gamma \qquad \mathrm{and} \qquad \underline{\vec{t}}_{\eta}(\vec{x}) = \int\limits_\Gamma \underline{N}(\vec{X})\; \vec{t}_{\eta}(\vec{x}) \; \mathrm{d}\Gamma * *Internal forces* .. math:: \underline{\vec{f}}(\vec{x}) = \int\limits_\Omega \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big] : \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega \qquad \mathrm{and} \qquad \underline{\vec{f}}(\vec{x}) = \int\limits_\Omega \big[\, \vec{\nabla} \underline{N}(\vec{X}) \,\big] : \bm{\sigma}_{\eta}(\dot{\vec{x}}) \; \mathrm{d}\Omega .. note:: In many problems it makes sense to assume the mass matrix constant, as any change of volume results in an equivalent change of the density, i.e. .. math:: \int\limits_{\Omega} \rho(\vec{x}) \;\mathrm{d}\Omega = \int\limits_{\Omega_0} \rho(\vec{X}) \;\mathrm{d}\Omega_0 This results in the following expression for the mass matrix: .. math:: \underline{\underline{M}}(\vec{X}) = \int\limits_{\Omega_0} \rho(\vec{X})\; \underline{N}(\vec{X})\; \underline{N}^\mathsf{T}(\vec{X}) \; \mathrm{d}\Omega_0 = \mathrm{constant} Time discretization ------------------- Here we will discuss several common time discretization steps. To simplify notation we will denote the velocity :math:`\vec{v} = \dot{\vec{x}}` and the acceleration :math:`\vec{a} = \ddot{\vec{x}}`. .. note:: Most time integration schemes result is some form like .. math:: \underline{\underline{M}}\; \underline{\vec{a}}_{n+1} = \underline{\vec{q}}_{n} where :math:`\underline{\vec{q}}_{n}` contains the boundary tractions and internal forces, including their damping equivalents. The subscript :math:`n` indicates that the variable is a known quantity, while :math:`n+1` indicates that it is an unknown quantity. To enhance computational efficiency, it may be a good option to approximate the mass matrix in such a way that it becomes diagonal. Consequently, no system has be solved to find :math:`\underline{\vec{a}}_{n+1}`. One only has to invert an array of scalars. Since in addition the mass matrix is almost often assumed constant, this factorization has to be performed only once for the entire simulation. Physically one can interpret this assumption as assuming the damping to be concentrated on the nodes. See: :ref:`fem_examples_dynamic_diagonal-mass`. .. note:: References `Syllabus of the course "Computational Physics (PY 502)" by Anders Sandvik, Department of Physics, Boston University `_. Velocity Verlet with damping ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1. Compute the position at :math:`t_{n+1} = t_{n} + \Delta_t`: .. math:: \vec{x}_{n+1} = \vec{x}_{n} + \Delta_t \vec{v}_{n} + \tfrac{1}{2} \Delta_t^2 \vec{a}_{n} 2. Estimate the velocity at :math:`t_{n+1} = t_{n} + \Delta_t`: .. math:: \hat{\vec{v}}_{n+1} = \vec{v}_{n} + \tfrac{1}{2} \Delta_t \Big[\, \vec{a}_{n} + \vec{a} ( \vec{x}_{n+1} , \vec{v}_{n} + \Delta_t \vec{a}_{n} , t_{n+1} ) \, \Big] 3. Correct :math:`\hat{\vec{v}}_{n+1}`: .. math:: \vec{v}_{n+1} = \vec{v}_{n} + \tfrac{1}{2} \Delta_t \Big[\, \vec{a}_{n} + \vec{a} ( \vec{x}_{n+1} , \hat{\vec{v}}_{n+1} , t_{n+1} ) \, \Big] Shape functions =============== In the Finite Element Method a geometry is discretized using nodes. The nodes are grouped in elements which define the domain :math:`\Omega^h_0`. The crux of the method is that nodal quantities, for example :math:`\vec{u}_i`, are extrapolated throughout the discretized domain :math:`\Omega^h_0` using shape functions :math:`N_i (\vec{X})`. Each shape function is globally supported, however, in such a way that :math:`N_i (\vec{X}) \neq 0` only in the elements containing node :math:`i`. It is, furthermore, imposed that :math:`N_i (\vec{X}_j) = \delta_{ij}`, i.e. it is one in the node :math:`i`, and zero in all other nodes. For a one-dimensional problem comprising four linear elements and five nodes the shape functions are sketched below (whereby the node numbers are in color, while the element numbers are in black, in between the nodes). .. image:: shape-functions-1d.svg :width: 600px :align: center From this it becomes obvious that :math:`N_i (\vec{X})` is polynomial through each of the nodes, and that :math:`\partial N_i / \partial \vec{X}` is discontinuous across element boundaries. Note once more that each of the shape functions :math:`N_i (X)` is globally supported, but zero outside the elements that contain the node :math:`i`. For node 2, the shape function is thus: .. image:: shape-functions-1d-node-2.svg :width: 600px :align: center As we can see, node 2 is only non-zero in elements 1 and 2, while it is zero everywhere else. To evaluate :math:`\vec{f}_2` we therefore only have to integrate on these elements (using `Isoparametric transformation and quadrature`_): .. math:: \vec{f}_2 = \int\limits_{\Omega^1} \big[\, \vec{\nabla} N^1_2(\vec{X}) \,\big] \cdot \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega + \int\limits_{\Omega^2} \big[\, \vec{\nabla} N^2_2(\vec{X}) \,\big] \cdot \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega By now it should be clear that the above allows us assemble :math:`\underline{f}` element-by-element. For this example, graphically this corresponds to the following sum: .. image:: shape-functions-1d-element-0.svg :width: 600px :align: center .. image:: shape-functions-1d-element-1.svg :width: 600px :align: center .. image:: shape-functions-1d-element-2.svg :width: 600px :align: center .. image:: shape-functions-1d-element-3.svg :width: 600px :align: center where the indices show that the *shape functions* are evaluated compared to some generic element definition (see `Isoparametric transformation and quadrature`_). Isoparametric transformation and quadrature =========================================== A very important concept in the Finite Element Method is the isoparametric transformation. It allows us to map an arbitrarily shaped element with volume :math:`\Omega^e` onto a generic *isoparametric element* of constant volume :math:`Q`. By using this mapping it is easy to perform numerical quadrature while even reusing an existing implementation (for example the one of `GooseFEM `_). .. image:: isoparametric-transform.svg :width: 600px :align: center The mapping between the generic domain :math:`Q` and the physical domain :math:`\Omega^e` is as follows .. math:: \vec{x} ( \vec{\xi} ) = \big[\, \underline{N}^{e} \,\big]^\mathsf{T} \underline{x}^e where the column :math:`\underline{x}^e` contains the real position vectors of the element nodes. In order to perform the quadrature on :math:`Q` we must also map the gradient operator: .. math:: \vec{\nabla}_{\xi}\, = \vec{e}_i \frac{\partial}{\partial \xi_i} = \vec{e}_i \frac{\partial x_j(\vec{\xi})}{\partial \xi_i} \frac{\partial}{\partial x_j} = \vec{e}_i \frac{\partial x_j(\vec{\xi})}{\partial \xi_i} \vec{e}_j \cdot \vec{e}_k \frac{\partial}{\partial x_k} = \big[\, \vec{\nabla}_{\xi}\, \vec{x}(\vec{\xi}) \,\big] \cdot \vec{\nabla} = \bm{J}(\vec{\xi}) \cdot \vec{\nabla} or .. math:: \vec{\nabla} = \bm{J}^{-1}(\vec{\xi}) \cdot \vec{\nabla}_{\xi}\, with .. math:: \bm{J}(\vec{\xi}) = \vec{\nabla}_{\xi}\, \vec{x}(\vec{\xi}) = \big[\, \vec{\nabla}_{\xi}\, \underline{N}^{e} \,\big]^\mathsf{T} \; \underline{x}^e Using the above: .. math:: \vec{\nabla} \underline{N}^{e} = \bm{J}^{-1}(\vec{\xi}) \cdot \big[\, \vec{\nabla}_{\xi}\, \underline{N}^{e} \,\big] We can now determine the mapping between the real and the master volume: .. math:: \mathrm{d} \Omega = \mathrm{d} \vec{x}_0 \times \mathrm{d} \vec{x}_1 \cdot \mathrm{d} \vec{x}_2 = \left[ \mathrm{d} \vec{x}_0 \cdot \bm{J}(\vec{\xi}) \right] \times \left[ \mathrm{d} \vec{x}_1 \cdot \bm{J}(\vec{\xi}) \right] \cdot \left[ \mathrm{d} \vec{x}_2 \cdot \bm{J}(\vec{\xi}) \right] = \det \big( \bm{J}(\vec{\xi}) \big)\, \mathrm{d} \vec{\xi}_0 \times \mathrm{d} \vec{\xi}_1 \cdot \mathrm{d} \vec{\xi}_2 = \det \big( \bm{J}(\vec{\xi}) \big)\, \mathrm{d} Q For example for the internal force this implies .. math:: \underline{\vec{f}^e} = \int\limits_{\Omega^e} \big[\, \vec{\nabla} \underline{N} \,\big] \cdot \bm{\sigma}(\vec{x}) \; \mathrm{d}\Omega = \int\limits_{Q} \big[\, \vec{\nabla} \underline{N} \,\big] \cdot \bm{\sigma}(\vec{x}) \; \det \big( \bm{J}(\vec{\xi}) \big) \; \mathrm{d}Q Numerical quadrature can be formulated (exactly) on the master element. It corresponds to taking the weighted sum of the integrand evaluated at specific *quadrature points* (or *integration-points*). Again, for our internal force: .. math:: \underline{\vec{f}^e} = \sum_{k}^{n_k} w_k \big[\, \vec{\nabla} \underline{N} \,\big]_{\vec{\xi} = \vec{\xi}_k} \cdot \bm{\sigma}\big(\vec{x}(\vec{\xi}_k)\big) \; \det \big( \bm{J}(\vec{\xi}_k) \big) \; .. note:: To obtain :math:`\vec{X}(\vec{\xi})`, :math:`\vec{\nabla}_0`, and :math:`\int\limits_{\Omega_0} . \;\mathrm{d}\Omega`, simply replace :math:`\underline{x}^e` with :math:`\underline{X}^e` in the first equation. For this reason the same element implementation (of for example `GooseFEM `_) can be used in small strain and finite strain (total Lagrange and updated Lagrange), proving either :math:`\underline{X}^e` or :math:`\underline{X}^e + \underline{u}^e` as input. .. note:: The details depend on the element type. Several standard elements types are implemented in `GooseFEM `_. diff --git a/docs/theory/readme.tex b/docs/theory/readme.tex new file mode 100644 index 0000000..f009cf5 --- /dev/null +++ b/docs/theory/readme.tex @@ -0,0 +1,34 @@ +%!TEX program = xelatex +\documentclass[garamond,namecite]{goose-article} + +\usepackage{lipsum,verbatim,mdframed} + +\title{% + The Finite Element Method -- some basic theory +} + +\author{T.W.J.~de~Geus} + +\contact{% + $^*$Contact: % + \href{http://goosefem.geus.me}{goosefem.geus.me}% + \hspace{1mm}--\hspace{1mm} % + \href{mailto:tom@geus.me}{tom@geus.me} % + \hspace{1mm}--\hspace{1mm} % + \href{http://www.geus.me}{www.geus.me}% +} + +\hypersetup{pdfauthor={T.W.J. de Geus}} + +\header{% + goose-article +} + +\begin{document} + +\maketitle + + +% \bibliography{example_refs} + +\end{document} diff --git a/docs/theory/unsrtnat.bst b/docs/theory/unsrtnat.bst new file mode 100644 index 0000000..1c7eddf --- /dev/null +++ b/docs/theory/unsrtnat.bst @@ -0,0 +1,1343 @@ +%% File: `unsrtnat.bst' +%% A modification of `unsrt.bst' for use with natbib package +%% +%% Copyright 1993-2007 Patrick W Daly +%% Max-Planck-Institut f\"ur Sonnensystemforschung +%% Max-Planck-Str. 2 +%% D-37191 Katlenburg-Lindau +%% Germany +%% E-mail: daly@mps.mpg.de +%% +%% This program can be redistributed and/or modified under the terms +%% of the LaTeX Project Public License Distributed from CTAN +%% archives in directory macros/latex/base/lppl.txt; either +%% version 1 of the License, or any later version. +%% + % Version and source file information: + % \ProvidesFile{natbst.mbs}[2007/11/26 1.93 (PWD)] + % + % BibTeX `plainnat' family + % version 0.99b for BibTeX versions 0.99a or later, + % for LaTeX versions 2.09 and 2e. + % + % For use with the `natbib.sty' package; emulates the corresponding + % member of the `plain' family, but with author-year citations. + % + % With version 6.0 of `natbib.sty', it may also be used for numerical + % citations, while retaining the commands \citeauthor, \citefullauthor, + % and \citeyear to print the corresponding information. + % + % For version 7.0 of `natbib.sty', the KEY field replaces missing + % authors/editors, and the date is left blank in \bibitem. + % + % Includes field EID for the sequence/citation number of electronic journals + % which is used instead of page numbers. + % + % Includes fields ISBN and ISSN. + % + % Includes field URL for Internet addresses. + % + % Includes field DOI for Digital Object Idenfifiers. + % + % Works best with the url.sty package of Donald Arseneau. + % + % Works with identical authors and year are further sorted by + % title text, as in the standard plain.bst etc. + % +ENTRY + { address + author + booktitle + chapter + doi + eprint + eid + edition + editor + howpublished + institution + isbn + issn + journal + key + month + note + number + organization + pages + publisher + school + series + title + type + url + volume + year + } + {} + { label extra.label sort.label short.list } + +INTEGERS { output.state before.all mid.sentence after.sentence after.block } + +FUNCTION {init.state.consts} +{ #0 'before.all := + #1 'mid.sentence := + #2 'after.sentence := + #3 'after.block := +} + +STRINGS { s t } + +FUNCTION {output.nonnull} +{ 's := + output.state mid.sentence = + { ", " * write$ } + { output.state after.block = + { add.period$ write$ + newline$ + "\newblock " write$ + } + { output.state before.all = + 'write$ + { add.period$ " " * write$ } + if$ + } + if$ + mid.sentence 'output.state := + } + if$ + s +} + +FUNCTION {output} +{ duplicate$ empty$ + 'pop$ + 'output.nonnull + if$ +} + +FUNCTION {output.check} +{ 't := + duplicate$ empty$ + { pop$ "empty " t * " in " * cite$ * warning$ } + 'output.nonnull + if$ +} + +FUNCTION {fin.entry} +{ add.period$ + write$ + newline$ +} + +FUNCTION {new.block} +{ output.state before.all = + 'skip$ + { after.block 'output.state := } + if$ +} + +FUNCTION {new.sentence} +{ output.state after.block = + 'skip$ + { output.state before.all = + 'skip$ + { after.sentence 'output.state := } + if$ + } + if$ +} + +FUNCTION {not} +{ { #0 } + { #1 } + if$ +} + +FUNCTION {and} +{ 'skip$ + { pop$ #0 } + if$ +} + +FUNCTION {or} +{ { pop$ #1 } + 'skip$ + if$ +} + +FUNCTION {new.block.checka} +{ empty$ + 'skip$ + 'new.block + if$ +} + +FUNCTION {new.block.checkb} +{ empty$ + swap$ empty$ + and + 'skip$ + 'new.block + if$ +} + +FUNCTION {new.sentence.checka} +{ empty$ + 'skip$ + 'new.sentence + if$ +} + +FUNCTION {new.sentence.checkb} +{ empty$ + swap$ empty$ + and + 'skip$ + 'new.sentence + if$ +} + +FUNCTION {field.or.null} +{ duplicate$ empty$ + { pop$ "" } + 'skip$ + if$ +} + +FUNCTION {emphasize} +{ duplicate$ empty$ + { pop$ "" } + { "\emph{" swap$ * "}" * } + if$ +} + +INTEGERS { nameptr namesleft numnames } + +FUNCTION {format.names} +{ 's := + #1 'nameptr := + s num.names$ 'numnames := + numnames 'namesleft := + { namesleft #0 > } + { s nameptr "{ff~}{vv~}{ll}{, jj}" format.name$ 't := + nameptr #1 > + { namesleft #1 > + { ", " * t * } + { numnames #2 > + { "," * } + 'skip$ + if$ + t "others" = + { " et~al." * } + { " and " * t * } + if$ + } + if$ + } + 't + if$ + nameptr #1 + 'nameptr := + namesleft #1 - 'namesleft := + } + while$ +} + +FUNCTION {format.key} +{ empty$ + { key field.or.null } + { "" } + if$ +} + +FUNCTION {format.authors} +{ author empty$ + { "" } + { author format.names } + if$ +} + +FUNCTION {format.editors} +{ editor empty$ + { "" } + { editor format.names + editor num.names$ #1 > + { ", editors" * } + { ", editor" * } + if$ + } + if$ +} + +FUNCTION {format.isbn} +{ isbn empty$ + { "" } + { new.block "ISBN " isbn * } + if$ +} + +FUNCTION {format.issn} +{ issn empty$ + { "" } + { new.block "ISSN " issn * } + if$ +} + +FUNCTION {format.url} +{ url empty$ + { "" } + { new.block "URL \url{" url * "}" * } + if$ +} + +FUNCTION {format.doi} +{ doi empty$ + { "" } + { new.block "\doi{" doi * "}" * } + if$ +} + +FUNCTION {format.eprint} +{ eprint empty$ + { "" } + { new.block "\eprint{" eprint * "}" * } + if$ +} + +FUNCTION {format.title} +{ title empty$ + { "" } + { title "t" change.case$ } + if$ +} + +FUNCTION {format.full.names} +{'s := + #1 'nameptr := + s num.names$ 'numnames := + numnames 'namesleft := + { namesleft #0 > } + { s nameptr + "{vv~}{ll}" format.name$ 't := + nameptr #1 > + { + namesleft #1 > + { ", " * t * } + { + numnames #2 > + { "," * } + 'skip$ + if$ + t "others" = + { " et~al." * } + { " and " * t * } + if$ + } + if$ + } + 't + if$ + nameptr #1 + 'nameptr := + namesleft #1 - 'namesleft := + } + while$ +} + +FUNCTION {author.editor.full} +{ author empty$ + { editor empty$ + { "" } + { editor format.full.names } + if$ + } + { author format.full.names } + if$ +} + +FUNCTION {author.full} +{ author empty$ + { "" } + { author format.full.names } + if$ +} + +FUNCTION {editor.full} +{ editor empty$ + { "" } + { editor format.full.names } + if$ +} + +FUNCTION {make.full.names} +{ type$ "book" = + type$ "inbook" = + or + 'author.editor.full + { type$ "proceedings" = + 'editor.full + 'author.full + if$ + } + if$ +} + +FUNCTION {output.bibitem} +{ newline$ + "\bibitem[" write$ + label write$ + ")" make.full.names duplicate$ short.list = + { pop$ } + { * } + if$ + "]{" * write$ + cite$ write$ + "}" write$ + newline$ + "" + before.all 'output.state := +} + +FUNCTION {n.dashify} +{ 't := + "" + { t empty$ not } + { t #1 #1 substring$ "-" = + { t #1 #2 substring$ "--" = not + { "--" * + t #2 global.max$ substring$ 't := + } + { { t #1 #1 substring$ "-" = } + { "-" * + t #2 global.max$ substring$ 't := + } + while$ + } + if$ + } + { t #1 #1 substring$ * + t #2 global.max$ substring$ 't := + } + if$ + } + while$ +} + +FUNCTION {format.date} +{ year duplicate$ empty$ + { "empty year in " cite$ * warning$ + pop$ "" } + 'skip$ + if$ + month empty$ + 'skip$ + { month + " " * swap$ * + } + if$ + extra.label * +} + +FUNCTION {format.btitle} +{ title emphasize +} + +FUNCTION {tie.or.space.connect} +{ duplicate$ text.length$ #3 < + { "~" } + { " " } + if$ + swap$ * * +} + +FUNCTION {either.or.check} +{ empty$ + 'pop$ + { "can't use both " swap$ * " fields in " * cite$ * warning$ } + if$ +} + +FUNCTION {format.bvolume} +{ volume empty$ + { "" } + { "volume" volume tie.or.space.connect + series empty$ + 'skip$ + { " of " * series emphasize * } + if$ + "volume and number" number either.or.check + } + if$ +} + +FUNCTION {format.number.series} +{ volume empty$ + { number empty$ + { series field.or.null } + { output.state mid.sentence = + { "number" } + { "Number" } + if$ + number tie.or.space.connect + series empty$ + { "there's a number but no series in " cite$ * warning$ } + { " in " * series * } + if$ + } + if$ + } + { "" } + if$ +} + +FUNCTION {format.edition} +{ edition empty$ + { "" } + { output.state mid.sentence = + { edition "l" change.case$ " edition" * } + { edition "t" change.case$ " edition" * } + if$ + } + if$ +} + +INTEGERS { multiresult } + +FUNCTION {multi.page.check} +{ 't := + #0 'multiresult := + { multiresult not + t empty$ not + and + } + { t #1 #1 substring$ + duplicate$ "-" = + swap$ duplicate$ "," = + swap$ "+" = + or or + { #1 'multiresult := } + { t #2 global.max$ substring$ 't := } + if$ + } + while$ + multiresult +} + +FUNCTION {format.pages} +{ pages empty$ + { "" } + { pages multi.page.check + { "pages" pages n.dashify tie.or.space.connect } + { "page" pages tie.or.space.connect } + if$ + } + if$ +} + +FUNCTION {format.eid} +{ eid empty$ + { "" } + { "art." eid tie.or.space.connect } + if$ +} + +FUNCTION {format.vol.num.pages} +{ volume field.or.null + number empty$ + 'skip$ + { "\penalty0 (" number * ")" * * + volume empty$ + { "there's a number but no volume in " cite$ * warning$ } + 'skip$ + if$ + } + if$ + pages empty$ + 'skip$ + { duplicate$ empty$ + { pop$ format.pages } + { ":\penalty0 " * pages n.dashify * } + if$ + } + if$ +} + +FUNCTION {format.vol.num.eid} +{ volume field.or.null + number empty$ + 'skip$ + { "\penalty0 (" number * ")" * * + volume empty$ + { "there's a number but no volume in " cite$ * warning$ } + 'skip$ + if$ + } + if$ + eid empty$ + 'skip$ + { duplicate$ empty$ + { pop$ format.eid } + { ":\penalty0 " * eid * } + if$ + } + if$ +} + +FUNCTION {format.chapter.pages} +{ chapter empty$ + 'format.pages + { type empty$ + { "chapter" } + { type "l" change.case$ } + if$ + chapter tie.or.space.connect + pages empty$ + 'skip$ + { ", " * format.pages * } + if$ + } + if$ +} + +FUNCTION {format.in.ed.booktitle} +{ booktitle empty$ + { "" } + { editor empty$ + { "In " booktitle emphasize * } + { "In " format.editors * ", " * booktitle emphasize * } + if$ + } + if$ +} + +FUNCTION {empty.misc.check} +{ author empty$ title empty$ howpublished empty$ + month empty$ year empty$ note empty$ + and and and and and + key empty$ not and + { "all relevant fields are empty in " cite$ * warning$ } + 'skip$ + if$ +} + +FUNCTION {format.thesis.type} +{ type empty$ + 'skip$ + { pop$ + type "t" change.case$ + } + if$ +} + +FUNCTION {format.tr.number} +{ type empty$ + { "Technical Report" } + 'type + if$ + number empty$ + { "t" change.case$ } + { number tie.or.space.connect } + if$ +} + +FUNCTION {format.article.crossref} +{ key empty$ + { journal empty$ + { "need key or journal for " cite$ * " to crossref " * crossref * + warning$ + "" + } + { "In \emph{" journal * "}" * } + if$ + } + { "In " } + if$ + " \citet{" * crossref * "}" * +} + +FUNCTION {format.book.crossref} +{ volume empty$ + { "empty volume in " cite$ * "'s crossref of " * crossref * warning$ + "In " + } + { "Volume" volume tie.or.space.connect + " of " * + } + if$ + editor empty$ + editor field.or.null author field.or.null = + or + { key empty$ + { series empty$ + { "need editor, key, or series for " cite$ * " to crossref " * + crossref * warning$ + "" * + } + { "\emph{" * series * "}" * } + if$ + } + 'skip$ + if$ + } + 'skip$ + if$ + " \citet{" * crossref * "}" * +} + +FUNCTION {format.incoll.inproc.crossref} +{ editor empty$ + editor field.or.null author field.or.null = + or + { key empty$ + { booktitle empty$ + { "need editor, key, or booktitle for " cite$ * " to crossref " * + crossref * warning$ + "" + } + { "In \emph{" booktitle * "}" * } + if$ + } + { "In " } + if$ + } + { "In " } + if$ + " \citet{" * crossref * "}" * +} + +FUNCTION {article} +{ output.bibitem + format.authors "author" output.check + author format.key output + new.block + format.title "title" output.check + new.block + crossref missing$ + { journal emphasize "journal" output.check + eid empty$ + { format.vol.num.pages output } + { format.vol.num.eid output } + if$ + format.date "year" output.check + } + { format.article.crossref output.nonnull + eid empty$ + { format.pages output } + { format.eid output } + if$ + } + if$ + format.issn output + format.doi output + format.eprint output + format.url output + new.block + note output + fin.entry +} + +FUNCTION {book} +{ output.bibitem + author empty$ + { format.editors "author and editor" output.check + editor format.key output + } + { format.authors output.nonnull + crossref missing$ + { "author and editor" editor either.or.check } + 'skip$ + if$ + } + if$ + new.block + format.btitle "title" output.check + crossref missing$ + { format.bvolume output + new.block + format.number.series output + new.sentence + publisher "publisher" output.check + address output + } + { new.block + format.book.crossref output.nonnull + } + if$ + format.edition output + format.date "year" output.check + format.isbn output + format.doi output + format.eprint output + format.url output + new.block + note output + fin.entry +} + +FUNCTION {booklet} +{ output.bibitem + format.authors output + author format.key output + new.block + format.title "title" output.check + howpublished address new.block.checkb + howpublished output + address output + format.date output + format.isbn output + format.doi output + format.eprint output + format.url output + new.block + note output + fin.entry +} + +FUNCTION {inbook} +{ output.bibitem + author empty$ + { format.editors "author and editor" output.check + editor format.key output + } + { format.authors output.nonnull + crossref missing$ + { "author and editor" editor either.or.check } + 'skip$ + if$ + } + if$ + new.block + format.btitle "title" output.check + crossref missing$ + { format.bvolume output + format.chapter.pages "chapter and pages" output.check + new.block + format.number.series output + new.sentence + publisher "publisher" output.check + address output + } + { format.chapter.pages "chapter and pages" output.check + new.block + format.book.crossref output.nonnull + } + if$ + format.edition output + format.date "year" output.check + format.isbn output + format.doi output + format.eprint output + format.url output + new.block + note output + fin.entry +} + +FUNCTION {incollection} +{ output.bibitem + format.authors "author" output.check + author format.key output + new.block + format.title "title" output.check + new.block + crossref missing$ + { format.in.ed.booktitle "booktitle" output.check + format.bvolume output + format.number.series output + format.chapter.pages output + new.sentence + publisher "publisher" output.check + address output + format.edition output + format.date "year" output.check + } + { format.incoll.inproc.crossref output.nonnull + format.chapter.pages output + } + if$ + format.isbn output + format.doi output + format.eprint output + format.url output + new.block + note output + fin.entry +} + +FUNCTION {inproceedings} +{ output.bibitem + format.authors "author" output.check + author format.key output + new.block + format.title "title" output.check + new.block + crossref missing$ + { format.in.ed.booktitle "booktitle" output.check + format.bvolume output + format.number.series output + format.pages output + address empty$ + { organization publisher new.sentence.checkb + organization output + publisher output + format.date "year" output.check + } + { address output.nonnull + format.date "year" output.check + new.sentence + organization output + publisher output + } + if$ + } + { format.incoll.inproc.crossref output.nonnull + format.pages output + } + if$ + format.isbn output + format.doi output + format.eprint output + format.url output + new.block + note output + fin.entry +} + +FUNCTION {conference} { inproceedings } + +FUNCTION {manual} +{ output.bibitem + format.authors output + author format.key output + new.block + format.btitle "title" output.check + organization address new.block.checkb + organization output + address output + format.edition output + format.date output + format.url output + new.block + note output + fin.entry +} + +FUNCTION {mastersthesis} +{ output.bibitem + format.authors "author" output.check + author format.key output + new.block + format.title "title" output.check + new.block + "Master's thesis" format.thesis.type output.nonnull + school "school" output.check + address output + format.date "year" output.check + format.url output + new.block + note output + fin.entry +} + +FUNCTION {misc} +{ output.bibitem + format.authors output + author format.key output + title howpublished new.block.checkb + format.title output + howpublished new.block.checka + howpublished output + format.date output + format.issn output + format.url output + new.block + note output + fin.entry + empty.misc.check +} + +FUNCTION {phdthesis} +{ output.bibitem + format.authors "author" output.check + author format.key output + new.block + format.btitle "title" output.check + new.block + "PhD thesis" format.thesis.type output.nonnull + school "school" output.check + address output + format.date "year" output.check + format.url output + new.block + note output + fin.entry +} + +FUNCTION {proceedings} +{ output.bibitem + format.editors output + editor format.key output + new.block + format.btitle "title" output.check + format.bvolume output + format.number.series output + address output + format.date "year" output.check + new.sentence + organization output + publisher output + format.isbn output + format.doi output + format.eprint output + format.url output + new.block + note output + fin.entry +} + +FUNCTION {techreport} +{ output.bibitem + format.authors "author" output.check + author format.key output + new.block + format.title "title" output.check + new.block + format.tr.number output.nonnull + institution "institution" output.check + address output + format.date "year" output.check + format.url output + new.block + note output + fin.entry +} + +FUNCTION {unpublished} +{ output.bibitem + format.authors "author" output.check + author format.key output + new.block + format.title "title" output.check + new.block + note "note" output.check + format.date output + format.url output + fin.entry +} + +FUNCTION {default.type} { misc } + + +MACRO {jan} {"January"} + +MACRO {feb} {"February"} + +MACRO {mar} {"March"} + +MACRO {apr} {"April"} + +MACRO {may} {"May"} + +MACRO {jun} {"June"} + +MACRO {jul} {"July"} + +MACRO {aug} {"August"} + +MACRO {sep} {"September"} + +MACRO {oct} {"October"} + +MACRO {nov} {"November"} + +MACRO {dec} {"December"} + + + +MACRO {acmcs} {"ACM Computing Surveys"} + +MACRO {acta} {"Acta Informatica"} + +MACRO {cacm} {"Communications of the ACM"} + +MACRO {ibmjrd} {"IBM Journal of Research and Development"} + +MACRO {ibmsj} {"IBM Systems Journal"} + +MACRO {ieeese} {"IEEE Transactions on Software Engineering"} + +MACRO {ieeetc} {"IEEE Transactions on Computers"} + +MACRO {ieeetcad} + {"IEEE Transactions on Computer-Aided Design of Integrated Circuits"} + +MACRO {ipl} {"Information Processing Letters"} + +MACRO {jacm} {"Journal of the ACM"} + +MACRO {jcss} {"Journal of Computer and System Sciences"} + +MACRO {scp} {"Science of Computer Programming"} + +MACRO {sicomp} {"SIAM Journal on Computing"} + +MACRO {tocs} {"ACM Transactions on Computer Systems"} + +MACRO {tods} {"ACM Transactions on Database Systems"} + +MACRO {tog} {"ACM Transactions on Graphics"} + +MACRO {toms} {"ACM Transactions on Mathematical Software"} + +MACRO {toois} {"ACM Transactions on Office Information Systems"} + +MACRO {toplas} {"ACM Transactions on Programming Languages and Systems"} + +MACRO {tcs} {"Theoretical Computer Science"} + + +READ + +FUNCTION {sortify} +{ purify$ + "l" change.case$ +} + +INTEGERS { len } + +FUNCTION {chop.word} +{ 's := + 'len := + s #1 len substring$ = + { s len #1 + global.max$ substring$ } + 's + if$ +} + +FUNCTION {format.lab.names} +{ 's := + s #1 "{vv~}{ll}" format.name$ + s num.names$ duplicate$ + #2 > + { pop$ " et~al." * } + { #2 < + 'skip$ + { s #2 "{ff }{vv }{ll}{ jj}" format.name$ "others" = + { " et~al." * } + { " and " * s #2 "{vv~}{ll}" format.name$ * } + if$ + } + if$ + } + if$ +} + +FUNCTION {author.key.label} +{ author empty$ + { key empty$ + { cite$ #1 #3 substring$ } + 'key + if$ + } + { author format.lab.names } + if$ +} + +FUNCTION {author.editor.key.label} +{ author empty$ + { editor empty$ + { key empty$ + { cite$ #1 #3 substring$ } + 'key + if$ + } + { editor format.lab.names } + if$ + } + { author format.lab.names } + if$ +} + +FUNCTION {author.key.organization.label} +{ author empty$ + { key empty$ + { organization empty$ + { cite$ #1 #3 substring$ } + { "The " #4 organization chop.word #3 text.prefix$ } + if$ + } + 'key + if$ + } + { author format.lab.names } + if$ +} + +FUNCTION {editor.key.organization.label} +{ editor empty$ + { key empty$ + { organization empty$ + { cite$ #1 #3 substring$ } + { "The " #4 organization chop.word #3 text.prefix$ } + if$ + } + 'key + if$ + } + { editor format.lab.names } + if$ +} + +FUNCTION {calc.short.authors} +{ type$ "book" = + type$ "inbook" = + or + 'author.editor.key.label + { type$ "proceedings" = + 'editor.key.organization.label + { type$ "manual" = + 'author.key.organization.label + 'author.key.label + if$ + } + if$ + } + if$ + 'short.list := +} + +FUNCTION {calc.label} +{ calc.short.authors + short.list + "(" + * + year duplicate$ empty$ + short.list key field.or.null = or + { pop$ "" } + 'skip$ + if$ + * + 'label := +} + +INTEGERS { seq.num } + +FUNCTION {init.seq} +{ #0 'seq.num :=} + +EXECUTE {init.seq} + +FUNCTION {int.to.fix} +{ "000000000" swap$ int.to.str$ * + #-1 #10 substring$ +} + + +FUNCTION {presort} +{ calc.label + label sortify + " " + * + seq.num #1 + 'seq.num := + seq.num int.to.fix + 'sort.label := + sort.label * + #1 entry.max$ substring$ + 'sort.key$ := +} + +ITERATE {presort} + +SORT + +STRINGS { longest.label last.label next.extra } + +INTEGERS { longest.label.width last.extra.num number.label } + +FUNCTION {initialize.longest.label} +{ "" 'longest.label := + #0 int.to.chr$ 'last.label := + "" 'next.extra := + #0 'longest.label.width := + #0 'last.extra.num := + #0 'number.label := +} + +FUNCTION {forward.pass} +{ last.label label = + { last.extra.num #1 + 'last.extra.num := + last.extra.num int.to.chr$ 'extra.label := + } + { "a" chr.to.int$ 'last.extra.num := + "" 'extra.label := + label 'last.label := + } + if$ + number.label #1 + 'number.label := +} + +FUNCTION {reverse.pass} +{ next.extra "b" = + { "a" 'extra.label := } + 'skip$ + if$ + extra.label 'next.extra := + extra.label + duplicate$ empty$ + 'skip$ + { "{\natexlab{" swap$ * "}}" * } + if$ + 'extra.label := + label extra.label * 'label := +} + +EXECUTE {initialize.longest.label} + +ITERATE {forward.pass} + +REVERSE {reverse.pass} + +FUNCTION {bib.sort.order} +{ sort.label 'sort.key$ := +} + +ITERATE {bib.sort.order} + +SORT + +FUNCTION {begin.bib} +{ preamble$ empty$ + 'skip$ + { preamble$ write$ newline$ } + if$ + "\begin{thebibliography}{" number.label int.to.str$ * "}" * + write$ newline$ + "\providecommand{\natexlab}[1]{#1}" + write$ newline$ + "\providecommand{\url}[1]{\texttt{#1}}" + write$ newline$ + "\expandafter\ifx\csname urlstyle\endcsname\relax" + write$ newline$ + " \providecommand{\eprint}[1]{eprint: #1}\else" + write$ newline$ + " \providecommand{\eprint}{eprint: \begingroup \urlstyle{rm}\Url}\fi" + write$ newline$ + "\expandafter\ifx\csname urlstyle\endcsname\relax" + write$ newline$ + " \providecommand{\doi}[1]{doi: #1}\else" + write$ newline$ + " \providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi" + write$ newline$ +} + +EXECUTE {begin.bib} + +EXECUTE {init.state.consts} + +ITERATE {call.type$} + +FUNCTION {end.bib} +{ newline$ + "\end{thebibliography}" write$ newline$ +} + +EXECUTE {end.bib} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ef0b952..3a3f709 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,109 +1,111 @@ # required to specify the c++ standard cmake_minimum_required(VERSION 3.0) # define project information project(GooseFEM) # automatically parse the version number file(READ "${CMAKE_CURRENT_SOURCE_DIR}/GooseFEM/GooseFEM.h" _goosefem_version_header) string(REGEX MATCH "define[ \t]+GOOSEFEM_WORLD_VERSION[ \t]+([0-9]+)" _goosefem_world_version_match "${_goosefem_version_header}") set(GOOSEFEM_WORLD_VERSION "${CMAKE_MATCH_1}") string(REGEX MATCH "define[ \t]+GOOSEFEM_MAJOR_VERSION[ \t]+([0-9]+)" _goosefem_major_version_match "${_goosefem_version_header}") set(GOOSEFEM_MAJOR_VERSION "${CMAKE_MATCH_1}") string(REGEX MATCH "define[ \t]+GOOSEFEM_MINOR_VERSION[ \t]+([0-9]+)" _goosefem_minor_version_match "${_goosefem_version_header}") set(GOOSEFEM_MINOR_VERSION "${CMAKE_MATCH_1}") set(GOOSEFEM_VERSION_NUMBER ${GOOSEFEM_WORLD_VERSION}.${GOOSEFEM_MAJOR_VERSION}.${GOOSEFEM_MINOR_VERSION}) # required for install include(CMakePackageConfigHelpers) include(GNUInstallDirs) # set c++ standard -> c++14 set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) if(NOT CMAKE_CXX_STANDARD OR CMAKE_CXX_STANDARD LESS 14) set(CMAKE_CXX_STANDARD 14) endif() # set paths where 'GooseFEM' will be installed # - root set(GOOSEFEM_ROOT_DIR ${CMAKE_INSTALL_PREFIX}) # - headers set(GOOSEFEM_INCLUDE_DIR "${CMAKE_INSTALL_PREFIX}/${INCLUDE_INSTALL_DIR}") # - headers set(INCLUDE_INSTALL_DIR "${CMAKE_INSTALL_INCLUDEDIR}/GooseFEM" CACHE PATH "Path, relative to CMAKE_PREFIX_PATH, where the GooseFEM header files are installed" ) # - CMake configuration set(CMAKEPACKAGE_INSTALL_DIR "${CMAKE_INSTALL_DATADIR}/GooseFEM" CACHE PATH "Path where GooseFEMConfig.cmake is installed" ) # - pkg-config set(PKGCONFIG_INSTALL_DIR "${CMAKE_INSTALL_DATADIR}/pkgconfig" CACHE PATH "Path where GooseFEM.pc is installed" ) # list all header files set(GOOSEFEM_HEADERS GooseFEM/GooseFEM.h GooseFEM/Mesh.h GooseFEM/Mesh.cpp GooseFEM/MeshTri3.h GooseFEM/MeshTri3.cpp GooseFEM/MeshQuad4.h GooseFEM/MeshQuad4.cpp GooseFEM/MeshHex8.h GooseFEM/MeshHex8.cpp GooseFEM/DynamicsDiagonal.h GooseFEM/DynamicsDiagonal.cpp + GooseFEM/OverdampedDynamicsDiagonal.h + GooseFEM/OverdampedDynamicsDiagonal.cpp ) # Disable pkg-config for native Windows builds if(NOT WIN32 OR NOT CMAKE_HOST_SYSTEM_NAME MATCHES Windows) option(GOOSEFEM_BUILD_PKGCONFIG "Build pkg-config .pc file for GooseFEM" ON) endif() # build pkg-config file -> fill in relevant values if(GOOSEFEM_BUILD_PKGCONFIG) configure_file(GooseFEM.pc.in GooseFEM.pc @ONLY) install(FILES ${CMAKE_CURRENT_BINARY_DIR}/GooseFEM.pc DESTINATION ${PKGCONFIG_INSTALL_DIR}) endif() # configure CMake file configure_package_config_file( ${CMAKE_CURRENT_SOURCE_DIR}/GooseFEMConfig.cmake.in ${CMAKE_CURRENT_BINARY_DIR}/GooseFEMConfig.cmake PATH_VARS GOOSEFEM_INCLUDE_DIR GOOSEFEM_ROOT_DIR INSTALL_DESTINATION ${CMAKEPACKAGE_INSTALL_DIR} NO_CHECK_REQUIRED_COMPONENTS_MACRO ) # install # - CMake-file install(FILES ${CMAKE_CURRENT_BINARY_DIR}/GooseFEMConfig.cmake DESTINATION ${CMAKEPACKAGE_INSTALL_DIR} ) # - headers install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/${GOOSEFEM_HEADERS} DESTINATION ${INCLUDE_INSTALL_DIR}/GooseFEM ) # print information to screen message(STATUS "") message(STATUS "+-------------------------------------------------------------------------------") message(STATUS "|") message(STATUS "| Use 'make install' to install in '${CMAKE_INSTALL_PREFIX}'") message(STATUS "| To specify a custom directory call") message(STATUS "| cmake .. -DCMAKE_INSTALL_PREFIX=yourprefix") message(STATUS "| For custom paths, add the following line to your '~/.bashrc'") message(STATUS "| export PKG_CONFIG_PATH=${CMAKE_INSTALL_PREFIX}/share/pkgconfig:$PKG_CONFIG_PATH") message(STATUS "|") message(STATUS "+-------------------------------------------------------------------------------") message(STATUS "") diff --git a/src/GooseFEM/DynamicsDiagonal.cpp b/src/GooseFEM/DynamicsDiagonal.cpp index 41cfc83..5fcb719 100644 --- a/src/GooseFEM/DynamicsDiagonal.cpp +++ b/src/GooseFEM/DynamicsDiagonal.cpp @@ -1,1111 +1,1134 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ -#ifndef GOOSEFEM_DYNAMICS_DIAGONAL_H -#define GOOSEFEM_DYNAMICS_DIAGONAL_H +#ifndef GOOSEFEM_DYNAMICS_DIAGONAL_CPP +#define GOOSEFEM_DYNAMICS_DIAGONAL_CPP // ------------------------------------------------------------------------------------------------- -#include "GooseFEM.h" +#include "DynamicsDiagonal.h" // ================================= GooseFEM::Dynamics::Diagonal ================================== namespace GooseFEM { namespace Dynamics { namespace Diagonal { // ===================================== Simulation - Periodic ===================================== template inline Periodic::Periodic( std::unique_ptr _elem, const MatD &_x, const MatS &_conn, const MatS &_dofs, double _dt ) +: elem(std::move(_elem)), conn(_conn), dofs(_dofs), x(_x), dt(_dt) { - // copy input - elem = std::move(_elem); - x = _x; - conn = _conn; - dofs = _dofs; - dt = _dt; - // compute sizes nnode = static_cast(x.rows()); ndim = static_cast(x.cols()); nelem = static_cast(conn.rows()); nne = static_cast(conn.cols()); ndof = dofs.maxCoeff()+1; // basic checks (mostly the user is 'trusted') assert( static_cast(dofs.size()) == nnode * ndim ); assert( ndof < nnode * ndim ); // allocate and zero-initialize nodal quantities u.conservativeResize(nnode,ndim); u.setZero(); v.conservativeResize(nnode,ndim); v.setZero(); a.conservativeResize(nnode,ndim); a.setZero(); // allocate and zero-initialize linear system (DOFs) Minv.conservativeResize(ndof); M .conservativeResize(ndof); D .conservativeResize(ndof); F .conservativeResize(ndof); V .conservativeResize(ndof); V_n .conservativeResize(ndof); V_n .setZero(); A .conservativeResize(ndof); A_n .conservativeResize(ndof); A_n .setZero(); // initialize all fields updated_x(); updated_u(true); updated_v(true); } // ------------------------------------------------------------------------------------------------- template inline void Periodic::velocityVerlet() { // (1) new positions (displacements) // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" u += dt * v + ( .5 * std::pow(dt,2.) ) * a; // - process update in displacements updated_u(); // (2a) estimate new velocities // - update velocities (DOFs) V.noalias() = V_n + dt * A; // - convert to nodal velocities (periodicity implies that several nodes depend on the same DOF) for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' // "v" == "V_n" -> new nodal velocities, and a 'back-up' // "u" -> new nodal displacements // The forces "F" correspond to this state of the system } // ------------------------------------------------------------------------------------------------- template inline void Periodic::Verlet() { // (1) new nodal positions (displacements) // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" u += dt * v + ( .5 * std::pow(dt,2.) ) * a; // - process update in displacements updated_u(); // - solve for accelerations (DOFs) A.noalias() = Minv.cwiseProduct( - F ); // - convert to nodal acceleration (periodicity implies that several nodes depend on the same DOF) for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' // "v" == "V_n" -> new nodal velocities, and a 'back-up' // "u" -> new nodal displacements // The forces "F" correspond to this state of the system } // ------------------------------------------------------------------------------------------------- +template +inline void Periodic::ground() +{ + // set accelerations and velocities zero (DOFs) + A.setZero(); + V.setZero(); + + // convert to nodal velocity/acceleration (several nodes may depend on the same DOF) + for ( size_t i=0; i inline void Periodic::updated_x() { // set the nodal positions of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->x(e,m,i) = x(conn(e,m),i); // signal update elem->updated_x(); } // ------------------------------------------------------------------------------------------------- template inline void Periodic::updated_u(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->u(e,m,i) = u(conn(e,m),i); // signal update elem->updated_u(); // update if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ------------------------------------------------------------------------------------------------- template inline void Periodic::updated_v(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->v(e,m,i) = v(conn(e,m),i); // signal update elem->updated_v(); // update if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ------------------------------------------------------------------------------------------------- template inline void Periodic::assemble_M() { // zero-initialize M.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - mass matrix, per thread ColD M_(ndof); M_.setZero(); // - assemble diagonal mass matrix, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) M_(dofs(conn(e,m),i)) += elem->M(e,m*ndim+i,m*ndim+i); // - reduce "M_" per thread to total "M" #pragma omp critical M += M_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); // compute inverse of the mass matrix Minv = M.cwiseInverse(); } // ------------------------------------------------------------------------------------------------- template inline void Periodic::assemble_D() { // zero-initialize D.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - damping matrix, per thread ColD D_(ndof); D_.setZero(); // - assemble diagonal damping matrix, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) D_(dofs(conn(e,m),i)) += elem->D(e,m*ndim+i,m*ndim+i); // - reduce "D_" per thread to total "D" #pragma omp critical D += D_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); } // ------------------------------------------------------------------------------------------------- template inline void Periodic::assemble_F() { // zero-initialize F.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - force, per thread ColD F_(ndof); F_.setZero(); // - assemble force, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) F_(dofs(conn(e,m),i)) += elem->f(e,m,i); // - reduce "F_" per thread to total "F" #pragma omp critical F += F_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); } // =================================== Simulation - SemiPeriodic =================================== template inline SemiPeriodic::SemiPeriodic( std::unique_ptr _elem, const MatD &_x, const MatS &_conn, const MatS &_dofs, const ColS &_fixedDofs, double _dt ) +: elem(std::move(_elem)), conn(_conn), dofs(_dofs), x(_x), fixedDofs(_fixedDofs), dt(_dt) { - // copy input - elem = std::move(_elem); - x = _x; - conn = _conn; - dofs = _dofs; - dt = _dt; - fixedDofs = _fixedDofs; - // compute sizes nfixed = static_cast(fixedDofs.size()); nnode = static_cast(x.rows()); ndim = static_cast(x.cols()); nelem = static_cast(conn.rows()); nne = static_cast(conn.cols()); ndof = dofs.maxCoeff()+1; // basic checks (mostly the user is 'trusted') - assert( dofs.size() == nnode * ndim ); - assert( ndof < nnode * ndim ); + assert( static_cast(dofs.size()) == nnode * ndim ); + assert( ndof < nnode * ndim ); #ifndef NDEBUG for ( size_t i = 0 ; i < nfixed ; ++i ) assert( fixedDofs(i) < ndof ); #endif // allocate and zero-initialize nodal quantities u.conservativeResize(nnode,ndim); u.setZero(); v.conservativeResize(nnode,ndim); v.setZero(); a.conservativeResize(nnode,ndim); a.setZero(); // allocate and zero-initialize linear system (DOFs) Minv.conservativeResize(ndof); M .conservativeResize(ndof); D .conservativeResize(ndof); F .conservativeResize(ndof); V .conservativeResize(ndof); V_n .conservativeResize(ndof); V_n .setZero(); A .conservativeResize(ndof); A_n .conservativeResize(ndof); A_n .setZero(); // fixed DOFs : default zero velocity and acceleration fixedV.conservativeResize(nfixed); fixedV.setZero(); fixedA.conservativeResize(nfixed); fixedA.setZero(); // initialize all fields updated_x(); updated_u(true); updated_v(true); } // ------------------------------------------------------------------------------------------------- template inline void SemiPeriodic::velocityVerlet() { // (1) new positions (displacements) // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" u += dt * v + ( .5 * std::pow(dt,2.) ) * a; // - process update in displacements updated_u(); // (2a) estimate new velocities // - update velocities (DOFs) V.noalias() = V_n + dt * A; // - apply the fixed velocities for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' // "v" == "V_n" -> new nodal velocities, and a 'back-up' // "u" -> new nodal displacements // The forces "F" correspond to this state of the system } // ------------------------------------------------------------------------------------------------- template inline void SemiPeriodic::Verlet() { // (1) new nodal positions (displacements) // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" u += dt * v + ( .5 * std::pow(dt,2.) ) * a; // - process update in displacements updated_u(); // - solve for accelerations (DOFs) A.noalias() = Minv.cwiseProduct( - F ); // - apply the fixed accelerations for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' // "v" == "V_n" -> new nodal velocities, and a 'back-up' // "u" -> new nodal displacements // The forces "F" correspond to this state of the system } // ------------------------------------------------------------------------------------------------- +template +inline void SemiPeriodic::ground() +{ + // set accelerations and velocities zero (DOFs) + A.setZero(); + V.setZero(); + + // convert to nodal velocity/acceleration (several nodes may depend on the same DOF) + for ( size_t i=0; i inline void SemiPeriodic::updated_x() { // set the nodal positions of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->x(e,m,i) = x(conn(e,m),i); // signal update elem->updated_x(); } // ------------------------------------------------------------------------------------------------- template inline void SemiPeriodic::updated_u(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->u(e,m,i) = u(conn(e,m),i); // signal update elem->updated_u(); // update if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ------------------------------------------------------------------------------------------------- template inline void SemiPeriodic::updated_v(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->v(e,m,i) = v(conn(e,m),i); // signal update elem->updated_v(); // update if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ------------------------------------------------------------------------------------------------- template inline void SemiPeriodic::assemble_M() { // zero-initialize M.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - mass matrix, per thread ColD M_(ndof); M_.setZero(); // - assemble diagonal mass matrix, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) M_(dofs(conn(e,m),i)) += elem->M(e,m*ndim+i,m*ndim+i); // - reduce "M_" per thread to total "M" #pragma omp critical M += M_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); // compute inverse of the mass matrix Minv = M.cwiseInverse(); } // ------------------------------------------------------------------------------------------------- template inline void SemiPeriodic::assemble_D() { // zero-initialize D.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - damping matrix, per thread ColD D_(ndof); D_.setZero(); // - assemble diagonal damping matrix, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) D_(dofs(conn(e,m),i)) += elem->D(e,m*ndim+i,m*ndim+i); // - reduce "D_" per thread to total "D" #pragma omp critical D += D_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); } // ------------------------------------------------------------------------------------------------- template inline void SemiPeriodic::assemble_F() { // zero-initialize F.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - force, per thread ColD F_(ndof); F_.setZero(); // - assemble force, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) F_(dofs(conn(e,m),i)) += elem->f(e,m,i); // - reduce "F_" per thread to total "F" #pragma omp critical F += F_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); } // ======================================== Element - Quad4 ======================================== namespace SmallStrain { // ------------------------------------------------------------------------------------------------- template inline 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}); // - vol .resize({nelem,nip}); vol_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 inline 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_, vol_; // 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)); vol_.map(&vol_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 vol_(k) = w_(k) * Jdet_; // - assemble element mass matrix // M(m+i,n+i) += N(m) * rho * vol * N(n); M_(k*2 ,k*2 ) = mat->rho(e,k) * vol_(k); M_(k*2+1,k*2+1) = mat->rho(e,k) * vol_(k); // - assemble element non-Galilean damping matrix // D(m+i,n+i) += N(m) * alpha * vol * N(n); D_(k*2 ,k*2 ) = mat->alpha(e,k) * vol_(k); D_(k*2+1,k*2+1) = mat->alpha(e,k) * vol_(k); } // Gaussian quadrature // ------------------- // pointer to element integration volume and weight vol_.map(&vol(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 vol_(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 inline 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 vol_; // 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)); vol_.map(&vol(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) * vol; for ( size_t m = 0 ; m < nne ; ++m ) { f_(m,0) += dNx_(m,0) * sig_(0,0) * vol_(k) + dNx_(m,1) * sig_(1,0) * vol_(k); f_(m,1) += dNx_(m,0) * sig_(0,1) * vol_(k) + dNx_(m,1) * sig_(1,1) * vol_(k); } } } // set signals changed_f = true; changed_M = false; changed_D = false; } } // ================================================================================================= template inline 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 vol_; // 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(&v (e)); vol_.map(&vol(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->epsdot(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) * vol; for ( size_t m = 0 ; m < nne ; ++m ) { f_(m,0) += dNx_(m,0) * sig_(0,0) * vol_(k) + dNx_(m,1) * sig_(1,0) * vol_(k); f_(m,1) += dNx_(m,0) * sig_(0,1) * vol_(k) + dNx_(m,1) * sig_(1,1) * vol_(k); } } } // set signals changed_f = true; changed_M = false; changed_D = false; } } // ================================================================================================= template inline cppmat::cartesian2d::tensor2s Quad4::mean_eps(size_t e) { cppmat::cartesian2d::tensor2s eps_, tot_eps_(0.0); cppmat::tiny::vector vol_; double tot_vol_ = 0.0; // pointer to integration volume vol_.map(&vol(e)); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to strain tensor (stored symmetric) eps_.map(&mat->eps(e,k)); // - add to average tot_eps_ += vol_(k) * eps_; tot_vol_ += vol_(k); } // return volume average return ( tot_eps_ / tot_vol_ ); } // ================================================================================================= template inline cppmat::cartesian2d::tensor2s Quad4::mean_sig(size_t e) { cppmat::cartesian2d::tensor2s sig_, tot_sig_(0.0); cppmat::tiny::vector vol_; double tot_vol_ = 0.0; // pointer to integration volume vol_.map(&vol(e)); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to strain tensor (stored symmetric) sig_.map(&mat->sig(e,k)); // - add to average tot_sig_ += vol_(k) * sig_; tot_vol_ += vol_(k); } // return volume average return ( tot_sig_ / tot_vol_ ); } // ================================================================================================= template inline cppmat::cartesian2d::tensor2s Quad4::mean_eps() { cppmat::cartesian2d::tensor2s eps_, tot_eps_(0.0); cppmat::tiny::vector vol_; double tot_vol_ = 0.0; // loop over all elements (in parallel) for ( size_t e = 0 ; e < nelem ; ++e ) { // pointer to integration volume vol_.map(&vol(e)); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to strain tensor (stored symmetric) eps_.map(&mat->eps(e,k)); // - add to average tot_eps_ += vol_(k) * eps_; tot_vol_ += vol_(k); } } // return volume average return ( tot_eps_ / tot_vol_ ); } // ================================================================================================= template inline cppmat::cartesian2d::tensor2s Quad4::mean_sig() { cppmat::cartesian2d::tensor2s sig_, tot_sig_(0.0); cppmat::tiny::vector vol_; double tot_vol_ = 0.0; // loop over all elements (in parallel) for ( size_t e = 0 ; e < nelem ; ++e ) { // pointer to integration volume vol_.map(&vol(e)); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to strain tensor (stored symmetric) sig_.map(&mat->sig(e,k)); // - add to average tot_sig_ += vol_(k) * sig_; tot_vol_ += vol_(k); } } // return volume average return ( tot_sig_ / tot_vol_ ); } // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/DynamicsDiagonal.h b/src/GooseFEM/DynamicsDiagonal.h index d8cffba..0e0fa5b 100644 --- a/src/GooseFEM/DynamicsDiagonal.h +++ b/src/GooseFEM/DynamicsDiagonal.h @@ -1,168 +1,170 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ -#ifndef GOOSEFEM_DYNAMICS_DIAGONAL_CPP -#define GOOSEFEM_DYNAMICS_DIAGONAL_CPP +#ifndef GOOSEFEM_DYNAMICS_DIAGONAL_H +#define GOOSEFEM_DYNAMICS_DIAGONAL_H // ------------------------------------------------------------------------------------------------- -#include "DynamicsDiagonal.h" +#include "GooseFEM.h" // ================================= GooseFEM::Dynamics::Diagonal ================================== namespace GooseFEM { namespace Dynamics { namespace Diagonal { // ===================================== Simulation - Periodic ===================================== template class Periodic { public: // variables // --------- // element/quadrature/material definition std::unique_ptr elem; // mesh: nodal position/displacement/velocity/acceleration, DOF-numbers, connectivity, dimensions size_t nnode, nelem, nne, ndim, ndof; MatS conn, dofs; MatD x, u, v, a; // linear system: columns (also "M" and "D" which are diagonal) ColD M, Minv, D, F, V, V_n, A, A_n; // time integration double dt, t=0.0; // constructor // ----------- Periodic(){}; Periodic( std::unique_ptr elem, const MatD &x, const MatS &conn, const MatS &dofs, double dt ); // functions // --------- - void velocityVerlet(); // one time step of time integrator - void Verlet(); // one time step of time integrator (no velocity dependence) + void velocityVerlet(); // one time step of the time integrator + void Verlet(); // one time step of the time integrator (no velocity dependence) + void ground(); // ground the system: set A = V = 0 (also for their history) void updated_x(); // process update in "x" void updated_u(bool init=false); // process update in "u", if init all possible updates are made void updated_v(bool init=false); // process update in "v", if init all possible updates are made void assemble_M(); // assemble the mass matrix from the element mass matrices void assemble_D(); // assemble the damping matrix from the element damping matrices void assemble_F(); // assemble the force from the element forces }; // =================================== Simulation - SemiPeriodic =================================== template class SemiPeriodic { public: // variables // --------- // element/quadrature/material definition std::unique_ptr elem; // mesh: nodal position/displacement/velocity/acceleration, DOF-numbers, connectivity, dimensions size_t nnode, nelem, nne, ndim, ndof; MatS conn, dofs; MatD x, u, v, a; // fixed DOFs: prescribed velocity/acceleration, DOF-numbers, dimensions size_t nfixed; ColS fixedDofs; ColD fixedV, fixedA; // linear system: columns (also "M" and "D" which are diagonal) ColD M, Minv, D, F, V, V_n, A, A_n; // time integration double dt, t=0.0; // constructor // ----------- SemiPeriodic(){}; SemiPeriodic( std::unique_ptr elem, const MatD &x, const MatS &conn, const MatS &dofs, const ColS &fixedDofs, double dt ); // functions // --------- - void velocityVerlet(); // one time step of time integrator - void Verlet(); // one time step of time integrator (no velocity dependence) + void velocityVerlet(); // one time step of the time integrator + void Verlet(); // one time step of the time integrator (no velocity dependence) + void ground(); // ground the system: set A = V = 0 (also for their history) void updated_x(); // process update in "x" void updated_u(bool init=false); // process update in "u", if init all possible updates are made void updated_v(bool init=false); // process update in "v", if init all possible updates are made void assemble_M(); // assemble the mass matrix from the element mass matrices void assemble_D(); // assemble the damping matrix from the element damping matrices void assemble_F(); // assemble the force from the element forces }; // ======================================== Element - Quad4 ======================================== namespace SmallStrain { 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, vol, dNxi_n, w_n, vol_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(); // return volume average stress and strain cppmat::cartesian2d::tensor2s mean_eps(size_t e); // of element "e" cppmat::cartesian2d::tensor2s mean_sig(size_t e); // of element "e" cppmat::cartesian2d::tensor2s mean_eps(); // of all elements cppmat::cartesian2d::tensor2s mean_sig(); // of all elements }; } // namespace ... // ================================================================================================= }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/GooseFEM.h b/src/GooseFEM/GooseFEM.h index e0e9d2d..6dc76a2 100644 --- a/src/GooseFEM/GooseFEM.h +++ b/src/GooseFEM/GooseFEM.h @@ -1,74 +1,76 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_H #define GOOSEFEM_H // ================================================================================================= #define _USE_MATH_DEFINES // to use "M_PI" from "math.h" #include #include #include #include #include #include #include #include #include #include // ================================================================================================= #define GOOSEFEM_WORLD_VERSION 0 #define GOOSEFEM_MAJOR_VERSION 0 #define GOOSEFEM_MINOR_VERSION 7 #define GOOSEFEM_VERSION_AT_LEAST(x,y,z) \ (GOOSEFEM_WORLD_VERSION>x || (GOOSEFEM_WORLD_VERSION>=x && \ (GOOSEFEM_MAJOR_VERSION>y || (GOOSEFEM_MAJOR_VERSION>=y && \ GOOSEFEM_MINOR_VERSION>=z)))) #define GOOSEFEM_VERSION(x,y,z) \ (GOOSEFEM_WORLD_VERSION==x && \ GOOSEFEM_MAJOR_VERSION==y && \ GOOSEFEM_MINOR_VERSION==z) // ================================================================================================= // dummy operation that can be use to suppress the "unused parameter" warnings #define UNUSED(p) ( (void)(p) ) // ================================================================================================= // alias Eigen types namespace GooseFEM { typedef Eigen::Matrix MatD; typedef Eigen::Matrix MatS; typedef Eigen::Matrix ColD; typedef Eigen::Matrix ColS; typedef Eigen::Matrix ColI; } // ================================================================================================= #include "Mesh.h" #include "MeshTri3.h" #include "MeshQuad4.h" #include "MeshHex8.h" #include "DynamicsDiagonal.h" +#include "OverdampedDynamicsDiagonal.h" #include "Mesh.cpp" #include "MeshTri3.cpp" #include "MeshQuad4.cpp" #include "MeshHex8.cpp" #include "DynamicsDiagonal.cpp" +#include "OverdampedDynamicsDiagonal.cpp" // ================================================================================================= #endif diff --git a/src/GooseFEM/DynamicsDiagonal.cpp b/src/GooseFEM/OverdampedDynamicsDiagonal.cpp similarity index 64% copy from src/GooseFEM/DynamicsDiagonal.cpp copy to src/GooseFEM/OverdampedDynamicsDiagonal.cpp index 41cfc83..e28b37d 100644 --- a/src/GooseFEM/DynamicsDiagonal.cpp +++ b/src/GooseFEM/OverdampedDynamicsDiagonal.cpp @@ -1,1111 +1,807 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ -#ifndef GOOSEFEM_DYNAMICS_DIAGONAL_H -#define GOOSEFEM_DYNAMICS_DIAGONAL_H +#ifndef GOOSEFEM_OVERDAMPEDDYNAMICS_DIAGONAL_CPP +#define GOOSEFEM_OVERDAMPEDDYNAMICS_DIAGONAL_CPP // ------------------------------------------------------------------------------------------------- -#include "GooseFEM.h" +#include "OverdampedDynamicsDiagonal.h" -// ================================= GooseFEM::Dynamics::Diagonal ================================== +// ============================ GooseFEM::OverdampedDynamics::Diagonal ============================= namespace GooseFEM { -namespace Dynamics { +namespace OverdampedDynamics { namespace Diagonal { // ===================================== Simulation - Periodic ===================================== template inline Periodic::Periodic( std::unique_ptr _elem, const MatD &_x, const MatS &_conn, const MatS &_dofs, double _dt ) +: elem(std::move(_elem)), conn(_conn), dofs(_dofs), x(_x), dt(_dt) { - // copy input - elem = std::move(_elem); - x = _x; - conn = _conn; - dofs = _dofs; - dt = _dt; - // compute sizes nnode = static_cast(x.rows()); ndim = static_cast(x.cols()); nelem = static_cast(conn.rows()); nne = static_cast(conn.cols()); ndof = dofs.maxCoeff()+1; // basic checks (mostly the user is 'trusted') assert( static_cast(dofs.size()) == nnode * ndim ); assert( ndof < nnode * ndim ); // allocate and zero-initialize nodal quantities u.conservativeResize(nnode,ndim); u.setZero(); v.conservativeResize(nnode,ndim); v.setZero(); - a.conservativeResize(nnode,ndim); a.setZero(); // allocate and zero-initialize linear system (DOFs) - Minv.conservativeResize(ndof); - M .conservativeResize(ndof); D .conservativeResize(ndof); + Dinv.conservativeResize(ndof); F .conservativeResize(ndof); V .conservativeResize(ndof); - V_n .conservativeResize(ndof); V_n .setZero(); - A .conservativeResize(ndof); - A_n .conservativeResize(ndof); A_n .setZero(); // initialize all fields updated_x(); updated_u(true); updated_v(true); } // ------------------------------------------------------------------------------------------------- template -inline void Periodic::velocityVerlet() +inline void Periodic::forwardEuler() { - // (1) new positions (displacements) - // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" - u += dt * v + ( .5 * std::pow(dt,2.) ) * a; - // - process update in displacements - updated_u(); - - // (2a) estimate new velocities - // - update velocities (DOFs) - V.noalias() = V_n + dt * A; + // (1) compute the velocities + // - solve for velocities (DOFs) + V.noalias() = Dinv.cwiseProduct( - F ); // - convert to nodal velocities (periodicity implies that several nodes depend on the same DOF) for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' - // "v" == "V_n" -> new nodal velocities, and a 'back-up' - // "u" -> new nodal displacements - // The forces "F" correspond to this state of the system -} - -// ------------------------------------------------------------------------------------------------- + // (2) update positions + u += dt * v; -template -inline void Periodic::Verlet() -{ - // (1) new nodal positions (displacements) - // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" - u += dt * v + ( .5 * std::pow(dt,2.) ) * a; - // - process update in displacements + // process update in displacements and velocities updated_u(); - // - solve for accelerations (DOFs) - A.noalias() = Minv.cwiseProduct( - F ); - // - convert to nodal acceleration (periodicity implies that several nodes depend on the same DOF) - for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' - // "v" == "V_n" -> new nodal velocities, and a 'back-up' - // "u" -> new nodal displacements - // The forces "F" correspond to this state of the system + // update time + t += dt; } // ------------------------------------------------------------------------------------------------- template inline void Periodic::updated_x() { // set the nodal positions of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->x(e,m,i) = x(conn(e,m),i); // signal update elem->updated_x(); } // ------------------------------------------------------------------------------------------------- template inline void Periodic::updated_u(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->u(e,m,i) = u(conn(e,m),i); // signal update elem->updated_u(); // update - if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ------------------------------------------------------------------------------------------------- template inline void Periodic::updated_v(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->v(e,m,i) = v(conn(e,m),i); // signal update elem->updated_v(); // update - if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ------------------------------------------------------------------------------------------------- -template -inline void Periodic::assemble_M() -{ - // zero-initialize - M.setZero(); - - // temporarily disable parallelization by Eigen - Eigen::setNbThreads(1); - // assemble - #pragma omp parallel - { - // - mass matrix, per thread - ColD M_(ndof); - M_.setZero(); - - // - assemble diagonal mass matrix, per thread - #pragma omp for - for ( size_t e = 0 ; e < nelem ; ++e ) - for ( size_t m = 0 ; m < nne ; ++m ) - for ( size_t i = 0 ; i < ndim ; ++i ) - M_(dofs(conn(e,m),i)) += elem->M(e,m*ndim+i,m*ndim+i); - - // - reduce "M_" per thread to total "M" - #pragma omp critical - M += M_; - } - // automatic parallelization by Eigen - Eigen::setNbThreads(0); - - // compute inverse of the mass matrix - Minv = M.cwiseInverse(); -} - -// ------------------------------------------------------------------------------------------------- - template inline void Periodic::assemble_D() { // zero-initialize D.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - damping matrix, per thread ColD D_(ndof); D_.setZero(); // - assemble diagonal damping matrix, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) D_(dofs(conn(e,m),i)) += elem->D(e,m*ndim+i,m*ndim+i); // - reduce "D_" per thread to total "D" #pragma omp critical D += D_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); + + // compute inverse + Dinv = D.cwiseInverse(); } // ------------------------------------------------------------------------------------------------- template inline void Periodic::assemble_F() { // zero-initialize F.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - force, per thread ColD F_(ndof); F_.setZero(); // - assemble force, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) F_(dofs(conn(e,m),i)) += elem->f(e,m,i); // - reduce "F_" per thread to total "F" #pragma omp critical F += F_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); } // =================================== Simulation - SemiPeriodic =================================== template inline SemiPeriodic::SemiPeriodic( std::unique_ptr _elem, const MatD &_x, const MatS &_conn, const MatS &_dofs, const ColS &_fixedDofs, double _dt ) +: elem(std::move(_elem)), conn(_conn), dofs(_dofs), x(_x), fixedDofs(_fixedDofs), dt(_dt) { - // copy input - elem = std::move(_elem); - x = _x; - conn = _conn; - dofs = _dofs; - dt = _dt; - fixedDofs = _fixedDofs; - // compute sizes nfixed = static_cast(fixedDofs.size()); nnode = static_cast(x.rows()); ndim = static_cast(x.cols()); nelem = static_cast(conn.rows()); nne = static_cast(conn.cols()); ndof = dofs.maxCoeff()+1; // basic checks (mostly the user is 'trusted') - assert( dofs.size() == nnode * ndim ); - assert( ndof < nnode * ndim ); + assert( static_cast(dofs.size()) == nnode * ndim ); + assert( ndof < nnode * ndim ); #ifndef NDEBUG for ( size_t i = 0 ; i < nfixed ; ++i ) assert( fixedDofs(i) < ndof ); #endif // allocate and zero-initialize nodal quantities u.conservativeResize(nnode,ndim); u.setZero(); v.conservativeResize(nnode,ndim); v.setZero(); - a.conservativeResize(nnode,ndim); a.setZero(); // allocate and zero-initialize linear system (DOFs) - Minv.conservativeResize(ndof); - M .conservativeResize(ndof); D .conservativeResize(ndof); + Dinv.conservativeResize(ndof); F .conservativeResize(ndof); V .conservativeResize(ndof); - V_n .conservativeResize(ndof); V_n .setZero(); - A .conservativeResize(ndof); - A_n .conservativeResize(ndof); A_n .setZero(); - // fixed DOFs : default zero velocity and acceleration + // fixed DOFs : default zero velocity fixedV.conservativeResize(nfixed); fixedV.setZero(); - fixedA.conservativeResize(nfixed); fixedA.setZero(); // initialize all fields updated_x(); updated_u(true); updated_v(true); } // ------------------------------------------------------------------------------------------------- template -inline void SemiPeriodic::velocityVerlet() +inline void SemiPeriodic::forwardEuler() { - // (1) new positions (displacements) - // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" - u += dt * v + ( .5 * std::pow(dt,2.) ) * a; - // - process update in displacements - updated_u(); - - // (2a) estimate new velocities - // - update velocities (DOFs) - V.noalias() = V_n + dt * A; - // - apply the fixed velocities - for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' - // "v" == "V_n" -> new nodal velocities, and a 'back-up' - // "u" -> new nodal displacements - // The forces "F" correspond to this state of the system -} - -// ------------------------------------------------------------------------------------------------- + // (2) update positions + u += dt * v; -template -inline void SemiPeriodic::Verlet() -{ - // (1) new nodal positions (displacements) - // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" - u += dt * v + ( .5 * std::pow(dt,2.) ) * a; - // - process update in displacements + // process update in displacements and velocities updated_u(); - // - solve for accelerations (DOFs) - A.noalias() = Minv.cwiseProduct( - F ); - // - apply the fixed accelerations - for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' - // "v" == "V_n" -> new nodal velocities, and a 'back-up' - // "u" -> new nodal displacements - // The forces "F" correspond to this state of the system + // update time + t += dt; } // ------------------------------------------------------------------------------------------------- template inline void SemiPeriodic::updated_x() { // set the nodal positions of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->x(e,m,i) = x(conn(e,m),i); // signal update elem->updated_x(); } // ------------------------------------------------------------------------------------------------- template inline void SemiPeriodic::updated_u(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->u(e,m,i) = u(conn(e,m),i); // signal update elem->updated_u(); // update - if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ------------------------------------------------------------------------------------------------- template inline void SemiPeriodic::updated_v(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->v(e,m,i) = v(conn(e,m),i); // signal update elem->updated_v(); // update - if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ------------------------------------------------------------------------------------------------- -template -inline void SemiPeriodic::assemble_M() -{ - // zero-initialize - M.setZero(); - - // temporarily disable parallelization by Eigen - Eigen::setNbThreads(1); - // assemble - #pragma omp parallel - { - // - mass matrix, per thread - ColD M_(ndof); - M_.setZero(); - - // - assemble diagonal mass matrix, per thread - #pragma omp for - for ( size_t e = 0 ; e < nelem ; ++e ) - for ( size_t m = 0 ; m < nne ; ++m ) - for ( size_t i = 0 ; i < ndim ; ++i ) - M_(dofs(conn(e,m),i)) += elem->M(e,m*ndim+i,m*ndim+i); - - // - reduce "M_" per thread to total "M" - #pragma omp critical - M += M_; - } - // automatic parallelization by Eigen - Eigen::setNbThreads(0); - - // compute inverse of the mass matrix - Minv = M.cwiseInverse(); -} - -// ------------------------------------------------------------------------------------------------- - template inline void SemiPeriodic::assemble_D() { // zero-initialize D.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - damping matrix, per thread ColD D_(ndof); D_.setZero(); // - assemble diagonal damping matrix, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) D_(dofs(conn(e,m),i)) += elem->D(e,m*ndim+i,m*ndim+i); // - reduce "D_" per thread to total "D" #pragma omp critical D += D_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); + + // compute inverse + Dinv = D.cwiseInverse(); } // ------------------------------------------------------------------------------------------------- template inline void SemiPeriodic::assemble_F() { // zero-initialize F.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - force, per thread ColD F_(ndof); F_.setZero(); // - assemble force, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) F_(dofs(conn(e,m),i)) += elem->f(e,m,i); // - reduce "F_" per thread to total "F" #pragma omp critical F += F_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); } // ======================================== Element - Quad4 ======================================== namespace SmallStrain { // ------------------------------------------------------------------------------------------------- template inline 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}); // - vol .resize({nelem,nip}); vol_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 inline 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_, vol_; // 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)); + // pointer to element damping matrix, nodal volume, and integration weight D_ .map(&D (e)); vol_.map(&vol_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 vol_(k) = w_(k) * Jdet_; - // - assemble element mass matrix - // M(m+i,n+i) += N(m) * rho * vol * N(n); - M_(k*2 ,k*2 ) = mat->rho(e,k) * vol_(k); - M_(k*2+1,k*2+1) = mat->rho(e,k) * vol_(k); - // - assemble element non-Galilean damping matrix // D(m+i,n+i) += N(m) * alpha * vol * N(n); D_(k*2 ,k*2 ) = mat->alpha(e,k) * vol_(k); D_(k*2+1,k*2+1) = mat->alpha(e,k) * vol_(k); } // Gaussian quadrature // ------------------- // pointer to element integration volume and weight vol_.map(&vol(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 vol_(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 inline 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 vol_; // 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)); vol_.map(&vol(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) * vol; for ( size_t m = 0 ; m < nne ; ++m ) { f_(m,0) += dNx_(m,0) * sig_(0,0) * vol_(k) + dNx_(m,1) * sig_(1,0) * vol_(k); f_(m,1) += dNx_(m,0) * sig_(0,1) * vol_(k) + dNx_(m,1) * sig_(1,1) * vol_(k); } } } // set signals changed_f = true; - changed_M = false; changed_D = false; } } // ================================================================================================= template inline 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 vol_; - - // 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(&v (e)); - vol_.map(&vol(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->epsdot(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) * vol; - for ( size_t m = 0 ; m < nne ; ++m ) - { - f_(m,0) += dNx_(m,0) * sig_(0,0) * vol_(k) + dNx_(m,1) * sig_(1,0) * vol_(k); - f_(m,1) += dNx_(m,0) * sig_(0,1) * vol_(k) + dNx_(m,1) * sig_(1,1) * vol_(k); - } - } - } - - // set signals changed_f = true; - changed_M = false; changed_D = false; } -} // ================================================================================================= template inline cppmat::cartesian2d::tensor2s Quad4::mean_eps(size_t e) { cppmat::cartesian2d::tensor2s eps_, tot_eps_(0.0); cppmat::tiny::vector vol_; double tot_vol_ = 0.0; // pointer to integration volume vol_.map(&vol(e)); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to strain tensor (stored symmetric) eps_.map(&mat->eps(e,k)); // - add to average tot_eps_ += vol_(k) * eps_; tot_vol_ += vol_(k); } // return volume average return ( tot_eps_ / tot_vol_ ); } // ================================================================================================= template inline cppmat::cartesian2d::tensor2s Quad4::mean_sig(size_t e) { cppmat::cartesian2d::tensor2s sig_, tot_sig_(0.0); cppmat::tiny::vector vol_; double tot_vol_ = 0.0; // pointer to integration volume vol_.map(&vol(e)); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to strain tensor (stored symmetric) sig_.map(&mat->sig(e,k)); // - add to average tot_sig_ += vol_(k) * sig_; tot_vol_ += vol_(k); } // return volume average return ( tot_sig_ / tot_vol_ ); } // ================================================================================================= template inline cppmat::cartesian2d::tensor2s Quad4::mean_eps() { cppmat::cartesian2d::tensor2s eps_, tot_eps_(0.0); cppmat::tiny::vector vol_; double tot_vol_ = 0.0; // loop over all elements (in parallel) for ( size_t e = 0 ; e < nelem ; ++e ) { // pointer to integration volume vol_.map(&vol(e)); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to strain tensor (stored symmetric) eps_.map(&mat->eps(e,k)); // - add to average tot_eps_ += vol_(k) * eps_; tot_vol_ += vol_(k); } } // return volume average return ( tot_eps_ / tot_vol_ ); } // ================================================================================================= template inline cppmat::cartesian2d::tensor2s Quad4::mean_sig() { cppmat::cartesian2d::tensor2s sig_, tot_sig_(0.0); cppmat::tiny::vector vol_; double tot_vol_ = 0.0; // loop over all elements (in parallel) for ( size_t e = 0 ; e < nelem ; ++e ) { // pointer to integration volume vol_.map(&vol(e)); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to strain tensor (stored symmetric) sig_.map(&mat->sig(e,k)); // - add to average tot_sig_ += vol_(k) * sig_; tot_vol_ += vol_(k); } } // return volume average return ( tot_sig_ / tot_vol_ ); } // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/DynamicsDiagonal.h b/src/GooseFEM/OverdampedDynamicsDiagonal.h similarity index 72% copy from src/GooseFEM/DynamicsDiagonal.h copy to src/GooseFEM/OverdampedDynamicsDiagonal.h index d8cffba..1242628 100644 --- a/src/GooseFEM/DynamicsDiagonal.h +++ b/src/GooseFEM/OverdampedDynamicsDiagonal.h @@ -1,168 +1,164 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ -#ifndef GOOSEFEM_DYNAMICS_DIAGONAL_CPP -#define GOOSEFEM_DYNAMICS_DIAGONAL_CPP +#ifndef GOOSEFEM_OVERDAMPEDDYNAMICS_DIAGONAL_H +#define GOOSEFEM_OVERDAMPEDDYNAMICS_DIAGONAL_H // ------------------------------------------------------------------------------------------------- -#include "DynamicsDiagonal.h" +#include "GooseFEM.h" -// ================================= GooseFEM::Dynamics::Diagonal ================================== +// ============================ GooseFEM::OverdampedDynamics::Diagonal ============================= namespace GooseFEM { -namespace Dynamics { +namespace OverdampedDynamics { namespace Diagonal { // ===================================== Simulation - Periodic ===================================== template class Periodic { public: // variables // --------- // element/quadrature/material definition std::unique_ptr elem; - // mesh: nodal position/displacement/velocity/acceleration, DOF-numbers, connectivity, dimensions + // mesh: nodal position/displacement/velocity, DOF-numbers, connectivity, dimensions size_t nnode, nelem, nne, ndim, ndof; MatS conn, dofs; - MatD x, u, v, a; + MatD x, u, v; - // linear system: columns (also "M" and "D" which are diagonal) - ColD M, Minv, D, F, V, V_n, A, A_n; + // linear system: columns (also "D" which is diagonal) + ColD D, Dinv, F, V; // time integration double dt, t=0.0; // constructor // ----------- Periodic(){}; Periodic( std::unique_ptr elem, const MatD &x, const MatS &conn, const MatS &dofs, double dt ); // functions // --------- - void velocityVerlet(); // one time step of time integrator - void Verlet(); // one time step of time integrator (no velocity dependence) + void forwardEuler(); // one time step of the time integrator void updated_x(); // process update in "x" void updated_u(bool init=false); // process update in "u", if init all possible updates are made void updated_v(bool init=false); // process update in "v", if init all possible updates are made - void assemble_M(); // assemble the mass matrix from the element mass matrices void assemble_D(); // assemble the damping matrix from the element damping matrices void assemble_F(); // assemble the force from the element forces }; // =================================== Simulation - SemiPeriodic =================================== template class SemiPeriodic { public: // variables // --------- // element/quadrature/material definition std::unique_ptr elem; - // mesh: nodal position/displacement/velocity/acceleration, DOF-numbers, connectivity, dimensions + // mesh: nodal position/displacement/velocity, DOF-numbers, connectivity, dimensions size_t nnode, nelem, nne, ndim, ndof; MatS conn, dofs; - MatD x, u, v, a; + MatD x, u, v; - // fixed DOFs: prescribed velocity/acceleration, DOF-numbers, dimensions + // fixed DOFs: prescribed velocity, DOF-numbers, dimensions size_t nfixed; ColS fixedDofs; - ColD fixedV, fixedA; + ColD fixedV; - // linear system: columns (also "M" and "D" which are diagonal) - ColD M, Minv, D, F, V, V_n, A, A_n; + // linear system: columns (also "D" which is diagonal) + ColD D, Dinv, F, V; // time integration double dt, t=0.0; // constructor // ----------- SemiPeriodic(){}; SemiPeriodic( std::unique_ptr elem, const MatD &x, const MatS &conn, const MatS &dofs, const ColS &fixedDofs, double dt ); // functions // --------- - void velocityVerlet(); // one time step of time integrator - void Verlet(); // one time step of time integrator (no velocity dependence) + void forwardEuler(); // one time step of the time integrator void updated_x(); // process update in "x" void updated_u(bool init=false); // process update in "u", if init all possible updates are made void updated_v(bool init=false); // process update in "v", if init all possible updates are made - void assemble_M(); // assemble the mass matrix from the element mass matrices void assemble_D(); // assemble the damping matrix from the element damping matrices void assemble_F(); // assemble the force from the element forces }; // ======================================== Element - Quad4 ======================================== namespace SmallStrain { 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, vol, dNxi_n, w_n, vol_n; + cppmat::matrix x, u, v, f, D, dNx, dNxi, w, vol, dNxi_n, w_n, vol_n; // dimensions size_t nelem, nip=4, nne=4, ndim=2; // signals for solvers - bool changed_f=false, changed_M=true, changed_D=true; + bool changed_f=false, 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(); // return volume average stress and strain cppmat::cartesian2d::tensor2s mean_eps(size_t e); // of element "e" cppmat::cartesian2d::tensor2s mean_sig(size_t e); // of element "e" cppmat::cartesian2d::tensor2s mean_eps(); // of all elements cppmat::cartesian2d::tensor2s mean_sig(); // of all elements }; } // namespace ... // ================================================================================================= }}} // namespace ... // ================================================================================================= #endif