diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0c185556..87d9d18f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -9,6 +9,7 @@ option(XDG_ENABLE_EMBREE  "Enable support for the Embree ray tracing library" ON
 option(XDG_ENABLE_GPRT    "Enable support for the GPRT ray tracing library"  OFF)
 option(XDG_BUILD_TESTS    "Enable C++ unit testing"                           ON)
 option(XDG_BUILD_TOOLS    "Enable tools and miniapps"                         ON)
+option(XDG_BUILD_EXAMPLES "Enable C++ examples"                               ON)
 
 # Set version numbers
 set(XDG_VERSION_MAJOR 0)
@@ -26,6 +27,7 @@ set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 # Set module path
 set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules)
 
+
 #===============================================================================
 # libMesh
 #===============================================================================
@@ -69,6 +71,13 @@ if (XDG_ENABLE_EMBREE)
   message(STATUS "Found Embree ${EMBREE_VERSION} at ${EMBREE_ROOT_DIR}")
 endif()
 
+#============================================================================
+# If we want to build examples
+#============================================================================
+if(XDG_BUILD_EXAMPLES)
+  add_subdirectory(examples)
+endif ()
+
 #===============================================================================
 # Update git submodules as needed (borrowed from OpenMC)
 #===============================================================================
@@ -323,7 +332,7 @@ if (XDG_ENABLE_GPRT)
   target_link_options(xdg PRIVATE -Wl,--unresolved-symbols=ignore-in-shared-libs)
 endif()
 
-target_link_libraries(xdg PRIVATE fmt::fmt)
+target_link_libraries(xdg PUBLIC fmt::fmt)
 
 # ==========================
 # Link ray tracing libraries
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
new file mode 100644
index 00000000..3cc87d28
--- /dev/null
+++ b/examples/CMakeLists.txt
@@ -0,0 +1,8 @@
+# =============================================================================
+# Examples for how to use XDG C++ APIs
+# =============================================================================
+
+add_executable(element_volume_estimation element_volume_estimation.cpp)
+target_link_libraries(element_volume_estimation PRIVATE xdg)
+
+
diff --git a/examples/element_volume_estimation.cpp b/examples/element_volume_estimation.cpp
new file mode 100644
index 00000000..af9d1003
--- /dev/null
+++ b/examples/element_volume_estimation.cpp
@@ -0,0 +1,52 @@
+#include <iostream>
+#include <iomanip>
+#include <random>
+
+#include "xdg/xdg.h"
+
+int main(int argc, char* argv[])
+{
+    if (argc < 2) {
+        std::cerr << "Usage: " << argv[0] << " <mesh_file>\n";
+        return 1;
+    }
+    const std::string mesh_file = argv[1];
+
+    std::shared_ptr<xdg::XDG> xdg = xdg::XDG::create();
+    const auto& mesh_manager = xdg->mesh_manager();
+    mesh_manager->load_file(mesh_file);
+    mesh_manager->init();
+    xdg->prepare_raytracer();
+
+    xdg::BoundingBox bbox = mesh_manager->global_bounding_box();
+    double bbox_volume = (bbox.max_x - bbox.min_x)
+                         * (bbox.max_y - bbox.min_y)
+                         * (bbox.max_z - bbox.min_z);
+
+    std::mt19937 gen(42);
+    std::uniform_real_distribution<double> x_dist(bbox.min_x, bbox.max_x);
+    std::uniform_real_distribution<double> y_dist(bbox.min_y, bbox.max_y);
+    std::uniform_real_distribution<double> z_dist(bbox.min_z, bbox.max_z);
+
+    constexpr size_t n_samples = 100000;
+
+    std::cout << std::setw(10) << "Volume ID"
+              << std::setw(20) << "MC Volume Estimate"
+              << std::setw(20) << "Hits / Samples\n";
+    std::cout << std::string(50, '-') << "\n";
+
+    for (xdg::MeshID volume : mesh_manager->volumes()) {
+        size_t hits = 0;
+        for (size_t i = 0; i < n_samples; ++i) {
+            xdg::Position sample = {x_dist(gen), y_dist(gen), z_dist(gen)};
+            if (xdg->point_in_volume(volume, sample))
+                hits++;
+        }
+        const double mc_volume = bbox_volume * static_cast<double>(hits) / n_samples;
+        std::cout << std::setw(10) << volume
+                  << std::setw(20) << std::fixed << std::setprecision(4) << mc_volume
+                  << std::setw(14) << hits << " / " << n_samples << "\n";
+    }
+
+    return 0;
+}
\ No newline at end of file