find_nesting_region.m 12.2 KB
Newer Older
1
function Nested = find_nesting_region(conf, Mobj)
Ricardo Torres's avatar
Ricardo Torres committed
2 3
% Creates a nesting structure array for direct/indirect or weighted nesting
%
4
% function Nested = find_nesting_region(conf, Mobj)
Ricardo Torres's avatar
Ricardo Torres committed
5 6
%
% DESCRIPTION:
7
%   Uses the Mesh object Mobj and the conf variable to search for the nodes
Ricardo Torres's avatar
Ricardo Torres committed
8 9 10 11 12 13 14 15 16 17 18 19
%   and elements originating from the open boundaries into de domain
%   interior for a specified number of levels.
%
% Optionally specify nesting type:
%   1/2: DIRECT/INDIRECT nesting:
%       - Full variables/no surface elevation respectively.
%   3:   RELAXATION nesting:
%       - Nesting with a relaxation method.
%
% INPUT:
%  conf     = struct whose field names are the model configuration options
%              as set by the user. These are generally assigned at the top of a
20 21 22 23 24 25 26 27 28 29
%              make_input.m type of script. Minimum fields are:
%              - Nested_type: Array of nesting types, one for each open boundary.
%              1 or 2 are direct nesting, 3 is weighted. If we're doing
%              type 3, we can specify the number of levels of nested
%              boundaries to use. The minimum valid value is 1. For
%              Indirect or direct nesting use 1
%              - levels: number of levels of nodes over which to relax the
%              boundary (if Nested_type is 3).
%              - power = determines drop of weights from 1. 0 is linear,
%              1-4 is 1/conf.levels.^conf.power.
Pierre Cazenave's avatar
Pierre Cazenave committed
30
%  Mobj     = Mesh object with the following fields:
31 32 33 34 35 36 37 38 39 40 41 42 43 44
%              - tri: Triangulation table as pupulated by read_sms_grid
%              - x: Node x coordinates (cartesian)
%              - y: Node y coordinates (cartesian)
%              - xc: element x coordinates (cartesian)
%              - yc: Element y coordinates (cartesian)
%              - nObcNodes: number of nodes as defined by SMS and read
%              through read_sms_grid
%              - read_obc_nodes = nodes indices at boundaries as read in
%              read_sms_grid.
%              - obc_type = Type of OBC as defined in mod_obcs.F [ values
%              of 1-10] and parsed to Mobj from conf by add_obc_nodes_list
%              - obc_nodes: matrix with node indices. Each row is a given
%              open boundary.
%              - nObs: number of open boundaries.
Ricardo Torres's avatar
Ricardo Torres committed
45 46
%
% OUTPUT:
Pierre Cazenave's avatar
Pierre Cazenave committed
47 48 49 50 51 52 53 54 55 56 57 58
%  Nested   = Mesh object with all the same fields as Mobj, plus the
%             following added and modified fields:
%              - read_obc_nodes: new nested boundary node IDs
%              - read_obc_elems: new nested boundary element IDs
%              - nObs: number of open boundaries (number of levels * number
%              - nObcNodes: number of nodes in each nested level
%              of original open boundaries)
%              - obc_type: the type for each nested boundary
%              - obc_nodes: the array-based nested boundary node IDs (for
%              backwards compatibility)
%              - weight_node: weights for each nested node boundary level
%              - weight_cell: weights for each nested element boundary level
Ricardo Torres's avatar
Ricardo Torres committed
59 60
%
% EXAMPLE USAGE:
Pierre Cazenave's avatar
Pierre Cazenave committed
61
%   conf.Nested_type = type of nesting [1, 2 == direct nesting, 3 == weighted]
Ricardo Torres's avatar
Ricardo Torres committed
62
%   conf.levels = number of boundary bands to use for weighted option
Pierre Cazenave's avatar
Pierre Cazenave committed
63 64
%   conf.power = determines drop of weights from 1 [0 is linear, anything
%   else is 1/conf.levels.^conf.power]
Ricardo Torres's avatar
Ricardo Torres committed
65 66 67 68 69 70 71 72 73 74 75
%   Mobj.tri = Triangulation table as pupulated by read_sms_grid
%   Mobj.x = Nodes x coordinate [we should make it possible to use lon instead]
%   Mobj.y = Nodes y coordinate [we should make it possible to use lat instead]
%   Mobj.nObcNodes = number of nodes as defined by SMS and read through read_sms_grid
%   Mobj.read_obc_nodes = nodes indices at boundaries as read in read_sms_grid
%   Mobj.obc_type = Type of OBC as defined in mod_obcs.F [ values of 1-10]
%       and parsed to Mobj from conf by add_obc_nodes_list
%   Mobj.xc = Nodes x coordinate [we should make it possible to use lonc instead]
%   Mobj.yc = Elements y coordinate [we should make it possible to use latc instead]
%   Mobj.obc_nodes = matrix with node indices. Each row is a boundary  level.
%   Mobj.nObs = total number of open boundary levels (I think this is set
Pierre Cazenave's avatar
Pierre Cazenave committed
76
%       in setup_metrics.m).
Ricardo Torres's avatar
Ricardo Torres committed
77 78 79 80
%
%   the global variable ftbverbose shows information at run time as well as
%   generating a figure with the location of the nesting nodes and elements
%
Pierre Cazenave's avatar
Pierre Cazenave committed
81
%   Nested = find_nesting_region(conf,Mobj)
Ricardo Torres's avatar
Ricardo Torres committed
82 83 84 85
%
%   Nested.nObcNodes = number of nodes in new nesting zone
%   Nested.read_obc_nodes = nodes indices in nesting zone
%   Nested.obc_type = Type of OBC as defined in mod_obcs.F [ values of 1-10]
86
%   Nested.obc_nodes = matrix with node indices. Each row is a boundary level.
Ricardo Torres's avatar
Ricardo Torres committed
87 88 89 90 91 92 93 94 95 96 97 98 99
%   Nested.nObs = total number of open boundary levels
%   Nested.weight_node = weights for nodes if using weighted type [0-1] see
%       FVCOM manual for further info (Chapter 6 section 4 in version 3.2)
%   Nested.weight_cell = weights for elements if using weighted type [0-1]

% Author(s):
%   Ricardo Torres (Plymouth Marine Laboratory)
%   Pierre Cazenave (Plymouth Marine Laboratory)
%   Darren Price (CH2MHill)
%   Hakeem Johnson (CH2MHill)
%
% Revision history:
%   2015-11-01 First version based on Hakeem and Darren code as provided to
100 101
%   Torres by Pierre.
%   2016-01-19 Updated to a stand alone function and general tidy up.
102 103
%   2016-12-14 Updated the help. Also disabled the plot to ease automated
%   runs.
104 105
%   2016-12-22 Fairly major rewrite to make things clearer and less prone
%   to subtle bugs.
106
%   2017-02-16 Fix for direct nesting (no weights needed).
107 108
%   2017-02-20 Fix non-critical bug which added empty cells for the nest
%   elements. Also make the node and element weight values match one another.
Ricardo Torres's avatar
Ricardo Torres committed
109 110 111
%
%==========================================================================

112
[~, subname] = fileparts(mfilename('fullpath'));
Ricardo Torres's avatar
Ricardo Torres committed
113 114

global ftbverbose
115 116 117
if isempty(ftbverbose)
    ftbverbose=false;
end
Ricardo Torres's avatar
Ricardo Torres committed
118 119 120 121
if ftbverbose
    fprintf('\nbegin : %s\n', subname)
end

122 123
TR = triangulation(Mobj.tri, [Mobj.x, Mobj.y]);

124
Nested = Mobj;
125 126 127 128 129 130
Nested.nObs = 0; % number of nodal levels is incremented for each level.

% Make cell arrays to store the element IDs for each nested level as well
% as for the weights on the nodes and elements.
Nested.read_obc_nodes = cell(0);
Nested.read_obc_elems = cell(0);
131
if any(conf.Nested_type == 3)
132 133 134
    Nested.weight_cell = cell(0);
    Nested.weight_node = cell(0);
end
135

Pierre Cazenave's avatar
Pierre Cazenave committed
136 137 138 139 140 141 142
% if ftbverbose
%     figure(1)
%     clf
%     triplot(Nested.tri, Nested.x, Nested.y)
%     axis('equal', 'tight')
%     hold on
% end
Ricardo Torres's avatar
Ricardo Torres committed
143

144 145 146 147 148 149 150 151
% Indices for the output cell arrays which are incremented for each nest
% level and each open boundary.
cumulative_node_idx = 1;
cumulative_elem_idx = 1;
for obc_idx = 1:Mobj.nObs
    % Generate the weights for the elements and nodes.
    if conf.Nested_type(obc_idx) == 3
        if conf.power == 0
152
            weights_nodes = fliplr((1:conf.levels(obc_idx) + 1)./(conf.levels(obc_idx) + 1));
Ricardo Torres's avatar
Ricardo Torres committed
153
        else
154
            weights_nodes = 1:conf.levels(obc_idx) + 1;
Ricardo Torres's avatar
Ricardo Torres committed
155 156
            weights_nodes = 1./weights_nodes.^conf.power;
        end
157 158 159
        % Use the same weights as the nodes, but drop the last level
        % (elements is always 1 smaller).
        weights_elems = weights_nodes(1:end - 1);
Ricardo Torres's avatar
Ricardo Torres committed
160
    end
161

162
    % Save the original open boundary nodes into the nested struct (Nested).
163 164 165 166 167 168 169
    Nested.read_obc_nodes{cumulative_node_idx} = Mobj.read_obc_nodes{obc_idx};

    % Given the current open boundary, find the elements connected to it
    % and give them some weights.
    ti = vertexAttachments(TR, double(Mobj.read_obc_nodes{obc_idx})');
    Nested.read_obc_elems{cumulative_elem_idx} = unique([ti{:}]);

170
    % Save the weights into the nested struct (Nested).
171 172 173 174
    if conf.Nested_type(obc_idx) ~= 1
        Nested.weight_node{cumulative_node_idx} = repmat(weights_nodes(1), 1, length(Nested.read_obc_nodes{cumulative_node_idx}));
        Nested.weight_cell{cumulative_elem_idx} = repmat(weights_elems(1), 1, length(Nested.read_obc_elems{cumulative_elem_idx}));
    end
175

176 177 178 179 180 181
    % Also save the type of open boundary we've got and update the open
    % boundary counter and number of open boundary nodes.
    Nested.nObcNodes(cumulative_node_idx) = length(Nested.read_obc_nodes{cumulative_node_idx});
    Nested.nObs = Nested.nObs + 1;
    Nested.obc_type(cumulative_node_idx) = conf.Nested_type(obc_idx);

182 183 184 185
    if ftbverbose && conf.Nested_type(obc_idx) ~= 1
%         figure(1)
%         scatter(Nested.x(Nested.read_obc_nodes{cumulative_node_idx}), Nested.y(Nested.read_obc_nodes{cumulative_node_idx}), 20, Nested.weight_node{cumulative_node_idx}, 'filled')
%         scatter(Nested.xc(Nested.read_obc_elems{cumulative_elem_idx}), Nested.yc(Nested.read_obc_elems{cumulative_elem_idx}), 20, Nested.weight_cell{cumulative_elem_idx}, 'filled')
186 187
        fprintf('Original open boundary %d\n', obc_idx)
    end
188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205

    % Now we have the original open boundary and the elements connected to
    % it we can move through the levels specified in conf.levels(obc_idx)
    % and repeat the process. Bump the cumulative counters accordingly.
    cumulative_node_idx = cumulative_node_idx + 1;
    cumulative_elem_idx = cumulative_elem_idx + 1;

    for lev = 1:conf.levels(obc_idx)
        % Find the nodes and elements for this level and assign their
        % weights. Use the most recent data in Nested.read_obc_nodes as the
        % anchor from which to work.
        Nested.read_obc_nodes{cumulative_node_idx} = int32(setdiff(unique(Nested.tri(Nested.read_obc_elems{cumulative_elem_idx - 1}, :)), ...
            [Nested.read_obc_nodes{1:cumulative_node_idx - 1}]))';
        ti = vertexAttachments(TR, double(Nested.read_obc_nodes{cumulative_node_idx})');
        Nested.nObs = Nested.nObs + 1;
        Nested.obc_type(cumulative_node_idx) = conf.Nested_type(obc_idx);
        Nested.nObcNodes(cumulative_node_idx) = length(Nested.read_obc_nodes{cumulative_node_idx});

206 207 208
        if conf.Nested_type(obc_idx) ~= 1
            Nested.weight_node{cumulative_node_idx} = repmat(weights_nodes(lev + 1), 1, length(Nested.read_obc_nodes{cumulative_node_idx}));
        end
209 210
        if lev ~= conf.levels(obc_idx)
            Nested.read_obc_elems{cumulative_elem_idx} = setdiff(unique([ti{:}]), [Nested.read_obc_elems{:}]);
211 212 213
            if conf.Nested_type(obc_idx) ~= 1
                Nested.weight_cell{cumulative_elem_idx} = repmat(weights_elems(lev + 1), 1, length(Nested.read_obc_elems{cumulative_elem_idx}));
            end
Ricardo Torres's avatar
Ricardo Torres committed
214
        end
215

216 217 218 219 220 221
        if ftbverbose && conf.Nested_type(obc_idx) ~= 1
%             figure(1)
%             scatter(Nested.x(Nested.read_obc_nodes{cumulative_node_idx}), Nested.y(Nested.read_obc_nodes{cumulative_node_idx}), 20, Nested.weight_node{cumulative_node_idx}, 'filled')
%             if lev ~= conf.levels(obc_idx)
%                 scatter(Nested.xc(Nested.read_obc_elems{cumulative_elem_idx}), Nested.yc(Nested.read_obc_elems{cumulative_elem_idx}), 20, Nested.weight_cell{cumulative_elem_idx}, 'filled')
%             end
222
            fprintf('Nested level %d\n', lev)
223
        end
224 225 226 227

        % Bump the node and element cumulative counters so the next loop
        % dumps everything into the right position in the cell arrays.
        cumulative_node_idx = cumulative_node_idx + 1;
228 229 230
        if lev ~= conf.levels(obc_idx)
            cumulative_elem_idx = cumulative_elem_idx + 1;
        end
Ricardo Torres's avatar
Ricardo Torres committed
231
    end
232 233 234 235 236 237 238 239
    % Check if all possible elements have been correctly identified
    % colapse all nodes into an array and extract all elements connected to
    % those nodes
    % find the elements that are attached to 3 nodes in the nesting region
    %      
    candidate= find(all(ismember(Mobj.tri,[Nested.read_obc_nodes{end}]),2));
    % test if it is different from existing elements
    [truecandidate]=setdiff(candidate,[Nested.read_obc_elems{:}]);
Ricardo Torres's avatar
Ricardo Torres committed
240
    % add to existing list. Catenate will ignore an empty result from
241
    % setdiff
Ricardo Torres's avatar
Ricardo Torres committed
242 243 244
    if ~isempty(truecandidate)
        for dd=1:length(truecandidate)
    Nested.read_obc_elems{end} = cat(2,Nested.read_obc_elems{end},truecandidate(dd));
245
    Nested.weight_cell{end} = cat(2,Nested.weight_cell{end},Nested.weight_cell{end}(end));
Ricardo Torres's avatar
Ricardo Torres committed
246 247
    end
    end
248 249 250
    if ftbverbose
        fprintf('\n')
    end
Ricardo Torres's avatar
Ricardo Torres committed
251 252
end

253
% Update the clunky obc_nodes array with the new node IDs from
254
% Nested.read_obc_nodes.
255 256 257 258 259
for nidx = 1:length(Nested.read_obc_nodes)
    Nested.obc_nodes(nidx, 1:length(Nested.read_obc_nodes{nidx})) = Nested.read_obc_nodes{nidx};
end

if ftbverbose
Pierre Cazenave's avatar
Pierre Cazenave committed
260 261 262
    % figure(1)
    % colorbar
    % title('Nest weights')
Ricardo Torres's avatar
Ricardo Torres committed
263 264
    fprintf('end   : %s \n', subname)
end