Commit b6cf40f7 authored by Pierre Cazenave's avatar Pierre Cazenave

Tool from the mesh2D toolbox to calculate unstructured grid connectivity...

Tool from the mesh2D toolbox to calculate unstructured grid connectivity (including finding edge nodes, which is useful for plotting the coastline)
parent 42c42b3b
function [e,te,e2t,bnd] = connectivity(p,t)
% CONNECTIVITY: Assemble connectivity data for a triangular mesh.
%
% The edge based connectivity is built for a triangular mesh and the
% boundary nodes identified. This data should be useful when implementing
% FE/FV methods using triangular meshes.
%
% [e,te,et2,bnd] = connectivity(p,t);
%
% p : Nx2 array of nodes coordinates, [x1,y1; x2,y2; etc]
% t : Mx3 array of triangles as indices, [n11,n12,n13; n21,n22,n23; etc]
%
% e : Kx2 array of unique mesh edges - [n11,n12; n21,n22; etc]
% te : Mx3 array of triangles as indices into E, [e11,e12,e13;
% e21,e22,e23; etc]
% e2t : Kx2 array of triangle neighbours for unique mesh edges -
% [t11,t12; t21,t22; etc]. Each row has two entries corresponding to
% the triangle numbers associated with each edge in E. Boundary
% edges have e2t(i,2)=0.
% bnd : Nx1 logical array identifying boundary nodes. P(i,:) is a boundary
% node if BND(i)=TRUE.
%
% See also MESH2D, REFINE
% Darren Engwirda - 2007
if nargin<2
error('Wrong number of inputs');
end
if nargout>4
error('Wrong number of outputs');
end
if numel(p)~=2*size(p,1)
error('P must be an Nx2 array');
end
if numel(t)~=3*size(t,1)
error('T must be an Mx3 array');
end
if any(t(:)<1) || max(t(:)>size(p,1))
error('Invalid T');
end
% Unique mesh edges as indices into P
numt = size(t,1);
vect = 1:numt; % Triangle indices
e = [t(:,[1,2]); t(:,[2,3]); t(:,[3,1])]; % Edges - not unique
[e,j,j] = unique(sort(e,2),'rows'); % Unique edges
te = [j(vect), j(vect+numt), j(vect+2*numt)]; % Unique edges in each triangle
% Edge-to-triangle connectivity
% Each row has two entries corresponding to the triangle numbers
% associated with each edge. Boundary edges have e2t(i,2)=0.
nume = size(e,1);
e2t = zeros(nume,2);
for k = 1:numt
j = 1;
while j<=3
ce = te(k,j);
if e2t(ce,1)==0
e2t(ce,1) = k;
else
e2t(ce,2) = k;
end
j = j+1;
end
end
% Flag boundary nodes
bnd = false(size(p,1),1);
bnd(e(e2t(:,2)==0,:)) = true; % True for bnd nodes
end % connectivity()
\ No newline at end of file
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