find_relaxation_boundary.m 3.86 KB
Newer Older
1
function Mobj = find_relaxation_boundary(Mobj)
2 3 4
% Find the elements which fall along the boundary for a nested
% configuration.
%
5
% Mobj = find_boundary_elements(Mobj)
6
%
7 8
% DESCRIPTION:
%   Find the elements which are bounded by the open boundaries described by
9 10 11
%   the nodes in Mobj.read_obc_nodes and which go to a depth of
%   Mobj.relax_bc_Nnodes.
%
12 13
% INPUT:
%   Mobj - required fields:
14 15 16 17 18 19
%           - read_obc_nodes - cell array of open boundary node IDs.
%           - tri - mesh triangulation
%           - relax_bc_Nnodes - array (length is Mobj.nObs) of number of
%           elements from the open boundary over which to relax the nested
%           inputs.
%
20
% OUTPUT:
21 22 23 24 25 26
%   Mobj - mesh object with the following new fields:
%           - relaxBC_nodes = node IDs for the relaxed region
%           - relaxBC_elems = element IDs for the relaxed region
%           - relaxnBC_nodes = number of nodes in the relaxed region
%           - relaxnBC_elems = number of elements in the relaxed region
%
27
% EXAMPLE USAGE:
28 29
%   Mobj = find_relaxation_boundary(Mobj)
%
30 31
% Author(s):
%   Pierre Cazenave (Plymouth Marine Laboratory)
32 33
%   Ricardo Torres (Plymouth Marine Laboratory)
%
34
% Revision history:
35 36 37 38 39 40 41 42 43 44 45
%   2016-12-15 Add support for varying depth of nested region over each
%   open boundary. Also update help to actually refer to what this function
%   does.
%
% Author(s):
%   Pierre Cazenave (Plymouth Marine Laboratory)
%   Ricardo Torres (Plymouth Marine Laboratory)
%
% Revision history:
%   2016-12-15 Update help to actually refer to what this function does.
%
46 47 48 49 50 51
%==========================================================================

subname = 'find_boundary_elements';

global ftbverbose
if ftbverbose
52
    fprintf('\nbegin : %s\n', subname)
53 54 55
end

nb = length(Mobj.read_obc_nodes); % number of boundaries
56
bc_width = Mobj.relax_bc_Nnodes;
57 58 59 60 61 62
obc_elems = cell(nb, 1);
nObcElements = nan(nb, 1);

for i = 1:nb
    
    % Do the current boundary's nodes
63 64 65 66 67 68 69 70
    nodeIDs = Mobj.read_obc_nodes{i};
    [C1,~] = ismember(Mobj.tri(:,1), nodeIDs(:), 'rows');
    [C2,~] = ismember(Mobj.tri(:,2), nodeIDs(:), 'rows');
    [C3,~] = ismember(Mobj.tri(:,3), nodeIDs(:), 'rows');
    obc_elems{i} = unique([find(C1); find(C2); find(C3)]);
    if iscolumn(obc_elems{i})
        obc_elems{i} = obc_elems{i}';
    end
71 72 73
    nObcElements(i) = numel(obc_elems{i}(:));
    
end
74 75
Mobj.relaxBC_nodes = {[Mobj.read_obc_nodes{:}]};
Mobj.relaxBC_elems = {[obc_elems{:}]};
76

77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
for bb = 2:bc_width
    nodeIDs = Mobj.tri(Mobj.relaxBC_elems{bb-1}, :);
    nodeIDs = unique(nodeIDs(:));
    % Find connected nodes.
    bc_nodes = Mobj.relaxBC_nodes{1:bb-1};
    if ~iscolumn(bc_nodes)
        bc_nodes = bc_nodes';
    end
    C1 = setdiff(nodeIDs(:), bc_nodes, 'rows');
    Mobj.relaxBC_nodes(bb) = {C1};
    % And now elements.
    [C1,~] = ismember(Mobj.tri(:,1), nodeIDs(:), 'rows');
    [C2,~] = ismember(Mobj.tri(:,2), nodeIDs(:), 'rows');
    [C3,~] = ismember(Mobj.tri(:,3), nodeIDs(:), 'rows');
    bc_elems = Mobj.relaxBC_elems{1:bb-1};
    if ~iscolumn(bc_elems)
        bc_elems = bc_elems';
    end
    C1 = setdiff(unique([find(C1); find(C2); find(C3)]), ...
        bc_elems, 'rows');
    Mobj.relaxBC_elems(bb) = {C1};
98 99
end

100 101 102 103 104
all_nest_elems = Mobj.relaxBC_elems{:};
all_nest_nodes = Mobj.relaxBC_nodes{:};
Mobj.relaxnBC_elems = length(all_nest_elems);
Mobj.relaxnBC_nodes = length(all_nest_nodes);

105 106 107 108
% Check it's worked for the first model boundary.
% xc = nodes2elems(Mobj.x, Mobj);
% yc = nodes2elems(Mobj.y, Mobj);
% figure
109
% clf
110
% triplot(Mobj.tri,Mobj.x,Mobj.y,'k');
111 112
% hold on
% 
113 114 115 116 117 118
% symbols = {'r.', 'k.', 'rx', 'kx', 'ro', 'ko'};
% for nn = 1:length(Mobj.relaxBC_nodes)
%     plot(Mobj.x(Mobj.relaxBC_nodes{nn}), ...
%         Mobj.y(Mobj.relaxBC_nodes{nn}), ...
%         symbols{mod(nn, length(Mobj.relaxBC_nodes)) + 1})
% end
119 120 121 122 123 124
% plot(xc(Mobj.relaxBC_elems{3}), yc(Mobj.relaxBC_elems{3}), 'kx')
% axis('equal', 'tight')

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