diff --git a/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h b/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h
index 270f4cdceb9..09737aca984 100644
--- a/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h
+++ b/Isosurfacing_3/doc/Isosurfacing_3/Concepts/IsosurfacingDomain.h
@@ -30,6 +30,11 @@ public:
     */
     typedef unspecified_type Point;
 
+    /*!
+    The vector type.
+    */
+    typedef unspecified_type Vector;
+
     /*!
     A handle to identify a vertex.
     */
@@ -78,10 +83,15 @@ public:
     Point position(const Vertex_handle& v) const;
 
     /*!
-    Returns the stored value of vertex v
+    Returns the value of vertex v
     */
     FT value(const Vertex_handle& v) const;
 
+    /*!
+    (Optional) Returns the gradient at the position p
+    */
+    Vector gradient(const Point& p) const;
+
     /*!
     Returns the two vertices incident to edge e
     */
diff --git a/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h b/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h
index e2184d61d86..a01b366f285 100644
--- a/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h
+++ b/Isosurfacing_3/include/CGAL/Cartesian_grid_domain.h
@@ -1,8 +1,20 @@
+// Copyright (c) 2022 INRIA Sophia-Antipolis (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Julian Stahl
+
 #ifndef CGAL_CARTESIAN_GRID_DOMAIN_H
 #define CGAL_CARTESIAN_GRID_DOMAIN_H
 
 #include <CGAL/Cartesian_grid_3.h>
 #include <CGAL/Cartesian_topology_base.h>
+#include <CGAL/Default_gradients.h>
 #include <CGAL/Isosurfacing_3/internal/Tables.h>
 
 #ifdef CGAL_LINKED_WITH_TBB
@@ -12,7 +24,7 @@
 namespace CGAL {
 namespace Isosurfacing {
 
-template <class GeomTraits>
+template <class GeomTraits, typename Gradient = Zero_gradient<GeomTraits>>
 class Cartesian_grid_domain : public Cartesian_topology_base {
 public:
     typedef GeomTraits Geom_traits;
@@ -20,9 +32,10 @@ public:
     typedef typename Geom_traits::Point_3 Point;
     typedef typename Geom_traits::Vector_3 Vector;
     typedef typename Geom_traits::Vector_3 Grid_spacing;
+    typedef Cartesian_grid_3<Geom_traits> Grid;
 
 public:
-    Cartesian_grid_domain(const Cartesian_grid_3<Geom_traits>& grid) : grid(&grid) {}
+    Cartesian_grid_domain(const Grid& grid, const Gradient& grad = Gradient()) : grid(&grid), grad(&grad) {}
 
     Point position(const Vertex_handle& v) const {
         const Bbox_3& bbox = grid->get_bbox();
@@ -32,8 +45,8 @@ public:
                      v[2] * spacing.z() + bbox.zmin());
     }
 
-    Vector gradient(const Vertex_handle& v) const {
-        return grid->gradient(v[0], v[1], v[2]);
+    Vector gradient(const Point& p) const {
+        return grad->operator()(p);
     }
 
     FT value(const Vertex_handle& v) const {
@@ -103,7 +116,7 @@ public:
         const std::size_t size_z = grid->zdim();
 
         //#pragma omp parallel for
-        //for (int x = 0; x < size_x - 1; x++) {
+        // for (int x = 0; x < size_x - 1; x++) {
         //    for (std::size_t y = 0; y < size_y - 1; y++) {
         //        for (std::size_t z = 0; z < size_z - 1; z++) {
         //            f({(std::size_t)x, y, z});
@@ -120,13 +133,15 @@ public:
                 }
             }
         };
-        
+
         tbb::parallel_for(tbb::blocked_range<std::size_t>(0, size_x - 1), iterator);
     }
 #endif  // CGAL_LINKED_WITH_TBB
 
 private:
-    const Cartesian_grid_3<Geom_traits>* grid;
+    const Grid* grid;
+
+    const Gradient* grad;
 };
 
 }  // namespace Isosurfacing
diff --git a/Isosurfacing_3/include/CGAL/Default_gradients.h b/Isosurfacing_3/include/CGAL/Default_gradients.h
index a159e385dcc..f241c780a07 100644
--- a/Isosurfacing_3/include/CGAL/Default_gradients.h
+++ b/Isosurfacing_3/include/CGAL/Default_gradients.h
@@ -1,3 +1,14 @@
+// Copyright (c) 2022 INRIA Sophia-Antipolis (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Julian Stahl
+
 #ifndef CGAL_DEFAULT_GRADIENT_H
 #define CGAL_DEFAULT_GRADIENT_H
 
@@ -28,7 +39,7 @@ public:
 public:
     Finite_difference_gradient(const Function& func, const FT delta = 0.001) : func(&func), delta(delta) {}
 
-    Vector operator()(const Point& point) const {
+    Vector operator()(const Point& point) const {  // TODO
         const Point p0 = point + Vector(delta, 0, 0);
         const Point p1 = point - Vector(delta, 0, 0);
         const Point p2 = point + Vector(0, delta, 0);
diff --git a/Isosurfacing_3/include/CGAL/Implicit_domain.h b/Isosurfacing_3/include/CGAL/Implicit_domain.h
index ba481812f17..35c0834bb2d 100644
--- a/Isosurfacing_3/include/CGAL/Implicit_domain.h
+++ b/Isosurfacing_3/include/CGAL/Implicit_domain.h
@@ -1,3 +1,14 @@
+// Copyright (c) 2022 INRIA Sophia-Antipolis (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Julian Stahl
+
 #ifndef CGAL_IMPLICIT_DOMAIN_H
 #define CGAL_IMPLICIT_DOMAIN_H
 
@@ -22,13 +33,13 @@ public:
     typedef typename Geom_traits::Vector_3 Grid_spacing;
 
 public:
-    Implicit_domain(const Bbox_3& domain, const Grid_spacing& spacing, const Function& func,
+    Implicit_domain(const Bbox_3& bbox, const Grid_spacing& spacing, const Function& func,
                     const Gradient& grad = Gradient())
-        : bbox(domain), spacing(spacing), func(&func), grad(&grad) {
+        : bbox(bbox), spacing(spacing), func(&func), grad(&grad) {
 
-        sizes[0] = domain.x_span() / spacing.x() + 1;
-        sizes[1] = domain.y_span() / spacing.y() + 1;
-        sizes[2] = domain.z_span() / spacing.z() + 1;
+        sizes[0] = bbox.x_span() / spacing.x() + 1;
+        sizes[1] = bbox.y_span() / spacing.y() + 1;
+        sizes[2] = bbox.z_span() / spacing.z() + 1;
     }
 
     Point position(const Vertex_handle& v) const {
@@ -36,8 +47,8 @@ public:
                      v[2] * spacing.z() + bbox.zmin());
     }
 
-    Vector gradient(const Vertex_handle& v) const {
-        return grad->operator()(position(v));
+    Vector gradient(const Point& p) const {
+        return grad->operator()(p);
     }
 
     FT value(const Vertex_handle& v) const {
diff --git a/Isosurfacing_3/include/CGAL/Octree_domain.h b/Isosurfacing_3/include/CGAL/Octree_domain.h
index 955b87fe647..9640ac5ce10 100644
--- a/Isosurfacing_3/include/CGAL/Octree_domain.h
+++ b/Isosurfacing_3/include/CGAL/Octree_domain.h
@@ -1,7 +1,19 @@
+// Copyright (c) 2022 INRIA Sophia-Antipolis (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s)     : Julian Stahl
+
 #ifndef CGAL_OCTREE_DOMAIN_H
 #define CGAL_OCTREE_DOMAIN_H
 
 #include <CGAL/Cell_type.h>
+#include <CGAL/Default_gradients.h>
 #include <CGAL/Octree_wrapper.h>
 
 #ifdef CGAL_LINKED_WITH_TBB
@@ -13,7 +25,7 @@
 namespace CGAL {
 namespace Isosurfacing {
 
-template <typename GeomTraits>
+template <typename GeomTraits, typename Gradient = Zero_gradient<GeomTraits>>
 class Octree_domain {
 public:
     typedef GeomTraits Geom_traits;
@@ -36,14 +48,14 @@ public:
     typedef std::array<Edge_handle, 12> Cell_edges;
 
 public:
-    Octree_domain(const Octree& octree) : octree_(&octree) {}
+    Octree_domain(const Octree& octree, const Gradient& grad = Gradient()) : octree_(&octree), grad(&grad) {}
 
     Point position(const Vertex_handle& v) const {
         return octree_->point(v);
     }
 
-    Vector gradient(const Vertex_handle& v) const {
-        return octree_->gradient(v);
+    Vector gradient(const Point& p) const {
+        return grad->operator()(p);
     }
 
     FT value(const Vertex_handle& v) const {
@@ -130,6 +142,8 @@ public:
 
 private:
     const Octree* octree_;
+
+    const Gradient* grad;
 };
 
 }  // namespace Isosurfacing