diff --git a/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp b/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp
index 5a84eac6cb2..c81aaeddb87 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp
@@ -27,19 +27,19 @@ int main(int, char**)
 {
   // create a Cartesian grid with 7^3 grid points and the bounding box [-1, 1]^3
   const CGAL::Bbox_3 bbox(-1.0, -1.0, -1.0,  1.0, 1.0, 1.0);
-  std::shared_ptr<Grid> grid = std::make_shared<Grid>(7, 7, 7, bbox);
+  Grid grid { 7, 7, 7, bbox };
 
   // calculate the value at all grid points
-  for(std::size_t x=0; x<grid->xdim(); ++x) {
-    for(std::size_t y=0; y<grid->ydim(); ++y) {
-      for(std::size_t z=0; z<grid->zdim(); ++z)
+  for(std::size_t x=0; x<grid.xdim(); ++x) {
+    for(std::size_t y=0; y<grid.ydim(); ++y) {
+      for(std::size_t z=0; z<grid.zdim(); ++z)
       {
-        const FT pos_x = x * grid->get_spacing()[0] + bbox.xmin();
-        const FT pos_y = y * grid->get_spacing()[1] + bbox.ymin();
-        const FT pos_z = z * grid->get_spacing()[2] + bbox.zmin();
+        const FT pos_x = x * grid.get_spacing()[0] + bbox.xmin();
+        const FT pos_y = y * grid.get_spacing()[1] + bbox.ymin();
+        const FT pos_z = z * grid.get_spacing()[2] + bbox.zmin();
 
         // L_inf distance to the origin
-        grid->value(x, y, z) = std::max({std::abs(pos_x), std::abs(pos_y), std::abs(pos_z)});
+        grid.value(x, y, z) = (std::max)({std::abs(pos_x), std::abs(pos_y), std::abs(pos_z)});
       }
     }
   }
@@ -48,7 +48,7 @@ int main(int, char**)
   auto cube_gradient = [](const Point& p)
   {
     // the normal depends on the side of the cube
-    const FT max_value = std::max({std::abs(p.x()), std::abs(p.y()), std::abs(p.z())});
+    const FT max_value = (std::max)({std::abs(p.x()), std::abs(p.y()), std::abs(p.z())});
 
     Vector g(0.0, 0.0, 0.0);
     if(max_value == std::abs(p.x()))
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 87b389bd5d6..94264709b13 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_cartesian_grid.cpp
@@ -21,22 +21,22 @@ int main(int, char**)
 {
   // create bounding box and grid
   const CGAL::Bbox_3 bbox(-1.0, -1.0, -1.0,  1.0, 1.0, 1.0);
-  std::shared_ptr<Grid> grid = std::make_shared<Grid>(30, 30, 30, bbox);
+  Grid grid { 30, 30, 30, bbox };
 
   // compute field values and gradients
-  for(std::size_t x=0; x<grid->xdim(); ++x) {
-    for(std::size_t y=0; y<grid->ydim(); ++y) {
-      for(std::size_t z=0; z<grid->zdim(); ++z)
+  for(std::size_t x=0; x<grid.xdim(); ++x) {
+    for(std::size_t y=0; y<grid.ydim(); ++y) {
+      for(std::size_t z=0; z<grid.zdim(); ++z)
       {
-        const FT pos_x = x * grid->get_spacing()[0] + bbox.xmin();
-        const FT pos_y = y * grid->get_spacing()[1] + bbox.ymin();
-        const FT pos_z = z * grid->get_spacing()[2] + bbox.zmin();
+        const FT pos_x = x * grid.get_spacing()[0] + bbox.xmin();
+        const FT pos_y = y * grid.get_spacing()[1] + bbox.ymin();
+        const FT pos_z = z * grid.get_spacing()[2] + bbox.zmin();
 
         const Vector direction(pos_x, pos_y, pos_z);
         const FT distance = CGAL::approximate_sqrt(direction.squared_length());
 
-        grid->value(x, y, z) = distance;
-        grid->gradient(x, y, z) = direction / distance; // @todo check division / 0
+        grid.value(x, y, z) = distance;
+        grid.gradient(x, y, z) = direction / distance; // @todo check division / 0
       }
     }
   }
@@ -45,7 +45,7 @@ int main(int, char**)
   CGAL::Isosurfacing::Explicit_cartesian_grid_gradient<Kernel> gradient(grid);
 
   // create domain from scalar and gradient fields
-  auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid, gradient);
+  auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain(grid, gradient);
 
   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 c56c22bd25d..e8e9f5f81f8 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_octree.cpp
@@ -19,7 +19,7 @@ using Point = typename Kernel::Point_3;
 using Point_range = std::vector<Point>;
 using Polygon_range = std::vector<std::vector<std::size_t> >;
 
-using Octree_wrapper_ = CGAL::Isosurfacing::internal::Octree_wrapper<Kernel>;
+using Octree_wrapper = CGAL::Isosurfacing::internal::Octree_wrapper<Kernel>;
 
 struct Refine_one_eighth
 {
@@ -28,11 +28,11 @@ struct Refine_one_eighth
 
   std::size_t octree_dim_;
 
-  Octree_wrapper_::Uniform_coords uniform_coordinates(const Octree_wrapper_::Octree::Node& node) const
+  Octree_wrapper::Uniform_coords uniform_coordinates(const Octree_wrapper::Octree::Node& node) const
   {
     auto coords = node.global_coordinates();
     const std::size_t depth_factor = std::size_t(1) << (max_depth_ - node.depth());
-    for(int i=0; i < Octree_wrapper_::Octree::Node::Dimension::value; ++i)
+    for(int i=0; i < Octree_wrapper::Octree::Node::Dimension::value; ++i)
       coords[i] *= uint32_t(depth_factor);
 
     return coords;
@@ -46,7 +46,7 @@ struct Refine_one_eighth
     octree_dim_ = std::size_t(1) << max_depth_;
   }
 
-  bool operator()(const Octree_wrapper_::Octree::Node& n) const
+  bool operator()(const Octree_wrapper::Octree::Node& n) const
   {
     // n.depth()
     if(n.depth() < min_depth_)
@@ -73,10 +73,10 @@ struct Refine_one_eighth
 int main(int, char**)
 {
   const CGAL::Bbox_3 bbox(-1., -1., -1., 1., 1., 1.);
-  std::shared_ptr<Octree_wrapper_> octree_wrap = std::make_shared<Octree_wrapper_>(bbox);
+  Octree_wrapper octree_wrap { bbox };
 
   Refine_one_eighth split_predicate(3, 4);
-  octree_wrap->refine(split_predicate);
+  octree_wrap.refine(split_predicate);
 
   auto sphere_function = [&](const Point& p)
   {
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 27106849c63..d929374579e 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
@@ -23,25 +23,25 @@ int main(int, char**)
 {
   // create a Cartesian grid with 100^3 grid points and the bounding box [-1, 1]^3
   const CGAL::Bbox_3 bbox(-1.0, -1.0, -1.0, 1.0, 1.0, 1.0);
-  std::shared_ptr<Grid> grid = std::make_shared<Grid>(50, 50, 50, bbox);
+  Grid grid { 50, 50, 50, bbox };
 
   // compute and store function values at all grid points
-  for(std::size_t x=0; x<grid->xdim(); ++x) {
-    for(std::size_t y=0; y<grid->ydim(); ++y) {
-      for(std::size_t z=0; z<grid->zdim(); ++z)
+  for(std::size_t x=0; x<grid.xdim(); ++x) {
+    for(std::size_t y=0; y<grid.ydim(); ++y) {
+      for(std::size_t z=0; z<grid.zdim(); ++z)
       {
-        const FT pos_x = x * grid->get_spacing()[0] + bbox.xmin();
-        const FT pos_y = y * grid->get_spacing()[1] + bbox.ymin();
-        const FT pos_z = z * grid->get_spacing()[2] + bbox.zmin();
+        const FT pos_x = x * grid.get_spacing()[0] + bbox.xmin();
+        const FT pos_y = y * grid.get_spacing()[1] + bbox.ymin();
+        const FT pos_z = z * grid.get_spacing()[2] + bbox.zmin();
 
         // Euclidean distance to the origin
-        grid->value(x, y, z) = std::sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z);
+        grid.value(x, y, z) = std::sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z);
       }
     }
   }
 
   // create a domain from the grid
-  auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid);
+  auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain(grid);
 
   // 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 96ae599a5cc..952f9bdda45 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp
@@ -29,10 +29,10 @@ int main(int, char**)
   }
 
   // convert image to a Cartesian grid
-  std::shared_ptr<Grid> grid = std::make_shared<Grid>(image);
+  Grid grid { image };
 
   // create a domain from the grid
-  auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid);
+  auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain(grid);
 
   // prepare collections for the output indexed mesh
   Point_range points;
diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_signed_mesh_offset.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_signed_mesh_offset.cpp
index 1e87e345cc8..e067c7310e6 100644
--- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_signed_mesh_offset.cpp
+++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_signed_mesh_offset.cpp
@@ -67,30 +67,30 @@ int main(int, char**)
   CGAL::Side_of_triangle_mesh<Mesh, CGAL::GetGeomTraits<Mesh>::type> sotm(mesh_input);
 
   // create grid
-  std::shared_ptr<Grid> grid = std::make_shared<Grid>(n_voxels, n_voxels, n_voxels, aabb_grid);
+  Grid grid { n_voxels, n_voxels, n_voxels, aabb_grid };
 
-  for(std::size_t z=0; z<grid->zdim(); ++z) {
-    for(std::size_t y=0; y<grid->ydim(); ++y) {
-      for(std::size_t x=0; x<grid->xdim(); ++x)
+  for(std::size_t z=0; z<grid.zdim(); ++z) {
+    for(std::size_t y=0; y<grid.ydim(); ++y) {
+      for(std::size_t x=0; x<grid.xdim(); ++x)
       {
-        const FT pos_x = x * grid->get_spacing()[0] + grid->get_bbox().xmin();
-        const FT pos_y = y * grid->get_spacing()[1] + grid->get_bbox().ymin();
-        const FT pos_z = z * grid->get_spacing()[2] + grid->get_bbox().zmin();
+        const FT pos_x = x * grid.get_spacing()[0] + grid.get_bbox().xmin();
+        const FT pos_y = y * grid.get_spacing()[1] + grid.get_bbox().ymin();
+        const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin();
         const Point p(pos_x, pos_y, pos_z);
 
         // compute unsigned distance to input mesh
-        grid->value(x, y, z) = distance_to_mesh(tree, p);
+        grid.value(x, y, z) = distance_to_mesh(tree, p);
 
         // sign distance so that it is negative inside the mesh
         const bool is_inside = (sotm(p) == CGAL::ON_BOUNDED_SIDE);
         if(is_inside)
-          grid->value(x, y, z) *= -1.0;
+          grid.value(x, y, z) *= -1.0;
       }
     }
   }
 
   // create domain from the grid
-  auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid);
+  auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain(grid);
 
   // containers for output indexed triangle soup
   Point_range points;
diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h
index dfc09df34d7..97de20f5076 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_cartesian_grid_domain.h
@@ -63,7 +63,7 @@ using Explicit_cartesian_grid_domain =
 template <typename GeomTraits,
           typename Gradient_ = Zero_gradient>
 Explicit_cartesian_grid_domain<GeomTraits, Gradient_>
-create_explicit_cartesian_grid_domain(const std::shared_ptr<Cartesian_grid_3<GeomTraits> > grid,
+create_explicit_cartesian_grid_domain(const Cartesian_grid_3<GeomTraits>& grid,
                                       const Gradient_& gradient = Gradient_())
 {
   using Domain = Explicit_cartesian_grid_domain<GeomTraits, Gradient_>;
@@ -73,21 +73,20 @@ create_explicit_cartesian_grid_domain(const std::shared_ptr<Cartesian_grid_3<Geo
   using Function = typename Domain::Function;
   using Gradient = typename Domain::Gradient;
 
-  const std::size_t size_i = grid->xdim();
-  const std::size_t size_j = grid->ydim();
-  const std::size_t size_k = grid->zdim();
+  const std::size_t size_i = grid.xdim();
+  const std::size_t size_j = grid.ydim();
+  const std::size_t size_k = grid.zdim();
 
-  const Bbox_3& bbox = grid->get_bbox();
+  const Bbox_3& bbox = grid.get_bbox();
   const typename GeomTraits::Vector_3 offset(bbox.xmin(), bbox.ymin(), bbox.zmin());
-  const typename GeomTraits::Vector_3 spacing = grid->get_spacing();
+  const typename GeomTraits::Vector_3 spacing = grid.get_spacing();
 
-  // create copies as shared_ptr for safe memory management
-  const Topology topo = std::make_shared<Topology::element_type>(size_i, size_j, size_k);
-  const Geometry geom = std::make_shared<Geometry::element_type>(offset, spacing);
-  const Function func = grid;
-  const Gradient grad = std::make_shared<Gradient::element_type>(gradient);
+  const Topology topo { size_i, size_j, size_k };
+  const Geometry geom { offset, spacing };
+  const Function func { grid };
+  const Gradient grad { gradient };
 
-  return Domain(topo, geom, func, grad);
+  return { topo, geom, func, grad };
 }
 
 } // namespace Isosurfacing
diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_cartesian_grid_gradient.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_cartesian_grid_gradient.h
index 60d709cc70b..231694e0775 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_cartesian_grid_gradient.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Explicit_cartesian_grid_gradient.h
@@ -41,7 +41,7 @@ public:
   using Point = typename Geom_traits::Point_3;
   using Vector = typename Geom_traits::Vector_3;
 
-  using Grid = std::shared_ptr<Cartesian_grid_3<Geom_traits> >;
+  using Grid = Cartesian_grid_3<Geom_traits>;
 
 public:
   /**
@@ -65,20 +65,20 @@ public:
   Vector operator()(const Point& point) const
   {
     // trilinear interpolation of stored gradients
-    const Bbox_3& bbox = grid->get_bbox();
-    const Vector& spacing = grid->get_spacing();
+    const Bbox_3& bbox = grid.get_bbox();
+    const Vector& spacing = grid.get_spacing();
 
     // calculate min index including border case
     std::size_t min_i = (point.x() - bbox.xmin()) / spacing.x();
     std::size_t min_j = (point.y() - bbox.ymin()) / spacing.y();
     std::size_t min_k = (point.z() - bbox.zmin()) / spacing.z();
-    if(min_i == grid->xdim() - 1)
+    if(min_i == grid.xdim() - 1)
       --min_i;
 
-    if(min_j == grid->ydim() - 1)
+    if(min_j == grid.ydim() - 1)
       --min_j;
 
-    if(min_k == grid->zdim() - 1)
+    if(min_k == grid.zdim() - 1)
       --min_k;
 
     // calculate coordinates of min index
@@ -92,14 +92,14 @@ public:
     const FT f_k = (point.z() - min_z) / spacing.z();
 
     // read the gradient at all 8 corner points
-    const Vector g000 = grid->gradient(min_i + 0, min_j + 0, min_k + 0);
-    const Vector g001 = grid->gradient(min_i + 0, min_j + 0, min_k + 1);
-    const Vector g010 = grid->gradient(min_i + 0, min_j + 1, min_k + 0);
-    const Vector g011 = grid->gradient(min_i + 0, min_j + 1, min_k + 1);
-    const Vector g100 = grid->gradient(min_i + 1, min_j + 0, min_k + 0);
-    const Vector g101 = grid->gradient(min_i + 1, min_j + 0, min_k + 1);
-    const Vector g110 = grid->gradient(min_i + 1, min_j + 1, min_k + 0);
-    const Vector g111 = grid->gradient(min_i + 1, min_j + 1, min_k + 1);
+    const Vector g000 = grid.gradient(min_i + 0, min_j + 0, min_k + 0);
+    const Vector g001 = grid.gradient(min_i + 0, min_j + 0, min_k + 1);
+    const Vector g010 = grid.gradient(min_i + 0, min_j + 1, min_k + 0);
+    const Vector g011 = grid.gradient(min_i + 0, min_j + 1, min_k + 1);
+    const Vector g100 = grid.gradient(min_i + 1, min_j + 0, min_k + 0);
+    const Vector g101 = grid.gradient(min_i + 1, min_j + 0, min_k + 1);
+    const Vector g110 = grid.gradient(min_i + 1, min_j + 1, min_k + 0);
+    const Vector g111 = grid.gradient(min_i + 1, min_j + 1, min_k + 1);
 
     // interpolate along all axes by weighting the corner points
     const Vector g0 = g000 * (1 - f_i) * (1 - f_j) * (1 - f_k);
@@ -116,7 +116,7 @@ public:
   }
 
 private:
-  const Grid grid;
+  const Grid& grid;
 };
 
 } // namespace Isosurfacing
diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_cartesian_grid_domain.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_cartesian_grid_domain.h
index ff06eccf050..77a146cd848 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_cartesian_grid_domain.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_cartesian_grid_domain.h
@@ -87,7 +87,7 @@ create_implicit_cartesian_grid_domain(const Bbox_3& bbox,
   using Geometry = typename Domain::Geometry;
   using Function = typename Domain::Function;
   using Gradient = typename Domain::Gradient;
-  using Point_function = typename Function::element_type::Point_function;
+  using Point_function = PointFunction;
 
   // calculate grid dimensions
   const std::size_t size_i = std::ceil(bbox.x_span() / spacing.x()) + 1;
@@ -96,14 +96,13 @@ create_implicit_cartesian_grid_domain(const Bbox_3& bbox,
 
   const typename GeomTraits::Vector_3 offset(bbox.xmin(), bbox.ymin(), bbox.zmin());
 
-  // create copies as shared_ptr for safe memory management
-  const Topology topo = std::make_shared<Topology::element_type>(size_i, size_j, size_k);
-  const Geometry geom = std::make_shared<Geometry::element_type>(offset, spacing);
-  const Point_function point_func = std::make_shared<Point_function::element_type>(point_function);
-  const Function func = std::make_shared<Function::element_type>(geom, point_func);
-  const Gradient grad = std::make_shared<Gradient::element_type>(gradient);
+  const Topology topo { size_i, size_j, size_k };
+  const Geometry geom { offset, spacing };
+  const Point_function point_func { point_function };
+  const Function func { geom, point_func };
+  const Gradient grad { gradient };
 
-  return Domain(topo, geom, func, grad);
+  return { topo, geom, func, grad };
 }
 
 } // namespace Isosurfacing
diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_octree_domain.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_octree_domain.h
index 1dc9d8efe43..5c2dfd0eed2 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_octree_domain.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/Implicit_octree_domain.h
@@ -73,7 +73,7 @@ template <typename GeomTraits,
           typename PointFunction,
           typename Gradient_ = Zero_gradient>
 Implicit_octree_domain<GeomTraits, PointFunction, Gradient_>
-create_implicit_octree_domain(const std::shared_ptr<internal::Octree_wrapper<GeomTraits> > octree,
+create_implicit_octree_domain(const internal::Octree_wrapper<GeomTraits>& octree,
                               const PointFunction& point_function,
                               const Gradient_& gradient = Gradient_())
 {
@@ -83,17 +83,16 @@ create_implicit_octree_domain(const std::shared_ptr<internal::Octree_wrapper<Geo
   using Geometry = typename Domain::Geometry;
   using Function = typename Domain::Function;
   using Gradient = typename Domain::Gradient;
-  using Point_function = typename Function::element_type::Point_function;
-  using Octree = typename Topology::element_type::Octree;
+  using Point_function = PointFunction;
+  using Octree = internal::Octree_wrapper<GeomTraits>;
 
-  const Octree oct = octree;
-  const Topology topo = std::make_shared<Topology::element_type>(oct);
-  const Geometry geom = std::make_shared<Geometry::element_type>(oct);
-  const Point_function point_func = std::make_shared<Point_function::element_type>(point_function);
-  const Function func = std::make_shared<Function::element_type>(geom, point_func);
-  const Gradient grad = std::make_shared<Gradient::element_type>(gradient);
+  const Topology topo { octree };
+  const Geometry geom { octree };
+  const Point_function point_func { point_function };
+  const Function func { geom, point_func };
+  const Gradient grad { gradient };
 
-  return Domain(topo, geom, func, grad);
+  return { topo, geom, func, grad };
 }
 
 } // namespace Isosurfacing
diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Base_domain.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Base_domain.h
index b62b930469b..769f9bab58f 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Base_domain.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Base_domain.h
@@ -38,7 +38,7 @@ public:
   using Point = typename Geom_traits::Point_3;
   using Vector = typename Geom_traits::Vector_3;
 
-  using Topology = std::shared_ptr<Topology_>;
+  using Topology = Topology_;
   using Vertex_descriptor = typename Topology_::Vertex_descriptor;
   using Edge_descriptor = typename Topology_::Edge_descriptor;
   using Cell_descriptor = typename Topology_::Cell_descriptor;
@@ -48,9 +48,9 @@ public:
   using Cell_vertices = typename Topology_::Cell_vertices;
   using Cell_edges = typename Topology_::Cell_edges;
 
-  using Geometry = std::shared_ptr<Geometry_>;
-  using Function = std::shared_ptr<Function_>;
-  using Gradient = std::shared_ptr<Gradient_>;
+  using Geometry = Geometry_;
+  using Function = Function_;
+  using Gradient = Gradient_;
 
   static constexpr Cell_type CELL_TYPE = Topology_::CELL_TYPE;
   static constexpr std::size_t VERTICES_PER_CELL = Topology_::VERTICES_PER_CELL;
@@ -71,64 +71,64 @@ public:
   // gets the position of vertex `v`
   Point position(const Vertex_descriptor& v) const
   {
-    return geom->operator()(v);
+    return geom.operator()(v);
   }
 
   // gets the value of the function at vertex `v`
   FT value(const Vertex_descriptor& v) const
   {
-    return func->operator()(v);
+    return func.operator()(v);
   }
 
   // gets the gradient at vertex `v`
   Vector gradient(const Point& p) const
   {
-    return grad->operator()(p);
+    return grad.operator()(p);
   }
 
   // gets a container with the two vertices incident to the edge `e`
   Vertices_incident_to_edge edge_vertices(const Edge_descriptor& e) const
   {
-    return topo->edge_vertices(e);
+    return topo.edge_vertices(e);
   }
 
   // gets a container with all cells incident to the edge `e`
   Cells_incident_to_edge cells_incident_to_edge(const Edge_descriptor& e) const
   {
-    return topo->cells_incident_to_edge(e);
+    return topo.cells_incident_to_edge(e);
   }
 
   // gets a container with all vertices of the cell `c`
   Cell_vertices cell_vertices(const Cell_descriptor& c) const
   {
-    return topo->cell_vertices(c);
+    return topo.cell_vertices(c);
   }
 
   // gets a container with all edges of the cell `c`
   Cell_edges cell_edges(const Cell_descriptor& c) const
   {
-    return topo->cell_edges(c);
+    return topo.cell_edges(c);
   }
 
   // iterates over all vertices `v`, calling `f(v)` on each of them
   template <typename Concurrency_tag, typename Functor>
   void iterate_vertices(Functor& f) const
   {
-    topo->iterate_vertices(f, Concurrency_tag());
+    topo.iterate_vertices(f, Concurrency_tag());
   }
 
   // iterates over all edges `e`, calling `f(e)` on each of them
   template <typename Concurrency_tag, typename Functor>
   void iterate_edges(Functor& f) const
   {
-    topo->iterate_edges(f, Concurrency_tag());
+    topo.iterate_edges(f, Concurrency_tag());
   }
 
   // iterates over all cells `c`, calling `f(c)` on each of them
   template <typename Concurrency_tag, typename Functor>
   void iterate_cells(Functor& f) const
   {
-    topo->iterate_cells(f, Concurrency_tag());
+    topo.iterate_cells(f, Concurrency_tag());
   }
 
 private:
diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology.h
index 186b59a7b4e..e89fd0da810 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology.h
@@ -18,12 +18,12 @@
 #include <CGAL/Isosurfacing_3/internal/Tables.h>
 #include <CGAL/tags.h>
 
-#include <array>
-
 #ifdef CGAL_LINKED_WITH_TBB
 #include <tbb/parallel_for.h>
 #endif // CGAL_LINKED_WITH_TBB
 
+#include <array>
+
 namespace CGAL {
 namespace Isosurfacing {
 namespace internal {
diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h
index 80b0b15e016..e43cca72d7d 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Implicit_function_with_geometry.h
@@ -32,8 +32,8 @@ public:
   using Geom_traits = GeomTraits;
   using FT = typename Geom_traits::FT;
 
-  using Geometry = std::shared_ptr<Geometry_>;
-  using Point_function = std::shared_ptr<PointFunction>;
+  using Geometry = Geometry_;
+  using Point_function = PointFunction;
 
 public:
   // creates a function that uses the geometry to evaluate the function at vertex positions.
@@ -47,7 +47,7 @@ public:
   template <typename VertexDescriptor>
   FT operator()(const VertexDescriptor& v) const
   {
-    return func->operator()(geom->operator()(v));
+    return func.operator()(geom.operator()(v));
   }
 
 private:
diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_geometry.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_geometry.h
index 17f7545803a..6483fe751eb 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_geometry.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_geometry.h
@@ -30,7 +30,7 @@ public:
   using Geom_traits = GeomTraits;
   using Point = typename Geom_traits::Point_3;
 
-  using Octree = std::shared_ptr<Octree_wrapper<Geom_traits> >;
+  using Octree = Octree_wrapper<Geom_traits>;
 
   using Vertex_descriptor = typename Octree_topology<Geom_traits>::Vertex_descriptor;
 
@@ -41,11 +41,11 @@ public:
 
   Point operator()(const Vertex_descriptor& v) const
   {
-    return octree->point(v);
+    return octree.point(v);
   }
 
 private:
-  const Octree octree;
+  const Octree& octree;
 };
 
 } // namespace internal
diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_topology.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_topology.h
index e02743b9c4b..cf58318c2d5 100644
--- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_topology.h
+++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Octree_topology.h
@@ -32,11 +32,10 @@ class Octree_topology
 {
 public:
   using Geom_traits = GeomTraits;
-  using Octree_ = Octree_wrapper<Geom_traits>;
-  using Octree = std::shared_ptr<Octree_wrapper<Geom_traits> >;
-  using Vertex_descriptor = typename Octree_::Vertex_handle;
-  using Edge_descriptor = typename Octree_::Edge_handle;
-  using Cell_descriptor = typename Octree_::Voxel_handle;
+  using Octree = Octree_wrapper<Geom_traits>;
+  using Vertex_descriptor = typename Octree::Vertex_handle;
+  using Edge_descriptor = typename Octree::Edge_handle;
+  using Cell_descriptor = typename Octree::Voxel_handle;
 
   static constexpr Cell_type CELL_TYPE = CUBICAL_CELL;
   static constexpr std::size_t VERTICES_PER_CELL = 8;
@@ -54,42 +53,42 @@ public:
 
   Vertices_incident_to_edge edge_vertices(const Edge_descriptor& e) const
   {
-    return octree->edge_vertices(e);
+    return octree.edge_vertices(e);
   }
 
   Cells_incident_to_edge cells_incident_to_edge(const Edge_descriptor& e) const
   {
-    return octree->edge_voxels(e);
+    return octree.edge_voxels(e);
   }
 
   Cell_vertices cell_vertices(const Cell_descriptor& c) const
   {
-    return octree->voxel_vertices(c);
+    return octree.voxel_vertices(c);
   }
 
   Cell_edges cell_edges(const Cell_descriptor& c) const
   {
-    return octree->voxel_edges(c);
+    return octree.voxel_edges(c);
   }
 
   template <typename Functor>
   void iterate_vertices(Functor& f, Sequential_tag) const
   {
-    for(const Vertex_descriptor& v : octree->leaf_vertices())
+    for(const Vertex_descriptor& v : octree.leaf_vertices())
       f(v);
   }
 
   template <typename Functor>
   void iterate_edges(Functor& f, Sequential_tag) const
   {
-    for(const Edge_descriptor& e : octree->leaf_edges())
+    for(const Edge_descriptor& e : octree.leaf_edges())
       f(e);
   }
 
   template <typename Functor>
   void iterate_cells(Functor& f, Sequential_tag) const
   {
-    for(const Cell_descriptor& v : octree->leaf_voxels())
+    for(const Cell_descriptor& v : octree.leaf_voxels())
       f(v);
   }
 
@@ -97,7 +96,7 @@ public:
   template <typename Functor>
   void iterate_vertices(Functor& f, Parallel_tag) const
   {
-    const auto& vertices = octree->leaf_vertices();
+    const auto& vertices = octree.leaf_vertices();
 
     auto iterator = [&](const tbb::blocked_range<std::size_t>& r)
     {
@@ -111,7 +110,7 @@ public:
   template <typename Functor>
   void iterate_edges(Functor& f, Parallel_tag) const
   {
-    const auto& edges = octree->leaf_edges();
+    const auto& edges = octree.leaf_edges();
 
     auto iterator = [&](const tbb::blocked_range<std::size_t>& r)
     {
@@ -125,7 +124,7 @@ public:
   template <typename Functor>
   void iterate_cells(Functor& f, Parallel_tag) const
   {
-    const auto& cells = octree->leaf_voxels();
+    const auto& cells = octree.leaf_voxels();
 
     auto iterator = [&](const tbb::blocked_range<std::size_t>& r)
     {
@@ -138,7 +137,7 @@ public:
 #endif // CGAL_LINKED_WITH_TBB
 
 private:
-  const Octree octree;
+  const Octree& octree;
 };
 
 } // namespace internal
diff --git a/Isosurfacing_3/test/Isosurfacing_3/test_dual_contouring.cpp b/Isosurfacing_3/test/Isosurfacing_3/test_dual_contouring.cpp
index 0b11f074512..fab57b0d6e9 100644
--- a/Isosurfacing_3/test/Isosurfacing_3/test_dual_contouring.cpp
+++ b/Isosurfacing_3/test/Isosurfacing_3/test_dual_contouring.cpp
@@ -44,17 +44,17 @@ int main(int, char**)
   const std::size_t ny = static_cast<std::size_t>(2.0 / spacing.y());
   const std::size_t nz = static_cast<std::size_t>(2.0 / spacing.z());
 
-  std::shared_ptr<Grid> grid = std::make_shared<Grid>(nx, ny, nz, bbox);
+  Grid grid { nx, ny, nz, bbox };
 
-  for(std::size_t x=0; x<grid->xdim(); ++x) {
-    for(std::size_t y=0; y<grid->ydim(); ++y) {
-      for(std::size_t z=0; z<grid->zdim(); ++z)
+  for(std::size_t x=0; x<grid.xdim(); ++x) {
+    for(std::size_t y=0; y<grid.ydim(); ++y) {
+      for(std::size_t z=0; z<grid.zdim(); ++z)
       {
         const Point pos(x * spacing.x() + bbox.xmin(),
                         y * spacing.y() + bbox.ymin(),
                         z * spacing.z() + bbox.zmin());
 
-        grid->value(x, y, z) = sphere_function(pos);
+        grid.value(x, y, z) = sphere_function(pos);
       }
     }
   }
diff --git a/Isosurfacing_3/test/Isosurfacing_3/test_marching_cubes.cpp b/Isosurfacing_3/test/Isosurfacing_3/test_marching_cubes.cpp
index aade1870a41..7021accb490 100644
--- a/Isosurfacing_3/test/Isosurfacing_3/test_marching_cubes.cpp
+++ b/Isosurfacing_3/test/Isosurfacing_3/test_marching_cubes.cpp
@@ -62,15 +62,15 @@ void test_grid_sphere(const std::size_t n)
 
   std::shared_ptr<Grid> grid = std::make_shared<Grid>(n, n, n, bbox);
 
-  for(std::size_t x=0; x<grid->xdim(); ++x) {
-    for(std::size_t y=0; y<grid->ydim(); ++y) {
-      for(std::size_t z=0; z<grid->zdim(); ++z)
+  for(std::size_t x=0; x<grid.xdim(); ++x) {
+    for(std::size_t y=0; y<grid.ydim(); ++y) {
+      for(std::size_t z=0; z<grid.zdim(); ++z)
       {
         const Point pos(x * spacing.x() + bbox.xmin(),
                         y * spacing.y() + bbox.ymin(),
                         z * spacing.z() + bbox.zmin());
 
-        grid->value(x, y, z) = sphere_function(pos);
+        grid.value(x, y, z) = sphere_function(pos);
       }
     }
   }