% Matrice de Masse locale et globale sur Omega
close all
clear all

meshfile=fc_oogmsh.gmsh.buildmesh2d('./geofile/square6dom4holes.geo',5);
Th=fc_simesh.siMesh(meshfile); % maillage global

C1=1;C2=2;
u1=@(x,y) C1+0*x; % fonction u sur Omega{10}
u2=@(x,y) C2+0*x; % fonction u sur Omega{2:2:8,20}
v=@(x,y) sin(x-y+1);
% integrale de u.v sur [a,b]x[c,d] avec Sage
Iexfun=@(a,b,c,d) -sin(-a+c-1) + sin(-a+d-1) + sin(-b+c-1) - sin(-b+d-1);
IexfunEh=@(Eh) Iexfun(Eh.bbox(1),Eh.bbox(2),Eh.bbox(3),Eh.bbox(4));

idx=Th.find(2);

%%% (1)
%%%   a. calcul integrale de u.v sur Omega(2)+Omega(20)
labels=[2,20];
Iex=0; 
for lab=labels % calcul exact : u.v = C1.v
  i=Th.find(2,lab);
  Iex=Iex+C1*IexfunEh(Th.sTh{i});
end

IL=0; % Calculs locaux
u=@(x,y) u1(x,y);
for lab=labels
  i=Th.find(2,lab);
  if lab==10,u=@(x,y) u1(x,y);else,u=@(x,y) u2(x,y);end
  ML=siMeshElt.Masse(Th.sTh{i}); % Matrice locale
  UL=Th.sTh{i}.eval(u);
  VL=Th.sTh{i}.eval(v);
  IL=IL+dot(ML*UL,VL);
end
fprintf('(1) a.  Iex=%.5e, |IL-Iex|=%.6e \n',Iex,abs(Iex-IL))

% Calculs globaux (methode 1)
IG=0; 
for lab=labels
  i=Th.find(2,lab);
  if lab==10,u=@(x,y) u1(x,y);else,u=@(x,y) u2(x,y);end
  MG=siMeshElt.Masse(Th.sTh{i},'local',false); % Matrice globale
  UG=zeros(Th.nq,1);UG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(u);
  VG=zeros(Th.nq,1);VG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(v);
  IG=IG+dot(MG*UG,VG);
end
fprintf('        |IL-IG|=%.6e  [IG, methode 1]\n',abs(IG-IL))

% Calculs globaux (methode 3)  u continue sur Omega(2)+Omega(20) 
MG=sparse(Th.nq,Th.nq);
UG=zeros(Th.nq,1);VG=zeros(Th.nq,1);
for lab=labels
  i=Th.find(2,lab);
  MG=MG+siMeshElt.Masse(Th.sTh{i},'local',false);
  if lab==10,u=@(x,y) u1(x,y);else,u=@(x,y) u2(x,y);end
  UG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(u);
  VG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(v);
end
IG=dot(MG*UG,VG);
fprintf('        |IL-IG|=%.6e  [IG, methode 3]\n',abs(IG-IL))

% Calculs globaux (methode 4)  u continue sur Omega(2)+Omega(20) 
IG=0;
for lab=labels
  i=Th.find(2,lab);
  MG=siMeshElt.Masse(Th.sTh{i},'local',false);
  if lab==10,u=@(x,y) u1(x,y);else,u=@(x,y) u2(x,y);end
  UG=zeros(Th.nq,1);VG=zeros(Th.nq,1);
  UG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(u);
  VG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(v);
  IG=IG+dot(MG*UG,VG);
end
fprintf('        |IL-IG|=%.6e  [IG, methode 4]\n',abs(IG-IL))


%%%   b. calcul integrale de u.v sur Omega(10+20)
labels=[10,20];
i10=Th.find(2,10);i20=Th.find(2,20);
Iex=C1*IexfunEh(Th.sTh{i10}); % u.v=C1.v sur Omega(10)
for lab=[1:8,20] % frontieres pour enlever les int. sur les carrees interieurs
  i=Th.find(1,lab);
  Iex=Iex-C1*IexfunEh(Th.sTh{i});  
end
Iex=Iex+C2*IexfunEh(Th.sTh{i20});  % u.v=C2.v sur Omega(20)

IL=0; % Calculs locaux
for lab=labels
  i=Th.find(2,lab);
  ML=siMeshElt.Masse(Th.sTh{i}); % Matrice locale
  if lab==10,u=@(x,y) u1(x,y);else,u=@(x,y) u2(x,y);end
  UL=Th.sTh{i}.eval(u);
  VL=Th.sTh{i}.eval(v);
  IL=IL+dot(ML*UL,VL);
end
fprintf('(1) b.  Iex=%.5e, |IL-Iex|=%.6e \n',Iex,abs(Iex-IL))

% Calculs globaux (methode 1)
IG=0; 
for lab=labels
  i=Th.find(2,lab);
  MG=siMeshElt.Masse(Th.sTh{i},'local',false); % Matrice globale
  if lab==10,u=@(x,y) u1(x,y);else,u=@(x,y) u2(x,y);end
  UG=zeros(Th.nq,1);UG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(u);
  VG=zeros(Th.nq,1);VG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(v);
  IG=IG+dot(MG*UG,VG);
end
fprintf('        |IL-IG|=%.6e  [IG, methode 1]\n',abs(IG-IL))

% Calculs globaux (methode 3) 
MG=sparse(Th.nq,Th.nq);
UG=zeros(Th.nq,1);VG=zeros(Th.nq,1);
for lab=labels
  i=Th.find(2,lab);
  MG=MG+siMeshElt.Masse(Th.sTh{i},'local',false);
  if lab==10,u=@(x,y) u1(x,y);else,u=@(x,y) u2(x,y);end
  UG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(u);
  VG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(v);
end
IG=dot(MG*UG,VG);
fprintf('        |IL-IG|=%.6e  [IG, methode 3]\n',abs(IG-IL))

% Calculs globaux (methode 4)  
IG=0;
for lab=labels
  i=Th.find(2,lab);
  MG=siMeshElt.Masse(Th.sTh{i},'local',false);
  if lab==10,u=@(x,y) u1(x,y);else,u=@(x,y) u2(x,y);end
  UG=zeros(Th.nq,1);VG=zeros(Th.nq,1);
  UG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(u);
  VG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(v);
  IG=IG+dot(MG*UG,VG);
end
fprintf('        |IL-IG|=%.6e  [IG, methode 4]\n',abs(IG-IL))


%%%   c. calcul integrale de u.v sur Omega
%  u=u1 sur Omega(10), u=u2 sinon  (donc u non continue sur Omega)

% calcul exact (par rapport aux bbox des bords)
i10=Th.find(2,10);
Iex=C1*IexfunEh(Th.sTh{i10}); % u.v=C1.v sur Omega(10)
for lab=[1:8,20] % label des frontieres (1-simplex), pour enlever les int. sur les carrees interieurs
  i=Th.find(1,lab);
  Iex=Iex-C1*IexfunEh(Th.sTh{i});  
end
for lab=[2:2:8,20] % label des domaines (2-simplex)
  i=Th.find(2,lab);
  Iex=Iex+C2*IexfunEh(Th.sTh{i});  % u.v=C2.v sur Omega(lab)
end

[idx,labels]=Th.find(2);
IL=0; % Calculs locaux
for i=idx
  ML=siMeshElt.Masse(Th.sTh{i}); % Matrice locale
  if Th.sTh{i}.label==10,u=@(x,y) u1(x,y);else,u=@(x,y) u2(x,y);end
  UL=Th.sTh{i}.eval(u);
  VL=Th.sTh{i}.eval(v);
  IL=IL+dot(ML*UL,VL);
end
fprintf('(1) c.  Iex=%.5e, |IL-Iex|=%.6e \n',Iex,abs(Iex-IL))

% Calculs globaux (methode 1)
IG=0; 
for i=idx
  MG=siMeshElt.Masse(Th.sTh{i},'local',false); % Matrice globale
  if Th.sTh{i}.label==10,u=@(x,y) u1(x,y);else,u=@(x,y) u2(x,y);end
  UG=zeros(Th.nq,1);UG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(u);
  VG=zeros(Th.nq,1);VG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(v);
  IG=IG+dot(MG*UG,VG);
end
fprintf('        |IL-IG|=%.6e  [IG, methode 1]\n',abs(IG-IL))

% Calculs globaux (methode 3) hyp. u non continue sur Omega et v continue sur Omega
%   u=u1 sur Omega(10), u=u2 sinon
MG=sparse(Th.nq,Th.nq);
UG=zeros(Th.nq,1);VG=zeros(Th.nq,1);
for lab=labels
  i=Th.find(2,lab);
  MG=MG+siMeshElt.Masse(Th.sTh{i},'local',false);
  if Th.sTh{i}.label==10,u=@(x,y) u1(x,y);else,u=@(x,y) u2(x,y);end
  UG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(u);
  VG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(v);
end
IG=dot(MG*UG,VG);
fprintf('        |IL-IG|=%.6e  [IG, methode 3]\n',abs(IG-IL))

