Merge ae108-meshing/main into '132-merge-ae108-meshing-into-ae108'

.gitignore

@@ -16,4 +16,5 @@
 # along with ae108. If not, see <https://www.gnu.org/licenses/>.
 
 build
-*.pyc
+*.pyc
\ No newline at end of file
+.vscode

CMakeLists.txt

@@ -27,7 +27,7 @@ include(cmake/ae108_install_library.cmake)
 include(cmake/ae108_add_cppptest.cmake)
 
 add_subdirectory(cppptest)
-foreach(AE108_LIBRARY elements cpppetsc cppslepc assembly solve cmdline)
+foreach(AE108_LIBRARY elements cpppetsc cppslepc assembly solve cmdline meshing)
     add_subdirectory(${AE108_LIBRARY})
 endforeach()
 add_subdirectory(examples)
@@ -50,6 +50,8 @@ install(FILES
         cmake/modules/FindAE108_PETSc.cmake
         cmake/modules/FindAE108_SLEPc.cmake
         cmake/modules/AE108_PETSc.cc
+        cmake/modules/FindVoro++.cmake
+        cmake/modules/FindGmsh.cmake
         DESTINATION "${CMAKE_INSTALL_LIBDIR}/cmake/${PROJECT_NAME}/modules"
 )
 
@@ -65,7 +67,7 @@ set(CPACK_GENERATOR "DEB")
 set(CPACK_PACKAGE_NAME "libae108-${AE108_PETSC_SCALAR_TYPE}-dev")
 set(CPACK_DEBIAN_PACKAGE_SECTION "science")
 set(CPACK_DEBIAN_PACKAGE_ARCHITECTURE "amd64")
-set(CPACK_DEBIAN_PACKAGE_DEPENDS "libboost-program-options-dev,libeigen3-dev,librange-v3-dev,libxrender1,libvtk9-dev,libpetsc-${AE108_PETSC_SCALAR_TYPE}-dev,libslepc-${AE108_PETSC_SCALAR_TYPE}3.15-dev")
+set(CPACK_DEBIAN_PACKAGE_DEPENDS "libboost-program-options-dev,libeigen3-dev,librange-v3-dev,libxrender1,libvtk9-dev,libpetsc-${AE108_PETSC_SCALAR_TYPE}-dev,libslepc-${AE108_PETSC_SCALAR_TYPE}3.15-dev,libgmsh-dev,voro++-dev")
 set(CPACK_PACKAGE_DESCRIPTION "C++ foundation for computational solid mechanics simulations using a variational framework, primarily focusing on the Finite Element Method (FEM)")
 set(CPACK_PACKAGE_CONTACT "Mechanics & Materials Lab at ETH Zurich <gerhard.braeunlich@id.ethz.ch>")
 

cmake/modules/FindGmsh.cmake

+# © 2021 ETH Zurich, Mechanics and Materials Lab
+#
+# Licensed under the GNU General Public License (GPL), 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
+#
+# https:#www.gnu.org/licenses/old-licenses/gpl-2.0.html
+#
+# 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.
+
+if (NOT TARGET external::gmsh)
+  include(FindPackageHandleStandardArgs)
+
+  if (NOT DEFINED ENV{GMSH_LIBRARIES} OR NOT DEFINED ENV{GMSH_INCLUDE_DIRS})
+    find_library(Gmsh_LIBRARIES NAMES gmsh)
+    find_path(Gmsh_INCLUDE_DIRS gmsh.h)
+  else()
+    find_library(Gmsh_LIBRARIES NAMES gmsh PATHS $ENV{GMSH_LIBRARIES} NO_DEFAULT_PATH)
+    find_path(Gmsh_INCLUDE_DIRS gmsh.h PATHS $ENV{GMSH_INCLUDE_DIRS} NO_DEFAULT_PATH)
+  endif()
+
+  find_package_handle_standard_args(Gmsh Gmsh_LIBRARIES Gmsh_INCLUDE_DIRS)
+
+  if (Gmsh_FOUND)
+    add_library(external::gmsh UNKNOWN IMPORTED)
+    set_property(TARGET external::gmsh PROPERTY IMPORTED_LOCATION ${Gmsh_LIBRARIES})
+    set_property(TARGET external::gmsh PROPERTY INTERFACE_INCLUDE_DIRECTORIES ${Gmsh_INCLUDE_DIRS})
+  endif()
+endif()
+
+

cmake/modules/FindVoro++.cmake

+# © 2021 ETH Zurich, Mechanics and Materials Lab
+#
+# This file is part of ae108.
+#
+# ae108 is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or any
+# later version.
+#
+# ae108 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+if (NOT TARGET external::voro++)
+  include(FindPackageHandleStandardArgs)
+  if (NOT DEFINED ENV{Voro++_LIBRARIES} OR NOT DEFINED ENV{Voro++_INCLUDE_DIRS})
+    find_library(Voro++_LIBRARIES NAMES voro++ PATHS $ENV{VOROPP_DIR}/lib)
+    find_path(Voro++_INCLUDE_DIRS voro++.hh PATH_SUFFIXES voro++ PATHS $ENV{VOROPP_DIR}/include/voro++)
+  else()
+    find_library(Voro++_LIBRARIES NAMES voro++ PATHS $ENV{Voro++_LIBRARIES} NO_DEFAULT_PATH)
+    find_path(Voro++_INCLUDE_DIRS voro++.hh PATH_SUFFIXES voro++ PATHS $ENV{Voro++_INCLUDE_DIRS} NO_DEFAULT_PATH)
+  endif()
+
+  find_package_handle_standard_args(Voro++ DEFAULT_MSG Voro++_LIBRARIES Voro++_INCLUDE_DIRS)
+
+  if (Voro++_FOUND)
+    add_library(external::voro++ UNKNOWN IMPORTED)
+    set_property(TARGET external::voro++ PROPERTY IMPORTED_LOCATION ${Voro++_LIBRARIES})
+    set_property(TARGET external::voro++ PROPERTY INTERFACE_INCLUDE_DIRECTORIES ${Voro++_INCLUDE_DIRS})
+  endif()
+endif()

docker-compose.yml

@@ -24,7 +24,7 @@ services:
       context: ./docker
       dockerfile: Dockerfile
     volumes:
-    - .:/mnt/io
+      - .:/mnt/io
     tty: true
     shm_size: '1gb'
     working_dir: /mnt/io
@@ -35,7 +35,7 @@ services:
       context: ./docker
       dockerfile: Dockerfile
     volumes:
-    - .:/mnt/io
+      - .:/mnt/io
     working_dir: /mnt/io
     tty: true
     shm_size: '1gb'
@@ -44,6 +44,6 @@ services:
       context: ./docker
       dockerfile: Dockerfile_docu
     volumes:
-    - .:/mnt/io
+      - .:/mnt/io
     working_dir: /mnt/io
-    tty: true
+    tty: true
\ No newline at end of file

docker/Dockerfile

@@ -38,7 +38,14 @@ RUN apt-get update && \
   libxrender1 \
   ninja-build \
   python3 python3-pip \
-  python3-scipy
+  python3-scipy 
+
+RUN apt-get update && \
+  DEBIAN_FRONTEND="noninteractive" \
+  TZ="Europe/Zurich" \
+  apt-get install -y \
+  libgmsh-dev \
+  voro++-dev
 
 RUN pip3 install \
   mypy==0.931 \

examples/CMakeLists.txt

@@ -156,3 +156,5 @@ if (NOT ${AE108_PETSC_COMPLEX})
                           PRIVATE ae108::elements
     )
 endif()
+
+add_subdirectory(meshing)
\ No newline at end of file

examples/meshing/CMakeLists.txt

+# © 2020 ETH Zurich, Mechanics and Materials Lab
+#
+# This file is part of ae108.
+#
+# ae108 is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or any
+# later version.
+#
+# ae108 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
+
+project(ae108-examples LANGUAGES C CXX)
+
+if (NOT TARGET ae108::meshing)
+    find_package(ae108 REQUIRED)
+endif()
+
+add_executable(${PROJECT_NAME}-SegmentMesh SegmentMesh.cc)
+target_link_libraries(${PROJECT_NAME}-SegmentMesh
+                        PRIVATE ae108::meshing
+)
+
+add_executable(${PROJECT_NAME}-PeriodicRVEMesh PeriodicRVEMesh.cc)
+target_link_libraries(${PROJECT_NAME}-PeriodicRVEMesh
+                        PRIVATE ae108::meshing
+)
+
+add_executable(${PROJECT_NAME}-TwoPhaseCube TwoPhaseCube.cc)
+target_link_libraries(${PROJECT_NAME}-TwoPhaseCube
+                        PRIVATE ae108::meshing
+)
+
+add_executable(${PROJECT_NAME}-PhysicalGroups PhysicalGroups.cc)
+target_link_libraries(${PROJECT_NAME}-PhysicalGroups
+                        PRIVATE ae108::meshing
+)
+add_executable(${PROJECT_NAME}-GradedLattice GradedLattice.cc)
+target_link_libraries(${PROJECT_NAME}-GradedLattice
+                        PRIVATE ae108::meshing
+)

examples/meshing/GradedLattice.cc

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// ae108 is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or any
+// later version.
+//
+// ae108 is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+#include "ae108/meshing/construct_rectilinear_grid.h"
+#include "ae108/meshing/cppgmsh/Context.h"
+#include "ae108/meshing/cppgmsh/construct_box.h"
+#include "ae108/meshing/cppgmsh/construct_rectangle.h"
+#include "ae108/meshing/cppgmsh/extract_mesh.h"
+#include "ae108/meshing/cppgmsh/extrude_surface.h"
+#include "ae108/meshing/cppgmsh/fuse_entities.h"
+#include "ae108/meshing/cppgmsh/generate_mesh.h"
+#include "ae108/meshing/cppgmsh/remove_entities.h"
+#include "ae108/meshing/cppgmsh/rotate_entities.h"
+#include "ae108/meshing/cppgmsh/set_granularity.h"
+#include "ae108/meshing/cppgmsh/synchronize.h"
+#include "ae108/meshing/cppgmsh/translate_entities.h"
+#include "ae108/meshing/cppgmsh/write_file.h"
+
+#include <gmsh.h>
+#include <iostream>
+#include <numeric>
+
+using namespace ae108::meshing;
+
+constexpr std::size_t dimension = 3;
+using Point = std::array<double, dimension>;
+
+int main(int argc, char **argv) {
+
+  const auto gmshContext = cppgmsh::Context(argc, argv, true, false);
+
+  const double gamma_start = 5.;
+  const double gamma_end = 10.;
+  const double length = 466.75;
+
+  const auto xcoords = [] {
+    std::vector<double> range(80);
+    std::iota(range.begin(), range.end(), -40);
+    std::transform(range.begin(), range.end(), range.begin(),
+                   [](auto &N) { return N * 5; });
+    return range;
+  }();
+
+  const auto ycoords = [&] {
+    std::vector<double> range(59);
+    std::iota(range.begin(), range.end(), -29);
+    std::transform(range.begin(), range.end(), range.begin(), [&](auto &N) {
+      return N * length * gamma_start /
+             (length - fabs(N) * gamma_end + fabs(N) * gamma_start);
+    });
+    return range;
+  }();
+
+  const auto crosssection =
+      cppgmsh::construct_rectangle({-.25, -.5, 0}, {.5, 1.});
+  cppgmsh::rotate_entities({crosssection}, {0, 0, 0}, {1, 0, 0}, M_PI_2);
+
+  std::vector<std::pair<int, int>> entities;
+  for (const auto &xcoord : xcoords)
+    entities.push_back(cppgmsh::extrude_surface(
+        crosssection.second, {xcoord, -200, 0}, {xcoord, 200, 0}, {0, 1, 0}));
+  for (const auto &ycoord : ycoords)
+    entities.push_back(cppgmsh::extrude_surface(
+        crosssection.second, {-200, ycoord, 0}, {200, ycoord, 0}, {0, 1, 0}));
+
+  entities.push_back(cppgmsh::construct_box({-225, -225, -.5}, {450, 25, 1}));
+  entities.push_back(cppgmsh::construct_box({-225, -225, -.5}, {25, 450, 1}));
+  entities.push_back(cppgmsh::construct_box({-225, 200, -.5}, {450, 25, 1}));
+  entities.push_back(cppgmsh::construct_box({200, -225, -.5}, {25, 450, 1}));
+
+  cppgmsh::remove_entities({crosssection}, true);
+
+  cppgmsh::fuse_entities(entities);
+
+  cppgmsh::synchronize();
+  cppgmsh::write_file("gradedLattice.step");
+
+  cppgmsh::generate_mesh();
+  cppgmsh::write_file("gradedLattice.vtk");
+
+  const auto [positions, connectivity, _, __] =
+      cppgmsh::extract_mesh<dimension>();
+
+  std::cout << "# points: " << positions.size() << std::endl
+            << "# elements: " << connectivity.size() << std::endl;
+}
\ No newline at end of file

examples/meshing/PeriodicRVEMesh.cc

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// ae108 is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or any
+// later version.
+//
+// ae108 is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+#include "ae108/meshing/bounding_box_of.h"
+#include "ae108/meshing/construct_periodic_point_cloud.h"
+#include "ae108/meshing/construct_voronoi_cell.h"
+#include "ae108/meshing/cppgmsh/Context.h"
+#include "ae108/meshing/cppgmsh/construct_cylinders.h"
+#include "ae108/meshing/cppgmsh/construct_polyhedron.h"
+#include "ae108/meshing/cppgmsh/extract_mesh.h"
+#include "ae108/meshing/cppgmsh/generate_mesh.h"
+#include "ae108/meshing/cppgmsh/get_entities_in.h"
+#include "ae108/meshing/cppgmsh/get_periodic_nodes_of.h"
+#include "ae108/meshing/cppgmsh/intersect_entities.h"
+#include "ae108/meshing/cppgmsh/set_domain_entities_periodic.h"
+#include "ae108/meshing/cppgmsh/set_granularity.h"
+#include "ae108/meshing/cppgmsh/synchronize.h"
+#include "ae108/meshing/cppgmsh/write_file.h"
+
+#include <iostream>
+
+using namespace ae108::meshing;
+
+// In this example we will generate a periodic mesh of a representative volume
+// element as it is often required for homogenization or Bloch wave analysis of
+// a periodic medium.
+
+constexpr std::size_t dimension = 3;
+using Point = std::array<double, dimension>;
+
+// A 3D periodic medium with translational symmetry is commonly described as a
+// Bravais lattice by three translation vectors. Three well-known lattices are
+// the simple cubic (sc), the body centered cubic (bcc), and the face centered
+// cubic (fcc) lattice. https://en.wikipedia.org/wiki/Bravais_lattice
+
+constexpr auto lattice_vectors =
+    std::array<std::array<double, dimension>, dimension>{
+        {{0.5, 0.5, 0.5}, {0.5, -0.5, 0.5}, {0.5, 0.5, -0.5}}}; // bcc
+//      {{1., 1., 0.}, {0., 1., 1.}, {1., 0., 1.}}}; //fcc
+//      {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}}; //sc
+
+int main(int argc, char **argv) {
+
+  // In a first step, we need to initialize the meshing engine "cppgmsh"
+
+  const auto gmshContext = cppgmsh::Context(argc, argv);
+
+  // Let us now define the domain of our representative volume element (rve).
+  // For a perodic medium, like above Bravais lattice, the smallest possible
+  // unit cell is a primitive unit cell: Two frequent choices are the
+  // parallelpiped or the Wigner-Seitz (a type of Voronoi cell) primitive cell.
+  // https://en.wikipedia.org/wiki/Unit_cell
+
+  std::vector<std::pair<std::size_t, std::size_t>> periodic_faces;
+  const auto domain = construct_voronoi_cell(
+      construct_periodic_point_cloud(lattice_vectors, {0, 0, 0}),
+      &periodic_faces);
+  const auto domain_entity = cppgmsh::construct_polyhedron(domain);
+
+  // Next, we define the solid within the primitive unit cell. We choose to
+  // place cylinders on the edges of the domain. The solid may also span several
+  // unit cells like in a finite lattice.
+
+  const auto solid_entity = cppgmsh::construct_cylinders(
+      domain.vertices, domain.edges,
+      std::vector<double>(domain.edges.size(), 0.05), true);
+
+  // Finally, we construct the rve by cropping all solid outside of the
+  // domain of the primitive unit cell.
+
+  const auto rve_entity =
+      cppgmsh::intersect_entities({solid_entity}, {domain_entity});
+
+  // We sync the constructive solid geometry engine and write a CAD file.
+
+  cppgmsh::synchronize();
+  //   cppgmsh::write_file("periodic.step");
+
+  // To obtain a periodic mesh (e.g. for homogenization or Bloch wave analysis),
+  // we define the periodic domain faces.
+
+  for (const auto &face_pair : periodic_faces)
+    cppgmsh::set_domain_entities_periodic(
+        bounding_box_of(domain.vertices_of_face(face_pair.first)),
+        bounding_box_of(domain.vertices_of_face(face_pair.second)), 2);
+
+  // Finally, we generate the periodic FE mesh.
+
+  cppgmsh::set_granularity(1. / 10);
+  cppgmsh::generate_mesh();
+
+  // We save the mesh directly as "periodic.vtk". Other formats are available.
+
+  cppgmsh::write_file("periodic.vtk");
+
+  // To use the mesh in an ae108 FE application, we extract the positions and
+  // connectivities of the mesh elements.
+
+  const auto [positions, connectivity, nodeTagToIndex, elementTagToIndex] =
+      cppgmsh::extract_mesh<dimension>();
+
+  for (const auto &face_pair : periodic_faces) {
+    std::pair<int, int> numper_of_nodes;
+    for (const auto &target_surface : cppgmsh::get_entities_in(
+             bounding_box_of(domain.vertices_of_face(face_pair.second)), 2)) {
+      const auto nodes = cppgmsh::get_periodic_nodes_of(target_surface);
+      assert(nodes.first.size() == nodes.second.size());
+      numper_of_nodes.first += nodes.first.size();
+      numper_of_nodes.second += nodes.second.size();
+    }
+    std::cout << "# target face (" << face_pair.second
+              << ") nodes: " << numper_of_nodes.second << std::endl
+              << "# source face (" << face_pair.first
+              << ") nodes: " << numper_of_nodes.first << std::endl
+              << std::endl;
+  }
+}

examples/meshing/PhysicalGroups.cc

+// © 2022 ETH Zurich, Mechanics and Materials Lab
+//
+// ae108 is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or any
+// later version.
+//
+// ae108 is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+#include "ae108/meshing/cppgmsh/Context.h"
+#include "ae108/meshing/cppgmsh/construct_box.h"
+#include "ae108/meshing/cppgmsh/extract_mesh.h"
+#include "ae108/meshing/cppgmsh/fragment_entities.h"
+#include "ae108/meshing/cppgmsh/generate_mesh.h"
+#include "ae108/meshing/cppgmsh/get_elements_in.h"
+#include "ae108/meshing/cppgmsh/get_entities_in.h"
+#include "ae108/meshing/cppgmsh/get_nodes_of.h"
+#include "ae108/meshing/cppgmsh/get_physical_groups.h"
+#include "ae108/meshing/cppgmsh/set_granularity.h"
+#include "ae108/meshing/cppgmsh/set_physical_group_of.h"
+#include "ae108/meshing/cppgmsh/synchronize.h"
+
+#include <array>
+#include <iostream>
+
+using namespace ae108::meshing;
+
+/*
+In this example we will construct and mesh a cube with unit side lengths, which
+consists of three layers (z < 0.3 -> layer 1, 0.3 < z < 0.7 -> layer 2, z > 0.7
+-> layer 3). Layers 1 and 3 are of one material, layer 2 of another. In order to
+account for this, physical groups are defined.
+
+     ________________    _
+    |                |  |
+    |   Material 1   |  | 0.3
+    |________________|  |_
+    |                |  |
+    |   Material 2   |  | 0.4
+    |                |  |
+    |________________|  |_
+    |                |  |
+    |   Material 1   |  | 0.3
+    |________________|  |_
+     ________________
+    |       1.0      |
+
+ z
+ ^
+ |
+    --> x,y
+
+*/
+
+int main(int argc, char **argv) {
+
+  // Initialize the gmsh context:
+  const auto gmshContext = cppgmsh::Context(argc, argv, true, false);
+
+  // Construct three boxes:
+  const auto boxes = std::vector<std::pair<int, int>>{
+      cppgmsh::construct_box({0., 0., 0.}, {1., 1., 0.3}),
+      cppgmsh::construct_box({0., 0., .3}, {1., 1., 0.4}),
+      cppgmsh::construct_box({0., 0., .7}, {1., 1., 0.3})};
+
+  std::cout << "Box 1: dim = " << boxes[0].first
+            << ", tag = " << boxes[0].second << std::endl
+            << "Box 2: dim = " << boxes[1].first
+            << ", tag = " << boxes[1].second << std::endl
+            << "Box 3: dim = " << boxes[2].first
+            << ", tag = " << boxes[2].second << std::endl;
+
+  // Fuse boxes, while keeping track of the original entities:
+  std::vector<std::pair<int, int>> fused_entitites =
+      cppgmsh::fragment_entities(boxes);
+  cppgmsh::synchronize();
+
+  // Define physical groups:
+  cppgmsh::set_physical_group_of({boxes[0], boxes[2]});
+  cppgmsh::set_physical_group_of({boxes[1]});
+
+  // Show that the physical groups have been defined as desired:
+  auto groups = cppgmsh::get_physical_groups();
+  std::cout << "Groups (dim, tag): ";
+  for (const auto &group : groups)
+    std::cout << "(" << group.first << "," << group.second << ") ";
+  std::cout << std::endl;
+
+  // Get the tags of the entities attributed two both groups:
+  for (const auto &group : groups) {
+    std::cout << "Group " << group.second << ", dim: " << group.first
+              << ", entities: ";
+    const auto entities = cppgmsh::get_entities_in(group);
+    for (const auto &entity : entities)
+      std::cout << entity.second << " ";
+    std::cout << std::endl;
+  }
+
+  // Create mesh:
+  cppgmsh::set_granularity(2.);
+  cppgmsh::generate_mesh(3, 1, 6);
+
+  const auto [positions, connectivity, nodeTagToIndex, elementTagToIndex] =
+      cppgmsh::extract_mesh<3>();
+
+  // Print element tags and indides, as well as node indides for the two groups:
+  for (const auto &group : groups) {
+    std::cout << "Elements in group " << group.second << ":\n";
+    std::cout << "\ttag\tindex\tnode indides\n";
+    const auto entities = cppgmsh::get_entities_in(group);
+    for (const auto &entity : entities) {
+      const auto elements = cppgmsh::get_elements_in(entity);
+      for (const auto &element : elements) {
+        const auto index = elementTagToIndex.at(element.second);
+        std::cout << "\t" << element.second << "\t" << index << "\t";
+        for (const auto &node : connectivity[index])
+          std::cout << node << "\t";
+        std::cout << std::endl;
+      }
+    }
+  }
+
+  // Print node tags, indides and coordinates for the two boxes:
+  for (const auto &group : groups) {
+    std::cout << "Nodes in group " << group.second << ":\n";
+    std::cout << "\ttag\tindex\tx\ty\tz\n";
+    const auto entities = cppgmsh::get_entities_in(group);
+    for (const auto &entity : entities) {
+      const auto nodes = cppgmsh::get_nodes_of<3>(entity);
+      for (const auto &node : nodes) {
+        std::cout << "\t" << node.id << "\t" << nodeTagToIndex.at(node.id)
+                  << "\t" << node.position[0] << "\t" << node.position[1]
+                  << "\t" << node.position[2] << std::endl;
+      }
+    }
+  }
+}
\ No newline at end of file

examples/meshing/SegmentMesh.cc

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// ae108 is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or any
+// later version.
+//
+// ae108 is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+#include "ae108/meshing/construct_rectilinear_grid.h"
+#include "ae108/meshing/cppgmsh/Context.h"
+#include "ae108/meshing/cppgmsh/construct_rectangle.h"
+#include "ae108/meshing/cppgmsh/extract_mesh.h"
+#include "ae108/meshing/cppgmsh/extrude_surface.h"
+#include "ae108/meshing/cppgmsh/fuse_entities.h"
+#include "ae108/meshing/cppgmsh/generate_mesh.h"
+#include "ae108/meshing/cppgmsh/remove_entities.h"
+#include "ae108/meshing/cppgmsh/set_granularity.h"
+#include "ae108/meshing/cppgmsh/synchronize.h"
+#include "ae108/meshing/cppgmsh/write_file.h"
+
+#include <gmsh.h>
+#include <iostream>
+
+using namespace ae108::meshing;
+
+constexpr std::size_t dimension = 3;
+using Point = std::array<double, dimension>;
+
+int main(int argc, char **argv) {
+
+  const auto gmshContext = cppgmsh::Context(argc, argv, true, false);
+
+  // construct a regular grid
+  //    *--*
+  //    |  |
+  //    *--*
+
+  const auto coordinates =
+      std::array<std::vector<double>, dimension>{{{0, 1}, {0, 1}, {0}}};
+  auto [nodes, segments] = construct_rectilinear_grid(coordinates);
+
+  const auto crosssection =
+      cppgmsh::construct_rectangle({-.05, -.05, 0}, {.1, .1});
+
+  std::vector<std::pair<int, int>> beams;
+  for (const auto &segment : segments)
+    beams.push_back(cppgmsh::extrude_surface(
+        crosssection.second, nodes[segment[0]], nodes[segment[1]]));
+
+  cppgmsh::remove_entities({crosssection}, true);
+
+  cppgmsh::fuse_entities(beams);
+
+  cppgmsh::synchronize();
+  //   cppgmsh::write_file("segmentMesh.step");
+
+  cppgmsh::set_granularity(0.25);
+  cppgmsh::generate_mesh();
+  cppgmsh::write_file("segmentMesh.vtk");
+
+  const auto [positions, connectivity, _, __] =
+      cppgmsh::extract_mesh<dimension>();
+
+  std::cout << "# points: " << positions.size() << std::endl
+            << "# elements: " << connectivity.size() << std::endl;
+}
\ No newline at end of file

examples/meshing/TwoPhaseCube.cc

+// © 2022 ETH Zurich, Mechanics and Materials Lab
+//
+// ae108 is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or any
+// later version.
+//
+// ae108 is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+#include "ae108/meshing/cppgmsh/Context.h"
+#include "ae108/meshing/cppgmsh/construct_box.h"
+#include "ae108/meshing/cppgmsh/extract_mesh.h"
+#include "ae108/meshing/cppgmsh/fragment_entities.h"
+#include "ae108/meshing/cppgmsh/generate_mesh.h"
+#include "ae108/meshing/cppgmsh/get_elements_in.h"
+#include "ae108/meshing/cppgmsh/get_nodes_of.h"
+#include "ae108/meshing/cppgmsh/set_granularity.h"
+#include "ae108/meshing/cppgmsh/synchronize.h"
+
+#include <array>
+#include <iostream>
+
+using namespace ae108::meshing;
+
+/*
+In this example we will construct and mesh a cube with unit side lengths, which
+consists of two different phases (z < 0.5 -> phase 1,z > 0.5 -> phase 2).
+     ________________    _
+    |                |  |
+    |     Phase 2    |  | 0.5
+    |________________|  |_
+    |                |  |
+    |     Phase 1    |  | 0.5
+    |________________|  |_
+     ________________
+    |       1.0      |
+ z
+ ^
+ |
+    --> x,y
+
+*/
+
+int main(int argc, char **argv) {
+
+  // Initialize the gmsh context:
+  const auto gmshContext = cppgmsh::Context(argc, argv, true, false);
+
+  // Construct two boxes:
+  const auto boxes = std::vector<std::pair<int, int>>{
+      cppgmsh::construct_box({0., 0., 0.}, {1., 1., 0.5}),
+      cppgmsh::construct_box({0., 0., .5}, {1., 1., 0.5})};
+
+  std::cout << "Box 1: dim = " << boxes[0].first
+            << ", tag = " << boxes[0].second << std::endl
+            << "Box 2: dim = " << boxes[1].first
+            << ", tag = " << boxes[1].second << std::endl;
+
+  // Fuse boxes, while keeping track of the original entities:
+  const auto fused_entitites = cppgmsh::fragment_entities(boxes);
+  cppgmsh::synchronize();
+
+  // Create mesh:
+  cppgmsh::set_granularity(2.);
+  cppgmsh::generate_mesh(3, 1, 6);
+
+  const auto [positions, connectivity, nodeTagToIndex, elementTagToIndex] =
+      cppgmsh::extract_mesh<3>();
+
+  // Print element tags and indides, as well as node indides for the two boxes:
+  for (const auto &box : boxes) {
+    std::cout << "Elements in box " << box.second << ":\n"
+              << "\ttag\tindex\tnode indices\n";
+
+    const auto elements = cppgmsh::get_elements_in(box);
+    for (const auto &element : elements) {
+      const auto index = elementTagToIndex.at(element.second);
+      std::cout << "\t" << element.second << "\t" << index << "\t";
+      for (const auto &node : connectivity[index])
+        std::cout << node << "\t";
+      std::cout << std::endl;
+    }
+  }
+
+  // Print node tags, indides and coordinates for the two boxes:
+  for (const auto &box : boxes) {
+    std::cout << "Nodes in box " << box.second << ":\n"
+              << "\ttag\tindex\tx\ty\tz\n";
+
+    const auto nodes = cppgmsh::get_nodes_of<3>(box);
+    for (const auto &node : nodes)
+      std::cout << "\t" << node.id << "\t" << nodeTagToIndex.at(node.id) << "\t"
+                << node.position[0] << "\t" << node.position[1] << "\t"
+                << node.position[2] << std::endl;
+  }
+}
\ No newline at end of file

meshing/CMakeLists.txt

+# © 2021 ETH Zurich, Mechanics and Materials Lab
+#
+# This file is part of ae108.
+#
+# ae108 is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or any
+# later version.
+#
+# ae108 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
+
+add_subdirectory(src)
+add_subdirectory(test)
\ No newline at end of file

meshing/src/BoundaryRepresentation.cc

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// ae108 is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or any
+// later version.
+//
+// ae108 is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+#include "ae108/meshing/BoundaryRepresentation.h"
\ No newline at end of file

meshing/src/BoundingBox.cc

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// ae108 is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or any
+// later version.
+//
+// ae108 is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+#include "ae108/meshing/BoundingBox.h"
\ No newline at end of file

meshing/src/CMakeLists.txt

+# © 2021 ETH Zurich, Mechanics and Materials Lab
+#
+# This file is part of ae108.
+#
+# ae108 is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or any
+# later version.
+#
+# ae108 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
+
+add_library(${PROJECT_NAME}-meshing
+    cppgmsh/Context.cc
+    cppgmsh/Node.cc
+    cppgmsh/as_affine_transform.cc
+    cppgmsh/construct_box.cc
+    cppgmsh/construct_cylinders.cc
+    cppgmsh/construct_polyhedron.cc
+    cppgmsh/construct_rectangle.cc
+    cppgmsh/copy_entities.cc
+    cppgmsh/extrude_entities.cc
+    cppgmsh/extrude_surface.cc
+    cppgmsh/fuse_entities.cc
+    cppgmsh/fragment_entities.cc
+    cppgmsh/intersect_entities.cc
+    cppgmsh/remove_entities.cc
+    cppgmsh/rotate_entities.cc
+    cppgmsh/translate_entities.cc
+    cppgmsh/get_boundary_of.cc
+    cppgmsh/get_centroid_of.cc
+    cppgmsh/get_coords_of.cc
+    cppgmsh/get_elements_in.cc
+    cppgmsh/get_entities_in.cc
+    cppgmsh/get_entities_of.cc
+    cppgmsh/get_nodes_of.cc
+    cppgmsh/get_normal_of.cc
+    cppgmsh/get_periodic_nodes_of.cc
+    cppgmsh/get_physical_groups.cc
+    cppgmsh/get_points_of.cc
+    cppgmsh/set_domain_entities_periodic.cc
+    cppgmsh/set_expert_mode.cc
+    cppgmsh/set_granularity.cc
+    cppgmsh/set_periodic.cc
+    cppgmsh/set_physical_group_of.cc
+    cppgmsh/heal_periodic_surfaces.cc
+    cppgmsh/heal_periodic_domains.cc
+    cppgmsh/import_shapes.cc
+    cppgmsh/read_file.cc
+    cppgmsh/write_file.cc
+    cppgmsh/extract_mesh.cc
+    cppgmsh/generate_mesh.cc
+    cppgmsh/synchronize.cc
+    BoundaryRepresentation.cc
+    BoundingBox.cc
+    bounding_box_of.cc
+    construct_rectilinear_grid.cc
+    construct_periodic_point_cloud.cc
+    construct_voronoi_cell.cc
+)
+
+target_include_directories(${PROJECT_NAME}-meshing PUBLIC
+                           $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
+)
+
+target_compile_features(${PROJECT_NAME}-meshing
+                        PUBLIC cxx_std_11
+)
+
+find_package(Gmsh 4.8.4 REQUIRED MODULE)
+find_package(Voro++ REQUIRED MODULE)
+find_package(Eigen3 3.3 CONFIG REQUIRED)
+find_package(range-v3 0.11.0 CONFIG REQUIRED)
+target_link_libraries(${PROJECT_NAME}-meshing 
+                            PUBLIC external::gmsh
+                            PUBLIC external::voro++
+                            PUBLIC Eigen3::Eigen
+                            PUBLIC range-v3::range-v3)
+
+add_library(${PROJECT_NAME}::meshing ALIAS ${PROJECT_NAME}-meshing)
+
+ae108_install_library(meshing)
\ No newline at end of file

meshing/src/bounding_box_of.cc

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// ae108 is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or any
+// later version.
+//
+// ae108 is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+#include "ae108/meshing/bounding_box_of.h"
\ No newline at end of file

meshing/src/construct_periodic_point_cloud.cc

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// ae108 is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or any
+// later version.
+//
+// ae108 is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+#include "ae108/meshing/construct_periodic_point_cloud.h"
+#include <Eigen/Dense>
+#include <range/v3/view/cartesian_product.hpp>
+#include <range/v3/view/zip.hpp>
+
+namespace ae108 {
+namespace meshing {
+
+constexpr auto get_array = [](auto &&...x) {
+  return std::array{std::forward<decltype(x)>(x)...};
+};
+
+template <std::size_t Dimension>
+const auto construct_point_cloud =
+    [](const std::array<std::array<double, Dimension>, Dimension> &translations,
+       const std::array<double, Dimension> &origin, const auto &permutations) {
+      std::vector<std::array<double, Dimension>> pointCloud(
+          permutations.size());
+      for (auto &&[point, tuple] :
+           ranges::views::zip(pointCloud, permutations)) {
+        const auto permutation = std::apply(get_array, tuple);
+        Eigen::Map<Eigen::Matrix<double, Dimension, 1>>(point.data()) =
+            Eigen::Map<const Eigen::Matrix<double, Dimension, 1>>(
+                origin.data()) +
+            (Eigen::Map<const Eigen::Matrix<double, Dimension, Dimension,
+                                            Eigen::ColMajor>>(
+                 translations.front().data()) *
+             Eigen::Map<const Eigen::Matrix<double, Dimension, 1>>(
+                 permutation.data()));
+      }
+      return pointCloud;
+    };
+
+template <>
+std::vector<std::array<double, 3>> construct_periodic_point_cloud(
+    const std::array<std::array<double, 3>, 3> &translations,
+    std::array<double, 3> origin, std::vector<double> set) noexcept {
+  return construct_point_cloud<3>(
+      translations, origin, ranges::views::cartesian_product(set, set, set));
+}
+
+template <>
+std::vector<std::array<double, 2>> construct_periodic_point_cloud(
+    const std::array<std::array<double, 2>, 2> &translations,
+    std::array<double, 2> origin, std::vector<double> set) noexcept {
+
+  return construct_point_cloud<2>(translations, origin,
+                                  ranges::views::cartesian_product(set, set));
+}
+
+template <>
+std::vector<std::array<double, 1>> construct_periodic_point_cloud(
+    const std::array<std::array<double, 1>, 1> &translations,
+    std::array<double, 1> origin, std::vector<double> set) noexcept {
+  return construct_point_cloud<1>(translations, origin,
+                                  ranges::views::cartesian_product(set));
+}
+
+} // namespace meshing
+} // namespace ae108
\ No newline at end of file

meshing/src/construct_rectilinear_grid.cc

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// ae108 is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or any
+// later version.
+//
+// ae108 is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with ae108. If not, see <https://www.gnu.org/licenses/>.
+
+#include "ae108/meshing/construct_rectilinear_grid.h"
+#include <Eigen/Dense>
+
+namespace ae108 {
+namespace meshing {
+
+template <>
+std::tuple<std::vector<std::array<double, 1>>,
+           std::vector<std::array<std::size_t, 2>>>
+construct_rectilinear_grid(
+    const std::array<std::vector<double>, 1> &instances) noexcept {
+
+  const auto generate_points = [&]() {
+    std::vector<std::array<double, 1>> points;
+    points.reserve(instances[0].size());
+
+    for (const auto &x : instances[0])
+      points.push_back({{x}});
+
+    return points;
+  };
+
+  const auto step_to_index = [&](const std::size_t step_x) { return step_x; };
+
+  const auto generate_edges = [&]() {
+    std::vector<std::array<std::size_t, 2>> edges;
+    edges.reserve(instances[0].size() - 1);
+    std::array<std::size_t, 1> step = {{0}};
+
+    // edges in x
+    for (step[0] = 0; step[0] < instances[0].size() - 1; ++step[0])
+      edges.push_back({{step_to_index(step[0]), step_to_index(step[0] + 1)}});
+
+    return edges;
+  };
+
+  return {generate_points(), generate_edges()};
+}
+
+template <>
+std::tuple<std::vector<std::array<double, 2>>,
+           std::vector<std::array<std::size_t, 2>>>
+construct_rectilinear_grid(
+    const std::array<std::vector<double>, 2> &instances) noexcept {
+
+  const auto generate_points = [&]() {
+    std::vector<std::array<double, 2>> points;
+    points.reserve(instances[0].size() * instances[1].size());
+
+    for (const auto &x : instances[0])
+      for (const auto &y : instances[1])
+        points.push_back({{x, y}});
+
+    return points;
+  };
+
+  const auto step_to_index = [&](const std::size_t step_x,
+                                 const std::size_t step_y) {
+    return step_x * instances[1].size() + step_y;
+  };
+
+  const auto generate_edges = [&]() {
+    std::vector<std::array<std::size_t, 2>> edges;
+    edges.reserve((instances[0].size() - 1) * instances[1].size() +
+                  instances[0].size() * (instances[1].size() - 1));
+
+    std::array<std::size_t, 2> step = {{0, 0}};
+
+    // edges in y
+    for (step[0] = 0; step[0] < instances[0].size(); ++step[0])
+      for (step[1] = 0; step[1] < instances[1].size() - 1; ++step[1])
+        edges.push_back({{step_to_index(step[0], step[1]),
+                          step_to_index(step[0], step[1] + 1)}});
+
+    // edges in x
+    for (step[0] = 0; step[0] < instances[0].size() - 1; ++step[0])
+      for (step[1] = 0; step[1] < instances[1].size(); ++step[1])
+        edges.push_back({{step_to_index(step[0], step[1]),
+                          step_to_index(step[0] + 1, step[1])}});
+    return edges;
+  };
+
+  return {generate_points(), generate_edges()};
+}
+
+template <>
+std::tuple<std::vector<std::array<double, 3>>,
+           std::vector<std::array<std::size_t, 2>>>
+construct_rectilinear_grid(
+    const std::array<std::vector<double>, 3> &instances) noexcept {
+
+  const auto generate_points = [&]() {
+    std::vector<std::array<double, 3>> points;
+    points.reserve(instances[0].size() * instances[1].size() *
+                   instances[2].size());
+
+    for (const auto &x : instances[0])
+      for (const auto &y : instances[1])
+        for (const auto &z : instances[2])
+          points.push_back({{x, y, z}});
+
+    return points;
+  };
+
+  const auto step_to_index = [&](const std::size_t step_x,
+                                 const std::size_t step_y,
+                                 const std::size_t step_z) {
+    return step_x * instances[1].size() * instances[2].size() +
+           step_y * instances[2].size() + step_z;
+  };
+
+  const auto generate_edges = [&]() {
+    std::vector<std::array<std::size_t, 2>> edges;
+    edges.reserve((instances[0].size() - 1) * (instances[1].size() - 1) *
+                      instances[2].size() +
+                  instances[0].size() * (instances[1].size() - 1) *
+                      (instances[2].size() - 1) +
+                  (instances[0].size() - 1) * instances[1].size() *
+                      (instances[2].size() - 1));
+
+    std::array<std::size_t, 3> step = {{0, 0, 0}};
+
+    // edges in z
+    for (step[0] = 0; step[0] < instances[0].size(); ++step[0])
+      for (step[1] = 0; step[1] < instances[1].size(); ++step[1])
+        for (step[2] = 0; step[2] < instances[2].size() - 1; ++step[2])
+          edges.push_back({{step_to_index(step[0], step[1], step[2]),
+                            step_to_index(step[0], step[1], step[2] + 1)}});
+
+    // edges in y
+    for (step[0] = 0; step[0] < instances[0].size(); ++step[0])
+      for (step[1] = 0; step[1] < instances[1].size() - 1; ++step[1])
+        for (step[2] = 0; step[2] < instances[2].size(); ++step[2])
+          edges.push_back({{step_to_index(step[0], step[1], step[2]),
+                            step_to_index(step[0], step[1] + 1, step[2])}});
+
+    // edges in x
+    for (step[0] = 0; step[0] < instances[0].size() - 1; ++step[0])
+      for (step[1] = 0; step[1] < instances[1].size(); ++step[1])
+        for (step[2] = 0; step[2] < instances[2].size(); ++step[2])
+          edges.push_back({{step_to_index(step[0], step[1], step[2]),
+                            step_to_index(step[0] + 1, step[1], step[2])}});
+
+    return edges;
+  };
+
+  return {generate_points(), generate_edges()};
+}
+
+} // namespace meshing
+} // namespace ae108
\ No newline at end of file