diff --git a/GooseFEM.pc.in b/GooseFEM.pc.in
index 4ba0cf4..0728150 100644
--- a/GooseFEM.pc.in
+++ b/GooseFEM.pc.in
@@ -1,7 +1,7 @@
 prefix=@CMAKE_INSTALL_PREFIX@
 includedir=${prefix}/include
 
 Name: @PROJECT_NAME@
 Description: Finite Element simulation, some simple meshes and operations.
-Version: @GOOSEFEM_VERSION@
+Version: @GooseFEM_VERSION@
 Cflags: -I${includedir}
diff --git a/include/GooseFEM/Mesh.h b/include/GooseFEM/Mesh.h
index 30e3c5b..46c3790 100644
--- a/include/GooseFEM/Mesh.h
+++ b/include/GooseFEM/Mesh.h
@@ -1,570 +1,590 @@
 /**
 Generic mesh operations.
 
 \file Mesh.h
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_MESH_H
 #define GOOSEFEM_MESH_H
 
 #include "config.h"
 
 namespace GooseFEM {
 namespace Mesh {
 
 
 /**
 Enumerator for element-types
 */
 enum class ElementType {
     Quad4, ///< Quadrilateral: 4-noded element in 2-d
     Hex8, ///< Hexahedron: 8-noded element in 3-d
     Tri3 ///< Triangle: 3-noded element in 2-d
 };
 
 /**
 Extract the element type based on the connectivity.
 
 \param coor Nodal coordinates.
 \param conn Connectivity.
 \return ElementType().
 */
 inline ElementType defaultElementType(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn);
 
 /**
 Find overlapping nodes. The output has the following structure:
 
     [[nodes_from_mesh_a],
      [nodes_from_mesh_b]]
 
 \param coor_a Nodal coordinates of mesh "a".
 \param coor_b Nodal coordinates of mesh "b".
 \param rtol Relative tolerance for position match.
 \param atol Absolute tolerance for position match.
 \return Overlapping nodes.
 */
 inline xt::xtensor<size_t, 2> overlapping(
     const xt::xtensor<double, 2>& coor_a,
     const xt::xtensor<double, 2>& coor_b,
     double rtol = 1e-5,
     double atol = 1e-8);
 
 /**
 Stitch two mesh objects, specifying overlapping nodes by hand.
 */
 class ManualStitch {
 public:
     ManualStitch() = default;
 
     /**
     \param coor_a Nodal coordinates of mesh "a".
     \param conn_a Connectivity of mesh "a".
     \param overlapping_nodes_a Node-numbers of mesh "a" that overlap with mesh "b".
     \param coor_b Nodal coordinates of mesh "b".
     \param conn_b Connectivity of mesh "b".
     \param overlapping_nodes_b Node-numbers of mesh "b" that overlap with mesh "a".
     \param check_position If ``true`` the nodes are checked for position overlap.
     \param rtol Relative tolerance for check on position overlap.
     \param atol Absolute tolerance for check on position overlap.
     */
     ManualStitch(
         const xt::xtensor<double, 2>& coor_a,
         const xt::xtensor<size_t, 2>& conn_a,
         const xt::xtensor<size_t, 1>& overlapping_nodes_a,
         const xt::xtensor<double, 2>& coor_b,
         const xt::xtensor<size_t, 2>& conn_b,
         const xt::xtensor<size_t, 1>& overlapping_nodes_b,
         bool check_position = true,
         double rtol = 1e-5,
         double atol = 1e-8);
 
     /**
     Number of sub meshes.
 
     \return 2.
     */
     size_t nmesh() const;
 
     /**
     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;
 
     /**
     Nodal coordinates.
 
     \return [#nnode, #ndim].
     */
     xt::xtensor<double, 2> coor() const;
 
     /**
     Connectivity.
 
     \return [#nelem, #nne].
     */
     xt::xtensor<size_t, 2> conn() const;
 
     /**
     DOF numbers for each node (numbered sequentially).
 
     \return [#nnode, #ndim].
     */
     xt::xtensor<size_t, 2> dofs() const;
 
+    /**
+    \return Node-map per sub-mesh.
+    */
+    std::vector<xt::xtensor<size_t, 1>> nodemap() const;
+
+    /**
+    \return Element-map per sub-mesh.
+    */
+    std::vector<xt::xtensor<size_t, 1>> elemmap() const;
+
     /**
     \param mesh_index Index of the mesh ("a" = 1, "b" = 1).
     \return Node-map for a given mesh.
     */
     xt::xtensor<size_t, 1> nodemap(size_t mesh_index) const;
 
     /**
     \param mesh_index Index of the mesh ("a" = 1, "b" = 1).
     \return Element-map for a given mesh.
     */
     xt::xtensor<size_t, 1> elemmap(size_t mesh_index) const;
 
     /**
     Convert set of node numbers for an original mesh to the stitched mesh.
 
     \param set List of node numbers.
     \param mesh_index Index of the mesh ("a" = 1, "b" = 1).
     \return List of node numbers for the stitched mesh.
     */
     xt::xtensor<size_t, 1> nodeset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const;
 
     /**
     Convert set of element numbers for an original mesh to the stitched mesh.
 
     \param set List of element numbers.
     \param mesh_index Index of the mesh ("a" = 1, "b" = 1).
     \return List of element numbers for the stitched mesh.
     */
     xt::xtensor<size_t, 1> elemset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const;
 
 private:
     xt::xtensor<double, 2> m_coor;
     xt::xtensor<size_t, 2> m_conn;
     xt::xtensor<size_t, 1> m_map_b;
     size_t m_nnd_a;
     size_t m_nel_a;
     size_t m_nel_b;
 };
 
 /**
 Stitch mesh objects, automatically searching for overlapping nodes.
 */
 class Stitch {
 public:
     /**
     \param rtol Relative tolerance for position match.
     \param atol Absolute tolerance for position match.
     */
     Stitch(double rtol = 1e-5, double atol = 1e-8);
 
     /**
     Add mesh to be stitched.
 
     \param coor Nodal coordinates.
     \param conn Connectivity.
     */
     void push_back(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn);
 
     /**
     Number of sub meshes.
 
     \return unsigned int
     */
     size_t nmesh() const;
 
     /**
     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;
 
     /**
     Nodal coordinates.
 
     \return [#nnode, #ndim].
     */
     xt::xtensor<double, 2> coor() const;
 
     /**
     Connectivity.
 
     \return [#nelem, #nne].
     */
     xt::xtensor<size_t, 2> conn() const;
 
     /**
     DOF numbers for each node (numbered sequentially).
 
     \return [#nnode, #ndim].
     */
     xt::xtensor<size_t, 2> dofs() const;
 
+    /**
+    \return Node-map per sub-mesh.
+    */
+    std::vector<xt::xtensor<size_t, 1>> nodemap() const;
+
+    /**
+    \return Element-map per sub-mesh.
+    */
+    std::vector<xt::xtensor<size_t, 1>> elemmap() const;
+
     /**
     The node numbers in the stitched mesh that are coming from a specific sub-mesh.
 
     \param mesh_index Index of the sub-mesh.
     \return List of node numbers.
     */
     xt::xtensor<size_t, 1> nodemap(size_t mesh_index) const;
 
     /**
     The element numbers in the stitched mesh that are coming from a specific sub-mesh.
 
     \param mesh_index Index of the sub-mesh.
     \return List of element numbers.
     */
     xt::xtensor<size_t, 1> elemmap(size_t mesh_index) const;
 
     /**
     Convert set of node-numbers for a sub-mesh to the stitched mesh.
 
     \param set List of node numbers.
     \param mesh_index Index of the sub-mesh.
     \return List of node numbers for the stitched mesh.
     */
     xt::xtensor<size_t, 1> nodeset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const;
 
     /**
     Convert set of element-numbers for a sub-mesh to the stitched mesh.
 
     \param set List of element numbers.
     \param mesh_index Index of the sub-mesh.
     \return List of element numbers for the stitched mesh.
     */
     xt::xtensor<size_t, 1> elemset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const;
 
     /**
     Combine set of node numbers for an original to the final mesh (removes duplicates).
 
     \param set List of node numbers per mesh.
     \return List of node numbers for the stitched mesh.
     */
     xt::xtensor<size_t, 1> nodeset(const std::vector<xt::xtensor<size_t, 1>>& set) const;
 
     /**
     Combine set of element numbers for an original to the final mesh.
 
     \param set List of element numbers per mesh.
     \return List of element numbers for the stitched mesh.
     */
     xt::xtensor<size_t, 1> elemset(const std::vector<xt::xtensor<size_t, 1>>& set) const;
 
 private:
     xt::xtensor<double, 2> m_coor;
     xt::xtensor<size_t, 2> m_conn;
     std::vector<xt::xtensor<size_t, 1>> m_map;
     std::vector<size_t> m_nel; ///< Number of elements per sub-mesh.
     std::vector<size_t> m_el_offset;
     double m_rtol;
     double m_atol;
 };
 
 /**
 \rst
 Renumber indices to lowest possible index. For example:
 
 .. math::
 
     \begin{bmatrix}
         0 & 1 \\
         5 & 4
     \end{bmatrix}
 
 is renumbered to
 
 .. math::
 
     \begin{bmatrix}
         0 & 1 \\
         3 & 2
     \end{bmatrix}
 
 Or, in pseudo-code, the result of this function is that:
 
 .. code-block:: python
 
     dofs = renumber(dofs)
 
     sort(unique(dofs[:])) == range(max(dofs+1))
 
 .. tip::
 
     One can use the wrapper function :cpp:func:`GooseFEM::Mesh::renumber`.
     This class gives more advanced features.
 \endrst
 */
 class Renumber {
 public:
     Renumber() = default;
 
     /**
     \param dofs DOF-numbers.
     */
     template <class T>
     Renumber(const T& dofs);
 
     /**
     Get renumbered DOFs (same as ``Renumber::apply(dofs)``).
 
     \param dofs List of (DOF-)numbers.
     \return Renumbered list of (DOF-)numbers.
     */
     [[deprecated]]
     xt::xtensor<size_t, 2> get(const xt::xtensor<size_t, 2>& dofs) const;
 
     /**
     Apply renumbering to other set.
 
     \param list List of (DOF-)numbers.
     \return Renumbered list of (DOF-)numbers.
     */
     template <class T>
     T apply(const T& list) const;
 
     /**
     Get the list needed to renumber, e.g.:
 
         dofs_renumbered(i, j) = index(dofs(i, j))
 
     \return Renumber-index.
     */
     xt::xtensor<size_t, 1> index() const;
 
 private:
     xt::xtensor<size_t, 1> m_renum;
 };
 
 /**
 Renumber to lowest possible index (see GooseFEM::Mesh::Renumber).
 
 \param dofs DOF-numbers.
 \return Renumbered DOF-numbers.
 */
 inline xt::xtensor<size_t, 2> renumber(const xt::xtensor<size_t, 2>& dofs);
 
 /**
 Reorder to lowest possible index, in specific order.
 
 For example for ``Reorder({iiu, iip})`` after reordering:
 
     iiu = xt::range<size_t>(nnu);
     iip = xt::range<size_t>(nnp) + nnu;
 */
 class Reorder {
 public:
     Reorder() = default;
 
     /**
     \param args List of (DOF-)numbers.
     */
     Reorder(const std::initializer_list<xt::xtensor<size_t, 1>> args);
 
     /**
     Get reordered DOFs (same as ``Reorder::apply(dofs)``).
 
     \param dofs List of (DOF-)numbers.
     \return Reordered list of (DOF-)numbers.
     */
     [[deprecated]]
     xt::xtensor<size_t, 2> get(const xt::xtensor<size_t, 2>& dofs) const;
 
     /**
     Apply reordering to other set.
 
     \param list List of (DOF-)numbers.
     \return Reordered list of (DOF-)numbers.
     */
     template <class T>
     T apply(const T& list) const;
 
     /**
     Get the list needed to reorder, e.g.:
 
         dofs_reordered(i, j) = index(dofs(i, j))
 
     \return Reorder-index.
     */
     xt::xtensor<size_t, 1> index() const;
 
 private:
     xt::xtensor<size_t, 1> m_renum;
 };
 
 /**
 List with DOF-numbers in sequential order.
 The output is a sequential list of DOF-numbers for each vector-component of each node.
 For example for 3 nodes in 2 dimensions the output is
 
 \rst
 .. math::
 
     \begin{bmatrix}
         0 & 1 \\
         2 & 3 \\
         4 & 5
     \end{bmatrix}
 \endrst
 
 \param nnode Number of nodes.
 \param ndim Number of dimensions.
 \return DOF-numbers.
 */
 inline xt::xtensor<size_t, 2> dofs(size_t nnode, size_t ndim);
 
 /**
 Number of elements connected to each node.
 
 \param conn Connectivity.
 \return Coordination per node.
 */
 inline xt::xtensor<size_t, 1> coordination(const xt::xtensor<size_t, 2>& conn);
 
 /**
 Elements connected to each node.
 
 \param conn Connectivity.
 \param sorted If ``true`` the output is sorted.
 \return Elements per node.
 */
 inline std::vector<std::vector<size_t>> elem2node(
     const xt::xtensor<size_t, 2>& conn,
     bool sorted = true);
 
 /**
 Return size of each element edge.
 
 \param coor Nodal coordinates.
 \param conn Connectivity.
 \param type ElementType.
 \return Edge-sizes per element.
 */
 inline xt::xtensor<double, 2> edgesize(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn,
     ElementType type);
 
 /**
 Return size of each element edge.
 The element-type is automatically determined, see defaultElementType().
 
 \param coor Nodal coordinates.
 \param conn Connectivity.
 \return Edge-sizes per element.
 */
 inline xt::xtensor<double, 2> edgesize(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn);
 
 /**
 Coordinates of the center of each element.
 
 \param coor Nodal coordinates.
 \param conn Connectivity.
 \param type ElementType.
 \return Center of each element.
 */
 inline xt::xtensor<double, 2> centers(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn,
     ElementType type);
 
 /**
 Coordinates of the center of each element.
 The element-type is automatically determined, see defaultElementType().
 
 \param coor Nodal coordinates.
 \param conn Connectivity.
 \return Center of each element.
 */
 inline xt::xtensor<double, 2> centers(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn);
 
 /**
 Convert an element-map to a node-map.
 
 \param elem_map Element-map such that ``new_elvar = elvar[elem_map]``.
 \param coor Nodal coordinates.
 \param conn Connectivity.
 \param type ElementType.
 \return Node-map such that ``new_nodevar = nodevar[node_map]``
 */
 inline xt::xtensor<size_t, 1> elemmap2nodemap(
     const xt::xtensor<size_t, 1>& elem_map,
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn,
     ElementType type);
 
 /**
 Convert an element-map to a node-map.
 The element-type is automatically determined, see defaultElementType().
 
 \param elem_map Element-map such that ``new_elvar = elvar[elem_map]``.
 \param coor Nodal coordinates.
 \param conn Connectivity.
 \return Node-map such that ``new_nodevar = nodevar[node_map]``
 */
 inline xt::xtensor<size_t, 1> elemmap2nodemap(
     const xt::xtensor<size_t, 1>& elem_map,
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn);
 
 } // namespace Mesh
 } // namespace GooseFEM
 
 #include "Mesh.hpp"
 
 #endif
diff --git a/include/GooseFEM/Mesh.hpp b/include/GooseFEM/Mesh.hpp
index e8f5cec..e01110e 100644
--- a/include/GooseFEM/Mesh.hpp
+++ b/include/GooseFEM/Mesh.hpp
@@ -1,620 +1,656 @@
 /**
 Implementation of Mesh.h
 
 \file Mesh.hpp
 \copyright Copyright 2017. Tom de Geus. All rights reserved.
 \license This project is released under the GNU Public License (GPLv3).
 */
 
 #ifndef GOOSEFEM_MESH_HPP
 #define GOOSEFEM_MESH_HPP
 
 #include "Mesh.h"
 
 namespace GooseFEM {
 namespace Mesh {
 
 inline ElementType defaultElementType(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn)
 {
     if (coor.shape(1) == 2ul && conn.shape(1) == 3ul) {
         return ElementType::Tri3;
     }
     if (coor.shape(1) == 2ul && conn.shape(1) == 4ul) {
         return ElementType::Quad4;
     }
     if (coor.shape(1) == 3ul && conn.shape(1) == 8ul) {
         return ElementType::Hex8;
     }
 
     throw std::runtime_error("Element-type not implemented");
 }
 
 namespace detail {
 
     template <class T, class R>
     inline T renum(const T& arg, const R& mapping)
     {
         T ret = T::from_shape(arg.shape());
 
         auto jt = ret.begin();
 
         for (auto it = arg.begin(); it != arg.end(); ++it, ++jt) {
             *jt = mapping(*it);
         }
 
         return ret;
     }
 
 } // namespace detail
 
 inline xt::xtensor<size_t, 2> overlapping(
     const xt::xtensor<double, 2>& coor_a,
     const xt::xtensor<double, 2>& coor_b,
     double rtol,
     double atol)
 {
     GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1));
 
     std::vector<size_t> ret_a;
     std::vector<size_t> ret_b;
 
     for (size_t i = 0; i < coor_a.shape(0); ++i) {
 
         auto idx = xt::flatten_indices(xt::argwhere(xt::prod(xt::isclose(
             coor_b, xt::view(coor_a, i, xt::all()), rtol, atol), 1)));
 
         for (auto& j : idx) {
             ret_a.push_back(i);
             ret_b.push_back(j);
         }
     }
 
     xt::xtensor<size_t, 2> ret = xt::empty<size_t>({size_t(2), ret_a.size()});
     for (size_t i = 0; i < ret_a.size(); ++i) {
         ret(0, i) = ret_a[i];
         ret(1, i) = ret_b[i];
     }
 
     return ret;
 }
 
 inline ManualStitch::ManualStitch(
     const xt::xtensor<double, 2>& coor_a,
     const xt::xtensor<size_t, 2>& conn_a,
     const xt::xtensor<size_t, 1>& overlapping_nodes_a,
     const xt::xtensor<double, 2>& coor_b,
     const xt::xtensor<size_t, 2>& conn_b,
     const xt::xtensor<size_t, 1>& overlapping_nodes_b,
     bool check_position,
     double rtol,
     double atol)
 {
     UNUSED(rtol);
     UNUSED(atol);
 
     GOOSEFEM_ASSERT(xt::has_shape(overlapping_nodes_a, overlapping_nodes_b.shape()));
     GOOSEFEM_ASSERT(coor_a.shape(1) == coor_b.shape(1));
     GOOSEFEM_ASSERT(conn_a.shape(1) == conn_b.shape(1));
 
     if (check_position) {
         GOOSEFEM_ASSERT(xt::allclose(
             xt::view(coor_a, xt::keep(overlapping_nodes_a), xt::all()),
             xt::view(coor_b, xt::keep(overlapping_nodes_b), xt::all()),
             rtol,
             atol));
     }
 
     size_t nnda = coor_a.shape(0);
     size_t nndb = coor_b.shape(0);
     size_t ndim = coor_a.shape(1);
     size_t nelim = overlapping_nodes_a.size();
 
     size_t nela = conn_a.shape(0);
     size_t nelb = conn_b.shape(0);
     size_t nne = conn_a.shape(1);
 
     m_nel_a = nela;
     m_nel_b = nelb;
     m_nnd_a = nnda;
 
     xt::xtensor<size_t, 1> keep_b = xt::setdiff1d(xt::arange<size_t>(nndb), overlapping_nodes_b);
 
     m_map_b = xt::empty<size_t>({nndb});
     xt::view(m_map_b, xt::keep(overlapping_nodes_b)) = overlapping_nodes_a;
     xt::view(m_map_b, xt::keep(keep_b)) = xt::arange<size_t>(keep_b.size()) + nnda;
 
     m_conn = xt::empty<size_t>({nela + nelb, nne});
     xt::view(m_conn, xt::range(0, nela), xt::all()) = conn_a;
     xt::view(m_conn, xt::range(nela, nela + nelb), xt::all()) = detail::renum(conn_b, m_map_b);
 
     m_coor = xt::empty<size_t>({nnda + nndb - nelim, ndim});
     xt::view(m_coor, xt::range(0, nnda), xt::all()) = coor_a;
     xt::view(m_coor, xt::range(nnda, nnda + nndb - nelim), xt::all()) =
         xt::view(coor_b, xt::keep(keep_b), xt::all());
 }
 
 inline xt::xtensor<double, 2> ManualStitch::coor() const
 {
     return m_coor;
 }
 
 inline xt::xtensor<size_t, 2> ManualStitch::conn() const
 {
     return m_conn;
 }
 
 inline size_t ManualStitch::nmesh() const
 {
     return 2;
 }
 
 inline size_t ManualStitch::nelem() const
 {
     return m_conn.shape(0);
 }
 
 inline size_t ManualStitch::nnode() const
 {
     return m_coor.shape(0);
 }
 
 inline size_t ManualStitch::nne() const
 {
     return m_conn.shape(1);
 }
 
 inline size_t ManualStitch::ndim() const
 {
     return m_coor.shape(1);
 }
 
 inline xt::xtensor<size_t, 2> ManualStitch::dofs() const
 {
     size_t nnode = this->nnode();
     size_t ndim = this->ndim();
     return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
 }
 
+inline std::vector<xt::xtensor<size_t, 1>> ManualStitch::nodemap() const
+{
+    std::vector<xt::xtensor<size_t, 1>> ret(this->nmesh());
+    for (size_t i = 0; i < this->nmesh(); ++i) {
+        ret[i] = this->nodemap(i);
+    }
+    return ret;
+}
+
+inline std::vector<xt::xtensor<size_t, 1>> ManualStitch::elemmap() const
+{
+    std::vector<xt::xtensor<size_t, 1>> ret(this->nmesh());
+    for (size_t i = 0; i < this->nmesh(); ++i) {
+        ret[i] = this->elemmap(i);
+    }
+    return ret;
+}
+
 inline xt::xtensor<size_t, 1> ManualStitch::nodemap(size_t mesh_index) const
 {
     GOOSEFEM_ASSERT(mesh_index <= 1);
 
     if (mesh_index == 0) {
         return xt::arange<size_t>(m_nnd_a);
     }
 
     return m_map_b;
 }
 
 inline xt::xtensor<size_t, 1> ManualStitch::elemmap(size_t mesh_index) const
 {
     GOOSEFEM_ASSERT(mesh_index <= 1);
 
     if (mesh_index == 0) {
         return xt::arange<size_t>(m_nel_a);
     }
 
     return xt::arange<size_t>(m_nel_b) + m_nel_a;
 }
 
 inline xt::xtensor<size_t, 1>
 ManualStitch::nodeset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const
 {
     GOOSEFEM_ASSERT(mesh_index <= 1);
 
     if (mesh_index == 0) {
         GOOSEFEM_ASSERT(xt::amax(set)() < m_nnd_a);
         return set;
     }
 
     GOOSEFEM_ASSERT(xt::amax(set)() < m_map_b.size());
     return detail::renum(set, m_map_b);
 }
 
 inline xt::xtensor<size_t, 1>
 ManualStitch::elemset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const
 {
     GOOSEFEM_ASSERT(mesh_index <= 1);
 
     if (mesh_index == 0) {
         GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_a);
         return set;
     }
 
     GOOSEFEM_ASSERT(xt::amax(set)() < m_nel_b);
     return set + m_nel_a;
 }
 
 inline Stitch::Stitch(double rtol, double atol)
 {
     m_rtol = rtol;
     m_atol = atol;
 }
 
 inline void Stitch::push_back(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn)
 {
     if (m_map.size() == 0) {
         m_coor = coor;
         m_conn = conn;
         m_map.push_back(xt::eval(xt::arange<size_t>(coor.shape(0))));
         m_nel.push_back(conn.shape(0));
         m_el_offset.push_back(0);
         return;
     }
 
     auto overlap = overlapping(m_coor, coor, m_rtol, m_atol);
     size_t index = m_map.size();
 
     ManualStitch stich(
         m_coor, m_conn, xt::view(overlap, 0, xt::all()),
         coor, conn, xt::view(overlap, 1, xt::all()),
         false);
 
     m_coor = stich.coor();
     m_conn = stich.conn();
     m_map.push_back(stich.nodemap(1));
     m_nel.push_back(conn.shape(0));
     m_el_offset.push_back(m_el_offset[index - 1] + m_nel[index - 1]);
 }
 
 inline size_t Stitch::nmesh() const
 {
     return m_map.size();
 }
 
 inline xt::xtensor<double, 2> Stitch::coor() const
 {
     return m_coor;
 }
 
 inline xt::xtensor<size_t, 2> Stitch::conn() const
 {
     return m_conn;
 }
 
 inline size_t Stitch::nelem() const
 {
     return m_conn.shape(0);
 }
 
 inline size_t Stitch::nnode() const
 {
     return m_coor.shape(0);
 }
 
 inline size_t Stitch::nne() const
 {
     return m_conn.shape(1);
 }
 
 inline size_t Stitch::ndim() const
 {
     return m_coor.shape(1);
 }
 
 inline xt::xtensor<size_t, 2> Stitch::dofs() const
 {
     size_t nnode = this->nnode();
     size_t ndim = this->ndim();
     return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
 }
 
+inline std::vector<xt::xtensor<size_t, 1>> Stitch::nodemap() const
+{
+    std::vector<xt::xtensor<size_t, 1>> ret(this->nmesh());
+    for (size_t i = 0; i < this->nmesh(); ++i) {
+        ret[i] = this->nodemap(i);
+    }
+    return ret;
+}
+
+inline std::vector<xt::xtensor<size_t, 1>> Stitch::elemmap() const
+{
+    std::vector<xt::xtensor<size_t, 1>> ret(this->nmesh());
+    for (size_t i = 0; i < this->nmesh(); ++i) {
+        ret[i] = this->elemmap(i);
+    }
+    return ret;
+}
+
 inline xt::xtensor<size_t, 1> Stitch::nodemap(size_t mesh_index) const
 {
     GOOSEFEM_ASSERT(mesh_index < m_map.size());
     return m_map[mesh_index];
 }
 
 inline xt::xtensor<size_t, 1> Stitch::elemmap(size_t mesh_index) const
 {
     GOOSEFEM_ASSERT(mesh_index < m_map.size());
     return xt::arange<size_t>(m_nel[mesh_index]) + m_el_offset[mesh_index];
 }
 
 inline xt::xtensor<size_t, 1> Stitch::nodeset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const
 {
     GOOSEFEM_ASSERT(mesh_index < m_map.size());
     GOOSEFEM_ASSERT(xt::amax(set)() < m_map[mesh_index].size());
     return detail::renum(set, m_map[mesh_index]);
 }
 
 inline xt::xtensor<size_t, 1> Stitch::elemset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const
 {
     GOOSEFEM_ASSERT(mesh_index < m_map.size());
     GOOSEFEM_ASSERT(xt::amax(set)() < m_nel[mesh_index]);
     return set + m_el_offset[mesh_index];
 }
 
 inline xt::xtensor<size_t, 1> Stitch::nodeset(const std::vector<xt::xtensor<size_t, 1>>& set) const
 {
     GOOSEFEM_ASSERT(set.size() == m_map.size());
 
     size_t n = 0;
 
     for (size_t i = 0; i < set.size(); ++i) {
         n += set[i].size();
     }
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({n});
 
     n = 0;
 
     for (size_t i = 0; i < set.size(); ++i) {
         xt::view(ret, xt::range(n, n + set[i].size())) = this->nodeset(set[i], i);
         n += set[i].size();
     }
 
     return xt::unique(ret);
 }
 
 inline xt::xtensor<size_t, 1> Stitch::elemset(const std::vector<xt::xtensor<size_t, 1>>& set) const
 {
     GOOSEFEM_ASSERT(set.size() == m_map.size());
 
     size_t n = 0;
 
     for (size_t i = 0; i < set.size(); ++i) {
         n += set[i].size();
     }
 
     xt::xtensor<size_t, 1> ret = xt::empty<size_t>({n});
 
     n = 0;
 
     for (size_t i = 0; i < set.size(); ++i) {
         xt::view(ret, xt::range(n, n + set[i].size())) = this->elemset(set[i], i);
         n += set[i].size();
     }
 
     return ret;
 }
 
 template <class T>
 inline Renumber::Renumber(const T& dofs)
 {
     size_t n = xt::amax(dofs)() + 1;
     size_t i = 0;
 
     xt::xtensor<size_t, 1> unique = xt::unique(dofs);
 
     m_renum = xt::empty<size_t>({n});
 
     for (auto& j : unique) {
         m_renum(j) = i;
         ++i;
     }
 }
 
 inline xt::xtensor<size_t, 2> Renumber::get(const xt::xtensor<size_t, 2>& dofs) const
 {
     GOOSEFEM_WARNING("Renumber::get is deprecated, use Renumber::apply");
     return this->apply(dofs);
 }
 
 template <class T>
 inline T Renumber::apply(const T& list) const
 {
     return detail::renum(list, m_renum);
 }
 
 inline xt::xtensor<size_t, 1> Renumber::index() const
 {
     return m_renum;
 }
 
 inline xt::xtensor<size_t, 2> renumber(const xt::xtensor<size_t, 2>& dofs)
 {
     return Renumber(dofs).apply(dofs);
 }
 
 inline Reorder::Reorder(const std::initializer_list<xt::xtensor<size_t, 1>> args)
 {
     size_t n = 0;
     size_t i = 0;
 
     for (auto& arg : args) {
         if (arg.size() == 0) {
             continue;
         }
         n = std::max(n, xt::amax(arg)() + 1);
     }
 
     #ifdef GOOSEFEM_ENABLE_ASSERT
     for (auto& arg : args) {
         GOOSEFEM_ASSERT(xt::unique(arg) == xt::sort(arg));
     }
     #endif
 
     m_renum = xt::empty<size_t>({n});
 
     for (auto& arg : args) {
         for (auto& j : arg) {
             m_renum(j) = i;
             ++i;
         }
     }
 }
 
 inline xt::xtensor<size_t, 2> Reorder::get(const xt::xtensor<size_t, 2>& dofs) const
 {
     GOOSEFEM_WARNING("Reorder::get is deprecated, use Reorder::apply");
     return this->apply(dofs);
 }
 
 template <class T>
 inline T Reorder::apply(const T& list) const
 {
     T ret = T::from_shape(list.shape());
 
     auto jt = ret.begin();
 
     for (auto it = list.begin(); it != list.end(); ++it, ++jt) {
         *jt = m_renum(*it);
     }
 
     return ret;
 }
 
 inline xt::xtensor<size_t, 1> Reorder::index() const
 {
     return m_renum;
 }
 
 inline xt::xtensor<size_t, 2> dofs(size_t nnode, size_t ndim)
 {
     return xt::reshape_view(xt::arange<size_t>(nnode * ndim), {nnode, ndim});
 }
 
 inline xt::xtensor<size_t, 1> coordination(const xt::xtensor<size_t, 2>& conn)
 {
     size_t nnode = xt::amax(conn)() + 1;
 
     xt::xtensor<size_t, 1> N = xt::zeros<size_t>({nnode});
 
     for (auto it = conn.begin(); it != conn.end(); ++it) {
         N(*it) += 1;
     }
 
     return N;
 }
 
 inline std::vector<std::vector<size_t>> elem2node(const xt::xtensor<size_t, 2>& conn, bool sorted)
 {
     auto N = coordination(conn);
     auto nnode = N.size();
 
     std::vector<std::vector<size_t>> ret(nnode);
     for (size_t i = 0; i < nnode; ++i) {
         ret[i].reserve(N(i));
     }
 
     for (size_t e = 0; e < conn.shape(0); ++e) {
         for (size_t m = 0; m < conn.shape(1); ++m) {
             ret[conn(e, m)].push_back(e);
         }
     }
 
     if (sorted) {
         for (auto& row : ret) {
             std::sort(row.begin(), row.end());
         }
     }
 
     return ret;
 }
 
 inline xt::xtensor<double, 2> edgesize(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn,
     ElementType type)
 {
     GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
 
     if (type == ElementType::Quad4) {
         GOOSEFEM_ASSERT(coor.shape(1) == 2ul);
         GOOSEFEM_ASSERT(conn.shape(1) == 4ul);
         xt::xtensor<size_t, 1> n0 = xt::view(conn, xt::all(), 0);
         xt::xtensor<size_t, 1> n1 = xt::view(conn, xt::all(), 1);
         xt::xtensor<size_t, 1> n2 = xt::view(conn, xt::all(), 2);
         xt::xtensor<size_t, 1> n3 = xt::view(conn, xt::all(), 3);
         xt::xtensor<double, 1> x0 = xt::view(coor, xt::keep(n0), 0);
         xt::xtensor<double, 1> x1 = xt::view(coor, xt::keep(n1), 0);
         xt::xtensor<double, 1> x2 = xt::view(coor, xt::keep(n2), 0);
         xt::xtensor<double, 1> x3 = xt::view(coor, xt::keep(n3), 0);
         xt::xtensor<double, 1> y0 = xt::view(coor, xt::keep(n0), 1);
         xt::xtensor<double, 1> y1 = xt::view(coor, xt::keep(n1), 1);
         xt::xtensor<double, 1> y2 = xt::view(coor, xt::keep(n2), 1);
         xt::xtensor<double, 1> y3 = xt::view(coor, xt::keep(n3), 1);
         xt::xtensor<double, 2> ret = xt::empty<double>(conn.shape());
         xt::view(ret, xt::all(), 0) = xt::sqrt(xt::pow(x1 - x0, 2.0) + xt::pow(y1 - y0, 2.0));
         xt::view(ret, xt::all(), 1) = xt::sqrt(xt::pow(x2 - x1, 2.0) + xt::pow(y2 - y1, 2.0));
         xt::view(ret, xt::all(), 2) = xt::sqrt(xt::pow(x3 - x2, 2.0) + xt::pow(y3 - y2, 2.0));
         xt::view(ret, xt::all(), 3) = xt::sqrt(xt::pow(x0 - x3, 2.0) + xt::pow(y0 - y3, 2.0));
         return ret;
     }
 
     throw std::runtime_error("Element-type not implemented");
 }
 
 inline xt::xtensor<double, 2> edgesize(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn)
 {
     return edgesize(coor, conn, defaultElementType(coor, conn));
 }
 
 inline xt::xtensor<double, 2> centers(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn,
     ElementType type)
 {
     GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
     xt::xtensor<double, 2> ret = xt::zeros<double>({conn.shape(0), coor.shape(1)});
 
     if (type == ElementType::Quad4) {
         GOOSEFEM_ASSERT(coor.shape(1) == 2);
         GOOSEFEM_ASSERT(conn.shape(1) == 4);
         for (size_t i = 0; i < 4; ++i) {
             auto n = xt::view(conn, xt::all(), i);
             ret += xt::view(coor, xt::keep(n), xt::all());
         }
         ret /= 4.0;
         return ret;
     }
 
     throw std::runtime_error("Element-type not implemented");
 }
 
 inline xt::xtensor<double, 2> centers(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn)
 {
     return centers(coor, conn, defaultElementType(coor, conn));
 }
 
 inline xt::xtensor<size_t, 1> elemmap2nodemap(
     const xt::xtensor<size_t, 1>& elem_map,
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn,
     ElementType type)
 {
     GOOSEFEM_ASSERT(xt::amax(conn)() < coor.shape(0));
     GOOSEFEM_ASSERT(elem_map.size() == conn.shape(0));
     size_t N = coor.shape(0);
 
     xt::xtensor<size_t, 1> ret = N * xt::ones<size_t>({N});
 
     if (type == ElementType::Quad4) {
         GOOSEFEM_ASSERT(coor.shape(1) == 2);
         GOOSEFEM_ASSERT(conn.shape(1) == 4);
 
         for (size_t i = 0; i < 4; ++i) {
             xt::xtensor<size_t, 1> t = N * xt::ones<size_t>({N});
             auto old_nd = xt::view(conn, xt::all(), i);
             auto new_nd = xt::view(conn, xt::keep(elem_map), i);
             xt::view(t, xt::keep(old_nd)) = new_nd;
             ret = xt::where(xt::equal(ret, N), t, ret);
         }
 
         return ret;
     }
 
     throw std::runtime_error("Element-type not implemented");
 }
 
 inline xt::xtensor<size_t, 1> elemmap2nodemap(
     const xt::xtensor<size_t, 1>& elem_map,
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn)
 {
     return elemmap2nodemap(elem_map, coor, conn, defaultElementType(coor, conn));
 }
 
 } // namespace Mesh
 } // namespace GooseFEM
 
 #endif
diff --git a/python/Mesh.hpp b/python/Mesh.hpp
index b5ffeb4..7be2623 100644
--- a/python/Mesh.hpp
+++ b/python/Mesh.hpp
@@ -1,404 +1,434 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #include <GooseFEM/GooseFEM.h>
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
 
 namespace py = pybind11;
 
 void init_Mesh(py::module& m)
 {
 
     py::enum_<GooseFEM::Mesh::ElementType>(m,
             "ElementType",
             "ElementType."
             "See :cpp:enum:`GooseFEM::Mesh::ElementType`.")
         .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,
         "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")
 
         .def(py::init<
                 const xt::xtensor<double, 2>&,
                 const xt::xtensor<size_t, 2>&,
                 const xt::xtensor<size_t, 1>&,
                 const xt::xtensor<double, 2>&,
                 const xt::xtensor<size_t, 2>&,
                 const xt::xtensor<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",
-             &GooseFEM::Mesh::ManualStitch::nodemap,
-             "Map to new node-numbers."
+             py::overload_cast<>(
+                &GooseFEM::Mesh::ManualStitch::nodemap, py::const_),
+             "Node-map for givall sub-meshes."
              "See :cpp:func:`GooseFEM::Mesh::ManualStitch::nodemap`.")
 
         .def("elemmap",
-             &GooseFEM::Mesh::ManualStitch::elemmap,
-             "Map to new element-numbers."
+             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,
              "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,
              "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")
 
         .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,
              "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",
-             &GooseFEM::Mesh::Stitch::nodemap,
-             "Node-map for given mesh."
+             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",
-             &GooseFEM::Mesh::Stitch::elemmap,
-             "Element-map for given mesh."
+             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",
              py::overload_cast<const xt::xtensor<size_t, 1>&, size_t>(
                  &GooseFEM::Mesh::Stitch::nodeset, py::const_),
              "Convert node-set to the stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::nodeset`.",
              py::arg("set"),
              py::arg("mesh_index"))
 
         .def("elemset",
              py::overload_cast<const xt::xtensor<size_t, 1>&, size_t>(
                  &GooseFEM::Mesh::Stitch::elemset, py::const_),
              "Convert element-set to the stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::elemset`.",
              py::arg("set"),
              py::arg("mesh_index"))
 
         .def("nodeset",
              py::overload_cast<const std::vector<xt::xtensor<size_t, 1>>&>(
                  &GooseFEM::Mesh::Stitch::nodeset, py::const_),
              "Convert node-set to the stitched mesh."
              "See :cpp:func:`GooseFEM::Mesh::Stitch::nodeset`.",
              py::arg("set"))
 
         .def("elemset",
              py::overload_cast<const std::vector<xt::xtensor<size_t, 1>>&>(
                  &GooseFEM::Mesh::Stitch::elemset, py::const_),
              "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::Renumber>(m, "Renumber")
 
         .def(
             py::init<const xt::xarray<size_t>&>(),
             "Renumber to lowest possible index."
             "See :cpp:class:`GooseFEM::Mesh::Renumber`.",
             py::arg("dofs"))
 
         .def(
             "get",
             &GooseFEM::Mesh::Renumber::apply<xt::xarray<size_t>>,
             "Get renumbered DOFs."
             "See :cpp:func:`GooseFEM::Mesh::Renumber::get`.")
 
         .def(
             "apply",
             &GooseFEM::Mesh::Renumber::apply<xt::xarray<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")
 
         .def(
             py::init([](xt::xtensor<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::xtensor<size_t, 1>& a, xt::xtensor<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::xtensor<size_t, 1>& a, xt::xtensor<size_t, 1>& b, xt::xtensor<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::xtensor<size_t, 1>& a,
                         xt::xtensor<size_t, 1>& b,
                         xt::xtensor<size_t, 1>& c,
                         xt::xtensor<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::xarray<size_t>>,
             "Reorder matrix (e.g. ``dofs``)."
             "See :cpp:func:`GooseFEM::Mesh::Reorder::get`.",
             py::arg("dofs"))
 
         .def(
             "apply",
             &GooseFEM::Mesh::Reorder::apply<xt::xarray<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,
         "Renumber to lowest possible indices."
         "See :cpp:func:`GooseFEM::Mesh::renumber`.",
         py::arg("dofs"));
 
     m.def(
         "coordination",
         &GooseFEM::Mesh::coordination,
         "Coordination number of each node."
         "See :cpp:func:`GooseFEM::Mesh::coordination`.",
         py::arg("conn"));
 
     m.def(
         "elem2node",
         &GooseFEM::Mesh::elem2node,
         "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::xtensor<double, 2>&, const xt::xtensor<size_t, 2>&>(
             &GooseFEM::Mesh::edgesize),
         "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::xtensor<double, 2>&,
             const xt::xtensor<size_t, 2>&,
             GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::edgesize),
         "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::xtensor<double, 2>&, const xt::xtensor<size_t, 2>&>(
             &GooseFEM::Mesh::centers),
         "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::xtensor<double, 2>&,
             const xt::xtensor<size_t, 2>&,
             GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::centers),
         "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::xtensor<size_t, 1>&,
             const xt::xtensor<double, 2>&,
             const xt::xtensor<size_t, 2>&>(&GooseFEM::Mesh::elemmap2nodemap),
         "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::xtensor<size_t, 1>&,
             const xt::xtensor<double, 2>&,
             const xt::xtensor<size_t, 2>&,
             GooseFEM::Mesh::ElementType>(&GooseFEM::Mesh::elemmap2nodemap),
         "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"));
 }
diff --git a/test/basic/Mesh.cpp b/test/basic/Mesh.cpp
index f6ca057..5d1d69f 100644
--- a/test/basic/Mesh.cpp
+++ b/test/basic/Mesh.cpp
@@ -1,368 +1,372 @@
 
 #include <catch2/catch.hpp>
 #include <xtensor/xrandom.hpp>
 #include <xtensor/xmath.hpp>
 #include <GooseFEM/GooseFEM.h>
 
 #define ISCLOSE(a,b) REQUIRE_THAT((a), Catch::WithinAbs((b), 1.e-12));
 
 TEST_CASE("GooseFEM::Mesh", "Mesh.h")
 {
     SECTION("overlapping")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(5, 5, 1.0);
 
         auto coor_a = mesh.coor();
         auto overlap_a = mesh.nodesTopEdge();
 
         auto coor_b = mesh.coor();
         auto overlap_b = mesh.nodesBottomEdge();
         xt::view(coor_b, xt::all(), 1) += 5.0;
 
         auto overlap = GooseFEM::Mesh::overlapping(coor_a, coor_b);
 
         REQUIRE(xt::all(xt::equal(xt::view(overlap, 0, xt::all()), overlap_a)));
         REQUIRE(xt::all(xt::equal(xt::view(overlap, 1, xt::all()), overlap_b)));
     }
 
     SECTION("ManualStitch")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(5, 5, 1.0);
 
         auto coor_a = mesh.coor();
         auto conn_a = mesh.conn();
         auto overlap_a = mesh.nodesTopEdge();
 
         auto coor_b = mesh.coor();
         auto conn_b = mesh.conn();
         auto overlap_b = mesh.nodesBottomEdge();
         xt::view(coor_b, xt::all(), 1) += 5.0;
 
         GooseFEM::Mesh::ManualStitch stitch(coor_a, conn_a, overlap_a, coor_b, conn_b, overlap_b);
 
         GooseFEM::Mesh::Quad4::Regular res(5, 10, 1.0);
 
         REQUIRE(xt::allclose(stitch.coor(), res.coor()));
         REQUIRE(xt::all(xt::equal(stitch.conn(), res.conn())));
 
+        REQUIRE(stitch.nodemap().size() == 2);
+        REQUIRE(stitch.elemmap().size() == 2);
         REQUIRE(xt::all(xt::equal(stitch.nodemap(0), xt::arange<size_t>(6 * 6))));
         REQUIRE(xt::all(xt::equal(stitch.nodemap(1), xt::arange<size_t>(6 * 6) + 5 * 6)));
         REQUIRE(xt::all(xt::equal(stitch.elemmap(0), xt::arange<size_t>(5 * 5))));
         REQUIRE(xt::all(xt::equal(stitch.elemmap(1), xt::arange<size_t>(5 * 5) + 5 * 5)));
 
         REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::arange<size_t>(6 * 6), 0), xt::arange<size_t>(6 * 6))));
         REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::arange<size_t>(6 * 6), 1), xt::arange<size_t>(6 * 6) + 5 * 6)));
         REQUIRE(xt::all(xt::equal(stitch.elemset(xt::arange<size_t>(5 * 5), 0), xt::arange<size_t>(5 * 5))));
         REQUIRE(xt::all(xt::equal(stitch.elemset(xt::arange<size_t>(5 * 5), 1), xt::arange<size_t>(5 * 5) + 5 * 5)));
 
         REQUIRE(xt::all(xt::equal(stitch.nodeset(overlap_a, 0), stitch.nodeset(overlap_b, 1))));
     }
 
     SECTION("Stitch")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(5, 5, 1.0);
 
         auto coor_a = mesh.coor();
         auto conn_a = mesh.conn();
         auto overlap_a = mesh.nodesTopEdge();
         auto nset = mesh.nodesLeftEdge();
         auto eset = xt::eval(xt::arange<size_t>(mesh.nelem()));
 
         auto coor_b = mesh.coor();
         auto conn_b = mesh.conn();
         auto overlap_b = mesh.nodesBottomEdge();
         xt::view(coor_b, xt::all(), 1) += 5.0;
 
         GooseFEM::Mesh::Stitch stitch;
         stitch.push_back(coor_a, conn_a);
         stitch.push_back(coor_b, conn_b);
 
         GooseFEM::Mesh::Quad4::Regular res(5, 10, 1.0);
 
         REQUIRE(xt::allclose(stitch.coor(), res.coor()));
         REQUIRE(xt::all(xt::equal(stitch.conn(), res.conn())));
 
+        REQUIRE(stitch.nodemap().size() == 2);
+        REQUIRE(stitch.elemmap().size() == 2);
         REQUIRE(xt::all(xt::equal(stitch.nodemap(0), xt::arange<size_t>(6 * 6))));
         REQUIRE(xt::all(xt::equal(stitch.nodemap(1), xt::arange<size_t>(6 * 6) + 5 * 6)));
         REQUIRE(xt::all(xt::equal(stitch.elemmap(0), xt::arange<size_t>(5 * 5))));
         REQUIRE(xt::all(xt::equal(stitch.elemmap(1), xt::arange<size_t>(5 * 5) + 5 * 5)));
 
         REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::arange<size_t>(6 * 6), 0), xt::arange<size_t>(6 * 6))));
         REQUIRE(xt::all(xt::equal(stitch.nodeset(xt::arange<size_t>(6 * 6), 1), xt::arange<size_t>(6 * 6) + 5 * 6)));
         REQUIRE(xt::all(xt::equal(stitch.elemset(xt::arange<size_t>(5 * 5), 0), xt::arange<size_t>(5 * 5))));
         REQUIRE(xt::all(xt::equal(stitch.elemset(xt::arange<size_t>(5 * 5), 1), xt::arange<size_t>(5 * 5) + 5 * 5)));
 
         REQUIRE(xt::all(xt::equal(stitch.nodeset(overlap_a, 0), stitch.nodeset(overlap_b, 1))));
 
         REQUIRE(xt::all(xt::equal(stitch.nodeset({nset, nset}), xt::arange<size_t>(0, 6 * 6 + 5 * 6, 6))));
         REQUIRE(xt::all(xt::equal(stitch.elemset({eset, eset}), xt::arange<size_t>(2 * 5 * 5))));
     }
 
     SECTION("edgesize")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(2, 2, 10.0);
         auto s = GooseFEM::Mesh::edgesize(mesh.coor(), mesh.conn());
         auto t = GooseFEM::Mesh::edgesize(mesh.coor(), mesh.conn(), mesh.getElementType());
         REQUIRE(xt::allclose(s, 10.0));
         REQUIRE(xt::allclose(t, 10.0));
     }
 
     SECTION("coordination")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
         auto N = GooseFEM::Mesh::coordination(mesh.conn());
         xt::xtensor<size_t, 1> ret = {1, 2, 2, 1, 2, 4, 4, 2, 2, 4, 4, 2, 1, 2, 2, 1};
         REQUIRE(xt::all(xt::equal(N, ret)));
     }
 
     SECTION("centers")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(2, 2, 2.0);
         xt::xtensor<double, 2> c = {
             {1.0, 1.0},
             {3.0, 1.0},
             {1.0, 3.0},
             {3.0, 3.0}};
 
         REQUIRE(xt::allclose(GooseFEM::Mesh::centers(mesh.coor(), mesh.conn()), c));
     }
 
     SECTION("elem2node")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
         auto tonode = GooseFEM::Mesh::elem2node(mesh.conn());
         REQUIRE(tonode.size() == 16);
         REQUIRE(tonode[0] == std::vector<size_t>{0});
         REQUIRE(tonode[1] == std::vector<size_t>{0, 1});
         REQUIRE(tonode[2] == std::vector<size_t>{1, 2});
         REQUIRE(tonode[3] == std::vector<size_t>{2});
         REQUIRE(tonode[4] == std::vector<size_t>{0, 3});
         REQUIRE(tonode[5] == std::vector<size_t>{0, 1, 3, 4});
         REQUIRE(tonode[6] == std::vector<size_t>{1, 2, 4, 5});
         REQUIRE(tonode[7] == std::vector<size_t>{2, 5});
         REQUIRE(tonode[8] == std::vector<size_t>{3, 6});
         REQUIRE(tonode[9] == std::vector<size_t>{3, 4, 6, 7});
         REQUIRE(tonode[10] == std::vector<size_t>{4, 5, 7, 8});
         REQUIRE(tonode[11] == std::vector<size_t>{5, 8});
         REQUIRE(tonode[12] == std::vector<size_t>{6});
         REQUIRE(tonode[13] == std::vector<size_t>{6, 7});
         REQUIRE(tonode[14] == std::vector<size_t>{7, 8});
         REQUIRE(tonode[15] == std::vector<size_t>{8});
     }
 
     SECTION("elemmap2nodemap")
     {
         GooseFEM::Mesh::Quad4::Regular mesh(3, 3);
 
         xt::xtensor<size_t, 1> elmap0 = {
             0, 1, 2,
             3, 4, 5,
             6, 7, 8
         };
         xt::xtensor<size_t, 1> elmap1 = {
             2, 0, 1,
             5, 3, 4,
             8, 6, 7
         };
         xt::xtensor<size_t, 1> elmap2 = {
             1, 2, 0,
             4, 5, 3,
             7, 8, 6
         };
 
         xt::xtensor<size_t, 1> nodemap0 = {
              0,  1,  2,  3,
              4,  5,  6,  7,
              8,  9, 10, 11,
             12, 13, 14, 15
         };
         xt::xtensor<size_t, 1> nodemap1 = {
               2,  0,  1,  2,
               6,  4,  5,  6,
              10,  8,  9, 10,
              14, 15, 13, 14
         };
         xt::xtensor<size_t, 1> nodemap2 = {
              1,  2,  0,  1,
              5,  6,  4,  5,
              9, 10,  8,  9,
             13, 14, 15, 13
         };
 
         REQUIRE(xt::all(xt::equal(GooseFEM::Mesh::elemmap2nodemap(elmap0, mesh.coor(), mesh.conn()), nodemap0)));
         REQUIRE(xt::all(xt::equal(GooseFEM::Mesh::elemmap2nodemap(elmap1, mesh.coor(), mesh.conn()), nodemap1)));
         REQUIRE(xt::all(xt::equal(GooseFEM::Mesh::elemmap2nodemap(elmap2, mesh.coor(), mesh.conn()), nodemap2)));
     }
 
     SECTION("elemmap2nodemap - example 1")
     {
         GooseFEM::Mesh::Quad4::FineLayer mesh(3, 3);
 
         xt::xtensor<int, 1> elemval = {
             1, 0, 0,
             1, 0, 0,
             1, 0, 0
         };
 
         xt::xtensor<int, 1> elemval_r1 = {
             0, 1, 0,
             0, 1, 0,
             0, 1, 0
         };
 
         xt::xtensor<int, 1> elemval_r2 = {
             0, 0, 1,
             0, 0, 1,
             0, 0, 1
         };
 
         xt::xtensor<int, 1> nodeval = {
             1, 0, 0, 1,
             1, 0, 0, 1,
             1, 0, 0, 1,
             1, 0, 0, 1
         };
 
         xt::xtensor<int, 1> nodeval_r1 = {
             0, 1, 0, 0,
             0, 1, 0, 0,
             0, 1, 0, 0,
             0, 1, 0, 0
         };
 
         xt::xtensor<int, 1> nodeval_r2 = {
             0, 0, 1, 0,
             0, 0, 1, 0,
             0, 0, 1, 0,
             0, 0, 1, 0
         };
 
         {
             auto elemmap = mesh.roll(0);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap)))));
         }
 
         {
             auto elemmap = mesh.roll(1);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap)))));
         }
 
         {
             auto elemmap = mesh.roll(2);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap)))));
         }
 
         {
             auto elemmap = mesh.roll(3);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap)))));
         }
 
         {
             auto elemmap = mesh.roll(4);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap)))));
         }
 
         {
             auto elemmap = mesh.roll(5);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap)))));
         }
     }
 
     SECTION("elemmap2nodemap - example 2")
     {
         GooseFEM::Mesh::Quad4::FineLayer mesh(3, 3);
 
         xt::xtensor<int, 1> elemval = {
             1, 0, 0,
             0, 1, 0,
             0, 0, 1
         };
 
         xt::xtensor<int, 1> elemval_r1 = {
             0, 1, 0,
             0, 0, 1,
             1, 0, 0
         };
 
         xt::xtensor<int, 1> elemval_r2 = {
             0, 0, 1,
             1, 0, 0,
             0, 1, 0
         };
 
         xt::xtensor<int, 1> nodeval = {
             1, 0, 0, 1,
             0, 1, 0, 0,
             0, 0, 1, 0,
             1, 0, 0, 1
         };
 
         xt::xtensor<int, 1> nodeval_r1 = {
             0, 1, 0, 0,
             0, 0, 1, 0,
             1, 0, 0, 1,
             0, 1, 0, 0
         };
 
         xt::xtensor<int, 1> nodeval_r2 = {
             0, 0, 1, 0,
             1, 0, 0, 1,
             0, 1, 0, 0,
             0, 0, 1, 0
         };
 
         {
             auto elemmap = mesh.roll(0);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap)))));
         }
 
         {
             auto elemmap = mesh.roll(1);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap)))));
         }
 
         {
             auto elemmap = mesh.roll(2);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap)))));
         }
 
         {
             auto elemmap = mesh.roll(3);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval, xt::view(nodeval, xt::keep(nodemap)))));
         }
 
         {
             auto elemmap = mesh.roll(4);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval_r1, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval_r1, xt::view(nodeval, xt::keep(nodemap)))));
         }
 
         {
             auto elemmap = mesh.roll(5);
             auto nodemap = GooseFEM::Mesh::elemmap2nodemap(elemmap, mesh.coor(), mesh.conn());
             REQUIRE(xt::all(xt::equal(elemval_r2, xt::view(elemval, xt::keep(elemmap)))));
             REQUIRE(xt::all(xt::equal(nodeval_r2, xt::view(nodeval, xt::keep(nodemap)))));
         }
     }
 }