diff --git a/Project.toml b/Project.toml
index dc15685..bee2040 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "3.0.2"
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
@@ -14,6 +15,7 @@ StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 ColorSchemes = "3"
 Colors = "0.12, 0.13"
 DocStringExtensions = "0.8, 0.9"
+GeometryBasics = "0.5.9"
 StaticArrays = "1.9.13"
 StaticArraysCore = "1.4"
 julia = "1.6"
diff --git a/src/GridVisualizeTools.jl b/src/GridVisualizeTools.jl
index a2c2563..a3fe542 100644
--- a/src/GridVisualizeTools.jl
+++ b/src/GridVisualizeTools.jl
@@ -11,6 +11,7 @@ import ColorSchemes
 using DocStringExtensions: SIGNATURES, TYPEDEF, TYPEDSIGNATURES
 using StaticArraysCore: SVector
 using StaticArrays: @MArray
+using GeometryBasics # replace by SArray!
 
 include("colors.jl")
 export region_cmap, bregion_cmap, rgbtuple, rgbcolor
@@ -18,9 +19,16 @@ export region_cmap, bregion_cmap, rgbtuple, rgbcolor
 include("extraction.jl")
 export extract_visible_cells3D, extract_visible_bfaces3D
 
+include("linearsimplex.jl")
+export LinearSimplex, LinearSimplexIterator
+export LinearEdge, LinearTriangle, LinearTetrahedron
+export testloop
+
 include("marching.jl")
+include("marching_iterators.jl")
 export marching_tetrahedra, marching_triangles
 
+
 include("markerpoints.jl")
 export markerpoints
 
diff --git a/src/colors.jl b/src/colors.jl
index 4461c7b..a05a6d3 100644
--- a/src/colors.jl
+++ b/src/colors.jl
@@ -50,6 +50,8 @@ function bregion_cmap(n)
 end
 
 """
+$(SIGNATURES)
+
 Create color tuple from  color description (e.g. string)
 
 ```jldoctest
diff --git a/src/experimental/gridsimplex.jl b/src/experimental/gridsimplex.jl
new file mode 100644
index 0000000..f8d086a
--- /dev/null
+++ b/src/experimental/gridsimplex.jl
@@ -0,0 +1,58 @@
+"""
+    GridSimplex
+Grid simplex with point coordinates, visible faces and regions.
+"""
+struct GridSimplex{D, N, Tv, Ti}
+    points::SVector{N, SVector{D, Tv}}
+    visibleFaces::SVector{N, Bool}
+    region::Ti
+end
+
+points(s::GridSimplex) = s.points
+visibleFaces(s::GridSimplex) = s.visibleFaces
+region(s::GridSimplex) = s.region
+
+function GridSimplex(::Type{Val{D}}, points::Union{Vector, Tuple}, visibleFaces::Union{Vector, Tuple}, region::Ti) where {D}
+    simplexPoints = SVector{D + 1, SVector{D, Float32}}(points...)
+    simplexVisibleFaces = SVector{D + 1, Bool}(visibleFaces...)
+    return GridSimplex(simplexPoints, simplexVisibleFaces, region)
+end
+
+function GridSimplex(::Type{Val{1}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti)
+    @views simplexPoints = SVector{2, SVector{1, Float32}}(points[:, 1], points[:, 2])
+    simplexVisibleFaces = SVector{2, Bool}(visibleFaces[1], visibleFaces[2])
+    return GridSimplex(simplexPoints, simplexVisibleFaces, region)
+end
+
+function GridSimplex(::Type{Val{2}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti)
+    @views simplexPoints = SVector{3, SVector{2, Float32}}(points[:, 1], points[:, 2], points[:, 3])
+    simplexVisibleFaces = SVector{3, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3])
+    return GridSimplex(simplexPoints, simplexVisibleFaces, region)
+end
+
+function GridSimplex(::Type{Val{3}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti)
+    @views simplexPoints = SVector{4, SVector{3, Float32}}(points[:, 1], points[:, 2], points[:, 3], points[:, 4])
+    simplexVisibleFaces = SVector{4, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3], visibleFaces[4])
+    return GridSimplex(simplexPoints, simplexVisibleFaces, region)
+end
+
+function GridSimplex(::Type{Val{1}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti, coordscale)
+    @views simplexPoints = SVector{2, SVector{1, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale)
+    simplexVisibleFaces = SVector{2, Bool}(visibleFaces[1], visibleFaces[2])
+    return GridSimplex(simplexPoints, simplexVisibleFaces, region)
+end
+
+function GridSimplex(::Type{Val{2}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti, coordscale)
+    @views simplexPoints = SVector{3, SVector{2, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale)
+    simplexVisibleFaces = SVector{3, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3])
+    return GridSimplex(simplexPoints, simplexVisibleFaces, region)
+end
+
+function GridSimplex(::Type{Val{3}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti, coordscale)
+    @views simplexPoints = SVector{4, SVector{3, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale, points[:, 4] * coordscale)
+    simplexVisibleFaces = SVector{4, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3], visibleFaces[4])
+    return GridSimplex(simplexPoints, simplexVisibleFaces, region)
+end
+
+GridEdge(points, visibleFaces, region) = GridSimplex(Val{1}, points, visibleFaces, region)
+GridTriangle(points, visibleFaces, region) = GridSimplex(Val)
diff --git a/src/experimental/linearsimplex.jl b/src/experimental/linearsimplex.jl
new file mode 100644
index 0000000..29f6167
--- /dev/null
+++ b/src/experimental/linearsimplex.jl
@@ -0,0 +1,99 @@
+"""
+    LinearSimplex
+Linear simplex with point coordinates and node values.
+"""
+struct LinearSimplex{D, N, Tv}
+    points::SVector{N, SVector{D, Tv}}
+    values::SVector{N, Tv}
+end
+
+values(s::LinearSimplex) = s.values
+points(s::LinearSimplex) = s.points
+
+
+function LinearSimplex(::Type{Val{D}}, points::Union{Vector, Tuple}, values::Union{Vector, Tuple}) where {D}
+    spoints = SVector{D + 1, SVector{D, Float32}}(points...)
+    svalues = SVector{D + 1, Float32}(values...)
+    return LinearSimplex(spoints, svalues)
+end
+
+function LinearSimplex(::Type{Val{1}}, points::AbstractMatrix, values::AbstractVector)
+    @views spoints = SVector{2, SVector{1, Float32}}(points[:, 1], points[:, 2])
+    svalues = SVector{2, Float32}(values[1], values[2])
+    return LinearSimplex(spoints, svalues)
+end
+
+function LinearSimplex(::Type{Val{2}}, points::AbstractMatrix, values::AbstractVector)
+    @views spoints = SVector{3, SVector{2, Float32}}(points[:, 1], points[:, 2], points[:, 3])
+    svalues = SVector{3, Float32}(values[1], values[2], values[3])
+    return LinearSimplex(spoints, svalues)
+end
+
+function LinearSimplex(::Type{Val{3}}, points::AbstractMatrix, values::AbstractVector)
+    @views spoints = SVector{4, SVector{3, Float32}}(points[:, 1], points[:, 2], points[:, 3], points[:, 4])
+    svalues = SVector{4, Float32}(values[1], values[2], values[3], values[4])
+    return LinearSimplex(spoints, svalues)
+end
+
+
+function LinearSimplex(::Type{Val{1}}, points::AbstractMatrix, values::AbstractVector, coordscale)
+    @views spoints = SVector{2, SVector{1, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale)
+    svalues = SVector{2, Float32}(values[1], values[2])
+    return LinearSimplex(spoints, svalues)
+end
+
+function LinearSimplex(::Type{Val{2}}, points::AbstractMatrix, values::AbstractVector, coordscale)
+    @views spoints = SVector{3, SVector{2, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale)
+    svalues = SVector{3, Float32}(values[1], values[2], values[3])
+    return LinearSimplex(spoints, svalues)
+end
+
+function LinearSimplex(::Type{Val{3}}, points::AbstractMatrix, values::AbstractVector, coordscale)
+    @views spoints = SVector{4, SVector{3, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale, points[:, 4] * coordscale)
+    svalues = SVector{4, Float32}(values[1], values[2], values[3], values[4])
+    return LinearSimplex(spoints, svalues)
+end
+
+LinearEdge(points, values) = LinearSimplex(Val{1}, points, values)
+LinearTriangle(points, values) = LinearSimplex(Val{2}, points, values)
+LinearTetrahedron(points, values) = LinearSimplex(Val{3}, points, values)
+
+
+"""
+    abstract type LinearSimplexIterator{D}
+
+Iterator over D-dimensional linear simplices. 
+
+Any subtype `TSub` should comply with the [iteration interface](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration)
+and implement
+```julia
+     Base.iterate(vs::TSub{D}, state):: (LinearSimplex{D}, state) where D
+```
+The optional methods (in particular, size, length) are not needed.
+"""
+abstract type LinearSimplexIterator{D} end
+
+"""
+    testloop(iterators::AbstractVector{T}) where T<:LinearSimplexIterator
+
+Sum up all values passed via the iterator. Designed for testing iterators.
+
+Useful for:
+- Consistency test - e.g. when passing constant ones, the sum can be calculated by 
+  other means and compared to the return value of testloop
+- Allocation test. `testloop` itself does allocate only a few (<1000) bytes in
+  very few (<10) allocations. So any allocations beyond this from a
+  call to `testloop` hint at possibilities to improve an iterator implementation.
+"""
+function testloop(iterators::AbstractVector{T}) where {T <: LinearSimplexIterator}
+    threads = map(iterators) do iterator
+        begin
+            local x = 0.0
+            for vt in iterator
+                x += sum(vt.values)
+            end
+            x
+        end
+    end
+    return sum(fetch.(threads))
+end
diff --git a/src/marching.jl b/src/marching.jl
index 85899ab..7acc106 100644
--- a/src/marching.jl
+++ b/src/marching.jl
@@ -78,7 +78,6 @@ function calculate_plane_tetrahedron_intersection!(
         )
         return 0
     end
-
     amount_intersections = 0
 
     @inbounds for n1 in 1:4
@@ -233,7 +232,7 @@ function marching_tetrahedra(
 
     function pushtris(ns, ixcoord, ixvalues)
         # number of intersection points can be 3 or 4
-        if ns >= 3
+        return if ns >= 3
             last_i = length(all_ixvalues)
             for is in 1:ns
                 @views push!(all_ixcoord, ixcoord[:, is])
@@ -356,7 +355,6 @@ function marching_triangles(
     points = Vector{Tp}(undef, 0)
     values = Vector{Tv}(undef, 0)
     adjacencies = Vector{SVector{2, Ti}}(undef, 0)
-
     for igrid in 1:length(coords)
         func = funcs[igrid]
         coord = coords[igrid]
@@ -416,7 +414,6 @@ function marching_triangles(
                 # connect last two points
                 push!(adjacencies, SVector{2, Ti}((length(points) - 1, length(points))))
             end
-
             return
         end
 
diff --git a/src/marching_iterators.jl b/src/marching_iterators.jl
new file mode 100644
index 0000000..b5f68e4
--- /dev/null
+++ b/src/marching_iterators.jl
@@ -0,0 +1,64 @@
+function pushintersection!(intersection_points, triangle::LinearSimplex{2}, levels)
+    f = values(triangle)
+    coord = points(triangle)
+
+    (n1, n2, n3) = (1, 2, 3)
+
+    f[1] <= f[2] ? (n1, n2) = (1, 2) : (n1, n2) = (2, 1)
+    f[n2] <= f[3] ? n3 = 3 : (n2, n3) = (3, n2)
+    f[n1] > f[n2] ? (n1, n2) = (n2, n1) : nothing
+
+    dx31 = coord[n3][1] - coord[n1][1]
+    dx21 = coord[n2][1] - coord[n1][1]
+    dx32 = coord[n3][1] - coord[n2][1]
+
+    dy31 = coord[n3][2] - coord[n1][2]
+    dy21 = coord[n2][2] - coord[n1][2]
+    dy32 = coord[n3][2] - coord[n2][2]
+
+    df31 = f[n3] != f[n1] ? 1 / (f[n3] - f[n1]) : 0.0
+    df21 = f[n2] != f[n1] ? 1 / (f[n2] - f[n1]) : 0.0
+    df32 = f[n3] != f[n2] ? 1 / (f[n3] - f[n2]) : 0.0
+
+    for level in levels
+        if (f[n1] <= level) && (level < f[n3])
+            α = (level - f[n1]) * df31
+            x1 = coord[n1][1] + α * dx31
+            y1 = coord[n1][2] + α * dy31
+
+            if (level < f[n2])
+                α = (level - f[n1]) * df21
+                x2 = coord[n1][1] + α * dx21
+                y2 = coord[n1][2] + α * dy21
+            else
+                α = (level - f[n2]) * df32
+                x2 = coord[n2][1] + α * dx32
+                y2 = coord[n2][2] + α * dy32
+            end
+            push!(intersection_points, Point{2, Float32}(x1, y1))
+            push!(intersection_points, Point{2, Float32}(x2, y2))
+        end
+    end
+    return
+end
+
+function intersections(triangles::T, levels) where {T <: LinearSimplexIterator{2}}
+    local intersection_points = Vector{Point{2, Float32}}(undef, 0)
+    for triangle in triangles
+        pushintersection!(intersection_points, triangle, levels)
+    end
+    return intersection_points
+end
+
+function marching_triangles(triangle_iterators::Vector{T}, levels) where {T <: LinearSimplexIterator{2}}
+    return if Threads.nthreads() == 1
+        map(triangle_iterators) do triangles
+            intersections(triangles, levels)
+        end
+    else
+        threads = map(triangle_iterators) do triangles
+            Threads.@spawn intersections(triangles, levels)
+        end
+        fetch.(threads)
+    end
+end