% Create a hexahedral block with a number of spherical holes in it.
% The holes may intersect each other but must do so with a clear margin.
% So this is like a Swiss cheese. 

% Number of holes
N=20

% Size of block
A=[2,3,4]'

% Average radius of holes
Rave=0.4;

% Epsilon to determine the "clear" cut
ceps=0.1;

% Variation (with an even distribution) in the size of the holes
cvar=0.3;

Reps=ceps*Rave;
Rvar=cvar*Rave;

%Total volume of the brick
Vtot=A(1)*A(2)*A(3);

%Volume of the holes neglecting cuts & overlaps
Vholes=0.0;

misses=0;
for i=1:N
  while 1 

% Random position and radius for a ball
    X=(A.*rand(3,1)+2*(Rave-Reps)*rand(3,1))-(Rave-Reps);
    R=Rave+(2*rand(1,1)-1)*Rvar;
    if i==1
      Xall=X;
      Rall=R;
      break;
    else
      cut=0;

% Check that the ball isn't too close to another ball that it would create too small elements
      for j=1:i-1
	dx=X-Xall(:,j);   
	s=sqrt(dx'*dx);
        r1=R;
        r2=Rall(j);
        dr=abs(r1-r2);
        sr=r1+r2;
        if (s<(r1+r2)/2) 
          cut=1;
          break;
        end
        if(abs(s-sr)<Reps) 
          cut=2;
          break;
        end 
      end

% Check that the ball isn't too close to the edge so that it would create too small elements
      for j=1,3
	if(abs(X(j)-R)<Reps | abs(A(j)-X(j)-R)<Reps) 
          cut=3;
          break;  
        end
      end 

% If there is no critical cut add the ball to the list of balls, otherwise continue search
      if(~cut) 
        Xall=[Xall X];
        Rall=[Rall R];
        Vholes=Vholes+(4*pi/3)*R^3;
        break;
      else
        misses=misses+1;
      end   
    end    
  end 
end 

Xall
Rall
misses
qholes=Vholes/Vtot

% write the file in .geo format for netgen
fid = fopen('cheese.geo','w');

fprintf(fid,'algebraic3d\n');
fprintf(fid,'solid rectangle = orthobrick(0,0,0;%6.3f,%6.3f,%6.3f);\n',A(1),A(2),A(3));
for i=1:N
  fprintf(fid,'solid ball%d = sphere(%6.3f,%6.3f,%6.3f; %6.3f);\n',i,Xall(1,i),Xall(2,i),Xall(3,i),Rall(i));
end 
fprintf(fid,'solid rest = rectangle');
for i=1:N
  fprintf(fid,' and not ball%d',i);
end 
fprintf(fid,';\n');
fprintf(fid,'tlo rest -col=[0,0,1];\n');
fclose(fid);
