clear all
close all

N=50;
meshfile=fc_oogmsh.gmsh.buildmesh2d('square4',N);
Th=fc_simesh.siMesh(meshfile,'trans',@(q) [-pi/3+q(1,:);q(2,:)]);

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);
du={@(x,y) -pi*sin(pi*(x-y)),@(x,y) pi*sin(pi*(x-y))}; % du{1}: fct du/dx, du{2}: fct du/dy,
Rlab=1:4; % Neumann labels
gR{1}=@(x,y) -du{1}(x,y);
gR{2}=@(x,y)  du{1}(x,y);
gR{3}=@(x,y) -du{2}(x,y);
gR{4}=@(x,y)  du{2}(x,y);

Uex=siMesh.eval(Th,u);

% Calcul de Ah et b
K=siMesh.Rigidite(Th);
M=siMesh.Masse(Th);
Ah=eta*M+K;
F=siMesh.eval(Th,f);
b=M*F;
for lab=Rlab
  b=b+siMesh.Masse(Th,'d',1,'labels',lab)*siMesh.eval(Th,gR{lab},'d',1,'labels',lab);
end
% Resolution
U=FEM.solveBVPsystem(Ah,b,[]);

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 


