% 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]=fe_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;
Coor=hstep*[0:nanr-1];

%A=spalloc(nall,nall,7*nall);
disp('Begin allocating memory...')
lengthA=nanr*nanr*7;
AI  = zeros(lengthA,1);
AJ  = zeros(lengthA,1);
AV  = zeros(lengthA,1);
nextA = 0;
disp('...end allocating memory.')%

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

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

AI = AI(1:nextA);
AJ = AJ(1:nextA);
AV = AV(1:nextA);

A  = sparse(AI, AJ, AV); 

return