function mesh_2d = flatten_mesh(mesh_3d, tol)


bb = boundingbox(mesh_3d);

if nargin < 2
    tol = (bb(2,3) - bb(1,3)) / 1e6;
end

[mesh_2d, elemind] = get_boundary(mesh_3d);
[mesh_2d, ~, keepind] = ...
    mesh_section(mesh_2d, [-inf -inf bb(1,3)-tol; inf inf bb(1,3)+tol]);

mesh_2d = drop_unused_nodes(mesh_2d);
mesh_2d = drop_mesh_IDs(mesh_2d);

mesh_2d.Nodes(:,4) = 0;

dataind = elemind(keepind);
mesh_2d.Elements(dataind,:) = mesh_2d.Elements;

mesh_2d.Nodes(:, 4) = 0;

end