diff --git a/CMakeLists.txt b/CMakeLists.txt
index cdfb12c0e74796918f0af197d5ca27c0ed6c7ecf..953f6ae254ed11e63a1c949ed2b46ed3d1183591 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,39 +1,49 @@
-cmake_minimum_required(VERSION 3.14)
-
-set(CMAKE_CXX_STANDARD 17)
-set(CMAKE_CXX_STANDARD_REQUIRED ON)
-set(CMAKE_CXX_EXTENSIONS OFF)
-
-project(MFIX-Parser
-   DESCRIPTION  "A parser for MFIX-Exa input files"
-   HOMEPAGE_URL "https://mfix.netl.doe.gov/gitlab/exa/mfix-parser"
-   LANGUAGES    CXX
-   )
-
-set( CMAKE_EXPORT_COMPILE_COMMANDS ON )
-set(USE_CCACHE "")
-find_program(CCACHE_FOUND ccache)
-if(CCACHE_FOUND)
-  set( CMAKE_CXX_COMPILER_LAUNCHER ccache )
+################################################################################
+# CSG Parser
+################################################################################
+add_library(csg)
+
+target_sources(csg PRIVATE
+  src/csg/impl/parser.cpp
+  src/csg/impl/levelset.cpp
+  src/csg/impl/csg.cpp
+  )
+
+target_include_directories(csg PRIVATE
+  ${CMAKE_CURRENT_SOURCE_DIR}/subprojects/PEGTL/include)
+
+target_link_libraries(csg PRIVATE stdc++fs)
+
+target_compile_features(csg PRIVATE cxx_std_17)
+
+SET(catch2_dir ${CMAKE_CURRENT_SOURCE_DIR}/subprojects/Catch2)
+
+if(NOT EXISTS "${catch2_dir}/CMakeLists.txt")
+  message(FATAL_ERROR "File ${catch2_dir}/CMakeLists.txt does not exist!  The submodules were not downloaded! GIT_SUBMODULE was turned off or failed. Please update submodules and try again.")
 endif()
 
-add_subdirectory(src)
-add_subdirectory(tests)
+add_subdirectory(${catch2_dir} ${CMAKE_CURRENT_BINARY_DIR}/catch_build)
+list(APPEND CMAKE_MODULE_PATH "${catch2_dir}/contrib/")
+
+include(CTest)
+include(Catch)
+
 
 find_program(CLANG_FMT NAMES clang-format-9 clang-format-8 clang-format-7 clang-format)
 
 if(CLANG_FMT)
-  get_target_property(_psrcs unit_tests SOURCES)
+  get_target_property(_csg_srcs csg SOURCES)
+  file(GLOB _csg_incs ${CMAKE_CURRENT_SOURCE_DIR}/*.hpp ${CMAKE_CURRENT_SOURCE_DIR}/impl/*.hpp)
+
+  add_custom_target(fmt_csg
+    COMMAND ${CLANG_FMT} -i ${_csg_srcs} ${_csg_incs}
+    WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+    )
   add_custom_target(fmt_test
-    COMMAND ${CLANG_FMT} -i ${_psrcs}
+    COMMAND ${CLANG_FMT} -i ${_ut_srcs} ${_ft_srcs}
     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/tests
     )
-  file(GLOB _includes ${PROJECT_SOURCE_DIR}/src/*.hpp)
 
   add_custom_target(fmt)
-  add_dependencies(fmt fmt_test)
+  add_dependencies(fmt fmt_test fmt_csg)
 endif()
-
-set(EXENAME "mfix-parser")
-add_executable(${EXENAME} src/main.cpp)
-target_link_libraries(${EXENAME} parser)
diff --git a/meson.build b/meson.build
index 0b8213fcf1d8e553d38bad757bcc33f0a8ab90d0..7c6e889c4b5ef8f2999bd0e19d1a3ffccb62ef19 100644
--- a/meson.build
+++ b/meson.build
@@ -4,7 +4,19 @@ project('mfix-parser', 'cpp',
 
 tao_inc = include_directories('subprojects/PEGTL/include')
 catch2_inc = include_directories('subprojects/Catch2/single_include')
-parser_inc = include_directories('src')
+parser_inc = include_directories(
+  'src/csg',
+  'src/csg/impl',
+  'src/inputs')
 
-subdir('src')
-subdir('tests')
+subdir('src/csg')
+subdir('src/inputs')
+
+executable(
+  'mfix-parser',
+  'src/main.cpp',
+  include_directories: parser_inc,
+  link_with: [
+    lib_csg_parser,
+    lib_inputs_parser,
+  ])
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
deleted file mode 100644
index b83493c6c62fc7fd428b877ef93ef04202bda6d9..0000000000000000000000000000000000000000
--- a/src/CMakeLists.txt
+++ /dev/null
@@ -1,13 +0,0 @@
-################################################################################
-# inputs Parser
-################################################################################
-
-add_library(parser
-  parser.cpp
-  solver.cpp
-  )
-
-target_include_directories(parser PRIVATE
-  ${CMAKE_SOURCE_DIR}/subprojects/PEGTL/include)
-
-target_link_libraries(parser PRIVATE stdc++fs)
diff --git a/src/csg/csg.hpp b/src/csg/csg.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1862b8305afc484ea004a513157274cc937b4454
--- /dev/null
+++ b/src/csg/csg.hpp
@@ -0,0 +1,52 @@
+#ifndef CSG_H_
+#define CSG_H_
+
+#if defined(CSG_USE_GPU) && !defined(CSG_GPU_HOST_DEVICE)
+#define CSG_GPU_HOST_DEVICE __host__ __device__
+#else
+#define CSG_GPU_HOST_DEVICE
+#endif
+
+#include <array>
+#include <memory>
+#include <tuple>
+#include <vector>
+
+namespace csg {
+
+struct GPUable {};
+template <class D, class Enable = void> struct IsGPUable : std::false_type {};
+template <class D>
+struct IsGPUable<
+    D, typename std::enable_if<std::is_base_of<GPUable, D>::value>::type>
+    : std::true_type {};
+
+struct Tree;
+std::shared_ptr<Tree> parse_csg(std::string);
+
+class CsgIF : GPUable {
+public:
+  CsgIF(std::shared_ptr<Tree> a_state);
+  CsgIF(const CsgIF &rhs);
+  CsgIF(CsgIF &&rhs) noexcept;
+  CsgIF &operator=(const CsgIF &rhs) = delete;
+  CsgIF &operator=(CsgIF &&rhs) = delete;
+  ~CsgIF();
+
+  CSG_GPU_HOST_DEVICE
+  double operator()(double xx, double yy, double zz) const noexcept;
+
+  inline double operator()(const std::array<double, 3> &p) const noexcept {
+    return this->operator()(p[0], p[1], p[2]);
+  }
+
+private:
+  std::shared_ptr<Tree> m_state;
+
+  class Impl;
+  std::unique_ptr<Impl> m_pimpl;
+};
+
+} // namespace csg
+
+#endif
diff --git a/src/csg/impl/csg.cpp b/src/csg/impl/csg.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1cb62eee5132ad1b2969f1f48d659fbb122fb937
--- /dev/null
+++ b/src/csg/impl/csg.cpp
@@ -0,0 +1,28 @@
+#include "../csg.hpp"
+
+#include "csg_types.hpp"
+
+namespace csg {
+
+class CsgIF::Impl {
+public:
+  double call_signed_distance(const std::shared_ptr<Tree> &a_tree, double xx,
+                              double yy, double zz) const {
+    return signed_distance(a_tree->top, xx, yy, zz);
+  }
+};
+
+CsgIF::CsgIF(const CsgIF &rhs)
+    : m_state(rhs.m_state), m_pimpl(std::make_unique<Impl>(*rhs.m_pimpl)) {}
+
+CsgIF::CsgIF(std::shared_ptr<Tree> a_state)
+    : m_state(a_state), m_pimpl(std::make_unique<Impl>()) {}
+
+CsgIF::CsgIF(CsgIF &&rhs) noexcept = default;
+CsgIF::~CsgIF() = default;
+
+double CsgIF::operator()(double xx, double yy, double zz) const noexcept {
+  return m_pimpl->call_signed_distance(m_state, xx, yy, zz);
+}
+
+} // namespace csg
diff --git a/src/csg/impl/csg_types.hpp b/src/csg/impl/csg_types.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e70f97bf6adb93e1b9122e1453d8fd1fd19e24be
--- /dev/null
+++ b/src/csg/impl/csg_types.hpp
@@ -0,0 +1,69 @@
+#ifndef CSG_TYPES_H_
+#define CSG_TYPES_H_
+
+#include <memory>
+#include <optional>
+#include <tuple>
+#include <variant>
+#include <vector>
+
+namespace csg {
+
+struct Sphere {
+  double radius;
+};
+
+struct Cube {
+  std::tuple<double, double, double> size;
+  bool center;
+};
+
+struct Cylinder {
+  double radius;
+  double height;
+  bool center;
+};
+
+struct Cone {
+  double radius1;
+  double radius2;
+  double height;
+  bool center;
+};
+
+struct Mulmatrix;
+struct Union;
+struct Intersection;
+struct Difference;
+
+using Type = std::variant<Sphere, Cube, Cylinder, Cone, Union, Intersection,
+                          Difference, Mulmatrix>;
+
+struct Union {
+  std::vector<Type> objs;
+};
+
+struct Intersection {
+  std::vector<Type> objs;
+};
+
+struct Difference {
+  std::shared_ptr<Type> first_obj;
+  Union next_objs;
+};
+
+struct Mulmatrix {
+  std::array<std::array<double, 3>, 3> rotation;
+  std::array<double, 3> translation;
+  Union group;
+};
+
+struct Tree {
+  Union top;
+};
+
+double signed_distance(const Union &, double, double, double);
+
+} // namespace csg
+
+#endif
diff --git a/src/csg/impl/levelset.cpp b/src/csg/impl/levelset.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0f23278424c5668f955fcc00265a18f256545b3c
--- /dev/null
+++ b/src/csg/impl/levelset.cpp
@@ -0,0 +1,135 @@
+#include "csg_types.hpp"
+
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <limits>
+
+namespace {
+
+const double EXTERNAL_FLOW = -1.0;
+
+template <class... Ts> struct overloaded : Ts... { using Ts::operator()...; };
+template <class... Ts>
+// clang-format off
+overloaded(Ts...) -> overloaded<Ts...>; // not needed as of C++20
+// clang-format on
+
+} // namespace
+
+namespace csg {
+
+double signed_distance(const Cone &, double, double, double);
+double signed_distance(const Cube &, double, double, double);
+double signed_distance(const Cylinder &, double, double, double);
+double signed_distance(const Difference &, double, double, double);
+double signed_distance(const Intersection &, double, double, double);
+double signed_distance(const Mulmatrix &, double, double, double);
+double signed_distance(const Sphere &, double, double, double);
+double signed_distance(const Type &, double, double, double);
+
+double signed_distance(const Union &group, double xx, double yy, double zz) {
+  auto sdist = -std::numeric_limits<double>::max();
+  for (const auto &member : group.objs) {
+    auto sd = signed_distance(member, xx, yy, zz);
+    sdist = std::max(sdist, sd);
+  };
+
+  return sdist;
+}
+
+double signed_distance(const Intersection &group, double xx, double yy,
+                       double zz) {
+  auto sdist = std::numeric_limits<double>::max();
+  for (const auto &member : group.objs) {
+    auto sd = signed_distance(member, xx, yy, zz);
+    sdist = std::min(sdist, sd);
+  };
+  return sdist;
+}
+
+double signed_distance(const Difference &group, double xx, double yy,
+                       double zz) {
+  auto sdist = signed_distance(*group.first_obj, xx, yy, zz);
+
+  for (const auto &member : group.next_objs.objs) {
+    auto sd = signed_distance(member, xx, yy, zz);
+    sdist = std::min(sdist, -sd);
+  };
+  return sdist;
+}
+
+double signed_distance(const Mulmatrix &mm, double xx, double yy, double zz) {
+  // TODO: Invert non-orthogonal matrices
+  auto XX = xx - mm.translation[0];
+  auto YY = yy - mm.translation[1];
+  auto ZZ = zz - mm.translation[2];
+  return signed_distance(
+      mm.group,
+      mm.rotation[0][0] * XX + mm.rotation[1][0] * YY + mm.rotation[2][0] * ZZ,
+      mm.rotation[0][1] * XX + mm.rotation[1][1] * YY + mm.rotation[2][1] * ZZ,
+      mm.rotation[0][2] * XX + mm.rotation[1][2] * YY + mm.rotation[2][2] * ZZ);
+}
+
+double signed_distance(const Cone &cone, double xx, double yy, double zz) {
+  double ZZ = cone.center ? zz + cone.height / 2 : zz;
+  double sign_z = (ZZ >= 0 && ZZ <= cone.height) ? -1.0 : 1.0;
+  auto dist_z = std::min(std::fabs(ZZ), std::fabs(ZZ - cone.height));
+
+  auto rr = cone.radius1 + (ZZ * (cone.radius2 - cone.radius1) / cone.height);
+  auto delta_r = xx * xx + yy * yy - rr * rr;
+  double sign_r = delta_r <= 0 ? -1.0 : 1.0;
+  auto dist_r = std::fabs(delta_r);
+
+  return EXTERNAL_FLOW * std::max(sign_z * dist_z, sign_r * dist_r);
+}
+
+double signed_distance(const Cube &cube, double xx, double yy, double zz) {
+
+  auto [Lx, Ly, Lz] = cube.size;
+  auto XX = cube.center ? xx + Lx / 2 : xx;
+  double sign_x = (XX >= 0 && XX <= Lx) ? -1.0 : 1.0;
+  auto dist_x = std::min(std::fabs(XX), std::fabs(XX - Lx));
+
+  auto YY = cube.center ? yy + Ly / 2 : yy;
+  double sign_y = (YY >= 0 && YY <= Ly) ? -1.0 : 1.0;
+  auto dist_y = std::min(std::fabs(YY), std::fabs(YY - Ly));
+
+  auto ZZ = cube.center ? zz + Lz / 2 : zz;
+  double sign_z = (ZZ >= 0 && ZZ <= Lz) ? -1.0 : 1.0;
+  auto dist_z = std::min(std::fabs(ZZ), std::fabs(ZZ - Lz));
+
+  return EXTERNAL_FLOW *
+         std::max({sign_x * dist_x, sign_y * dist_y, sign_z * dist_z});
+}
+
+double signed_distance(const Sphere &sph, double xx, double yy, double zz) {
+  auto delta_r = xx * xx + yy * yy + zz * zz - sph.radius * sph.radius;
+  double sign = delta_r <= 0 ? -1.0 : 1.0;
+  auto dist = std::fabs(delta_r);
+  return EXTERNAL_FLOW * sign * dist;
+}
+
+double signed_distance(const Cylinder &cyl, double xx, double yy, double zz) {
+  auto ZZ = cyl.center ? zz + cyl.height / 2 : zz;
+  double sign_z = (ZZ >= 0 && ZZ <= cyl.height) ? -1.0 : 1.0;
+  auto dist_z = std::min(std::fabs(ZZ), std::fabs(ZZ - cyl.height));
+
+  auto rr = cyl.radius;
+  auto delta_r = xx * xx + yy * yy - rr * rr;
+  double sign_r = delta_r <= 0 ? -1.0 : 1.0;
+  auto dist_r = std::fabs(delta_r);
+
+  return EXTERNAL_FLOW * std::max(sign_z * dist_z, sign_r * dist_r);
+}
+
+double signed_distance(const Type &obj, double xx, double yy, double zz) {
+
+  return std::visit(
+      overloaded{
+          [xx, yy, zz](auto &&arg) { return signed_distance(arg, xx, yy, zz); },
+      },
+      obj);
+}
+
+} // namespace csg
diff --git a/src/csg/impl/parser.cpp b/src/csg/impl/parser.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a6cd4b1bfb4067f878c29f434a871523bf279716
--- /dev/null
+++ b/src/csg/impl/parser.cpp
@@ -0,0 +1,403 @@
+// standard includes
+#include <assert.h>
+#include <iostream>
+#include <memory>
+
+// subproject includes
+#include <tao/pegtl.hpp>
+
+// includes
+#include "csg_types.hpp"
+
+using namespace tao::pegtl;
+
+namespace {
+
+struct parser_state {
+  std::vector<std::vector<csg::Type>> current_objs;
+  std::vector<csg::Type> current_group;
+  std::string current_name;
+  std::vector<double> current_vec;
+  std::vector<std::vector<double>> current_matrix;
+  std::vector<std::vector<std::vector<double>>> current_matrices;
+  std::map<std::string,
+           std::variant<double, bool, std::vector<double>, std::string>>
+      curr_attr;
+};
+
+} // namespace
+
+namespace {
+
+struct escaped_x : seq<one<'x'>, rep<2, must<xdigit>>> {};
+struct escaped_u : seq<one<'u'>, rep<4, must<xdigit>>> {};
+struct escaped_U : seq<one<'U'>, rep<8, must<xdigit>>> {};
+struct escaped_c
+    : one<'\'', '"', '?', '\\', 'a', 'b', 'f', 'n', 'r', 't', 'v'> {};
+
+struct escaped : sor<escaped_x, escaped_u, escaped_U, escaped_c> {};
+
+struct character
+    : if_must_else<one<'\\'>, escaped, utf8::range<0x20, 0x10FFFF>> {};
+struct string_literal : if_must<one<'"'>, until<one<'"'>, character>> {};
+
+struct plus_minus : opt<one<'+', '-'>> {};
+struct dot : one<'.'> {};
+
+struct decimal
+    : if_then_else<dot, plus<digit>, seq<plus<digit>, opt<dot, star<digit>>>> {
+};
+struct e : one<'e', 'E'> {};
+struct p : one<'p', 'P'> {};
+struct exponent : seq<plus_minus, plus<digit>> {};
+struct double_ : seq<plus_minus, decimal, opt<e, exponent>> {};
+
+struct padded_double : pad<double_, space> {};
+
+struct L_FUN : pad<string<'('>, space> {};
+struct R_FUN : pad<string<')'>, space> {};
+struct L_ARR : pad<string<'['>, space> {};
+struct R_ARR : pad<string<']'>, space> {};
+struct L_BLK : pad<string<'{'>, space> {};
+struct R_BLK : pad<string<'}'>, space> {};
+struct S_CLN : pad<string<';'>, space> {};
+
+struct true_literal : string<'t', 'r', 'u', 'e'> {};
+struct false_literal : string<'f', 'a', 'l', 's', 'e'> {};
+
+struct boolean_literal : sor<true_literal, false_literal> {};
+
+struct special_identifier : seq<string<'$'>, identifier> {};
+
+struct name : sor<identifier, special_identifier> {};
+
+struct vector_cell : padded_double {};
+
+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 keyval : seq<pad<name, space>, string<'='>, pad<value, space>> {};
+
+struct attr_list : list<keyval, one<','>> {};
+
+struct row : vector {};
+
+struct matrix : seq<L_ARR, list<row, one<','>>, R_ARR> {};
+
+struct sphere
+    : seq<string<'s', 'p', 'h', 'e', 'r', 'e'>, L_FUN, attr_list, R_FUN> {};
+
+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 shape : seq<sor<cube, cylinder, sphere>, opt<S_CLN>> {};
+
+struct obj_list;
+
+struct bool_union : seq<string<'u', 'n', 'i', 'o', 'n'>, L_FUN, R_FUN, L_BLK,
+                        obj_list, R_BLK> {};
+
+struct bool_intersection
+    : seq<string<'i', 'n', 't', 'e', 'r', 's', 'e', 'c', 't', 'i', 'o', 'n'>,
+          L_FUN, R_FUN, L_BLK, obj_list, R_BLK> {};
+
+struct bool_diff : seq<string<'d', 'i', 'f', 'f', 'e', 'r', 'e', 'n', 'c', 'e'>,
+                       L_FUN, R_FUN, L_BLK, obj_list, R_BLK> {};
+
+struct bool_exp : sor<bool_union, bool_intersection, bool_diff> {};
+
+struct mulmat : seq<string<'m', 'u', 'l', 't', 'm', 'a', 't', 'r', 'i', 'x'>,
+                    L_FUN, matrix, R_FUN, L_BLK, obj_list, R_BLK> {};
+
+struct group : seq<string<'g', 'r', 'o', 'u', 'p'>, L_FUN, R_FUN,
+                   opt<seq<L_BLK, obj_list, R_BLK>>, opt<S_CLN>> {};
+
+struct render : seq<string<'r', 'e', 'n', 'd', 'e', 'r'>, L_FUN, attr_list,
+                    R_FUN, opt<seq<L_BLK, obj_list, R_BLK>>, opt<S_CLN>> {};
+
+struct colorgroup : seq<string<'c', 'o', 'l', 'o', 'r'>, L_FUN, vector, R_FUN,
+                        opt<seq<L_BLK, obj_list, R_BLK>>, opt<S_CLN>> {};
+
+struct csg_obj : sor<shape, mulmat, bool_exp, group, render, colorgroup> {};
+
+struct obj_list : plus<pad<csg_obj, space>> {};
+
+struct grammar : seq<obj_list, eof> {};
+
+template <typename Rule> struct action {};
+
+template <> struct action<L_BLK> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    // std::cout << " { " << std::endl;
+    std::vector<csg::Type> new_group;
+    st.current_objs.push_back(new_group);
+  }
+};
+
+template <> struct action<R_BLK> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    // std::cout << " } " << std::endl;
+    st.current_group.clear();
+    for (auto obj : st.current_objs.back()) {
+      st.current_group.push_back(obj);
+    }
+    st.current_objs.pop_back();
+  }
+};
+
+template <> struct action<vector_cell> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    std::stringstream ss(in.string());
+    double cell;
+    ss >> cell;
+    st.current_vec.push_back(cell);
+  }
+};
+
+template <> struct action<matrix> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    st.current_matrices.push_back(st.current_matrix);
+    st.current_matrix.clear();
+  }
+};
+
+template <> struct action<row> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    st.current_matrix.push_back(st.current_vec);
+    st.current_vec.clear();
+  }
+};
+
+void add_group(parser_state &st) {
+  csg::Union group;
+  for (const auto &curr_obj : st.current_group) {
+    group.objs.push_back(curr_obj);
+  }
+  st.current_objs.back().push_back(group);
+}
+
+template <> struct action<group> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    add_group(st);
+  }
+};
+
+template <> struct action<render> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    add_group(st);
+  }
+};
+
+template <> struct action<colorgroup> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    add_group(st);
+  }
+};
+
+template <> struct action<bool_union> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    add_group(st);
+  }
+};
+
+template <> struct action<bool_intersection> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    auto csg_in = csg::Intersection();
+    for (const auto &curr_obj : st.current_group) {
+      csg_in.objs.push_back(curr_obj);
+    }
+    st.current_objs.back().push_back(csg_in);
+  }
+};
+
+template <> struct action<bool_diff> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    auto csg_diff = csg::Difference();
+    for (auto it = st.current_group.begin(); it != st.current_group.end();
+         ++it) {
+      if (it == st.current_group.begin()) {
+        csg_diff.first_obj = std::make_shared<csg::Type>(*it);
+      } else {
+        csg_diff.next_objs.objs.push_back(*it);
+      }
+    }
+    st.current_objs.back().push_back(csg_diff);
+  }
+};
+
+template <> struct action<mulmat> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    auto mulmat = csg::Mulmatrix();
+    auto mat = st.current_matrices.back();
+    // clang-format off
+    mulmat.rotation[0] = {mat[0][0], mat[0][1], mat[0][2]};
+    mulmat.rotation[1] = {mat[1][0], mat[1][1], mat[1][2]};
+    mulmat.rotation[2] = {mat[2][0], mat[2][1], mat[2][2]};
+    // clang-format on
+
+    mulmat.translation = {
+        mat[0][3],
+        mat[1][3],
+        mat[2][3],
+    };
+
+    for (const auto &curr_obj : st.current_group) {
+      mulmat.group.objs.push_back(curr_obj);
+    }
+    st.current_objs.back().push_back(mulmat);
+    st.current_matrices.pop_back();
+  }
+};
+
+template <> struct action<cube> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    csg::Cube cub;
+    auto size = std::get<std::vector<double>>(st.curr_attr["size"]);
+    cub.size = {size[0], size[1], size[2]};
+    cub.center = std::get<bool>(st.curr_attr["center"]);
+
+    st.current_objs.back().push_back(cub);
+    st.curr_attr.clear();
+  }
+};
+
+template <> struct action<cylinder> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    if (st.curr_attr.count("r")) {
+      // proper cylinder
+      if (st.curr_attr.count("r1") || st.curr_attr.count("r2")) {
+        std::cout << " ERROR: cannot specify both r and (r1 or r2); ambiguous "
+                  << std::endl;
+      }
+      csg::Cylinder cyl;
+      cyl.center = std::get<bool>(st.curr_attr["center"]);
+      cyl.height = std::get<double>(st.curr_attr["h"]);
+      cyl.radius = std::get<double>(st.curr_attr["r"]);
+      st.current_objs.back().push_back(cyl);
+
+    } else {
+      // conic "cylinder"
+      csg::Cone cone;
+      cone.center = std::get<bool>(st.curr_attr["center"]);
+      cone.height = std::get<double>(st.curr_attr["h"]);
+      cone.radius1 = std::get<double>(st.curr_attr["r1"]);
+      cone.radius2 = std::get<double>(st.curr_attr["r2"]);
+      st.current_objs.back().push_back(cone);
+    }
+
+    st.curr_attr.clear();
+  }
+};
+
+template <> struct action<sphere> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    csg::Sphere sph;
+    sph.radius = std::get<double>(st.curr_attr["r"]);
+
+    st.current_objs.back().push_back(sph);
+    st.curr_attr.clear();
+  }
+};
+
+template <> struct action<name> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    std::stringstream ss(in.string());
+    ss >> st.current_name;
+  }
+};
+
+template <> struct action<string_literal> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    std::stringstream ss(in.string());
+    std::string v;
+    ss >> v;
+    st.curr_attr[st.current_name] = v;
+  }
+};
+
+template <> struct action<double_> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    std::stringstream ss(in.string());
+    double v;
+    ss >> v;
+    st.curr_attr[st.current_name] = v;
+  }
+};
+
+template <> struct action<true_literal> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    st.curr_attr[st.current_name] = true;
+  }
+};
+
+template <> struct action<false_literal> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    st.curr_attr[st.current_name] = false;
+  }
+};
+
+template <> struct action<vector_attr> {
+  template <typename Input>
+  static void apply(const Input &in, parser_state &st) {
+    st.curr_attr[st.current_name] = st.current_vec;
+    st.current_vec.clear();
+  }
+};
+
+std::optional<parser_state> do_parse(std::string str) {
+  parser_state st;
+  std::vector<csg::Type> new_group;
+  st.current_objs.push_back(new_group);
+  memory_input in(str, "std::cin");
+  if (!parse<grammar, action>(in, st)) {
+    return std::nullopt;
+  }
+  return st;
+}
+
+} // namespace
+
+namespace csg {
+
+std::shared_ptr<Tree> parse_csg(std::string str) {
+  auto maybe_state = do_parse(str);
+  if (!maybe_state.has_value()) {
+    return nullptr;
+  }
+  auto st = maybe_state.value();
+
+  assert(st.current_objs.size() == 1);
+  Tree tree;
+  for (auto obj : st.current_objs.back()) {
+    tree.top.objs.push_back(obj);
+  }
+  assert(tree.top.objs.size() > 0); // Disallow empty .csg file
+  return std::make_shared<Tree>(tree);
+}
+} // namespace csg
diff --git a/src/csg/meson.build b/src/csg/meson.build
new file mode 100644
index 0000000000000000000000000000000000000000..ca6bcd515dad2ba4616fd2c7cf0e32908a7ca2e8
--- /dev/null
+++ b/src/csg/meson.build
@@ -0,0 +1,9 @@
+lib_csg_parser = static_library(
+  'csg-parser',
+  'impl/csg.cpp',
+  'impl/levelset.cpp',
+  'impl/parser.cpp',
+  include_directories: tao_inc,
+  install : true)
+
+subdir('tests')
diff --git a/src/csg/tests/levelset/boolean.t.cpp b/src/csg/tests/levelset/boolean.t.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..461d254ccec3188da2dae66e0ca19ea72c632cd8
--- /dev/null
+++ b/src/csg/tests/levelset/boolean.t.cpp
@@ -0,0 +1,92 @@
+#include "catch2/catch.hpp"
+
+#include <csg.hpp>
+#include <csg_types.hpp>
+
+#include <memory>
+
+namespace {
+
+TEST_CASE("Union", "[Levelset Boolean]") {
+  // Create a union of a cube of size 12 and a sphere of radius 8
+  csg::Cube my_cub{.size = {12, 12, 12}, .center = true};
+  csg::Sphere my_sph{.radius = 8};
+
+  auto my_union = csg::Union();
+  my_union.objs.push_back(my_cub);
+  my_union.objs.push_back(my_sph);
+
+  auto my_tree = std::make_shared<csg::Tree>();
+  my_tree->top.objs.push_back(my_union);
+  csg::CsgIF my_levelset(my_tree);
+
+  CHECK(0 < my_levelset(0, 7, 0));
+  CHECK(Approx(8 * 8 - 7 * 7) == my_levelset(0, 7, 0));
+
+  CHECK(0 < my_levelset(0, 0, -7.5));
+  CHECK(Approx(8 * 8 - 7.5 * 7.5) == my_levelset(0, 0, -7.5));
+
+  CHECK_FALSE(0 < my_levelset(0, -9, 0));
+  CHECK_FALSE(Approx(-1) == my_levelset(0, -9, 0));
+
+  CHECK_FALSE(0 < my_levelset(9, 0, 0));
+  CHECK_FALSE(Approx(-1) == my_levelset(9, 0, 0));
+}
+
+TEST_CASE("Intersection", "[Levelset Boolean]") {
+  // Create an of a cube of size 12 and a sphere of radius 8
+  csg::Cube my_cub{.size = {12, 12, 12}, .center = true};
+  csg::Sphere my_sph{.radius = 8};
+
+  auto my_in = csg::Intersection();
+  my_in.objs.push_back(my_cub);
+  my_in.objs.push_back(my_sph);
+
+  auto my_tree = std::make_shared<csg::Tree>();
+  my_tree->top.objs.push_back(my_in);
+  csg::CsgIF my_levelset(my_tree);
+
+  CHECK_FALSE(0 < my_levelset(0, 7, 0));
+  CHECK(Approx(-1) == my_levelset(0, 7, 0));
+
+  CHECK_FALSE(0 < my_levelset(0, 0, -7.5));
+  CHECK(Approx(-1.5) == my_levelset(0, 0, -7.5));
+
+  CHECK(0 < my_levelset(0, 5, 0));
+  CHECK(Approx(1) == my_levelset(0, 5, 0));
+
+  CHECK(0 < my_levelset(-5, 0, 0));
+  CHECK(Approx(1) == my_levelset(-5, 0, 0));
+}
+
+TEST_CASE("Difference", "[Levelset Boolean]") {
+  // Create an of a cube of size 12 and remove a sphere of radius 8
+  csg::Cube my_cub{.size = {12, 12, 12}, .center = true};
+  csg::Sphere my_sph{.radius = 8};
+
+  csg::Union my_union;
+  my_union.objs.push_back(my_sph);
+
+  auto my_diff = csg::Difference({
+      std::make_shared<csg::Type>(my_cub),
+      csg::Union(my_union),
+  });
+
+  auto my_tree = std::make_shared<csg::Tree>();
+  my_tree->top.objs.push_back(my_diff);
+  csg::CsgIF my_levelset(my_tree);
+
+  CHECK_FALSE(0 < my_levelset(0, 0, 0));
+  CHECK(Approx(-(8 * 8)) == my_levelset(0, 0, 0));
+
+  CHECK_FALSE(0 < my_levelset(-4, 0, 0));
+  CHECK(Approx(-(8 * 8 - 4 * 4)) == my_levelset(-4, 0, 0));
+
+  CHECK(0 < my_levelset(5.99, 5.99, 0));
+  CHECK(Approx(0.01) == my_levelset(5.99, 5.99, 0));
+
+  CHECK(0 < my_levelset(-5.99, 0, 5.99));
+  CHECK(Approx(0.01) == my_levelset(-5.99, 0, 5.99));
+}
+
+} // namespace
diff --git a/src/csg/tests/levelset/primitives.t.cpp b/src/csg/tests/levelset/primitives.t.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5a02af4771405728f5f05ef8c1168a8b2291a547
--- /dev/null
+++ b/src/csg/tests/levelset/primitives.t.cpp
@@ -0,0 +1,539 @@
+#include "catch2/catch.hpp"
+
+#include <csg.hpp>
+#include <csg_types.hpp>
+
+namespace {
+
+TEST_CASE("Cylinder", "[Levelset Primitives]") {
+
+  csg::Cylinder my_cyl{.radius = 1, .height = 2, .center = false};
+
+  SECTION("Not Centered") {
+    my_cyl.center = false;
+
+    auto my_tree = std::make_shared<csg::Tree>();
+    my_tree->top.objs.push_back(my_cyl);
+    csg::CsgIF my_levelset(my_tree);
+
+    SECTION("Axis") {
+      CHECK_FALSE(0 < my_levelset(0, 0, -1));
+      CHECK(Approx(-1) == my_levelset(0, 0, -1));
+
+      CHECK(0 < my_levelset(0, 0, 1));
+      CHECK(Approx(1) == my_levelset(0, 0, 1));
+
+      CHECK_FALSE(0 < my_levelset(0, 0, 3));
+      CHECK(Approx(-1) == my_levelset(0, 0, 3));
+    }
+
+    SECTION("Cross section") {
+      CHECK_FALSE(0 < my_levelset(-0.8, 0.8, 1));
+      CHECK(Approx(-0.28) == my_levelset(-0.8, 0.8, 1));
+
+      CHECK_FALSE(0 < my_levelset(0.8, -0.8, 1));
+      CHECK(Approx(-0.28) == my_levelset(0.8, -0.8, 1));
+
+      CHECK(0 < my_levelset(0, 0.8, 1));
+      CHECK(Approx(0.36) == my_levelset(0, 0.8, 1));
+
+      CHECK(0 < my_levelset(0.8, 0, 1));
+      CHECK(Approx(0.36) == my_levelset(0.8, 0, 1));
+    }
+  }
+  SECTION("Centered") {
+    my_cyl.center = true;
+
+    auto my_tree = std::make_shared<csg::Tree>();
+    my_tree->top.objs.push_back(my_cyl);
+    csg::CsgIF my_levelset(my_tree);
+
+    SECTION("Axis") {
+      CHECK_FALSE(0 < my_levelset(0, 0, -2));
+      CHECK(Approx(-1) == my_levelset(0, 0, -2));
+
+      CHECK(0 < my_levelset(0, 0, 0));
+      CHECK(Approx(1) == my_levelset(0, 0, 0));
+
+      CHECK_FALSE(0 < my_levelset(0, 0, 2));
+      CHECK(Approx(-1) == my_levelset(0, 0, 2));
+    }
+
+    SECTION("Cross section") {
+      CHECK_FALSE(0 < my_levelset(-0.8, 0.8, 0));
+      CHECK(Approx(-0.28) == my_levelset(-0.8, 0.8, 0));
+
+      CHECK_FALSE(0 < my_levelset(0.8, -0.8, 0));
+      CHECK(Approx(-0.28) == my_levelset(0.8, -0.8, 0));
+
+      CHECK(0 < my_levelset(0, 0.8, 0));
+      CHECK(Approx(0.36) == my_levelset(0, 0.8, 0));
+
+      CHECK(0 < my_levelset(0.8, 0, 0));
+      CHECK(Approx(0.36) == my_levelset(0.8, 0, 0));
+    }
+  }
+}
+
+TEST_CASE("Cube", "[Levelset Primitives]") {
+
+  csg::Cube my_cub{.size = {2, 3, 5}, .center = false};
+
+  SECTION("Not Centered") {
+    my_cub.center = false;
+
+    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, 1, -1));
+      CHECK(Approx(-1) == my_levelset(1, 1, -1));
+
+      CHECK_FALSE(0 < my_levelset(1, -1, 1));
+      CHECK(Approx(-1) == my_levelset(1, -1, 1));
+
+      CHECK_FALSE(0 < my_levelset(-1, 1, 1));
+      CHECK(Approx(-1) == my_levelset(-1, 1, 1));
+
+      CHECK_FALSE(0 < my_levelset(3, 1, 1));
+      CHECK(Approx(-1) == my_levelset(3, 1, 1));
+
+      CHECK_FALSE(0 < my_levelset(1, 4, 1));
+      CHECK(Approx(-1) == my_levelset(1, 4, 1));
+
+      CHECK_FALSE(0 < my_levelset(1, 1, 6));
+      CHECK(Approx(-1) == my_levelset(1, 1, 6));
+    }
+    SECTION("Inside") {
+      CHECK(0 < my_levelset(1, 1, 4));
+      CHECK(Approx(1) == my_levelset(1, 1, 4));
+
+      CHECK(0 < my_levelset(1, 2, 1));
+      CHECK(Approx(1) == my_levelset(1, 2, 1));
+
+      CHECK(0 < my_levelset(1, 1, 1));
+      CHECK(Approx(1) == my_levelset(1, 1, 1));
+    }
+  }
+
+  SECTION("Centered") {
+    my_cub.center = true;
+
+    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(-0.5, 1, -3));
+      CHECK(Approx(-0.5) == my_levelset(-0.5, 1, -3));
+
+      CHECK_FALSE(0 < my_levelset(-0.5, 2, -2));
+      CHECK(Approx(-0.5) == my_levelset(-0.5, 2, -2));
+
+      CHECK_FALSE(0 < my_levelset(-1.5, 1, -2));
+      CHECK(Approx(-0.5) == my_levelset(-1.5, 1, -2));
+
+      CHECK_FALSE(0 < my_levelset(0.5, -1, 3));
+      CHECK(Approx(-0.5) == my_levelset(0.5, -1, 3));
+
+      CHECK_FALSE(0 < my_levelset(0.5, -2, 2));
+      CHECK(Approx(-0.5) == my_levelset(0.5, -2, 2));
+
+      CHECK_FALSE(0 < my_levelset(1.5, -1, 2));
+      CHECK(Approx(-0.5) == my_levelset(1.5, -1, 2));
+    }
+    SECTION("Inside") {
+      CHECK(0 < my_levelset(0.5, -1, 2));
+      CHECK(Approx(0.5) == my_levelset(0.5, -1, 2));
+
+      CHECK(0 < my_levelset(-0.5, 1, -2));
+      CHECK(Approx(0.5) == my_levelset(-0.5, 1, -2));
+    }
+  }
+}
+
+TEST_CASE("Sphere", "[Levelset Primitives]") {
+
+  csg::Sphere my_sph{.radius = 10};
+
+  auto my_tree = std::make_shared<csg::Tree>();
+  my_tree->top.objs.push_back(my_sph);
+  csg::CsgIF my_levelset(my_tree);
+
+  SECTION("Outside") {
+    CHECK_FALSE(0 < my_levelset(0, 8, 8));
+    CHECK(Approx(-28) == my_levelset(0, 8, 8));
+
+    CHECK_FALSE(0 < my_levelset(-8, 0, 8));
+    CHECK(Approx(-28) == my_levelset(-8, 0, 8));
+
+    CHECK_FALSE(0 < my_levelset(0, 8, -8));
+    CHECK(Approx(-28) == my_levelset(0, 8, -8));
+  }
+  SECTION("Inside") {
+    CHECK(0 < my_levelset(0, 0, 0));
+    CHECK(Approx(100) == my_levelset(0, 0, 0));
+
+    CHECK(0 < my_levelset(0, 7, 7));
+    CHECK(Approx(2) == my_levelset(0, 7, 7));
+
+    CHECK(0 < my_levelset(-7, 0, 7));
+    CHECK(Approx(2) == my_levelset(-7, 0, 7));
+
+    CHECK(0 < my_levelset(0, 7, -7));
+    CHECK(Approx(2) == my_levelset(0, 7, -7));
+  }
+}
+TEST_CASE("Cone", "[Levelset Primitives]") {
+
+  SECTION("radius1 < radius 2; regular)") {
+
+    csg::Cone my_cone{.radius1 = 0, .radius2 = 1, .height = 2, .center = false};
+
+    SECTION("Untruncated") {
+      my_cone.radius1 = 0;
+
+      SECTION("Not Centered") {
+        my_cone.center = false;
+
+        auto my_tree = std::make_shared<csg::Tree>();
+        my_tree->top.objs.push_back(my_cone);
+        csg::CsgIF my_levelset(my_tree);
+
+        SECTION("Outside") {
+          CHECK_FALSE(0 < my_levelset(0.15, -0.15, 0.1));
+          CHECK(
+              Approx((0.5 * 0.1) * (0.5 * 0.1) - (0.15 * 0.15 + 0.15 * .15)) ==
+              my_levelset(0.15, -0.15, 0.1));
+
+          CHECK_FALSE(0 < my_levelset(-0.15, 0.15, 0.1));
+          CHECK(
+              Approx((0.5 * 0.1) * (0.5 * 0.1) - (0.15 * 0.15 + 0.15 * .15)) ==
+              my_levelset(-0.15, 0.15, 0.1));
+
+          CHECK_FALSE(0 < my_levelset(0, 0.99, 0.9));
+          CHECK(Approx((0.5 * 0.9) * (0.5 * 0.9) - 0.99 * 0.99) ==
+                my_levelset(0, 0.99, 0.9));
+
+          CHECK_FALSE(0 < my_levelset(-0.99, 0, 0.9));
+          CHECK(Approx((0.5 * 0.9) * (0.5 * 0.9) - 0.99 * 0.99) ==
+                my_levelset(-0.99, 0, 0.9));
+
+          CHECK_FALSE(0 < my_levelset(0, 0, -0.1));
+          CHECK(Approx(-0.1) == my_levelset(0, 0, -0.1));
+
+          CHECK_FALSE(0 < my_levelset(0, 0, 2.1));
+          CHECK(Approx(-0.1) == my_levelset(0, 0, 2.1));
+        }
+
+        SECTION("Inside") {
+          CHECK(0 < my_levelset(0, 0, 0.1));
+          CHECK(Approx((0.5 * 0.1) * (0.5 * 0.1)) == my_levelset(0, 0, 0.1));
+
+          CHECK(0 < my_levelset(0, 0, 1.9));
+          CHECK(Approx(0.1) == my_levelset(0, 0, 1.9));
+
+          CHECK(0 < my_levelset(0.01, 0, 0.1));
+          CHECK(Approx((0.5 * 0.1) * (0.5 * 0.1) - 0.01 * 0.01) ==
+                my_levelset(0.01, 0, 0.1));
+
+          CHECK(0 < my_levelset(0, -0.01, 0.1));
+          CHECK(Approx((0.5 * 0.1) * (0.5 * 0.1) - 0.01 * 0.01) ==
+                my_levelset(0, -0.01, 0.1));
+
+          CHECK(0 < my_levelset(0, 0.4, 0.99));
+          CHECK(Approx((0.5 * 0.99) * (0.5 * 0.99) - 0.4 * 0.4) ==
+                my_levelset(0, 0.4, 0.99));
+
+          CHECK(0 < my_levelset(-0.4, 0, 0.9));
+          CHECK(Approx((0.5 * 0.9) * (0.5 * 0.9) - 0.4 * 0.4) ==
+                my_levelset(-0.4, 0, 0.9));
+        }
+      }
+
+      SECTION("Centered") {
+        my_cone.center = true;
+
+        auto my_tree = std::make_shared<csg::Tree>();
+        my_tree->top.objs.push_back(my_cone);
+        csg::CsgIF my_levelset(my_tree);
+
+        SECTION("Outside") {
+          CHECK_FALSE(0 < my_levelset(0.15, -0.15, -0.9));
+          CHECK(
+              Approx((0.5 * 0.1) * (0.5 * 0.1) - (0.15 * 0.15 + 0.15 * 0.15)) ==
+              my_levelset(0.15, -0.15, -0.9));
+
+          CHECK_FALSE(0 < my_levelset(-0.15, 0.15, -0.9));
+          CHECK(
+              Approx((0.5 * 0.1) * (0.5 * 0.1) - (0.15 * 0.15 + 0.15 * 0.15)) ==
+              my_levelset(-0.15, 0.15, -0.9));
+
+          CHECK_FALSE(0 < my_levelset(0, 0.99, -0.1));
+          CHECK(Approx((0.5 * 0.9) * (0.5 * 0.9) - 0.99 * 0.99) ==
+                my_levelset(0, 0.99, -0.1));
+
+          CHECK_FALSE(0 < my_levelset(-0.99, 0, -0.1));
+          CHECK(Approx((0.5 * 0.9) * (0.5 * 0.9) - 0.99 * 0.99) ==
+                my_levelset(-0.99, 0, -0.1));
+
+          CHECK_FALSE(0 < my_levelset(0, 0, -1.1));
+          CHECK(Approx(-0.1) == my_levelset(0, 0, -1.1));
+
+          CHECK_FALSE(0 < my_levelset(0, 0, 1.1));
+          CHECK(Approx(-0.1) == my_levelset(0, 0, 1.1));
+        }
+
+        SECTION("Inside") {
+          CHECK(0 < my_levelset(0, 0, -0.9));
+          CHECK(Approx((0.5 * 0.1) * (0.5 * 0.1)) == my_levelset(0, 0, -0.9));
+
+          CHECK(0 < my_levelset(0, 0, 0.9));
+          CHECK(Approx(0.1) == my_levelset(0, 0, 0.9));
+
+          CHECK(0 < my_levelset(0.01, 0, -0.9));
+          CHECK(Approx((0.5 * 0.1) * (0.5 * 0.1) - 0.01 * 0.01) ==
+                my_levelset(0.01, 0, -0.9));
+
+          CHECK(0 < my_levelset(0, -0.01, -0.9));
+          CHECK(Approx((0.5 * 0.1) * (0.5 * 0.1) - 0.01 * 0.01) ==
+                my_levelset(0, -0.01, -0.9));
+
+          CHECK(0 < my_levelset(0, 0.4, -0.01));
+          CHECK(Approx((0.5 * 0.99) * (0.5 * 0.99) - 0.4 * 0.4) ==
+                my_levelset(0, 0.4, -0.01));
+
+          CHECK(0 < my_levelset(-0.4, 0, -0.1));
+          CHECK(Approx((0.5 * 0.9) * (0.5 * 0.9) - 0.4 * 0.4) ==
+                my_levelset(-0.4, 0, -0.1));
+        }
+      }
+    }
+
+    SECTION("Truncated") {
+      my_cone.radius1 = 0.6;
+
+      SECTION("Not Centered") {
+        my_cone.center = false;
+
+        auto my_tree = std::make_shared<csg::Tree>();
+        my_tree->top.objs.push_back(my_cone);
+        csg::CsgIF my_levelset(my_tree);
+
+        SECTION("smaller end") {
+          CHECK_FALSE(0 < my_levelset(0, 0, -0.01));
+          CHECK(Approx(-0.01) == my_levelset(0, 0, -0.01));
+
+          CHECK(0 < my_levelset(0, 0, 0.01));
+          CHECK(Approx(0.01) == my_levelset(0, 0, 0.01));
+        }
+        SECTION("curved surface near smaller end") {
+          CHECK_FALSE(0 < my_levelset(0.61, 0, 0.01));
+          CHECK(Approx((0.6 + 0.2 * 0.01) * (0.6 + 0.2 * 0.01) - 0.61 * 0.61) ==
+                my_levelset(0.61, 0, 0.01));
+
+          CHECK(0 < my_levelset(0.59, 0, 0.01));
+          CHECK(Approx(0.01) == my_levelset(0.59, 0, 0.01));
+
+          CHECK_FALSE(0 < my_levelset(-0.5, 0.5, 0.4));
+          CHECK(Approx((0.6 + 0.2 * 0.4) * (0.6 + 0.2 * 0.4) -
+                       (0.5 * 0.5 + 0.5 * 0.5)) == my_levelset(-0.5, 0.5, 0.4));
+
+          CHECK(0 < my_levelset(-0.3, 0.5, 0.4));
+          CHECK(Approx((0.6 + 0.2 * 0.4) * (0.6 + 0.2 * 0.4) -
+                       (0.3 * 0.3 + 0.5 * 0.5)) == my_levelset(-0.3, 0.5, 0.4));
+        }
+        SECTION("curved surface near the middle") {
+          CHECK_FALSE(0 < my_levelset(0.73, -0.3, 0.8));
+          CHECK(Approx((0.6 + 0.2 * 0.8) * (0.6 + 0.2 * 0.8) -
+                       (0.73 * 0.73 + 0.3 * 0.3)) ==
+                my_levelset(0.73, -0.3, 0.8));
+
+          CHECK(0 < my_levelset(0.63, -0.3, 0.8));
+          CHECK(Approx((0.6 + 0.2 * 0.8) * (0.6 + 0.2 * 0.8) -
+                       (0.63 * 0.63 + 0.3 * 0.3)) ==
+                my_levelset(0.63, -0.3, 0.8));
+
+          CHECK_FALSE(0 < my_levelset(-0.9, 0.3, 1.4));
+          CHECK(Approx((0.6 + 0.2 * 1.4) * (0.6 + 0.2 * 1.4) -
+                       (0.9 * 0.9 + 0.3 * 0.3)) == my_levelset(-0.9, 0.3, 1.4));
+
+          CHECK(0 < my_levelset(-0.8, 0.3, 1.4));
+          CHECK(Approx((0.6 + 0.2 * 1.4) * (0.6 + 0.2 * 1.4) -
+                       (0.8 * 0.8 + 0.3 * 0.3)) == my_levelset(-0.8, 0.3, 1.4));
+        }
+        SECTION("larger end") {
+          CHECK_FALSE(0 < my_levelset(0, 0, 2.01));
+          CHECK(Approx(-0.01) == my_levelset(0, 0, 2.01));
+
+          CHECK(0 < my_levelset(0, 0, 1.99));
+          CHECK(Approx(0.01) == my_levelset(0, 0, 1.99));
+        }
+        SECTION("curved surface near larger end") {
+          CHECK_FALSE(0 < my_levelset(1.1, 0, 1.9));
+          CHECK(Approx((0.6 + 0.2 * 1.9) * (0.6 + 0.2 * 1.9) - 1.1 * 1.1) ==
+                my_levelset(1.1, 0, 1.9));
+
+          CHECK(0 < my_levelset(0.9, 0, 1.9));
+          CHECK(Approx(0.1) == my_levelset(0.9, 0, 1.9));
+
+          CHECK_FALSE(0 < my_levelset(1, -0.3, 1.9));
+          CHECK(Approx((0.6 + 0.2 * 1.9) * (0.6 + 0.2 * 1.9) -
+                       (0.3 * 0.3 + 1 * 1)) == my_levelset(1, -0.3, 1.9));
+
+          CHECK(0 < my_levelset(0.9, -0.3, 1.9));
+          CHECK(Approx((0.6 + 0.2 * 1.9) * (0.6 + 0.2 * 1.9) -
+                       (0.3 * 0.3 + 0.9 * 0.9)) == my_levelset(0.9, -0.3, 1.9));
+        }
+      }
+
+      SECTION("Centered") {
+        my_cone.center = true;
+
+        auto my_tree = std::make_shared<csg::Tree>();
+        my_tree->top.objs.push_back(my_cone);
+        csg::CsgIF my_levelset(my_tree);
+
+        SECTION("smaller end") {
+          CHECK_FALSE(0 < my_levelset(0, 0, -1.01));
+          CHECK(Approx(-0.01) == my_levelset(0, 0, -1.01));
+
+          CHECK(0 < my_levelset(0, 0, -0.99));
+          CHECK(Approx(0.01) == my_levelset(0, 0, -0.99));
+
+          CHECK(0 < my_levelset(0, 0, -0.99));
+          CHECK(Approx(0.01) == my_levelset(0, 0, -0.99));
+        }
+        SECTION("curved surface near smaller end") {
+          CHECK_FALSE(0 < my_levelset(0.61, 0, -0.99));
+          CHECK(Approx((0.6 + 0.2 * 0.01) * (0.6 + 0.2 * 0.01) - 0.61 * 0.61) ==
+                my_levelset(0.61, 0, -0.99));
+
+          CHECK(0 < my_levelset(0.5, 0, -0.99));
+          CHECK(Approx(0.01) == my_levelset(0.5, 0, -0.99));
+
+          CHECK_FALSE(0 < my_levelset(0.65, -0.1, -0.9));
+          CHECK(Approx((0.6 + 0.2 * 0.1) * (0.6 + 0.2 * 0.1) -
+                       (0.65 * 0.65 + 0.1 * 0.1)) ==
+                my_levelset(0.65, -0.1, -0.9));
+
+          CHECK(0 < my_levelset(0.57, -0.1, -0.9));
+          CHECK(Approx((0.6 + 0.2 * 0.1) * (0.6 + 0.2 * 0.1) -
+                       (0.57 * 0.57 + 0.1 * 0.1)) ==
+                my_levelset(0.57, -0.1, -0.9));
+        }
+        SECTION("curved surface near middle") {
+          CHECK_FALSE(0 < my_levelset(-0.9, 0.3, 0.2));
+          CHECK(Approx((0.6 + 0.2 * 1.2) * (0.6 + 0.2 * 1.2) -
+                       (0.9 * 0.9 + 0.3 * 0.3)) == my_levelset(-0.9, 0.3, 0.2));
+
+          CHECK(0 < my_levelset(-0.7, 0.3, 0.2));
+          CHECK(Approx((0.6 + 0.2 * 1.2) * (0.6 + 0.2 * 1.2) -
+                       (0.7 * 0.7 + 0.3 * 0.3)) == my_levelset(-0.7, 0.3, 0.2));
+
+          CHECK_FALSE(0 < my_levelset(0.75, -0.3, -0.2));
+          CHECK(Approx((0.6 + 0.2 * 0.8) * (0.6 + 0.2 * 0.8) -
+                       (0.75 * 0.75 + 0.3 * 0.3)) ==
+                my_levelset(0.75, -0.3, -0.2));
+
+          CHECK(0 < my_levelset(0.6, -0.3, -0.2));
+          CHECK(Approx((0.6 + 0.2 * 0.8) * (0.6 + 0.2 * 0.8) -
+                       (0.6 * 0.6 + 0.3 * 0.3)) ==
+                my_levelset(0.6, -0.3, -0.2));
+        }
+        SECTION("larger end") {
+          CHECK_FALSE(0 < my_levelset(0, 0, 1.01));
+          CHECK(Approx(-0.01) == my_levelset(0, 0, 1.01));
+
+          CHECK(0 < my_levelset(0, 0, 0.99));
+          CHECK(Approx(0.01) == my_levelset(0, 0, 0.99));
+        }
+        SECTION("curved surface near the larger end") {
+          CHECK_FALSE(0 < my_levelset(1.1, 0, 0.95));
+          CHECK(Approx((0.6 + 0.2 * 1.95) * (0.6 + 0.2 * 1.95) - 1.1 * 1.1) ==
+                my_levelset(1.1, 0, 0.95));
+
+          CHECK(0 < my_levelset(0.8, 0, 0.95));
+          CHECK(Approx(0.05) == my_levelset(0.8, 0, 0.95));
+
+          CHECK_FALSE(0 < my_levelset(1.0, -0.3, 0.8));
+          CHECK(Approx((0.6 + 0.2 * 1.8) * (0.6 + 0.2 * 1.8) -
+                       (0.3 * 0.3 + 1 * 1)) == my_levelset(1.0, -0.3, 0.8));
+
+          CHECK(0 < my_levelset(0.73, -0.3, 0.8));
+          CHECK(Approx(0.2) == my_levelset(0.73, -0.3, 0.8));
+        }
+      }
+    }
+  }
+
+  SECTION("radius1 > radius2; inverted") {
+    csg::Cone my_cone{.radius1 = 1, .radius2 = 0, .height = 2, .center = false};
+
+    SECTION("Untruncated") {
+      my_cone.radius2 = 0;
+
+      SECTION("Not Centered") {
+        my_cone.center = false;
+
+        auto my_tree = std::make_shared<csg::Tree>();
+        my_tree->top.objs.push_back(my_cone);
+        csg::CsgIF my_levelset(my_tree);
+
+        SECTION("Outside") {
+          CHECK_FALSE(0 < my_levelset(0.15, -0.15, 1.9));
+          CHECK(Approx((1 - 0.5 * 1.9) * (1 - 0.5 * 1.9) -
+                       (0.15 * 0.15 + 0.15 * .15)) ==
+                my_levelset(0.15, -0.15, 1.9));
+
+          CHECK_FALSE(0 < my_levelset(-0.15, 0.15, 1.9));
+          CHECK(Approx((1 - 0.5 * 1.9) * (1 - 0.5 * 1.9) -
+                       (0.15 * 0.15 + 0.15 * .15)) ==
+                my_levelset(-0.15, 0.15, 1.9));
+
+          CHECK_FALSE(0 < my_levelset(0, 0.99, 1.1));
+          CHECK(Approx((1 - 0.5 * 1.1) * (1 - 0.5 * 1.1) - 0.99 * 0.99) ==
+                my_levelset(0, 0.99, 1.1));
+
+          CHECK_FALSE(0 < my_levelset(-0.99, 0, 1.1));
+          CHECK(Approx((1 - 0.5 * 1.1) * (1 - 0.5 * 1.1) - 0.99 * 0.99) ==
+                my_levelset(-0.99, 0, 1.1));
+
+          CHECK_FALSE(0 < my_levelset(0, 0, -0.1));
+          CHECK(Approx(-0.1) == my_levelset(0, 0, -0.1));
+
+          CHECK_FALSE(0 < my_levelset(0, 0, 2.1));
+          CHECK(Approx(-0.1) == my_levelset(0, 0, 2.1));
+        }
+
+        SECTION("Inside") {
+          CHECK(0 < my_levelset(0, 0, 1.9));
+          CHECK(Approx((1 - 0.5 * 1.9) * (1 - 0.5 * 1.9)) ==
+                my_levelset(0, 0, 1.9));
+
+          CHECK(0 < my_levelset(0, 0, 0.1));
+          CHECK(Approx(0.1) == my_levelset(0, 0, 0.1));
+
+          CHECK(0 < my_levelset(0.01, 0, 1.9));
+          CHECK(Approx((1 - 0.5 * 1.9) * (1 - 0.5 * 1.9) - 0.01 * 0.01) ==
+                my_levelset(0.01, 0, 1.9));
+
+          CHECK(0 < my_levelset(0, -0.01, 1.9));
+          CHECK(Approx((1 - 0.5 * 1.9) * (1 - 0.5 * 1.9) - 0.01 * 0.01) ==
+                my_levelset(0, -0.01, 1.9));
+
+          CHECK(0 < my_levelset(0, 0.4, 1.01));
+          CHECK(Approx((1 - 0.5 * 1.01) * (1 - 0.5 * 1.01) - 0.4 * 0.4) ==
+                my_levelset(0, 0.4, 1.01));
+
+          CHECK(0 < my_levelset(-0.4, 0, 1.1));
+          CHECK(Approx((1 - 0.5 * 1.1) * (1 - 0.5 * 1.1) - 0.4 * 0.4) ==
+                my_levelset(-0.4, 0, 1.1));
+        }
+      }
+    }
+  }
+}
+
+} // namespace
diff --git a/src/csg/tests/levelset/transform.t.cpp b/src/csg/tests/levelset/transform.t.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..055bf4fc462c82b97a46c7811746be93e7752f83
--- /dev/null
+++ b/src/csg/tests/levelset/transform.t.cpp
@@ -0,0 +1,140 @@
+#include <memory>
+
+#include "catch2/catch.hpp"
+
+#include <csg.hpp>
+#include <csg_types.hpp>
+
+namespace {
+
+TEST_CASE("Translation", "[Levelset Transform]") {
+
+  csg::Cylinder my_cyl{.radius = 2, .height = 20, .center = true};
+
+  auto my_mat = csg::Mulmatrix();
+  my_mat.rotation[0] = std::array<double, 3>{{1, 0, 0}};
+  my_mat.rotation[1] = std::array<double, 3>{{0, 1, 0}};
+  my_mat.rotation[2] = std::array<double, 3>{{0, 0, 1}};
+  my_mat.translation = {100, 200, 300};
+  my_mat.group.objs.push_back(my_cyl);
+
+  auto my_tree = std::make_shared<csg::Tree>();
+  my_tree->top.objs.push_back(my_mat);
+  csg::CsgIF my_levelset(my_tree);
+
+  SECTION("Outside") {
+    CHECK_FALSE(0 < my_levelset(100, 200, 280));
+    CHECK_FALSE(0 < my_levelset(100, 200, 320));
+    CHECK_FALSE(0 < my_levelset(103, 200, 300));
+    CHECK_FALSE(0 < my_levelset(97, 200, 300));
+    CHECK_FALSE(0 < my_levelset(100, 203, 300));
+    CHECK_FALSE(0 < my_levelset(100, 197, 300));
+  }
+
+  SECTION("Inside") {
+    CHECK(0 < my_levelset(100, 200, 295));
+    CHECK(0 < my_levelset(100, 200, 305));
+    CHECK(0 < my_levelset(101, 200, 300));
+    CHECK(0 < my_levelset(99, 200, 300));
+    CHECK(0 < my_levelset(100, 201, 300));
+    CHECK(0 < my_levelset(100, 199, 300));
+  }
+}
+
+TEST_CASE("90° Rotation around y-axis, rotating z-axis into x-axis",
+          "[Levelset Transform]") {
+
+  csg::Cylinder my_cyl{.radius = 2, .height = 20, .center = false};
+
+  auto my_mat = csg::Mulmatrix();
+  my_mat.rotation[0] = std::array<double, 3>{{0, 0, 1}};
+  my_mat.rotation[1] = std::array<double, 3>{{0, 1, 0}};
+  my_mat.rotation[2] = std::array<double, 3>{{-1, 0, 0}};
+  my_mat.translation = {0, 0, 0};
+  my_mat.group.objs.push_back(my_cyl);
+
+  auto my_tree = std::make_shared<csg::Tree>();
+  my_tree->top.objs.push_back(my_mat);
+  csg::CsgIF my_levelset(my_tree);
+
+  SECTION("near origin end") {
+    CHECK_FALSE(0 < my_levelset(-1, 0, 0));
+    CHECK(0 < my_levelset(1, 0, 0));
+
+    CHECK_FALSE(0 < my_levelset(3, 0, 2.1));
+    CHECK(0 < my_levelset(3, 0, 1.9));
+
+    CHECK_FALSE(0 < my_levelset(3, 0, -2.1));
+    CHECK(0 < my_levelset(3, 0, -1.9));
+  }
+  SECTION("near middle") {
+    CHECK_FALSE(0 < my_levelset(10, 2.1, 0));
+    CHECK(0 < my_levelset(10, 1.9, 0));
+
+    CHECK_FALSE(0 < my_levelset(10, -2.1, 0));
+    CHECK(0 < my_levelset(10, -1.9, 0));
+  }
+  SECTION("near other end") {
+    CHECK_FALSE(0 < my_levelset(21, 0, 0));
+    CHECK(0 < my_levelset(19, 0, 0));
+
+    CHECK_FALSE(0 < my_levelset(18, 0, 2.1));
+    CHECK(0 < my_levelset(18, 0, 1.9));
+
+    CHECK_FALSE(0 < my_levelset(18, 0, -2.1));
+    CHECK(0 < my_levelset(18, 0, -1.9));
+  }
+}
+
+TEST_CASE("Orthogonal rotation + translation of cylinder",
+          "[Levelset Transform]") {
+
+  double radius = 4.5, height = 10;
+  double Cx = 10, Cy = 5, Cz = 5;
+
+  csg::Cylinder my_cyl{.radius = radius, .height = height, .center = true};
+  auto my_mat = csg::Mulmatrix();
+  my_mat.rotation[0] = std::array<double, 3>{{0, 0, 1}};
+  my_mat.rotation[1] = std::array<double, 3>{{0, 1, 0}};
+  my_mat.rotation[2] = std::array<double, 3>{{-1, 0, 0}};
+  my_mat.translation = {Cx, Cy, Cz};
+  my_mat.group.objs.push_back(my_cyl);
+
+  auto my_tree = std::make_shared<csg::Tree>();
+  my_tree->top.objs.push_back(my_mat);
+  csg::CsgIF my_levelset(my_tree);
+
+  SECTION("near origin end") {
+    CHECK_FALSE(0 < my_levelset(Cx - (1.1 * height / 2), Cy, Cz));
+    CHECK(0 < my_levelset(Cx - (0.9 * height / 2), Cy, Cz));
+
+    CHECK_FALSE(0 <
+                my_levelset(Cx - (0.9 * height / 2), Cy, Cz + 1.1 * radius));
+    CHECK(0 < my_levelset(Cx - (0.9 * height / 2), Cy, 0.9 * radius));
+
+    CHECK_FALSE(0 <
+                my_levelset(Cx - (0.9 * height / 2), Cy, Cz - 1.1 * radius));
+    CHECK(0 < my_levelset(Cx - (0.9 * height / 2), Cy, Cz - 0.9 * radius));
+  }
+  SECTION("near middle") {
+    CHECK_FALSE(0 < my_levelset(Cx, Cy + 1.1 * radius, Cz));
+    CHECK(0 < my_levelset(Cx, Cy + 0.9 * radius, Cz));
+
+    CHECK_FALSE(0 < my_levelset(Cx, Cy - 1.1 * radius, Cz));
+    CHECK(0 < my_levelset(Cx, Cy - 0.9 * radius, Cz));
+  }
+  SECTION("near other end") {
+    CHECK_FALSE(0 < my_levelset(Cx + (1.1 * height / 2), Cy, Cz));
+    CHECK(0 < my_levelset(Cx + (0.9 * height / 2), Cy, Cz));
+
+    CHECK_FALSE(0 <
+                my_levelset(Cx + (0.9 * height / 2), Cy, Cz + 1.1 * radius));
+    CHECK(0 < my_levelset(Cx + (0.9 * height / 2), Cy, Cz + 0.9 * radius));
+
+    CHECK_FALSE(0 <
+                my_levelset(Cx + (0.9 * height / 2), Cy, Cz - 1.1 * radius));
+    CHECK(0 < my_levelset(Cx + (0.9 * height / 2), Cy, Cz - 0.9 * radius));
+  }
+}
+
+} // namespace
diff --git a/src/csg/tests/meson.build b/src/csg/tests/meson.build
new file mode 100644
index 0000000000000000000000000000000000000000..4f068140acca0929e72f0ace30caaae0c29d36b6
--- /dev/null
+++ b/src/csg/tests/meson.build
@@ -0,0 +1,15 @@
+unit_exe = executable(
+  'unit_tests',
+  'levelset/boolean.t.cpp',
+  'levelset/primitives.t.cpp',
+  'levelset/transform.t.cpp',
+  'parser/boolean.t.cpp',
+  'parser/main.cpp',
+  'parser/nest.cpp',
+  'parser/other.t.cpp',
+  'parser/primitives.t.cpp',
+  'parser/transform.t.cpp',
+  include_directories: [parser_inc, catch2_inc],
+  link_with: lib_csg_parser
+)
+test('unit_test', unit_exe)
diff --git a/src/csg/tests/parser/boolean.t.cpp b/src/csg/tests/parser/boolean.t.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..45fa06290b0f5c9fcb216641c4e726a8cead1a8a
--- /dev/null
+++ b/src/csg/tests/parser/boolean.t.cpp
@@ -0,0 +1,86 @@
+#include "catch2/catch.hpp"
+
+#include <csg.hpp>
+#include <csg_types.hpp>
+
+// Tests for CSG union(), intersection(), and difference()
+
+TEST_CASE("two shapes", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+sphere(r = 10);
+cube(size = [1,2,3], center=true);
+)");
+  auto sph = std::get<csg::Sphere>(st.top.objs.at(0));
+  CHECK(sph.radius == 10);
+  auto cub = std::get<csg::Cube>(st.top.objs.at(1));
+
+  auto [XX, YY, ZZ] = cub.size;
+  CHECK(XX == 1);
+  CHECK(YY == 2);
+  CHECK(ZZ == 3);
+}
+
+TEST_CASE("union", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+union() {
+  cube(size = [12, 14, 15], center = true);
+  sphere(r = 8);
+}
+)");
+
+  auto un = std::get<csg::Union>(st.top.objs.back());
+  CHECK(st.top.objs.size() == 1);
+
+  auto cub = std::get<csg::Cube>(un.objs.at(0));
+  auto sph = std::get<csg::Sphere>(un.objs.at(1));
+
+  auto [XX, YY, ZZ] = cub.size;
+  CHECK(XX == 12);
+  CHECK(YY == 14);
+  CHECK(ZZ == 15);
+  CHECK(sph.radius == 8);
+}
+
+TEST_CASE("intersection", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+intersection() {
+  cube(size = [15, 15, 15], center = true);
+  sphere(r = 10);
+}
+)");
+
+  CHECK(st.top.objs.size() == 1);
+  auto ints = std::get<csg::Intersection>(st.top.objs.back());
+  CHECK(ints.objs.size() == 2);
+
+  auto cub = std::get<csg::Cube>(ints.objs.at(0));
+
+  auto sph = std::get<csg::Sphere>(ints.objs.at(1));
+
+  auto [XX, YY, ZZ] = cub.size;
+  CHECK(XX == 15);
+  CHECK(YY == 15);
+  CHECK(ZZ == 15);
+  CHECK(sph.radius == 10);
+}
+
+TEST_CASE("difference", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+difference() {
+  cube(size = [12, 12, 12], center = true);
+  sphere(r = 8);
+}
+)");
+
+  auto diff = std::get<csg::Difference>(st.top.objs.back());
+  CHECK(st.top.objs.size() == 1);
+
+  auto cub = std::get<csg::Cube>(*diff.first_obj);
+  auto sph = std::get<csg::Sphere>(diff.next_objs.objs.at(0));
+
+  auto [XX, YY, ZZ] = cub.size;
+  CHECK(XX == 12);
+  CHECK(YY == 12);
+  CHECK(ZZ == 12);
+  CHECK(sph.radius == 8);
+}
diff --git a/tests/main.cpp b/src/csg/tests/parser/main.cpp
similarity index 100%
rename from tests/main.cpp
rename to src/csg/tests/parser/main.cpp
diff --git a/src/csg/tests/parser/nest.cpp b/src/csg/tests/parser/nest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..60a28a3965d838e38755447664e21b136b072971
--- /dev/null
+++ b/src/csg/tests/parser/nest.cpp
@@ -0,0 +1,35 @@
+#include "catch2/catch.hpp"
+
+#include <csg.hpp>
+#include <csg_types.hpp>
+
+TEST_CASE("matmul in matmul", "[csg]") {
+  auto maybe_st = csg::parse_csg(R"(
+multmatrix(
+[
+[1, 0, 0, 1],
+[0, 1, 0, 2],
+[0, 0, 1, 3],
+[0, 0, 0, 1]
+]
+) {
+    cylinder(h = 2, r = 10, center=true);
+multmatrix(
+[
+[1, 0, 0, 4],
+[0, 1, 0, 5],
+[0, 0, 1, 6],
+[0, 0, 0, 1]
+]
+) {
+    cylinder(h = 2, r = 10, center=true);
+}}
+)");
+  CHECK(maybe_st != nullptr);
+  auto st = maybe_st.get();
+  auto mat = std::get<csg::Mulmatrix>(st->top.objs.back());
+  auto mat2 = std::get<csg::Mulmatrix>(mat.group.objs.back());
+  auto cyl = std::get<csg::Cylinder>(mat2.group.objs.at(0));
+  CHECK(cyl.height == 2);
+  CHECK(cyl.radius == 10);
+}
diff --git a/src/csg/tests/parser/other.t.cpp b/src/csg/tests/parser/other.t.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fee47796991b5acae4dffb147204aac151c739d3
--- /dev/null
+++ b/src/csg/tests/parser/other.t.cpp
@@ -0,0 +1,49 @@
+#include "catch2/catch.hpp"
+
+#include <csg.hpp>
+#include <csg_types.hpp>
+
+// Tests that don't fit elsewhere
+
+//  group() same as union()
+
+TEST_CASE("group", "[csg]") {
+  CHECK(csg::parse_csg(R"(
+group() {
+   cylinder($fn = 0, $fa = 5, $fs = 0.1, h = 20, r1 = 5, r2 = 5, center = true);
+}
+)"));
+}
+
+//  render() should be accepted, but ignored
+TEST_CASE("render", "[csg]") {
+  CHECK(csg::parse_csg(R"(
+render(convexity=2) {
+   cylinder($fn = 0, $fa = 5, $fs = 0.1, h = 20, r1 = 5, r2 = 5, center = true);
+}
+)"));
+}
+
+//  First Basic Example when running OpenSCAD
+
+TEST_CASE("OpenSCAD Basic Example", "[csg]") {
+  CHECK(csg::parse_csg(R"(
+multmatrix([[1, 0, 0, -24], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) {
+  union() {
+    cube(size = [15, 15, 15], center = true);
+    sphere($fn = 0, $fa = 12, $fs = 2, r = 10);
+  }
+}
+intersection() {
+  cube(size = [15, 15, 15], center = true);
+  sphere($fn = 0, $fa = 12, $fs = 2, r = 10);
+}
+multmatrix([[1, 0, 0, 24], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) {
+  difference() {
+    cube(size = [15, 15, 15], center = true);
+    sphere(r = 10);
+  }
+}
+group();
+)"));
+}
diff --git a/src/csg/tests/parser/primitives.t.cpp b/src/csg/tests/parser/primitives.t.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..03493ba7192b3b30ebebd56bd8e039a19ea2c654
--- /dev/null
+++ b/src/csg/tests/parser/primitives.t.cpp
@@ -0,0 +1,44 @@
+#include "catch2/catch.hpp"
+
+#include <csg.hpp>
+#include <csg_types.hpp>
+
+// Tests for the primitives on their own
+
+TEST_CASE("cylinder", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+cylinder(h = 2, r = 10, center=true);
+)");
+  auto cyl = std::get<csg::Cylinder>(st.top.objs.at(0));
+  CHECK(cyl.radius == 10);
+  CHECK(cyl.height == 2);
+}
+
+TEST_CASE("cube", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+cube(size = [1,2,3], center=true);
+)");
+  auto cub = std::get<csg::Cube>(st.top.objs.at(0));
+  auto [XX, YY, ZZ] = cub.size;
+  CHECK(XX == 1);
+  CHECK(YY == 2);
+  CHECK(ZZ == 3);
+}
+
+TEST_CASE("sphere", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+sphere(r = 10);
+)");
+  auto sph = std::get<csg::Sphere>(st.top.objs.at(0));
+  CHECK(sph.radius == 10);
+}
+
+TEST_CASE("semicolon optional", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+sphere(r = 10)
+)");
+  auto sph = std::get<csg::Sphere>(st.top.objs.at(0));
+  CHECK(sph.radius == 10);
+}
+
+// TODO: polyhedron()
diff --git a/src/csg/tests/parser/transform.t.cpp b/src/csg/tests/parser/transform.t.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d1fce41662b1376860b710e26547e0d92f743b2b
--- /dev/null
+++ b/src/csg/tests/parser/transform.t.cpp
@@ -0,0 +1,133 @@
+#include "catch2/catch.hpp"
+
+#include <csg.hpp>
+#include <csg_types.hpp>
+
+//  color() should be accepted, but ignored
+TEST_CASE("color", "[csg]") {
+  CHECK(csg::parse_csg(R"(
+color([0, 1, 0, 1]) {
+ cylinder($fn = 0, $fa = 5, $fs = 0.1, h = 20, r1 = 5, r2 = 5, center = true);
+}
+)"));
+}
+
+TEST_CASE("one shape in matmul", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+multmatrix(
+[
+[1, 0, 0, 0.0020],
+[0, 1, 0, 0.0005],
+[0, 0, 1, 0.0005],
+[0, 0, 0, 1]
+]
+) {
+    cylinder(h = 2, r = 10, center=true);
+}
+)");
+  auto mat = std::get<csg::Mulmatrix>(st.top.objs.back());
+  CHECK(mat.group.objs.size() == 1);
+  auto cyl = std::get<csg::Cylinder>(mat.group.objs.at(0));
+  CHECK(cyl.height == 2);
+  CHECK(cyl.radius == 10);
+}
+
+TEST_CASE("two shapes in matmul", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+multmatrix(
+[
+[1, 0, 0, 0.0020],
+[0, 1, 0, 0.0005],
+[0, 0, 1, 0.0005],
+[0, 0, 0, 1]
+]
+) {
+    cylinder(h = 2, r = 10, center=true);
+    sphere(r = 10);
+}
+)");
+  auto mat = std::get<csg::Mulmatrix>(st.top.objs.back());
+  auto cyl = std::get<csg::Cylinder>(mat.group.objs.at(0));
+  auto sph = std::get<csg::Sphere>(mat.group.objs.at(1));
+  CHECK(cyl.height == 2);
+  CHECK(cyl.radius == 10);
+  CHECK(sph.radius == 10);
+}
+
+TEST_CASE("two matmuls", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+multmatrix(
+[
+[1, 0, 0, 0.0020],
+[0, 1, 0, 0.0005],
+[0, 0, 1, 0.0005],
+[0, 0, 0, 1]
+]
+) {
+    cylinder(h = 2, r = 10, center=true);
+    sphere(r = 11);
+}
+multmatrix(
+[
+[1, 0, 0, 0.0020],
+[0, 1, 0, 0.0005],
+[0, 0, 1, 0.0005],
+[0, 0, 0, 1]
+]
+) {
+    cube(size = [1,2,3], center=true);
+    cylinder(h=4, r1=1, r2=2, center=true);
+}
+)");
+  auto mat = std::get<csg::Mulmatrix>(st.top.objs.at(0));
+  auto cyl = std::get<csg::Cylinder>(mat.group.objs.at(0));
+  auto sph = std::get<csg::Sphere>(mat.group.objs.at(1));
+  CHECK(cyl.height == 2);
+  CHECK(cyl.radius == 10);
+  CHECK(sph.radius == 11);
+
+  auto mat2 = std::get<csg::Mulmatrix>(st.top.objs.at(1));
+  auto cube = std::get<csg::Cube>(mat2.group.objs.at(0));
+  auto cone = std::get<csg::Cone>(mat2.group.objs.at(1));
+  auto [XX, YY, ZZ] = cube.size;
+  CHECK(XX == 1);
+  CHECK(YY == 2);
+  CHECK(ZZ == 3);
+  CHECK(cone.height == 4);
+  CHECK(cone.radius1 == 1);
+  CHECK(cone.radius2 == 2);
+}
+
+TEST_CASE("two shapes one matmul", "[csg]") {
+  auto st = *csg::parse_csg(R"(
+cylinder(h = 2, r = 10, center=true);
+sphere(r = 11);
+multmatrix(
+[
+[1, 0, 0, 0.0020],
+[0, 1, 0, 0.0005],
+[0, 0, 1, 0.0005],
+[0, 0, 0, 1]
+]
+) {
+    cube(size = [1,2,3], center=true);
+    cylinder(h=4, r1=1, r2=2, center=true);
+}
+)");
+  auto cyl = std::get<csg::Cylinder>(st.top.objs.at(0));
+  auto sph = std::get<csg::Sphere>(st.top.objs.at(1));
+  auto mat = std::get<csg::Mulmatrix>(st.top.objs.at(2));
+  CHECK(cyl.height == 2);
+  CHECK(cyl.radius == 10);
+  CHECK(sph.radius == 11);
+
+  auto cube = std::get<csg::Cube>(mat.group.objs.at(0));
+  auto cone = std::get<csg::Cone>(mat.group.objs.at(1));
+  auto [XX, YY, ZZ] = cube.size;
+  CHECK(XX == 1);
+  CHECK(YY == 2);
+  CHECK(ZZ == 3);
+  CHECK(cone.height == 4);
+  CHECK(cone.radius1 == 1);
+  CHECK(cone.radius2 == 2);
+}
diff --git a/src/geometry.cpp b/src/inputs/geometry.cpp
similarity index 100%
rename from src/geometry.cpp
rename to src/inputs/geometry.cpp
diff --git a/src/inputs.hpp b/src/inputs/inputs.hpp
similarity index 100%
rename from src/inputs.hpp
rename to src/inputs/inputs.hpp
diff --git a/src/mesh.cpp b/src/inputs/mesh.cpp
similarity index 100%
rename from src/mesh.cpp
rename to src/inputs/mesh.cpp
diff --git a/src/inputs/meson.build b/src/inputs/meson.build
new file mode 100644
index 0000000000000000000000000000000000000000..a435c0ac4fd146e7ecf97a3632d292691b8c2bf4
--- /dev/null
+++ b/src/inputs/meson.build
@@ -0,0 +1,11 @@
+lib_inputs_parser = static_library(
+  'inputs-parser',
+  'geometry.cpp',
+  'mesh.cpp',
+  'parser.cpp',
+  'solver.cpp',
+  'time.cpp',
+  include_directories: tao_inc,
+  install : true)
+
+subdir('tests')
diff --git a/src/parser.cpp b/src/inputs/parser.cpp
similarity index 100%
rename from src/parser.cpp
rename to src/inputs/parser.cpp
diff --git a/src/solver.cpp b/src/inputs/solver.cpp
similarity index 100%
rename from src/solver.cpp
rename to src/inputs/solver.cpp
diff --git a/src/solver.hpp b/src/inputs/solver.hpp
similarity index 100%
rename from src/solver.hpp
rename to src/inputs/solver.hpp
diff --git a/src/inputs/tests/main.cpp b/src/inputs/tests/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fc72584f1b1bd7c8bc9414ccee5fd417b547cae5
--- /dev/null
+++ b/src/inputs/tests/main.cpp
@@ -0,0 +1,2 @@
+#define CATCH_CONFIG_MAIN // Catch main() - only add to one cpp file
+#include "catch2/catch.hpp"
diff --git a/src/inputs/tests/meson.build b/src/inputs/tests/meson.build
new file mode 100644
index 0000000000000000000000000000000000000000..74f0aa91a640b4e78f172f2575977308cf691c48
--- /dev/null
+++ b/src/inputs/tests/meson.build
@@ -0,0 +1,9 @@
+unit_exe = executable(
+  'unit_tests',
+  'parser.t.cpp',
+  'solver.t.cpp',
+  'main.cpp',
+  include_directories: [parser_inc, catch2_inc],
+  link_with: lib_inputs_parser
+)
+test('unit_test', unit_exe)
diff --git a/tests/parser.t.cpp b/src/inputs/tests/parser.t.cpp
similarity index 100%
rename from tests/parser.t.cpp
rename to src/inputs/tests/parser.t.cpp
diff --git a/tests/solver.t.cpp b/src/inputs/tests/solver.t.cpp
similarity index 100%
rename from tests/solver.t.cpp
rename to src/inputs/tests/solver.t.cpp
diff --git a/src/time.cpp b/src/inputs/time.cpp
similarity index 100%
rename from src/time.cpp
rename to src/inputs/time.cpp
diff --git a/src/meson.build b/src/meson.build
deleted file mode 100644
index abe80e4e70b76dec37f7a8103d161d2e5620548e..0000000000000000000000000000000000000000
--- a/src/meson.build
+++ /dev/null
@@ -1,12 +0,0 @@
-lib_parser = static_library('mfix-parser',
-                            'geometry.cpp',
-                            'mesh.cpp',
-                            'parser.cpp',
-                            'solver.cpp',
-                            'time.cpp',
-                            include_directories: tao_inc,
-                            install : true)
-
-executable('mfix-parser',
-           'main.cpp',
-           link_with: lib_parser)
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
deleted file mode 100644
index ddb31a995803f0961c9ac2f4cc641f8bcdea1e2c..0000000000000000000000000000000000000000
--- a/tests/CMakeLists.txt
+++ /dev/null
@@ -1,46 +0,0 @@
-################################################################################
-#  Parser Tests
-################################################################################
-
-add_executable(unit_tests EXCLUDE_FROM_ALL
-  inputs.t.cpp
-  parser.t.cpp
-  main.cpp
-  solver.t.cpp
-  )
-target_link_libraries(unit_tests parser)
-
-target_include_directories(unit_tests PRIVATE
-  ${CMAKE_SOURCE_DIR}/subprojects/PEGTL/include
-  ${CMAKE_SOURCE_DIR}/src
-  )
-
-SET(catch2_dir ${PROJECT_SOURCE_DIR}/subprojects/Catch2)
-
-find_package(Git QUIET)
-if(GIT_FOUND AND EXISTS "${PROJECT_SOURCE_DIR}/.git")
-  # Update submodules as needed
-  option(GIT_SUBMODULE "Check submodules during build" ON)
-  if(GIT_SUBMODULE)
-    message(STATUS "Submodule update")
-    execute_process(COMMAND ${GIT_EXECUTABLE} submodule update --init
-      WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
-      RESULT_VARIABLE GIT_SUBMOD_RESULT)
-    if(NOT GIT_SUBMOD_RESULT EQUAL "0")
-      message(FATAL_ERROR "git submodule update --init failed with ${GIT_SUBMOD_RESULT}, please checkout submodules")
-    endif()
-  endif()
-endif()
-
-if(NOT EXISTS "${catch2_dir}/CMakeLists.txt")
-  message(FATAL_ERROR "The submodules were not downloaded! GIT_SUBMODULE was turned off or failed. Please update submodules and try again.")
-endif()
-
-add_subdirectory(${catch2_dir} ${CMAKE_CURRENT_BINARY_DIR}/catch_build)
-list(APPEND CMAKE_MODULE_PATH "${catch2_dir}/contrib/")
-
-target_link_libraries(unit_tests Catch2::Catch2)
-
-include(CTest)
-include(Catch)
-catch_discover_tests(unit_tests)
diff --git a/tests/meson.build b/tests/meson.build
deleted file mode 100644
index 2b302c396ee2de167025908845427c7ff4cd7da1..0000000000000000000000000000000000000000
--- a/tests/meson.build
+++ /dev/null
@@ -1,8 +0,0 @@
-unit_exe = executable('unit_tests',
-                      'parser.t.cpp',
-                      'solver.t.cpp',
-                      'main.cpp',
-                      include_directories: [parser_inc, catch2_inc],
-                      link_with: lib_parser
-                     )
-test('unit_test', unit_exe)