diff --git a/.gitignore b/.gitignore
index c182d408..6e2974ee 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,4 +18,5 @@ _generate/
 *.py[cod]
 *.egg-info
 *env*
-.mypy_cache/
\ No newline at end of file
+.mypy_cache/
+.eggs/
\ No newline at end of file
diff --git a/Python/pyWrap.cpp b/Python/pyWrap.cpp
index 512f4d81..3f8c1e93 100644
--- a/Python/pyWrap.cpp
+++ b/Python/pyWrap.cpp
@@ -27,6 +27,7 @@
 #include <lsDetectFeatures.hpp>
 #include <lsDomain.hpp>
 #include <lsExpand.hpp>
+#include <lsExtrude.hpp>
 #include <lsFileFormats.hpp>
 #include <lsFromMesh.hpp>
 #include <lsFromSurfaceMesh.hpp>
@@ -509,6 +510,39 @@ PYBIND11_MODULE(VIENNALS_MODULE_NAME, module) {
       .def("setWidth", &lsExpand<T, D>::setWidth, "Set the width to expand to.")
       .def("apply", &lsExpand<T, D>::apply, "Perform expansion.");
 
+  // lsExtrude
+  // Does not work in current implementation, because one can not import both 2D
+  // and 3D ViennaLS libraries in Python in the same file
+  //   pybind11::class_<lsExtrude<T>, lsSmartPointer<lsExtrude<T>>>(module,
+  //                                                                "lsExtrude")
+  //       // constructors
+  //       .def(pybind11::init(&lsSmartPointer<lsExtrude<T>>::New<>))
+  //       .def(
+  //           pybind11::init(&lsSmartPointer<lsExtrude<T>>::New<
+  //                          lsSmartPointer<lsDomain<T, 2>> &,
+  //                          lsSmartPointer<lsDomain<T, 3>> &, std::array<T,
+  //                          2>, const int,
+  //                          std::array<lsBoundaryConditionEnum<3>, 3>>))
+  //       // methods
+  //       .def("setInputLevelSet", &lsExtrude<T>::setInputLevelSet,
+  //            "Set 2D input Level Set")
+  //       .def("setOutputLevelSet", &lsExtrude<T>::setOutputLevelSet,
+  //            "Set 3D output Level Set")
+  //       .def("setExtent", &lsExtrude<T>::setExtent,
+  //            "Set the extent in the extruded dimension")
+  //       .def("setExtrudeDimension", &lsExtrude<T>::setExtrudeDimension,
+  //            "Set the dimension which should be extruded")
+  //       .def("setBoundaryConditions",
+  //            pybind11::overload_cast<std::array<lsBoundaryConditionEnum<3>,
+  //            3>>(
+  //                &lsExtrude<T>::setBoundaryConditions),
+  //            "Set the boundary conditions in the 3D extruded domain.")
+  //       .def("setBoundaryConditions",
+  //            pybind11::overload_cast<lsBoundaryConditionEnum<3> *>(
+  //                &lsExtrude<T>::setBoundaryConditions),
+  //            "Set the boundary conditions in the 3D extruded domain.")
+  //       .def("apply", &lsExtrude<T>::apply, "Perform extrusion.");
+
   // lsFileFormats
   pybind11::enum_<lsFileFormatEnum>(module, "lsFileFormatEnum")
       .value("VTK_LEGACY", lsFileFormatEnum::VTK_LEGACY)
diff --git a/Tests/Extrude/CMakeLists.txt b/Tests/Extrude/CMakeLists.txt
new file mode 100644
index 00000000..4ed0d7b1
--- /dev/null
+++ b/Tests/Extrude/CMakeLists.txt
@@ -0,0 +1,10 @@
+cmake_minimum_required(VERSION 3.14)
+
+project("Extrude")
+
+add_executable(${PROJECT_NAME} ${PROJECT_NAME}.cpp)
+target_include_directories(${PROJECT_NAME} PUBLIC ${VIENNALS_INCLUDE_DIRS})
+target_link_libraries(${PROJECT_NAME} PRIVATE ${VIENNALS_LIBRARIES})
+
+add_dependencies(buildTests ${PROJECT_NAME})
+add_test(NAME ${PROJECT_NAME} COMMAND $<TARGET_FILE:${PROJECT_NAME}>)
diff --git a/Tests/Extrude/Extrude.cpp b/Tests/Extrude/Extrude.cpp
new file mode 100644
index 00000000..e4dfdd05
--- /dev/null
+++ b/Tests/Extrude/Extrude.cpp
@@ -0,0 +1,77 @@
+#include <iostream>
+
+#include <lsBooleanOperation.hpp>
+#include <lsDomain.hpp>
+#include <lsExtrude.hpp>
+#include <lsMakeGeometry.hpp>
+#include <lsMesh.hpp>
+#include <lsToMesh.hpp>
+#include <lsVTKWriter.hpp>
+
+/**
+  Simple example showing how to use lsExtrude in order
+  to transform a 2D trench into a 3D level set.
+  \example Extrude.cpp
+*/
+
+int main() {
+
+  omp_set_num_threads(4);
+
+  double extent = 15;
+  double gridDelta = 0.5;
+
+  // 2D domain boundaries
+  double bounds[2 * 2] = {-extent, extent, -extent, extent};
+  lsDomain<double, 2>::BoundaryType boundaryCons[2];
+  boundaryCons[0] = lsDomain<double, 2>::BoundaryType::REFLECTIVE_BOUNDARY;
+  boundaryCons[1] = lsDomain<double, 2>::BoundaryType::INFINITE_BOUNDARY;
+
+  auto trench =
+      lsSmartPointer<lsDomain<double, 2>>::New(bounds, boundaryCons, gridDelta);
+
+  {
+    double origin[2] = {0., 0.};
+    double normal[2] = {0., 1.};
+
+    lsMakeGeometry<double, 2>(
+        trench, lsSmartPointer<lsPlane<double, 2>>::New(origin, normal))
+        .apply();
+
+    auto cutOut = lsSmartPointer<lsDomain<double, 2>>::New(bounds, boundaryCons,
+                                                           gridDelta);
+
+    double minPoint[2] = {-5., -5.};
+    double maxPoint[2] = {5., gridDelta};
+
+    lsMakeGeometry<double, 2>(
+        cutOut, lsSmartPointer<lsBox<double, 2>>::New(minPoint, maxPoint))
+        .apply();
+
+    lsBooleanOperation<double, 2>(trench, cutOut,
+                                  lsBooleanOperationEnum::RELATIVE_COMPLEMENT)
+        .apply();
+  }
+
+  {
+    auto mesh = lsSmartPointer<lsMesh<>>::New();
+    lsToMesh<double, 2>(trench, mesh).apply();
+    lsVTKWriter<double>(mesh, "trench_initial.vtp").apply();
+  }
+
+  std::array<double, 2> extrudeExtent = {-5., 5.};
+  std::array<lsBoundaryConditionEnum<3>, 3> boundaryConds{
+      lsBoundaryConditionEnum<3>::REFLECTIVE_BOUNDARY,
+      lsBoundaryConditionEnum<3>::INFINITE_BOUNDARY,
+      lsBoundaryConditionEnum<3>::REFLECTIVE_BOUNDARY};
+  auto trench_3D = lsSmartPointer<lsDomain<double, 3>>::New();
+  lsExtrude<double>(trench, trench_3D, extrudeExtent, 2, boundaryConds).apply();
+
+  {
+    auto mesh = lsSmartPointer<lsMesh<>>::New();
+    lsToMesh<double, 3>(trench_3D, mesh).apply();
+    lsVTKWriter<double>(mesh, "trench_extrude.vtp").apply();
+  }
+
+  return 0;
+}
diff --git a/include/lsExtrude.hpp b/include/lsExtrude.hpp
new file mode 100644
index 00000000..9c83977e
--- /dev/null
+++ b/include/lsExtrude.hpp
@@ -0,0 +1,145 @@
+#ifndef LS_EXTRUDE_HPP
+#define LS_EXTRUDE_HPP
+
+#include <lsDomain.hpp>
+#include <lsFromSurfaceMesh.hpp>
+#include <lsToSurfaceMesh.hpp>
+
+/// Extrudes a 2D Level Set into a 3D domain. The axis in which should be
+/// extruded can be set and boundary conditions in the 3D domain must be
+/// specified.
+template <class T> class lsExtrude {
+  lsSmartPointer<lsDomain<T, 2>> inputLevelSet = nullptr;
+  lsSmartPointer<lsDomain<T, 3>> outputLevelSet = nullptr;
+  std::array<T, 2> extent = {0., 0.};
+  int extrudeDim = 0;
+  std::array<lsBoundaryConditionEnum<3>, 3> boundaryConds;
+
+public:
+  lsExtrude() {}
+  lsExtrude(lsSmartPointer<lsDomain<T, 2>> passedInputLS,
+            lsSmartPointer<lsDomain<T, 3>> passedOutputLS,
+            std::array<T, 2> passedExtent, const int passedExtrudeDim,
+            std::array<lsBoundaryConditionEnum<3>, 3> passedBoundaryConds)
+      : inputLevelSet(passedInputLS), outputLevelSet(passedOutputLS),
+        extent(passedExtent), extrudeDim(passedExtrudeDim),
+        boundaryConds(passedBoundaryConds) {}
+
+  void setInputLevelSet(lsSmartPointer<lsDomain<T, 2>> passedInputLS) {
+    inputLevelSet = passedInputLS;
+  }
+
+  // The 3D output LS will be overwritten by the extruded LS
+  void setOutputLevelSet(lsSmartPointer<lsDomain<T, 3>> &passedOutputLS) {
+    outputLevelSet = passedOutputLS;
+  }
+
+  // Set the min and max extent in the extruded dimension
+  void setExtent(std::array<T, 2> passedExtent) { extent = passedExtent; }
+
+  // Set which index of the added dimension (x: 0, y: 1, z: 2)
+  void setExtrudeDimension(const int passedExtrudeDim) {
+    extrudeDim = passedExtrudeDim;
+  }
+
+  void setBoundaryConditions(
+      std::array<lsBoundaryConditionEnum<3>, 3> passedBoundaryConds) {
+    boundaryConds = passedBoundaryConds;
+  }
+
+  void
+  setBoundaryConditions(lsBoundaryConditionEnum<3> passedBoundaryConds[3]) {
+    for (int i = 0; i < 3; i++)
+      boundaryConds[i] = passedBoundaryConds[i];
+  }
+
+  void apply() {
+    if (inputLevelSet == nullptr) {
+      lsMessage::getInstance()
+          .addWarning(
+              "No input Level Set supplied to lsExtrude! Not converting.")
+          .print();
+    }
+    if (outputLevelSet == nullptr) {
+      lsMessage::getInstance()
+          .addWarning(
+              "No output Level Set supplied to lsExtrude! Not converting.")
+          .print();
+      return;
+    }
+
+    // x and y of the input LS get transformed to these indices
+    const auto extrudeDims = getExtrudeDims();
+
+    // create new domain based on 2D extent
+    {
+      const T gridDelta = inputLevelSet->getGrid().getGridDelta();
+      auto minBounds = inputLevelSet->getGrid().getMinBounds();
+      auto maxBounds = inputLevelSet->getGrid().getMaxBounds();
+
+      double domainBounds[2 * 3];
+      domainBounds[2 * extrudeDim] = extent[0];
+      domainBounds[2 * extrudeDim + 1] = extent[1];
+      domainBounds[2 * extrudeDims[0]] = gridDelta * minBounds[0];
+      domainBounds[2 * extrudeDims[0] + 1] = gridDelta * maxBounds[0];
+      domainBounds[2 * extrudeDims[1]] = gridDelta * minBounds[1];
+      domainBounds[2 * extrudeDims[1] + 1] = gridDelta * maxBounds[1];
+
+      auto tmpLevelSet = lsSmartPointer<lsDomain<T, 3>>::New(
+          domainBounds, boundaryConds.data(), gridDelta);
+      outputLevelSet->deepCopy(tmpLevelSet);
+    }
+
+    auto surface = lsSmartPointer<lsMesh<T>>::New();
+    lsToSurfaceMesh<T, 2>(inputLevelSet, surface).apply();
+
+    auto &lines = surface->template getElements<2>();
+    auto &nodes = surface->getNodes();
+    const unsigned numNodes = nodes.size();
+
+    // add new nodes shifted by the extent
+    for (unsigned i = 0; i < numNodes; i++) {
+      nodes[i][extrudeDims[1]] = nodes[i][1];
+      nodes[i][extrudeDims[0]] = nodes[i][0];
+      nodes[i][extrudeDim] = extent[1];
+
+      nodes.push_back(nodes[i]);
+
+      nodes[i][extrudeDim] = extent[0];
+    }
+
+    // add triangles in places of lines
+    for (unsigned i = 0; i < lines.size(); i++) {
+      std::array<unsigned, 3> triangle = {lines[i][1], lines[i][0],
+                                          lines[i][0] + numNodes};
+      if (extrudeDim == 1) {
+        std::swap(triangle[0], triangle[2]);
+      }
+      surface->insertNextTriangle(triangle);
+      triangle[0] = lines[i][0] + numNodes;
+      triangle[1] = lines[i][1] + numNodes;
+      triangle[2] = lines[i][1];
+      if (extrudeDim == 1) {
+        std::swap(triangle[0], triangle[2]);
+      }
+      surface->insertNextTriangle(triangle);
+    }
+    surface->template getElements<2>().clear(); // remove lines
+
+    lsFromSurfaceMesh<T, 3>(outputLevelSet, surface).apply();
+  }
+
+private:
+  inline std::array<unsigned, 2> getExtrudeDims() const {
+    assert(extrudeDim < 3);
+    if (extrudeDim == 0) {
+      return {1, 2};
+    } else if (extrudeDim == 1) {
+      return {0, 2};
+    } else {
+      return {0, 1};
+    }
+  }
+};
+
+#endif // LS_EXTRUDE_HPP
\ No newline at end of file