diff --git a/docs/details/figures/MeshQuad4/FineLayer/stitch.py b/docs/details/figures/MeshQuad4/FineLayer/stitch.py
index ce3dd3c..f3297f2 100644
--- a/docs/details/figures/MeshQuad4/FineLayer/stitch.py
+++ b/docs/details/figures/MeshQuad4/FineLayer/stitch.py
@@ -1,44 +1,44 @@
 
 import GooseFEM
 import matplotlib.pyplot as plt
 import GooseMPL as gplt
 import numpy as np
 
 plt.style.use(['goose'])
 
 fig, ax = plt.subplots()
 
 mesh = GooseFEM.Mesh.Quad4.FineLayer(6 * 9, 51)
 
 coor_a = mesh.coor()
 conn_a = mesh.conn()
 
 coor_b = mesh.coor()
 conn_b = mesh.conn()
 
 coor_c = mesh.coor()
 conn_c = mesh.conn()
 
 coor_b[:, 1] += np.max(coor_a[:, 1])
 coor_c[:, 1] += np.max(coor_a[:, 1]) * 2
 
 stitch = GooseFEM.Mesh.Stitch()
 stitch.push_back(coor_a, conn_a)
 stitch.push_back(coor_b, conn_b)
 stitch.push_back(coor_c, conn_c)
 
 coor = stitch.coor()
 conn = stitch.conn()
 
 cindex = np.zeros(conn.shape[0])
-cindex[stitch.elementset(np.arange(mesh.nelem()), 0)] = 1
-cindex[stitch.elementset(np.arange(mesh.nelem()), 1)] = 2
-cindex[stitch.elementset(np.arange(mesh.nelem()), 2)] = 3
+cindex[stitch.elemset(np.arange(mesh.nelem()), 0)] = 1
+cindex[stitch.elemset(np.arange(mesh.nelem()), 1)] = 2
+cindex[stitch.elemset(np.arange(mesh.nelem()), 2)] = 3
 
 im = gplt.patch(coor=coor, conn=conn, cindex=cindex, cmap='jet', axis=ax)
 
 ax.set_aspect('equal')
 ax.get_xaxis().set_visible(False)
 ax.get_yaxis().set_visible(False)
 
 plt.show()
diff --git a/include/GooseFEM/Mesh.h b/include/GooseFEM/Mesh.h
index 1bc73b1..2213a2f 100644
--- a/include/GooseFEM/Mesh.h
+++ b/include/GooseFEM/Mesh.h
@@ -1,206 +1,209 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #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
 
 inline ElementType defaultElementType(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn);
 
 // Stitch two mesh objects, specifying overlapping nodes by hand
 
 class ManualStitch {
 public:
     ManualStitch() = default;
 
     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);
 
     xt::xtensor<double, 2> coor() const;
     xt::xtensor<size_t, 2> conn() const;
     xt::xtensor<size_t, 1> nodemap(size_t mesh_index) const;
+    xt::xtensor<size_t, 1> elemmap(size_t mesh_index) const;
 
     // convert set of of node/element-numbers for an original mesh to the final mesh
     xt::xtensor<size_t, 1> nodeset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const;
-    xt::xtensor<size_t, 1> elementset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const;
+    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, searching for overlapping nodes
 
 class Stitch {
 public:
     Stitch() = default;
 
     Stitch(double rtol, double atol);
 
     void push_back(const xt::xtensor<double, 2>& coor, const xt::xtensor<size_t, 2>& conn);
 
-    // return connectivity
     xt::xtensor<double, 2> coor() const;
     xt::xtensor<size_t, 2> conn() const;
     xt::xtensor<size_t, 1> nodemap(size_t mesh_index) const;
+    xt::xtensor<size_t, 1> elemmap(size_t mesh_index) const;
 
     // convert set of node/element-numbers for an original mesh to the final mesh
     xt::xtensor<size_t, 1> nodeset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const;
-    xt::xtensor<size_t, 1> elementset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const;
+    xt::xtensor<size_t, 1> elemset(const xt::xtensor<size_t, 1>& set, size_t mesh_index) const;
 
     // combine set of node/element-numbers for an original to the final mesh (removes duplicates)
     xt::xtensor<size_t, 1> nodeset(const std::vector<xt::xtensor<size_t, 1>>& set) const;
-    xt::xtensor<size_t, 1> elementset(const std::vector<xt::xtensor<size_t, 1>>& set) const;
+    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;
+    std::vector<size_t> m_el_offset;
     double m_rtol = 1e-5;
     double m_atol = 1e-8;
 };
 
 // Renumber to lowest possible index. For example [0,3,4,2] -> [0,2,3,1]
 
 class Renumber {
 public:
     // constructors
     Renumber() = default;
     Renumber(const xt::xarray<size_t>& dofs);
 
     // get renumbered DOFs (same as "Renumber::apply(dofs)")
     xt::xtensor<size_t, 2> get(const xt::xtensor<size_t, 2>& dofs) const;
 
     // apply renumbering to other set
     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))
     xt::xtensor<size_t, 1> index() const;
 
 private:
     xt::xtensor<size_t, 1> m_renum;
 };
 
 // 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:
     // constructors
     Reorder() = default;
     Reorder(const std::initializer_list<xt::xtensor<size_t, 1>> args);
 
     // get reordered DOFs (same as "Reorder::apply(dofs)")
     xt::xtensor<size_t, 2> get(const xt::xtensor<size_t, 2>& dofs) const;
 
     // apply renumbering to other set
     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))
     xt::xtensor<size_t, 1> index() const;
 
 private:
     xt::xtensor<size_t, 1> m_renum;
 };
 
 // list with DOF-numbers in sequential order
 inline xt::xtensor<size_t, 2> dofs(size_t nnode, size_t ndim);
 
 // renumber to lowest possible index (see "GooseFEM::Mesh::Renumber")
 inline xt::xtensor<size_t, 2> renumber(const xt::xtensor<size_t, 2>& dofs);
 
 // number of elements connected to each node
 inline xt::xtensor<size_t, 1> coordination(const xt::xtensor<size_t, 2>& conn);
 
 // elements connected to each node
 inline std::vector<std::vector<size_t>> elem2node(
     const xt::xtensor<size_t, 2>& conn,
     bool sorted=true); // ensure the output to be sorted
 
 // return size of each element edge
 inline xt::xtensor<double, 2> edgesize(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn,
     ElementType type);
 
 inline xt::xtensor<double, 2> edgesize(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn); // extract element-type based on shape of "conn"
 
 // Coordinates of the 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);
 
 inline xt::xtensor<double, 2> centers(
     const xt::xtensor<double, 2>& coor,
     const xt::xtensor<size_t, 2>& conn); // extract element-type based on shape of "conn"
 
 // Convert an element-map to a node-map.
 // Input: new_elvar = elvar[elem_map]
 // Return: 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); // extract element-type based on shape of "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);
 
 // Find overlapping nodes
 // Return: [[nodes_from_mesh_a],
 //          [nodes_from_mesh_b]]
 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);
 
 } // namespace Mesh
 } // namespace GooseFEM
 
 #include "Mesh.hpp"
 
 #endif
diff --git a/include/GooseFEM/Mesh.hpp b/include/GooseFEM/Mesh.hpp
index c4ba2fb..cc97568 100644
--- a/include/GooseFEM/Mesh.hpp
+++ b/include/GooseFEM/Mesh.hpp
@@ -1,527 +1,542 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #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 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 xt::xtensor<size_t, 1> ManualStitch::nodemap(size_t index) const
 {
     GOOSEFEM_ASSERT(index <= 1);
 
     if (index == 0) {
         return xt::arange<size_t>(m_nnd_a);
     }
 
     return m_map_b;
 }
 
+inline xt::xtensor<size_t, 1> ManualStitch::elemmap(size_t index) const
+{
+    GOOSEFEM_ASSERT(index <= 1);
+
+    if (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 index) const
 {
     GOOSEFEM_ASSERT(index <= 1);
 
     if (index == 0) {
         return set;
     }
 
     return detail::renum(set, m_map_b);
 }
 
 inline xt::xtensor<size_t, 1>
-ManualStitch::elementset(const xt::xtensor<size_t, 1>& set, size_t index) const
+ManualStitch::elemset(const xt::xtensor<size_t, 1>& set, size_t index) const
 {
     GOOSEFEM_ASSERT(index <= 1);
 
     if (index == 0) {
         return set;
     }
 
     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] + conn.shape(0));
 }
 
 inline xt::xtensor<double, 2> Stitch::coor() const
 {
     return m_coor;
 }
 
 inline xt::xtensor<size_t, 2> Stitch::conn() const
 {
     return m_conn;
 }
 
 inline xt::xtensor<size_t, 1> Stitch::nodemap(size_t index) const
 {
     GOOSEFEM_ASSERT(index < m_map.size());
     return m_map[index];
 }
 
-inline xt::xtensor<size_t, 1>
-Stitch::nodeset(const xt::xtensor<size_t, 1>& set, size_t index) const
+inline xt::xtensor<size_t, 1> Stitch::elemmap(size_t index) const
+{
+    GOOSEFEM_ASSERT(index < m_map.size());
+    return xt::arange<size_t>(m_nel[index]) + m_el_offset[index];
+}
+
+inline xt::xtensor<size_t, 1> Stitch::nodeset(const xt::xtensor<size_t, 1>& set, size_t index) const
 {
     GOOSEFEM_ASSERT(index < m_map.size());
     return detail::renum(set, m_map[index]);
 }
 
-inline xt::xtensor<size_t, 1>
-Stitch::elementset(const xt::xtensor<size_t, 1>& set, size_t index) const
+inline xt::xtensor<size_t, 1> Stitch::elemset(const xt::xtensor<size_t, 1>& set, size_t index) const
 {
     GOOSEFEM_ASSERT(index < m_map.size());
-    size_t offset = 0;
-    for (size_t i = 0; i < index; ++i) {
-        offset += m_nel[i];
-    }
-    return set + offset;
+    return set + m_el_offset[index];
 }
 
-xt::xtensor<size_t, 1> Stitch::nodeset(const std::vector<xt::xtensor<size_t, 1>>& set) const
+inline xt::xtensor<size_t, 1> Stitch::nodeset(const std::vector<xt::xtensor<size_t, 1>>& set) const
 {
     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);
 }
 
-xt::xtensor<size_t, 1> Stitch::elementset(const std::vector<xt::xtensor<size_t, 1>>& set) const
+inline xt::xtensor<size_t, 1> Stitch::elemset(const std::vector<xt::xtensor<size_t, 1>>& set) const
 {
     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->elementset(set[i], i);
+        xt::view(ret, xt::range(n, n + set[i].size())) = this->elemset(set[i], i);
         n += set[i].size();
     }
 
     return ret;
 }
 
 inline Renumber::Renumber(const xt::xarray<size_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;
     }
 }
 
 // ret(i,j) = renum(list(i,j))
 template <class T>
 T Renumber::apply(const T& list) const
 {
     return detail::renum(list, m_renum);
 }
 
 inline xt::xtensor<size_t, 2> Renumber::get(const xt::xtensor<size_t, 2>& dofs) const
 {
     return this->apply(dofs);
 }
 
 inline xt::xtensor<size_t, 1> Renumber::index() const
 {
     return m_renum;
 }
 
 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
 {
     return this->apply(dofs);
 }
 
 inline xt::xtensor<size_t, 1> Reorder::index() const
 {
     return m_renum;
 }
 
 // apply renumbering, e.g. for a matrix:
 //
 //   ret(i,j) = renum(list(i,j))
 
 template <class T>
 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, 2> renumber(const xt::xtensor<size_t, 2>& dofs)
 {
     return Renumber(dofs).get(dofs);
 }
 
 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));
 }
 
 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;
 }
 
 } // namespace Mesh
 } // namespace GooseFEM
 
 #endif
diff --git a/python/Mesh.hpp b/python/Mesh.hpp
index 12d0e41..30841c5 100644
--- a/python/Mesh.hpp
+++ b/python/Mesh.hpp
@@ -1,264 +1,266 @@
 /* =================================================================================================
 
 (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")
         .value("Tri3", GooseFEM::Mesh::ElementType::Tri3)
         .value("Quad4", GooseFEM::Mesh::ElementType::Quad4)
         .value("Hex8", GooseFEM::Mesh::ElementType::Hex8)
         .export_values();
 
     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>(),
             "ManualStitch meshes",
             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("coor", &GooseFEM::Mesh::ManualStitch::coor, "Return coordinates of stitched mesh")
         .def("conn", &GooseFEM::Mesh::ManualStitch::conn, "Return connectivity of stitched mesh")
         .def("nodemap", &GooseFEM::Mesh::ManualStitch::nodemap, "Map to new node-numbers")
+        .def("elemmap", &GooseFEM::Mesh::ManualStitch::elemmap, "Map to new element-numbers")
 
         .def(
             "nodeset",
             &GooseFEM::Mesh::ManualStitch::nodeset,
             "Transfer nodeset to the stitched mesh",
-            py::arg("arg"),
+            py::arg("set"),
             py::arg("mesh_index"))
 
         .def(
-            "elementset",
-            &GooseFEM::Mesh::ManualStitch::elementset,
+            "elemset",
+            &GooseFEM::Mesh::ManualStitch::elemset,
             "Transfer elementset to the stitched mesh",
-            py::arg("arg"),
+            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",
             py::arg("rtol") = 1e-5,
             py::arg("atol") = 1e-8)
 
         .def(
             "push_back",
             &GooseFEM::Mesh::Stitch::push_back,
             "Add mesh",
             py::arg("coor"),
             py::arg("conn"))
 
         .def("coor", &GooseFEM::Mesh::Stitch::coor, "Return coordinates of stitched mesh")
         .def("conn", &GooseFEM::Mesh::Stitch::conn, "Return connectivity of stitched mesh")
         .def("nodemap", &GooseFEM::Mesh::Stitch::nodemap, "Map to new node-numbers")
+        .def("elemmap", &GooseFEM::Mesh::Stitch::elemmap, "Map to new element-numbers")
 
         .def(
             "nodeset",
             py::overload_cast<const xt::xtensor<size_t, 1>&, size_t>(
                 &GooseFEM::Mesh::Stitch::nodeset, py::const_),
             "Transfer nodeset to the stitched mesh",
-            py::arg("arg"),
+            py::arg("set"),
             py::arg("mesh_index"))
 
         .def(
-            "elementset",
+            "elemset",
             py::overload_cast<const xt::xtensor<size_t, 1>&, size_t>(
-                &GooseFEM::Mesh::Stitch::elementset, py::const_),
+                &GooseFEM::Mesh::Stitch::elemset, py::const_),
             "Transfer elementset to the stitched mesh",
-            py::arg("arg"),
+            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_),
             "Transfer nodeset to the stitched mesh",
-            py::arg("arg"))
+            py::arg("set"))
 
         .def(
-            "elementset",
+            "elemset",
             py::overload_cast<const std::vector<xt::xtensor<size_t, 1>>&>(
-                &GooseFEM::Mesh::Stitch::elementset, py::const_),
+                &GooseFEM::Mesh::Stitch::elemset, py::const_),
             "Transfer elementset to the stitched mesh",
-            py::arg("arg"))
+            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>&>(),
             "Class to renumber to the lowest possible indices. Use ``Renumber(...).get()`` to get "
             "the renumbered result",
             py::arg("dofs"))
 
         .def("get", &GooseFEM::Mesh::Renumber::get, "Get result of renumbering")
 
         .def(
             "index",
             &GooseFEM::Mesh::Renumber::index,
             "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``")
 
         .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}); }))
 
         .def(py::init([](xt::xtensor<size_t, 1>& a, xt::xtensor<size_t, 1>& b) {
             return new GooseFEM::Mesh::Reorder({a, b});
         }))
 
         .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});
             }))
 
         .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});
         }))
 
         .def("get", &GooseFEM::Mesh::Reorder::get, "Reorder matrix (e.g. ``dofs``)")
 
         .def(
             "index",
             &GooseFEM::Mesh::Reorder::index,
             "Get index list to apply renumbering. Apply renumbering using ``index[dofs]``")
 
         .def("__repr__", [](const GooseFEM::Mesh::Reorder&) { return "<GooseFEM.Mesh.Reorder>"; });
 
     m.def(
         "dofs",
         &GooseFEM::Mesh::dofs,
         "List with DOF-numbers (in sequential order)",
         py::arg("nnode"),
         py::arg("ndim"));
 
     m.def(
         "renumber",
         &GooseFEM::Mesh::renumber,
         "Renumber to lowest possible indices. Use ``GooseFEM.Mesh.Renumber`` for advanced "
         "functionality",
         py::arg("dofs"));
 
     m.def(
         "coordination",
         &GooseFEM::Mesh::coordination,
         "Coordination number of each node (number of elements connected to each node)",
         py::arg("conn"));
 
     m.def(
         "elem2node",
         &GooseFEM::Mesh::elem2node,
         "Element-numbers connected to each node",
         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",
         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",
         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",
         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",
         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",
         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",
         py::arg("elem_map"),
         py::arg("coor"),
         py::arg("conn"),
         py::arg("type"));
 
     m.def(
         "overlapping",
         &GooseFEM::Mesh::overlapping,
         "Find overlapping nodes",
         py::arg("coor_a"),
         py::arg("coor_b"),
         py::arg("rtol") = 1e-5,
         py::arg("atol") = 1e-8);
 }
diff --git a/test/basic/Mesh.cpp b/test/basic/Mesh.cpp
index f3e1ee6..f6ca057 100644
--- a/test/basic/Mesh.cpp
+++ b/test/basic/Mesh.cpp
@@ -1,345 +1,368 @@
 
 #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(res.coor(), stitch.coor()));
-        REQUIRE(xt::all(xt::equal(res.conn(), stitch.conn())));
+        REQUIRE(xt::allclose(stitch.coor(), res.coor()));
+        REQUIRE(xt::all(xt::equal(stitch.conn(), res.conn())));
+
+        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(res.coor(), stitch.coor()));
-        REQUIRE(xt::all(xt::equal(res.conn(), stitch.conn())));
+        REQUIRE(xt::allclose(stitch.coor(), res.coor()));
+        REQUIRE(xt::all(xt::equal(stitch.conn(), res.conn())));
+
+        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.elementset({eset, eset}), xt::arange<size_t>(2 * 5 * 5))));
+        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)))));
         }
     }
 }