%% SCRIPTS TO TEST THE p1-FEM

f  = @(x,y) 0*x;
uD = @(x,y) x+2*y;
uN = @(x,y) [1+0*x 2+0*y];
ux = @(x,y) 1+0*x;
uy = @(x,y) 2+0*y;

load fourTriang
Tcollect={myT1 myT2 myT3 myT4};
for j=1:length(Tcollect)
    T=Tcollect{j};
    [u,S,M]=P1FEM(T,f,uD,uN,0);
    x=T.coordinates(:,1); y=T.coordinates(:,2);
    [errorL2(j),errorH1(j)]=computeError(T,u,uD,ux,uy);
    subplot(2,2,j)
    trisurf(T.elements,x,y,u), shading flat
end

%% 

uD = @(x,y) y.*x.^2+3*y.^2.*x;
c=0;
f  = @(x,y) -2*y-6*x+c*uD(x,y);
uN = @(x,y) [2*x.*y+3*y.^2  x.^2+6*x.*y];
ux = @(x,y) 2*x.*y+3*y.^2;
uy = @(x,y) x.^2+6*x.*y;

load fourTriang
Tcollect={myT1 myT2 myT3 myT4};
for j=1:length(Tcollect)
    T=Tcollect{j};
    [u,S,M]=P1FEM(T,f,uD,uN,c);
    [errorL2(j),errorH1(j)]=computeError(T,u,uD,ux,uy);
end

h=[1 1/2 1/4 1/8];
loglog(h,errorL2,'-o'), hold on
loglog(h,errorH1,'-or')

%%

uD = @(x,y) sin(4*x-y);
c=1;
f  = @(x,y) 17*sin(4*x-y)+c*uD(x,y);
uN = @(x,y) [4*cos(4*x-y) -cos(4*x-y)];
ux = @(x,y) 4*cos(4*x-y);
uy = @(x,y) -cos(4*x-y);

load fourTriang
Tcollect={myT1 myT2 myT3 myT4};
for j=1:length(Tcollect)
    T=Tcollect{j};
    [u,S,M]=P1FEM(T,f,uD,uN,c);
    [errorL2(j),errorH1(j)]=computeError(T,u,uD,ux,uy);
end

h=[1 1/2 1/4 1/8];
loglog(h,errorL2,'-o'), hold on
loglog(h,errorH1,'-or')

%%