A technique to improve quality of elements on the boundary #305
Description
- Using the functions below, one can discretize a polygon so the spacing of its nodes follows a sizing function.
- Along each boundary segment of the polygon, you can place a new node at a distance from the segment's midpoint so a equilateral triangle is formed.
- During mesh generation, these points associated with the boundary are fixed so they cannot move. In other words you can pass these points as fixed. Note that you should remove any existing initialized points in the vicinity of these fixed points. Second note: in areas where the boundary folds in on itself more work is needed to avoid degenerate triangles.
The result is dramatically improved mesh quality along the boundary
function [pout, t, converged] = mesh1d(poly, edgeFuncHandles, minEdgeLength, fixedPoints, boundaryBox, boxNum, varargin)
% Mesh Generator using Distance Functions.
%
% Inputs:
% poly - Nx2 array of polygon vertices.
% edgeFuncHandles - Cell array of functions for edge length determination.
% minEdgeLength - Minimum edge length.
% fixedPoints - Indices of fixed polygon vertices.
% boundaryBox - Boundary box for mesh.
% boxNum - Index of the boundary box for initial edge length.
%
% Outputs:
% pout - Nx2 array of resampled node positions.
% t - Edge indices (NTx2).
% converged - Boolean flag indicating convergence.
% Initialize variables
X = poly(:,1);
Y = poly(:,2);
xd = diff(X);
yd = diff(Y);
distances = sqrt(xd.^2 + yd.^2);
cumulativeDist = [0; cumsum(distances)];
totalDist = cumulativeDist(end);
% Calculate number of points and parameterize
numPoints = totalDist / min(minEdgeLength);
t = linspace(0, max(cumulativeDist), ceil(numPoints));
% Interpolate new points
newX = interp1(cumulativeDist, X, t);
newY = interp1(cumulativeDist, Y, t);
% Prepare for mesh generation
p = [newX', newY']; % Initial point distribution
fixedIndices = unique([1, length(p)]); % Ensure the start and end points are fixed
p(fixedIndices, :) = poly(fixedPoints, :); % Apply fixed points
% Distance function for the domain
fdist = @(x, y) sqrt(min((x - boundaryBox(:,1)).^2 + (y - boundaryBox(:,2)).^2));
% Edge length function (simplified for demonstration)
fh = @(x, y) minEdgeLength;
% Main mesh generation loop
[p, t, converged] = generateMesh(p, fdist, fh, fixedIndices, totalDist, minEdgeLength);
% Map the final points back to the original coordinates
pout = [interp1(cumulativeDist, X, p(:,1)), interp1(cumulativeDist, Y, p(:,2))];
if any(isnan(pout))
warning('NaNs detected in 1D mesh');
end
end
function [p, t, converged] = generateMesh(p, fdist, fh, fixedIndices, totalDist, minEdgeLength)
% Initialize variables
converged = false;
maxIterations = 100;
iteration = 0;
ptol = 1e-3; % Tolerance for point movement
% Start the iteration
while ~converged && iteration < maxIterations
iteration = iteration + 1;
% Calculate edge lengths and desired adjustments
edges = diff(p, 1, 1); % Differences between points
edgeLengths = sqrt(sum(edges.^2, 2));
desiredLengths = fh(p(:,1), p(:,2)); % Assuming fh returns a scalar or vector the same length as p
% Calculate movement for each point based on edge length discrepancy
moveDist = (desiredLengths - edgeLengths) .* (edges ./ edgeLengths);
moveDist = [0, 0; moveDist]; % No movement for the first point
% Apply movements, excluding fixed points
isFixed = false(size(p, 1), 1);
isFixed(fixedIndices) = true;
p(~isFixed, :) = p(~isFixed, :) + moveDist(~isFixed, :);
% Check for convergence
maxMove = max(sqrt(sum(moveDist.^2, 2)));
if maxMove < ptol
converged = true;
end
end
% Generate edge indices (for visualization or further processing)
t = [(1:size(p, 1)-1)', (2:size(p, 1))'];
end
function [initial_points, layer_of_points] = form_initial_points_near_and_along_boundary(outer_polygon, ...
bounding_polygon, ...
minimum_resolution_in_meters, ...
mesh_sizing_fx)
% Form a layer of points one row along the boundary rim of the domain
% so that when triangulated it forms perfectly equilateral triangles.
minimum_resolution_in_degrees = minimum_resolution_in_meters / 111e3;
% ensure the first point is not repeated etc.
outer_polygon = unique(outer_polygon,'rows','stable');
%% Resample the outer polygon to match the point spacing following mesh_sizing_fx
[outer_polygon_rs,~,~]=mesh1d(outer_polygon, ...
mesh_sizing_fx, ...
minimum_resolution_in_degrees, ...
[], ...
bounding_polygon, ...
1, ... % box_number (assumed 1).
[] ...
);
layer_of_points = [];
%figure;
for i = 1 : length(outer_polygon_rs)-1
% 1) Compute the midpoint
p1 = outer_polygon_rs(i,:);
p2 = outer_polygon_rs(i+1,:);
midpoint = (p1 + p2) / 2.0;
% 2) Compute the tangent and then the unit normal vector pointing inwards
tangent = p2 - p1; % Tangent vector
tangent_length = norm(tangent); % Length of the tangent vector
unit_tangent = tangent / tangent_length; % Normalize to get unit tangent
unit_normal = [-unit_tangent(2), unit_tangent(1)]; % Rotate by 90 degrees to get the unit normal
% 3) Calculate the distance for the equilateral triangle and compute p3
local_size_in_degrees = mesh_sizing_fx{1}(midpoint) / 111e3;
% The distance from the midpoint to the third point to form an equilateral triangle
% Assuming local_size_in_degrees gives the desired edge length of the equilateral triangle
height_of_triangle = (sqrt(3) / 2) * local_size_in_degrees; % Height of an equilateral triangle
p3 = midpoint + height_of_triangle * unit_normal; % New point forming equilateral triangle
% For debugging plot
% hold on;
% plot([p1(1), p2(1)], [p1(2), p2(2)], 'k-'); % Plot edge
% plot(midpoint(1), midpoint(2), 'k.'); % Plot midpoint
% quiver(midpoint(1), midpoint(2), unit_normal(1), unit_normal(2), 0.05, 'r'); % Plot unit normal
% plot(p3(1), p3(2), 'bs'); % Plot new point
layer_of_points = [layer_of_points; p3];
end
% Combine initial resampled outer_polygon with layer_of_points
initial_points = [outer_polygon_rs; layer_of_points];
end