% Matrice de Masse locale et globale sur Omega
% en utilisant siMeshElt.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:
% ------- sage ---------
% var('x,y,a,b,c,d'); 
% u(x,y)=cos(x+y-pi/3);v(x,y)=sin(x-y+1);
% integrate(integrate(u(x,y)*v(x,y),(x,a,b)),(y,c,d))
% ----------------------
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)
i=Th.find(2,20);
ML=siMeshElt.Masse(Th.sTh{i}); % Matrice locale
UL=Th.sTh{i}.eval(u);
VL=Th.sTh{i}.eval(v);
IL=dot(ML*UL,VL);
Iex=IexfunEh(Th.sTh{i});
fprintf('(1) a.  Iex=%.5e, |IL-Iex|=%.6e \n',Iex,abs(Iex-IL))

MG=siMeshElt.Masse(Th.sTh{i},'local',false); % Matrice globale
UG=siMeshElt.relevement(Th.sTh{i},UL,'funinit',@rand);
VG=siMeshElt.relevement(Th.sTh{i},VL,'funinit',@rand);
%UG=zeros(Th.nq,1);UG(Th.sTh{i}.toGlobal)=UL; 
%VG=zeros(Th.nq,1);VG(Th.sTh{i}.toGlobal)=VL; 
IG=dot(MG*UG,VG);


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

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

IL=0; % Calculs locaux
for lab=labels
  i=Th.find(2,lab);
  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) 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
  UG=siMeshElt.relevement(Th.sTh{i},Th.sTh{i}.eval(u),'funinit',@rand);
  VG=siMeshElt.relevement(Th.sTh{i},Th.sTh{i}.eval(v),'funinit',@rand);
  %   ou
  % 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);
  %   ou
  % UG=Th.eval(u); VG=Th.eval(v);
  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=sparse(Th.nq,Th.nq);
for lab=labels
  i=Th.find(2,lab);
  MG=MG+siMeshElt.Masse(Th.sTh{i},'local',false); 
end
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=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);
  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))


%%%   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);
IL=0; % Calculs locaux
for i=idx
  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) 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
  UG=zeros(Th.nq,1);UG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(u); % ou utiliser siMeshElt.relevement
  VG=zeros(Th.nq,1);VG(Th.sTh{i}.toGlobal)=Th.sTh{i}.eval(v); % ou utiliser siMeshElt.relevement
  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=sparse(Th.nq,Th.nq);
for lab=labels
  i=Th.find(2,lab);
  MG=MG+siMeshElt.Masse(Th.sTh{i},'local',false); 
end
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=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);
  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))

