diff --git a/include/GooseFEM/Mesh.h b/include/GooseFEM/Mesh.h
index 65cf851..5f99eff 100644
--- a/include/GooseFEM/Mesh.h
+++ b/include/GooseFEM/Mesh.h
@@ -1,98 +1,114 @@
 /*
 
 (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 };
 
+inline ElementType defaultElementType(
+    const xt::xtensor<double, 2>& coor,
+    const xt::xtensor<size_t, 2>& conn);
+
 // 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);
+    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);
+    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"
 
 } // namespace Mesh
 } // namespace GooseFEM
 
 #include "Mesh.hpp"
 
 #endif
diff --git a/include/GooseFEM/Mesh.hpp b/include/GooseFEM/Mesh.hpp
index f8eaf77..43c516f 100644
--- a/include/GooseFEM/Mesh.hpp
+++ b/include/GooseFEM/Mesh.hpp
@@ -1,209 +1,247 @@
 /*
 
 (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");
+}
+
 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
 {
     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::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)
+    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);
+    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));
+}
+
+
 } // namespace Mesh
 } // namespace GooseFEM
 
 #endif
diff --git a/test/basic/Mesh.cpp b/test/basic/Mesh.cpp
index 42121f6..b6f126b 100644
--- a/test/basic/Mesh.cpp
+++ b/test/basic/Mesh.cpp
@@ -1,50 +1,62 @@
 
 #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("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});
     }
 }