Commit 2b826e58 authored by conmccoid's avatar conmccoid
Browse files

Homepage: updated Intersect.m

parent 136f9e22
......@@ -8,21 +8,19 @@ function [P,n,M]=Intersect(X,Y)
% element Y. The numerical challenges are handled by including
% points on the boundary and removing duplicates at the end.
[P1,Q,R,n] = TriIntersect(X,Y);
[P1,Q,R,n]=TriIntersect(X,Y);
if size(P1,2)>1 % if two or more interior points
n=[1 1 1]; % the triangle is candidate for all
n=[1 1 1]; % the triangle is candidate for all
end % neighbors
P=[P1 Q R];
P=SortAndRemoveDoubles(P); % sort counter clock wise
M=zeros(3,3);
if size(P,2)>0
for j=2:size(P,2)-1 % compute interface matrix
M=M+MortarInt(P(:,[1 j j+1]),X,Y);
end
patch(P(1,:),P(2,:),'m') % draw intersection for illustration
% H=line([P(1,:) P(1,1)],[P(2,:),P(2,1)]);
% set(H,'LineWidth',3,'Color','m');
pause(0)
for j=2:size(P,2)-1 % compute interface matrix
M=M+MortarInt(P(:,[1 j j+1]),X,Y);
end
patch(P(1,:),P(2,:),'m') % draw intersection for illustration
pause(0)
end
end
......@@ -31,84 +29,53 @@ function [P,Q,R,n]=TriIntersect(U,V)
% [P,Q,R,n]=TriIntersect(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
X=[v1,v2]\(U-V(:,1)); % Step 1: change of coordinates
Q=zeros(2,9); indQ=zeros(1,9); n=zeros(1,3); indP=ones(1,3);
indR=n; % 1-(0,0); 2-(1,0); 3-(0,1)
for i=1:3 % Step 2: select ref. line
[Qtemp,ntemp,indRtemp,sp]=EdgeIntersect(X,i);
ind=(1:3) + 3*(i-1); % indices of Q for this ref. line
Q(:,ind)=Qtemp; indQ(ind)=ntemp; % update Q and its non-empty index
n=max(n,ntemp); indR=max(indR,indRtemp); % update other indices
indP=indP.*sp; % Step 3: vertices of X in Y
v1=V(:,2)-V(:,1); % vector between V1 and V2
v2=V(:,3)-V(:,1); % vector between V1 and V3
X=[v1,v2]\(U-V(:,1)); % Step 1: change of coordinates
p=[X([2,1],:);1-X(1,:)-X(2,:)]; % values of the parameters p
sp=sign(p)>=0; % sign of the p values
D=cross([X(1,:);1-X(1,:);X(1,:)],[X(2,:);X(2,:);1-X(2,:)],2); % num.s of q0 and 1-q0
Q=zeros(2,9); indQ=zeros(1,9); n=zeros(1,3); indR=ones(3);
for i=1:3 % Step 2: select ref. line
ind=(1:3) + 3*(i-1); % indices of Q for this ref. line
[Q(:,ind),indQ(ind),indR(i,:)]=EdgeIntersect(p,sp,i,D(:,[3,1,2]));
n=max(n,indQ(ind)); % update nhbr index
end
P=U(:,indP==1); % vertices of U in V
Q=V(:,1) + [v1,v2]*Q(:,indQ==1); % Step 4: reverse change of coordinates
R=V(:,indR==1); % vertices of V in U
P=U(:,prod(sp)==1); % vertices of U in V
Q=V(:,1) + [v1,v2]*Q(:,indQ==1);% Step 4: reverse change of coordinates
R=V(:,max(indR==1)); % vertices of V in U
end
function [Q,ind,indR,sp]=EdgeIntersect(X,edge)
function [Q,ind,indR]=EdgeIntersect(p,sp,edge,D)
% EDGEINTERSECT intersection of edges of X with reference line of Y
% [Q,ind,indR,sp]=EdgeIntersect(X,edge) calculates the intersections
% [Q,ind,indR]=EdgeIntersect(p,sp,edge,D) calculates the intersections
% between all edges of X and the reference line of Y indicated by edge
% and stores the values in Q. It also calculates Boolean indicators for
% neighbours of X intersecting Y (ind), vertices of Y inside X (indR),
% and vertices of X inside Y (sp).
if edge==1 % v1
p=X(2,:); q=X(1,:); % parametrization of edge
a=[1;0]; b=[0;0]; % inversion of parametrization
R1=1; R2=2; % vertices of V on edge
elseif edge==2 % v2
p=X(1,:); q=X(2,:); a=[0;1]; b=[0;0]; R1=1; R2=3;
elseif edge==3 % line connecting v1 and v2
p=1-X(1,:)-X(2,:); q=0.5*(1-X(1,:)+X(2,:)); a=[-1;1]; b=[1;0]; R1=2; R2=3;
% neighbours of X intersecting Y (ind) and vertices of Y inside X (indR).
if edge==1 % v1
a=[1;0]; b=[0;0]; % inversion of parametrization
R1=1; R2=2; % vertices of V on edge
D1=D(R1,:); D2=D(R2,:); % relevant numerators
E0=2; E1=3; % adjacent ref lines
elseif edge==2 % v2
a=[0;1]; b=[0;0]; R1=1; R2=3; D1=-D(R1,:); D2=-D(R2,:); E0=1; E1=3;
elseif edge==3 % line connecting v1 and v2
a=[-1;1]; b=[1;0]; R1=2; R2=3; D1=-D(R1,:); D2=D(R2,:); E0=1; E1=2;
end
Q=zeros(2,3); ind=zeros(1,3); indR=ind; sp=sign(p)>=0; qold=NaN;
Q=zeros(2,3); ind=zeros(1,3); indR=ind; D2=sign(D2)>=0;
for i=1:3
j=mod(i,3)+1;
if sp(i)~=sp(j) % vertices on opposite sides of edge?
q0=(q(i)*p(j)-q(j)*p(i))/(p(j)-p(i)); % calculate intersection
if q0>=0 && q0<=1 % intersection lies on edge of Y?
Q(:,i)=q0*a+b;
ind(i)=1;
end
if isnan(qold) % 1st intersection already calculated?
qold=q0;
else
indR(R1)=indR(R1)|sign(qold)~=sign(q0);
indR(R2)=indR(R2)|sign(qold-1)~=sign(q0-1);
end
end
end
end
function P=SortAndRemoveDoubles(P)
% SORTANDREMOVEDOUBLES sort points and remove duplicates
% P=SortAndRemoveDoubles(P); orders polygon corners in P counter
% clock wise and removes duplicates
ep=10*eps; % tolerance for identical nodes
m=size(P,2);
if m>0
c=sum(P,2)/m; % order polygon corners counter
ao=zeros(1,m);
for i=1:m % clockwise
d=P(:,i)-c; ao(i)=angle(d(1)+1i*d(2));
end
[~,id]=sort(ao);
P=P(:,id);
i=1;j=2; % remove duplicates
while j<=m
if norm(P(:,i)-P(:,j))>ep
i=i+1;P(:,i)=P(:,j);j=j+1;
else
j=j+1;
end
j=mod(i,3)+1;
if sp(edge,i)~=sp(edge,j) % vertices on opposite sides of edge?
q0=D1(i)/(p(edge,j)-p(edge,i)); % calculate intersection, eqn (2)
if (q0<0 && sp(E0,i)~=sp(E0,j)) || (~sp(E0,i) && ~sp(E0,j)) % see Sec. 4.1
indR([R1,R2])=indR([R1,R2])+[0.25,0.75];
elseif (D2(i)~=sp(edge,j) && sp(E1,i)~=sp(E1,j)) || (~sp(E1,i) && ~sp(E1,j))
indR([R1,R2])=indR([R1,R2])+[0.75,0.25];
else
indR([R1,R2])=indR([R1,R2])+[0.75,0.75];
Q(:,i)=q0*a+b; ind(i)=1;
end
P=P(:,1:i);
end
end
end
......
Supports Markdown
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