diff --git a/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp b/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp
index 79ac89e3182..3e84bdd8cb6 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp
@@ -1,6 +1,6 @@
 #include <CGAL/Cartesian_grid_3.h>
-#include <CGAL/Cartesian_grid_domain.h>
 #include <CGAL/Dual_contouring_3.h>
+#include <CGAL/Isosurfacing_domains.h>
 #include <CGAL/Marching_cubes_3.h>
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/boost/graph/IO/OFF.h>
@@ -60,7 +60,7 @@ int main() {
     };
 
     // create a domain from the grid
-    CGAL::Isosurfacing::Cartesian_grid_domain<Kernel, decltype(cube_gradient)> domain(grid, cube_gradient);
+    auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid, cube_gradient);
 
     // prepare collections for the results
     Point_range points_mc, points_tmc, points_dc;
diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp
index c4ad83a107a..dc06704f2cf 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp
@@ -1,6 +1,7 @@
 #include <CGAL/Cartesian_grid_3.h>
-#include <CGAL/Cartesian_grid_domain.h>
+#include <CGAL/Default_gradients.h>
 #include <CGAL/Dual_contouring_3.h>
+#include <CGAL/Isosurfacing_domains.h>
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/boost/graph/IO/OFF.h>
 
@@ -35,7 +36,9 @@ int main() {
         }
     }
 
-    CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
+    CGAL::Isosurfacing::Explicit_cartesian_grid_gradient<Kernel> gradient(grid);
+
+    auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid, gradient);
 
     Point_range points;
     Polygon_range polygons;
diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp
index f40b5ade29a..9e7118f3286 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp
@@ -1,5 +1,6 @@
+#include <CGAL/Bbox_3.h>
 #include <CGAL/Dual_contouring_3.h>
-#include <CGAL/Implicit_domain.h>
+#include <CGAL/Isosurfacing_domains.h>
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/boost/graph/IO/OFF.h>
 
@@ -36,9 +37,12 @@ int main() {
         return Vector(gx, gy, gz);
     };
 
+    const CGAL::Bbox_3 bbox{-1, -1, -1, 1, 1, 1};
+    const Vector spacing(0.02f, 0.02f, 0.02f);
+
     // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
-    CGAL::Isosurfacing::Implicit_domain<Kernel, decltype(iwp_value), decltype(iwp_gradient)> domain(
-        {-1, -1, -1, 1, 1, 1}, Vector(0.02f, 0.02f, 0.02f), iwp_value, iwp_gradient);
+    auto domain =
+        CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(bbox, spacing, iwp_value, iwp_gradient);
 
     // prepare collections for the result
     Point_range points;
diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp
index 866dd6a9b34..742cb4e9efc 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp
@@ -1,8 +1,9 @@
 #include <CGAL/AABB_face_graph_triangle_primitive.h>
 #include <CGAL/AABB_traits.h>
 #include <CGAL/AABB_tree.h>
+#include <CGAL/Bbox_3.h>
 #include <CGAL/Dual_contouring_3.h>
-#include <CGAL/Implicit_domain.h>
+#include <CGAL/Isosurfacing_domains.h>
 #include <CGAL/Side_of_triangle_mesh.h>
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Surface_mesh.h>
@@ -25,9 +26,9 @@ typedef std::vector<Point> Point_range;
 typedef std::vector<std::vector<std::size_t>> Polygon_range;
 
 int main() {
-    const std::string input_name = CGAL::data_file_path("bunny.off");
-    const Vector grid_spacing(0.005, 0.005, 0.005);
-    const FT offset_value = 0.01;
+    const std::string input_name = CGAL::data_file_path("meshes/triceratops.off");
+    const Vector grid_spacing(0.5, 0.5, 0.5);
+    const FT offset_value = 0.1;
 
     Mesh mesh_input;
     if (!CGAL::IO::read_OFF(input_name, mesh_input)) {
@@ -59,8 +60,8 @@ int main() {
     };
 
     // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
-    CGAL::Isosurfacing::Implicit_domain<Kernel, decltype(mesh_distance), decltype(mesh_normal)> domain(
-        aabb_grid, grid_spacing, mesh_distance, mesh_normal);
+    auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(aabb_grid, grid_spacing,
+                                                                                    mesh_distance, mesh_normal);
 
     Point_range points;
     Polygon_range polygons;
diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp
index 457839aa2c8..6481cad0c01 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp
@@ -1,5 +1,5 @@
 #include <CGAL/Dual_contouring_3.h>
-#include <CGAL/Octree_domain.h>
+#include <CGAL/Isosurfacing_domains.h>
 #include <CGAL/Octree_wrapper.h>
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/boost/graph/IO/OFF.h>
@@ -9,18 +9,13 @@
 
 typedef CGAL::Simple_cartesian<double> Kernel;
 typedef typename Kernel::FT FT;
-typedef typename Kernel::Vector_3 Vector_3;
-typedef typename Kernel::Point_3 Point_3;
+typedef typename Kernel::Vector_3 Vector;
+typedef typename Kernel::Point_3 Point;
 
-typedef std::vector<Point_3> Point_range;
+typedef std::vector<Point> Point_range;
 typedef std::vector<std::vector<std::size_t>> Polygon_range;
 
 typedef CGAL::Isosurfacing::Octree_wrapper<Kernel> Octree_wrapper_;
-typedef CGAL::Isosurfacing::Octree_domain<Kernel> Octree_domain_;
-
-Kernel::FT sphere_function(const Point_3& point) {
-    return std::sqrt(point.x() * point.x() + point.y() * point.y() + point.z() * point.z());
-}
 
 struct Refine_one_eighth {
     std::size_t min_depth_;
@@ -72,22 +67,19 @@ int main() {
     Refine_one_eighth split_predicate(4, 6);
     octree_wrap.refine(split_predicate);
 
-    Octree_domain_ octree_domain(octree_wrap);
+    auto sphere_function = [&](const Point& p) { return std::sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z()); };
 
-    auto lam = [&](const Octree_domain_::Vertex_handle& v) {
-        Point_3 p = octree_domain.position(v);
-        const auto val = sphere_function(p);
-        Vector_3 gradient = p - CGAL::ORIGIN;
-        gradient = gradient / std::sqrt(gradient.squared_length());
-        octree_wrap.value(v) = val;
-        octree_wrap.gradient(v) = gradient;
+    auto sphere_gradient = [&](const Point& p) {
+        const Vector g = p - CGAL::ORIGIN;
+        return g / std::sqrt(g.squared_length());
     };
-    octree_domain.iterate_vertices(lam, CGAL::Sequential_tag());
+
+    auto domain = CGAL::Isosurfacing::create_implicit_octree_domain(octree_wrap, sphere_function, sphere_gradient);
 
     Point_range points;
     Polygon_range polygons;
 
-    CGAL::Isosurfacing::dual_contouring(octree_domain, 0.8, points, polygons);
+    CGAL::Isosurfacing::dual_contouring(domain, 0.8, points, polygons);
 
     CGAL::IO::write_OFF("result.off", points, polygons);
 }
diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp
index a43fb409f48..8a363703a6e 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp
@@ -1,6 +1,6 @@
 
 #include <CGAL/Cartesian_grid_3.h>
-#include <CGAL/Cartesian_grid_domain.h>
+#include <CGAL/Isosurfacing_domains.h>
 #include <CGAL/Marching_cubes_3.h>
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/boost/graph/IO/OFF.h>
@@ -34,7 +34,7 @@ int main() {
     }
 
     // create a domain from the grid
-    CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
+    auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid);
 
     // prepare collections for the result
     Point_range points;
diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp
index 28ded19e538..9d9584907f5 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp
@@ -1,5 +1,6 @@
 
-#include <CGAL/Implicit_domain.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/Isosurfacing_domains.h>
 #include <CGAL/Marching_cubes_3.h>
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/boost/graph/IO/OFF.h>
@@ -12,14 +13,13 @@ typedef std::vector<Point> Point_range;
 typedef std::vector<std::vector<std::size_t>> Polygon_range;
 
 int main() {
-    // distance to the origin
-    auto sphere_function = [](const Point& point) {
-        return std::sqrt(point.x() * point.x() + point.y() * point.y() + point.z() * point.z());
-    };
+    const CGAL::Bbox_3 bbox{-1, -1, -1, 1, 1, 1};
+    const Vector spacing(0.02f, 0.02f, 0.02f);
+
+    auto sphere_function = [&](const Point& p) { return std::sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z()); };
 
     // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
-    CGAL::Isosurfacing::Implicit_domain<Kernel, decltype(sphere_function)> domain(
-        {-1, -1, -1, 1, 1, 1}, Vector(0.02f, 0.02f, 0.02f), sphere_function);  // TODO: this is ugly
+    auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(bbox, spacing, sphere_function);
 
     // prepare collections for the result
     Point_range points;
diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp
index e67b269a17d..0a5e14907a8 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp
@@ -1,5 +1,5 @@
 #include <CGAL/Cartesian_grid_3.h>
-#include <CGAL/Cartesian_grid_domain.h>
+#include <CGAL/Isosurfacing_domains.h>
 #include <CGAL/Marching_cubes_3.h>
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/boost/graph/IO/OFF.h>
@@ -27,7 +27,7 @@ int main() {
     const Grid grid(image);
 
     // create a domain from the grid
-    CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
+    auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid);
 
     // prepare collections for the result
     Point_range points;
diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp
index 851f85b125f..c428546342f 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp
@@ -2,7 +2,7 @@
 #include <CGAL/AABB_traits.h>
 #include <CGAL/AABB_tree.h>
 #include <CGAL/Cartesian_grid_3.h>
-#include <CGAL/Cartesian_grid_domain.h>
+#include <CGAL/Isosurfacing_domains.h>
 #include <CGAL/Marching_cubes_3.h>
 #include <CGAL/Side_of_triangle_mesh.h>
 #include <CGAL/Simple_cartesian.h>
@@ -82,7 +82,7 @@ int main() {
     }
 
     // create a domain from the grid
-    CGAL::Isosurfacing::Cartesian_grid_domain<Kernel> domain(grid);
+    auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid);
 
     // prepare collections for the result
     Point_range points;