diff --git a/include/GooseFEM/MeshQuad4.h b/include/GooseFEM/MeshQuad4.h
index 37695bf..139c2df 100644
--- a/include/GooseFEM/MeshQuad4.h
+++ b/include/GooseFEM/MeshQuad4.h
@@ -1,527 +1,559 @@
 /**
 Generate mesh with 4-noded quadrilateral elements.
 
 \file MeshQuad4.h
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_MESHQUAD4_H
 #define GOOSEFEM_MESHQUAD4_H
 
 #include "config.h"
 
 namespace GooseFEM {
 namespace Mesh {
 namespace Quad4 {
 
 // pre-allocation
 
 namespace Map {
     class FineLayer2Regular;
 }
 
 /**
 Regular mesh: equi-sized elements.
 */
 class Regular {
 public:
 
     Regular() = default;
 
     /**
     Constructor.
 
     \param nelx Number of elements in horizontal (x) direction.
     \param nely Number of elements in vertical (y) direction.
     \param h Edge size (width == height).
     */
     Regular(size_t nelx, size_t nely, double h = 1.0);
 
     /**
     Number of elements.
 
     \return unsigned int.
     */
     size_t nelem() const;
 
     /**
     Number of nodes.
 
     \return unsigned int.
     */
     size_t nnode() const;
 
     /**
     Number of nodes-per-element.
 
     \return unsigned int.
     */
     size_t nne() const;
 
     /**
     Number of dimensions.
 
     \return unsigned int.
     */
     size_t ndim() const;
 
     /**
     Number of elements in x-direction == width of the mesh in units of h().
 
     \return unsigned int.
     */
     size_t nelx() const;
 
     /**
     Number of elements in y-direction == height of the mesh, in units of h(),
 
     \return unsigned int.
     */
     size_t nely() const;
 
     /**
     Edge size of one element.
 
     \return double.
     */
     double h() const;
 
     /**
     Element type.
 
     \return GooseFEM::Mesh::ElementType().
     */
     ElementType getElementType() const;
 
     /**
     Nodal coordinates.
 
     \return ``[nnode, ndim]``.
     */
     xt::xtensor<double, 2> coor() const;
 
     /**
     Connectivity.
 
     \return ``[nelem, nne]``.
     */
     xt::xtensor<size_t, 2> conn() const;
 
     // boundary nodes: edges
     xt::xtensor<size_t, 1> nodesBottomEdge() const;
     xt::xtensor<size_t, 1> nodesTopEdge() const;
     xt::xtensor<size_t, 1> nodesLeftEdge() const;
     xt::xtensor<size_t, 1> nodesRightEdge() const;
 
     // boundary nodes: edges, without corners
     xt::xtensor<size_t, 1> nodesBottomOpenEdge() const;
     xt::xtensor<size_t, 1> nodesTopOpenEdge() const;
     xt::xtensor<size_t, 1> nodesLeftOpenEdge() const;
     xt::xtensor<size_t, 1> nodesRightOpenEdge() const;
 
     // boundary nodes: corners (including aliases)
     size_t nodesBottomLeftCorner() const;
     size_t nodesBottomRightCorner() const;
     size_t nodesTopLeftCorner() const;
     size_t nodesTopRightCorner() const;
     size_t nodesLeftBottomCorner() const;
     size_t nodesLeftTopCorner() const;
     size_t nodesRightBottomCorner() const;
     size_t nodesRightTopCorner() const;
 
     // DOF-numbers for each component of each node (sequential)
     xt::xtensor<size_t, 2> dofs() const;
 
     // DOF-numbers for the case that the periodicity if fully eliminated
     xt::xtensor<size_t, 2> dofsPeriodic() const;
 
     // periodic node pairs [:,2]: (independent, dependent)
     xt::xtensor<size_t, 2> nodesPeriodic() const;
 
     // front-bottom-left node, used as reference for periodicity
     size_t nodesOrigin() const;
 
     // element numbers as matrix
     xt::xtensor<size_t, 2> elementgrid() const;
 
 private:
     double m_h;                     // elementary element edge-size (in all directions)
     size_t m_nelx;                  // number of elements in x-direction (length == "m_nelx * m_h")
     size_t m_nely;                  // number of elements in y-direction (length == "m_nely * m_h")
     size_t m_nelem;                 // number of elements
     size_t m_nnode;                 // number of nodes
     static const size_t m_nne = 4;  // number of nodes-per-element
     static const size_t m_ndim = 2; // number of dimensions
 };
 
 // Mesh with fine middle layer, and coarser elements towards the top and bottom
 
 class FineLayer {
 public:
 
     FineLayer() = default;
 
     FineLayer(size_t nelx, size_t nely, double h = 1.0, size_t nfine = 1);
 
     /**
     Reconstruct class for given coordinates / connectivity.
 
     \param coor Nodal coordinates ``[nnode, ndim]`` with ``ndim == 2``.
     \param conn Connectivity ``[nne, nne]`` with ``nne == 4``.
     \throw GOOSEFEM_CHECK()
     */
     FineLayer(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn);
 
     /**
     Number of elements.
 
     \return unsigned int.
     */
     size_t nelem() const;
 
     /**
     Number of nodes.
 
     \return unsigned int.
     */
     size_t nnode() const;
 
     /**
     Number of nodes-per-element.
 
     \return unsigned int.
     */
     size_t nne() const;
 
     /**
     Number of dimensions.
 
     \return unsigned int.
     */
     size_t ndim() const;
 
     /**
     Number of elements in x-direction along the middle layer == width of the mesh in units of h().
 
     \return unsigned int.
     */
     size_t nelx() const;
 
     /**
     Height of the mesh, in units of h()
 
     \return unsigned int.
     */
     size_t nely() const;
 
     /**
     Edge size of the smallest elements (along the middle layer).
 
     \return double.
     */
     double h() const;
 
     // edge size, per row of elements (in units of "h")
     xt::xtensor<size_t, 1> elemrow_nhx() const;
     xt::xtensor<size_t, 1> elemrow_nhy() const;
     xt::xtensor<size_t, 1> elemrow_nelem() const;
 
     /**
     Element type.
 
     \return GooseFEM::Mesh::ElementType().
     */
     ElementType getElementType() const;
 
     /**
     Nodal coordinates.
 
     \return ``[nnode, ndim]``.
     */
     xt::xtensor<double, 2> coor() const;
 
     /**
     Connectivity.
 
     \return ``[nelem, nne]``.
     */
     xt::xtensor<size_t, 2> conn() const;
 
     // elements in the middle (fine) layer
     xt::xtensor<size_t, 1> elementsMiddleLayer() const;
 
     // extract elements along a layer
     xt::xtensor<size_t, 1> elementsLayer(size_t layer) const;
 
     // select region of elements from 'matrix' of element numbers
     xt::xtensor<size_t, 1> elementgrid_ravel(
         std::vector<size_t> rows_start_stop,
         std::vector<size_t> cols_start_stop) const;
 
     /**
     Select region of elements from 'matrix' of element numbers around an element:
     square box with edge-size ``(2 * size + 1) * h``, around ``element``.
 
     \param element The element around which to select elements.
     \param size Edge size of the square box encapsulating the selected element.
     \param periodic Assume the mesh periodic.
     \returns List of elements.
     */
     xt::xtensor<size_t, 1> elementgrid_around_ravel(
         size_t element,
         size_t size,
         bool periodic = true);
 
     /**
     Select region of elements from 'matrix' of element numbers around an element:
     left/right from ``element`` (on the same layer).
 
     \param element The element around which to select elements.
     \param left Number of elements to select to the left.
     \param right Number of elements to select to the right.
     \param periodic Assume the mesh periodic.
     \returns List of elements.
     */
     // -
     xt::xtensor<size_t, 1> elementgrid_leftright(
         size_t element,
         size_t left,
         size_t right,
         bool periodic = true);
 
     // boundary nodes: edges
     xt::xtensor<size_t, 1> nodesBottomEdge() const;
     xt::xtensor<size_t, 1> nodesTopEdge() const;
     xt::xtensor<size_t, 1> nodesLeftEdge() const;
     xt::xtensor<size_t, 1> nodesRightEdge() const;
 
     // boundary nodes: edges, without corners
     xt::xtensor<size_t, 1> nodesBottomOpenEdge() const;
     xt::xtensor<size_t, 1> nodesTopOpenEdge() const;
     xt::xtensor<size_t, 1> nodesLeftOpenEdge() const;
     xt::xtensor<size_t, 1> nodesRightOpenEdge() const;
 
     // boundary nodes: corners (including aliases)
     size_t nodesBottomLeftCorner() const;
     size_t nodesBottomRightCorner() const;
     size_t nodesTopLeftCorner() const;
     size_t nodesTopRightCorner() const;
     size_t nodesLeftBottomCorner() const;
     size_t nodesLeftTopCorner() const;
     size_t nodesRightBottomCorner() const;
     size_t nodesRightTopCorner() const;
 
     // DOF-numbers for each component of each node (sequential)
     xt::xtensor<size_t, 2> dofs() const;
 
     // DOF-numbers for the case that the periodicity if fully eliminated
     xt::xtensor<size_t, 2> dofsPeriodic() const;
 
     // periodic node pairs [:,2]: (independent, dependent)
     xt::xtensor<size_t, 2> nodesPeriodic() const;
 
     // front-bottom-left node, used as reference for periodicity
     size_t nodesOrigin() const;
 
     // mapping to 'roll' periodically in the x-direction,
     // returns element mapping, such that: new_elemvar = elemvar[elem_map]
     xt::xtensor<size_t, 1> roll(size_t n);
 
 private:
     double m_h;                         // elementary element edge-size (in all directions)
     double m_Lx;                        // mesh size in "x"
     size_t m_nelem;                     // number of elements
     size_t m_nnode;                     // number of nodes
     static const size_t m_nne = 4;      // number of nodes-per-element
     static const size_t m_ndim = 2;     // number of dimensions
     xt::xtensor<size_t, 1> m_nelx;      // number of elements in "x" (*)
     xt::xtensor<size_t, 1> m_nnd;       // total number of nodes in the main node layer (**)
     xt::xtensor<size_t, 1> m_nhx;       // element size in x-direction (*)
     xt::xtensor<size_t, 1> m_nhy;       // element size in y-direction (*)
     xt::xtensor<int, 1> m_refine;       // refine direction (-1:no refine, 0:"x" (*)
     xt::xtensor<size_t, 1> m_startElem; // start element (*)
     xt::xtensor<size_t, 1> m_startNode; // start node (**)
     // (*) per element layer in "y"
     // (**) per node layer in "y"
 
     void init(size_t nelx, size_t nely, double h, size_t nfine = 1);
     void map(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn);
 
     friend class GooseFEM::Mesh::Quad4::Map::FineLayer2Regular;
 };
 
 // Mesh mappings
 
 namespace Map {
 
     // Return "FineLayer"-class responsible for generating a connectivity
     // Throws if conversion is not possible
 
     [[ deprecated ]]
     GooseFEM::Mesh::Quad4::FineLayer FineLayer(
         const xt::xtensor<double, 2>& coor,
         const xt::xtensor<size_t, 2>& conn);
 
     /**
     Refine a Regular mesh: subdivide elements in several smaller elements.
     */
     class RefineRegular {
     public:
-        // Constructors
+
         RefineRegular() = default;
+
+        /**
+        Constructor.
+
+        \param mesh the coarse mesh.
+        \param nx for each coarse element: number of fine elements in x-direction.
+        \param ny for each coarse element: number of fine elements in y-direction.
+        */
         RefineRegular(const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny);
 
         /**
         For each coarse element: number of fine elements in x-direction.
 
         \return unsigned int (same as used in constructor)
         */
         size_t nx() const;
 
         /**
         For each coarse element: number of fine elements in y-direction.
 
         \return unsigned int (same as used in constructor)
         */
         size_t ny() const;
 
-        // return the coarse or the fine mesh objects
+        /**
+        Obtain the coarse mesh (copy of the mesh passed to the constructor).
+
+        \return mesh
+        */
         GooseFEM::Mesh::Quad4::Regular getCoarseMesh() const;
+
+        /**
+        Obtain the fine mesh.
+
+        \return mesh
+        */
         GooseFEM::Mesh::Quad4::Regular getFineMesh() const;
 
-        // elements of the fine mesh per element of the coarse mesh
+        /**
+        Get element-mapping: elements of the fine mesh per element of the coarse mesh.
+
+        \return [nelem_coarse, nx() * ny()]
+        */
         xt::xtensor<size_t, 2> getMap() const;
 
         // map field
-
         [[ deprecated ]]
         xt::xtensor<double, 2> mapToCoarse(const xt::xtensor<double, 1>& data) const; // scalar per el
 
         [[ deprecated ]]
         xt::xtensor<double, 2> mapToCoarse(const xt::xtensor<double, 2>& data) const; // scalar per intpnt
 
         [[ deprecated ]]
         xt::xtensor<double, 4> mapToCoarse(const xt::xtensor<double, 4>& data) const; // tensor per intpnt
 
         /**
         Compute the mean of the quantity define on the fine mesh when mapped on the coarse mesh.
 
+        \tparam T type of the data (e.g. ``double``).
+        \tparam rank rank of the data.
         \param data the data [nelem_fine, ...]
         \return the average data of the coarse mesh [nelem_coarse, ...]
         */
         template <class T, size_t rank>
         xt::xtensor<T, rank> meanToCoarse(const xt::xtensor<T, rank>& data) const;
 
         /**
         Compute the average of the quantity define on the fine mesh when mapped on the coarse mesh.
 
+        \tparam T type of the data (e.g. ``double``).
+        \tparam rank rank of the data.
+        \tparam S type of the weights (e.g. ``double``).
         \param data the data [nelem_fine, ...]
         \param weights the weights [nelem_fine, ...]
         \return the average data of the coarse mesh [nelem_coarse, ...]
         */
         template <class T, size_t rank, class S>
         xt::xtensor<T, rank> averageToCoarse(
             const xt::xtensor<T, rank>& data,
             const xt::xtensor<S, rank>& weights) const;
 
         /**
         Map element quantities to the fine mesh.
         The mapping is a bit simplistic: no interpolation is involved.
         The mapping is such that::
 
-            ret[e_fine, ...] <- arg[e_coarse, ...]
+            ret[e_fine, ...] <- data[e_coarse, ...]
 
         \tparam T type of the data (e.g. ``double``).
         \tparam rank rank of the data.
-        \param arg data.
+        \param data the data.
         \return mapped data.
         */
         template <class T, size_t rank>
         xt::xtensor<T, rank> mapToFine(const xt::xtensor<T, rank>& data) const;
 
     private:
 
         GooseFEM::Mesh::Quad4::Regular m_coarse; ///< the coarse mesh
         GooseFEM::Mesh::Quad4::Regular m_fine; ///< the fine mesh
-        size_t m_nx; ///< for each coarse element: number of fine elements in x-direction
-        size_t m_ny; ///< for each coarse element: number of fine elements in y-direction
-
-        /**
-        Fine elements for each coarse elements [nelem_coarse, m_nx * m_ny]
-        */
-        xt::xtensor<size_t, 2> m_coarse2fine;
+        size_t m_nx; ///< see nx()
+        size_t m_ny; ///< see ny()
+        xt::xtensor<size_t, 2> m_coarse2fine; ///< see getMap()
     };
 
     /**
     Map a FineLayer mesh to a Regular mesh.
-    The element size is based on the element size of the FineLayer mesh.
+    The element size of the Regular corresponds to the smallest elements of the FineLayer mesh
+    (along the middle layer).
     */
     class FineLayer2Regular {
     public:
 
         FineLayer2Regular() = default;
 
         /**
         Constructors.
 
         \param mesh The FineLayer mesh.
         */
         FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh);
 
         /**
-        Obtain Regular mesh.
+        Obtain the Regular mesh.
 
         \return mesh.
         */
         GooseFEM::Mesh::Quad4::Regular getRegularMesh() const;
 
         /**
-        Obtain FineLayer mesh (copy of the mesh passed to the constructor).
+        Obtain the FineLayer mesh (copy of the mesh passed to the constructor).
 
         \return mesh.
         */
         GooseFEM::Mesh::Quad4::FineLayer getFineLayerMesh() const;
 
         // elements of the Regular mesh per element of the FineLayer mesh
         // and the fraction by which the overlap is
+
+        /**
+        Get element-mapping: elements of the Regular mesh per element of the FineLayer mesh.
+        The number of Regular elements varies between elements of the FineLayer mesh.
+
+        \return [nelem_finelayer, ?]
+        */
         std::vector<std::vector<size_t>> getMap() const;
+
+        /**
+        To overlap fraction for each item in the mapping in getMap().
+
+        \return [nelem_finelayer, ?]
+        */
         std::vector<std::vector<double>> getMapFraction() const;
 
         /**
         Map element quantities to Regular.
         The mapping is a bit simplistic: no interpolation is involved, the function just
         accounts the fraction of overlap between the FineLayer element and the Regular element.
         The mapping is such that::
 
             ret[e_regular, ...] <- arg[e_finelayer, ...]
 
         \tparam T type of the data (e.g. ``double``).
         \tparam rank rank of the data.
         \param arg data.
         \return mapped data.
         */
         template <class T, size_t rank>
         xt::xtensor<T, rank> mapToRegular(const xt::xtensor<T, rank>& arg) const;
 
     private:
-        // the "FineLayer" mesh to map
-        GooseFEM::Mesh::Quad4::FineLayer m_finelayer;
-
-        // the new "Regular" mesh to which to map
-        GooseFEM::Mesh::Quad4::Regular m_regular;
 
-        // mapping
-        std::vector<std::vector<size_t>> m_elem_regular;
-        std::vector<std::vector<double>> m_frac_regular;
+        GooseFEM::Mesh::Quad4::FineLayer m_finelayer; ///< the FineLayer mesh to map
+        GooseFEM::Mesh::Quad4::Regular m_regular; ///< the new Regular mesh to which to map
+        std::vector<std::vector<size_t>> m_elem_regular; ///< see getMap()
+        std::vector<std::vector<double>> m_frac_regular; ///< see getMapFraction()
     };
 
 } // namespace Map
 
 } // namespace Quad4
 } // namespace Mesh
 } // namespace GooseFEM
 
 #include "MeshQuad4.hpp"
 
 #endif
diff --git a/include/GooseFEM/MeshQuad4.hpp b/include/GooseFEM/MeshQuad4.hpp
index 2826760..01a5262 100644
--- a/include/GooseFEM/MeshQuad4.hpp
+++ b/include/GooseFEM/MeshQuad4.hpp
@@ -1,1653 +1,1653 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_MESHQUAD4_HPP
 #define GOOSEFEM_MESHQUAD4_HPP
 
 #include "MeshQuad4.h"
 
 namespace GooseFEM {
 namespace Mesh {
 namespace Quad4 {
 
 inline Regular::Regular(size_t nelx, size_t nely, double h) : m_h(h), m_nelx(nelx), m_nely(nely)
 {
     GOOSEFEM_ASSERT(m_nelx >= 1ul);
     GOOSEFEM_ASSERT(m_nely >= 1ul);
 
     m_nnode = (m_nelx + 1) * (m_nely + 1);
     m_nelem = m_nelx * m_nely;
 }
 
 inline size_t Regular::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t Regular::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t Regular::nne() const
 {
     return m_nne;
 }
 
 inline size_t Regular::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t Regular::nelx() const
 {
     return m_nelx;
 }
 
 inline size_t Regular::nely() const
 {
     return m_nely;
 }
 
 inline double Regular::h() const
 {
     return m_h;
 }
 
 inline ElementType Regular::getElementType() const
 {
     return ElementType::Quad4;
 }
 
 inline xt::xtensor<double, 2> Regular::coor() const
 {
     xt::xtensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});
 
     xt::xtensor<double, 1> x =
         xt::linspace<double>(0.0, m_h * static_cast<double>(m_nelx), m_nelx + 1);
     xt::xtensor<double, 1> y =
         xt::linspace<double>(0.0, m_h * static_cast<double>(m_nely), m_nely + 1);
 
     size_t inode = 0;
 
     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);
             ++inode;
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 2> Regular::conn() const
 {
     xt::xtensor<size_t, 2> ret = xt::empty<size_t>({m_nelem, m_nne});
 
     size_t ielem = 0;
 
     for (size_t iy = 0; iy < m_nely; ++iy) {
         for (size_t ix = 0; ix < m_nelx; ++ix) {
             ret(ielem, 0) = (iy) * (m_nelx + 1) + (ix);
             ret(ielem, 1) = (iy) * (m_nelx + 1) + (ix + 1);
             ret(ielem, 3) = (iy + 1) * (m_nelx + 1) + (ix);
             ret(ielem, 2) = (iy + 1) * (m_nelx + 1) + (ix + 1);
             ++ielem;
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBottomEdge() const
 {
     return xt::arange<size_t>(m_nelx + 1);
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesTopEdge() const
 {
     return xt::arange<size_t>(m_nelx + 1) + m_nely * (m_nelx + 1);
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesLeftEdge() const
 {
     return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1);
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesRightEdge() const
 {
     return xt::arange<size_t>(m_nely + 1) * (m_nelx + 1) + m_nelx;
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesBottomOpenEdge() const
 {
     return xt::arange<size_t>(1, m_nelx);
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesTopOpenEdge() const
 {
     return xt::arange<size_t>(1, m_nelx) + m_nely * (m_nelx + 1);
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesLeftOpenEdge() const
 {
     return xt::arange<size_t>(1, m_nely) * (m_nelx + 1);
 }
 
 inline xt::xtensor<size_t, 1> Regular::nodesRightOpenEdge() const
 {
     return xt::arange<size_t>(1, m_nely) * (m_nelx + 1) + m_nelx;
 }
 
 inline size_t Regular::nodesBottomLeftCorner() const
 {
     return 0;
 }
 
 inline size_t Regular::nodesBottomRightCorner() const
 {
     return m_nelx;
 }
 
 inline size_t Regular::nodesTopLeftCorner() const
 {
     return m_nely * (m_nelx + 1);
 }
 
 inline size_t Regular::nodesTopRightCorner() const
 {
     return m_nely * (m_nelx + 1) + m_nelx;
 }
 
 inline size_t Regular::nodesLeftBottomCorner() const
 {
     return nodesBottomLeftCorner();
 }
 
 inline size_t Regular::nodesLeftTopCorner() const
 {
     return nodesTopLeftCorner();
 }
 
 inline size_t Regular::nodesRightBottomCorner() const
 {
     return nodesBottomRightCorner();
 }
 
 inline size_t Regular::nodesRightTopCorner() const
 {
     return nodesTopRightCorner();
 }
 
 inline xt::xtensor<size_t, 2> Regular::nodesPeriodic() const
 {
     xt::xtensor<size_t, 1> bot = nodesBottomOpenEdge();
     xt::xtensor<size_t, 1> top = nodesTopOpenEdge();
     xt::xtensor<size_t, 1> lft = nodesLeftOpenEdge();
     xt::xtensor<size_t, 1> rgt = nodesRightOpenEdge();
     std::array<size_t, 2> shape = {bot.size() + lft.size() + 3ul, 2ul};
     xt::xtensor<size_t, 2> ret = xt::empty<size_t>(shape);
 
     ret(0, 0) = nodesBottomLeftCorner();
     ret(0, 1) = nodesBottomRightCorner();
 
     ret(1, 0) = nodesBottomLeftCorner();
     ret(1, 1) = nodesTopRightCorner();
 
     ret(2, 0) = nodesBottomLeftCorner();
     ret(2, 1) = nodesTopLeftCorner();
 
     size_t i = 3;
 
     xt::view(ret, xt::range(i, i + bot.size()), 0) = bot;
     xt::view(ret, xt::range(i, i + bot.size()), 1) = top;
 
     i += bot.size();
 
     xt::view(ret, xt::range(i, i + lft.size()), 0) = lft;
     xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt;
 
     return ret;
 }
 
 inline size_t Regular::nodesOrigin() const
 {
     return nodesBottomLeftCorner();
 }
 
 inline xt::xtensor<size_t, 2> Regular::dofs() const
 {
     return GooseFEM::Mesh::dofs(m_nnode, m_ndim);
 }
 
 inline xt::xtensor<size_t, 2> Regular::dofsPeriodic() const
 {
     xt::xtensor<size_t, 2> ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim);
     xt::xtensor<size_t, 2> nodePer = nodesPeriodic();
     xt::xtensor<size_t, 1> independent = xt::view(nodePer, xt::all(), 0);
     xt::xtensor<size_t, 1> dependent = xt::view(nodePer, xt::all(), 1);
 
     for (size_t j = 0; j < m_ndim; ++j) {
         xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j);
     }
 
     return GooseFEM::Mesh::renumber(ret);
 }
 
 inline xt::xtensor<size_t, 2> Regular::elementgrid() const
 {
     return xt::arange<size_t>(m_nelem).reshape({m_nely, m_nelx});
 }
 
 inline FineLayer::FineLayer(size_t nelx, size_t nely, double h, size_t nfine)
 {
     this->init(nelx, nely, h, nfine);
 }
 
 inline FineLayer::FineLayer(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn)
 {
     this->map(coor, conn);
 }
 
 inline void FineLayer::init(size_t nelx, size_t nely, double h, size_t nfine)
 {
     GOOSEFEM_ASSERT(nelx >= 1ul);
     GOOSEFEM_ASSERT(nely >= 1ul);
 
     m_h = h;
     m_Lx = m_h * static_cast<double>(nelx);
 
     // 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<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-direction should fit the total number of elements in x-direction
     // (3) a certain number of layers have the minimum size "1" (are fine)
     for (size_t iy = nfine;;) {
         // initialize current size in y-direction
         if (iy == nfine) {
             ntot = nfine;
         }
         // check to stop
         if (iy >= nely || ntot >= nmin) {
             nely = iy;
             break;
         }
 
         // rules (1,2) satisfied: coarsen in x-direction
         if (3 * nhy(iy) <= ntot && nelx % (3 * nhx(iy)) == 0 && ntot + nhy(iy) < nmin) {
             refine(iy) = 0;
             nhy(iy) *= 2;
             auto vnhy = xt::view(nhy, xt::range(iy + 1, _));
             auto vnhx = xt::view(nhx, xt::range(iy, _));
             vnhy *= 3;
             vnhx *= 3;
         }
 
         // update the number of elements in y-direction
         ntot += nhy(iy);
         // proceed to next element layer in y-direction
         ++iy;
         // check to stop
         if (iy >= nely || 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_refine = xt::empty<int>({nely * 2 - 1});
     m_nelx = 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_refine(iy) = refine(nely - iy - 1);
     }
     // - upper half
     for (size_t iy = 0; iy < nely - 1; ++iy) {
         m_nhx(iy + nely) = nhx(iy + 1);
         m_nhy(iy + nely) = nhy(iy + 1);
         m_refine(iy + nely) = refine(iy + 1);
     }
 
     // update size
     nely = m_nhx.size();
 
     // compute the number of elements per element layer in y-direction
     for (size_t iy = 0; iy < nely; ++iy) {
         m_nelx(iy) = nelx / m_nhx(iy);
     }
 
     // compute the number of nodes per node layer in y-direction
     for (size_t iy = 0; iy < (nely + 1) / 2; ++iy) {
         m_nnd(iy) = m_nelx(iy) + 1;
     }
     for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
         m_nnd(iy + 1) = m_nelx(iy) + 1;
     }
 
     // compute mesh dimensions
 
     // initialize
     m_nnode = 0;
     m_nelem = 0;
     m_startNode(0) = 0;
 
     // loop over element layers (bottom -> middle, elements become finer)
     for (size_t i = 0; i < (nely - 1) / 2; ++i) {
         // - store the first element of the layer
         m_startElem(i) = m_nelem;
         // - add the nodes of this layer
         if (m_refine(i) == 0) {
             m_nnode += (3 * m_nelx(i) + 1);
         }
         else {
             m_nnode += (m_nelx(i) + 1);
         }
         // - add the elements of this layer
         if (m_refine(i) == 0) {
             m_nelem += (4 * m_nelx(i));
         }
         else {
             m_nelem += (m_nelx(i));
         }
         // - store the starting node of the next layer
         m_startNode(i + 1) = m_nnode;
     }
 
     // loop over element layers (middle -> top, elements become coarser)
     for (size_t i = (nely - 1) / 2; i < nely; ++i) {
         // - store the first element of the layer
         m_startElem(i) = m_nelem;
         // - add the nodes of this layer
         if (m_refine(i) == 0) {
             m_nnode += (5 * m_nelx(i) + 1);
         }
         else {
             m_nnode += (m_nelx(i) + 1);
         }
         // - add the elements of this layer
         if (m_refine(i) == 0) {
             m_nelem += (4 * m_nelx(i));
         }
         else {
             m_nelem += (m_nelx(i));
         }
         // - store the starting node of the next layer
         m_startNode(i + 1) = m_nnode;
     }
     // - add the top row of nodes
     m_nnode += m_nelx(nely - 1) + 1;
 }
 
 inline size_t FineLayer::nelem() const
 {
     return m_nelem;
 }
 
 inline size_t FineLayer::nnode() const
 {
     return m_nnode;
 }
 
 inline size_t FineLayer::nne() const
 {
     return m_nne;
 }
 
 inline size_t FineLayer::ndim() const
 {
     return m_ndim;
 }
 
 inline size_t FineLayer::nelx() const
 {
     return xt::amax(m_nelx)();
 }
 
 inline size_t FineLayer::nely() const
 {
     return xt::sum(m_nhy)();
 }
 
 inline double FineLayer::h() const
 {
     return m_h;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::elemrow_nhx() const
 {
     return m_nhx;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::elemrow_nhy() const
 {
     return m_nhy;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::elemrow_nelem() const
 {
     return m_nelx;
 }
 
 inline ElementType FineLayer::getElementType() const
 {
     return ElementType::Quad4;
 }
 
 inline xt::xtensor<double, 2> FineLayer::coor() const
 {
     // allocate output
     xt::xtensor<double, 2> ret = xt::empty<double>({m_nnode, m_ndim});
 
     // current node, number of element layers
     size_t inode = 0;
     size_t nely = static_cast<size_t>(m_nhy.size());
 
     // y-position of each main node layer (i.e. excluding node layers for refinement/coarsening)
     // - allocate
     xt::xtensor<double, 1> y = xt::empty<double>({nely + 1});
     // - initialize
     y(0) = 0.0;
     // - compute
     for (size_t iy = 1; iy < nely + 1; ++iy) {
         y(iy) = y(iy - 1) + m_nhy(iy - 1) * m_h;
     }
 
     // loop over element layers (bottom -> middle) : add bottom layer (+ refinement layer) of nodes
 
     for (size_t iy = 0;; ++iy) {
         // get positions along the x- and z-axis
         xt::xtensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_nelx(iy) + 1);
 
         // add nodes of the bottom layer of this element
         for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
             ret(inode, 0) = x(ix);
             ret(inode, 1) = y(iy);
             ++inode;
         }
 
         // stop at middle layer
         if (iy == (nely - 1) / 2) {
             break;
         }
 
         // add extra nodes of the intermediate layer, for refinement in x-direction
         if (m_refine(iy) == 0) {
             // - get position offset in x- and y-direction
             double dx = m_h * static_cast<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 ix = 0; ix < m_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;
                     ++inode;
                 }
             }
         }
     }
 
     // loop over element layers (middle -> top) : add (refinement layer +) top layer of nodes
 
     for (size_t iy = (nely - 1) / 2; iy < nely; ++iy) {
         // get positions along the x- and z-axis
         xt::xtensor<double, 1> x = xt::linspace<double>(0.0, m_Lx, m_nelx(iy) + 1);
 
         // 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 ix = 0; ix < m_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;
                     ++inode;
                 }
             }
         }
 
         // add nodes of the top layer of this element
         for (size_t ix = 0; ix < m_nelx(iy) + 1; ++ix) {
             ret(inode, 0) = x(ix);
             ret(inode, 1) = y(iy + 1);
             ++inode;
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 2> FineLayer::conn() const
 {
     // allocate output
     xt::xtensor<size_t, 2> ret = xt::empty<size_t>({m_nelem, m_nne});
 
     // current element, number of element layers, starting nodes of each node layer
     size_t ielem = 0;
     size_t nely = static_cast<size_t>(m_nhy.size());
     size_t bot, mid, top;
 
     // loop over all element layers
     for (size_t iy = 0; iy < nely; ++iy) {
         // - get: starting nodes of bottom(, middle) and top layer
         bot = m_startNode(iy);
         mid = m_startNode(iy) + m_nnd(iy);
         top = m_startNode(iy + 1);
 
         // - define connectivity: no coarsening/refinement
         if (m_refine(iy) == -1) {
             for (size_t ix = 0; ix < m_nelx(iy); ++ix) {
                 ret(ielem, 0) = bot + (ix);
                 ret(ielem, 1) = bot + (ix + 1);
                 ret(ielem, 2) = top + (ix + 1);
                 ret(ielem, 3) = top + (ix);
                 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 ix = 0; ix < m_nelx(iy); ++ix) {
                 // -- bottom element
                 ret(ielem, 0) = bot + (ix);
                 ret(ielem, 1) = bot + (ix + 1);
                 ret(ielem, 2) = mid + (2 * ix + 1);
                 ret(ielem, 3) = mid + (2 * ix);
                 ielem++;
                 // -- top-right element
                 ret(ielem, 0) = bot + (ix + 1);
                 ret(ielem, 1) = top + (3 * ix + 3);
                 ret(ielem, 2) = top + (3 * ix + 2);
                 ret(ielem, 3) = mid + (2 * ix + 1);
                 ielem++;
                 // -- top-center element
                 ret(ielem, 0) = mid + (2 * ix);
                 ret(ielem, 1) = mid + (2 * ix + 1);
                 ret(ielem, 2) = top + (3 * ix + 2);
                 ret(ielem, 3) = top + (3 * ix + 1);
                 ielem++;
                 // -- top-left element
                 ret(ielem, 0) = bot + (ix);
                 ret(ielem, 1) = mid + (2 * ix);
                 ret(ielem, 2) = top + (3 * ix + 1);
                 ret(ielem, 3) = top + (3 * ix);
                 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 ix = 0; ix < m_nelx(iy); ++ix) {
                 // -- lower-left element
                 ret(ielem, 0) = bot + (3 * ix);
                 ret(ielem, 1) = bot + (3 * ix + 1);
                 ret(ielem, 2) = mid + (2 * ix);
                 ret(ielem, 3) = top + (ix);
                 ielem++;
                 // -- lower-center element
                 ret(ielem, 0) = bot + (3 * ix + 1);
                 ret(ielem, 1) = bot + (3 * ix + 2);
                 ret(ielem, 2) = mid + (2 * ix + 1);
                 ret(ielem, 3) = mid + (2 * ix);
                 ielem++;
                 // -- lower-right element
                 ret(ielem, 0) = bot + (3 * ix + 2);
                 ret(ielem, 1) = bot + (3 * ix + 3);
                 ret(ielem, 2) = top + (ix + 1);
                 ret(ielem, 3) = mid + (2 * ix + 1);
                 ielem++;
                 // -- upper element
                 ret(ielem, 0) = mid + (2 * ix);
                 ret(ielem, 1) = mid + (2 * ix + 1);
                 ret(ielem, 2) = top + (ix + 1);
                 ret(ielem, 3) = top + (ix);
                 ielem++;
             }
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::elementsMiddleLayer() const
 {
     size_t nely = m_nhy.size();
     size_t iy = (nely - 1) / 2;
     return m_startElem(iy) + xt::arange<size_t>(m_nelx(iy));
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::elementsLayer(size_t iy) const
 {
     GOOSEFEM_ASSERT(iy < m_nelx.size());
     size_t n = m_nelx(iy);
     if (m_refine(iy) != -1) {
         n *= 4;
     }
     return m_startElem(iy) + xt::arange<size_t>(n);
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::elementgrid_ravel(
     std::vector<size_t> start_stop_rows,
     std::vector<size_t> start_stop_cols) const
 {
     GOOSEFEM_ASSERT(start_stop_rows.size() == 0 || start_stop_rows.size() == 2);
     GOOSEFEM_ASSERT(start_stop_cols.size() == 0 || start_stop_cols.size() == 2);
 
     std::array<size_t, 2> rows;
     std::array<size_t, 2> cols;
 
     if (start_stop_rows.size() == 2) {
         std::copy(start_stop_rows.begin(), start_stop_rows.end(), rows.begin());
         GOOSEFEM_ASSERT(rows[0] <= this->nely());
         GOOSEFEM_ASSERT(rows[1] <= this->nely());
     }
     else {
         rows[0] = 0;
         rows[1] = this->nely();
     }
 
     if (start_stop_cols.size() == 2) {
         std::copy(start_stop_cols.begin(), start_stop_cols.end(), cols.begin());
         GOOSEFEM_ASSERT(cols[0] <= this->nelx());
         GOOSEFEM_ASSERT(cols[1] <= this->nelx());
     }
     else {
         cols[0] = 0;
         cols[1] = this->nelx();
     }
 
     if (rows[0] == rows[1] || cols[0] == cols[1]) {
         xt::xtensor<size_t, 1> ret = xt::empty<size_t>({0});
         return ret;
     }
 
     // Compute dimensions
 
     auto H = xt::cumsum(m_nhy);
     size_t yl = 0;
     if (rows[0] > 0) {
         yl = xt::argmax(H > rows[0])();
     }
     size_t yu = xt::argmax(H >= rows[1])();
     size_t hx = std::max(m_nhx(yl), m_nhx(yu));
     size_t xl = (cols[0] - cols[0] % hx) / hx;
     size_t xu = (cols[1] - cols[1] % hx) / hx;
 
     // Allocate output
 
     size_t N = 0;
 
     for (size_t i = yl; i <= yu; ++i) {
         // no refinement
         size_t n = (xu - xl) * hx / m_nhx(i);
         // refinement
         if (m_refine(i) != -1) {
             n *= 4;
         }
         N += n;
     }
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({N});
 
     // Write output
 
     N = 0;
 
     for (size_t i = yl; i <= yu; ++i) {
         // no refinement
         size_t n = (xu - xl) * hx / m_nhx(i);
         size_t h = hx;
         // refinement
         if (m_refine(i) != -1) {
             n *= 4;
             h *= 4;
         }
         xt::xtensor<size_t, 1> e = m_startElem(i) + xl * h / m_nhx(i) + xt::arange<size_t>(n);
         xt::view(ret, xt::range(N, N + n)) = e;
         N += n;
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::elementgrid_around_ravel(
     size_t e,
     size_t size,
     bool periodic)
 {
     GOOSEFEM_WIP_ASSERT(periodic == true);
 
     size_t iy = xt::argmin(m_startElem <= e)() - 1;
     size_t nel = m_nelx(iy);
 
     GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2);
 
     auto H = xt::cumsum(m_nhy);
 
     if (2 * size >= H(H.size() - 1)) {
         return xt::arange<size_t>(this->nelem());
     }
 
     size_t hy = H(iy);
     size_t l = xt::argmax(H > (hy - size - 1))();
     size_t u = xt::argmax(H >= (hy + size))();
     size_t lh = 0;
     if (l > 0) {
         lh = H(l - 1);
     }
     size_t uh = H(u);
 
     size_t step = xt::amax(m_nhx)();
     size_t relx = (e - m_startElem(iy)) % step;
     size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx;
     size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step;
     size_t dx = m_nhx(u);
     size_t xl = 0;
     size_t xu = nel;
     if (mid >= size) {
         xl = mid - size;
     }
     if (mid + size < nel) {
         xu = mid + size + 1;
     }
     xl = xl - xl % dx;
     xu = xu - xu % dx;
     if (mid - xl < size) {
         if (xl < dx) {
             xl = 0;
         }
         else {
             xl -= dx;
         }
     }
     if (xu - mid < size) {
         if (xu > nel - dx) {
             xu = nel;
         }
         else {
             xu += dx;
         }
     }
 
     auto ret = this->elementgrid_ravel({lh, uh}, {xl, xu});
     auto map = this->roll(nroll);
     return xt::view(map, xt::keep(ret));
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::elementgrid_leftright(
     size_t e,
     size_t left,
     size_t right,
     bool periodic)
 {
     GOOSEFEM_WIP_ASSERT(periodic == true);
 
     size_t iy = xt::argmin(m_startElem <= e)() - 1;
     size_t nel = m_nelx(iy);
 
     GOOSEFEM_WIP_ASSERT(iy == (m_nhy.size() - 1) / 2);
 
     size_t step = xt::amax(m_nhx)();
     size_t relx = (e - m_startElem(iy)) % step;
     size_t mid = (nel / step - (nel / step) % 2) / 2 * step + relx;
     size_t nroll = (nel - (nel - mid + e - m_startElem(iy)) % nel) / step;
     size_t dx = m_nhx(iy);
     size_t xl = 0;
     size_t xu = nel;
     if (mid >= left) {
         xl = mid - left;
     }
     if (mid + right < nel) {
         xu = mid + right + 1;
     }
     xl = xl - xl % dx;
     xu = xu - xu % dx;
     if (mid - xl < left) {
         if (xl < dx) {
             xl = 0;
         }
         else {
             xl -= dx;
         }
     }
     if (xu - mid < right) {
         if (xu > nel - dx) {
             xu = nel;
         }
         else {
             xu += dx;
         }
     }
 
     auto H = xt::cumsum(m_nhy);
     size_t yl = 0;
     if (iy > 0) {
         yl = H(iy - 1);
     }
     auto ret = this->elementgrid_ravel({yl, H(iy)}, {xl, xu});
     auto map = this->roll(nroll);
     return xt::view(map, xt::keep(ret));
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBottomEdge() const
 {
     return m_startNode(0) + xt::arange<size_t>(m_nelx(0) + 1);
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesTopEdge() const
 {
     size_t nely = m_nhy.size();
     return m_startNode(nely) + xt::arange<size_t>(m_nelx(nely - 1) + 1);
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesLeftEdge() const
 {
     size_t nely = m_nhy.size();
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
 
     size_t i = 0;
     size_t j = (nely + 1) / 2;
     size_t k = (nely - 1) / 2;
     size_t l = nely;
 
     xt::view(ret, xt::range(i, j)) = xt::view(m_startNode, xt::range(i, j));
     xt::view(ret, xt::range(k + 1, l + 1)) = xt::view(m_startNode, xt::range(k + 1, l + 1));
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesRightEdge() const
 {
     size_t nely = m_nhy.size();
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely + 1});
 
     size_t i = 0;
     size_t j = (nely + 1) / 2;
     size_t k = (nely - 1) / 2;
     size_t l = nely;
 
     xt::view(ret, xt::range(i, j)) =
         xt::view(m_startNode, xt::range(i, j)) + xt::view(m_nelx, xt::range(i, j));
 
     xt::view(ret, xt::range(k + 1, l + 1)) =
         xt::view(m_startNode, xt::range(k + 1, l + 1)) + xt::view(m_nelx, xt::range(k, l));
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesBottomOpenEdge() const
 {
     return m_startNode(0) + xt::arange<size_t>(1, m_nelx(0));
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesTopOpenEdge() const
 {
     size_t nely = m_nhy.size();
 
     return m_startNode(nely) + xt::arange<size_t>(1, m_nelx(nely - 1));
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesLeftOpenEdge() const
 {
     size_t nely = m_nhy.size();
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
 
     size_t i = 0;
     size_t j = (nely + 1) / 2;
     size_t k = (nely - 1) / 2;
     size_t l = nely;
 
     xt::view(ret, xt::range(i, j - 1)) = xt::view(m_startNode, xt::range(i + 1, j));
     xt::view(ret, xt::range(k, l - 1)) = xt::view(m_startNode, xt::range(k + 1, l));
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::nodesRightOpenEdge() const
 {
     size_t nely = m_nhy.size();
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({nely - 1});
 
     size_t i = 0;
     size_t j = (nely + 1) / 2;
     size_t k = (nely - 1) / 2;
     size_t l = nely;
 
     xt::view(ret, xt::range(i, j - 1)) =
         xt::view(m_startNode, xt::range(i + 1, j)) + xt::view(m_nelx, xt::range(i + 1, j));
 
     xt::view(ret, xt::range(k, l - 1)) =
         xt::view(m_startNode, xt::range(k + 1, l)) + xt::view(m_nelx, xt::range(k, l - 1));
 
     return ret;
 }
 
 inline size_t FineLayer::nodesBottomLeftCorner() const
 {
     return m_startNode(0);
 }
 
 inline size_t FineLayer::nodesBottomRightCorner() const
 {
     return m_startNode(0) + m_nelx(0);
 }
 
 inline size_t FineLayer::nodesTopLeftCorner() const
 {
     size_t nely = m_nhy.size();
 
     return m_startNode(nely);
 }
 
 inline size_t FineLayer::nodesTopRightCorner() const
 {
     size_t nely = m_nhy.size();
 
     return m_startNode(nely) + m_nelx(nely - 1);
 }
 
 inline size_t FineLayer::nodesLeftBottomCorner() const
 {
     return nodesBottomLeftCorner();
 }
 
 inline size_t FineLayer::nodesRightBottomCorner() const
 {
     return nodesBottomRightCorner();
 }
 
 inline size_t FineLayer::nodesLeftTopCorner() const
 {
     return nodesTopLeftCorner();
 }
 
 inline size_t FineLayer::nodesRightTopCorner() const
 {
     return nodesTopRightCorner();
 }
 
 inline xt::xtensor<size_t, 2> FineLayer::nodesPeriodic() const
 {
     xt::xtensor<size_t, 1> bot = nodesBottomOpenEdge();
     xt::xtensor<size_t, 1> top = nodesTopOpenEdge();
     xt::xtensor<size_t, 1> lft = nodesLeftOpenEdge();
     xt::xtensor<size_t, 1> rgt = nodesRightOpenEdge();
     std::array<size_t, 2> shape = {bot.size() + lft.size() + 3ul, 2ul};
     xt::xtensor<size_t, 2> ret = xt::empty<size_t>(shape);
 
     ret(0, 0) = nodesBottomLeftCorner();
     ret(0, 1) = nodesBottomRightCorner();
 
     ret(1, 0) = nodesBottomLeftCorner();
     ret(1, 1) = nodesTopRightCorner();
 
     ret(2, 0) = nodesBottomLeftCorner();
     ret(2, 1) = nodesTopLeftCorner();
 
     size_t i = 3;
 
     xt::view(ret, xt::range(i, i + bot.size()), 0) = bot;
     xt::view(ret, xt::range(i, i + bot.size()), 1) = top;
 
     i += bot.size();
 
     xt::view(ret, xt::range(i, i + lft.size()), 0) = lft;
     xt::view(ret, xt::range(i, i + lft.size()), 1) = rgt;
 
     return ret;
 }
 
 inline size_t FineLayer::nodesOrigin() const
 {
     return nodesBottomLeftCorner();
 }
 
 inline xt::xtensor<size_t, 2> FineLayer::dofs() const
 {
     return GooseFEM::Mesh::dofs(m_nnode, m_ndim);
 }
 
 inline xt::xtensor<size_t, 2> FineLayer::dofsPeriodic() const
 {
     xt::xtensor<size_t, 2> ret = GooseFEM::Mesh::dofs(m_nnode, m_ndim);
     xt::xtensor<size_t, 2> nodePer = nodesPeriodic();
     xt::xtensor<size_t, 1> independent = xt::view(nodePer, xt::all(), 0);
     xt::xtensor<size_t, 1> dependent = xt::view(nodePer, xt::all(), 1);
 
     for (size_t j = 0; j < m_ndim; ++j) {
         xt::view(ret, xt::keep(dependent), j) = xt::view(ret, xt::keep(independent), j);
     }
 
     return GooseFEM::Mesh::renumber(ret);
 }
 
 inline xt::xtensor<size_t, 1> FineLayer::roll(size_t n)
 {
     auto conn = this->conn();
     size_t nely = static_cast<size_t>(m_nhy.size());
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({m_nelem});
 
     // loop over all element layers
     for (size_t iy = 0; iy < nely; ++iy) {
 
         // no refinement
         size_t shift = n * (m_nelx(iy) / m_nelx(0));
         size_t nel = m_nelx(iy);
 
         // refinement
         if (m_refine(iy) != -1) {
             shift = n * (m_nelx(iy) / m_nelx(0)) * 4;
             nel = m_nelx(iy) * 4;
         }
 
         // element numbers of the layer, and roll them
         auto e = m_startElem(iy) + xt::arange<size_t>(nel);
         xt::view(ret, xt::range(m_startElem(iy), m_startElem(iy) + nel)) = xt::roll(e, shift);
     }
 
     return ret;
 }
 
 inline void FineLayer::map(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn)
 {
     GOOSEFEM_ASSERT(coor.shape(1) == 2);
     GOOSEFEM_ASSERT(conn.shape(1) == 4);
     GOOSEFEM_ASSERT(conn.shape(0) > 0);
     GOOSEFEM_ASSERT(coor.shape(0) >= 4);
     GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
 
     if (conn.shape(0) == 1) {
         this->init(1, 1, coor(conn(0, 1), 0) - coor(conn(0, 0), 0));
         GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn)));
         GOOSEFEM_CHECK(xt::allclose(this->coor(), coor));
         return;
     }
 
     // Identify the middle layer
 
     size_t emid = (conn.shape(0) - conn.shape(0) % 2) / 2;
     size_t eleft = emid;
     size_t eright = emid;
 
     while (conn(eleft, 0) == conn(eleft - 1, 1) && eleft > 0) {
         eleft--;
     }
 
     while (conn(eright, 1) == conn(eright + 1, 0) && eright < conn.shape(0) - 1) {
         eright++;
     }
 
     GOOSEFEM_CHECK(xt::allclose(coor(conn(eleft, 0), 0), 0.0));
 
     // Get element sizes along the middle layer
 
     auto n0 = xt::view(conn, xt::range(eleft, eright + 1), 0);
     auto n1 = xt::view(conn, xt::range(eleft, eright + 1), 1);
     auto n2 = xt::view(conn, xt::range(eleft, eright + 1), 2);
     auto dx = xt::view(coor, xt::keep(n1), 0) - xt::view(coor, xt::keep(n0), 0);
     auto dy = xt::view(coor, xt::keep(n2), 1) - xt::view(coor, xt::keep(n1), 1);
     auto hx = xt::amin(dx)();
     auto hy = xt::amin(dy)();
 
     GOOSEFEM_CHECK(xt::allclose(hx, hy));
     GOOSEFEM_CHECK(xt::allclose(dx, hx));
     GOOSEFEM_CHECK(xt::allclose(dy, hy));
 
     // Extract shape and initialise
 
     size_t nelx = eright - eleft + 1;
     size_t nely = static_cast<size_t>((coor(coor.shape(0) - 1, 1) - coor(0, 1)) / hx);
     this->init(nelx, nely, hx);
     GOOSEFEM_CHECK(xt::all(xt::equal(this->conn(), conn)));
     GOOSEFEM_CHECK(xt::allclose(this->coor(), coor));
     GOOSEFEM_CHECK(xt::all(xt::equal(this->elementsMiddleLayer(), eleft + xt::arange<size_t>(nelx))));
 }
 
 namespace Map {
 
 inline RefineRegular::RefineRegular(
     const GooseFEM::Mesh::Quad4::Regular& mesh, size_t nx, size_t ny)
     : m_coarse(mesh), m_nx(nx), m_ny(ny)
 {
     m_fine = Regular(nx * m_coarse.nelx(), ny * m_coarse.nely(), m_coarse.h());
 
     xt::xtensor<size_t, 2> elmat_coarse = m_coarse.elementgrid();
     xt::xtensor<size_t, 2> elmat_fine = m_fine.elementgrid();
 
     m_coarse2fine = xt::empty<size_t>({m_coarse.nelem(), nx * ny});
 
     for (size_t i = 0; i < elmat_coarse.shape(0); ++i) {
         for (size_t j = 0; j < elmat_coarse.shape(1); ++j) {
             xt::view(m_coarse2fine, elmat_coarse(i, j), xt::all()) = xt::flatten(xt::view(
                 elmat_fine, xt::range(i * ny, (i + 1) * ny), xt::range(j * nx, (j + 1) * nx)));
         }
     }
 }
 
 inline size_t RefineRegular::nx() const
 {
     return m_nx;
 }
 
 inline size_t RefineRegular::ny() const
 {
     return m_ny;
 }
 
 inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getCoarseMesh() const
 {
     return m_coarse;
 }
 
 inline GooseFEM::Mesh::Quad4::Regular RefineRegular::getFineMesh() const
 {
     return m_fine;
 }
 
 inline xt::xtensor<size_t, 2> RefineRegular::getMap() const
 {
     return m_coarse2fine;
 }
 
 template <class T, size_t rank>
-xt::xtensor<T, rank> RefineRegular::meanToCoarse(const xt::xtensor<T, rank>& data) const
+inline xt::xtensor<T, rank> RefineRegular::meanToCoarse(const xt::xtensor<T, rank>& data) const
 {
     GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
 
     std::array<size_t, rank> shape;
     std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]);
     shape[0] = m_coarse2fine.shape(0);
 
     xt::xtensor<T, rank> ret = xt::empty<T>(shape);
 
     for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
         auto e = xt::view(m_coarse2fine, i, xt::all());
         auto d = xt::view(data, xt::keep(e));
         xt::view(ret, i) = xt::mean(d, 0);
     }
 
     return ret;
 }
 
 template <class T, size_t rank, class S>
-xt::xtensor<T, rank> RefineRegular::averageToCoarse(
+inline xt::xtensor<T, rank> RefineRegular::averageToCoarse(
     const xt::xtensor<T, rank>& data,
     const xt::xtensor<S, rank>& weights) const
 {
     GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
 
     std::array<size_t, rank> shape;
     std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]);
     shape[0] = m_coarse2fine.shape(0);
 
     xt::xtensor<T, rank> ret = xt::empty<T>(shape);
 
     for (size_t i = 0; i < m_coarse2fine.shape(0); ++i) {
         auto e = xt::view(m_coarse2fine, i, xt::all());
         xt::xtensor<T, rank> d = xt::view(data, xt::keep(e));
         xt::xtensor<T, rank> w = xt::view(weights, xt::keep(e));
         xt::view(ret, i) = xt::average(d, w, {0});
     }
 
     return ret;
 }
 
 inline xt::xtensor<double, 2> RefineRegular::mapToCoarse(const xt::xtensor<double, 1>& data) const
 {
     std::cout << "mapToCoarse is deprecated, use meanToCoarse or averageToCoarse" << std::endl;
 
     GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
 
     size_t m = m_coarse2fine.shape(0);
     size_t n = m_coarse2fine.shape(1);
 
     xt::xtensor<double, 2> ret = xt::empty<double>({m, n});
 
     for (size_t i = 0; i < m; ++i) {
         auto e = xt::view(m_coarse2fine, i, xt::all());
         xt::view(ret, i) = xt::view(data, xt::keep(e));
     }
 
     return ret;
 }
 
 inline xt::xtensor<double, 2> RefineRegular::mapToCoarse(const xt::xtensor<double, 2>& data) const
 {
     std::cout << "mapToCoarse is deprecated, use meanToCoarse or averageToCoarse" << std::endl;
 
     GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
 
     size_t m = m_coarse2fine.shape(0);
     size_t n = m_coarse2fine.shape(1);
     size_t N = data.shape(1);
 
     xt::xtensor<double, 2> ret = xt::empty<double>({m, n * N});
 
     for (size_t i = 0; i < m; ++i) {
         auto e = xt::view(m_coarse2fine, i, xt::all());
         for (size_t q = 0; q < data.shape(1); ++q) {
             xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q);
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<double, 4> RefineRegular::mapToCoarse(const xt::xtensor<double, 4>& data) const
 {
     std::cout << "mapToCoarse is deprecated, use meanToCoarse or averageToCoarse" << std::endl;
 
     GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.size());
 
     size_t m = m_coarse2fine.shape(0);
     size_t n = m_coarse2fine.shape(1);
     size_t N = data.shape(1);
 
     xt::xtensor<double, 4> ret = xt::empty<double>({m, n * N, data.shape(2), data.shape(3)});
 
     for (size_t i = 0; i < m; ++i) {
         auto e = xt::view(m_coarse2fine, i, xt::all());
         for (size_t q = 0; q < data.shape(1); ++q) {
             xt::view(ret, i, xt::range(q + 0, q + n * N, N)) = xt::view(data, xt::keep(e), q);
         }
     }
 
     return ret;
 }
 
 template <class T, size_t rank>
 inline xt::xtensor<T, rank> RefineRegular::mapToFine(const xt::xtensor<T, rank>& data) const
 {
     GOOSEFEM_ASSERT(data.shape(0) == m_coarse2fine.shape(0));
 
     std::array<size_t, rank> shape;
     std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]);
     shape[0] = m_coarse2fine.size();
 
     xt::xtensor<T, rank> ret = xt::empty<T>(shape);
 
     for (size_t e = 0; e < m_coarse2fine.shape(0); ++e) {
         for (size_t i = 0; i < m_coarse2fine.shape(1); ++i) {
             xt::view(ret, m_coarse2fine(e, i)) = xt::view(data, e);
         }
     }
 
     return ret;
 }
 
 inline FineLayer2Regular::FineLayer2Regular(const GooseFEM::Mesh::Quad4::FineLayer& mesh)
     : m_finelayer(mesh)
 {
     // ------------
     // Regular-mesh
     // ------------
 
     m_regular = GooseFEM::Mesh::Quad4::Regular(
         xt::amax(m_finelayer.m_nelx)(), xt::sum(m_finelayer.m_nhy)(), m_finelayer.m_h);
 
     // -------
     // mapping
     // -------
 
     // allocate mapping
     m_elem_regular.resize(m_finelayer.m_nelem);
     m_frac_regular.resize(m_finelayer.m_nelem);
 
     // alias
     xt::xtensor<size_t, 1> nhx = m_finelayer.m_nhx;
     xt::xtensor<size_t, 1> nhy = m_finelayer.m_nhy;
     xt::xtensor<size_t, 1> nelx = m_finelayer.m_nelx;
     xt::xtensor<size_t, 1> start = m_finelayer.m_startElem;
 
     // 'matrix' of element numbers of the Regular-mesh
     xt::xtensor<size_t, 2> elementgrid = m_regular.elementgrid();
 
     // cumulative number of element-rows of the Regular-mesh per layer of the FineLayer-mesh
     xt::xtensor<size_t, 1> cum_nhy =
         xt::concatenate(xt::xtuple(xt::zeros<size_t>({1}), xt::cumsum(nhy)));
 
     // number of element layers in y-direction of the FineLayer-mesh
     size_t nely = nhy.size();
 
     // loop over layers of the FineLayer-mesh
     for (size_t iy = 0; iy < nely; ++iy) {
         // element numbers of the Regular-mesh along this layer of the FineLayer-mesh
         auto el_new = xt::view(elementgrid, xt::range(cum_nhy(iy), cum_nhy(iy + 1)), xt::all());
 
         // no coarsening/refinement
         // ------------------------
 
         if (m_finelayer.m_refine(iy) == -1) {
             // element numbers of the FineLayer-mesh along this layer
             xt::xtensor<size_t, 1> el_old = start(iy) + xt::arange<size_t>(nelx(iy));
 
             // loop along this layer of the FineLayer-mesh
             for (size_t ix = 0; ix < nelx(iy); ++ix) {
                 // get the element numbers of the Regular-mesh for this element of the
                 // FineLayer-mesh
                 auto block =
                     xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
 
                 // write to mapping
                 for (auto& i : block) {
                     m_elem_regular[el_old(ix)].push_back(i);
                     m_frac_regular[el_old(ix)].push_back(1.0);
                 }
             }
         }
 
         // refinement along the x-direction (below the middle layer)
         // ---------------------------------------------------------
 
         else if (m_finelayer.m_refine(iy) == 0 && iy <= (nely - 1) / 2) {
             // element numbers of the FineLayer-mesh along this layer
             // rows: coarse block, columns element numbers per block
             xt::xtensor<size_t, 2> el_old =
                 start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});
 
             // loop along this layer of the FineLayer-mesh
             for (size_t ix = 0; ix < nelx(iy); ++ix) {
                 // get the element numbers of the Regular-mesh for this block of the FineLayer-mesh
                 auto block =
                     xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
 
                 // bottom: wide-to-narrow
                 {
                     for (size_t j = 0; j < nhy(iy) / 2; ++j) {
                         auto e = xt::view(block, j, xt::range(j, nhx(iy) - j));
 
                         m_elem_regular[el_old(ix, 0)].push_back(e(0));
                         m_frac_regular[el_old(ix, 0)].push_back(0.5);
 
                         for (size_t k = 1; k < e.size() - 1; ++k) {
                             m_elem_regular[el_old(ix, 0)].push_back(e(k));
                             m_frac_regular[el_old(ix, 0)].push_back(1.0);
                         }
 
                         m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1));
                         m_frac_regular[el_old(ix, 0)].push_back(0.5);
                     }
                 }
 
                 // top: regular small element
                 {
                     auto e = xt::view(
                         block,
                         xt::range(nhy(iy) / 2, nhy(iy)),
                         xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));
 
                     for (auto& i : e) {
                         m_elem_regular[el_old(ix, 2)].push_back(i);
                         m_frac_regular[el_old(ix, 2)].push_back(1.0);
                     }
                 }
 
                 // left
                 {
                     // left-bottom: narrow-to-wide
                     for (size_t j = 0; j < nhy(iy) / 2; ++j) {
                         auto e = xt::view(block, j, xt::range(0, j + 1));
 
                         for (size_t k = 0; k < e.size() - 1; ++k) {
                             m_elem_regular[el_old(ix, 3)].push_back(e(k));
                             m_frac_regular[el_old(ix, 3)].push_back(1.0);
                         }
 
                         m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1));
                         m_frac_regular[el_old(ix, 3)].push_back(0.5);
                     }
 
                     // left-top: regular
                     {
                         auto e = xt::view(
                             block,
                             xt::range(nhy(iy) / 2, nhy(iy)),
                             xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));
 
                         for (auto& i : e) {
                             m_elem_regular[el_old(ix, 3)].push_back(i);
                             m_frac_regular[el_old(ix, 3)].push_back(1.0);
                         }
                     }
                 }
 
                 // right
                 {
                     // right-bottom: narrow-to-wide
                     for (size_t j = 0; j < nhy(iy) / 2; ++j) {
                         auto e = xt::view(block, j, xt::range(nhx(iy) - j - 1, nhx(iy)));
 
                         m_elem_regular[el_old(ix, 1)].push_back(e(0));
                         m_frac_regular[el_old(ix, 1)].push_back(0.5);
 
                         for (size_t k = 1; k < e.size(); ++k) {
                             m_elem_regular[el_old(ix, 1)].push_back(e(k));
                             m_frac_regular[el_old(ix, 1)].push_back(1.0);
                         }
                     }
 
                     // right-top: regular
                     {
                         auto e = xt::view(
                             block,
                             xt::range(nhy(iy) / 2, nhy(iy)),
                             xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));
 
                         for (auto& i : e) {
                             m_elem_regular[el_old(ix, 1)].push_back(i);
                             m_frac_regular[el_old(ix, 1)].push_back(1.0);
                         }
                     }
                 }
             }
         }
 
         // coarsening along the x-direction (above the middle layer)
         else if (m_finelayer.m_refine(iy) == 0 && iy > (nely - 1) / 2) {
             // element numbers of the FineLayer-mesh along this layer
             // rows: coarse block, columns element numbers per block
             xt::xtensor<size_t, 2> el_old =
                 start(iy) + xt::arange<size_t>(nelx(iy) * 4ul).reshape({-1, 4});
 
             // loop along this layer of the FineLayer-mesh
             for (size_t ix = 0; ix < nelx(iy); ++ix) {
                 // get the element numbers of the Regular-mesh for this block of the FineLayer-mesh
                 auto block =
                     xt::view(el_new, xt::all(), xt::range(ix * nhx(iy), (ix + 1) * nhx(iy)));
 
                 // top: narrow-to-wide
                 {
                     for (size_t j = 0; j < nhy(iy) / 2; ++j) {
                         auto e = xt::view(
                             block,
                             nhy(iy) / 2 + j,
                             xt::range(1 * nhx(iy) / 3 - j - 1, 2 * nhx(iy) / 3 + j + 1));
 
                         m_elem_regular[el_old(ix, 3)].push_back(e(0));
                         m_frac_regular[el_old(ix, 3)].push_back(0.5);
 
                         for (size_t k = 1; k < e.size() - 1; ++k) {
                             m_elem_regular[el_old(ix, 3)].push_back(e(k));
                             m_frac_regular[el_old(ix, 3)].push_back(1.0);
                         }
 
                         m_elem_regular[el_old(ix, 3)].push_back(e(e.size() - 1));
                         m_frac_regular[el_old(ix, 3)].push_back(0.5);
                     }
                 }
 
                 // bottom: regular small element
                 {
                     auto e = xt::view(
                         block,
                         xt::range(0, nhy(iy) / 2),
                         xt::range(1 * nhx(iy) / 3, 2 * nhx(iy) / 3));
 
                     for (auto& i : e) {
                         m_elem_regular[el_old(ix, 1)].push_back(i);
                         m_frac_regular[el_old(ix, 1)].push_back(1.0);
                     }
                 }
 
                 // left
                 {
                     // left-bottom: regular
                     {
                         auto e = xt::view(
                             block,
                             xt::range(0, nhy(iy) / 2),
                             xt::range(0 * nhx(iy) / 3, 1 * nhx(iy) / 3));
 
                         for (auto& i : e) {
                             m_elem_regular[el_old(ix, 0)].push_back(i);
                             m_frac_regular[el_old(ix, 0)].push_back(1.0);
                         }
                     }
 
                     // left-top: narrow-to-wide
                     for (size_t j = 0; j < nhy(iy) / 2; ++j) {
                         auto e =
                             xt::view(block, nhy(iy) / 2 + j, xt::range(0, 1 * nhx(iy) / 3 - j));
 
                         for (size_t k = 0; k < e.size() - 1; ++k) {
                             m_elem_regular[el_old(ix, 0)].push_back(e(k));
                             m_frac_regular[el_old(ix, 0)].push_back(1.0);
                         }
 
                         m_elem_regular[el_old(ix, 0)].push_back(e(e.size() - 1));
                         m_frac_regular[el_old(ix, 0)].push_back(0.5);
                     }
                 }
 
                 // right
                 {
                     // right-bottom: regular
                     {
                         auto e = xt::view(
                             block,
                             xt::range(0, nhy(iy) / 2),
                             xt::range(2 * nhx(iy) / 3, 3 * nhx(iy) / 3));
 
                         for (auto& i : e) {
                             m_elem_regular[el_old(ix, 2)].push_back(i);
                             m_frac_regular[el_old(ix, 2)].push_back(1.0);
                         }
                     }
 
                     // right-top: narrow-to-wide
                     for (size_t j = 0; j < nhy(iy) / 2; ++j) {
                         auto e = xt::view(
                             block, nhy(iy) / 2 + j, xt::range(2 * nhx(iy) / 3 + j, nhx(iy)));
 
                         m_elem_regular[el_old(ix, 2)].push_back(e(0));
                         m_frac_regular[el_old(ix, 2)].push_back(0.5);
 
                         for (size_t k = 1; k < e.size(); ++k) {
                             m_elem_regular[el_old(ix, 2)].push_back(e(k));
                             m_frac_regular[el_old(ix, 2)].push_back(1.0);
                         }
                     }
                 }
             }
         }
     }
 }
 
 inline GooseFEM::Mesh::Quad4::Regular FineLayer2Regular::getRegularMesh() const
 {
     return m_regular;
 }
 
 inline GooseFEM::Mesh::Quad4::FineLayer FineLayer2Regular::getFineLayerMesh() const
 {
     return m_finelayer;
 }
 
 inline std::vector<std::vector<size_t>> FineLayer2Regular::getMap() const
 {
     return m_elem_regular;
 }
 
 inline std::vector<std::vector<double>> FineLayer2Regular::getMapFraction() const
 {
     return m_frac_regular;
 }
 
 template <class T, size_t rank>
 inline xt::xtensor<T, rank> FineLayer2Regular::mapToRegular(const xt::xtensor<T, rank>& data) const
 {
     GOOSEFEM_ASSERT(data.shape(0) == m_finelayer.nelem());
 
     std::array<size_t, rank> shape;
     std::copy(data.shape().cbegin(), data.shape().cend(), &shape[0]);
     shape[0] = m_regular.nelem();
 
     xt::xtensor<T, rank> ret = xt::zeros<T>(shape);
 
     for (size_t e = 0; e < m_finelayer.nelem(); ++e) {
         for (size_t i = 0; i < m_elem_regular[e].size(); ++i) {
             xt::view(ret, m_elem_regular[e][i]) += m_frac_regular[e][i] * xt::view(data, e);
         }
     }
 
     return ret;
 }
 
 } // namespace Map
 
 } // namespace Quad4
 } // namespace Mesh
 } // namespace GooseFEM
 
 #endif