diff --git a/CHANGELOG.md b/CHANGELOG.md
index a3afbfe..4fe035d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,14 @@
 
 All notable changes to this project will be documented in this file.
 
+## [3.0.0] - 2025-03-03
+
+### Changed
+
+- `marching_triangles` takes now an additional `lines` argument to compute intersections with given lines and
+   returns also values and adjacency information to be able to construct a 1D grid
+   To restore the old behavior, use `lines = []` and neglect the second and third return value
+
 ## [1.1.1] - 2024-11-02
 - Update links after move to WIAS-PDELib org
 
diff --git a/Project.toml b/Project.toml
index 507a1ba..0645e39 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "GridVisualizeTools"
 uuid = "5573ae12-3b76-41d9-b48c-81d0b6e61cc5"
-authors = ["Jürgen Fuhrmann <juergen-fuhrmann@web.de>"]
-version = "2.0.0"
+authors = ["Jürgen Fuhrmann <juergen-fuhrmann@web.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de"]
+version = "2.1.0"
 
 [deps]
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
diff --git a/src/marching.jl b/src/marching.jl
index 49fba57..5686e94 100644
--- a/src/marching.jl
+++ b/src/marching.jl
@@ -259,43 +259,81 @@ end
 """
     $(SIGNATURES)
 
-Collect isoline snippets on triangles ready for linesegments!
+March through the given grid and extract points and values for given iso-line levels and/or given intersection lines.
+From the returned point list and value list a line plot can be created.
+
+Input:
+    coord: matrix storing the coordinates of the grid
+    cellnodes: connectivity matrix
+    func: function on the grid nodes to be evaluated
+    lines: vector of line definitions [a,b,c], s.t., ax + by + c = 0 defines a line
+    levels: vector of levels for the iso-surface
+    Tc: scalar type of coordinates
+    Tp: vector type of coordinates
+    Tv: scalar type of function values
+
+Output:
+    points: vector of 2D points of the intersections of the grid with the iso-surfaces or lines
+    adjacencies: vector of 2D vectors storing connected points in the grid
+    value: interpolated values of `func` at the intersection points
+
+Note that passing both nonempty `lines` and `levels` will create a result with both types of points mixed.
 """
 function marching_triangles(
-        coord::Matrix{Tv},
+        coord::Matrix{T},
         cellnodes::Matrix{Ti},
         func,
+        lines,
         levels;
-        Tc = Float32,
-        Tp = SVector{2, Tc}
-    ) where {Tv <: Number, Ti <: Number}
-    return marching_triangles([coord], [cellnodes], [func], levels; Tc, Tp)
+        Tc = T,
+        Tp = SVector{2, Tc},
+        Tv = Float64
+    ) where {T <: Number, Ti <: Number}
+    return marching_triangles([coord], [cellnodes], [func], lines, levels; Tc, Tp, Tv)
 end
 
+
+"""
+    $(SIGNATURES)
+
+
+Variant of `marching_triangles` with multiple grid input
+"""
 function marching_triangles(
-        coords::Vector{Matrix{Tv}},
+        coords::Vector{Matrix{T}},
         cellnodes::Vector{Matrix{Ti}},
         funcs,
+        lines,
         levels;
-        Tc = Float32,
-        Tp = SVector{2, Tc}
-    ) where {Tv <: Number, Ti <: Number}
+        Tc = T,
+        Tp = SVector{2, Tc},
+        Tv = Float64
+    ) where {T <: Number, Ti <: Number}
     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]
 
-        function isect(nodes)
+        # pre-allcate memory for triangle values (3 nodes per triangle)
+        objective_values = Vector{Tv}(undef, 3)
+
+        # the objective_func is used to determine the intersection (line equation or iso levels)
+        # the value_func is used to interpolate values at the intersections
+        function isect(tri_nodes, objective_func, value_func)
             (i1, i2, i3) = (1, 2, 3)
 
-            f = (func[nodes[1]], func[nodes[2]], func[nodes[3]])
+            # 3 values of the objective function
+            f = objective_func
 
+            # sort f[i1] ≤ f[i2] ≤ f[i3]
             f[1] <= f[2] ? (i1, i2) = (1, 2) : (i1, i2) = (2, 1)
             f[i2] <= f[3] ? i3 = 3 : (i2, i3) = (3, i2)
             f[i1] > f[i2] ? (i1, i2) = (i2, i1) : nothing
 
-            (n1, n2, n3) = (nodes[i1], nodes[i2], nodes[i3])
+            (n1, n2, n3) = (tri_nodes[i1], tri_nodes[i2], tri_nodes[i3])
 
             dx31 = coord[1, n3] - coord[1, n1]
             dx21 = coord[1, n2] - coord[1, n1]
@@ -305,35 +343,67 @@ function marching_triangles(
             dy21 = coord[2, n2] - coord[2, n1]
             dy32 = coord[2, n3] - coord[2, n2]
 
-            df31 = f[i3] != f[i1] ? 1 / (f[i3] - f[i1]) : 0.0
-            df21 = f[i2] != f[i1] ? 1 / (f[i2] - f[i1]) : 0.0
-            df32 = f[i3] != f[i2] ? 1 / (f[i3] - f[i2]) : 0.0
-
-            for level in levels
-                if (f[i1] <= level) && (level < f[i3])
-                    α = (level - f[i1]) * df31
-                    x1 = coord[1, n1] + α * dx31
-                    y1 = coord[2, n1] + α * dy31
-
-                    if (level < f[i2])
-                        α = (level - f[i1]) * df21
-                        x2 = coord[1, n1] + α * dx21
-                        y2 = coord[2, n1] + α * dy21
-                    else
-                        α = (level - f[i2]) * df32
-                        x2 = coord[1, n2] + α * dx32
-                        y2 = coord[2, n2] + α * dy32
-                    end
-                    push!(points, SVector{2, Tc}((x1, y1)))
-                    push!(points, SVector{2, Tc}((x2, y2)))
+            df31 = f[i3] != f[i1] ? 1 / (f[i1] - f[i3]) : 0.0
+            df21 = f[i2] != f[i1] ? 1 / (f[i1] - f[i2]) : 0.0
+            df32 = f[i3] != f[i2] ? 1 / (f[i2] - f[i3]) : 0.0
+
+            if (f[i1] <= 0) && (0 < f[i3])
+                α = f[i1] * df31
+                x1 = coord[1, n1] + α * dx31
+                y1 = coord[2, n1] + α * dy31
+                value1 = value_func[n1] + α * (value_func[n3] - value_func[n1])
+
+                if (0 < f[i2])
+                    α = f[i1] * df21
+                    x2 = coord[1, n1] + α * dx21
+                    y2 = coord[2, n1] + α * dy21
+                    value2 = value_func[n1] + α * (value_func[n2] - value_func[n1])
+                else
+                    α = f[i2] * df32
+                    x2 = coord[1, n2] + α * dx32
+                    y2 = coord[2, n2] + α * dy32
+                    value2 = value_func[n2] + α * (value_func[n3] - value_func[n2])
                 end
+
+                push!(points, SVector{2, Tc}((x1, y1)))
+                push!(points, SVector{2, Tc}((x2, y2)))
+                push!(values, value1)
+                push!(values, value2)
+                # connect last two points
+                push!(adjacencies, SVector{2, Ti}((length(points) - 1, length(points))))
             end
+
             return
         end
 
         for itri in 1:size(cellnodes[igrid], 2)
-            @views isect(cellnodes[igrid][:, itri])
+
+            # nodes of the current triangle
+            tri_nodes = @views cellnodes[igrid][:, itri]
+
+            for level in levels
+                # objective func is iso-level equation
+                @views @fastmath map!(
+                    inode -> (func[inode] - level),
+                    objective_values,
+                    tri_nodes
+                )
+                @views isect(tri_nodes, objective_values, func)
+            end
+
+            for line in lines
+                @fastmath line_equation(line, coord) = coord[1] * line[1] + coord[2] * line[2] + line[3]
+
+                # objective func is iso-level equation
+                @views @fastmath map!(
+                    inode -> (line_equation(line, coord[:, inode])),
+                    objective_values,
+                    tri_nodes
+                )
+                @views isect(tri_nodes, objective_values, func)
+            end
+
         end
     end
-    return points
+    return points, adjacencies, values
 end