diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0c6a348..70abc4e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,105 +1,106 @@ # 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/Macros.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/DynamicsDiagonalLinearStrainQuad4.h + GooseFEM/DynamicsDiagonalSemiPeriodic.h GooseFEM/DynamicsDiagonalPeriodic.h GooseFEM/GooseFEM.h GooseFEM/Macros.h GooseFEM/Mesh.h GooseFEM/MeshQuad4.h GooseFEM/MeshTri3.h ) # 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/DynamicsDiagonalPeriodic.h b/src/GooseFEM/DynamicsDiagonalPeriodic.h index 6924631..720f8bb 100644 --- a/src/GooseFEM/DynamicsDiagonalPeriodic.h +++ b/src/GooseFEM/DynamicsDiagonalPeriodic.h @@ -1,351 +1,351 @@ /* ========================================== DESCRIPTION ========================================== (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_DYNAMICS_DIAGONAL_PERIODIC_H #define GOOSEFEM_DYNAMICS_DIAGONAL_PERIODIC_H #include "Macros.h" #include // ------------------------------------------------------------------------------------------------- namespace GooseFEM { namespace Dynamics { namespace Diagonal { // =========================================== OVERVIEW ============================================ template class Periodic { public: // variables // --------- // element/quadrature/material definition (should match mesh) std::shared_ptr elem; // mesh : nodal quantities and connectivity - MatS dofs; // DOF-numbers of each node (column of 'vectors') - MatS conn; // node numbers of each element (column of elements) - MatD x0; // initial positions of each node (column of vectors) - MatD u; // displacements of each node (column of vectors) - MatD v; // velocities of each node (column of vectors) - MatD a; // accelerations of each node (column of vectors) - size_t nnode; // number of nodes - size_t nelem; // number of elements - size_t nne; // number of nodes-per-element - size_t ndim; // number of dimensions - size_t ndof; // number of DOFs (after eliminating periodic dependencies) + MatS dofs; // DOF-numbers of each node (column of 'vectors') + MatS conn; // node numbers of each element (column of elements) + MatD x0; // initial positions of each node (column of vectors) + MatD u; // displacements of each node (column of vectors) + MatD v; // velocities of each node (column of vectors) + MatD a; // accelerations of each node (column of vectors) + size_t nnode; // number of nodes + size_t nelem; // number of elements + size_t nne; // number of nodes-per-element + size_t ndim; // number of dimensions + size_t ndof; // number of DOFs (after eliminating periodic dependencies) // linear system - ColD Minv; // inverse of mass matrix (constructed diagonal -> inverse == diagonal) - ColD Fu; // force that depends on displacement (column) - ColD Fv; // force that depends on velocity (and displacement) (column) - ColD V; // velocity (column) - ColD V_n; // velocity (column, last increment) - ColD A; // acceleration (column) - ColD A_n; // acceleration (column, last increment) + ColD Minv; // inverse of mass matrix (constructed diagonal -> inverse == diagonal) + ColD Fu; // force that depends on displacement (column) + ColD Fv; // force that depends on velocity (and displacement) (column) + ColD V; // velocity (column) + ColD V_n; // velocity (column, last increment) + ColD A; // acceleration (column) + ColD A_n; // acceleration (column, last increment) // time integration - double dt; // time step - double t = 0.0; // current time - double alpha; // non-Galilean damping coefficient + double dt; // time step + double t = 0.0; // current time + double alpha; // non-Galilean damping coefficient // constructor // ----------- Periodic(std::shared_ptr elem, const MatD &x0, const MatS &conn, const MatS &dofs, double dt, double alpha=0.0 ); // functions // --------- void velocityVerlet(); // one time step of time integrator void Verlet(); // one time step of time integrator (Fv == 0) void computeMinv(); // element loop: inverse mass matrix <- "elem->computeM" void computeFu(); // element loop: displacement dependent forces <- "elem->computeFu" void computeFv(); // element loop: velocity dependent forces <- "elem->computeFv" void post(); // element loop: post-process <- "elem->post" }; // ========================================== SOURCE CODE ========================================== template Periodic::Periodic( std::shared_ptr _elem, const MatD &_x0, const MatS &_conn, const MatS &_dofs, double _dt, double _alpha ) { // problem specific settings elem = _elem; // mesh definition x0 = _x0; conn = _conn; dofs = _dofs; // time integration dt = _dt; alpha = _alpha; // extract mesh size from definition nnode = static_cast(x0 .rows()); ndim = static_cast(x0 .cols()); nelem = static_cast(conn.rows()); nne = static_cast(conn.cols()); ndof = dofs.maxCoeff()+1; // basic check (mostly the user is 'trusted') assert( dofs.size() == nnode * ndim ); assert( ndof < nnode * ndim ); // nodal quantities u.conservativeResize(nnode,ndim); u.setZero(); v.conservativeResize(nnode,ndim); v.setZero(); a.conservativeResize(nnode,ndim); a.setZero(); // linear system (DOFs) Minv.conservativeResize(ndof); Minv.setZero(); Fu .conservativeResize(ndof); Fu .setZero(); Fv .conservativeResize(ndof); Fv .setZero(); V .conservativeResize(ndof); V .setZero(); V_n .conservativeResize(ndof); V_n .setZero(); A .conservativeResize(ndof); A .setZero(); A_n .conservativeResize(ndof); A_n .setZero(); // compute inverse mass matrix : assumed constant in time computeMinv(); } // ================================================================================================= template void Periodic::velocityVerlet() { // (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; // - compute forces that are dependent on the displacement, but not on the velocity computeFu(); // (2a) estimate nodal velocities // - update velocities (DOFs) V.noalias() = V_n + dt * A; // - convert to nodal velocity (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 "Fu" and "Fv" correspond to this state of the system } // ================================================================================================= template 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; // - compute forces that are dependent on the displacement, but not on the velocity computeFu(); // - solve for accelerations (DOFs) A.noalias() = Minv.cwiseProduct( - Fu ); // - 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 "Fu" correspond to this state of the system } // ================================================================================================= template void Periodic::computeMinv() { // array with DOF numbers of each node cppmat::matrix2 dofe(nne,ndim); // zero-initialize mass matrix (only the inverse is stored) ColD M(ndof); M.setZero(); // loop over all elements for ( size_t e = 0 ; e < nelem ; ++e ) { // - nodal positions, displacement, and DOF numbers for element "e" for ( size_t m = 0 ; m < nne ; ++m ) { for ( size_t i = 0 ; i < ndim ; ++i ) { elem->xe(m,i) = x0 (conn(e,m),i); elem->ue(m,i) = u (conn(e,m),i); dofe (m,i) = dofs(conn(e,m),i); } } // - compute element mass matrix using problem specific "Element" class elem->computeM(e); // - assemble element mass matrix : take only the diagonal for ( size_t i = 0 ; i < nne*ndim ; ++i ) M(dofe[i]) += elem->M(i,i); // - check that the user provided a diagonal mass matrix #ifndef NDEBUG for ( size_t i = 0 ; i < nne*ndim ; ++i ) { for ( size_t j = 0 ; j < nne*ndim ; ++j ) { if ( i != j ) assert( ! elem->M(i,j) ); else assert( elem->M(i,i) ); } } #endif } // compute inverse of the mass matrix Minv = M.cwiseInverse(); } // ================================================================================================= template void Periodic::computeFu() { // array with DOF numbers of each node cppmat::matrix2 dofe(nne,ndim); // zero-initialize displacement dependent force Fu.setZero(); // loop over all elements for ( size_t e = 0 ; e < nelem ; ++e ) { // - nodal positions, displacement, and DOF numbers for element "e" for ( size_t m = 0 ; m < nne ; ++m ) { for ( size_t i = 0 ; i < ndim ; ++i ) { elem->xe(m,i) = x0 (conn(e,m),i); elem->ue(m,i) = u (conn(e,m),i); dofe (m,i) = dofs(conn(e,m),i); } } // - compute element force using problem specific "Element" class elem->computeFu(e); // - assemble force to global system for ( size_t i = 0 ; i < nne*ndim ; ++i ) Fu(dofe[i]) += elem->fu(i); } } // ================================================================================================= template void Periodic::computeFv() { // array with DOF numbers of each node cppmat::matrix2 dofe(nne,ndim); // zero-initialize velocity dependent force Fv.setZero(); // loop over all elements for ( size_t e = 0 ; e < nelem ; ++e ) { // - nodal positions, displacement, velocity, and DOF numbers for element "e" for ( size_t m = 0 ; m < nne ; ++m ) { for ( size_t i = 0 ; i < ndim ; ++i ) { elem->xe(m,i) = x0 (conn(e,m),i); elem->ue(m,i) = u (conn(e,m),i); elem->ve(m,i) = v (conn(e,m),i); dofe (m,i) = dofs(conn(e,m),i); } } // - compute element force using problem specific "Element" class elem->computeFv(e); // - assemble force to global system for ( size_t i = 0 ; i < nne*ndim ; ++i ) Fv(dofe[i]) += elem->fv(i); } } // ================================================================================================= template void Periodic::post() { // array with DOF numbers of each node cppmat::matrix2 dofe(nne,ndim); // loop over all elements for ( size_t e = 0 ; e < nelem ; ++e ) { // - nodal positions, displacement, velocity, and DOF numbers for element "e" for ( size_t m = 0 ; m < nne ; ++m ) { for ( size_t i = 0 ; i < ndim ; ++i ) { elem->xe(m,i) = x0 (conn(e,m),i); elem->ue(m,i) = u (conn(e,m),i); elem->ve(m,i) = v (conn(e,m),i); dofe (m,i) = dofs(conn(e,m),i); } } // - run post-process function of the problem specific "Element" class (output is handled there) elem->post(e); } } // ================================================================================================= }}} // namespace GooseFEM::Dynamics::Diagonal #endif diff --git a/src/GooseFEM/DynamicsDiagonalPeriodic.h b/src/GooseFEM/DynamicsDiagonalSemiPeriodic.h similarity index 78% copy from src/GooseFEM/DynamicsDiagonalPeriodic.h copy to src/GooseFEM/DynamicsDiagonalSemiPeriodic.h index 6924631..0708763 100644 --- a/src/GooseFEM/DynamicsDiagonalPeriodic.h +++ b/src/GooseFEM/DynamicsDiagonalSemiPeriodic.h @@ -1,351 +1,377 @@ /* ========================================== DESCRIPTION ========================================== (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ -#ifndef GOOSEFEM_DYNAMICS_DIAGONAL_PERIODIC_H -#define GOOSEFEM_DYNAMICS_DIAGONAL_PERIODIC_H +#ifndef GOOSEFEM_DYNAMICS_DIAGONAL_SEMIPERIODIC_H +#define GOOSEFEM_DYNAMICS_DIAGONAL_SEMIPERIODIC_H #include "Macros.h" #include // ------------------------------------------------------------------------------------------------- namespace GooseFEM { namespace Dynamics { namespace Diagonal { // =========================================== OVERVIEW ============================================ template -class Periodic +class SemiPeriodic { public: // variables // --------- // element/quadrature/material definition (should match mesh) std::shared_ptr elem; // mesh : nodal quantities and connectivity - MatS dofs; // DOF-numbers of each node (column of 'vectors') - MatS conn; // node numbers of each element (column of elements) - MatD x0; // initial positions of each node (column of vectors) - MatD u; // displacements of each node (column of vectors) - MatD v; // velocities of each node (column of vectors) - MatD a; // accelerations of each node (column of vectors) - size_t nnode; // number of nodes - size_t nelem; // number of elements - size_t nne; // number of nodes-per-element - size_t ndim; // number of dimensions - size_t ndof; // number of DOFs (after eliminating periodic dependencies) + MatS dofs; // DOF-numbers of each node (column of 'vectors') + MatS conn; // node numbers of each element (column of elements) + MatD x0; // initial positions of each node (column of vectors) + MatD u; // displacements of each node (column of vectors) + MatD v; // velocities of each node (column of vectors) + MatD a; // accelerations of each node (column of vectors) + size_t nnode; // number of nodes + size_t nelem; // number of elements + size_t nne; // number of nodes-per-element + size_t ndim; // number of dimensions + size_t ndof; // number of DOFs (after eliminating periodic dependencies) + + // fixed DOFs + size_t nfixed; // number of fixed DOFs + ColS fixedDofs; // DOF-number that are prescribed + ColD fixedV; // prescribed velocity for the fixed DOFs + ColD fixedA; // prescribed acceleration for the fixed DOFs // linear system - ColD Minv; // inverse of mass matrix (constructed diagonal -> inverse == diagonal) - ColD Fu; // force that depends on displacement (column) - ColD Fv; // force that depends on velocity (and displacement) (column) - ColD V; // velocity (column) - ColD V_n; // velocity (column, last increment) - ColD A; // acceleration (column) - ColD A_n; // acceleration (column, last increment) + ColD Minv; // inverse of mass matrix (constructed diagonal -> inverse == diagonal) + ColD Fu; // force that depends on displacement (column) + ColD Fv; // force that depends on velocity (and displacement) (column) + ColD V; // velocity (column) + ColD V_n; // velocity (column, last increment) + ColD A; // acceleration (column) + ColD A_n; // acceleration (column, last increment) // time integration - double dt; // time step - double t = 0.0; // current time - double alpha; // non-Galilean damping coefficient + double dt; // time step + double t = 0.0; // current time + double alpha; // non-Galilean damping coefficient // constructor // ----------- - Periodic(std::shared_ptr elem, const MatD &x0, const MatS &conn, const MatS &dofs, - double dt, double alpha=0.0 ); + SemiPeriodic(std::shared_ptr elem, const MatD &x0, const MatS &conn, const MatS &dofs, + const ColS &fixedDofs, double dt, double alpha=0.0 ); // functions // --------- void velocityVerlet(); // one time step of time integrator void Verlet(); // one time step of time integrator (Fv == 0) void computeMinv(); // element loop: inverse mass matrix <- "elem->computeM" void computeFu(); // element loop: displacement dependent forces <- "elem->computeFu" void computeFv(); // element loop: velocity dependent forces <- "elem->computeFv" void post(); // element loop: post-process <- "elem->post" }; // ========================================== SOURCE CODE ========================================== template -Periodic::Periodic( +SemiPeriodic::SemiPeriodic( std::shared_ptr _elem, const MatD &_x0, const MatS &_conn, const MatS &_dofs, - double _dt, double _alpha + const ColS &_fixedDofs, double _dt, double _alpha ) { // problem specific settings elem = _elem; // mesh definition x0 = _x0; conn = _conn; dofs = _dofs; // time integration dt = _dt; alpha = _alpha; + // fixed DOFs + fixedDofs = _fixedDofs; + nfixed = static_cast(fixedDofs.size()); + // extract mesh size from definition nnode = static_cast(x0 .rows()); ndim = static_cast(x0 .cols()); nelem = static_cast(conn.rows()); nne = static_cast(conn.cols()); ndof = dofs.maxCoeff()+1; // basic check (mostly the user is 'trusted') assert( dofs.size() == nnode * ndim ); assert( ndof < nnode * ndim ); // nodal quantities u.conservativeResize(nnode,ndim); u.setZero(); v.conservativeResize(nnode,ndim); v.setZero(); a.conservativeResize(nnode,ndim); a.setZero(); // linear system (DOFs) Minv.conservativeResize(ndof); Minv.setZero(); Fu .conservativeResize(ndof); Fu .setZero(); Fv .conservativeResize(ndof); Fv .setZero(); V .conservativeResize(ndof); V .setZero(); V_n .conservativeResize(ndof); V_n .setZero(); A .conservativeResize(ndof); A .setZero(); A_n .conservativeResize(ndof); A_n .setZero(); + // fixed DOFs : default zero velocity and acceleration + fixedV.conservativeResize(nfixed); fixedV.setZero(); + fixedA.conservativeResize(nfixed); fixedA.setZero(); + // compute inverse mass matrix : assumed constant in time computeMinv(); } // ================================================================================================= template -void Periodic::velocityVerlet() +void SemiPeriodic::velocityVerlet() { // (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; // - compute forces that are dependent on the displacement, but not on the velocity computeFu(); // (2a) estimate nodal velocities // - update velocities (DOFs) V.noalias() = V_n + dt * A; // - convert to nodal velocity (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 "Fu" and "Fv" correspond to this state of the system } // ================================================================================================= template -void Periodic::Verlet() +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; // - compute forces that are dependent on the displacement, but not on the velocity computeFu(); // - solve for accelerations (DOFs) A.noalias() = Minv.cwiseProduct( - Fu ); // - 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 "Fu" correspond to this state of the system } // ================================================================================================= template -void Periodic::computeMinv() +void SemiPeriodic::computeMinv() { // array with DOF numbers of each node cppmat::matrix2 dofe(nne,ndim); // zero-initialize mass matrix (only the inverse is stored) ColD M(ndof); M.setZero(); // loop over all elements for ( size_t e = 0 ; e < nelem ; ++e ) { // - nodal positions, displacement, and DOF numbers for element "e" for ( size_t m = 0 ; m < nne ; ++m ) { for ( size_t i = 0 ; i < ndim ; ++i ) { elem->xe(m,i) = x0 (conn(e,m),i); elem->ue(m,i) = u (conn(e,m),i); dofe (m,i) = dofs(conn(e,m),i); } } // - compute element mass matrix using problem specific "Element" class elem->computeM(e); // - assemble element mass matrix : take only the diagonal for ( size_t i = 0 ; i < nne*ndim ; ++i ) M(dofe[i]) += elem->M(i,i); // - check that the user provided a diagonal mass matrix #ifndef NDEBUG for ( size_t i = 0 ; i < nne*ndim ; ++i ) { for ( size_t j = 0 ; j < nne*ndim ; ++j ) { if ( i != j ) assert( ! elem->M(i,j) ); else assert( elem->M(i,i) ); } } #endif } // compute inverse of the mass matrix Minv = M.cwiseInverse(); } // ================================================================================================= template -void Periodic::computeFu() +void SemiPeriodic::computeFu() { // array with DOF numbers of each node cppmat::matrix2 dofe(nne,ndim); // zero-initialize displacement dependent force Fu.setZero(); // loop over all elements for ( size_t e = 0 ; e < nelem ; ++e ) { // - nodal positions, displacement, and DOF numbers for element "e" for ( size_t m = 0 ; m < nne ; ++m ) { for ( size_t i = 0 ; i < ndim ; ++i ) { elem->xe(m,i) = x0 (conn(e,m),i); elem->ue(m,i) = u (conn(e,m),i); dofe (m,i) = dofs(conn(e,m),i); } } // - compute element force using problem specific "Element" class elem->computeFu(e); // - assemble force to global system for ( size_t i = 0 ; i < nne*ndim ; ++i ) Fu(dofe[i]) += elem->fu(i); } } // ================================================================================================= template -void Periodic::computeFv() +void SemiPeriodic::computeFv() { // array with DOF numbers of each node cppmat::matrix2 dofe(nne,ndim); // zero-initialize velocity dependent force Fv.setZero(); // loop over all elements for ( size_t e = 0 ; e < nelem ; ++e ) { // - nodal positions, displacement, velocity, and DOF numbers for element "e" for ( size_t m = 0 ; m < nne ; ++m ) { for ( size_t i = 0 ; i < ndim ; ++i ) { elem->xe(m,i) = x0 (conn(e,m),i); elem->ue(m,i) = u (conn(e,m),i); elem->ve(m,i) = v (conn(e,m),i); dofe (m,i) = dofs(conn(e,m),i); } } // - compute element force using problem specific "Element" class elem->computeFv(e); // - assemble force to global system for ( size_t i = 0 ; i < nne*ndim ; ++i ) Fv(dofe[i]) += elem->fv(i); } } // ================================================================================================= template -void Periodic::post() +void SemiPeriodic::post() { // array with DOF numbers of each node cppmat::matrix2 dofe(nne,ndim); // loop over all elements for ( size_t e = 0 ; e < nelem ; ++e ) { // - nodal positions, displacement, velocity, and DOF numbers for element "e" for ( size_t m = 0 ; m < nne ; ++m ) { for ( size_t i = 0 ; i < ndim ; ++i ) { elem->xe(m,i) = x0 (conn(e,m),i); elem->ue(m,i) = u (conn(e,m),i); elem->ve(m,i) = v (conn(e,m),i); dofe (m,i) = dofs(conn(e,m),i); } } // - run post-process function of the problem specific "Element" class (output is handled there) elem->post(e); } } // ================================================================================================= }}} // namespace GooseFEM::Dynamics::Diagonal #endif diff --git a/src/GooseFEM/GooseFEM.h b/src/GooseFEM/GooseFEM.h index a7938a7..7f8fa9a 100644 --- a/src/GooseFEM/GooseFEM.h +++ b/src/GooseFEM/GooseFEM.h @@ -1,19 +1,20 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_H #define GOOSEFEM_H #include "Macros.h" #include "DynamicsDiagonalPeriodic.h" +#include "DynamicsDiagonalSemiPeriodic.h" #include "DynamicsDiagonalLinearStrainQuad4.h" #include "Mesh.h" #include "MeshTri3.h" #include "MeshQuad4.h" #endif