Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F107212719
MeshHex8.hpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Apr 6, 00:35
Size
95 KB
Mime Type
text/x-c
Expires
Tue, Apr 8, 00:35 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
25344857
Attached To
rGOOSEFEM GooseFEM
MeshHex8.hpp
View Options
/* =================================================================================================
(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
================================================================================================= */
#ifndef XGOOSEFEM_MESHHEX8_CPP
#define XGOOSEFEM_MESHHEX8_CPP
// -------------------------------------------------------------------------------------------------
#include "MeshHex8.h"
// ===================================== GooseFEM::Mesh::Hex8 ======================================
namespace xGooseFEM {
namespace Mesh {
namespace Hex8 {
// ------------------------------------------ constructor ------------------------------------------
inline Regular::Regular(size_t nelx, size_t nely, size_t nelz, double h):
m_h(h), m_nelx(nelx), m_nely(nely), m_nelz(nelz)
{
assert( m_nelx >= 1 );
assert( m_nely >= 1 );
assert( m_nelz >= 1 );
m_nnode = (m_nelx+1) * (m_nely+1) * (m_nelz+1);
m_nelem = m_nelx * m_nely * m_nelz ;
}
// -------------------------------------- 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 xt::xtensor<double,2> Regular::coor() const
{
xt::xtensor<double,2> out = xt::empty<double>({m_nnode, m_ndim});
xt::xtensor<double,1> x = xt::linspace<double>(0.0, m_h*static_cast<double>(m_nelx), m_nelx+1);
xt::xtensor<double,1> y = xt::linspace<double>(0.0, m_h*static_cast<double>(m_nely), m_nely+1);
xt::xtensor<double,1> z = xt::linspace<double>(0.0, m_h*static_cast<double>(m_nelz), m_nelz+1);
size_t inode = 0;
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz ) {
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);
out(inode,2) = z(iz);
++inode;
}
}
}
return out;
}
// ---------------------------- connectivity (node-numbers per element) ----------------------------
inline xt::xtensor<size_t,2> Regular::conn() const
{
xt::xtensor<size_t,2> out = xt::empty<size_t>({m_nelem,m_nne});
size_t ielem = 0;
for ( size_t iz = 0 ; iz < m_nelz ; ++iz ) {
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 ) + (iz )*(m_nely+1)*(m_nelx+1);
out(ielem,1) = (iy )*(m_nelx+1) + (ix+1) + (iz )*(m_nely+1)*(m_nelx+1);
out(ielem,3) = (iy+1)*(m_nelx+1) + (ix ) + (iz )*(m_nely+1)*(m_nelx+1);
out(ielem,2) = (iy+1)*(m_nelx+1) + (ix+1) + (iz )*(m_nely+1)*(m_nelx+1);
out(ielem,4) = (iy )*(m_nelx+1) + (ix ) + (iz+1)*(m_nely+1)*(m_nelx+1);
out(ielem,5) = (iy )*(m_nelx+1) + (ix+1) + (iz+1)*(m_nely+1)*(m_nelx+1);
out(ielem,7) = (iy+1)*(m_nelx+1) + (ix ) + (iz+1)*(m_nely+1)*(m_nelx+1);
out(ielem,6) = (iy+1)*(m_nelx+1) + (ix+1) + (iz+1)*(m_nely+1)*(m_nelx+1);
++ielem;
}
}
}
return out;
}
// ------------------------------ node-numbers along the front plane -------------------------------
inline xt::xtensor<size_t,1> Regular::nodesFront() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nelx+1)*(m_nely+1)});
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(iy*(m_nelx+1)+ix) = iy*(m_nelx+1) + ix;
return out;
}
// ------------------------------- node-numbers along the back plane -------------------------------
inline xt::xtensor<size_t,1> Regular::nodesBack() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nelx+1)*(m_nely+1)});
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(iy*(m_nelx+1)+ix) = iy*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1);
return out;
}
// ------------------------------- node-numbers along the left plane -------------------------------
inline xt::xtensor<size_t,1> Regular::nodesLeft() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nely+1)*(m_nelz+1)});
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iz*(m_nely+1)+iy) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ------------------------------ node-numbers along the right plane -------------------------------
inline xt::xtensor<size_t,1> Regular::nodesRight() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nely+1)*(m_nelz+1)});
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iz*(m_nely+1)+iy) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// ------------------------------ node-numbers along the bottom plane ------------------------------
inline xt::xtensor<size_t,1> Regular::nodesBottom() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nelx+1)*(m_nelz+1)});
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(iz*(m_nelx+1)+ix) = ix + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ------------------------------- node-numbers along the top plane --------------------------------
inline xt::xtensor<size_t,1> Regular::nodesTop() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nelx+1)*(m_nelz+1)});
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(iz*(m_nelx+1)+ix) = ix + m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ------------------------------- node-numbers along the front face -------------------------------
inline xt::xtensor<size_t,1> Regular::nodesFrontFace() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nelx-1)*(m_nely-1)});
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out((iy-1)*(m_nelx-1)+(ix-1)) = iy*(m_nelx+1) + ix;
return out;
}
// ------------------------------- node-numbers along the back face --------------------------------
inline xt::xtensor<size_t,1> Regular::nodesBackFace() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nelx-1)*(m_nely-1)});
for ( size_t iy = 1 ; iy < m_nely ; ++iy ) {
for ( size_t ix = 1 ; ix < m_nelx ; ++ix ) {
out((iy-1)*(m_nelx-1)+(ix-1)) = iy*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1);
}
}
return out;
}
// ------------------------------- node-numbers along the left face --------------------------------
inline xt::xtensor<size_t,1> Regular::nodesLeftFace() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nely-1)*(m_nelz-1)});
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out((iz-1)*(m_nely-1)+(iy-1)) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ------------------------------- node-numbers along the right face -------------------------------
inline xt::xtensor<size_t,1> Regular::nodesRightFace() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nely-1)*(m_nelz-1)});
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out((iz-1)*(m_nely-1)+(iy-1)) = iy*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// ------------------------------ node-numbers along the bottom face -------------------------------
inline xt::xtensor<size_t,1> Regular::nodesBottomFace() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nelx-1)*(m_nelz-1)});
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out((iz-1)*(m_nelx-1)+(ix-1)) = ix + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// -------------------------------- node-numbers along the top face --------------------------------
inline xt::xtensor<size_t,1> Regular::nodesTopFace() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nelx-1)*(m_nelz-1)});
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out((iz-1)*(m_nelx-1)+(ix-1)) = ix + m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// --------------------------- node-numbers along the front-bottom edge ----------------------------
inline xt::xtensor<size_t,1> Regular::nodesFrontBottomEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelx+1});
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(ix) = ix;
return out;
}
// ----------------------------- node-numbers along the front-top edge -----------------------------
inline xt::xtensor<size_t,1> Regular::nodesFrontTopEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 front-left edge -----------------------------
inline xt::xtensor<size_t,1> Regular::nodesFrontLeftEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 front-right edge ----------------------------
inline xt::xtensor<size_t,1> Regular::nodesFrontRightEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 back-bottom edge ----------------------------
inline xt::xtensor<size_t,1> Regular::nodesBackBottomEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelx+1});
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(ix) = ix + m_nelz*(m_nely+1)*(m_nelx+1);
return out;
}
// ----------------------------- node-numbers along the back-top edge ------------------------------
inline xt::xtensor<size_t,1> Regular::nodesBackTopEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelx+1});
for ( size_t ix = 0 ; ix < m_nelx+1 ; ++ix )
out(ix) = m_nely*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1);
return out;
}
// ----------------------------- node-numbers along the back-left edge -----------------------------
inline xt::xtensor<size_t,1> Regular::nodesBackLeftEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nely+1});
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iy) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1);
return out;
}
// ---------------------------- node-numbers along the back-right edge -----------------------------
inline xt::xtensor<size_t,1> Regular::nodesBackRightEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nely+1});
for ( size_t iy = 0 ; iy < m_nely+1 ; ++iy )
out(iy) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// ---------------------------- node-numbers along the bottom-left edge ----------------------------
inline xt::xtensor<size_t,1> Regular::nodesBottomLeftEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz+1});
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
out(iz) = iz*(m_nelx+1)*(m_nely+1);
return out;
}
// --------------------------- node-numbers along the bottom-right edge ----------------------------
inline xt::xtensor<size_t,1> Regular::nodesBottomRightEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz+1});
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
out(iz) = iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// ----------------------------- node-numbers along the top-left edge ------------------------------
inline xt::xtensor<size_t,1> Regular::nodesTopLeftEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz+1});
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
out(iz) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ----------------------------- node-numbers along the top-right edge -----------------------------
inline xt::xtensor<size_t,1> Regular::nodesTopRightEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz+1});
for ( size_t iz = 0 ; iz < m_nelz+1 ; ++iz )
out(iz) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// -------------------------------------------- aliases --------------------------------------------
inline xt::xtensor<size_t,1> Regular::nodesBottomFrontEdge() const { return nodesFrontBottomEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesBottomBackEdge() const { return nodesBackBottomEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesTopFrontEdge() const { return nodesFrontTopEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesTopBackEdge() const { return nodesBackTopEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesLeftBottomEdge() const { return nodesBottomLeftEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesLeftFrontEdge() const { return nodesFrontLeftEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesLeftBackEdge() const { return nodesBackLeftEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesLeftTopEdge() const { return nodesTopLeftEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesRightBottomEdge() const { return nodesBottomRightEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesRightTopEdge() const { return nodesTopRightEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesRightFrontEdge() const { return nodesFrontRightEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesRightBackEdge() const { return nodesBackRightEdge(); }
// ------------------- node-numbers along the front-bottom edge, without corners -------------------
inline xt::xtensor<size_t,1> Regular::nodesFrontBottomOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelx-1});
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out(ix-1) = ix;
return out;
}
// -------------------- node-numbers along the front-top edge, without corners ---------------------
inline xt::xtensor<size_t,1> Regular::nodesFrontTopOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 front-left edge, without corners --------------------
inline xt::xtensor<size_t,1> Regular::nodesFrontLeftOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 front-right edge, without corners --------------------
inline xt::xtensor<size_t,1> Regular::nodesFrontRightOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nely-1});
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out(iy-1) = iy*(m_nelx+1) + m_nelx;
return out;
}
// ------------------- node-numbers along the back-bottom edge, without corners --------------------
inline xt::xtensor<size_t,1> Regular::nodesBackBottomOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelx-1});
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out(ix-1) = ix + m_nelz*(m_nely+1)*(m_nelx+1);
return out;
}
// --------------------- node-numbers along the back-top edge, without corners ---------------------
inline xt::xtensor<size_t,1> Regular::nodesBackTopOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelx-1});
for ( size_t ix = 1 ; ix < m_nelx ; ++ix )
out(ix-1) = m_nely*(m_nelx+1) + ix + m_nelz*(m_nely+1)*(m_nelx+1);
return out;
}
// -------------------- node-numbers along the back-left edge, without corners ---------------------
inline xt::xtensor<size_t,1> Regular::nodesBackLeftOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nely-1});
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out(iy-1) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1);
return out;
}
// -------------------- node-numbers along the back-right edge, without corners --------------------
inline xt::xtensor<size_t,1> Regular::nodesBackRightOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nely-1});
for ( size_t iy = 1 ; iy < m_nely ; ++iy )
out(iy-1) = iy*(m_nelx+1) + m_nelz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// ------------------- node-numbers along the bottom-left edge, without corners --------------------
inline xt::xtensor<size_t,1> Regular::nodesBottomLeftOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz-1});
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
out(iz-1) = iz*(m_nelx+1)*(m_nely+1);
return out;
}
// ------------------- node-numbers along the bottom-right edge, without corners -------------------
inline xt::xtensor<size_t,1> Regular::nodesBottomRightOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz-1});
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
out(iz-1) = iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// --------------------- node-numbers along the top-left edge, without corners ---------------------
inline xt::xtensor<size_t,1> Regular::nodesTopLeftOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz-1});
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
out(iz-1) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1);
return out;
}
// -------------------- node-numbers along the top-right edge, without corners ---------------------
inline xt::xtensor<size_t,1> Regular::nodesTopRightOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz-1});
for ( size_t iz = 1 ; iz < m_nelz ; ++iz )
out(iz-1) = m_nely*(m_nelx+1) + iz*(m_nelx+1)*(m_nely+1) + m_nelx;
return out;
}
// -------------------------------------------- aliases --------------------------------------------
inline xt::xtensor<size_t,1> Regular::nodesBottomFrontOpenEdge() const { return nodesFrontBottomOpenEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesBottomBackOpenEdge() const { return nodesBackBottomOpenEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesTopFrontOpenEdge() const { return nodesFrontTopOpenEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesTopBackOpenEdge() const { return nodesBackTopOpenEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesLeftBottomOpenEdge() const { return nodesBottomLeftOpenEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesLeftFrontOpenEdge() const { return nodesFrontLeftOpenEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesLeftBackOpenEdge() const { return nodesBackLeftOpenEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesLeftTopOpenEdge() const { return nodesTopLeftOpenEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesRightBottomOpenEdge() const { return nodesBottomRightOpenEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesRightTopOpenEdge() const { return nodesTopRightOpenEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesRightFrontOpenEdge() const { return nodesFrontRightOpenEdge(); }
inline xt::xtensor<size_t,1> Regular::nodesRightBackOpenEdge() const { return nodesBackRightOpenEdge(); }
// -------------------------- node-number of the front-bottom-left corner --------------------------
inline size_t Regular::nodesFrontBottomLeftCorner() const
{
return 0;
}
// ------------------------- node-number of the front-bottom-right corner --------------------------
inline size_t Regular::nodesFrontBottomRightCorner() const
{
return m_nelx;
}
// --------------------------- node-number of the front-top-left corner ----------------------------
inline size_t Regular::nodesFrontTopLeftCorner() const
{
return m_nely*(m_nelx+1);
}
// --------------------------- node-number of the front-top-right corner ---------------------------
inline size_t Regular::nodesFrontTopRightCorner() const
{
return m_nely*(m_nelx+1) + m_nelx;
}
// -------------------------- node-number of the back-bottom-left corner ---------------------------
inline size_t Regular::nodesBackBottomLeftCorner() const
{
return m_nelz*(m_nely+1)*(m_nelx+1);
}
// -------------------------- node-number of the back-bottom-right corner --------------------------
inline size_t Regular::nodesBackBottomRightCorner() const
{
return m_nelx + m_nelz*(m_nely+1)*(m_nelx+1);
}
// ---------------------------- node-number of the back-top-left corner ----------------------------
inline size_t Regular::nodesBackTopLeftCorner() const
{
return m_nely*(m_nelx+1) + m_nelz*(m_nely+1)*(m_nelx+1);
}
// --------------------------- node-number of the back-top-right corner ----------------------------
inline size_t Regular::nodesBackTopRightCorner() const
{
return m_nely*(m_nelx+1) + m_nelx + m_nelz*(m_nely+1)*(m_nelx+1);
}
// -------------------------------------------- aliases --------------------------------------------
inline size_t Regular::nodesFrontLeftBottomCorner() const { return nodesFrontBottomLeftCorner(); }
inline size_t Regular::nodesBottomFrontLeftCorner() const { return nodesFrontBottomLeftCorner(); }
inline size_t Regular::nodesBottomLeftFrontCorner() const { return nodesFrontBottomLeftCorner(); }
inline size_t Regular::nodesLeftFrontBottomCorner() const { return nodesFrontBottomLeftCorner(); }
inline size_t Regular::nodesLeftBottomFrontCorner() const { return nodesFrontBottomLeftCorner(); }
inline size_t Regular::nodesFrontRightBottomCorner()const { return nodesFrontBottomRightCorner(); }
inline size_t Regular::nodesBottomFrontRightCorner()const { return nodesFrontBottomRightCorner(); }
inline size_t Regular::nodesBottomRightFrontCorner()const { return nodesFrontBottomRightCorner(); }
inline size_t Regular::nodesRightFrontBottomCorner()const { return nodesFrontBottomRightCorner(); }
inline size_t Regular::nodesRightBottomFrontCorner()const { return nodesFrontBottomRightCorner(); }
inline size_t Regular::nodesFrontLeftTopCorner() const { return nodesFrontTopLeftCorner(); }
inline size_t Regular::nodesTopFrontLeftCorner() const { return nodesFrontTopLeftCorner(); }
inline size_t Regular::nodesTopLeftFrontCorner() const { return nodesFrontTopLeftCorner(); }
inline size_t Regular::nodesLeftFrontTopCorner() const { return nodesFrontTopLeftCorner(); }
inline size_t Regular::nodesLeftTopFrontCorner() const { return nodesFrontTopLeftCorner(); }
inline size_t Regular::nodesFrontRightTopCorner() const { return nodesFrontTopRightCorner(); }
inline size_t Regular::nodesTopFrontRightCorner() const { return nodesFrontTopRightCorner(); }
inline size_t Regular::nodesTopRightFrontCorner() const { return nodesFrontTopRightCorner(); }
inline size_t Regular::nodesRightFrontTopCorner() const { return nodesFrontTopRightCorner(); }
inline size_t Regular::nodesRightTopFrontCorner() const { return nodesFrontTopRightCorner(); }
inline size_t Regular::nodesBackLeftBottomCorner() const { return nodesBackBottomLeftCorner(); }
inline size_t Regular::nodesBottomBackLeftCorner() const { return nodesBackBottomLeftCorner(); }
inline size_t Regular::nodesBottomLeftBackCorner() const { return nodesBackBottomLeftCorner(); }
inline size_t Regular::nodesLeftBackBottomCorner() const { return nodesBackBottomLeftCorner(); }
inline size_t Regular::nodesLeftBottomBackCorner() const { return nodesBackBottomLeftCorner(); }
inline size_t Regular::nodesBackRightBottomCorner() const { return nodesBackBottomRightCorner(); }
inline size_t Regular::nodesBottomBackRightCorner() const { return nodesBackBottomRightCorner(); }
inline size_t Regular::nodesBottomRightBackCorner() const { return nodesBackBottomRightCorner(); }
inline size_t Regular::nodesRightBackBottomCorner() const { return nodesBackBottomRightCorner(); }
inline size_t Regular::nodesRightBottomBackCorner() const { return nodesBackBottomRightCorner(); }
inline size_t Regular::nodesBackLeftTopCorner() const { return nodesBackTopLeftCorner(); }
inline size_t Regular::nodesTopBackLeftCorner() const { return nodesBackTopLeftCorner(); }
inline size_t Regular::nodesTopLeftBackCorner() const { return nodesBackTopLeftCorner(); }
inline size_t Regular::nodesLeftBackTopCorner() const { return nodesBackTopLeftCorner(); }
inline size_t Regular::nodesLeftTopBackCorner() const { return nodesBackTopLeftCorner(); }
inline size_t Regular::nodesBackRightTopCorner() const { return nodesBackTopRightCorner(); }
inline size_t Regular::nodesTopBackRightCorner() const { return nodesBackTopRightCorner(); }
inline size_t Regular::nodesTopRightBackCorner() const { return nodesBackTopRightCorner(); }
inline size_t Regular::nodesRightBackTopCorner() const { return nodesBackTopRightCorner(); }
inline size_t Regular::nodesRightTopBackCorner() const { return nodesBackTopRightCorner(); }
// ------------------------------ node-numbers of periodic node-pairs ------------------------------
inline xt::xtensor<size_t,2> Regular::nodesPeriodic() const
{
// faces
xt::xtensor<size_t,1> fro = nodesFrontFace();
xt::xtensor<size_t,1> bck = nodesBackFace();
xt::xtensor<size_t,1> lft = nodesLeftFace();
xt::xtensor<size_t,1> rgt = nodesRightFace();
xt::xtensor<size_t,1> bot = nodesBottomFace();
xt::xtensor<size_t,1> top = nodesTopFace();
// edges
xt::xtensor<size_t,1> froBot = nodesFrontBottomOpenEdge();
xt::xtensor<size_t,1> froTop = nodesFrontTopOpenEdge();
xt::xtensor<size_t,1> froLft = nodesFrontLeftOpenEdge();
xt::xtensor<size_t,1> froRgt = nodesFrontRightOpenEdge();
xt::xtensor<size_t,1> bckBot = nodesBackBottomOpenEdge();
xt::xtensor<size_t,1> bckTop = nodesBackTopOpenEdge();
xt::xtensor<size_t,1> bckLft = nodesBackLeftOpenEdge();
xt::xtensor<size_t,1> bckRgt = nodesBackRightOpenEdge();
xt::xtensor<size_t,1> botLft = nodesBottomLeftOpenEdge();
xt::xtensor<size_t,1> botRgt = nodesBottomRightOpenEdge();
xt::xtensor<size_t,1> topLft = nodesTopLeftOpenEdge();
xt::xtensor<size_t,1> topRgt = nodesTopRightOpenEdge();
// allocate nodal ties
// - number of tying per category
size_t tface = fro.size() + lft.size() + bot.size();
size_t tedge = 3*froBot.size() + 3*froLft.size() + 3*botLft.size();
size_t tnode = 7;
// - allocate
xt::xtensor<size_t,2> out = xt::empty<size_t>({tface+tedge+tnode, std::size_t(2)});
// counter
size_t i = 0;
// tie all corners to one corner
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontBottomRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomLeftCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopLeftCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopLeftCorner(); ++i;
// tie all corresponding edges to each other (exclude corners)
for ( size_t j = 0 ; j<froBot.size() ; ++j ){ out(i,0) = froBot(j); out(i,1) = bckBot(j); ++i; }
for ( size_t j = 0 ; j<froBot.size() ; ++j ){ out(i,0) = froBot(j); out(i,1) = bckTop(j); ++i; }
for ( size_t j = 0 ; j<froBot.size() ; ++j ){ out(i,0) = froBot(j); out(i,1) = froTop(j); ++i; }
for ( size_t j = 0 ; j<botLft.size() ; ++j ){ out(i,0) = botLft(j); out(i,1) = botRgt(j); ++i; }
for ( size_t j = 0 ; j<botLft.size() ; ++j ){ out(i,0) = botLft(j); out(i,1) = topRgt(j); ++i; }
for ( size_t j = 0 ; j<botLft.size() ; ++j ){ out(i,0) = botLft(j); out(i,1) = topLft(j); ++i; }
for ( size_t j = 0 ; j<froLft.size() ; ++j ){ out(i,0) = froLft(j); out(i,1) = froRgt(j); ++i; }
for ( size_t j = 0 ; j<froLft.size() ; ++j ){ out(i,0) = froLft(j); out(i,1) = bckRgt(j); ++i; }
for ( size_t j = 0 ; j<froLft.size() ; ++j ){ out(i,0) = froLft(j); out(i,1) = bckLft(j); ++i; }
// tie faces to each-other
for ( size_t j = 0 ; j<fro.size() ; ++j ){ out(i,0) = fro(j); out(i,1) = bck(j); ++i; }
for ( size_t j = 0 ; j<lft.size() ; ++j ){ out(i,0) = lft(j); out(i,1) = rgt(j); ++i; }
for ( size_t j = 0 ; j<bot.size() ; ++j ){ out(i,0) = bot(j); out(i,1) = top(j); ++i; }
return out;
}
// ------------------------------ node-number that lies in the origin ------------------------------
inline size_t Regular::nodesOrigin() const
{
return nodesFrontBottomLeftCorner();
}
// ------------------------- DOF numbers per node (sequentially numbered) --------------------------
inline xt::xtensor<size_t,2> Regular::dofs() const
{
return xGooseFEM::Mesh::dofs(m_nnode,m_ndim);
}
// ------------------------ DOP-numbers with periodic dependencies removed -------------------------
inline xt::xtensor<size_t,2> Regular::dofsPeriodic() const
{
// DOF-numbers for each component of each node (sequential)
xt::xtensor<size_t,2> out = xGooseFEM::Mesh::dofs(m_nnode,m_ndim);
// periodic node-pairs
xt::xtensor<size_t,2> 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 xGooseFEM::Mesh::renumber(out);
}
// ------------------------------------------ constructor ------------------------------------------
inline FineLayer::FineLayer(size_t nelx, size_t nely, size_t nelz, double h, size_t nfine):
m_h(h)
{
// basic assumptions
assert( nelx >= 1 );
assert( nely >= 1 );
assert( nelz >= 1 );
// store basic info
m_Lx = m_h * static_cast<double>(nelx);
m_Lz = m_h * static_cast<double>(nelz);
// compute element size in y-direction (use symmetry, compute upper half)
// ----------------------------------------------------------------------
// temporary variables
size_t nmin, ntot;
xt::xtensor<size_t,1> nhx = xt::ones<size_t>({nely});
xt::xtensor<size_t,1> nhy = xt::ones<size_t>({nely});
xt::xtensor<size_t,1> nhz = xt::ones<size_t>({nely});
xt::xtensor<int ,1> refine = -1 * xt::ones<int> ({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-(z-)direction should fit the total number of elements in x-(z-)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 (and z-direction)
if ( 3*nhy(iy) <= ntot and nelx%(3*nhx(iy)) == 0 and ntot+nhy(iy) < nmin )
{
// - process refinement in x-direction
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;
// - rule (2) satisfied: coarsen next element layer in z-direction
if ( iy+1 < nely and ntot+2*nhy(iy) < nmin )
{
if ( nelz%(3*nhz(iy+1)) == 0 )
{
// - update the number of elements in y-direction
ntot += nhy(iy);
// - proceed to next element layer in y-direction
++iy;
// - process refinement in z-direction
refine(iy) = 2;
nhy (iy) = nhy(iy-1);
auto vnhz = xt::view(nhz, xt::range(iy, _));
vnhz *= 3;
}
}
}
// rules (1,2) satisfied: coarse in z-direction
else if ( 3*nhy(iy) <= ntot and nelz%(3*nhz(iy)) == 0 and ntot+nhy(iy) < nmin )
{
// - process refinement in z-direction
refine(iy) = 2;
nhy (iy) *= 2;
auto vnhy = xt::view(nhy, xt::range(iy+1, _));
auto vnhz = xt::view(nhz, xt::range(iy , _));
vnhy *= 3;
vnhz *= 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<size_t>({nely*2-1});
m_nhy = xt::empty<size_t>({nely*2-1});
m_nhz = xt::empty<size_t>({nely*2-1});
m_refine = xt::empty<int> ({nely*2-1});
m_nelx = xt::empty<size_t>({nely*2-1});
m_nelz = xt::empty<size_t>({nely*2-1});
m_nnd = xt::empty<size_t>({nely*2 });
m_startElem = xt::empty<size_t>({nely*2-1});
m_startNode = xt::empty<size_t>({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_nhz (iy) = nhz (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_nhz (iy+nely) = nhz (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);
m_nelz(iy) = nelz / m_nhz(iy);
}
// compute the number of nodes per node layer in y-direction
// - bottom half
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
m_nnd(iy) = (m_nelx(iy)+1) * (m_nelz(iy)+1);
// - top half
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
m_nnd(iy+1) = (m_nelx(iy)+1) * (m_nelz(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) * ( m_nelz(i)+1); }
else if ( m_refine(i) == 2 ) { m_nnode += ( m_nelx(i)+1) * (3*m_nelz(i)+1); }
else { m_nnode += ( m_nelx(i)+1) * ( m_nelz(i)+1); }
// - add the elements of this layer
if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ) * ( m_nelz(i) ); }
else if ( m_refine(i) == 2 ) { m_nelem += ( m_nelx(i) ) * (4*m_nelz(i) ); }
else { m_nelem += ( m_nelx(i) ) * ( m_nelz(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) * ( m_nelz(i)+1); }
else if ( m_refine(i) == 2 ) { m_nnode += ( m_nelx(i)+1) * (5*m_nelz(i)+1); }
else { m_nnode += ( m_nelx(i)+1) * ( m_nelz(i)+1); }
// - add the elements of this layer
if ( m_refine(i) == 0 ) { m_nelem += (4*m_nelx(i) ) * ( m_nelz(i) ); }
else if ( m_refine(i) == 2 ) { m_nelem += ( m_nelx(i) ) * (4*m_nelz(i) ); }
else { m_nelem += ( m_nelx(i) ) * ( m_nelz(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) * (m_nelz(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 <= 2 );
if ( i == 0 ) return xt::amax(m_nelx)[0];
else if ( i == 2 ) return xt::amax(m_nelz)[0];
else return xt::sum (m_nhy )[0];
}
// --------------------------------- coordinates (nodal positions) ---------------------------------
inline xt::xtensor<double,2> FineLayer::coor() const
{
// allocate output
xt::xtensor<double,2> out = xt::empty<double>({m_nnode, m_ndim});
// current node, number of element layers
size_t inode = 0;
size_t nely = static_cast<size_t>(m_nhy.size());
// y-position of each main node layer (i.e. excluding node layers for refinement/coarsening)
// - allocate
xt::xtensor<double,1> y = xt::empty<double>({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<double,1> x = xt::linspace<double>(0.0, m_Lx, m_nelx(iy)+1);
xt::xtensor<double,1> z = xt::linspace<double>(0.0, m_Lz, m_nelz(iy)+1);
// add nodes of the bottom layer of this element
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(inode,0) = x(ix);
out(inode,1) = y(iy);
out(inode,2) = z(iz);
++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<double>(m_nhx(iy)/3);
double dy = m_h * static_cast<double>(m_nhy(iy)/2);
// - add nodes of the intermediate layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
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<double>(j+1);
out(inode,1) = y(iy) + dy;
out(inode,2) = z(iz);
++inode;
}
}
}
}
// add extra nodes of the intermediate layer, for refinement in z-direction
else if ( m_refine(iy) == 2 )
{
// - get position offset in y- and z-direction
double dz = m_h * static_cast<double>(m_nhz(iy)/3);
double dy = m_h * static_cast<double>(m_nhy(iy)/2);
// - add nodes of the intermediate layer
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t j = 0 ; j < 2 ; ++j ) {
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(inode,0) = x(ix);
out(inode,1) = y(iy) + dy;
out(inode,2) = z(iz) + dz * static_cast<double>(j+1);
++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<double,1> x = xt::linspace<double>(0.0, m_Lx, m_nelx(iy)+1);
xt::xtensor<double,1> z = xt::linspace<double>(0.0, m_Lz, m_nelz(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<double>(m_nhx(iy)/3);
double dy = m_h * static_cast<double>(m_nhy(iy)/2);
// - add nodes of the intermediate layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
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<double>(j+1);
out(inode,1) = y(iy) + dy;
out(inode,2) = z(iz);
++inode;
}
}
}
}
// add extra nodes of the intermediate layer, for refinement in z-direction
else if ( m_refine(iy) == 2 )
{
// - get position offset in y- and z-direction
double dz = m_h * static_cast<double>(m_nhz(iy)/3);
double dy = m_h * static_cast<double>(m_nhy(iy)/2);
// - add nodes of the intermediate layer
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t j = 0 ; j < 2 ; ++j ) {
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(inode,0) = x(ix);
out(inode,1) = y(iy) + dy;
out(inode,2) = z(iz) + dz * static_cast<double>(j+1);
++inode;
}
}
}
}
// add nodes of the top layer of this element
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(inode,0) = x(ix );
out(inode,1) = y(iy+1);
out(inode,2) = z(iz );
++inode;
}
}
}
return out;
}
// ---------------------------- connectivity (node-numbers per element) ----------------------------
inline xt::xtensor<size_t,2> FineLayer::conn() const
{
// allocate output
xt::xtensor<size_t,2> out = xt::empty<size_t>({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<size_t>(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 iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
out(ielem,0) = bot + (ix ) + (iz ) * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix+1) + (iz ) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + (iz ) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + (iz ) * (m_nelx(iy)+1);
out(ielem,4) = bot + (ix ) + (iz+1) * (m_nelx(iy)+1);
out(ielem,5) = bot + (ix+1) + (iz+1) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + (iz+1) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + (iz+1) * (m_nelx(iy)+1);
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 iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
// -- bottom element
out(ielem,0) = bot + ( ix ) + (iz ) * ( m_nelx(iy)+1);
out(ielem,1) = bot + ( ix+1) + (iz ) * ( m_nelx(iy)+1);
out(ielem,2) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,3) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,4) = bot + ( ix ) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,5) = bot + ( ix+1) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,6) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
out(ielem,7) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
ielem++;
// -- top-right element
out(ielem,0) = bot + ( ix+1) + (iz ) * ( m_nelx(iy)+1);
out(ielem,1) = top + (3*ix+3) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,2) = top + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,3) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,4) = bot + ( ix+1) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,5) = top + (3*ix+3) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,6) = top + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,7) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
ielem++;
// -- top-center element
out(ielem,0) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,1) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,2) = top + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,3) = top + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,4) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
out(ielem,5) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
out(ielem,6) = top + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,7) = top + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1);
ielem++;
// -- top-left element
out(ielem,0) = bot + ( ix ) + (iz ) * ( m_nelx(iy)+1);
out(ielem,1) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,2) = top + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,3) = top + (3*ix ) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,4) = bot + ( ix ) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,5) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
out(ielem,6) = top + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,7) = top + (3*ix ) + (iz+1) * (3*m_nelx(iy)+1);
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 iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
// -- lower-left element
out(ielem,0) = bot + (3*ix ) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,1) = bot + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,2) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,3) = top + ( ix ) + (iz ) * ( m_nelx(iy)+1);
out(ielem,4) = bot + (3*ix ) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,5) = bot + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,6) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
out(ielem,7) = top + ( ix ) + (iz+1) * ( m_nelx(iy)+1);
ielem++;
// -- lower-center element
out(ielem,0) = bot + (3*ix+1) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,1) = bot + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,2) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,3) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,4) = bot + (3*ix+1) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,5) = bot + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,6) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
out(ielem,7) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
ielem++;
// -- lower-right element
out(ielem,0) = bot + (3*ix+2) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,1) = bot + (3*ix+3) + (iz ) * (3*m_nelx(iy)+1);
out(ielem,2) = top + ( ix+1) + (iz ) * ( m_nelx(iy)+1);
out(ielem,3) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,4) = bot + (3*ix+2) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,5) = bot + (3*ix+3) + (iz+1) * (3*m_nelx(iy)+1);
out(ielem,6) = top + ( ix+1) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,7) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
ielem++;
// -- upper element
out(ielem,0) = mid + (2*ix ) + (iz ) * (2*m_nelx(iy) );
out(ielem,1) = mid + (2*ix+1) + (iz ) * (2*m_nelx(iy) );
out(ielem,2) = top + ( ix+1) + (iz ) * ( m_nelx(iy)+1);
out(ielem,3) = top + ( ix ) + (iz ) * ( m_nelx(iy)+1);
out(ielem,4) = mid + (2*ix ) + (iz+1) * (2*m_nelx(iy) );
out(ielem,5) = mid + (2*ix+1) + (iz+1) * (2*m_nelx(iy) );
out(ielem,6) = top + ( ix+1) + (iz+1) * ( m_nelx(iy)+1);
out(ielem,7) = top + ( ix ) + (iz+1) * ( m_nelx(iy)+1);
ielem++;
}
}
}
// - define connectivity: refinement along the z-direction (below the middle layer)
else if ( m_refine(iy) == 2 and iy <= (nely-1)/2 )
{
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
// -- bottom element
out(ielem,0) = bot + (ix ) + iz * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix ) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,2) = bot + (ix+1) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,3) = bot + (ix+1) + iz * (m_nelx(iy)+1);
out(ielem,4) = mid + (ix ) + 2*iz * (m_nelx(iy)+1);
out(ielem,5) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,6) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,7) = mid + (ix+1) + 2*iz * (m_nelx(iy)+1);
ielem++;
// -- top-back element
out(ielem,0) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,1) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,4) = bot + (ix ) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,5) = bot + (ix+1) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + (3*iz+3) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + (3*iz+3) * (m_nelx(iy)+1);
ielem++;
// -- top-center element
out(ielem,0) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,1) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,4) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,5) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + (3*iz+2) * (m_nelx(iy)+1);
ielem++;
// -- top-front element
out(ielem,0) = bot + (ix ) + ( iz ) * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix+1) + ( iz ) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + (3*iz ) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + (3*iz ) * (m_nelx(iy)+1);
out(ielem,4) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,5) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + (3*iz+1) * (m_nelx(iy)+1);
ielem++;
}
}
}
// - define connectivity: coarsening along the z-direction (above the middle layer)
else if ( m_refine(iy) == 2 and iy > (nely-1)/2 )
{
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz ) {
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix ) {
// -- bottom-front element
out(ielem,0) = bot + (ix ) + (3*iz ) * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix+1) + (3*iz ) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + ( iz ) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + ( iz ) * (m_nelx(iy)+1);
out(ielem,4) = bot + (ix ) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,5) = bot + (ix+1) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,6) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,7) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1);
ielem++;
// -- bottom-center element
out(ielem,0) = bot + (ix ) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix+1) + (3*iz+1) * (m_nelx(iy)+1);
out(ielem,2) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,3) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,4) = bot + (ix ) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,5) = bot + (ix+1) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,6) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,7) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
ielem++;
// -- bottom-back element
out(ielem,0) = bot + (ix ) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,1) = bot + (ix+1) + (3*iz+2) * (m_nelx(iy)+1);
out(ielem,2) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,3) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,4) = bot + (ix ) + (3*iz+3) * (m_nelx(iy)+1);
out(ielem,5) = bot + (ix+1) + (3*iz+3) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + ( iz+1) * (m_nelx(iy)+1);
ielem++;
// -- top element
out(ielem,0) = mid + (ix ) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,1) = mid + (ix+1) + (2*iz ) * (m_nelx(iy)+1);
out(ielem,2) = top + (ix+1) + ( iz ) * (m_nelx(iy)+1);
out(ielem,3) = top + (ix ) + ( iz ) * (m_nelx(iy)+1);
out(ielem,4) = mid + (ix ) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,5) = mid + (ix+1) + (2*iz+1) * (m_nelx(iy)+1);
out(ielem,6) = top + (ix+1) + ( iz+1) * (m_nelx(iy)+1);
out(ielem,7) = top + (ix ) + ( iz+1) * (m_nelx(iy)+1);
ielem++;
}
}
}
}
return out;
}
// ------------------------------ element numbers of the middle layer ------------------------------
inline xt::xtensor<size_t,1> FineLayer::elementsMiddleLayer() const
{
// number of element layers in y-direction, the index of the middle layer
size_t nely = static_cast<size_t>(m_nhy.size());
size_t iy = (nely-1)/2;
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelx(iy)*m_nelz(iy)});
for ( size_t ix = 0 ; ix < m_nelx(iy) ; ++ix )
for ( size_t iz = 0 ; iz < m_nelz(iy) ; ++iz )
out(ix+iz*m_nelx(iy)) = m_startElem(iy) + ix + iz*m_nelx(iy);
return out;
}
// ------------------------------ node-numbers along the front plane -------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesFront() const
{
// number of element layers in y-direction
size_t nely = static_cast<size_t>(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1;
else n += m_nelx(iy) + 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1;
else n += m_nelx(iy) + 1;
}
// allocate node-list
xt::xtensor<size_t,1> out = xt::empty<size_t>({n});
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(j) = m_startNode(iy) + ix;
++j;
}
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy);
++j;
}
}
// -- top node layer
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(j) = m_startNode(iy+1) + ix;
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the back plane -------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBack() const
{
// number of element layers in y-direction
size_t nely = static_cast<size_t>(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1;
else n += m_nelx(iy) + 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 + 1;
else n += m_nelx(iy) + 1;
}
// allocate node-list
xt::xtensor<size_t,1> out = xt::empty<size_t>({n});
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(j) = m_startNode(iy) + ix + (m_nelx(iy)+1)*m_nelz(iy);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy);
++j;
}
}
// -- top node layer
for ( size_t ix = 0 ; ix < m_nelx(iy)+1 ; ++ix ) {
out(j) = m_startNode(iy+1) + ix + (m_nelx(iy)+1)*m_nelz(iy);
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the left plane -------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesLeft() const
{
// number of element layers in y-direction
size_t nely = static_cast<size_t>(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1;
else n += m_nelz(iy) + 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1;
else n += m_nelz(iy) + 1;
}
// allocate node-list
xt::xtensor<size_t,1> out = xt::empty<size_t>({n});
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy);
++j;
}
}
// -- top node layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1);
++j;
}
}
return out;
}
// ------------------------------ node-numbers along the right plane -------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesRight() const
{
// number of element layers in y-direction
size_t nely = static_cast<size_t>(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1;
else n += m_nelz(iy) + 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 + 1;
else n += m_nelz(iy) + 1;
}
// allocate node-list
xt::xtensor<size_t,1> out = xt::empty<size_t>({n});
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
// -- top node layer
for ( size_t iz = 0 ; iz < m_nelz(iy)+1 ; ++iz ) {
out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
return out;
}
// ------------------------------ node-numbers along the bottom plane ------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBottom() const
{
// number of element layers in y-direction
size_t nely = static_cast<size_t>(m_nhy.size());
// allocate node list
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nnd(nely)});
// counter
size_t j = 0;
// fill node list
for ( size_t ix = 0 ; ix < m_nelx(0)+1 ; ++ix ) {
for ( size_t iz = 0 ; iz < m_nelz(0)+1 ; ++iz ) {
out(j) = m_startNode(0) + ix + iz * (m_nelx(0)+1);
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the top plane --------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesTop() const
{
// number of element layers in y-direction
size_t nely = static_cast<size_t>(m_nhy.size());
// allocate node list
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nnd(nely)});
// counter
size_t j = 0;
// fill node list
for ( size_t ix = 0 ; ix < m_nelx(nely-1)+1 ; ++ix ) {
for ( size_t iz = 0 ; iz < m_nelz(nely-1)+1 ; ++iz ) {
out(j) = m_startNode(nely) + ix + iz * (m_nelx(nely-1)+1);
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the front face -------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesFrontFace() const
{
// number of element layers in y-direction
size_t nely = static_cast<size_t>(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1;
else n += m_nelx(iy) - 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1;
else n += m_nelx(iy) - 1;
}
// allocate node-list
xt::xtensor<size_t,1> out = xt::empty<size_t>({n});
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix;
++j;
}
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy);
++j;
}
}
// -- top node layer
for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy+1) + ix;
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the back face --------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBackFace() const
{
// number of element layers in y-direction
size_t nely = static_cast<size_t>(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1;
else n += m_nelx(iy) - 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
if ( m_refine(iy) == 0 ) n += m_nelx(iy) * 3 - 1;
else n += m_nelx(iy) - 1;
}
// allocate node-list
xt::xtensor<size_t,1> out = xt::empty<size_t>({n});
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + (m_nelx(iy)+1)*m_nelz(iy);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 0 ) {
for ( size_t ix = 0 ; ix < 2*m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy) + ix + m_nnd(iy) + 2*m_nelx(iy)*m_nelz(iy);
++j;
}
}
// -- top node layer
for ( size_t ix = 1 ; ix < m_nelx(iy) ; ++ix ) {
out(j) = m_startNode(iy+1) + ix + (m_nelx(iy)+1)*m_nelz(iy);
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the left face --------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesLeftFace() const
{
// number of element layers in y-direction
size_t nely = static_cast<size_t>(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1;
else n += m_nelz(iy) - 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1;
else n += m_nelz(iy) - 1;
}
// allocate node-list
xt::xtensor<size_t,1> out = xt::empty<size_t>({n});
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nnd(iy);
++j;
}
}
// -- top node layer
for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1);
++j;
}
}
return out;
}
// ------------------------------- node-numbers along the right face -------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesRightFace() const
{
// number of element layers in y-direction
size_t nely = static_cast<size_t>(m_nhy.size());
// number of boundary nodes
// - initialize
size_t n = 0;
// - bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1;
else n += m_nelz(iy) - 1;
}
// - top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
if ( m_refine(iy) == 2 ) n += m_nelz(iy) * 3 - 1;
else n += m_nelz(iy) - 1;
}
// allocate node-list
xt::xtensor<size_t,1> out = xt::empty<size_t>({n});
// initialize counter: current index in the node-list "out"
size_t j = 0;
// bottom half: bottom node layer (+ middle node layer)
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
{
// -- bottom node layer
for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
}
// top half: (middle node layer +) top node layer
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
{
// -- refinement layer
if ( m_refine(iy) == 2 ) {
for ( size_t iz = 0 ; iz < 2*m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
// -- top node layer
for ( size_t iz = 1 ; iz < m_nelz(iy) ; ++iz ) {
out(j) = m_startNode(iy+1) + iz * (m_nelx(iy)+1) + m_nelx(iy);
++j;
}
}
return out;
}
// ------------------------------ node-numbers along the bottom face -------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBottomFace() const
{
// allocate node list
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nelx(0)-1)*(m_nelz(0)-1)});
// counter
size_t j = 0;
// fill node list
for ( size_t ix = 1 ; ix < m_nelx(0) ; ++ix ) {
for ( size_t iz = 1 ; iz < m_nelz(0) ; ++iz ) {
out(j) = m_startNode(0) + ix + iz * (m_nelx(0)+1);
++j;
}
}
return out;
}
// -------------------------------- node-numbers along the top face --------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesTopFace() const
{
// number of element layers in y-direction
size_t nely = static_cast<size_t>(m_nhy.size());
// allocate node list
xt::xtensor<size_t,1> out = xt::empty<size_t>({(m_nelx(nely-1)-1)*(m_nelz(nely-1)-1)});
// counter
size_t j = 0;
// fill node list
for ( size_t ix = 1 ; ix < m_nelx(nely-1) ; ++ix ) {
for ( size_t iz = 1 ; iz < m_nelz(nely-1) ; ++iz ) {
out(j) = m_startNode(nely) + ix + iz * (m_nelx(nely-1)+1);
++j;
}
}
return out;
}
// --------------------------- node-numbers along the front-bottom edge ----------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesFrontBottomEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 front-top edge -----------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesFrontTopEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 front-left edge -----------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesFrontLeftEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 front-right edge ----------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesFrontRightEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 back-bottom edge ----------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBackBottomEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelx(0)+1});
for ( size_t ix = 0 ; ix < m_nelx(0)+1 ; ++ix )
out(ix) = m_startNode(0) + ix + (m_nelx(0)+1)*(m_nelz(0));
return out;
}
// ----------------------------- node-numbers along the back-top edge ------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBackTopEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelx(nely-1)+1});
for ( size_t ix = 0 ; ix < m_nelx(nely-1)+1 ; ++ix )
out(ix) = m_startNode(nely) + ix + (m_nelx(nely-1)+1)*(m_nelz(nely-1));
return out;
}
// ----------------------------- node-numbers along the back-left edge -----------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBackLeftEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({nely+1});
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
out(iy) = m_startNode(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
out(iy+1) = m_startNode(iy+1) + (m_nelx(iy)+1)*(m_nelz(iy));
return out;
}
// ---------------------------- node-numbers along the back-right edge -----------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBackRightEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({nely+1});
for ( size_t iy = 0 ; iy < (nely+1)/2 ; ++iy )
out(iy) = m_startNode(iy) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
for ( size_t iy = (nely-1)/2 ; iy < nely ; ++iy )
out(iy+1) = m_startNode(iy+1) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
return out;
}
// ---------------------------- node-numbers along the bottom-left edge ----------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBottomLeftEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz(0)+1});
for ( size_t iz = 0 ; iz < m_nelz(0)+1 ; ++iz )
out(iz) = m_startNode(0) + iz * (m_nelx(0)+1);
return out;
}
// --------------------------- node-numbers along the bottom-right edge ----------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBottomRightEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz(0)+1});
for ( size_t iz = 0 ; iz < m_nelz(0)+1 ; ++iz )
out(iz) = m_startNode(0) + m_nelx(0) + iz * (m_nelx(0)+1);
return out;
}
// ----------------------------- node-numbers along the top-left edge ------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesTopLeftEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz(nely-1)+1});
for ( size_t iz = 0 ; iz < m_nelz(nely-1)+1 ; ++iz )
out(iz) = m_startNode(nely) + iz * (m_nelx(nely-1)+1);
return out;
}
// ----------------------------- node-numbers along the top-right edge -----------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesTopRightEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz(nely-1)+1});
for ( size_t iz = 0 ; iz < m_nelz(nely-1)+1 ; ++iz )
out(iz) = m_startNode(nely) + m_nelx(nely-1) + iz * (m_nelx(nely-1)+1);
return out;
}
// -------------------------------------------- aliases --------------------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBottomFrontEdge() const { return nodesFrontBottomEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesBottomBackEdge() const { return nodesBackBottomEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesTopFrontEdge() const { return nodesFrontTopEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesTopBackEdge() const { return nodesBackTopEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesLeftBottomEdge() const { return nodesBottomLeftEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesLeftFrontEdge() const { return nodesFrontLeftEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesLeftBackEdge() const { return nodesBackLeftEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesLeftTopEdge() const { return nodesTopLeftEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesRightBottomEdge() const { return nodesBottomRightEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesRightTopEdge() const { return nodesTopRightEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesRightFrontEdge() const { return nodesFrontRightEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesRightBackEdge() const { return nodesBackRightEdge(); }
// ------------------- node-numbers along the front-bottom edge, without corners -------------------
inline xt::xtensor<size_t,1> FineLayer::nodesFrontBottomOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 front-top edge, without corners ---------------------
inline xt::xtensor<size_t,1> FineLayer::nodesFrontTopOpenEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 front-left edge, without corners --------------------
inline xt::xtensor<size_t,1> FineLayer::nodesFrontLeftOpenEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({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 front-right edge, without corners --------------------
inline xt::xtensor<size_t,1> FineLayer::nodesFrontRightOpenEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({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-numbers along the back-bottom edge, without corners --------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBackBottomOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelx(0)-1});
for ( size_t ix = 1 ; ix < m_nelx(0) ; ++ix )
out(ix-1) = m_startNode(0) + ix + (m_nelx(0)+1)*(m_nelz(0));
return out;
}
// --------------------- node-numbers along the back-top edge, without corners ---------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBackTopOpenEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelx(nely-1)-1});
for ( size_t ix = 1 ; ix < m_nelx(nely-1) ; ++ix )
out(ix-1) = m_startNode(nely) + ix + (m_nelx(nely-1)+1)*(m_nelz(nely-1));
return out;
}
// -------------------- node-numbers along the back-left edge, without corners ---------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBackLeftOpenEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({nely-1});
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
out(iy-1) = m_startNode(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
out(iy) = m_startNode(iy+1) + (m_nelx(iy)+1)*(m_nelz(iy));
return out;
}
// -------------------- node-numbers along the back-right edge, without corners --------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBackRightOpenEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({nely-1});
for ( size_t iy = 1 ; iy < (nely+1)/2 ; ++iy )
out(iy-1) = m_startNode(iy) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
for ( size_t iy = (nely-1)/2 ; iy < nely-1 ; ++iy )
out(iy) = m_startNode(iy+1) + m_nelx(iy) + (m_nelx(iy)+1)*(m_nelz(iy));
return out;
}
// ------------------- node-numbers along the bottom-left edge, without corners --------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBottomLeftOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz(0)-1});
for ( size_t iz = 1 ; iz < m_nelz(0) ; ++iz )
out(iz-1) = m_startNode(0) + iz * (m_nelx(0)+1);
return out;
}
// ------------------- node-numbers along the bottom-right edge, without corners -------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBottomRightOpenEdge() const
{
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz(0)-1});
for ( size_t iz = 1 ; iz < m_nelz(0) ; ++iz )
out(iz-1) = m_startNode(0) + m_nelx(0) + iz * (m_nelx(0)+1);
return out;
}
// --------------------- node-numbers along the top-left edge, without corners ---------------------
inline xt::xtensor<size_t,1> FineLayer::nodesTopLeftOpenEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz(nely-1)-1});
for ( size_t iz = 1 ; iz < m_nelz(nely-1) ; ++iz )
out(iz-1) = m_startNode(nely) + iz * (m_nelx(nely-1)+1);
return out;
}
// -------------------- node-numbers along the top-right edge, without corners ---------------------
inline xt::xtensor<size_t,1> FineLayer::nodesTopRightOpenEdge() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
xt::xtensor<size_t,1> out = xt::empty<size_t>({m_nelz(nely-1)-1});
for ( size_t iz = 1 ; iz < m_nelz(nely-1) ; ++iz )
out(iz-1) = m_startNode(nely) + m_nelx(nely-1) + iz * (m_nelx(nely-1)+1);
return out;
}
// -------------------------------------------- aliases --------------------------------------------
inline xt::xtensor<size_t,1> FineLayer::nodesBottomFrontOpenEdge() const { return nodesFrontBottomOpenEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesBottomBackOpenEdge() const { return nodesBackBottomOpenEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesTopFrontOpenEdge() const { return nodesFrontTopOpenEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesTopBackOpenEdge() const { return nodesBackTopOpenEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesLeftBottomOpenEdge() const { return nodesBottomLeftOpenEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesLeftFrontOpenEdge() const { return nodesFrontLeftOpenEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesLeftBackOpenEdge() const { return nodesBackLeftOpenEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesLeftTopOpenEdge() const { return nodesTopLeftOpenEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesRightBottomOpenEdge() const { return nodesBottomRightOpenEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesRightTopOpenEdge() const { return nodesTopRightOpenEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesRightFrontOpenEdge() const { return nodesFrontRightOpenEdge(); }
inline xt::xtensor<size_t,1> FineLayer::nodesRightBackOpenEdge() const { return nodesBackRightOpenEdge(); }
// -------------------------- node-number of the front-bottom-left corner --------------------------
inline size_t FineLayer::nodesFrontBottomLeftCorner() const
{
return m_startNode(0);
}
// ------------------------- node-number of the front-bottom-right corner --------------------------
inline size_t FineLayer::nodesFrontBottomRightCorner() const
{
return m_startNode(0) + m_nelx(0);
}
// --------------------------- node-number of the front-top-left corner ----------------------------
inline size_t FineLayer::nodesFrontTopLeftCorner() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
return m_startNode(nely);
}
// --------------------------- node-number of the front-top-right corner ---------------------------
inline size_t FineLayer::nodesFrontTopRightCorner() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
return m_startNode(nely) + m_nelx(nely-1);
}
// -------------------------- node-number of the back-bottom-left corner ---------------------------
inline size_t FineLayer::nodesBackBottomLeftCorner() const
{
return m_startNode(0) + (m_nelx(0)+1)*(m_nelz(0));
}
// -------------------------- node-number of the back-bottom-right corner --------------------------
inline size_t FineLayer::nodesBackBottomRightCorner() const
{
return m_startNode(0) + m_nelx(0) + (m_nelx(0)+1)*(m_nelz(0));
}
// ---------------------------- node-number of the back-top-left corner ----------------------------
inline size_t FineLayer::nodesBackTopLeftCorner() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
return m_startNode(nely) + (m_nelx(nely-1)+1)*(m_nelz(nely-1));
}
// --------------------------- node-number of the back-top-right corner ----------------------------
inline size_t FineLayer::nodesBackTopRightCorner() const
{
size_t nely = static_cast<size_t>(m_nhy.size());
return m_startNode(nely) + m_nelx(nely-1) + (m_nelx(nely-1)+1)*(m_nelz(nely-1));
}
// -------------------------------------------- aliases --------------------------------------------
inline size_t FineLayer::nodesFrontLeftBottomCorner() const { return nodesFrontBottomLeftCorner(); }
inline size_t FineLayer::nodesBottomFrontLeftCorner() const { return nodesFrontBottomLeftCorner(); }
inline size_t FineLayer::nodesBottomLeftFrontCorner() const { return nodesFrontBottomLeftCorner(); }
inline size_t FineLayer::nodesLeftFrontBottomCorner() const { return nodesFrontBottomLeftCorner(); }
inline size_t FineLayer::nodesLeftBottomFrontCorner() const { return nodesFrontBottomLeftCorner(); }
inline size_t FineLayer::nodesFrontRightBottomCorner() const { return nodesFrontBottomRightCorner(); }
inline size_t FineLayer::nodesBottomFrontRightCorner() const { return nodesFrontBottomRightCorner(); }
inline size_t FineLayer::nodesBottomRightFrontCorner() const { return nodesFrontBottomRightCorner(); }
inline size_t FineLayer::nodesRightFrontBottomCorner() const { return nodesFrontBottomRightCorner(); }
inline size_t FineLayer::nodesRightBottomFrontCorner() const { return nodesFrontBottomRightCorner(); }
inline size_t FineLayer::nodesFrontLeftTopCorner() const { return nodesFrontTopLeftCorner(); }
inline size_t FineLayer::nodesTopFrontLeftCorner() const { return nodesFrontTopLeftCorner(); }
inline size_t FineLayer::nodesTopLeftFrontCorner() const { return nodesFrontTopLeftCorner(); }
inline size_t FineLayer::nodesLeftFrontTopCorner() const { return nodesFrontTopLeftCorner(); }
inline size_t FineLayer::nodesLeftTopFrontCorner() const { return nodesFrontTopLeftCorner(); }
inline size_t FineLayer::nodesFrontRightTopCorner() const { return nodesFrontTopRightCorner(); }
inline size_t FineLayer::nodesTopFrontRightCorner() const { return nodesFrontTopRightCorner(); }
inline size_t FineLayer::nodesTopRightFrontCorner() const { return nodesFrontTopRightCorner(); }
inline size_t FineLayer::nodesRightFrontTopCorner() const { return nodesFrontTopRightCorner(); }
inline size_t FineLayer::nodesRightTopFrontCorner() const { return nodesFrontTopRightCorner(); }
inline size_t FineLayer::nodesBackLeftBottomCorner() const { return nodesBackBottomLeftCorner(); }
inline size_t FineLayer::nodesBottomBackLeftCorner() const { return nodesBackBottomLeftCorner(); }
inline size_t FineLayer::nodesBottomLeftBackCorner() const { return nodesBackBottomLeftCorner(); }
inline size_t FineLayer::nodesLeftBackBottomCorner() const { return nodesBackBottomLeftCorner(); }
inline size_t FineLayer::nodesLeftBottomBackCorner() const { return nodesBackBottomLeftCorner(); }
inline size_t FineLayer::nodesBackRightBottomCorner() const { return nodesBackBottomRightCorner(); }
inline size_t FineLayer::nodesBottomBackRightCorner() const { return nodesBackBottomRightCorner(); }
inline size_t FineLayer::nodesBottomRightBackCorner() const { return nodesBackBottomRightCorner(); }
inline size_t FineLayer::nodesRightBackBottomCorner() const { return nodesBackBottomRightCorner(); }
inline size_t FineLayer::nodesRightBottomBackCorner() const { return nodesBackBottomRightCorner(); }
inline size_t FineLayer::nodesBackLeftTopCorner() const { return nodesBackTopLeftCorner(); }
inline size_t FineLayer::nodesTopBackLeftCorner() const { return nodesBackTopLeftCorner(); }
inline size_t FineLayer::nodesTopLeftBackCorner() const { return nodesBackTopLeftCorner(); }
inline size_t FineLayer::nodesLeftBackTopCorner() const { return nodesBackTopLeftCorner(); }
inline size_t FineLayer::nodesLeftTopBackCorner() const { return nodesBackTopLeftCorner(); }
inline size_t FineLayer::nodesBackRightTopCorner() const { return nodesBackTopRightCorner(); }
inline size_t FineLayer::nodesTopBackRightCorner() const { return nodesBackTopRightCorner(); }
inline size_t FineLayer::nodesTopRightBackCorner() const { return nodesBackTopRightCorner(); }
inline size_t FineLayer::nodesRightBackTopCorner() const { return nodesBackTopRightCorner(); }
inline size_t FineLayer::nodesRightTopBackCorner() const { return nodesBackTopRightCorner(); }
// ------------------------------ node-numbers of periodic node-pairs ------------------------------
inline xt::xtensor<size_t,2> FineLayer::nodesPeriodic() const
{
// faces
xt::xtensor<size_t,1> fro = nodesFrontFace();
xt::xtensor<size_t,1> bck = nodesBackFace();
xt::xtensor<size_t,1> lft = nodesLeftFace();
xt::xtensor<size_t,1> rgt = nodesRightFace();
xt::xtensor<size_t,1> bot = nodesBottomFace();
xt::xtensor<size_t,1> top = nodesTopFace();
// edges
xt::xtensor<size_t,1> froBot = nodesFrontBottomOpenEdge();
xt::xtensor<size_t,1> froTop = nodesFrontTopOpenEdge();
xt::xtensor<size_t,1> froLft = nodesFrontLeftOpenEdge();
xt::xtensor<size_t,1> froRgt = nodesFrontRightOpenEdge();
xt::xtensor<size_t,1> bckBot = nodesBackBottomOpenEdge();
xt::xtensor<size_t,1> bckTop = nodesBackTopOpenEdge();
xt::xtensor<size_t,1> bckLft = nodesBackLeftOpenEdge();
xt::xtensor<size_t,1> bckRgt = nodesBackRightOpenEdge();
xt::xtensor<size_t,1> botLft = nodesBottomLeftOpenEdge();
xt::xtensor<size_t,1> botRgt = nodesBottomRightOpenEdge();
xt::xtensor<size_t,1> topLft = nodesTopLeftOpenEdge();
xt::xtensor<size_t,1> topRgt = nodesTopRightOpenEdge();
// allocate nodal ties
// - number of tying per category
size_t tface = fro.size() + lft.size() + bot.size();
size_t tedge = 3*froBot.size() + 3*froLft.size() + 3*botLft.size();
size_t tnode = 7;
// - allocate
xt::xtensor<size_t,2> out = xt::empty<size_t>({tface+tedge+tnode, std::size_t(2)});
// counter
size_t i = 0;
// tie all corners to one corner
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontBottomRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackBottomLeftCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopLeftCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesFrontTopRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopRightCorner(); ++i;
out(i,0) = nodesFrontBottomLeftCorner(); out(i,1) = nodesBackTopLeftCorner(); ++i;
// tie all corresponding edges to each other (exclude corners)
for ( size_t j = 0 ; j<froBot.size() ; ++j ){ out(i,0) = froBot(j); out(i,1) = bckBot(j); ++i; }
for ( size_t j = 0 ; j<froBot.size() ; ++j ){ out(i,0) = froBot(j); out(i,1) = bckTop(j); ++i; }
for ( size_t j = 0 ; j<froBot.size() ; ++j ){ out(i,0) = froBot(j); out(i,1) = froTop(j); ++i; }
for ( size_t j = 0 ; j<botLft.size() ; ++j ){ out(i,0) = botLft(j); out(i,1) = botRgt(j); ++i; }
for ( size_t j = 0 ; j<botLft.size() ; ++j ){ out(i,0) = botLft(j); out(i,1) = topRgt(j); ++i; }
for ( size_t j = 0 ; j<botLft.size() ; ++j ){ out(i,0) = botLft(j); out(i,1) = topLft(j); ++i; }
for ( size_t j = 0 ; j<froLft.size() ; ++j ){ out(i,0) = froLft(j); out(i,1) = froRgt(j); ++i; }
for ( size_t j = 0 ; j<froLft.size() ; ++j ){ out(i,0) = froLft(j); out(i,1) = bckRgt(j); ++i; }
for ( size_t j = 0 ; j<froLft.size() ; ++j ){ out(i,0) = froLft(j); out(i,1) = bckLft(j); ++i; }
// tie faces to each-other
for ( size_t j = 0 ; j<fro.size() ; ++j ){ out(i,0) = fro(j); out(i,1) = bck(j); ++i; }
for ( size_t j = 0 ; j<lft.size() ; ++j ){ out(i,0) = lft(j); out(i,1) = rgt(j); ++i; }
for ( size_t j = 0 ; j<bot.size() ; ++j ){ out(i,0) = bot(j); out(i,1) = top(j); ++i; }
return out;
}
// ------------------------------ node-number that lies in the origin ------------------------------
inline size_t FineLayer::nodesOrigin() const
{
return nodesFrontBottomLeftCorner();
}
// ------------------------- DOF numbers per node (sequentially numbered) --------------------------
inline xt::xtensor<size_t,2> FineLayer::dofs() const
{
return xGooseFEM::Mesh::dofs(m_nnode,m_ndim);
}
// ------------------------ DOP-numbers with periodic dependencies removed -------------------------
inline xt::xtensor<size_t,2> FineLayer::dofsPeriodic() const
{
// DOF-numbers for each component of each node (sequential)
xt::xtensor<size_t,2> out = xGooseFEM::Mesh::dofs(m_nnode,m_ndim);
// periodic node-pairs
xt::xtensor<size_t,2> 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 xGooseFEM::Mesh::renumber(out);
}
// -------------------------------------------------------------------------------------------------
}}} // namespace ...
// =================================================================================================
#endif
Event Timeline
Log In to Comment