Cleanup: remove redundant doxygen \file argument
[blender.git] / source / blender / blenkernel / intern / particle_distribute.c
1 /*
2  * This program is free software; you can redistribute it and/or
3  * modify it under the terms of the GNU General Public License
4  * as published by the Free Software Foundation; either version 2
5  * of the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program; if not, write to the Free Software Foundation,
14  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15  *
16  * The Original Code is Copyright (C) 2007 by Janne Karhu.
17  * All rights reserved.
18  */
19
20 /** \file \ingroup bke
21  */
22
23 #include <string.h>
24
25 #include "MEM_guardedalloc.h"
26
27 #include "BLI_utildefines.h"
28 #include "BLI_jitter_2d.h"
29 #include "BLI_kdtree.h"
30 #include "BLI_math.h"
31 #include "BLI_math_geom.h"
32 #include "BLI_rand.h"
33 #include "BLI_sort.h"
34 #include "BLI_task.h"
35
36 #include "DNA_mesh_types.h"
37 #include "DNA_meshdata_types.h"
38 #include "DNA_modifier_types.h"
39 #include "DNA_particle_types.h"
40 #include "DNA_scene_types.h"
41
42 #include "BKE_global.h"
43 #include "BKE_library.h"
44 #include "BKE_mesh.h"
45 #include "BKE_object.h"
46 #include "BKE_particle.h"
47
48 #include "DEG_depsgraph_query.h"
49
50 static void alloc_child_particles(ParticleSystem *psys, int tot)
51 {
52         if (psys->child) {
53                 /* only re-allocate if we have to */
54                 if (psys->part->childtype && psys->totchild == tot) {
55                         memset(psys->child, 0, tot*sizeof(ChildParticle));
56                         return;
57                 }
58
59                 MEM_freeN(psys->child);
60                 psys->child=NULL;
61                 psys->totchild=0;
62         }
63
64         if (psys->part->childtype) {
65                 psys->totchild= tot;
66                 if (psys->totchild)
67                         psys->child= MEM_callocN(psys->totchild*sizeof(ChildParticle), "child_particles");
68         }
69 }
70
71 static void distribute_simple_children(Scene *scene, Object *ob, Mesh *final_mesh, Mesh *deform_mesh, ParticleSystem *psys, const bool use_render_params)
72 {
73         ChildParticle *cpa = NULL;
74         int i, p;
75         int child_nbr= psys_get_child_number(scene, psys, use_render_params);
76         int totpart= psys_get_tot_child(scene, psys, use_render_params);
77         RNG *rng = BLI_rng_new_srandom(31415926 + psys->seed + psys->child_seed);
78
79         alloc_child_particles(psys, totpart);
80
81         cpa = psys->child;
82         for (i=0; i<child_nbr; i++) {
83                 for (p=0; p<psys->totpart; p++,cpa++) {
84                         float length=2.0;
85                         cpa->parent=p;
86
87                         /* create even spherical distribution inside unit sphere */
88                         while (length>=1.0f) {
89                                 cpa->fuv[0]=2.0f*BLI_rng_get_float(rng)-1.0f;
90                                 cpa->fuv[1]=2.0f*BLI_rng_get_float(rng)-1.0f;
91                                 cpa->fuv[2]=2.0f*BLI_rng_get_float(rng)-1.0f;
92                                 length=len_v3(cpa->fuv);
93                         }
94
95                         cpa->num=-1;
96                 }
97         }
98         /* dmcache must be updated for parent particles if children from faces is used */
99         psys_calc_dmcache(ob, final_mesh, deform_mesh, psys);
100
101         BLI_rng_free(rng);
102 }
103 static void distribute_grid(Mesh *mesh, ParticleSystem *psys)
104 {
105         ParticleData *pa=NULL;
106         float min[3], max[3], delta[3], d;
107         MVert *mv, *mvert = mesh->mvert;
108         int totvert=mesh->totvert, from=psys->part->from;
109         int i, j, k, p, res=psys->part->grid_res, size[3], axis;
110
111         /* find bounding box of dm */
112         if (totvert > 0) {
113                 mv=mvert;
114                 copy_v3_v3(min, mv->co);
115                 copy_v3_v3(max, mv->co);
116                 mv++;
117                 for (i = 1; i < totvert; i++, mv++) {
118                         minmax_v3v3_v3(min, max, mv->co);
119                 }
120         }
121         else {
122                 zero_v3(min);
123                 zero_v3(max);
124         }
125
126         sub_v3_v3v3(delta, max, min);
127
128         /* determine major axis */
129         axis = axis_dominant_v3_single(delta);
130
131         d = delta[axis]/(float)res;
132
133         size[axis] = res;
134         size[(axis+1)%3] = (int)ceil(delta[(axis+1)%3]/d);
135         size[(axis+2)%3] = (int)ceil(delta[(axis+2)%3]/d);
136
137         /* float errors grrr.. */
138         size[(axis+1)%3] = MIN2(size[(axis+1)%3],res);
139         size[(axis+2)%3] = MIN2(size[(axis+2)%3],res);
140
141         size[0] = MAX2(size[0], 1);
142         size[1] = MAX2(size[1], 1);
143         size[2] = MAX2(size[2], 1);
144
145         /* no full offset for flat/thin objects */
146         min[0]+= d < delta[0] ? d/2.f : delta[0]/2.f;
147         min[1]+= d < delta[1] ? d/2.f : delta[1]/2.f;
148         min[2]+= d < delta[2] ? d/2.f : delta[2]/2.f;
149
150         for (i=0,p=0,pa=psys->particles; i<res; i++) {
151                 for (j=0; j<res; j++) {
152                         for (k=0; k<res; k++,p++,pa++) {
153                                 pa->fuv[0] = min[0] + (float)i*d;
154                                 pa->fuv[1] = min[1] + (float)j*d;
155                                 pa->fuv[2] = min[2] + (float)k*d;
156                                 pa->flag |= PARS_UNEXIST;
157                                 pa->hair_index = 0; /* abused in volume calculation */
158                         }
159                 }
160         }
161
162         /* enable particles near verts/edges/faces/inside surface */
163         if (from==PART_FROM_VERT) {
164                 float vec[3];
165
166                 pa=psys->particles;
167
168                 min[0] -= d/2.0f;
169                 min[1] -= d/2.0f;
170                 min[2] -= d/2.0f;
171
172                 for (i=0,mv=mvert; i<totvert; i++,mv++) {
173                         sub_v3_v3v3(vec,mv->co,min);
174                         vec[0]/=delta[0];
175                         vec[1]/=delta[1];
176                         vec[2]/=delta[2];
177                         pa[((int)(vec[0] * (size[0] - 1))  * res +
178                             (int)(vec[1] * (size[1] - 1))) * res +
179                             (int)(vec[2] * (size[2] - 1))].flag &= ~PARS_UNEXIST;
180                 }
181         }
182         else if (ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
183                 float co1[3], co2[3];
184
185                 MFace *mface= NULL, *mface_array;
186                 float v1[3], v2[3], v3[3], v4[4], lambda;
187                 int a, a1, a2, a0mul, a1mul, a2mul, totface;
188                 int amax= from==PART_FROM_FACE ? 3 : 1;
189
190                 totface = mesh->totface;
191                 mface = mface_array = mesh->mface;
192
193                 for (a=0; a<amax; a++) {
194                         if (a==0) { a0mul=res*res; a1mul=res; a2mul=1; }
195                         else if (a==1) { a0mul=res; a1mul=1; a2mul=res*res; }
196                         else { a0mul=1; a1mul=res*res; a2mul=res; }
197
198                         for (a1=0; a1<size[(a+1)%3]; a1++) {
199                                 for (a2=0; a2<size[(a+2)%3]; a2++) {
200                                         mface= mface_array;
201
202                                         pa = psys->particles + a1*a1mul + a2*a2mul;
203                                         copy_v3_v3(co1, pa->fuv);
204                                         co1[a] -= d < delta[a] ? d/2.f : delta[a]/2.f;
205                                         copy_v3_v3(co2, co1);
206                                         co2[a] += delta[a] + 0.001f*d;
207                                         co1[a] -= 0.001f*d;
208
209                                         struct IsectRayPrecalc isect_precalc;
210                                         float ray_direction[3];
211                                         sub_v3_v3v3(ray_direction, co2, co1);
212                                         isect_ray_tri_watertight_v3_precalc(&isect_precalc, ray_direction);
213
214                                         /* lets intersect the faces */
215                                         for (i=0; i<totface; i++,mface++) {
216                                                 copy_v3_v3(v1, mvert[mface->v1].co);
217                                                 copy_v3_v3(v2, mvert[mface->v2].co);
218                                                 copy_v3_v3(v3, mvert[mface->v3].co);
219
220                                                 bool intersects_tri = isect_ray_tri_watertight_v3(co1,
221                                                                                                   &isect_precalc,
222                                                                                                   v1, v2, v3,
223                                                                                                   &lambda, NULL);
224                                                 if (intersects_tri) {
225                                                         if (from==PART_FROM_FACE)
226                                                                 (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
227                                                         else /* store number of intersections */
228                                                                 (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
229                                                 }
230
231                                                 if (mface->v4 && (!intersects_tri || from==PART_FROM_VOLUME)) {
232                                                         copy_v3_v3(v4, mvert[mface->v4].co);
233
234                                                         if (isect_ray_tri_watertight_v3(
235                                                                     co1,
236                                                                     &isect_precalc,
237                                                                     v1, v3, v4,
238                                                                     &lambda, NULL))
239                                                         {
240                                                                 if (from==PART_FROM_FACE)
241                                                                         (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
242                                                                 else
243                                                                         (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
244                                                         }
245                                                 }
246                                         }
247
248                                         if (from==PART_FROM_VOLUME) {
249                                                 int in=pa->hair_index%2;
250                                                 if (in) pa->hair_index++;
251                                                 for (i=0; i<size[0]; i++) {
252                                                         if (in || (pa+i*a0mul)->hair_index%2)
253                                                                 (pa+i*a0mul)->flag &= ~PARS_UNEXIST;
254                                                         /* odd intersections == in->out / out->in */
255                                                         /* even intersections -> in stays same */
256                                                         in=(in + (pa+i*a0mul)->hair_index) % 2;
257                                                 }
258                                         }
259                                 }
260                         }
261                 }
262         }
263
264         if (psys->part->flag & PART_GRID_HEXAGONAL) {
265                 for (i=0,p=0,pa=psys->particles; i<res; i++) {
266                         for (j=0; j<res; j++) {
267                                 for (k=0; k<res; k++,p++,pa++) {
268                                         if (j%2)
269                                                 pa->fuv[0] += d/2.f;
270
271                                         if (k%2) {
272                                                 pa->fuv[0] += d/2.f;
273                                                 pa->fuv[1] += d/2.f;
274                                         }
275                                 }
276                         }
277                 }
278         }
279
280         if (psys->part->flag & PART_GRID_INVERT) {
281                 for (i=0; i<size[0]; i++) {
282                         for (j=0; j<size[1]; j++) {
283                                 pa=psys->particles + res*(i*res + j);
284                                 for (k=0; k<size[2]; k++, pa++) {
285                                         pa->flag ^= PARS_UNEXIST;
286                                 }
287                         }
288                 }
289         }
290
291         if (psys->part->grid_rand > 0.f) {
292                 float rfac = d * psys->part->grid_rand;
293                 for (p=0,pa=psys->particles; p<psys->totpart; p++,pa++) {
294                         if (pa->flag & PARS_UNEXIST)
295                                 continue;
296
297                         pa->fuv[0] += rfac * (psys_frand(psys, p + 31) - 0.5f);
298                         pa->fuv[1] += rfac * (psys_frand(psys, p + 32) - 0.5f);
299                         pa->fuv[2] += rfac * (psys_frand(psys, p + 33) - 0.5f);
300                 }
301         }
302 }
303
304 /* modified copy from rayshade.c */
305 static void hammersley_create(float *out, int n, int seed, float amount)
306 {
307         RNG *rng;
308
309         double offs[2], t;
310
311         rng = BLI_rng_new(31415926 + n + seed);
312         offs[0] = BLI_rng_get_double(rng) + (double)amount;
313         offs[1] = BLI_rng_get_double(rng) + (double)amount;
314         BLI_rng_free(rng);
315
316         for (int k = 0; k < n; k++) {
317                 BLI_hammersley_1D(k, &t);
318
319                 out[2*k + 0] = fmod((double)k/(double)n + offs[0], 1.0);
320                 out[2*k + 1] = fmod(t + offs[1], 1.0);
321         }
322 }
323
324 /* almost exact copy of BLI_jitter_init */
325 static void init_mv_jit(float *jit, int num, int seed2, float amount)
326 {
327         RNG *rng;
328         float *jit2, x, rad1, rad2, rad3;
329         int i, num2;
330
331         if (num==0) return;
332
333         rad1= (float)(1.0f/sqrtf((float)num));
334         rad2= (float)(1.0f/((float)num));
335         rad3= (float)sqrtf((float)num)/((float)num);
336
337         rng = BLI_rng_new(31415926 + num + seed2);
338         x= 0;
339         num2 = 2 * num;
340         for (i=0; i<num2; i+=2) {
341
342                 jit[i] = x + amount*rad1*(0.5f - BLI_rng_get_float(rng));
343                 jit[i+1] = i/(2.0f*num) + amount*rad1*(0.5f - BLI_rng_get_float(rng));
344
345                 jit[i]-= (float)floor(jit[i]);
346                 jit[i+1]-= (float)floor(jit[i+1]);
347
348                 x+= rad3;
349                 x -= (float)floor(x);
350         }
351
352         jit2= MEM_mallocN(12 + 2*sizeof(float)*num, "initjit");
353
354         for (i=0 ; i<4 ; i++) {
355                 BLI_jitterate1((float (*)[2])jit, (float (*)[2])jit2, num, rad1);
356                 BLI_jitterate1((float (*)[2])jit, (float (*)[2])jit2, num, rad1);
357                 BLI_jitterate2((float (*)[2])jit, (float (*)[2])jit2, num, rad2);
358         }
359         MEM_freeN(jit2);
360         BLI_rng_free(rng);
361 }
362
363 static void psys_uv_to_w(float u, float v, int quad, float *w)
364 {
365         float vert[4][3], co[3];
366
367         if (!quad) {
368                 if (u+v > 1.0f)
369                         v= 1.0f-v;
370                 else
371                         u= 1.0f-u;
372         }
373
374         vert[0][0] = 0.0f; vert[0][1] = 0.0f; vert[0][2] = 0.0f;
375         vert[1][0] = 1.0f; vert[1][1] = 0.0f; vert[1][2] = 0.0f;
376         vert[2][0] = 1.0f; vert[2][1] = 1.0f; vert[2][2] = 0.0f;
377
378         co[0] = u;
379         co[1] = v;
380         co[2] = 0.0f;
381
382         if (quad) {
383                 vert[3][0] = 0.0f; vert[3][1] = 1.0f; vert[3][2] = 0.0f;
384                 interp_weights_poly_v3( w,vert, 4, co);
385         }
386         else {
387                 interp_weights_poly_v3( w,vert, 3, co);
388                 w[3] = 0.0f;
389         }
390 }
391
392 /* Find the index in "sum" array before "value" is crossed. */
393 static int distribute_binary_search(float *sum, int n, float value)
394 {
395         int mid, low = 0, high = n - 1;
396
397         if (high == low)
398                 return low;
399
400         if (sum[low] >= value)
401                 return low;
402
403         if (sum[high - 1] < value)
404                 return high;
405
406         while (low < high) {
407                 mid = (low + high) / 2;
408
409                 if ((sum[mid] >= value) && (sum[mid - 1] < value))
410                         return mid;
411
412                 if (sum[mid] > value) {
413                         high = mid - 1;
414                 }
415                 else {
416                         low = mid + 1;
417                 }
418         }
419
420         return low;
421 }
422
423 /* the max number if calls to rng_* funcs within psys_thread_distribute_particle
424  * be sure to keep up to date if this changes */
425 #define PSYS_RND_DIST_SKIP 3
426
427 /* note: this function must be thread safe, for from == PART_FROM_CHILD */
428 #define ONLY_WORKING_WITH_PA_VERTS 0
429 static void distribute_from_verts_exec(ParticleTask *thread, ParticleData *pa, int p)
430 {
431         ParticleThreadContext *ctx= thread->ctx;
432         MFace *mface;
433
434         mface = ctx->mesh->mface;
435
436         int rng_skip_tot = PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
437
438         /* TODO_PARTICLE - use original index */
439         pa->num = ctx->index[p];
440
441         zero_v4(pa->fuv);
442
443         if (pa->num != DMCACHE_NOTFOUND && pa->num < ctx->mesh->totvert) {
444
445                 /* This finds the first face to contain the emitting vertex,
446                  * this is not ideal, but is mostly fine as UV seams generally
447                  * map to equal-colored parts of a texture */
448                 for (int i = 0; i < ctx->mesh->totface; i++, mface++) {
449                         if (ELEM(pa->num, mface->v1, mface->v2, mface->v3, mface->v4)) {
450                                 unsigned int *vert = &mface->v1;
451
452                                 for (int j = 0; j < 4; j++, vert++) {
453                                         if (*vert == pa->num) {
454                                                 pa->fuv[j] = 1.0f;
455                                                 break;
456                                         }
457                                 }
458
459                                 break;
460                         }
461                 }
462         }
463
464 #if ONLY_WORKING_WITH_PA_VERTS
465         if (ctx->tree) {
466                 KDTreeNearest ptn[3];
467                 int w, maxw;
468
469                 psys_particle_on_dm(ctx->mesh,from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co1,0,0,0,orco1,0);
470                 BKE_mesh_orco_verts_transform(ob->data, &orco1, 1, 1);
471                 maxw = BLI_kdtree_find_nearest_n(ctx->tree,orco1,ptn,3);
472
473                 for (w=0; w<maxw; w++) {
474                         pa->verts[w]=ptn->num;
475                 }
476         }
477 #endif
478
479         BLI_assert(rng_skip_tot >= 0);  /* should never be below zero */
480         if (rng_skip_tot > 0) {
481                 BLI_rng_skip(thread->rng, rng_skip_tot);
482         }
483 }
484
485 static void distribute_from_faces_exec(ParticleTask *thread, ParticleData *pa, int p) {
486         ParticleThreadContext *ctx= thread->ctx;
487         Mesh *mesh = ctx->mesh;
488         float randu, randv;
489         int distr= ctx->distr;
490         int i;
491         int rng_skip_tot= PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
492
493         MFace *mface;
494
495         pa->num = i = ctx->index[p];
496         mface = &mesh->mface[i];
497
498         switch (distr) {
499                 case PART_DISTR_JIT:
500                         if (ctx->jitlevel == 1) {
501                                 if (mface->v4)
502                                         psys_uv_to_w(0.5f, 0.5f, mface->v4, pa->fuv);
503                                 else
504                                         psys_uv_to_w(1.0f / 3.0f, 1.0f / 3.0f, mface->v4, pa->fuv);
505                         }
506                         else {
507                                 float offset = fmod(ctx->jitoff[i] + (float)p, (float)ctx->jitlevel);
508                                 if (!isnan(offset)) {
509                                         psys_uv_to_w(ctx->jit[2*(int)offset], ctx->jit[2*(int)offset+1], mface->v4, pa->fuv);
510                                 }
511                         }
512                         break;
513                 case PART_DISTR_RAND:
514                         randu= BLI_rng_get_float(thread->rng);
515                         randv= BLI_rng_get_float(thread->rng);
516                         rng_skip_tot -= 2;
517
518                         psys_uv_to_w(randu, randv, mface->v4, pa->fuv);
519                         break;
520         }
521         pa->foffset= 0.0f;
522
523         BLI_assert(rng_skip_tot >= 0);  /* should never be below zero */
524         if (rng_skip_tot > 0) {
525                 BLI_rng_skip(thread->rng, rng_skip_tot);
526         }
527 }
528
529 static void distribute_from_volume_exec(ParticleTask *thread, ParticleData *pa, int p) {
530         ParticleThreadContext *ctx= thread->ctx;
531         Mesh *mesh = ctx->mesh;
532         float *v1, *v2, *v3, *v4, nor[3], co[3];
533         float cur_d, min_d, randu, randv;
534         int distr= ctx->distr;
535         int i, intersect, tot;
536         int rng_skip_tot= PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
537
538         MFace *mface;
539         MVert *mvert = mesh->mvert;
540
541         pa->num = i = ctx->index[p];
542         mface = &mesh->mface[i];
543
544         switch (distr) {
545                 case PART_DISTR_JIT:
546                         if (ctx->jitlevel == 1) {
547                                 if (mface->v4)
548                                         psys_uv_to_w(0.5f, 0.5f, mface->v4, pa->fuv);
549                                 else
550                                         psys_uv_to_w(1.0f / 3.0f, 1.0f / 3.0f, mface->v4, pa->fuv);
551                         }
552                         else {
553                                 float offset = fmod(ctx->jitoff[i] + (float)p, (float)ctx->jitlevel);
554                                 if (!isnan(offset)) {
555                                         psys_uv_to_w(ctx->jit[2*(int)offset], ctx->jit[2*(int)offset+1], mface->v4, pa->fuv);
556                                 }
557                         }
558                         break;
559                 case PART_DISTR_RAND:
560                         randu= BLI_rng_get_float(thread->rng);
561                         randv= BLI_rng_get_float(thread->rng);
562                         rng_skip_tot -= 2;
563
564                         psys_uv_to_w(randu, randv, mface->v4, pa->fuv);
565                         break;
566         }
567         pa->foffset= 0.0f;
568
569         /* experimental */
570         tot = mesh->totface;
571
572         psys_interpolate_face(mvert,mface,0,0,pa->fuv,co,nor,0,0,0);
573
574         normalize_v3(nor);
575         negate_v3(nor);
576
577         min_d=FLT_MAX;
578         intersect=0;
579
580         for (i=0, mface=mesh->mface; i<tot; i++,mface++) {
581                 if (i==pa->num) continue;
582
583                 v1=mvert[mface->v1].co;
584                 v2=mvert[mface->v2].co;
585                 v3=mvert[mface->v3].co;
586
587                 if (isect_ray_tri_v3(co, nor, v2, v3, v1, &cur_d, NULL)) {
588                         if (cur_d<min_d) {
589                                 min_d=cur_d;
590                                 pa->foffset=cur_d*0.5f; /* to the middle of volume */
591                                 intersect=1;
592                         }
593                 }
594                 if (mface->v4) {
595                         v4=mvert[mface->v4].co;
596
597                         if (isect_ray_tri_v3(co, nor, v4, v1, v3, &cur_d, NULL)) {
598                                 if (cur_d<min_d) {
599                                         min_d=cur_d;
600                                         pa->foffset=cur_d*0.5f; /* to the middle of volume */
601                                         intersect=1;
602                                 }
603                         }
604                 }
605         }
606         if (intersect==0)
607                 pa->foffset=0.0;
608         else {
609                 switch (distr) {
610                         case PART_DISTR_JIT:
611                                 pa->foffset *= ctx->jit[p % (2 * ctx->jitlevel)];
612                                 break;
613                         case PART_DISTR_RAND:
614                                 pa->foffset *= BLI_rng_get_float(thread->rng);
615                                 rng_skip_tot--;
616                                 break;
617                 }
618         }
619
620         BLI_assert(rng_skip_tot >= 0); /* should never be below zero */
621         if (rng_skip_tot > 0) {
622                 BLI_rng_skip(thread->rng, rng_skip_tot);
623         }
624 }
625
626 static void distribute_children_exec(ParticleTask *thread, ChildParticle *cpa, int p) {
627         ParticleThreadContext *ctx= thread->ctx;
628         Object *ob= ctx->sim.ob;
629         Mesh *mesh = ctx->mesh;
630         float orco1[3], co1[3], nor1[3];
631         float randu, randv;
632         int cfrom= ctx->cfrom;
633         int i;
634         int rng_skip_tot= PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
635
636         MFace *mf;
637
638         if (ctx->index[p] < 0) {
639                 cpa->num=0;
640                 cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]=0.0f;
641                 cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
642                 return;
643         }
644
645         mf = &mesh->mface[ctx->index[p]];
646
647         randu= BLI_rng_get_float(thread->rng);
648         randv= BLI_rng_get_float(thread->rng);
649         rng_skip_tot -= 2;
650
651         psys_uv_to_w(randu, randv, mf->v4, cpa->fuv);
652
653         cpa->num = ctx->index[p];
654
655         if (ctx->tree) {
656                 KDTreeNearest ptn[10];
657                 int w,maxw;//, do_seams;
658                 float maxd /*, mind,dd */, totw= 0.0f;
659                 int parent[10];
660                 float pweight[10];
661
662                 psys_particle_on_dm(mesh,cfrom,cpa->num,DMCACHE_ISCHILD,cpa->fuv,cpa->foffset,co1,nor1,NULL,NULL,orco1);
663                 BKE_mesh_orco_verts_transform(ob->data, &orco1, 1, 1);
664                 maxw = BLI_kdtree_find_nearest_n(ctx->tree,orco1,ptn,3);
665
666                 maxd=ptn[maxw-1].dist;
667                 /* mind=ptn[0].dist; */ /* UNUSED */
668
669                 /* the weights here could be done better */
670                 for (w=0; w<maxw; w++) {
671                         parent[w]=ptn[w].index;
672                         pweight[w]=(float)pow(2.0,(double)(-6.0f*ptn[w].dist/maxd));
673                 }
674                 for (;w<10; w++) {
675                         parent[w]=-1;
676                         pweight[w]=0.0f;
677                 }
678
679                 for (w=0,i=0; w<maxw && i<4; w++) {
680                         if (parent[w]>=0) {
681                                 cpa->pa[i]=parent[w];
682                                 cpa->w[i]=pweight[w];
683                                 totw+=pweight[w];
684                                 i++;
685                         }
686                 }
687                 for (;i<4; i++) {
688                         cpa->pa[i]=-1;
689                         cpa->w[i]=0.0f;
690                 }
691
692                 if (totw > 0.0f) {
693                         for (w = 0; w < 4; w++) {
694                                 cpa->w[w] /= totw;
695                         }
696                 }
697
698                 cpa->parent=cpa->pa[0];
699         }
700
701         if (rng_skip_tot > 0) /* should never be below zero */
702                 BLI_rng_skip(thread->rng, rng_skip_tot);
703 }
704
705 static void exec_distribute_parent(TaskPool * __restrict UNUSED(pool), void *taskdata, int UNUSED(threadid))
706 {
707         ParticleTask *task = taskdata;
708         ParticleSystem *psys= task->ctx->sim.psys;
709         ParticleData *pa;
710         int p;
711
712         BLI_rng_skip(task->rng, PSYS_RND_DIST_SKIP * task->begin);
713
714         pa= psys->particles + task->begin;
715         switch (psys->part->from) {
716                 case PART_FROM_FACE:
717                         for (p = task->begin; p < task->end; ++p, ++pa)
718                                 distribute_from_faces_exec(task, pa, p);
719                         break;
720                 case PART_FROM_VOLUME:
721                         for (p = task->begin; p < task->end; ++p, ++pa)
722                                 distribute_from_volume_exec(task, pa, p);
723                         break;
724                 case PART_FROM_VERT:
725                         for (p = task->begin; p < task->end; ++p, ++pa)
726                                 distribute_from_verts_exec(task, pa, p);
727                         break;
728         }
729 }
730
731 static void exec_distribute_child(TaskPool * __restrict UNUSED(pool), void *taskdata, int UNUSED(threadid))
732 {
733         ParticleTask *task = taskdata;
734         ParticleSystem *psys = task->ctx->sim.psys;
735         ChildParticle *cpa;
736         int p;
737
738         /* RNG skipping at the beginning */
739         cpa = psys->child;
740         for (p = 0; p < task->begin; ++p, ++cpa) {
741                 BLI_rng_skip(task->rng, PSYS_RND_DIST_SKIP);
742         }
743
744         for (; p < task->end; ++p, ++cpa) {
745                 distribute_children_exec(task, cpa, p);
746         }
747 }
748
749 static int distribute_compare_orig_index(const void *p1, const void *p2, void *user_data)
750 {
751         int *orig_index = (int *) user_data;
752         int index1 = orig_index[*(const int *)p1];
753         int index2 = orig_index[*(const int *)p2];
754
755         if (index1 < index2)
756                 return -1;
757         else if (index1 == index2) {
758                 /* this pointer comparison appears to make qsort stable for glibc,
759                  * and apparently on solaris too, makes the renders reproducible */
760                 if (p1 < p2)
761                         return -1;
762                 else if (p1 == p2)
763                         return 0;
764                 else
765                         return 1;
766         }
767         else
768                 return 1;
769 }
770
771 static void distribute_invalid(ParticleSimulationData *sim, int from)
772 {
773         Scene *scene = sim->scene;
774         ParticleSystem *psys = sim->psys;
775         const bool use_render_params = (DEG_get_mode(sim->depsgraph) == DAG_EVAL_RENDER);
776
777         if (from == PART_FROM_CHILD) {
778                 ChildParticle *cpa;
779                 int p, totchild = psys_get_tot_child(scene, psys, use_render_params);
780
781                 if (psys->child && totchild) {
782                         for (p=0,cpa=psys->child; p<totchild; p++,cpa++) {
783                                 cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3] = 0.0;
784                                 cpa->foffset= 0.0f;
785                                 cpa->parent=0;
786                                 cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
787                                 cpa->num= -1;
788                         }
789                 }
790         }
791         else {
792                 PARTICLE_P;
793                 LOOP_PARTICLES {
794                         pa->fuv[0] = pa->fuv[1] = pa->fuv[2] = pa->fuv[3] = 0.0;
795                         pa->foffset= 0.0f;
796                         pa->num= -1;
797                 }
798         }
799 }
800
801 /* Creates a distribution of coordinates on a Mesh */
802 static int psys_thread_context_init_distribute(ParticleThreadContext *ctx, ParticleSimulationData *sim, int from)
803 {
804         Scene *scene = sim->scene;
805         Mesh *final_mesh = sim->psmd->mesh_final;
806         Object *ob = sim->ob;
807         ParticleSystem *psys= sim->psys;
808         ParticleData *pa=0, *tpars= 0;
809         ParticleSettings *part;
810         ParticleSeam *seams= 0;
811         KDTree *tree=0;
812         Mesh *mesh = NULL;
813         float *jit= NULL;
814         int i, p=0;
815         int cfrom=0;
816         int totelem=0, totpart, *particle_element=0, children=0, totseam=0;
817         int jitlevel= 1, distr;
818         float *element_weight=NULL,*jitter_offset=NULL, *vweight=NULL;
819         float cur, maxweight=0.0, tweight, totweight, inv_totweight, co[3], nor[3], orco[3];
820         RNG *rng = NULL;
821
822         if (ELEM(NULL, ob, psys, psys->part))
823                 return 0;
824
825         part=psys->part;
826         totpart=psys->totpart;
827         if (totpart==0)
828                 return 0;
829
830         if (!final_mesh->runtime.deformed_only && !CustomData_get_layer(&final_mesh->fdata, CD_ORIGINDEX)) {
831                 printf("Can't create particles with the current modifier stack, disable destructive modifiers\n");
832 // XXX          error("Can't paint with the current modifier stack, disable destructive modifiers");
833                 return 0;
834         }
835
836         /* XXX This distribution code is totally broken in case from == PART_FROM_CHILD, it's always using finaldm
837          *     even if use_modifier_stack is unset... But making things consistent here break all existing edited
838          *     hair systems, so better wait for complete rewrite.
839          */
840
841         psys_thread_context_init(ctx, sim);
842
843         const bool use_render_params = (DEG_get_mode(sim->depsgraph) == DAG_EVAL_RENDER);
844
845         /* First handle special cases */
846         if (from == PART_FROM_CHILD) {
847                 /* Simple children */
848                 if (part->childtype != PART_CHILD_FACES) {
849                         distribute_simple_children(scene, ob, final_mesh, sim->psmd->mesh_original, psys, use_render_params);
850                         return 0;
851                 }
852         }
853         else {
854                 /* Grid distribution */
855                 if (part->distr==PART_DISTR_GRID && from != PART_FROM_VERT) {
856                         if (psys->part->use_modifier_stack) {
857                                 mesh = final_mesh;
858                         }
859                         else {
860                                 BKE_id_copy_ex(NULL, ob->data, (ID **)&mesh, LIB_ID_COPY_LOCALIZE);
861                         }
862                         BKE_mesh_tessface_ensure(mesh);
863
864                         distribute_grid(mesh,psys);
865
866                         if (mesh != final_mesh) {
867                                 BKE_id_free(NULL, mesh);
868                         }
869
870                         return 0;
871                 }
872         }
873
874         /* Create trees and original coordinates if needed */
875         if (from == PART_FROM_CHILD) {
876                 distr = PART_DISTR_RAND;
877                 rng = BLI_rng_new_srandom(31415926 + psys->seed + psys->child_seed);
878                 mesh= final_mesh;
879
880                 /* BMESH ONLY */
881                 BKE_mesh_tessface_ensure(mesh);
882
883                 children=1;
884
885                 tree=BLI_kdtree_new(totpart);
886
887                 for (p=0,pa=psys->particles; p<totpart; p++,pa++) {
888                         psys_particle_on_dm(mesh,part->from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co,nor,0,0,orco);
889                         BKE_mesh_orco_verts_transform(ob->data, &orco, 1, 1);
890                         BLI_kdtree_insert(tree, p, orco);
891                 }
892
893                 BLI_kdtree_balance(tree);
894
895                 totpart = psys_get_tot_child(scene, psys, use_render_params);
896                 cfrom = from = PART_FROM_FACE;
897         }
898         else {
899                 distr = part->distr;
900
901                 rng = BLI_rng_new_srandom(31415926 + psys->seed);
902
903                 if (psys->part->use_modifier_stack)
904                         mesh = final_mesh;
905                 else
906                         BKE_id_copy_ex(NULL, ob->data, (ID **)&mesh, LIB_ID_COPY_LOCALIZE);
907
908                 BKE_mesh_tessface_ensure(mesh);
909
910                 /* we need orco for consistent distributions */
911                 if (!CustomData_has_layer(&mesh->vdata, CD_ORCO))
912                         CustomData_add_layer(&mesh->vdata, CD_ORCO, CD_ASSIGN, BKE_mesh_orco_verts_get(ob), mesh->totvert);
913
914                 if (from == PART_FROM_VERT) {
915                         MVert *mv = mesh->mvert;
916                         float (*orcodata)[3] = CustomData_get_layer(&mesh->vdata, CD_ORCO);
917                         int totvert = mesh->totvert;
918
919                         tree=BLI_kdtree_new(totvert);
920
921                         for (p=0; p<totvert; p++) {
922                                 if (orcodata) {
923                                         copy_v3_v3(co,orcodata[p]);
924                                         BKE_mesh_orco_verts_transform(ob->data, &co, 1, 1);
925                                 }
926                                 else
927                                         copy_v3_v3(co,mv[p].co);
928                                 BLI_kdtree_insert(tree, p, co);
929                         }
930
931                         BLI_kdtree_balance(tree);
932                 }
933         }
934
935         /* Get total number of emission elements and allocate needed arrays */
936         totelem = (from == PART_FROM_VERT) ? mesh->totvert : mesh->totface;
937
938         if (totelem == 0) {
939                 distribute_invalid(sim, children ? PART_FROM_CHILD : 0);
940
941                 if (G.debug & G_DEBUG)
942                         fprintf(stderr,"Particle distribution error: Nothing to emit from!\n");
943
944                 if (mesh != final_mesh) BKE_id_free(NULL, mesh);
945
946                 BLI_kdtree_free(tree);
947                 BLI_rng_free(rng);
948
949                 return 0;
950         }
951
952         element_weight   = MEM_callocN(sizeof(float) * totelem, "particle_distribution_weights");
953         particle_element = MEM_callocN(sizeof(int) * totpart, "particle_distribution_indexes");
954         jitter_offset    = MEM_callocN(sizeof(float) * totelem, "particle_distribution_jitoff");
955
956         /* Calculate weights from face areas */
957         if ((part->flag & PART_EDISTR || children) && from != PART_FROM_VERT) {
958                 MVert *v1, *v2, *v3, *v4;
959                 float totarea=0.f, co1[3], co2[3], co3[3], co4[3];
960                 float (*orcodata)[3];
961
962                 orcodata = CustomData_get_layer(&mesh->vdata, CD_ORCO);
963
964                 for (i=0; i<totelem; i++) {
965                         MFace *mf = &mesh->mface[i];
966
967                         if (orcodata) {
968                                 copy_v3_v3(co1, orcodata[mf->v1]);
969                                 copy_v3_v3(co2, orcodata[mf->v2]);
970                                 copy_v3_v3(co3, orcodata[mf->v3]);
971                                 BKE_mesh_orco_verts_transform(ob->data, &co1, 1, 1);
972                                 BKE_mesh_orco_verts_transform(ob->data, &co2, 1, 1);
973                                 BKE_mesh_orco_verts_transform(ob->data, &co3, 1, 1);
974                                 if (mf->v4) {
975                                         copy_v3_v3(co4, orcodata[mf->v4]);
976                                         BKE_mesh_orco_verts_transform(ob->data, &co4, 1, 1);
977                                 }
978                         }
979                         else {
980                                 v1 = &mesh->mvert[mf->v1];
981                                 v2 = &mesh->mvert[mf->v2];
982                                 v3 = &mesh->mvert[mf->v3];
983                                 copy_v3_v3(co1, v1->co);
984                                 copy_v3_v3(co2, v2->co);
985                                 copy_v3_v3(co3, v3->co);
986                                 if (mf->v4) {
987                                         v4 = &mesh->mvert[mf->v4];
988                                         copy_v3_v3(co4, v4->co);
989                                 }
990                         }
991
992                         cur = mf->v4 ? area_quad_v3(co1, co2, co3, co4) : area_tri_v3(co1, co2, co3);
993
994                         if (cur > maxweight)
995                                 maxweight = cur;
996
997                         element_weight[i] = cur;
998                         totarea += cur;
999                 }
1000
1001                 for (i=0; i<totelem; i++)
1002                         element_weight[i] /= totarea;
1003
1004                 maxweight /= totarea;
1005         }
1006         else {
1007                 float min=1.0f/(float)(MIN2(totelem,totpart));
1008                 for (i=0; i<totelem; i++)
1009                         element_weight[i]=min;
1010                 maxweight=min;
1011         }
1012
1013         /* Calculate weights from vgroup */
1014         vweight = psys_cache_vgroup(mesh,psys,PSYS_VG_DENSITY);
1015
1016         if (vweight) {
1017                 if (from==PART_FROM_VERT) {
1018                         for (i=0;i<totelem; i++)
1019                                 element_weight[i]*=vweight[i];
1020                 }
1021                 else { /* PART_FROM_FACE / PART_FROM_VOLUME */
1022                         for (i=0;i<totelem; i++) {
1023                                 MFace *mf = &mesh->mface[i];
1024                                 tweight = vweight[mf->v1] + vweight[mf->v2] + vweight[mf->v3];
1025
1026                                 if (mf->v4) {
1027                                         tweight += vweight[mf->v4];
1028                                         tweight /= 4.0f;
1029                                 }
1030                                 else {
1031                                         tweight /= 3.0f;
1032                                 }
1033
1034                                 element_weight[i]*=tweight;
1035                         }
1036                 }
1037                 MEM_freeN(vweight);
1038         }
1039
1040         /* Calculate total weight of all elements */
1041         int totmapped = 0;
1042         totweight = 0.0f;
1043         for (i = 0; i < totelem; i++) {
1044                 if (element_weight[i] > 0.0f) {
1045                         totmapped++;
1046                         totweight += element_weight[i];
1047                 }
1048         }
1049
1050         if (totmapped == 0) {
1051                 /* We are not allowed to distribute particles anywhere... */
1052                 return 0;
1053         }
1054
1055         inv_totweight = 1.0f / totweight;
1056
1057         /* Calculate cumulative weights.
1058          * We remove all null-weighted elements from element_sum, and create a new mapping
1059          * 'activ'_elem_index -> orig_elem_index.
1060          * This simplifies greatly the filtering of zero-weighted items - and can be much more efficient
1061          * especially in random case (reducing a lot the size of binary-searched array)...
1062          */
1063         float *element_sum = MEM_mallocN(sizeof(*element_sum) * totmapped, __func__);
1064         int *element_map = MEM_mallocN(sizeof(*element_map) * totmapped, __func__);
1065         int i_mapped = 0;
1066
1067         for (i = 0; i < totelem && element_weight[i] == 0.0f; i++);
1068         element_sum[i_mapped] = element_weight[i] * inv_totweight;
1069         element_map[i_mapped] = i;
1070         i_mapped++;
1071         for (i++; i < totelem; i++) {
1072                 if (element_weight[i] > 0.0f) {
1073                         element_sum[i_mapped] = element_sum[i_mapped - 1] + element_weight[i] * inv_totweight;
1074                         /* Skip elements which weight is so small that it does not affect the sum. */
1075                         if (element_sum[i_mapped] > element_sum[i_mapped - 1]) {
1076                                 element_map[i_mapped] = i;
1077                                 i_mapped++;
1078                         }
1079                 }
1080         }
1081         totmapped = i_mapped;
1082
1083         /* Finally assign elements to particles */
1084         if (part->flag & PART_TRAND) {
1085                 for (p = 0; p < totpart; p++) {
1086                         /* In theory element_sum[totmapped - 1] should be 1.0,
1087                          * but due to float errors this is not necessarily always true, so scale pos accordingly. */
1088                         const float pos = BLI_rng_get_float(rng) * element_sum[totmapped - 1];
1089                         const int eidx = distribute_binary_search(element_sum, totmapped, pos);
1090                         particle_element[p] = element_map[eidx];
1091                         BLI_assert(pos <= element_sum[eidx]);
1092                         BLI_assert(eidx ? (pos > element_sum[eidx - 1]) : (pos >= 0.0f));
1093                         jitter_offset[particle_element[p]] = pos;
1094                 }
1095         }
1096         else {
1097                 double step, pos;
1098
1099                 step = (totpart < 2) ? 0.5 : 1.0 / (double)totpart;
1100                 /* This is to address tricky issues with vertex-emitting when user tries (and expects) exact 1-1 vert/part
1101                  * distribution (see T47983 and its two example files). It allows us to consider pos as
1102                  * 'midpoint between v and v+1' (or 'p and p+1', depending whether we have more vertices than particles or not),
1103                  * and avoid stumbling over float impression in element_sum.
1104                  * Note: moved face and volume distribution to this as well (instead of starting at zero),
1105                  * for the same reasons, see T52682. */
1106                 pos = (totpart < totmapped) ? 0.5 / (double)totmapped : step * 0.5;  /* We choose the smaller step. */
1107
1108                 for (i = 0, p = 0; p < totpart; p++, pos += step) {
1109                         for ( ; (i < totmapped - 1) && (pos > (double)element_sum[i]); i++);
1110
1111                         particle_element[p] = element_map[i];
1112
1113                         jitter_offset[particle_element[p]] = pos;
1114                 }
1115         }
1116
1117         MEM_freeN(element_sum);
1118         MEM_freeN(element_map);
1119
1120         /* For hair, sort by origindex (allows optimization's in rendering), */
1121         /* however with virtual parents the children need to be in random order. */
1122         if (part->type == PART_HAIR && !(part->childtype==PART_CHILD_FACES && part->parents != 0.0f)) {
1123                 int *orig_index = NULL;
1124
1125                 if (from == PART_FROM_VERT) {
1126                         if (mesh->totvert)
1127                                 orig_index = CustomData_get_layer(&mesh->vdata, CD_ORIGINDEX);
1128                 }
1129                 else {
1130                         if (mesh->totface)
1131                                 orig_index = CustomData_get_layer(&mesh->fdata, CD_ORIGINDEX);
1132                 }
1133
1134                 if (orig_index) {
1135                         BLI_qsort_r(particle_element, totpart, sizeof(int), distribute_compare_orig_index, orig_index);
1136                 }
1137         }
1138
1139         /* Create jittering if needed */
1140         if (distr==PART_DISTR_JIT && ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
1141                 jitlevel= part->userjit;
1142
1143                 if (jitlevel == 0) {
1144                         jitlevel= totpart/totelem;
1145                         if (part->flag & PART_EDISTR) jitlevel*= 2;     /* looks better in general, not very scientific */
1146                         if (jitlevel<3) jitlevel= 3;
1147                 }
1148
1149                 jit= MEM_callocN((2+ jitlevel*2)*sizeof(float), "jit");
1150
1151                 /* for small amounts of particles we use regular jitter since it looks
1152                  * a bit better, for larger amounts we switch to hammersley sequence
1153                  * because it is much faster */
1154                 if (jitlevel < 25)
1155                         init_mv_jit(jit, jitlevel, psys->seed, part->jitfac);
1156                 else
1157                         hammersley_create(jit, jitlevel+1, psys->seed, part->jitfac);
1158                 BLI_array_randomize(jit, 2*sizeof(float), jitlevel, psys->seed); /* for custom jit or even distribution */
1159         }
1160
1161         /* Setup things for threaded distribution */
1162         ctx->tree= tree;
1163         ctx->seams= seams;
1164         ctx->totseam= totseam;
1165         ctx->sim.psys= psys;
1166         ctx->index= particle_element;
1167         ctx->jit= jit;
1168         ctx->jitlevel= jitlevel;
1169         ctx->jitoff= jitter_offset;
1170         ctx->weight= element_weight;
1171         ctx->maxweight= maxweight;
1172         ctx->cfrom= cfrom;
1173         ctx->distr= distr;
1174         ctx->mesh= mesh;
1175         ctx->tpars= tpars;
1176
1177         if (children) {
1178                 alloc_child_particles(psys, totpart);
1179         }
1180
1181         BLI_rng_free(rng);
1182
1183         return 1;
1184 }
1185
1186 static void psys_task_init_distribute(ParticleTask *task, ParticleSimulationData *sim)
1187 {
1188         /* init random number generator */
1189         int seed = 31415926 + sim->psys->seed;
1190
1191         task->rng = BLI_rng_new(seed);
1192 }
1193
1194 static void distribute_particles_on_dm(ParticleSimulationData *sim, int from)
1195 {
1196         TaskScheduler *task_scheduler;
1197         TaskPool *task_pool;
1198         ParticleThreadContext ctx;
1199         ParticleTask *tasks;
1200         Mesh *final_mesh = sim->psmd->mesh_final;
1201         int i, totpart, numtasks;
1202
1203         /* create a task pool for distribution tasks */
1204         if (!psys_thread_context_init_distribute(&ctx, sim, from))
1205                 return;
1206
1207         task_scheduler = BLI_task_scheduler_get();
1208         task_pool = BLI_task_pool_create(task_scheduler, &ctx);
1209
1210         totpart = (from == PART_FROM_CHILD ? sim->psys->totchild : sim->psys->totpart);
1211         psys_tasks_create(&ctx, 0, totpart, &tasks, &numtasks);
1212         for (i = 0; i < numtasks; ++i) {
1213                 ParticleTask *task = &tasks[i];
1214
1215                 psys_task_init_distribute(task, sim);
1216                 if (from == PART_FROM_CHILD)
1217                         BLI_task_pool_push(task_pool, exec_distribute_child, task, false, TASK_PRIORITY_LOW);
1218                 else
1219                         BLI_task_pool_push(task_pool, exec_distribute_parent, task, false, TASK_PRIORITY_LOW);
1220         }
1221         BLI_task_pool_work_and_wait(task_pool);
1222
1223         BLI_task_pool_free(task_pool);
1224
1225         psys_calc_dmcache(sim->ob, final_mesh, sim->psmd->mesh_original, sim->psys);
1226
1227         if (ctx.mesh != final_mesh)
1228                 BKE_id_free(NULL, ctx.mesh);
1229
1230         psys_tasks_free(tasks, numtasks);
1231
1232         psys_thread_context_free(&ctx);
1233 }
1234
1235 /* ready for future use, to emit particles without geometry */
1236 static void distribute_particles_on_shape(ParticleSimulationData *sim, int UNUSED(from))
1237 {
1238         distribute_invalid(sim, 0);
1239
1240         fprintf(stderr,"Shape emission not yet possible!\n");
1241 }
1242
1243 void distribute_particles(ParticleSimulationData *sim, int from)
1244 {
1245         PARTICLE_PSMD;
1246         int distr_error=0;
1247
1248         if (psmd) {
1249                 if (psmd->mesh_final)
1250                         distribute_particles_on_dm(sim, from);
1251                 else
1252                         distr_error=1;
1253         }
1254         else
1255                 distribute_particles_on_shape(sim, from);
1256
1257         if (distr_error) {
1258                 distribute_invalid(sim, from);
1259
1260                 fprintf(stderr,"Particle distribution error!\n");
1261         }
1262 }