clear all
close all

N=50;
meshfile=fc_oogmsh.gmsh.buildmesh2d('square4',N);
Th=fc_simesh.siMesh(meshfile);

u=@(x,y) cos(pi*(x-y));
eta=1/2;
f=@(x,y) 2*pi^2*cos(pi*(x - y))+eta*u(x,y);
gD=@(x,y) u(x,y);

Uex=siMesh.eval(Th,u);

% Calcul de Ah et b (sans conditions de Dirichlet)
F=siMesh.eval(Th,f);
K=siMesh.Rigidite(Th);
M=siMesh.Masse(Th);
Ah=eta*M+K;
b=M*F;

Dlab=1:4; % Dirichlet labels
% Calcul de ID, indices globaux bords Dirichlet 
ID=siMesh.get_toGlobal(Th,'d',1,'labels',Dlab);
% Relevement
R=siMesh.eval(Th,gD,'d',1,'labels',Dlab);
% Resolution
U=FEM.solveBVPsystem(Ah,b,ID,R,2);

figure(1)
Th.plot(Uex)
shading interp
colorbar
axis off;axis image 

figure(2)
Th.plot(U)
shading interp
colorbar
axis off;axis image 

figure(3)
colormap jet
Th.plot(abs(U-Uex))
shading interp
colorbar
axis off;axis image 
