clear
% 2D
syms x y
u = x^2-y^2; %x^3-3*x*y^2                 % "manufactured solution"
dux=diff(u,'x')
duy=diff(u,'y')
diff(dux,'x')+diff(duy,'y')
% generacja siatki punktow w kwadracie

h=.1;
[X,Y]=meshgrid(0:h:1,0:h:1);
U = double(subs(u,{x,y},{X,Y}));
figure(1), surf(X,Y,U), axis equal, colorbar

% boundary elements
[M,M]=size(X);
m=M-1;
Nelb=0;

ii=1;
for jj=1:m
    Nelb=Nelb+1;
    xk=h/2+(jj-1)*h;
    ieltop(Nelb,1:4)=[ii,jj,ii,jj+1];
    celnel(Nelb,1:4)=[xk, 0, 0,  -1];
    ub(Nelb)=subs(u,{x,y},{xk,0});
    tb(Nelb)=subs(dux,{x,y},{xk,0})*0+subs(duy,{x,y},{xk,0})*(-1);
end

jj=M;
for ii=1:m
    Nelb=Nelb+1;
    yk=h/2+(ii-1)*h;
    ieltop(Nelb,1:4)=[ii,jj,ii+1,jj];
    celnel(Nelb,1:4)=[1, yk, 1, 0];
    ub(Nelb)=subs(u,{x,y},{1,yk});
    tb(Nelb)=subs(dux,{x,y},{1,yk})*1+subs(duy,{x,y},{1,yk})*0;
end

ii=M;
for jj=M:(-1):2
    Nelb=Nelb+1;
    xk=h/2+(jj-2)*h;
    ieltop(Nelb,1:4)=[ii,jj,ii,jj-1];
    celnel(Nelb,1:4)=[xk, 1, 0, 1];
    ub(Nelb)=subs(u,{x,y},{xk,1});
    tb(Nelb)=subs(dux,{x,y},{xk,1})*0+subs(duy,{x,y},{xk,1})*1;
end

jj=1;
for ii=M:(-1):2
    Nelb=Nelb+1;
    yk=h/2+(ii-2)*h;
    ieltop(Nelb,1:4)=[ii,jj,ii-1,jj];
    celnel(Nelb,1:4)=[0,yk , -1, 0];
    ub(Nelb)=subs(u,{x,y},{0,yk});
    tb(Nelb)=subs(dux,{x,y},{0,yk})*(-1)+subs(duy,{x,y},{0,yk})*0;
end
figure(2), plot(ub), hold on
figure(3), plot(tb), hold on

% boundary values approximation
for iel=1:Nelb
    xx=celnel(iel,1); yy=celnel(iel,2);
    xn=celnel(iel,3); yn=celnel(iel,4);
    ub(iel)=subs(u,{x,y},{xx,yy});
    qb(iel)=subs(dux,{x,y},{xx,yy})*xn+subs(duy,{x,y},{xx,yy})*yn;
end

% algebraic equations
for ii=1:Nelb
    xk=celnel(ii,1);
    yk=celnel(ii,2);
    for jj=1:Nelb
        [vgitg,tgitg]=elemb(ii,jj,xk,yk,h,celnel,ieltop,X,Y);
        G(ii,jj)=vgitg;
        H(ii,jj)=tgitg;
    end
    H(ii,ii)=-sum(H(ii,:));
end
TB=G\(H*ub');
plot(TB,'r')

% postprocessing along the boundary
Uh=U*0;
for iel=1:Nelb
    ii=ieltop(iel,1);
    jj=ieltop(iel,2);
    Uh(ii,jj)=Uh(ii,jj)+ub(iel)/2;
    ii=ieltop(iel,3);
    jj=ieltop(iel,4);
    Uh(ii,jj)=Uh(ii,jj)+ub(iel)/2;
end

% postprocessing in interior
for ii=2:M-1
    for jj=2:M-1
        x0=X(ii,jj);
        y0=Y(ii,jj);
        for kk=1:Nelb
            [vgitg,tgitg]=elemb(0,kk,x0,y0,h,celnel,ieltop,X,Y);
            GG(kk)=vgitg;
            HH(kk)=tgitg;
        end
        Uh(ii,jj)=GG*tb'-HH*ub';
    end
end
figure(4), surf(X,Y,Uh), axis equal, colorbar

% ASSIGNMENT
% Compute exact an dapproximate solutions at the interior point
x0=3/4;
y0=1/4;

