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

meshfile=fc_oogmsh.gmsh.buildmesh2d('condenser11',10);
Th=fc_simesh.siMesh(meshfile); % maillage global

u=@(x,y) cos(pi*x).*sin(pi/3*y)+1;

%%% 2-simplex in dim=2
i=Th.find(2,20); % get index of the 2-simplices elementary mesh with label=20
Eh=Th.sTh{i};

Ug=Th.eval(u); % Ug: dimension Th.nq
Ul=Eh.eval(u); % Ul: dimension Eh.nq
Ig=Eh.toGlobal;
Igc=setdiff(1:Th.nq,Ig); % complementaire de Ig

% prolongement par 0 dans Th (maillage global)
P=zeros(Th.nq,1);
P(Eh.toGlobal)=Ul;  

% Matrices locales
Ml=FEM.Masse(Eh.nq,Eh.me,Eh.vols);  
Kl=FEM.Rigidite(Eh.nq,Eh.me,Eh.vols,Eh.gradBaCo);

% Matrices globales
Mg=sparse(Th.nq,Th.nq);
Mg(Ig,Ig)=Ml;
Kg=sparse(Th.nq,Th.nq);
Kg(Ig,Ig)=Kl;

Ug=rand(Th.nq,1); % vecteur global
Ul=Ug(Ig); % vecteur local, restriction de Ug sur Eh

% Avec matrice de Masse
Vl=Ml*Ul; % calcul local
Vg=Mg*Ug; % calcul global

norm(Vl-Vg(Ig),Inf)
norm(Vg(Igc),Inf)


% Avec matrice de Rigidite
Wl=Kl*Ul; % calcul local
Wg=Kg*Ug; % calcul global

norm(Wl-Wg(Ig),Inf)
norm(Wg(Igc),Inf)
