fix_inside_boundary.m 7.29 KB
Newer Older
1
function [adjx, adjy] = fix_inside_boundary(x, y, node_ids, ang_thresh)
2
% Fix unstructured grid points inside the given boundary.
3 4
%
% [adjx, adjy] = fix_inside_boundary(x, y, node_ids, ang_thresh)
5 6 7 8 9 10 11 12
% 
% 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.
13
%
14 15 16
%   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.
17
%
18 19 20 21
% 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.
22 23 24
%   ang_thresh - [optional] Specify a minimum angle in degrees between the
%                two adjacent nodal positions to deal with corners better.
%
25 26 27 28 29
% 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.
30
%
31
% EXAMPLE USAGE:
32 33
%   [adjx, adjy] = fix_inside_boundary(Mobj.x, Mobj.y, Mobj.read_obc_nodes{1}, 90)
%
34 35 36
% NOTES:
%   This works best with cartesian coordinates but will work with spherical
%   too, although the angles for large elements will be incorrect.
37 38 39 40
%   Secondly, this will sometimes place put points outside the model domain
%   (though I'm not yet sure why). The net result is that you have to
%   re-edit the grid SMS by deleting that particular node and recreating
%   the triangulation manually.
41
%
42 43
% Author(s):
%   Pierre Cazenave (Plymouth Marine Laboratory)
44
%
45 46
% Revision history:
%   2013-03-11 First version.
47 48
%   2013-03-19 Add optional minimum angle support.
%
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
%==========================================================================

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

66 67 68 69 70 71
if nargin == 4
    minAng = true;
else
    minAng = false;
end

72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
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
112 113 114 115 116 117 118 119 120

    % Check if we've been given a threshold minimum angle and skip this
    % open boundary point if we have.
    if minAng
        if ang1 < ang_thresh
            continue
        end
    end

121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
    % 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')