1 //
2 // Copyright 2013 Pixar
3 //
4 // Licensed under the Apache License, Version 2.0 (the "Apache License")
5 // with the following modification; you may not use this file except in
6 // compliance with the Apache License and the following modification to it:
7 // Section 6. Trademarks. is deleted and replaced with:
8 //
9 // 6. Trademarks. This License does not grant permission to use the trade
10 // names, trademarks, service marks, or product names of the Licensor
11 // and its affiliates, except as required to comply with Section 4(c) of
12 // the License and to reproduce the content of the NOTICE file.
13 //
14 // You may obtain a copy of the Apache License at
15 //
16 // http://www.apache.org/licenses/LICENSE-2.0
17 //
18 // Unless required by applicable law or agreed to in writing, software
19 // distributed under the Apache License with the above modification is
20 // distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
21 // KIND, either express or implied. See the Apache License for the specific
22 // language governing permissions and limitations under the Apache License.
23 //
24
25 #include "particles.h"
26
27 #include <opensubdiv/far/ptexIndices.h>
28 #include <opensubdiv/far/patchMap.h>
29 #include <opensubdiv/sdc/types.h>
30
31 #include <cassert>
32 #include <cmath>
33
34 using namespace OpenSubdiv;
35
36 void
37 UpdateParticle(float speed,
38 STParticles::Position *p,
39 float *dp,
40 Osd::PatchCoord *patchCoord,
41 int regFaceSize,
42 std::vector<STParticles::FaceInfo> const &adjacency,
43 Far::PatchMap const *patchMap);
44
45 #ifdef OPENSUBDIV_HAS_TBB
46 #include <tbb/parallel_for.h>
47 class TbbUpdateKernel {
48 public:
TbbUpdateKernel(float speed,STParticles::Position * positions,float * velocities,Osd::PatchCoord * patchCoords,int regFaceSize,std::vector<STParticles::FaceInfo> const & adjacency,Far::PatchMap const * patchMap)49 TbbUpdateKernel(float speed,
50 STParticles::Position *positions,
51 float *velocities,
52 Osd::PatchCoord *patchCoords,
53 int regFaceSize,
54 std::vector<STParticles::FaceInfo> const &adjacency,
55 Far::PatchMap const *patchMap) :
56 _speed(speed), _positions(positions), _velocities(velocities),
57 _patchCoords(patchCoords),
58 _regFaceSize(regFaceSize), _adjacency(adjacency), _patchMap(patchMap) {
59 }
60
operator ()(tbb::blocked_range<int> const & r) const61 void operator () (tbb::blocked_range<int> const &r) const {
62 for (int i = r.begin(); i < r.end(); ++i) {
63 STParticles::Position * p = _positions + i;
64 float *dp = _velocities + i*2;
65 Osd::PatchCoord *patchCoord = &_patchCoords[i];
66
67 UpdateParticle(_speed, p, dp, patchCoord, _regFaceSize, _adjacency, _patchMap);
68 }
69 }
70 private:
71 float _speed;
72 STParticles::Position *_positions;
73 float *_velocities;
74 Osd::PatchCoord *_patchCoords;
75 int _regFaceSize;
76 std::vector<STParticles::FaceInfo> const &_adjacency;
77 Far::PatchMap const *_patchMap;
78 };
79 #endif
80
STParticles(Refiner const & refiner,PatchTable const * patchTable,int nParticles,bool centered)81 STParticles::STParticles(Refiner const & refiner,
82 PatchTable const *patchTable,
83 int nParticles, bool centered)
84 : _speed(1.0f)
85 , _regFaceSize(
86 Sdc::SchemeTypeTraits::GetRegularFaceSize(refiner.GetSchemeType())) {
87
88 Far::PtexIndices ptexIndices(refiner);
89
90 // Create a far patch map
91 _patchMap = new Far::PatchMap(*patchTable);
92
93 int nPtexFaces = ptexIndices.GetNumFaces();
94
95 srand(static_cast<int>(2147483647));
96
97 { // initialize positions
98
99 _positions.resize(nParticles);
100 Position * pos = &_positions[0];
101
102 for (int i = 0; i < nParticles; ++i) {
103 pos->ptexIndex = std::min(
104 (int)(((float)rand()/(float)RAND_MAX) * nPtexFaces), nPtexFaces-1);
105 if (_regFaceSize==3) {
106 pos->s = centered ? 1.0f/3.0f : (float)rand()/(float)RAND_MAX;
107 pos->t = centered ? 1.0f/3.0f : (float)rand()/(float)RAND_MAX;
108 // Keep locations within the triangular parametric domain
109 if ((pos->s+pos->t) >= 1.0f) {
110 pos->s = 1.0f - pos->s;
111 pos->t = 1.0f - pos->t;
112 }
113 } else {
114 pos->s = centered ? 0.5f : (float)rand()/(float)RAND_MAX;
115 pos->t = centered ? 0.5f : (float)rand()/(float)RAND_MAX;
116 }
117 ++pos;
118 }
119 }
120
121 { // initialize velocities
122 _velocities.resize(nParticles * 2);
123
124 for (int i = 0; i < nParticles; ++i) {
125 // initialize normalized random directions
126 float s = 2.0f*(float)rand()/(float)RAND_MAX - 1.0f,
127 t = 2.0f*(float)rand()/(float)RAND_MAX - 1.0f,
128 l = sqrtf(s*s+t*t);
129
130 _velocities[2*i ] = s / l;
131 _velocities[2*i+1] = t / l;
132 }
133 }
134
135 if (_regFaceSize == 4) { // initialize topology adjacency
136 _adjacency.resize(nPtexFaces);
137
138 Far::TopologyLevel const & refBaseLevel = refiner.GetLevel(0);
139
140 int nfaces = refBaseLevel.GetNumFaces(),
141 adjfaces[4],
142 adjedges[4];
143
144 for (int face=0, ptexface=0; face<nfaces; ++face) {
145
146 Far::ConstIndexArray fverts = refBaseLevel.GetFaceVertices(face);
147
148 if (fverts.size()==_regFaceSize) {
149 ptexIndices.GetAdjacency(refiner, face, 0, adjfaces, adjedges);
150 _adjacency[ptexface] = FaceInfo(adjfaces, adjedges, false);
151 ++ptexface;
152 } else {
153 for (int vert=0; vert<fverts.size(); ++vert) {
154 ptexIndices.GetAdjacency(refiner, face, vert, adjfaces, adjedges);
155 _adjacency[ptexface+vert] =
156 FaceInfo(adjfaces, adjedges, true);
157 }
158 ptexface+=fverts.size();
159 }
160 }
161 }
162 //std::cout << *this;
163 }
164
~STParticles()165 STParticles::~STParticles() {
166 delete _patchMap;
167 }
168
169 inline void
FlipS(STParticles::Position * p,float * dp)170 FlipS(STParticles::Position * p, float * dp) {
171 p->s = 1.0f-p->s;
172 dp[0] = - dp[0];
173 }
174
175 inline void
FlipT(STParticles::Position * p,float * dp)176 FlipT(STParticles::Position * p, float * dp) {
177 p->t = 1.0f-p->t;
178 dp[1] = -dp[1];
179 }
180
181 inline void
SwapST(STParticles::Position * p,float * dp)182 SwapST(STParticles::Position * p, float * dp) {
183 std::swap(p->s, p->t);
184 std::swap(dp[0], dp[1]);
185 }
186
187 inline void
RotateQuad(int rot,STParticles::Position * p,float * dp)188 RotateQuad(int rot, STParticles::Position * p, float * dp) {
189 switch (rot & 3) {
190 default: return;
191 case 1: FlipS(p, dp); SwapST(p, dp); break;
192 case 2: FlipS(p, dp); FlipT(p, dp); break;
193 case 3: FlipT(p, dp); SwapST(p, dp); break;
194 }
195 assert((p->s>=0.0f) && (p->s<=1.0f) && (p->t>=0.0f) && (p->t<=1.0f));
196 }
197
198 inline void
TrimQuad(STParticles::Position * p)199 TrimQuad(STParticles::Position * p) {
200 if (p->s <0.0f) p->s = 1.0f + p->s;
201 if (p->s>=1.0f) p->s = p->s - 1.0f;
202 if (p->t <0.0f) p->t = 1.0f + p->t;
203 if (p->t>=1.0f) p->t = p->t - 1.0f;
204 assert((p->s>=0.0f) && (p->s<=1.0f) && (p->t>=0.0f) && (p->t<=1.0f));
205 }
206
207 inline void
ClampQuad(STParticles::Position * p)208 ClampQuad(STParticles::Position * p) {
209 if (p->s<0.0f) {
210 p->s=0.0f;
211 } else if (p->s>1.0f) {
212 p->s=1.0f;
213 }
214 if (p->t<0.0f) {
215 p->t=0.0f;
216 } else if (p->t>1.0f) {
217 p->t=1.0f;
218 }
219 }
220
221 inline void
BounceQuad(int edge,STParticles::Position * p,float * dp)222 BounceQuad(int edge, STParticles::Position * p, float * dp) {
223 switch (edge) {
224 case 0: assert(p->t<=0.0f); p->t = -p->t; dp[1] = -dp[1]; break;
225 case 1: assert(p->s>=1.0f); p->s = 2.0f - p->s; dp[0] = -dp[0]; break;
226 case 2: assert(p->t>=1.0f); p->t = 2.0f - p->t; dp[1] = -dp[1]; break;
227 case 3: assert(p->s<=0.0f); p->s = -p->s; dp[0] = -dp[0]; break;
228 }
229
230 // because 'diagonal' cases aren't handled, stick particles to edges when
231 // if they cross 2 boundaries
232 ClampQuad(p);
233 assert((p->s>=0.0f) && (p->s<=1.0f) && (p->t>=0.0f) && (p->t<=1.0f));
234 }
235
236 void
WarpQuad(std::vector<STParticles::FaceInfo> const & adjacency,int edge,STParticles::Position * p,float * dp)237 WarpQuad(std::vector<STParticles::FaceInfo> const &adjacency,
238 int edge, STParticles::Position * p, float * dp) {
239 assert(p->ptexIndex<(int)adjacency.size() && (edge>=0 && edge<4));
240
241 STParticles::FaceInfo const & f = adjacency[p->ptexIndex];
242
243 int afid = f.adjface(edge),
244 aeid = f.adjedge(edge);
245
246 if (afid==-1) {
247 // boundary detected: bounce the particle
248 BounceQuad(edge, p, dp);
249 } else {
250 STParticles::FaceInfo const & af = adjacency[afid];
251 int rot = edge - aeid + 2;
252
253 bool fIsSubface = f.isSubface(),
254 afIsSubface = af.isSubface();
255
256 if (fIsSubface != afIsSubface) {
257 // XXXX manuelk domain should be split properly
258 BounceQuad(edge, p, dp);
259 } else {
260 TrimQuad(p);
261 RotateQuad(rot, p, dp);
262 p->ptexIndex = afid; // move particle to adjacent face
263 }
264 }
265 assert((p->s>=0.0f) && (p->s<=1.0f) && (p->t>=0.0f) && (p->t<=1.0f));
266 }
267
268 void
ConstrainQuad(STParticles::Position * p,float * dp,std::vector<STParticles::FaceInfo> const & adjacency)269 ConstrainQuad(STParticles::Position *p,
270 float *dp,
271 std::vector<STParticles::FaceInfo> const &adjacency) {
272
273 // make sure particles can't skip more than 1 face boundary at a time
274 assert((p->s>-2.0f) && (p->s<2.0f) && (p->t>-2.0f) && (p->t<2.0f));
275
276 // check if the particle is jumping a boundary
277 // note: a particle can jump 2 edges at a time (a "diagonal" jump)
278 // this is not treated here.
279 int edge = -1;
280 if (p->s >= 1.0f) edge = 1;
281 if (p->s <= 0.0f) edge = 3;
282 if (p->t >= 1.0f) edge = 2;
283 if (p->t <= 0.0f) edge = 0;
284
285 if (edge>=0) {
286 // warp the particle to the other side of the boundary
287 WarpQuad(adjacency, edge, p, dp);
288 }
289 assert((p->s>=0.0f) && (p->s<=1.0f) && (p->t>=0.0f) && (p->t<=1.0f));
290 }
291
292 inline void
ClampTri(STParticles::Position * p)293 ClampTri(STParticles::Position * p) {
294 if (p->s<0.0f) {
295 p->s=0.0f;
296 } else if (p->s>1.0f) {
297 p->s=1.0f;
298 }
299 if (p->t<0.0f) {
300 p->t=0.0f;
301 } else if (p->t>1.0f) {
302 p->t=1.0f;
303 }
304 if ((p->s+p->t)>=1.0f) {
305 p->s = 1.0f-p->t;
306 p->t = 1.0f-p->s;
307 }
308 }
309
310 inline void
BounceTri(int edge,STParticles::Position * p,float * dp)311 BounceTri(int edge, STParticles::Position * p, float * dp) {
312 switch (edge) {
313 case 0:
314 assert(p->t<=0.0f);
315 p->t = -p->t; dp[1] = -dp[1];
316 break;
317 case 1:
318 assert((p->s+p->t)>=1.0f);
319 p->s = 1.0f-p->s; dp[0] = -dp[0];
320 p->t = 1.0f-p->t; dp[1] = -dp[1];
321 break;
322 case 2:
323 assert(p->s<=0.0f);
324 p->s = -p->s; dp[0] = -dp[0];
325 break;
326 }
327
328 // because 'diagonal' cases aren't handled, stick particles to edges when
329 // if they cross 2 boundaries
330 ClampTri(p);
331 assert((p->s>=0.0f) && (p->s<=1.0f) && (p->t>=0.0f) && (p->t<=1.0f) &&
332 ((p->s+p->t)<=1.0f));
333 }
334
335 void
WarpTri(std::vector<STParticles::FaceInfo> const &,int edge,STParticles::Position * p,float * dp)336 WarpTri(std::vector<STParticles::FaceInfo> const &,
337 int edge, STParticles::Position * p, float * dp) {
338
339 // For now, particles on triangle meshes just bounce.
340 BounceTri(edge, p, dp);
341 assert((p->s>=0.0f) && (p->s<=1.0f) && (p->t>=0.0f) && (p->t<=1.0f) &&
342 ((p->s+p->t)<=1.0f));
343 }
344
345 void
ConstrainTri(STParticles::Position * p,float * dp,std::vector<STParticles::FaceInfo> const & adjacency)346 ConstrainTri(STParticles::Position *p,
347 float *dp,
348 std::vector<STParticles::FaceInfo> const &adjacency) {
349
350 // make sure particles can't skip more than 1 face boundary at a time
351 assert((p->s>-2.0f) && (p->s<2.0f) && (p->t>-2.0f) && (p->t<2.0f) &&
352 ((p->s+p->t)>-2.0f) && ((p->s+p->t)<2.0f));
353
354 // check if the particle is jumping a boundary
355 // note: a particle can jump 2 edges at a time (a "diagonal" jump)
356 // this is not treated here.
357 int edge = -1;
358 if (p->t <= 0.0f) edge = 0;
359 if (p->s <= 0.0f) edge = 2;
360 if ((p->s+p->t) >= 1.0f) edge = 1;
361
362 if (edge>=0) {
363 // warp the particle to the other side of the boundary
364 WarpTri(adjacency, edge, p, dp);
365 }
366
367 assert((p->s>-2.0f) && (p->s<2.0f) && (p->t>-2.0f) && (p->t<2.0f) &&
368 ((p->s+p->t)>-2.0f) && ((p->s+p->t)<2.0f));
369 }
370
371 void
UpdateParticle(float speed,STParticles::Position * p,float * dp,Osd::PatchCoord * patchCoord,int regFaceSize,std::vector<STParticles::FaceInfo> const & adjacency,Far::PatchMap const * patchMap)372 UpdateParticle(float speed,
373 STParticles::Position *p,
374 float *dp,
375 Osd::PatchCoord *patchCoord,
376 int regFaceSize,
377 std::vector<STParticles::FaceInfo> const &adjacency,
378 Far::PatchMap const *patchMap) {
379
380 // apply velocity
381 p->s += dp[0] * speed;
382 p->t += dp[1] * speed;
383
384 if (regFaceSize == 3) {
385 ConstrainTri(p, dp, adjacency);
386 } else {
387 ConstrainQuad(p, dp, adjacency);
388 }
389
390 // resolve particle positions into patch handles
391 Far::PatchTable::PatchHandle const *handle =
392 patchMap->FindPatch(p->ptexIndex, p->s, p->t);
393 if (handle) {
394 *patchCoord = Osd::PatchCoord(*handle, p->s, p->t);
395 }
396 }
397
398 void
Update(float deltaTime)399 STParticles::Update(float deltaTime) {
400
401 if (deltaTime == 0) return;
402 float speed = GetSpeed() * std::max(0.001f, std::min(deltaTime, 0.5f));
403
404 _patchCoords.resize(GetNumParticles());
405
406 #ifdef OPENSUBDIV_HAS_TBB
407 TbbUpdateKernel kernel(speed, &_positions[0], &_velocities[0],
408 &_patchCoords[0],
409 _regFaceSize, _adjacency, _patchMap);
410 tbb::blocked_range<int> range(0, GetNumParticles(), 256);
411 tbb::parallel_for(range, kernel);
412 #else
413 for (int i=0; i<GetNumParticles(); ++i) {
414 Position * p = &_positions[i];
415 float * dp = &_velocities[i*2];
416 Osd::PatchCoord *patchCoord = &_patchCoords[i];
417
418 UpdateParticle(speed, p, dp, patchCoord, _regFaceSize, _adjacency, _patchMap);
419 }
420 #endif
421 }
422
423 // Dump adjacency info
operator <<(std::ostream & os,STParticles::FaceInfo const & f)424 std::ostream & operator << (std::ostream & os,
425 STParticles::FaceInfo const & f) {
426
427 os << " adjface: " << f.adjfaces[0] << ' '
428 << f.adjfaces[1] << ' '
429 << f.adjfaces[2] << ' '
430 << f.adjfaces[3]
431 << " adjedge: " << f.adjedge(0) << ' '
432 << f.adjedge(1) << ' '
433 << f.adjedge(2) << ' '
434 << f.adjedge(3)
435 << " flags:";
436
437 if (f.flags == 0) {
438 os << " (none)";
439 } else {
440 if (f.isSubface()) {
441 std::cout << " subface";
442 }
443 }
444 os << std::endl;
445
446 return os;
447 }
448
operator <<(std::ostream & os,STParticles const & particles)449 std::ostream & operator << (std::ostream & os,
450 STParticles const & particles) {
451
452 for (int i=0; i<(int)particles._adjacency.size(); ++i) {
453 os << particles._adjacency[i];
454 }
455
456 return os;
457 }
458
459
460