1
// © 2009 Václav Šmilauer <eudoxos@arcig.cz>
3
#include<yade/pkg/dem/SpherePack.hpp>
5
#include<yade/core/Omega.hpp>
6
#include<yade/core/Scene.hpp>
7
#include<yade/pkg/common/Sphere.hpp>
8
#include<yade/pkg/dem/Shop.hpp>
10
#include <boost/random/linear_congruential.hpp>
11
#include <boost/random/uniform_real.hpp>
12
#include <boost/random/variate_generator.hpp>
14
#include<yade/core/Timing.hpp>
16
// not a serializable in the sense of YADE_PLUGIN
18
CREATE_LOGGER(SpherePack);
21
using namespace boost;
22
namespace py=boost::python;
25
void SpherePack::fromList(const py::list& l){
27
size_t len=py::len(l);
28
for(size_t i=0; i<len; i++){
29
const py::tuple& t=py::extract<py::tuple>(l[i]);
30
py::extract<Vector3r> vec(t[0]);
31
if(vec.check()) { pack.push_back(Sph(vec(),py::extract<double>(t[1]),(py::len(t)>2?py::extract<int>(t[2]):-1))); continue; }
32
PyErr_SetString(PyExc_TypeError, "List elements must be (Vector3, float) or (Vector3, float, int)!");
33
py::throw_error_already_set();
37
void SpherePack::fromLists(const vector<Vector3r>& centers, const vector<Real>& radii){
39
if(centers.size()!=radii.size()) throw std::invalid_argument(("The same number of centers and radii must be given (is "+lexical_cast<string>(centers.size())+", "+lexical_cast<string>(radii.size())+")").c_str());
40
size_t l=centers.size();
41
for(size_t i=0; i<l; i++){
42
add(centers[i],radii[i]);
44
cellSize=Vector3r::Zero();
47
py::list SpherePack::toList() const {
49
FOREACH(const Sph& s, pack) ret.append(s.asTuple());
53
void SpherePack::fromFile(string file) {
54
typedef pair<Vector3r,Real> pairVector3rReal;
55
vector<pairVector3rReal> ss;
57
ss=Shop::loadSpheresFromFile(file,mn,mx,&cellSize);
59
FOREACH(const pairVector3rReal& s, ss) pack.push_back(Sph(s.first,s.second));
62
void SpherePack::toFile(const string fname) const {
63
ofstream f(fname.c_str());
64
if(!f.good()) throw runtime_error("Unable to open file `"+fname+"'");
65
if(cellSize!=Vector3r::Zero()){ f<<"##PERIODIC:: "<<cellSize[0]<<" "<<cellSize[1]<<" "<<cellSize[2]<<endl; }
66
FOREACH(const Sph& s, pack){
67
if(s.clumpId>=0) throw std::invalid_argument("SpherePack with clumps cannot be (currently) exported to a text file.");
68
f<<s.c[0]<<" "<<s.c[1]<<" "<<s.c[2]<<" "<<s.r<<endl;
73
void SpherePack::fromSimulation() {
75
Scene* scene=Omega::instance().getScene().get();
76
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
78
shared_ptr<Sphere> intSph=dynamic_pointer_cast<Sphere>(b->shape);
80
pack.push_back(Sph(b->state->pos,intSph->radius,(b->isClumpMember()?b->clumpId:-1)));
82
if(scene->isPeriodic) {
83
cellSize=scene->cell->getSize();
88
long SpherePack::makeCloud(Vector3r mn, Vector3r mx, Real rMean, Real rRelFuzz, int num, bool periodic, Real porosity, const vector<Real>& psdSizes, const vector<Real>& psdCumm, bool distributeMass, int seed, Matrix3r hSize){
89
isPeriodic = periodic;
90
static boost::minstd_rand randGen(seed!=0?seed:(int)TimingInfo::getNow(/* get the number even if timing is disabled globally */ true));
91
static boost::variate_generator<boost::minstd_rand&, boost::uniform_real<Real> > rnd(randGen, boost::uniform_real<Real>(0,1));
92
vector<Real> psdRadii; // holds plain radii (rather than diameters), scaled down in some situations to get the target number
93
vector<Real> psdCumm2; // psdCumm but dimensionally transformed to match mass distribution
95
bool hSizeFound =(hSize!=Matrix3r::Zero());//is hSize passed to the function?
97
if (!hSizeFound) {hSize=size.asDiagonal();}
98
if (hSizeFound && !periodic) LOG_WARN("hSize can be defined only for periodic cells.");
99
Real volume=hSize.determinant();
100
Matrix3r invHsize =hSize.inverse();
101
Real area=abs(size[0]*size[2]+size[0]*size[1]+size[1]*size[2]);//2 terms will be null if one coordinate is 0, the other is the area
103
if (hSizeFound || !periodic) throw invalid_argument("The box as null volume, this is not supported. One exception is for flat box defined by min-max in periodic conditions, if hSize argument defines a non-null volume (or if the hSize argument is left undefined).");
104
else LOG_WARN("The volume of the min-max box is null, we will assume that the packing is 2D. If it is not what you want then you defined wrong input values; check that min and max corners are defined correctly.");}
105
int mode=-1; bool err=false;
106
// determine the way we generate radii
107
if(porosity<=0) {LOG_WARN("porosity must be >0, changing it for you. It will be ineffective if rMean>0."); porosity=0.5;}
108
//If rMean is not defined, then in will be defined in RDIST_NUM
109
if(rMean>0) mode=RDIST_RMEAN;
110
else if(num>0 && psdSizes.size()==0) {
112
// the term (1+rRelFuzz²) comes from the mean volume for uniform distribution : Vmean = 4/3*pi*Rmean*(1+rRelFuzz²)
113
if (volume) rMean=pow(volume*(1-porosity)/(Mathr::PI*(4/3.)*(1+rRelFuzz*rRelFuzz)*num),1/3.);
114
else {//The volume is null, we will generate a 2D packing with the following rMean
115
if (!area) throw invalid_argument("The box defined has null volume AND null surface. Define at least maxCorner of the box, or hSize if periodic.");
116
rMean=pow(area*(1-porosity)/(Mathr::PI*(1+rRelFuzz*rRelFuzz)*num),0.5);}
118
// transform sizes and cummulated fractions values in something convenient for the generation process
119
if(psdSizes.size()>0){
120
err=(mode>=0); mode=RDIST_PSD;
121
if(psdSizes.size()!=psdCumm.size()) throw invalid_argument(("SpherePack.makeCloud: psdSizes and psdCumm must have same dimensions ("+lexical_cast<string>(psdSizes.size())+"!="+lexical_cast<string>(psdCumm.size())).c_str());
122
if(psdSizes.size()<=1) throw invalid_argument("SpherePack.makeCloud: psdSizes must have at least 2 items");
123
if((*psdCumm.begin())!=0. && (*psdCumm.rbegin())!=1.) throw invalid_argument("SpherePack.makeCloud: first and last items of psdCumm *must* be exactly 0 and 1.");
124
psdRadii.reserve(psdSizes.size());
125
for(size_t i=0; i<psdSizes.size(); i++) {
126
psdRadii.push_back(/* radius, not diameter */ .5*psdSizes[i]);
128
//psdCumm2 is first obtained by integrating the number of particles over the volumic PSD (dN/dSize = totV*(dPassing/dSize)*1/(4/3πr³)). The total cumulated number will be the number of spheres in volume*(1-porosity), it is used to decide if the PSD will be scaled down. psdCumm2 is normalized below in order to fit in [0,1]. (Bruno C.)
129
if (i==0) psdCumm2.push_back(0);
130
else psdCumm2.push_back(psdCumm2[i-1] + 3.0* (volume?volume:(area*psdSizes[psdSizes.size()-1])) *(1-porosity)/Mathr::PI*(psdCumm[i]-psdCumm[i-1])/(psdSizes[i]-psdSizes[i-1])*(pow(psdSizes[i-1],-2)-pow(psdSizes[i],-2)));
132
LOG_DEBUG(i<<". "<<psdRadii[i]<<", cdf="<<psdCumm[i]<<", cdf2="<<(distributeMass?lexical_cast<string>(psdCumm2[i]):string("--")));
133
// check monotonicity
134
if(i>0 && (psdSizes[i-1]>psdSizes[i] || psdCumm[i-1]>psdCumm[i])) throw invalid_argument("SpherePack:makeCloud: psdSizes and psdCumm must be both non-decreasing.");
136
// check the consistency between sizes, num, and poro if all three are imposed. If target number will not fit in (1-poro)*volume, scale down particles sizes
140
if (psdCumm2[psdSizes.size()-1]<num) appliedPsdScaling=pow(psdCumm2[psdSizes.size()-1]/num,1./3.);
143
for(size_t i=1; i<psdSizes.size(); i++) totVol+= 4/3*Mathr::PI*(psdCumm[i]-psdCumm[i-1])*num*
144
pow(0.5*(psdSizes[i]+psdSizes[i-1]),3)*(1+pow(0.5*(psdSizes[i]-psdSizes[i-1]),2));
145
Real volumeRatio = totVol/((1-porosity)*(volume?volume:(area*psdSizes[psdSizes.size()-1])));
146
if (volumeRatio>1) appliedPsdScaling=pow(volumeRatio,-1./3.);
148
if (appliedPsdScaling<1) for(size_t i=0; i<psdSizes.size(); i++) psdRadii[i]*=appliedPsdScaling;
150
//Normalize psdCumm2 so it's between 0 and 1
151
if(distributeMass) for(size_t i=1; i<psdSizes.size(); i++) psdCumm2[i]/=psdCumm2[psdSizes.size()-1];
153
if(err || mode<0) throw invalid_argument("SpherePack.makeCloud: at least one of rMean, porosity, psdSizes & psdCumm arguments must be specified. rMean can't be combined with psdSizes.");
154
// adjust uniform distribution parameters with distributeMass; rMean has the meaning (dimensionally) of _volume_
155
const int maxTry=1000;
156
if(periodic && volume && !hSizeFound)(cellSize=size);
158
for(int i=0; (i<num) || (num<0); i++) {
160
//Determine radius of the next sphere that will be placed in space. If (num>0), generate radii the deterministic way, in decreasing order, else radii are stochastic since we don't know what the final number will be
161
if (num>0) rand = ((Real)num-(Real)i+0.5)/((Real)num+1.);
166
//FIXME : r is never defined, it will be zero at first iteration, but it will have values in the next ones. Some magic?
168
if(distributeMass) r=pow3Interp(rand,rMean*(1-rRelFuzz),rMean*(1+rRelFuzz));
169
else r=rMean*(2*(rand-.5)*rRelFuzz+1); // uniform distribution in rMean*(1±rRelFuzz)
173
int piece=psdGetPiece(rand,psdCumm2,norm);
174
r=pow3Interp(norm,psdRadii[piece],psdRadii[piece+1]);
176
int piece=psdGetPiece(rand,psdCumm,norm);
177
r=psdRadii[piece]+norm*(psdRadii[piece+1]-psdRadii[piece]);}
179
// try to put the sphere into a free spot
180
for(t=0; t<maxTry; ++t){
182
if(!periodic) { for(int axis=0; axis<3; axis++) c[axis]=mn[axis]+(size[axis]?(size[axis]-2*r)*rnd()+r:0);}//we handle 2D with the special case size[axis]==0
183
else { for(int axis=0; axis<3; axis++) c[axis]=rnd();//coordinates in [0,1]
184
c=mn+hSize*c;}//coordinates in reference frame (inside the base cell)
185
size_t packSize=pack.size(); bool overlap=false;
186
if(!periodic) for(size_t j=0;j<packSize;j++) {if(pow(pack[j].r+r,2)>=(pack[j].c-c).squaredNorm()) {overlap=true; break;}}
188
for(size_t j=0; j<packSize; j++){
189
Vector3r dr=Vector3r::Zero();
190
if (!hSizeFound) {//The box is axis-aligned, use the wrap methods
191
for(int axis=0; axis<3; axis++) dr[axis]=size[axis]? min(cellWrapRel(c[axis],pack[j].c[axis],pack[j].c[axis]+size[axis]),cellWrapRel(pack[j].c[axis],c[axis],c[axis]+size[axis])) : 0;
192
} else {//not aligned, find closest neighbor in a cube of size 1, then transform distance to cartesian coordinates
193
Vector3r c1c2=invHsize*(pack[j].c-c);
194
for(int axis=0; axis<3; axis++){
195
if (abs(c1c2[axis])<abs(c1c2[axis] - Mathr::Sign(c1c2[axis]))) dr[axis]=c1c2[axis];
196
else dr[axis] = c1c2[axis] - Mathr::Sign(c1c2[axis]);}
197
dr=hSize*dr;//now in cartesian coordinates
199
if(pow(pack[j].r+r,2)>= dr.squaredNorm()){ overlap=true; break; }
202
if(!overlap) { pack.push_back(Sph(c,r)); break; }
206
if (mode!=RDIST_RMEAN) {
207
//if rMean is not imposed, then we call makeCloud recursively, scaling the PSD down until the target num is obtained
208
Real nextPoro = porosity+(1-porosity)/10.;
209
LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<num<<". Trying again with porosity "<<nextPoro<<". The size distribution is being scaled down");
211
return makeCloud(mn, mx, -1., rRelFuzz, num, periodic, nextPoro, psdSizes, psdCumm, distributeMass,seed,hSizeFound?hSize:Matrix3r::Zero());}
212
else LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<num<<".");
216
if (appliedPsdScaling<1) LOG_WARN("The size distribution has been scaled down by a factor pack.appliedPsdScaling="<<appliedPsdScaling);
220
void SpherePack::cellFill(Vector3r vol){
222
for(int i=0; i<3; i++) count[i]=(int)(ceil(vol[i]/cellSize[i]));
223
LOG_DEBUG("Filling volume "<<vol<<" with cell "<<cellSize<<", repeat counts are "<<count);
227
void SpherePack::cellRepeat(Vector3i count){
228
if(cellSize==Vector3r::Zero()){ throw std::runtime_error("cellRepeat cannot be used on non-periodic packing."); }
229
if(count[0]<=0 || count[1]<=0 || count[2]<=0){ throw std::invalid_argument("Repeat count components must be positive."); }
230
size_t origSize=pack.size();
231
pack.reserve(origSize*count[0]*count[1]*count[2]);
232
for(int i=0; i<count[0]; i++){
233
for(int j=0; j<count[1]; j++){
234
for(int k=0; k<count[2]; k++){
235
if((i==0) && (j==0) && (k==0)) continue; // original cell
236
Vector3r off(cellSize[0]*i,cellSize[1]*j,cellSize[2]*k);
237
for(size_t l=0; l<origSize; l++){
238
const Sph& s=pack[l]; pack.push_back(Sph(s.c+off,s.r));
243
cellSize=Vector3r(cellSize[0]*count[0],cellSize[1]*count[1],cellSize[2]*count[2]);
246
int SpherePack::psdGetPiece(Real x, const vector<Real>& cumm, Real& norm){
247
int sz=cumm.size(); int i=0;
248
while(i<sz && cumm[i]<=x) i++; // upper interval limit index
249
if((i==sz-1) && cumm[i]<=x){ i=sz-2; norm=1.; return i;}
250
i--; // lower interval limit intex
251
norm=(x-cumm[i])/(cumm[i+1]-cumm[i]);
252
//LOG_TRACE("count="<<sz<<", x="<<x<<", piece="<<i<<" in "<<cumm[i]<<"…"<<cumm[i+1]<<", norm="<<norm);
256
py::tuple SpherePack::psd(int bins, bool mass) const {
257
if(pack.size()==0) return py::make_tuple(py::list(),py::list()); // empty packing
259
Real minD=std::numeric_limits<Real>::infinity(); Real maxD=-minD;
260
// volume, but divided by π*4/3
261
Real vol=0; long N=pack.size();
262
FOREACH(const Sph& s, pack){ maxD=max(2*s.r,maxD); minD=min(2*s.r,minD); vol+=pow(s.r,3); }
263
if(minD==maxD){ minD-=.5; maxD+=.5; } // emulates what numpy.histogram does
264
// create bins and bin edges
265
vector<Real> hist(bins,0); vector<Real> cumm(bins+1,0); /* cummulative values compute from hist at the end */
266
vector<Real> edges(bins+1); for(int i=0; i<=bins; i++){ edges[i]=minD+i*(maxD-minD)/bins; }
267
// weight each grain by its "volume" relative to overall "volume"
268
FOREACH(const Sph& s, pack){
269
int bin=int(bins*(2*s.r-minD)/(maxD-minD)); bin=min(bin,bins-1); // to make sure
270
if (mass) hist[bin]+=pow(s.r,3)/vol; else hist[bin]+=1./N;
272
for(int i=0; i<bins; i++) cumm[i+1]=min(1.,cumm[i]+hist[i]); // cumm[i+1] is OK, cumm.size()==bins+1
273
return py::make_tuple(edges,cumm);
276
/* possible enhacement: proportions parameter, so that the domain is not cube, but box with sides having given proportions */
277
long SpherePack::particleSD2(const vector<Real>& radii, const vector<Real>& passing, int numSph, bool periodic, Real cloudPorosity, int seed){
278
typedef Eigen::Matrix<Real,Eigen::Dynamic,Eigen::Dynamic> MatrixXr;
279
typedef Eigen::Matrix<Real,Eigen::Dynamic,1> VectorXr;
281
int dim=radii.size()+1;
282
if(passing.size()!=radii.size()) throw std::invalid_argument("SpherePack.particleSD2: radii and passing must have the same length.");
283
MatrixXr M=MatrixXr::Zero(dim,dim);
284
VectorXr rhs=VectorXr::Zero(dim);
287
We know percentages for each fraction (Δpi) and their radii (ri), and want to find
288
the number of sphere for each fraction Ni and total solid volume Vs. For each fraction,
289
we know that the volume is equal to Ni*(4/3*πri³), which must be equal to Vs*Δpi (Δpi is
290
relative solid volume of the i-th fraction).
292
The last equation says that total number of particles (sum of fractions) is equal to N,
293
which is the total number of particles requested by the user.
303
for(int i=0; i<dim-1; i++){
304
M(i,i)=(4/3.)*Mathr::PI*pow(radii[i],3);
305
M(i,dim-1)=-(passing[i]-(i>0?passing[i-1]:0))/100.;
309
// NumsVs=M^-1*rhs: number of spheres and volume of solids
310
VectorXr NumsVs(dim); NumsVs=M.inverse()*rhs;
311
Real Vs=NumsVs[dim-1]; // total volume of solids
312
Real Vtot=Vs/(1-cloudPorosity); // total volume of cell containing the packing
313
Vector3r cellSize=pow(Vtot,1/3.)*Vector3r().Ones(); // for now, assume always cubic sample
314
Real rMean=pow(Vs/(numSph*(4/3.)*Mathr::PI),1/3.); // make rMean such that particleSD will compute the right Vs (called a bit confusingly Vtot anyway) inversely
315
// cerr<<"Vs="<<Vs<<", Vtot="<<Vtot<<", rMean="<<rMean<<endl;
316
// cerr<<"cellSize="<<cellSize<<", rMean="<<rMean<<", numSph="<<numSph<<endl;
317
return particleSD(Vector3r::Zero(),cellSize,rMean,periodic,"",numSph,radii,passing,false);
320
// TODO: header, python wrapper, default params
322
// Discrete particle size distribution
323
long SpherePack::particleSD(Vector3r mn, Vector3r mx, Real rMean, bool periodic, string name, int numSph, const vector<Real>& radii, const vector<Real>& passing, bool passingIsNotPercentageButCount, int seed){
324
vector<Real> numbers;
325
if(!passingIsNotPercentageButCount){
326
Real Vtot=numSph*4./3.*Mathr::PI*pow(rMean,3.); // total volume of the packing (computed with rMean)
328
// calculate number of spheres necessary per each radius to match the wanted psd
329
// passing has to contain increasing values
330
for (size_t i=0; i<radii.size(); i++){
331
Real volS=4./3.*Mathr::PI*pow(radii[i],3.);
332
if (i==0) {numbers.push_back(passing[i]/100.*Vtot/volS);}
333
else {numbers.push_back((passing[i]-passing[i-1])/100.*Vtot/volS);} //
334
cout<<"fraction #"<<i<<" ("<<passing[i]<<"%, r="<<radii[i]<<"): "<<numbers[i]<<" spheres, fraction/cloud volumes "<<volS<<"/"<<Vtot<<endl;
337
FOREACH(Real p, passing) numbers.push_back(p);
340
static boost::minstd_rand randGen(seed!=0?seed:(int)TimingInfo::getNow(true));
341
static boost::variate_generator<boost::minstd_rand&, boost::uniform_real<> > rnd(randGen, boost::uniform_real<>(0,1));
343
const int maxTry=1000;
345
if(periodic)(cellSize=size);
346
for (int ii=(int)radii.size()-1; ii>=0; ii--){
347
Real r=radii[ii]; // select radius
348
for(int i=0; i<numbers[ii]; i++) { // place as many spheres as required by the psd for the selected radius into the free spot
350
for(t=0; t<maxTry; ++t){
352
if(!periodic) { for(int axis=0; axis<3; axis++) c[axis]=mn[axis]+r+(size[axis]-2*r)*rnd(); }
353
else { for(int axis=0; axis<3; axis++) c[axis]=mn[axis]+size[axis]*rnd(); }
354
size_t packSize=pack.size(); bool overlap=false;
356
for(size_t j=0; j<packSize; j++){ if(pow(pack[j].r+r,2) >= (pack[j].c-c).squaredNorm()) { overlap=true; break; } }
358
for(size_t j=0; j<packSize; j++){
360
for(int axis=0; axis<3; axis++) dr[axis]=min(cellWrapRel(c[axis],pack[j].c[axis],pack[j].c[axis]+size[axis]),cellWrapRel(pack[j].c[axis],c[axis],c[axis]+size[axis]));
361
if(pow(pack[j].r+r,2)>= dr.squaredNorm()){ overlap=true; break; }
364
if(!overlap) { pack.push_back(Sph(c,r)); break; }
367
if(numbers[ii]>0) LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<numbers[ii]<<" with radius "<<radii[ii]);
376
long SpherePack::particleSD_2d(Vector2r mn, Vector2r mx, Real rMean, bool periodic, string name, int numSph, const vector<Real>& radii, const vector<Real>& passing, bool passingIsNotPercentageButCount, int seed){
377
vector<Real> numbers;
378
if(!passingIsNotPercentageButCount){
379
Real Vtot=numSph*4./3.*Mathr::PI*pow(rMean,3.); // total volume of the packing (computed with rMean)
381
// calculate number of spheres necessary per each radius to match the wanted psd
382
// passing has to contain increasing values
383
for (size_t i=0; i<radii.size(); i++){
384
Real volS=4./3.*Mathr::PI*pow(radii[i],3.);
385
if (i==0) {numbers.push_back(passing[i]/100.*Vtot/volS);}
386
else {numbers.push_back((passing[i]-passing[i-1])/100.*Vtot/volS);} //
387
cout<<"fraction #"<<i<<" ("<<passing[i]<<"%, r="<<radii[i]<<"): "<<numbers[i]<<" spheres, fraction/cloud volumes "<<volS<<"/"<<Vtot<<endl;
390
FOREACH(Real p, passing) numbers.push_back(p);
393
static boost::minstd_rand randGen(seed!=0?seed:(int)TimingInfo::getNow(true));
394
static boost::variate_generator<boost::minstd_rand&, boost::uniform_real<> > rnd(randGen, boost::uniform_real<>(0,1));
396
const int maxTry=1000;
398
//if(periodic)(cellSize=size); in this case, it must be defined in py script as cell needs the third dimension
399
for (int ii=(int)radii.size()-1; ii>=0; ii--){
400
Real r=radii[ii]; // select radius
401
for(int i=0; i<numbers[ii]; i++) { // place as many spheres as required by the psd for the selected radius into the free spot
403
for(t=0; t<maxTry; ++t){
405
if(!periodic) { for(int axis=0; axis<2; axis++) c[axis]=mn[axis]+r+(size[axis]-2*r)*rnd(); }
406
else { for(int axis=0; axis<2; axis++) c[axis]=mn[axis]+size[axis]*rnd(); }
407
size_t packSize=pack.size(); bool overlap=false;
409
for(size_t j=0; j<packSize; j++){ if(pow(pack[j].r+r,2) >= (pack[j].c-c).squaredNorm()) { overlap=true; break; } }
411
for(size_t j=0; j<packSize; j++){
413
for(int axis=0; axis<2; axis++) dr[axis]=min(cellWrapRel(c[axis],pack[j].c[axis],pack[j].c[axis]+size[axis]),cellWrapRel(pack[j].c[axis],c[axis],c[axis]+size[axis]));
414
if(pow(pack[j].r+r,2)>= dr.squaredNorm()){ overlap=true; break; }
417
if(!overlap) { pack.push_back(Sph(c,r)); break; }
420
if(numbers[ii]>0) LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<numbers[ii]<<" with radius "<<radii[ii]);
428
long SpherePack::makeClumpCloud(const Vector3r& mn, const Vector3r& mx, const vector<shared_ptr<SpherePack> >& _clumps, bool periodic, int num){
429
// recenter given clumps and compute their margins
430
vector<SpherePack> clumps; /* vector<Vector3r> margins; */ Vector3r boxMargins(Vector3r::Zero()); Real maxR=0;
431
vector<Real> boundRad; // squared radii of bounding sphere for each clump
432
FOREACH(const shared_ptr<SpherePack>& c, _clumps){
434
c2.translate(c2.midPt()); //recenter
435
clumps.push_back(c2);
437
FOREACH(const Sph& s, c2.pack) r=max(r,s.c.norm()+s.r);
438
boundRad.push_back(r);
439
Vector3r cMn,cMx; c2.aabb(cMn,cMx); // centered at zero now, this gives margin
440
//margins.push_back(periodic?cMx:Vector3r::Zero());
441
//boxMargins=boxMargins.cwise().max(cMx);
442
FOREACH(const Sph& s, c2.pack) maxR=max(maxR,s.r); // keep track of maximum sphere radius
444
std::list<ClumpInfo> clumpInfos;
446
if(periodic)(cellSize=size);
447
const int maxTry=200;
448
int nGen=0; // number of clumps generated
449
// random point coordinate generator, with non-zero margins if aperiodic
450
static boost::minstd_rand randGen(TimingInfo::getNow(true));
451
typedef boost::variate_generator<boost::minstd_rand&, boost::uniform_real<> > UniRandGen;
452
static UniRandGen rndX(randGen,boost::uniform_real<>(mn[0],mx[0]));
453
static UniRandGen rndY(randGen,boost::uniform_real<>(mn[1],mx[1]));
454
static UniRandGen rndZ(randGen,boost::uniform_real<>(mn[2],mx[2]));
455
static UniRandGen rndUnit(randGen,boost::uniform_real<>(0,1));
456
while(nGen<num || num<0){
457
int clumpChoice=rand()%clumps.size();
459
while(true){ // check for tries at the end
460
Vector3r pos(rndX(),rndY(),rndZ()); // random point
461
// TODO: check this random orientation is homogeneously distributed
462
Quaternionr ori(rndUnit(),rndUnit(),rndUnit(),rndUnit()); ori.normalize();
463
// copy the packing and rotate
464
SpherePack C(clumps[clumpChoice]); C.rotateAroundOrigin(ori); C.translate(pos);
465
const Real& rad(boundRad[clumpChoice]);
466
ClumpInfo ci; // to be used later, but must be here because of goto's
469
// check against bounding spheres of other clumps, and only check individual spheres if there is overlap
471
// check overlap with box margins first
472
if((pos+rad*Vector3r::Ones()).cwise().max(mx)!=mx || (pos-rad*Vector3r::Ones()).cwise().min(mn)!=mn){ FOREACH(const Sph& s, C.pack) if((s.c+s.r*Vector3r::Ones()).cwise().max(mx)!=mx || (s.c-s.r*Vector3r::Ones()).cwise().min(mn)!=mn) goto overlap; }
473
// check overlaps with bounding spheres of other clumps
474
FOREACH(const ClumpInfo& cInfo, clumpInfos){
475
bool detailedCheck=false;
476
// check overlaps between individual spheres and bounding sphere of the other clump
477
if((pos-cInfo.center).squaredNorm()<pow(rad+cInfo.rad,2)){ FOREACH(const Sph& s, C.pack) if(pow(s.r+cInfo.rad,2)>(s.c-cInfo.center).squaredNorm()){ detailedCheck=true; break; }}
478
// check sphere-by-sphere, since bounding spheres did overlap
479
if(detailedCheck){ FOREACH(const Sph& s, C.pack) for(int id=cInfo.minId; id<=cInfo.maxId; id++) if((s.c-pack[id].c).squaredNorm()<pow(s.r+pack[id].r,2)) goto overlap;}
482
FOREACH(const ClumpInfo& cInfo, clumpInfos){
483
// bounding spheres overlap (in the periodic space)
484
if(periPtDistSq(pos,cInfo.center)<pow(rad+cInfo.rad,2)){
485
bool detailedCheck=false;
486
// check spheres with bounding sphere of the other clump
487
FOREACH(const Sph& s, C.pack) if(pow(s.r+cInfo.rad,2)>periPtDistSq(s.c,cInfo.center)){ detailedCheck=true; break; }
488
// check sphere-by-sphere
489
if(detailedCheck){ FOREACH(const Sph& s, C.pack) for(int id=cInfo.minId; id<=cInfo.maxId; id++) if(periPtDistSq(s.c,pack[id].c)<pow(s.r+pack[id].r,2)) goto overlap; }
495
// crude algorithm: check all spheres against all other spheres (slow!!)
496
// use vtkPointLocator, add all existing points and check distance of r+maxRadius, then refine
497
// for periodicity, duplicate points close than boxMargins to the respective boundary
499
for(size_t i=0; i<C.pack.size(); i++){
500
for(size_t j=0; j<pack.size(); j++){
501
const Vector3r& c(C.pack[i].c); const Real& r(C.pack[i].r);
502
if(pow(r+pack[j].r,2)>=(c-pack[j].c).squaredNorm()) goto overlap;
503
// check that we are not over the box boundary
504
// this could be handled by adjusting the position random interval (by taking off the smallest radius in the clump)
505
// but usually the margin band is relatively small and this does not make the code as hairy
506
if((c+r*Vector3r::Ones()).cwise().max(mx)!=mx || (c-r*Vector3r::Ones()).cwise().min(mn)!=mn) goto overlap;
510
for(size_t i=0; i<C.pack.size(); i++){
511
for(size_t j=0; j<pack.size(); j++){
512
const Vector3r& c(C.pack[i].c); const Real& r(C.pack[i].r);
514
for(int axis=0; axis<3; axis++) dr[axis]=min(cellWrapRel(c[axis],pack[j].c[axis],pack[j].c[axis]+size[axis]),cellWrapRel(pack[j].c[axis],c[axis],c[axis]+size[axis]));
515
if(pow(pack[j].r+r,2)>= dr.squaredNorm()) goto overlap;
521
// add the clump, if no collisions
522
/*number clumps consecutively*/ ci.clumpId=nGen; ci.center=pos; ci.rad=rad; ci.minId=pack.size(); ci.maxId=pack.size()+C.pack.size();
523
FOREACH(const Sph& s, C.pack){
524
pack.push_back(Sph(s.c,s.r,ci.clumpId));
526
clumpInfos.push_back(ci);
529
break; // break away from the try-loop
533
if(tries++==maxTry){ // last loop
534
if(num>0) LOG_WARN("Exceeded "<<maxTry<<" attempts to place non-overlapping clump. Only "<<nGen<<" clumps were added, although you requested "<<num);
542
bool SpherePack::hasClumps() const { FOREACH(const Sph& s, pack){ if(s.clumpId>=0) return true; } return false; }
543
python::tuple SpherePack::getClumps() const{
544
std::map<int,py::list> clumps;
545
py::list standalone; size_t packSize=pack.size();
546
for(size_t i=0; i<packSize; i++){
547
const Sph& s(pack[i]);
548
if(s.clumpId<0) { standalone.append(i); continue; }
549
if(clumps.count(s.clumpId)==0) clumps[s.clumpId]=py::list();
550
clumps[s.clumpId].append(i);
553
typedef std::pair<int,py::list> intListPair;
554
FOREACH(const intListPair& c, clumps) clumpList.append(c.second);
555
return python::make_tuple(standalone,clumpList);
added
removed
pkg/dem/SpherePack.cpp
