convert.m

解决麦克斯韦的matlab源文件

function mesh = convert(x,y,tri)
% convert: Conversion of simple mesh data structure to 
% an extended data format describing an unstructured mesh
% Usage:
%          mesh = convert(x,y,tri)
%
% x -> vector of length N containing the x-coordinates of the vertices
% y -> vector of length N containing the y-coordinates of the vertices
% tri -> Nx3 matrix whose rows contain the vertex indices [1,...,N] of
%        each triangle's vertices.
%
% Data structure for describing a mesh
%
% mesh.Nv       -> # of vertices
% mesh.Ne       -> # of edges
% mesh.Nt       -> # of triangles
% mesh.vt(Nv,2) -> (x,y)-coordinates of vertices
% mesh.vl(Nv,2) -> boundary flag for vertex
%                  vl(k,1) > 0 -> vertex k is located on a vertical boundary
%                  vl(k,2) > 0 -> vertex k is located on a horizontal boundary
% mesh.ep(Ne,2) -> index numbers of endpoints of edges
% mesh.ebfl(Ne,1) -> boundary flags for edges
%                    ebfl(k) = 1 -> edge is located on a vertical boundary
%                    ebfl(k) = 2 -> edge is located on a horizontal boundary
% mesh.etr(Ne,2) -> index numbers of triangles adjacent to an edge
% mesh.trv(Nt,3) -> index numbers of vertices of triangles
% mesh.tre(Nt,3) -> index numbers of edges of triangles
%                    tre(k,i) refers to the edge opposite the i-th vertex 
%                    of the k-th triangle 
% mesh.treo(Nt,3) -> relative orientations of edges of triangles
%                   (all vertices of triangles have to be arranged in
%                    counterclockwise fashion)

N = length(x);
if (length(y) ~= N), error('Size mismatch of x and y'); end;
if (size(tri,2) ~= 3), error('Three vertices required for each triangle'); end

mesh.Nv = N; % Total number of vertices in the fl_mesh
mesh.Nt = size(tri,1); % Total number of triangles in the fl_mesh

% Copy vertex coordinates

mesh.vt = zeros(N,1);
for i=1:N
  mesh.vt(i,1) = x(i);
  mesh.vt(i,2) = y(i);
end

% Copy vertex numbers of triangles
% make sure that triangles have uniform orientation

mesh.trv = [];

for i = tri'
  if (det(mesh.vt([i(2) i(3)],:) - mesh.vt([i(1) i(1)],:)) > 0)
    mesh.trv = [mesh.trv; i'];
  else
    mesh.trv = [mesh.trv; i([1 3 2])'];
  end
end

% Retrieve edges from incidence information for triangles and vertices

EL = [sort([mesh.trv(:,1) mesh.trv(:,2)]')' (1:mesh.Nt)' 3*ones(mesh.Nt,1); ...
      sort([mesh.trv(:,1) mesh.trv(:,3)]')' (1:mesh.Nt)' 2*ones(mesh.Nt,1); ...
      sort([mesh.trv(:,2) mesh.trv(:,3)]')' (1:mesh.Nt)' 1*ones(mesh.Nt,1)];

% mesh.Nt, size(EL)

[B,i1,j1] = unique(EL(:,1:2),'rows');
[B,i2,j3] = unique(flipud(EL(:,1:2)),'rows');
EM = [EL(i1,1:4) EL((3*mesh.Nt)-i2+1,3:4)];

% EM(k,1) -> first endpoint of edge k
% EM(k,2) -> second endpoint of edge k
% EM(k,3) -> index of first adjacent triangle of edge k
% EM(k,4) -> local index of edge k in triangle EM(k,3)
% EM(k,5) -> index of second adjacent triangle of edge k
%            (if edge is on the boundary EM(k,3) = EM(k,5) !)
% EM(k,6) -> local index of edge k in triangle EM(k,5)

mesh.ep = EM(:,1:2); mesh.etr = EM(:,[3 5]); 
mesh.Ne = size(EM,1);

% Set boundary flags for edges

mesh.ebfl = zeros(mesh.Ne,1);
mesh.ebfl(find(EM(:,3) == EM(:,5))) = 1;

% Run through the boundary edges and set boundary flags for vertices

mesh.vl = zeros(mesh.Nv,2);
for i= find(mesh.ebfl(:,1) == 1)'
  edge_vertices = mesh.ep(i,:)';
  edge_direction = (mesh.vt(edge_vertices(2),:) - mesh.vt(edge_vertices(1),:))';
  if (abs(edge_direction(1)) < eps)
    mesh.vl(edge_vertices,1) = mesh.vl(edge_vertices,1)+1;
  elseif (abs(edge_direction(2)) < eps)
    mesh.vl(edge_vertices,2) = mesh.vl(edge_vertices,2)+1;
  else error('Code can only cope with either horizontal or vertical edges'); end
end

% Set edges for triangles

mesh.tre = zeros(mesh.Nt,3);
for i=1:mesh.Ne
  mesh.tre(EM(i,3),EM(i,4)) = i;
  mesh.tre(EM(i,5),EM(i,6)) = i;
end

% Determine orientation of edges w.r.t. to triangle

mesh.treo = zeros(mesh.Nt,3);

for i=1:mesh.Nt
  ed1 = [mesh.trv(i,2) mesh.trv(i,3)];
  ed2 = [mesh.trv(i,3) mesh.trv(i,1)];
  ed3 = [mesh.trv(i,1) mesh.trv(i,2)];
  if (ed1 == mesh.ep(mesh.tre(i,1),:)), mesh.treo(i,1) = 1;
  elseif (fliplr(ed1) == mesh.ep(mesh.tre(i,1),:)), mesh.treo(i,1) = -1;
  else fprintf('ERROR Edge 1 of triangle %d!',i); ed1, mesh.ep(mesh.tre(i,1),:)
    error('Inconsistent edge');
  end
  ed2 = [mesh.trv(i,3) mesh.trv(i,1)];
  if (ed2 == mesh.ep(mesh.tre(i,2),:)), mesh.treo(i,2) = 1;
  elseif (fliplr(ed2) == mesh.ep(mesh.tre(i,2),:)), mesh.treo(i,2) = -1;
  else fprintf('ERROR Edge 2 of triangle %d!',i); ed2, mesh.ep(mesh.tre(i,2),:)
    error('Inconsistent edge');
  end
  ed3 = [mesh.trv(i,1) mesh.trv(i,2)];
  if (ed3 == mesh.ep(mesh.tre(i,3),:)), mesh.treo(i,3) = 1;
  elseif (fliplr(ed3) == mesh.ep(mesh.tre(i,3),:)), mesh.treo(i,3) = -1;
  else fprintf('ERROR Edge 3 of triangle %d!',i); ed3, mesh.ep(mesh.tre(i,3),:)
    error('Inconsistent edge');
  end
end