%  clear all
%  close all

N=30;
meshfile=fc_oogmsh.gmsh.buildmesh2d('condenser11',N);
Th=fc_simesh.siMesh(meshfile);

% Calcul de Ah et b (sans conditions de Dirichlet)
Ah=siMesh.Rigidite(Th);
b=zeros(Th.nq,1);
% Prise en compte des conditions de Dirichlet (comme inconnues, methode 1)
Dlab=[1:2:7,101,102,103,104]; % Dirichlet labels
ID=siMesh.get_toGlobal(Th,'d',1,'labels',Dlab);
%IDc=setdiff(1:Th.nq,ID);

% Les bords 1,3,5 et 7 sont disjoints
% et 0 sur les autres bords (donc continue)
% Relevement
R=siMesh.eval(Th,+12,'d',1,'labels',[1,5])+siMesh.eval(Th,-12,'d',1,'labels',[3,7]);

MySolvefun=@(A,b) gmres(A,b,20,1e-8);
[U,flag,relres,iter]=FEM.solveBVPsystemV2(Ah,b,ID,R,'method',2,'solvefun',MySolvefun);

% Using pcg function:
% MySolvefun=@(A,b) pcg(A,b,1e-12,1000);
% [U,flag,relres,iter]=FEM.solveBVPsystemV2(Ah,b,ID,R,'method',2,'solvefun',MySolvefun);

% Using bicg function:
% U=FEM.solveBVPsystemV2(Ah,b,ID,R,'method',2,'solvefun',@(A,b) bicg(A,b,1e-8,300));

Th.plotmeshinfo()

figure()
Th.plot(U)
shading interp
colorbar
axis off;axis image
hold on
Th.plotiso(U,'color','w')
