Commit 67b09248 authored by Pierre Cazenave's avatar Pierre Cazenave

Change the code which identifies nodes just inside the model open boundaries...

Change the code which identifies nodes just inside the model open boundaries to adjust their positions instead of creating a new array of 'ideal' positions. This means the tool now performs an optimisation on the unstructured grid rather than being a tool which simply provides a set of points which can be used in SMS to try and guide the triangulation (which I discovered does not work)
parent e3409693
function [x2, y2] = find_inside_boundary(x, y)
% Find points inside the given boundary.
%
% [x2, y2] = find_inside_boundary(x, y)
%
% DESCRIPTION:
% Find the coordinates of points which are normal to the open boundary
% described by the coordinates (x, y). The distance from the boundary is
% determined from the mean of the length of the two adjacent boundary
% element lengths.
%
% INPUT:
% x, y - coordinate pairs for the open boundary.
%
% OUTPUT:
% x2, y2 - coordinate pairs for the points as normal to the open boundary
% as possible (i.e. bisecting the angle between the two
% adjacent open boundary element faces).
%
% EXAMPLE USAGE:
% [x2, y2] = find_inside_boundary(x, y)
%
% NOTES:
% This works best with cartesian coordinates but will work with spherical
% too, although the angles for large elements will be incorrect (in an
% absolute sense).
%
% Author(s):
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history:
% 2013-03-11 First version.
%
%==========================================================================
subname = 'find_inside_boundary';
global ftbverbose
if ftbverbose
fprintf('\n'); fprintf(['begin : ' subname '\n']);
end
% Check the inputs
if length(x) ~= length(y)
error('Size of inputs (x, y) do not match.')
else
% Set the number of points in the boundary
np = length(x);
end
x2 = nan(np, 1);
y2 = nan(np, 1);
% Order the boundary points in clockwise order (for the polygon centroid
% calculation).
[x1, y1] = poly2cw(x, y);
[cx, cy] = centroid(x1(:), y1(:));
% For each node, find the two closest nodes to it and use those as the
% adjacent vectors. Doing this is slower than just assuming the nodes are
% provided in a sorted order, but means we don't have to worry too much
% about any irregularities from the sorting above (poly2cw).
for pp = 1:np
[~, idx] = sort(sqrt((x1(pp) - x1).^2 + (y1(pp) - y1).^2));
% Get the coordinates of the two nearest points (skip the first closest
% because that's the current point).
[px1, py1] = deal(x1(idx(2)), y1(idx(2)));
[px2, py2] = deal(x1(idx(3)), y1(idx(3)));
% Find the length of the edges of the triangle formed from the three
% points we're currently considering. 1 and 2 are those adjacent and 3
% is the remaining side.
ln1 = sqrt((x1(pp) - px1)^2 + (y1(pp) - py1)^2);
ln2 = sqrt((x1(pp) - px2)^2 + (y1(pp) - py2)^2);
ln3 = sqrt((px1 - px2)^2 + (py1 - py2)^2);
% Find the angle between the two element edges and the current node
% (cosine rule). Use the real component only for cases where the three
% points lie on a straight line (in which case the angle should be
% 180 degrees).
ang1 = real(acosd((ln1^2 + ln2^2 - ln3^2) / (2 * ln1 * ln2)));
ang1b = ang1 / 2; % bisect the angle
% Find the angle to the horizontal for the current node and one of the
% other points.
ang2 = atan2((py1 - y1(pp)), (px1 - x1(pp))) * (180 / pi);
% Find the difference between the two.
ang3 = ang2 - ang1b;
% Now get the mean length of the two closest element edges and use that
% to create the new point inside the boundary.
ml = mean([ln1, ln2]);
dy = ml * sind(ang3);
dx = ml * cosd(ang3);
% Add the offsets to the current node to get the new node's position.
% This is where things get a bit hairy: we'll assume that the boundary
% is approximately circular in nature. This means we can use its
% centroid as a tool to find the points inside the current boundary and
% those outside.
[xx(1), yy(1)] = deal(x1(pp) + dx, y1(pp) + dy);
[xx(2), yy(2)] = deal(x1(pp) - dx, y1(pp) - dy);
% Find the distances from the centroid.
dist = sqrt((cx - xx).^2 + (cy - yy).^2);
if dist(1) < dist(2)
[x2(pp, 1), y2(pp, 1)] = deal(xx(1), yy(1));
% [x2(pp, 2), y2(pp, 2)] = deal(xx(2), yy(2));
else
[x2(pp, 1), y2(pp, 1)] = deal(xx(2), yy(2));
% [x2(pp, 2), y2(pp, 2)] = deal(xx(1), yy(1));
end
end
if ftbverbose
fprintf('end : %s \n', subname)
end
function [adjx, adjy] = fix_inside_boundary(x, y, node_ids)
% Fix unstructured grid points inside the given boundary.
%
% [adjx, adjy] = fix_inside_boundary(Mobj, node_ids)
%
% DESCRIPTION:
% Find the coordinates of points which are normal to the open boundary
% described by the coordinates (x(node_ids), y(node_ids)). The distance
% from the boundary is determined from the mean of the length of the two
% adjacent boundary element lengths. Once the 'ideal' position has been
% identified, find the closest existing nodal position and change it to
% the 'ideal' position.
%
% The resulting x2 and y2 coordinates can be exported to 2dm to be
% checked in SMS for mesh quality with the fvcom-toolbox function
% write_SMS_2dm.
%
% INPUT:
% x, y - Unstructured grid coordinates.
% node_ids - List of IDs of the nodes within the grid which are on the
% open boundary of interest.
%
% OUTPUT:
% adjx, adjy - New unstructured grid coordinate pairs in which the points
% just inside the open boundary have been adjusted so as to
% bisect the angle formed by the two adjacent boundary
% faces.
%
% EXAMPLE USAGE:
% [adjx, adjy] = fix_inside_boundary(x, y, node_ids)
%
% NOTES:
% This works best with cartesian coordinates but will work with spherical
% too, although the angles for large elements will be incorrect.
%
% Author(s):
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history:
% 2013-03-11 First version.
%
%==========================================================================
subname = 'fix_inside_boundary';
global ftbverbose
if ftbverbose
fprintf('\n'); fprintf(['begin : ' subname '\n']);
end
% Check the inputs
if length(x) ~= length(y)
error('Size of inputs (x, y) do not match.')
else
% Set the number of points in the boundary
np = length(x(node_ids));
end
normx = nan(np, 1);
normy = nan(np, 1);
% Order the boundary points in clockwise order.
[boundx, boundy] = poly2cw(x(node_ids), y(node_ids));
% Create an array for the list of indices which have been moved which are
% adjacent to the open boundary.
seen = nan(np, 1);
% Create the output arrays from the input data. These values will be
% adjusted for the nodes adjacent to the open boundary defined by node_ids.
adjx = x;
adjy = y;
for pp = 1:np
% For each node, find the two closest nodes to it and use those as the
% adjacent vectors. Doing this is slower than just assuming the nodes
% are provided in a sorted order, but means we don't have to worry too
% much about any irregularities from the sorting above (poly2cw).
[~, idx] = sort(sqrt((boundx(pp) - boundx).^2 + (boundy(pp) - boundy).^2));
% Get the coordinates of the two nearest points (skip the first closest
% because that's the current point).
[px1, py1] = deal(boundx(idx(2)), boundy(idx(2)));
[px2, py2] = deal(boundx(idx(3)), boundy(idx(3)));
% Find the length of the edges of the triangle formed from the three
% points we're currently considering. 1 and 2 are those adjacent and 3
% is the remaining side.
ln1 = sqrt((boundx(pp) - px1)^2 + (boundy(pp) - py1)^2);
ln2 = sqrt((boundx(pp) - px2)^2 + (boundy(pp) - py2)^2);
ln3 = sqrt((px1 - px2)^2 + (py1 - py2)^2);
% Find the angle between the two element edges and the current node
% (cosine rule). Use the real component only for cases where the three
% points lie on a straight line (in which case the angle should be
% 180 degrees).
ang1 = real(acosd((ln1^2 + ln2^2 - ln3^2) / (2 * ln1 * ln2)));
ang1b = ang1 / 2; % bisect the angle
% Find the angle to the horizontal for the current node and one of the
% other points.
ang2 = atan2((py1 - boundy(pp)), (px1 - boundx(pp))) * (180 / pi);
% Find the difference between the two.
ang3 = ang2 - ang1b;
% Now get the mean length of the two closest element edges and use that
% to create the new point inside the boundary. Scale it to 90% of the
% value to make a cleaner transition when we're shifting points below.
ml = 0.9 * mean([ln1, ln2]);
dy = ml * sind(ang3);
dx = ml * cosd(ang3);
% Add the offsets to the current node to get the new node's position.
[xx(1), yy(1)] = deal(boundx(pp) + dx, boundy(pp) + dy);
[xx(2), yy(2)] = deal(boundx(pp) - dx, boundy(pp) - dy);
% Check which of the two sets above is inside the polygon defined by
% the open boundary.
if inpolygon(xx(1), yy(1), boundx, boundy) == 1
[normx(pp, 1), normy(pp, 1)] = deal(xx(1), yy(1));
elseif inpolygon(xx(2), yy(2), boundx, boundy) == 1
[normx(pp, 1), normy(pp, 1)] = deal(xx(2), yy(2));
else
warning('Both versions of the calculated point are outside the model domain. Skipping.')
continue
end
% OK, so now we have a new point approximately orthogonal to the
% current node, we can use its position to find the nearest existing
% node inside the domain and replace its coordinates with the new
% node's.
[~, idx2] = sort(sqrt((normx(pp, 1) - x).^2 + (normy(pp, 1) - y).^2));
% We need to check we haven't seen this node before (in 'seen') and
% also that it's not an open boundary point.
c = 1;
while true
if ismember(idx2(c), node_ids) || ismember(idx2(c), seen)
% Keep going until we find one that's not been used before.
c = c + 1;
else
break
end
end
% Append to the list of coordinates we've seen so we don't move the
% same node twice.
seen(pp) = idx2(c);
% Replace the coordinates.
adjx(idx2(c)) = normx(pp, 1);
adjy(idx2(c)) = normy(pp, 1);
end
if ftbverbose
fprintf('end : %s\n', subname)
end
% Do a figure to see the effect
% close all
% h = figure(1);
%
% % Original node positions
% plot(x, y, 'o', 'MarkerEdgeColor', [0.9, 0.1, 0.1], ...
% 'MarkerFaceColor', [1, 0.1, 0.1])
% hold on
% axis('equal', 'tight')
% % Adjusted node positions
% plot(adjx, adjy, 'o', 'MarkerEdgeColor', [0.1, 0.1, 0.8], ...
% 'MarkerFaceColor', [0.1, 0.1, 0.9])
%
% % Original triangulation
% patch('Vertices', [x, y], 'Faces', Mobj.tri, 'CData', [], ...
% 'edgecolor', [0.9, 0.1, 0.1], 'FaceColor', 'w', 'LineWidth', 1);
% % New triangulation
% patch('Vertices', [adjx, adjy], 'Faces', Mobj.tri, 'CData', [], ...
% 'edgecolor', [0.1, 0.1, 0.8], 'FaceColor', 'w', 'LineWidth', 1);
%
% % Open boundary nodes
% % plot(x(node_ids), y(node_ids), 'ko', 'MarkerFaceColor', 'k')
%
% legend('Original nodes', 'Adjusted nodes')
% legend('BoxOff')
%
% xlabel('Eastings')
% ylabel('Northings')
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment