diff --git a/src/visualization/types.jl b/src/visualization/types.jl
index be7bc78a075..fa7e77fc5bf 100644
--- a/src/visualization/types.jl
+++ b/src/visualization/types.jl
@@ -102,6 +102,19 @@ struct PlotData2DTriangulated{DataType, NodeType, FaceNodeType, FaceDataType, Va
   face_data::FaceDataType
   variable_names::VariableNames
 end
+        
+struct PlotData3DTriangulated{DataType, NodeType, FaceNodeType, FaceDataType, VariableNames, PlottingTriangulation}
+  x::NodeType # physical nodal coordinates, size (num_plotting_nodes x num_elements)
+  y::NodeType
+  z::NodeType
+  data::DataType
+  t::PlottingTriangulation
+  x_face::FaceNodeType
+  y_face::FaceNodeType
+  z_face::FaceNodeType
+  face_data::FaceDataType
+  variable_names::VariableNames
+end
 
 # Show only a truncated output for convenience (the full data does not make sense)
 function Base.show(io::IO, pd::PlotData2DTriangulated)
diff --git a/src/visualization/utilities.jl b/src/visualization/utilities.jl
index 9eeb25cdaa7..3d28b1c722d 100644
--- a/src/visualization/utilities.jl
+++ b/src/visualization/utilities.jl
@@ -1276,6 +1276,45 @@ function global_plotting_triangulation_makie(pds::PlotDataSeries{<:PlotData2DTri
   plotting_mesh = merge([trimesh...]) # merge meshes on each element into one large mesh
   return plotting_mesh
 end
+        
+# 3D isosurface case for global_plotting_triangulation_makie function
+function global_plotting_triangulation_makie(plot_data, level)
+
+  xp = plot_data.x
+  yp = plot_data.y
+  zp = plot_data.z
+  func = plot_data.data
+  connectivity = plot_data.t
+
+  plotting_coordinates = zeros(3, size(xp, 1))
+  num_elements = size(xp, 2)
+  list_of_meshes = Vector{GeometryBasics.Mesh{3, Float32}}(undef, num_elements)
+  sk = 1
+  planes = []
+
+  for e in 1:size(func, 2) # for each column
+    plotting_coordinates[1, :] .= xp[:, e]
+    plotting_coordinates[2, :] .= yp[:, e]
+    plotting_coordinates[3, :] .= zp[:, e]
+
+    pts, trngls, fvals = GridVisualize.marching_tetrahedra(plotting_coordinates,
+                                                           connectivity,
+                                                           func[:, e], planes, level)
+
+    #output mesh that encompasses isosurface
+    if length(pts) > 0
+      makie_triangles = Makie.to_triangles(hcat((getindex.(trngls,i) for i in 1:3)...))
+      iso_mesh = GeometryBasics.normal_mesh(Makie.to_vertices(pts), makie_triangles)
+
+      # add newly found mesh to list of meshes
+      list_of_meshes[sk] = iso_mesh
+      sk += 1
+    end    
+  end
+
+  plotting_mesh = merge(list_of_meshes[1:sk-1])
+  return plotting_mesh
+end
 
 # Returns a list of `Makie.Point`s which can be used to plot the mesh, or a solution "wireframe"
 # (e.g., a plot of the mesh lines but with the z-coordinate equal to the value of the solution).