% Assembly of the global stiffness matrix 
% for the discontinuous coefficients Laplace problem
% Default discontinuity region: 0.25<=x<=0.75; 0.25<=y<=0.75
%
%    ___  __
%   |\  |\  |
%   |2\3| \ |
%   |__\|6_\| <- j      Dirichlet BC on {x=1, y=1}
%   |\ 1|\  |           Neumann   BC on {x=0, y=0}
%   | \ |4\5|
%   |__\|__\|
%
%       i
%
function [A]=matgen_disco(times,Vdsr)

Ndsr     = 1;   % one discontinuity region
dsr(1,1) = 0.25d0;
dsr(1,2) = 0.75d0;
dsr(1,3) = 0.25d0;
dsr(1,4) = 0.75d0;
%off      Vdsr(1)   = 1.0d-4;
%off      Vdsr(1)   = 1.0d0;

disp(['Size of discontinuity: ',num2str(Vdsr)])

%%%%%%%%%%%%%%%%%% Generate the matrix A %%%%%%%%%%%%%%
init=2;
%off times=input('How many times to refine (in one direction)?<2^n> ');
nanr = init^times;  hstep=1/nanr; h2=hstep*hstep; nall=nanr^2;
A=spalloc(nall,nall,7*nall);

%
for i=1:nanr,
 for j=1:nanr,
    row=(i-1)*nanr+j;
%________________________________________________
   if i-1~=0
      worka=disco_stiff(i,j,1,2,2,dsr,Vdsr,Ndsr,hstep);
      if j~=1,
        worka=worka+disco_stiff(i,j,2,1,1,dsr,Vdsr,Ndsr,hstep);
      end
        if worka~=0,
           cola=row-nanr;
           A(row,cola)=worka;
        end
   end
%________________________________________________
    if (i-1~=0) & (j~=nanr)
       worka=disco_stiff(i,j,3,2,2,dsr,Vdsr,Ndsr,hstep) ...
            +disco_stiff(i,j,2,3,3,dsr,Vdsr,Ndsr,hstep);
        if worka~=0,
           cola=row-nanr+1;
           A(row,cola)=worka;
        end
   end
%________________________________________________
   if j-1~=0,
        worka=disco_stiff(i,j,1,3,4,dsr,Vdsr,Ndsr,hstep);
        if i~=1,
           worka=worka+disco_stiff(i,j,3,1,1,dsr,Vdsr,Ndsr,hstep);
	end
        if worka~=0,
           cola=row-1;
           A(row,cola)=worka;
        end
   end
%________________________________________________ the central node
	cola=row;
if (i~=1)&(j~=1),
A(row,cola)= disco_stiff(i,j,1,1,1,dsr,Vdsr,Ndsr,hstep) ...
           + disco_stiff(i,j,2,2,2,dsr,Vdsr,Ndsr,hstep) ...
           + disco_stiff(i,j,3,3,3,dsr,Vdsr,Ndsr,hstep) ...
           + disco_stiff(i,j,1,1,6,dsr,Vdsr,Ndsr,hstep) ...
           + disco_stiff(i,j,3,3,4,dsr,Vdsr,Ndsr,hstep) ...
           + disco_stiff(i,j,2,2,5,dsr,Vdsr,Ndsr,hstep);
end
if (i~=1)&(j==1),
A(row,cola) = disco_stiff(i,j,1,1,6,dsr,Vdsr,Ndsr,hstep) ...
            + disco_stiff(i,j,2,2,2,dsr,Vdsr,Ndsr,hstep) ...
            + disco_stiff(i,j,3,3,3,dsr,Vdsr,Ndsr,hstep);
end
if (i==1)&(j~=1),
A(row,cola) = disco_stiff(i,j,1,1,6,dsr,Vdsr,Ndsr,hstep) ...
            + disco_stiff(i,j,2,2,5,dsr,Vdsr,Ndsr,hstep) ...
            + disco_stiff(i,j,3,3,4,dsr,Vdsr,Ndsr,hstep);
end
if (i==1)&(j==1),
  A(row,cola) = disco_stiff(i,j,1,1,6,dsr,Vdsr,Ndsr,hstep);
end

%________________________________________________
   if j~=nanr,
        worka=disco_stiff(i,j,3,1,6,dsr,Vdsr,Ndsr,hstep);
        if i~=1,
           worka=worka+disco_stiff(i,j,1,3,3,dsr,Vdsr,Ndsr,hstep);
	end
        if worka~=0,
           cola=row+1;
           A(row,cola)=worka;
        end
   end
%________________________________________________
    if (i~=nanr) & (j~=1) ,
       worka = disco_stiff(i,j,2,3,4,dsr,Vdsr,Ndsr,hstep) ...
             + disco_stiff(i,j,3,2,5,dsr,Vdsr,Ndsr,hstep);
        if worka~=0,
           cola=row+nanr-1;
           A(row,cola)=worka;
        end
   end
%________________________________________________
   if i~=nanr,
        worka=disco_stiff(i,j,2,1,6,dsr,Vdsr,Ndsr,hstep);
        if j~=1,
           worka=worka+disco_stiff(i,j,1,2,5,dsr,Vdsr,Ndsr,hstep);
	end
        if worka~=0,
           cola=row+nanr;
           A(row,cola)=worka;
        end		
   end
%________________________________________________
end
end

return