### Updates to InterfaceMatrixNew.m, putting it in line with founding paper

parent e438033f
 ... ... @@ -79,13 +79,12 @@ end end function [P,Q,R,n] = Geometry(U,V) % GEOMETRY figures out the intersection of the triangles X and Y. It does % 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 Y, decides if the points of X lie within or outside Y, then % calculates the intersections of the lines of X with those of Y and again % deciding if they are inside or outside Y. % 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. ind=[0,0,0]; v1=V(:,2)-V(:,1); % vector between V1 and V2 v2=V(:,3)-V(:,1); % vector between V1 and V3 d00=v1'*v1; % norm of v1 ... ... @@ -93,7 +92,7 @@ d01=v1'*v2; d11=v2'*v2; % norm of v2 id =d00*d11 - d01*d01; % this value is used to normalize quantities %---Interior points---% %---Change of Coordinates (sec.2)---% X = zeros(2,3); for i=1:3 r=U(:,i)-V(:,1); % vector between ith vertex of triangle U and V1 ... ... @@ -102,63 +101,65 @@ for i=1:3 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 % if u and v are positive and their sum less than 1 then v2 points % inside the simplex created by the vectors v0 and v1 if X(1,i)>=0 && X(2,i)>=0 && sum(X(:,i))<=1 % also include nodes on the boundary ind(i) = 1; end end P = U(:,ind==1); %---Intersections---% %---Intersections (sec.3)---% if nargout>1 Q = zeros(2,9); Q = zeros(2,9); indQ = zeros(1,9); n = [0,0,0]; n = [0,0,0]; indR = n; % 1-(0,0); 2-(1,0); 3-(0,1) indP = ones(1,3); for i=1:3 [Qtemp,ntemp,indRtemp] = EdgeIntersect(X,i); ind = (1:3) + 3*(i-1); [Qtemp,ntemp,indRtemp,sp] = EdgeIntersect(X,i); ind = (1:3) + 3*(i-1); Q(:,ind) = Qtemp; indQ(ind)= ntemp; n = max(n,ntemp); indR = max(indR,indRtemp); n = max(n,ntemp); indR = max(indR,indRtemp); indP = indP.*sp; end P = U(:,indP==1); Q = V(:,1) + [v1,v2]*Q(:,indQ==1); R = V(:,indR==1); end end function [Q,ind,indR]=EdgeIntersect(X,edge) 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. % Particulars for each edge of reference triangle if edge==1 % v0 if edge==1 % v1 Param1 = @(x,y) y; Param2 = @(x,y) x; Retrieve = @(q) [q;0]; R1 = 1; R2 = 2; elseif edge==2 % v1 elseif edge==2 % v2 Param1 = @(x,y) x; Param2 = @(x,y) y; Retrieve = @(q) [0;q]; R1 = 1; R2 = 3; elseif edge==3 % line connecting v0 and v1 elseif edge==3 % line connecting v1 and v2 Param1 = @(x,y) 1-x-y; Param2 = @(x,y) 0.5*(1-x+y); Retrieve = @(q) [1-q;q]; R1 = 2; R2 = 3; end Q = zeros(2,3); ind = zeros(1,3); indR = ind; p = Param1(X(1,:),X(2,:)); 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) qi=Param2(X(1,i),X(2,i)); qj=Param2(X(1,j),X(2,j)); q0=(qi*p(j) - qj*p(i))/(p(j)-p(i)); if isnan(q(i)) 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 Q(:,i) = Retrieve(q0); ind(i) = 1; ... ...