fileID = fopen('N60_12.dat','r');
formatSpec = '%f %f %f';
sizeA = [3 Inf]; % column order, [x y radius]

A = fscanf(fileID,formatSpec,sizeA);
A = A';
fclose(fileID);

% since the inclusions are generated for DEM, they are all in contact
% scale their radius to separate them from each other

factor = 0.9;

A(:,3) = A(:,3)*factor;

numberOfInclusions = size(A,1);

th = 0:pi/50:2*pi;

% plot the inclusions 

figure
hold on

for i=1:numberOfInclusions
    xunit = A(i,3) * cos(th) + A(i,1);
    yunit = A(i,3) * sin(th) + A(i,2);
    plot(xunit, yunit);
end
axis equal

% plot the outer box
dL = 0.1;
L  = 1;

plot([0-dL L+dL],[-dL -dL]);
plot([L+dL L+dL],[-dL L+dL]);
plot([L+dL -dL],[L+dL L+dL]);
plot([-dL -dL],[L+dL -dL]);

%%%%%%%%%%%%%%
% make a Gmsh geo file

gmshFile = fopen('N60.geo','w');

h = 0.03;
fprintf(gmshFile, 'h=%f;\n',h);
fprintf(gmshFile, 'h2=%f;\n',h*2);

for i=1:numberOfInclusions
    x = A(i,1);
    y = A(i,2);
    r = A(i,3);
    
    fprintf(gmshFile, 'Point(%i)={%22.14E,%22.14E,0,h2};\n',10*i,  x,y  );
    fprintf(gmshFile, 'Point(%i)={%22.14E,%22.14E,0,h};\n',10*i+1,x+r,y);
    fprintf(gmshFile, 'Point(%i)={%22.14E,%22.14E,0,h};\n',10*i+2,x,y+r);
    fprintf(gmshFile, 'Point(%i)={%22.14E,%22.14E,0,h};\n',10*i+3,x-r,y);
    fprintf(gmshFile, 'Point(%i)={%22.14E,%22.14E,0,h};\n',10*i+4,x,y-r);
    
    fprintf(gmshFile, 'Circle(%i)={%i,%i,%i};\n',10*i,  10*i+1,10*i,10*i+2);
    fprintf(gmshFile, 'Circle(%i)={%i,%i,%i};\n',10*i+1,10*i+2,10*i,10*i+3);
    fprintf(gmshFile, 'Circle(%i)={%i,%i,%i};\n',10*i+2,10*i+3,10*i,10*i+4);
    fprintf(gmshFile, 'Circle(%i)={%i,%i,%i};\n',10*i+3,10*i+4,10*i,10*i+1);
    
    fprintf(gmshFile, 'Line Loop(%i)={%i,%i,%i,%i};\n',i,10*i,10*i+1,10*i+2,10*i+3);
    
    fprintf(gmshFile, 'Plane Surface(%i)={%i};\n',i,i);
end

fprintf(gmshFile, '\n');

% use an array to faciliate Plane definition later

fprintf(gmshFile, 'For i In {1:%i}\n',numberOfInclusions);
fprintf(gmshFile, 'inclusions[i]=i;\n');
fprintf(gmshFile, 'fields[i]=2*i;\n');
fprintf(gmshFile, 'EndFor\n\n');

% put the inclusions in a box (1,2,3,4)

fprintf(gmshFile, 'Point(%i)={%22.4E,%22.4E,0,h};\n',1,  -dL,-dL  );
fprintf(gmshFile, 'Point(%i)={%22.4E,%22.4E,0,h};\n',2,  L+dL,-dL );
fprintf(gmshFile, 'Point(%i)={%22.4E,%22.4E,0,h};\n',3,  L+dL,L+dL);
fprintf(gmshFile, 'Point(%i)={%22.4E,%22.4E,0,h};\n',4,  -dL,L+dL );

fprintf(gmshFile, 'Line(%i)={%i,%i};\n',1,1,2);
fprintf(gmshFile, 'Line(%i)={%i,%i};\n',2,2,3);
fprintf(gmshFile, 'Line(%i)={%i,%i};\n',3,3,4);
fprintf(gmshFile, 'Line(%i)={%i,%i};\n',4,4,1);

fprintf(gmshFile, 'Line Loop(%i)={%i,%i,%i,%i};\n',1000,1,2,3,4);

fprintf(gmshFile, 'Plane Surface(1000)={1000,inclusions[{1:%i}]};\n',60);

% physical quantities

fprintf(gmshFile, 'Physical Line(%i)={%i};\n',1,1); % lower edge
fprintf(gmshFile, 'Physical Line(%i)={%i};\n',2,2); % right edge
fprintf(gmshFile, 'Physical Line(%i)={%i};\n',3,3); % upper edge
fprintf(gmshFile, 'Physical Line(%i)={%i};\n',4,4); % left edge

fprintf(gmshFile, 'Physical Surface(%i)={%i};\n',1,1000); % matrix
fprintf(gmshFile, 'Physical Surface(%i)={inclusions[{1:%i}]};\n\n',2,numberOfInclusions);

% inclusions: coarse elements inside


% for i=1:numberOfInclusions
%     r = A(i,3);
%     
%     fprintf(gmshFile, 'Field[%i]=Attractor;\n',2*i-1);
%     fprintf(gmshFile, 'Field[%i].EdgesList={%i,%i,%i,%i};\n',2*i-1,10*i,10*i+1,10*i+2,10*i+3);
%     fprintf(gmshFile, 'Field[%i]=Threshold;\n',2*i);
%     fprintf(gmshFile, 'Field[%i].IField=1;\n',2*i);
%     fprintf(gmshFile, 'Field[%i].LcMin=h;\n',2*i);
%     fprintf(gmshFile, 'Field[%i].LcMax=h2;\n',2*i);
%     fprintf(gmshFile, 'Field[%i].DistMin=%f;\n',2*i,0.1*r);
%     fprintf(gmshFile, 'Field[%i].DistMax=%f;\n',2*i,0.3*r);
% end
% 
% fprintf(gmshFile, 'Field[%i]=Min;\n',2*numberOfInclusions+1);
% fprintf(gmshFile, 'Field[%i].FieldsList={fields[{1:%i}]};\n',2*numberOfInclusions+1,numberOfInclusions);
% fprintf(gmshFile, 'Background Field=%i;\n',2*numberOfInclusions+1);
% fprintf(gmshFile, 'Mesh.CharacteristicLengthExtendFromBoundary=0;\n');

% another solution

for i=1:numberOfInclusions
    fprintf(gmshFile, 'Point{%i} In Surface {%i};\n',10*i,i);
end

