find_boundary_elements.m 2.54 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
function Mobj = find_boundary_elements(Mobj)
% Find the elements which fall along the boundary.
% 
% Mobj = find_boundary_elements(Mobj)
% 
% DESCRIPTION:
%   Find the elements which are bounded by the open boundaries described by
%   the nodes in Mobj.read_obc_nodes.
% 
% INPUT:
%   Mobj - required fields:
%           - read_obc_nodes
%           - obc_nodes
%           - tri
% 
% OUTPUT:
%   Mobj - new field of a cell array read_obc_elements which contains the
18 19 20
%          IDs of the elements which fall on the model open boundaries and
%          nObcElements which is the total number of boundary elements
%          along each boundary.
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
% 
% NOTES:
%   This will be pretty slow if your unstructured grid has an enormous
%   number of elements in it (it loops through every element and compares
%   against the boundary nodes). I'm sure there's a quicker way, so feel
%   free to have at it.
% 
% EXAMPLE USAGE:
%   Mobj = find_boundary_elements(Mobj)
% 
% Author(s):
%   Pierre Cazenave (Plymouth Marine Laboratory)
% 
% Revision history:
%   2013-02-26 First version.
36 37
%   2013-02-28 Add new field to the output (total number of boundary
%   elements as nObcElements).
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
% 
%==========================================================================

subname = 'find_boundary_elements';

global ftbverbose
if ftbverbose
    fprintf('\n')
    fprintf('begin : %s\n', subname)
end

ne = length(Mobj.tri); % number of elements
nb = length(Mobj.read_obc_nodes); % number of boundaries

obc_elems = cell(nb, 1);
53
nObcElements = nan(nb, 1);
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74

for i = 1:nb

    % Do the current boundary's nodes
    nodeIDs = Mobj.obc_nodes(i, Mobj.obc_nodes(i, :) ~= 0);

    f = 0;
    
    for ttt = 1:ne
        tri = Mobj.tri(ttt, :);
        C = intersect(tri, nodeIDs);
        % Only those with a face along the boundary count (i.e. two nodes
        % on the boundary), particularly for the mean flow.
        if numel(C) == 2
            f = f + 1; % increment the found counter
            obc_elems{i}(f) = ttt;
            if ftbverbose
                fprintf('Found boundary element ID %i\n', ttt)
            end
        end
    end
75
    nObcElements(i) = numel(obc_elems{i}(:));
76 77 78
end

Mobj.read_obc_elements = obc_elems;
79
Mobj.nObcElements = nObcElements;
80 81 82 83 84 85 86 87 88 89 90 91 92 93

% Check it's worked for the first model boundary.
% xc = nodes2elems(Mobj.x, Mobj);
% yc = nodes2elems(Mobj.y, Mobj);
% figure(1)
% clf
% plot(Mobj.x, Mobj.y, 'r.', xc, yc, 'ko')
% hold on
% plot(xc(Mobj.read_obc_elements{1}), yc(Mobj.read_obc_elements{1}), 'gx')
% axis('equal', 'tight')

if ftbverbose
    fprintf('end   : %s \n', subname)
end