diff --git a/cpp/LabelledTetrahedralMesh.cpp b/cpp/LabelledTetrahedralMesh.cpp
index e7ab4d2ec7cfaf467f77615084409c63a22a8ba2..a4655c76348005ee5ae4b41eb1bbea673a8f72a9 100644
--- a/cpp/LabelledTetrahedralMesh.cpp
+++ b/cpp/LabelledTetrahedralMesh.cpp
@@ -17,3 +17,30 @@ void LabelledTetrahedralMesh::addPoint(Eigen::Array3d point)
     std::vector<double> _point = std::vector<double> {point(0, 0), point(0, 1), point(0, 2)};
     this->_vertex_points.push_back(_point);
 }
+
+Eigen::Array<unsigned short, Eigen::Dynamic, 4> LabelledTetrahedralMesh::tetraVertices(void) const
+{
+    Eigen::Array<unsigned short, Eigen::Dynamic, 4> tetra_vertices(this->_tetra_vertices.size(), 4);
+    for (int i=0; i<this->_tetra_vertices.size(); i++) {
+        for (int v=0; v<4; v++) {
+            tetra_vertices(i, v) = this->_tetra_vertices[i][v];
+        }
+    }
+    return tetra_vertices;
+}
+
+Eigen::Array<unsigned short, Eigen::Dynamic, 1> LabelledTetrahedralMesh::tetraLabels(void) const
+{
+    Eigen::Array<unsigned short, Eigen::Dynamic, 1> tetra_labels(this->_tetra_labels.size(), 1);
+    for (int i=0; i<this->_tetra_labels.size(); i++) {
+        tetra_labels(i, 0) = this->_tetra_labels[i];
+    }
+    return tetra_labels;
+}
+
+void LabelledTetrahedralMesh::addTetra(Eigen::Array<unsigned short, 4, 1>  tetra, unsigned short label)
+{
+    std::vector<unsigned short> _tetra = std::vector<unsigned short> {tetra(0, 0), tetra(0, 1), tetra(0, 2), tetra(0, 3)};
+    this->_tetra_vertices.push_back(_tetra);
+    this->_tetra_labels.push_back(label);
+}
diff --git a/cpp/LabelledTetrahedralMesh.h b/cpp/LabelledTetrahedralMesh.h
index c6be404731067ed8f319b01d4d270bad909b7c33..1305737338c21d832c86bb78a4da5e4cf04608ba 100644
--- a/cpp/LabelledTetrahedralMesh.h
+++ b/cpp/LabelledTetrahedralMesh.h
@@ -13,6 +13,12 @@ public:
     Eigen::ArrayX3d vertexPoints(void) const;
     void addPoint(Eigen::Array3d point);
 
+    Eigen::Array<unsigned short, Eigen::Dynamic, 4> tetraVertices(void) const;
+    Eigen::Array<unsigned short, Eigen::Dynamic, 1> tetraLabels(void) const;
+    void addTetra(Eigen::Array<unsigned short, 4, 1>  tetra, unsigned short label);
+
 private:
     std::vector< std::vector<double> > _vertex_points;
+    std::vector< std::vector<unsigned short> > _tetra_vertices;
+    std::vector< unsigned short > _tetra_labels;
 };
\ No newline at end of file
diff --git a/tests/test_tetrahedral_mesh.py b/tests/test_tetrahedral_mesh.py
index 3ae1ffd962f7e0e0392a5e40869eebedc3fc25cb..778cf7f3c25f9c478e75f1cc6c801707a1615638 100644
--- a/tests/test_tetrahedral_mesh.py
+++ b/tests/test_tetrahedral_mesh.py
@@ -14,3 +14,20 @@ class TestTetrahedralMesh(unittest.TestCase):
         mesh.add_point([0, 0, 0])
         points = mesh.vertex_points()
         assert points.shape == (1, 3)
+
+    def test_tetrahedral_mesh_tetras(self):
+        mesh = LabelledTetrahedralMesh()
+        mesh.add_point([0, 0, 0])
+        mesh.add_point([1, 0, 0])
+        mesh.add_point([0, 1, 0])
+        mesh.add_point([0, 0, 1])
+        points = mesh.vertex_points()
+        assert points.shape == (4, 3)
+
+        mesh.add_tetra([0, 1, 2, 3], 0)
+        tetras = mesh.tetra_vertices()
+        assert tetras.shape == (1, 4)
+
+        tetra_labels = mesh.tetra_labels()
+        assert tetra_labels.shape == (1,)
+        assert tetra_labels[0] == 0
diff --git a/wrp/LabelledTetrahedralMeshWrapper.cpp b/wrp/LabelledTetrahedralMeshWrapper.cpp
index e07adf85fb2e390cb106ff81d5fb56ec37f44426..b38c2cea5959ab58f6abf907b6e8d47026d3f59b 100644
--- a/wrp/LabelledTetrahedralMeshWrapper.cpp
+++ b/wrp/LabelledTetrahedralMeshWrapper.cpp
@@ -10,5 +10,8 @@ PYBIND11_MODULE(_wrp, m) {
 py::class_<LabelledTetrahedralMesh>(m, "LabelledTetrahedralMesh")
     .def(py::init<>())
     .def("vertex_points", &LabelledTetrahedralMesh::vertexPoints)
-    .def("add_point", &LabelledTetrahedralMesh::addPoint, py::arg("point"));
+    .def("add_point", &LabelledTetrahedralMesh::addPoint, py::arg("point"))
+    .def("tetra_vertices", &LabelledTetrahedralMesh::tetraVertices)
+    .def("tetra_labels", &LabelledTetrahedralMesh::tetraLabels)
+    .def("add_tetra", &LabelledTetrahedralMesh::addTetra, py::arg("tetra"), py::arg("label"));
 }
\ No newline at end of file