%MakeMeshData A script that builds two circular meshes for 2D EIT Package demos 
% MakeMeshData is a script that builds two circular  
% meshes for 2D EIT Package demos. Denser mesh (mesh 2) is for forward computations and
% a coarse mesh (mesh 1) for inverse computations. See also 
% meshgen_eit2d.m , MakeElement and  MakeNode. 

clear
%S=2.5; % Length of the electrode.
%N=16;  % Number of the electrodes.
%r=14;  % Radius of the circle.
%style='s'; % 's' for the structured mesh and 'u' for the unstructured mesh

%[H1,g1,H2,g2,E1,E2,Ind2,Indb1,Indb2] = meshgen_eit2d(S,N,r,style);
%[Element1,Nodelist1]=MakeElement(H1,Indb1,E1);
%[Node1]=MakeNode(Element1,Nodelist1,g1);
%[Element2,Nodelist2]=MakeElement(H2,Indb2,E2);
%[Node2]=MakeNode(Element2,Nodelist2,g2);



%r=[14,12,10,8,5,3,0]; % Radii of rings in course mesh
r=[10,9.13,8.25,7.38,6.5,5.64,4.76,3.89,3,2.14,1.27,0];

eI=[1,2];% Number of elements in the electrode and between electrodes

%N=[[32,26,20,16,12,8]/2*3,1]; % Number of nodes in each ring
N=[[24,22,20,18,15,13,11,9,7,5,3]*sum(eI),1];


[g1,gp,H1,E1]=cirgrid_eit(r,N,eI); % make the course circular grid
% G 1 is the coods of nodes H1 is a list of nodes of each triangle
% gp is the polar coordinates (not needed)
% Each row of E1 is a list of elements under an electrode
[g2,H2,Ind2] = RefineMesh(g1,H1); % Make a finer mesh for the potential
% Mesh 1 is for resistivity Mesh 2 is for potential
% Ind2 which elements of mesh 2 are in which element of mesh 1
%     a sparse boolean array

%Now calculate the electrode information for the fine mesh
E2=[];
for ii=1:size(E1,1)
 Eii=E1(ii,:);
 E2ii=[];
  for jj=1:size(Eii,2)
   Ind=find(Ind2(:,Eii(jj)));
   E2ii=[E2ii,Ind(1:2)']; 
  end
 E2=[E2;E2ii];
end 

Indb1=findboundary(g1,H1); % Index of boundary
Indb2=findboundary(g2,H2);

[Element1,Nodelist1]=MakeElement(H1,Indb1,E1);
[Node1]=MakeNode(Element1,Nodelist1,g1);
[Element2,Nodelist2]=MakeElement(H2,Indb2,E2);
[Node2]=MakeNode(Element2,Nodelist2,g2);

save('hsmeshdata24_3','Node1','Element1','Node2','Element2','Ind2')
% Element1(k).Topology is H(1,:) etc
% Element1(k).Face(:,1) is the three eage of this element
% Element1(k).Face(:,2) is neighbouring elements 0 if no neighbour
% Element1(k).Face(:,3) is electrode number ,0 if not on an electrode
% Node1(k).Coordinate is coordinates kth node