randomly generated spheres in a volume in Matlab
104 views (last 30 days)
Show older comments
Hey,
I want to generate n spheres of defined radius with in a defined volume at random locations. Following the condition of not intersecting each other.
Can someone help me with this.
Regards
Om
Accepted Answer
Jan
on 13 May 2019
Edited: Jan
on 13 May 2019
function P = GetRandomSpheres(nWant, Width, Radius)
% INPUT:
% nWant: Number of spheres
% Width: Dimension of 3d box as [1 x 3] double vector
% Radius: Radius of spheres
% OUTPUT:
% P: [nWant x 3] matrix, centers
P = zeros(nWant, 3);
R2 = (2 * Radius) ^ 2; % Squared once instead of SQRT each time
W = Width - 2 * Radius; % Avoid interesction with borders
iLoop = 1; % Security break to avoid infinite loop
nValid = 0;
while nValid < nWant && iLoop < 1e6
newP = rand(1, 3) .* W + Radius;
% Auto-expanding, need Matlab >= R2016b. For earlier versions:
% Dist2 = sum(bsxfun(@minus, P(1:nValid, :), newP) .^ 2, 2);
Dist2 = sum((P(1:nValid, :) - newP) .^ 2, 2);
if all(Dist2 > R2)
% Success: The new point does not touch existing sheres:
nValid = nValid + 1; % Append this point
P(nValid, :) = newP;
end
iLoop = iLoop + 1;
end
% Stop if too few values have been found:
if nValid < nWant
error('Cannot find wanted number of points in %d iterations.', iLoop)
end
end
And a demo code:
n = 10;
R = 10;
P = GetRandomSpheres(n, [100, 100, 100], R);
figure
axes('NextPlot', 'add', ...
'XLim', [0, 100], 'YLim', [0, 100], 'ZLim', [0, 100]);
view(3);
[X, Y, Z] = sphere();
for k = 1:n
surf(X * R + P(k, 1), Y * R + P(k, 2), Z * R + P(k, 3));
end
28 Comments
Russell Jones
on 22 Apr 2023
I am looking into a similar problem where my volume fraction is 0.59, and from testing @Jan's modified code, it unfortunately seems to only work for volume fractions of 0.31 or lower. I have increased the iterations to 1e9 and I still can't get above 0.31. From what I can tell, it would require a different approach. Did you have any success with a different method / does anyone know what method I should attempt for a packing ratio of 0.59?
Walter Roberson
on 22 Apr 2023
@Jan wrote code for the requirement that the spheres be non-intersecting. To get to the volume fractions you are talking about, you would need intersecting (touching) spheres, and you will need to use some sort of compaction algorithm. As the theoretical limit is only 63.4% you will need to do a fair amount of compaction to get 0.59. Or you could deliberately generate one of the regular packings.
More Answers (3)
Mehdi Chouchane
on 4 Mar 2020
Here is a code I'm using which is a custom version of a code taken from Matlab Community.
It generates the desired volume fraction of non-overlapping spheres with radii varying between rmin and rmax.
function [centers, rads] = sampleSpheres( Fract,dims,rmin,rmax,verbosity)
% main function which is to be called for adding multiple spheres in a cube
% dims is assumed to be a row vector of size 1-by-ndim
% For demo take dims = [ 2 2 2 ] and n = 50
% preallocate
V = Fract * dims(1) * dims(2) * dims(3) ;
v = 0 ;
ndim = numel(dims);
n = round(1.5*V / ((4*pi*((rmax+rmin)/2)^3)/3)) ;
centers = zeros( n, ndim );
rads = zeros( n, 1 );
ii = 1;
n = 0 ;
while v < V
[centers(ii,:), rads(ii) ] = randomSphere2( dims,rmin,rmax);
if nonOver2( centers(1:ii,:), rads(1:ii))
n = n + 1 ;
v = v + (4*pi*rads(ii)^3)/(3) ;
ii = ii + 1; % accept and move on
if verbosity == 1
100*v/V
end
end
end
centers = centers(~all(centers == 0, 2),:);
rads = rads(~all(rads == 0, 2),:);
end
function [ c, r ] = randomSphere2( dims,rmin,rmax)
% creating one sphere at random inside [0..dims(1)]x[0..dims(2)]x...
% In general cube dimension would be 10mm*10mm*10mm
% radius and center coordinates are sampled from a uniform distribution
% over the relevant domain.
% output: c - center of sphere (vector cx, cy,... )
% r - radius of sphere (scalar)
r = rmin + ( rmax - rmin) .* rand(1);
c = bsxfun(@times,(dims - r) , rand(1,3)) + r; % to make sure sphere is placed inside the cube
end
function not_ovlp = nonOver2( centers, rads) % to make sure spheres do not overlap
if numel( rads ) == 1
not_ovlp = true;
return; % nothing to check for a single or the first sphere
end
center_dist = sqrt(sum(bsxfun(@minus,centers(1:end-1,:),centers(end,:)).^2,2));
radsum = rads(end) + rads(1:end-1);
not_ovlp = all(center_dist >= radsum);
return;
end
6 Comments
Russell Jones
on 25 Apr 2023
I tried using your code and it worked better to achieve a higher volumetric fraction. However, I am unable to get the code to work for volumetric fractions greater than 0.5. As I need 0.59, I am struggling to figure out how to add to the code to get to this fraction. Do you have any tips?
Walter Roberson
on 25 Apr 2023
As I already discussed with you, this kind of code is for creating non-overlapping spheres that might have any amount of distance between them (though the larger the gap the more likely that another sphere will be randomly placed in the gap volume.)
In order to get the kinds of packings you need, you need code that actively compacts.
For example you might want to see if you can find the discussion from 5-ish years ago from the person who needed to model accreation, which was handled by "dropping" particles from the top and letting them "fall" until they were in a stable situation. That kind of code might still not get you close enough for your purposes, but it would be better than the random-placement-in-a-volume approaches.
sadedbouta
on 21 Apr 2020
Hey, are you help me for this lines in script matlab with Comsol:
n = 100;
Vsum = 0;
% Generate position vectors
Pos = zeros(n,3) ; % XYZ
R = zeros(n,3);
idx = 1; % index for sphere
flag = 0;
while (Vsum < Vsq * vf)
r = abs( normrnd(miu,sigma) );
% generates a random number from the normal distribution with mean parameter mu and
% standard deviation parameter sigma.
pos = [blk_size * rand(1,1) blk_size * rand(1,1) blk_size * rand(1,1)];
for k = 1:idx %Check the distance between the randomly generated sphere and all existing spheres.
Distance = sqrt((pos(1)-Pos(k,1))^2+(pos(2)-Pos(k,2))^2+(pos(3)-Pos(k,3))^2);
rsum = r
%rsum = r+R(k);
if Distance < 2*rsum
flag = 1;
break;
end
end
Can you idea me corrected the error in this script?
1 Comment
Serigne Saliou Sene
on 31 Jan 2022
@JanHello Jan I hope you are well, I would like you to improve my script I have already randomly generated several ellipsoides in the cube I want a plot (volumeCube +aggregate+fiber) Thanks in advance
clear ; clc; close all;
%input
%------------------
H=150; %Cube Height (mm)
W=150; %Cube Width (mm)
L=150; %Cube Length (mm)
x=[0 L L 0]; %Specimen dimensions
y=[0 W]; %Specimen dimensions
z=[0 0 H H]; %Specimen dimensions
Classes_diameters=[19 12.70 9.5 4.75 2.36]; %Particles classes diameters (descendingly)
Alpha=0.45; %Fuller's curve exponent
m=3; %Particles shape distribution factor
Particle_ratio=0.50; %Particles ratio of the volume
er=0.05; %Spacing factor between particles
dist=W/2; %Cutting distance for ellipses
r_min=2.36; %Minimum ellipse radius involved
L=3; %Length of fibers
N=1200; %Number of fibers
DFiber=3; %Diameter of fibers
Orientation=[]; %Orientation: Orientation of fibers as [l; m; n] column vector where l, m, and n are direction cosines of orientation in x, y, and z directions, respectively.
Ndiv=1; %Number of fiber mesh divisions
%------------------------
Classes=Particles_Generation(x,y,z,Classes_diameters,Alpha,m,Particle_ratio);
Plot_Sieve(Classes,x,y,z,Classes_diameters,Alpha,Particle_ratio);
Ellipsoids=Particles_Distribution(Classes,x,y,z,er);
Plot_Ellipsoids(Ellipsoids,x,y,z);
Ellipses=Ellipsoids_to_Ellipses(Ellipsoids,dist,r_min);
Plot_Ellipses(Ellipses,x,z);
[Nodes_Fibers, Fibers]=Generate_Fiber(x,y,z,L,N,DFiber,Orientation,Ndiv,Ellipsoids);
Plot_Fiber(x,y,z,Nodes_Fibers,Fibers,DFiber);
Plot_Ellipsoids_Fiber(Ellipsoids,x,y,z,Nodes_Fibers,Fibers,DFiber);
3 Comments
Serigne Saliou Sene
on 2 Feb 2022
Yes this is the function I used.
I must add a 6th plot in which I would have the aggregates and the fibers in the packaging (Volume) of the cube clearly visible in order to form the ITZ, as in the photo. Thanks in advance
Yahriel Salinas-Reyes
on 14 Aug 2023
Hello,
I am running into the similar problem and trying to plot the Interfacial Transition Zone (ITZ) all in one 3d plot. Where you able to find a solution or have any help to offer?
See Also
Categories
Find more on Descriptive Statistics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!