diff --git a/CMakeLists.txt b/CMakeLists.txt index 6c3be32..b4a3a8e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,137 +1,145 @@ #==================================================================================================# # # # (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM # # # #==================================================================================================# # initialization # -------------- # required to specify the c++ standard cmake_minimum_required(VERSION 3.0) # required for install include(CMakePackageConfigHelpers) include(GNUInstallDirs) # options option(PKGCONFIG "Build pkg-config" ON) option(BUILD_TESTS "${PROJECT_NAME} test suite" OFF) # project settings # ---------------- # name project(xGooseFEM) # file that contains the version information set(parse_version include/xGooseFEM/GooseFEM.h) # header files set(headers include/xGooseFEM/GooseFEM.h include/xGooseFEM/Mesh.hpp include/xGooseFEM/Mesh.h include/xGooseFEM/MeshTri3.hpp include/xGooseFEM/MeshTri3.h include/xGooseFEM/MeshQuad4.hpp include/xGooseFEM/MeshQuad4.h include/xGooseFEM/MeshHex8.hpp include/xGooseFEM/MeshHex8.h include/xGooseFEM/Element.hpp include/xGooseFEM/Element.h include/xGooseFEM/ElementQuad4.hpp include/xGooseFEM/ElementQuad4.h + include/xGooseFEM/ElementQuad4Planar.hpp + include/xGooseFEM/ElementQuad4Planar.h include/xGooseFEM/ElementHex8.hpp include/xGooseFEM/ElementHex8.h include/xGooseFEM/Vector.hpp include/xGooseFEM/Vector.h + include/xGooseFEM/VectorPartitioned.hpp + include/xGooseFEM/VectorPartitioned.h + include/xGooseFEM/MatrixPartitioned.hpp + include/xGooseFEM/MatrixPartitioned.h include/xGooseFEM/MatrixDiagonal.hpp include/xGooseFEM/MatrixDiagonal.h + include/xGooseFEM/MatrixDiagonalPartitioned.hpp + include/xGooseFEM/MatrixDiagonalPartitioned.h include/xGooseFEM/Iterate.hpp include/xGooseFEM/Iterate.h include/xGooseFEM/Dynamics.hpp include/xGooseFEM/Dynamics.h ) # automatically parse the version number file(READ "${parse_version}" version) string(REGEX MATCH "define[ \t]+GOOSEFEM_WORLD_VERSION[ \t]+([0-9]+)" _ "${version}") set(world_version "${CMAKE_MATCH_1}") string(REGEX MATCH "define[ \t]+GOOSEFEM_MAJOR_VERSION[ \t]+([0-9]+)" _ "${version}") set(major_version "${CMAKE_MATCH_1}") string(REGEX MATCH "define[ \t]+GOOSEFEM_MINOR_VERSION[ \t]+([0-9]+)" _ "${version}") set(minor_version "${CMAKE_MATCH_1}") set(GOOSEFEM_VERSION ${world_version}.${major_version}.${minor_version}) # print information to screen message(STATUS "Building ${PROJECT_NAME} v${GOOSEFEM_VERSION}") # paths # ----- set(GOOSEFEM_ROOT_DIR "${CMAKE_INSTALL_PREFIX}") set(GOOSEFEM_INCLUDE_DIR "${CMAKE_INSTALL_PREFIX}/${INCLUDE_INSTALL_DIR}") set(INCLUDE_INSTALL_DIR "${CMAKE_INSTALL_INCLUDEDIR}/${PROJECT_NAME}") set(CMAKEPACKAGE_INSTALL_DIR "${CMAKE_INSTALL_DATADIR}/${PROJECT_NAME}") set(PKGCONFIG_INSTALL_DIR "${CMAKE_INSTALL_DATADIR}/pkgconfig") set(fcmake "GooseFEMConfig.cmake") set(fpkg "GooseFEM.pc") # configure CMake # --------------- 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/build # ------------- # build test cases if(BUILD_TESTS) add_subdirectory(test) endif() # disable pkg-config for native Windows builds if(WIN32 OR CMAKE_HOST_SYSTEM_NAME MATCHES Windows) option(PKGCONFIG "Build pkg-config ${fpkg} file" OFF) endif() # pkg-config if(PKGCONFIG) configure_file(${fpkg}.in ${fpkg} @ONLY) install(FILES ${CMAKE_CURRENT_BINARY_DIR}/${fpkg} DESTINATION ${PKGCONFIG_INSTALL_DIR}) endif() # CMake install(FILES ${CMAKE_CURRENT_BINARY_DIR}/${fcmake} DESTINATION ${CMAKEPACKAGE_INSTALL_DIR}) # header files install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/${headers} DESTINATION ${INCLUDE_INSTALL_DIR}) # print information to screen # --------------------------- message(STATUS "") message(STATUS "+---------------------------------------------------------------------------------") message(STATUS "|") message(STATUS "| Use 'make install' to install in") message(STATUS "| ${CMAKE_INSTALL_PREFIX}") if(BUILD_TESTS) message(STATUS "|") message(STATUS "| Use 'make' and './test' to run the tests") endif() message(STATUS "|") message(STATUS "| To specify a custom directory call") message(STATUS "| cmake /path/to/${PROJECT_NAME} -DCMAKE_INSTALL_PREFIX=yourprefix") message(STATUS "|") 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/docs/examples/Statics/LinearElasticity.cpp b/docs/examples/Statics/LinearElasticity.cpp new file mode 100644 index 0000000..eae2885 --- /dev/null +++ b/docs/examples/Statics/LinearElasticity.cpp @@ -0,0 +1,71 @@ +#include +#include + +int main() +{ + // mesh + // ---- + + xGooseFEM::Mesh::Quad4::Regular mesh(5,5); + + xt::xtensor coor = mesh.coor(); + xt::xtensor conn = mesh.conn(); + xt::xtensor dofs = mesh.dofs(); + xt::xtensor disp = xt::zeros(coor.shape()); + + xt::xtensor nodesLeft = mesh.nodesLeftEdge(); + xt::xtensor nodesRight = mesh.nodesRightEdge(); + xt::xtensor nodesTop = mesh.nodesTopEdge(); + xt::xtensor nodesBottom = mesh.nodesBottomEdge(); + + // boundary displacements + // ---------------------- + + xt::xtensor iip = xt::empty({2*nodesLeft.size()+2*nodesBottom.size()}); + xt::xtensor u_p = xt::zeros(iip.shape()); + + size_t i = 0; + for ( auto &n : nodesRight ) { iip(i) = dofs(n,0); u_p(i) = 0.1; ++i; } + for ( auto &n : nodesTop ) { iip(i) = dofs(n,1); u_p(i) = -0.1; ++i; } + for ( auto &n : nodesLeft ) { iip(i) = dofs(n,0); u_p(i) = 0.0; ++i; } + for ( auto &n : nodesBottom ) { iip(i) = dofs(n,1); u_p(i) = 0.0; ++i; } + + // element definition + // ------------------ + + xGooseFEM::VectorPartitioned vector(conn, dofs, iip); + xGooseFEM::MatrixPartitioned K (conn, dofs, iip); + + xGooseFEM::Element::Quad4::QuadraturePlanar elem(vector.asElement(coor)); + + GMatLinearElastic::Cartesian3d::Material mat(mesh.nelem(), elem.nip(), 1., 1.); + + // solve + // ----- + + // strain + xt::xtensor Eps = elem.symGradN_vector(vector.asElement(disp)); + + // stress + xt::xtensor Sig = mat.Sig(Eps); + + // tangent + xt::xtensor C = mat.Tangent(); + + // internal force + xt::xtensor fint = vector.asNode(elem.int_gradN_dot_tensor2_dV(Sig)); + + // stiffness matrix + K.assemble(elem.int_gradN_dot_tensor4_dot_gradNT_dV(C)); + + // solve + xt::xtensor u_u = K.solve(vector.asDofs_u(fint), u_p); + + // assemble to nodal vector + disp = vector.asNode(u_u, u_p); + + // print result + std::cout << disp << std::endl; + + return 0; +} diff --git a/docs/examples/Statics/test b/docs/examples/Statics/test new file mode 100755 index 0000000..5bb7e18 Binary files /dev/null and b/docs/examples/Statics/test differ diff --git a/include/xGooseFEM/ElementQuad4.h b/include/xGooseFEM/ElementQuad4.h index 4311538..79ad9ad 100644 --- a/include/xGooseFEM/ElementQuad4.h +++ b/include/xGooseFEM/ElementQuad4.h @@ -1,161 +1,161 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_ELEMENTQUAD4_H #define XGOOSEFEM_ELEMENTQUAD4_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // ================================================================================================= namespace xGooseFEM { namespace Element { namespace Quad4 { // ================================================================================================= using T2 = xt::xtensor_fixed>; inline double inv(const T2 &A, T2 &Ainv); // ================================================================================================= namespace Gauss { inline size_t nip(); // number of integration points inline xt::xtensor xi(); // integration point coordinates (local coordinates) inline xt::xtensor w(); // integration point weights } // ================================================================================================= namespace Nodal { inline size_t nip(); // number of integration points inline xt::xtensor xi(); // integration point coordinates (local coordinates) inline xt::xtensor w(); // integration point weights } // ================================================================================================= -// ------------------------------------------ quadrature ------------------------------------------- +// ------------------------------------------------------------------------------------------------- class Quadrature { public: // fixed dimensions: // ndim = 2 - number of dimensions // nne = 4 - number of nodes per element // // naming convention: // "elemmat" - matrices stored per element - [nelem, nne*ndim, nne*ndim] // "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] // "qtensor" - integration point tensor - [nelem, nip, ndim, ndim] // "qscalar" - integration point scalar - [nelem, nip] // constructor: integration point coordinates and weights are optional (default: Gauss) Quadrature() = default; Quadrature(const xt::xtensor &x) : Quadrature(x, Gauss::xi(), Gauss::w()){}; Quadrature(const xt::xtensor &x, const xt::xtensor &xi, const xt::xtensor &w); // update the nodal positions (shape of "x" should match the earlier definition) void update_x(const xt::xtensor &x); // return dimensions size_t nelem() const { return m_nelem; }; // number of elements size_t nne() const { return m_nne; }; // number of nodes per element size_t ndim() const { return m_ndim; }; // number of dimension size_t nip() const { return m_nip; }; // number of integration points // return integration volume void dV(xt::xtensor &qscalar) const; void dV(xt::xtensor &qtensor) const; // same volume for all tensor components // dyadic product "qtensor(i,j) += dNdx(m,i) * elemvec(m,j)", its transpose and its symmetric part void gradN_vector(const xt::xtensor &elemvec, xt::xtensor &qtensor) const; void gradN_vector_T(const xt::xtensor &elemvec, xt::xtensor &qtensor) const; void symGradN_vector(const xt::xtensor &elemvec, xt::xtensor &qtensor) const; // integral of the scalar product "elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV" void int_N_scalar_NT_dV(const xt::xtensor &qscalar, xt::xtensor &elemmat) const; // integral of the dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" void int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor, xt::xtensor &elemvec) const; // integral of the dot product "elemmat(m*2+j, n*2+k) += dNdx(m,i) * qtensor(i,j,k,l) * dNdx(n,l) * dV" void int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor &qtensor, xt::xtensor &elemmat) const; // auto allocation of the functions above xt::xtensor dV() const; xt::xtensor dVtensor() const; xt::xtensor gradN_vector(const xt::xtensor &elemvec) const; xt::xtensor gradN_vector_T(const xt::xtensor &elemvec) const; xt::xtensor symGradN_vector(const xt::xtensor &elemvec) const; xt::xtensor int_N_scalar_NT_dV(const xt::xtensor &qscalar) const; xt::xtensor int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor) const; xt::xtensor int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor &qtensor) const; private: // compute "vol" and "dNdx" based on current "x" void compute_dN(); private: // dimensions (flexible) size_t m_nelem; // number of elements size_t m_nip; // number of integration points // dimensions (fixed for this element type) static const size_t m_nne=4; // number of nodes per element static const size_t m_ndim=2; // number of dimensions // data arrays xt::xtensor m_x; // nodal positions stored per element [nelem, nne, ndim] xt::xtensor m_w; // weight of each integration point [nip] xt::xtensor m_xi; // local coordinate of each integration point [nip, ndim] xt::xtensor m_N; // shape functions [nip, nne] xt::xtensor m_dNxi; // shape function gradients w.r.t. local coordinate [nip, nne, ndim] xt::xtensor m_dNx; // shape function gradients w.r.t. global coordinate [nelem, nip, nne, ndim] xt::xtensor m_vol; // integration point volume [nelem, nip] }; // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/ElementQuad4.h b/include/xGooseFEM/ElementQuad4Planar.h similarity index 77% copy from include/xGooseFEM/ElementQuad4.h copy to include/xGooseFEM/ElementQuad4Planar.h index 4311538..02c8868 100644 --- a/include/xGooseFEM/ElementQuad4.h +++ b/include/xGooseFEM/ElementQuad4Planar.h @@ -1,161 +1,142 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ -#ifndef XGOOSEFEM_ELEMENTQUAD4_H -#define XGOOSEFEM_ELEMENTQUAD4_H +#ifndef XGOOSEFEM_ELEMENTQUAD4PLANAR_H +#define XGOOSEFEM_ELEMENTQUAD4PLANAR_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // ================================================================================================= namespace xGooseFEM { namespace Element { namespace Quad4 { -// ================================================================================================= - -using T2 = xt::xtensor_fixed>; - -inline double inv(const T2 &A, T2 &Ainv); - -// ================================================================================================= - -namespace Gauss { -inline size_t nip(); // number of integration points -inline xt::xtensor xi(); // integration point coordinates (local coordinates) -inline xt::xtensor w(); // integration point weights -} - -// ================================================================================================= - -namespace Nodal { -inline size_t nip(); // number of integration points -inline xt::xtensor xi(); // integration point coordinates (local coordinates) -inline xt::xtensor w(); // integration point weights -} - -// ================================================================================================= - -// ------------------------------------------ quadrature ------------------------------------------- +// ------------------------------------------------------------------------------------------------- -class Quadrature +class QuadraturePlanar { public: // fixed dimensions: // ndim = 2 - number of dimensions // nne = 4 - number of nodes per element + // tdim = 3 - number of dimensions of tensors // // naming convention: // "elemmat" - matrices stored per element - [nelem, nne*ndim, nne*ndim] // "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] - // "qtensor" - integration point tensor - [nelem, nip, ndim, ndim] + // "qtensor" - integration point tensor - [nelem, nip, tdim, tdim] // "qscalar" - integration point scalar - [nelem, nip] // constructor: integration point coordinates and weights are optional (default: Gauss) - Quadrature() = default; + QuadraturePlanar() = default; - Quadrature(const xt::xtensor &x) : Quadrature(x, Gauss::xi(), Gauss::w()){}; + QuadraturePlanar(const xt::xtensor &x, double thick=1.) : QuadraturePlanar(x, Gauss::xi(), Gauss::w()){}; - Quadrature(const xt::xtensor &x, const xt::xtensor &xi, const xt::xtensor &w); + QuadraturePlanar(const xt::xtensor &x, const xt::xtensor &xi, const xt::xtensor &w, double thick=1.); // update the nodal positions (shape of "x" should match the earlier definition) void update_x(const xt::xtensor &x); // return dimensions size_t nelem() const { return m_nelem; }; // number of elements size_t nne() const { return m_nne; }; // number of nodes per element size_t ndim() const { return m_ndim; }; // number of dimension size_t nip() const { return m_nip; }; // number of integration points // return integration volume void dV(xt::xtensor &qscalar) const; void dV(xt::xtensor &qtensor) const; // same volume for all tensor components // dyadic product "qtensor(i,j) += dNdx(m,i) * elemvec(m,j)", its transpose and its symmetric part void gradN_vector(const xt::xtensor &elemvec, xt::xtensor &qtensor) const; void gradN_vector_T(const xt::xtensor &elemvec, xt::xtensor &qtensor) const; void symGradN_vector(const xt::xtensor &elemvec, xt::xtensor &qtensor) const; // integral of the scalar product "elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV" void int_N_scalar_NT_dV(const xt::xtensor &qscalar, xt::xtensor &elemmat) const; // integral of the dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" void int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor, xt::xtensor &elemvec) const; // integral of the dot product "elemmat(m*2+j, n*2+k) += dNdx(m,i) * qtensor(i,j,k,l) * dNdx(n,l) * dV" void int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor &qtensor, xt::xtensor &elemmat) const; // auto allocation of the functions above xt::xtensor dV() const; xt::xtensor dVtensor() const; xt::xtensor gradN_vector(const xt::xtensor &elemvec) const; xt::xtensor gradN_vector_T(const xt::xtensor &elemvec) const; xt::xtensor symGradN_vector(const xt::xtensor &elemvec) const; xt::xtensor int_N_scalar_NT_dV(const xt::xtensor &qscalar) const; xt::xtensor int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor) const; xt::xtensor int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor &qtensor) const; private: // compute "vol" and "dNdx" based on current "x" void compute_dN(); private: // dimensions (flexible) size_t m_nelem; // number of elements size_t m_nip; // number of integration points // dimensions (fixed for this element type) static const size_t m_nne=4; // number of nodes per element static const size_t m_ndim=2; // number of dimensions + static const size_t m_tdim=3; // number of dimensions of tensors // data arrays xt::xtensor m_x; // nodal positions stored per element [nelem, nne, ndim] xt::xtensor m_w; // weight of each integration point [nip] xt::xtensor m_xi; // local coordinate of each integration point [nip, ndim] xt::xtensor m_N; // shape functions [nip, nne] xt::xtensor m_dNxi; // shape function gradients w.r.t. local coordinate [nip, nne, ndim] xt::xtensor m_dNx; // shape function gradients w.r.t. global coordinate [nelem, nip, nne, ndim] xt::xtensor m_vol; // integration point volume [nelem, nip] + // thickness + double m_thick; + }; // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif diff --git a/include/xGooseFEM/ElementQuad4Planar.hpp b/include/xGooseFEM/ElementQuad4Planar.hpp new file mode 100644 index 0000000..9e0dbbc --- /dev/null +++ b/include/xGooseFEM/ElementQuad4Planar.hpp @@ -0,0 +1,514 @@ +/* ================================================================================================= + +(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM + +================================================================================================= */ + +#ifndef XGOOSEFEM_ELEMENTQUAD4PLANAR_CPP +#define XGOOSEFEM_ELEMENTQUAD4PLANAR_CPP + +// ------------------------------------------------------------------------------------------------- + +#include "GooseFEM.h" + +// ================================================================================================= + +namespace xGooseFEM { +namespace Element { +namespace Quad4 { + +// ================================================================================================= + +inline QuadraturePlanar::QuadraturePlanar(const xt::xtensor &x, + const xt::xtensor &xi, const xt::xtensor &w, double thick) : + m_x(x), m_w(w), m_xi(xi), m_thick(thick) +{ + // check input + assert( m_x.shape()[1] == m_nne ); + assert( m_x.shape()[2] == m_ndim ); + + // extract number of elements and number of integration points + m_nelem = m_x.shape()[0]; + m_nip = m_w.size(); + + // check input + assert( m_xi.shape()[0] == m_nip ); + assert( m_xi.shape()[1] == m_ndim ); + assert( m_w .size() == m_nip ); + + // allocate arrays + m_N = xt::empty({ m_nip, m_nne }); + m_dNxi = xt::empty({ m_nip, m_nne, m_ndim}); + m_dNx = xt::empty({m_nelem, m_nip, m_nne, m_ndim}); + m_vol = xt::empty({m_nelem, m_nip }); + + // shape functions + for ( size_t q = 0 ; q < m_nip ; ++q ) + { + m_N(q,0) = .25 * (1.-m_xi(q,0)) * (1.-m_xi(q,1)); + m_N(q,1) = .25 * (1.+m_xi(q,0)) * (1.-m_xi(q,1)); + m_N(q,2) = .25 * (1.+m_xi(q,0)) * (1.+m_xi(q,1)); + m_N(q,3) = .25 * (1.-m_xi(q,0)) * (1.+m_xi(q,1)); + } + + // shape function gradients in local coordinates + for ( size_t q = 0 ; q < m_nip ; ++q ) + { + // - dN / dxi_0 + m_dNxi(q,0,0) = -.25*(1.-m_xi(q,1)); + m_dNxi(q,1,0) = +.25*(1.-m_xi(q,1)); + m_dNxi(q,2,0) = +.25*(1.+m_xi(q,1)); + m_dNxi(q,3,0) = -.25*(1.+m_xi(q,1)); + // - dN / dxi_1 + m_dNxi(q,0,1) = -.25*(1.-m_xi(q,0)); + m_dNxi(q,1,1) = -.25*(1.+m_xi(q,0)); + m_dNxi(q,2,1) = +.25*(1.+m_xi(q,0)); + m_dNxi(q,3,1) = +.25*(1.-m_xi(q,0)); + } + + // compute the shape function gradients, based on "x" + compute_dN(); +} + +// ------------------------------------------------------------------------------------------------- + +inline void QuadraturePlanar::dV(xt::xtensor &qscalar) const +{ + assert( qscalar.shape()[0] == m_nelem ); + assert( qscalar.shape()[1] == m_nip ); + + #pragma omp parallel for + for ( size_t e = 0 ; e < m_nelem ; ++e ) + for ( size_t q = 0 ; q < m_nip ; ++q ) + qscalar(e,q) = m_vol(e,q); +} + +// ------------------------------------------------------------------------------------------------- + +inline void QuadraturePlanar::dV(xt::xtensor &qtensor) const +{ + assert( qtensor.shape()[0] == m_nelem ); + assert( qtensor.shape()[1] == m_nne ); + assert( qtensor.shape()[2] == m_tdim ); + assert( qtensor.shape()[3] == m_tdim ); + + #pragma omp parallel for + for ( size_t e = 0 ; e < m_nelem ; ++e ) + for ( size_t q = 0 ; q < m_nip ; ++q ) + for ( size_t i = 0 ; i < m_tdim ; ++i ) + for ( size_t j = 0 ; j < m_tdim ; ++j ) + qtensor(e,q,i,j) = m_vol(e,q); +} + +// ------------------------------------------------------------------------------------------------- + +inline void QuadraturePlanar::update_x(const xt::xtensor &x) +{ + assert( x.shape()[0] == m_nelem ); + assert( x.shape()[1] == m_nne ); + assert( x.shape()[2] == m_ndim ); + assert( x.size() == m_x.size() ); + + // update positions + xt::noalias(m_x) = x; + + // update the shape function gradients for the new "x" + compute_dN(); +} + +// ------------------------------------------------------------------------------------------------- + +inline void QuadraturePlanar::compute_dN() +{ + // loop over all elements (in parallel) + #pragma omp parallel + { + // allocate local variables + T2 J, Jinv; + + #pragma omp for + for ( size_t e = 0 ; e < m_nelem ; ++e ) + { + // alias nodal positions + auto x = xt::adapt(&m_x(e,0,0), xt::xshape()); + + // loop over integration points + for ( size_t q = 0 ; q < m_nip ; ++q ) + { + // - alias + auto dNxi = xt::adapt(&m_dNxi( q,0,0), xt::xshape()); + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + + // - Jacobian (loops unrolled for efficiency) + // J(i,j) += dNxi(m,i) * x(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 + double Jdet = inv(J, Jinv); + + // - shape function gradients wrt global coordinates (loops partly unrolled for efficiency) + // dNx(m,i) += Jinv(i,j) * dNxi(m,j); + for ( size_t m = 0 ; m < 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); + } + + // - integration point volume + m_vol(e,q) = m_w(q) * Jdet * m_thick; + } + } + } +} + +// ------------------------------------------------------------------------------------------------- + +inline void QuadraturePlanar::gradN_vector( + const xt::xtensor &elemvec, xt::xtensor &qtensor) const +{ + assert( elemvec.shape()[0] == m_nelem ); + assert( elemvec.shape()[1] == m_nne ); + assert( elemvec.shape()[2] == m_ndim ); + assert( qtensor.shape()[0] == m_nelem ); + assert( qtensor.shape()[1] == m_nne ); + assert( qtensor.shape()[2] == m_tdim ); + assert( qtensor.shape()[3] == m_tdim ); + + // zero-initialize output: matrix of tensors + qtensor.fill(0.0); + + // loop over all elements (in parallel) + #pragma omp parallel for + for ( size_t e = 0 ; e < m_nelem ; ++e ) + { + // alias element vector (e.g. nodal displacements) + auto u = xt::adapt(&elemvec(e,0,0), xt::xshape()); + + // loop over all integration points in element "e" + for ( size_t q = 0 ; q < m_nip ; ++q ) + { + // - alias + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto gradu = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); + + // - evaluate dyadic product (loops unrolled for efficiency) + // gradu(i,j) += dNx(m,i) * u(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); + } + } +} + +// ------------------------------------------------------------------------------------------------- + +inline void QuadraturePlanar::gradN_vector_T( + const xt::xtensor &elemvec, xt::xtensor &qtensor) const +{ + assert( elemvec.shape()[0] == m_nelem ); + assert( elemvec.shape()[1] == m_nne ); + assert( elemvec.shape()[2] == m_ndim ); + assert( qtensor.shape()[0] == m_nelem ); + assert( qtensor.shape()[1] == m_nne ); + assert( qtensor.shape()[2] == m_tdim ); + assert( qtensor.shape()[3] == m_tdim ); + + // zero-initialize output: matrix of tensors + qtensor.fill(0.0); + + // loop over all elements (in parallel) + #pragma omp parallel for + for ( size_t e = 0 ; e < m_nelem ; ++e ) + { + // alias element vector (e.g. nodal displacements) + auto u = xt::adapt(&elemvec(e,0,0), xt::xshape()); + + // loop over all integration points in element "e" + for ( size_t q = 0 ; q < m_nip ; ++q ) + { + // - alias + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto gradu = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); + + // - evaluate transpose of dyadic product (loops unrolled for efficiency) + // gradu(j,i) += dNx(m,i) * u(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(1,0) = 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(0,1) = 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); + } + } +} + +// ------------------------------------------------------------------------------------------------- + +inline void QuadraturePlanar::symGradN_vector( + const xt::xtensor &elemvec, xt::xtensor &qtensor) const +{ + assert( elemvec.shape()[0] == m_nelem ); + assert( elemvec.shape()[1] == m_nne ); + assert( elemvec.shape()[2] == m_ndim ); + assert( qtensor.shape()[0] == m_nelem ); + assert( qtensor.shape()[1] == m_nne ); + assert( qtensor.shape()[2] == m_tdim ); + assert( qtensor.shape()[3] == m_tdim ); + + // zero-initialize output: matrix of tensors + qtensor.fill(0.0); + + // loop over all elements (in parallel) + #pragma omp parallel for + for ( size_t e = 0 ; e < m_nelem ; ++e ) + { + // alias element vector (e.g. nodal displacements) + auto u = xt::adapt(&elemvec(e,0,0), xt::xshape()); + + // loop over all integration points in element "e" + for ( size_t q = 0 ; q < m_nip ; ++q ) + { + // - alias + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto eps = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); + + // - evaluate symmetrized dyadic product (loops unrolled for efficiency) + // grad(i,j) += dNx(m,i) * u(m,j) + // eps (j,i) = 0.5 * ( grad(i,j) + grad(j,i) ) + eps(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); + eps(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); + eps(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) + + dNx(0,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,0) ) * 0.5; + eps(1,0) = eps(0,1); + } + } +} + +// ------------------------------------------------------------------------------------------------- + +inline void QuadraturePlanar::int_N_scalar_NT_dV( + const xt::xtensor &qscalar, xt::xtensor &elemmat) const +{ + assert( qscalar.shape()[0] == m_nelem ); + assert( qscalar.shape()[1] == m_nip ); + assert( elemmat.shape()[0] == m_nelem ); + assert( elemmat.shape()[1] == m_nne*m_ndim ); + assert( elemmat.shape()[2] == m_nne*m_ndim ); + + // zero-initialize: matrix of matrices + elemmat.fill(0.0); + + // loop over all elements (in parallel) + #pragma omp parallel for + for ( size_t e = 0 ; e < m_nelem ; ++e ) + { + // alias (e.g. mass matrix) + auto M = xt::adapt(&elemmat(e,0,0), xt::xshape()); + + // loop over all integration points in element "e" + for ( size_t q = 0 ; q < m_nip ; ++q ) + { + // - alias + auto N = xt::adapt(&m_N(q,0), xt::xshape()); + auto& vol = m_vol (e,q); + auto& rho = qscalar(e,q); + + // - evaluate scalar product, for all dimensions, and assemble + // M(m*ndim+i,n*ndim+i) += N(m) * scalar * N(n) * dV + for ( size_t m = 0 ; m < m_nne ; ++m ) { + for ( size_t n = 0 ; n < m_nne ; ++n ) { + M(m*m_ndim+0, n*m_ndim+0) += N(m) * rho * N(n) * vol; + M(m*m_ndim+1, n*m_ndim+1) += N(m) * rho * N(n) * vol; + } + } + } + } +} + +// ------------------------------------------------------------------------------------------------- + +inline void QuadraturePlanar::int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor, + xt::xtensor &elemvec) const +{ + assert( qtensor.shape()[0] == m_nelem ); + assert( qtensor.shape()[1] == m_nip ); + assert( qtensor.shape()[2] == m_tdim ); + assert( qtensor.shape()[3] == m_tdim ); + assert( elemvec.shape()[0] == m_nelem ); + assert( elemvec.shape()[1] == m_nne ); + assert( elemvec.shape()[2] == m_ndim ); + + // zero-initialize output: matrix of vectors + elemvec.fill(0.0); + + // loop over all elements (in parallel) + #pragma omp parallel for + for ( size_t e = 0 ; e < m_nelem ; ++e ) + { + // alias (e.g. nodal force) + auto f = xt::adapt(&elemvec(e,0,0), xt::xshape()); + + // loop over all integration points in element "e" + for ( size_t q = 0 ; q < m_nip ; ++q ) + { + // - alias + auto dNx = xt::adapt(&m_dNx (e,q,0,0), xt::xshape()); + auto sig = xt::adapt(&qtensor(e,q,0,0), xt::xshape()); + auto& vol = m_vol(e,q); + + // - evaluate dot product, and assemble + for ( size_t m = 0 ; m < m_nne ; ++m ) + { + f(m,0) += ( dNx(m,0) * sig(0,0) + dNx(m,1) * sig(1,0) ) * vol; + f(m,1) += ( dNx(m,0) * sig(0,1) + dNx(m,1) * sig(1,1) ) * vol; + } + } + } +} + +// ------------------------------------------------------------------------------------------------- + +inline void QuadraturePlanar::int_gradN_dot_tensor4_dot_gradNT_dV(const xt::xtensor &qtensor, + xt::xtensor &elemmat) const +{ + assert( qtensor.shape()[0] == m_nelem ); + assert( qtensor.shape()[1] == m_nip ); + assert( qtensor.shape()[2] == m_tdim ); + assert( qtensor.shape()[3] == m_tdim ); + assert( qtensor.shape()[4] == m_tdim ); + assert( qtensor.shape()[5] == m_tdim ); + + assert( elemmat.shape()[0] == m_nelem ); + assert( elemmat.shape()[1] == m_nne*m_ndim ); + assert( elemmat.shape()[2] == m_nne*m_ndim ); + + // zero-initialize output: matrix of vector + elemmat.fill(0.0); + + // loop over all elements (in parallel) + #pragma omp parallel for + for ( size_t e = 0 ; e < m_nelem ; ++e ) + { + // alias (e.g. nodal force) + auto K = xt::adapt(&elemmat(e,0,0), xt::xshape()); + + // loop over all integration points in element "e" + for ( size_t q = 0 ; q < m_nip ; ++q ){ + + // - alias + auto dNx = xt::adapt(&m_dNx(e,q,0,0), xt::xshape()); + auto C = xt::adapt(&qtensor(e,q,0,0,0,0), xt::xshape()); + auto& vol = m_vol(e,q); + + // - evaluate dot product, and assemble + for ( size_t m = 0 ; m < m_nne ; ++m ) + for ( size_t n = 0 ; n < m_nne ; ++n ) + for ( size_t i = 0 ; i < m_ndim ; ++i ) + for ( size_t j = 0 ; j < m_ndim ; ++j ) + for ( size_t k = 0 ; k < m_ndim ; ++k ) + for ( size_t l = 0 ; l < m_ndim ; ++l ) + K(m*m_ndim+j, n*m_ndim+k) += dNx(m,i) * C(i,j,k,l) * dNx(n,l) * vol; + } + } +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor QuadraturePlanar::dV() const +{ + xt::xtensor out = xt::empty({m_nelem, m_nip}); + + this->dV(out); + + return out; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor QuadraturePlanar::dVtensor() const +{ + xt::xtensor out = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); + + this->dV(out); + + return out; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor QuadraturePlanar::gradN_vector(const xt::xtensor &elemvec) const +{ + xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); + + this->gradN_vector(elemvec, qtensor); + + return qtensor; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor QuadraturePlanar::gradN_vector_T(const xt::xtensor &elemvec) const +{ + xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); + + this->gradN_vector_T(elemvec, qtensor); + + return qtensor; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor QuadraturePlanar::symGradN_vector(const xt::xtensor &elemvec) const +{ + xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_tdim, m_tdim}); + + this->symGradN_vector(elemvec, qtensor); + + return qtensor; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor QuadraturePlanar::int_N_scalar_NT_dV( + const xt::xtensor &qscalar) const +{ + xt::xtensor elemmat = xt::empty({m_nelem, m_nne*m_ndim, m_nne*m_ndim}); + + this->int_N_scalar_NT_dV(qscalar, elemmat); + + return elemmat; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor QuadraturePlanar::int_gradN_dot_tensor2_dV( + const xt::xtensor &qtensor) const +{ + xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); + + this->int_gradN_dot_tensor2_dV(qtensor, elemvec); + + return elemvec; +} + +// ------------------------------------------------------------------------------------------------- + +inline xt::xtensor QuadraturePlanar::int_gradN_dot_tensor4_dot_gradNT_dV( + const xt::xtensor &qtensor) const + { + xt::xtensor elemmat = xt::empty({m_nelem, m_ndim*m_nne, m_ndim*m_nne}); + + this->int_gradN_dot_tensor4_dot_gradNT_dV(qtensor, elemmat); + + return elemmat; + } + +// ------------------------------------------------------------------------------------------------- + +}}} // namespace ... + +// ================================================================================================= + +#endif diff --git a/include/xGooseFEM/GooseFEM.h b/include/xGooseFEM/GooseFEM.h index d885316..4c1f4ee 100644 --- a/include/xGooseFEM/GooseFEM.h +++ b/include/xGooseFEM/GooseFEM.h @@ -1,111 +1,112 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef XGOOSEFEM_H #define XGOOSEFEM_H // ================================================================================================= #define _USE_MATH_DEFINES // to use "M_PI" from "math.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace xt::placeholders; // ================================================================================================= #define GOOSEFEM_WORLD_VERSION 0 #define GOOSEFEM_MAJOR_VERSION 1 #define GOOSEFEM_MINOR_VERSION 0 #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)))) + 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 sparse matrices namespace xGooseFEM { typedef Eigen::SparseMatrix SpMatD; typedef Eigen::SparseMatrix SpMatS; typedef Eigen::Triplet TripD; } // ================================================================================================= #include "Mesh.h" -#include "MeshTri3.h" -#include "MeshQuad4.h" -#include "MeshHex8.h" -#include "Element.h" -#include "ElementQuad4.h" -#include "ElementHex8.h" -#include "Vector.h" -#include "VectorPartitioned.h" -#include "MatrixPartitioned.h" -#include "MatrixDiagonal.h" -#include "MatrixDiagonalPartitioned.h" -#include "Iterate.h" -#include "Dynamics.h" - #include "Mesh.hpp" +#include "MeshTri3.h" #include "MeshTri3.hpp" +#include "MeshQuad4.h" #include "MeshQuad4.hpp" +#include "MeshHex8.h" #include "MeshHex8.hpp" +#include "Element.h" #include "Element.hpp" +#include "ElementQuad4.h" #include "ElementQuad4.hpp" +#include "ElementQuad4Planar.h" +#include "ElementQuad4Planar.hpp" +#include "ElementHex8.h" #include "ElementHex8.hpp" +#include "Vector.h" #include "Vector.hpp" +#include "VectorPartitioned.h" #include "VectorPartitioned.hpp" +#include "MatrixPartitioned.h" #include "MatrixPartitioned.hpp" +#include "MatrixDiagonal.h" #include "MatrixDiagonal.hpp" +#include "MatrixDiagonalPartitioned.h" #include "MatrixDiagonalPartitioned.hpp" +#include "Iterate.h" #include "Iterate.hpp" +#include "Dynamics.h" #include "Dynamics.hpp" // ================================================================================================= #endif