diff --git a/examples/__init__.py b/examples/__init__.py
new file mode 100644
index 00000000..0a0428f7
--- /dev/null
+++ b/examples/__init__.py
@@ -0,0 +1 @@
+"""scikit-gmsh example for 3D mesh generation."""
diff --git a/examples/cylinder.py b/examples/cylinder.py
new file mode 100644
index 00000000..6ec77576
--- /dev/null
+++ b/examples/cylinder.py
@@ -0,0 +1,13 @@
+"""scikit-gmsh package for 3D mesh generation test."""
+
+from __future__ import annotations
+
+import pyvista as pv
+
+import skgmsh as sg
+
+edge_source = pv.Cylinder(resolution=16)
+edge_source.merge(pv.PolyData(edge_source.points), merge_points=True, inplace=True)
+edge_source.plot(show_edges=True)
+delaunay_3d = sg.Delaunay3D(edge_source)
+delaunay_3d.mesh.shrink(0.9).plot(show_edges=True)
diff --git a/src/skgmsh/__init__.py b/src/skgmsh/__init__.py
index c436f0eb..e8e01d89 100644
--- a/src/skgmsh/__init__.py
+++ b/src/skgmsh/__init__.py
@@ -16,6 +16,7 @@
 INITIAL_MESH_ONLY_2D = 3
 FRONTAL_DELAUNAY_2D = 6
 DELAUNAY_3D = 1
+INITIAL_MESH_ONLY_3D = 3
 
 SILENT = 0
 
@@ -125,12 +126,21 @@ def delaunay_3d(
 
     """
     points = edge_source.points
-    faces = edge_source.regular_faces
+    faces = edge_source.irregular_faces
 
     gmsh.initialize()
-    gmsh.option.set_number("Mesh.Algorithm3D", DELAUNAY_3D)
+    if target_sizes is None:
+        gmsh.option.set_number("Mesh.Algorithm", INITIAL_MESH_ONLY_3D)
+        gmsh.option.set_number("Mesh.MeshSizeExtendFromBoundary", 0)
+        gmsh.option.set_number("Mesh.MeshSizeFromPoints", 0)
+        gmsh.option.set_number("Mesh.MeshSizeFromCurvature", 0)
+    else:
+        gmsh.option.set_number("Mesh.Algorithm3D", DELAUNAY_3D)
     gmsh.option.set_number("General.Verbosity", SILENT)
 
+    if target_sizes is None:
+        target_sizes = 0.0
+
     if isinstance(target_sizes, float):
         target_sizes = [target_sizes] * edge_source.number_of_points
 
@@ -140,11 +150,13 @@ def delaunay_3d(
 
     surface_loop = []
     for i, face in enumerate(faces):
-        gmsh.model.geo.add_line(face[0] + 1, face[1] + 1, i * 4 + 0)
-        gmsh.model.geo.add_line(face[1] + 1, face[2] + 1, i * 4 + 1)
-        gmsh.model.geo.add_line(face[2] + 1, face[3] + 1, i * 4 + 2)
-        gmsh.model.geo.add_line(face[3] + 1, face[0] + 1, i * 4 + 3)
-        gmsh.model.geo.add_curve_loop([i * 4 + 0, i * 4 + 1, i * 4 + 2, i * 4 + 3], i + 1)
+        curve_tags = []
+        for j, _ in enumerate(face):
+            start_tag = face[j - 1] + 1
+            end_tag = face[j] + 1
+            curve_tag = gmsh.model.geo.add_line(start_tag, end_tag)
+            curve_tags.append(curve_tag)
+        gmsh.model.geo.add_curve_loop(curve_tags, i + 1)
         gmsh.model.geo.add_plane_surface([i + 1], i + 1)
         gmsh.model.geo.remove_all_duplicates()
         gmsh.model.geo.synchronize()
@@ -291,3 +303,42 @@ def edge_source(self: Delaunay2D) -> pv.PolyData:
     def mesh(self: Delaunay2D) -> pv.UnstructuredGrid:
         """Get the mesh."""
         return self._mesh
+
+
+class Delaunay3D:
+    """
+    Delaunay 3D mesh algorithm.
+
+    Parameters
+    ----------
+    edge_source : pyvista.PolyData
+        Specify the source object used to specify constrained
+        edges and loops. If set, and lines/polygons are defined, a
+        constrained triangulation is created. The lines/polygons
+        are assumed to reference points in the input point set
+        (i.e. point ids are identical in the input and
+        source).
+
+    Notes
+    -----
+    .. versionadded:: 0.2.0
+
+    """
+
+    def __init__(
+        self: Delaunay3D,
+        edge_source: pv.PolyData,
+    ) -> None:
+        """Initialize the Delaunay3D class."""
+        self._edge_source = edge_source
+        self._mesh = delaunay_3d(edge_source)
+
+    @property
+    def edge_source(self: Delaunay3D) -> pv.PolyData:
+        """Get the edge source."""
+        return self._edge_source
+
+    @property
+    def mesh(self: Delaunay3D) -> pv.UnstructuredGrid:
+        """Get the mesh."""
+        return self._mesh
diff --git a/tests/test_skgmsh.py b/tests/test_skgmsh.py
index 3d8b1e67..b463baa3 100644
--- a/tests/test_skgmsh.py
+++ b/tests/test_skgmsh.py
@@ -65,3 +65,12 @@ def test_delaunay_3d(target_sizes: float | Sequence[float]) -> None:
     assert np.allclose(mesh.volume, edge_source.volume)
     # TODO @tkoyama010: Compare cell type. # noqa: FIX002
     # https://github.com/pyvista/scikit-gmsh/pull/125
+
+
+@pytest.mark.parametrize("edge_source", [pv.Cube(), pv.Cylinder()])
+def test_delaunay_3d_default(edge_source: pv.PolyData) -> None:
+    """Delaunay 3D mesh algorithm test code."""
+    edge_source.merge(pv.PolyData(edge_source.points), merge_points=True, inplace=True)
+    delaunay_3d = sg.Delaunay3D(edge_source)
+    mesh = delaunay_3d.mesh
+    assert np.allclose(mesh.volume, edge_source.volume)