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