diff --git a/setup.py b/setup.py index 84e4b30..ed32092 100644 --- a/setup.py +++ b/setup.py @@ -1,49 +1,51 @@ desc = ''' GooseFEM is a C++ module, wrapped in Python, that provides several predefined finite element meshes. The original C++ module also includes element definitions and several standard finite element simulations. ''' from setuptools import setup, Extension import sys, re import setuptools import pybind11 -import cppmat +import pyxtensor header = open('src/GooseFEM/GooseFEM.h','r').read() world = re.split(r'(.*)(\#define GOOSEFEM_WORLD_VERSION\ )([0-9]+)(.*)',header)[3] major = re.split(r'(.*)(\#define GOOSEFEM_MAJOR_VERSION\ )([0-9]+)(.*)',header)[3] minor = re.split(r'(.*)(\#define GOOSEFEM_MINOR_VERSION\ )([0-9]+)(.*)',header)[3] __version__ = '.'.join([world,major,minor]) ext_modules = [ Extension( 'GooseFEM', ['src/GooseFEM/python.cpp'], include_dirs=[ - pybind11.get_include(False), - pybind11.get_include(True ), - cppmat .get_include(False), - cppmat .get_include(True ), - cppmat .find_eigen() + pybind11 .get_include(False), + pybind11 .get_include(True ), + pyxtensor.get_include(False), + pyxtensor.get_include(True ), + pyxtensor.find_xtensor(), + pyxtensor.find_xtl(), + pyxtensor.find_eigen(), ], language='c++' ), ] setup( name = 'GooseFEM', description = 'Finite element meshes, quadrature, and assembly tools', long_description = desc, version = __version__, license = 'GPLv3', author = 'Tom de Geus', author_email = 'tom@geus.me', url = 'https://github.com/tdegeus/GooseFEM', ext_modules = ext_modules, - install_requires = ['pybind11>=2.2.0','cppmat>=1.0.6'], - cmdclass = {'build_ext': cppmat.BuildExt}, + install_requires = ['pybind11>=2.2.0','pyxtensor>=0.0.1'], + cmdclass = {'build_ext': pyxtensor.BuildExt}, zip_safe = False, ) diff --git a/src/GooseFEM/ElementHex8.h b/src/GooseFEM/ElementHex8.h index 203e682..34b2a30 100644 --- a/src/GooseFEM/ElementHex8.h +++ b/src/GooseFEM/ElementHex8.h @@ -1,134 +1,133 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_ELEMENTHEX8_H #define GOOSEFEM_ELEMENTHEX8_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // ==================================== GooseFEM::Element::Hex8 ==================================== namespace GooseFEM { namespace Element { namespace Hex8 { // ======================================== tensor algebra ========================================= static const size_t ndim = 3; using T2 = xt::xtensor_fixed>; -inline double det(const T2 &A); -inline T2 inv(const T2 &A); +inline double inv(const T2 &A, T2 &Ainv); // ================================ GooseFEM::Element::Hex8::Gauss ================================= 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 } // ================================ GooseFEM::Element::Hex8::Nodal ================================= 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 { 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=8; // number of nodes per element static const size_t m_ndim=3; // 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] private: // compute "vol" and "dNdx" based on current "x" void compute_dN(); public: // convention: // "elemmat" - matrices stored per element - ArrD - [nelem, nne*ndim, nne*ndim] // "elemvec" - nodal vectors stored per element - ArrD - [nelem, nne, ndim] // "qtensor" - integration point tensor - ArrD - [nelem, nip, ndim, ndim] // "qscalar" - integration point scalar - ArrD - [nelem, nip] // constructor: integration point coordinates and weights are optional (default: Gauss) Quadrature() = default; Quadrature(const xt::xtensor &x); 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; // number of elements size_t nne() const; // number of nodes per element size_t ndim() const; // number of dimension size_t nip() const; // number of integration points // return integration volume // - in-place void dV(xt::xtensor &qscalar) const; void dV(xt::xtensor &qtensor) const; // same volume for all tensor components // - return qscalar xt::xtensor dV() const; // dyadic product "qtensor(i,j) += dNdx(m,i) * elemvec(m,j)", its transpose and its symmetric part // - in-place 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; // - return qtensor 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; // integral of the scalar product "elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV" // - in-place void int_N_scalar_NT_dV(const xt::xtensor &qscalar, xt::xtensor &elemmat) const; // - return elemmat xt::xtensor int_N_scalar_NT_dV(const xt::xtensor &qscalar) const; // integral of the dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" // - in-place void int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor, xt::xtensor &elemvec) const; // - return elemvec xt::xtensor int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor) const; }; // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/ElementHex8.hpp b/src/GooseFEM/ElementHex8.hpp index f663e88..4a209bb 100644 --- a/src/GooseFEM/ElementHex8.hpp +++ b/src/GooseFEM/ElementHex8.hpp @@ -1,719 +1,714 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_ELEMENTHEX8_CPP #define GOOSEFEM_ELEMENTHEX8_CPP // ------------------------------------------------------------------------------------------------- #include "ElementHex8.h" // ==================================== GooseFEM::Element::Hex8 ==================================== namespace GooseFEM { namespace Element { namespace Hex8 { // ======================================== tensor algebra ========================================= -inline double det(const T2 &A) -{ - return ( A[0] * A[4] * A[8] + - A[1] * A[5] * A[6] + - A[2] * A[3] * A[7] ) - - ( A[2] * A[4] * A[6] + - A[1] * A[3] * A[8] + - A[0] * A[5] * A[7] ); -} - -// ------------------------------------------------------------------------------------------------- - -inline double inv(const T2 &A, T2 &C) +inline double inv(const T2 &A, T2 &Ainv) { // compute determinant double det = ( A[0] * A[4] * A[8] + A[1] * A[5] * A[6] + A[2] * A[3] * A[7] ) - ( A[2] * A[4] * A[6] + A[1] * A[3] * A[8] + A[0] * A[5] * A[7] ); // compute inverse - C[0] = (A[4]*A[8]-A[5]*A[7]) / det; - C[1] = (A[2]*A[7]-A[1]*A[8]) / det; - C[2] = (A[1]*A[5]-A[2]*A[4]) / det; - C[3] = (A[5]*A[6]-A[3]*A[8]) / det; - C[4] = (A[0]*A[8]-A[2]*A[6]) / det; - C[5] = (A[2]*A[3]-A[0]*A[5]) / det; - C[6] = (A[3]*A[7]-A[4]*A[6]) / det; - C[7] = (A[1]*A[6]-A[0]*A[7]) / det; - C[8] = (A[0]*A[4]-A[1]*A[3]) / det; + Ainv[0] = (A[4]*A[8]-A[5]*A[7]) / det; + Ainv[1] = (A[2]*A[7]-A[1]*A[8]) / det; + Ainv[2] = (A[1]*A[5]-A[2]*A[4]) / det; + Ainv[3] = (A[5]*A[6]-A[3]*A[8]) / det; + Ainv[4] = (A[0]*A[8]-A[2]*A[6]) / det; + Ainv[5] = (A[2]*A[3]-A[0]*A[5]) / det; + Ainv[6] = (A[3]*A[7]-A[4]*A[6]) / det; + Ainv[7] = (A[1]*A[6]-A[0]*A[7]) / det; + Ainv[8] = (A[0]*A[4]-A[1]*A[3]) / det; return det; } // ================================ GooseFEM::Element::Hex8::Gauss ================================= namespace Gauss { // --------------------------------- number of integration points ---------------------------------- inline size_t nip() { return 8; } // ----------------------- integration point coordinates (local coordinates) ----------------------- inline xt::xtensor xi() { static const size_t nip = 8; static const size_t ndim = 3; xt::xtensor xi = xt::empty({nip,ndim}); xi(0,0) = -1./std::sqrt(3.); xi(0,1) = -1./std::sqrt(3.); xi(0,2) = -1./std::sqrt(3.); xi(1,0) = +1./std::sqrt(3.); xi(1,1) = -1./std::sqrt(3.); xi(1,2) = -1./std::sqrt(3.); xi(2,0) = +1./std::sqrt(3.); xi(2,1) = +1./std::sqrt(3.); xi(2,2) = -1./std::sqrt(3.); xi(3,0) = -1./std::sqrt(3.); xi(3,1) = +1./std::sqrt(3.); xi(3,2) = -1./std::sqrt(3.); xi(4,0) = -1./std::sqrt(3.); xi(4,1) = -1./std::sqrt(3.); xi(4,2) = +1./std::sqrt(3.); xi(5,0) = +1./std::sqrt(3.); xi(5,1) = -1./std::sqrt(3.); xi(5,2) = +1./std::sqrt(3.); xi(6,0) = +1./std::sqrt(3.); xi(6,1) = +1./std::sqrt(3.); xi(6,2) = +1./std::sqrt(3.); xi(7,0) = -1./std::sqrt(3.); xi(7,1) = +1./std::sqrt(3.); xi(7,2) = +1./std::sqrt(3.); return xi; } // ----------------------------------- integration point weights ----------------------------------- inline xt::xtensor w() { static const size_t nip = 8; xt::xtensor w = xt::empty({nip}); w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; w(4) = 1.; w(5) = 1.; w(6) = 1.; w(7) = 1.; return w; } // ------------------------------------------------------------------------------------------------- } // ================================ GooseFEM::Element::Hex8::Nodal ================================ namespace Nodal { // --------------------------------- number of integration points ---------------------------------- inline size_t nip() { return 8; } // ----------------------- integration point coordinates (local coordinates) ----------------------- inline xt::xtensor xi() { static const size_t nip = 8; static const size_t ndim = 3; xt::xtensor xi = xt::empty({nip,ndim}); xi(0,0) = -1.; xi(0,1) = -1.; xi(0,2) = -1.; xi(1,0) = +1.; xi(1,1) = -1.; xi(1,2) = -1.; xi(2,0) = +1.; xi(2,1) = +1.; xi(2,2) = -1.; xi(3,0) = -1.; xi(3,1) = +1.; xi(3,2) = -1.; xi(4,0) = -1.; xi(4,1) = -1.; xi(4,2) = +1.; xi(5,0) = +1.; xi(5,1) = -1.; xi(5,2) = +1.; xi(6,0) = +1.; xi(6,1) = +1.; xi(6,2) = +1.; xi(7,0) = -1.; xi(7,1) = +1.; xi(7,2) = +1.; return xi; } // ----------------------------------- integration point weights ----------------------------------- inline xt::xtensor w() { static const size_t nip = 8; xt::xtensor w = xt::empty({nip}); w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; w(4) = 1.; w(5) = 1.; w(6) = 1.; w(7) = 1.; return w; } // ------------------------------------------------------------------------------------------------- } // ================================================================================================= // ------------------------------------------ constructor ------------------------------------------ inline Quadrature::Quadrature(const xt::xtensor &x) : m_x(x) { assert( m_x.shape()[1] == m_nne ); assert( m_x.shape()[2] == m_ndim ); // set integration scheme m_xi = Gauss::xi(); m_w = Gauss::w (); // extract number of elements m_nelem = m_x.shape()[0]; m_nip = m_w.size(); // allocate arrays // - shape functions m_N = xt::empty({m_nip,m_nne}); // - shape function gradients in local coordinates m_dNxi = xt::empty({m_nip,m_nne,m_ndim}); // - shape function gradients in global coordinates m_dNx = xt::empty({m_nelem,m_nip,m_nne,m_ndim}); // - integration point volume m_vol = xt::empty({m_nelem,m_nip}); // shape functions for ( size_t k = 0 ; k < m_nip ; ++k ) { m_N(k,0) = .125 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)) * (1.-m_xi(k,2)); m_N(k,1) = .125 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)) * (1.-m_xi(k,2)); m_N(k,2) = .125 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)) * (1.-m_xi(k,2)); m_N(k,3) = .125 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)) * (1.-m_xi(k,2)); m_N(k,4) = .125 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)) * (1.+m_xi(k,2)); m_N(k,5) = .125 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)) * (1.+m_xi(k,2)); m_N(k,6) = .125 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)) * (1.+m_xi(k,2)); m_N(k,7) = .125 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)) * (1.+m_xi(k,2)); } // shape function gradients in local coordinates for ( size_t k = 0 ; k < m_nip ; ++k ) { // - dN / dxi_0 m_dNxi(k,0,0) = -.125*(1.-m_xi(k,1))*(1.-m_xi(k,2)); m_dNxi(k,1,0) = +.125*(1.-m_xi(k,1))*(1.-m_xi(k,2)); m_dNxi(k,2,0) = +.125*(1.+m_xi(k,1))*(1.-m_xi(k,2)); m_dNxi(k,3,0) = -.125*(1.+m_xi(k,1))*(1.-m_xi(k,2)); m_dNxi(k,4,0) = -.125*(1.-m_xi(k,1))*(1.+m_xi(k,2)); m_dNxi(k,5,0) = +.125*(1.-m_xi(k,1))*(1.+m_xi(k,2)); m_dNxi(k,6,0) = +.125*(1.+m_xi(k,1))*(1.+m_xi(k,2)); m_dNxi(k,7,0) = -.125*(1.+m_xi(k,1))*(1.+m_xi(k,2)); // - dN / dxi_1 m_dNxi(k,0,1) = -.125*(1.-m_xi(k,0))*(1.-m_xi(k,2)); m_dNxi(k,1,1) = -.125*(1.+m_xi(k,0))*(1.-m_xi(k,2)); m_dNxi(k,2,1) = +.125*(1.+m_xi(k,0))*(1.-m_xi(k,2)); m_dNxi(k,3,1) = +.125*(1.-m_xi(k,0))*(1.-m_xi(k,2)); m_dNxi(k,4,1) = -.125*(1.-m_xi(k,0))*(1.+m_xi(k,2)); m_dNxi(k,5,1) = -.125*(1.+m_xi(k,0))*(1.+m_xi(k,2)); m_dNxi(k,6,1) = +.125*(1.+m_xi(k,0))*(1.+m_xi(k,2)); m_dNxi(k,7,1) = +.125*(1.-m_xi(k,0))*(1.+m_xi(k,2)); // - dN / dxi_2 m_dNxi(k,0,2) = -.125*(1.-m_xi(k,0))*(1.-m_xi(k,1)); m_dNxi(k,1,2) = -.125*(1.+m_xi(k,0))*(1.-m_xi(k,1)); m_dNxi(k,2,2) = -.125*(1.+m_xi(k,0))*(1.+m_xi(k,1)); m_dNxi(k,3,2) = -.125*(1.-m_xi(k,0))*(1.+m_xi(k,1)); m_dNxi(k,4,2) = +.125*(1.-m_xi(k,0))*(1.-m_xi(k,1)); m_dNxi(k,5,2) = +.125*(1.+m_xi(k,0))*(1.-m_xi(k,1)); m_dNxi(k,6,2) = +.125*(1.+m_xi(k,0))*(1.+m_xi(k,1)); m_dNxi(k,7,2) = +.125*(1.-m_xi(k,0))*(1.+m_xi(k,1)); } // compute the shape function gradients, based on "x" compute_dN(); } // ------------------------------------------ constructor ------------------------------------------ inline Quadrature::Quadrature(const xt::xtensor &x, const xt::xtensor &xi, const xt::xtensor &w) : m_x(x), m_w(w), m_xi(xi) { assert( m_x.shape()[1] == m_nne ); assert( m_x.shape()[2] == m_ndim ); // extract shape m_nelem = m_x.shape()[0]; m_nip = m_w.size(); assert( m_xi.shape()[0] == m_nip ); assert( m_xi.shape()[1] == m_ndim ); assert( m_w .size() == m_nip ); // allocate arrays // - shape functions m_N = xt::empty({m_nip,m_nne}); // - shape function gradients in local coordinates m_dNxi = xt::empty({m_nip,m_nne,m_ndim}); // - shape function gradients in global coordinates m_dNx = xt::empty({m_nelem,m_nip,m_nne,m_ndim}); // - integration point volume m_vol = xt::empty({m_nelem,m_nip}); // shape functions for ( size_t k = 0 ; k < m_nip ; ++k ) { m_N(k,0) = .125 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)) * (1.-m_xi(k,2)); m_N(k,1) = .125 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)) * (1.-m_xi(k,2)); m_N(k,2) = .125 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)) * (1.-m_xi(k,2)); m_N(k,3) = .125 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)) * (1.-m_xi(k,2)); m_N(k,4) = .125 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)) * (1.+m_xi(k,2)); m_N(k,5) = .125 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)) * (1.+m_xi(k,2)); m_N(k,6) = .125 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)) * (1.+m_xi(k,2)); m_N(k,7) = .125 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)) * (1.+m_xi(k,2)); } // shape function gradients in local coordinates for ( size_t k = 0 ; k < m_nip ; ++k ) { // - dN / dxi_0 m_dNxi(k,0,0) = -.125*(1.-m_xi(k,1))*(1.-m_xi(k,2)); m_dNxi(k,1,0) = +.125*(1.-m_xi(k,1))*(1.-m_xi(k,2)); m_dNxi(k,2,0) = +.125*(1.+m_xi(k,1))*(1.-m_xi(k,2)); m_dNxi(k,3,0) = -.125*(1.+m_xi(k,1))*(1.-m_xi(k,2)); m_dNxi(k,4,0) = -.125*(1.-m_xi(k,1))*(1.+m_xi(k,2)); m_dNxi(k,5,0) = +.125*(1.-m_xi(k,1))*(1.+m_xi(k,2)); m_dNxi(k,6,0) = +.125*(1.+m_xi(k,1))*(1.+m_xi(k,2)); m_dNxi(k,7,0) = -.125*(1.+m_xi(k,1))*(1.+m_xi(k,2)); // - dN / dxi_1 m_dNxi(k,0,1) = -.125*(1.-m_xi(k,0))*(1.-m_xi(k,2)); m_dNxi(k,1,1) = -.125*(1.+m_xi(k,0))*(1.-m_xi(k,2)); m_dNxi(k,2,1) = +.125*(1.+m_xi(k,0))*(1.-m_xi(k,2)); m_dNxi(k,3,1) = +.125*(1.-m_xi(k,0))*(1.-m_xi(k,2)); m_dNxi(k,4,1) = -.125*(1.-m_xi(k,0))*(1.+m_xi(k,2)); m_dNxi(k,5,1) = -.125*(1.+m_xi(k,0))*(1.+m_xi(k,2)); m_dNxi(k,6,1) = +.125*(1.+m_xi(k,0))*(1.+m_xi(k,2)); m_dNxi(k,7,1) = +.125*(1.-m_xi(k,0))*(1.+m_xi(k,2)); // - dN / dxi_2 m_dNxi(k,0,2) = -.125*(1.-m_xi(k,0))*(1.-m_xi(k,1)); m_dNxi(k,1,2) = -.125*(1.+m_xi(k,0))*(1.-m_xi(k,1)); m_dNxi(k,2,2) = -.125*(1.+m_xi(k,0))*(1.+m_xi(k,1)); m_dNxi(k,3,2) = -.125*(1.-m_xi(k,0))*(1.+m_xi(k,1)); m_dNxi(k,4,2) = +.125*(1.-m_xi(k,0))*(1.-m_xi(k,1)); m_dNxi(k,5,2) = +.125*(1.+m_xi(k,0))*(1.-m_xi(k,1)); m_dNxi(k,6,2) = +.125*(1.+m_xi(k,0))*(1.+m_xi(k,1)); m_dNxi(k,7,2) = +.125*(1.-m_xi(k,0))*(1.+m_xi(k,1)); } // compute the shape function gradients, based on "x" compute_dN(); } // --------------------------- integration volume (per tensor-component) --------------------------- inline void Quadrature::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 k = 0 ; k < m_nip ; ++k ) qscalar(e,k) = m_vol(e,k); } // ------------------------------------------------------------------------------------------------- inline void Quadrature::dV(xt::xtensor &qtensor) const { assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nne ); assert( qtensor.shape()[2] >= m_ndim ); assert( qtensor.shape()[3] >= m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t k = 0 ; k < m_nip ; ++k ) for ( size_t i = 0 ; i < qtensor.shape()[2] ; ++i ) for ( size_t j = 0 ; j < qtensor.shape()[3] ; ++j ) qtensor(e,k,i,j) = m_vol(e,k); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::dV() const { xt::xtensor out = xt::empty({m_nelem, m_nip}); this->dV(out); return out; } // -------------------------------------- number of elements --------------------------------------- inline size_t Quadrature::nelem() const { return m_nelem; } // ---------------------------------- number of nodes per element ---------------------------------- inline size_t Quadrature::nne() const { return m_nne; } // ------------------------------------- number of dimensions -------------------------------------- inline size_t Quadrature::ndim() const { return m_ndim; } // --------------------------------- number of integration points ---------------------------------- inline size_t Quadrature::nip() const { return m_nip; } // --------------------------------------- update positions ---------------------------------------- inline void Quadrature::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 m_x = x; // update the shape function gradients for the new "x" compute_dN(); } // ------------------------ shape function gradients in global coordinates ------------------------- inline void Quadrature::compute_dN() { // loop over all elements (in parallel) #pragma omp parallel { // - allocate T2 J; T2 Jinv; #pragma omp for for ( size_t e = 0 ; e < m_nelem ; ++e ) { // alias nodal positions auto x = xt::view(m_x, e, xt::all(), xt::all()); // loop over integration points for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias auto dNxi = xt::view(m_dNxi, k, xt::all(), xt::all()); auto dNx = xt::view(m_dNx , e, k, xt::all(), xt::all()); // - zero-initialize J *= 0.0; // - Jacobian for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) J(i,j) += dNxi(m,i) * x(m,j); // - determinant and inverse of the Jacobian double Jdet = inv(J, Jinv); - // - shape function gradients wrt global coordinates - for ( size_t m = 0 ; m < m_nne ; ++m ) { + // - 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) + Jinv(0,2) * dNxi(m,2); dNx(m,1) = Jinv(1,0) * dNxi(m,0) + Jinv(1,1) * dNxi(m,1) + Jinv(1,2) * dNxi(m,2); dNx(m,2) = Jinv(2,0) * dNxi(m,0) + Jinv(2,1) * dNxi(m,1) + Jinv(2,2) * dNxi(m,2); } // - copy to matrix: integration point volume m_vol(e,k) = m_w(k) * Jdet; } } } } // ------------------- dyadic product "qtensor(i,j) = dNdx(m,i) * elemvec(m,j)" -------------------- inline void Quadrature::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_ndim ); assert( qtensor.shape()[3] >= m_ndim ); // zero-initialize output: matrix of tensors qtensor *= 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::view(elemvec, e, xt::all(), xt::all()); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias auto dNx = xt::view(m_dNx , e, k, xt::all() , xt::all() ); auto gradu = xt::view(qtensor, e, k, xt::range(0,m_ndim), xt::range(0,m_ndim)); // - evaluate dyadic product for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) gradu(i,j) += dNx(m,i) * u(m,j); } } } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::gradN_vector(const xt::xtensor &elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); this->gradN_vector(elemvec, qtensor); return qtensor; } // ---------------------------------- transpose of "gradN_vector" ---------------------------------- inline void Quadrature::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_ndim ); assert( qtensor.shape()[3] >= m_ndim ); // zero-initialize output: matrix of tensors qtensor *= 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::view(elemvec, e, xt::all(), xt::all()); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias auto dNx = xt::view(m_dNx , e, k, xt::all() , xt::all() ); auto gradu = xt::view(qtensor, e, k, xt::range(0,m_ndim), xt::range(0,m_ndim)); - // - evaluate dyadic product + // - evaluate transpose of dyadic product for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) gradu(j,i) += dNx(m,i) * u(m,j); } } } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::gradN_vector_T(const xt::xtensor &elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); this->gradN_vector_T(elemvec, qtensor); return qtensor; } // ------------------------------- symmetric part of "gradN_vector" -------------------------------- inline void Quadrature::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_ndim ); assert( qtensor.shape()[3] >= m_ndim ); // zero-initialize output: matrix of tensors qtensor *= 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::view(elemvec, e, xt::all(), xt::all()); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias shape function gradients (local coordinates) auto dNx = xt::view(m_dNx , e, k, xt::all() , xt::all() ); auto eps = xt::view(qtensor, e, k, xt::range(0,m_ndim), xt::range(0,m_ndim)); - // - evaluate dyadic product + // - evaluate symmetrized dyadic product for ( size_t m = 0 ; m < m_nne ; ++m ) { for ( size_t i = 0 ; i < m_ndim ; ++i ) { for ( size_t j = 0 ; j < m_ndim ; ++j ) { eps(i,j) += dNx(m,i) * u(m,j) / 2.; eps(j,i) += dNx(m,i) * u(m,j) / 2.; } } } } } } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::symGradN_vector(const xt::xtensor &elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); this->symGradN_vector(elemvec, qtensor); return qtensor; } // ------- scalar product "elemmat(m*ndim+i,n*ndim+i) = N(m) * qscalar * N(n)"; for all "i" -------- inline void Quadrature::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 *= 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::view(elemmat, e, xt::all(), xt::all()); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias shape functions auto N = xt::view(m_N, k, xt::all()); // - alias double vol = m_vol (e,k); // integration point volume double rho = qscalar(e,k); // integration point scalar (e.g. density) // - 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 ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - M(m*m_ndim+i, n*m_ndim+i) += N(m) * rho * N(n) * vol; + 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; + M(m*m_ndim+2, n*m_ndim+2) += N(m) * rho * N(n) * vol; + } + } } } } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::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; } // ------------ integral of dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" ------------ inline void Quadrature::int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor, xt::xtensor &elemvec) const { assert( qtensor.shape()[0] == m_nelem ); // number of elements assert( qtensor.shape()[1] == m_nip ); // number of integration points assert( qtensor.shape()[2] >= m_ndim ); // number of dimensions assert( qtensor.shape()[3] >= m_ndim ); // number of dimensions assert( elemvec.shape()[0] == m_nelem ); // number of elements assert( elemvec.shape()[1] == m_nne ); // number of nodes per element assert( elemvec.shape()[2] == m_ndim ); // number of dimensions // zero-initialize output: matrix of vectors elemvec *= 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::view(elemvec, e, xt::all(), xt::all()); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias auto dNx = xt::view(m_dNx , e, k, xt::all(), xt::all()); auto sig = xt::view(qtensor, e, k, xt::range(0,m_ndim), xt::range(0,m_ndim)); double vol = m_vol(e,k); // - evaluate dot product, and assemble for ( size_t m = 0 ; m < m_nne ; ++m ) - for ( size_t i = 0 ; i < m_ndim ; ++i ) - for ( size_t j = 0 ; j < m_ndim ; ++j ) - f(m,j) += dNx(m,i) * sig(i,j) * vol; + { + f(m,0) += ( dNx(m,0) * sig(0,0) + dNx(m,1) * sig(1,0) + dNx(m,2) * sig(2,0) ) * vol; + f(m,1) += ( dNx(m,0) * sig(0,1) + dNx(m,1) * sig(1,1) + dNx(m,2) * sig(2,1) ) * vol; + f(m,2) += ( dNx(m,0) * sig(0,2) + dNx(m,1) * sig(1,2) + dNx(m,2) * sig(2,2) ) * vol; + } } } } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::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; } // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/ElementQuad4.h b/src/GooseFEM/ElementQuad4.h index 3f096a6..03b61f9 100644 --- a/src/GooseFEM/ElementQuad4.h +++ b/src/GooseFEM/ElementQuad4.h @@ -1,134 +1,133 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_ELEMENTQUAD4_H #define GOOSEFEM_ELEMENTQUAD4_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // =================================== GooseFEM::Element::Quad4 ==================================== namespace GooseFEM { namespace Element { namespace Quad4 { // ======================================== tensor algebra ========================================= static const size_t ndim = 2; using T2 = xt::xtensor_fixed>; -inline double det(const T2 &A); -inline T2 inv(const T2 &A); +inline double inv(const T2 &A, T2 &Ainv); // ================================ GooseFEM::Element::Quad4::Gauss ================================ 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 } // ================================ GooseFEM::Element::Quad4::Nodal ================================ 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 { 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] private: // compute "vol" and "dNdx" based on current "x" void compute_dN(); public: // convention: // "elemmat" - matrices stored per element - ArrD - [nelem, nne*ndim, nne*ndim] // "elemvec" - nodal vectors stored per element - ArrD - [nelem, nne, ndim] // "qtensor" - integration point tensor - ArrD - [nelem, nip, ndim, ndim] // "qscalar" - integration point scalar - ArrD - [nelem, nip] // constructor: integration point coordinates and weights are optional (default: Gauss) Quadrature() = default; Quadrature(const xt::xtensor &x); 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; // number of elements size_t nne() const; // number of nodes per element size_t ndim() const; // number of dimension size_t nip() const; // number of integration points // return integration volume // - in-place void dV(xt::xtensor &qscalar) const; void dV(xt::xtensor &qtensor) const; // same volume for all tensor components // - return qscalar xt::xtensor dV() const; // dyadic product "qtensor(i,j) += dNdx(m,i) * elemvec(m,j)", its transpose and its symmetric part // - in-place 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; // - return qtensor 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; // integral of the scalar product "elemmat(m*ndim+i,n*ndim+i) += N(m) * qscalar * N(n) * dV" // - in-place void int_N_scalar_NT_dV(const xt::xtensor &qscalar, xt::xtensor &elemmat) const; // - return elemmat xt::xtensor int_N_scalar_NT_dV(const xt::xtensor &qscalar) const; // integral of the dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" // - in-place void int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor, xt::xtensor &elemvec) const; // - return elemvec xt::xtensor int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor) const; }; // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/ElementQuad4.hpp b/src/GooseFEM/ElementQuad4.hpp index e203d72..b04354f 100644 --- a/src/GooseFEM/ElementQuad4.hpp +++ b/src/GooseFEM/ElementQuad4.hpp @@ -1,647 +1,642 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_ELEMENTQUAD4_CPP #define GOOSEFEM_ELEMENTQUAD4_CPP // ------------------------------------------------------------------------------------------------- #include "ElementQuad4.h" // =================================== GooseFEM::Element::Quad4 ==================================== namespace GooseFEM { namespace Element { namespace Quad4 { // ======================================== tensor algebra ========================================= -inline double det(const T2 &A) -{ - return A[0] * A[3] - A[1] * A[2]; -} - -// ------------------------------------------------------------------------------------------------- - -inline T2 inv(const T2 &A) +inline double inv(const T2 &A, T2 &Ainv) { // compute determinant - double D = det(A); - - // allocate result - T2 C = xt::empty({ndim,ndim}); + double det = A[0] * A[3] - A[1] * A[2]; // compute inverse - C[0] = A[3] / D; - C[1] = -1. * A[1] / D; - C[2] = -1. * A[2] / D; - C[3] = A[0] / D; + Ainv[0] = A[3] / det; + Ainv[1] = -1. * A[1] / det; + Ainv[2] = -1. * A[2] / det; + Ainv[3] = A[0] / det; - return C; + return det; } // ================================ GooseFEM::Element::Quad4::Gauss ================================ namespace Gauss { // --------------------------------- number of integration points ---------------------------------- inline size_t nip() { return 4; } // ----------------------- integration point coordinates (local coordinates) ----------------------- inline xt::xtensor xi() { static const size_t nip = 4; static const size_t ndim = 2; xt::xtensor xi = xt::empty({nip,ndim}); xi(0,0) = -1./std::sqrt(3.); xi(0,1) = -1./std::sqrt(3.); xi(1,0) = +1./std::sqrt(3.); xi(1,1) = -1./std::sqrt(3.); xi(2,0) = +1./std::sqrt(3.); xi(2,1) = +1./std::sqrt(3.); xi(3,0) = -1./std::sqrt(3.); xi(3,1) = +1./std::sqrt(3.); return xi; } // ----------------------------------- integration point weights ----------------------------------- inline xt::xtensor w() { static const size_t nip = 4; xt::xtensor w = xt::empty({nip}); w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; return w; } // ------------------------------------------------------------------------------------------------- } // ================================ GooseFEM::Element::Quad4::Nodal ================================ namespace Nodal { // --------------------------------- number of integration points ---------------------------------- inline size_t nip() { return 4; } // ----------------------- integration point coordinates (local coordinates) ----------------------- inline xt::xtensor xi() { static const size_t nip = 4; static const size_t ndim = 2; xt::xtensor xi = xt::empty({nip,ndim}); xi(0,0) = -1.; xi(0,1) = -1.; xi(1,0) = +1.; xi(1,1) = -1.; xi(2,0) = +1.; xi(2,1) = +1.; xi(3,0) = -1.; xi(3,1) = +1.; return xi; } // ----------------------------------- integration point weights ----------------------------------- inline xt::xtensor w() { static const size_t nip = 4; xt::xtensor w = xt::empty({nip}); w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; return w; } // ------------------------------------------------------------------------------------------------- } // ================================================================================================= // ------------------------------------------ constructor ------------------------------------------ inline Quadrature::Quadrature(const xt::xtensor &x) : m_x(x) { assert( m_x.shape()[1] == m_nne ); assert( m_x.shape()[2] == m_ndim ); // set integration scheme m_xi = Gauss::xi(); m_w = Gauss::w (); // extract number of elements m_nelem = m_x.shape()[0]; m_nip = m_w.size(); // allocate arrays // - shape functions m_N = xt::empty({m_nip,m_nne}); // - shape function gradients in local coordinates m_dNxi = xt::empty({m_nip,m_nne,m_ndim}); // - shape function gradients in global coordinates m_dNx = xt::empty({m_nelem,m_nip,m_nne,m_ndim}); // - integration point volume m_vol = xt::empty({m_nelem,m_nip}); // shape functions for ( size_t k = 0 ; k < m_nip ; ++k ) { m_N(k,0) = .25 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)); m_N(k,1) = .25 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)); m_N(k,2) = .25 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)); m_N(k,3) = .25 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)); } // shape function gradients in local coordinates for ( size_t k = 0 ; k < m_nip ; ++k ) { // - dN / dxi_0 m_dNxi(k,0,0) = -.25*(1.-m_xi(k,1)); m_dNxi(k,1,0) = +.25*(1.-m_xi(k,1)); m_dNxi(k,2,0) = +.25*(1.+m_xi(k,1)); m_dNxi(k,3,0) = -.25*(1.+m_xi(k,1)); // - dN / dxi_1 m_dNxi(k,0,1) = -.25*(1.-m_xi(k,0)); m_dNxi(k,1,1) = -.25*(1.+m_xi(k,0)); m_dNxi(k,2,1) = +.25*(1.+m_xi(k,0)); m_dNxi(k,3,1) = +.25*(1.-m_xi(k,0)); } // compute the shape function gradients, based on "x" compute_dN(); } // ------------------------------------------ constructor ------------------------------------------ inline Quadrature::Quadrature(const xt::xtensor &x, const xt::xtensor &xi, const xt::xtensor &w) : m_x(x), m_w(w), m_xi(xi) { assert( m_x.shape()[1] == m_nne ); assert( m_x.shape()[2] == m_ndim ); // extract shape m_nelem = m_x.shape()[0]; m_nip = m_w.size(); assert( m_xi.shape()[0] == m_nip ); assert( m_xi.shape()[1] == m_ndim ); assert( m_w .size() == m_nip ); // allocate arrays // - shape functions m_N = xt::empty({m_nip,m_nne}); // - shape function gradients in local coordinates m_dNxi = xt::empty({m_nip,m_nne,m_ndim}); // - shape function gradients in global coordinates m_dNx = xt::empty({m_nelem,m_nip,m_nne,m_ndim}); // - integration point volume m_vol = xt::empty({m_nelem,m_nip}); // shape functions for ( size_t k = 0 ; k < m_nip ; ++k ) { m_N(k,0) = .25 * (1.-m_xi(k,0)) * (1.-m_xi(k,1)); m_N(k,1) = .25 * (1.+m_xi(k,0)) * (1.-m_xi(k,1)); m_N(k,2) = .25 * (1.+m_xi(k,0)) * (1.+m_xi(k,1)); m_N(k,3) = .25 * (1.-m_xi(k,0)) * (1.+m_xi(k,1)); } // shape function gradients in local coordinates for ( size_t k = 0 ; k < m_nip ; ++k ) { // - dN / dxi_0 m_dNxi(k,0,0) = -.25*(1.-m_xi(k,1)); m_dNxi(k,1,0) = +.25*(1.-m_xi(k,1)); m_dNxi(k,2,0) = +.25*(1.+m_xi(k,1)); m_dNxi(k,3,0) = -.25*(1.+m_xi(k,1)); // - dN / dxi_1 m_dNxi(k,0,1) = -.25*(1.-m_xi(k,0)); m_dNxi(k,1,1) = -.25*(1.+m_xi(k,0)); m_dNxi(k,2,1) = +.25*(1.+m_xi(k,0)); m_dNxi(k,3,1) = +.25*(1.-m_xi(k,0)); } // compute the shape function gradients, based on "x" compute_dN(); } // --------------------------- integration volume (per tensor-component) --------------------------- inline void Quadrature::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 k = 0 ; k < m_nip ; ++k ) qscalar(e,k) = m_vol(e,k); } // ------------------------------------------------------------------------------------------------- inline void Quadrature::dV(xt::xtensor &qtensor) const { assert( qtensor.shape()[0] == m_nelem ); assert( qtensor.shape()[1] == m_nne ); assert( qtensor.shape()[2] >= m_ndim ); assert( qtensor.shape()[3] >= m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t k = 0 ; k < m_nip ; ++k ) for ( size_t i = 0 ; i < qtensor.shape()[2] ; ++i ) for ( size_t j = 0 ; j < qtensor.shape()[3] ; ++j ) qtensor(e,k,i,j) = m_vol(e,k); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::dV() const { xt::xtensor out = xt::empty({m_nelem, m_nip}); this->dV(out); return out; } // -------------------------------------- number of elements --------------------------------------- inline size_t Quadrature::nelem() const { return m_nelem; } // ---------------------------------- number of nodes per element ---------------------------------- inline size_t Quadrature::nne() const { return m_nne; } // ------------------------------------- number of dimensions -------------------------------------- inline size_t Quadrature::ndim() const { return m_ndim; } // --------------------------------- number of integration points ---------------------------------- inline size_t Quadrature::nip() const { return m_nip; } // --------------------------------------- update positions ---------------------------------------- inline void Quadrature::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 m_x = x; // update the shape function gradients for the new "x" compute_dN(); } // ------------------------ shape function gradients in global coordinates ------------------------- inline void Quadrature::compute_dN() { // loop over all elements (in parallel) - #pragma omp parallel for - for ( size_t e = 0 ; e < m_nelem ; ++e ) + #pragma omp parallel { - // alias nodal positions - auto x = xt::view(m_x, e, xt::all(), xt::all()); + // - allocate + T2 J; + T2 Jinv; - // loop over integration points - for ( size_t k = 0 ; k < m_nip ; ++k ) + #pragma omp for + for ( size_t e = 0 ; e < m_nelem ; ++e ) { - // - alias - auto dNxi = xt::view(m_dNxi, k, xt::all(), xt::all()); - auto dNx = xt::view(m_dNx , e, k, xt::all(), xt::all()); - - // - Jacobian (loops unrolled for efficiency) - // J(i,j) += dNxi(m,i) * x(m,j); - T2 J = xt::empty({m_ndim,m_ndim}); - 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 = det(J); - T2 Jinv = inv(J); - - // - 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 ) + // alias nodal positions + auto x = xt::view(m_x, e, xt::all(), xt::all()); + + // loop over integration points + for ( size_t k = 0 ; k < m_nip ; ++k ) { - 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); - } + // - alias + auto dNxi = xt::view(m_dNxi, k, xt::all(), xt::all()); + auto dNx = xt::view(m_dNx , e, k, xt::all(), xt::all()); + + // - 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); + } - // - copy to matrix: integration point volume - m_vol(e,k) = m_w(k) * Jdet; + // - copy to matrix: integration point volume + m_vol(e,k) = m_w(k) * Jdet; + } } } } // ------------------- dyadic product "qtensor(i,j) = dNdx(m,i) * elemvec(m,j)" -------------------- inline void Quadrature::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_ndim ); assert( qtensor.shape()[3] >= m_ndim ); // zero-initialize output: matrix of tensors qtensor *= 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::view(elemvec, e, xt::all(), xt::all()); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias auto dNx = xt::view(m_dNx , e, k, xt::all() , xt::all() ); auto gradu = xt::view(qtensor, e, k, xt::range(0,m_ndim), xt::range(0,m_ndim)); // - 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 xt::xtensor Quadrature::gradN_vector(const xt::xtensor &elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); this->gradN_vector(elemvec, qtensor); return qtensor; } // ---------------------------------- transpose of "gradN_vector" ---------------------------------- inline void Quadrature::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_ndim ); assert( qtensor.shape()[3] >= m_ndim ); // zero-initialize output: matrix of tensors qtensor *= 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::view(elemvec, e, xt::all(), xt::all()); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias auto dNx = xt::view(m_dNx , e, k, xt::all() , xt::all() ); auto gradu = xt::view(qtensor, e, k, xt::range(0,m_ndim), xt::range(0,m_ndim)); - // - evaluate dyadic product (loops unrolled for efficiency) + // - 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 xt::xtensor Quadrature::gradN_vector_T(const xt::xtensor &elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); this->gradN_vector_T(elemvec, qtensor); return qtensor; } // ------------------------------- symmetric part of "gradN_vector" -------------------------------- inline void Quadrature::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_ndim ); assert( qtensor.shape()[3] >= m_ndim ); // zero-initialize output: matrix of tensors qtensor *= 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::view(elemvec, e, xt::all(), xt::all()); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias shape function gradients (local coordinates) auto dNx = xt::view(m_dNx , e, k, xt::all() , xt::all() ); auto eps = xt::view(qtensor, e, k, xt::range(0,m_ndim), xt::range(0,m_ndim)); - // - evaluate dyadic product (loops unrolled for efficiency) + // - 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,1)*u(0,0) + dNx(1,1)*u(1,0) + dNx(2,1)*u(2,0) + dNx(3,1)*u(3,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) ) / 2.; eps(1,0) = eps(0,1); } } } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Quadrature::symGradN_vector(const xt::xtensor &elemvec) const { xt::xtensor qtensor = xt::empty({m_nelem, m_nip, m_ndim, m_ndim}); this->symGradN_vector(elemvec, qtensor); return qtensor; } // ------- scalar product "elemmat(m*ndim+i,n*ndim+i) = N(m) * qscalar * N(n)"; for all "i" -------- inline void Quadrature::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 *= 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::view(elemmat, e, xt::all(), xt::all()); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias shape functions auto N = xt::view(m_N, k, xt::all()); // - alias double vol = m_vol (e,k); // integration point volume double rho = qscalar(e,k); // integration point scalar (e.g. density) // - 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 xt::xtensor Quadrature::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; } // ------------ integral of dot product "elemvec(m,j) += dNdx(m,i) * qtensor(i,j) * dV" ------------ inline void Quadrature::int_gradN_dot_tensor2_dV(const xt::xtensor &qtensor, xt::xtensor &elemvec) const { assert( qtensor.shape()[0] == m_nelem ); // number of elements assert( qtensor.shape()[1] == m_nip ); // number of integration points assert( qtensor.shape()[2] >= m_ndim ); // number of dimensions assert( qtensor.shape()[3] >= m_ndim ); // number of dimensions assert( elemvec.shape()[0] == m_nelem ); // number of elements assert( elemvec.shape()[1] == m_nne ); // number of nodes per element assert( elemvec.shape()[2] == m_ndim ); // number of dimensions // zero-initialize output: matrix of vectors elemvec *= 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::view(elemvec, e, xt::all(), xt::all()); // loop over all integration points in element "e" for ( size_t k = 0 ; k < m_nip ; ++k ) { // - alias auto dNx = xt::view(m_dNx , e, k, xt::all(), xt::all()); auto sig = xt::view(qtensor, e, k, xt::range(0,m_ndim), xt::range(0,m_ndim)); double vol = m_vol(e,k); // - 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 xt::xtensor Quadrature::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; } // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/Vector.h b/src/GooseFEM/Vector.h index ee59e46..9f4fb20 100644 --- a/src/GooseFEM/Vector.h +++ b/src/GooseFEM/Vector.h @@ -1,100 +1,124 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_VECTOR_H #define GOOSEFEM_VECTOR_H // ------------------------------------------------------------------------------------------------- #include "GooseFEM.h" // =========================================== GooseFEM ============================================ namespace GooseFEM { // ------------------------------------------------------------------------------------------------- class Vector { private: // information xt::xtensor m_conn; // connectivity [nelem, nne ] xt::xtensor m_dofs; // DOF-numbers per node [nnode, ndim] xt::xtensor m_part; // DOF-numbers per node, after partitioning [nnode, ndim] xt::xtensor m_iiu; // DOF-numbers that are unknown [nnu] xt::xtensor m_iip; // DOF-numbers that are prescribed [nnp] // dimensions size_t m_nelem; // number of elements size_t m_nne; // number of nodes per element size_t m_nnode; // number of nodes size_t m_ndim; // number of dimensions size_t m_ndof; // number of DOFs size_t m_nnu; // number of unknown DOFs size_t m_nnp; // number of prescribed DOFs public: // notation: // "nodevec" - nodal vectors - [nnode, ndim] // "elemvec" - nodal vectors stored per element - [nelem, nne, ndim] // "dofval" - DOF values - [ndof] // "dofval_u" - DOF values (Unknown) - [nnu] // "dofval_p" - DOF values (Prescribed) - [nnp] // constructor Vector(){}; Vector(const xt::xtensor &conn, const xt::xtensor &dofs); Vector(const xt::xtensor &conn, const xt::xtensor &dofs, const xt::xtensor &iip); // dimensions size_t nelem() const; // number of elements size_t nne() const; // number of nodes per element size_t nnode() const; // number of nodes size_t ndim() const; // number of dimensions size_t ndof() const; // number of DOFs size_t nnu() const; // number of unknown DOFs size_t nnp() const; // number of prescribed DOFs // DOF lists xt::xtensor dofs() const; // DOFs xt::xtensor iiu() const; // unknown DOFs xt::xtensor iip() const; // prescribed DOFs + // convert vectors (overwrite entries that occur more that once) -- no allocation + void asDofs (const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &dofval ) const; + void asDofs (const xt::xtensor &nodevec , xt::xtensor &dofval ) const; + void asDofs (const xt::xtensor &elemvec , xt::xtensor &dofval ) const; + void asDofs_u (const xt::xtensor &nodevec , xt::xtensor &dofval ) const; + void asDofs_u (const xt::xtensor &elemvec , xt::xtensor &dofval ) const; + void asDofs_p (const xt::xtensor &nodevec , xt::xtensor &dofval ) const; + void asDofs_p (const xt::xtensor &elemvec , xt::xtensor &dofval ) const; + void asNode (const xt::xtensor &dofval , xt::xtensor &nodevec) const; + void asNode (const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &nodevec) const; + void asNode (const xt::xtensor &elemvec , xt::xtensor &nodevec) const; + void asElement(const xt::xtensor &dofval , xt::xtensor &elemvec) const; + void asElement(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &elemvec) const; + void asElement(const xt::xtensor &nodevec , xt::xtensor &elemvec) const; + + // assemble vectors (adds entries that occur more that once) -- no allocation + void assembleDofs (const xt::xtensor &nodevec, xt::xtensor &dofval ) const; + void assembleDofs (const xt::xtensor &elemvec, xt::xtensor &dofval ) const; + void assembleDofs_u(const xt::xtensor &nodevec, xt::xtensor &dofval ) const; + void assembleDofs_u(const xt::xtensor &elemvec, xt::xtensor &dofval ) const; + void assembleDofs_p(const xt::xtensor &nodevec, xt::xtensor &dofval ) const; + void assembleDofs_p(const xt::xtensor &elemvec, xt::xtensor &dofval ) const; + void assembleNode (const xt::xtensor &elemvec, xt::xtensor &nodevec) const; + // convert vectors (overwrite entries that occur more that once) xt::xtensor asDofs (const xt::xtensor &dofval_u, const xt::xtensor &dofval_p) const; xt::xtensor asDofs (const xt::xtensor &nodevec ) const; xt::xtensor asDofs (const xt::xtensor &elemvec ) const; xt::xtensor asDofs_u (const xt::xtensor &nodevec ) const; xt::xtensor asDofs_u (const xt::xtensor &elemvec ) const; xt::xtensor asDofs_p (const xt::xtensor &nodevec ) const; xt::xtensor asDofs_p (const xt::xtensor &elemvec ) const; xt::xtensor asNode (const xt::xtensor &dofval ) const; xt::xtensor asNode (const xt::xtensor &dofval_u, const xt::xtensor &dofval_p) const; xt::xtensor asNode (const xt::xtensor &elemvec ) const; xt::xtensor asElement(const xt::xtensor &dofval ) const; xt::xtensor asElement(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p) const; xt::xtensor asElement(const xt::xtensor &nodevec ) const; // assemble vectors (adds entries that occur more that once) xt::xtensor assembleDofs (const xt::xtensor &nodevec) const; xt::xtensor assembleDofs (const xt::xtensor &elemvec) const; xt::xtensor assembleDofs_u(const xt::xtensor &nodevec) const; xt::xtensor assembleDofs_u(const xt::xtensor &elemvec) const; xt::xtensor assembleDofs_p(const xt::xtensor &nodevec) const; xt::xtensor assembleDofs_p(const xt::xtensor &elemvec) const; xt::xtensor assembleNode (const xt::xtensor &elemvec) const; }; // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/Vector.hpp b/src/GooseFEM/Vector.hpp index a9ec46d..db9ec8c 100644 --- a/src/GooseFEM/Vector.hpp +++ b/src/GooseFEM/Vector.hpp @@ -1,726 +1,719 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_VECTOR_CPP #define GOOSEFEM_VECTOR_CPP // ------------------------------------------------------------------------------------------------- #include "Vector.h" // =========================================== GooseFEM ============================================ namespace GooseFEM { // ------------------------------------------ constructor ------------------------------------------ inline Vector::Vector(const xt::xtensor &conn, const xt::xtensor &dofs) : m_conn(conn), m_dofs(dofs) { m_nelem = m_conn.shape()[0]; m_nne = m_conn.shape()[1]; m_nnode = m_dofs.shape()[0]; m_ndim = m_dofs.shape()[1]; m_ndof = xt::amax(m_dofs)[0] + 1; m_nnp = 0; m_nnu = m_ndof; m_iiu = xt::arange(m_ndof); m_iip = xt::empty({0}); // check consistency assert( xt::amax(m_conn)[0] + 1 == m_nnode ); assert( m_ndof <= m_nnode * m_ndim ); } // ------------------------------------------ constructor ------------------------------------------ inline Vector::Vector(const xt::xtensor &conn, const xt::xtensor &dofs, const xt::xtensor &iip) : m_conn(conn), m_dofs(dofs), m_iip(iip) { m_nelem = m_conn.shape()[0]; m_nne = m_conn.shape()[1]; m_nnode = m_dofs.shape()[0]; m_ndim = m_dofs.shape()[1]; m_ndof = xt::amax(m_dofs)[0] + 1; m_nnp = m_iip.size(); m_nnu = m_ndof - m_nnp; m_iiu = xt::empty({m_nnu}); // check consistency assert( xt::amax(m_conn)[0] + 1 == m_nnode ); assert( m_ndof <= m_nnode * m_ndim ); // reorder DOFs such that they can be used for partitioning; renumber such that // "iiu" -> beginning // "iip" -> end // (otherwise preserving the order) // this array can be used to assemble to/from partitioned arrays m_part = Mesh::reorder(m_dofs, m_iip, "end"); // extract unknown DOFs #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(n,i) < m_nnu ) m_iiu(m_part(n,i)) = m_dofs(n,i); } // -------------------------------------- number of elements --------------------------------------- inline size_t Vector::nelem() const { return m_nelem; } // ---------------------------------- number of nodes per element ---------------------------------- inline size_t Vector::nne() const { return m_nne; } // ---------------------------------------- number of nodes ---------------------------------------- inline size_t Vector::nnode() const { return m_nnode; } // ------------------------------------- number of dimensions -------------------------------------- inline size_t Vector::ndim() const { return m_ndim; } // ---------------------------------------- number of DOFs ----------------------------------------- inline size_t Vector::ndof() const { return m_ndof; } // ------------------------------------ number of unknown DOFs ------------------------------------- inline size_t Vector::nnu() const { return m_nnu; } // ----------------------------------- number of prescribed DOFs ----------------------------------- inline size_t Vector::nnp() const { return m_nnp; } // ------------------------------------------ return DOFs ------------------------------------------ inline xt::xtensor Vector::dofs() const { return m_dofs; } // -------------------------------------- return unknown DOFs -------------------------------------- inline xt::xtensor Vector::iiu() const { return m_iiu; } // ------------------------------------ return prescribed DOFs ------------------------------------- inline xt::xtensor Vector::iip() const { return m_iip; } // --------------------------------------- dofval -> dofval ---------------------------------------- inline void Vector::asDofs(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &dofval) const { assert( dofval_u.size() == m_nnu ); assert( dofval_p.size() == m_nnp ); assert( dofval.size() == m_ndof ); #pragma omp parallel for for ( size_t i = 0 ; i < m_nnu ; ++i ) dofval(m_iiu(i)) = dofval_u(i); #pragma omp parallel for for ( size_t i = 0 ; i < m_nnp ; ++i ) dofval(m_iip(i)) = dofval_p(i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asDofs(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p) const { xt::xtensor dofval = xt::empty({m_ndof}); this->asDofs(dofval_u, dofval_p, dofval); return dofval; } // --------------------------------------- nodevec -> dofval --------------------------------------- inline void Vector::asDofs(const xt::xtensor &nodevec, xt::xtensor &dofval) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval.size() == m_ndof ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) dofval(m_dofs(n,i)) = nodevec(n,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asDofs(const xt::xtensor &nodevec) const { xt::xtensor dofval = xt::empty({m_ndof}); this->asDofs(nodevec, dofval); return dofval; } // --------------------------------------- nodevec -> dofval --------------------------------------- inline void Vector::asDofs_u(const xt::xtensor &nodevec, xt::xtensor &dofval) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval.size() == m_nnu ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(n,i) < m_nnu ) dofval(m_part(n,i)) = nodevec(n,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asDofs_u(const xt::xtensor &nodevec) const { xt::xtensor dofval = xt::empty({m_nnu}); this->asDofs_u(nodevec, dofval); return dofval; } // --------------------------------------- nodevec -> dofval --------------------------------------- inline void Vector::asDofs_p(const xt::xtensor &nodevec, xt::xtensor &dofval) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval.size() == m_nnp ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(n,i) >= m_nnu ) dofval(m_part(n,i)-m_nnu) = nodevec(n,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asDofs_p(const xt::xtensor &nodevec) const { xt::xtensor dofval = xt::empty({m_nnp}); this->asDofs_p(nodevec, dofval); return dofval; } // --------------------------------------- elemvec -> dofval --------------------------------------- inline void Vector::asDofs(const xt::xtensor &elemvec, xt::xtensor &dofval) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval.size() == m_ndof ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) dofval(m_dofs(m_conn(e,m),i)) = elemvec(e,m,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asDofs(const xt::xtensor &elemvec) const { xt::xtensor dofval = xt::empty({m_ndof}); this->asDofs(elemvec, dofval); return dofval; } // --------------------------------------- elemvec -> dofval --------------------------------------- inline void Vector::asDofs_u(const xt::xtensor &elemvec, xt::xtensor &dofval) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval.size() == m_nnu ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(m_conn(e,m),i) < m_nnu ) dofval(m_part(m_conn(e,m),i)) = elemvec(e,m,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asDofs_u(const xt::xtensor &elemvec) const { xt::xtensor dofval = xt::empty({m_nnu}); this->asDofs_u(elemvec, dofval); return dofval; } // --------------------------------------- elemvec -> dofval --------------------------------------- inline void Vector::asDofs_p(const xt::xtensor &elemvec, xt::xtensor &dofval) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval.size() == m_nnp ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(m_conn(e,m),i) >= m_nnu ) dofval(m_part(m_conn(e,m),i)-m_nnu) = elemvec(e,m,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asDofs_p(const xt::xtensor &elemvec) const { xt::xtensor dofval = xt::empty({m_nnp}); this->asDofs_p(elemvec, dofval); return dofval; } // --------------------------------------- dofval -> nodevec --------------------------------------- inline void Vector::asNode(const xt::xtensor &dofval, xt::xtensor &nodevec) const { assert( dofval.size() == m_ndof ); assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) nodevec(n,i) = dofval(m_dofs(n,i)); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asNode(const xt::xtensor &dofval) const { xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); this->asNode(dofval, nodevec); return nodevec; } // --------------------------------------- dofval -> nodevec --------------------------------------- inline void Vector::asNode(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &nodevec) const { assert( dofval_u.size() == m_nnu ); assert( dofval_p.size() == m_nnp ); assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) { for ( size_t i = 0 ; i < m_ndim ; ++i ) { if ( m_part(n,i) < m_nnu ) nodevec(n,i) = dofval_u(m_part(n,i) ); else nodevec(n,i) = dofval_p(m_part(n,i)-m_nnu); } } } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asNode(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p) const { xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); this->asNode(dofval_u, dofval_p, nodevec); return nodevec; } // --------------------------------------- elemvec -> nodevec --------------------------------------- inline void Vector::asNode(const xt::xtensor &elemvec, xt::xtensor &nodevec) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) nodevec(m_conn(e,m),i) = elemvec(e,m,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asNode(const xt::xtensor &elemvec) const { xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); this->asNode(elemvec, nodevec); return nodevec; } // --------------------------------------- dofval -> elemvec --------------------------------------- inline void Vector::asElement(const xt::xtensor &dofval, xt::xtensor &elemvec) const { assert( dofval.size() == m_ndof ); assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) elemvec(e,m,i) = dofval(m_dofs(m_conn(e,m),i)); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asElement(const xt::xtensor &dofval) const { xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); this->asElement(dofval, elemvec); return elemvec; } // --------------------------------------- dofval -> elemvec --------------------------------------- inline void Vector::asElement(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p, xt::xtensor &elemvec) const { assert( dofval_u.size() == m_nnu ); assert( dofval_p.size() == m_nnp ); assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) { for ( size_t m = 0 ; m < m_nne ; ++m ) { for ( size_t i = 0 ; i < m_ndim ; ++i ) { if ( m_part(m_conn(e,m),i) Vector::asElement(const xt::xtensor &dofval_u, const xt::xtensor &dofval_p) const { xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); this->asElement(dofval_u, dofval_p, elemvec); return elemvec; } // -------------------------------------- nodevec -> elemvec --------------------------------------- inline void Vector::asElement(const xt::xtensor &nodevec, xt::xtensor &elemvec) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) elemvec(e,m,i) = nodevec(m_conn(e,m),i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::asElement(const xt::xtensor &nodevec) const { xt::xtensor elemvec = xt::empty({m_nelem, m_nne, m_ndim}); this->asElement(nodevec, elemvec); return elemvec; } // --------------------------------------- nodevec -> dofval --------------------------------------- inline void Vector::assembleDofs(const xt::xtensor &nodevec, xt::xtensor &dofval) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval.size() == m_ndof ); dofval *= 0.0; - #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) dofval(m_dofs(n,i)) += nodevec(n,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::assembleDofs(const xt::xtensor &nodevec) const { xt::xtensor dofval = xt::empty({m_ndof}); this->assembleDofs(nodevec, dofval); return dofval; } // --------------------------------------- nodevec -> dofval --------------------------------------- inline void Vector::assembleDofs_u(const xt::xtensor &nodevec, xt::xtensor &dofval) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval.size() == m_nnu ); dofval *= 0.0; - #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(n,i) < m_nnu ) dofval(m_part(n,i)) += nodevec(n,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::assembleDofs_u(const xt::xtensor &nodevec) const { xt::xtensor dofval = xt::empty({m_nnu}); this->assembleDofs_u(nodevec, dofval); return dofval; } // --------------------------------------- nodevec -> dofval --------------------------------------- inline void Vector::assembleDofs_p(const xt::xtensor &nodevec, xt::xtensor &dofval) const { assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); assert( dofval.size() == m_nnp ); dofval *= 0.0; - #pragma omp parallel for for ( size_t n = 0 ; n < m_nnode ; ++n ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(n,i) >= m_nnu ) dofval(m_part(n,i)-m_nnu) += nodevec(n,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::assembleDofs_p(const xt::xtensor &nodevec) const { xt::xtensor dofval = xt::empty({m_nnp}); this->assembleDofs_p(nodevec, dofval); return dofval; } // --------------------------------------- elemvec -> dofval --------------------------------------- inline void Vector::assembleDofs(const xt::xtensor &elemvec, xt::xtensor &dofval) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval.size() == m_ndof ); dofval *= 0.0; - #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) dofval(m_dofs(m_conn(e,m),i)) += elemvec(e,m,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::assembleDofs(const xt::xtensor &elemvec) const { xt::xtensor dofval = xt::empty({m_ndof}); this->assembleDofs(elemvec, dofval); return dofval; } // --------------------------------------- elemvec -> dofval --------------------------------------- inline void Vector::assembleDofs_u(const xt::xtensor &elemvec, xt::xtensor &dofval) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval.size() == m_nnu ); dofval *= 0.0; - #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(m_conn(e,m),i) < m_nnu ) dofval(m_dofs(m_conn(e,m),i)) += elemvec(e,m,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::assembleDofs_u(const xt::xtensor &elemvec) const { xt::xtensor dofval = xt::empty({m_nnu}); this->assembleDofs_u(elemvec, dofval); return dofval; } // --------------------------------------- elemvec -> dofval --------------------------------------- inline void Vector::assembleDofs_p(const xt::xtensor &elemvec, xt::xtensor &dofval) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( dofval.size() == m_nnp ); dofval *= 0.0; - #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) if ( m_part(m_conn(e,m),i) >= m_nnu ) dofval(m_dofs(m_conn(e,m),i)-m_nnu) += elemvec(e,m,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::assembleDofs_p(const xt::xtensor &elemvec) const { xt::xtensor dofval = xt::empty({m_nnp}); this->assembleDofs_p(elemvec, dofval); return dofval; } // --------------------------------------- elemvec -> nodevec --------------------------------------- inline void Vector::assembleNode(const xt::xtensor &elemvec, xt::xtensor &nodevec) const { assert( elemvec.shape()[0] == m_nelem ); assert( elemvec.shape()[1] == m_nne ); assert( elemvec.shape()[2] == m_ndim ); assert( nodevec.shape()[0] == m_nnode ); assert( nodevec.shape()[1] == m_ndim ); nodevec *= 0.0; - #pragma omp parallel for for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t m = 0 ; m < m_nne ; ++m ) for ( size_t i = 0 ; i < m_ndim ; ++i ) nodevec(m_conn(e,m),i) += elemvec(e,m,i); } // ------------------------------------------------------------------------------------------------- inline xt::xtensor Vector::assembleNode(const xt::xtensor &elemvec) const { xt::xtensor nodevec = xt::empty({m_nnode, m_ndim}); this->assembleNode(elemvec, nodevec); return nodevec; } // ------------------------------------------------------------------------------------------------- } // namespace ... // ================================================================================================= #endif diff --git a/src/GooseFEM/python.cpp b/src/GooseFEM/python.cpp index f5cd57d..d919410 100644 --- a/src/GooseFEM/python.cpp +++ b/src/GooseFEM/python.cpp @@ -1,926 +1,937 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #include #include #include #include #include +#include + #include "GooseFEM.h" // ================================================================================================= // abbreviate name-space namespace py = pybind11; namespace M = GooseFEM; -// abbreviate types -typedef GooseFEM::ColD ColD; -typedef GooseFEM::ColS ColS; -typedef GooseFEM::MatD MatD; -typedef GooseFEM::MatS MatS; -typedef GooseFEM::ArrD ArrD; - -// abbreviate const types -typedef const GooseFEM::ColD cColD; -typedef const GooseFEM::ColS cColS; -typedef const GooseFEM::MatD cMatD; -typedef const GooseFEM::MatS cMatS; -typedef const GooseFEM::ArrD cArrD; - // ================================================================================================= class PyGeometry : public GooseFEM::Dynamics::Geometry { public: // inherit the constructors using GooseFEM::Dynamics::Geometry::Geometry; + using Arr1 = xt::xtensor; + using Arr2 = xt::xtensor; + using Arr3 = xt::xtensor; // trampoline - ColD solve_A() override { PYBIND11_OVERLOAD_PURE( ColD, GooseFEM::Dynamics::Geometry, solve_A ); } - ColD solve_V() override { PYBIND11_OVERLOAD_PURE( ColD, GooseFEM::Dynamics::Geometry, solve_V ); } - MatD u() const override { PYBIND11_OVERLOAD_PURE( MatD, GooseFEM::Dynamics::Geometry, u ); } - MatD v() const override { PYBIND11_OVERLOAD_PURE( MatD, GooseFEM::Dynamics::Geometry, v ); } - MatD a() const override { PYBIND11_OVERLOAD_PURE( MatD, GooseFEM::Dynamics::Geometry, a ); } - ColD dofs_u() const override { PYBIND11_OVERLOAD_PURE( ColD, GooseFEM::Dynamics::Geometry, dofs_u ); } - ColD dofs_v() const override { PYBIND11_OVERLOAD_PURE( ColD, GooseFEM::Dynamics::Geometry, dofs_v ); } - ColD dofs_a() const override { PYBIND11_OVERLOAD_PURE( ColD, GooseFEM::Dynamics::Geometry, dofs_a ); } - void set_u(const MatD &nodevec) override { PYBIND11_OVERLOAD_PURE( void, GooseFEM::Dynamics::Geometry, set_u , nodevec); } - void set_u(const ColD &dofval ) override { PYBIND11_OVERLOAD_PURE( void, GooseFEM::Dynamics::Geometry, set_u , dofval ); } - void set_v(const ColD &dofval ) override { PYBIND11_OVERLOAD_PURE( void, GooseFEM::Dynamics::Geometry, set_v , dofval ); } - void set_a(const ColD &dofval ) override { PYBIND11_OVERLOAD_PURE( void, GooseFEM::Dynamics::Geometry, set_a , dofval ); } + xt::xtensor solve_A() override { PYBIND11_OVERLOAD_PURE( Arr1, GooseFEM::Dynamics::Geometry, solve_A ); } + xt::xtensor solve_V() override { PYBIND11_OVERLOAD_PURE( Arr1, GooseFEM::Dynamics::Geometry, solve_V ); } + xt::xtensor u() const override { PYBIND11_OVERLOAD_PURE( Arr2, GooseFEM::Dynamics::Geometry, u ); } + xt::xtensor v() const override { PYBIND11_OVERLOAD_PURE( Arr2, GooseFEM::Dynamics::Geometry, v ); } + xt::xtensor a() const override { PYBIND11_OVERLOAD_PURE( Arr2, GooseFEM::Dynamics::Geometry, a ); } + xt::xtensor dofs_u() const override { PYBIND11_OVERLOAD_PURE( Arr1, GooseFEM::Dynamics::Geometry, dofs_u ); } + xt::xtensor dofs_v() const override { PYBIND11_OVERLOAD_PURE( Arr1, GooseFEM::Dynamics::Geometry, dofs_v ); } + xt::xtensor dofs_a() const override { PYBIND11_OVERLOAD_PURE( Arr1, GooseFEM::Dynamics::Geometry, dofs_a ); } + void set_u(const xt::xtensor &nodevec) override { PYBIND11_OVERLOAD_PURE( void, GooseFEM::Dynamics::Geometry, set_u , nodevec); } + void set_u(const xt::xtensor &dofval ) override { PYBIND11_OVERLOAD_PURE( void, GooseFEM::Dynamics::Geometry, set_u , dofval ); } + void set_v(const xt::xtensor &dofval ) override { PYBIND11_OVERLOAD_PURE( void, GooseFEM::Dynamics::Geometry, set_v , dofval ); } + void set_a(const xt::xtensor &dofval ) override { PYBIND11_OVERLOAD_PURE( void, GooseFEM::Dynamics::Geometry, set_a , dofval ); } }; // =========================================== GooseFEM ============================================ PYBIND11_MODULE(GooseFEM, m) { m.doc() = "Some simple finite element meshes and operations"; // ================================= GooseFEM - GooseFEM/Vector.h ================================== py::class_(m, "Vector") // constructor .def( - py::init(), + py::init &, const xt::xtensor &>(), + "Class to switch between DOF/nodal/element views of vectors", + py::arg("conn"), + py::arg("dofs") + ) + // constructor + .def( + py::init &, const xt::xtensor &, const xt::xtensor &>(), "Class to switch between DOF/nodal/element views of vectors", py::arg("conn"), py::arg("dofs"), - py::arg("iip")=GooseFEM::ColS() + py::arg("iip") ) // methods .def("nelem", &M::Vector::nelem) .def("nne" , &M::Vector::nne ) .def("nnode", &M::Vector::nnode) .def("ndim" , &M::Vector::ndim ) .def("ndof" , &M::Vector::ndof ) .def("nnu" , &M::Vector::nnu ) .def("nnp" , &M::Vector::nnp ) // - .def("iiu" , &M::Vector::iiu ) .def("iip" , &M::Vector::iip ) // - - .def("asDofs" , py::overload_cast(&M::Vector::asDofs , py::const_)) - .def("asDofs" , py::overload_cast(&M::Vector::asDofs , py::const_)) - .def("asDofs" , py::overload_cast(&M::Vector::asDofs , py::const_)) - .def("asDofs_u" , py::overload_cast(&M::Vector::asDofs_u , py::const_)) - .def("asDofs_u" , py::overload_cast(&M::Vector::asDofs_u , py::const_)) - .def("asDofs_p" , py::overload_cast(&M::Vector::asDofs_p , py::const_)) - .def("asDofs_p" , py::overload_cast(&M::Vector::asDofs_p , py::const_)) - .def("asNode" , py::overload_cast(&M::Vector::asNode , py::const_)) - .def("asNode" , py::overload_cast(&M::Vector::asNode , py::const_)) - .def("asNode" , py::overload_cast(&M::Vector::asNode , py::const_)) - .def("asElement" , py::overload_cast(&M::Vector::asElement , py::const_)) - .def("asElement" , py::overload_cast(&M::Vector::asElement , py::const_)) - .def("asElement" , py::overload_cast(&M::Vector::asElement , py::const_)) + .def("asDofs" , py::overload_cast&,const xt::xtensor&>(&M::Vector::asDofs , py::const_)) + .def("asDofs" , py::overload_cast& >(&M::Vector::asDofs , py::const_)) + .def("asDofs" , py::overload_cast& >(&M::Vector::asDofs , py::const_)) + .def("asDofs_u" , py::overload_cast& >(&M::Vector::asDofs_u , py::const_)) + .def("asDofs_u" , py::overload_cast& >(&M::Vector::asDofs_u , py::const_)) + .def("asDofs_p" , py::overload_cast& >(&M::Vector::asDofs_p , py::const_)) + .def("asDofs_p" , py::overload_cast& >(&M::Vector::asDofs_p , py::const_)) + .def("asNode" , py::overload_cast& >(&M::Vector::asNode , py::const_)) + .def("asNode" , py::overload_cast&,const xt::xtensor&>(&M::Vector::asNode , py::const_)) + .def("asNode" , py::overload_cast& >(&M::Vector::asNode , py::const_)) + .def("asElement" , py::overload_cast& >(&M::Vector::asElement , py::const_)) + .def("asElement" , py::overload_cast&,const xt::xtensor&>(&M::Vector::asElement , py::const_)) + .def("asElement" , py::overload_cast& >(&M::Vector::asElement , py::const_)) // - - .def("assembleDofs" , py::overload_cast(&M::Vector::assembleDofs , py::const_)) - .def("assembleDofs" , py::overload_cast(&M::Vector::assembleDofs , py::const_)) - .def("assembleDofs_u", py::overload_cast(&M::Vector::assembleDofs_u, py::const_)) - .def("assembleDofs_u", py::overload_cast(&M::Vector::assembleDofs_u, py::const_)) - .def("assembleDofs_p", py::overload_cast(&M::Vector::assembleDofs_p, py::const_)) - .def("assembleDofs_p", py::overload_cast(&M::Vector::assembleDofs_p, py::const_)) - .def("assembleNode" , py::overload_cast(&M::Vector::assembleNode , py::const_)) + .def("assembleDofs" , py::overload_cast&>(&M::Vector::assembleDofs , py::const_)) + .def("assembleDofs" , py::overload_cast&>(&M::Vector::assembleDofs , py::const_)) + .def("assembleDofs_u", py::overload_cast&>(&M::Vector::assembleDofs_u, py::const_)) + .def("assembleDofs_u", py::overload_cast&>(&M::Vector::assembleDofs_u, py::const_)) + .def("assembleDofs_p", py::overload_cast&>(&M::Vector::assembleDofs_p, py::const_)) + .def("assembleDofs_p", py::overload_cast&>(&M::Vector::assembleDofs_p, py::const_)) + .def("assembleNode" , py::overload_cast&>(&M::Vector::assembleNode , py::const_)) // print to screen .def("__repr__", - [](const GooseFEM::Vector &a){ return ""; } + [](const GooseFEM::Vector &){ return ""; } ); // ============================= GooseFEM - GooseFEM/MatrixDiagonal.h ============================== py::class_(m, "MatrixDiagonal") +// constructor + .def( + py::init &, const xt::xtensor &>(), + "Class to switch between DOF/nodal/element views of vectors", + py::arg("conn"), + py::arg("dofs") + ) // constructor .def( - py::init(), + py::init &, const xt::xtensor &, const xt::xtensor &>(), "Class to switch between DOF/nodal/element views of vectors", py::arg("conn"), py::arg("dofs"), - py::arg("iip")=GooseFEM::ColS() + py::arg("iip") ) // methods .def("nelem", &M::MatrixDiagonal::nelem) .def("nne" , &M::MatrixDiagonal::nne ) .def("nnode", &M::MatrixDiagonal::nnode) .def("ndim" , &M::MatrixDiagonal::ndim ) .def("ndof" , &M::MatrixDiagonal::ndof ) .def("nnu" , &M::MatrixDiagonal::nnu ) .def("nnp" , &M::MatrixDiagonal::nnp ) // - .def("iiu" , &M::MatrixDiagonal::iiu ) .def("iip" , &M::MatrixDiagonal::iip ) // - - .def("dot" , &M::MatrixDiagonal::dot ) - .def("dot_u", py::overload_cast(&M::MatrixDiagonal::dot_u, py::const_)) - .def("dot_u", py::overload_cast(&M::MatrixDiagonal::dot_u, py::const_)) - .def("dot_p", py::overload_cast(&M::MatrixDiagonal::dot_p, py::const_)) - .def("dot_p", py::overload_cast(&M::MatrixDiagonal::dot_p, py::const_)) + .def("dot" , &M::MatrixDiagonal::dot ) + .def("dot_u", py::overload_cast& >(&M::MatrixDiagonal::dot_u, py::const_)) + .def("dot_u", py::overload_cast&,const xt::xtensor&>(&M::MatrixDiagonal::dot_u, py::const_)) + .def("dot_p", py::overload_cast& >(&M::MatrixDiagonal::dot_p, py::const_)) + .def("dot_p", py::overload_cast&,const xt::xtensor&>(&M::MatrixDiagonal::dot_p, py::const_)) // - .def("check_diagonal", &M::MatrixDiagonal::check_diagonal) .def("assemble" , &M::MatrixDiagonal::assemble ) // .def("set" , &M::MatrixDiagonal::set ) // .def("set_uu" , &M::MatrixDiagonal::set_uu ) // .def("set_pp" , &M::MatrixDiagonal::set_pp ) - .def("solve" , &M::MatrixDiagonal::solve, "Solve", py::arg("rhs"), py::arg("u_p")=ColD()) - .def("solve_u" , &M::MatrixDiagonal::solve_u ) + .def("solve" , py::overload_cast& >(&M::MatrixDiagonal::solve ), "Solve", py::arg("rhs" ) ) + .def("solve" , py::overload_cast&, const xt::xtensor&>(&M::MatrixDiagonal::solve ), "Solve", py::arg("rhs" ), py::arg("u_p")) + .def("solve_u" , py::overload_cast& >(&M::MatrixDiagonal::solve_u), "Solve", py::arg("rhs_u") ) + .def("solve_u" , py::overload_cast&, const xt::xtensor&>(&M::MatrixDiagonal::solve_u), "Solve", py::arg("rhs_u"), py::arg("u_p")) // .def("rhs_p" , &M::MatrixDiagonal::rhs_p ) .def("asDiagonal" , &M::MatrixDiagonal::asDiagonal ) // .def("asDiagonal_uu" , &M::MatrixDiagonal::asDiagonal_uu ) // .def("asDiagonal_pp" , &M::MatrixDiagonal::asDiagonal_pp ) // .def("asSparse" , &M::MatrixDiagonal::asSparse ) // .def("asSparse_uu" , &M::MatrixDiagonal::asSparse_uu ) // .def("asSparse_up" , &M::MatrixDiagonal::asSparse_up ) // .def("asSparse_pu" , &M::MatrixDiagonal::asSparse_pu ) // .def("asSparse_pp" , &M::MatrixDiagonal::asSparse_pp ) // .def("asDense" , &M::MatrixDiagonal::asDense ) // .def("asDense_uu" , &M::MatrixDiagonal::asDense_uu ) // .def("asDense_up" , &M::MatrixDiagonal::asDense_up ) // .def("asDense_pu" , &M::MatrixDiagonal::asDense_pu ) // .def("asDense_pp" , &M::MatrixDiagonal::asDense_pp ) // print to screen .def("__repr__", - [](const GooseFEM::MatrixDiagonal &a){ return ""; } + [](const GooseFEM::MatrixDiagonal &){ return ""; } ); // =========================== GooseFEM::Dynamics - GooseFEM/Dynamics.h ============================ py::module mDynamics = m.def_submodule("Dynamics", "Solve routines for dynamic FEM"); // ------------------------------------------------------------------------------------------------- mDynamics.def("Verlet", &GooseFEM::Dynamics::Verlet, "Verlet time integration", py::arg("geometry"), py::arg("dt"), py::arg("nstep")=1 ); // ------------------------------------------------------------------------------------------------- mDynamics.def("velocityVerlet", &GooseFEM::Dynamics::velocityVerlet, "Velocity-Verlet time integration", py::arg("geometry"), py::arg("dt"), py::arg("nstep")=1 ); -// ------------------------------------------------------------------------------------------------- +// // ------------------------------------------------------------------------------------------------- py::class_(mDynamics, "Geometry") // constructor .def(py::init<>()) // methods .def("solve_A" , &GooseFEM::Dynamics::Geometry::solve_A) .def("solve_V" , &GooseFEM::Dynamics::Geometry::solve_V) .def("u" , &GooseFEM::Dynamics::Geometry::u ) .def("v" , &GooseFEM::Dynamics::Geometry::v ) .def("a" , &GooseFEM::Dynamics::Geometry::a ) .def("dofs_u" , &GooseFEM::Dynamics::Geometry::dofs_u ) .def("dofs_v" , &GooseFEM::Dynamics::Geometry::dofs_v ) .def("dofs_a" , &GooseFEM::Dynamics::Geometry::dofs_a ) .def("set_v" , &GooseFEM::Dynamics::Geometry::set_v ) .def("set_a" , &GooseFEM::Dynamics::Geometry::set_a ) - .def("set_u" , py::overload_cast(&GooseFEM::Dynamics::Geometry::set_u)) - .def("set_u" , py::overload_cast(&GooseFEM::Dynamics::Geometry::set_u)) + .def("set_u" , py::overload_cast &>(&GooseFEM::Dynamics::Geometry::set_u)) + .def("set_u" , py::overload_cast &>(&GooseFEM::Dynamics::Geometry::set_u)) // print to screen .def("__repr__", - [](const GooseFEM::Dynamics::Geometry &a){ return ""; } + [](const GooseFEM::Dynamics::Geometry &){ return ""; } ); // ============================ GooseFEM::Element - GooseFEM/Element.h ============================= py::module mElement = m.def_submodule("Element", "Generic element routines"); // ------------------------------------------------------------------------------------------------- mElement.def("asElementVector", &GooseFEM::Element::asElementVector, "convert nodal vector [nnode, ndim] to nodal vector stored per element [nelem, nne, ndim]", py::arg("conn"), py::arg("nodevec") ); // ------------------------------------------------------------------------------------------------- mElement.def("assembleElementVector", &GooseFEM::Element::assembleNodeVector, "assemble nodal vector stored per element [nelem, nne, ndim] to nodal vector [nnode, ndim]", py::arg("conn"), py::arg("elemvec") ); // ====================== GooseFEM::Element::Quad4 - GooseFEM/ElementQuad4.h ======================= { // create sub-module py::module sm = mElement.def_submodule("Quad4", "Linear quadrilateral elements (2D)"); // abbreviate name-space namespace SM = GooseFEM::Element::Quad4; -using T2 = cppmat::tiny::cartesian::tensor2 ; -using T2s = cppmat::tiny::cartesian::tensor2s; - // ------------------------------------------------------------------------------------------------- py::class_(sm, "Quadrature") // constructor .def( - py::init(), + py::init &>(), + "Quadrature", + py::arg("x") + ) + // constructor + .def( + py::init &, const xt::xtensor &, const xt::xtensor &>(), "Quadrature", py::arg("x"), - py::arg("xi")=ArrD::Zero({0}), - py::arg("w")=ArrD::Zero({0}) + py::arg("xi"), + py::arg("w") ) // sizes .def("nelem" , &SM::Quadrature::nelem) .def("nne" , &SM::Quadrature::nne) .def("ndim" , &SM::Quadrature::ndim) .def("nip" , &SM::Quadrature::nip) - .def("dV" , &SM::Quadrature::dV, py::arg("ncomp")=0) - .def("gradN_vector" , &SM::Quadrature::gradN_vector) - .def("gradN_vector_T" , &SM::Quadrature::gradN_vector_T) - .def("symGradN_vector" , &SM::Quadrature::symGradN_vector) - .def("int_N_scalar_NT_dV" , &SM::Quadrature::int_N_scalar_NT_dV) - .def("int_gradN_dot_tensor2_dV" , &SM::Quadrature::int_gradN_dot_tensor2_dV) - .def("int_gradN_dot_tensor2s_dV", &SM::Quadrature::int_gradN_dot_tensor2_dV) + .def("dV" , py::overload_cast<>(&SM::Quadrature::dV, py::const_)) + .def("gradN_vector" , py::overload_cast &>(&SM::Quadrature::gradN_vector, py::const_)) + .def("gradN_vector_T" , py::overload_cast &>(&SM::Quadrature::gradN_vector_T, py::const_)) + .def("symGradN_vector" , py::overload_cast &>(&SM::Quadrature::symGradN_vector, py::const_)) + .def("int_N_scalar_NT_dV" , py::overload_cast &>(&SM::Quadrature::int_N_scalar_NT_dV, py::const_)) + .def("int_gradN_dot_tensor2_dV" , py::overload_cast &>(&SM::Quadrature::int_gradN_dot_tensor2_dV, py::const_)) // print to screen .def("__repr__", - [](const SM::Quadrature &a){ return ""; } + [](const SM::Quadrature &){ return ""; } ); // ------------------------------------------------------------------------------------------------- { py::module ssm = sm.def_submodule("Gauss", "Gauss quadrature"); namespace SSM = GooseFEM::Element::Quad4::Gauss; ssm.def("nip", &SSM::nip); ssm.def("xi" , &SSM::xi); ssm.def("w" , &SSM::w); } // ------------------------------------------------------------------------------------------------- { py::module ssm = sm.def_submodule("Nodal", "Nodal quadrature"); namespace SSM = GooseFEM::Element::Quad4::Nodal; ssm.def("nip", &SSM::nip); ssm.def("xi" , &SSM::xi); ssm.def("w" , &SSM::w); } // ------------------------------------------------------------------------------------------------- } // ====================== GooseFEM::Element::Hex8 - GooseFEM/ElementHex8.h ======================= { // create sub-module py::module sm = mElement.def_submodule("Hex8", "Linear hexahedron (brick) elements (3D)"); // abbreviate name-space namespace SM = GooseFEM::Element::Hex8; -using T2 = cppmat::tiny::cartesian::tensor2 ; -using T2s = cppmat::tiny::cartesian::tensor2s; - // ------------------------------------------------------------------------------------------------- py::class_(sm, "Quadrature") // constructor .def( - py::init(), + py::init &>(), + "Quadrature", + py::arg("x") + ) + // constructor + .def( + py::init &, const xt::xtensor &, const xt::xtensor &>(), "Quadrature", py::arg("x"), - py::arg("xi")=ArrD::Zero({0}), - py::arg("w")=ArrD::Zero({0}) + py::arg("xi"), + py::arg("w") ) // sizes .def("nelem" , &SM::Quadrature::nelem) .def("nne" , &SM::Quadrature::nne) .def("ndim" , &SM::Quadrature::ndim) .def("nip" , &SM::Quadrature::nip) - .def("dV" , &SM::Quadrature::dV, py::arg("ncomp")=0) - .def("gradN_vector" , &SM::Quadrature::gradN_vector) - .def("gradN_vector_T" , &SM::Quadrature::gradN_vector_T) - .def("symGradN_vector" , &SM::Quadrature::symGradN_vector) - .def("int_N_scalar_NT_dV" , &SM::Quadrature::int_N_scalar_NT_dV) - .def("int_gradN_dot_tensor2_dV" , &SM::Quadrature::int_gradN_dot_tensor2_dV) - .def("int_gradN_dot_tensor2s_dV", &SM::Quadrature::int_gradN_dot_tensor2_dV) + .def("dV" , py::overload_cast<>(&SM::Quadrature::dV, py::const_)) + .def("gradN_vector" , py::overload_cast &>(&SM::Quadrature::gradN_vector, py::const_)) + .def("gradN_vector_T" , py::overload_cast &>(&SM::Quadrature::gradN_vector_T, py::const_)) + .def("symGradN_vector" , py::overload_cast &>(&SM::Quadrature::symGradN_vector, py::const_)) + .def("int_N_scalar_NT_dV" , py::overload_cast &>(&SM::Quadrature::int_N_scalar_NT_dV, py::const_)) + .def("int_gradN_dot_tensor2_dV" , py::overload_cast &>(&SM::Quadrature::int_gradN_dot_tensor2_dV, py::const_)) // print to screen .def("__repr__", - [](const SM::Quadrature &a){ return ""; } + [](const SM::Quadrature &){ return ""; } ); // ------------------------------------------------------------------------------------------------- { py::module ssm = sm.def_submodule("Gauss", "Gauss quadrature"); namespace SSM = GooseFEM::Element::Hex8::Gauss; ssm.def("nip", &SSM::nip); ssm.def("xi" , &SSM::xi); ssm.def("w" , &SSM::w); } // ------------------------------------------------------------------------------------------------- { py::module ssm = sm.def_submodule("Nodal", "Nodal quadrature"); namespace SSM = GooseFEM::Element::Hex8::Nodal; ssm.def("nip", &SSM::nip); ssm.def("xi" , &SSM::xi); ssm.def("w" , &SSM::w); } // ------------------------------------------------------------------------------------------------- } // =============================== GooseFEM::Mesh - GooseFEM/Mesh.h ================================ py::module mMesh = m.def_submodule("Mesh", "Generic mesh routines"); // ------------------------------------------------------------------------------------------------- mMesh.def("elem2node", &GooseFEM::Mesh::elem2node, "Elements connect to each node: [ number of elements , element numbers ]", py::arg("conn") ); // ------------------------------------------------------------------------------------------------- mMesh.def("coordination", &GooseFEM::Mesh::coordination, "Get the coordination number of each node: the number of elements connected to it", py::arg("conn") ); // ------------------------------------------------------------------------------------------------- mMesh.def("dofs", &GooseFEM::Mesh::dofs, "List with DOF-numbers (in sequential order)", py::arg("nnode"), py::arg("ndim") ); // ------------------------------------------------------------------------------------------------- -using renumber = GooseFEM::MatS(cMatS &); +using renumber = xt::xtensor(const xt::xtensor &); mMesh.def("renumber", - py::overload_cast((renumber*)&GooseFEM::Mesh::renumber), + py::overload_cast &>((renumber*)&GooseFEM::Mesh::renumber), "Renumber DOF-list to use the lowest possible index", py::arg("dofs") ); // ------------------------------------------------------------------------------------------------- -using reorder = GooseFEM::MatS(cMatS &, cColS&, std::string); +using reorder = xt::xtensor(const xt::xtensor &, const xt::xtensor&, std::string); mMesh.def("reorder", - py::overload_cast( + py::overload_cast&,const xt::xtensor&,std::string>( (reorder*)&GooseFEM::Mesh::reorder), "Renumber DOF-list to begin or end with 'idx'", py::arg("dofs"), py::arg("idx"), py::arg("location")="end" ); // ========================== GooseFEM::Mesh::Hex8 - GooseFEM/MeshHex8.h =========================== { // create sub-module py::module sm = mMesh.def_submodule("Hex8", "Linear hexahedron (brick) elements (3D)"); // abbreviate name-space namespace SM = GooseFEM::Mesh::Hex8; // ------------------------------------------------------------------------------------------------- py::class_(sm, "Regular") // constructor .def( py::init(), "mesh with nx*ny*nz 'pixels' and edge size h", py::arg("nx"), py::arg("ny"), py::arg("nz"), py::arg("h")=1. ) // sizes .def("nelem" , &SM::Regular::nelem ) .def("nnode" , &SM::Regular::nnode ) .def("nne" , &SM::Regular::nne ) .def("ndim" , &SM::Regular::ndim ) // mesh .def("coor" , &SM::Regular::coor ) .def("conn" , &SM::Regular::conn ) // boundary nodes: planes .def("nodesFront" , &SM::Regular::nodesFront ) .def("nodesBack" , &SM::Regular::nodesBack ) .def("nodesLeft" , &SM::Regular::nodesLeft ) .def("nodesRight" , &SM::Regular::nodesRight ) .def("nodesBottom" , &SM::Regular::nodesBottom ) .def("nodesTop" , &SM::Regular::nodesTop ) // boundary nodes: faces .def("nodesFrontFace" , &SM::Regular::nodesFrontFace ) .def("nodesBackFace" , &SM::Regular::nodesBackFace ) .def("nodesLeftFace" , &SM::Regular::nodesLeftFace ) .def("nodesRightFace" , &SM::Regular::nodesRightFace ) .def("nodesBottomFace" , &SM::Regular::nodesBottomFace ) .def("nodesTopFace" , &SM::Regular::nodesTopFace ) // boundary nodes: edges .def("nodesFrontBottomEdge" , &SM::Regular::nodesFrontBottomEdge ) .def("nodesFrontTopEdge" , &SM::Regular::nodesFrontTopEdge ) .def("nodesFrontLeftEdge" , &SM::Regular::nodesFrontLeftEdge ) .def("nodesFrontRightEdge" , &SM::Regular::nodesFrontRightEdge ) .def("nodesBackBottomEdge" , &SM::Regular::nodesBackBottomEdge ) .def("nodesBackTopEdge" , &SM::Regular::nodesBackTopEdge ) .def("nodesBackLeftEdge" , &SM::Regular::nodesBackLeftEdge ) .def("nodesBackRightEdge" , &SM::Regular::nodesBackRightEdge ) .def("nodesBottomLeftEdge" , &SM::Regular::nodesBottomLeftEdge ) .def("nodesBottomRightEdge" , &SM::Regular::nodesBottomRightEdge ) .def("nodesTopLeftEdge" , &SM::Regular::nodesTopLeftEdge ) .def("nodesTopRightEdge" , &SM::Regular::nodesTopRightEdge ) // boundary nodes: faces (aliases) .def("nodesBottomFrontEdge" , &SM::Regular::nodesBottomFrontEdge ) .def("nodesBottomBackEdge" , &SM::Regular::nodesBottomBackEdge ) .def("nodesTopFrontEdge" , &SM::Regular::nodesTopFrontEdge ) .def("nodesTopBackEdge" , &SM::Regular::nodesTopBackEdge ) .def("nodesLeftBottomEdge" , &SM::Regular::nodesLeftBottomEdge ) .def("nodesLeftFrontEdge" , &SM::Regular::nodesLeftFrontEdge ) .def("nodesLeftBackEdge" , &SM::Regular::nodesLeftBackEdge ) .def("nodesLeftTopEdge" , &SM::Regular::nodesLeftTopEdge ) .def("nodesRightBottomEdge" , &SM::Regular::nodesRightBottomEdge ) .def("nodesRightTopEdge" , &SM::Regular::nodesRightTopEdge ) .def("nodesRightFrontEdge" , &SM::Regular::nodesRightFrontEdge ) .def("nodesRightBackEdge" , &SM::Regular::nodesRightBackEdge ) // boundary nodes: edges, without corners .def("nodesFrontBottomOpenEdge" , &SM::Regular::nodesFrontBottomOpenEdge ) .def("nodesFrontTopOpenEdge" , &SM::Regular::nodesFrontTopOpenEdge ) .def("nodesFrontLeftOpenEdge" , &SM::Regular::nodesFrontLeftOpenEdge ) .def("nodesFrontRightOpenEdge" , &SM::Regular::nodesFrontRightOpenEdge ) .def("nodesBackBottomOpenEdge" , &SM::Regular::nodesBackBottomOpenEdge ) .def("nodesBackTopOpenEdge" , &SM::Regular::nodesBackTopOpenEdge ) .def("nodesBackLeftOpenEdge" , &SM::Regular::nodesBackLeftOpenEdge ) .def("nodesBackRightOpenEdge" , &SM::Regular::nodesBackRightOpenEdge ) .def("nodesBottomLeftOpenEdge" , &SM::Regular::nodesBottomLeftOpenEdge ) .def("nodesBottomRightOpenEdge" , &SM::Regular::nodesBottomRightOpenEdge ) .def("nodesTopLeftOpenEdge" , &SM::Regular::nodesTopLeftOpenEdge ) .def("nodesTopRightOpenEdge" , &SM::Regular::nodesTopRightOpenEdge ) // boundary nodes: edges, without corners (aliases) .def("nodesBottomFrontOpenEdge" , &SM::Regular::nodesBottomFrontOpenEdge ) .def("nodesBottomBackOpenEdge" , &SM::Regular::nodesBottomBackOpenEdge ) .def("nodesTopFrontOpenEdge" , &SM::Regular::nodesTopFrontOpenEdge ) .def("nodesTopBackOpenEdge" , &SM::Regular::nodesTopBackOpenEdge ) .def("nodesLeftBottomOpenEdge" , &SM::Regular::nodesLeftBottomOpenEdge ) .def("nodesLeftFrontOpenEdge" , &SM::Regular::nodesLeftFrontOpenEdge ) .def("nodesLeftBackOpenEdge" , &SM::Regular::nodesLeftBackOpenEdge ) .def("nodesLeftTopOpenEdge" , &SM::Regular::nodesLeftTopOpenEdge ) .def("nodesRightBottomOpenEdge" , &SM::Regular::nodesRightBottomOpenEdge ) .def("nodesRightTopOpenEdge" , &SM::Regular::nodesRightTopOpenEdge ) .def("nodesRightFrontOpenEdge" , &SM::Regular::nodesRightFrontOpenEdge ) .def("nodesRightBackOpenEdge" , &SM::Regular::nodesRightBackOpenEdge ) // boundary nodes: corners .def("nodesFrontBottomLeftCorner" , &SM::Regular::nodesFrontBottomLeftCorner ) .def("nodesFrontBottomRightCorner", &SM::Regular::nodesFrontBottomRightCorner) .def("nodesFrontTopLeftCorner" , &SM::Regular::nodesFrontTopLeftCorner ) .def("nodesFrontTopRightCorner" , &SM::Regular::nodesFrontTopRightCorner ) .def("nodesBackBottomLeftCorner" , &SM::Regular::nodesBackBottomLeftCorner ) .def("nodesBackBottomRightCorner" , &SM::Regular::nodesBackBottomRightCorner ) .def("nodesBackTopLeftCorner" , &SM::Regular::nodesBackTopLeftCorner ) .def("nodesBackTopRightCorner" , &SM::Regular::nodesBackTopRightCorner ) // boundary nodes: corners (aliases) .def("nodesFrontLeftBottomCorner" , &SM::Regular::nodesFrontLeftBottomCorner ) .def("nodesBottomFrontLeftCorner" , &SM::Regular::nodesBottomFrontLeftCorner ) .def("nodesBottomLeftFrontCorner" , &SM::Regular::nodesBottomLeftFrontCorner ) .def("nodesLeftFrontBottomCorner" , &SM::Regular::nodesLeftFrontBottomCorner ) .def("nodesLeftBottomFrontCorner" , &SM::Regular::nodesLeftBottomFrontCorner ) .def("nodesFrontRightBottomCorner", &SM::Regular::nodesFrontRightBottomCorner) .def("nodesBottomFrontRightCorner", &SM::Regular::nodesBottomFrontRightCorner) .def("nodesBottomRightFrontCorner", &SM::Regular::nodesBottomRightFrontCorner) .def("nodesRightFrontBottomCorner", &SM::Regular::nodesRightFrontBottomCorner) .def("nodesRightBottomFrontCorner", &SM::Regular::nodesRightBottomFrontCorner) .def("nodesFrontLeftTopCorner" , &SM::Regular::nodesFrontLeftTopCorner ) .def("nodesTopFrontLeftCorner" , &SM::Regular::nodesTopFrontLeftCorner ) .def("nodesTopLeftFrontCorner" , &SM::Regular::nodesTopLeftFrontCorner ) .def("nodesLeftFrontTopCorner" , &SM::Regular::nodesLeftFrontTopCorner ) .def("nodesLeftTopFrontCorner" , &SM::Regular::nodesLeftTopFrontCorner ) .def("nodesFrontRightTopCorner" , &SM::Regular::nodesFrontRightTopCorner ) .def("nodesTopFrontRightCorner" , &SM::Regular::nodesTopFrontRightCorner ) .def("nodesTopRightFrontCorner" , &SM::Regular::nodesTopRightFrontCorner ) .def("nodesRightFrontTopCorner" , &SM::Regular::nodesRightFrontTopCorner ) .def("nodesRightTopFrontCorner" , &SM::Regular::nodesRightTopFrontCorner ) .def("nodesBackLeftBottomCorner" , &SM::Regular::nodesBackLeftBottomCorner ) .def("nodesBottomBackLeftCorner" , &SM::Regular::nodesBottomBackLeftCorner ) .def("nodesBottomLeftBackCorner" , &SM::Regular::nodesBottomLeftBackCorner ) .def("nodesLeftBackBottomCorner" , &SM::Regular::nodesLeftBackBottomCorner ) .def("nodesLeftBottomBackCorner" , &SM::Regular::nodesLeftBottomBackCorner ) .def("nodesBackRightBottomCorner" , &SM::Regular::nodesBackRightBottomCorner ) .def("nodesBottomBackRightCorner" , &SM::Regular::nodesBottomBackRightCorner ) .def("nodesBottomRightBackCorner" , &SM::Regular::nodesBottomRightBackCorner ) .def("nodesRightBackBottomCorner" , &SM::Regular::nodesRightBackBottomCorner ) .def("nodesRightBottomBackCorner" , &SM::Regular::nodesRightBottomBackCorner ) .def("nodesBackLeftTopCorner" , &SM::Regular::nodesBackLeftTopCorner ) .def("nodesTopBackLeftCorner" , &SM::Regular::nodesTopBackLeftCorner ) .def("nodesTopLeftBackCorner" , &SM::Regular::nodesTopLeftBackCorner ) .def("nodesLeftBackTopCorner" , &SM::Regular::nodesLeftBackTopCorner ) .def("nodesLeftTopBackCorner" , &SM::Regular::nodesLeftTopBackCorner ) .def("nodesBackRightTopCorner" , &SM::Regular::nodesBackRightTopCorner ) .def("nodesTopBackRightCorner" , &SM::Regular::nodesTopBackRightCorner ) .def("nodesTopRightBackCorner" , &SM::Regular::nodesTopRightBackCorner ) .def("nodesRightBackTopCorner" , &SM::Regular::nodesRightBackTopCorner ) .def("nodesRightTopBackCorner" , &SM::Regular::nodesRightTopBackCorner ) // periodicity .def("nodesPeriodic" , &SM::Regular::nodesPeriodic ) .def("nodesOrigin" , &SM::Regular::nodesOrigin ) .def("dofs" , &SM::Regular::dofs ) .def("dofsPeriodic" , &SM::Regular::dofsPeriodic ) // print to screen .def("__repr__", - [](const SM::Regular &a){ return ""; } + [](const SM::Regular &){ return ""; } ); // ------------------------------------------------------------------------------------------------- py::class_(sm, "FineLayer") // constructor .def( py::init(), "mesh with nx*ny*nz 'pixels' and edge size h", py::arg("nx"), py::arg("ny"), py::arg("nz"), py::arg("h")=1., py::arg("nfine")=1 ) // sizes .def("nelem" , &SM::FineLayer::nelem ) .def("nnode" , &SM::FineLayer::nnode ) .def("nne" , &SM::FineLayer::nne ) .def("ndim" , &SM::FineLayer::ndim ) .def("shape" , &SM::FineLayer::shape ) // mesh .def("coor" , &SM::FineLayer::coor ) .def("conn" , &SM::FineLayer::conn ) // element sets .def("elementsMiddleLayer" , &SM::FineLayer::elementsMiddleLayer ) // boundary nodes: planes .def("nodesFront" , &SM::FineLayer::nodesFront ) .def("nodesBack" , &SM::FineLayer::nodesBack ) .def("nodesLeft" , &SM::FineLayer::nodesLeft ) .def("nodesRight" , &SM::FineLayer::nodesRight ) .def("nodesBottom" , &SM::FineLayer::nodesBottom ) .def("nodesTop" , &SM::FineLayer::nodesTop ) // boundary nodes: faces .def("nodesFrontFace" , &SM::FineLayer::nodesFrontFace ) .def("nodesBackFace" , &SM::FineLayer::nodesBackFace ) .def("nodesLeftFace" , &SM::FineLayer::nodesLeftFace ) .def("nodesRightFace" , &SM::FineLayer::nodesRightFace ) .def("nodesBottomFace" , &SM::FineLayer::nodesBottomFace ) .def("nodesTopFace" , &SM::FineLayer::nodesTopFace ) // boundary nodes: edges .def("nodesFrontBottomEdge" , &SM::FineLayer::nodesFrontBottomEdge ) .def("nodesFrontTopEdge" , &SM::FineLayer::nodesFrontTopEdge ) .def("nodesFrontLeftEdge" , &SM::FineLayer::nodesFrontLeftEdge ) .def("nodesFrontRightEdge" , &SM::FineLayer::nodesFrontRightEdge ) .def("nodesBackBottomEdge" , &SM::FineLayer::nodesBackBottomEdge ) .def("nodesBackTopEdge" , &SM::FineLayer::nodesBackTopEdge ) .def("nodesBackLeftEdge" , &SM::FineLayer::nodesBackLeftEdge ) .def("nodesBackRightEdge" , &SM::FineLayer::nodesBackRightEdge ) .def("nodesBottomLeftEdge" , &SM::FineLayer::nodesBottomLeftEdge ) .def("nodesBottomRightEdge" , &SM::FineLayer::nodesBottomRightEdge ) .def("nodesTopLeftEdge" , &SM::FineLayer::nodesTopLeftEdge ) .def("nodesTopRightEdge" , &SM::FineLayer::nodesTopRightEdge ) // boundary nodes: faces (aliases) .def("nodesBottomFrontEdge" , &SM::FineLayer::nodesBottomFrontEdge ) .def("nodesBottomBackEdge" , &SM::FineLayer::nodesBottomBackEdge ) .def("nodesTopFrontEdge" , &SM::FineLayer::nodesTopFrontEdge ) .def("nodesTopBackEdge" , &SM::FineLayer::nodesTopBackEdge ) .def("nodesLeftBottomEdge" , &SM::FineLayer::nodesLeftBottomEdge ) .def("nodesLeftFrontEdge" , &SM::FineLayer::nodesLeftFrontEdge ) .def("nodesLeftBackEdge" , &SM::FineLayer::nodesLeftBackEdge ) .def("nodesLeftTopEdge" , &SM::FineLayer::nodesLeftTopEdge ) .def("nodesRightBottomEdge" , &SM::FineLayer::nodesRightBottomEdge ) .def("nodesRightTopEdge" , &SM::FineLayer::nodesRightTopEdge ) .def("nodesRightFrontEdge" , &SM::FineLayer::nodesRightFrontEdge ) .def("nodesRightBackEdge" , &SM::FineLayer::nodesRightBackEdge ) // boundary nodes: edges, without corners .def("nodesFrontBottomOpenEdge" , &SM::FineLayer::nodesFrontBottomOpenEdge ) .def("nodesFrontTopOpenEdge" , &SM::FineLayer::nodesFrontTopOpenEdge ) .def("nodesFrontLeftOpenEdge" , &SM::FineLayer::nodesFrontLeftOpenEdge ) .def("nodesFrontRightOpenEdge" , &SM::FineLayer::nodesFrontRightOpenEdge ) .def("nodesBackBottomOpenEdge" , &SM::FineLayer::nodesBackBottomOpenEdge ) .def("nodesBackTopOpenEdge" , &SM::FineLayer::nodesBackTopOpenEdge ) .def("nodesBackLeftOpenEdge" , &SM::FineLayer::nodesBackLeftOpenEdge ) .def("nodesBackRightOpenEdge" , &SM::FineLayer::nodesBackRightOpenEdge ) .def("nodesBottomLeftOpenEdge" , &SM::FineLayer::nodesBottomLeftOpenEdge ) .def("nodesBottomRightOpenEdge" , &SM::FineLayer::nodesBottomRightOpenEdge ) .def("nodesTopLeftOpenEdge" , &SM::FineLayer::nodesTopLeftOpenEdge ) .def("nodesTopRightOpenEdge" , &SM::FineLayer::nodesTopRightOpenEdge ) // boundary nodes: edges, without corners (aliases) .def("nodesBottomFrontOpenEdge" , &SM::FineLayer::nodesBottomFrontOpenEdge ) .def("nodesBottomBackOpenEdge" , &SM::FineLayer::nodesBottomBackOpenEdge ) .def("nodesTopFrontOpenEdge" , &SM::FineLayer::nodesTopFrontOpenEdge ) .def("nodesTopBackOpenEdge" , &SM::FineLayer::nodesTopBackOpenEdge ) .def("nodesLeftBottomOpenEdge" , &SM::FineLayer::nodesLeftBottomOpenEdge ) .def("nodesLeftFrontOpenEdge" , &SM::FineLayer::nodesLeftFrontOpenEdge ) .def("nodesLeftBackOpenEdge" , &SM::FineLayer::nodesLeftBackOpenEdge ) .def("nodesLeftTopOpenEdge" , &SM::FineLayer::nodesLeftTopOpenEdge ) .def("nodesRightBottomOpenEdge" , &SM::FineLayer::nodesRightBottomOpenEdge ) .def("nodesRightTopOpenEdge" , &SM::FineLayer::nodesRightTopOpenEdge ) .def("nodesRightFrontOpenEdge" , &SM::FineLayer::nodesRightFrontOpenEdge ) .def("nodesRightBackOpenEdge" , &SM::FineLayer::nodesRightBackOpenEdge ) // boundary nodes: corners .def("nodesFrontBottomLeftCorner" , &SM::FineLayer::nodesFrontBottomLeftCorner ) .def("nodesFrontBottomRightCorner", &SM::FineLayer::nodesFrontBottomRightCorner) .def("nodesFrontTopLeftCorner" , &SM::FineLayer::nodesFrontTopLeftCorner ) .def("nodesFrontTopRightCorner" , &SM::FineLayer::nodesFrontTopRightCorner ) .def("nodesBackBottomLeftCorner" , &SM::FineLayer::nodesBackBottomLeftCorner ) .def("nodesBackBottomRightCorner" , &SM::FineLayer::nodesBackBottomRightCorner ) .def("nodesBackTopLeftCorner" , &SM::FineLayer::nodesBackTopLeftCorner ) .def("nodesBackTopRightCorner" , &SM::FineLayer::nodesBackTopRightCorner ) // boundary nodes: corners (aliases) .def("nodesFrontLeftBottomCorner" , &SM::FineLayer::nodesFrontLeftBottomCorner ) .def("nodesBottomFrontLeftCorner" , &SM::FineLayer::nodesBottomFrontLeftCorner ) .def("nodesBottomLeftFrontCorner" , &SM::FineLayer::nodesBottomLeftFrontCorner ) .def("nodesLeftFrontBottomCorner" , &SM::FineLayer::nodesLeftFrontBottomCorner ) .def("nodesLeftBottomFrontCorner" , &SM::FineLayer::nodesLeftBottomFrontCorner ) .def("nodesFrontRightBottomCorner", &SM::FineLayer::nodesFrontRightBottomCorner) .def("nodesBottomFrontRightCorner", &SM::FineLayer::nodesBottomFrontRightCorner) .def("nodesBottomRightFrontCorner", &SM::FineLayer::nodesBottomRightFrontCorner) .def("nodesRightFrontBottomCorner", &SM::FineLayer::nodesRightFrontBottomCorner) .def("nodesRightBottomFrontCorner", &SM::FineLayer::nodesRightBottomFrontCorner) .def("nodesFrontLeftTopCorner" , &SM::FineLayer::nodesFrontLeftTopCorner ) .def("nodesTopFrontLeftCorner" , &SM::FineLayer::nodesTopFrontLeftCorner ) .def("nodesTopLeftFrontCorner" , &SM::FineLayer::nodesTopLeftFrontCorner ) .def("nodesLeftFrontTopCorner" , &SM::FineLayer::nodesLeftFrontTopCorner ) .def("nodesLeftTopFrontCorner" , &SM::FineLayer::nodesLeftTopFrontCorner ) .def("nodesFrontRightTopCorner" , &SM::FineLayer::nodesFrontRightTopCorner ) .def("nodesTopFrontRightCorner" , &SM::FineLayer::nodesTopFrontRightCorner ) .def("nodesTopRightFrontCorner" , &SM::FineLayer::nodesTopRightFrontCorner ) .def("nodesRightFrontTopCorner" , &SM::FineLayer::nodesRightFrontTopCorner ) .def("nodesRightTopFrontCorner" , &SM::FineLayer::nodesRightTopFrontCorner ) .def("nodesBackLeftBottomCorner" , &SM::FineLayer::nodesBackLeftBottomCorner ) .def("nodesBottomBackLeftCorner" , &SM::FineLayer::nodesBottomBackLeftCorner ) .def("nodesBottomLeftBackCorner" , &SM::FineLayer::nodesBottomLeftBackCorner ) .def("nodesLeftBackBottomCorner" , &SM::FineLayer::nodesLeftBackBottomCorner ) .def("nodesLeftBottomBackCorner" , &SM::FineLayer::nodesLeftBottomBackCorner ) .def("nodesBackRightBottomCorner" , &SM::FineLayer::nodesBackRightBottomCorner ) .def("nodesBottomBackRightCorner" , &SM::FineLayer::nodesBottomBackRightCorner ) .def("nodesBottomRightBackCorner" , &SM::FineLayer::nodesBottomRightBackCorner ) .def("nodesRightBackBottomCorner" , &SM::FineLayer::nodesRightBackBottomCorner ) .def("nodesRightBottomBackCorner" , &SM::FineLayer::nodesRightBottomBackCorner ) .def("nodesBackLeftTopCorner" , &SM::FineLayer::nodesBackLeftTopCorner ) .def("nodesTopBackLeftCorner" , &SM::FineLayer::nodesTopBackLeftCorner ) .def("nodesTopLeftBackCorner" , &SM::FineLayer::nodesTopLeftBackCorner ) .def("nodesLeftBackTopCorner" , &SM::FineLayer::nodesLeftBackTopCorner ) .def("nodesLeftTopBackCorner" , &SM::FineLayer::nodesLeftTopBackCorner ) .def("nodesBackRightTopCorner" , &SM::FineLayer::nodesBackRightTopCorner ) .def("nodesTopBackRightCorner" , &SM::FineLayer::nodesTopBackRightCorner ) .def("nodesTopRightBackCorner" , &SM::FineLayer::nodesTopRightBackCorner ) .def("nodesRightBackTopCorner" , &SM::FineLayer::nodesRightBackTopCorner ) .def("nodesRightTopBackCorner" , &SM::FineLayer::nodesRightTopBackCorner ) // periodicity .def("nodesPeriodic" , &SM::FineLayer::nodesPeriodic ) .def("nodesOrigin" , &SM::FineLayer::nodesOrigin ) .def("dofs" , &SM::FineLayer::dofs ) .def("dofsPeriodic" , &SM::FineLayer::dofsPeriodic ) // print to screen .def("__repr__", - [](const SM::FineLayer &a){ return ""; } + [](const SM::FineLayer &){ return ""; } ); // ------------------------------------------------------------------------------------------------- } // ========================= GooseFEM::Mesh::Quad4 - GooseFEM/MeshQuad4.h ========================== { // create sub-module py::module sm = mMesh.def_submodule("Quad4", "Linear quadrilateral elements (2D)"); // abbreviate name-space namespace SM = GooseFEM::Mesh::Quad4; // ------------------------------------------------------------------------------------------------- py::class_(sm, "Regular") .def( py::init(), "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge size 'h'", py::arg("nx"), py::arg("ny"), py::arg("h")=1. ) .def("coor" , &SM::Regular::coor ) .def("conn" , &SM::Regular::conn ) .def("nelem" , &SM::Regular::nelem ) .def("nnode" , &SM::Regular::nnode ) .def("nne" , &SM::Regular::nne ) .def("ndim" , &SM::Regular::ndim ) .def("nodesBottomEdge" , &SM::Regular::nodesBottomEdge ) .def("nodesTopEdge" , &SM::Regular::nodesTopEdge ) .def("nodesLeftEdge" , &SM::Regular::nodesLeftEdge ) .def("nodesRightEdge" , &SM::Regular::nodesRightEdge ) .def("nodesBottomOpenEdge" , &SM::Regular::nodesBottomOpenEdge ) .def("nodesTopOpenEdge" , &SM::Regular::nodesTopOpenEdge ) .def("nodesLeftOpenEdge" , &SM::Regular::nodesLeftOpenEdge ) .def("nodesRightOpenEdge" , &SM::Regular::nodesRightOpenEdge ) .def("nodesBottomLeftCorner" , &SM::Regular::nodesBottomLeftCorner ) .def("nodesBottomRightCorner", &SM::Regular::nodesBottomRightCorner) .def("nodesTopLeftCorner" , &SM::Regular::nodesTopLeftCorner ) .def("nodesTopRightCorner" , &SM::Regular::nodesTopRightCorner ) .def("nodesLeftBottomCorner" , &SM::Regular::nodesLeftBottomCorner ) .def("nodesLeftTopCorner" , &SM::Regular::nodesLeftTopCorner ) .def("nodesRightBottomCorner", &SM::Regular::nodesRightBottomCorner) .def("nodesRightTopCorner" , &SM::Regular::nodesRightTopCorner ) .def("nodesPeriodic" , &SM::Regular::nodesPeriodic ) .def("nodesOrigin" , &SM::Regular::nodesOrigin ) .def("dofs" , &SM::Regular::dofs ) .def("dofsPeriodic" , &SM::Regular::dofsPeriodic ) .def("__repr__", - [](const SM::Regular &a){ return ""; } + [](const SM::Regular &){ return ""; } ); // ------------------------------------------------------------------------------------------------- py::class_(sm, "FineLayer") .def( py::init(), "FineLayer mesh: 'nx' pixels in horizontal direction (length 'Lx'), idem in vertical direction", py::arg("nx"), py::arg("ny"), py::arg("h")=1., py::arg("nfine")=1 ) .def("shape" , &SM::FineLayer::shape ) .def("coor" , &SM::FineLayer::coor ) .def("conn" , &SM::FineLayer::conn ) .def("nelem" , &SM::FineLayer::nelem ) .def("nnode" , &SM::FineLayer::nnode ) .def("nne" , &SM::FineLayer::nne ) .def("ndim" , &SM::FineLayer::ndim ) .def("elementsMiddleLayer" , &SM::FineLayer::elementsMiddleLayer ) .def("nodesBottomEdge" , &SM::FineLayer::nodesBottomEdge ) .def("nodesTopEdge" , &SM::FineLayer::nodesTopEdge ) .def("nodesLeftEdge" , &SM::FineLayer::nodesLeftEdge ) .def("nodesRightEdge" , &SM::FineLayer::nodesRightEdge ) .def("nodesBottomOpenEdge" , &SM::FineLayer::nodesBottomOpenEdge ) .def("nodesTopOpenEdge" , &SM::FineLayer::nodesTopOpenEdge ) .def("nodesLeftOpenEdge" , &SM::FineLayer::nodesLeftOpenEdge ) .def("nodesRightOpenEdge" , &SM::FineLayer::nodesRightOpenEdge ) .def("nodesBottomLeftCorner" , &SM::FineLayer::nodesBottomLeftCorner ) .def("nodesBottomRightCorner", &SM::FineLayer::nodesBottomRightCorner) .def("nodesTopLeftCorner" , &SM::FineLayer::nodesTopLeftCorner ) .def("nodesTopRightCorner" , &SM::FineLayer::nodesTopRightCorner ) .def("nodesLeftBottomCorner" , &SM::FineLayer::nodesLeftBottomCorner ) .def("nodesLeftTopCorner" , &SM::FineLayer::nodesLeftTopCorner ) .def("nodesRightBottomCorner", &SM::FineLayer::nodesRightBottomCorner) .def("nodesRightTopCorner" , &SM::FineLayer::nodesRightTopCorner ) .def("nodesPeriodic" , &SM::FineLayer::nodesPeriodic ) .def("nodesOrigin" , &SM::FineLayer::nodesOrigin ) .def("dofs" , &SM::FineLayer::dofs ) .def("dofsPeriodic" , &SM::FineLayer::dofsPeriodic ) .def("__repr__", - [](const SM::FineLayer &a){ return ""; } + [](const SM::FineLayer &){ return ""; } ); // ------------------------------------------------------------------------------------------------- } // ========================== GooseFEM::Mesh::Tri3 - GooseFEM/MeshTri3.h =========================== { // create sub-module py::module sm = mMesh.def_submodule("Tri3" , "Linear triangular elements (2D)"); // abbreviate name-space namespace SM = GooseFEM::Mesh::Tri3; // ------------------------------------------------------------------------------------------------- py::class_(sm, "Regular") .def( py::init(), "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge size 'h'", py::arg("nx"), py::arg("ny"), py::arg("h")=1. ) .def("coor" , &SM::Regular::coor ) .def("conn" , &SM::Regular::conn ) .def("nelem" , &SM::Regular::nelem ) .def("nnode" , &SM::Regular::nnode ) .def("nne" , &SM::Regular::nne ) .def("ndim" , &SM::Regular::ndim ) .def("nodesBottomEdge" , &SM::Regular::nodesBottomEdge ) .def("nodesTopEdge" , &SM::Regular::nodesTopEdge ) .def("nodesLeftEdge" , &SM::Regular::nodesLeftEdge ) .def("nodesRightEdge" , &SM::Regular::nodesRightEdge ) .def("nodesBottomOpenEdge" , &SM::Regular::nodesBottomOpenEdge ) .def("nodesTopOpenEdge" , &SM::Regular::nodesTopOpenEdge ) .def("nodesLeftOpenEdge" , &SM::Regular::nodesLeftOpenEdge ) .def("nodesRightOpenEdge" , &SM::Regular::nodesRightOpenEdge ) .def("nodesBottomLeftCorner" , &SM::Regular::nodesBottomLeftCorner ) .def("nodesBottomRightCorner", &SM::Regular::nodesBottomRightCorner) .def("nodesTopLeftCorner" , &SM::Regular::nodesTopLeftCorner ) .def("nodesTopRightCorner" , &SM::Regular::nodesTopRightCorner ) .def("nodesLeftBottomCorner" , &SM::Regular::nodesLeftBottomCorner ) .def("nodesLeftTopCorner" , &SM::Regular::nodesLeftTopCorner ) .def("nodesRightBottomCorner", &SM::Regular::nodesRightBottomCorner) .def("nodesRightTopCorner" , &SM::Regular::nodesRightTopCorner ) .def("nodesPeriodic" , &SM::Regular::nodesPeriodic ) .def("nodesOrigin" , &SM::Regular::nodesOrigin ) .def("dofs" , &SM::Regular::dofs ) .def("dofsPeriodic" , &SM::Regular::dofsPeriodic ) .def("__repr__", - [](const SM::Regular &a){ return ""; } + [](const SM::Regular &){ return ""; } ); // ------------------------------------------------------------------------------------------------- sm.def("getOrientation",&SM::getOrientation, "Get the orientation of each element", py::arg("coor"), py::arg("conn") ); // ------------------------------------------------------------------------------------------------- sm.def("retriangulate",&SM::retriangulate, "Re-triangulate existing mesh", py::arg("coor"), py::arg("conn"), py::arg("orientation")=-1 ); // ------------------------------------------------------------------------------------------------- } // ================================================================================================= }