Commit cd3bef1e authored by Ricardo Torres's avatar Ricardo Torres 💬

update of nesting tools

parent 593cca9b
% create a list of the nesting boundary nodes
% subname = 'create_nesting_nodes'
%
% DESCRIPTION:
% READ in the following data:
% - FVCOM OBC file
% - FVCOM Grid file (connectivity + nodes)
% Process data:
% - Extract all nodes belonging to elements at the boundary
% Write output
% - Write Node nesting file
%
% INPUT :
% FVCOM OBC file and Grid file
%
% OUTPUT:
% Node nesting file
%
% EXAMPLE USAGE
% create_nesting_nodes
%
% Author(s):
% Hakeem Johnson (CH2MHILL, Warrington, UK)
% Darren Price(CH2MHILL, Warrington, UK)
%
% Revision history
% Jan 2014: Beta version including the correct nesting format
%==============================================================================
clear;
subname = 'create_nesting_nodes';
global ftbverbose;
if(ftbverbose);
fprintf('\n')
fprintf(['begin : ' subname '\n'])
end
%---------------------------------------------------------
% user specify input data & parameters ...
%---------------------------------------------------------
% work & data folders ...
workDir='\\hand-fs-01\maritime\Projects\Scottish Waters\Calcs\Models\WLLS\mesh10\nesting\';
dataDir='\\hand-fs-01\maritime\Projects\Scottish Waters\Calcs\Models\WLLS\mesh10\mesh\';
basename = 'WLLS_v3smth';
% input file names
meshFile=[basename '_grd.dat'];
obcFile=[basename '_obc.dat'];
% output file
NestNodesFile = [basename '_node_nest.dat'];
%---------------------------------------------------------
% set full file names ...
%---------------------------------------------------------
meshFile = [dataDir meshFile];
obcFile = [dataDir obcFile];
NestNodesFile = [workDir NestNodesFile];
checkfile(meshFile);
checkfile(obcFile);
%---------------------------------------------------------
% Read data into matlab
% 1) grid mesh; 2)boundary nodes; ...
%---------------------------------------------------------
% get fvcom grid file in as mesh object ...
FV_Mesh = read_fvcom_mesh_lonlat(meshFile); % ... from matlab toolbox
% get wave nesting nodes and associated lat/lon coord file
FV_OBC = get_HD_OBC(obcFile);
%---------------------------------------------------------------------
% get elements at boundary & element centres
%---------------------------------------------------------------------
% get elements at boundary
T = FV_Mesh.tri;
X1 = FV_Mesh.lon;
Y1 = FV_Mesh.lat;
P = [X1,Y1];
TR = triangulation(T,P);
vi = FV_OBC.nnodesID;
ti = vertexAttachments(TR,vi);
% arrange as column data; since elements overlap, get unique elements
temp1 = [ti{:}]';
% bdcell = unique(temp1,'stable');
bdcell = unique(temp1);
nCells = numel(bdcell);
nElems = nCells;
bdElem = bdcell;
% get element nodes & arrange in required format
nv= zeros(3*nElems,1);
k = 1;
for i=1:nElems
ielem = bdElem(i);
nv(k) = FV_Mesh.tri(ielem,1);
nv(k+1) = FV_Mesh.tri(ielem,2);
nv(k+2) = FV_Mesh.tri(ielem,3);
k = k+3;
end
% I don't think preserving the order is the problem - remove next line
% nestNodeID= unique(nv,'stable');
nestNodeID= unique(nv);
nestNodes = numel(nestNodeID);
out = ones(nestNodes,3);
for i=1:nestNodes
out(i,1) = i;
out(i,2) = nestNodeID(i);
end
%------------------------------------------------------------------------------
% Dump the file
%------------------------------------------------------------------------------
filename = NestNodesFile;
if(ftbverbose); fprintf('writing FVCOM obc %s\n',filename); end;
%------------------------------------------------------------------------------
% Check output file
%------------------------------------------------------------------------------
fid = fopen(filename,'w');
fprintf(fid,'Nest_Node Number = %12d\n',nestNodes);
fprintf(fid,'No. Node_ID Node_Type\n');
for i=1:nestNodes
fprintf(fid,'%12d %12d %12d\n',out(i,:));
end
fclose(fid);
disp('finished creating nesting file');
\ No newline at end of file
function [Nested]=find_nesting_region(conf,Mobj)
% Creates a nesting structure array for direct/indirect or weighted nesting
%
% function [Nested]=find_nesting_region(conf,M)
%
% DESCRIPTION:
% Uses the Mesh object M and the conf variable to search for the nodes
% 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
% make_input.m type of script. Minimum fields are
% Nesting type (1, 2 == direct nesting, 3 == weighted)
% conf.Nested_type = 3;
% 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
% conf.levels = 5;
% conf.power = determines drop of weights from 1 [0 is linear, 1-4 is 1/conf.levels.^conf.power]
% M = Mesh object
%
% OUTPUT:
% Nested Mesh object with added nesting variables .
%
% EXAMPLE USAGE:
% conf.Nested_type = type of nesting [1, 2 == direct nesting, 3 == weighted]
% conf.levels = number of boundary bands to use for weighted option
% conf.power = determines drop of weights from 1 [0 is linear, 1-4 is 1/conf.levels.^conf.power]
% 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
% in setup_metrics.m
%
% the global variable ftbverbose shows information at run time as well as
% generating a figure with the location of the nesting nodes and elements
%
% [Nested]=find_nesting_region(conf,Mobj)
%
% 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]
% Nested.obc_nodes = matrix with node indices. Each row is a boundary level.
% 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
% Torres by Pierre
% script.
% 2016-01-19 Updated to a stand alone function and general tidy up
% general tidy up.
%
%==========================================================================
dump = dbstack;
subname = dump.name;
clear dump
global ftbverbose
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
M=Mobj;
M.nObcNodes = M.nObcNodes.*(~isnan(conf.levels./conf.levels));
nBC = sum(~isnan(conf.levels./conf.levels));
TR = triangulation(M.tri, [M.x, M.y]);
et = cell(sum(conf.levels), 1);
M.weight_cell = cell(1,nBC*sum(conf.levels));
for oo=1:length(M.read_obc_nodes)
M.weight_node{oo} = ones(1, M.nObcNodes(oo));
end
nn=0;
if ftbverbose
figure,cla
patch('Vertices', [M.x, M.y], 'Faces', M.tri, ...
'Cdata', -M.h, 'edgecolor', 'b', 'facecolor', 'interp');
end
%% The initial OB are first in the order. The nested OB
% follow
etdx=0;
for oo=1:length(Mobj.read_obc_nodes)
if conf.Nested_type(oo) == 3
% probably best to have them as a linear decrease?
if conf.power ==0
% one should be the outermost boundary level which is why the vectors
% are flipped
weights_nodes = fliplr((1:conf.levels(oo))./conf.levels(oo));
weights_nodes(end+1) =0;
weights_elems = fliplr((1:conf.levels(oo)-1)./conf.levels(oo)-1);
weights_elems(end+1) =0;
else
weights_nodes = 1:conf.levels(oo)+1;
weights_nodes = 1./weights_nodes.^conf.power;
weights_elems = 1:conf.levels(oo);
weights_elems = 1./weights_elems.^conf.power;
end
M.weight_node{oo} = repmat(weights_nodes(1), ...
1, M.nObcNodes(oo));
end
for n = 1:conf.levels(oo)
nn=nn+1;etdx=etdx+1;
if (oo>1 && n == 1)
M.read_obc_nodes{nn+1} = Mobj.read_obc_nodes{oo};
M.obc_type(nn+1)=conf.Nested_type(oo);
M.nObcNodes(nn+1) = length(Mobj.read_obc_nodes{oo});
M.obc_nodes(nn+1, 1:Mobj.nObcNodes(oo)) = Mobj.read_obc_nodes{oo};
ti = vertexAttachments(TR, double(M.read_obc_nodes{nn+1})');
et{etdx} = setdiff(unique([ti{:}]),[et{1:end}]);
nn=nn+1;
M.read_obc_nodes{nn+1} = int32(setdiff(unique(M.tri(et{etdx}, :)), ...
[M.read_obc_nodes{1:nn}]))';
M.obc_type(nn+1)=conf.Nested_type(nn);
M.nObs = M.nObs+1;
M.nObcNodes(nn+1) = length(M.read_obc_nodes{nn+1});
M.obc_nodes(nn+1, 1:M.nObcNodes(nn+1)) = M.read_obc_nodes{nn+1};
else
ti = vertexAttachments(TR, double(M.read_obc_nodes{nn})');
et{etdx} = setdiff(unique([ti{:}]),[et{1:end}]);
M.read_obc_nodes{nn+1} = int32(setdiff(unique(M.tri(et{etdx}, :)), ...
[M.read_obc_nodes{1:nn}]))';
M.obc_type(nn+1)=M.obc_type(nn);
M.nObs = M.nObs+1;
M.nObcNodes(nn+1) = length(M.read_obc_nodes{nn+1});
M.obc_nodes(nn+1, 1:M.nObcNodes(nn+1)) = M.read_obc_nodes{nn+1};
end
% Plot nesting region.
if ftbverbose
if n == 1 | nn==1
axis('equal', 'tight')
colormap('gray')
hold on
plot(M.x(M.read_obc_nodes{nn}), M.y(M.read_obc_nodes{nn}), 'wo')
plot(M.x(M.read_obc_nodes{nn+1}), M.y(M.read_obc_nodes{nn+1}), 'ro')
plot(M.xc(et{etdx}), M.yc(et{etdx}), 'wx')
else
plot(M.x(M.read_obc_nodes{nn+1}), M.y(M.read_obc_nodes{nn+1}), 'ro')
plot(M.xc(et{etdx}), M.yc(et{etdx}), 'wx')
end
end
if conf.Nested_type(oo) == 3
M.weight_node{nn+1} = repmat(weights_nodes(n+1), ...
1, M.nObcNodes(nn+1));
M.weight_cell{etdx} = repmat(weights_elems(n), ...
1, length(et{etdx}), 1);
end
end
end
%%
Nested=M;
Nested.read_obc_elems=et;
if ftbverbose
fprintf('end : %s \n', subname)
end
return
This diff is collapsed.
function Mobj=distance_along_BC(Mobj,BCnodes,conf)
% Calculates the distance from coast along the open boundary nodes
%
% Mobj=distance_along_BC(Mobj,BCnodes,conf)
%
% DESCRIPTION:
% Calculates distance from coast along open boundary nodes and
% store them in a matlab mesh object
%
% INPUT
% Mobj = Mesh object structure variable
% BCnodes = indices of nodes located at the boundary
% conf = configuration structure variable with the
% directory where the HJB_Solver_Package is installed
%
% OUTPUT:
% Mobj = matlab structure containing distance data
%
% EXAMPLE USAGE
% Mobj=distance_along_BC(Mobj,BCnodes,conf)
% This function needs the HJB_solver package by Shawn Walker and can be downloaded from Matlab central
% http://www.mathworks.com/matlabcentral/fileexchange/24827-hamilton-jacobi-solver-on-unstructured-triangular-grids
% Author(s):
% Ricardo Torres (Plymouth Marine Laboratory)
%
% Revision history
%
% 2015-11-20 First version
%
%==============================================================================
dump = dbstack;
subname = dump.name;
clear dump
global ftbverbose;
if ftbverbose
fprintf('\nbegin : %s \n', subname)
end
CD=pwd;
% setup HPJ solver to calculate the distance function for the SMS mesh
[~,~,~,bnd] = connectivity([Mobj.x,Mobj.y],Mobj.tri);
% remove nodestring from coast.
% BCnodes=[Mobj.read_obc_nodes{:}];
coast_ind=(BCnodes);
% % calculate distance function
myParam.Max_Tri_per_Star = 20;
myParam.NumGaussSeidel = 40;
myParam.INF_VAL = 1000000;
myParam.TOL = 1e-12;
% in this case, we will assume the standard Euclidean metric
myMetric = [];
myTM.Vtx=[Mobj.x(:),Mobj.y(:)];
myTM.DoFmap=[Mobj.tri];
myTM.NegMask=false(size(bnd));
myBdy.Nodes=coast_ind(:);
myBdy.Data=zeros(size(myBdy.Nodes));
%
cd (conf.HJB_Solver_Package)
%
%
SEmex = SolveEikonalmex(myTM,myBdy,myParam,myMetric);
tic
Mobj.distOB = SEmex.Compute_Soln;
cd(CD)
if ftbverbose
fprintf('end : %s \n', subname)
end
return
function Mobj=distance_to_coast(Mobj)
function Mobj=distance_to_coast(Mobj,conf)
% Calculates the distance from the coast in all mesh nodes
%
% Mobj=distance_to_coast(Mobj,conf)
%
% DESCRIPTION:
% Calculates distance from coast within the domain mesh and
% stores it in a matlab mesh object
%
% INPUT
% Mobj = Mesh object structure variable
% conf = configuration structure variable with the
% directory where the HJB_Solver_Package is installed
%
% OUTPUT:
% Mobj = matlab structure containing distance data
%
% EXAMPLE USAGE
% Mobj=distance_to_coast(Mobj,conf)
% This function needs the HJB_solver package by Shawn Walker and can be downloaded from Matlab central
% http://www.mathworks.com/matlabcentral/fileexchange/24827-hamilton-jacobi-solver-on-unstructured-triangular-grids
% Author(s):
% Ricardo Torres (Plymouth Marine Laboratory)
%
% Revision history
%
% 2015-11-20 First version
%
%==============================================================================
dump = dbstack;
subname = dump.name;
clear dump
global ftbverbose;
if ftbverbose
fprintf('\nbegin : %s \n', subname)
end
CD=pwd;
% setup HPJ solver to calculate the distance function for the SMS mesh
[~,~,~,bnd] = connectivity([Mobj.x,Mobj.y],Mobj.tri);
pv = [Mobj.x(find(bnd)) Mobj.y(find(bnd));Mobj.x(find(bnd(1))) Mobj.y(find(bnd(1)))];
% remove nodestring from coast.
BCnodes=Mobj.read_obc_nodes{1};
BCnodes=[Mobj.read_obc_nodes{:}];
bnd(BCnodes)=0;
coast_ind=find(bnd);
coast_ind(1:max(BCnodes))=[];
% % calculate distance function
myParam.Max_Tri_per_Star = 20;
......@@ -21,12 +59,15 @@ myTM.NegMask=false(size(bnd));
myBdy.Nodes=coast_ind(:);
myBdy.Data=zeros(size(myBdy.Nodes));
%
cd /users/modellers/rito/matlab/HJB_Solver_Package/HJB_Solver_Package
cd (conf.HJB_Solver_Package)
%
%
SEmex = SolveEikonalmex(myTM,myBdy,myParam,myMetric);
tic
Mobj.dist = SEmex.Compute_Soln;
Mobj.dist =SEmex.Compute_Soln;
cd(CD)
if ftbverbose
fprintf('end : %s \n', subname)
end
return
\ No newline at end of file
......@@ -60,11 +60,12 @@ for i = 1:nb
[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
nObcElements(i) = numel(obc_elems{i}(:));
end
Mobj.relaxBC_nodes={reshape([Mobj.read_obc_nodes{:}],[],1)};
Mobj.relaxBC_elems={reshape([obc_elems{:}],[],1)};
Mobj.relaxBC_nodes={[Mobj.read_obc_nodes{:}]};
Mobj.relaxBC_elems={[obc_elems{:}]};
for bb=2:bc_width
nodeIDs = Mobj.tri(Mobj.relaxBC_elems{bb-1},:);
......
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