Commit 79f8e72a authored by Conor McCoid's avatar Conor McCoid

Added comments/help files to CircleMesh and InterfaceMatrixNew

parent a4c7787b
function [N,T,Np] = CircleMesh(n,error)
% CIRCLEMESH constructs a circular mesh
% [N,T] = CircleMesh(n) constructs a grid of n+1 points N,
% containing the origin and n points evenly distributed on the unit
% circle, and a triangulation of these points T such that the unit circle
% is divided into n equal triangles.
%
% [N,T,Np] = CircleMesh(n,error) additionally constructs a second grid of
% points Np with the central point shifted randomly by O(error) and the
% points along the unit circle rotated by a random O(error) number of
% radians. Np uses the same triangulation, T, as N.
ind = 0:(n-1);
N = [cos(ind*2*pi/n), 0; sin(ind*2*pi/n), 0];
T = [(n+1)*ones(n,1), mod(ind'+1,n), mod(ind'+2,n), mod(ind',n), (n+1)*ones(n,1), mod(ind'+2,n)];
T(T==0)=n;
N = [cos(ind*2*pi/n), 0; sin(ind*2*pi/n), 0]; % grid points
T = [(n+1)*ones(n,1), mod(ind'+1,n), mod(ind'+2,n),... % triangle vertices
mod(ind',n), (n+1)*ones(n,1), mod(ind'+2,n)]; % triangle neighbours
% neighbours are listed in the order of the vertices they share:
% T(i,4) shares vertices T(i,1) and T(i,2); T(i,5)--T(i,2) and T(i,3);
% T(i,6)--T(i,3) and T(i,1).
T(T==0)=n; % periodicity: the last grid point/triangle neighbours the first
if nargin==2 && nargout==3
r = error*rand(1)*2*pi/n;
Np= [cos(ind*2*pi/n + r), error*rand(1)-0.5*error; sin(ind*2*pi/n + r), error*rand(1)-0.5*error];
r = error*rand(1)*2*pi/n; % radians rotated
Np= [cos(ind*2*pi/n + r), error*rand(1)-0.5*error;... % x-shifts
sin(ind*2*pi/n + r), error*rand(1)-0.5*error]; % y-shifts
end
\ No newline at end of file
......@@ -81,11 +81,10 @@ end
end
function [P,Q,R,n] = Geometry(U,V)
% GEOMETRY figures out the intersection of the triangles U and V. It does
% so by first putting all variables in terms of two of the lines of the
% triangle V, decides if the points of U lie within or outside V, then
% calculates the intersections of the lines of U with those of V and again
% deciding if they are inside or outside V.
% GEOMETRY intersection of two triangles
% [P,Q,R,n] = Geometry(U,V) gives the vertices of U inside V (P),
% intersections between the edges of U and those of V (Q), the vertices
% of V inside U (R), and the neighbours of U that also intersect V (n).
v1=V(:,2)-V(:,1); % vector between V1 and V2
v2=V(:,3)-V(:,1); % vector between V1 and V3
......@@ -102,7 +101,7 @@ for i=1:3
d12=v2'*r; % projection of r in v2 direction
X(1,i)=(d11*d02-d01*d12)/id; % amount of r in v1 direction
X(2,i)=(d00*d12-d01*d02)/id; % amount of r in v2 direction
% r = u*v0 + v*v1
% end result: r = X(1,i)*v1 + X(2,i)*v2
end
%---Intersections (sec.3)---%
......@@ -130,15 +129,19 @@ end
end
function [Q,ind,indR,sp]=EdgeIntersect(X,edge)
% EdgeIntersect computes the intersection of an edge of triangle U with the
% reference triangle. The specific reference line intersected is defined
% by edge.
% EDGEINTERSECT intersection of edges of X with reference line of Y
% [Q,ind,indR,sp] = EdgeIntersect(X,edge) calculates the intersections
% between all edges of X and the reference line of Y indicated by edge
% and stores the values in Q.
% Boolean indicators for neighbours of U intersecting V, vertices of V
% inside U, and vertices of U inside V are stored in ind, indR and sp,
% respectively.
% Particulars for each edge of reference triangle
if edge==1 % v1
Param1 = @(x,y) y; Param2 = @(x,y) x;
Retrieve = @(q) [q;0];
R1 = 1; R2 = 2;
Param1 = @(x,y) y; Param2 = @(x,y) x; % parametrization of edge
Retrieve = @(q) [q;0]; % inversion of parametrization
R1 = 1; R2 = 2; % vertices of V on edge
elseif edge==2 % v2
Param1 = @(x,y) x; Param2 = @(x,y) y;
Retrieve = @(q) [0;q];
......@@ -154,19 +157,19 @@ p = Param1(X(1,:),X(2,:)); q = NaN*ones(1,3);
sp= sign(p)>=0;
for i = 1:3
j = mod(i,3) + 1;
if sp(i) ~= sp(j)
if isnan(q(i))
if sp(i) ~= sp(j) % test if vertices of X lie on opposite sides of edge
if isnan(q(i))% test if q-coord. of vertex already calculated
q(i) = Param2(X(1,i),X(2,i));
end
if isnan(q(j))
q(j) = Param2(X(1,j),X(2,j));
end
q0=(q(i)*p(j) - q(j)*p(i))/(p(j)-p(i));
if q0>=0 && q0<=1
q0=(q(i)*p(j) - q(j)*p(i))/(p(j)-p(i)); % calculate intersection
if q0>=0 && q0<=1 % test if intersection lies on edge of Y
Q(:,i) = Retrieve(q0);
ind(i) = 1;
end
if exist('qold','var')
if exist('qold','var') % test if 1st intersect. already calculated
if sign(qold)~=sign(q0)
indR(R1) = 1;
end
......
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