function webassemb
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% assembles PDE for 2D using WEB-method
% n                  : degree
% xx,yy, x,y(gridded): position of grid cells
% D                  : all b-splines supporting the region 
% bb                 : inner(2) and outer(1) b-splines
% gi,gj              : inner&outer b-splines
% TH                 : the direction angle of wave
% k2                 : square of wave number
% K,M,F,DNq,DNg      : Matrices for assembly
% Created by G.Apaydin in June 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global name n bb h D C bound KM col row K M F DNq DNg TH k2 NC step DNNq
%----------------------------INITILIAZE MATRIX-----------------------------
ld=length(D);F=sparse(ld,1);%K=zeros(ld);M=(ld,ld);
ik=reshape(1:ld,size(bb));ik=ik(1:end-n,1:end-n);ik=ik(:);
row=zeros(length(ik),(n+1)^4);col=row;KM=row;
lboun=cell(1,NC);DNg=cell(1,NC);DNq=cell(1,NC);BNq=cell(1,NC);BNg=cell(1,NC);
Nx=cell(1,NC);Ny=cell(1,NC);Nxx=cell(1,NC);Nyy=cell(1,NC);w=cell(1,NC);
bNxy1=cell(1,NC);bNxy2=cell(1,NC);
for i=1:NC
    if bound{i}.con==1
        DNg{i}=sparse(ld,1);DNq{i}=KM;
        BNq{i}=sqrt(-k2);BNg{i}=1;
    end    
end
%----------------------------THE COEFFICIENTS------------------------------
MM=-k2;FK=0;
%--------------------------------------------------------------------------

[Kinti,Minti,Finti,BSPL,DBSPLX,DBSPLY,BXYval,DBXval,DBYval,BXYrval,DBXrval,DBYrval,wwr]=inn_out_integ(MM,FK);
%assembly
[mm,nn]=size(bb);
for ii=1:length(ik)
    k=ik(ii);
   
    if bb(k)>0,
        Dx1=D(k,3)+n*h;Dx2=Dx1+h;Dy1=D(k,4)+n*h;Dy2=Dy1+h;
        cc=isinpoly([Dx1 Dx2 Dx1 Dx2],[Dy1 Dy1 Dy2 Dy2],C(1,:),C(2,:));
        if NC>1,for tt=2:NC,cc=cc.*(1-isinpoly([Dx1 Dx2 Dx1 Dx2],[Dy1 Dy1 Dy2 Dy2],C(2*tt-1,:),C(2*tt,:)));end,end
       
        kkk=[];
        if max(cc)>0,
            for q=0:n,kkk=[kkk k+q*mm+(0:n)];end,
            
            row(ii,:)=kron(kkk,ones(1,(n+1)^2));col(ii,:)=kron(ones(1,(n+1)^2),kkk);
            con2=(bb(row(ii,:))~=0);con1=(con2.*bb(col(ii,:)))~=0;
            if sum(cc)==4,
                KM(ii,:)=con1.*(Kinti+Minti);
                F(row(ii,:))=F(row(ii,:))+(con2.*Finti).';
            else    
                step=51;
                for i=1:NC,
                    if bound{i}.cont==1,
                        [Nx{i},Ny{i},w{i}]=bounellipse(Dx1,Dx2,Dy1,Dy2,bound{i});lboun{i}=1;
                        Nxx{i}=(Nx{i}-Dx1)/h;Nyy{i}=(Ny{i}-Dy1)/h;
                        [X,Y]=meshgrid(linspace(Dx1,Dx2,step),linspace(Dy1,Dy2,step));cond=bound{i}.fun(X,Y)>=0;
                    elseif bound{i}.cont==0,
                        xmin=max(min(C(1,:)),Dx1);xmax=min(max(C(1,:)),Dx2);
                        ymin=max(min(C(2,:)),Dy1);ymax=min(max(C(2,:)),Dy2);
                        [Nx{i},Ny{i},w{i}]=bounrect(Dx1,Dx2,Dy1,Dy2,bound{i});lboun{i}=1;
                        Nxx{i}=(Nx{i}-Dx1)/h;Nyy{i}=(Ny{i}-Dy1)/h;
                    elseif bound{i}.cont==2,
                        [Nx{i},Ny{i},w{i}]=boununarbit(Dx1,Dx2,Dy1,Dy2,bound{i});lboun{i}=1;
                        Nxx{i}=(Nx{i}-Dx1)/h;Nyy{i}=(Ny{i}-Dy1)/h;
                    end
                    NSUMq{i}=zeros((n+1)^2);
                end
                KMint=zeros((n+1)^2);
                for q=1:(n+1)^2
                    kkkk=kkk(q);
                    Fint=0;NSUM=0;
                    if bb(kkkk)>0,
                        BXY1=BSPL{q};DBX1=DBSPLX{q};DBY1=DBSPLY{q};
                        BXY1val=BXYval{q};DBX1val=DBXval{q};DBY1val=DBYval{q};
                        BXY1rval=BXYrval{q};DBX1rval=DBXrval{q};DBY1rval=DBYrval{q};
                        for i=1:NC,
                            if lboun{i}>0
                                bNxy1{i}=BXY1(Nxx{i},Nyy{i});
                            
                                if BNg{i}~=0,
                                    
                                    if bound{i}.cont==1,
                                        TH=0;NSUM=1/h*BNg{i}*BNq{i}*sum(w{i}.*bNxy1{i}.*exp(BNq{i}*Nx{i}).*(Nx{i}/1+1));
                                    
                                    elseif bound{i}.cont==0,
                                        NSUM=1/h*BNg{i}*BNq{i}*sum(w{i}.*bNxy1{i}.*exp(BNq{i}*(cos(TH)*Nx{i}+sin(TH)*Ny{i})).*(...
                                            ((Nx{i}>=xmin)&(Nx{i}<=xmax)&(Ny{i}==ymin)).*(1-sin(TH))+...
                                            ((Nx{i}==xmax)&(Ny{i}<ymax)&(Ny{i}>ymin)).*(1+cos(TH))+...
                                            ((Nx{i}>=xmin)&(Nx{i}<=xmax)&(Ny{i}==ymax)).*(1+sin(TH))+...
                                            ((Nx{i}==xmin)&(Ny{i}<ymax)&(Ny{i}>ymin)).*(1-cos(TH))    ));
                                       
                                    end
                                    DNg{i}(kkkk)=DNg{i}(kkkk)+NSUM;
                                end
                            
                            end
                        
                        end
                        for r=q:(n+1)^2,
                            llll=kkk(r);
                            if bb(llll)>0,
                                BXY2=BSPL{r};DBX2=DBSPLX{r};DBY2=DBSPLY{r};
                                BXY2val=BXYval{r};DBX2val=DBXval{r};DBY2val=DBYval{r};
                                BXY2rval=BXYrval{r};DBX2rval=DBXrval{r};DBY2rval=DBYrval{r};
                                
                                for i=1:NC,
                                    if lboun{i}>0&&BNq{i}~=0,
                                        bNxy2{i}=BXY2(Nxx{i},Nyy{i});
                                        %integral(BN1q u v) for the curve
                                        NSUMq{i}(q,r)=BNq{i}/h/h*sum(w{i}.*bNxy1{i}.*bNxy2{i});
                                        NSUMq{i}(r,q)=NSUMq{i}(q,r);
                                    end
                                
                                    if bound{i}.cont==1,
                                        [KMint(q,r),Fint]=outsplnintegellipse(cond,BXY1val,BXY2val,DBX1val,DBX2val,DBY1val,DBY2val,MM,FK);
                                        KMint(r,q)=KMint(q,r);
                                    elseif bound{i}.cont==0,
                                        [KMint(q,r),Fint]=outsplnintegrect(xmin,xmax,ymin,ymax,BXY1rval,BXY2rval,DBX1rval,DBX2rval,DBY1rval,DBY2rval,MM,FK,wwr);
                                        KMint(r,q)=KMint(q,r); 
                                    elseif bound{i}.cont==2,
                                        %[Nx{i},Ny{i},w{i}]=boununarbit(Dx1,Dx2,Dy1,Dy2,bound{i});lboun{i}=1;
                                        KMint(p,r)=0;Fint=0;KMint(r,q)=KMint(q,r); 
                                    end
                                end
                            
                            end
                        end
                        F(kkkk)=F(kkkk)+Fint;
                    end                
                end
                KMint=KMint(:);
                KM(ii,:)=KMint.';
                for i=1:NC,
                    if lboun{i}>0
                        DNq{i}(ii,:)=NSUMq{i}(:);
                    end
                end
            end
        else %display('Empty grid cell');
        end
    end
end
index=find(sum(row,2)==0);row(index,:)=[];col(index,:)=[];KM(index,:)=[];
K=sparse(row,col,KM,ld,ld);
for i=1:NC,
    if lboun{i}>0
        DNq{i}(index,:)=[];DNNq{i}=sparse(row,col,DNq{i},ld,ld);
    end
end