diff --git a/get_elements2d.m b/get_elements2d.m
new file mode 100644
index 0000000..c6233fa
--- /dev/null
+++ b/get_elements2d.m
@@ -0,0 +1,9 @@
+function [e] = get_elements2d(mesh)
+    % function GET_ELEMENTS2D - 2D elements of a 2D or 3D mesh.
+    
+    if isa(mesh, 'mesh2d')
+        e = mesh.elements;
+    else
+        e = mesh.elements2d;
+    end
+end
diff --git a/get_x2d.m b/get_x2d.m
new file mode 100644
index 0000000..5d363a7
--- /dev/null
+++ b/get_x2d.m
@@ -0,0 +1,8 @@
+function [e] = get_x2d(mesh)
+    % function GET_X2D - 2D x coordinates of a 2D or 3D mesh.
+    if isa(mesh, 'mesh2d')
+        e = mesh.x;
+    else
+        e = mesh.x2d;
+    end
+end
diff --git a/get_y2d.m b/get_y2d.m
new file mode 100644
index 0000000..63c92c3
--- /dev/null
+++ b/get_y2d.m
@@ -0,0 +1,9 @@
+function [e] = get_y2d(mesh)
+    % function GET_Y2D - 2D y coordinates of a 2d or 3d mesh
+    
+    if isa(mesh, 'mesh2d')
+        e = mesh.y;
+    else
+        e = mesh.y2d;
+    end
+end
diff --git a/levelset_area.m b/levelset_area.m
new file mode 100644
index 0000000..2e69d27
--- /dev/null
+++ b/levelset_area.m
@@ -0,0 +1,73 @@
+function [sum_area] = levelset_area(md, levelset, poly)
+    %function LEVELSET_AREA - area of the region where the levelset is negative
+    %
+    % Sum of the area of elements of the mesh that are inside the region where the levelset is negative.
+    % Partial area is computed for elements that are partially in the region by linear interpolation of the levelset field.
+    % If a polygon is passed as the optional third argument, only the area inside that polygon is computed.
+    % Partial area is computed for elements that are partially within the polygon by intersecting the element with the polygon. 
+    % Assumes a mesh of (extruded) triangles.
+    % 
+    % md: ISSM model
+    % levelset: levelset field that defines the area
+    % poly: optional, polygon ([x1 y1; x2, y2; ...]) that defines the polygon
+    %
+    % e.g. ice_area = levelset_area(md, md.mask.ice_levelset)
+    %      ice_area_northeast = levelset_area(md, md.mask.ice_levelset, [0 0; 1e6 0; 1e6 1e6; 0, 1e6])
+
+    if nargin == 2
+        % no polygon given, use polygon that encloses the whole model
+        [minx, maxx] = bounds(md.mesh.x);
+        [miny, maxy] = bounds(md.mesh.y);
+        lx = maxx - minx;
+        ly = maxy - miny;
+        minx = minx - lx;
+        maxx = maxx + lx;
+        miny = miny - ly;
+        maxy = maxy + ly;
+        poly = [minx, miny; maxx, miny; maxx, maxy; minx, maxy];
+    end
+
+    [front_vertices, front_element_indices] = levelset_interface(md, levelset);
+    el2d = get_elements2d(md.mesh);
+    sum_area = 0;
+
+    for i = 1:size(el2d, 1)
+        % area of each element
+        vertex_is_inside = levelset(el2d(i, :)) <= 0;
+        if sum(vertex_is_inside) == 0
+            continue; % no vertices in the region of negative levelset, ignore element
+        end
+
+        element_front_vertex_idx = find(front_element_indices == i);
+
+        % find the polygon that describes the part of the element that is inside the negative levelset 
+        if isempty(element_front_vertex_idx)
+            % no front vertices -> whole element inside the front
+            p = [md.mesh.x(el2d(i, :)), md.mesh.y(el2d(i, :))];
+        else
+            % one or two points inside the front ... 
+            p = [md.mesh.x(el2d(i, find(vertex_is_inside))), md.mesh.y(el2d(i, find(vertex_is_inside)))];
+            % ... plus the front segment form a triangle or quadrilateral.
+            p = [p; flipud(squeeze(front_vertices(element_front_vertex_idx, :, :)))];
+            % for quadrilateral, the order of the front vertices needs to be flipped, otherwise the polygon self-intersects.
+            % (depends on the order of vertices delivered by the `levelset_interface` function)
+            % for triangle, order doesn't matter.
+        end
+
+        %intersect with polygon
+        vertex_is_inside = inpolygon(p(:, 1), p(:, 2), poly(:, 1), poly(:, 2));
+        if ~any(vertex_is_inside)
+            % all vertices outside the polygon, ignore element
+            continue;
+        elseif all(vertex_is_inside)
+            % all vertices inside the polygon
+            % cheap path, area of a polygon is quick
+            sum_area = sum_area + polyarea(p(:, 1), p(:, 2)); 
+        else
+            % some vertices inside the polygon, use polygon intersection
+            % expensive path, matlab functions `polyshape` and `intersect` are slow, but hopefully only a few elements.
+            % probably could be replaced with a faster implementations.
+            sum_area = sum_area + area(intersect(polyshape(p), polyshape(poly)));
+        end
+    end
+end
diff --git a/levelset_interface.m b/levelset_interface.m
new file mode 100644
index 0000000..0ce25de
--- /dev/null
+++ b/levelset_interface.m
@@ -0,0 +1,54 @@
+function [contour, element_indices] = levelset_interface(md, levelset)
+    % Function LEVELSET_INTERFACE - zero isocontour of the levelset
+    %
+    % The contour where the levelset field is zero as an unordered list of vertex pairs.
+    % Each vertex pair is the points where the contour intersects an element of the mesh, forming a line segment of the contour.
+    % Intersections are computed by linear interpolation of the levelset field along the edges of an element.
+    % Returns an N x 2 x 2 matrix where N is the number of elements that are intersected by the contour.
+    % contour(i, 1, 1) and contour(i, 1, 2) are the x and y coordinates of the first vertex of the i-th contour segment.
+    % contour(i, 2, 1) and contour(i, 2, 2) are the x and y coordinates of the second vertex of the i-th contour segment.
+    %
+    % md: ISSM model, 2d or 3d, (extruded) triangle mesh.
+    % levelset: levelset field to find the zero isocontour for.
+    %
+
+    % which elements are intersected by the zero isocontour
+    mesh = md.mesh;
+    elements = get_elements2d(mesh);
+    n_in = sum(levelset(elements) <= 0, 2);
+    is_interface = n_in > 0 & n_in < 3;
+    interface_elements = elements(is_interface, :);
+    
+    % compute points where each element is intersected
+    num_elements = size(interface_elements, 1);
+    contour = zeros(num_elements, 2, 2);
+    if nargout > 1
+        element_indices = find(is_interface);
+    end
+    for i = 1:num_elements
+        current_element = interface_elements(i, :);
+
+        %find the edges that are intersected by the interface
+        vert_in = current_element(levelset(current_element) <= 0);
+        vert_out = current_element(levelset(current_element) > 0);
+        [I, O] = meshgrid(vert_in, vert_out);
+        interface_edges = [I(:), O(:)];
+        
+        %pt 1
+        edge1 = interface_edges(1, :);
+        v1 = abs(levelset(edge1(1)));
+        v2 = abs(levelset(edge1(2)));
+        w = v1 / (v1 + v2);
+        pt1 = w * [mesh.x(edge1(2)), mesh.y(edge1(2))] + (1-w) * [mesh.x(edge1(1)), mesh.y(edge1(1))];
+
+        %pt 1
+        edge2 = interface_edges(2, :);
+        v1 = abs(levelset(edge2(1)));
+        v2 = abs(levelset(edge2(2)));
+        w = v1 / (v1 + v2);
+        pt2 = w * [mesh.x(edge2(2)), mesh.y(edge2(2))] + (1-w) * [mesh.x(edge2(1)), mesh.y(edge2(1))];
+        
+        contour(i, 1, :) = pt1;
+        contour(i, 2, :) = pt2;
+    end
+end