diff --git a/elements/src/CMakeLists.txt b/elements/src/CMakeLists.txt
index 4360252728b0a176a982801249b299ebc5f9c5f6..abc32eb4f651af4bdc3589aaa4fb41d0af45045d 100644
--- a/elements/src/CMakeLists.txt
+++ b/elements/src/CMakeLists.txt
@@ -58,6 +58,7 @@ add_library(${PROJECT_NAME}-elements
             materialmodels/compute_stress.cc
             materialmodels/compute_tangent_matrix.cc
             mesh/Mesh.cc
+            mesh/generate_cuboid_mesh.cc
             mesh/generate_quadratic_tetrahedron_mesh.cc
             mesh/generate_quadratic_triangle_mesh.cc
             mesh/generate_tetrahedron_mesh.cc
diff --git a/elements/src/include/ae108/elements/mesh/generate_cuboid_mesh.h b/elements/src/include/ae108/elements/mesh/generate_cuboid_mesh.h
new file mode 100644
index 0000000000000000000000000000000000000000..1685d650bc87592c1f6832b000a43423d15bbf1e
--- /dev/null
+++ b/elements/src/include/ae108/elements/mesh/generate_cuboid_mesh.h
@@ -0,0 +1,39 @@
+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#pragma once
+
+#include "ae108/elements/mesh/Mesh.h"
+#include "ae108/elements/tensor/Tensor.h"
+#include <cstddef>
+#include <vector>
+
+namespace ae108 {
+namespace elements {
+namespace mesh {
+
+/**
+ * @brief Creates a simple mesh of cuboid shape by splitting
+ * into smaller cuboids as defined by `granularity`.
+ *
+ * @param size The vertex of the cuboid that is opposite to [0., 0., 0.].
+ * @param granularity The number of cuboids along the three axes.
+ */
+Mesh<tensor::Tensor<double, 3>> generate_cuboid_mesh(
+    const tensor::Tensor<double, 3> &size,
+    const tensor::Tensor<std::size_t, 3> &granularity) noexcept;
+
+} // namespace mesh
+} // namespace elements
+} // namespace ae108
\ No newline at end of file
diff --git a/elements/src/mesh/generate_cuboid_mesh.cc b/elements/src/mesh/generate_cuboid_mesh.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c67c9b187c4a9e319174a0ebd148d2c79daeb087
--- /dev/null
+++ b/elements/src/mesh/generate_cuboid_mesh.cc
@@ -0,0 +1,141 @@
+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ae108/elements/mesh/generate_cuboid_mesh.h"
+#include <array>
+#include <numeric>
+
+namespace ae108 {
+namespace elements {
+namespace mesh {
+
+namespace {
+
+using Point = Mesh<tensor::Tensor<double, 3>>::Point;
+using size_type = Mesh<Point>::size_type;
+using value_type = Point::value_type;
+using Positions = Mesh<Point>::Positions;
+using Connectivity = Mesh<Point>::Connectivity;
+using Index = tensor::Tensor<size_type, 3>;
+
+constexpr auto points_per_cuboid = size_type{8};
+
+/**
+ * @brief Computes the number of smaller cuboids the cuboid
+ * will be split into.
+ */
+size_type number_of_cuboids(const Index &granularity) noexcept {
+  return std::accumulate(granularity.begin(), granularity.end(), size_type{1},
+                         std::multiplies<size_type>{});
+}
+
+/**
+ * @brief Converts the steps in x, y, and z
+ * direction to the ID of the point at that location.
+ */
+size_type steps_to_id(const Index &steps, const Index &granularity) noexcept {
+  return steps[2] * (granularity[1] + size_type{1}) *
+             (granularity[0] + size_type{1}) +
+         steps[1] * (granularity[0] + size_type{1}) + steps[0];
+}
+
+/**
+ * @brief Converts the steps in x, y, and z
+ * direction to the position of the point at that location.
+ */
+Point steps_to_position(const Index &steps, const Point &size,
+                        const Index &granularity) noexcept {
+  return {{
+      static_cast<value_type>(steps[0]) * static_cast<value_type>(size[0]) /
+          static_cast<value_type>(granularity[0]),
+      static_cast<value_type>(steps[1]) * static_cast<value_type>(size[1]) /
+          static_cast<value_type>(granularity[1]),
+      static_cast<value_type>(steps[2]) * static_cast<value_type>(size[2]) /
+          static_cast<value_type>(granularity[2]),
+  }};
+}
+
+/**
+ * @brief Generates the positions of the points in the mesh.
+ */
+Positions generate_positions(const Point &size,
+                             const Index &granularity) noexcept {
+  auto positions = Positions();
+  positions.reserve(number_of_cuboids(granularity));
+
+  auto steps = Index();
+  for (steps[2] = 0; steps[2] <= granularity[2]; ++steps[2])
+    for (steps[1] = 0; steps[1] <= granularity[1]; ++steps[1])
+      for (steps[0] = 0; steps[0] <= granularity[0]; ++steps[0])
+        positions.emplace_back(steps_to_position(steps, size, granularity));
+
+  return positions;
+}
+
+/**
+ * @brief Generates the connectivity of the mesh.
+ */
+Connectivity generate_connectivity(const Point &size,
+                                   const Index &granularity) noexcept {
+  auto connectivity = Connectivity();
+  connectivity.reserve(number_of_cuboids(granularity));
+
+  constexpr std::array<Index, points_per_cuboid> steps = {{
+      {{0, 0, 0}},
+      {{1, 0, 0}},
+      {{1, 1, 0}},
+      {{0, 1, 0}},
+      {{0, 0, 1}},
+      {{1, 0, 1}},
+      {{1, 1, 1}},
+      {{0, 1, 1}},
+  }};
+
+  const auto add_cuboid = [&](const Index &offset) {
+    connectivity.emplace_back();
+    connectivity.back().reserve(points_per_cuboid);
+    for (auto &&step : steps) {
+      connectivity.back().push_back(steps_to_id({{
+                                                    step[0] + offset[0],
+                                                    step[1] + offset[1],
+                                                    step[2] + offset[2],
+                                                }},
+                                                granularity));
+    }
+  };
+
+  auto offset = Index();
+  for (offset[2] = 0; offset[2] < granularity[2]; ++offset[2])
+    for (offset[1] = 0; offset[1] < granularity[1]; ++offset[1])
+      for (offset[0] = 0; offset[0] < granularity[0]; ++offset[0])
+        add_cuboid(offset);
+
+  return connectivity;
+}
+
+} // namespace
+
+Mesh<Point> generate_cuboid_mesh(const Point &size,
+                                 const Index &granularity) noexcept {
+  if (number_of_cuboids(granularity) == 0) {
+    return Mesh<Point>({}, {});
+  }
+
+  return Mesh<Point>(generate_connectivity(size, granularity),
+                     generate_positions(size, granularity));
+}
+
+} // namespace mesh
+} // namespace elements
+} // namespace ae108
\ No newline at end of file
diff --git a/elements/test/CMakeLists.txt b/elements/test/CMakeLists.txt
index 7d925ade8b2d3b2c3e5c25240750a0d1bdbd2a00..49076723e46a3e2adc9e8e0910bddfe19089ace3 100644
--- a/elements/test/CMakeLists.txt
+++ b/elements/test/CMakeLists.txt
@@ -24,6 +24,7 @@ add_executable(${PROJECT_NAME}-ElementsTests
                materialmodels/Hookean_Test.cc
                materialmodels/Minimal_Test.cc
                mesh/Mesh_Test.cc
+               mesh/generate_cuboid_mesh_Test.cc
                mesh/generate_quadratic_tetrahedron_mesh_Test.cc
                mesh/generate_quadratic_triangle_mesh_Test.cc
                mesh/generate_tetrahedron_mesh_Test.cc
diff --git a/elements/test/mesh/generate_cuboid_mesh_Test.cc b/elements/test/mesh/generate_cuboid_mesh_Test.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f8967bbd2b61ee97e0d8d2f5fb9e831a7859ca1f
--- /dev/null
+++ b/elements/test/mesh/generate_cuboid_mesh_Test.cc
@@ -0,0 +1,286 @@
+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ae108/elements/mesh/generate_cuboid_mesh.h"
+#include <gmock/gmock.h>
+
+using testing::DoubleEq;
+using testing::ElementsAre;
+using testing::Eq;
+using testing::IsEmpty;
+using testing::Test;
+
+namespace ae108 {
+namespace elements {
+namespace mesh {
+namespace {
+
+struct generate_cuboid_mesh_Test : Test {};
+
+TEST_F(generate_cuboid_mesh_Test, returns_empty_mesh_for_zero_x_granularity) {
+  const auto mesh = generate_cuboid_mesh({{1., 1., 1.}}, {{0, 1, 1}});
+
+  EXPECT_THAT(mesh.number_of_positions(), Eq(0));
+  EXPECT_THAT(mesh.connectivity(), IsEmpty());
+}
+
+TEST_F(generate_cuboid_mesh_Test, returns_empty_mesh_for_zero_y_granularity) {
+  const auto mesh = generate_cuboid_mesh({{1., 1., 1.}}, {{1, 0, 1}});
+
+  EXPECT_THAT(mesh.number_of_positions(), Eq(0));
+  EXPECT_THAT(mesh.connectivity(), IsEmpty());
+}
+
+TEST_F(generate_cuboid_mesh_Test, returns_empty_mesh_for_zero_z_granularity) {
+  const auto mesh = generate_cuboid_mesh({{1., 1., 1.}}, {{1, 1, 0}});
+
+  EXPECT_THAT(mesh.number_of_positions(), Eq(0));
+  EXPECT_THAT(mesh.connectivity(), IsEmpty());
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_positions_for_cube) {
+  const auto mesh = generate_cuboid_mesh({{1., 1., 1.}}, {{1, 1, 1}});
+
+  ASSERT_THAT(mesh.number_of_positions(), Eq(8));
+  EXPECT_THAT(mesh.position_of_vertex(0),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(1),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(2),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(3),
+              ElementsAre(DoubleEq(1.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(4),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(5),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(6),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(7),
+              ElementsAre(DoubleEq(1.), DoubleEq(1.), DoubleEq(1.)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_connectivity_for_cube) {
+  const auto mesh = generate_cuboid_mesh({{1., 1., 1.}}, {{1, 1, 1}});
+
+  EXPECT_THAT(mesh.connectivity(),
+              ElementsAre(ElementsAre(0, 1, 3, 2, 4, 5, 7, 6)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_positions_for_x_scaled_cube) {
+  const auto mesh = generate_cuboid_mesh({{2., 1., 1.}}, {{1, 1, 1}});
+
+  ASSERT_THAT(mesh.number_of_positions(), Eq(8));
+  EXPECT_THAT(mesh.position_of_vertex(0),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(1),
+              ElementsAre(DoubleEq(2.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(2),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(3),
+              ElementsAre(DoubleEq(2.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(4),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(5),
+              ElementsAre(DoubleEq(2.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(6),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(7),
+              ElementsAre(DoubleEq(2.), DoubleEq(1.), DoubleEq(1.)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_positions_for_y_scaled_cube) {
+  const auto mesh = generate_cuboid_mesh({{1., 2., 1.}}, {{1, 1, 1}});
+
+  ASSERT_THAT(mesh.number_of_positions(), Eq(8));
+  EXPECT_THAT(mesh.position_of_vertex(0),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(1),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(2),
+              ElementsAre(DoubleEq(0.), DoubleEq(2.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(3),
+              ElementsAre(DoubleEq(1.), DoubleEq(2.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(4),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(5),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(6),
+              ElementsAre(DoubleEq(0.), DoubleEq(2.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(7),
+              ElementsAre(DoubleEq(1.), DoubleEq(2.), DoubleEq(1.)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_positions_for_z_scaled_cube) {
+  const auto mesh = generate_cuboid_mesh({{1., 1., 2.}}, {{1, 1, 1}});
+
+  ASSERT_THAT(mesh.number_of_positions(), Eq(8));
+  EXPECT_THAT(mesh.position_of_vertex(0),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(1),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(2),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(3),
+              ElementsAre(DoubleEq(1.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(4),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(2.)));
+  EXPECT_THAT(mesh.position_of_vertex(5),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(2.)));
+  EXPECT_THAT(mesh.position_of_vertex(6),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(2.)));
+  EXPECT_THAT(mesh.position_of_vertex(7),
+              ElementsAre(DoubleEq(1.), DoubleEq(1.), DoubleEq(2.)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_connectivity_for_x_scaled_cube) {
+  const auto mesh = generate_cuboid_mesh({{2., 1., 1.}}, {{1, 1, 1}});
+
+  EXPECT_THAT(mesh.connectivity(),
+              ElementsAre(ElementsAre(0, 1, 3, 2, 4, 5, 7, 6)));
+  ;
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_connectivity_for_y_scaled_cube) {
+  const auto mesh = generate_cuboid_mesh({{1., 2., 1.}}, {{1, 1, 1}});
+
+  EXPECT_THAT(mesh.connectivity(),
+              ElementsAre(ElementsAre(0, 1, 3, 2, 4, 5, 7, 6)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_connectivity_for_z_scaled_cube) {
+  const auto mesh = generate_cuboid_mesh({{1., 1., 2.}}, {{1, 1, 1}});
+
+  EXPECT_THAT(mesh.connectivity(),
+              ElementsAre(ElementsAre(0, 1, 3, 2, 4, 5, 7, 6)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_positions_for_x_split_cuboid) {
+  const auto mesh = generate_cuboid_mesh({{2., 1., 1.}}, {{2, 1, 1}});
+
+  ASSERT_THAT(mesh.number_of_positions(), Eq(12));
+  EXPECT_THAT(mesh.position_of_vertex(0),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(1),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(2),
+              ElementsAre(DoubleEq(2.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(3),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(4),
+              ElementsAre(DoubleEq(1.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(5),
+              ElementsAre(DoubleEq(2.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(6),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(7),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(8),
+              ElementsAre(DoubleEq(2.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(9),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(10),
+              ElementsAre(DoubleEq(1.), DoubleEq(1.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(11),
+              ElementsAre(DoubleEq(2.), DoubleEq(1.), DoubleEq(1.)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_positions_for_y_split_cuboid) {
+  const auto mesh = generate_cuboid_mesh({{1., 2., 1.}}, {{1, 2, 1}});
+
+  ASSERT_THAT(mesh.number_of_positions(), Eq(12));
+  EXPECT_THAT(mesh.position_of_vertex(0),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(1),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(2),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(3),
+              ElementsAre(DoubleEq(1.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(4),
+              ElementsAre(DoubleEq(0.), DoubleEq(2.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(5),
+              ElementsAre(DoubleEq(1.), DoubleEq(2.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(6),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(7),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(8),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(9),
+              ElementsAre(DoubleEq(1.), DoubleEq(1.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(10),
+              ElementsAre(DoubleEq(0.), DoubleEq(2.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(11),
+              ElementsAre(DoubleEq(1.), DoubleEq(2.), DoubleEq(1.)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_positions_for_z_split_cuboid) {
+  const auto mesh = generate_cuboid_mesh({{1., 1., 2.}}, {{1, 1, 2}});
+
+  ASSERT_THAT(mesh.number_of_positions(), Eq(12));
+  EXPECT_THAT(mesh.position_of_vertex(0),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(1),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(2),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(3),
+              ElementsAre(DoubleEq(1.), DoubleEq(1.), DoubleEq(0.)));
+  EXPECT_THAT(mesh.position_of_vertex(4),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(5),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(6),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(7),
+              ElementsAre(DoubleEq(1.), DoubleEq(1.), DoubleEq(1.)));
+  EXPECT_THAT(mesh.position_of_vertex(8),
+              ElementsAre(DoubleEq(0.), DoubleEq(0.), DoubleEq(2.)));
+  EXPECT_THAT(mesh.position_of_vertex(9),
+              ElementsAre(DoubleEq(1.), DoubleEq(0.), DoubleEq(2.)));
+  EXPECT_THAT(mesh.position_of_vertex(10),
+              ElementsAre(DoubleEq(0.), DoubleEq(1.), DoubleEq(2.)));
+  EXPECT_THAT(mesh.position_of_vertex(11),
+              ElementsAre(DoubleEq(1.), DoubleEq(1.), DoubleEq(2.)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_connectivity_for_x_split_cuboid) {
+  const auto mesh = generate_cuboid_mesh({{2., 1., 1.}}, {{2, 1, 1}});
+
+  EXPECT_THAT(mesh.connectivity(),
+              ElementsAre(ElementsAre(0, 1, 4, 3, 6, 7, 10, 9),
+                          ElementsAre(1, 2, 5, 4, 7, 8, 11, 10)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_connectivity_for_y_split_cuboid) {
+  const auto mesh = generate_cuboid_mesh({{1., 2., 1.}}, {{1, 2, 1}});
+
+  EXPECT_THAT(mesh.connectivity(),
+              ElementsAre(ElementsAre(0, 1, 3, 2, 6, 7, 9, 8),
+                          ElementsAre(2, 3, 5, 4, 8, 9, 11, 10)));
+}
+
+TEST_F(generate_cuboid_mesh_Test, has_correct_connectivity_for_z_split_cuboid) {
+  const auto mesh = generate_cuboid_mesh({{1., 1., 2.}}, {{1, 1, 2}});
+
+  EXPECT_THAT(mesh.connectivity(),
+              ElementsAre(ElementsAre(0, 1, 3, 2, 4, 5, 7, 6),
+                          ElementsAre(4, 5, 7, 6, 8, 9, 11, 10)));
+}
+
+} // namespace
+} // namespace mesh
+} // namespace elements
+} // namespace ae108
\ No newline at end of file