diff --git a/include/GooseFEM/MeshHex8.h b/include/GooseFEM/MeshHex8.h
index ec90df7..1857fbf 100644
--- a/include/GooseFEM/MeshHex8.h
+++ b/include/GooseFEM/MeshHex8.h
@@ -1,210 +1,213 @@
 /**
 Generate simple meshes of 8-noded hexahedral elements in 3d (GooseFEM::Mesh::ElementType::Hex8).
 
 \file MeshHex8.h
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_MESHHEX8_H
 #define GOOSEFEM_MESHHEX8_H
 
 #include "config.h"
 #include "Mesh.h"
 
 namespace GooseFEM {
 namespace Mesh {
 
 /**
 Simple meshes of 8-noded hexahedral elements in 3d (ElementType::Hex8).
 */
 namespace Hex8 {
 
 /**
 Regular mesh: equi-sized elements.
 */
 class Regular : public RegularBase3d<Regular> {
 public:
 
     /**
     Constructor.
 
     \param nelx Number of elements in horizontal (x) direction.
     \param nely Number of elements in vertical (y) direction.
     \param nelz Number of elements in vertical (z) direction.
     \param h Edge size (width == height == depth).
     */
     Regular(size_t nelx, size_t nely, size_t nelz, double h = 1.0);
 
 private:
     friend class RegularBase<Regular>;
     friend class RegularBase3d<Regular>;
 
+    size_t nelx_impl() const;
+    size_t nely_impl() const;
+    size_t nelz_impl() const;
     ElementType getElementType_impl() const;
     xt::xtensor<double, 2> coor_impl() const;
     xt::xtensor<size_t, 2> conn_impl() const;
     xt::xtensor<size_t, 1> nodesFront_impl() const;
     xt::xtensor<size_t, 1> nodesBack_impl() const;
     xt::xtensor<size_t, 1> nodesLeft_impl() const;
     xt::xtensor<size_t, 1> nodesRight_impl() const;
     xt::xtensor<size_t, 1> nodesBottom_impl() const;
     xt::xtensor<size_t, 1> nodesTop_impl() const;
     xt::xtensor<size_t, 1> nodesFrontFace_impl() const;
     xt::xtensor<size_t, 1> nodesBackFace_impl() const;
     xt::xtensor<size_t, 1> nodesLeftFace_impl() const;
     xt::xtensor<size_t, 1> nodesRightFace_impl() const;
     xt::xtensor<size_t, 1> nodesBottomFace_impl() const;
     xt::xtensor<size_t, 1> nodesTopFace_impl() const;
     xt::xtensor<size_t, 1> nodesFrontBottomEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontTopEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontLeftEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontRightEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackBottomEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackTopEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackLeftEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackRightEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBottomLeftEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBottomRightEdge_impl() const;
     xt::xtensor<size_t, 1> nodesTopLeftEdge_impl() const;
     xt::xtensor<size_t, 1> nodesTopRightEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontBottomOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontTopOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontLeftOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontRightOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackBottomOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackTopOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackLeftOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackRightOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBottomLeftOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBottomRightOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesTopLeftOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesTopRightOpenEdge_impl() const;
     size_t nodesFrontBottomLeftCorner_impl() const;
     size_t nodesFrontBottomRightCorner_impl() const;
     size_t nodesFrontTopLeftCorner_impl() const;
     size_t nodesFrontTopRightCorner_impl() const;
     size_t nodesBackBottomLeftCorner_impl() const;
     size_t nodesBackBottomRightCorner_impl() const;
     size_t nodesBackTopLeftCorner_impl() const;
     size_t nodesBackTopRightCorner_impl() const;
 
     double m_h;     ///< See h()
     size_t m_nelx;  ///< See nelx()
     size_t m_nely;  ///< See nely()
     size_t m_nelz;  ///< See nely()
     size_t m_nelem; ///< See nelem()
     size_t m_nnode; ///< See nnode()
     size_t m_nne;   ///< See nne()
     size_t m_ndim;  ///< See ndim()
 };
 
 /**
 Mesh with fine middle layer, and coarser elements towards the top and bottom.
 */
 class FineLayer : public RegularBase3d<FineLayer> {
 public:
 
     /**
     Constructor.
 
     \param nelx Number of elements (along the middle layer) in horizontal (x) direction.
     \param nely Approximate equivalent number of elements in vertical (y) direction.
     \param nelz Number of elements (along the middle layer) in depth (z) direction.
     \param h Edge size (width == height == depth) of elements along the weak layer.
 
     \param nfine
         Extra number of fine layers around the middle layer.
         By default the element size is kept smaller than the distance to the middle layer.
     */
     FineLayer(size_t nelx, size_t nely, size_t nelz, double h = 1.0, size_t nfine = 1);
 
     /**
     Elements in the middle (fine) layer.
 
     \return List of element numbers.
     */
     xt::xtensor<size_t, 1> elementsMiddleLayer() const;
 
 private:
     friend class RegularBase<FineLayer>;
     friend class RegularBase3d<FineLayer>;
 
     size_t nelx_impl() const;
     size_t nely_impl() const;
     size_t nelz_impl() const;
     ElementType getElementType_impl() const;
     xt::xtensor<double, 2> coor_impl() const;
     xt::xtensor<size_t, 2> conn_impl() const;
     xt::xtensor<size_t, 1> nodesFront_impl() const;
     xt::xtensor<size_t, 1> nodesBack_impl() const;
     xt::xtensor<size_t, 1> nodesLeft_impl() const;
     xt::xtensor<size_t, 1> nodesRight_impl() const;
     xt::xtensor<size_t, 1> nodesBottom_impl() const;
     xt::xtensor<size_t, 1> nodesTop_impl() const;
     xt::xtensor<size_t, 1> nodesFrontFace_impl() const;
     xt::xtensor<size_t, 1> nodesBackFace_impl() const;
     xt::xtensor<size_t, 1> nodesLeftFace_impl() const;
     xt::xtensor<size_t, 1> nodesRightFace_impl() const;
     xt::xtensor<size_t, 1> nodesBottomFace_impl() const;
     xt::xtensor<size_t, 1> nodesTopFace_impl() const;
     xt::xtensor<size_t, 1> nodesFrontBottomEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontTopEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontLeftEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontRightEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackBottomEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackTopEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackLeftEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackRightEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBottomLeftEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBottomRightEdge_impl() const;
     xt::xtensor<size_t, 1> nodesTopLeftEdge_impl() const;
     xt::xtensor<size_t, 1> nodesTopRightEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontBottomOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontTopOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontLeftOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesFrontRightOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackBottomOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackTopOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackLeftOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBackRightOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBottomLeftOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesBottomRightOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesTopLeftOpenEdge_impl() const;
     xt::xtensor<size_t, 1> nodesTopRightOpenEdge_impl() const;
     size_t nodesFrontBottomLeftCorner_impl() const;
     size_t nodesFrontBottomRightCorner_impl() const;
     size_t nodesFrontTopLeftCorner_impl() const;
     size_t nodesFrontTopRightCorner_impl() const;
     size_t nodesBackBottomLeftCorner_impl() const;
     size_t nodesBackBottomRightCorner_impl() const;
     size_t nodesBackTopLeftCorner_impl() const;
     size_t nodesBackTopRightCorner_impl() const;
 
     double m_h;                          ///< See h()
     size_t m_nelx;                       ///< See nelx()
     size_t m_nely;                       ///< See nely()
     size_t m_nelz;                       ///< See nely()
     size_t m_nelem;                      ///< See nelem()
     size_t m_nnode;                      ///< See nnode()
     size_t m_nne;                        ///< See nne()
     size_t m_ndim;                       ///< See ndim()
     double m_Lx;                         ///< mesh size in "x"
     double m_Lz;                         ///< mesh size in "z"
     xt::xtensor<size_t, 1> m_layer_nelx; ///< number of elements in "x" per element layer in "y"
     xt::xtensor<size_t, 1> m_layer_nelz; ///< number of elements in "z" per element layer in "y"
     xt::xtensor<size_t, 1> m_nnd;        ///< number of nodes in the main node layer per node layer in "y"
     xt::xtensor<size_t, 1> m_nhx;        ///< element size in x-direction per element layer in "y"
     xt::xtensor<size_t, 1> m_nhy;        ///< element size in y-direction per element layer in "y"
     xt::xtensor<size_t, 1> m_nhz;        ///< element size in z-direction per element layer in "y"
     xt::xtensor<int, 1> m_refine;        ///< refine direction (-1:no refine, 0:"x", 2:"z") per element layer in "y"
     xt::xtensor<size_t, 1> m_startElem;  ///< start element per element layer in "y"
     xt::xtensor<size_t, 1> m_startNode;  ///< start node per node layer in "y"
 };
 
 } // namespace Hex8
 } // namespace Mesh
 } // namespace GooseFEM
 
 #include "MeshHex8.hpp"
 
 #endif
diff --git a/include/GooseFEM/MeshHex8.hpp b/include/GooseFEM/MeshHex8.hpp
index 4334a5c..0dad188 100644
--- a/include/GooseFEM/MeshHex8.hpp
+++ b/include/GooseFEM/MeshHex8.hpp
@@ -1,2189 +1,2200 @@
 /**
 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;
     m_nne = 8;
     m_ndim = 3;
 
     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 ElementType Regular::getElementType_impl() const
 {
     return ElementType::Hex8;
 }
 
+inline size_t Regular::nelx_impl() const
+{
+    return m_nelx;
+}
+
+inline size_t Regular::nely_impl() const
+{
+    return m_nely;
+}
+
+inline size_t Regular::nelz_impl() const
+{
+    return m_nelz;
+}
+
 inline xt::xtensor<double, 2> Regular::coor_impl() 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_impl() 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, 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);
+                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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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::nodesFrontBottomOpenEdge_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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_impl() 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 size_t Regular::nodesFrontBottomLeftCorner_impl() const
 {
     return 0;
 }
 
 inline size_t Regular::nodesFrontBottomRightCorner_impl() const
 {
     return m_nelx;
 }
 
 inline size_t Regular::nodesFrontTopLeftCorner_impl() const
 {
     return m_nely * (m_nelx + 1);
 }
 
 inline size_t Regular::nodesFrontTopRightCorner_impl() const
 {
     return m_nely * (m_nelx + 1) + m_nelx;
 }
 
 inline size_t Regular::nodesBackBottomLeftCorner_impl() const
 {
     return m_nelz * (m_nely + 1) * (m_nelx + 1);
 }
 
 inline size_t Regular::nodesBackBottomRightCorner_impl() const
 {
     return m_nelx + m_nelz * (m_nely + 1) * (m_nelx + 1);
 }
 
 inline size_t Regular::nodesBackTopLeftCorner_impl() const
 {
     return m_nely * (m_nelx + 1) + m_nelz * (m_nely + 1) * (m_nelx + 1);
 }
 
 inline size_t Regular::nodesBackTopRightCorner_impl() const
 {
     return m_nely * (m_nelx + 1) + m_nelx + m_nelz * (m_nely + 1) * (m_nelx + 1);
 }
 
 inline FineLayer::FineLayer(size_t nelx, size_t nely, size_t nelz, double h, size_t nfine)
 {
     m_h = h;
     m_nne = 8;
     m_ndim = 3;
 
     // 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_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_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_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_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_layer_nelx(i) + 1) * (m_layer_nelz(i) + 1);
         }
         else if (m_refine(i) == 2) {
             m_nnode += (m_layer_nelx(i) + 1) * (3 * m_layer_nelz(i) + 1);
         }
         else {
             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_layer_nelx(i)) * (m_layer_nelz(i));
         }
         else if (m_refine(i) == 2) {
             m_nelem += (m_layer_nelx(i)) * (4 * m_layer_nelz(i));
         }
         else {
             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_layer_nelx(i) + 1) * (m_layer_nelz(i) + 1);
         }
         else if (m_refine(i) == 2) {
             m_nnode += (m_layer_nelx(i) + 1) * (5 * m_layer_nelz(i) + 1);
         }
         else {
             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_layer_nelx(i)) * (m_layer_nelz(i));
         }
         else if (m_refine(i) == 2) {
             m_nelem += (m_layer_nelx(i)) * (4 * m_layer_nelz(i));
         }
         else {
             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_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1) + 1);
 }
 
 inline size_t FineLayer::nelx_impl() const
 {
     return xt::amax(m_layer_nelx)();
 }
 
 inline size_t FineLayer::nely_impl() const
 {
     return xt::sum(m_nhy)();
 }
 
 inline size_t FineLayer::nelz_impl() const
 {
     return xt::amax(m_layer_nelz)();
 }
 
 inline ElementType FineLayer::getElementType_impl() const
 {
     return ElementType::Hex8;
 }
 
 inline xt::xtensor<double, 2> FineLayer::coor_impl() 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_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_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_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_layer_nelz(iy); ++iz) {
                 for (size_t j = 0; j < 2; ++j) {
                     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_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_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_layer_nelz(iy); ++iz) {
                 for (size_t j = 0; j < 2; ++j) {
                     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_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_impl() 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_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_layer_nelz(iy); ++iz) {
                 for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
                     // -- bottom element
                     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_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_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_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_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_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_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_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_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_layer_nelz(iy); ++iz) {
                 for (size_t ix = 0; ix < m_layer_nelx(iy); ++ix) {
                     // -- bottom element
                     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_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_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_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_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_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_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_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_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_layer_nelx(iy) * m_layer_nelz(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_impl() 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_layer_nelx(iy) * 3 + 1;
         }
         else {
             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_layer_nelx(iy) * 3 + 1;
         }
         else {
             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_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_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_layer_nelx(iy); ++ix) {
                 ret(j) = m_startNode(iy) + ix + m_nnd(iy);
                 ++j;
             }
         }
         // -- top node layer
         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_impl() 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_layer_nelx(iy) * 3 + 1;
         }
         else {
             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_layer_nelx(iy) * 3 + 1;
         }
         else {
             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_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_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_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_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_impl() 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_layer_nelz(iy) * 3 + 1;
         }
         else {
             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_layer_nelz(iy) * 3 + 1;
         }
         else {
             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_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_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_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_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_impl() 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_layer_nelz(iy) * 3 + 1;
         else
             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_layer_nelz(iy) * 3 + 1;
         else
             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_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_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_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_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_impl() 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_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_impl() 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_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_impl() 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_layer_nelx(iy) * 3 - 1;
         }
         else {
             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_layer_nelx(iy) * 3 - 1;
         }
         else {
             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_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_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_layer_nelx(iy); ++ix) {
                 ret(j) = m_startNode(iy) + ix + m_nnd(iy);
                 ++j;
             }
         }
         // -- top node layer
         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_impl() 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_layer_nelx(iy) * 3 - 1;
         }
         else {
             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_layer_nelx(iy) * 3 - 1;
         }
         else {
             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_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_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_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_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_impl() 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_layer_nelz(iy) * 3 - 1;
         }
         else {
             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_layer_nelz(iy) * 3 - 1;
         }
         else {
             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_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_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_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_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_impl() 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_layer_nelz(iy) * 3 - 1;
         }
         else {
             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_layer_nelz(iy) * 3 - 1;
         }
         else {
             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_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_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_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_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_impl() const
 {
     // allocate node list
     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_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_impl() 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_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_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_impl() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(0) + 1});
 
     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_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(nely - 1) + 1});
 
     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_impl() 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_impl() 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_layer_nelx(iy);
     }
 
     for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
         ret(iy + 1) = m_startNode(iy + 1) + m_layer_nelx(iy);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackBottomEdge_impl() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(0) + 1});
 
     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_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(nely - 1) + 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_impl() 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_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_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackRightEdge_impl() 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_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_layer_nelx(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBottomLeftEdge_impl() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(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_impl() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(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_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(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_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(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::nodesFrontBottomOpenEdge_impl() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(0) - 1});
 
     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_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(nely - 1) - 1});
 
     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_impl() 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_impl() 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_layer_nelx(iy);
     }
 
     for (size_t iy = (nely - 1) / 2; iy < nely - 1; ++iy) {
         ret(iy) = m_startNode(iy + 1) + m_layer_nelx(iy);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackBottomOpenEdge_impl() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(0) - 1});
 
     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_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelx(nely - 1) - 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_impl() 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_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_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBackRightOpenEdge_impl() 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_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_layer_nelx(iy) + (m_layer_nelx(iy) + 1) * (m_layer_nelz(iy));
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBottomLeftOpenEdge_impl() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(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_impl() const
 {
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(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_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(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_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_layer_nelz(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 size_t FineLayer::nodesFrontBottomLeftCorner_impl() const
 {
     return m_startNode(0);
 }
 
 inline size_t FineLayer::nodesFrontBottomRightCorner_impl() const
 {
     return m_startNode(0) + m_layer_nelx(0);
 }
 
 inline size_t FineLayer::nodesFrontTopLeftCorner_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
     return m_startNode(nely);
 }
 
 inline size_t FineLayer::nodesFrontTopRightCorner_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
     return m_startNode(nely) + m_layer_nelx(nely - 1);
 }
 
 inline size_t FineLayer::nodesBackBottomLeftCorner_impl() const
 {
     return m_startNode(0) + (m_layer_nelx(0) + 1) * (m_layer_nelz(0));
 }
 
 inline size_t FineLayer::nodesBackBottomRightCorner_impl() const
 {
     return m_startNode(0) + m_layer_nelx(0) + (m_layer_nelx(0) + 1) * (m_layer_nelz(0));
 }
 
 inline size_t FineLayer::nodesBackTopLeftCorner_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
     return m_startNode(nely) + (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1));
 }
 
 inline size_t FineLayer::nodesBackTopRightCorner_impl() const
 {
     size_t nely = static_cast<size_t>(m_nhy.size());
     return m_startNode(nely) + m_layer_nelx(nely - 1) + (m_layer_nelx(nely - 1) + 1) * (m_layer_nelz(nely - 1));
 }
 
 } // namespace Hex8
 } // namespace Mesh
 } // namespace GooseFEM
 
 #endif
diff --git a/python/Element.hpp b/python/Element.hpp
index 355dfae..567b06a 100644
--- a/python/Element.hpp
+++ b/python/Element.hpp
@@ -1,194 +1,194 @@
 /**
 \file
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef PYGOOSEFEM_ELEMENT_H
 #define PYGOOSEFEM_ELEMENT_H
 
-#include <GooseFEM/GooseFEM.h>
+#include <GooseFEM/Element.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 #include <xtensor-python/pyarray.hpp>
 #include <xtensor-python/pytensor.hpp>
 #include <pyxtensor/pyxtensor.hpp>
 
 namespace py = pybind11;
 
 template <class C, class P>
 void register_Element_QuadratureBase(P& cls)
 {
     cls.def("nelem", &C::nelem, "Number of elements");
     cls.def("nne", &C::nne, "Number of nodes per element");
     cls.def("ndim", &C::ndim, "Number of dimensions");
     cls.def("tdim", &C::tdim, "Number of dimensions of tensors");
     cls.def("nip", &C::nip, "Number of integration points");
 
     // todo: https://github.com/xtensor-stack/xtensor-python/issues/265
     // cls.def("asTensor",
     //         static_cast<void (C::*)(const xt::pyarray<double>&, xt::pyarray<double>&) const>(&C::asTensor),
     //         "Convert 'qscalar' to 'qtensor' of certain rank",
     //         py::arg("qscalar"),
     //         py::arg("qtensor"));
 
     cls.def("AsTensor",
             static_cast<xt::xarray<double> (C::*)(size_t, const xt::pyarray<double>&) const>(&C::AsTensor),
             "Convert 'qscalar' to 'qtensor' of certain rank",
             py::arg("rank"),
             py::arg("qscalar"));
 
     cls.def("shape_elemvec",
             static_cast<std::array<size_t, 3> (C::*)() const>(&C::shape_elemvec),
             "Shape of 'elemvec'");
 
     cls.def("shape_elemvec",
             static_cast<std::array<size_t, 3> (C::*)(size_t) const>(&C::shape_elemvec),
             "Shape of 'elemvec'",
             py::arg("tdim"));
 
     cls.def("shape_elemmat",
             &C::shape_elemmat,
             "Shape of 'elemmat'");
 
     cls.def("shape_qtensor",
             static_cast<std::vector<size_t> (C::*)(size_t) const>(&C::shape_qtensor),
             "Shape of 'qtensor'",
             py::arg("rank"));
 
     cls.def("shape_qtensor",
             static_cast<std::vector<size_t> (C::*)(size_t, size_t) const>(&C::shape_qtensor),
             "Shape of 'qtensor'",
             py::arg("rank"),
             py::arg("tdim"));
 
     cls.def("shape_qscalar",
             &C::shape_qscalar,
             "Shape of 'qscalar'");
 
     cls.def("shape_qvector",
             static_cast<std::array<size_t, 3> (C::*)() const>(&C::shape_qvector),
             "Shape of 'qvector'");
 
     cls.def("shape_qvector",
             static_cast<std::array<size_t, 3> (C::*)(size_t) const>(&C::shape_qvector),
             "Shape of 'qvector'",
             py::arg("tdim"));
 }
 
 template <class C, class P>
 void register_Element_QuadratureBaseCartesian(P& cls)
 {
     cls.def("update_x",
             &C::template update_x<xt::pytensor<double, 3>>,
             "Update the nodal positions",
             py::arg("x"));
 
     cls.def("dV", &C::dV, "Integration point volume (qscalar)");
 
     cls.def("InterpQuad_vector",
             &C::template InterpQuad_vector<xt::pytensor<double, 3>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::InterpQuad_vector`.",
             py::arg("elemvec"));
 
     cls.def("interpQuad_vector",
             &C::template interpQuad_vector<xt::pytensor<double, 3>, xt::pytensor<double, 3>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::interpQuad_vector`.",
             py::arg("elemvec"),
             py::arg("qvector"));
 
     cls.def("GradN_vector",
             &C::template GradN_vector<xt::pytensor<double, 3>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::GradN_vector`.",
             py::arg("elemvec"));
 
     cls.def("gradN_vector",
             &C::template gradN_vector<xt::pytensor<double, 3>, xt::pytensor<double, 4>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::gradN_vector`.",
             py::arg("elemvec"),
             py::arg("qtensor"));
 
     cls.def("GradN_vector_T",
             &C::template GradN_vector_T<xt::pytensor<double, 3>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::GradN_vector_T`.",
             py::arg("elemvec"));
 
     cls.def("gradN_vector_T",
             &C::template gradN_vector_T<xt::pytensor<double, 3>, xt::pytensor<double, 4>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::gradN_vector_T`.",
             py::arg("elemvec"),
             py::arg("qtensor"));
 
     cls.def("SymGradN_vector",
             &C::template SymGradN_vector<xt::pytensor<double, 3>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::SymGradN_vector`.",
             py::arg("elemvec"));
 
     cls.def("symGradN_vector",
             &C::template symGradN_vector<xt::pytensor<double, 3>, xt::pytensor<double, 4>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::symGradN_vector`.",
             py::arg("elemvec"),
             py::arg("qtensor"));
 
     cls.def("Int_N_vector_dV",
             &C::template Int_N_vector_dV<xt::pytensor<double, 3>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::Int_N_vector_dV`.",
             py::arg("qvector"));
 
     cls.def("int_N_vector_dV",
             &C::template int_N_vector_dV<xt::pytensor<double, 3>, xt::pytensor<double, 3>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::int_N_vector_dV`.",
             py::arg("qvector"),
             py::arg("elemvec"));
 
     cls.def("Int_N_scalar_NT_dV",
             &C::template Int_N_scalar_NT_dV<xt::pytensor<double, 2>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::Int_N_scalar_NT_dV`.",
             py::arg("qscalar"));
 
     cls.def("int_N_scalar_NT_dV",
             &C::template int_N_scalar_NT_dV<xt::pytensor<double, 2>, xt::pytensor<double, 3>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::int_N_scalar_NT_dV`.",
             py::arg("qscalar"),
             py::arg("elemmat"));
 
     cls.def("Int_gradN_dot_tensor2_dV",
             &C::template Int_gradN_dot_tensor2_dV<xt::pytensor<double, 4>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::Int_gradN_dot_tensor2_dV`.",
             py::arg("qtensor"));
 
     cls.def("int_gradN_dot_tensor2_dV",
             &C::template int_gradN_dot_tensor2_dV<xt::pytensor<double, 4>, xt::pytensor<double, 3>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::int_gradN_dot_tensor2_dV`.",
             py::arg("qtensor"),
             py::arg("elemvec"));
 
     cls.def("Int_gradN_dot_tensor4_dot_gradNT_dV",
             &C::template Int_gradN_dot_tensor4_dot_gradNT_dV<xt::pytensor<double, 6>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::Int_gradN_dot_tensor4_dot_gradNT_dV`.",
             py::arg("qtensor"));
 
     cls.def("int_gradN_dot_tensor4_dot_gradNT_dV",
             &C::template int_gradN_dot_tensor4_dot_gradNT_dV<xt::pytensor<double, 6>, xt::pytensor<double, 3>>,
             "See :cpp:class:`GooseFEM::Quad4::Quadrature::int_gradN_dot_tensor4_dot_gradNT_dV`.",
             py::arg("qtensor"),
             py::arg("elemmat"));
 }
 
 void init_Element(py::module& m)
 {
 
     m.def("asElementVector",
         &GooseFEM::Element::asElementVector,
         "Convert nodal vector [nnode, ndim] to nodal vector stored per element [nelem, nne, ndim]",
         py::arg("conn"),
         py::arg("nodevec"));
 
     m.def("assembleElementVector",
         &GooseFEM::Element::assembleNodeVector,
         "Assemble nodal vector stored per element [nelem, nne, ndim] to nodal vector [nnode, ndim]",
         py::arg("conn"),
         py::arg("elemvec"));
 }
 
 #endif
diff --git a/python/Mesh.hpp b/python/Mesh.hpp
index 09e44af..58ecf38 100644
--- a/python/Mesh.hpp
+++ b/python/Mesh.hpp
@@ -1,436 +1,594 @@
-/* =================================================================================================
+/**
+\file
+\copyright Copyright 2017. Tom de Geus. All rights reserved.
+\license This project is released under the GNU Public License (GPLv3).
+*/
 
-(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
-
-================================================================================================= */
+#ifndef PYGOOSEFEM_MESH_H
+#define PYGOOSEFEM_MESH_H
 
 #include <GooseFEM/Mesh.h>
 #include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
 #include <xtensor-python/pyarray.hpp>
 #include <xtensor-python/pytensor.hpp>
 
 namespace py = pybind11;
 
-void init_Mesh(py::module& m)
+template <class C, class P>
+void register_Element_RegularBase(P& cls)
+{
+    cls.def("nelem", &C::nelem);
+    cls.def("nnode", &C::nnode);
+    cls.def("nne", &C::nne);
+    cls.def("ndim", &C::ndim);
+    cls.def("nelx", &C::nelx);
+    cls.def("nely", &C::nely);
+    cls.def("h", &C::h);
+    cls.def("getElementType", &C::getElementType);
+    cls.def("coor", &C::coor);
+    cls.def("conn", &C::conn);
+    cls.def("dofs", &C::dofs);
+    cls.def("dofsPeriodic", &C::dofsPeriodic);
+    cls.def("nodesPeriodic", &C::nodesPeriodic);
+    cls.def("nodesOrigin", &C::nodesOrigin);
+}
+
+template <class C, class P>
+void register_Element_RegularBase2d(P& cls)
+{
+     cls.def("nodesBottomEdge", &C::nodesBottomEdge);
+     cls.def("nodesTopEdge", &C::nodesTopEdge);
+     cls.def("nodesLeftEdge", &C::nodesLeftEdge);
+     cls.def("nodesRightEdge", &C::nodesRightEdge);
+     cls.def("nodesBottomOpenEdge", &C::nodesBottomOpenEdge);
+     cls.def("nodesTopOpenEdge", &C::nodesTopOpenEdge);
+     cls.def("nodesLeftOpenEdge", &C::nodesLeftOpenEdge);
+     cls.def("nodesRightOpenEdge", &C::nodesRightOpenEdge);
+     cls.def("nodesBottomLeftCorner", &C::nodesBottomLeftCorner);
+     cls.def("nodesBottomRightCorner", &C::nodesBottomRightCorner);
+     cls.def("nodesTopLeftCorner", &C::nodesTopLeftCorner);
+     cls.def("nodesTopRightCorner", &C::nodesTopRightCorner);
+     cls.def("nodesLeftBottomCorner", &C::nodesLeftBottomCorner);
+     cls.def("nodesLeftTopCorner", &C::nodesLeftTopCorner);
+     cls.def("nodesRightBottomCorner", &C::nodesRightBottomCorner);
+     cls.def("nodesRightTopCorner", &C::nodesRightTopCorner);
+}
+
+template <class C, class P>
+void register_Element_RegularBase3d(P& cls)
 {
+     cls.def("nodesFront", &C::nodesFront);
+     cls.def("nodesBack", &C::nodesBack);
+     cls.def("nodesLeft", &C::nodesLeft);
+     cls.def("nodesRight", &C::nodesRight);
+     cls.def("nodesBottom", &C::nodesBottom);
+     cls.def("nodesTop", &C::nodesTop);
+     cls.def("nodesFrontFace", &C::nodesFrontFace);
+     cls.def("nodesBackFace", &C::nodesBackFace);
+     cls.def("nodesLeftFace", &C::nodesLeftFace);
+     cls.def("nodesRightFace", &C::nodesRightFace);
+     cls.def("nodesBottomFace", &C::nodesBottomFace);
+     cls.def("nodesTopFace", &C::nodesTopFace);
+     cls.def("nodesFrontBottomEdge", &C::nodesFrontBottomEdge);
+     cls.def("nodesFrontTopEdge", &C::nodesFrontTopEdge);
+     cls.def("nodesFrontLeftEdge", &C::nodesFrontLeftEdge);
+     cls.def("nodesFrontRightEdge", &C::nodesFrontRightEdge);
+     cls.def("nodesBackBottomEdge", &C::nodesBackBottomEdge);
+     cls.def("nodesBackTopEdge", &C::nodesBackTopEdge);
+     cls.def("nodesBackLeftEdge", &C::nodesBackLeftEdge);
+     cls.def("nodesBackRightEdge", &C::nodesBackRightEdge);
+     cls.def("nodesBottomLeftEdge", &C::nodesBottomLeftEdge);
+     cls.def("nodesBottomRightEdge", &C::nodesBottomRightEdge);
+     cls.def("nodesTopLeftEdge", &C::nodesTopLeftEdge);
+     cls.def("nodesTopRightEdge", &C::nodesTopRightEdge);
+     cls.def("nodesBottomFrontEdge", &C::nodesBottomFrontEdge);
+     cls.def("nodesBottomBackEdge", &C::nodesBottomBackEdge);
+     cls.def("nodesTopFrontEdge", &C::nodesTopFrontEdge);
+     cls.def("nodesTopBackEdge", &C::nodesTopBackEdge);
+     cls.def("nodesLeftBottomEdge", &C::nodesLeftBottomEdge);
+     cls.def("nodesLeftFrontEdge", &C::nodesLeftFrontEdge);
+     cls.def("nodesLeftBackEdge", &C::nodesLeftBackEdge);
+     cls.def("nodesLeftTopEdge", &C::nodesLeftTopEdge);
+     cls.def("nodesRightBottomEdge", &C::nodesRightBottomEdge);
+     cls.def("nodesRightTopEdge", &C::nodesRightTopEdge);
+     cls.def("nodesRightFrontEdge", &C::nodesRightFrontEdge);
+     cls.def("nodesRightBackEdge", &C::nodesRightBackEdge);
+     cls.def("nodesFrontBottomOpenEdge", &C::nodesFrontBottomOpenEdge);
+     cls.def("nodesFrontTopOpenEdge", &C::nodesFrontTopOpenEdge);
+     cls.def("nodesFrontLeftOpenEdge", &C::nodesFrontLeftOpenEdge);
+     cls.def("nodesFrontRightOpenEdge", &C::nodesFrontRightOpenEdge);
+     cls.def("nodesBackBottomOpenEdge", &C::nodesBackBottomOpenEdge);
+     cls.def("nodesBackTopOpenEdge", &C::nodesBackTopOpenEdge);
+     cls.def("nodesBackLeftOpenEdge", &C::nodesBackLeftOpenEdge);
+     cls.def("nodesBackRightOpenEdge", &C::nodesBackRightOpenEdge);
+     cls.def("nodesBottomLeftOpenEdge", &C::nodesBottomLeftOpenEdge);
+     cls.def("nodesBottomRightOpenEdge", &C::nodesBottomRightOpenEdge);
+     cls.def("nodesTopLeftOpenEdge", &C::nodesTopLeftOpenEdge);
+     cls.def("nodesTopRightOpenEdge", &C::nodesTopRightOpenEdge);
+     cls.def("nodesBottomFrontOpenEdge", &C::nodesBottomFrontOpenEdge);
+     cls.def("nodesBottomBackOpenEdge", &C::nodesBottomBackOpenEdge);
+     cls.def("nodesTopFrontOpenEdge", &C::nodesTopFrontOpenEdge);
+     cls.def("nodesTopBackOpenEdge", &C::nodesTopBackOpenEdge);
+     cls.def("nodesLeftBottomOpenEdge", &C::nodesLeftBottomOpenEdge);
+     cls.def("nodesLeftFrontOpenEdge", &C::nodesLeftFrontOpenEdge);
+     cls.def("nodesLeftBackOpenEdge", &C::nodesLeftBackOpenEdge);
+     cls.def("nodesLeftTopOpenEdge", &C::nodesLeftTopOpenEdge);
+     cls.def("nodesRightBottomOpenEdge", &C::nodesRightBottomOpenEdge);
+     cls.def("nodesRightTopOpenEdge", &C::nodesRightTopOpenEdge);
+     cls.def("nodesRightFrontOpenEdge", &C::nodesRightFrontOpenEdge);
+     cls.def("nodesRightBackOpenEdge", &C::nodesRightBackOpenEdge);
+     cls.def("nodesFrontBottomLeftCorner", &C::nodesFrontBottomLeftCorner);
+     cls.def("nodesFrontBottomRightCorner", &C::nodesFrontBottomRightCorner);
+     cls.def("nodesFrontTopLeftCorner", &C::nodesFrontTopLeftCorner);
+     cls.def("nodesFrontTopRightCorner", &C::nodesFrontTopRightCorner);
+     cls.def("nodesBackBottomLeftCorner", &C::nodesBackBottomLeftCorner);
+     cls.def("nodesBackBottomRightCorner", &C::nodesBackBottomRightCorner);
+     cls.def("nodesBackTopLeftCorner", &C::nodesBackTopLeftCorner);
+     cls.def("nodesBackTopRightCorner", &C::nodesBackTopRightCorner);
+     cls.def("nodesFrontLeftBottomCorner", &C::nodesFrontLeftBottomCorner);
+     cls.def("nodesBottomFrontLeftCorner", &C::nodesBottomFrontLeftCorner);
+     cls.def("nodesBottomLeftFrontCorner", &C::nodesBottomLeftFrontCorner);
+     cls.def("nodesLeftFrontBottomCorner", &C::nodesLeftFrontBottomCorner);
+     cls.def("nodesLeftBottomFrontCorner", &C::nodesLeftBottomFrontCorner);
+     cls.def("nodesFrontRightBottomCorner", &C::nodesFrontRightBottomCorner);
+     cls.def("nodesBottomFrontRightCorner", &C::nodesBottomFrontRightCorner);
+     cls.def("nodesBottomRightFrontCorner", &C::nodesBottomRightFrontCorner);
+     cls.def("nodesRightFrontBottomCorner", &C::nodesRightFrontBottomCorner);
+     cls.def("nodesRightBottomFrontCorner", &C::nodesRightBottomFrontCorner);
+     cls.def("nodesFrontLeftTopCorner", &C::nodesFrontLeftTopCorner);
+     cls.def("nodesTopFrontLeftCorner", &C::nodesTopFrontLeftCorner);
+     cls.def("nodesTopLeftFrontCorner", &C::nodesTopLeftFrontCorner);
+     cls.def("nodesLeftFrontTopCorner", &C::nodesLeftFrontTopCorner);
+     cls.def("nodesLeftTopFrontCorner", &C::nodesLeftTopFrontCorner);
+     cls.def("nodesFrontRightTopCorner", &C::nodesFrontRightTopCorner);
+     cls.def("nodesTopFrontRightCorner", &C::nodesTopFrontRightCorner);
+     cls.def("nodesTopRightFrontCorner", &C::nodesTopRightFrontCorner);
+     cls.def("nodesRightFrontTopCorner", &C::nodesRightFrontTopCorner);
+     cls.def("nodesRightTopFrontCorner", &C::nodesRightTopFrontCorner);
+     cls.def("nodesBackLeftBottomCorner", &C::nodesBackLeftBottomCorner);
+     cls.def("nodesBottomBackLeftCorner", &C::nodesBottomBackLeftCorner);
+     cls.def("nodesBottomLeftBackCorner", &C::nodesBottomLeftBackCorner);
+     cls.def("nodesLeftBackBottomCorner", &C::nodesLeftBackBottomCorner);
+     cls.def("nodesLeftBottomBackCorner", &C::nodesLeftBottomBackCorner);
+     cls.def("nodesBackRightBottomCorner", &C::nodesBackRightBottomCorner);
+     cls.def("nodesBottomBackRightCorner", &C::nodesBottomBackRightCorner);
+     cls.def("nodesBottomRightBackCorner", &C::nodesBottomRightBackCorner);
+     cls.def("nodesRightBackBottomCorner", &C::nodesRightBackBottomCorner);
+     cls.def("nodesRightBottomBackCorner", &C::nodesRightBottomBackCorner);
+     cls.def("nodesBackLeftTopCorner", &C::nodesBackLeftTopCorner);
+     cls.def("nodesTopBackLeftCorner", &C::nodesTopBackLeftCorner);
+     cls.def("nodesTopLeftBackCorner", &C::nodesTopLeftBackCorner);
+     cls.def("nodesLeftBackTopCorner", &C::nodesLeftBackTopCorner);
+     cls.def("nodesLeftTopBackCorner", &C::nodesLeftTopBackCorner);
+     cls.def("nodesBackRightTopCorner", &C::nodesBackRightTopCorner);
+     cls.def("nodesTopBackRightCorner", &C::nodesTopBackRightCorner);
+     cls.def("nodesTopRightBackCorner", &C::nodesTopRightBackCorner);
+     cls.def("nodesRightBackTopCorner", &C::nodesRightBackTopCorner);
+     cls.def("nodesRightTopBackCorner", &C::nodesRightTopBackCorner);
+}
 
-    py::enum_<GooseFEM::Mesh::ElementType>(m,
+void init_Mesh(py::module& mod)
+{
+    py::enum_<GooseFEM::Mesh::ElementType>(mod,
             "ElementType",
             "ElementType."
             "See :cpp:enum:`GooseFEM::Mesh::ElementType`.")
         .value("Unknown", GooseFEM::Mesh::ElementType::Unknown)
         .value("Tri3", GooseFEM::Mesh::ElementType::Tri3)
         .value("Quad4", GooseFEM::Mesh::ElementType::Quad4)
         .value("Hex8", GooseFEM::Mesh::ElementType::Hex8)
         .export_values();
 
-    m.def("overlapping",
-          &GooseFEM::Mesh::overlapping<xt::pytensor<double, 2>, xt::pytensor<double, 2>>,
-          "Find overlapping nodes."
-          "See :cpp:func:`GooseFEM::Mesh::overlapping`.",
-          py::arg("coor_a"),
-          py::arg("coor_b"),
-          py::arg("rtol") = 1e-5,
-          py::arg("atol") = 1e-8);
+    mod.def("overlapping",
+            &GooseFEM::Mesh::overlapping<xt::pytensor<double, 2>, xt::pytensor<double, 2>>,
+            "Find overlapping nodes."
+            "See :cpp:func:`GooseFEM::Mesh::overlapping`.",
+            py::arg("coor_a"),
+            py::arg("coor_b"),
+            py::arg("rtol") = 1e-5,
+            py::arg("atol") = 1e-8);
 
-    py::class_<GooseFEM::Mesh::ManualStitch>(m, "ManualStitch")
+    py::class_<GooseFEM::Mesh::ManualStitch>(mod, "ManualStitch")
 
         .def(py::init<
                 const xt::pytensor<double, 2>&,
                 const xt::pytensor<size_t, 2>&,
                 const xt::pytensor<size_t, 1>&,
                 const xt::pytensor<double, 2>&,
                 const xt::pytensor<size_t, 2>&,
                 const xt::pytensor<size_t, 1>&,
                 bool,
                 double,
                 double>(),
              "Manually stitch meshes."
              "See :cpp:class:`GooseFEM::Mesh::ManualStitch`.",
              py::arg("coor_a"),
              py::arg("conn_a"),
              py::arg("overlapping_nodes_a"),
              py::arg("coor_b"),
              py::arg("conn_b"),
              py::arg("overlapping_nodes_b"),
              py::arg("check_position") = true,
              py::arg("rtol") = 1e-5,
              py::arg("atol") = 1e-8)
 
         .def("nmesh",
              &GooseFEM::Mesh::ManualStitch::nmesh,
              "Number of sub meshes."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nmesh`.")
 
         .def("nnode",
              &GooseFEM::Mesh::ManualStitch::nnode,
              "Number of nodes of stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nnode`.")
 
         .def("nelem",
              &GooseFEM::Mesh::ManualStitch::nelem,
              "Number of elements of stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nelem`.")
 
         .def("ndim",
              &GooseFEM::Mesh::ManualStitch::ndim,
              "Number of dimensions (of stitched mesh)."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::ndim`.")
 
         .def("nne",
              &GooseFEM::Mesh::ManualStitch::nne,
              "Number of nodes per element (of stitched mesh)."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nne`.")
 
         .def("coor",
              &GooseFEM::Mesh::ManualStitch::coor,
              "Coordinates of stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::coor`.")
 
         .def("conn",
              &GooseFEM::Mesh::ManualStitch::conn,
              "Connectivity of stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::conn`.")
 
         .def("dofs",
              &GooseFEM::Mesh::ManualStitch::dofs,
              "DOF numbers per node."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::dofs`.")
 
         .def("nodemap",
              py::overload_cast<>(
                 &GooseFEM::Mesh::ManualStitch::nodemap, py::const_),
              "Node-map for givall sub-meshes."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nodemap`.")
 
         .def("elemmap",
              py::overload_cast<>(
                 &GooseFEM::Mesh::ManualStitch::elemmap, py::const_),
              "Element-map for all sub-meshes."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::elemmap`.")
 
         .def("nodemap",
              py::overload_cast<size_t>(
                 &GooseFEM::Mesh::ManualStitch::nodemap, py::const_),
              "Node-map for given sub-mesh."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nodemap`.",
              py::arg("mesh_index"))
 
         .def("elemmap",
              py::overload_cast<size_t>(
                 &GooseFEM::Mesh::ManualStitch::elemmap, py::const_),
              "Element-map for given sub-mesh."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::elemmap`.",
              py::arg("mesh_index"))
 
         .def("nodeset",
              &GooseFEM::Mesh::ManualStitch::nodeset<xt::pytensor<size_t, 1>>,
              "Convert node-set to the stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nodeset`.",
              py::arg("set"),
              py::arg("mesh_index"))
 
         .def("elemset",
              &GooseFEM::Mesh::ManualStitch::elemset<xt::pytensor<size_t, 1>>,
              "Convert element-set to the stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::elemset`.",
              py::arg("set"),
              py::arg("mesh_index"))
 
         .def("__repr__",
             [](const GooseFEM::Mesh::ManualStitch&) { return "<GooseFEM.Mesh.ManualStitch>"; });
 
-    py::class_<GooseFEM::Mesh::Stitch>(m, "Stitch")
+    py::class_<GooseFEM::Mesh::Stitch>(mod, "Stitch")
 
         .def(py::init<double, double>(),
              "Stitch meshes."
              "See :cpp:class:`GooseFEM::Mesh::Stitch`.",
              py::arg("rtol") = 1e-5,
              py::arg("atol") = 1e-8)
 
         .def("push_back",
              &GooseFEM::Mesh::Stitch::push_back<xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>,
              "Add mesh to be stitched."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::push_back`.",
              py::arg("coor"),
              py::arg("conn"))
 
         .def("nmesh",
              &GooseFEM::Mesh::Stitch::nmesh,
              "Number of sub meshes."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::nmesh`.")
 
         .def("nnode",
              &GooseFEM::Mesh::Stitch::nnode,
              "Number of nodes of stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::nnode`.")
 
         .def("nelem",
              &GooseFEM::Mesh::Stitch::nelem,
              "Number of elements of stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::nelem`.")
 
         .def("ndim",
              &GooseFEM::Mesh::Stitch::ndim,
              "Number of dimensions (of stitched mesh)."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::ndim`.")
 
         .def("nne",
              &GooseFEM::Mesh::Stitch::nne,
              "Number of nodes per element (of stitched mesh)."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::nne`.")
 
         .def("coor",
              &GooseFEM::Mesh::Stitch::coor,
              "Coordinates of stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::coor`.")
 
         .def("conn",
              &GooseFEM::Mesh::Stitch::conn,
              "Connectivity of stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::conn`.")
 
         .def("dofs",
              &GooseFEM::Mesh::Stitch::dofs,
              "DOF numbers per node."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::dofs`.")
 
         .def("nodemap",
              py::overload_cast<>(
                 &GooseFEM::Mesh::Stitch::nodemap, py::const_),
              "Node-map for givall sub-meshes."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::nodemap`.")
 
         .def("elemmap",
              py::overload_cast<>(
                 &GooseFEM::Mesh::Stitch::elemmap, py::const_),
              "Element-map for all sub-meshes."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::elemmap`.")
 
         .def("nodemap",
              py::overload_cast<size_t>(
                 &GooseFEM::Mesh::Stitch::nodemap, py::const_),
              "Node-map for given sub-mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::nodemap`.",
              py::arg("mesh_index"))
 
         .def("elemmap",
              py::overload_cast<size_t>(
                 &GooseFEM::Mesh::Stitch::elemmap, py::const_),
              "Element-map for given sub-mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::elemmap`.",
              py::arg("mesh_index"))
 
         .def("nodeset",
              static_cast<xt::pytensor<size_t, 1> (GooseFEM::Mesh::Stitch::*)(
                 const xt::pytensor<size_t, 1>&, size_t) const>(&GooseFEM::Mesh::Stitch::nodeset),
              "Convert node-set to the stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::nodeset`.",
              py::arg("set"),
              py::arg("mesh_index"))
 
         .def("elemset",
              static_cast<xt::pytensor<size_t, 1> (GooseFEM::Mesh::Stitch::*)(
                 const xt::pytensor<size_t, 1>&, size_t) const>(&GooseFEM::Mesh::Stitch::elemset),
              "Convert element-set to the stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::elemset`.",
              py::arg("set"),
              py::arg("mesh_index"))
 
         .def("nodeset",
              static_cast<xt::pytensor<size_t, 1> (GooseFEM::Mesh::Stitch::*)(
                 const std::vector<xt::pytensor<size_t, 1>>&) const>(&GooseFEM::Mesh::Stitch::nodeset),
              "Convert node-set to the stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::nodeset`.",
              py::arg("set"))
 
         .def("elemset",
              static_cast<xt::pytensor<size_t, 1> (GooseFEM::Mesh::Stitch::*)(
                 const std::vector<xt::pytensor<size_t, 1>>&) const>(&GooseFEM::Mesh::Stitch::elemset),
              "Convert element-set to the stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::elemset`.",
              py::arg("set"))
 
         .def("__repr__", [](const GooseFEM::Mesh::Stitch&) { return "<GooseFEM.Mesh.Stitch>"; });
 
-    py::class_<GooseFEM::Mesh::Vstack, GooseFEM::Mesh::Stitch>(m, "Vstack")
+    py::class_<GooseFEM::Mesh::Vstack, GooseFEM::Mesh::Stitch>(mod, "Vstack")
 
         .def(py::init<bool, double, double>(),
              "Vstack meshes."
              "See :cpp:class:`GooseFEM::Mesh::Vstack`.",
              py::arg("check_overlap") = true,
              py::arg("rtol") = 1e-5,
              py::arg("atol") = 1e-8)
 
         .def("push_back",
              &GooseFEM::Mesh::Vstack::push_back<xt::pytensor<double, 2>, xt::pytensor<size_t, 2>, xt::pytensor<size_t, 1>>,
              "Add mesh to be stitched."
              "See :cpp:func:`GooseFEM::Mesh::Vstack::push_back`.",
              py::arg("coor"),
              py::arg("conn"),
              py::arg("nodes_bottom"),
              py::arg("nodes_top"))
 
         .def("__repr__", [](const GooseFEM::Mesh::Vstack&) { return "<GooseFEM.Mesh.Vstack>"; });
 
-    py::class_<GooseFEM::Mesh::Renumber>(m, "Renumber")
+    py::class_<GooseFEM::Mesh::Renumber>(mod, "Renumber")
 
         .def(py::init<const xt::pyarray<size_t>&>(),
              "Renumber to lowest possible index."
              "See :cpp:class:`GooseFEM::Mesh::Renumber`.",
              py::arg("dofs"))
 
         .def("get",
              &GooseFEM::Mesh::Renumber::apply<xt::pyarray<size_t>>,
              "Get renumbered DOFs."
              "See :cpp:func:`GooseFEM::Mesh::Renumber::get`.")
 
         .def("apply",
              &GooseFEM::Mesh::Renumber::apply<xt::pyarray<size_t>>,
              "Get renumbered list."
              "See :cpp:func:`GooseFEM::Mesh::Renumber::apply`.")
 
         .def("index",
              &GooseFEM::Mesh::Renumber::index,
              "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``."
              "See :cpp:func:`GooseFEM::Mesh::Renumber::index`.")
 
         .def("__repr__",
              [](const GooseFEM::Mesh::Renumber&) { return "<GooseFEM.Mesh.Renumber>"; });
 
-    py::class_<GooseFEM::Mesh::Reorder>(m, "Reorder")
+    py::class_<GooseFEM::Mesh::Reorder>(mod, "Reorder")
 
         .def(py::init([](xt::pytensor<size_t, 1>& a) { return new GooseFEM::Mesh::Reorder({a}); }),
              "Reorder to lowest possible index."
              "See :cpp:class:`GooseFEM::Mesh::Reorder`.")
 
         .def(py::init([](xt::pytensor<size_t, 1>& a, xt::pytensor<size_t, 1>& b) {
                 return new GooseFEM::Mesh::Reorder({a, b});
              }),
              "Reorder to lowest possible index."
              "See :cpp:class:`GooseFEM::Mesh::Reorder`.")
 
         .def(py::init(
                 [](xt::pytensor<size_t, 1>& a, xt::pytensor<size_t, 1>& b, xt::pytensor<size_t, 1>& c) {
                     return new GooseFEM::Mesh::Reorder({a, b, c});
                 }),
              "Reorder to lowest possible index."
              "See :cpp:class:`GooseFEM::Mesh::Reorder`.")
 
         .def(py::init([](xt::pytensor<size_t, 1>& a,
                          xt::pytensor<size_t, 1>& b,
                          xt::pytensor<size_t, 1>& c,
                          xt::pytensor<size_t, 1>& d) {
                 return new GooseFEM::Mesh::Reorder({a, b, c, d});
                 }),
              "Reorder to lowest possible index."
              "See :cpp:class:`GooseFEM::Mesh::Reorder`.")
 
         .def("get",
              &GooseFEM::Mesh::Reorder::apply<xt::pyarray<size_t>>,
              "Reorder matrix (e.g. ``dofs``)."
              "See :cpp:func:`GooseFEM::Mesh::Reorder::get`.",
              py::arg("dofs"))
 
         .def("apply",
              &GooseFEM::Mesh::Reorder::apply<xt::pyarray<size_t>>,
              "Get reordered list."
              "See :cpp:func:`GooseFEM::Mesh::Reorder::apply`.")
 
         .def("index",
              &GooseFEM::Mesh::Reorder::index,
              "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``."
              "See :cpp:func:`GooseFEM::Mesh::Reorder::index`.")
 
         .def("__repr__", [](const GooseFEM::Mesh::Reorder&) { return "<GooseFEM.Mesh.Reorder>"; });
 
-    m.def("dofs",
-          &GooseFEM::Mesh::dofs,
-          "List with DOF-numbers in sequential order."
-          "See :cpp:func:`GooseFEM::Mesh::dofs`.",
-          py::arg("nnode"),
-          py::arg("ndim"));
-
-    m.def("renumber",
-          &GooseFEM::Mesh::renumber<xt::pytensor<size_t, 2>>,
-          "Renumber to lowest possible indices."
-          "See :cpp:func:`GooseFEM::Mesh::renumber`.",
-          py::arg("dofs"));
-
-    m.def("coordination",
-          &GooseFEM::Mesh::coordination<xt::pytensor<size_t, 2>>,
-          "Coordination number of each node."
-          "See :cpp:func:`GooseFEM::Mesh::coordination`.",
-          py::arg("conn"));
-
-    m.def("elem2node",
-          &GooseFEM::Mesh::elem2node<xt::pytensor<size_t, 2>>,
-          "Element-numbers connected to each node."
-          "See :cpp:func:`GooseFEM::Mesh::elem2node`.",
-          py::arg("conn"),
-          py::arg("sorted") = true);
-
-    m.def("edgesize",
-          py::overload_cast<const xt::pytensor<double, 2>&, const xt::pytensor<size_t, 2>&>(
-            &GooseFEM::Mesh::edgesize<xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
-          "Get the edge size of all elements."
-          "See :cpp:func:`GooseFEM::Mesh::edgesize`.",
-          py::arg("coor"),
-          py::arg("conn"));
-
-    m.def("edgesize",
-          py::overload_cast<
-            const xt::pytensor<double, 2>&,
-            const xt::pytensor<size_t, 2>&,
-            GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::edgesize<xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
-          "Get the edge size of all elements."
-          "See :cpp:func:`GooseFEM::Mesh::edgesize`.",
-          py::arg("coor"),
-          py::arg("conn"),
-          py::arg("type"));
-
-    m.def("centers",
-          py::overload_cast<const xt::pytensor<double, 2>&, const xt::pytensor<size_t, 2>&>(
-            &GooseFEM::Mesh::centers<xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
-          "Coordinates of the center of each element."
-          "See :cpp:func:`GooseFEM::Mesh::centers`.",
-          py::arg("coor"),
-          py::arg("conn"));
-
-    m.def("centers",
-          py::overload_cast<
-            const xt::pytensor<double, 2>&,
-            const xt::pytensor<size_t, 2>&,
-            GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::centers<xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
-          "Coordinates of the center of each element."
-          "See :cpp:func:`GooseFEM::Mesh::centers`.",
-          py::arg("coor"),
-          py::arg("conn"),
-          py::arg("type"));
-
-    m.def("elemmap2nodemap",
-          py::overload_cast<
-            const xt::pytensor<size_t, 1>&,
-            const xt::pytensor<double, 2>&,
-            const xt::pytensor<size_t, 2>&>(&GooseFEM::Mesh::elemmap2nodemap<
-                xt::pytensor<size_t, 1>, xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
-          "Convert an element-map to a node-map."
-          "See :cpp:func:`GooseFEM::Mesh::elemmap2nodemap`.",
-          py::arg("elem_map"),
-          py::arg("coor"),
-          py::arg("conn"));
-
-    m.def("elemmap2nodemap",
-          py::overload_cast<
-            const xt::pytensor<size_t, 1>&,
-            const xt::pytensor<double, 2>&,
-            const xt::pytensor<size_t, 2>&,
-            GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::elemmap2nodemap<
-                xt::pytensor<size_t, 1>, xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
-          "Convert an element-map to a node-map."
-          "See :cpp:func:`GooseFEM::Mesh::elemmap2nodemap`.",
-          py::arg("elem_map"),
-          py::arg("coor"),
-          py::arg("conn"),
-          py::arg("type"));
+    mod.def("dofs",
+            &GooseFEM::Mesh::dofs,
+            "List with DOF-numbers in sequential order."
+            "See :cpp:func:`GooseFEM::Mesh::dofs`.",
+            py::arg("nnode"),
+            py::arg("ndim"));
+
+    mod.def("renumber",
+            &GooseFEM::Mesh::renumber<xt::pytensor<size_t, 2>>,
+            "Renumber to lowest possible indices."
+            "See :cpp:func:`GooseFEM::Mesh::renumber`.",
+            py::arg("dofs"));
+
+    mod.def("coordination",
+            &GooseFEM::Mesh::coordination<xt::pytensor<size_t, 2>>,
+            "Coordination number of each node."
+            "See :cpp:func:`GooseFEM::Mesh::coordination`.",
+            py::arg("conn"));
+
+    mod.def("elem2node",
+            &GooseFEM::Mesh::elem2node<xt::pytensor<size_t, 2>>,
+            "Element-numbers connected to each node."
+            "See :cpp:func:`GooseFEM::Mesh::elem2node`.",
+            py::arg("conn"),
+            py::arg("sorted") = true);
+
+    mod.def("edgesize",
+            py::overload_cast<const xt::pytensor<double, 2>&, const xt::pytensor<size_t, 2>&>(
+                &GooseFEM::Mesh::edgesize<xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
+            "Get the edge size of all elements."
+            "See :cpp:func:`GooseFEM::Mesh::edgesize`.",
+            py::arg("coor"),
+            py::arg("conn"));
+
+    mod.def("edgesize",
+            py::overload_cast<
+                const xt::pytensor<double, 2>&,
+                const xt::pytensor<size_t, 2>&,
+                GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::edgesize<xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
+            "Get the edge size of all elements."
+            "See :cpp:func:`GooseFEM::Mesh::edgesize`.",
+            py::arg("coor"),
+            py::arg("conn"),
+            py::arg("type"));
+
+    mod.def("centers",
+            py::overload_cast<const xt::pytensor<double, 2>&, const xt::pytensor<size_t, 2>&>(
+                &GooseFEM::Mesh::centers<xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
+            "Coordinates of the center of each element."
+            "See :cpp:func:`GooseFEM::Mesh::centers`.",
+            py::arg("coor"),
+            py::arg("conn"));
+
+    mod.def("centers",
+            py::overload_cast<
+                const xt::pytensor<double, 2>&,
+                const xt::pytensor<size_t, 2>&,
+                GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::centers<xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
+            "Coordinates of the center of each element."
+            "See :cpp:func:`GooseFEM::Mesh::centers`.",
+            py::arg("coor"),
+            py::arg("conn"),
+            py::arg("type"));
+
+    mod.def("elemmap2nodemap",
+            py::overload_cast<
+                const xt::pytensor<size_t, 1>&,
+                const xt::pytensor<double, 2>&,
+                const xt::pytensor<size_t, 2>&>(&GooseFEM::Mesh::elemmap2nodemap<
+                    xt::pytensor<size_t, 1>, xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
+            "Convert an element-map to a node-map."
+            "See :cpp:func:`GooseFEM::Mesh::elemmap2nodemap`.",
+            py::arg("elem_map"),
+            py::arg("coor"),
+            py::arg("conn"));
+
+    mod.def("elemmap2nodemap",
+            py::overload_cast<
+                const xt::pytensor<size_t, 1>&,
+                const xt::pytensor<double, 2>&,
+                const xt::pytensor<size_t, 2>&,
+                GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::elemmap2nodemap<
+                    xt::pytensor<size_t, 1>, xt::pytensor<double, 2>, xt::pytensor<size_t, 2>>),
+            "Convert an element-map to a node-map."
+            "See :cpp:func:`GooseFEM::Mesh::elemmap2nodemap`.",
+            py::arg("elem_map"),
+            py::arg("coor"),
+            py::arg("conn"),
+            py::arg("type"));
 }
+
+#endif
diff --git a/python/MeshHex8.hpp b/python/MeshHex8.hpp
index 7337bf7..a5eed0c 100644
--- a/python/MeshHex8.hpp
+++ b/python/MeshHex8.hpp
@@ -1,614 +1,60 @@
-/* =================================================================================================
+/**
+\file
+\copyright Copyright 2017. Tom de Geus. All rights reserved.
+\license This project is released under the GNU Public License (GPLv3).
+*/
 
-(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
+#ifndef PYGOOSEFEM_MESHHEX8_H
+#define PYGOOSEFEM_MESHHEX8_H
 
-================================================================================================= */
-
-#include <GooseFEM/GooseFEM.h>
+#include <GooseFEM/MeshHex8.h>
 #include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
 #include <xtensor-python/pyarray.hpp>
 #include <xtensor-python/pytensor.hpp>
 #include <pyxtensor/pyxtensor.hpp>
 
+#include "Mesh.hpp"
+
 namespace py = pybind11;
 
-void init_MeshHex8(py::module& m)
+void init_MeshHex8(py::module& mod)
 {
+    {
+        py::class_<GooseFEM::Mesh::Hex8::Regular> cls(mod, "Regular");
 
-    py::class_<GooseFEM::Mesh::Hex8::Regular>(m, "Regular")
-
-        .def(
-            py::init<size_t, size_t, size_t, double>(),
-            "Mesh with nx*ny*nz 'pixels' and edge size h",
-            py::arg("nx"),
-            py::arg("ny"),
-            py::arg("nz"),
-            py::arg("h") = 1.)
-
-        .def("nelem", &GooseFEM::Mesh::Hex8::Regular::nelem)
-
-        .def("nnode", &GooseFEM::Mesh::Hex8::Regular::nnode)
-
-        .def("nne", &GooseFEM::Mesh::Hex8::Regular::nne)
-
-        .def("ndim", &GooseFEM::Mesh::Hex8::Regular::ndim)
-
-        .def("coor", &GooseFEM::Mesh::Hex8::Regular::coor)
-
-        .def("conn", &GooseFEM::Mesh::Hex8::Regular::conn)
-
-        .def("getElementType", &GooseFEM::Mesh::Hex8::Regular::getElementType)
-
-        .def("nodesFront", &GooseFEM::Mesh::Hex8::Regular::nodesFront)
-
-        .def("nodesBack", &GooseFEM::Mesh::Hex8::Regular::nodesBack)
-
-        .def("nodesLeft", &GooseFEM::Mesh::Hex8::Regular::nodesLeft)
-
-        .def("nodesRight", &GooseFEM::Mesh::Hex8::Regular::nodesRight)
-
-        .def("nodesBottom", &GooseFEM::Mesh::Hex8::Regular::nodesBottom)
-
-        .def("nodesTop", &GooseFEM::Mesh::Hex8::Regular::nodesTop)
-
-        .def("nodesFrontFace", &GooseFEM::Mesh::Hex8::Regular::nodesFrontFace)
-
-        .def("nodesBackFace", &GooseFEM::Mesh::Hex8::Regular::nodesBackFace)
-
-        .def("nodesLeftFace", &GooseFEM::Mesh::Hex8::Regular::nodesLeftFace)
-
-        .def("nodesRightFace", &GooseFEM::Mesh::Hex8::Regular::nodesRightFace)
-
-        .def("nodesBottomFace", &GooseFEM::Mesh::Hex8::Regular::nodesBottomFace)
-
-        .def("nodesTopFace", &GooseFEM::Mesh::Hex8::Regular::nodesTopFace)
-
-        .def("nodesFrontBottomEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomEdge)
-
-        .def("nodesFrontTopEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopEdge)
-
-        .def("nodesFrontLeftEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftEdge)
-
-        .def("nodesFrontRightEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightEdge)
-
-        .def("nodesBackBottomEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomEdge)
-
-        .def("nodesBackTopEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackTopEdge)
-
-        .def("nodesBackLeftEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftEdge)
-
-        .def("nodesBackRightEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackRightEdge)
-
-        .def("nodesBottomLeftEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftEdge)
-
-        .def("nodesBottomRightEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightEdge)
-
-        .def("nodesTopLeftEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftEdge)
-
-        .def("nodesTopRightEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopRightEdge)
-
-        .def("nodesBottomFrontEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontEdge)
-
-        .def("nodesBottomBackEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackEdge)
-
-        .def("nodesTopFrontEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontEdge)
-
-        .def("nodesTopBackEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopBackEdge)
-
-        .def("nodesLeftBottomEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomEdge)
-
-        .def("nodesLeftFrontEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontEdge)
-
-        .def("nodesLeftBackEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackEdge)
-
-        .def("nodesLeftTopEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopEdge)
-
-        .def("nodesRightBottomEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomEdge)
-
-        .def("nodesRightTopEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightTopEdge)
-
-        .def("nodesRightFrontEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontEdge)
-
-        .def("nodesRightBackEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightBackEdge)
-
-        .def("nodesFrontBottomOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomOpenEdge)
-
-        .def("nodesFrontTopOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopOpenEdge)
-
-        .def("nodesFrontLeftOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftOpenEdge)
-
-        .def("nodesFrontRightOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightOpenEdge)
-
-        .def("nodesBackBottomOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomOpenEdge)
-
-        .def("nodesBackTopOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackTopOpenEdge)
-
-        .def("nodesBackLeftOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftOpenEdge)
-
-        .def("nodesBackRightOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBackRightOpenEdge)
-
-        .def("nodesBottomLeftOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftOpenEdge)
-
-        .def("nodesBottomRightOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightOpenEdge)
-
-        .def("nodesTopLeftOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftOpenEdge)
-
-        .def("nodesTopRightOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopRightOpenEdge)
-
-        .def("nodesBottomFrontOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontOpenEdge)
-
-        .def("nodesBottomBackOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackOpenEdge)
-
-        .def("nodesTopFrontOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontOpenEdge)
+        cls.def(py::init<size_t, size_t, size_t, double>(),
+                "Mesh with nx*ny*nz 'pixels' and edge size h",
+                py::arg("nx"),
+                py::arg("ny"),
+                py::arg("nz"),
+                py::arg("h") = 1.);
 
-        .def("nodesTopBackOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesTopBackOpenEdge)
+        register_Element_RegularBase<GooseFEM::Mesh::Hex8::Regular>(cls);
+        register_Element_RegularBase3d<GooseFEM::Mesh::Hex8::Regular>(cls);
 
-        .def("nodesLeftBottomOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomOpenEdge)
-
-        .def("nodesLeftFrontOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontOpenEdge)
-
-        .def("nodesLeftBackOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackOpenEdge)
-
-        .def("nodesLeftTopOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopOpenEdge)
-
-        .def("nodesRightBottomOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomOpenEdge)
-
-        .def("nodesRightTopOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightTopOpenEdge)
-
-        .def("nodesRightFrontOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontOpenEdge)
-
-        .def("nodesRightBackOpenEdge", &GooseFEM::Mesh::Hex8::Regular::nodesRightBackOpenEdge)
-
-        .def(
-            "nodesFrontBottomLeftCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomLeftCorner)
-
-        .def(
-            "nodesFrontBottomRightCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesFrontBottomRightCorner)
-
-        .def("nodesFrontTopLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopLeftCorner)
-
-        .def("nodesFrontTopRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontTopRightCorner)
-
-        .def("nodesBackBottomLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomLeftCorner)
-
-        .def(
-            "nodesBackBottomRightCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesBackBottomRightCorner)
-
-        .def("nodesBackTopLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackTopLeftCorner)
-
-        .def("nodesBackTopRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackTopRightCorner)
-
-        .def(
-            "nodesFrontLeftBottomCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftBottomCorner)
-
-        .def(
-            "nodesBottomFrontLeftCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontLeftCorner)
-
-        .def(
-            "nodesBottomLeftFrontCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftFrontCorner)
-
-        .def(
-            "nodesLeftFrontBottomCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontBottomCorner)
-
-        .def(
-            "nodesLeftBottomFrontCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomFrontCorner)
-
-        .def(
-            "nodesFrontRightBottomCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightBottomCorner)
-
-        .def(
-            "nodesBottomFrontRightCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesBottomFrontRightCorner)
-
-        .def(
-            "nodesBottomRightFrontCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightFrontCorner)
-
-        .def(
-            "nodesRightFrontBottomCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontBottomCorner)
-
-        .def(
-            "nodesRightBottomFrontCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomFrontCorner)
-
-        .def("nodesFrontLeftTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontLeftTopCorner)
-
-        .def("nodesTopFrontLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontLeftCorner)
-
-        .def("nodesTopLeftFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftFrontCorner)
-
-        .def("nodesLeftFrontTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftFrontTopCorner)
-
-        .def("nodesLeftTopFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopFrontCorner)
-
-        .def("nodesFrontRightTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesFrontRightTopCorner)
-
-        .def("nodesTopFrontRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopFrontRightCorner)
-
-        .def("nodesTopRightFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopRightFrontCorner)
-
-        .def("nodesRightFrontTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightFrontTopCorner)
-
-        .def("nodesRightTopFrontCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightTopFrontCorner)
-
-        .def("nodesBackLeftBottomCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftBottomCorner)
-
-        .def("nodesBottomBackLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackLeftCorner)
-
-        .def("nodesBottomLeftBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBottomLeftBackCorner)
-
-        .def("nodesLeftBackBottomCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackBottomCorner)
-
-        .def("nodesLeftBottomBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBottomBackCorner)
-
-        .def(
-            "nodesBackRightBottomCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesBackRightBottomCorner)
-
-        .def(
-            "nodesBottomBackRightCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesBottomBackRightCorner)
-
-        .def(
-            "nodesBottomRightBackCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesBottomRightBackCorner)
-
-        .def(
-            "nodesRightBackBottomCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesRightBackBottomCorner)
-
-        .def(
-            "nodesRightBottomBackCorner",
-            &GooseFEM::Mesh::Hex8::Regular::nodesRightBottomBackCorner)
-
-        .def("nodesBackLeftTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackLeftTopCorner)
-
-        .def("nodesTopBackLeftCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopBackLeftCorner)
-
-        .def("nodesTopLeftBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopLeftBackCorner)
-
-        .def("nodesLeftBackTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftBackTopCorner)
-
-        .def("nodesLeftTopBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesLeftTopBackCorner)
-
-        .def("nodesBackRightTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesBackRightTopCorner)
-
-        .def("nodesTopBackRightCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopBackRightCorner)
-
-        .def("nodesTopRightBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesTopRightBackCorner)
-
-        .def("nodesRightBackTopCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightBackTopCorner)
-
-        .def("nodesRightTopBackCorner", &GooseFEM::Mesh::Hex8::Regular::nodesRightTopBackCorner)
-
-        .def("nodesPeriodic", &GooseFEM::Mesh::Hex8::Regular::nodesPeriodic)
-
-        .def("nodesOrigin", &GooseFEM::Mesh::Hex8::Regular::nodesOrigin)
-
-        .def("dofs", &GooseFEM::Mesh::Hex8::Regular::dofs)
-
-        .def("dofsPeriodic", &GooseFEM::Mesh::Hex8::Regular::dofsPeriodic)
-
-        .def("__repr__", [](const GooseFEM::Mesh::Hex8::Regular&) {
+        cls.def("__repr__", [](const GooseFEM::Mesh::Hex8::Regular&) {
             return "<GooseFEM.Mesh.Hex8.Regular>";
         });
+    }
 
-    py::class_<GooseFEM::Mesh::Hex8::FineLayer>(m, "FineLayer")
-
-        .def(
-            py::init<size_t, size_t, size_t, double, size_t>(),
-            "Mesh with nx*ny*nz 'pixels' and edge size h",
-            py::arg("nx"),
-            py::arg("ny"),
-            py::arg("nz"),
-            py::arg("h") = 1.0,
-            py::arg("nfine") = 1)
+    {
+        py::class_<GooseFEM::Mesh::Hex8::FineLayer> cls(mod, "FineLayer");
 
-        .def("nelem", &GooseFEM::Mesh::Hex8::FineLayer::nelem)
+        cls.def(py::init<size_t, size_t, size_t, double, size_t>(),
+                "Mesh with nx*ny*nz 'pixels' and edge size h",
+                py::arg("nx"),
+                py::arg("ny"),
+                py::arg("nz"),
+                py::arg("h") = 1.0,
+                py::arg("nfine") = 1);
 
-        .def("nnode", &GooseFEM::Mesh::Hex8::FineLayer::nnode)
+        cls.def("elementsMiddleLayer", &GooseFEM::Mesh::Hex8::FineLayer::elementsMiddleLayer);
 
-        .def("nne", &GooseFEM::Mesh::Hex8::FineLayer::nne)
-
-        .def("ndim", &GooseFEM::Mesh::Hex8::FineLayer::ndim)
-
-        .def("nelx", &GooseFEM::Mesh::Hex8::FineLayer::nelx)
-
-        .def("nely", &GooseFEM::Mesh::Hex8::FineLayer::nely)
-
-        .def("nelz", &GooseFEM::Mesh::Hex8::FineLayer::nelz)
-
-        .def("coor", &GooseFEM::Mesh::Hex8::FineLayer::coor)
-
-        .def("conn", &GooseFEM::Mesh::Hex8::FineLayer::conn)
-
-        .def("getElementType", &GooseFEM::Mesh::Hex8::FineLayer::getElementType)
-
-        .def("elementsMiddleLayer", &GooseFEM::Mesh::Hex8::FineLayer::elementsMiddleLayer)
-
-        .def("nodesFront", &GooseFEM::Mesh::Hex8::FineLayer::nodesFront)
-
-        .def("nodesBack", &GooseFEM::Mesh::Hex8::FineLayer::nodesBack)
-
-        .def("nodesLeft", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeft)
-
-        .def("nodesRight", &GooseFEM::Mesh::Hex8::FineLayer::nodesRight)
-
-        .def("nodesBottom", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottom)
-
-        .def("nodesTop", &GooseFEM::Mesh::Hex8::FineLayer::nodesTop)
-
-        .def("nodesFrontFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontFace)
-
-        .def("nodesBackFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackFace)
-
-        .def("nodesLeftFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFace)
-
-        .def("nodesRightFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFace)
-
-        .def("nodesBottomFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFace)
-
-        .def("nodesTopFace", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFace)
-
-        .def("nodesFrontBottomEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomEdge)
-
-        .def("nodesFrontTopEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopEdge)
-
-        .def("nodesFrontLeftEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftEdge)
-
-        .def("nodesFrontRightEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightEdge)
-
-        .def("nodesBackBottomEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomEdge)
-
-        .def("nodesBackTopEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopEdge)
-
-        .def("nodesBackLeftEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftEdge)
-
-        .def("nodesBackRightEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightEdge)
-
-        .def("nodesBottomLeftEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftEdge)
-
-        .def("nodesBottomRightEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightEdge)
-
-        .def("nodesTopLeftEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftEdge)
-
-        .def("nodesTopRightEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightEdge)
-
-        .def("nodesBottomFrontEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontEdge)
-
-        .def("nodesBottomBackEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackEdge)
-
-        .def("nodesTopFrontEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontEdge)
-
-        .def("nodesTopBackEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackEdge)
-
-        .def("nodesLeftBottomEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomEdge)
-
-        .def("nodesLeftFrontEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontEdge)
-
-        .def("nodesLeftBackEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackEdge)
-
-        .def("nodesLeftTopEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopEdge)
-
-        .def("nodesRightBottomEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomEdge)
-
-        .def("nodesRightTopEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopEdge)
-
-        .def("nodesRightFrontEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontEdge)
-
-        .def("nodesRightBackEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackEdge)
-
-        .def("nodesFrontBottomOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomOpenEdge)
-
-        .def("nodesFrontTopOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopOpenEdge)
-
-        .def("nodesFrontLeftOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftOpenEdge)
-
-        .def("nodesFrontRightOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightOpenEdge)
-
-        .def("nodesBackBottomOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomOpenEdge)
-
-        .def("nodesBackTopOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopOpenEdge)
-
-        .def("nodesBackLeftOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftOpenEdge)
-
-        .def("nodesBackRightOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightOpenEdge)
-
-        .def("nodesBottomLeftOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftOpenEdge)
-
-        .def("nodesBottomRightOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightOpenEdge)
-
-        .def("nodesTopLeftOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftOpenEdge)
-
-        .def("nodesTopRightOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightOpenEdge)
-
-        .def("nodesBottomFrontOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontOpenEdge)
-
-        .def("nodesBottomBackOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackOpenEdge)
-
-        .def("nodesTopFrontOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontOpenEdge)
-
-        .def("nodesTopBackOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackOpenEdge)
-
-        .def("nodesLeftBottomOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomOpenEdge)
-
-        .def("nodesLeftFrontOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontOpenEdge)
-
-        .def("nodesLeftBackOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackOpenEdge)
-
-        .def("nodesLeftTopOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopOpenEdge)
-
-        .def("nodesRightBottomOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomOpenEdge)
-
-        .def("nodesRightTopOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopOpenEdge)
-
-        .def("nodesRightFrontOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontOpenEdge)
-
-        .def("nodesRightBackOpenEdge", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackOpenEdge)
-
-        .def(
-            "nodesFrontBottomLeftCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomLeftCorner)
-
-        .def(
-            "nodesFrontBottomRightCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontBottomRightCorner)
-
-        .def("nodesFrontTopLeftCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopLeftCorner)
-
-        .def("nodesFrontTopRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontTopRightCorner)
-
-        .def(
-            "nodesBackBottomLeftCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomLeftCorner)
-
-        .def(
-            "nodesBackBottomRightCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBackBottomRightCorner)
-
-        .def("nodesBackTopLeftCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopLeftCorner)
-
-        .def("nodesBackTopRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackTopRightCorner)
-
-        .def(
-            "nodesFrontLeftBottomCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftBottomCorner)
-
-        .def(
-            "nodesBottomFrontLeftCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontLeftCorner)
-
-        .def(
-            "nodesBottomLeftFrontCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftFrontCorner)
-
-        .def(
-            "nodesLeftFrontBottomCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontBottomCorner)
-
-        .def(
-            "nodesLeftBottomFrontCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomFrontCorner)
-
-        .def(
-            "nodesFrontRightBottomCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightBottomCorner)
-
-        .def(
-            "nodesBottomFrontRightCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomFrontRightCorner)
-
-        .def(
-            "nodesBottomRightFrontCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightFrontCorner)
-
-        .def(
-            "nodesRightFrontBottomCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontBottomCorner)
-
-        .def(
-            "nodesRightBottomFrontCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomFrontCorner)
-
-        .def("nodesFrontLeftTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontLeftTopCorner)
-
-        .def("nodesTopFrontLeftCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontLeftCorner)
-
-        .def("nodesTopLeftFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftFrontCorner)
-
-        .def("nodesLeftFrontTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftFrontTopCorner)
-
-        .def("nodesLeftTopFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopFrontCorner)
-
-        .def("nodesFrontRightTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesFrontRightTopCorner)
-
-        .def("nodesTopFrontRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopFrontRightCorner)
-
-        .def("nodesTopRightFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightFrontCorner)
-
-        .def("nodesRightFrontTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightFrontTopCorner)
-
-        .def("nodesRightTopFrontCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopFrontCorner)
-
-        .def(
-            "nodesBackLeftBottomCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftBottomCorner)
-
-        .def(
-            "nodesBottomBackLeftCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackLeftCorner)
-
-        .def(
-            "nodesBottomLeftBackCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomLeftBackCorner)
-
-        .def(
-            "nodesLeftBackBottomCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackBottomCorner)
-
-        .def(
-            "nodesLeftBottomBackCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBottomBackCorner)
-
-        .def(
-            "nodesBackRightBottomCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightBottomCorner)
-
-        .def(
-            "nodesBottomBackRightCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomBackRightCorner)
-
-        .def(
-            "nodesBottomRightBackCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesBottomRightBackCorner)
-
-        .def(
-            "nodesRightBackBottomCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackBottomCorner)
-
-        .def(
-            "nodesRightBottomBackCorner",
-            &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBottomBackCorner)
-
-        .def("nodesBackLeftTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackLeftTopCorner)
-
-        .def("nodesTopBackLeftCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackLeftCorner)
-
-        .def("nodesTopLeftBackCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopLeftBackCorner)
-
-        .def("nodesLeftBackTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftBackTopCorner)
-
-        .def("nodesLeftTopBackCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesLeftTopBackCorner)
-
-        .def("nodesBackRightTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesBackRightTopCorner)
-
-        .def("nodesTopBackRightCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopBackRightCorner)
-
-        .def("nodesTopRightBackCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesTopRightBackCorner)
-
-        .def("nodesRightBackTopCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightBackTopCorner)
-
-        .def("nodesRightTopBackCorner", &GooseFEM::Mesh::Hex8::FineLayer::nodesRightTopBackCorner)
-
-        .def("nodesPeriodic", &GooseFEM::Mesh::Hex8::FineLayer::nodesPeriodic)
-
-        .def("nodesOrigin", &GooseFEM::Mesh::Hex8::FineLayer::nodesOrigin)
-
-        .def("dofs", &GooseFEM::Mesh::Hex8::FineLayer::dofs)
-
-        .def("dofsPeriodic", &GooseFEM::Mesh::Hex8::FineLayer::dofsPeriodic)
-
-        .def("__repr__", [](const GooseFEM::Mesh::Hex8::FineLayer&) {
+        cls.def("__repr__", [](const GooseFEM::Mesh::Hex8::FineLayer&) {
             return "<GooseFEM.Mesh.Hex8.FineLayer>";
         });
+    }
 }
+
+#endif
diff --git a/python/MeshQuad4.hpp b/python/MeshQuad4.hpp
index 04de340..158ff49 100644
--- a/python/MeshQuad4.hpp
+++ b/python/MeshQuad4.hpp
@@ -1,250 +1,199 @@
-/* =================================================================================================
+/**
+\file
+\copyright Copyright 2017. Tom de Geus. All rights reserved.
+\license This project is released under the GNU Public License (GPLv3).
+*/
 
-(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
+#ifndef PYGOOSEFEM_MESHQUAD4_H
+#define PYGOOSEFEM_MESHQUAD4_H
 
-================================================================================================= */
-
-#include <GooseFEM/GooseFEM.h>
+#include <GooseFEM/MeshQuad4.h>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 #include <xtensor-python/pyarray.hpp>
 #include <xtensor-python/pytensor.hpp>
 #include <pyxtensor/pyxtensor.hpp>
 
+#include "Mesh.hpp"
+
 namespace py = pybind11;
 
 void init_MeshQuad4(py::module& m)
 {
+    {
+        py::class_<GooseFEM::Mesh::Quad4::Regular> cls(m, "Regular");
 
-    py::class_<GooseFEM::Mesh::Quad4::Regular>(m, "Regular")
-
-        .def(
-            py::init<size_t, size_t, double>(),
-            "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge "
-            "size 'h'",
-            py::arg("nx"),
-            py::arg("ny"),
-            py::arg("h") = 1.)
-
-        .def("coor", &GooseFEM::Mesh::Quad4::Regular::coor)
-        .def("conn", &GooseFEM::Mesh::Quad4::Regular::conn)
-        .def("nelem", &GooseFEM::Mesh::Quad4::Regular::nelem)
-        .def("nnode", &GooseFEM::Mesh::Quad4::Regular::nnode)
-        .def("nne", &GooseFEM::Mesh::Quad4::Regular::nne)
-        .def("ndim", &GooseFEM::Mesh::Quad4::Regular::ndim)
-        .def("nelx", &GooseFEM::Mesh::Quad4::Regular::nelx)
-        .def("nely", &GooseFEM::Mesh::Quad4::Regular::nely)
-        .def("h", &GooseFEM::Mesh::Quad4::Regular::h)
-        .def("getElementType", &GooseFEM::Mesh::Quad4::Regular::getElementType)
-        .def("nodesBottomEdge", &GooseFEM::Mesh::Quad4::Regular::nodesBottomEdge)
-        .def("nodesTopEdge", &GooseFEM::Mesh::Quad4::Regular::nodesTopEdge)
-        .def("nodesLeftEdge", &GooseFEM::Mesh::Quad4::Regular::nodesLeftEdge)
-        .def("nodesRightEdge", &GooseFEM::Mesh::Quad4::Regular::nodesRightEdge)
-        .def("nodesBottomOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesBottomOpenEdge)
-        .def("nodesTopOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesTopOpenEdge)
-        .def("nodesLeftOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesLeftOpenEdge)
-        .def("nodesRightOpenEdge", &GooseFEM::Mesh::Quad4::Regular::nodesRightOpenEdge)
-        .def("nodesBottomLeftCorner", &GooseFEM::Mesh::Quad4::Regular::nodesBottomLeftCorner)
-        .def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::Regular::nodesBottomRightCorner)
-        .def("nodesTopLeftCorner", &GooseFEM::Mesh::Quad4::Regular::nodesTopLeftCorner)
-        .def("nodesTopRightCorner", &GooseFEM::Mesh::Quad4::Regular::nodesTopRightCorner)
-        .def("nodesLeftBottomCorner", &GooseFEM::Mesh::Quad4::Regular::nodesLeftBottomCorner)
-        .def("nodesLeftTopCorner", &GooseFEM::Mesh::Quad4::Regular::nodesLeftTopCorner)
-        .def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::Regular::nodesRightBottomCorner)
-        .def("nodesRightTopCorner", &GooseFEM::Mesh::Quad4::Regular::nodesRightTopCorner)
-        .def("dofs", &GooseFEM::Mesh::Quad4::Regular::dofs)
-        .def("nodesPeriodic", &GooseFEM::Mesh::Quad4::Regular::nodesPeriodic)
-        .def("nodesOrigin", &GooseFEM::Mesh::Quad4::Regular::nodesOrigin)
-        .def("dofsPeriodic", &GooseFEM::Mesh::Quad4::Regular::dofsPeriodic)
-        .def("elementgrid", &GooseFEM::Mesh::Quad4::Regular::elementgrid)
-
-        .def("__repr__", [](const GooseFEM::Mesh::Quad4::Regular&) {
-            return "<GooseFEM.Mesh.Quad4.Regular>";
-        });
+        cls.def(py::init<size_t, size_t, double>(),
+                "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge "
+                "size 'h'",
+                py::arg("nx"),
+                py::arg("ny"),
+                py::arg("h") = 1.0);
 
-    py::class_<GooseFEM::Mesh::Quad4::FineLayer>(m, "FineLayer")
+        register_Element_RegularBase<GooseFEM::Mesh::Quad4::Regular>(cls);
+        register_Element_RegularBase2d<GooseFEM::Mesh::Quad4::Regular>(cls);
 
-        .def(
-            py::init<size_t, size_t, double, size_t>(),
-            "FineLayer mesh: 'nx' pixels in horizontal direction (length 'Lx'), idem in vertical "
-            "direction",
-            py::arg("nx"),
-            py::arg("ny"),
-            py::arg("h") = 1.,
-            py::arg("nfine") = 1)
+        cls.def("elementgrid", &GooseFEM::Mesh::Quad4::Regular::elementgrid);
 
-        .def(
-            py::init<const xt::pytensor<double, 2>&, const xt::pytensor<size_t, 2>&>(),
-            "Map connectivity to generating FineLayer-object.",
-            py::arg("coor"),
-            py::arg("conn"))
-
-        .def("coor", &GooseFEM::Mesh::Quad4::FineLayer::coor)
-        .def("conn", &GooseFEM::Mesh::Quad4::FineLayer::conn)
-        .def("nelem", &GooseFEM::Mesh::Quad4::FineLayer::nelem)
-        .def("nnode", &GooseFEM::Mesh::Quad4::FineLayer::nnode)
-        .def("nne", &GooseFEM::Mesh::Quad4::FineLayer::nne)
-        .def("ndim", &GooseFEM::Mesh::Quad4::FineLayer::ndim)
-        .def("nelx", &GooseFEM::Mesh::Quad4::FineLayer::nelx)
-        .def("nely", &GooseFEM::Mesh::Quad4::FineLayer::nely)
-        .def("h", &GooseFEM::Mesh::Quad4::FineLayer::h)
-        .def("elemrow_nhx", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nhx)
-        .def("elemrow_nhy", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nhy)
-        .def("elemrow_nelem", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nelem)
-        .def("getElementType", &GooseFEM::Mesh::Quad4::FineLayer::getElementType)
-        .def("elementsMiddleLayer", &GooseFEM::Mesh::Quad4::FineLayer::elementsMiddleLayer)
-        .def("elementsLayer", &GooseFEM::Mesh::Quad4::FineLayer::elementsLayer)
-
-        .def(
-            "elementgrid_ravel",
-            &GooseFEM::Mesh::Quad4::FineLayer::elementgrid_ravel,
-            py::arg("rows_range"),
-            py::arg("cols_range"))
-
-        .def(
-            "elementgrid_around_ravel",
-            &GooseFEM::Mesh::Quad4::FineLayer::elementgrid_around_ravel,
-            py::arg("element"),
-            py::arg("size"),
-            py::arg("periodic") = true)
-
-        .def(
-            "elementgrid_leftright",
-            &GooseFEM::Mesh::Quad4::FineLayer::elementgrid_leftright,
-            py::arg("element"),
-            py::arg("left"),
-            py::arg("right"),
-            py::arg("periodic") = true)
-
-        .def("nodesBottomEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomEdge)
-        .def("nodesTopEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopEdge)
-        .def("nodesLeftEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftEdge)
-        .def("nodesRightEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightEdge)
-        .def("nodesBottomOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomOpenEdge)
-        .def("nodesTopOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopOpenEdge)
-        .def("nodesLeftOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftOpenEdge)
-        .def("nodesRightOpenEdge", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightOpenEdge)
-        .def("nodesBottomLeftCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomLeftCorner)
-        .def("nodesBottomRightCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesBottomRightCorner)
-        .def("nodesTopLeftCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopLeftCorner)
-        .def("nodesTopRightCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesTopRightCorner)
-        .def("nodesLeftBottomCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftBottomCorner)
-        .def("nodesLeftTopCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesLeftTopCorner)
-        .def("nodesRightBottomCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightBottomCorner)
-        .def("nodesRightTopCorner", &GooseFEM::Mesh::Quad4::FineLayer::nodesRightTopCorner)
-        .def("dofs", &GooseFEM::Mesh::Quad4::FineLayer::dofs)
-        .def("nodesPeriodic", &GooseFEM::Mesh::Quad4::FineLayer::nodesPeriodic)
-        .def("nodesOrigin", &GooseFEM::Mesh::Quad4::FineLayer::nodesOrigin)
-        .def("dofsPeriodic", &GooseFEM::Mesh::Quad4::FineLayer::dofsPeriodic)
-        .def("roll", &GooseFEM::Mesh::Quad4::FineLayer::roll)
-
-        .def("__repr__", [](const GooseFEM::Mesh::Quad4::FineLayer&) {
+        cls.def("__repr__", [](const GooseFEM::Mesh::Quad4::Regular&) {
+            return "<GooseFEM.Mesh.Quad4.Regular>";
+        });
+    }
+
+    {
+        py::class_<GooseFEM::Mesh::Quad4::FineLayer> cls(m, "FineLayer");
+
+        cls.def(py::init<size_t, size_t, double, size_t>(),
+                "FineLayer mesh: 'nx' pixels in horizontal direction (length 'Lx'), idem in vertical direction",
+                py::arg("nx"),
+                py::arg("ny"),
+                py::arg("h") = 1.,
+                py::arg("nfine") = 1);
+
+        cls.def(py::init<const xt::pytensor<double, 2>&, const xt::pytensor<size_t, 2>&>(),
+                "Map connectivity to generating FineLayer-object.",
+                py::arg("coor"),
+                py::arg("conn"));
+
+        register_Element_RegularBase<GooseFEM::Mesh::Quad4::FineLayer>(cls);
+        register_Element_RegularBase2d<GooseFEM::Mesh::Quad4::FineLayer>(cls);
+
+        cls.def("elemrow_nhx", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nhx);
+        cls.def("elemrow_nhy", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nhy);
+        cls.def("elemrow_nelem", &GooseFEM::Mesh::Quad4::FineLayer::elemrow_nelem);
+        cls.def("elementsMiddleLayer", &GooseFEM::Mesh::Quad4::FineLayer::elementsMiddleLayer);
+        cls.def("elementsLayer", &GooseFEM::Mesh::Quad4::FineLayer::elementsLayer);
+
+        cls.def("elementgrid_ravel",
+                &GooseFEM::Mesh::Quad4::FineLayer::elementgrid_ravel,
+                py::arg("rows_range"),
+                py::arg("cols_range"));
+
+        cls.def("elementgrid_around_ravel",
+                &GooseFEM::Mesh::Quad4::FineLayer::elementgrid_around_ravel,
+                py::arg("element"),
+                py::arg("size"),
+                py::arg("periodic") = true);
+
+        cls.def("elementgrid_leftright",
+                &GooseFEM::Mesh::Quad4::FineLayer::elementgrid_leftright,
+                py::arg("element"),
+                py::arg("left"),
+                py::arg("right"),
+                py::arg("periodic") = true);
+
+        cls.def("roll", &GooseFEM::Mesh::Quad4::FineLayer::roll);
+
+        cls.def("__repr__", [](const GooseFEM::Mesh::Quad4::FineLayer&) {
             return "<GooseFEM.Mesh.Quad4.FineLayer>";
         });
+    }
 }
 
 void init_MeshQuad4Map(py::module& m)
 {
 
     py::class_<GooseFEM::Mesh::Quad4::Map::RefineRegular>(m, "RefineRegular")
 
         .def(
             py::init<const GooseFEM::Mesh::Quad4::Regular&, size_t, size_t>(),
             "Refine a regular mesh",
             py::arg("mesh"),
             py::arg("nx"),
             py::arg("ny"))
 
         .def("getCoarseMesh", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getCoarseMesh)
 
         .def("getFineMesh", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getFineMesh)
 
         .def("getMap", &GooseFEM::Mesh::Quad4::Map::RefineRegular::getMap)
 
         .def("meanToCoarse",
              py::overload_cast<const xt::xtensor<double, 1>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::meanToCoarse<double, 1>, py::const_))
 
         .def("meanToCoarse",
              py::overload_cast<const xt::xtensor<double, 2>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::meanToCoarse<double, 2>, py::const_))
 
         .def("meanToCoarse",
              py::overload_cast<const xt::xtensor<double, 3>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::meanToCoarse<double, 3>, py::const_))
 
         .def("meanToCoarse",
              py::overload_cast<const xt::xtensor<double, 4>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::meanToCoarse<double, 4>, py::const_))
 
         .def("averageToCoarse",
              py::overload_cast<const xt::xtensor<double, 1>&, const xt::xtensor<double, 1>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::averageToCoarse<double, 1, double>, py::const_))
 
         .def("averageToCoarse",
              py::overload_cast<const xt::xtensor<double, 2>&, const xt::xtensor<double, 2>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::averageToCoarse<double, 2, double>, py::const_))
 
         .def("averageToCoarse",
              py::overload_cast<const xt::xtensor<double, 3>&, const xt::xtensor<double, 3>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::averageToCoarse<double, 3, double>, py::const_))
 
         .def("averageToCoarse",
              py::overload_cast<const xt::xtensor<double, 4>&, const xt::xtensor<double, 4>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::averageToCoarse<double, 4, double>, py::const_))
 
         .def("mapToFine",
              py::overload_cast<const xt::xtensor<double, 1>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine<double, 1>, py::const_))
 
         .def("mapToFine",
              py::overload_cast<const xt::xtensor<double, 2>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine<double, 2>, py::const_))
 
         .def("mapToFine",
              py::overload_cast<const xt::xtensor<double, 3>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine<double, 3>, py::const_))
 
         .def("mapToFine",
              py::overload_cast<const xt::xtensor<double, 4>&>(
                 &GooseFEM::Mesh::Quad4::Map::RefineRegular::mapToFine<double, 4>, py::const_))
 
         .def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::RefineRegular&) {
             return "<GooseFEM.Mesh.Quad4.Map.RefineRegular>";
         });
 
     py::class_<GooseFEM::Mesh::Quad4::Map::FineLayer2Regular>(m, "FineLayer2Regular")
 
         .def(
             py::init<const GooseFEM::Mesh::Quad4::FineLayer&>(),
             "Map a FineLayer mesh to a Regular mesh",
             py::arg("mesh"))
 
         .def("getRegularMesh", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getRegularMesh)
 
         .def("getFineLayerMesh", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getFineLayerMesh)
 
         .def("getMap", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMap)
 
         .def("getMapFraction", &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::getMapFraction)
 
         .def("mapToRegular",
              py::overload_cast<const xt::xtensor<double, 1>&>(
                 &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular<double, 1>, py::const_))
 
         .def("mapToRegular",
              py::overload_cast<const xt::xtensor<double, 2>&>(
                 &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular<double, 2>, py::const_))
 
         .def("mapToRegular",
              py::overload_cast<const xt::xtensor<double, 3>&>(
                 &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular<double, 3>, py::const_))
 
         .def("mapToRegular",
              py::overload_cast<const xt::xtensor<double, 4>&>(
                 &GooseFEM::Mesh::Quad4::Map::FineLayer2Regular::mapToRegular<double, 4>, py::const_))
 
         .def("__repr__", [](const GooseFEM::Mesh::Quad4::Map::FineLayer2Regular&) {
             return "<GooseFEM.Mesh.Quad4.Map.FineLayer2Regular>";
         });
 }
+
+#endif
diff --git a/python/MeshTri3.hpp b/python/MeshTri3.hpp
index b58213c..93e9655 100644
--- a/python/MeshTri3.hpp
+++ b/python/MeshTri3.hpp
@@ -1,114 +1,65 @@
-/* =================================================================================================
+/**
+\file
+\copyright Copyright 2017. Tom de Geus. All rights reserved.
+\license This project is released under the GNU Public License (GPLv3).
+*/
 
-(c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
+#ifndef PYGOOSEFEM_MESHTRI3_H
+#define PYGOOSEFEM_MESHTRI3_H
 
-================================================================================================= */
-
-#include <GooseFEM/GooseFEM.h>
+#include <GooseFEM/MeshTri3.h>
 #include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
 #include <xtensor-python/pyarray.hpp>
 #include <xtensor-python/pytensor.hpp>
 #include <pyxtensor/pyxtensor.hpp>
 
+#include "Mesh.hpp"
+
 namespace py = pybind11;
 
-void init_MeshTri3(py::module& m)
+void init_MeshTri3(py::module& mod)
 {
+    py::class_<GooseFEM::Mesh::Tri3::Regular> cls(mod, "Regular");
 
-    py::class_<GooseFEM::Mesh::Tri3::Regular>(m, "Regular")
-
-        .def(
-            py::init<size_t, size_t, double>(),
-            "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge "
-            "size 'h'",
+    cls.def(py::init<size_t, size_t, double>(),
+            "Regular mesh: 'nx' pixels in horizontal direction, 'ny' in vertical direction, edge size 'h'",
             py::arg("nx"),
             py::arg("ny"),
-            py::arg("h") = 1.)
-
-        .def("coor", &GooseFEM::Mesh::Tri3::Regular::coor)
-
-        .def("conn", &GooseFEM::Mesh::Tri3::Regular::conn)
-
-        .def("nelem", &GooseFEM::Mesh::Tri3::Regular::nelem)
-
-        .def("nnode", &GooseFEM::Mesh::Tri3::Regular::nnode)
-
-        .def("nne", &GooseFEM::Mesh::Tri3::Regular::nne)
-
-        .def("ndim", &GooseFEM::Mesh::Tri3::Regular::ndim)
-
-        .def("getElementType", &GooseFEM::Mesh::Tri3::Regular::getElementType)
-
-        .def("nodesBottomEdge", &GooseFEM::Mesh::Tri3::Regular::nodesBottomEdge)
-
-        .def("nodesTopEdge", &GooseFEM::Mesh::Tri3::Regular::nodesTopEdge)
-
-        .def("nodesLeftEdge", &GooseFEM::Mesh::Tri3::Regular::nodesLeftEdge)
-
-        .def("nodesRightEdge", &GooseFEM::Mesh::Tri3::Regular::nodesRightEdge)
-
-        .def("nodesBottomOpenEdge", &GooseFEM::Mesh::Tri3::Regular::nodesBottomOpenEdge)
-
-        .def("nodesTopOpenEdge", &GooseFEM::Mesh::Tri3::Regular::nodesTopOpenEdge)
-
-        .def("nodesLeftOpenEdge", &GooseFEM::Mesh::Tri3::Regular::nodesLeftOpenEdge)
-
-        .def("nodesRightOpenEdge", &GooseFEM::Mesh::Tri3::Regular::nodesRightOpenEdge)
-
-        .def("nodesBottomLeftCorner", &GooseFEM::Mesh::Tri3::Regular::nodesBottomLeftCorner)
-
-        .def("nodesBottomRightCorner", &GooseFEM::Mesh::Tri3::Regular::nodesBottomRightCorner)
-
-        .def("nodesTopLeftCorner", &GooseFEM::Mesh::Tri3::Regular::nodesTopLeftCorner)
-
-        .def("nodesTopRightCorner", &GooseFEM::Mesh::Tri3::Regular::nodesTopRightCorner)
-
-        .def("nodesLeftBottomCorner", &GooseFEM::Mesh::Tri3::Regular::nodesLeftBottomCorner)
-
-        .def("nodesLeftTopCorner", &GooseFEM::Mesh::Tri3::Regular::nodesLeftTopCorner)
-
-        .def("nodesRightBottomCorner", &GooseFEM::Mesh::Tri3::Regular::nodesRightBottomCorner)
-
-        .def("nodesRightTopCorner", &GooseFEM::Mesh::Tri3::Regular::nodesRightTopCorner)
-
-        .def("nodesPeriodic", &GooseFEM::Mesh::Tri3::Regular::nodesPeriodic)
-
-        .def("nodesOrigin", &GooseFEM::Mesh::Tri3::Regular::nodesOrigin)
-
-        .def("dofs", &GooseFEM::Mesh::Tri3::Regular::dofs)
-
-        .def("dofsPeriodic", &GooseFEM::Mesh::Tri3::Regular::dofsPeriodic)
-
-        .def("__repr__", [](const GooseFEM::Mesh::Tri3::Regular&) {
-            return "<GooseFEM.Mesh.Tri3.Regular>";
-        });
-
-    m.def(
-        "getOrientation",
-        &GooseFEM::Mesh::Tri3::getOrientation,
-        "Get the orientation of each element",
-        py::arg("coor"),
-        py::arg("conn"));
-
-    m.def(
-        "setOrientation",
-        py::overload_cast<const xt::xtensor<double, 2>&, const xt::xtensor<size_t, 2>&, int>(
-            &GooseFEM::Mesh::Tri3::setOrientation),
-        "Set the orientation of each element",
-        py::arg("coor"),
-        py::arg("conn"),
-        py::arg("orientation"));
-
-    m.def(
-        "setOrientation",
-        py::overload_cast<
-            const xt::xtensor<double, 2>&,
-            const xt::xtensor<size_t, 2>&,
-            const xt::xtensor<int, 1>&,
-            int>(&GooseFEM::Mesh::Tri3::setOrientation),
-        "Set the orientation of each element",
-        py::arg("coor"),
-        py::arg("conn"),
-        py::arg("val"),
-        py::arg("orientation"));
+            py::arg("h") = 1.);
+
+    register_Element_RegularBase<GooseFEM::Mesh::Tri3::Regular>(cls);
+    register_Element_RegularBase2d<GooseFEM::Mesh::Tri3::Regular>(cls);
+
+    cls.def("__repr__", [](const GooseFEM::Mesh::Tri3::Regular&) {
+        return "<GooseFEM.Mesh.Tri3.Regular>";
+    });
+
+    mod.def("getOrientation",
+            &GooseFEM::Mesh::Tri3::getOrientation,
+            "Get the orientation of each element",
+            py::arg("coor"),
+            py::arg("conn"));
+
+    mod.def("setOrientation",
+            py::overload_cast<const xt::xtensor<double, 2>&, const xt::xtensor<size_t, 2>&, int>(
+                &GooseFEM::Mesh::Tri3::setOrientation),
+            "Set the orientation of each element",
+            py::arg("coor"),
+            py::arg("conn"),
+            py::arg("orientation"));
+
+    mod.def("setOrientation",
+            py::overload_cast<
+                const xt::xtensor<double, 2>&,
+                const xt::xtensor<size_t, 2>&,
+                const xt::xtensor<int, 1>&,
+                int>(&GooseFEM::Mesh::Tri3::setOrientation),
+            "Set the orientation of each element",
+            py::arg("coor"),
+            py::arg("conn"),
+            py::arg("val"),
+            py::arg("orientation"));
 }
+
+#endif