diff --git a/include/GooseFEM/MeshHex8.hpp b/include/GooseFEM/MeshHex8.hpp
index 68f7b14..0e340ea 100644
--- a/include/GooseFEM/MeshHex8.hpp
+++ b/include/GooseFEM/MeshHex8.hpp
@@ -1,3017 +1,3017 @@
 /**
 Implementation of MeshHex8.h
 
 \file MeshHex8.hpp
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_MESHHEX8_HPP
 #define GOOSEFEM_MESHHEX8_HPP
 
 #include "MeshHex8.h"
 
 namespace GooseFEM {
 namespace Mesh {
 namespace Hex8 {
 
 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)
 {
     GOOSEFEM_ASSERT(m_nelx >= 1ul);
     GOOSEFEM_ASSERT(m_nely >= 1ul);
     GOOSEFEM_ASSERT(m_nelz >= 1ul);
 
     m_nnode = (m_nelx + 1) * (m_nely + 1) * (m_nelz + 1);
     m_nelem = m_nelx * m_nely * m_nelz;
 }
 
 inline size_t Regular::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t Regular::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t Regular::nne() const
 {
     return m_nne;
 }
 
 inline size_t Regular::ndim() const
 {
     return m_ndim;
 }
 
 inline ElementType Regular::getElementType() const
 {
     return ElementType::Hex8;
 }
 
 inline xt::xtensor<double, 2> Regular::coor() const
 {
     xt::xtensor<double, 2> ret = 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) {
                 ret(inode, 0) = x(ix);
                 ret(inode, 1) = y(iy);
                 ret(inode, 2) = z(iz);
                 ++inode;
             }
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 2> Regular::conn() const
 {
     xt::xtensor<size_t, 2> ret = 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) {
                 ret(ielem, 0) = iy * (m_nelx + 1) + ix + iz * (m_nely + 1) * (m_nelx + 1);
                 ret(ielem, 1) = iy * (m_nelx + 1) + (ix + 1) + iz * (m_nely + 1) * (m_nelx + 1);
                 ret(ielem, 3) = (iy + 1) * (m_nelx + 1) + ix + iz * (m_nely + 1) * (m_nelx + 1);
                 ret(ielem, 2) =
                     (iy + 1) * (m_nelx + 1) + (ix + 1) + iz * (m_nely + 1) * (m_nelx + 1);
                 ret(ielem, 4) = iy * (m_nelx + 1) + ix + (iz + 1) * (m_nely + 1) * (m_nelx + 1);
                 ret(ielem, 5) =
                     (iy) * (m_nelx + 1) + (ix + 1) + (iz + 1) * (m_nely + 1) * (m_nelx + 1);
                 ret(ielem, 7) =
                     (iy + 1) * (m_nelx + 1) + ix + (iz + 1) * (m_nely + 1) * (m_nelx + 1);
                 ret(ielem, 6) =
                     (iy + 1) * (m_nelx + 1) + (ix + 1) + (iz + 1) * (m_nely + 1) * (m_nelx + 1);
                 ++ielem;
             }
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesFront() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret(iy * (m_nelx + 1) + ix) = iy * (m_nelx + 1) + ix;
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBack() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret(iy * (m_nelx + 1) + ix) =
                 iy * (m_nelx + 1) + ix + m_nelz * (m_nely + 1) * (m_nelx + 1);
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesLeft() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret(iz * (m_nely + 1) + iy) = iy * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1);
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesRight() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret(iz * (m_nely + 1) + iy) =
                 iy * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1) + m_nelx;
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBottom() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret(iz * (m_nelx + 1) + ix) = ix + iz * (m_nelx + 1) * (m_nely + 1);
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesTop() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret(iz * (m_nelx + 1) + ix) =
                 ix + m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1);
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesFrontFace() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret((iy - 1) * (m_nelx - 1) + (ix - 1)) = iy * (m_nelx + 1) + ix;
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBackFace() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret((iy - 1) * (m_nelx - 1) + (ix - 1)) =
                 iy * (m_nelx + 1) + ix + m_nelz * (m_nely + 1) * (m_nelx + 1);
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesLeftFace() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret((iz - 1) * (m_nely - 1) + (iy - 1)) =
                 iy * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1);
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesRightFace() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret((iz - 1) * (m_nely - 1) + (iy - 1)) =
                 iy * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1) + m_nelx;
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBottomFace() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret((iz - 1) * (m_nelx - 1) + (ix - 1)) = ix + iz * (m_nelx + 1) * (m_nely + 1);
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesTopFace() const
 {
     xt::xtensor<size_t, 1> ret = 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) {
             ret((iz - 1) * (m_nelx - 1) + (ix - 1)) =
                 ix + m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1);
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesFrontBottomEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx + 1});
 
     for (size_t ix = 0; ix < m_nelx + 1; ++ix) {
         ret(ix) = ix;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesFrontTopEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx + 1});
 
     for (size_t ix = 0; ix < m_nelx + 1; ++ix) {
         ret(ix) = ix + m_nely * (m_nelx + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesFrontLeftEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nely + 1});
 
     for (size_t iy = 0; iy < m_nely + 1; ++iy) {
         ret(iy) = iy * (m_nelx + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesFrontRightEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nely + 1});
 
     for (size_t iy = 0; iy < m_nely + 1; ++iy) {
         ret(iy) = iy * (m_nelx + 1) + m_nelx;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBackBottomEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx + 1});
 
     for (size_t ix = 0; ix < m_nelx + 1; ++ix) {
         ret(ix) = ix + m_nelz * (m_nely + 1) * (m_nelx + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBackTopEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx + 1});
 
     for (size_t ix = 0; ix < m_nelx + 1; ++ix) {
         ret(ix) = m_nely * (m_nelx + 1) + ix + m_nelz * (m_nely + 1) * (m_nelx + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBackLeftEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nely + 1});
 
     for (size_t iy = 0; iy < m_nely + 1; ++iy) {
         ret(iy) = iy * (m_nelx + 1) + m_nelz * (m_nelx + 1) * (m_nely + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBackRightEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nely + 1});
 
     for (size_t iy = 0; iy < m_nely + 1; ++iy) {
         ret(iy) = iy * (m_nelx + 1) + m_nelz * (m_nelx + 1) * (m_nely + 1) + m_nelx;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBottomLeftEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz + 1});
 
     for (size_t iz = 0; iz < m_nelz + 1; ++iz) {
         ret(iz) = iz * (m_nelx + 1) * (m_nely + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBottomRightEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz + 1});
 
     for (size_t iz = 0; iz < m_nelz + 1; ++iz) {
         ret(iz) = iz * (m_nelx + 1) * (m_nely + 1) + m_nelx;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesTopLeftEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz + 1});
 
     for (size_t iz = 0; iz < m_nelz + 1; ++iz) {
         ret(iz) = m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesTopRightEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz + 1});
 
     for (size_t iz = 0; iz < m_nelz + 1; ++iz) {
         ret(iz) = m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1) + m_nelx;
     }
 
     return ret;
 }
 
 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();
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesFrontBottomOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx - 1});
 
     for (size_t ix = 1; ix < m_nelx; ++ix) {
         ret(ix - 1) = ix;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesFrontTopOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx - 1});
 
     for (size_t ix = 1; ix < m_nelx; ++ix) {
         ret(ix - 1) = ix + m_nely * (m_nelx + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesFrontLeftOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nely - 1});
 
     for (size_t iy = 1; iy < m_nely; ++iy) {
         ret(iy - 1) = iy * (m_nelx + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesFrontRightOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nely - 1});
 
     for (size_t iy = 1; iy < m_nely; ++iy) {
         ret(iy - 1) = iy * (m_nelx + 1) + m_nelx;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBackBottomOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx - 1});
 
     for (size_t ix = 1; ix < m_nelx; ++ix) {
         ret(ix - 1) = ix + m_nelz * (m_nely + 1) * (m_nelx + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBackTopOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx - 1});
 
     for (size_t ix = 1; ix < m_nelx; ++ix) {
         ret(ix - 1) = m_nely * (m_nelx + 1) + ix + m_nelz * (m_nely + 1) * (m_nelx + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBackLeftOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nely - 1});
 
     for (size_t iy = 1; iy < m_nely; ++iy) {
         ret(iy - 1) = iy * (m_nelx + 1) + m_nelz * (m_nelx + 1) * (m_nely + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBackRightOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nely - 1});
 
     for (size_t iy = 1; iy < m_nely; ++iy) {
         ret(iy - 1) = iy * (m_nelx + 1) + m_nelz * (m_nelx + 1) * (m_nely + 1) + m_nelx;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBottomLeftOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz - 1});
 
     for (size_t iz = 1; iz < m_nelz; ++iz) {
         ret(iz - 1) = iz * (m_nelx + 1) * (m_nely + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBottomRightOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz - 1});
 
     for (size_t iz = 1; iz < m_nelz; ++iz) {
         ret(iz - 1) = iz * (m_nelx + 1) * (m_nely + 1) + m_nelx;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesTopLeftOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz - 1});
 
     for (size_t iz = 1; iz < m_nelz; ++iz) {
         ret(iz - 1) = m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesTopRightOpenEdge() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz - 1});
 
     for (size_t iz = 1; iz < m_nelz; ++iz) {
         ret(iz - 1) = m_nely * (m_nelx + 1) + iz * (m_nelx + 1) * (m_nely + 1) + m_nelx;
     }
 
     return ret;
 }
 
 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();
 }
 
 inline size_t Regular::nodesFrontBottomLeftCorner() const
 {
     return 0;
 }
 
 inline size_t Regular::nodesFrontBottomRightCorner() const
 {
     return m_nelx;
 }
 
 inline size_t Regular::nodesFrontTopLeftCorner() const
 {
     return m_nely * (m_nelx + 1);
 }
 
 inline size_t Regular::nodesFrontTopRightCorner() const
 {
     return m_nely * (m_nelx + 1) + m_nelx;
 }
 
 inline size_t Regular::nodesBackBottomLeftCorner() const
 {
     return m_nelz * (m_nely + 1) * (m_nelx + 1);
 }
 
 inline size_t Regular::nodesBackBottomRightCorner() const
 {
     return m_nelx + m_nelz * (m_nely + 1) * (m_nelx + 1);
 }
 
 inline size_t Regular::nodesBackTopLeftCorner() const
 {
     return m_nely * (m_nelx + 1) + m_nelz * (m_nely + 1) * (m_nelx + 1);
 }
 
 inline size_t Regular::nodesBackTopRightCorner() const
 {
     return m_nely * (m_nelx + 1) + m_nelx + m_nelz * (m_nely + 1) * (m_nelx + 1);
 }
 
 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();
 }
 
 inline xt::xtensor<size_t, 2> Regular::nodesPeriodic() const
 {
     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();
 
     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();
 
     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;
     xt::xtensor<size_t, 2> ret = xt::empty<size_t>({tface + tedge + tnode, std::size_t(2)});
 
     size_t i = 0;
 
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesFrontBottomRightCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesBackBottomRightCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesBackBottomLeftCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesFrontTopLeftCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesFrontTopRightCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesBackTopRightCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesBackTopLeftCorner();
     ++i;
 
     for (size_t j = 0; j < froBot.size(); ++j) {
         ret(i, 0) = froBot(j);
         ret(i, 1) = bckBot(j);
         ++i;
     }
     for (size_t j = 0; j < froBot.size(); ++j) {
         ret(i, 0) = froBot(j);
         ret(i, 1) = bckTop(j);
         ++i;
     }
     for (size_t j = 0; j < froBot.size(); ++j) {
         ret(i, 0) = froBot(j);
         ret(i, 1) = froTop(j);
         ++i;
     }
     for (size_t j = 0; j < botLft.size(); ++j) {
         ret(i, 0) = botLft(j);
         ret(i, 1) = botRgt(j);
         ++i;
     }
     for (size_t j = 0; j < botLft.size(); ++j) {
         ret(i, 0) = botLft(j);
         ret(i, 1) = topRgt(j);
         ++i;
     }
     for (size_t j = 0; j < botLft.size(); ++j) {
         ret(i, 0) = botLft(j);
         ret(i, 1) = topLft(j);
         ++i;
     }
     for (size_t j = 0; j < froLft.size(); ++j) {
         ret(i, 0) = froLft(j);
         ret(i, 1) = froRgt(j);
         ++i;
     }
     for (size_t j = 0; j < froLft.size(); ++j) {
         ret(i, 0) = froLft(j);
         ret(i, 1) = bckRgt(j);
         ++i;
     }
     for (size_t j = 0; j < froLft.size(); ++j) {
         ret(i, 0) = froLft(j);
         ret(i, 1) = bckLft(j);
         ++i;
     }
 
     for (size_t j = 0; j < fro.size(); ++j) {
         ret(i, 0) = fro(j);
         ret(i, 1) = bck(j);
         ++i;
     }
     for (size_t j = 0; j < lft.size(); ++j) {
         ret(i, 0) = lft(j);
         ret(i, 1) = rgt(j);
         ++i;
     }
     for (size_t j = 0; j < bot.size(); ++j) {
         ret(i, 0) = bot(j);
         ret(i, 1) = top(j);
         ++i;
     }
 
     return ret;
 }
 
 inline size_t Regular::nodesOrigin() const
 {
     return nodesFrontBottomLeftCorner();
 }
 
 inline xt::xtensor<size_t, 2> Regular::dofs() const
 {
     return GooseFEM::Mesh::dofs(m_nnode, m_ndim);
 }
 
 inline xt::xtensor<size_t, 2> Regular::dofsPeriodic() const
 {
     xt::xtensor<size_t, 2> ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim);
     xt::xtensor<size_t, 2> nodePer = nodesPeriodic();
 
     for (size_t i = 0; i < nodePer.shape(0); ++i) {
         for (size_t j = 0; j < m_ndim; ++j) {
             ret(nodePer(i, 1), j) = ret(nodePer(i, 0), j);
         }
     }
 
     return GooseFEM::Mesh::renumber(ret);
 }
 
 inline FineLayer::FineLayer(size_t nelx, size_t nely, size_t nelz, double h, size_t nfine) : m_h(h)
 {
     // basic assumptions
     GOOSEFEM_ASSERT(nelx >= 1ul);
     GOOSEFEM_ASSERT(nely >= 1ul);
     GOOSEFEM_ASSERT(nelz >= 1ul);
 
     // 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 || ntot >= nmin) {
             nely = iy;
             break;
         }
 
         // rules (1,2) satisfied: coarsen in x-direction (and z-direction)
         if (3 * nhy(iy) <= ntot && nelx % (3 * nhx(iy)) == 0 && 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 && 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 && nelz % (3 * nhz(iy)) == 0 && 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 || 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_layer_nelx = xt::empty<size_t>({nely * 2 - 1});
+    m_layer_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);
+        m_layer_nelx(iy) = nelx / m_nhx(iy);
+        m_layer_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);
+        m_nnd(iy) = (m_layer_nelx(iy) + 1) * (m_layer_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);
+        m_nnd(iy + 1) = (m_layer_nelx(iy) + 1) * (m_layer_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);
+            m_nnode += (3 * m_layer_nelx(i) + 1) * (m_layer_nelz(i) + 1);
         }
         else if (m_refine(i) == 2) {
-            m_nnode += (m_nelx(i) + 1) * (3 * m_nelz(i) + 1);
+            m_nnode += (m_layer_nelx(i) + 1) * (3 * m_layer_nelz(i) + 1);
         }
         else {
-            m_nnode += (m_nelx(i) + 1) * (m_nelz(i) + 1);
+            m_nnode += (m_layer_nelx(i) + 1) * (m_layer_nelz(i) + 1);
         }
         // - add the elements of this layer
         if (m_refine(i) == 0) {
-            m_nelem += (4 * m_nelx(i)) * (m_nelz(i));
+            m_nelem += (4 * m_layer_nelx(i)) * (m_layer_nelz(i));
         }
         else if (m_refine(i) == 2) {
-            m_nelem += (m_nelx(i)) * (4 * m_nelz(i));
+            m_nelem += (m_layer_nelx(i)) * (4 * m_layer_nelz(i));
         }
         else {
-            m_nelem += (m_nelx(i)) * (m_nelz(i));
+            m_nelem += (m_layer_nelx(i)) * (m_layer_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);
+            m_nnode += (5 * m_layer_nelx(i) + 1) * (m_layer_nelz(i) + 1);
         }
         else if (m_refine(i) == 2) {
-            m_nnode += (m_nelx(i) + 1) * (5 * m_nelz(i) + 1);
+            m_nnode += (m_layer_nelx(i) + 1) * (5 * m_layer_nelz(i) + 1);
         }
         else {
-            m_nnode += (m_nelx(i) + 1) * (m_nelz(i) + 1);
+            m_nnode += (m_layer_nelx(i) + 1) * (m_layer_nelz(i) + 1);
         }
         // - add the elements of this layer
         if (m_refine(i) == 0) {
-            m_nelem += (4 * m_nelx(i)) * (m_nelz(i));
+            m_nelem += (4 * m_layer_nelx(i)) * (m_layer_nelz(i));
         }
         else if (m_refine(i) == 2) {
-            m_nelem += (m_nelx(i)) * (4 * m_nelz(i));
+            m_nelem += (m_layer_nelx(i)) * (4 * m_layer_nelz(i));
         }
         else {
-            m_nelem += (m_nelx(i)) * (m_nelz(i));
+            m_nelem += (m_layer_nelx(i)) * (m_layer_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);
+    m_nnode += (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1) + 1);
 }
 
 inline size_t FineLayer::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t FineLayer::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t FineLayer::nne() const
 {
     return m_nne;
 }
 
 inline size_t FineLayer::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t FineLayer::nelx() const
 {
-    return xt::amax(m_nelx)();
+    return xt::amax(m_layer_nelx)();
 }
 
 inline size_t FineLayer::nely() const
 {
     return xt::sum(m_nhy)();
 }
 
 inline size_t FineLayer::nelz() const
 {
-    return xt::amax(m_nelz)();
+    return xt::amax(m_layer_nelz)();
 }
 
 inline ElementType FineLayer::getElementType() const
 {
     return ElementType::Hex8;
 }
 
 inline xt::xtensor<double, 2> FineLayer::coor() const
 {
     // allocate output
     xt::xtensor<double, 2> ret = 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);
+        xt::xtensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_layer_nelx(iy) + 1);
+        xt::xtensor<double, 1> z = xt::linspace<double>(0.0, m_Lz, m_layer_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) {
+        for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) {
+            for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) {
                 ret(inode, 0) = x(ix);
                 ret(inode, 1) = y(iy);
                 ret(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 iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) {
+                for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
                     for (size_t j = 0; j < 2; ++j) {
                         ret(inode, 0) = x(ix) + dx * static_cast<double>(j + 1);
                         ret(inode, 1) = y(iy) + dy;
                         ret(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 iz = 0; iz < m_layer_nelz(iy); ++iz) {
                 for (size_t j = 0; j < 2; ++j) {
-                    for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
+                    for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) {
                         ret(inode, 0) = x(ix);
                         ret(inode, 1) = y(iy) + dy;
                         ret(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);
+        xt::xtensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_layer_nelx(iy) + 1);
+        xt::xtensor<double, 1> z = xt::linspace<double>(0.0, m_Lz, m_layer_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 iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) {
+                for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
                     for (size_t j = 0; j < 2; ++j) {
                         ret(inode, 0) = x(ix) + dx * static_cast<double>(j + 1);
                         ret(inode, 1) = y(iy) + dy;
                         ret(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 iz = 0; iz < m_layer_nelz(iy); ++iz) {
                 for (size_t j = 0; j < 2; ++j) {
-                    for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
+                    for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) {
                         ret(inode, 0) = x(ix);
                         ret(inode, 1) = y(iy) + dy;
                         ret(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) {
+        for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) {
+            for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) {
                 ret(inode, 0) = x(ix);
                 ret(inode, 1) = y(iy + 1);
                 ret(inode, 2) = z(iz);
                 ++inode;
             }
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 2> FineLayer::conn() const
 {
     // allocate output
     xt::xtensor<size_t, 2> ret = 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) {
-                    ret(ielem, 0) = bot + ix + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = bot + (ix + 1) + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 2) = top + (ix + 1) + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 3) = top + ix + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 4) = bot + ix + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = bot + (ix + 1) + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 7) = top + ix + (iz + 1) * (m_nelx(iy) + 1);
+            for (size_t iz = 0; iz < m_layer_nelz(iy); ++iz) {
+                for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
+                    ret(ielem, 0) = bot + ix + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = bot + (ix + 1) + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = top + (ix + 1) + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = top + ix + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = bot + ix + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = bot + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = top + ix + (iz + 1) * (m_layer_nelx(iy) + 1);
                     ielem++;
                 }
             }
         }
 
         // - define connectivity: refinement along the x-direction (below the middle layer)
         else if (m_refine(iy) == 0 && iy <= (nely - 1) / 2) {
-            for (size_t iz = 0; iz < m_nelz(iy); ++iz) {
-                for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
+            for (size_t iz = 0; iz < m_layer_nelz(iy); ++iz) {
+                for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
                     // -- bottom element
-                    ret(ielem, 0) = bot + ix + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = bot + (ix + 1) + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 2) = mid + (2 * ix + 1) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 3) = mid + (2 * ix) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 4) = bot + ix + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = bot + (ix + 1) + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 6) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_nelx(iy));
-                    ret(ielem, 7) = mid + (2 * ix) + (iz + 1) * (2 * m_nelx(iy));
+                    ret(ielem, 0) = bot + ix + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = bot + (ix + 1) + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 3) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 4) = bot + ix + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = bot + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy));
+                    ret(ielem, 7) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy));
                     ielem++;
                     // -- top-right element
-                    ret(ielem, 0) = bot + (ix + 1) + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = top + (3 * ix + 3) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 2) = top + (3 * ix + 2) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 3) = mid + (2 * ix + 1) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 4) = bot + (ix + 1) + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = top + (3 * ix + 3) + (iz + 1) * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 6) = top + (3 * ix + 2) + (iz + 1) * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 7) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_nelx(iy));
+                    ret(ielem, 0) = bot + (ix + 1) + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = top + (3 * ix + 3) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = top + (3 * ix + 2) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 4) = bot + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = top + (3 * ix + 3) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = top + (3 * ix + 2) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy));
                     ielem++;
                     // -- top-center element
-                    ret(ielem, 0) = mid + (2 * ix) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 1) = mid + (2 * ix + 1) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 2) = top + (3 * ix + 2) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 3) = top + (3 * ix + 1) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 4) = mid + (2 * ix) + (iz + 1) * (2 * m_nelx(iy));
-                    ret(ielem, 5) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_nelx(iy));
-                    ret(ielem, 6) = top + (3 * ix + 2) + (iz + 1) * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 7) = top + (3 * ix + 1) + (iz + 1) * (3 * m_nelx(iy) + 1);
+                    ret(ielem, 0) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 1) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 2) = top + (3 * ix + 2) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = top + (3 * ix + 1) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy));
+                    ret(ielem, 5) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy));
+                    ret(ielem, 6) = top + (3 * ix + 2) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = top + (3 * ix + 1) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
                     ielem++;
                     // -- top-left element
-                    ret(ielem, 0) = bot + ix + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = mid + (2 * ix) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 2) = top + (3 * ix + 1) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 3) = top + (3 * ix) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 4) = bot + ix + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = mid + (2 * ix) + (iz + 1) * (2 * m_nelx(iy));
-                    ret(ielem, 6) = top + (3 * ix + 1) + (iz + 1) * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 7) = top + (3 * ix) + (iz + 1) * (3 * m_nelx(iy) + 1);
+                    ret(ielem, 0) = bot + ix + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 2) = top + (3 * ix + 1) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = top + (3 * ix) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = bot + ix + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy));
+                    ret(ielem, 6) = top + (3 * ix + 1) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = top + (3 * ix) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
                     ielem++;
                 }
             }
         }
 
         // - define connectivity: coarsening along the x-direction (above the middle layer)
         else if (m_refine(iy) == 0 && iy > (nely - 1) / 2) {
-            for (size_t iz = 0; iz < m_nelz(iy); ++iz) {
-                for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
+            for (size_t iz = 0; iz < m_layer_nelz(iy); ++iz) {
+                for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
                     // -- lower-left element
-                    ret(ielem, 0) = bot + (3 * ix) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 1) = bot + (3 * ix + 1) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 2) = mid + (2 * ix) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 3) = top + ix + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 4) = bot + (3 * ix) + (iz + 1) * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 5) = bot + (3 * ix + 1) + (iz + 1) * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 6) = mid + (2 * ix) + (iz + 1) * (2 * m_nelx(iy));
-                    ret(ielem, 7) = top + ix + (iz + 1) * (m_nelx(iy) + 1);
+                    ret(ielem, 0) = bot + (3 * ix) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = bot + (3 * ix + 1) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 3) = top + ix + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = bot + (3 * ix) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = bot + (3 * ix + 1) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy));
+                    ret(ielem, 7) = top + ix + (iz + 1) * (m_layer_nelx(iy) + 1);
                     ielem++;
                     // -- lower-center element
-                    ret(ielem, 0) = bot + (3 * ix + 1) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 1) = bot + (3 * ix + 2) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 2) = mid + (2 * ix + 1) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 3) = mid + (2 * ix) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 4) = bot + (3 * ix + 1) + (iz + 1) * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 5) = bot + (3 * ix + 2) + (iz + 1) * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 6) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_nelx(iy));
-                    ret(ielem, 7) = mid + (2 * ix) + (iz + 1) * (2 * m_nelx(iy));
+                    ret(ielem, 0) = bot + (3 * ix + 1) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = bot + (3 * ix + 2) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 3) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 4) = bot + (3 * ix + 1) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = bot + (3 * ix + 2) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy));
+                    ret(ielem, 7) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy));
                     ielem++;
                     // -- lower-right element
-                    ret(ielem, 0) = bot + (3 * ix + 2) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 1) = bot + (3 * ix + 3) + iz * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 2) = top + (ix + 1) + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 3) = mid + (2 * ix + 1) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 4) = bot + (3 * ix + 2) + (iz + 1) * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 5) = bot + (3 * ix + 3) + (iz + 1) * (3 * m_nelx(iy) + 1);
-                    ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 7) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_nelx(iy));
+                    ret(ielem, 0) = bot + (3 * ix + 2) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = bot + (3 * ix + 3) + iz * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = top + (ix + 1) + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 4) = bot + (3 * ix + 2) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = bot + (3 * ix + 3) + (iz + 1) * (3 * m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy));
                     ielem++;
                     // -- upper element
-                    ret(ielem, 0) = mid + (2 * ix) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 1) = mid + (2 * ix + 1) + iz * (2 * m_nelx(iy));
-                    ret(ielem, 2) = top + (ix + 1) + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 3) = top + ix + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 4) = mid + (2 * ix) + (iz + 1) * (2 * m_nelx(iy));
-                    ret(ielem, 5) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_nelx(iy));
-                    ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 7) = top + ix + (iz + 1) * (m_nelx(iy) + 1);
+                    ret(ielem, 0) = mid + (2 * ix) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 1) = mid + (2 * ix + 1) + iz * (2 * m_layer_nelx(iy));
+                    ret(ielem, 2) = top + (ix + 1) + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = top + ix + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = mid + (2 * ix) + (iz + 1) * (2 * m_layer_nelx(iy));
+                    ret(ielem, 5) = mid + (2 * ix + 1) + (iz + 1) * (2 * m_layer_nelx(iy));
+                    ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = top + ix + (iz + 1) * (m_layer_nelx(iy) + 1);
                     ielem++;
                 }
             }
         }
 
         // - define connectivity: refinement along the z-direction (below the middle layer)
         else if (m_refine(iy) == 2 && iy <= (nely - 1) / 2) {
-            for (size_t iz = 0; iz < m_nelz(iy); ++iz) {
-                for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
+            for (size_t iz = 0; iz < m_layer_nelz(iy); ++iz) {
+                for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
                     // -- bottom element
-                    ret(ielem, 0) = bot + ix + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = bot + ix + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 2) = bot + (ix + 1) + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 3) = bot + (ix + 1) + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 4) = mid + ix + 2 * iz * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = mid + ix + (2 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 6) = mid + (ix + 1) + (2 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 7) = mid + (ix + 1) + 2 * iz * (m_nelx(iy) + 1);
+                    ret(ielem, 0) = bot + ix + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = bot + ix + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = bot + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = bot + (ix + 1) + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = mid + ix + 2 * iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = mid + (ix + 1) + 2 * iz * (m_layer_nelx(iy) + 1);
                     ielem++;
                     // -- top-back element
-                    ret(ielem, 0) = mid + ix + (2 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = mid + (ix + 1) + (2 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 2) = top + (ix + 1) + (3 * iz + 2) * (m_nelx(iy) + 1);
-                    ret(ielem, 3) = top + ix + (3 * iz + 2) * (m_nelx(iy) + 1);
-                    ret(ielem, 4) = bot + ix + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = bot + (ix + 1) + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 6) = top + (ix + 1) + (3 * iz + 3) * (m_nelx(iy) + 1);
-                    ret(ielem, 7) = top + ix + (3 * iz + 3) * (m_nelx(iy) + 1);
+                    ret(ielem, 0) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = top + (ix + 1) + (3 * iz + 2) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = top + ix + (3 * iz + 2) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = bot + ix + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = bot + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = top + (ix + 1) + (3 * iz + 3) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = top + ix + (3 * iz + 3) * (m_layer_nelx(iy) + 1);
                     ielem++;
                     // -- top-center element
-                    ret(ielem, 0) = mid + ix + (2 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = mid + (ix + 1) + (2 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 2) = top + (ix + 1) + (3 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 3) = top + ix + (3 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 4) = mid + ix + (2 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = mid + (ix + 1) + (2 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 6) = top + (ix + 1) + (3 * iz + 2) * (m_nelx(iy) + 1);
-                    ret(ielem, 7) = top + ix + (3 * iz + 2) * (m_nelx(iy) + 1);
+                    ret(ielem, 0) = mid + ix + (2 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = mid + (ix + 1) + (2 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = top + (ix + 1) + (3 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = top + ix + (3 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = top + (ix + 1) + (3 * iz + 2) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = top + ix + (3 * iz + 2) * (m_layer_nelx(iy) + 1);
                     ielem++;
                     // -- top-front element
-                    ret(ielem, 0) = bot + ix + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = bot + (ix + 1) + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 2) = top + (ix + 1) + (3 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 3) = top + ix + (3 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 4) = mid + ix + (2 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = mid + (ix + 1) + (2 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 6) = top + (ix + 1) + (3 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 7) = top + ix + (3 * iz + 1) * (m_nelx(iy) + 1);
+                    ret(ielem, 0) = bot + ix + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = bot + (ix + 1) + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = top + (ix + 1) + (3 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = top + ix + (3 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = mid + ix + (2 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = mid + (ix + 1) + (2 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = top + (ix + 1) + (3 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = top + ix + (3 * iz + 1) * (m_layer_nelx(iy) + 1);
                     ielem++;
                 }
             }
         }
 
         // - define connectivity: coarsening along the z-direction (above the middle layer)
         else if (m_refine(iy) == 2 && iy > (nely - 1) / 2) {
-            for (size_t iz = 0; iz < m_nelz(iy); ++iz) {
-                for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
+            for (size_t iz = 0; iz < m_layer_nelz(iy); ++iz) {
+                for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
                     // -- bottom-front element
-                    ret(ielem, 0) = bot + ix + (3 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = bot + (ix + 1) + (3 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 2) = top + (ix + 1) + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 3) = top + ix + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 4) = bot + ix + (3 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = bot + (ix + 1) + (3 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 6) = mid + (ix + 1) + (2 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 7) = mid + ix + (2 * iz) * (m_nelx(iy) + 1);
+                    ret(ielem, 0) = bot + ix + (3 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = bot + (ix + 1) + (3 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = top + (ix + 1) + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = top + ix + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = bot + ix + (3 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = bot + (ix + 1) + (3 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = mid + (ix + 1) + (2 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = mid + ix + (2 * iz) * (m_layer_nelx(iy) + 1);
                     ielem++;
                     // -- bottom-center element
-                    ret(ielem, 0) = bot + ix + (3 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = bot + (ix + 1) + (3 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 2) = mid + (ix + 1) + (2 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 3) = mid + ix + (2 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 4) = bot + ix + (3 * iz + 2) * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = bot + (ix + 1) + (3 * iz + 2) * (m_nelx(iy) + 1);
-                    ret(ielem, 6) = mid + (ix + 1) + (2 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 7) = mid + ix + (2 * iz + 1) * (m_nelx(iy) + 1);
+                    ret(ielem, 0) = bot + ix + (3 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = bot + (ix + 1) + (3 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = mid + (ix + 1) + (2 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = mid + ix + (2 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = bot + ix + (3 * iz + 2) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = bot + (ix + 1) + (3 * iz + 2) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
                     ielem++;
                     // -- bottom-back element
-                    ret(ielem, 0) = bot + ix + (3 * iz + 2) * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = bot + (ix + 1) + (3 * iz + 2) * (m_nelx(iy) + 1);
-                    ret(ielem, 2) = mid + (ix + 1) + (2 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 3) = mid + ix + (2 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 4) = bot + ix + (3 * iz + 3) * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = bot + (ix + 1) + (3 * iz + 3) * (m_nelx(iy) + 1);
-                    ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 7) = top + ix + (iz + 1) * (m_nelx(iy) + 1);
+                    ret(ielem, 0) = bot + ix + (3 * iz + 2) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = bot + (ix + 1) + (3 * iz + 2) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = bot + ix + (3 * iz + 3) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = bot + (ix + 1) + (3 * iz + 3) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = top + ix + (iz + 1) * (m_layer_nelx(iy) + 1);
                     ielem++;
                     // -- top element
-                    ret(ielem, 0) = mid + ix + (2 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 1) = mid + (ix + 1) + (2 * iz) * (m_nelx(iy) + 1);
-                    ret(ielem, 2) = top + (ix + 1) + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 3) = top + ix + iz * (m_nelx(iy) + 1);
-                    ret(ielem, 4) = mid + ix + (2 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 5) = mid + (ix + 1) + (2 * iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_nelx(iy) + 1);
-                    ret(ielem, 7) = top + ix + (iz + 1) * (m_nelx(iy) + 1);
+                    ret(ielem, 0) = mid + ix + (2 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 1) = mid + (ix + 1) + (2 * iz) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 2) = top + (ix + 1) + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 3) = top + ix + iz * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 4) = mid + ix + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 5) = mid + (ix + 1) + (2 * iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 6) = top + (ix + 1) + (iz + 1) * (m_layer_nelx(iy) + 1);
+                    ret(ielem, 7) = top + ix + (iz + 1) * (m_layer_nelx(iy) + 1);
                     ielem++;
                 }
             }
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::elementsMiddleLayer() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
     size_t iy = (nely - 1) / 2;
 
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx(iy) * m_nelz(iy)});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(iy) * m_layer_nelz(iy)});
 
-    for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
-        for (size_t iz = 0; iz < m_nelz(iy); ++iz) {
-            ret(ix + iz * m_nelx(iy)) = m_startElem(iy) + ix + iz * m_nelx(iy);
+    for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
+        for (size_t iz = 0; iz < m_layer_nelz(iy); ++iz) {
+            ret(ix + iz * m_layer_nelx(iy)) = m_startElem(iy) + ix + iz * m_layer_nelx(iy);
         }
     }
 
     return ret;
 }
 
 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;
+            n += m_layer_nelx(iy) * 3 + 1;
         }
         else {
-            n += m_nelx(iy) + 1;
+            n += m_layer_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;
+            n += m_layer_nelx(iy) * 3 + 1;
         }
         else {
-            n += m_nelx(iy) + 1;
+            n += m_layer_nelx(iy) + 1;
         }
     }
 
     // allocate node-list
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({n});
 
     // initialize counter: current index in the node-list "ret"
     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) {
+        for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) {
             ret(j) = m_startNode(iy) + ix;
             ++j;
         }
         // -- refinement layer
         if (m_refine(iy) == 0) {
-            for (size_t ix = 0; ix < 2 * m_nelx(iy); ++ix) {
+            for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) {
                 ret(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) {
+            for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) {
                 ret(j) = m_startNode(iy) + ix + m_nnd(iy);
                 ++j;
             }
         }
         // -- top node layer
-        for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
+        for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) {
             ret(j) = m_startNode(iy + 1) + ix;
             ++j;
         }
     }
 
     return ret;
 }
 
 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;
+            n += m_layer_nelx(iy) * 3 + 1;
         }
         else {
-            n += m_nelx(iy) + 1;
+            n += m_layer_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;
+            n += m_layer_nelx(iy) * 3 + 1;
         }
         else {
-            n += m_nelx(iy) + 1;
+            n += m_layer_nelx(iy) + 1;
         }
     }
 
     // allocate node-list
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({n});
 
     // initialize counter: current index in the node-list "ret"
     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) {
-            ret(j) = m_startNode(iy) + ix + (m_nelx(iy) + 1) * m_nelz(iy);
+        for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) {
+            ret(j) = m_startNode(iy) + ix + (m_layer_nelx(iy) + 1) * m_layer_nelz(iy);
             ++j;
         }
         // -- refinement layer
         if (m_refine(iy) == 0) {
-            for (size_t ix = 0; ix < 2 * m_nelx(iy); ++ix) {
-                ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_nelx(iy) * m_nelz(iy);
+            for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) {
+                ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_layer_nelx(iy) * m_layer_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) {
-                ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_nelx(iy) * m_nelz(iy);
+            for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) {
+                ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_layer_nelx(iy) * m_layer_nelz(iy);
                 ++j;
             }
         }
         // -- top node layer
-        for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
-            ret(j) = m_startNode(iy + 1) + ix + (m_nelx(iy) + 1) * m_nelz(iy);
+        for (size_t ix = 0; ix < m_layer_nelx(iy) + 1; ++ix) {
+            ret(j) = m_startNode(iy + 1) + ix + (m_layer_nelx(iy) + 1) * m_layer_nelz(iy);
             ++j;
         }
     }
 
     return ret;
 }
 
 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;
+            n += m_layer_nelz(iy) * 3 + 1;
         }
         else {
-            n += m_nelz(iy) + 1;
+            n += m_layer_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;
+            n += m_layer_nelz(iy) * 3 + 1;
         }
         else {
-            n += m_nelz(iy) + 1;
+            n += m_layer_nelz(iy) + 1;
         }
     }
 
     // allocate node-list
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({n});
 
     // initialize counter: current index in the node-list "ret"
     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) {
-            ret(j) = m_startNode(iy) + iz * (m_nelx(iy) + 1);
+        for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) {
+            ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1);
             ++j;
         }
         // -- refinement layer
         if (m_refine(iy) == 2) {
-            for (size_t iz = 0; iz < 2 * m_nelz(iy); ++iz) {
-                ret(j) = m_startNode(iy) + iz * (m_nelx(iy) + 1) + m_nnd(iy);
+            for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) {
+                ret(j) = m_startNode(iy) + iz * (m_layer_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) {
-                ret(j) = m_startNode(iy) + iz * (m_nelx(iy) + 1) + m_nnd(iy);
+            for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) {
+                ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1) + m_nnd(iy);
                 ++j;
             }
         }
         // -- top node layer
-        for (size_t iz = 0; iz < m_nelz(iy) + 1; ++iz) {
-            ret(j) = m_startNode(iy + 1) + iz * (m_nelx(iy) + 1);
+        for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) {
+            ret(j) = m_startNode(iy + 1) + iz * (m_layer_nelx(iy) + 1);
             ++j;
         }
     }
 
     return ret;
 }
 
 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;
+            n += m_layer_nelz(iy) * 3 + 1;
         else
-            n += m_nelz(iy) + 1;
+            n += m_layer_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;
+            n += m_layer_nelz(iy) * 3 + 1;
         else
-            n += m_nelz(iy) + 1;
+            n += m_layer_nelz(iy) + 1;
     }
 
     // allocate node-list
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({n});
 
     // initialize counter: current index in the node-list "ret"
     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) {
-            ret(j) = m_startNode(iy) + iz * (m_nelx(iy) + 1) + m_nelx(iy);
+        for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) {
+            ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy);
             ++j;
         }
         // -- refinement layer
         if (m_refine(iy) == 2) {
-            for (size_t iz = 0; iz < 2 * m_nelz(iy); ++iz) {
-                ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy) + 1) + m_nelx(iy);
+            for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) {
+                ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_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) {
-                ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy) + 1) + m_nelx(iy);
+            for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) {
+                ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy);
                 ++j;
             }
         }
         // -- top node layer
-        for (size_t iz = 0; iz < m_nelz(iy) + 1; ++iz) {
-            ret(j) = m_startNode(iy + 1) + iz * (m_nelx(iy) + 1) + m_nelx(iy);
+        for (size_t iz = 0; iz < m_layer_nelz(iy) + 1; ++iz) {
+            ret(j) = m_startNode(iy + 1) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy);
             ++j;
         }
     }
 
     return ret;
 }
 
 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> ret = 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) {
-            ret(j) = m_startNode(0) + ix + iz * (m_nelx(0) + 1);
+    for (size_t ix = 0; ix < m_layer_nelx(0) + 1; ++ix) {
+        for (size_t iz = 0; iz < m_layer_nelz(0) + 1; ++iz) {
+            ret(j) = m_startNode(0) + ix + iz * (m_layer_nelx(0) + 1);
             ++j;
         }
     }
 
     return ret;
 }
 
 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> ret = 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) {
-            ret(j) = m_startNode(nely) + ix + iz * (m_nelx(nely - 1) + 1);
+    for (size_t ix = 0; ix < m_layer_nelx(nely - 1) + 1; ++ix) {
+        for (size_t iz = 0; iz < m_layer_nelz(nely - 1) + 1; ++iz) {
+            ret(j) = m_startNode(nely) + ix + iz * (m_layer_nelx(nely - 1) + 1);
             ++j;
         }
     }
 
     return ret;
 }
 
 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;
+            n += m_layer_nelx(iy) * 3 - 1;
         }
         else {
-            n += m_nelx(iy) - 1;
+            n += m_layer_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;
+            n += m_layer_nelx(iy) * 3 - 1;
         }
         else {
-            n += m_nelx(iy) - 1;
+            n += m_layer_nelx(iy) - 1;
         }
     }
 
     // allocate node-list
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({n});
 
     // initialize counter: current index in the node-list "ret"
     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) {
+        for (size_t ix = 1; ix < m_layer_nelx(iy); ++ix) {
             ret(j) = m_startNode(iy) + ix;
             ++j;
         }
         // -- refinement layer
         if (m_refine(iy) == 0) {
-            for (size_t ix = 0; ix < 2 * m_nelx(iy); ++ix) {
+            for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) {
                 ret(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) {
+            for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) {
                 ret(j) = m_startNode(iy) + ix + m_nnd(iy);
                 ++j;
             }
         }
         // -- top node layer
-        for (size_t ix = 1; ix < m_nelx(iy); ++ix) {
+        for (size_t ix = 1; ix < m_layer_nelx(iy); ++ix) {
             ret(j) = m_startNode(iy + 1) + ix;
             ++j;
         }
     }
 
     return ret;
 }
 
 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;
+            n += m_layer_nelx(iy) * 3 - 1;
         }
         else {
-            n += m_nelx(iy) - 1;
+            n += m_layer_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;
+            n += m_layer_nelx(iy) * 3 - 1;
         }
         else {
-            n += m_nelx(iy) - 1;
+            n += m_layer_nelx(iy) - 1;
         }
     }
 
     // allocate node-list
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({n});
 
     // initialize counter: current index in the node-list "ret"
     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) {
-            ret(j) = m_startNode(iy) + ix + (m_nelx(iy) + 1) * m_nelz(iy);
+        for (size_t ix = 1; ix < m_layer_nelx(iy); ++ix) {
+            ret(j) = m_startNode(iy) + ix + (m_layer_nelx(iy) + 1) * m_layer_nelz(iy);
             ++j;
         }
         // -- refinement layer
         if (m_refine(iy) == 0) {
-            for (size_t ix = 0; ix < 2 * m_nelx(iy); ++ix) {
-                ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_nelx(iy) * m_nelz(iy);
+            for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) {
+                ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_layer_nelx(iy) * m_layer_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) {
-                ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_nelx(iy) * m_nelz(iy);
+            for (size_t ix = 0; ix < 2 * m_layer_nelx(iy); ++ix) {
+                ret(j) = m_startNode(iy) + ix + m_nnd(iy) + 2 * m_layer_nelx(iy) * m_layer_nelz(iy);
                 ++j;
             }
         }
         // -- top node layer
-        for (size_t ix = 1; ix < m_nelx(iy); ++ix) {
-            ret(j) = m_startNode(iy + 1) + ix + (m_nelx(iy) + 1) * m_nelz(iy);
+        for (size_t ix = 1; ix < m_layer_nelx(iy); ++ix) {
+            ret(j) = m_startNode(iy + 1) + ix + (m_layer_nelx(iy) + 1) * m_layer_nelz(iy);
             ++j;
         }
     }
 
     return ret;
 }
 
 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;
+            n += m_layer_nelz(iy) * 3 - 1;
         }
         else {
-            n += m_nelz(iy) - 1;
+            n += m_layer_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;
+            n += m_layer_nelz(iy) * 3 - 1;
         }
         else {
-            n += m_nelz(iy) - 1;
+            n += m_layer_nelz(iy) - 1;
         }
     }
 
     // allocate node-list
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({n});
 
     // initialize counter: current index in the node-list "ret"
     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) {
-            ret(j) = m_startNode(iy) + iz * (m_nelx(iy) + 1);
+        for (size_t iz = 1; iz < m_layer_nelz(iy); ++iz) {
+            ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1);
             ++j;
         }
         // -- refinement layer
         if (m_refine(iy) == 2) {
-            for (size_t iz = 0; iz < 2 * m_nelz(iy); ++iz) {
-                ret(j) = m_startNode(iy) + iz * (m_nelx(iy) + 1) + m_nnd(iy);
+            for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) {
+                ret(j) = m_startNode(iy) + iz * (m_layer_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) {
-                ret(j) = m_startNode(iy) + iz * (m_nelx(iy) + 1) + m_nnd(iy);
+            for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) {
+                ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1) + m_nnd(iy);
                 ++j;
             }
         }
         // -- top node layer
-        for (size_t iz = 1; iz < m_nelz(iy); ++iz) {
-            ret(j) = m_startNode(iy + 1) + iz * (m_nelx(iy) + 1);
+        for (size_t iz = 1; iz < m_layer_nelz(iy); ++iz) {
+            ret(j) = m_startNode(iy + 1) + iz * (m_layer_nelx(iy) + 1);
             ++j;
         }
     }
 
     return ret;
 }
 
 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;
+            n += m_layer_nelz(iy) * 3 - 1;
         }
         else {
-            n += m_nelz(iy) - 1;
+            n += m_layer_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;
+            n += m_layer_nelz(iy) * 3 - 1;
         }
         else {
-            n += m_nelz(iy) - 1;
+            n += m_layer_nelz(iy) - 1;
         }
     }
 
     // allocate node-list
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({n});
 
     // initialize counter: current index in the node-list "ret"
     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) {
-            ret(j) = m_startNode(iy) + iz * (m_nelx(iy) + 1) + m_nelx(iy);
+        for (size_t iz = 1; iz < m_layer_nelz(iy); ++iz) {
+            ret(j) = m_startNode(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy);
             ++j;
         }
         // -- refinement layer
         if (m_refine(iy) == 2) {
-            for (size_t iz = 0; iz < 2 * m_nelz(iy); ++iz) {
-                ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy) + 1) + m_nelx(iy);
+            for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) {
+                ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_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) {
-                ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_nelx(iy) + 1) + m_nelx(iy);
+            for (size_t iz = 0; iz < 2 * m_layer_nelz(iy); ++iz) {
+                ret(j) = m_startNode(iy) + m_nnd(iy) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy);
                 ++j;
             }
         }
         // -- top node layer
-        for (size_t iz = 1; iz < m_nelz(iy); ++iz) {
-            ret(j) = m_startNode(iy + 1) + iz * (m_nelx(iy) + 1) + m_nelx(iy);
+        for (size_t iz = 1; iz < m_layer_nelz(iy); ++iz) {
+            ret(j) = m_startNode(iy + 1) + iz * (m_layer_nelx(iy) + 1) + m_layer_nelx(iy);
             ++j;
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBottomFace() const
 {
     // allocate node list
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({(m_nelx(0) - 1) * (m_nelz(0) - 1)});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({(m_layer_nelx(0) - 1) * (m_layer_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) {
-            ret(j) = m_startNode(0) + ix + iz * (m_nelx(0) + 1);
+    for (size_t ix = 1; ix < m_layer_nelx(0); ++ix) {
+        for (size_t iz = 1; iz < m_layer_nelz(0); ++iz) {
+            ret(j) = m_startNode(0) + ix + iz * (m_layer_nelx(0) + 1);
             ++j;
         }
     }
 
     return ret;
 }
 
 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> ret =
-        xt::empty<size_t>({(m_nelx(nely - 1) - 1) * (m_nelz(nely - 1) - 1)});
+        xt::empty<size_t>({(m_layer_nelx(nely - 1) - 1) * (m_layer_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) {
-            ret(j) = m_startNode(nely) + ix + iz * (m_nelx(nely - 1) + 1);
+    for (size_t ix = 1; ix < m_layer_nelx(nely - 1); ++ix) {
+        for (size_t iz = 1; iz < m_layer_nelz(nely - 1); ++iz) {
+            ret(j) = m_startNode(nely) + ix + iz * (m_layer_nelx(nely - 1) + 1);
             ++j;
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesFrontBottomEdge() const
 {
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx(0) + 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(0) + 1});
 
-    for (size_t ix = 0; ix < m_nelx(0) + 1; ++ix) {
+    for (size_t ix = 0; ix < m_layer_nelx(0) + 1; ++ix) {
         ret(ix) = m_startNode(0) + ix;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesFrontTopEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx(nely - 1) + 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(nely - 1) + 1});
 
-    for (size_t ix = 0; ix < m_nelx(nely - 1) + 1; ++ix) {
+    for (size_t ix = 0; ix < m_layer_nelx(nely - 1) + 1; ++ix) {
         ret(ix) = m_startNode(nely) + ix;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesFrontLeftEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
 
     for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) {
         ret(iy) = m_startNode(iy);
     }
 
     for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
         ret(iy + 1) = m_startNode(iy + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesFrontRightEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
 
     for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) {
-        ret(iy) = m_startNode(iy) + m_nelx(iy);
+        ret(iy) = m_startNode(iy) + m_layer_nelx(iy);
     }
 
     for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
-        ret(iy + 1) = m_startNode(iy + 1) + m_nelx(iy);
+        ret(iy + 1) = m_startNode(iy + 1) + m_layer_nelx(iy);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackBottomEdge() const
 {
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx(0) + 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(0) + 1});
 
-    for (size_t ix = 0; ix < m_nelx(0) + 1; ++ix) {
-        ret(ix) = m_startNode(0) + ix + (m_nelx(0) + 1) * (m_nelz(0));
+    for (size_t ix = 0; ix < m_layer_nelx(0) + 1; ++ix) {
+        ret(ix) = m_startNode(0) + ix + (m_layer_nelx(0) + 1) * (m_layer_nelz(0));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackTopEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx(nely - 1) + 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(nely - 1) + 1});
 
-    for (size_t ix = 0; ix < m_nelx(nely - 1) + 1; ++ix) {
-        ret(ix) = m_startNode(nely) + ix + (m_nelx(nely - 1) + 1) * (m_nelz(nely - 1));
+    for (size_t ix = 0; ix < m_layer_nelx(nely - 1) + 1; ++ix) {
+        ret(ix) = m_startNode(nely) + ix + (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackLeftEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
 
     for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) {
-        ret(iy) = m_startNode(iy) + (m_nelx(iy) + 1) * (m_nelz(iy));
+        ret(iy) = m_startNode(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
-        ret(iy + 1) = m_startNode(iy + 1) + (m_nelx(iy) + 1) * (m_nelz(iy));
+        ret(iy + 1) = m_startNode(iy + 1) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackRightEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
 
     for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) {
-        ret(iy) = m_startNode(iy) + m_nelx(iy) + (m_nelx(iy) + 1) * (m_nelz(iy));
+        ret(iy) = m_startNode(iy) + m_layer_nelx(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
-        ret(iy + 1) = m_startNode(iy + 1) + m_nelx(iy) + (m_nelx(iy) + 1) * (m_nelz(iy));
+        ret(iy + 1) = m_startNode(iy + 1) + m_layer_nelx(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBottomLeftEdge() const
 {
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz(0) + 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(0) + 1});
 
-    for (size_t iz = 0; iz < m_nelz(0) + 1; ++iz) {
-        ret(iz) = m_startNode(0) + iz * (m_nelx(0) + 1);
+    for (size_t iz = 0; iz < m_layer_nelz(0) + 1; ++iz) {
+        ret(iz) = m_startNode(0) + iz * (m_layer_nelx(0) + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBottomRightEdge() const
 {
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz(0) + 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(0) + 1});
 
-    for (size_t iz = 0; iz < m_nelz(0) + 1; ++iz) {
-        ret(iz) = m_startNode(0) + m_nelx(0) + iz * (m_nelx(0) + 1);
+    for (size_t iz = 0; iz < m_layer_nelz(0) + 1; ++iz) {
+        ret(iz) = m_startNode(0) + m_layer_nelx(0) + iz * (m_layer_nelx(0) + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesTopLeftEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz(nely - 1) + 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(nely - 1) + 1});
 
-    for (size_t iz = 0; iz < m_nelz(nely - 1) + 1; ++iz) {
-        ret(iz) = m_startNode(nely) + iz * (m_nelx(nely - 1) + 1);
+    for (size_t iz = 0; iz < m_layer_nelz(nely - 1) + 1; ++iz) {
+        ret(iz) = m_startNode(nely) + iz * (m_layer_nelx(nely - 1) + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesTopRightEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz(nely - 1) + 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(nely - 1) + 1});
 
-    for (size_t iz = 0; iz < m_nelz(nely - 1) + 1; ++iz) {
-        ret(iz) = m_startNode(nely) + m_nelx(nely - 1) + iz * (m_nelx(nely - 1) + 1);
+    for (size_t iz = 0; iz < m_layer_nelz(nely - 1) + 1; ++iz) {
+        ret(iz) = m_startNode(nely) + m_layer_nelx(nely - 1) + iz * (m_layer_nelx(nely - 1) + 1);
     }
 
     return ret;
 }
 
 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();
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesFrontBottomOpenEdge() const
 {
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx(0) - 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(0) - 1});
 
-    for (size_t ix = 1; ix < m_nelx(0); ++ix) {
+    for (size_t ix = 1; ix < m_layer_nelx(0); ++ix) {
         ret(ix - 1) = m_startNode(0) + ix;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesFrontTopOpenEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx(nely - 1) - 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(nely - 1) - 1});
 
-    for (size_t ix = 1; ix < m_nelx(nely - 1); ++ix) {
+    for (size_t ix = 1; ix < m_layer_nelx(nely - 1); ++ix) {
         ret(ix - 1) = m_startNode(nely) + ix;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesFrontLeftOpenEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
 
     for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) {
         ret(iy - 1) = m_startNode(iy);
     }
 
     for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) {
         ret(iy) = m_startNode(iy + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesFrontRightOpenEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
 
     for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) {
-        ret(iy - 1) = m_startNode(iy) + m_nelx(iy);
+        ret(iy - 1) = m_startNode(iy) + m_layer_nelx(iy);
     }
 
     for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) {
-        ret(iy) = m_startNode(iy + 1) + m_nelx(iy);
+        ret(iy) = m_startNode(iy + 1) + m_layer_nelx(iy);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackBottomOpenEdge() const
 {
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx(0) - 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(0) - 1});
 
-    for (size_t ix = 1; ix < m_nelx(0); ++ix) {
-        ret(ix - 1) = m_startNode(0) + ix + (m_nelx(0) + 1) * (m_nelz(0));
+    for (size_t ix = 1; ix < m_layer_nelx(0); ++ix) {
+        ret(ix - 1) = m_startNode(0) + ix + (m_layer_nelx(0) + 1) * (m_layer_nelz(0));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackTopOpenEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelx(nely - 1) - 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(nely - 1) - 1});
 
-    for (size_t ix = 1; ix < m_nelx(nely - 1); ++ix) {
-        ret(ix - 1) = m_startNode(nely) + ix + (m_nelx(nely - 1) + 1) * (m_nelz(nely - 1));
+    for (size_t ix = 1; ix < m_layer_nelx(nely - 1); ++ix) {
+        ret(ix - 1) = m_startNode(nely) + ix + (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackLeftOpenEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
 
     for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) {
-        ret(iy - 1) = m_startNode(iy) + (m_nelx(iy) + 1) * (m_nelz(iy));
+        ret(iy - 1) = m_startNode(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) {
-        ret(iy) = m_startNode(iy + 1) + (m_nelx(iy) + 1) * (m_nelz(iy));
+        ret(iy) = m_startNode(iy + 1) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackRightOpenEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
 
     for (size_t iy = 1; iy < (nely + 1) / 2; ++iy) {
-        ret(iy - 1) = m_startNode(iy) + m_nelx(iy) + (m_nelx(iy) + 1) * (m_nelz(iy));
+        ret(iy - 1) = m_startNode(iy) + m_layer_nelx(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) {
-        ret(iy) = m_startNode(iy + 1) + m_nelx(iy) + (m_nelx(iy) + 1) * (m_nelz(iy));
+        ret(iy) = m_startNode(iy + 1) + m_layer_nelx(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBottomLeftOpenEdge() const
 {
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz(0) - 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(0) - 1});
 
-    for (size_t iz = 1; iz < m_nelz(0); ++iz) {
-        ret(iz - 1) = m_startNode(0) + iz * (m_nelx(0) + 1);
+    for (size_t iz = 1; iz < m_layer_nelz(0); ++iz) {
+        ret(iz - 1) = m_startNode(0) + iz * (m_layer_nelx(0) + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBottomRightOpenEdge() const
 {
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz(0) - 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(0) - 1});
 
-    for (size_t iz = 1; iz < m_nelz(0); ++iz) {
-        ret(iz - 1) = m_startNode(0) + m_nelx(0) + iz * (m_nelx(0) + 1);
+    for (size_t iz = 1; iz < m_layer_nelz(0); ++iz) {
+        ret(iz - 1) = m_startNode(0) + m_layer_nelx(0) + iz * (m_layer_nelx(0) + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesTopLeftOpenEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz(nely - 1) - 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(nely - 1) - 1});
 
-    for (size_t iz = 1; iz < m_nelz(nely - 1); ++iz) {
-        ret(iz - 1) = m_startNode(nely) + iz * (m_nelx(nely - 1) + 1);
+    for (size_t iz = 1; iz < m_layer_nelz(nely - 1); ++iz) {
+        ret(iz - 1) = m_startNode(nely) + iz * (m_layer_nelx(nely - 1) + 1);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesTopRightOpenEdge() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
-    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelz(nely - 1) - 1});
+    xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(nely - 1) - 1});
 
-    for (size_t iz = 1; iz < m_nelz(nely - 1); ++iz) {
-        ret(iz - 1) = m_startNode(nely) + m_nelx(nely - 1) + iz * (m_nelx(nely - 1) + 1);
+    for (size_t iz = 1; iz < m_layer_nelz(nely - 1); ++iz) {
+        ret(iz - 1) = m_startNode(nely) + m_layer_nelx(nely - 1) + iz * (m_layer_nelx(nely - 1) + 1);
     }
 
     return ret;
 }
 
 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();
 }
 
 inline size_t FineLayer::nodesFrontBottomLeftCorner() const
 {
     return m_startNode(0);
 }
 
 inline size_t FineLayer::nodesFrontBottomRightCorner() const
 {
-    return m_startNode(0) + m_nelx(0);
+    return m_startNode(0) + m_layer_nelx(0);
 }
 
 inline size_t FineLayer::nodesFrontTopLeftCorner() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
     return m_startNode(nely);
 }
 
 inline size_t FineLayer::nodesFrontTopRightCorner() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
-    return m_startNode(nely) + m_nelx(nely - 1);
+    return m_startNode(nely) + m_layer_nelx(nely - 1);
 }
 
 inline size_t FineLayer::nodesBackBottomLeftCorner() const
 {
-    return m_startNode(0) + (m_nelx(0) + 1) * (m_nelz(0));
+    return m_startNode(0) + (m_layer_nelx(0) + 1) * (m_layer_nelz(0));
 }
 
 inline size_t FineLayer::nodesBackBottomRightCorner() const
 {
-    return m_startNode(0) + m_nelx(0) + (m_nelx(0) + 1) * (m_nelz(0));
+    return m_startNode(0) + m_layer_nelx(0) + (m_layer_nelx(0) + 1) * (m_layer_nelz(0));
 }
 
 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));
+    return m_startNode(nely) + (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1));
 }
 
 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));
+    return m_startNode(nely) + m_layer_nelx(nely - 1) + (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1));
 }
 
 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();
 }
 
 inline xt::xtensor<size_t, 2> FineLayer::nodesPeriodic() const
 {
     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();
 
     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();
 
     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;
     xt::xtensor<size_t, 2> ret = xt::empty<size_t>({tface + tedge + tnode, std::size_t(2)});
 
     size_t i = 0;
 
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesFrontBottomRightCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesBackBottomRightCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesBackBottomLeftCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesFrontTopLeftCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesFrontTopRightCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesBackTopRightCorner();
     ++i;
     ret(i, 0) = nodesFrontBottomLeftCorner();
     ret(i, 1) = nodesBackTopLeftCorner();
     ++i;
 
     for (size_t j = 0; j < froBot.size(); ++j) {
         ret(i, 0) = froBot(j);
         ret(i, 1) = bckBot(j);
         ++i;
     }
     for (size_t j = 0; j < froBot.size(); ++j) {
         ret(i, 0) = froBot(j);
         ret(i, 1) = bckTop(j);
         ++i;
     }
     for (size_t j = 0; j < froBot.size(); ++j) {
         ret(i, 0) = froBot(j);
         ret(i, 1) = froTop(j);
         ++i;
     }
     for (size_t j = 0; j < botLft.size(); ++j) {
         ret(i, 0) = botLft(j);
         ret(i, 1) = botRgt(j);
         ++i;
     }
     for (size_t j = 0; j < botLft.size(); ++j) {
         ret(i, 0) = botLft(j);
         ret(i, 1) = topRgt(j);
         ++i;
     }
     for (size_t j = 0; j < botLft.size(); ++j) {
         ret(i, 0) = botLft(j);
         ret(i, 1) = topLft(j);
         ++i;
     }
     for (size_t j = 0; j < froLft.size(); ++j) {
         ret(i, 0) = froLft(j);
         ret(i, 1) = froRgt(j);
         ++i;
     }
     for (size_t j = 0; j < froLft.size(); ++j) {
         ret(i, 0) = froLft(j);
         ret(i, 1) = bckRgt(j);
         ++i;
     }
     for (size_t j = 0; j < froLft.size(); ++j) {
         ret(i, 0) = froLft(j);
         ret(i, 1) = bckLft(j);
         ++i;
     }
 
     for (size_t j = 0; j < fro.size(); ++j) {
         ret(i, 0) = fro(j);
         ret(i, 1) = bck(j);
         ++i;
     }
     for (size_t j = 0; j < lft.size(); ++j) {
         ret(i, 0) = lft(j);
         ret(i, 1) = rgt(j);
         ++i;
     }
     for (size_t j = 0; j < bot.size(); ++j) {
         ret(i, 0) = bot(j);
         ret(i, 1) = top(j);
         ++i;
     }
 
     return ret;
 }
 
 inline size_t FineLayer::nodesOrigin() const
 {
     return nodesFrontBottomLeftCorner();
 }
 
 inline xt::xtensor<size_t, 2> FineLayer::dofs() const
 {
     return GooseFEM::Mesh::dofs(m_nnode, m_ndim);
 }
 
 inline xt::xtensor<size_t, 2> FineLayer::dofsPeriodic() const
 {
     xt::xtensor<size_t, 2> ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim);
     xt::xtensor<size_t, 2> nodePer = nodesPeriodic();
 
     for (size_t i = 0; i < nodePer.shape(0); ++i) {
         for (size_t j = 0; j < m_ndim; ++j) {
             ret(nodePer(i, 1), j) = ret(nodePer(i, 0), j);
         }
     }
 
     return GooseFEM::Mesh::renumber(ret);
 }
 
 } // namespace Hex8
 } // namespace Mesh
 } // namespace GooseFEM
 
 #endif