function webassemb
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% assemble PDE for 1D 
% n             : degree
% xx            : 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
% k             : square of wave number
% K,M,F,DNq,DNg : Matrices for assembly
%
% Created by G.Apaydin in June 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global D bb n bound xi xe h K F weight weightt NP px
%-------------------------INITILIAZE MATRICES------------------------------
ld=length(D);F=sparse(ld,1);
ik=reshape(1:ld,size(bb));ik=ik(1:end-n);ik=ik(:);
row=zeros(length(ik),(n+1)^2);K=row;M=row;
for i=1:length(bound.BC)
    if bound.BC{i}==1,DNg{i}=sparse(ld,1);DNq{i}=K;NSUMq{i}=zeros(n+1);end
end
%----------------B-SPLINES AND THEIR DERIVATIVE FUNCTIONS------------------
[BSPL,DBSPL]=bsplfunc;
%------------------------------ASSEMBLY------------------------------------                   
for ii=1:length(ik)
    k=ik(ii);
    
    if bb(k)>0
        Dx1=D(k,2)+n*h;Dx2=Dx1+h; %grid coordinates
        if (Dx1>=xi&&Dx1<xe)||(Dx2>xi&&Dx2<=xe)
            kkk=k+[0:n];
            row(ii,:)=kron(kkk,ones(1,n+1));
            col(ii,:)=kron(ones(1,n+1),kkk);
            if Dx1<=xi,[X,W]=gnodewt(NP,xi,Dx2);boun=1;
            elseif Dx2>=xe,[X,W]=gnodewt(NP,Dx1,xe);boun=1;
            else [X,W]=gnodewt(NP,Dx1,Dx2);boun=0;
            end
            XX=(X-Dx1)/h;          
            
            ww=weight(X);dww=weightt(X);
           
            Kint=zeros(n+1);Mint=zeros(n+1);
            for q=1:(n+1)  
 
                kkkk=kkk(q);
                
                Fint=0;
                if bb(kkkk)>0,
                    BX1=BSPL{q};DB1=DBSPL{q};
                    %curve for function v
                    if boun==1
                        for i=1:length(bound.BC)
                            if bound.BC{i}==1&&Dx1<=bound.x{i}&&bound.x{i}<=Dx2
                                %NEUMANN CONDITION
                                Nx{i}=bound.x{i};Nxx{i}=(bound.x{i}-Dx1)/h;
                                bNx1{i}=weight(Nx{i})*BX1(Nxx{i})/sqrt(h);
                                NSUM=bound.BNg{i}*bNx1{i}*px(Nx{i});
                                DNg{i}(kkkk)=DNg{i}(kkkk)+NSUM;
                            end
                        end
                    end

                    for r=q:(n+1),
                        llll=kkk(r);

                       
                
                        if bb(llll)>0,
                            BX2=BSPL{r};DB2=DBSPL{r};
                            [Kint(q,r),Mint(q,r),Fint]=splninteg(ww.*BX1(XX)/sqrt(h),ww.*BX2(XX)/sqrt(h),dww.*BX1(XX)/sqrt(h)+ww.*DB1(XX)/h/sqrt(h),dww.*BX2(XX)/sqrt(h)+ww.*DB2(XX)/h/sqrt(h),W,X);
                            Kint(r,q)=Kint(q,r).';Mint(r,q)=Mint(q,r).';
                            
                            if boun==1
                                for i=1:length(bound.BC)
                                    if bound.BC{i}==1&&Dx1<=bound.x{i}&&bound.x{i}<=Dx2
                                        %NEUMANN integral(BN1q u v) for the curve
                                        bNx2{i}=weight(Nx{i})*BX2(Nxx{i})/sqrt(h);
                                        NSUMq{i}(q,r)=bound.BNq{i}*bNx1{i}*bNx2{i}*px(Nx{i});
                                        NSUMq{i}(r,q)=NSUMq{i}(q,r);                                        
                                    end
                                end
                            end
                        end
                    end
                    F(kkkk)=F(kkkk)+Fint;
                end
            end
            Kint=Kint(:);Mint=Mint(:);
            K(ii,:)=Kint.';M(ii,:)=Mint.';
            for i=1:length(bound.BC)
                if bound.BC{i}==1
                    DNq{i}(ii,:)=NSUMq{i}(:);
                    NSUMq{i}=zeros(n+1);
                end
            end
            
        else
            display('Disarida');
            
        end
    end
end
index=find(sum(row,2)==0);row(index,:)=[];col(index,:)=[];
K(index,:)=[];M(index,:)=[];
K=sparse(row,col,K,ld,ld);M=sparse(row,col,M,ld,ld);

K=K+M;
for i=1:length(bound.BC)
    if bound.BC{i}==1
        DNq{i}(index,:)=[];DNNq{i}=sparse(row,col,DNq{i},ld,ld);
        F=F+DNg{i};K=K+DNNq{i};
    end
end