diff --git a/src/GooseFEM/MeshQuad4.hpp b/src/GooseFEM/MeshQuad4.hpp index 3a699d5..054a3c6 100644 --- a/src/GooseFEM/MeshQuad4.hpp +++ b/src/GooseFEM/MeshQuad4.hpp @@ -1,922 +1,921 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MESHQUAD4_CPP #define GOOSEFEM_MESHQUAD4_CPP // ------------------------------------------------------------------------------------------------- #include "MeshQuad4.h" // ===================================== GooseFEM::Mesh::Quad4 ===================================== namespace GooseFEM { namespace Mesh { namespace Quad4 { // ------------------------------------------ constructor ------------------------------------------ inline Regular::Regular(size_t nelx, size_t nely, double h): m_h(h), m_nelx(nelx), m_nely(nely) { assert( m_nelx >= 1 ); assert( m_nely >= 1 ); m_nnode = (m_nelx+1) * (m_nely+1); m_nelem = m_nelx * m_nely ; } // -------------------------------------- number of elements --------------------------------------- inline size_t Regular::nelem() const { return m_nelem; } // ---------------------------------------- number of nodes ---------------------------------------- inline size_t Regular::nnode() const { return m_nnode; } // ---------------------------------- number of nodes per element ---------------------------------- inline size_t Regular::nne() const { return m_nne; } // ------------------------------------- number of dimensions -------------------------------------- inline size_t Regular::ndim() const { return m_ndim; } // ------------------------ number of nodes, after eliminating periodicity ------------------------- inline size_t Regular::nnodePeriodic() const { return (m_nelx+1) * (m_nely+1) - (m_nely+1) - (m_nelx); } // --------------------------------- coordinates (nodal positions) --------------------------------- inline xt::xtensor Regular::coor() const { xt::xtensor out = xt::empty({m_nnode, m_ndim}); xt::xtensor x = xt::linspace(0.0, m_h*static_cast(m_nelx), m_nelx+1); xt::xtensor y = xt::linspace(0.0, m_h*static_cast(m_nely), m_nely+1); size_t inode = 0; for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) { - for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) { out(inode,0) = x(ix); out(inode,1) = y(iy); ++inode; } } return out; } // ---------------------------- connectivity (node-numbers per element) ---------------------------- inline xt::xtensor Regular::conn() const { xt::xtensor out = xt::empty({m_nelem,m_nne}); size_t ielem = 0; for ( size_t iy = 0 ; iy < m_nely ; ++iy ) { for ( size_t ix = 0 ; ix < m_nelx ; ++ix ) { out(ielem,0) = (iy )*(m_nelx+1) + (ix ); out(ielem,1) = (iy )*(m_nelx+1) + (ix+1); out(ielem,3) = (iy+1)*(m_nelx+1) + (ix ); out(ielem,2) = (iy+1)*(m_nelx+1) + (ix+1); ++ielem; } } return out; } // ------------------------------ node-numbers along the bottom edge ------------------------------- inline xt::xtensor Regular::nodesBottomEdge() const { xt::xtensor out = xt::empty({m_nelx+1}); for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) out(ix) = ix; return out; } // -------------------------------- node-numbers along the top edge -------------------------------- inline xt::xtensor Regular::nodesTopEdge() const { xt::xtensor out = xt::empty({m_nelx+1}); for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) out(ix) = ix + m_nely*(m_nelx+1); return out; } // ------------------------------- node-numbers along the left edge -------------------------------- inline xt::xtensor Regular::nodesLeftEdge() const { xt::xtensor out = xt::empty({m_nely+1}); for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) out(iy) = iy*(m_nelx+1); return out; } // ------------------------------- node-numbers along the right edge ------------------------------- inline xt::xtensor Regular::nodesRightEdge() const { xt::xtensor out = xt::empty({m_nely+1}); for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) out(iy) = iy*(m_nelx+1) + m_nelx; return out; } // ---------------------- node-numbers along the bottom edge, without corners ---------------------- inline xt::xtensor Regular::nodesBottomOpenEdge() const { xt::xtensor out = xt::empty({m_nelx-1}); for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) out(ix-1) = ix; return out; } // ----------------------- node-numbers along the top edge, without corners ------------------------ inline xt::xtensor Regular::nodesTopOpenEdge() const { xt::xtensor out = xt::empty({m_nelx-1}); for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) out(ix-1) = ix + m_nely*(m_nelx+1); return out; } // ----------------------- node-numbers along the left edge, without corners ----------------------- inline xt::xtensor Regular::nodesLeftOpenEdge() const { xt::xtensor out = xt::empty({m_nely-1}); for ( size_t iy = 1 ; iy < m_nely ; ++iy ) out(iy-1) = iy*(m_nelx+1); return out; } // ---------------------- node-numbers along the right edge, without corners ----------------------- inline xt::xtensor Regular::nodesRightOpenEdge() const { xt::xtensor out = xt::empty({m_nely-1}); for ( size_t iy = 1 ; iy < m_nely ; ++iy ) out(iy-1) = iy*(m_nelx+1) + m_nelx; return out; } // ----------------------------- node-number of the bottom-left corner ----------------------------- inline size_t Regular::nodesBottomLeftCorner() const { return 0; } // ---------------------------- node-number of the bottom-right corner ----------------------------- inline size_t Regular::nodesBottomRightCorner() const { return m_nelx; } // ------------------------------ node-number of the top-left corner ------------------------------- inline size_t Regular::nodesTopLeftCorner() const { return m_nely*(m_nelx+1); } // ------------------------------ node-number of the top-right corner ------------------------------ inline size_t Regular::nodesTopRightCorner() const { return m_nely*(m_nelx+1) + m_nelx; } // ----------------------------- node-number of the corners (aliases) ------------------------------ inline size_t Regular::nodesLeftBottomCorner() const { return nodesBottomLeftCorner(); } inline size_t Regular::nodesLeftTopCorner() const { return nodesTopLeftCorner(); } inline size_t Regular::nodesRightBottomCorner() const { return nodesBottomRightCorner(); } inline size_t Regular::nodesRightTopCorner() const { return nodesTopRightCorner(); } // ------------------------------ node-numbers of periodic node-pairs ------------------------------ inline xt::xtensor Regular::nodesPeriodic() const { // edges (without corners) xt::xtensor bot = nodesBottomOpenEdge(); xt::xtensor top = nodesTopOpenEdge(); xt::xtensor lft = nodesLeftOpenEdge(); xt::xtensor rgt = nodesRightOpenEdge(); // allocate nodal ties // - number of tying per category size_t tedge = bot.size() + lft.size(); size_t tnode = 3; // - allocate xt::xtensor out = xt::empty({tedge+tnode, std::size_t(2)}); // counter size_t i = 0; // tie all corners to one corner out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesBottomRightCorner(); ++i; out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopRightCorner(); ++i; out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopLeftCorner(); ++i; // tie all corresponding edges to each other for ( size_t j = 0 ; j Regular::dofs() const { return GooseFEM::Mesh::dofs(m_nnode,m_ndim); } // ------------------------ DOP-numbers with periodic dependencies removed ------------------------- inline xt::xtensor Regular::dofsPeriodic() const { // DOF-numbers for each component of each node (sequential) xt::xtensor out = GooseFEM::Mesh::dofs(m_nnode,m_ndim); // periodic node-pairs xt::xtensor nodePer = nodesPeriodic(); // eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs for ( size_t i = 0 ; i < nodePer.shape()[0] ; ++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); } // ------------------------------------------ constructor ------------------------------------------ inline FineLayer::FineLayer(size_t nelx, size_t nely, double h, size_t nfine): m_h(h) { // basic assumptions assert( nelx >= 1 ); assert( nely >= 1 ); // store basic info m_Lx = m_h * static_cast(nelx); // compute element size in y-direction (use symmetry, compute upper half) // ---------------------------------------------------------------------- // temporary variables size_t nmin, ntot; xt::xtensor nhx = xt::ones({nely}); xt::xtensor nhy = xt::ones({nely}); xt::xtensor refine = -1 * xt::ones ({nely}); // minimum height in y-direction (half of the height because of symmetry) if ( nely % 2 == 0 ) nmin = nely /2; else nmin = (nely +1)/2; // minimum number of fine layers in y-direction (minimum 1, middle layer part of this half) if ( nfine % 2 == 0 ) nfine = nfine /2 + 1; else nfine = (nfine+1)/2; if ( nfine < 1 ) nfine = 1; if ( nfine > nmin ) nfine = nmin; // loop over element layers in y-direction, try to coarsen using these rules: // (1) element size in y-direction <= distance to origin in y-direction // (2) element size in x-direction should fit the total number of elements in x-direction // (3) a certain number of layers have the minimum size "1" (are fine) for ( size_t iy = nfine ; ; ) { // initialize current size in y-direction if ( iy == nfine ) ntot = nfine; // check to stop if ( iy >= nely or ntot >= nmin ) { nely = iy; break; } // rules (1,2) satisfied: coarsen in x-direction if ( 3*nhy(iy) <= ntot and nelx%(3*nhx(iy)) == 0 and ntot+nhy(iy) < nmin ) { refine(iy) = 0; nhy (iy) *= 2; auto vnhy = xt::view(nhy, xt::range(iy+1, _)); auto vnhx = xt::view(nhx, xt::range(iy , _)); vnhy *= 3; vnhx *= 3; } // update the number of elements in y-direction ntot += nhy(iy); // proceed to next element layer in y-direction ++iy; // check to stop if ( iy >= nely or ntot >= nmin ) { nely = iy; break; } } // symmetrize, compute full information // ------------------------------------ // allocate mesh constructor parameters m_nhx = xt::empty({nely*2-1}); m_nhy = xt::empty({nely*2-1}); m_refine = xt::empty ({nely*2-1}); m_nelx = xt::empty({nely*2-1}); m_nnd = xt::empty({nely*2 }); m_startElem = xt::empty({nely*2-1}); m_startNode = xt::empty({nely*2 }); // fill // - lower half for ( size_t iy = 0 ; iy < nely ; ++iy ) { m_nhx (iy) = nhx (nely-iy-1); m_nhy (iy) = nhy (nely-iy-1); m_refine(iy) = refine(nely-iy-1); } // - upper half for ( size_t iy = 0 ; iy < nely-1 ; ++iy ) { m_nhx (iy+nely) = nhx (iy+1); m_nhy (iy+nely) = nhy (iy+1); m_refine(iy+nely) = refine(iy+1); } // update size nely = m_nhx.size(); // compute the number of elements per element layer in y-direction for ( size_t iy = 0 ; iy < nely ; ++iy ) m_nelx(iy) = nelx / m_nhx(iy); // compute the number of nodes per node layer in y-direction for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) m_nnd(iy ) = m_nelx(iy)+1; for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) m_nnd(iy+1) = m_nelx(iy)+1; // compute mesh dimensions // ----------------------- // initialize m_nnode = 0; m_nelem = 0; m_startNode(0) = 0; // loop over element layers (bottom -> middle, elements become finer) for ( size_t i = 0 ; i < (nely-1)/2 ; ++i ) { // - store the first element of the layer m_startElem(i) = m_nelem; // - add the nodes of this layer if ( m_refine(i) == 0 ) { m_nnode += (3*m_nelx(i)+1); } else { m_nnode += ( m_nelx(i)+1); } // - add the elements of this layer if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ); } else { m_nelem += ( m_nelx(i) ); } // - store the starting node of the next layer m_startNode(i+1) = m_nnode; } // loop over element layers (middle -> top, elements become coarser) for ( size_t i = (nely-1)/2 ; i < nely ; ++i ) { // - store the first element of the layer m_startElem(i) = m_nelem; // - add the nodes of this layer if ( m_refine(i) == 0 ) { m_nnode += (5*m_nelx(i)+1); } else { m_nnode += ( m_nelx(i)+1); } // - add the elements of this layer if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ); } else { m_nelem += ( m_nelx(i) ); } // - store the starting node of the next layer m_startNode(i+1) = m_nnode; } // - add the top row of nodes m_nnode += m_nelx(nely-1)+1; } // -------------------------------------- number of elements --------------------------------------- inline size_t FineLayer::nelem() const { return m_nelem; } // ---------------------------------------- number of nodes ---------------------------------------- inline size_t FineLayer::nnode() const { return m_nnode; } // ---------------------------------- number of nodes per element ---------------------------------- inline size_t FineLayer::nne() const { return m_nne; } // ------------------------------------- number of dimensions -------------------------------------- inline size_t FineLayer::ndim() const { return m_ndim; } // ---------------------------- actual number of nodes in one direction ---------------------------- inline size_t FineLayer::shape(size_t i) const { assert( i >= 0 and i <= 1 ); if ( i == 0 ) return xt::amax(m_nelx)[0]; else return xt::sum (m_nhy )[0]; } // --------------------------------- coordinates (nodal positions) --------------------------------- inline xt::xtensor FineLayer::coor() const { // allocate output xt::xtensor out = xt::empty({m_nnode, m_ndim}); // current node, number of element layers size_t inode = 0; size_t nely = static_cast(m_nhy.size()); // y-position of each main node layer (i.e. excluding node layers for refinement/coarsening) // - allocate xt::xtensor y = xt::empty({nely+1}); // - initialize y(0) = 0.0; // - compute for ( size_t iy = 1 ; iy < nely+1 ; ++iy ) y(iy) = y(iy-1) + m_nhy(iy-1) * m_h; // loop over element layers (bottom -> middle) : add bottom layer (+ refinement layer) of nodes // -------------------------------------------------------------------------------------------- for ( size_t iy = 0 ; ; ++iy ) { // get positions along the x- and z-axis xt::xtensor x = xt::linspace(0.0, m_Lx, m_nelx(iy)+1); // add nodes of the bottom layer of this element for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { out(inode,0) = x(ix); out(inode,1) = y(iy); ++inode; } // stop at middle layer if ( iy == (nely-1)/2 ) break; // add extra nodes of the intermediate layer, for refinement in x-direction if ( m_refine(iy) == 0 ) { // - get position offset in x- and y-direction double dx = m_h * static_cast(m_nhx(iy)/3); double dy = m_h * static_cast(m_nhy(iy)/2); // - add nodes of the intermediate layer for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { for ( size_t j = 0 ; j < 2 ; ++j ) { out(inode,0) = x(ix) + dx * static_cast(j+1); out(inode,1) = y(iy) + dy; ++inode; } } } } // loop over element layers (middle -> top) : add (refinement layer +) top layer of nodes // -------------------------------------------------------------------------------------- for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) { // get positions along the x- and z-axis xt::xtensor x = xt::linspace(0.0, m_Lx, m_nelx(iy)+1); // add extra nodes of the intermediate layer, for refinement in x-direction if ( m_refine(iy) == 0 ) { // - get position offset in x- and y-direction double dx = m_h * static_cast(m_nhx(iy)/3); double dy = m_h * static_cast(m_nhy(iy)/2); // - add nodes of the intermediate layer for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { for ( size_t j = 0 ; j < 2 ; ++j ) { out(inode,0) = x(ix) + dx * static_cast(j+1); out(inode,1) = y(iy) + dy; ++inode; } } } // add nodes of the top layer of this element for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) { out(inode,0) = x(ix ); out(inode,1) = y(iy+1); ++inode; } } return out; } // ---------------------------- connectivity (node-numbers per element) ---------------------------- inline xt::xtensor FineLayer::conn() const { // allocate output xt::xtensor out = xt::empty({m_nelem, m_nne}); // current element, number of element layers, starting nodes of each node layer size_t ielem = 0; size_t nely = static_cast(m_nhy.size()); size_t bot,mid,top; // loop over all element layers for ( size_t iy = 0 ; iy < nely ; ++iy ) { // - get: starting nodes of bottom(, middle) and top layer bot = m_startNode(iy ); mid = m_startNode(iy ) + m_nnd(iy); top = m_startNode(iy+1); // - define connectivity: no coarsening/refinement if ( m_refine(iy) == -1 ) { for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { out(ielem,0) = bot + (ix ); out(ielem,1) = bot + (ix+1); out(ielem,2) = top + (ix+1); out(ielem,3) = top + (ix ); ielem++; } } // - define connectivity: refinement along the x-direction (below the middle layer) else if ( m_refine(iy) == 0 and iy <= (nely-1)/2 ) { for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { // -- bottom element out(ielem,0) = bot + ( ix ); out(ielem,1) = bot + ( ix+1); out(ielem,2) = mid + (2*ix+1); out(ielem,3) = mid + (2*ix ); ielem++; // -- top-right element out(ielem,0) = bot + ( ix+1); out(ielem,1) = top + (3*ix+3); out(ielem,2) = top + (3*ix+2); out(ielem,3) = mid + (2*ix+1); ielem++; // -- top-center element out(ielem,0) = mid + (2*ix ); out(ielem,1) = mid + (2*ix+1); out(ielem,2) = top + (3*ix+2); out(ielem,3) = top + (3*ix+1); ielem++; // -- top-left element out(ielem,0) = bot + ( ix ); out(ielem,1) = mid + (2*ix ); out(ielem,2) = top + (3*ix+1); out(ielem,3) = top + (3*ix ); ielem++; } } // - define connectivity: coarsening along the x-direction (above the middle layer) else if ( m_refine(iy) == 0 and iy > (nely-1)/2 ) { for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) { // -- lower-left element out(ielem,0) = bot + (3*ix ); out(ielem,1) = bot + (3*ix+1); out(ielem,2) = mid + (2*ix ); out(ielem,3) = top + ( ix ); ielem++; // -- lower-center element out(ielem,0) = bot + (3*ix+1); out(ielem,1) = bot + (3*ix+2); out(ielem,2) = mid + (2*ix+1); out(ielem,3) = mid + (2*ix ); ielem++; // -- lower-right element out(ielem,0) = bot + (3*ix+2); out(ielem,1) = bot + (3*ix+3); out(ielem,2) = top + ( ix+1); out(ielem,3) = mid + (2*ix+1); ielem++; // -- upper element out(ielem,0) = mid + (2*ix ); out(ielem,1) = mid + (2*ix+1); out(ielem,2) = top + ( ix+1); out(ielem,3) = top + ( ix ); ielem++; } } } return out; } // ------------------------------ element numbers of the middle layer ------------------------------ inline xt::xtensor FineLayer::elementsMiddleLayer() const { // number of element layers in y-direction, the index of the middle layer size_t nely = static_cast(m_nhy.size()); size_t iy = (nely-1)/2; xt::xtensor out = xt::empty({m_nelx(iy)}); for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) out(ix) = m_startElem(iy) + ix; return out; } // ------------------------------ node-numbers along the bottom edge ------------------------------- inline xt::xtensor FineLayer::nodesBottomEdge() const { xt::xtensor out = xt::empty({m_nelx(0)+1}); for ( size_t ix = 0 ; ix < m_nelx(0)+1 ; ++ix ) out(ix) = m_startNode(0) + ix; return out; } // -------------------------------- node-numbers along the top edge -------------------------------- inline xt::xtensor FineLayer::nodesTopEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor out = xt::empty({m_nelx(nely-1)+1}); for ( size_t ix = 0 ; ix < m_nelx(nely-1)+1 ; ++ix ) out(ix) = m_startNode(nely) + ix; return out; } // ------------------------------- node-numbers along the left edge -------------------------------- inline xt::xtensor FineLayer::nodesLeftEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor out = xt::empty({nely+1}); for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) out(iy) = m_startNode(iy); for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) out(iy+1) = m_startNode(iy+1); return out; } // ------------------------------- node-numbers along the right edge ------------------------------- inline xt::xtensor FineLayer::nodesRightEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor out = xt::empty({nely+1}); for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy ) out(iy) = m_startNode(iy) + m_nelx(iy); for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy ) out(iy+1) = m_startNode(iy+1) + m_nelx(iy); return out; } // ---------------------- node-numbers along the bottom edge, without corners ---------------------- inline xt::xtensor FineLayer::nodesBottomOpenEdge() const { xt::xtensor out = xt::empty({m_nelx(0)-1}); for ( size_t ix = 1 ; ix < m_nelx(0) ; ++ix ) out(ix-1) = m_startNode(0) + ix; return out; } // ----------------------- node-numbers along the top edge, without corners ------------------------ inline xt::xtensor FineLayer::nodesTopOpenEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor out = xt::empty({m_nelx(nely-1)-1}); for ( size_t ix = 1 ; ix < m_nelx(nely-1) ; ++ix ) out(ix-1) = m_startNode(nely) + ix; return out; } // ----------------------- node-numbers along the left edge, without corners ----------------------- inline xt::xtensor FineLayer::nodesLeftOpenEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor out = xt::empty({nely-1}); for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) out(iy-1) = m_startNode(iy); for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) out(iy) = m_startNode(iy+1); return out; } // ---------------------- node-numbers along the right edge, without corners ----------------------- inline xt::xtensor FineLayer::nodesRightOpenEdge() const { size_t nely = static_cast(m_nhy.size()); xt::xtensor out = xt::empty({nely-1}); for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy ) out(iy-1) = m_startNode(iy) + m_nelx(iy); for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy ) out(iy) = m_startNode(iy+1) + m_nelx(iy); return out; } // ----------------------------- node-number of the bottom-left corner ----------------------------- inline size_t FineLayer::nodesBottomLeftCorner() const { return m_startNode(0); } // ---------------------------- node-number of the bottom-right corner ----------------------------- inline size_t FineLayer::nodesBottomRightCorner() const { return m_startNode(0) + m_nelx(0); } // ------------------------------ node-number of the top-left corner ------------------------------- inline size_t FineLayer::nodesTopLeftCorner() const { size_t nely = static_cast(m_nhy.size()); return m_startNode(nely); } // ------------------------------ node-number of the top-right corner ------------------------------ inline size_t FineLayer::nodesTopRightCorner() const { size_t nely = static_cast(m_nhy.size()); return m_startNode(nely) + m_nelx(nely-1); } // -------------------------------------------- aliases -------------------------------------------- inline size_t FineLayer::nodesLeftBottomCorner() const { return nodesBottomLeftCorner(); } inline size_t FineLayer::nodesRightBottomCorner() const { return nodesBottomRightCorner(); } inline size_t FineLayer::nodesLeftTopCorner() const { return nodesTopLeftCorner(); } inline size_t FineLayer::nodesRightTopCorner() const { return nodesTopRightCorner(); } // ------------------------------ node-numbers of periodic node-pairs ------------------------------ inline xt::xtensor FineLayer::nodesPeriodic() const { // edges (without corners) xt::xtensor bot = nodesBottomOpenEdge(); xt::xtensor top = nodesTopOpenEdge(); xt::xtensor lft = nodesLeftOpenEdge(); xt::xtensor rgt = nodesRightOpenEdge(); // allocate nodal ties // - number of tying per category size_t tedge = bot.size() + lft.size(); size_t tnode = 3; // - allocate xt::xtensor out = xt::empty({tedge+tnode, std::size_t(2)}); // counter size_t i = 0; // tie all corners to one corner out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesBottomRightCorner(); ++i; out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopRightCorner(); ++i; out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopLeftCorner(); ++i; // tie all corresponding edges to each other for ( size_t j = 0 ; j FineLayer::dofs() const { return GooseFEM::Mesh::dofs(m_nnode,m_ndim); } // ------------------------ DOP-numbers with periodic dependencies removed ------------------------- inline xt::xtensor FineLayer::dofsPeriodic() const { // DOF-numbers for each component of each node (sequential) xt::xtensor out = GooseFEM::Mesh::dofs(m_nnode,m_ndim); // periodic node-pairs xt::xtensor nodePer = nodesPeriodic(); // eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs for ( size_t i = 0 ; i < nodePer.shape()[0] ; ++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 ... // ================================================================================================= #endif diff --git a/src/GooseFEM/MeshTri3.h b/src/GooseFEM/MeshTri3.h index c1b41cb..419b19d 100644 --- a/src/GooseFEM/MeshTri3.h +++ b/src/GooseFEM/MeshTri3.h @@ -1,155 +1,154 @@ /* ================================================================================================= (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 "Mesh.h" // ===================================== GooseFEM::Mesh::Tri3 ====================================== namespace GooseFEM { namespace Mesh { namespace Tri3 { // ------------------------------------------ regular mesh ----------------------------------------- class Regular { private: - double m_h; // elementary element edge-size (in both directions) - size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h") - size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h") - 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 + double m_h; // elementary element edge-size (in both directions) + size_t m_nelx; // number of elements in x-direction (length == "m_nelx * m_h") + size_t m_nely; // number of elements in y-direction (length == "m_nely * m_h") + size_t m_nelem; // number of elements + size_t m_nnode; // number of nodes + static const size_t m_nne=3; // number of nodes-per-element + static const size_t m_ndim=2; // number of dimensions public: // mesh with "2*nelx*nely" 'elements' of edge size "h" Regular(size_t nelx, size_t nely, double h=1.); // sizes - size_t nelem() const; // number of elements - size_t nnode() const; // number of nodes - size_t nne() const; // number of nodes-per-element - size_t ndim() const; // number of dimensions + size_t nelem() const; // number of elements + size_t nnode() const; // number of nodes + size_t nne() const; // number of nodes-per-element + size_t ndim() const; // number of dimensions // mesh - MatD coor() const; // nodal positions [nnode ,ndim] - MatS conn() const; // connectivity [nelem ,nne ] + xt::xtensor coor() const; // nodal positions [nnode ,ndim] + xt::xtensor conn() const; // connectivity [nelem ,nne ] // boundary nodes: edges - ColS nodesBottomEdge() const; // node-numbers along the bottom edge - ColS nodesTopEdge() const; // node-numbers along the top edge - ColS nodesLeftEdge() const; // node-numbers along the left edge - ColS nodesRightEdge() const; // node-numbers along the right edge + xt::xtensor nodesBottomEdge() const; // node-numbers along the bottom edge + xt::xtensor nodesTopEdge() const; // node-numbers along the top edge + xt::xtensor nodesLeftEdge() const; // node-numbers along the left edge + xt::xtensor nodesRightEdge() const; // node-numbers along the right edge // boundary nodes: edges, without corners - ColS nodesBottomOpenEdge() const; // node-numbers along the bottom edge - ColS nodesTopOpenEdge() const; // node-numbers along the top edge - ColS nodesLeftOpenEdge() const; // node-numbers along the left edge - ColS nodesRightOpenEdge() const; // node-numbers along the right edge + xt::xtensor nodesBottomOpenEdge() const; // node-numbers along the bottom edge + xt::xtensor nodesTopOpenEdge() const; // node-numbers along the top edge + xt::xtensor nodesLeftOpenEdge() const; // node-numbers along the left edge + xt::xtensor nodesRightOpenEdge() const; // node-numbers along the right edge // boundary nodes: corners size_t nodesBottomLeftCorner() const; // node-number of the bottom - left corner size_t nodesBottomRightCorner() const; // node-number of the bottom - right corner size_t nodesTopLeftCorner() const; // node-number of the top - left corner size_t nodesTopRightCorner() const; // node-number of the top - right corner // boundary nodes: corners (aliases) size_t nodesLeftBottomCorner() const; // alias, see above: nodesBottomLeftCorner size_t nodesLeftTopCorner() const; // alias, see above: nodesBottomRightCorner size_t nodesRightBottomCorner() const; // alias, see above: nodesTopLeftCorner size_t nodesRightTopCorner() const; // alias, see above: nodesTopRightCorner // periodicity - MatS nodesPeriodic() const; // periodic node pairs [:,2]: (independent, dependent) - size_t nodesOrigin() const; // bottom-left node, used as reference for periodicity - MatS dofs() const; // DOF-numbers for each component of each node (sequential) - MatS dofsPeriodic() const; // ,, for the case that the periodicity if fully eliminated + xt::xtensor nodesPeriodic() const; // periodic node pairs [:,2]: (independent, dependent) + size_t nodesOrigin() const; // bottom-left node, used as reference for periodicity + xt::xtensor dofs() const; // DOF-numbers for each component of each node (sequential) + xt::xtensor dofsPeriodic() const; // ,, for the case that the periodicity if fully eliminated }; // ----------------------------------------- mesh analysis ----------------------------------------- // read / set the orientation (-1 / +1) of all triangles -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); +inline xt::xtensor getOrientation(const xt::xtensor &coor, const xt::xtensor &conn ); +inline xt::xtensor setOrientation(const xt::xtensor &coor, const xt::xtensor &conn, int orientation=-1); +inline xt::xtensor setOrientation(const xt::xtensor &coor, const xt::xtensor &conn, const xt::xtensor &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" -inline MatS retriangulate(const MatD &coor, const MatS &conn, int orientation=-1); - +inline xt::xtensor retriangulate(const xt::xtensor &coor, const xt::xtensor &conn, int orientation=-1); // ================================= GooseFEM::Mesh::Tri3::Private ================================= namespace Private { // ------------------------- support class - update existing triangulation ------------------------- // 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) + size_t m_nelem; // #elements + size_t m_nnode; // #nodes + size_t m_nne; // #nodes-per-element + size_t m_ndim; // #dimensions + xt::xtensor m_edge; // the element that neighbors along each edge (m_nelem: no neighbor) + xt::xtensor m_conn; // connectivity (updated) + xt::xtensor m_coor; // nodal positions (does not change) + xt::xtensor_fixed> m_elem; // the two elements in the last element change + xt::xtensor_fixed> m_node; // the four nodes in the last element change // 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); + TriUpdate(const xt::xtensor &coor, const xt::xtensor &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 + bool eval (); // re-triangulate the full mesh (returns "true" if changed) + bool increment(); // one re-triangulation step (returns "true" if changed) + xt::xtensor conn() { return m_conn; }; // (new) connectivity + xt::xtensor ch_elem() { return m_elem; }; // element involved in last element change + xt::xtensor ch_node() { return m_node; }; // nodes involved in last element change }; // -------------------------------- support class - edge definition -------------------------------- // 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); }; // ------------------------------------------------------------------------------------------------- // compare edges inline bool Edge_cmp (Edge a, Edge b); // check equality inline bool Edge_sort(Edge a, Edge b); // check if "a" is smaller than "b" in terms of node-numbers // ------------------------------------------------------------------------------------------------- } // namespace Private // ------------------------------------------------------------------------------------------------- }}} // namespace ... #endif diff --git a/src/GooseFEM/MeshTri3.hpp b/src/GooseFEM/MeshTri3.hpp index 1c242ce..185a2ae 100644 --- a/src/GooseFEM/MeshTri3.hpp +++ b/src/GooseFEM/MeshTri3.hpp @@ -1,606 +1,600 @@ /* ================================================================================================= (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM ================================================================================================= */ #ifndef GOOSEFEM_MESHTRI3_CPP #define GOOSEFEM_MESHTRI3_CPP // ------------------------------------------------------------------------------------------------- #include "MeshTri3.h" // ===================================== GooseFEM::Mesh::Tri3 ====================================== namespace GooseFEM { namespace Mesh { namespace Tri3 { // ------------------------------------------ constructor ------------------------------------------ inline Regular::Regular(size_t nelx, size_t nely, double h): m_h(h), m_nelx(nelx), m_nely(nely) { assert( m_nelx >= 1 ); assert( m_nely >= 1 ); - m_nnode = (m_nelx+1) * (m_nely+1) ; - m_nelem = m_nelx * m_nely * 2; + m_nnode = (m_nelx+1) * (m_nely+1); + m_nelem = m_nelx * m_nely * 2; } // -------------------------------------- number of elements --------------------------------------- inline size_t Regular::nelem() const { return m_nelem; } // ---------------------------------------- number of nodes ---------------------------------------- inline size_t Regular::nnode() const { return m_nnode; } // ---------------------------------- number of nodes per element ---------------------------------- inline size_t Regular::nne() const { return m_nne; } // ------------------------------------- number of dimensions -------------------------------------- inline size_t Regular::ndim() const { return m_ndim; } // --------------------------------- coordinates (nodal positions) --------------------------------- -inline MatD Regular::coor() const +inline xt::xtensor Regular::coor() const { - MatD out(m_nnode , m_ndim); + xt::xtensor out = xt::empty({m_nnode, m_ndim}); - ColD x = ColD::LinSpaced( m_nelx+1 , 0.0 , m_h * static_cast(m_nelx) ); - ColD y = ColD::LinSpaced( m_nely+1 , 0.0 , m_h * static_cast(m_nely) ); + xt::xtensor x = xt::linspace(0.0, m_h*static_cast(m_nelx), m_nelx+1); + xt::xtensor y = xt::linspace(0.0, m_h*static_cast(m_nely), m_nely+1); size_t inode = 0; for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) { for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) { out(inode,0) = x(ix); out(inode,1) = y(iy); ++inode; } } return out; } // ---------------------------- connectivity (node-numbers per element) ---------------------------- -inline MatS Regular::conn() const +inline xt::xtensor Regular::conn() const { - MatS out(m_nelem , m_nne); + xt::xtensor out = xt::empty({m_nelem,m_nne}); size_t ielem = 0; for ( size_t iy = 0 ; iy < m_nely ; ++iy ) { for ( size_t ix = 0 ; ix < m_nelx ; ++ix ) { out(ielem,0) = (iy )*(m_nelx+1) + (ix ); out(ielem,1) = (iy )*(m_nelx+1) + (ix+1); out(ielem,2) = (iy+1)*(m_nelx+1) + (ix ); ++ielem; out(ielem,0) = (iy )*(m_nelx+1) + (ix+1); out(ielem,1) = (iy+1)*(m_nelx+1) + (ix+1); out(ielem,2) = (iy+1)*(m_nelx+1) + (ix ); ++ielem; } } return out; } // ------------------------------ node-numbers along the bottom edge ------------------------------- -inline ColS Regular::nodesBottomEdge() const +inline xt::xtensor Regular::nodesBottomEdge() const { - ColS out(m_nelx+1); + xt::xtensor out = xt::empty({m_nelx+1}); for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) out(ix) = ix; return out; } // -------------------------------- node-numbers along the top edge -------------------------------- -inline ColS Regular::nodesTopEdge() const +inline xt::xtensor Regular::nodesTopEdge() const { - ColS out(m_nelx+1); + xt::xtensor out = xt::empty({m_nelx+1}); for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix ) out(ix) = ix + m_nely*(m_nelx+1); return out; } // ------------------------------- node-numbers along the left edge -------------------------------- -inline ColS Regular::nodesLeftEdge() const +inline xt::xtensor Regular::nodesLeftEdge() const { - ColS out(m_nely+1); + xt::xtensor out = xt::empty({m_nely+1}); for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) out(iy) = iy*(m_nelx+1); return out; } // ------------------------------- node-numbers along the right edge ------------------------------- -inline ColS Regular::nodesRightEdge() const +inline xt::xtensor Regular::nodesRightEdge() const { - ColS out(m_nely+1); + xt::xtensor out = xt::empty({m_nely+1}); for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy ) out(iy) = iy*(m_nelx+1) + m_nelx; return out; } // ---------------------- node-numbers along the bottom edge, without corners ---------------------- -inline ColS Regular::nodesBottomOpenEdge() const +inline xt::xtensor Regular::nodesBottomOpenEdge() const { - ColS out(m_nelx-1); + xt::xtensor out = xt::empty({m_nelx-1}); for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) out(ix-1) = ix; return out; } // ----------------------- node-numbers along the top edge, without corners ------------------------ -inline ColS Regular::nodesTopOpenEdge() const +inline xt::xtensor Regular::nodesTopOpenEdge() const { - ColS out(m_nelx-1); + xt::xtensor out = xt::empty({m_nelx-1}); for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) out(ix-1) = ix + m_nely*(m_nelx+1); return out; } // ----------------------- node-numbers along the left edge, without corners ----------------------- -inline ColS Regular::nodesLeftOpenEdge() const +inline xt::xtensor Regular::nodesLeftOpenEdge() const { - ColS out(m_nely-1); + xt::xtensor out = xt::empty({m_nely-1}); for ( size_t iy = 1 ; iy < m_nely ; ++iy ) out(iy-1) = iy*(m_nelx+1); return out; } // ---------------------- node-numbers along the right edge, without corners ----------------------- -inline ColS Regular::nodesRightOpenEdge() const +inline xt::xtensor Regular::nodesRightOpenEdge() const { - ColS out(m_nely-1); + xt::xtensor out = xt::empty({m_nely-1}); for ( size_t iy = 1 ; iy < m_nely ; ++iy ) out(iy-1) = iy*(m_nelx+1) + m_nelx; return out; } // ----------------------------- node-number of the bottom-left corner ----------------------------- inline size_t Regular::nodesBottomLeftCorner() const { return 0; } // ---------------------------- node-number of the bottom-right corner ----------------------------- inline size_t Regular::nodesBottomRightCorner() const { return m_nelx; } // ------------------------------ node-number of the top-left corner ------------------------------- inline size_t Regular::nodesTopLeftCorner() const { return m_nely*(m_nelx+1); } // ------------------------------ node-number of the top-right corner ------------------------------ inline size_t Regular::nodesTopRightCorner() const { return m_nely*(m_nelx+1) + m_nelx; } // ----------------------------- node-number of the corners (aliases) ------------------------------ -inline size_t Regular::nodesLeftBottomCorner() const { return nodesBottomLeftCorner(); } -inline size_t Regular::nodesLeftTopCorner() const { return nodesTopLeftCorner(); } +inline size_t Regular::nodesLeftBottomCorner() const { return nodesBottomLeftCorner(); } +inline size_t Regular::nodesLeftTopCorner() const { return nodesTopLeftCorner(); } inline size_t Regular::nodesRightBottomCorner() const { return nodesBottomRightCorner(); } -inline size_t Regular::nodesRightTopCorner() const { return nodesTopRightCorner(); } +inline size_t Regular::nodesRightTopCorner() const { return nodesTopRightCorner(); } // ------------------------------ node-numbers of periodic node-pairs ------------------------------ -inline MatS Regular::nodesPeriodic() const +inline xt::xtensor Regular::nodesPeriodic() const { // edges (without corners) - ColS bot = nodesBottomOpenEdge(); - ColS top = nodesTopOpenEdge(); - ColS lft = nodesLeftOpenEdge(); - ColS rgt = nodesRightOpenEdge(); + xt::xtensor bot = nodesBottomOpenEdge(); + xt::xtensor top = nodesTopOpenEdge(); + xt::xtensor lft = nodesLeftOpenEdge(); + xt::xtensor rgt = nodesRightOpenEdge(); // allocate nodal ties // - number of tying per category size_t tedge = bot.size() + lft.size(); size_t tnode = 3; // - allocate - MatS out(tedge+tnode, 2); + xt::xtensor out = xt::empty({tedge+tnode, std::size_t(2)}); // counter size_t i = 0; // tie all corners to one corner out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesBottomRightCorner(); ++i; out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopRightCorner(); ++i; out(i,0) = nodesBottomLeftCorner(); out(i,1) = nodesTopLeftCorner(); ++i; // tie all corresponding edges to each other - for ( auto j = 0 ; j Regular::dofs() const { return GooseFEM::Mesh::dofs(m_nnode,m_ndim); } // ------------------------ DOP-numbers with periodic dependencies removed ------------------------- -inline MatS Regular::dofsPeriodic() const +inline xt::xtensor Regular::dofsPeriodic() const { // DOF-numbers for each component of each node (sequential) - MatS out = GooseFEM::Mesh::dofs(m_nnode,m_ndim); + xt::xtensor out = GooseFEM::Mesh::dofs(m_nnode,m_ndim); // periodic node-pairs - MatS nodePer = nodesPeriodic(); - size_t nper = static_cast(nodePer.rows()); + xt::xtensor nodePer = nodesPeriodic(); // eliminate 'dependent' DOFs; renumber "out" to be sequential for the remaining DOFs - for ( size_t i = 0 ; i < nper ; ++i ) + for ( size_t i = 0 ; i < nodePer.shape()[0] ; ++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); } // ------------------------------ get the orientation of each element ------------------------------ -inline ColI getOrientation(const MatD &coor, const MatS &conn) +inline xt::xtensor getOrientation(const xt::xtensor &coor, const xt::xtensor &conn) { - assert( conn.cols() == 3 ); - assert( coor.cols() == 2 ); + assert( conn.shape()[1] == 3 ); + assert( coor.shape()[1] == 2 ); - Eigen::Vector2d v1,v2; double k; - size_t nelem = static_cast( conn.rows() ); + size_t nelem = conn.shape()[0]; - ColI out( nelem ); + xt::xtensor out = xt::empty({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) ); + auto v1 = xt::view(coor, conn(ielem,0), xt::all()) - xt::view(coor, conn(ielem,1), xt::all()); + auto v2 = xt::view(coor, conn(ielem,2), xt::all()) - xt::view(coor, conn(ielem,1), xt::all()); - k = v1(0) * v2(1) - v2(0) * v1(1); + k = v1(0) * v2(1) - v2(0) * v1(1); - if ( k < 0 ) out( ielem ) = -1; - else out( ielem ) = +1; + if ( k < 0 ) out(ielem) = -1; + else out(ielem) = +1; } return out; } // ------------------------------ set the orientation of each element ------------------------------ -inline MatS setOrientation(const MatD &coor, const MatS &conn, int orientation) +inline xt::xtensor setOrientation(const xt::xtensor &coor, const xt::xtensor &conn, int orientation) { - assert( conn.cols() == 3 ); - assert( coor.cols() == 2 ); + assert( conn.shape()[1] == 3 ); + assert( coor.shape()[1] == 2 ); assert( orientation == -1 || orientation == +1 ); - ColI val = getOrientation(coor, conn); + xt::xtensor val = getOrientation(coor, conn); - return setOrientation( coor, conn, val, orientation ); + return setOrientation(coor, conn, val, orientation); } // -------------------- set the orientation of each element to a certain value --------------------- -inline MatS setOrientation(const MatD &coor, const MatS &conn, const ColI &val, int orientation) +inline xt::xtensor setOrientation(const xt::xtensor &coor, const xt::xtensor &conn, const xt::xtensor &val, int orientation) { - assert( conn.cols() == 3 ); - assert( coor.cols() == 2 ); - assert( conn.rows() == val.size() ); + assert( conn.shape()[1] == 3 ); + assert( coor.shape()[1] == 2 ); + assert( conn.shape()[0] == val.size() ); assert( orientation == -1 || orientation == +1 ); // avoid compiler warning UNUSED(coor); - size_t nelem = static_cast(conn.rows()); - MatS out = conn; + size_t nelem = conn.shape()[0]; + xt::xtensor 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 API -------------------------------------- -inline MatS retriangulate(const MatD &coor, const MatS &conn, int orientation) +inline xt::xtensor retriangulate(const xt::xtensor &coor, const xt::xtensor &conn, int orientation) { // get the orientation of all elements - ColI dir = getOrientation(coor, conn); + xt::xtensor dir = getOrientation(coor, conn); // check the orientation - bool eq = std::abs(dir.sum()) == conn.rows(); + bool eq = static_cast(std::abs(xt::sum(dir)[0])) == conn.shape()[0]; // new connectivity - MatS out; + xt::xtensor out; // perform re-triangulation // - use "TriUpdate" if ( eq ) { Private::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); } // ================================= GooseFEM::Mesh::Tri3::Private ================================= namespace Private { // ------------------------------------------ constructor ------------------------------------------ -inline TriUpdate::TriUpdate(const MatD &coor, const MatS &conn): m_conn(conn), m_coor(coor) +inline TriUpdate::TriUpdate(const xt::xtensor &coor, const xt::xtensor &conn): m_conn(conn), m_coor(coor) { - assert( conn.cols() == 3 ); - assert( coor.cols() == 2 ); + assert( conn.shape()[1] == 3 ); + assert( coor.shape()[1] == 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); + m_nnode = coor.shape()[0]; + m_ndim = coor.shape()[1]; + m_nelem = conn.shape()[0]; + m_nne = conn.shape()[1]; // set default to out-of-bounds, to make clear that nothing happened yet - m_elem.setConstant(m_nelem); - m_node.setConstant(m_nnode); + m_elem = m_nelem * xt::ones({2}); + m_node = m_nnode * xt::ones({4}); edge(); } // -------------------------- compute neighbors per edge of all elements --------------------------- inline void TriUpdate::edge() { - m_edge.resize(m_nelem , m_nne); - m_edge.setConstant(m_nelem); // signal that nothing has been set + // signal that nothing has been set + m_edge = m_nelem * xt::ones({m_nelem , m_nne}); 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; } } } // ---------------------------- update edges around renumbered elements ---------------------------- 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; } // --------------------------------- re-triangulate the full mesh ---------------------------------- inline bool TriUpdate::eval() { bool change = false; while ( increment() ) { change = true; } return change; } // ----------------------------------- one re-triangulation step ----------------------------------- inline bool TriUpdate::increment() { size_t ielem,jelem,iedge,jedge; double phi1,phi2; - ColS c(4); - ColS n(4); + xt::xtensor_fixed> c = xt::empty({4}); + xt::xtensor_fixed> n = xt::empty({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; } // ------------------------------------------ constructor ------------------------------------------ inline Edge::Edge(size_t i, size_t j, size_t el, size_t ed, bool sort): n1(i), n2(j), elem(el), edge(ed) { if ( sort && n1>n2 ) std::swap(n1,n2); } // --------------------------------------- compare two edges --------------------------------------- inline bool Edge_cmp(Edge a, Edge b) { if ( a.n1 == b.n1 and a.n2 == b.n2 ) return true; return false; } // ----------------------- sort edges by comparing the first and second node ----------------------- inline bool Edge_sort(Edge a, Edge b) { if ( a.n1 < b.n1 or a.n2 < b.n2 ) return true; return false; } // ------------------------------------------------------------------------------------------------- } // namespace Private // ------------------------------------------------------------------------------------------------- }}} // namespace ... // ================================================================================================= #endif