compute ice area in quadrant

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

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

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

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

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