diff --git a/src/csg/CMakeLists.txt b/src/csg/CMakeLists.txt
index 0efaa7b49a67ea2f6f632165c3096af6afe23913..ed0bd36038694bf9f1edd1f82b1560abef65dae3 100644
--- a/src/csg/CMakeLists.txt
+++ b/src/csg/CMakeLists.txt
@@ -7,6 +7,7 @@ add_library(csg
   levelset_3d.cpp
   matrix_functions.cpp
   parser.cpp
+  cgal_helper.cpp
   )
 
 target_link_libraries(csg PRIVATE
diff --git a/src/csg/cgal_helper.cpp b/src/csg/cgal_helper.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..72378cc3869651a1a5a28844c4e45241a4706c62
--- /dev/null
+++ b/src/csg/cgal_helper.cpp
@@ -0,0 +1,90 @@
+#include "csg_cgal_helper.hpp"
+
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+#include <CGAL/Polyhedron_incremental_builder_3.h>
+#include <CGAL/Side_of_triangle_mesh.h>
+#include <CGAL/algorithm.h>
+#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
+
+namespace cgal_helper {
+
+namespace {
+
+typedef Polyhedron::HalfedgeDS HalfedgeDS;
+typedef typename HalfedgeDS::Vertex Vertex;
+typedef typename Vertex::Point Point;
+typedef CGAL::AABB_face_graph_triangle_primitive<Polyhedron> Primitive;
+typedef CGAL::AABB_traits<Kernel, Primitive> Traits;
+typedef CGAL::AABB_tree<Traits> Tree;
+typedef CGAL::Side_of_triangle_mesh<Polyhedron, Kernel> Point_inside;
+
+// A modifier creating a polyhedron with the incremental builder.
+template <class HDS> class PolyhedronBuilder : public CGAL::Modifier_base<HDS> {
+private:
+  const std::vector<std::tuple<double, double, double>> &m_points;
+  const std::vector<std::vector<unsigned int>> &m_faces;
+
+public:
+  PolyhedronBuilder(
+      const std::vector<std::tuple<double, double, double>> &points,
+      const std::vector<std::vector<unsigned int>> &faces)
+      : m_points(points), m_faces(faces) {}
+
+  void operator()(HDS &hds) {
+    CGAL::Polyhedron_incremental_builder_3<HDS> B(hds, true);
+    B.begin_surface(m_points.size(), m_faces.size());
+
+    // Add all the vertices first
+    for (const auto &p : m_points) {
+      auto [px, py, pz] = p;
+      B.add_vertex(Point(px, py, pz));
+    }
+
+    // Add facets next
+    for (const auto &face : m_faces) {
+      B.begin_facet();
+      for (const auto &p_index : face) {
+        B.add_vertex_to_facet(p_index);
+      }
+      B.end_facet();
+    }
+
+    B.end_surface();
+  }
+};
+
+} // namespace
+
+Polyhedron
+create_polyhedron(const std::vector<std::tuple<double, double, double>> &points,
+                  const std::vector<std::vector<unsigned int>> &faces) {
+  Polyhedron p;
+
+  // Build incrementally
+  PolyhedronBuilder<HalfedgeDS> bp(points, faces);
+  p.delegate(bp);
+  CGAL_assertion(p.is_valid());
+
+  // Triangulate faces - needed for levelset
+  CGAL::Polygon_mesh_processing::triangulate_faces(p);
+  CGAL_assertion(p.is_valid());
+
+  return p;
+}
+
+bool inside(const Polyhedron &polyhedron, double xx, double yy, double zz) {
+  Kernel::Point_3 pt(xx, yy, zz);
+  // Construct AABB tree with a KdTree
+  Tree tree(faces(polyhedron).first, faces(polyhedron).second, polyhedron);
+  tree.accelerate_distance_queries();
+  // Initialize the point-in-polyhedron tester
+  Point_inside inside_tester(tree);
+
+  // Determine the side and return true if inside!
+  return inside_tester(pt) == CGAL::ON_BOUNDED_SIDE;
+}
+
+} // namespace cgal_helper
diff --git a/src/csg/levelset_3d.cpp b/src/csg/levelset_3d.cpp
index 1eee5aac2f79b256896eb6a98e088d4e9de5d487..3e7ce17be26ca97f27fc760b1a99626cf8f8ab17 100644
--- a/src/csg/levelset_3d.cpp
+++ b/src/csg/levelset_3d.cpp
@@ -1,3 +1,4 @@
+#include "csg_cgal_helper.hpp"
 #include "csg_types.hpp"
 
 #include <algorithm>
@@ -27,6 +28,7 @@ double signed_distance_3d(const Sphere &, double, double, double);
 double signed_distance_3d(const Type3D &, double, double, double);
 double signed_distance_3d(const LinearExtrude &, double, double, double);
 double signed_distance_3d(const RotateExtrude &, double, double, double);
+double signed_distance_3d(const Polyhedron &, double, double, double);
 
 double signed_distance_2d(const Union2D &, double, double);
 
@@ -149,6 +151,13 @@ double signed_distance_3d(const RotateExtrude &rot_ext, double xx, double yy,
   return signed_distance_2d(rot_ext.group, XX, YY);
 }
 
+double signed_distance_3d(const Polyhedron &polyhedron, double xx, double yy,
+                          double zz) {
+  // TODO: support signed distance instead of -1.0/1.0
+  return cgal_helper::inside(polyhedron.cgal_polyhedron(), xx, yy, zz) ? 1.0
+                                                                       : -1.0;
+}
+
 double signed_distance_3d(const Type3D &obj, double xx, double yy, double zz) {
 
   return std::visit(
diff --git a/src/csg/meson.build b/src/csg/meson.build
index c9f7140ed7259a289720df49a5fceb3109da3885..dafffc17315de7500a66a35b91efa38696722698 100644
--- a/src/csg/meson.build
+++ b/src/csg/meson.build
@@ -5,6 +5,7 @@ lib_csg_parser = static_library(
   'levelset_2d.cpp',
   'parser.cpp',
   'matrix_functions.cpp',
+  'cgal_helper.cpp',
   include_directories: parser_inc,
   dependencies: [pegtl, cgal],
   install : true)
diff --git a/src/csg/parser.cpp b/src/csg/parser.cpp
index 3e1429bd0d49a76a302b3996514f60640080ce1f..ec6209a29bd4e81119c4d51620622cf7c3fc490e 100644
--- a/src/csg/parser.cpp
+++ b/src/csg/parser.cpp
@@ -13,7 +13,8 @@ using namespace tao::pegtl;
 
 namespace {
 
-using Attr = std::variant<double, bool, std::vector<double>, std::string>;
+using Attr = std::variant<double, bool, std::vector<double>, std::string,
+                          std::vector<std::vector<double>>>;
 using AttrMap = std::map<std::string, Attr>;
 
 struct parser_state {
@@ -83,16 +84,18 @@ struct vector : seq<L_ARR, list<vector_cell, one<','>>, R_ARR> {};
 
 struct vector_attr : vector {};
 
-struct value
-    : sor<padded_double, boolean_literal, string_literal, vector_attr> {};
+struct row : vector {};
 
-struct keyval : seq<pad<name, space>, string<'='>, pad<value, space>> {};
+struct matrix : seq<L_ARR, list<row, one<','>>, R_ARR> {};
 
-struct attr_list : list<keyval, one<','>> {};
+struct matrix_attr : matrix {};
 
-struct row : vector {};
+struct value : sor<padded_double, boolean_literal, string_literal, vector_attr,
+                   matrix_attr> {};
 
-struct matrix : seq<L_ARR, list<row, one<','>>, R_ARR> {};
+struct keyval : seq<pad<name, space>, string<'='>, pad<value, space>> {};
+
+struct attr_list : list<keyval, one<','>> {};
 
 struct sphere
     : seq<string<'s', 'p', 'h', 'e', 'r', 'e'>, L_FUN, attr_list, R_FUN> {};
@@ -102,6 +105,10 @@ struct cube : seq<string<'c', 'u', 'b', 'e'>, L_FUN, attr_list, R_FUN> {};
 struct cylinder : seq<string<'c', 'y', 'l', 'i', 'n', 'd', 'e', 'r'>, L_FUN,
                       attr_list, R_FUN> {};
 
+struct polyhedron
+    : seq<string<'p', 'o', 'l', 'y', 'h', 'e', 'd', 'r', 'o', 'n'>, L_FUN,
+          attr_list, R_FUN> {};
+
 struct circle
     : seq<string<'c', 'i', 'r', 'c', 'l', 'e'>, L_FUN, attr_list, R_FUN> {};
 
@@ -111,7 +118,8 @@ struct square
 template <Dimension dim> struct shape;
 
 template <>
-struct shape<Dimension::D3> : seq<sor<cube, cylinder, sphere>, opt<S_CLN>> {};
+struct shape<Dimension::D3>
+    : seq<sor<cube, cylinder, sphere, polyhedron>, opt<S_CLN>> {};
 
 template <>
 struct shape<Dimension::D2> : seq<sor<circle, square>, opt<S_CLN>> {};
@@ -562,6 +570,24 @@ template <> struct action<square> {
   }
 };
 
+template <> struct action<polyhedron> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    std::stringstream ss(in.string());
+    auto &curr_attr = st.curr_attrs.back();
+
+    auto points =
+        std::get<std::vector<std::vector<double>>>(curr_attr["points"]);
+    auto faces = std::get<std::vector<std::vector<double>>>(curr_attr["faces"]);
+
+    csg::Polyhedron polyh(points, faces);
+    polyh.name = get_name(curr_attr);
+
+    st.current_3d_objs.back().push_back(polyh);
+    st.curr_attrs.pop_back();
+  }
+};
+
 template <> struct action<name> {
   template <typename Input>
   static void apply(const Input &in, parser_state &st) {
@@ -615,6 +641,15 @@ template <> struct action<vector_attr> {
   }
 };
 
+template <> struct action<matrix_attr> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    std::stringstream ss(in.string());
+    st.curr_attr[st.current_name] = st.current_matrix;
+    st.current_matrix.clear();
+  }
+};
+
 std::optional<parser_state> do_parse(std::string str) {
   parser_state st;
   std::vector<csg::Type3D> new_3d_group;
diff --git a/src/csg/tests/CMakeLists.txt b/src/csg/tests/CMakeLists.txt
index f3fcb5536414e722da53f3c3d599a7ba5d15003b..a5ee6f4baf3bda05e4ed1a1f4c5712782c6c16d7 100644
--- a/src/csg/tests/CMakeLists.txt
+++ b/src/csg/tests/CMakeLists.txt
@@ -23,6 +23,7 @@ target_include_directories(unit_tests_csg
 target_link_libraries(unit_tests_csg
   csg
   CONAN_PKG::catch2
+  CONAN_PKG::cgal
   )
 
 catch_discover_tests(unit_tests_csg)
diff --git a/src/csg/tests/levelset/primitives.t.cpp b/src/csg/tests/levelset/primitives.t.cpp
index b2b3b0152f9aa08138c9f5219ceab694f60426a1..b950f8285275bf2f9fb6ff1f5df8489e766ef7de 100644
--- a/src/csg/tests/levelset/primitives.t.cpp
+++ b/src/csg/tests/levelset/primitives.t.cpp
@@ -545,4 +545,44 @@ TEST_CASE("Cone", "[Levelset Primitives]") {
   }
 }
 
+TEST_CASE("Polyhedron", "[Levelset Primitives]") {
+  double Lx = 10.0, Ly = 7.0, Lz = 5.0;
+  // A 10 x 7 x 5 cuboid in positive quandrant formed using polyhedron
+  csg::Polyhedron my_cub({{0, 0, 0},
+                          {Lx, 0, 0},
+                          {Lx, Ly, 0},
+                          {0, Ly, 0},
+                          {0, 0, Lz},
+                          {Lx, 0, Lz},
+                          {Lx, Ly, Lz},
+                          {0, Ly, Lz}},
+                         {{0, 1, 2, 3},
+                          {4, 5, 1, 0},
+                          {7, 6, 5, 4},
+                          {5, 6, 2, 1},
+                          {6, 7, 3, 2},
+                          {7, 4, 0, 3}});
+
+  auto my_tree = std::make_shared<csg::Tree>();
+  my_tree->top.objs.push_back(my_cub);
+  csg::CsgIF my_levelset(my_tree);
+
+  SECTION("Outside") {
+    CHECK_FALSE(0 < my_levelset(1.01 * Lx, 0.1 * Ly, 0.1 * Lz));
+    CHECK_FALSE(0 < my_levelset(-0.01 * Lx, 0.5 * Ly, 0.5 * Lz));
+    CHECK_FALSE(0 < my_levelset(0.3 * Lx, 1.01 * Ly, 0.2 * Lz));
+    CHECK_FALSE(0 < my_levelset(0.7 * Lx, -0.01 * Ly, 0.7 * Lz));
+    CHECK_FALSE(0 < my_levelset(0.3 * Lx, 0.9 * Ly, 1.01 * Lz));
+    CHECK_FALSE(0 < my_levelset(0.7 * Lx, 0.1 * Ly, -0.01 * Lz));
+  }
+  SECTION("Inside") {
+    CHECK(0 < my_levelset(0.99 * Lx, 0.1 * Ly, 0.1 * Lz));
+    CHECK(0 < my_levelset(0.01 * Lx, 0.5 * Ly, 0.5 * Lz));
+    CHECK(0 < my_levelset(0.3 * Lx, 0.99 * Ly, 0.2 * Lz));
+    CHECK(0 < my_levelset(0.7 * Lx, 0.01 * Ly, 0.7 * Lz));
+    CHECK(0 < my_levelset(0.3 * Lx, 0.9 * Ly, 0.99 * Lz));
+    CHECK(0 < my_levelset(0.7 * Lx, 0.1 * Ly, 0.01 * Lz));
+  }
+}
+
 } // namespace
diff --git a/src/csg/tests/parser/primitives.t.cpp b/src/csg/tests/parser/primitives.t.cpp
index ca2175cd07cb8ae83f6f22a5f0d6dedf538bc427..7847718ffdf6c44f23d6e30189b174a23e2a39a8 100644
--- a/src/csg/tests/parser/primitives.t.cpp
+++ b/src/csg/tests/parser/primitives.t.cpp
@@ -54,4 +54,52 @@ sphere(r = 10)
   CHECK(sph.radius == 10);
 }
 
-// TODO: polyhedron()
+TEST_CASE("square pyramid polyhedron", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+polyhedron(
+points = [[10, 10, 0], [10, -10, 0], [-10, -10, 0], [-10, 10, 0], [0, 0, 10]], 
+faces = [[0, 1, 4], [1, 2, 4], [2, 3, 4], [3, 0, 4], [1, 0, 3], [2, 1, 3]], 
+$name="my_polyhedron");
+)");
+  auto polyh = std::get<csg::Polyhedron>(st.top.objs.at(0));
+  CHECK(polyh.name == "my_polyhedron");
+  CHECK(polyh.cgal_polyhedron().size_of_vertices() == 5);
+  CHECK(polyh.cgal_polyhedron().size_of_facets() == 6);
+}
+
+TEST_CASE("bad vs good polyhedron openscad example", "[csg]") {
+  // https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/Primitive_Solids#polyhedron
+  CHECK_THROWS(csg::parse_csg(R"(
+polyhedron(
+points = [[0, -10, 60], [0, 10, 60], [0, 10, 0], 
+[0, -10, 0], [60, -10, 60], [60, 10, 60], [10, -10, 50], 
+[10, 10, 50], [10, 10, 30], [10, -10, 30], [30, -10, 50], 
+[30, 10, 50]], 
+
+faces = [[0, 2, 3], [0, 1, 2], [0, 4, 5], [0, 5, 1], 
+[5, 4, 2], [2, 4, 3], [6, 8, 9], [6, 7, 8], [6, 10, 11], 
+[6, 11, 7], [10, 8, 11], [10, 9, 8], [0, 3, 9], [9, 0, 6], 
+[10, 6, 0], [0, 4, 10], [3, 9, 10], [3, 10, 4], [1, 7, 11], 
+[1, 11, 5], [1, 7, 8], [1, 8, 2], [2, 8, 11], [2, 11, 5]], 
+
+convexity = 1);
+)"));
+
+  auto st = *csg::parse_csg(R"(
+polyhedron(
+points = [[0, -10, 60], [0, 10, 60], [0, 10, 0], 
+[0, -10, 0], [60, -10, 60], [60, 10, 60], [10, -10, 50], 
+[10, 10, 50], [10, 10, 30], [10, -10, 30], [30, -10, 50], 
+[30, 10, 50]], 
+
+faces = [[0, 3, 2], [0, 2, 1], [4, 0, 5], [5, 0, 1], 
+[5, 2, 4], [4, 2, 3], [6, 8, 9], [6, 7, 8], [6, 10, 11], 
+[6, 11, 7], [10, 8, 11], [10, 9, 8], [3, 0, 9], [9, 0, 6], 
+[10, 6, 0], [0, 4, 10], [3, 9, 10], [3, 10, 4], [1, 7, 11], 
+[1, 11, 5], [1, 8, 7], [2, 8, 1], [8, 2, 11], [5, 11, 2]], 
+
+convexity = 1);
+)");
+  auto polyh = std::get<csg::Polyhedron>(st.top.objs.at(0));
+  CHECK(polyh.cgal_polyhedron().size_of_vertices() == 12);
+}
diff --git a/src/csg_cgal_helper.hpp b/src/csg_cgal_helper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..52d9aed823128607a308aec6ab7f31445485115a
--- /dev/null
+++ b/src/csg_cgal_helper.hpp
@@ -0,0 +1,21 @@
+#ifndef CGAL_HELPER_H_
+#define CGAL_HELPER_H_
+
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Simple_cartesian.h>
+
+namespace cgal_helper {
+
+typedef CGAL::Simple_cartesian<double> Kernel;
+typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
+
+Polyhedron
+create_polyhedron(const std::vector<std::tuple<double, double, double>> &points,
+                  const std::vector<std::vector<unsigned int>> &faces);
+
+bool inside(const Polyhedron &polyhedron, double xx, double yy,
+                       double zz);
+
+} // namespace cgal_helper
+
+#endif
diff --git a/src/csg_types.hpp b/src/csg_types.hpp
index 5077100942fe8e0a4c5185230a3e51804334700a..0ffe415c64b41b7e889c752a4b63fa6e59f2c186 100644
--- a/src/csg_types.hpp
+++ b/src/csg_types.hpp
@@ -1,6 +1,7 @@
 #ifndef CSG_TYPES_H_
 #define CSG_TYPES_H_
 
+#include "csg_cgal_helper.hpp"
 #include "csg_matrix_functions.hpp"
 
 #include <array>
@@ -49,6 +50,31 @@ struct Cone {
   bool center;
 };
 
+struct Polyhedron {
+private:
+  std::vector<std::tuple<double, double, double>> m_points;
+  std::vector<std::vector<unsigned int>> m_faces;
+  cgal_helper::Polyhedron m_cgal_polyhedron;
+
+public:
+  std::optional<std::string> name;
+  Polyhedron(const std::vector<std::vector<double>> &points,
+             const std::vector<std::vector<double>> &faces) {
+    for (const auto &pt : points) {
+      assert(pt.size() == 3);
+      m_points.push_back({pt[0], pt[1], pt[2]});
+    }
+    for (const auto &f : faces) {
+      m_faces.push_back(std::vector<unsigned int>(f.begin(), f.end()));
+    }
+    m_cgal_polyhedron = cgal_helper::create_polyhedron(m_points, m_faces);
+  }
+
+  const cgal_helper::Polyhedron &cgal_polyhedron() const {
+    return m_cgal_polyhedron;
+  }
+};
+
 enum Dimension { D2, D3 };
 
 template <Dimension dim> struct Mulmatrix;
@@ -68,10 +94,10 @@ template <> struct TypeHelper<Dimension::D2> {
 };
 
 template <> struct TypeHelper<Dimension::D3> {
-  using Type =
-      std::variant<Sphere, Cube, Cylinder, Cone, Union<Dimension::D3>,
-                   Intersection<Dimension::D3>, Difference<Dimension::D3>,
-                   Mulmatrix<Dimension::D3>, LinearExtrude, RotateExtrude>;
+  using Type = std::variant<Sphere, Cube, Cylinder, Cone, Union<Dimension::D3>,
+                            Intersection<Dimension::D3>,
+                            Difference<Dimension::D3>, Mulmatrix<Dimension::D3>,
+                            LinearExtrude, RotateExtrude, Polyhedron>;
 };
 
 template <Dimension dim> struct Union {