% Matrice de Masse locale et globale sur Omega
% en utilisant siMesh.Masse
close all
clear all

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


u=@(x,y) cos(x+y-pi/3);
v=@(x,y) sin(x-y+1);
% integrale de u.v sur [a,b]x[c,d] avec Sage
Iexfun=@(a,b,c,d) -c*(cos(1-pi/3+2*a)-cos(1-pi/3+2*b))/4 + d*(cos(1-pi/3+2*a)-cos(1-pi/3+2*b))/4 + (a-b)*cos(2*c-1-pi/3)/4 - (a-b)*cos(2*d-1-pi/3)/4;
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(20)
lab=20;
ML=siMesh.Masse(Th,'labels',lab,'local',true);
UL=siMesh.eval(Th,u,'labels',lab,'local',true);
VL=siMesh.eval(Th,v,'labels',lab,'local',true);
IL=dot(ML*UL,VL);
i=Th.find(2,lab);
Iex=IexfunEh(Th.sTh{i});
fprintf('(1) a.  Iex=%.5e, |IL-Iex|=%.6e \n',Iex,abs(Iex-IL))

MG=siMesh.Masse(Th,'labels',lab,'local',false); % Matrice globale
UG=siMesh.eval(Th,u,'labels',lab,'local',false);
VG=siMesh.eval(Th,v,'labels',lab,'local',false);
IG=dot(MG*UG,VG);

fprintf('        |IL-IG|=%.6e\n',abs(IL-IL))

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


ML=siMesh.Masse(Th,'labels',labels,'local',true);
UL=siMesh.eval(Th,u,'labels',labels,'local',true);
VL=siMesh.eval(Th,v,'labels',labels,'local',true);
IL=0; % Calculs locaux
for i=1:length(labels)
  IL=IL+dot(ML{i}*UL{i},VL{i});
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=siMesh.Masse(Th,'labels',lab,'local',false);
  UG=siMesh.eval(Th,u,'labels',lab,'local',false);
  VG=siMesh.eval(Th,v,'labels',lab,'local',false);
  IG=IG+dot(MG*UG,VG);
end
fprintf('        |IL-IG|=%.6e  [IG, methode 1]\n',abs(IG-IL))

% Calculs globaux (methode 2) hyp. u et v continues sur Omega
MG=siMesh.Masse(Th,'labels',labels,'local',false);
UG=Th.eval(u);
VG=Th.eval(v);
IG=dot(MG*UG,VG);
fprintf('        |IL-IG|=%.6e  [IG, methode 2]\n',abs(IG-IL))

% Calculs globaux (methode 3) hyp. u et v continues sur Omega
MG=siMesh.Masse(Th,'labels',labels,'local',false);
UG=siMesh.eval(Th,u,'labels',labels,'local',false);
VG=siMesh.eval(Th,v,'labels',labels,'local',false);
IG=dot(MG*UG,VG);
fprintf('        |IL-IG|=%.6e  [IG, methode 3]\n',abs(IG-IL))

%%%   c. calcul integrale de u.v sur Omega

% calcul exact (par rapport aux bbox des bords)
i=Th.find(1,10);
Iex=IexfunEh(Th.sTh{i}); 
for lab=[1,3,5,7] 
  i=Th.find(1,lab);
  Iex=Iex-IexfunEh(Th.sTh{i});
end

[idx,labels]=Th.find(2);
ML=siMesh.Masse(Th,'labels',labels,'local',true);
UL=siMesh.eval(Th,u,'labels',labels,'local',true);
VL=siMesh.eval(Th,v,'labels',labels,'local',true);
IL=0; % Calculs locaux
for i=1:length(labels)
  IL=IL+dot(ML{i}*UL{i},VL{i});
end
fprintf('(1) c.  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=siMesh.Masse(Th,'labels',lab,'local',false);
  UG=siMesh.eval(Th,u,'labels',lab,'local',false);
  VG=siMesh.eval(Th,v,'labels',lab,'local',false);
  IG=IG+dot(MG*UG,VG);
end
fprintf('        |IL-IG|=%.6e  [IG, methode 1]\n',abs(IG-IL))

% Calculs globaux (methode 2) hyp. u et v continues sur Omega
MG=siMesh.Masse(Th,'labels',labels,'local',false);
UG=Th.eval(u);
VG=Th.eval(v);
IG=dot(MG*UG,VG);
fprintf('        |IL-IG|=%.6e  [IG, methode 2]\n',abs(IG-IL))

% Calculs globaux (methode 3) hyp. u et v continues sur Omega
MG=siMesh.Masse(Th,'labels',labels,'local',false);
UG=siMesh.eval(Th,u,'labels',labels,'local',false);
VG=siMesh.eval(Th,v,'labels',labels,'local',false);
IG=dot(MG*UG,VG);
fprintf('        |IL-IG|=%.6e  [IG, methode 3]\n',abs(IG-IL))

