%This is software to paper
%A 3D Crouzeix-Raviart mortar finite element
%by authors Leszek Marcinkowski, Talal Rahman, Jan Valdman
%send comments to Jan Valdman (Jan.Valdman@gmail.com)
clear all
level_max=2; %set the finest mesh level
%X1_coarse=[0 0 0; 1 0 0; 0 1 0; 1 1 0; 0 0 1; 1 0 1; 0 1 1; 1 1 1];
%T1_coarse=[1 3 5 6; 1 2 3 6; 2 3 4 6; 3 5 6 7; 3 6 7 8; 3 4 6 8];
%X2_coarse=X1_coarse-[zeros(size(X1_coarse,1),2) ones(size(X1_coarse,1),1)];
%T2_coarse=[1 3 5 6; 1 2 3 6; 2 3 4 6; 3 5 6 7; 3 6 7 8; 3 4 6 8];
dx = 0:1:1; dy = 0:1:1; dz = 0:1:1;
[x,y,z] = meshgrid(dx,dy,dz);
X1_coarse = [x(:) y(:) z(:)]; T1_coarse = delaunay3(x,y,z,{'QJ','Pp'}); %master geometry
dx = 0:1:1; dy = 0:1:1; dz = -1:1:0;
[x,y,z] = meshgrid(dx,dy,dz);
X2_coarse = [x(:) y(:) z(:)]; T2_coarse = delaunay3(x,y,z,{'QJ','Pp'}); %slave geometry
visualization=0;
for level=0:level_max
X1=X1_coarse; T1=T1_coarse;
X2=X2_coarse; T2=T2_coarse;
for i=1:level+1 %decides how fine the master triangulation is.
[X1,T1]=tetrarefine3(X1,T1);
end
for i=1:level %decides how fine the slave triangulation is.
[X2,T2]=tetrarefine3(X2,T2);
end
np1 = size(X1,1); nt1 = size(T1,1);
np2 = size(X2,1); nt2 = size(T2,1);
%visualization of the coarsest meshes only
if (level < 2)
figure(level+1);
subplot(2,1,1);
show_mesh(T1,X1); text(2,0,0.5,'master subdomain');
subplot(2,1,2); show_mesh(T2,X2); text(2,0,-0.5,'slave subdomain');
visualization=visualization+1;
end
%vertices rotated so that the last vertex is lying of the intersection
[T1, T1_Plane_index]=rotate_intersection_elements(T1,X1);
[T2, T2_Plane_index]=rotate_intersection_elements(T2,X2);
%face numbers and face middle points computed
[T2F1, F1]=getFaces(T1); nF1=size(F1,1); middle_points1=face_middlepoint(X1,F1);
[T2F2, F2]=getFaces(T2); nF2=size(F2,1); middle_points2=face_middlepoint(X2,F2);
%MASS matrices only useful for L2 norm error computation
MASS1=stima_massCR_3D(T1,X1,T2F1);
MASS2=stima_massCR_3D(T2,X2,T2F2);
%solution of subproblems
A1=stima_gradgradCR_3D(T1,X1,T2F1,F1);
A2=stima_gradgradCR_3D(T2,X2,T2F2,F2);
u1_exact=exact_solution(middle_points1);
u2_exact=exact_solution(middle_points2);
b1=rhs_CR_3D(T1,X1,T2F1,F1);
b2=rhs_CR_3D(T2,X2,T2F2,F2);
%find internal faces for subproblems
T2F1_sorted=sort(T2F1(:));
internal_faces1=T2F1_sorted(find(T2F1_sorted-[0; T2F1_sorted(1:end-1)]==0));
T2F2_sorted=sort(T2F2(:));
internal_faces2=T2F2_sorted(find(T2F2_sorted-[0; T2F2_sorted(1:end-1)]==0));
fprintf('level %d:',level);
fprintf(' %d ',size(A1,2)+size(A2,2));
fprintf('(%d',size(A1,2)); fprintf('%d) ',size(A2,2));
fprintf(' DOF total (master|slave subdomain) \n');
%master subproblem solution
%u1=zeros(size(A1,1),1);
%[solution1,flag]=pcg(A1(internal_faces1,internal_faces1),b1(internal_faces1),1e-6,10000);
%u1(internal_faces1)=solution1;
%error1_H1=sqrt((u1-u1_exact)'*A1*(u1-u1_exact));
%slave subproblem solution
%u2=zeros(size(A2,1),1);
%[solution2, flag]=pcg(A2(internal_faces2,internal_faces2),b2(internal_faces2),1e-6,10000);
%u2(internal_faces2)=solution2;
%error2_H1=sqrt((u2-u2_exact)'*A2*(u2-u2_exact));
%intersection information
T1_onIntersection=T1(T1_Plane_index,:); %tetrahedras on intersection from master side
T2_onIntersection=T2(T2_Plane_index,:); %tetrahedras on intersection from slave side
T2F1_onIntersection=T2F1(T1_Plane_index,:); %faces of tetrahedras on intersection from master side
T2F2_onIntersection=T2F2(T2_Plane_index,:); %faces of tetrahedras on intersection from slave side
FAC1=T2F1_onIntersection(:,4); %faces on intersection from master side
FAC2=T2F2_onIntersection(:,4); %faces on intersection from slave side
TRI1=T1_onIntersection(:,1:3); %triangles on intersection from master side
TRI2=T2_onIntersection(:,1:3); %triangles on intersection from slave side
[TRI2E1,E1]=getEdges(TRI1); %edges on intersection from the master side
E2TRI1=entryInWhichRows(TRI2E1); %edges to triangles mapping
%[element2edges, edge2nodes]=getEdges(elements);
%edge2elements=entryInWhichRows(element2edges);
%figure(3); show_vertices(unique(TRI1)',X1);
%figure(4); show_vertices(unique(TRI2),X2);
%figure(4); show_mesh(T1_onIntersection(:,1:3),X1(:,1:2));
%text(middle_points1(F1_onIntersection(:,4),1),middle_points1(F1_onIntersection(:,4),2),int2str(F1_onIntersection(:,4)));
%figure(5); show_mesh(T2_onIntersection(:,1:3),X2(:,1:2));
%text(middle_points2(F2_onIntersection(:,4),1),middle_points2(F2_onIntersection(:,4),2),int2str(F2_onIntersection(:,4)));
%triangle neighbouring with a triangle on and intersection and its face numbers
TRI1_neighboursTRI=zeros(size(TRI1,1),3);
TRI1_neighboursFAC=zeros(size(TRI1,1),3);
for i=1:size(TRI1,1);
neighbours=setdiff(unique(E2TRI1(TRI2E1(i,:),:)),[0 i]);
TRI1_neighboursTRI(i,1:size(neighbours,2))=neighbours;
TRI1_neighboursFAC(i,1:size(neighbours,2))=FAC1(neighbours);
end
%figure(20+level);clf;
%show_vertices(TRI1(2,:),X1);
%figure(30+level);clf;
%show_vertices(unique(TRI1(setdiff(TRI1_neighboursTRI(2,:),0),:)),X1);
row_index=0;
IX=zeros(size(T1_onIntersection,1)*size(T2_onIntersection,1),4); IY=IX; IZ=IX;
for i=1:size(T1_onIntersection,1)
triangle1=TRI1(i,:); %triangle on the master side
%preparation of reduction of internal master degrees of freedom
if (TRI1_neighboursFAC(i,3))
middle_point=middle_points1(FAC1(i),1:2);
middle_point_neighbour1=middle_points1(TRI1_neighboursFAC(i,1),1:2);
middle_point_neighbour2=middle_points1(TRI1_neighboursFAC(i,2),1:2);
middle_point_neighbour3=middle_points1(TRI1_neighboursFAC(i,3),1:2);
face=T2F1_onIntersection(i,4);
face_neighbour1=TRI1_neighboursFAC(i,1);
face_neighbour2=TRI1_neighboursFAC(i,2);
face_neighbour3=TRI1_neighboursFAC(i,3);
else
ii=TRI1_neighboursTRI(i,1);
if (TRI1_neighboursFAC(ii,3)==0)
ii=TRI1_neighboursTRI(i,2);
if (ii==0)
fprintf('Geometry problem, program stops now! \n'); %only problem if the master triangulation is very coarse!
pause
end
end
middle_point=middle_points1(FAC1(i),1:2);
middle_point_neighbour1=middle_points1(TRI1_neighboursFAC(ii,1),1:2);
middle_point_neighbour2=middle_points1(TRI1_neighboursFAC(ii,2),1:2);
middle_point_neighbour3=middle_points1(TRI1_neighboursFAC(ii,3),1:2);
face=T2F1_onIntersection(i,4);
face_neighbour1=TRI1_neighboursFAC(ii,1);
face_neighbour2=TRI1_neighboursFAC(ii,2);
face_neighbour3=TRI1_neighboursFAC(ii,3);
end
%z component is not needed since the intersection satisfies z=0
x1=X1(triangle1,1); x1_closed=[x1; x1(1)];
y1=X1(triangle1,2); y1_closed=[y1; y1(1)];
for j=1:size(T2_onIntersection,1)
triangle2=TRI2(j,:); %triangle on the slave side
x2=X2(triangle2,1); x2_closed=[x2; x2(1)];
y2=X2(triangle2,2); y2_closed=[y2; y2(1)];
%the use of matlab functions inpolygon and polyxpoly
%slows down this code, one should consider to replace
%by own code which finds and intersection of two triangles
%triangle1 and triangle2
%start of slow part
list=[];
if inpolygon(x1(1),y1(1),x2,y2)
list=[list; [x1(1) y1(1)]];
end
if inpolygon(x1(2),y1(2),x2,y2)
list=[list; [x1(2) y1(2)]];
end
if inpolygon(x1(3),y1(3),x2,y2)
list=[list; [x1(3) y1(3)]];
end
if inpolygon(x2(1),y2(1),x1,y1)
list=[list; [x2(1) y2(1)]];
end
if inpolygon(x2(2),y2(2),x1,y1)
list=[list; [x2(2) y2(2)]];
end
if inpolygon(x2(3),y2(3),x1,y1)
list=[list; [x2(3) y2(3)]];
end
[XI,YI] =polyxpoly(x1_closed,y1_closed,x2_closed,y2_closed,'unique');
intersection_points=[list; [XI,YI]];
intersection_points=unique(intersection_points,'rows');
%end of slow part
if (size(intersection_points,1)>2)
row_index=row_index+1;
[intersection_points_x,intersection_points_y]=poly2cw(intersection_points(:,1),intersection_points(:,2));
ispolycw(intersection_points_x,intersection_points_y);
area_intersection=polyarea(intersection_points_x,intersection_points_y);
[center_intersection_x, center_intersection_y] =center_polygon(intersection_points_x,intersection_points_y);
center_intersection=[center_intersection_x, center_intersection_y 0];
%reducing internal master degrees of freedom using
%new interpolation operator
IX1(row_index,:)=[T2F2_onIntersection(j,4) T2F2_onIntersection(j,4) T2F2_onIntersection(j,4) T2F2_onIntersection(j,4)];
IY1(row_index,:)=[face_neighbour1 face_neighbour2 face_neighbour3 face];
inverse=inv([middle_point_neighbour1 1 0; middle_point_neighbour2 1 0; middle_point_neighbour3 1 0; 0 0 0 1]);
IZ1(row_index,:)=area_intersection*[center_intersection(1:2)-middle_point 0 1]*inverse;
%without interpolation operator, standart mortar method
IX2(row_index,:)=[T2F2_onIntersection(j,4) T2F2_onIntersection(j,4) T2F2_onIntersection(j,4) T2F2_onIntersection(j,4)];
IY2(row_index,:)=T2F1_onIntersection(i,:);
IZ2(row_index,:)=area_intersection*[center_intersection 1]*inv([middle_points1(T2F1_onIntersection(i,:)',:) [1; 1; 1; 1]]);
end
end
end
IX1=IX1(1:row_index,:); IY1=IY1(1:row_index,:); IZ1=IZ1(1:row_index,:); M1=sparse(IX1,IY1,IZ1,nF2,nF1); %new interpolation operator
IX2=IX2(1:row_index,:); IY2=IY2(1:row_index,:); IZ2=IZ2(1:row_index,:); M2=sparse(IX2,IY2,IZ2,nF2,nF1); %standard mortar method
clear IX1 IY1 IZ1 IX2 IY2 IZ2
%slave mass matrix, same for both methods
areas_intersection2=zeros(size(T2_onIntersection,1),1);
for j=1:size(T2_onIntersection,1)
triangle2=T2_onIntersection(j,1:3);
x2=X2(triangle2,1);
y2=X2(triangle2,2);
[dummy,dummy,areas_intersection2(j)]=center_polygon(x2,y2);
end
areas_intersection2_inverse=areas_intersection2.^(-1);
S=sparse(T2F2_onIntersection(:,4),T2F2_onIntersection(:,4),areas_intersection2,nF2,nF2);
%testing
%M1_intersection= M1(T2F2_onIntersection(:,4),unique(T2F1_onIntersection(:)));
%S_inverse=sparse(T2F2_onIntersection(:,4),T2F2_onIntersection(:,4),areas_intersection2_inverse,nF2,nF2);
%S_inverse_intersection=S_inverse(T2F2_onIntersection(:,4),T2F2_onIntersection(:,4));
%(S_inverse_intersection*M1_intersection*ones(size(M1_intersection,2),1
%))'; %should give vector of ones
boundary_faces2=setdiff((1:nF2)',union(internal_faces2,T2F2_onIntersection(:,4)));
internal_faces=union(internal_faces1,nF1+internal_faces2);
internal_faces=union(internal_faces,FAC1);
internal_faces=union(internal_faces,nF1+FAC2);
internal_faces=union(internal_faces,nF1+nF2+1:nF1+nF2+size(M1,1));
%solution new interpolation operator
A=[A1, sparse(size(A1,1),size(A2,2)), M1'; ...
sparse(size(A2,1),size(A1,2)), A2, -S';...
M1, -S, sparse(size(S,2),size(S,1))];
b=[b1; b2; zeros(size(M1,1),1)];
solution=zeros(size(A,1),1);
[solution12,flag]=minres(A(internal_faces,internal_faces),b(internal_faces),1e-6,10000);
clear A;
solution(internal_faces)=solution12;
u1_mortar=solution(1:nF1);
u2_mortar=solution(nF1+1:nF1+nF2);
%lambda_mortar=solution(nF1+nF2+1:end); %this Lagrange multiplier is not needed
%error computation new interpolation operator
error_H1=sqrt((u1_mortar-u1_exact)'*A1*(u1_mortar-u1_exact)+(u2_mortar-u2_exact)'*A2*(u2_mortar-u2_exact));
error_L2=sqrt((u1_mortar-u1_exact)'*MASS1*(u1_mortar-u1_exact)+(u2_mortar-u2_exact)'*MASS2*(u2_mortar-u2_exact));
fprintf(' new interpolation operator method - error: %f (H1 norm)',error_H1); fprintf(' %f (L2 norm) \n',error_L2);
%solution standard mortar method
A=[A1, sparse(size(A1,1),size(A2,2)), M2'; ...
sparse(size(A2,1),size(A1,2)), A2, -S';...
M2, -S, sparse(size(S,2),size(S,1))];
b=[b1; b2; zeros(size(M2,1),1)];
solution=zeros(size(A,1),1);
[solution12,flag]=minres(A(internal_faces,internal_faces),b(internal_faces),1e-6,10000);
clear A;
solution(internal_faces)=solution12;
u1_mortar=solution(1:nF1);
u2_mortar=solution(nF1+1:nF1+nF2);
lambda_mortar=solution(nF1+nF2+1:end);
%error computation standard mortar method
error_H1=sqrt((u1_mortar-u1_exact)'*A1*(u1_mortar-u1_exact)+(u2_mortar-u2_exact)'*A2*(u2_mortar-u2_exact));
error_L2=sqrt((u1_mortar-u1_exact)'*MASS1*(u1_mortar-u1_exact)+(u2_mortar-u2_exact)'*MASS2*(u2_mortar-u2_exact));
fprintf(' standard mortar method - error: %f (H1 norm)',error_H1); fprintf(' %f (L2 norm) \n',error_L2);
end
%u1_onIntersection=u1(FAC1); u2_onIntersection=u2(FAC2);
%figure(10)
%plot3(middle_points1(FAC1,1),middle_points1(FAC1,2),u1_onIntersection,'x');
%figure(11)
%plot3(middle_points2(FAC2,1),middle_points2(FAC2,2),u2_onIntersection,
%'x')
