diff --git a/src/GooseFEM/DynamicsDiagonalPeriodic.h b/src/GooseFEM/DynamicsDiagonalPeriodic.h index 92e6445..b77bbc3 100644 --- a/src/GooseFEM/DynamicsDiagonalPeriodic.h +++ b/src/GooseFEM/DynamicsDiagonalPeriodic.h @@ -1,360 +1,360 @@ /* ========================================== DESCRIPTION ========================================== (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_DYNAMICS_DIAGONAL_PERIODIC_H #define GOOSEFEM_DYNAMICS_DIAGONAL_PERIODIC_H #include "Macros.h" #include // ------------------------------------------------------------------------------------------------- namespace GooseFEM { namespace Dynamics { namespace Diagonal { // =========================================== OVERVIEW ============================================ -template +template class Periodic { public: // variables // --------- // element/quadrature/material definition std::unique_ptr elem; // mesh: nodal position/displacement/velocity/acceleration, DOF-numbers, connectivity, dimensions size_t nnode, nelem, nne, ndim, ndof; MatS conn, dofs; MatD x, u, v, a; // linear system: columns (also "M" and "D" which are diagonal) ColD M, Minv, D, F, V, V_n, A, A_n; // time integration double dt, t=0.0; // constructor // ----------- Periodic(){}; Periodic( std::unique_ptr elem, const MatD &x, const MatS &conn, const MatS &dofs, double dt ); // functions // --------- void velocityVerlet(); // one time step of time integrator void Verlet(); // one time step of time integrator (no velocity dependence) void updated_x(); // process update in "x" void updated_u(bool init=false); // process update in "u", if init all possible updates are made void updated_v(bool init=false); // process update in "v", if init all possible updates are made void assemble_M(); // assemble the mass matrix from the element mass matrices void assemble_D(); // assemble the damping matrix from the element damping matrices void assemble_F(); // assemble the force from the element forces }; // ========================================== SOURCE CODE ========================================== template -Periodic::Periodic( +inline Periodic::Periodic( std::unique_ptr _elem, const MatD &_x, const MatS &_conn, const MatS &_dofs, double _dt ) { // copy input elem = std::move(_elem); x = _x; conn = _conn; dofs = _dofs; dt = _dt; // compute sizes nnode = static_cast(x.rows()); ndim = static_cast(x.cols()); nelem = static_cast(conn.rows()); nne = static_cast(conn.cols()); ndof = dofs.maxCoeff()+1; // basic checks (mostly the user is 'trusted') assert( dofs.size() == nnode * ndim ); assert( ndof < nnode * ndim ); // allocate and zero-initialize nodal quantities u.conservativeResize(nnode,ndim); u.setZero(); v.conservativeResize(nnode,ndim); v.setZero(); a.conservativeResize(nnode,ndim); a.setZero(); // allocate and zero-initialize linear system (DOFs) Minv.conservativeResize(ndof); M .conservativeResize(ndof); D .conservativeResize(ndof); F .conservativeResize(ndof); V .conservativeResize(ndof); V_n .conservativeResize(ndof); V_n .setZero(); A .conservativeResize(ndof); A_n .conservativeResize(ndof); A_n .setZero(); // initialize all fields updated_x(); updated_u(true); updated_v(true); } // ================================================================================================= template -void Periodic::velocityVerlet() +inline void Periodic::velocityVerlet() { // (1) new positions (displacements) // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" u += dt * v + ( .5 * std::pow(dt,2.) ) * a; // - process update in displacements updated_u(); // (2a) estimate new velocities // - update velocities (DOFs) V.noalias() = V_n + dt * A; // - convert to nodal velocities (periodicity implies that several nodes depend on the same DOF) for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' // "v" == "V_n" -> new nodal velocities, and a 'back-up' // "u" -> new nodal displacements // The forces "F" correspond to this state of the system } // ================================================================================================= template -void Periodic::Verlet() +inline void Periodic::Verlet() { // (1) new nodal positions (displacements) // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" u += dt * v + ( .5 * std::pow(dt,2.) ) * a; // - process update in displacements updated_u(); // - solve for accelerations (DOFs) A.noalias() = Minv.cwiseProduct( - F ); // - convert to nodal acceleration (periodicity implies that several nodes depend on the same DOF) for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' // "v" == "V_n" -> new nodal velocities, and a 'back-up' // "u" -> new nodal displacements // The forces "F" correspond to this state of the system } // ================================================================================================= template -void Periodic::updated_x() +inline void Periodic::updated_x() { // set the nodal positions of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->x(e,m,i) = x(conn(e,m),i); // signal update elem->updated_x(); } // ================================================================================================= template -void Periodic::updated_u(bool init) +inline void Periodic::updated_u(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->u(e,m,i) = u(conn(e,m),i); // signal update elem->updated_u(); // update if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ================================================================================================= template -void Periodic::updated_v(bool init) +inline void Periodic::updated_v(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->v(e,m,i) = v(conn(e,m),i); // signal update elem->updated_v(); // update if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ================================================================================================= template -void Periodic::assemble_M() +inline void Periodic::assemble_M() { // zero-initialize M.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - mass matrix, per thread ColD M_(ndof); M_.setZero(); // - assemble diagonal mass matrix, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) M_(dofs(conn(e,m),i)) += elem->M(e,m*ndim+i,m*ndim+i); // - reduce "M_" per thread to total "M" #pragma omp critical M += M_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); // compute inverse of the mass matrix Minv = M.cwiseInverse(); } // ================================================================================================= template -void Periodic::assemble_D() +inline void Periodic::assemble_D() { // zero-initialize D.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - damping matrix, per thread ColD D_(ndof); D_.setZero(); // - assemble diagonal damping matrix, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) D_(dofs(conn(e,m),i)) += elem->D(e,m*ndim+i,m*ndim+i); // - reduce "D_" per thread to total "D" #pragma omp critical D += D_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); } // ================================================================================================= template -void Periodic::assemble_F() +inline void Periodic::assemble_F() { // zero-initialize F.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - force, per thread ColD F_(ndof); F_.setZero(); // - assemble force, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) F_(dofs(conn(e,m),i)) += elem->f(e,m,i); // - reduce "F_" per thread to total "F" #pragma omp critical F += F_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); } // ================================================================================================= }}} // namespace GooseFEM::Dynamics::Diagonal #endif diff --git a/src/GooseFEM/DynamicsDiagonalSemiPeriodic.h b/src/GooseFEM/DynamicsDiagonalSemiPeriodic.h index 6ccbf08..93edbab 100644 --- a/src/GooseFEM/DynamicsDiagonalSemiPeriodic.h +++ b/src/GooseFEM/DynamicsDiagonalSemiPeriodic.h @@ -1,389 +1,389 @@ /* ========================================== DESCRIPTION ========================================== (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_DYNAMICS_DIAGONAL_SEMIPERIODIC_H #define GOOSEFEM_DYNAMICS_DIAGONAL_SEMIPERIODIC_H #include "Macros.h" #include // ------------------------------------------------------------------------------------------------- namespace GooseFEM { namespace Dynamics { namespace Diagonal { // =========================================== OVERVIEW ============================================ -template +template class SemiPeriodic { public: // variables // --------- // element/quadrature/material definition std::unique_ptr elem; // mesh: nodal position/displacement/velocity/acceleration, DOF-numbers, connectivity, dimensions size_t nnode, nelem, nne, ndim, ndof; MatS conn, dofs; MatD x, u, v, a; // fixed DOFs: prescribed velocity/acceleration, DOF-numbers, dimensions size_t nfixed; ColS fixedDofs; ColD fixedV, fixedA; // linear system: columns (also "M" and "D" which are diagonal) ColD M, Minv, D, F, V, V_n, A, A_n; // time integration double dt, t=0.0; // constructor // ----------- SemiPeriodic(){}; SemiPeriodic( std::unique_ptr elem, const MatD &x, const MatS &conn, const MatS &dofs, const ColS &fixedDofs, double dt ); // functions // --------- void velocityVerlet(); // one time step of time integrator void Verlet(); // one time step of time integrator (no velocity dependence) void updated_x(); // process update in "x" void updated_u(bool init=false); // process update in "u", if init all possible updates are made void updated_v(bool init=false); // process update in "v", if init all possible updates are made void assemble_M(); // assemble the mass matrix from the element mass matrices void assemble_D(); // assemble the damping matrix from the element damping matrices void assemble_F(); // assemble the force from the element forces }; // ========================================== SOURCE CODE ========================================== template -SemiPeriodic::SemiPeriodic( +inline SemiPeriodic::SemiPeriodic( std::unique_ptr _elem, const MatD &_x, const MatS &_conn, const MatS &_dofs, const ColS &_fixedDofs, double _dt ) { // copy input elem = std::move(_elem); x = _x; conn = _conn; dofs = _dofs; dt = _dt; fixedDofs = _fixedDofs; // compute sizes nfixed = static_cast(fixedDofs.size()); nnode = static_cast(x.rows()); ndim = static_cast(x.cols()); nelem = static_cast(conn.rows()); nne = static_cast(conn.cols()); ndof = dofs.maxCoeff()+1; // basic checks (mostly the user is 'trusted') assert( dofs.size() == nnode * ndim ); assert( ndof < nnode * ndim ); #ifndef NDEBUG for ( size_t i = 0 ; i < nfixed ; ++i ) assert( fixedDofs(i) < ndof ); #endif // allocate and zero-initialize nodal quantities u.conservativeResize(nnode,ndim); u.setZero(); v.conservativeResize(nnode,ndim); v.setZero(); a.conservativeResize(nnode,ndim); a.setZero(); // allocate and zero-initialize linear system (DOFs) Minv.conservativeResize(ndof); M .conservativeResize(ndof); D .conservativeResize(ndof); F .conservativeResize(ndof); V .conservativeResize(ndof); V_n .conservativeResize(ndof); V_n .setZero(); A .conservativeResize(ndof); A_n .conservativeResize(ndof); A_n .setZero(); // fixed DOFs : default zero velocity and acceleration fixedV.conservativeResize(nfixed); fixedV.setZero(); fixedA.conservativeResize(nfixed); fixedA.setZero(); // initialize all fields updated_x(); updated_u(true); updated_v(true); } // ================================================================================================= template -void SemiPeriodic::velocityVerlet() +inline void SemiPeriodic::velocityVerlet() { // (1) new positions (displacements) // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" u += dt * v + ( .5 * std::pow(dt,2.) ) * a; // - process update in displacements updated_u(); // (2a) estimate new velocities // - update velocities (DOFs) V.noalias() = V_n + dt * A; // - apply the fixed velocities for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' // "v" == "V_n" -> new nodal velocities, and a 'back-up' // "u" -> new nodal displacements // The forces "F" correspond to this state of the system } // ================================================================================================= template -void SemiPeriodic::Verlet() +inline void SemiPeriodic::Verlet() { // (1) new nodal positions (displacements) // - apply update (nodal) : x_{n+1} = x_n + dt * v_n + .5 * dt^2 * a_n" u += dt * v + ( .5 * std::pow(dt,2.) ) * a; // - process update in displacements updated_u(); // - solve for accelerations (DOFs) A.noalias() = Minv.cwiseProduct( - F ); // - apply the fixed accelerations for ( size_t i=0; i new nodal accelerations, their DOF equivalents, and a 'back-up' // "v" == "V_n" -> new nodal velocities, and a 'back-up' // "u" -> new nodal displacements // The forces "F" correspond to this state of the system } // ================================================================================================= template -void SemiPeriodic::updated_x() +inline void SemiPeriodic::updated_x() { // set the nodal positions of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->x(e,m,i) = x(conn(e,m),i); // signal update elem->updated_x(); } // ================================================================================================= template -void SemiPeriodic::updated_u(bool init) +inline void SemiPeriodic::updated_u(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->u(e,m,i) = u(conn(e,m),i); // signal update elem->updated_u(); // update if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ================================================================================================= template -void SemiPeriodic::updated_v(bool init) +inline void SemiPeriodic::updated_v(bool init) { // set the nodal displacements of all elements (in parallel) #pragma omp parallel for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) elem->v(e,m,i) = v(conn(e,m),i); // signal update elem->updated_v(); // update if ( elem->changed_M or init ) assemble_M(); if ( elem->changed_D or init ) assemble_D(); if ( elem->changed_f or init ) assemble_F(); } // ================================================================================================= template -void SemiPeriodic::assemble_M() +inline void SemiPeriodic::assemble_M() { // zero-initialize M.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - mass matrix, per thread ColD M_(ndof); M_.setZero(); // - assemble diagonal mass matrix, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) M_(dofs(conn(e,m),i)) += elem->M(e,m*ndim+i,m*ndim+i); // - reduce "M_" per thread to total "M" #pragma omp critical M += M_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); // compute inverse of the mass matrix Minv = M.cwiseInverse(); } // ================================================================================================= template -void SemiPeriodic::assemble_D() +inline void SemiPeriodic::assemble_D() { // zero-initialize D.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - damping matrix, per thread ColD D_(ndof); D_.setZero(); // - assemble diagonal damping matrix, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) D_(dofs(conn(e,m),i)) += elem->D(e,m*ndim+i,m*ndim+i); // - reduce "D_" per thread to total "D" #pragma omp critical D += D_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); } // ================================================================================================= template -void SemiPeriodic::assemble_F() +inline void SemiPeriodic::assemble_F() { // zero-initialize F.setZero(); // temporarily disable parallelization by Eigen Eigen::setNbThreads(1); // assemble #pragma omp parallel { // - force, per thread ColD F_(ndof); F_.setZero(); // - assemble force, per thread #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) for ( size_t m = 0 ; m < nne ; ++m ) for ( size_t i = 0 ; i < ndim ; ++i ) F_(dofs(conn(e,m),i)) += elem->f(e,m,i); // - reduce "F_" per thread to total "F" #pragma omp critical F += F_; } // automatic parallelization by Eigen Eigen::setNbThreads(0); } // ================================================================================================= }}} // namespace GooseFEM::Dynamics::Diagonal #endif diff --git a/src/GooseFEM/DynamicsDiagonalSmallStrainQuad4.h b/src/GooseFEM/DynamicsDiagonalSmallStrainQuad4.h index 5cb4a4d..838505b 100644 --- a/src/GooseFEM/DynamicsDiagonalSmallStrainQuad4.h +++ b/src/GooseFEM/DynamicsDiagonalSmallStrainQuad4.h @@ -1,416 +1,416 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_DYNAMICS_DIAGONAL_SMALLSTRAIN_QUAD4_H #define GOOSEFEM_DYNAMICS_DIAGONAL_SMALLSTRAIN_QUAD4_H // ------------------------------------------------------------------------------------------------- #include "Macros.h" #include // ------------------------------------------------------------------------------------------------- namespace GooseFEM { namespace Element { namespace Diagonal { namespace SmallStrain { // ========================== N.B. most loops are unrolled for efficiency ========================== -template +template class Quad4 { public: // class variables // --------------- // constitutive response per integration point, per element std::unique_ptr mat; // matrices to store the element data cppmat::matrix x, u, v, f, M, D, dNx, dNxi, w, V, dNxi_n, w_n, V_n; // dimensions size_t nelem, nip=4, nne=4, ndim=2; // signals for solvers bool changed_f=false, changed_M=true, changed_D=true; // class functions // --------------- // constructor Quad4(std::unique_ptr mat, size_t nelem); // recompute relevant quantities after "x", "u", or "v" have been externally updated void updated_x(); void updated_u(); void updated_v(); }; // ================================================================================================= template -Quad4::Quad4(std::unique_ptr _mat, size_t _nelem) +inline Quad4::Quad4(std::unique_ptr _mat, size_t _nelem) { // copy from input nelem = _nelem; mat = std::move(_mat); // allocate matrices // - x .resize({nelem,nne,ndim}); u .resize({nelem,nne,ndim}); v .resize({nelem,nne,ndim}); f .resize({nelem,nne,ndim}); // - M .resize({nelem,nne*ndim,nne*ndim}); D .resize({nelem,nne*ndim,nne*ndim}); // - dNx .resize({nelem,nip,nne,ndim}); // - V .resize({nelem,nip}); V_n .resize({nelem,nne}); // - dNxi .resize({nip,nne,ndim}); dNxi_n.resize({nne,nne,ndim}); // - w .resize({nip}); w_n .resize({nne}); // zero initialize matrices (only diagonal is written, no new zero-initialization necessary) M.setZero(); D.setZero(); // shape function gradient at all Gauss points, in local coordinates // - k == 0 dNxi(0,0,0) = -.25*(1.+1./std::sqrt(3.)); dNxi(0,0,1) = -.25*(1.+1./std::sqrt(3.)); dNxi(0,1,0) = +.25*(1.+1./std::sqrt(3.)); dNxi(0,1,1) = -.25*(1.-1./std::sqrt(3.)); dNxi(0,2,0) = +.25*(1.-1./std::sqrt(3.)); dNxi(0,2,1) = +.25*(1.-1./std::sqrt(3.)); dNxi(0,3,0) = -.25*(1.-1./std::sqrt(3.)); dNxi(0,3,1) = +.25*(1.+1./std::sqrt(3.)); // - k == 1 dNxi(1,0,0) = -.25*(1.+1./std::sqrt(3.)); dNxi(1,0,1) = -.25*(1.-1./std::sqrt(3.)); dNxi(1,1,0) = +.25*(1.+1./std::sqrt(3.)); dNxi(1,1,1) = -.25*(1.+1./std::sqrt(3.)); dNxi(1,2,0) = +.25*(1.-1./std::sqrt(3.)); dNxi(1,2,1) = +.25*(1.+1./std::sqrt(3.)); dNxi(1,3,0) = -.25*(1.-1./std::sqrt(3.)); dNxi(1,3,1) = +.25*(1.-1./std::sqrt(3.)); // - k == 2 dNxi(2,0,0) = -.25*(1.-1./std::sqrt(3.)); dNxi(2,0,1) = -.25*(1.-1./std::sqrt(3.)); dNxi(2,1,0) = +.25*(1.-1./std::sqrt(3.)); dNxi(2,1,1) = -.25*(1.+1./std::sqrt(3.)); dNxi(2,2,0) = +.25*(1.+1./std::sqrt(3.)); dNxi(2,2,1) = +.25*(1.+1./std::sqrt(3.)); dNxi(2,3,0) = -.25*(1.+1./std::sqrt(3.)); dNxi(2,3,1) = +.25*(1.-1./std::sqrt(3.)); // - k == 3 dNxi(3,0,0) = -.25*(1.-1./std::sqrt(3.)); dNxi(3,0,1) = -.25*(1.+1./std::sqrt(3.)); dNxi(3,1,0) = +.25*(1.-1./std::sqrt(3.)); dNxi(3,1,1) = -.25*(1.-1./std::sqrt(3.)); dNxi(3,2,0) = +.25*(1.+1./std::sqrt(3.)); dNxi(3,2,1) = +.25*(1.-1./std::sqrt(3.)); dNxi(3,3,0) = -.25*(1.+1./std::sqrt(3.)); dNxi(3,3,1) = +.25*(1.+1./std::sqrt(3.)); // integration point weight at all Gauss points w(0) = 1.; w(1) = 1.; w(2) = 1.; w(3) = 1.; // shape function gradient at all nodes, in local coordinates // - k == 0 dNxi_n(0,0,0) = -0.5; dNxi_n(0,0,1) = -0.5; dNxi_n(0,1,0) = +0.5; dNxi_n(0,1,1) = 0.0; dNxi_n(0,2,0) = 0.0; dNxi_n(0,2,1) = 0.0; dNxi_n(0,3,0) = 0.0; dNxi_n(0,3,1) = +0.5; // - k == 1 dNxi_n(1,0,0) = -0.5; dNxi_n(1,0,1) = 0.0; dNxi_n(1,1,0) = +0.5; dNxi_n(1,1,1) = -0.5; dNxi_n(1,2,0) = 0.0; dNxi_n(1,2,1) = +0.5; dNxi_n(1,3,0) = 0.0; dNxi_n(1,3,1) = 0.0; // - k == 2 dNxi_n(2,0,0) = 0.0; dNxi_n(2,0,1) = 0.0; dNxi_n(2,1,0) = 0.0; dNxi_n(2,1,1) = -0.5; dNxi_n(2,2,0) = +0.5; dNxi_n(2,2,1) = +0.5; dNxi_n(2,3,0) = -0.5; dNxi_n(2,3,1) = 0.0; // - k == 3 dNxi_n(3,0,0) = 0.0; dNxi_n(3,0,1) = -0.5; dNxi_n(3,1,0) = 0.0; dNxi_n(3,1,1) = 0.0; dNxi_n(3,2,0) = +0.5; dNxi_n(3,2,1) = 0.0; dNxi_n(3,3,0) = -0.5; dNxi_n(3,3,1) = +0.5; // integration point weight at all nodes w_n(0) = 1.; w_n(1) = 1.; w_n(2) = 1.; w_n(3) = 1.; // Note, the above is a specialization of the following: // // - Shape function gradients // // dNxi(0,0) = -.25*(1.-xi(k,1)); dNxi(0,1) = -.25*(1.-xi(k,0)); // dNxi(1,0) = +.25*(1.-xi(k,1)); dNxi(1,1) = -.25*(1.+xi(k,0)); // dNxi(2,0) = +.25*(1.+xi(k,1)); dNxi(2,1) = +.25*(1.+xi(k,0)); // dNxi(3,0) = -.25*(1.+xi(k,1)); dNxi(3,1) = +.25*(1.-xi(k,0)); // // - Gauss point coordinates and weights // // xi(0,0) = -1./std::sqrt(3.); xi(0,1) = -1./std::sqrt(3.); w(0) = 1.; // xi(1,0) = +1./std::sqrt(3.); xi(1,1) = -1./std::sqrt(3.); w(1) = 1.; // xi(2,0) = +1./std::sqrt(3.); xi(2,1) = +1./std::sqrt(3.); w(2) = 1.; // xi(3,0) = -1./std::sqrt(3.); xi(3,1) = +1./std::sqrt(3.); w(3) = 1.; // // - Nodal coordinates and weights // // xi(0,0) = -1.; xi(0,1) = -1.; w(0) = 1.; // xi(1,0) = +1.; xi(1,1) = -1.; w(1) = 1.; // xi(2,0) = +1.; xi(2,1) = +1.; w(2) = 1.; // xi(3,0) = -1.; xi(3,1) = +1.; w(3) = 1.; } // ================================================================================================= template -void Quad4::updated_x() +inline void Quad4::updated_x() { #pragma omp parallel { // intermediate quantities cppmat::cartesian2d::tensor2 J_, Jinv_; double Jdet_; // local views of the global arrays (speeds up indexing, and increases readability) cppmat::tiny::matrix2 M_, D_; cppmat::tiny::matrix2 dNxi_, dNx_, x_; cppmat::tiny::vector w_, V_; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) { // pointer to element positions x_.map(&x(e)); // nodal quadrature // ---------------- // pointer to element mass/damping matrix, nodal volume, and integration weight M_.map(&M (e)); D_.map(&D (e)); V_.map(&V_n(e)); w_.map(&w_n(0)); // loop over nodes, i.e. the integration points for ( size_t k = 0 ; k < nne ; ++k ) { // - pointer to the shape function gradients (local coordinates) dNxi_.map(&dNxi_n(k)); // - Jacobian // J(i,j) += dNxi(m,i) * xe(m,j) J_(0,0) = dNxi_(0,0)*x_(0,0) + dNxi_(1,0)*x_(1,0) + dNxi_(2,0)*x_(2,0) + dNxi_(3,0)*x_(3,0); J_(0,1) = dNxi_(0,0)*x_(0,1) + dNxi_(1,0)*x_(1,1) + dNxi_(2,0)*x_(2,1) + dNxi_(3,0)*x_(3,1); J_(1,0) = dNxi_(0,1)*x_(0,0) + dNxi_(1,1)*x_(1,0) + dNxi_(2,1)*x_(2,0) + dNxi_(3,1)*x_(3,0); J_(1,1) = dNxi_(0,1)*x_(0,1) + dNxi_(1,1)*x_(1,1) + dNxi_(2,1)*x_(2,1) + dNxi_(3,1)*x_(3,1); // - determinant of the Jacobian Jdet_ = J_.det(); // - integration point volume V_(k) = w_(k) * Jdet_; // - assemble element mass matrix // M(m+i,n+i) += N(m) * rho * V * N(n); M_(k*2 ,k*2 ) = mat->rho(e,k) * V_(k); M_(k*2+1,k*2+1) = mat->rho(e,k) * V_(k); // - assemble element non-Galilean damping matrix // D(m+i,n+i) += N(m) * alpha * V * N(n); D_(k*2 ,k*2 ) = mat->alpha(e,k) * V_(k); D_(k*2+1,k*2+1) = mat->alpha(e,k) * V_(k); } // Gaussian quadrature // ------------------- // pointer to element integration volume and weight V_.map(&V(e)); w_.map(&w(0)); // loop over Gauss points for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to the shape function gradients (local/global coordinates) dNxi_.map(&dNxi( k)); dNx_ .map(&dNx (e,k)); // - Jacobian // J(i,j) += dNxi(m,i) * xe(m,j) J_(0,0) = dNxi_(0,0)*x_(0,0) + dNxi_(1,0)*x_(1,0) + dNxi_(2,0)*x_(2,0) + dNxi_(3,0)*x_(3,0); J_(0,1) = dNxi_(0,0)*x_(0,1) + dNxi_(1,0)*x_(1,1) + dNxi_(2,0)*x_(2,1) + dNxi_(3,0)*x_(3,1); J_(1,0) = dNxi_(0,1)*x_(0,0) + dNxi_(1,1)*x_(1,0) + dNxi_(2,1)*x_(2,0) + dNxi_(3,1)*x_(3,0); J_(1,1) = dNxi_(0,1)*x_(0,1) + dNxi_(1,1)*x_(1,1) + dNxi_(2,1)*x_(2,1) + dNxi_(3,1)*x_(3,1); // - determinant and inverse of the Jacobian Jdet_ = J_.det(); Jinv_ = J_.inv(); // - integration point volume V_(k) = w_(k) * Jdet_; // - shape function gradients (global coordinates) // dNx(m,i) += Jinv(i,j) * dNxi(m,j) for ( size_t m = 0 ; m < nne ; ++m ) { dNx_(m,0) = Jinv_(0,0) * dNxi_(m,0) + Jinv_(0,1) * dNxi_(m,1); dNx_(m,1) = Jinv_(1,0) * dNxi_(m,0) + Jinv_(1,1) * dNxi_(m,1); } } } // set signals changed_f = false; changed_M = true; changed_D = true; } // #pragma omp parallel } // ================================================================================================= template -void Quad4::updated_u() +inline void Quad4::updated_u() { #pragma omp parallel { // intermediate quantities cppmat::cartesian2d::tensor2 gradu_; cppmat::cartesian2d::tensor2s eps_, sig_; // local views of the global arrays (speeds up indexing, and increases readability) cppmat::tiny::matrix2 dNx_, u_, f_; cppmat::tiny::vector V_; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) { // pointer to element forces, displacements, and integration volume f_.map(&f(e)); u_.map(&u(e)); V_.map(&V(e)); // zero initialize forces f_.setZero(); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to the shape function gradients, strain and stress tensor (stored symmetric) dNx_.map(&dNx (e,k)); eps_.map(&mat->eps(e,k)); sig_.map(&mat->sig(e,k)); // - displacement gradient // gradu_(i,j) += dNx(m,i) * ue(m,j) gradu_(0,0) = dNx_(0,0)*u_(0,0) + dNx_(1,0)*u_(1,0) + dNx_(2,0)*u_(2,0) + dNx_(3,0)*u_(3,0); gradu_(0,1) = dNx_(0,0)*u_(0,1) + dNx_(1,0)*u_(1,1) + dNx_(2,0)*u_(2,1) + dNx_(3,0)*u_(3,1); gradu_(1,0) = dNx_(0,1)*u_(0,0) + dNx_(1,1)*u_(1,0) + dNx_(2,1)*u_(2,0) + dNx_(3,1)*u_(3,0); gradu_(1,1) = dNx_(0,1)*u_(0,1) + dNx_(1,1)*u_(1,1) + dNx_(2,1)*u_(2,1) + dNx_(3,1)*u_(3,1); // - strain (stored symmetric) // eps(i,j) = .5 * ( gradu_(i,j) + gradu_(j,i) ) eps_(0,0) = gradu_(0,0); eps_(0,1) = .5 * ( gradu_(0,1) + gradu_(1,0) ); eps_(1,1) = gradu_(1,1); // - constitutive response mat->updated_eps(e,k); // - assemble to element force // f(m,j) += dNx(m,i) * sig(i,j) * V; for ( size_t m = 0 ; m < nne ; ++m ) { f_(m,0) += dNx_(m,0) * sig_(0,0) * V_(k) + dNx_(m,1) * sig_(1,0) * V_(k); f_(m,1) += dNx_(m,0) * sig_(0,1) * V_(k) + dNx_(m,1) * sig_(1,1) * V_(k); } } } // set signals changed_f = true; changed_M = false; changed_D = false; } } // ================================================================================================= template -void Quad4::updated_v() +inline void Quad4::updated_v() { #pragma omp parallel { // intermediate quantities cppmat::cartesian2d::tensor2 gradv_; cppmat::cartesian2d::tensor2s epsdot_, sig_; // local views of the global arrays (speeds up indexing, and increases readability) cppmat::tiny::matrix2 dNx_, v_, f_; cppmat::tiny::vector V_; // loop over all elements (in parallel) #pragma omp for for ( size_t e = 0 ; e < nelem ; ++e ) { // pointer to element forces, displacements, and integration volume f_.map(&f(e)); v_.map(&v(e)); V_.map(&V(e)); // zero initialize forces f_.setZero(); // loop over all integration points in element "e" for ( size_t k = 0 ; k < nip ; ++k ) { // - pointer to the shape function gradients, strain-rate and stress tensor (stored symmetric) dNx_ .map(&dNx (e,k)); epsdot_.map(&mat->epsdot(e,k)); sig_ .map(&mat->sig (e,k)); // - displacement gradient // gradv_(i,j) += dNx(m,i) * ue(m,j) gradv_(0,0) = dNx_(0,0)*v_(0,0) + dNx_(1,0)*v_(1,0) + dNx_(2,0)*v_(2,0) + dNx_(3,0)*v_(3,0); gradv_(0,1) = dNx_(0,0)*v_(0,1) + dNx_(1,0)*v_(1,1) + dNx_(2,0)*v_(2,1) + dNx_(3,0)*v_(3,1); gradv_(1,0) = dNx_(0,1)*v_(0,0) + dNx_(1,1)*v_(1,0) + dNx_(2,1)*v_(2,0) + dNx_(3,1)*v_(3,0); gradv_(1,1) = dNx_(0,1)*v_(0,1) + dNx_(1,1)*v_(1,1) + dNx_(2,1)*v_(2,1) + dNx_(3,1)*v_(3,1); // - strain (stored symmetric) // epsdot(i,j) = .5 * ( gradv_(i,j) + gradv_(j,i) ) epsdot_(0,0) = gradv_(0,0); epsdot_(0,1) = .5 * ( gradv_(0,1) + gradv_(1,0) ); epsdot_(1,1) = gradv_(1,1); // - constitutive response mat->updated_epsdot(e,k); // - assemble to element force // f(m,j) += dNdx(m,i) * sig(i,j) * V; for ( size_t m = 0 ; m < nne ; ++m ) { f_(m,0) += dNx_(m,0) * sig_(0,0) * V_(k) + dNx_(m,1) * sig_(1,0) * V_(k); f_(m,1) += dNx_(m,0) * sig_(0,1) * V_(k) + dNx_(m,1) * sig_(1,1) * V_(k); } } } // set signals changed_f = true; changed_M = false; changed_D = false; } } // ================================================================================================= }}}} // namespace ... #endif diff --git a/src/GooseFEM/Mesh.h b/src/GooseFEM/Mesh.h index cb9c53d..67903db 100644 --- a/src/GooseFEM/Mesh.h +++ b/src/GooseFEM/Mesh.h @@ -1,160 +1,160 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MESH_H #define GOOSEFEM_MESH_H #include "Macros.h" // ------------------------------------------------------------------------------------------------- namespace GooseFEM { namespace Mesh { // ====================================== EXTRACT INFORMATION ====================================== // elements connected to each node: // out[ : , 0 ] = number of elements connected to each node // out[ j , i+1 ] = "i"th element connected to node "j" -MatS elem2node ( const MatS &conn ); +inline MatS elem2node ( const MatS &conn ); // list with DOF-numbers in sequential order -MatS dofs ( size_t nnode , size_t ndim ); +inline MatS dofs ( size_t nnode , size_t ndim ); // renumber list of indices // - renumber to lowest possible numbers (e.g. [0,3,4,2] -> [0,2,3,1]) -MatS renumber ( const MatS &in ); +inline MatS renumber ( const MatS &in ); // - renumber to begin [ idx , ... ] or end [ ... , idx ] (e.g. [0,1,3,2] -> [3,0,2,1]; with idx=0) -MatS renumber ( const MatS &in , const ColS &idx , std::string location="end" ); +inline MatS renumber ( const MatS &in , const ColS &idx , std::string location="end" ); // ========================================== SOURCE CODE ========================================== -MatS elem2node ( const MatS &conn ) +inline MatS elem2node ( const MatS &conn ) { size_t nnode = conn.maxCoeff() + 1; ColS N( nnode ); N *= 0; for ( auto e = 0 ; e < conn.rows() ; ++e ) for ( auto m = 0 ; m < conn.cols() ; ++m ) N( conn(e,m) ) += 1; MatS out( nnode , N.maxCoeff()+1 ); out *= 0; for ( auto e = 0 ; e < conn.rows() ; ++e ) { for ( auto m = 0 ; m < conn.cols() ; ++m ) { auto nd = conn(e,m); out( nd, 0 ) += 1; out( nd, out(nd,0) ) = e; } } return out; } // ------------------------------------------------------------------------------------------------- -MatS dofs ( size_t nnode , size_t ndim ) +inline MatS dofs ( size_t nnode , size_t ndim ) { ColS dofs_vec = ColS::LinSpaced( nnode*ndim , 0 , nnode*ndim ); Eigen::Map dofs( dofs_vec.data() , nnode , ndim ); return dofs; } // ------------------------------------------------------------------------------------------------- -MatS renumber ( const MatS &in ) +inline MatS renumber ( const MatS &in ) { // size parameters size_t N = in.size(); // number of items in "in" (and in "out") size_t M = in.maxCoeff() + 1; // size of renumber-list // allocate output MatS out ( in.rows() , in.cols() ); // allocate renumber-list ColS renum ( M ); // set all items as being 'excluded' ( == 0 ) renum.setZero(); // selectively set items that are 'included' ( == 1 ) for ( size_t i = 0 ; i < N ; ++i ) renum ( in(i) ) = 1; // convert to indices, step 1: correct such that first index will start at 0 renum[0] -= 1; // convert to indices, step 2: cumulative sum for ( size_t i = 1 ; i < M ; ++i ) renum ( i ) += renum ( i-1 ); // renumber DOFs for ( size_t i = 0 ; i < N ; ++i ) out ( i ) = renum ( in(i) ); return out; } // ------------------------------------------------------------------------------------------------- -MatS renumber ( const MatS &in , const ColS &idx , std::string loc ) +inline MatS renumber ( const MatS &in , const ColS &idx , std::string loc ) { // copy input to make modifications ColS tmp = idx; MatS out = in; ColS dat = ColS::LinSpaced(in.maxCoeff()+1, 0, in.maxCoeff()+1); // store size size_t ndx = static_cast( idx.size() ); size_t N = static_cast( out.size() ); // sort input std::sort ( tmp.data(), tmp.data()+tmp.size() ); // find "jdx = setdiff ( dat , idx )" // - allocate ColS jdx ( dat.size() - ndx ); // - compute std::set_difference ( dat.data() , dat.data()+dat.size() , tmp.data() , tmp.data()+tmp.size() , jdx.data() ); // store size size_t mdx = static_cast( jdx.size() ); // renumber-list // - allocate ColS renum ( dat.size() ); // - fill if ( loc == "end" ) { // -- order [ jdx , idx ] for ( size_t i = 0 ; i < mdx ; ++i ) renum ( jdx(i) ) = i; for ( size_t i = 0 ; i < ndx ; ++i ) renum ( idx(i) ) = i + mdx; } else if ( loc == "begin" ) { // -- order [ jdx, idx ] for ( size_t i = 0 ; i < ndx ; ++i ) renum ( idx(i) ) = i; for ( size_t i = 0 ; i < mdx ; ++i ) renum ( jdx(i) ) = i + ndx; } // renumber for ( size_t i = 0 ; i < N ; ++i ) out ( i ) = renum ( out (i) ); return out; } // ------------------------------------------------------------------------------------------------- } // namespace Mesh } // namespace GooseFEM #endif diff --git a/src/GooseFEM/MeshQuad4.h b/src/GooseFEM/MeshQuad4.h index 6cf7c54..f779dde 100644 --- a/src/GooseFEM/MeshQuad4.h +++ b/src/GooseFEM/MeshQuad4.h @@ -1,764 +1,764 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MESHQUAD4_H #define GOOSEFEM_MESHQUAD4_H #include "Macros.h" #include "Mesh.h" // ------------------------------------------------------------------------------------------------- namespace GooseFEM { namespace Mesh { namespace Quad4 { // ========================================== REGULAR MESH ========================================= class Regular { private: size_t m_nx; // number of 'pixels' horizontal direction (length == "m_nx * m_h") size_t m_ny; // number of 'pixels' vertical direction (length == "m_ny * m_h") double m_h; // size of the element edge (equal in both directions) size_t m_nelem; // number of elements size_t m_nnode; // number of nodes size_t m_nne=4; // number of nodes-per-element size_t m_ndim=2; // number of dimensions public: // mesh with "nx" pixels in horizontal direction, "ny" in vertical direction and "h" the edge size Regular(size_t nx, size_t ny, double h=1.); size_t nelem() { return m_nelem;}; // number of elements size_t nnode() { return m_nnode;}; // number of nodes size_t nne () { return m_nne; }; // number of nodes-per-element size_t ndim () { return m_ndim; }; // number of dimensions MatD coor (); // nodal positions [ nnode , ndim ] MatS conn (); // connectivity [ nelem , nne ] ColS nodesBottom (); // nodes along the bottom edge ColS nodesTop (); // nodes along the top edge ColS nodesLeft (); // nodes along the left edge ColS nodesRight (); // nodes along the right edge MatS nodesPeriodic(); // periodic node pairs [ : , 2 ]: ( independent , dependent ) size_t nodeOrigin (); // bottom-left node, to be used as reference for periodicity MatS dofs (); // DOF-numbers for each component of each node (sequential) MatS dofsPeriodic (); // DOF-numbers for each component of each node (sequential) }; // ====================== MESH WITH A FINE LAYER THAT EXPONENTIALLY COARSENS ======================= class FineLayer { private: double m_h; // base size of the element edge (equal in both directions) size_t m_nx; // number of elements in vertical direction ColS m_nh; // element size in vertical direction (number of time "h") ColS m_startNode; // start node of each row ColS m_startElem; // start element of each row size_t m_nelem; // number of elements size_t m_nnode; // number of nodes size_t m_nne=4; // number of nodes-per-element size_t m_ndim=2; // number of dimensions public: // mesh with "nx" pixels in horizontal direction, "ny" in vertical direction and "h" the edge size FineLayer(size_t nx, size_t ny, double h=1., size_t nfine=0, size_t nskip=0); size_t nelem() { return m_nelem;}; // number of elements size_t nnode() { return m_nnode;}; // number of nodes size_t nne () { return m_nne; }; // number of nodes-per-element size_t ndim () { return m_ndim; }; // number of dimensions size_t shape(size_t i); // actual shape in horizontal and vertical direction MatD coor (); // nodal positions [ nnode , ndim ] MatS conn (); // connectivity [ nelem , nne ] ColS elementsFine (); // elements in the middle, fine, layer ColS nodesBottom (); // nodes along the bottom edge ColS nodesTop (); // nodes along the top edge ColS nodesLeft (); // nodes along the left edge ColS nodesRight (); // nodes along the right edge MatS nodesPeriodic(); // periodic node pairs [ : , 2 ]: ( independent , dependent ) size_t nodeOrigin (); // bottom-left node, to be used as reference for periodicity MatS dofs (); // DOF-numbers for each component of each node (sequential) MatS dofsPeriodic (); // DOF-numbers for each component of each node (sequential) }; // ===================================== SOURCE CODE : REGULAR ===================================== -Regular::Regular(size_t nx, size_t ny, double h): m_nx(nx), m_ny(ny), m_h(h) +inline Regular::Regular(size_t nx, size_t ny, double h): m_nx(nx), m_ny(ny), m_h(h) { assert( m_nx >= 1 ); assert( m_ny >= 1 ); m_nnode = (m_nx+1) * (m_ny+1); m_nelem = m_nx * m_ny ; } // ------------------------------------------------------------------------------------------------- -MatD Regular::coor() +inline MatD Regular::coor() { MatD coor( m_nnode , m_ndim ); ColD x = ColD::LinSpaced( m_nx+1 , 0.0 , m_h * static_cast(m_nx) ); ColD y = ColD::LinSpaced( m_ny+1 , 0.0 , m_h * static_cast(m_ny) ); size_t inode = 0; for ( size_t row = 0 ; row < m_ny+1 ; ++row ) { for ( size_t col = 0 ; col < m_nx+1 ; ++col ) { coor(inode,0) = x(col); coor(inode,1) = y(row); ++inode; } } return coor; } // ------------------------------------------------------------------------------------------------- -MatS Regular::conn() +inline MatS Regular::conn() { MatS conn( m_nelem , m_nne ); size_t ielem = 0; for ( size_t row = 0 ; row < m_ny ; ++row ) { for ( size_t col = 0 ; col < m_nx ; ++col ) { conn(ielem,0) = (row+0)*(m_nx+1)+col+0; conn(ielem,1) = (row+0)*(m_nx+1)+col+1; conn(ielem,3) = (row+1)*(m_nx+1)+col+0; conn(ielem,2) = (row+1)*(m_nx+1)+col+1; ++ielem; } } return conn; } // ------------------------------------------------------------------------------------------------- -ColS Regular::nodesBottom() +inline ColS Regular::nodesBottom() { ColS nodes(m_nx+1); for ( size_t col = 0 ; col < m_nx+1 ; ++col ) nodes(col) = col; return nodes; } // ------------------------------------------------------------------------------------------------- -ColS Regular::nodesTop() +inline ColS Regular::nodesTop() { ColS nodes(m_nx+1); for ( size_t col = 0 ; col < m_nx+1 ; ++col ) nodes(col) = col+m_ny*(m_nx+1); return nodes; } // ------------------------------------------------------------------------------------------------- -ColS Regular::nodesLeft() +inline ColS Regular::nodesLeft() { ColS nodes(m_ny+1); for ( size_t row = 0 ; row < m_ny+1 ; ++row ) nodes(row) = row*(m_nx+1); return nodes; } // ------------------------------------------------------------------------------------------------- -ColS Regular::nodesRight() +inline ColS Regular::nodesRight() { ColS nodes(m_ny+1); for ( size_t row = 0 ; row < m_ny+1 ; ++row ) nodes(row) = row*(m_nx+1)+m_nx; return nodes; } // ------------------------------------------------------------------------------------------------- -MatS Regular::nodesPeriodic() +inline MatS Regular::nodesPeriodic() { ColS bot = nodesBottom(); ColS top = nodesTop (); ColS rgt = nodesRight (); ColS lft = nodesLeft (); MatS nodes( bot.size()-2 + lft.size()-2 + 3 , 2 ); size_t i = 0; nodes(i,0) = bot(0); nodes(i,1) = bot(bot.size()-1); ++i; nodes(i,0) = bot(0); nodes(i,1) = top(top.size()-1); ++i; nodes(i,0) = bot(0); nodes(i,1) = top(0 ); ++i; for ( auto j = 1 ; j < bot.size()-1 ; ++j ) { nodes(i,0) = bot(j); nodes(i,1) = top(j); ++i; } for ( auto j = 1 ; j < lft.size()-1 ; ++j ) { nodes(i,0) = lft(j); nodes(i,1) = rgt(j); ++i; } return nodes; } // ------------------------------------------------------------------------------------------------- -size_t Regular::nodeOrigin() +inline size_t Regular::nodeOrigin() { return 0; } // ------------------------------------------------------------------------------------------------- -MatS Regular::dofs() +inline MatS Regular::dofs() { return GooseFEM::Mesh::dofs(m_nnode,m_ndim); } // ------------------------------------------------------------------------------------------------- -MatS Regular::dofsPeriodic() +inline MatS Regular::dofsPeriodic() { // DOF-numbers for each component of each node (sequential) MatS out = GooseFEM::Mesh::dofs(m_nnode,m_ndim); // periodic node-pairs MatS nodePer = nodesPeriodic(); size_t nper = static_cast(nodePer.rows()); // eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs for ( size_t i = 0 ; i < nper ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) out(nodePer(i,1),j) = out(nodePer(i,0),j); // renumber "out" to be sequential return GooseFEM::Mesh::renumber(out); } // ==================================== SOURCE CODE : FINELAYER ==================================== -FineLayer::FineLayer(size_t nx, size_t ny, double h, size_t nfine, size_t nskip) +inline FineLayer::FineLayer(size_t nx, size_t ny, double h, size_t nfine, size_t nskip) { assert( nx >= 1 ); assert( ny >= 1 ); m_nx = nx; m_h = h; // compute the element size : based on symmetric half // -------------------------------------------------- // convert height to half of the height if ( ny % 2 == 0 ) ny = ny /2; else ny = (ny+1)/2; // check the number of fine layers from the center assert( nfine <= ny ); // define arrays to determine to coarsen ColS n (ny+1); // size of the element (number of times "h") ColS nsum(ny+1); // current size of the mesh in y-direction (in number of times "nsum") // initialize all elements of size "1" size_t next = 1; n.setOnes(); // size of the mesh in y-direction after "i" element // - initialize nsum(0) = 1; // - compute based on the initially equi-sized elements for ( size_t i = 1 ; i < ny+1 ; ++i ) nsum(i) = nsum(i-1) + n(i-1); // loop in vertical direction and check to coarsen; rules: // - the size of the element cannot be greater than the distance // - the size of the element should match the horizontal direction // - skip a certain number of elements size_t min = 1; if ( nfine > min ) min = nfine; for ( size_t i = min ; i < ny+1 ; ++i ) { if ( 3*n(i-1) <= nsum(i-1) and m_nx % ( 3*next ) == 0 ) { n (i) = n (i-1) * 2; nsum(i) = nsum(i-1) + n(i); next *= 3; } else { n (i) = next; nsum(i) = nsum(i-1) + n(i); } } // skip the last "nskip" coarsening steps // - counter : number of steps skipped size_t iskip = 0; // - loop to skip for ( size_t i = ny+1 ; i-- > 1 ; ) { // -- check -> quit if the number of steps is reached if ( iskip >= nskip ) break; // -- if coarsening step found -> remove if ( n(i) % 2 == 0 ) { for ( size_t j = i ; j < ny+1 ; ++j ) n(j) = n(i-1); iskip++; } } // - update the height after "i" elements for ( size_t i = 1 ; i < ny+1 ; ++i ) nsum(i) = nsum(i-1) + n(i); // truncate such that the height does not exceed "ny" size_t N; for ( N = 0 ; N < ny+1 ; ++N ) if ( nsum(N) > ny+1 ) break; // truncate n.conservativeResize(N); // compute mesh dimensions : based on the entire mesh // -------------------------------------------------- // symmetrize // - get size N = static_cast(n.size()); // - allocate m_nh.conservativeResize(2*N-1); // - store symmetrized for ( size_t i = 0 ; i < N ; i++ ) m_nh( i ) = n(N-1-i); for ( size_t i = 1 ; i < N ; i++ ) m_nh(N+i-1) = n( i); // allocate counters m_startElem.conservativeResize(m_nh.size() ); m_startNode.conservativeResize(m_nh.size()+1); // compute the dimensions of the mesh // - initialize m_nnode = 0; m_nelem = 0; N = static_cast(m_nh.size()); m_startNode(0) = 0; // loop from bottom to middle : elements become finer for ( size_t i = 0 ; i < (N-1)/2 ; ++i ) { // - store the first element of the row m_startElem(i) = m_nelem; // - add the nodes of this row if ( m_nh(i)%2 == 0 ) { m_nnode += 3 * m_nx/(3*m_nh(i)/2) + 1; } else { m_nnode += m_nx/ m_nh(i) + 1; } // - add the elements of this row if ( m_nh(i)%2 == 0 ) { m_nelem += 4 * m_nx/(3*m_nh(i)/2); } else { m_nelem += m_nx/ m_nh(i) ; } // - store the starting node of the next row m_startNode(i+1) = m_nnode; } // loop from middle to top : elements become coarser for ( size_t i = (N-1)/2 ; i < N ; ++i ) { // - store the first element of the row m_startElem(i) = m_nelem; // - add the nodes of this row if ( m_nh(i)%2 == 0 ) { m_nnode += 1 + m_nx/(m_nh(i)/2) + 2 * m_nx/(3*m_nh(i)/2); } else { m_nnode += 1 + m_nx/ m_nh(i); } // - add the elements of this row if ( m_nh(i)%2 == 0 ) { m_nelem += 4 * m_nx/(3*m_nh(i)/2); } else { m_nelem += m_nx/ m_nh(i) ; } // - store the starting node of the next row m_startNode(i+1) = m_nnode; } // add the top row of nodes to the number of rows (starting node already stored) if ( m_nh(N-1)%2 == 0 ) { m_nnode += 1 + m_nx/(3*m_nh(N-1)/2); } else { m_nnode += 1 + m_nx/ m_nh(N-1); } } // ------------------------------------------------------------------------------------------------- -size_t FineLayer::shape(size_t i) +inline size_t FineLayer::shape(size_t i) { // check the index assert( i >= 0 and i < 2 ); // horizontal direction if ( i == 0 ) return m_nx; // vertical direction // - zero-initialize size_t n = 0; // - cumulative sum of the size per row for ( size_t i = 0 ; i < static_cast(m_nh.size()) ; ++i ) n += m_nh(i); return n; } // ------------------------------------------------------------------------------------------------- -MatD FineLayer::coor() +inline MatD FineLayer::coor() { // allocate output MatD out( m_nnode , m_ndim ); // initialize position in horizontal direction (of the finest nodes) ColD x = ColD::LinSpaced( m_nx+1 , 0.0 , m_h * static_cast(m_nx) ); // zero-initialize current node and height: loop from top to bottom; number of rows double h = 0; size_t inode = 0; size_t N = static_cast(m_nh.size()); // lower half for ( size_t i = 0 ; i < (N-1)/2 ; ++i ) { // - refinement row if ( m_nh(i)%2 == 0 ) { // -- bottom row: coarse nodes for ( size_t j = 0 ; j < m_nx+1 ; j += 3*m_nh(i)/2 ) { out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; } h += static_cast(m_nh(i)/2); // -- middle row: part the fine nodes for ( size_t j = 0 ; j < m_nx ; j += 3*m_nh(i)/2 ) { out(inode,0) = x(j+m_nh(i)/2); out(inode,1) = h*m_h; inode++; out(inode,0) = x(j+m_nh(i) ); out(inode,1) = h*m_h; inode++; } h += static_cast(m_nh(i)/2); } // - normal row else { for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(i) ) { out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; } h += static_cast(m_nh(i)); } } // top half for ( size_t i = (N-1)/2 ; i < N ; ++i ) { // - coarsening row if ( m_nh(i)%2 == 0 ) { // -- bottom row: fine nodes for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(i)/2 ) { out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; } h += static_cast(m_nh(i)/2); // -- middle row: part the fine nodes for ( size_t j = 0 ; j < m_nx ; j += 3*m_nh(i)/2 ) { out(inode,0) = x(j+m_nh(i)/2); out(inode,1) = h*m_h; inode++; out(inode,0) = x(j+m_nh(i) ); out(inode,1) = h*m_h; inode++; } h += static_cast(m_nh(i)/2); } // - normal row else { for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(i) ) { out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; } h += static_cast(m_nh(i)); } } // top row if ( m_nh(N-1)%2 == 0 ) { for ( size_t j = 0 ; j < m_nx+1 ; j += 3*m_nh(N-1)/2 ) { out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; } } else { for ( size_t j = 0 ; j < m_nx+1 ; j += m_nh(N-1) ) { out(inode,0) = x(j); out(inode,1) = h*m_h; inode++; } } return out; } // // ------------------------------------------------------------------------------------------------- -MatS FineLayer::conn() +inline MatS FineLayer::conn() { // allocate output MatS out( m_nelem , m_nne ); // current element, number of rows, beginning nodes size_t ielem = 0; size_t N = static_cast(m_nh.size()); size_t bot,mid,top; // bottom half for ( size_t i = 0 ; i < (N-1)/2 ; ++i ) { // - starting nodes of bottom and top layer (and middle layer, only relevant for even shape) if ( true ) bot = m_startNode(i ); if ( m_nh(i)%2 == 0 ) mid = m_startNode(i )+m_nx/(3*m_nh(i)/2)+1; if ( true ) top = m_startNode(i+1); // - even sized shape: refinement layer if ( m_nh(i)%2 == 0 ) { for ( size_t j = 0 ; j < m_nx/(3*m_nh(i)/2) ; ++j ) { // -- lower out(ielem,0) = j +bot; out(ielem,1) = j +1+bot; out(ielem,2) = j*2+1+mid; out(ielem,3) = j*2 +mid; ielem++; // -- upper right out(ielem,0) = j +1+bot; out(ielem,1) = j*3+3+top; out(ielem,2) = j*3+2+top; out(ielem,3) = j*2+1+mid; ielem++; // -- upper middle out(ielem,0) = j*2 +mid; out(ielem,1) = j*2+1+mid; out(ielem,2) = j*3+2+top; out(ielem,3) = j*3+1+top; ielem++; // -- upper left out(ielem,0) = j +bot; out(ielem,1) = j*2 +mid; out(ielem,2) = j*3+1+top; out(ielem,3) = j*3 +top; ielem++; } } // - odd sized shape: normal layer else { for ( size_t j = 0 ; j < m_nx/m_nh(i) ; ++j ) { out(ielem,0) = j +bot; out(ielem,1) = j+1+bot; out(ielem,2) = j+1+top; out(ielem,3) = j +top; ielem++; } } } // top half for ( size_t i = (N-1)/2 ; i < N ; ++i ) { // - starting nodes of bottom and top layer (and middle layer, only relevant for even shape) if ( true ) bot = m_startNode(i ); if ( m_nh(i)%2 == 0 ) mid = m_startNode(i )+m_nx/(m_nh(i)/2)+1; if ( true ) top = m_startNode(i+1); // - even sized shape: refinement layer if ( m_nh(i)%2 == 0 ) { for ( size_t j = 0 ; j < m_nx/(3*m_nh(i)/2) ; ++j ) { // -- lower left out(ielem,0) = j*3 +bot; out(ielem,1) = j*3+1+bot; out(ielem,2) = j*2 +mid; out(ielem,3) = j +top; ielem++; // -- lower middle out(ielem,0) = j*3+1+bot; out(ielem,1) = j*3+2+bot; out(ielem,2) = j*2+1+mid; out(ielem,3) = j*2 +mid; ielem++; // -- lower right out(ielem,0) = j*3+2+bot; out(ielem,1) = j*3+3+bot; out(ielem,2) = j +1+top; out(ielem,3) = j*2+1+mid; ielem++; // -- upper out(ielem,0) = j +top; out(ielem,1) = j*2 +mid; out(ielem,2) = j*2+1+mid; out(ielem,3) = j +1+top; ielem++; } } // - odd sized shape: normal layer else { for ( size_t j = 0 ; j < m_nx/m_nh(i) ; ++j ) { out(ielem,0) = j +bot; out(ielem,1) = j+1+bot; out(ielem,2) = j+1+top; out(ielem,3) = j +top; ielem++; } } } return out; } // ------------------------------------------------------------------------------------------------- -ColS FineLayer::elementsFine() +inline ColS FineLayer::elementsFine() { ColS out(m_nx); size_t j = m_startElem((m_startElem.size()-1)/2); for ( size_t i = 0 ; i < m_nx ; ++i ) out(i) = i+j; return out; } // ------------------------------------------------------------------------------------------------- -ColS FineLayer::nodesBottom() +inline ColS FineLayer::nodesBottom() { ColS nodes; if ( m_nh(0)%2 == 0 ) nodes.conservativeResize(m_nx/(3*m_nh(0)/2)+1); else nodes.conservativeResize(m_nx/ m_nh(0) +1); for ( size_t i = 0 ; i < static_cast(nodes.size()) ; ++i ) nodes(i) = i; return nodes; } // ------------------------------------------------------------------------------------------------- -ColS FineLayer::nodesTop() +inline ColS FineLayer::nodesTop() { ColS nodes; if ( m_nh(0)%2 == 0 ) nodes.conservativeResize(m_nx/(3*m_nh(0)/2)+1); else nodes.conservativeResize(m_nx/ m_nh(0) +1); size_t j = m_startNode(m_startNode.size()-1); for ( size_t i = 0 ; i < static_cast(nodes.size()) ; ++i ) nodes(i) = i+j; return nodes; } // ------------------------------------------------------------------------------------------------- -ColS FineLayer::nodesLeft() +inline ColS FineLayer::nodesLeft() { return m_startNode; } // ------------------------------------------------------------------------------------------------- -ColS FineLayer::nodesRight() +inline ColS FineLayer::nodesRight() { size_t N = static_cast(m_nh.size()); ColS nodes(N+1); // bottom half for ( size_t i = 0 ; i < (N-1)/2 ; ++i ) { if ( m_nh(i)%2 == 0 ) nodes(i) = m_startNode(i) + m_nx/(3*m_nh(i)/2); else nodes(i) = m_startNode(i) + m_nx/ m_nh(i) ; } // top half for ( size_t i = (N-1)/2 ; i < N ; ++i ) { if ( m_nh(i)%2 == 0 ) nodes(i) = m_startNode(i) + m_nx/(m_nh(i)/2); else nodes(i) = m_startNode(i) + m_nx/ m_nh(i) ; } // last node nodes(N) = m_nnode-1; return nodes; } // ------------------------------------------------------------------------------------------------- -MatS FineLayer::nodesPeriodic() +inline MatS FineLayer::nodesPeriodic() { ColS bot = nodesBottom(); ColS top = nodesTop (); ColS rgt = nodesRight (); ColS lft = nodesLeft (); MatS nodes( bot.size()-2 + lft.size()-2 + 3 , 2 ); size_t i = 0; nodes(i,0) = bot(0); nodes(i,1) = bot(bot.size()-1); ++i; nodes(i,0) = bot(0); nodes(i,1) = top(top.size()-1); ++i; nodes(i,0) = bot(0); nodes(i,1) = top(0 ); ++i; for ( auto j = 1 ; j < bot.size()-1 ; ++j ) { nodes(i,0) = bot(j); nodes(i,1) = top(j); ++i; } for ( auto j = 1 ; j < lft.size()-1 ; ++j ) { nodes(i,0) = lft(j); nodes(i,1) = rgt(j); ++i; } return nodes; } // ------------------------------------------------------------------------------------------------- -size_t FineLayer::nodeOrigin() +inline size_t FineLayer::nodeOrigin() { return 0; } // ------------------------------------------------------------------------------------------------- -MatS FineLayer::dofs() +inline MatS FineLayer::dofs() { return GooseFEM::Mesh::dofs(m_nnode,m_ndim); } // ------------------------------------------------------------------------------------------------- -MatS FineLayer::dofsPeriodic() +inline MatS FineLayer::dofsPeriodic() { // DOF-numbers for each component of each node (sequential) MatS out = GooseFEM::Mesh::dofs(m_nnode,m_ndim); // periodic node-pairs MatS nodePer = nodesPeriodic(); size_t nper = static_cast(nodePer.rows()); // eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs for ( size_t i = 0 ; i < nper ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) out(nodePer(i,1),j) = out(nodePer(i,0),j); // renumber "out" to be sequential return GooseFEM::Mesh::renumber(out); } // ================================================================================================= }}} // namespace GooseFEM::Mesh::Quad4 #endif diff --git a/src/GooseFEM/MeshTri3.h b/src/GooseFEM/MeshTri3.h index b8cf1ab..7b7ae4d 100644 --- a/src/GooseFEM/MeshTri3.h +++ b/src/GooseFEM/MeshTri3.h @@ -1,565 +1,565 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MESHTRI3_H #define GOOSEFEM_MESHTRI3_H #include "Macros.h" #include "Mesh.h" // ------------------------------------------------------------------------------------------------- namespace GooseFEM { namespace Mesh { namespace Tri3 { // ========================================== REGULAR MESH ========================================= class Regular { private: size_t m_nx; // number of 'pixels' horizontal direction (length == "m_nx * m_h") size_t m_ny; // number of 'pixels' vertical direction (length == "m_ny * m_h") double m_h; // size of the element edge (equal in both directions) size_t m_nelem; // number of elements size_t m_nnode; // number of nodes size_t m_nne=3; // number of nodes-per-element size_t m_ndim=2; // number of dimensions public: // mesh with "nx" pixels in horizontal direction, "ny" in vertical direction and "h" the edge size Regular(size_t nx, size_t ny, double h=1.); size_t nelem() { return m_nelem;}; // number of elements size_t nnode() { return m_nnode;}; // number of nodes size_t nne () { return m_nne; }; // number of nodes-per-element size_t ndim () { return m_ndim; }; // number of dimensions MatD coor (); // nodal positions [ nnode , ndim ] MatS conn (); // connectivity [ nelem , nne ] ColS nodesBottom (); // nodes along the bottom edge ColS nodesTop (); // nodes along the top edge ColS nodesLeft (); // nodes along the left edge ColS nodesRight (); // nodes along the right edge MatS nodesPeriodic(); // periodic node pairs [ : , 2 ]: ( independent , dependent ) size_t nodeOrigin (); // bottom-left node, to be used as reference for periodicity MatS dofs (); // DOF-numbers for each component of each node (sequential) MatS dofsPeriodic (); // DOF-numbers for each component of each node (sequential) }; // ========================================= MESH ANALYSIS ========================================= // read / set the orientation (-1 / +1) of all triangles -ColI getOrientation ( const MatD &coor, const MatS &conn ); -MatS setOrientation ( const MatD &coor, const MatS &conn, int orientation=-1 ); -MatS setOrientation ( const MatD &coor, const MatS &conn, const ColI &val, int orientation=-1 ); +inline ColI getOrientation ( const MatD &coor, const MatS &conn ); +inline MatS setOrientation ( const MatD &coor, const MatS &conn, int orientation=-1 ); +inline MatS setOrientation ( const MatD &coor, const MatS &conn, const ColI &val, int orientation=-1 ); // ======================================= RE-TRIANGULATION ======================================== // simple interface to compute the full re-triangulation; it uses, depending on the input mesh: // (1) the minimal evasive "TriUpdate" // (2) the more rigorous "TriCompute" -MatS retriangulate ( const MatD &coor, const MatS &conn, int orientation=-1 ); +inline MatS retriangulate ( const MatD &coor, const MatS &conn, int orientation=-1 ); // ------------------------------------------------------------------------------------------------- // minimal evasive re-triangulation which only flips edges of the existing connectivity class TriUpdate { private: MatS m_edge; // the element that neighbors along each edge (m_nelem: no neighbor) MatS m_conn; // connectivity (updated) MatD m_coor; // nodal positions (does not change) size_t m_nelem; // #elements size_t m_nnode; // #nodes size_t m_nne; // #nodes-per-element size_t m_ndim; // #dimensions ColS m_elem; // the two elements involved in the last element change (see below) ColS m_node; // the four nodes involved in the last element change (see below) // old: m_elem(0) = [ m_node(0) , m_node(1) , m_node(2) ] // m_elem(1) = [ m_node(1) , m_node(3) , m_node(2) ] // new: m_elem(0) = [ m_node(0) , m_node(3) , m_node(2) ] // m_elem(1) = [ m_node(0) , m_node(1) , m_node(3) ] // compute neighbors per edge of all elements void edge(); // update edges around renumbered elements void chedge(size_t edge, size_t old_elem, size_t new_elem); public: TriUpdate(){}; TriUpdate(const MatD &coor, const MatS &conn); bool eval (); // re-triangulate the full mesh (returns "true" if changed) bool increment(); // one re-triangulation step (returns "true" if changed) MatS conn () { return m_conn; }; // return (new) connectivity MatS ch_elem () { return m_elem; }; // return element involved in last element change MatS ch_node () { return m_node; }; // return nodes involved in last element change }; // ================================== REGULAR MESH - SOURCE CODE =================================== -Regular::Regular(size_t nx, size_t ny, double h): m_nx(nx), m_ny(ny), m_h(h) +inline Regular::Regular(size_t nx, size_t ny, double h): m_nx(nx), m_ny(ny), m_h(h) { assert( m_nx >= 1 ); assert( m_ny >= 1 ); m_nnode = (m_nx+1) * (m_ny+1) ; m_nelem = m_nx * m_ny * 2; } // ------------------------------------------------------------------------------------------------- -MatD Regular::coor() +inline MatD Regular::coor() { MatD coor( m_nnode , m_ndim ); ColD x = ColD::LinSpaced( m_nx+1 , 0.0 , m_h * static_cast(m_nx) ); ColD y = ColD::LinSpaced( m_ny+1 , 0.0 , m_h * static_cast(m_ny) ); size_t inode = 0; for ( size_t row = 0 ; row < m_ny+1 ; ++row ) { for ( size_t col = 0 ; col < m_nx+1 ; ++col ) { coor(inode,0) = x(col); coor(inode,1) = y(row); ++inode; } } return coor; } // ------------------------------------------------------------------------------------------------- -MatS Regular::conn() +inline MatS Regular::conn() { MatS conn( m_nelem , m_nne ); size_t ielem = 0; for ( size_t row = 0 ; row < m_ny ; ++row ) { for ( size_t col = 0 ; col < m_nx ; ++col ) { conn(ielem,0) = (row+0)*(m_nx+1)+col+0; conn(ielem,1) = (row+0)*(m_nx+1)+col+1; conn(ielem,2) = (row+1)*(m_nx+1)+col+0; ++ielem; conn(ielem,0) = (row+0)*(m_nx+1)+col+1; conn(ielem,1) = (row+1)*(m_nx+1)+col+1; conn(ielem,2) = (row+1)*(m_nx+1)+col+0; ++ielem; } } return conn; } // ------------------------------------------------------------------------------------------------- -ColS Regular::nodesBottom() +inline ColS Regular::nodesBottom() { ColS nodes(m_nx+1); for ( size_t col = 0 ; col < m_nx+1 ; ++col ) nodes(col) = col; return nodes; } // ------------------------------------------------------------------------------------------------- -ColS Regular::nodesTop() +inline ColS Regular::nodesTop() { ColS nodes(m_nx+1); for ( size_t col = 0 ; col < m_nx+1 ; ++col ) nodes(col) = col+m_ny*(m_nx+1); return nodes; } // ------------------------------------------------------------------------------------------------- -ColS Regular::nodesLeft() +inline ColS Regular::nodesLeft() { ColS nodes(m_ny+1); for ( size_t row = 0 ; row < m_ny+1 ; ++row ) nodes(row) = row*(m_nx+1); return nodes; } // ------------------------------------------------------------------------------------------------- -ColS Regular::nodesRight() +inline ColS Regular::nodesRight() { ColS nodes(m_ny+1); for ( size_t row = 0 ; row < m_ny+1 ; ++row ) nodes(row) = row*(m_nx+1)+m_nx; return nodes; } // ------------------------------------------------------------------------------------------------- -MatS Regular::nodesPeriodic() +inline MatS Regular::nodesPeriodic() { ColS bot = nodesBottom(); ColS top = nodesTop (); ColS rgt = nodesRight (); ColS lft = nodesLeft (); MatS nodes( bot.size()-2 + lft.size()-2 + 3 , 2 ); size_t i = 0; nodes(i,0) = bot(0); nodes(i,1) = bot(bot.size()-1); ++i; nodes(i,0) = bot(0); nodes(i,1) = top(top.size()-1); ++i; nodes(i,0) = bot(0); nodes(i,1) = top(0 ); ++i; for ( auto j = 1 ; j < bot.size()-1 ; ++j ) { nodes(i,0) = bot(j); nodes(i,1) = top(j); ++i; } for ( auto j = 1 ; j < lft.size()-1 ; ++j ) { nodes(i,0) = lft(j); nodes(i,1) = rgt(j); ++i; } return nodes; } // ------------------------------------------------------------------------------------------------- -size_t Regular::nodeOrigin() +inline size_t Regular::nodeOrigin() { return 0; } // ------------------------------------------------------------------------------------------------- -MatS Regular::dofs() +inline MatS Regular::dofs() { return GooseFEM::Mesh::dofs(m_nnode,m_ndim); } // ------------------------------------------------------------------------------------------------- -MatS Regular::dofsPeriodic() +inline MatS Regular::dofsPeriodic() { // DOF-numbers for each component of each node (sequential) MatS out = GooseFEM::Mesh::dofs(m_nnode,m_ndim); // periodic node-pairs MatS nodePer = nodesPeriodic(); size_t nper = static_cast(nodePer.rows()); // eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs for ( size_t i = 0 ; i < nper ; ++i ) for ( size_t j = 0 ; j < m_ndim ; ++j ) out(nodePer(i,1),j) = out(nodePer(i,0),j); // renumber "out" to be sequential return GooseFEM::Mesh::renumber(out); } // ================================== MESH ANALYSIS - SOURCE CODE ================================== -ColI getOrientation ( const MatD &coor, const MatS &conn ) +inline ColI getOrientation ( const MatD &coor, const MatS &conn ) { assert( conn.cols() == 3 ); assert( coor.cols() == 2 ); Eigen::Vector2d v1,v2; double k; size_t nelem = static_cast( conn.rows() ); ColI out( nelem ); for ( size_t ielem = 0 ; ielem < nelem ; ++ielem ) { v1 = coor.row( conn(ielem,0) ) - coor.row( conn(ielem,1) ); v2 = coor.row( conn(ielem,2) ) - coor.row( conn(ielem,1) ); k = v1(0) * v2(1) - v2(0) * v1(1); if ( k < 0 ) out( ielem ) = -1; else out( ielem ) = +1; } return out; } // ------------------------------------------------------------------------------------------------- -MatS setOrientation ( const MatD &coor, const MatS &conn, int orientation ) +inline MatS setOrientation ( const MatD &coor, const MatS &conn, int orientation ) { assert( conn.cols() == 3 ); assert( coor.cols() == 2 ); assert( orientation == -1 || orientation == +1 ); ColI val = getOrientation(coor, conn); return setOrientation( coor, conn, val, orientation ); } // ------------------------------------------------------------------------------------------------- -MatS setOrientation ( const MatD &coor, const MatS &conn, const ColI &val, int orientation ) +inline MatS setOrientation ( const MatD &coor, const MatS &conn, const ColI &val, int orientation ) { assert( conn.cols() == 3 ); assert( coor.cols() == 2 ); assert( conn.rows() == val.size() ); assert( orientation == -1 || orientation == +1 ); // avoid compiler warning UNUSED(coor); size_t nelem = static_cast(conn.rows()); MatS out = conn; for ( size_t ielem = 0 ; ielem < nelem ; ++ielem ) if ( ( orientation == -1 and val(ielem) > 0 ) or ( orientation == +1 and val(ielem) < 0 ) ) std::swap( out(ielem,2) , out(ielem,1) ); return out; } // ================================ RE-TRIANGULATION - SOURCE CODE ================================= -MatS retriangulate ( const MatD &coor, const MatS &conn, int orientation ) +inline MatS retriangulate ( const MatD &coor, const MatS &conn, int orientation ) { // get the orientation of all elements ColI dir = getOrientation( coor, conn ); // check the orientation bool eq = std::abs( dir.sum() ) == conn.rows(); // new connectivity MatS out; // perform re-triangulation // - use "TriUpdate" if ( eq ) { TriUpdate obj(coor,conn); obj.eval(); out = obj.conn(); } // - using TriCompute else { throw std::runtime_error("Work-in-progress, has to be re-triangulated using 'TriCompute'"); } return setOrientation(coor,out,orientation); } // ================================================================================================= // support class to allow the storage of a list of edges class Edge { public: size_t n1 ; // node 1 (edge from node 1-2) size_t n2 ; // node 2 (edge from node 1-2) size_t elem; // element to which the edge belong size_t edge; // edge index within the element (e.g. edge==1 -> n1=conn(0,elem), n2=conn(1,elem)) Edge(){} Edge(size_t i, size_t j, size_t el, size_t ed, bool sort=false): n1(i), n2(j), elem(el), edge(ed) { if ( sort && n1>n2 ) { std::swap(n1,n2); } } }; // ------------------------------------------------------------------------------------------------- -bool Edge_cmp( Edge a , Edge b ) +inline bool Edge_cmp( Edge a , Edge b ) { if ( a.n1 == b.n1 and a.n2 == b.n2 ) return true; return false; } // ------------------------------------------------------------------------------------------------- -bool Edge_sort( Edge a , Edge b ) +inline bool Edge_sort( Edge a , Edge b ) { if ( a.n1 < b.n1 or a.n2 < b.n2 ) return true; return false; } // ================================================================================================= -TriUpdate::TriUpdate(const MatD &coor, const MatS &conn): m_conn(conn), m_coor(coor) +inline TriUpdate::TriUpdate(const MatD &coor, const MatS &conn): m_conn(conn), m_coor(coor) { assert( conn.cols() == 3 ); assert( coor.cols() == 2 ); // store shapes m_nnode = static_cast( coor.rows() ); m_ndim = static_cast( coor.cols() ); m_nelem = static_cast( conn.rows() ); m_nne = static_cast( conn.cols() ); // resize internal arrays m_elem.resize(2); m_node.resize(4); // set default to out-of-bounds, to make clear that nothing happened yet m_elem.setConstant(m_nelem); m_node.setConstant(m_nnode); edge(); } // ------------------------------------------------------------------------------------------------- -void TriUpdate::edge() +inline void TriUpdate::edge() { m_edge.resize( m_nelem , m_nne ); m_edge.setConstant( m_nelem ); // signal that nothing has been set std::vector idx = {0,1,2}; // lists to convert connectivity -> edges std::vector jdx = {1,2,0}; // lists to convert connectivity -> edges std::vector edge; edge.reserve( m_nelem*idx.size() ); for ( size_t e = 0 ; e < m_nelem ; ++e ) for ( size_t i = 0 ; i < m_nne ; ++i ) edge.push_back( Edge( m_conn(e,idx[i]), m_conn(e,jdx[i]) , e , i , true ) ); std::sort( edge.begin() , edge.end() , Edge_sort ); for ( size_t i = 0 ; i < edge.size()-1 ; ++i ) { if ( edge[i].n1 == edge[i+1].n1 and edge[i].n2 == edge[i+1].n2 ) { m_edge( edge[i ].elem , edge[i ].edge ) = edge[i+1].elem; m_edge( edge[i+1].elem , edge[i+1].edge ) = edge[i ].elem; } } } // ------------------------------------------------------------------------------------------------- -void TriUpdate::chedge(size_t edge, size_t old_elem, size_t new_elem) +inline void TriUpdate::chedge(size_t edge, size_t old_elem, size_t new_elem) { size_t m; size_t neigh = m_edge( old_elem , edge ); if ( neigh >= m_nelem ) return; for ( m = 0 ; m < m_nne ; ++m ) if ( m_edge( neigh , m ) == old_elem ) break; m_edge( neigh , m ) = new_elem; } // ------------------------------------------------------------------------------------------------- -bool TriUpdate::increment() +inline bool TriUpdate::increment() { size_t ielem,jelem,iedge,jedge; double phi1,phi2; ColS c(4); ColS n(4); // loop over all elements for ( ielem = 0 ; ielem < m_nelem ; ++ielem ) { // loop over all edges for ( iedge = 0 ; iedge < m_nne ; ++iedge ) { // only proceed if the edge is shared with another element if ( m_edge(ielem,iedge) >= m_nelem ) continue; // read "jelem" jelem = m_edge(ielem,iedge); // find the edge involved for "jelem" for ( jedge=0; jedge M_PI ) { // update connectivity m_conn(ielem,0) = c(0); m_conn(ielem,1) = c(3); m_conn(ielem,2) = c(2); m_conn(jelem,0) = c(0); m_conn(jelem,1) = c(1); m_conn(jelem,2) = c(3); // change list with neighbors for the elements around (only two neighbors change) if ( iedge==0 ) { chedge(2,ielem,jelem); } else if ( iedge==1 ) { chedge(0,ielem,jelem); } else if ( iedge==2 ) { chedge(1,ielem,jelem); } if ( jedge==0 ) { chedge(2,jelem,ielem); } else if ( jedge==1 ) { chedge(0,jelem,ielem); } else if ( jedge==2 ) { chedge(1,jelem,ielem); } // convert to four static nodes if ( iedge==0 ) { n(0)=m_edge(ielem,2); n(3)=m_edge(ielem,1); } else if ( iedge==1 ) { n(0)=m_edge(ielem,0); n(3)=m_edge(ielem,2); } else if ( iedge==2 ) { n(0)=m_edge(ielem,1); n(3)=m_edge(ielem,0); } if ( jedge==0 ) { n(1)=m_edge(jelem,1); n(2)=m_edge(jelem,2); } else if ( jedge==1 ) { n(1)=m_edge(jelem,2); n(2)=m_edge(jelem,0); } else if ( jedge==2 ) { n(1)=m_edge(jelem,0); n(2)=m_edge(jelem,1); } // store the neighbors for the changed elements m_edge(ielem,0) = jelem; m_edge(jelem,0) = n(0) ; m_edge(ielem,1) = n(2) ; m_edge(jelem,1) = n(1) ; m_edge(ielem,2) = n(3) ; m_edge(jelem,2) = ielem; // store information for transfer algorithm m_node = c; m_elem(0) = ielem; m_elem(1) = jelem; return true; } } } return false; } // ------------------------------------------------------------------------------------------------- -bool TriUpdate::eval() +inline bool TriUpdate::eval() { bool change = false; while ( increment() ) { change = true; } return change; } // ================================================================================================= } // namespace Tri3 } // namespace Mesh } // namespace GooseFEM #endif