diff --git a/include/GooseFEM/Mesh.h b/include/GooseFEM/Mesh.h
index 14b21bc..878b79a 100644
--- a/include/GooseFEM/Mesh.h
+++ b/include/GooseFEM/Mesh.h
@@ -1,88 +1,99 @@
 /*
 
 (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 {
 
 enum class ElementType {
     Quad4,
     Hex8,
     Tri3 };
 
 // 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);
 
+// 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);
+
+// return size of each element edge: extract element-type based on shape of "conn"
+inline xt::xtensor<double,2> edgesize(
+    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 d2612e5..3d33899 100644
--- a/include/GooseFEM/Mesh.hpp
+++ b/include/GooseFEM/Mesh.hpp
@@ -1,156 +1,206 @@
 /*
 
 (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 Renumber::Renumber(const xt::xarray<size_t>& dofs)
 {
     size_t n = xt::amax(dofs)[0] + 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;
     }
 }
 
 // out(i,j) = renum(list(i,j))
 template <class T>
 T Renumber::apply(const T& list) const
 {
     T out = T::from_shape(list.shape());
 
     auto jt = out.begin();
 
     for (auto it = list.begin(); it != list.end(); ++it, ++jt) {
         *jt = m_renum(*it);
     }
 
     return out;
 }
 
 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)[0] + 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:
 //
 //   out(i,j) = renum(list(i,j))
 
 template <class T>
 T Reorder::apply(const T& list) const
 {
     T out = T::from_shape(list.shape());
 
     auto jt = out.begin();
 
     for (auto it = list.begin(); it != list.end(); ++it, ++jt) {
         *jt = m_renum(*it);
     }
 
     return out;
 }
 
 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)[0] + 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)
 {
     auto N = coordination(conn);
     auto nnode = N.size();
 
     std::vector<std::vector<size_t>> out;
     out.resize(nnode);
     for (size_t i = 0; i < nnode; ++i) {
         out[i].reserve(N(i));
     }
 
     for (size_t e = 0; e < conn.shape(0); ++e) {
         for (size_t m = 0; m < conn.shape(1); ++m) {
             out[conn(e, m)].push_back(e);
         }
     }
 
     return out;
 }
 
+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> out = xt::empty<double>(conn.shape());
+        xt::view(out, xt::all(), 0) = xt::sqrt(xt::pow(x1 - x0, 2.0) + xt::pow(y1 - y0, 2.0));
+        xt::view(out, xt::all(), 1) = xt::sqrt(xt::pow(x2 - x1, 2.0) + xt::pow(y2 - y1, 2.0));
+        xt::view(out, xt::all(), 2) = xt::sqrt(xt::pow(x3 - x2, 2.0) + xt::pow(y3 - y2, 2.0));
+        xt::view(out, xt::all(), 3) = xt::sqrt(xt::pow(x0 - x3, 2.0) + xt::pow(y0 - y3, 2.0));
+        return out;
+    }
+
+    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)
+{
+    if (coor.shape(1) == 2ul && conn.shape(1) == 3ul) {
+        return edgesize(coor, conn, ElementType::Tri3);
+    }
+    if (coor.shape(1) == 2ul && conn.shape(1) == 4ul) {
+        return edgesize(coor, conn, ElementType::Quad4);
+    }
+    if (coor.shape(1) == 3ul && conn.shape(1) == 8ul) {
+        return edgesize(coor, conn, ElementType::Hex8);
+    }
+
+    throw std::runtime_error("Element-type not implemented");
+}
+
 } // namespace Mesh
 } // namespace GooseFEM
 
 #endif
diff --git a/include/GooseFEM/config.h b/include/GooseFEM/config.h
index 86a77ea..41e9741 100644
--- a/include/GooseFEM/config.h
+++ b/include/GooseFEM/config.h
@@ -1,76 +1,76 @@
 /*
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 */
 
 #ifndef GOOSEFEM_CONFIG_H
 #define GOOSEFEM_CONFIG_H
 
 #define _USE_MATH_DEFINES // to use "M_PI" from "math.h"
 
 #include <assert.h>
 #include <cstdlib>
 #include <vector>
 #include <string>
 #include <memory>
 #include <iostream>
 #include <iomanip>
 #include <numeric>
 #include <limits>
 #include <algorithm>
 #include <math.h>
 
 #include <xtensor/xarray.hpp>
 #include <xtensor/xtensor.hpp>
 #include <xtensor/xadapt.hpp>
 #include <xtensor/xfixed.hpp>
 #include <xtensor/xinfo.hpp>
 #include <xtensor/xio.hpp>
 #include <xtensor/xlayout.hpp>
 #include <xtensor/xmath.hpp>
 #include <xtensor/xnoalias.hpp>
 #include <xtensor/xsort.hpp>
 #include <xtensor/xstrided_view.hpp>
 #include <xtensor/xutils.hpp>
 #include <xtensor/xview.hpp>
 
 using namespace xt::placeholders;
 
 #define UNUSED(p) ((void)(p))
 
 #ifdef GOOSEFEM_ENABLE_ASSERT
     #define GOOSEFEM_ASSERT(expr) GOOSEFEM_ASSERT_IMPL(expr, __FILE__, __LINE__)
     #define GOOSEFEM_ASSERT_IMPL(expr, file, line) \
         if (!(expr)) { \
             throw std::runtime_error( \
                 std::string(file) + ':' + std::to_string(line) + \
                 ": assertion failed (" #expr ") \n\t"); \
         }
 #else
     #define GOOSEFEM_ASSERT(expr)
 #endif
 
 #define GOOSEFEM_CHECK(expr) GOOSEFEM_CHECK_IMPL(expr, __FILE__, __LINE__)
 #define GOOSEFEM_CHECK_IMPL(expr, file, line) \
     if (!(expr)) { \
         throw std::runtime_error( \
             std::string(file) + ':' + std::to_string(line) + \
             ": assertion failed (" #expr ") \n\t"); \
     }
 
 #define GOOSEFEM_VERSION_MAJOR 0
 #define GOOSEFEM_VERSION_MINOR 3
-#define GOOSEFEM_VERSION_PATCH 0
+#define GOOSEFEM_VERSION_PATCH 1
 
 #define GOOSEFEM_VERSION_AT_LEAST(x, y, z) \
     (GOOSEFEM_VERSION_MAJOR > x || (GOOSEFEM_VERSION_MAJOR >= x && \
     (GOOSEFEM_VERSION_MINOR > y || (GOOSEFEM_VERSION_MINOR >= y && \
                                     GOOSEFEM_VERSION_PATCH >= z))))
 
 #define GOOSEFEM_VERSION(x, y, z) \
     (GOOSEFEM_VERSION_MAJOR == x && \
      GOOSEFEM_VERSION_MINOR == y && \
      GOOSEFEM_VERSION_PATCH == z)
 
 #endif
diff --git a/python/Mesh.hpp b/python/Mesh.hpp
index eae3d68..0e2d9c5 100644
--- a/python/Mesh.hpp
+++ b/python/Mesh.hpp
@@ -1,100 +1,117 @@
 /* =================================================================================================
 
 (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
 
 ================================================================================================= */
 
 #include <pybind11/pybind11.h>
 #include <pyxtensor/pyxtensor.hpp>
 #include "../include/GooseFEM/GooseFEM.h"
 
 // =================================================================================================
 
 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::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"));
 
+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"));
+
 }
 
 // =================================================================================================
 
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 9d3ae8b..3e02f10 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,42 +1,43 @@
 
 cmake_minimum_required(VERSION 3.0)
 
 if(CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR)
     project(GooseFEM-test)
     find_package(GooseFEM REQUIRED CONFIG)
 endif()
 
 option(WARNINGS "Show build warnings" ON)
 option(SIMD "Enable xsimd" ON)
 option(ASSERT "Enable assertions" ON)
 option(DEBUG "Enable all assertions" ON)
 
 find_package(Catch2 REQUIRED)
 
 add_executable(main
     main.cpp
     ElementHex8.cpp
     ElementQuad4.cpp
     Iterate.cpp
     MatrixDiagonal.cpp
+    Mesh.cpp
     MeshQuad4.cpp
     Vector.cpp
     VectorPartitioned.cpp)
 
 target_link_libraries(main PRIVATE Catch2::Catch2 GooseFEM)
 
 if (SIMD)
     target_link_libraries(main PRIVATE xtensor::optimize xtensor::use_xsimd)
 endif()
 
 if (WARNINGS)
     target_link_libraries(main PRIVATE GooseFEM::compiler_warnings)
 endif()
 
 if (ASSERT)
     target_link_libraries(main PRIVATE GooseFEM::assert)
 endif()
 
 if (DEBUG)
     target_link_libraries(main PRIVATE GooseFEM::debug)
 endif()
diff --git a/test/Mesh.cpp b/test/Mesh.cpp
new file mode 100644
index 0000000..a922289
--- /dev/null
+++ b/test/Mesh.cpp
@@ -0,0 +1,15 @@
+#include "support.h"
+
+TEST_CASE("GooseFEM::Mesh", "Mesh.h")
+{
+
+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));
+}
+
+}