Blender  V3.3
implicit_blender.c
Go to the documentation of this file.
1 /* SPDX-License-Identifier: GPL-2.0-or-later
2  * Copyright Blender Foundation. All rights reserved. */
3 
8 #include "implicit.h"
9 
10 #ifdef IMPLICIT_SOLVER_BLENDER
11 
12 # include "MEM_guardedalloc.h"
13 
14 # include "DNA_meshdata_types.h"
15 # include "DNA_object_force_types.h"
16 # include "DNA_object_types.h"
17 # include "DNA_scene_types.h"
18 # include "DNA_texture_types.h"
19 
20 # include "BLI_math.h"
21 # include "BLI_utildefines.h"
22 
23 # include "BKE_cloth.h"
24 # include "BKE_collision.h"
25 # include "BKE_effect.h"
26 
27 # include "SIM_mass_spring.h"
28 
29 # ifdef __GNUC__
30 # pragma GCC diagnostic ignored "-Wtype-limits"
31 # endif
32 
33 # ifdef _OPENMP
34 # define CLOTH_OPENMP_LIMIT 512
35 # endif
36 
37 //#define DEBUG_TIME
38 
39 # ifdef DEBUG_TIME
40 # include "PIL_time.h"
41 # endif
42 
43 static float I[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
44 static float ZERO[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
45 
46 # if 0
47 # define C99
48 # ifdef C99
49 # defineDO_INLINE inline
50 # else
51 # defineDO_INLINE static
52 # endif
53 # endif /* if 0 */
54 
55 struct Cloth;
56 
58 /* fast vector / matrix library, enhancements are welcome :) -dg */
60 
61 /* DEFINITIONS */
62 typedef float lfVector[3];
63 typedef struct fmatrix3x3 {
64  float m[3][3]; /* 3x3 matrix */
65  unsigned int c, r; /* column and row number */
66  // int pinned; /* is this vertex allowed to move? */
67  float n1, n2, n3; /* three normal vectors for collision constrains */
68  unsigned int vcount; /* vertex count */
69  unsigned int scount; /* spring count */
71 
73 /* float[3] vector */
75 /* simple vector code */
76 /* STATUS: verified */
77 DO_INLINE void mul_fvector_S(float to[3], const float from[3], float scalar)
78 {
79  to[0] = from[0] * scalar;
80  to[1] = from[1] * scalar;
81  to[2] = from[2] * scalar;
82 }
83 /* simple v^T * v product ("outer product") */
84 /* STATUS: HAS TO BE verified (*should* work) */
85 DO_INLINE void mul_fvectorT_fvector(float to[3][3], const float vectorA[3], const float vectorB[3])
86 {
87  mul_fvector_S(to[0], vectorB, vectorA[0]);
88  mul_fvector_S(to[1], vectorB, vectorA[1]);
89  mul_fvector_S(to[2], vectorB, vectorA[2]);
90 }
91 /* simple v^T * v product with scalar ("outer product") */
92 /* STATUS: HAS TO BE verified (*should* work) */
93 DO_INLINE void mul_fvectorT_fvectorS(float to[3][3], float vectorA[3], float vectorB[3], float aS)
94 {
95  mul_fvectorT_fvector(to, vectorA, vectorB);
96 
97  mul_fvector_S(to[0], to[0], aS);
98  mul_fvector_S(to[1], to[1], aS);
99  mul_fvector_S(to[2], to[2], aS);
100 }
101 
102 # if 0
103 /* printf vector[3] on console: for debug output */
104 static void print_fvector(float m3[3])
105 {
106  printf("%f\n%f\n%f\n\n", m3[0], m3[1], m3[2]);
107 }
108 
110 /* long float vector float (*)[3] */
112 /* print long vector on console: for debug output */
113 DO_INLINE void print_lfvector(float (*fLongVector)[3], unsigned int verts)
114 {
115  unsigned int i = 0;
116  for (i = 0; i < verts; i++) {
117  print_fvector(fLongVector[i]);
118  }
119 }
120 # endif
121 
122 /* create long vector */
124 {
125  /* TODO: check if memory allocation was successful */
126  return (lfVector *)MEM_callocN(verts * sizeof(lfVector), "cloth_implicit_alloc_vector");
127  // return (lfVector *)cloth_aligned_malloc(&MEMORY_BASE, verts * sizeof(lfVector));
128 }
129 /* delete long vector */
130 DO_INLINE void del_lfvector(float (*fLongVector)[3])
131 {
132  if (fLongVector != NULL) {
133  MEM_freeN(fLongVector);
134  // cloth_aligned_free(&MEMORY_BASE, fLongVector);
135  }
136 }
137 /* copy long vector */
138 DO_INLINE void cp_lfvector(float (*to)[3], float (*from)[3], unsigned int verts)
139 {
140  memcpy(to, from, verts * sizeof(lfVector));
141 }
142 /* init long vector with float[3] */
143 DO_INLINE void init_lfvector(float (*fLongVector)[3], const float vector[3], unsigned int verts)
144 {
145  unsigned int i = 0;
146  for (i = 0; i < verts; i++) {
147  copy_v3_v3(fLongVector[i], vector);
148  }
149 }
150 /* zero long vector with float[3] */
151 DO_INLINE void zero_lfvector(float (*to)[3], unsigned int verts)
152 {
153  memset(to, 0.0f, verts * sizeof(lfVector));
154 }
155 /* Multiply long vector with scalar. */
156 DO_INLINE void mul_lfvectorS(float (*to)[3],
157  float (*fLongVector)[3],
158  float scalar,
159  unsigned int verts)
160 {
161  unsigned int i = 0;
162 
163  for (i = 0; i < verts; i++) {
164  mul_fvector_S(to[i], fLongVector[i], scalar);
165  }
166 }
167 /* Multiply long vector with scalar.
168  * `A -= B * float` */
169 DO_INLINE void submul_lfvectorS(float (*to)[3],
170  float (*fLongVector)[3],
171  float scalar,
172  unsigned int verts)
173 {
174  unsigned int i = 0;
175  for (i = 0; i < verts; i++) {
176  VECSUBMUL(to[i], fLongVector[i], scalar);
177  }
178 }
179 /* dot product for big vector */
180 DO_INLINE float dot_lfvector(float (*fLongVectorA)[3],
181  float (*fLongVectorB)[3],
182  unsigned int verts)
183 {
184  long i = 0;
185  float temp = 0.0;
186  /* XXX brecht, disabled this for now (first schedule line was already disabled),
187  * due to non-commutative nature of floating point ops this makes the sim give
188  * different results each time you run it!
189  * schedule(guided, 2) */
190  //#pragma omp parallel for reduction(+: temp) if (verts > CLOTH_OPENMP_LIMIT)
191  for (i = 0; i < (long)verts; i++) {
192  temp += dot_v3v3(fLongVectorA[i], fLongVectorB[i]);
193  }
194  return temp;
195 }
196 /* `A = B + C` -> for big vector. */
197 DO_INLINE void add_lfvector_lfvector(float (*to)[3],
198  float (*fLongVectorA)[3],
199  float (*fLongVectorB)[3],
200  unsigned int verts)
201 {
202  unsigned int i = 0;
203 
204  for (i = 0; i < verts; i++) {
205  add_v3_v3v3(to[i], fLongVectorA[i], fLongVectorB[i]);
206  }
207 }
208 /* `A = B + C * float` -> for big vector. */
209 DO_INLINE void add_lfvector_lfvectorS(float (*to)[3],
210  float (*fLongVectorA)[3],
211  float (*fLongVectorB)[3],
212  float bS,
213  unsigned int verts)
214 {
215  unsigned int i = 0;
216 
217  for (i = 0; i < verts; i++) {
218  VECADDS(to[i], fLongVectorA[i], fLongVectorB[i], bS);
219  }
220 }
221 /* `A = B * float + C * float` -> for big vector */
222 DO_INLINE void add_lfvectorS_lfvectorS(float (*to)[3],
223  float (*fLongVectorA)[3],
224  float aS,
225  float (*fLongVectorB)[3],
226  float bS,
227  unsigned int verts)
228 {
229  unsigned int i = 0;
230 
231  for (i = 0; i < verts; i++) {
232  VECADDSS(to[i], fLongVectorA[i], aS, fLongVectorB[i], bS);
233  }
234 }
235 /* `A = B - C * float` -> for big vector. */
236 DO_INLINE void sub_lfvector_lfvectorS(float (*to)[3],
237  float (*fLongVectorA)[3],
238  float (*fLongVectorB)[3],
239  float bS,
240  unsigned int verts)
241 {
242  unsigned int i = 0;
243  for (i = 0; i < verts; i++) {
244  VECSUBS(to[i], fLongVectorA[i], fLongVectorB[i], bS);
245  }
246 }
247 /* `A = B - C` -> for big vector. */
248 DO_INLINE void sub_lfvector_lfvector(float (*to)[3],
249  float (*fLongVectorA)[3],
250  float (*fLongVectorB)[3],
251  unsigned int verts)
252 {
253  unsigned int i = 0;
254 
255  for (i = 0; i < verts; i++) {
256  sub_v3_v3v3(to[i], fLongVectorA[i], fLongVectorB[i]);
257  }
258 }
260 // 3x3 matrix
262 # if 0
263 /* printf 3x3 matrix on console: for debug output */
264 static void print_fmatrix(float m3[3][3])
265 {
266  printf("%f\t%f\t%f\n", m3[0][0], m3[0][1], m3[0][2]);
267  printf("%f\t%f\t%f\n", m3[1][0], m3[1][1], m3[1][2]);
268  printf("%f\t%f\t%f\n\n", m3[2][0], m3[2][1], m3[2][2]);
269 }
270 
271 static void print_sparse_matrix(fmatrix3x3 *m)
272 {
273  if (m) {
274  unsigned int i;
275  for (i = 0; i < m[0].vcount + m[0].scount; i++) {
276  printf("%d:\n", i);
277  print_fmatrix(m[i].m);
278  }
279  }
280 }
281 # endif
282 
283 # if 0
284 static void print_lvector(lfVector *v, int numverts)
285 {
286  int i;
287  for (i = 0; i < numverts; i++) {
288  if (i > 0) {
289  printf("\n");
290  }
291 
292  printf("%f,\n", v[i][0]);
293  printf("%f,\n", v[i][1]);
294  printf("%f,\n", v[i][2]);
295  }
296 }
297 # endif
298 
299 # if 0
300 static void print_bfmatrix(fmatrix3x3 *m)
301 {
302  int tot = m[0].vcount + m[0].scount;
303  int size = m[0].vcount * 3;
304  float *t = MEM_callocN(sizeof(float) * size * size, "bfmatrix");
305  int q, i, j;
306 
307  for (q = 0; q < tot; q++) {
308  int k = 3 * m[q].r;
309  int l = 3 * m[q].c;
310 
311  for (j = 0; j < 3; j++) {
312  for (i = 0; i < 3; i++) {
313 # if 0
314  if (t[k + i + (l + j) * size] != 0.0f) {
315  printf("warning: overwriting value at %d, %d\n", m[q].r, m[q].c);
316  }
317 # endif
318  if (k == l) {
319  t[k + i + (k + j) * size] += m[q].m[i][j];
320  }
321  else {
322  t[k + i + (l + j) * size] += m[q].m[i][j];
323  t[l + j + (k + i) * size] += m[q].m[j][i];
324  }
325  }
326  }
327  }
328 
329  for (j = 0; j < size; j++) {
330  if (j > 0 && j % 3 == 0) {
331  printf("\n");
332  }
333 
334  for (i = 0; i < size; i++) {
335  if (i > 0 && i % 3 == 0) {
336  printf(" ");
337  }
338 
340  }
341  printf("\n");
342  }
343 
344  MEM_freeN(t);
345 }
346 # endif
347 
348 /* copy 3x3 matrix */
349 DO_INLINE void cp_fmatrix(float to[3][3], const float from[3][3])
350 {
351  // memcpy(to, from, sizeof(float[3][3]));
352  copy_v3_v3(to[0], from[0]);
353  copy_v3_v3(to[1], from[1]);
354  copy_v3_v3(to[2], from[2]);
355 }
356 
357 /* copy 3x3 matrix */
358 DO_INLINE void initdiag_fmatrixS(float to[3][3], float aS)
359 {
360  cp_fmatrix(to, ZERO);
361 
362  to[0][0] = aS;
363  to[1][1] = aS;
364  to[2][2] = aS;
365 }
366 
367 # if 0
368 /* calculate determinant of 3x3 matrix */
369 DO_INLINE float det_fmatrix(float m[3][3])
370 {
371  return m[0][0] * m[1][1] * m[2][2] + m[1][0] * m[2][1] * m[0][2] + m[0][1] * m[1][2] * m[2][0] -
372  m[0][0] * m[1][2] * m[2][1] - m[0][1] * m[1][0] * m[2][2] - m[2][0] * m[1][1] * m[0][2];
373 }
374 
375 DO_INLINE void inverse_fmatrix(float to[3][3], float from[3][3])
376 {
377  unsigned int i, j;
378  float d;
379 
380  if ((d = det_fmatrix(from)) == 0) {
381  printf("can't build inverse");
382  exit(0);
383  }
384  for (i = 0; i < 3; i++) {
385  for (j = 0; j < 3; j++) {
386  int i1 = (i + 1) % 3;
387  int i2 = (i + 2) % 3;
388  int j1 = (j + 1) % 3;
389  int j2 = (j + 2) % 3;
391  to[j][i] = (from[i1][j1] * from[i2][j2] - from[i1][j2] * from[i2][j1]) / d;
402  }
403  }
404 }
405 # endif
406 
407 /* 3x3 matrix multiplied by a scalar */
408 /* STATUS: verified */
409 DO_INLINE void mul_fmatrix_S(float matrix[3][3], float scalar)
410 {
411  mul_fvector_S(matrix[0], matrix[0], scalar);
412  mul_fvector_S(matrix[1], matrix[1], scalar);
413  mul_fvector_S(matrix[2], matrix[2], scalar);
414 }
415 
416 /* a vector multiplied by a 3x3 matrix */
417 /* STATUS: verified */
418 DO_INLINE void mul_fvector_fmatrix(float *to, const float *from, const float matrix[3][3])
419 {
420  to[0] = matrix[0][0] * from[0] + matrix[1][0] * from[1] + matrix[2][0] * from[2];
421  to[1] = matrix[0][1] * from[0] + matrix[1][1] * from[1] + matrix[2][1] * from[2];
422  to[2] = matrix[0][2] * from[0] + matrix[1][2] * from[1] + matrix[2][2] * from[2];
423 }
424 
425 /* 3x3 matrix multiplied by a vector */
426 /* STATUS: verified */
427 DO_INLINE void mul_fmatrix_fvector(float *to, const float matrix[3][3], const float from[3])
428 {
429  to[0] = dot_v3v3(matrix[0], from);
430  to[1] = dot_v3v3(matrix[1], from);
431  to[2] = dot_v3v3(matrix[2], from);
432 }
433 /* 3x3 matrix addition with 3x3 matrix */
434 DO_INLINE void add_fmatrix_fmatrix(float to[3][3],
435  const float matrixA[3][3],
436  const float matrixB[3][3])
437 {
438  add_v3_v3v3(to[0], matrixA[0], matrixB[0]);
439  add_v3_v3v3(to[1], matrixA[1], matrixB[1]);
440  add_v3_v3v3(to[2], matrixA[2], matrixB[2]);
441 }
442 /* `A -= B*x + (C * y)` (3x3 matrix sub-addition with 3x3 matrix). */
444  float to[3][3], const float matrixA[3][3], float aS, const float matrixB[3][3], float bS)
445 {
446  VECSUBADDSS(to[0], matrixA[0], aS, matrixB[0], bS);
447  VECSUBADDSS(to[1], matrixA[1], aS, matrixB[1], bS);
448  VECSUBADDSS(to[2], matrixA[2], aS, matrixB[2], bS);
449 }
450 /* `A = B - C` (3x3 matrix subtraction with 3x3 matrix). */
451 DO_INLINE void sub_fmatrix_fmatrix(float to[3][3],
452  const float matrixA[3][3],
453  const float matrixB[3][3])
454 {
455  sub_v3_v3v3(to[0], matrixA[0], matrixB[0]);
456  sub_v3_v3v3(to[1], matrixA[1], matrixB[1]);
457  sub_v3_v3v3(to[2], matrixA[2], matrixB[2]);
458 }
460 /* special functions */
462 /* 3x3 matrix multiplied+added by a vector */
463 /* STATUS: verified */
464 DO_INLINE void muladd_fmatrix_fvector(float to[3], const float matrix[3][3], const float from[3])
465 {
466  to[0] += dot_v3v3(matrix[0], from);
467  to[1] += dot_v3v3(matrix[1], from);
468  to[2] += dot_v3v3(matrix[2], from);
469 }
470 
471 DO_INLINE void muladd_fmatrixT_fvector(float to[3], const float matrix[3][3], const float from[3])
472 {
473  to[0] += matrix[0][0] * from[0] + matrix[1][0] * from[1] + matrix[2][0] * from[2];
474  to[1] += matrix[0][1] * from[0] + matrix[1][1] * from[1] + matrix[2][1] * from[2];
475  to[2] += matrix[0][2] * from[0] + matrix[1][2] * from[1] + matrix[2][2] * from[2];
476 }
477 
478 BLI_INLINE void outerproduct(float r[3][3], const float a[3], const float b[3])
479 {
480  mul_v3_v3fl(r[0], a, b[0]);
481  mul_v3_v3fl(r[1], a, b[1]);
482  mul_v3_v3fl(r[2], a, b[2]);
483 }
484 
485 BLI_INLINE void cross_m3_v3m3(float r[3][3], const float v[3], const float m[3][3])
486 {
487  cross_v3_v3v3(r[0], v, m[0]);
488  cross_v3_v3v3(r[1], v, m[1]);
489  cross_v3_v3v3(r[2], v, m[2]);
490 }
491 
492 BLI_INLINE void cross_v3_identity(float r[3][3], const float v[3])
493 {
494  r[0][0] = 0.0f;
495  r[1][0] = v[2];
496  r[2][0] = -v[1];
497  r[0][1] = -v[2];
498  r[1][1] = 0.0f;
499  r[2][1] = v[0];
500  r[0][2] = v[1];
501  r[1][2] = -v[0];
502  r[2][2] = 0.0f;
503 }
504 
505 BLI_INLINE void madd_m3_m3fl(float r[3][3], const float m[3][3], float f)
506 {
507  r[0][0] += m[0][0] * f;
508  r[0][1] += m[0][1] * f;
509  r[0][2] += m[0][2] * f;
510  r[1][0] += m[1][0] * f;
511  r[1][1] += m[1][1] * f;
512  r[1][2] += m[1][2] * f;
513  r[2][0] += m[2][0] * f;
514  r[2][1] += m[2][1] * f;
515  r[2][2] += m[2][2] * f;
516 }
517 
519 
521 /* SPARSE SYMMETRIC big matrix with 3x3 matrix entries */
523 /* printf a big matrix on console: for debug output */
524 # if 0
525 static void print_bfmatrix(fmatrix3x3 *m3)
526 {
527  unsigned int i = 0;
528 
529  for (i = 0; i < m3[0].vcount + m3[0].scount; i++) {
530  print_fmatrix(m3[i].m);
531  }
532 }
533 # endif
534 
535 BLI_INLINE void init_fmatrix(fmatrix3x3 *matrix, int r, int c)
536 {
537  matrix->r = r;
538  matrix->c = c;
539 }
540 
541 /* create big matrix */
542 DO_INLINE fmatrix3x3 *create_bfmatrix(unsigned int verts, unsigned int springs)
543 {
544  /* TODO: check if memory allocation was successful */
545  fmatrix3x3 *temp = (fmatrix3x3 *)MEM_callocN(sizeof(fmatrix3x3) * (verts + springs),
546  "cloth_implicit_alloc_matrix");
547  int i;
548 
549  temp[0].vcount = verts;
550  temp[0].scount = springs;
551 
552  /* vertex part of the matrix is diagonal blocks */
553  for (i = 0; i < verts; i++) {
554  init_fmatrix(temp + i, i, i);
555  }
556 
557  return temp;
558 }
559 /* delete big matrix */
561 {
562  if (matrix != NULL) {
563  MEM_freeN(matrix);
564  }
565 }
566 
567 /* copy big matrix */
569 {
570  /* TODO: bounds checking. */
571  memcpy(to, from, sizeof(fmatrix3x3) * (from[0].vcount + from[0].scount));
572 }
573 
574 /* init big matrix */
575 /* slow in parallel */
576 DO_INLINE void init_bfmatrix(fmatrix3x3 *matrix, float m3[3][3])
577 {
578  unsigned int i;
579 
580  for (i = 0; i < matrix[0].vcount + matrix[0].scount; i++) {
581  cp_fmatrix(matrix[i].m, m3);
582  }
583 }
584 
585 /* init the diagonal of big matrix */
586 /* slow in parallel */
587 DO_INLINE void initdiag_bfmatrix(fmatrix3x3 *matrix, float m3[3][3])
588 {
589  unsigned int i, j;
590  float tmatrix[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
591 
592  for (i = 0; i < matrix[0].vcount; i++) {
593  cp_fmatrix(matrix[i].m, m3);
594  }
595  for (j = matrix[0].vcount; j < matrix[0].vcount + matrix[0].scount; j++) {
596  cp_fmatrix(matrix[j].m, tmatrix);
597  }
598 }
599 
600 /* SPARSE SYMMETRIC multiply big matrix with long vector. */
601 /* STATUS: verified */
602 DO_INLINE void mul_bfmatrix_lfvector(float (*to)[3], fmatrix3x3 *from, lfVector *fLongVector)
603 {
604  unsigned int vcount = from[0].vcount;
605  lfVector *temp = create_lfvector(vcount);
606 
607  zero_lfvector(to, vcount);
608 
609 # pragma omp parallel sections if (vcount > CLOTH_OPENMP_LIMIT)
610  {
611 # pragma omp section
612  {
613  for (unsigned int i = from[0].vcount; i < from[0].vcount + from[0].scount; i++) {
614  /* This is the lower triangle of the sparse matrix,
615  * therefore multiplication occurs with transposed submatrices. */
616  muladd_fmatrixT_fvector(to[from[i].c], from[i].m, fLongVector[from[i].r]);
617  }
618  }
619 # pragma omp section
620  {
621  for (unsigned int i = 0; i < from[0].vcount + from[0].scount; i++) {
622  muladd_fmatrix_fvector(temp[from[i].r], from[i].m, fLongVector[from[i].c]);
623  }
624  }
625  }
626  add_lfvector_lfvector(to, to, temp, from[0].vcount);
627 
628  del_lfvector(temp);
629 }
630 
631 /* SPARSE SYMMETRIC sub big matrix with big matrix. */
632 /* A -= B * float + C * float --> for big matrix */
633 /* VERIFIED */
635  fmatrix3x3 *to, fmatrix3x3 *from, float aS, fmatrix3x3 *matrix, float bS)
636 {
637  unsigned int i = 0;
638 
639  /* process diagonal elements */
640  for (i = 0; i < matrix[0].vcount + matrix[0].scount; i++) {
641  subadd_fmatrixS_fmatrixS(to[i].m, from[i].m, aS, matrix[i].m, bS);
642  }
643 }
644 
646 /* simulator start */
648 
649 typedef struct Implicit_Data {
650  /* inputs */
651  fmatrix3x3 *bigI; /* identity (constant) */
652  fmatrix3x3 *tfm; /* local coordinate transform */
653  fmatrix3x3 *M; /* masses */
654  lfVector *F; /* forces */
655  fmatrix3x3 *dFdV, *dFdX; /* force jacobians */
656  int num_blocks; /* number of off-diagonal blocks (springs) */
657 
658  /* motion state data */
659  lfVector *X, *Xnew; /* positions */
660  lfVector *V, *Vnew; /* velocities */
661 
662  /* internal solver data */
663  lfVector *B; /* B for A*dV = B */
664  fmatrix3x3 *A; /* A for A*dV = B */
665 
666  lfVector *dV; /* velocity change (solution of A*dV = B) */
667  lfVector *z; /* target velocity in constrained directions */
668  fmatrix3x3 *S; /* filtering matrix for constraints */
669  fmatrix3x3 *P, *Pinv; /* pre-conditioning matrix */
671 
672 Implicit_Data *SIM_mass_spring_solver_create(int numverts, int numsprings)
673 {
674  Implicit_Data *id = (Implicit_Data *)MEM_callocN(sizeof(Implicit_Data), "implicit vecmat");
675 
676  /* process diagonal elements */
677  id->tfm = create_bfmatrix(numverts, 0);
678  id->A = create_bfmatrix(numverts, numsprings);
679  id->dFdV = create_bfmatrix(numverts, numsprings);
680  id->dFdX = create_bfmatrix(numverts, numsprings);
681  id->S = create_bfmatrix(numverts, 0);
682  id->Pinv = create_bfmatrix(numverts, numsprings);
683  id->P = create_bfmatrix(numverts, numsprings);
684  id->bigI = create_bfmatrix(numverts, numsprings); /* TODO: 0 springs. */
685  id->M = create_bfmatrix(numverts, numsprings);
686  id->X = create_lfvector(numverts);
687  id->Xnew = create_lfvector(numverts);
688  id->V = create_lfvector(numverts);
689  id->Vnew = create_lfvector(numverts);
690  id->F = create_lfvector(numverts);
691  id->B = create_lfvector(numverts);
692  id->dV = create_lfvector(numverts);
693  id->z = create_lfvector(numverts);
694 
695  initdiag_bfmatrix(id->bigI, I);
696 
697  return id;
698 }
699 
701 {
702  del_bfmatrix(id->tfm);
703  del_bfmatrix(id->A);
704  del_bfmatrix(id->dFdV);
705  del_bfmatrix(id->dFdX);
706  del_bfmatrix(id->S);
707  del_bfmatrix(id->P);
708  del_bfmatrix(id->Pinv);
709  del_bfmatrix(id->bigI);
710  del_bfmatrix(id->M);
711 
712  del_lfvector(id->X);
713  del_lfvector(id->Xnew);
714  del_lfvector(id->V);
715  del_lfvector(id->Vnew);
716  del_lfvector(id->F);
717  del_lfvector(id->B);
718  del_lfvector(id->dV);
719  del_lfvector(id->z);
720 
721  MEM_freeN(id);
722 }
723 
724 /* ==== Transformation from/to root reference frames ==== */
725 
726 BLI_INLINE void world_to_root_v3(Implicit_Data *data, int index, float r[3], const float v[3])
727 {
728  copy_v3_v3(r, v);
729  mul_transposed_m3_v3(data->tfm[index].m, r);
730 }
731 
732 BLI_INLINE void root_to_world_v3(Implicit_Data *data, int index, float r[3], const float v[3])
733 {
734  mul_v3_m3v3(r, data->tfm[index].m, v);
735 }
736 
738  int index,
739  float r[3][3],
740  const float m[3][3])
741 {
742  float trot[3][3];
743  copy_m3_m3(trot, data->tfm[index].m);
744  transpose_m3(trot);
745  mul_m3_m3m3(r, trot, m);
746 }
747 
749  int index,
750  float r[3][3],
751  const float m[3][3])
752 {
753  mul_m3_m3m3(r, data->tfm[index].m, m);
754 }
755 
756 /* ================================ */
757 
759 {
760  unsigned int i = 0;
761 
762  for (i = 0; i < S[0].vcount; i++) {
763  mul_m3_v3(S[i].m, V[S[i].r]);
764  }
765 }
766 
767 /* this version of the CG algorithm does not work very well with partial constraints
768  * (where S has non-zero elements). */
769 # if 0
770 static int cg_filtered(lfVector *ldV, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fmatrix3x3 *S)
771 {
772  /* Solves for unknown X in equation AX=B */
773  unsigned int conjgrad_loopcount = 0, conjgrad_looplimit = 100;
774  float conjgrad_epsilon = 0.0001f /* , conjgrad_lasterror=0 */ /* UNUSED */;
775  lfVector *q, *d, *tmp, *r;
776  float s, starget, a, s_prev;
777  unsigned int numverts = lA[0].vcount;
778  q = create_lfvector(numverts);
779  d = create_lfvector(numverts);
780  tmp = create_lfvector(numverts);
781  r = create_lfvector(numverts);
782 
783  // zero_lfvector(ldV, CLOTHPARTICLES);
784  filter(ldV, S);
785 
786  add_lfvector_lfvector(ldV, ldV, z, numverts);
787 
788  // r = B - Mul(tmp, A, X); /* just use B if X known to be zero. */
789  cp_lfvector(r, lB, numverts);
790  mul_bfmatrix_lfvector(tmp, lA, ldV);
791  sub_lfvector_lfvector(r, r, tmp, numverts);
792 
793  filter(r, S);
794 
795  cp_lfvector(d, r, numverts);
796 
797  s = dot_lfvector(r, r, numverts);
798  starget = s * sqrtf(conjgrad_epsilon);
799 
800  while (s > starget && conjgrad_loopcount < conjgrad_looplimit) {
801  // Mul(q, A, d); /* q = A*d; */
802  mul_bfmatrix_lfvector(q, lA, d);
803 
804  filter(q, S);
805 
806  a = s / dot_lfvector(d, q, numverts);
807 
808  /* X = X + d*a; */
809  add_lfvector_lfvectorS(ldV, ldV, d, a, numverts);
810 
811  /* r = r - q*a; */
812  sub_lfvector_lfvectorS(r, r, q, a, numverts);
813 
814  s_prev = s;
815  s = dot_lfvector(r, r, numverts);
816 
817  // d = r+d*(s/s_prev);
818  add_lfvector_lfvectorS(d, r, d, (s / s_prev), numverts);
819 
820  filter(d, S);
821 
822  conjgrad_loopcount++;
823  }
824  /* conjgrad_lasterror = s; */ /* UNUSED */
825 
826  del_lfvector(q);
827  del_lfvector(d);
828  del_lfvector(tmp);
829  del_lfvector(r);
830  // printf("W/O conjgrad_loopcount: %d\n", conjgrad_loopcount);
831 
832  return conjgrad_loopcount <
833  conjgrad_looplimit; /* true means we reached desired accuracy in given time - ie stable */
834 }
835 # endif
836 
837 static int cg_filtered(lfVector *ldV,
838  fmatrix3x3 *lA,
839  lfVector *lB,
840  lfVector *z,
841  fmatrix3x3 *S,
843 {
844  /* Solves for unknown X in equation AX=B */
845  unsigned int conjgrad_loopcount = 0, conjgrad_looplimit = 100;
846  float conjgrad_epsilon = 0.01f;
847 
848  unsigned int numverts = lA[0].vcount;
849  lfVector *fB = create_lfvector(numverts);
850  lfVector *AdV = create_lfvector(numverts);
851  lfVector *r = create_lfvector(numverts);
852  lfVector *c = create_lfvector(numverts);
853  lfVector *q = create_lfvector(numverts);
854  lfVector *s = create_lfvector(numverts);
855  float bnorm2, delta_new, delta_old, delta_target, alpha;
856 
857  cp_lfvector(ldV, z, numverts);
858 
859  /* d0 = filter(B)^T * P * filter(B) */
860  cp_lfvector(fB, lB, numverts);
861  filter(fB, S);
862  bnorm2 = dot_lfvector(fB, fB, numverts);
863  delta_target = conjgrad_epsilon * conjgrad_epsilon * bnorm2;
864 
865  /* r = filter(B - A * dV) */
866  mul_bfmatrix_lfvector(AdV, lA, ldV);
867  sub_lfvector_lfvector(r, lB, AdV, numverts);
868  filter(r, S);
869 
870  /* c = filter(P^-1 * r) */
871  cp_lfvector(c, r, numverts);
872  filter(c, S);
873 
874  /* delta = r^T * c */
875  delta_new = dot_lfvector(r, c, numverts);
876 
877 # ifdef IMPLICIT_PRINT_SOLVER_INPUT_OUTPUT
878  printf("==== A ====\n");
879  print_bfmatrix(lA);
880  printf("==== z ====\n");
881  print_lvector(z, numverts);
882  printf("==== B ====\n");
883  print_lvector(lB, numverts);
884  printf("==== S ====\n");
885  print_bfmatrix(S);
886 # endif
887 
888  while (delta_new > delta_target && conjgrad_loopcount < conjgrad_looplimit) {
889  mul_bfmatrix_lfvector(q, lA, c);
890  filter(q, S);
891 
892  alpha = delta_new / dot_lfvector(c, q, numverts);
893 
894  add_lfvector_lfvectorS(ldV, ldV, c, alpha, numverts);
895 
896  add_lfvector_lfvectorS(r, r, q, -alpha, numverts);
897 
898  /* s = P^-1 * r */
899  cp_lfvector(s, r, numverts);
900  delta_old = delta_new;
901  delta_new = dot_lfvector(r, s, numverts);
902 
903  add_lfvector_lfvectorS(c, s, c, delta_new / delta_old, numverts);
904  filter(c, S);
905 
906  conjgrad_loopcount++;
907  }
908 
909 # ifdef IMPLICIT_PRINT_SOLVER_INPUT_OUTPUT
910  printf("==== dV ====\n");
911  print_lvector(ldV, numverts);
912  printf("========\n");
913 # endif
914 
915  del_lfvector(fB);
916  del_lfvector(AdV);
917  del_lfvector(r);
918  del_lfvector(c);
919  del_lfvector(q);
920  del_lfvector(s);
921  // printf("W/O conjgrad_loopcount: %d\n", conjgrad_loopcount);
922 
923  result->status = conjgrad_loopcount < conjgrad_looplimit ? SIM_SOLVER_SUCCESS :
925  result->iterations = conjgrad_loopcount;
926  result->error = bnorm2 > 0.0f ? sqrtf(delta_new / bnorm2) : 0.0f;
927 
928  return conjgrad_loopcount <
929  conjgrad_looplimit; /* true means we reached desired accuracy in given time - ie stable */
930 }
931 
932 # if 0
933 /* block diagonalizer */
934 DO_INLINE void BuildPPinv(fmatrix3x3 *lA, fmatrix3x3 *P, fmatrix3x3 *Pinv)
935 {
936  unsigned int i = 0;
937 
938  /* Take only the diagonal blocks of A */
939  // #pragma omp parallel for private(i) if (lA[0].vcount > CLOTH_OPENMP_LIMIT)
940  for (i = 0; i < lA[0].vcount; i++) {
941  /* block diagonalizer */
942  cp_fmatrix(P[i].m, lA[i].m);
943  inverse_fmatrix(Pinv[i].m, P[i].m);
944  }
945 }
946 
947 # if 0
948 /* version 1.3 */
949 static int cg_filtered_pre(lfVector *dv,
950  fmatrix3x3 *lA,
951  lfVector *lB,
952  lfVector *z,
953  fmatrix3x3 *S,
954  fmatrix3x3 *P,
955  fmatrix3x3 *Pinv)
956 {
957  unsigned int numverts = lA[0].vcount, iterations = 0, conjgrad_looplimit = 100;
958  float delta0 = 0, deltaNew = 0, deltaOld = 0, alpha = 0;
959  float conjgrad_epsilon = 0.0001; /* 0.2 is dt for steps=5 */
960  lfVector *r = create_lfvector(numverts);
961  lfVector *p = create_lfvector(numverts);
962  lfVector *s = create_lfvector(numverts);
963  lfVector *h = create_lfvector(numverts);
964 
965  BuildPPinv(lA, P, Pinv);
966 
967  filter(dv, S);
968  add_lfvector_lfvector(dv, dv, z, numverts);
969 
970  mul_bfmatrix_lfvector(r, lA, dv);
971  sub_lfvector_lfvector(r, lB, r, numverts);
972  filter(r, S);
973 
974  mul_prevfmatrix_lfvector(p, Pinv, r);
975  filter(p, S);
976 
977  deltaNew = dot_lfvector(r, p, numverts);
978 
979  delta0 = deltaNew * sqrt(conjgrad_epsilon);
980 
981 # ifdef DEBUG_TIME
982  double start = PIL_check_seconds_timer();
983 # endif
984 
985  while ((deltaNew > delta0) && (iterations < conjgrad_looplimit)) {
986  iterations++;
987 
988  mul_bfmatrix_lfvector(s, lA, p);
989  filter(s, S);
990 
991  alpha = deltaNew / dot_lfvector(p, s, numverts);
992 
993  add_lfvector_lfvectorS(dv, dv, p, alpha, numverts);
994 
995  add_lfvector_lfvectorS(r, r, s, -alpha, numverts);
996 
997  mul_prevfmatrix_lfvector(h, Pinv, r);
998  filter(h, S);
999 
1000  deltaOld = deltaNew;
1001 
1002  deltaNew = dot_lfvector(r, h, numverts);
1003 
1004  add_lfvector_lfvectorS(p, h, p, deltaNew / deltaOld, numverts);
1005 
1006  filter(p, S);
1007  }
1008 
1009 # ifdef DEBUG_TIME
1010  double end = PIL_check_seconds_timer();
1011  printf("cg_filtered_pre time: %f\n", (float)(end - start));
1012 # endif
1013 
1014  del_lfvector(h);
1015  del_lfvector(s);
1016  del_lfvector(p);
1017  del_lfvector(r);
1018 
1019  printf("iterations: %d\n", iterations);
1020 
1021  return iterations < conjgrad_looplimit;
1022 }
1023 # endif
1024 
1025 /* version 1.4 */
1026 static int cg_filtered_pre(lfVector *dv,
1027  fmatrix3x3 *lA,
1028  lfVector *lB,
1029  lfVector *z,
1030  fmatrix3x3 *S,
1031  fmatrix3x3 *P,
1032  fmatrix3x3 *Pinv,
1033  fmatrix3x3 *bigI)
1034 {
1035  unsigned int numverts = lA[0].vcount, iterations = 0, conjgrad_looplimit = 100;
1036  float delta0 = 0, deltaNew = 0, deltaOld = 0, alpha = 0, tol = 0;
1037  lfVector *r = create_lfvector(numverts);
1038  lfVector *p = create_lfvector(numverts);
1039  lfVector *s = create_lfvector(numverts);
1040  lfVector *h = create_lfvector(numverts);
1041  lfVector *bhat = create_lfvector(numverts);
1042  lfVector *btemp = create_lfvector(numverts);
1043 
1044  BuildPPinv(lA, P, Pinv);
1045 
1046  initdiag_bfmatrix(bigI, I);
1047  sub_bfmatrix_Smatrix(bigI, bigI, S);
1048 
1049  /* x = Sx_0+(I-S)z */
1050  filter(dv, S);
1051  add_lfvector_lfvector(dv, dv, z, numverts);
1052 
1053  /* b_hat = S(b-A(I-S)z) */
1054  mul_bfmatrix_lfvector(r, lA, z);
1055  mul_bfmatrix_lfvector(bhat, bigI, r);
1056  sub_lfvector_lfvector(bhat, lB, bhat, numverts);
1057 
1058  /* r = S(b-Ax) */
1059  mul_bfmatrix_lfvector(r, lA, dv);
1060  sub_lfvector_lfvector(r, lB, r, numverts);
1061  filter(r, S);
1062 
1063  /* p = SP^-1r */
1064  mul_prevfmatrix_lfvector(p, Pinv, r);
1065  filter(p, S);
1066 
1067  /* delta0 = bhat^TP^-1bhat */
1068  mul_prevfmatrix_lfvector(btemp, Pinv, bhat);
1069  delta0 = dot_lfvector(bhat, btemp, numverts);
1070 
1071  /* deltaNew = r^TP */
1072  deltaNew = dot_lfvector(r, p, numverts);
1073 
1074 # if 0
1075  filter(dv, S);
1076  add_lfvector_lfvector(dv, dv, z, numverts);
1077 
1078  mul_bfmatrix_lfvector(r, lA, dv);
1079  sub_lfvector_lfvector(r, lB, r, numverts);
1080  filter(r, S);
1081 
1082  mul_prevfmatrix_lfvector(p, Pinv, r);
1083  filter(p, S);
1084 
1085  deltaNew = dot_lfvector(r, p, numverts);
1086 
1087  delta0 = deltaNew * sqrt(conjgrad_epsilon);
1088 # endif
1089 
1090 # ifdef DEBUG_TIME
1091  double start = PIL_check_seconds_timer();
1092 # endif
1093 
1094  tol = (0.01 * 0.2);
1095 
1096  while ((deltaNew > delta0 * tol * tol) && (iterations < conjgrad_looplimit)) {
1097  iterations++;
1098 
1099  mul_bfmatrix_lfvector(s, lA, p);
1100  filter(s, S);
1101 
1102  alpha = deltaNew / dot_lfvector(p, s, numverts);
1103 
1104  add_lfvector_lfvectorS(dv, dv, p, alpha, numverts);
1105 
1106  add_lfvector_lfvectorS(r, r, s, -alpha, numverts);
1107 
1108  mul_prevfmatrix_lfvector(h, Pinv, r);
1109  filter(h, S);
1110 
1111  deltaOld = deltaNew;
1112 
1113  deltaNew = dot_lfvector(r, h, numverts);
1114 
1115  add_lfvector_lfvectorS(p, h, p, deltaNew / deltaOld, numverts);
1116 
1117  filter(p, S);
1118  }
1119 
1120 # ifdef DEBUG_TIME
1121  double end = PIL_check_seconds_timer();
1122  printf("cg_filtered_pre time: %f\n", (float)(end - start));
1123 # endif
1124 
1125  del_lfvector(btemp);
1126  del_lfvector(bhat);
1127  del_lfvector(h);
1128  del_lfvector(s);
1129  del_lfvector(p);
1130  del_lfvector(r);
1131 
1132  // printf("iterations: %d\n", iterations);
1133 
1134  return iterations < conjgrad_looplimit;
1135 }
1136 # endif
1137 
1139 {
1140  unsigned int numverts = data->dFdV[0].vcount;
1141 
1142  lfVector *dFdXmV = create_lfvector(numverts);
1143  zero_lfvector(data->dV, numverts);
1144 
1145  cp_bfmatrix(data->A, data->M);
1146 
1147  subadd_bfmatrixS_bfmatrixS(data->A, data->dFdV, dt, data->dFdX, (dt * dt));
1148 
1149  mul_bfmatrix_lfvector(dFdXmV, data->dFdX, data->V);
1150 
1151  add_lfvectorS_lfvectorS(data->B, data->F, dt, dFdXmV, (dt * dt), numverts);
1152 
1153 # ifdef DEBUG_TIME
1154  double start = PIL_check_seconds_timer();
1155 # endif
1156 
1157  /* Conjugate gradient algorithm to solve Ax=b. */
1158  cg_filtered(data->dV, data->A, data->B, data->z, data->S, result);
1159 
1160  // cg_filtered_pre(id->dV, id->A, id->B, id->z, id->S, id->P, id->Pinv, id->bigI);
1161 
1162 # ifdef DEBUG_TIME
1163  double end = PIL_check_seconds_timer();
1164  printf("cg_filtered calc time: %f\n", (float)(end - start));
1165 # endif
1166 
1167  /* advance velocities */
1168  add_lfvector_lfvector(data->Vnew, data->V, data->dV, numverts);
1169 
1170  del_lfvector(dFdXmV);
1171 
1172  return result->status == SIM_SOLVER_SUCCESS;
1173 }
1174 
1176 {
1177  int numverts = data->M[0].vcount;
1178 
1179  /* advance positions */
1180  add_lfvector_lfvectorS(data->Xnew, data->X, data->Vnew, dt, numverts);
1181 
1182  return true;
1183 }
1184 
1186 {
1187  int numverts = data->M[0].vcount;
1188  cp_lfvector(data->X, data->Xnew, numverts);
1189  cp_lfvector(data->V, data->Vnew, numverts);
1190 }
1191 
1193 {
1194  unit_m3(data->M[index].m);
1195  mul_m3_fl(data->M[index].m, mass);
1196 }
1197 
1198 void SIM_mass_spring_set_rest_transform(Implicit_Data *data, int index, float tfm[3][3])
1199 {
1200 # ifdef CLOTH_ROOT_FRAME
1201  copy_m3_m3(data->tfm[index].m, tfm);
1202 # else
1203  unit_m3(data->tfm[index].m);
1204  (void)tfm;
1205 # endif
1206 }
1207 
1209  int index,
1210  const float x[3],
1211  const float v[3])
1212 {
1213  world_to_root_v3(data, index, data->X[index], x);
1214  world_to_root_v3(data, index, data->V[index], v);
1215 }
1216 
1217 void SIM_mass_spring_set_position(Implicit_Data *data, int index, const float x[3])
1218 {
1219  world_to_root_v3(data, index, data->X[index], x);
1220 }
1221 
1222 void SIM_mass_spring_set_velocity(Implicit_Data *data, int index, const float v[3])
1223 {
1224  world_to_root_v3(data, index, data->V[index], v);
1225 }
1226 
1228  int index,
1229  float x[3],
1230  float v[3])
1231 {
1232  if (x) {
1233  root_to_world_v3(data, index, x, data->X[index]);
1234  }
1235  if (v) {
1236  root_to_world_v3(data, index, v, data->V[index]);
1237  }
1238 }
1239 
1240 void SIM_mass_spring_get_position(struct Implicit_Data *data, int index, float x[3])
1241 {
1242  root_to_world_v3(data, index, x, data->X[index]);
1243 }
1244 
1245 void SIM_mass_spring_get_velocity(struct Implicit_Data *data, int index, float v[3])
1246 {
1247  root_to_world_v3(data, index, v, data->V[index]);
1248 }
1249 
1250 void SIM_mass_spring_get_new_position(struct Implicit_Data *data, int index, float x[3])
1251 {
1252  root_to_world_v3(data, index, x, data->Xnew[index]);
1253 }
1254 
1255 void SIM_mass_spring_set_new_position(struct Implicit_Data *data, int index, const float x[3])
1256 {
1257  world_to_root_v3(data, index, data->Xnew[index], x);
1258 }
1259 
1260 void SIM_mass_spring_get_new_velocity(struct Implicit_Data *data, int index, float v[3])
1261 {
1262  root_to_world_v3(data, index, v, data->Vnew[index]);
1263 }
1264 
1265 void SIM_mass_spring_set_new_velocity(struct Implicit_Data *data, int index, const float v[3])
1266 {
1267  world_to_root_v3(data, index, data->Vnew[index], v);
1268 }
1269 
1270 /* -------------------------------- */
1271 
1273 {
1274  int s = data->M[0].vcount + data->num_blocks; /* index from array start */
1275  BLI_assert(s < data->M[0].vcount + data->M[0].scount);
1276  ++data->num_blocks;
1277 
1278  /* tfm and S don't have spring entries (diagonal blocks only) */
1279  init_fmatrix(data->bigI + s, v1, v2);
1280  init_fmatrix(data->M + s, v1, v2);
1281  init_fmatrix(data->dFdX + s, v1, v2);
1282  init_fmatrix(data->dFdV + s, v1, v2);
1283  init_fmatrix(data->A + s, v1, v2);
1284  init_fmatrix(data->P + s, v1, v2);
1285  init_fmatrix(data->Pinv + s, v1, v2);
1286 
1287  return s;
1288 }
1289 
1291 {
1292  int i, numverts = data->S[0].vcount;
1293  for (i = 0; i < numverts; i++) {
1294  unit_m3(data->S[i].m);
1295  zero_v3(data->z[i]);
1296  }
1297 }
1298 
1299 void SIM_mass_spring_add_constraint_ndof0(Implicit_Data *data, int index, const float dV[3])
1300 {
1301  zero_m3(data->S[index].m);
1302 
1303  world_to_root_v3(data, index, data->z[index], dV);
1304 }
1305 
1307  Implicit_Data *data, int index, const float c1[3], const float c2[3], const float dV[3])
1308 {
1309  float m[3][3], p[3], q[3], u[3], cmat[3][3];
1310 
1311  world_to_root_v3(data, index, p, c1);
1312  mul_fvectorT_fvector(cmat, p, p);
1313  sub_m3_m3m3(m, I, cmat);
1314 
1315  world_to_root_v3(data, index, q, c2);
1316  mul_fvectorT_fvector(cmat, q, q);
1317  sub_m3_m3m3(m, m, cmat);
1318 
1319  /* XXX not sure but multiplication should work here */
1320  copy_m3_m3(data->S[index].m, m);
1321  // mul_m3_m3m3(data->S[index].m, data->S[index].m, m);
1322 
1323  world_to_root_v3(data, index, u, dV);
1324  add_v3_v3(data->z[index], u);
1325 }
1326 
1328  int index,
1329  const float c1[3],
1330  const float dV[3])
1331 {
1332  float m[3][3], p[3], u[3], cmat[3][3];
1333 
1334  world_to_root_v3(data, index, p, c1);
1335  mul_fvectorT_fvector(cmat, p, p);
1336  sub_m3_m3m3(m, I, cmat);
1337 
1338  copy_m3_m3(data->S[index].m, m);
1339  // mul_m3_m3m3(data->S[index].m, data->S[index].m, m);
1340 
1341  world_to_root_v3(data, index, u, dV);
1342  add_v3_v3(data->z[index], u);
1343 }
1344 
1346 {
1347  int numverts = data->M[0].vcount;
1348  zero_lfvector(data->F, numverts);
1349  init_bfmatrix(data->dFdX, ZERO);
1350  init_bfmatrix(data->dFdV, ZERO);
1351 
1352  data->num_blocks = 0;
1353 }
1354 
1356  int index,
1357  const float acceleration[3],
1358  const float omega[3],
1359  const float domega_dt[3],
1360  float mass)
1361 {
1362 # ifdef CLOTH_ROOT_FRAME
1363  float acc[3], w[3], dwdt[3];
1364  float f[3], dfdx[3][3], dfdv[3][3];
1365  float euler[3], coriolis[3], centrifugal[3], rotvel[3];
1366  float deuler[3][3], dcoriolis[3][3], dcentrifugal[3][3], drotvel[3][3];
1367 
1368  world_to_root_v3(data, index, acc, acceleration);
1369  world_to_root_v3(data, index, w, omega);
1370  world_to_root_v3(data, index, dwdt, domega_dt);
1371 
1372  cross_v3_v3v3(euler, dwdt, data->X[index]);
1373  cross_v3_v3v3(coriolis, w, data->V[index]);
1374  mul_v3_fl(coriolis, 2.0f);
1375  cross_v3_v3v3(rotvel, w, data->X[index]);
1376  cross_v3_v3v3(centrifugal, w, rotvel);
1377 
1378  sub_v3_v3v3(f, acc, euler);
1379  sub_v3_v3(f, coriolis);
1380  sub_v3_v3(f, centrifugal);
1381 
1382  mul_v3_fl(f, mass); /* F = m * a */
1383 
1384  cross_v3_identity(deuler, dwdt);
1385  cross_v3_identity(dcoriolis, w);
1386  mul_m3_fl(dcoriolis, 2.0f);
1387  cross_v3_identity(drotvel, w);
1388  cross_m3_v3m3(dcentrifugal, w, drotvel);
1389 
1390  add_m3_m3m3(dfdx, deuler, dcentrifugal);
1391  negate_m3(dfdx);
1392  mul_m3_fl(dfdx, mass);
1393 
1394  copy_m3_m3(dfdv, dcoriolis);
1395  negate_m3(dfdv);
1396  mul_m3_fl(dfdv, mass);
1397 
1398  add_v3_v3(data->F[index], f);
1399  add_m3_m3m3(data->dFdX[index].m, data->dFdX[index].m, dfdx);
1400  add_m3_m3m3(data->dFdV[index].m, data->dFdV[index].m, dfdv);
1401 # else
1402  (void)data;
1403  (void)index;
1404  (void)acceleration;
1405  (void)omega;
1406  (void)domega_dt;
1407 # endif
1408 }
1409 
1410 void SIM_mass_spring_force_gravity(Implicit_Data *data, int index, float mass, const float g[3])
1411 {
1412  /* force = mass * acceleration (in this case: gravity) */
1413  float f[3];
1414  world_to_root_v3(data, index, f, g);
1415  mul_v3_fl(f, mass);
1416 
1417  add_v3_v3(data->F[index], f);
1418 }
1419 
1421 {
1422  int i, numverts = data->M[0].vcount;
1423  for (i = 0; i < numverts; i++) {
1424  float tmp[3][3];
1425 
1426  /* NOTE: Uses root space velocity, no need to transform. */
1427  madd_v3_v3fl(data->F[i], data->V[i], -drag);
1428 
1429  copy_m3_m3(tmp, I);
1430  mul_m3_fl(tmp, -drag);
1431  add_m3_m3m3(data->dFdV[i].m, data->dFdV[i].m, tmp);
1432  }
1433 }
1434 
1436  struct Implicit_Data *data, int i, const float f[3], float dfdx[3][3], float dfdv[3][3])
1437 {
1438  float tf[3], tdfdx[3][3], tdfdv[3][3];
1439  world_to_root_v3(data, i, tf, f);
1440  world_to_root_m3(data, i, tdfdx, dfdx);
1441  world_to_root_m3(data, i, tdfdv, dfdv);
1442 
1443  add_v3_v3(data->F[i], tf);
1444  add_m3_m3m3(data->dFdX[i].m, data->dFdX[i].m, tdfdx);
1445  add_m3_m3m3(data->dFdV[i].m, data->dFdV[i].m, tdfdv);
1446 }
1447 
1448 static float calc_nor_area_tri(float nor[3],
1449  const float v1[3],
1450  const float v2[3],
1451  const float v3[3])
1452 {
1453  float n1[3], n2[3];
1454 
1455  sub_v3_v3v3(n1, v1, v2);
1456  sub_v3_v3v3(n2, v2, v3);
1457 
1458  cross_v3_v3v3(nor, n1, n2);
1459  return normalize_v3(nor) / 2.0f;
1460 }
1461 
1463  Implicit_Data *data, int v1, int v2, int v3, const float (*winvec)[3])
1464 {
1465  /* XXX does not support force jacobians yet,
1466  * since the effector system does not provide them either. */
1467 
1468  const float effector_scale = 0.02f;
1469  const int vs[3] = {v1, v2, v3};
1470  float win[3], nor[3], area;
1471  float factor, base_force;
1472  float force[3];
1473 
1474  /* calculate face normal and area */
1475  area = calc_nor_area_tri(nor, data->X[v1], data->X[v2], data->X[v3]);
1476  /* The force is calculated and split up evenly for each of the three face verts */
1477  factor = effector_scale * area / 3.0f;
1478 
1479  /* Calculate wind pressure at each vertex by projecting the wind field on the normal. */
1480  for (int i = 0; i < 3; i++) {
1481  world_to_root_v3(data, vs[i], win, winvec[vs[i]]);
1482 
1483  force[i] = dot_v3v3(win, nor);
1484  }
1485 
1486  /* Compute per-vertex force values from local pressures.
1487  * From integrating the pressure over the triangle and deriving
1488  * equivalent vertex forces, it follows that:
1489  *
1490  * force[idx] = (sum(pressure) + pressure[idx]) * area / 12
1491  *
1492  * Effectively, 1/4 of the pressure acts just on its vertex,
1493  * while 3/4 is split evenly over all three.
1494  */
1495  mul_v3_fl(force, factor / 4.0f);
1496 
1497  base_force = force[0] + force[1] + force[2];
1498 
1499  /* add pressure to each of the face verts */
1500  madd_v3_v3fl(data->F[v1], nor, base_force + force[0]);
1501  madd_v3_v3fl(data->F[v2], nor, base_force + force[1]);
1502  madd_v3_v3fl(data->F[v3], nor, base_force + force[2]);
1503 }
1504 
1506  Implicit_Data *data, int v1, int v2, int v3, const float (*forcevec)[3])
1507 {
1508  const float effector_scale = 0.02f;
1509  const int vs[3] = {v1, v2, v3};
1510  float nor[3], area;
1511  float factor, base_force[3];
1512  float force[3][3];
1513 
1514  /* calculate face normal and area */
1515  area = calc_nor_area_tri(nor, data->X[v1], data->X[v2], data->X[v3]);
1516  /* The force is calculated and split up evenly for each of the three face verts */
1517  factor = effector_scale * area / 3.0f;
1518 
1519  /* Compute common and per-vertex force vectors from the original inputs. */
1520  zero_v3(base_force);
1521 
1522  for (int i = 0; i < 3; i++) {
1523  world_to_root_v3(data, vs[i], force[i], forcevec[vs[i]]);
1524 
1525  mul_v3_fl(force[i], factor / 4.0f);
1526  add_v3_v3(base_force, force[i]);
1527  }
1528 
1529  /* Apply the common and vertex components to all vertices. */
1530  for (int i = 0; i < 3; i++) {
1531  add_v3_v3(force[i], base_force);
1532  add_v3_v3(data->F[vs[i]], force[i]);
1533  }
1534 }
1535 
1537 {
1538  /* The result will be 6x the volume */
1539  return volume_tri_tetrahedron_signed_v3_6x(data->X[v1], data->X[v2], data->X[v3]);
1540 }
1541 
1542 float SIM_tri_area(struct Implicit_Data *data, int v1, int v2, int v3)
1543 {
1544  float nor[3];
1545 
1546  return calc_nor_area_tri(nor, data->X[v1], data->X[v2], data->X[v3]);
1547 }
1548 
1550  int v1,
1551  int v2,
1552  int v3,
1553  float common_pressure,
1554  const float *vertex_pressure,
1555  const float weights[3])
1556 {
1557  float nor[3], area;
1558  float factor, base_force;
1559  float force[3];
1560 
1561  /* calculate face normal and area */
1562  area = calc_nor_area_tri(nor, data->X[v1], data->X[v2], data->X[v3]);
1563  /* The force is calculated and split up evenly for each of the three face verts */
1564  factor = area / 3.0f;
1565  base_force = common_pressure * factor;
1566 
1567  /* Compute per-vertex force values from local pressures.
1568  * From integrating the pressure over the triangle and deriving
1569  * equivalent vertex forces, it follows that:
1570  *
1571  * force[idx] = (sum(pressure) + pressure[idx]) * area / 12
1572  *
1573  * Effectively, 1/4 of the pressure acts just on its vertex,
1574  * while 3/4 is split evenly over all three.
1575  */
1576  if (vertex_pressure) {
1577  copy_v3_fl3(force, vertex_pressure[v1], vertex_pressure[v2], vertex_pressure[v3]);
1578  mul_v3_fl(force, factor / 4.0f);
1579 
1580  base_force += force[0] + force[1] + force[2];
1581  }
1582  else {
1583  zero_v3(force);
1584  }
1585 
1586  /* add pressure to each of the face verts */
1587  madd_v3_v3fl(data->F[v1], nor, (base_force + force[0]) * weights[0]);
1588  madd_v3_v3fl(data->F[v2], nor, (base_force + force[1]) * weights[1]);
1589  madd_v3_v3fl(data->F[v3], nor, (base_force + force[2]) * weights[2]);
1590 }
1591 
1592 static void edge_wind_vertex(const float dir[3],
1593  float length,
1594  float radius,
1595  const float wind[3],
1596  float f[3],
1597  float UNUSED(dfdx[3][3]),
1598  float UNUSED(dfdv[3][3]))
1599 {
1600  const float density = 0.01f; /* XXX arbitrary value, corresponds to effect of air density */
1601  float cos_alpha, sin_alpha, cross_section;
1602  float windlen = len_v3(wind);
1603 
1604  if (windlen == 0.0f) {
1605  zero_v3(f);
1606  return;
1607  }
1608 
1609  /* angle of wind direction to edge */
1610  cos_alpha = dot_v3v3(wind, dir) / windlen;
1611  sin_alpha = sqrtf(1.0f - cos_alpha * cos_alpha);
1612  cross_section = radius * ((float)M_PI * radius * sin_alpha + length * cos_alpha);
1613 
1614  mul_v3_v3fl(f, wind, density * cross_section);
1615 }
1616 
1618  Implicit_Data *data, int v1, int v2, float radius1, float radius2, const float (*winvec)[3])
1619 {
1620  float win[3], dir[3], length;
1621  float f[3], dfdx[3][3], dfdv[3][3];
1622 
1623  sub_v3_v3v3(dir, data->X[v1], data->X[v2]);
1624  length = normalize_v3(dir);
1625 
1626  world_to_root_v3(data, v1, win, winvec[v1]);
1627  edge_wind_vertex(dir, length, radius1, win, f, dfdx, dfdv);
1628  add_v3_v3(data->F[v1], f);
1629 
1630  world_to_root_v3(data, v2, win, winvec[v2]);
1631  edge_wind_vertex(dir, length, radius2, win, f, dfdx, dfdv);
1632  add_v3_v3(data->F[v2], f);
1633 }
1634 
1636  int v,
1637  float UNUSED(radius),
1638  const float (*winvec)[3])
1639 {
1640  const float density = 0.01f; /* XXX arbitrary value, corresponds to effect of air density */
1641 
1642  float wind[3];
1643  float f[3];
1644 
1645  world_to_root_v3(data, v, wind, winvec[v]);
1646  mul_v3_v3fl(f, wind, density);
1647  add_v3_v3(data->F[v], f);
1648 }
1649 
1650 BLI_INLINE void dfdx_spring(float to[3][3], const float dir[3], float length, float L, float k)
1651 {
1652  /* dir is unit length direction, rest is spring's restlength, k is spring constant. */
1653  // return ( (I-outerprod(dir, dir))*Min(1.0f, rest/length) - I) * -k;
1654  outerproduct(to, dir, dir);
1655  sub_m3_m3m3(to, I, to);
1656 
1657  mul_m3_fl(to, (L / length));
1658  sub_m3_m3m3(to, to, I);
1659  mul_m3_fl(to, k);
1660 }
1661 
1662 /* unused */
1663 # if 0
1664 BLI_INLINE void dfdx_damp(float to[3][3],
1665  const float dir[3],
1666  float length,
1667  const float vel[3],
1668  float rest,
1669  float damping)
1670 {
1671  /* Inner spring damping `vel` is the relative velocity of the endpoints. */
1672  // return (I - outerprod(dir, dir)) * (-damping * -(dot(dir, vel) / Max(length, rest)));
1673  mul_fvectorT_fvector(to, dir, dir);
1674  sub_fmatrix_fmatrix(to, I, to);
1675  mul_fmatrix_S(to, (-damping * -(dot_v3v3(dir, vel) / MAX2(length, rest))));
1676 }
1677 # endif
1678 
1679 BLI_INLINE void dfdv_damp(float to[3][3], const float dir[3], float damping)
1680 {
1681  /* Derivative of force with regards to velocity. */
1682  outerproduct(to, dir, dir);
1683  mul_m3_fl(to, -damping);
1684 }
1685 
1686 BLI_INLINE float fb(float length, float L)
1687 {
1688  float x = length / L;
1689  float xx = x * x;
1690  float xxx = xx * x;
1691  float xxxx = xxx * x;
1692  return (-11.541f * xxxx + 34.193f * xxx - 39.083f * xx + 23.116f * x - 9.713f);
1693 }
1694 
1695 BLI_INLINE float fbderiv(float length, float L)
1696 {
1697  float x = length / L;
1698  float xx = x * x;
1699  float xxx = xx * x;
1700  return (-46.164f * xxx + 102.579f * xx - 78.166f * x + 23.116f);
1701 }
1702 
1703 BLI_INLINE float fbstar(float length, float L, float kb, float cb)
1704 {
1705  float tempfb_fl = kb * fb(length, L);
1706  float fbstar_fl = cb * (length - L);
1707 
1708  if (tempfb_fl < fbstar_fl) {
1709  return fbstar_fl;
1710  }
1711 
1712  return tempfb_fl;
1713 }
1714 
1715 /* Function to calculate bending spring force (taken from Choi & Co). */
1716 BLI_INLINE float fbstar_jacobi(float length, float L, float kb, float cb)
1717 {
1718  float tempfb_fl = kb * fb(length, L);
1719  float fbstar_fl = cb * (length - L);
1720 
1721  if (tempfb_fl < fbstar_fl) {
1722  return -cb;
1723  }
1724 
1725  return -kb * fbderiv(length, L);
1726 }
1727 
1728 /* calculate elongation */
1730  int i,
1731  int j,
1732  float r_extent[3],
1733  float r_dir[3],
1734  float *r_length,
1735  float r_vel[3])
1736 {
1737  sub_v3_v3v3(r_extent, data->X[j], data->X[i]);
1738  sub_v3_v3v3(r_vel, data->V[j], data->V[i]);
1739  *r_length = len_v3(r_extent);
1740 
1741  if (*r_length > ALMOST_ZERO) {
1742 # if 0
1743  if (length > L) {
1744  if ((clmd->sim_parms->flags & CSIMSETT_FLAG_TEARING_ENABLED) &&
1745  (((length - L) * 100.0f / L) > clmd->sim_parms->maxspringlen)) {
1746  /* cut spring! */
1747  s->flags |= CSPRING_FLAG_DEACTIVATE;
1748  return false;
1749  }
1750  }
1751 # endif
1752  mul_v3_v3fl(r_dir, r_extent, 1.0f / (*r_length));
1753  }
1754  else {
1755  zero_v3(r_dir);
1756  }
1757 
1758  return true;
1759 }
1760 
1762  int i,
1763  int j,
1764  const float f[3],
1765  const float dfdx[3][3],
1766  const float dfdv[3][3])
1767 {
1768  int block_ij = SIM_mass_spring_add_block(data, i, j);
1769 
1770  add_v3_v3(data->F[i], f);
1771  sub_v3_v3(data->F[j], f);
1772 
1773  add_m3_m3m3(data->dFdX[i].m, data->dFdX[i].m, dfdx);
1774  add_m3_m3m3(data->dFdX[j].m, data->dFdX[j].m, dfdx);
1775  sub_m3_m3m3(data->dFdX[block_ij].m, data->dFdX[block_ij].m, dfdx);
1776 
1777  add_m3_m3m3(data->dFdV[i].m, data->dFdV[i].m, dfdv);
1778  add_m3_m3m3(data->dFdV[j].m, data->dFdV[j].m, dfdv);
1779  sub_m3_m3m3(data->dFdV[block_ij].m, data->dFdV[block_ij].m, dfdv);
1780 }
1781 
1783  int i,
1784  int j,
1785  float restlen,
1786  float stiffness_tension,
1787  float damping_tension,
1788  float stiffness_compression,
1789  float damping_compression,
1790  bool resist_compress,
1791  bool new_compress,
1792  float clamp_force)
1793 {
1794  float extent[3], length, dir[3], vel[3];
1795  float f[3], dfdx[3][3], dfdv[3][3];
1796  float damping = 0;
1797 
1798  /* calculate elongation */
1799  spring_length(data, i, j, extent, dir, &length, vel);
1800 
1801  /* This code computes not only the force, but also its derivative.
1802  * Zero derivative effectively disables the spring for the implicit solver.
1803  * Thus length > restlen makes cloth unconstrained at the start of simulation. */
1804  if ((length >= restlen && length > 0) || resist_compress) {
1805  float stretch_force;
1806 
1807  damping = damping_tension;
1808 
1809  stretch_force = stiffness_tension * (length - restlen);
1810  if (clamp_force > 0.0f && stretch_force > clamp_force) {
1811  stretch_force = clamp_force;
1812  }
1813  mul_v3_v3fl(f, dir, stretch_force);
1814 
1815  dfdx_spring(dfdx, dir, length, restlen, stiffness_tension);
1816  }
1817  else if (new_compress) {
1818  /* This is based on the Choi and Ko bending model,
1819  * which works surprisingly well for compression. */
1820  float kb = stiffness_compression;
1821  float cb = kb; /* cb equal to kb seems to work, but a factor can be added if necessary */
1822 
1823  damping = damping_compression;
1824 
1825  mul_v3_v3fl(f, dir, fbstar(length, restlen, kb, cb));
1826 
1827  outerproduct(dfdx, dir, dir);
1828  mul_m3_fl(dfdx, fbstar_jacobi(length, restlen, kb, cb));
1829  }
1830  else {
1831  return false;
1832  }
1833 
1834  madd_v3_v3fl(f, dir, damping * dot_v3v3(vel, dir));
1835  dfdv_damp(dfdv, dir, damping);
1836 
1837  apply_spring(data, i, j, f, dfdx, dfdv);
1838 
1839  return true;
1840 }
1841 
1843  Implicit_Data *data, int i, int j, float restlen, float kb, float cb)
1844 {
1845  /* See "Stable but Responsive Cloth" (Choi, Ko 2005). */
1846 
1847  float extent[3], length, dir[3], vel[3];
1848 
1849  /* calculate elongation */
1850  spring_length(data, i, j, extent, dir, &length, vel);
1851 
1852  if (length < restlen) {
1853  float f[3], dfdx[3][3], dfdv[3][3];
1854 
1855  mul_v3_v3fl(f, dir, fbstar(length, restlen, kb, cb));
1856 
1857  outerproduct(dfdx, dir, dir);
1858  mul_m3_fl(dfdx, fbstar_jacobi(length, restlen, kb, cb));
1859 
1860  /* XXX damping not supported */
1861  zero_m3(dfdv);
1862 
1863  apply_spring(data, i, j, f, dfdx, dfdv);
1864 
1865  return true;
1866  }
1867 
1868  return false;
1869 }
1870 
1871 BLI_INLINE void poly_avg(lfVector *data, const int *inds, int len, float r_avg[3])
1872 {
1873  float fact = 1.0f / (float)len;
1874 
1875  zero_v3(r_avg);
1876 
1877  for (int i = 0; i < len; i++) {
1878  madd_v3_v3fl(r_avg, data[inds[i]], fact);
1879  }
1880 }
1881 
1882 BLI_INLINE void poly_norm(lfVector *data, int i, int j, int *inds, int len, float r_dir[3])
1883 {
1884  float mid[3];
1885 
1886  poly_avg(data, inds, len, mid);
1887 
1888  normal_tri_v3(r_dir, data[i], data[j], mid);
1889 }
1890 
1891 BLI_INLINE void edge_avg(lfVector *data, int i, int j, float r_avg[3])
1892 {
1893  r_avg[0] = (data[i][0] + data[j][0]) * 0.5f;
1894  r_avg[1] = (data[i][1] + data[j][1]) * 0.5f;
1895  r_avg[2] = (data[i][2] + data[j][2]) * 0.5f;
1896 }
1897 
1898 BLI_INLINE void edge_norm(lfVector *data, int i, int j, float r_dir[3])
1899 {
1900  sub_v3_v3v3(r_dir, data[i], data[j]);
1901  normalize_v3(r_dir);
1902 }
1903 
1904 BLI_INLINE float bend_angle(const float dir_a[3], const float dir_b[3], const float dir_e[3])
1905 {
1906  float cos, sin;
1907  float tmp[3];
1908 
1909  cos = dot_v3v3(dir_a, dir_b);
1910 
1911  cross_v3_v3v3(tmp, dir_a, dir_b);
1912  sin = dot_v3v3(tmp, dir_e);
1913 
1914  return atan2f(sin, cos);
1915 }
1916 
1918  int i,
1919  int j,
1920  int *i_a,
1921  int *i_b,
1922  int len_a,
1923  int len_b,
1924  float r_dir_a[3],
1925  float r_dir_b[3],
1926  float *r_angle,
1927  float r_vel_a[3],
1928  float r_vel_b[3])
1929 {
1930  float dir_e[3], vel_e[3];
1931 
1932  poly_norm(data->X, j, i, i_a, len_a, r_dir_a);
1933  poly_norm(data->X, i, j, i_b, len_b, r_dir_b);
1934 
1935  edge_norm(data->X, i, j, dir_e);
1936 
1937  *r_angle = bend_angle(r_dir_a, r_dir_b, dir_e);
1938 
1939  poly_avg(data->V, i_a, len_a, r_vel_a);
1940  poly_avg(data->V, i_b, len_b, r_vel_b);
1941 
1942  edge_avg(data->V, i, j, vel_e);
1943 
1944  sub_v3_v3(r_vel_a, vel_e);
1945  sub_v3_v3(r_vel_b, vel_e);
1946 }
1947 
1949  int i,
1950  int j,
1951  int *i_a,
1952  int *i_b,
1953  int len_a,
1954  int len_b,
1955  float restang,
1956  float stiffness,
1957  float damping)
1958 {
1959  float angle, dir_a[3], dir_b[3], vel_a[3], vel_b[3];
1960  float f_a[3], f_b[3], f_e[3];
1961  float force;
1962  int x;
1963 
1964  spring_angle(data, i, j, i_a, i_b, len_a, len_b, dir_a, dir_b, &angle, vel_a, vel_b);
1965 
1966  /* spring force */
1967  force = stiffness * (angle - restang);
1968 
1969  /* damping force */
1970  force += -damping * (dot_v3v3(vel_a, dir_a) + dot_v3v3(vel_b, dir_b));
1971 
1972  mul_v3_v3fl(f_a, dir_a, force / len_a);
1973  mul_v3_v3fl(f_b, dir_b, force / len_b);
1974 
1975  for (x = 0; x < len_a; x++) {
1976  add_v3_v3(data->F[i_a[x]], f_a);
1977  }
1978 
1979  for (x = 0; x < len_b; x++) {
1980  add_v3_v3(data->F[i_b[x]], f_b);
1981  }
1982 
1983  mul_v3_v3fl(f_a, dir_a, force * 0.5f);
1984  mul_v3_v3fl(f_b, dir_b, force * 0.5f);
1985 
1986  add_v3_v3v3(f_e, f_a, f_b);
1987 
1988  sub_v3_v3(data->F[i], f_e);
1989  sub_v3_v3(data->F[j], f_e);
1990 
1991  return true;
1992 }
1993 
1994 /* Jacobian of a direction vector.
1995  * Basically the part of the differential orthogonal to the direction,
1996  * inversely proportional to the length of the edge.
1997  *
1998  * dD_ij/dx_i = -dD_ij/dx_j = (D_ij * D_ij^T - I) / len_ij
1999  */
2001  Implicit_Data *data, int i, int j, float edge[3], float dir[3], float grad_dir[3][3])
2002 {
2003  float length;
2004 
2005  sub_v3_v3v3(edge, data->X[j], data->X[i]);
2006  length = normalize_v3_v3(dir, edge);
2007 
2008  if (length > ALMOST_ZERO) {
2009  outerproduct(grad_dir, dir, dir);
2010  sub_m3_m3m3(grad_dir, I, grad_dir);
2011  mul_m3_fl(grad_dir, 1.0f / length);
2012  }
2013  else {
2014  zero_m3(grad_dir);
2015  }
2016 }
2017 
2019  int i,
2020  int j,
2021  int k,
2022  const float goal[3],
2023  float stiffness,
2024  float damping,
2025  int q,
2026  const float dx[3],
2027  const float dv[3],
2028  float r_f[3])
2029 {
2030  float edge_ij[3], dir_ij[3];
2031  float edge_jk[3], dir_jk[3];
2032  float vel_ij[3], vel_jk[3], vel_ortho[3];
2033  float f_bend[3], f_damp[3];
2034  float fk[3];
2035  float dist[3];
2036 
2037  zero_v3(fk);
2038 
2039  sub_v3_v3v3(edge_ij, data->X[j], data->X[i]);
2040  if (q == i) {
2041  sub_v3_v3(edge_ij, dx);
2042  }
2043  if (q == j) {
2044  add_v3_v3(edge_ij, dx);
2045  }
2046  normalize_v3_v3(dir_ij, edge_ij);
2047 
2048  sub_v3_v3v3(edge_jk, data->X[k], data->X[j]);
2049  if (q == j) {
2050  sub_v3_v3(edge_jk, dx);
2051  }
2052  if (q == k) {
2053  add_v3_v3(edge_jk, dx);
2054  }
2055  normalize_v3_v3(dir_jk, edge_jk);
2056 
2057  sub_v3_v3v3(vel_ij, data->V[j], data->V[i]);
2058  if (q == i) {
2059  sub_v3_v3(vel_ij, dv);
2060  }
2061  if (q == j) {
2062  add_v3_v3(vel_ij, dv);
2063  }
2064 
2065  sub_v3_v3v3(vel_jk, data->V[k], data->V[j]);
2066  if (q == j) {
2067  sub_v3_v3(vel_jk, dv);
2068  }
2069  if (q == k) {
2070  add_v3_v3(vel_jk, dv);
2071  }
2072 
2073  /* bending force */
2074  sub_v3_v3v3(dist, goal, edge_jk);
2075  mul_v3_v3fl(f_bend, dist, stiffness);
2076 
2077  add_v3_v3(fk, f_bend);
2078 
2079  /* damping force */
2080  madd_v3_v3v3fl(vel_ortho, vel_jk, dir_jk, -dot_v3v3(vel_jk, dir_jk));
2081  mul_v3_v3fl(f_damp, vel_ortho, damping);
2082 
2083  sub_v3_v3(fk, f_damp);
2084 
2085  copy_v3_v3(r_f, fk);
2086 }
2087 
2088 /* Finite Differences method for estimating the jacobian of the force */
2090  int i,
2091  int j,
2092  int k,
2093  const float goal[3],
2094  float stiffness,
2095  float damping,
2096  int q,
2097  float dfdx[3][3])
2098 {
2099  const float delta = 0.00001f; /* TODO: find a good heuristic for this. */
2100  float dvec_null[3][3], dvec_pos[3][3], dvec_neg[3][3];
2101  float f[3];
2102  int a, b;
2103 
2104  zero_m3(dvec_null);
2105  unit_m3(dvec_pos);
2106  mul_m3_fl(dvec_pos, delta * 0.5f);
2107  copy_m3_m3(dvec_neg, dvec_pos);
2108  negate_m3(dvec_neg);
2109 
2110  /* XXX TODO: offset targets to account for position dependency. */
2111 
2112  for (a = 0; a < 3; a++) {
2114  data, i, j, k, goal, stiffness, damping, q, dvec_pos[a], dvec_null[a], f);
2115  copy_v3_v3(dfdx[a], f);
2116 
2118  data, i, j, k, goal, stiffness, damping, q, dvec_neg[a], dvec_null[a], f);
2119  sub_v3_v3(dfdx[a], f);
2120 
2121  for (b = 0; b < 3; b++) {
2122  dfdx[a][b] /= delta;
2123  }
2124  }
2125 }
2126 
2127 /* Finite Differences method for estimating the jacobian of the force */
2129  int i,
2130  int j,
2131  int k,
2132  const float goal[3],
2133  float stiffness,
2134  float damping,
2135  int q,
2136  float dfdv[3][3])
2137 {
2138  const float delta = 0.00001f; /* TODO: find a good heuristic for this. */
2139  float dvec_null[3][3], dvec_pos[3][3], dvec_neg[3][3];
2140  float f[3];
2141  int a, b;
2142 
2143  zero_m3(dvec_null);
2144  unit_m3(dvec_pos);
2145  mul_m3_fl(dvec_pos, delta * 0.5f);
2146  copy_m3_m3(dvec_neg, dvec_pos);
2147  negate_m3(dvec_neg);
2148 
2149  /* XXX TODO: offset targets to account for position dependency. */
2150 
2151  for (a = 0; a < 3; a++) {
2153  data, i, j, k, goal, stiffness, damping, q, dvec_null[a], dvec_pos[a], f);
2154  copy_v3_v3(dfdv[a], f);
2155 
2157  data, i, j, k, goal, stiffness, damping, q, dvec_null[a], dvec_neg[a], f);
2158  sub_v3_v3(dfdv[a], f);
2159 
2160  for (b = 0; b < 3; b++) {
2161  dfdv[a][b] /= delta;
2162  }
2163  }
2164 }
2165 
2167  int i,
2168  int j,
2169  int k,
2170  const float target[3],
2171  float stiffness,
2172  float damping)
2173 {
2174  /* Angular springs roughly based on the bending model proposed by Baraff and Witkin in
2175  * "Large Steps in Cloth Simulation". */
2176 
2177  float goal[3];
2178  float fj[3], fk[3];
2179  float dfj_dxi[3][3], dfj_dxj[3][3], dfk_dxi[3][3], dfk_dxj[3][3], dfk_dxk[3][3];
2180  float dfj_dvi[3][3], dfj_dvj[3][3], dfk_dvi[3][3], dfk_dvj[3][3], dfk_dvk[3][3];
2181 
2182  const float vecnull[3] = {0.0f, 0.0f, 0.0f};
2183 
2184  int block_ij = SIM_mass_spring_add_block(data, i, j);
2185  int block_jk = SIM_mass_spring_add_block(data, j, k);
2186  int block_ik = SIM_mass_spring_add_block(data, i, k);
2187 
2188  world_to_root_v3(data, j, goal, target);
2189 
2190  spring_hairbend_forces(data, i, j, k, goal, stiffness, damping, k, vecnull, vecnull, fk);
2191  negate_v3_v3(fj, fk); /* Counter-force. */
2192 
2193  spring_hairbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, i, dfk_dxi);
2194  spring_hairbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, j, dfk_dxj);
2195  spring_hairbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, k, dfk_dxk);
2196  copy_m3_m3(dfj_dxi, dfk_dxi);
2197  negate_m3(dfj_dxi);
2198  copy_m3_m3(dfj_dxj, dfk_dxj);
2199  negate_m3(dfj_dxj);
2200 
2201  spring_hairbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, i, dfk_dvi);
2202  spring_hairbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, j, dfk_dvj);
2203  spring_hairbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, k, dfk_dvk);
2204  copy_m3_m3(dfj_dvi, dfk_dvi);
2205  negate_m3(dfj_dvi);
2206  copy_m3_m3(dfj_dvj, dfk_dvj);
2207  negate_m3(dfj_dvj);
2208 
2209  /* add forces and jacobians to the solver data */
2210 
2211  add_v3_v3(data->F[j], fj);
2212  add_v3_v3(data->F[k], fk);
2213 
2214  add_m3_m3m3(data->dFdX[j].m, data->dFdX[j].m, dfj_dxj);
2215  add_m3_m3m3(data->dFdX[k].m, data->dFdX[k].m, dfk_dxk);
2216 
2217  add_m3_m3m3(data->dFdX[block_ij].m, data->dFdX[block_ij].m, dfj_dxi);
2218  add_m3_m3m3(data->dFdX[block_jk].m, data->dFdX[block_jk].m, dfk_dxj);
2219  add_m3_m3m3(data->dFdX[block_ik].m, data->dFdX[block_ik].m, dfk_dxi);
2220 
2221  add_m3_m3m3(data->dFdV[j].m, data->dFdV[j].m, dfj_dvj);
2222  add_m3_m3m3(data->dFdV[k].m, data->dFdV[k].m, dfk_dvk);
2223 
2224  add_m3_m3m3(data->dFdV[block_ij].m, data->dFdV[block_ij].m, dfj_dvi);
2225  add_m3_m3m3(data->dFdV[block_jk].m, data->dFdV[block_jk].m, dfk_dvj);
2226  add_m3_m3m3(data->dFdV[block_ik].m, data->dFdV[block_ik].m, dfk_dvi);
2227 
2228  /* XXX analytical calculation of derivatives below is incorrect.
2229  * This proved to be difficult, but for now just using the finite difference method for
2230  * estimating the jacobians should be sufficient.
2231  */
2232 # if 0
2233  float edge_ij[3], dir_ij[3], grad_dir_ij[3][3];
2234  float edge_jk[3], dir_jk[3], grad_dir_jk[3][3];
2235  float dist[3], vel_jk[3], vel_jk_ortho[3], projvel[3];
2236  float target[3];
2237  float tmp[3][3];
2238  float fi[3], fj[3], fk[3];
2239  float dfi_dxi[3][3], dfj_dxi[3][3], dfj_dxj[3][3], dfk_dxi[3][3], dfk_dxj[3][3], dfk_dxk[3][3];
2240  float dfdvi[3][3];
2241 
2242  /* TESTING */
2243  damping = 0.0f;
2244 
2245  zero_v3(fi);
2246  zero_v3(fj);
2247  zero_v3(fk);
2248  zero_m3(dfi_dxi);
2249  zero_m3(dfj_dxi);
2250  zero_m3(dfk_dxi);
2251  zero_m3(dfk_dxj);
2252  zero_m3(dfk_dxk);
2253 
2254  /* jacobian of direction vectors */
2255  spring_grad_dir(data, i, j, edge_ij, dir_ij, grad_dir_ij);
2256  spring_grad_dir(data, j, k, edge_jk, dir_jk, grad_dir_jk);
2257 
2258  sub_v3_v3v3(vel_jk, data->V[k], data->V[j]);
2259 
2260  /* bending force */
2261  mul_v3_v3fl(target, dir_ij, restlen);
2262  sub_v3_v3v3(dist, target, edge_jk);
2263  mul_v3_v3fl(fk, dist, stiffness);
2264 
2265  /* damping force */
2266  madd_v3_v3v3fl(vel_jk_ortho, vel_jk, dir_jk, -dot_v3v3(vel_jk, dir_jk));
2267  madd_v3_v3fl(fk, vel_jk_ortho, damping);
2268 
2269  /* XXX this only holds true as long as we assume straight rest shape!
2270  * eventually will become a bit more involved since the opposite segment
2271  * gets its own target, under condition of having equal torque on both sides.
2272  */
2273  copy_v3_v3(fi, fk);
2274 
2275  /* counterforce on the middle point */
2276  sub_v3_v3(fj, fi);
2277  sub_v3_v3(fj, fk);
2278 
2279  /* === derivatives === */
2280 
2281  madd_m3_m3fl(dfk_dxi, grad_dir_ij, stiffness * restlen);
2282 
2283  madd_m3_m3fl(dfk_dxj, grad_dir_ij, -stiffness * restlen);
2284  madd_m3_m3fl(dfk_dxj, I, stiffness);
2285 
2286  madd_m3_m3fl(dfk_dxk, I, -stiffness);
2287 
2288  copy_m3_m3(dfi_dxi, dfk_dxk);
2289  negate_m3(dfi_dxi);
2290 
2291  /* dfj_dfi == dfi_dfj due to symmetry,
2292  * dfi_dfj == dfk_dfj due to fi == fk
2293  * XXX see comment above on future bent rest shapes
2294  */
2295  copy_m3_m3(dfj_dxi, dfk_dxj);
2296 
2297  /* dfj_dxj == -(dfi_dxj + dfk_dxj) due to fj == -(fi + fk) */
2298  sub_m3_m3m3(dfj_dxj, dfj_dxj, dfj_dxi);
2299  sub_m3_m3m3(dfj_dxj, dfj_dxj, dfk_dxj);
2300 
2301  /* add forces and jacobians to the solver data */
2302  add_v3_v3(data->F[i], fi);
2303  add_v3_v3(data->F[j], fj);
2304  add_v3_v3(data->F[k], fk);
2305 
2306  add_m3_m3m3(data->dFdX[i].m, data->dFdX[i].m, dfi_dxi);
2307  add_m3_m3m3(data->dFdX[j].m, data->dFdX[j].m, dfj_dxj);
2308  add_m3_m3m3(data->dFdX[k].m, data->dFdX[k].m, dfk_dxk);
2309 
2310  add_m3_m3m3(data->dFdX[block_ij].m, data->dFdX[block_ij].m, dfj_dxi);
2311  add_m3_m3m3(data->dFdX[block_jk].m, data->dFdX[block_jk].m, dfk_dxj);
2312  add_m3_m3m3(data->dFdX[block_ik].m, data->dFdX[block_ik].m, dfk_dxi);
2313 # endif
2314 
2315  return true;
2316 }
2317 
2319  int i,
2320  const float goal_x[3],
2321  const float goal_v[3],
2322  float stiffness,
2323  float damping)
2324 {
2325  float root_goal_x[3], root_goal_v[3], extent[3], length, dir[3], vel[3];
2326  float f[3], dfdx[3][3], dfdv[3][3];
2327 
2328  /* goal is in world space */
2329  world_to_root_v3(data, i, root_goal_x, goal_x);
2330  world_to_root_v3(data, i, root_goal_v, goal_v);
2331 
2332  sub_v3_v3v3(extent, root_goal_x, data->X[i]);
2333  sub_v3_v3v3(vel, root_goal_v, data->V[i]);
2334  length = normalize_v3_v3(dir, extent);
2335 
2336  if (length > ALMOST_ZERO) {
2337  mul_v3_v3fl(f, dir, stiffness * length);
2338 
2339  /* Ascher & Boxman, p.21: Damping only during elongation
2340  * something wrong with it. */
2341  madd_v3_v3fl(f, dir, damping * dot_v3v3(vel, dir));
2342 
2343  dfdx_spring(dfdx, dir, length, 0.0f, stiffness);
2344  dfdv_damp(dfdv, dir, damping);
2345 
2346  add_v3_v3(data->F[i], f);
2347  add_m3_m3m3(data->dFdX[i].m, data->dFdX[i].m, dfdx);
2348  add_m3_m3m3(data->dFdV[i].m, data->dFdV[i].m, dfdv);
2349 
2350  return true;
2351  }
2352 
2353  return false;
2354 }
2355 
2356 #endif /* IMPLICIT_SOLVER_BLENDER */
typedef float(TangentPoint)[2]
#define VECADDS(v1, v2, v3, bS)
Definition: BKE_cloth.h:152
#define DO_INLINE
Definition: BKE_cloth.h:24
#define VECSUBS(v1, v2, v3, bS)
Definition: BKE_cloth.h:166
#define VECADDSS(v1, v2, aS, v3, bS)
Definition: BKE_cloth.h:145
#define VECSUBADDSS(v1, v2, aS, v3, bS)
Definition: BKE_cloth.h:138
#define VECSUBMUL(v1, v2, aS)
Definition: BKE_cloth.h:159
#define ALMOST_ZERO
Definition: BKE_cloth.h:31
#define BLI_assert(a)
Definition: BLI_assert.h:46
#define BLI_INLINE
sqrt(x)+1/max(0
#define M_PI
Definition: BLI_math_base.h:20
float volume_tri_tetrahedron_signed_v3_6x(const float v1[3], const float v2[3], const float v3[3])
Definition: math_geom.c:255
float normal_tri_v3(float n[3], const float v1[3], const float v2[3], const float v3[3])
Definition: math_geom.c:33
void sub_m3_m3m3(float R[3][3], const float A[3][3], const float B[3][3])
Definition: math_matrix.c:1076
void negate_m3(float R[3][3])
Definition: math_matrix.c:989
void mul_m3_v3(const float M[3][3], float r[3])
Definition: math_matrix.c:926
void copy_m3_m3(float m1[3][3], const float m2[3][3])
Definition: math_matrix.c:71
void unit_m3(float m[3][3])
Definition: math_matrix.c:40
void mul_m3_fl(float R[3][3], float f)
Definition: math_matrix.c:956
void add_m3_m3m3(float R[3][3], const float A[3][3], const float B[3][3])
Definition: math_matrix.c:1032
void zero_m3(float m[3][3])
Definition: math_matrix.c:23
void mul_v3_m3v3(float r[3], const float M[3][3], const float a[3])
Definition: math_matrix.c:897
void transpose_m3(float R[3][3])
Definition: math_matrix.c:1332
void mul_transposed_m3_v3(const float M[3][3], float r[3])
Definition: math_matrix.c:936
void mul_m3_m3m3(float R[3][3], const float A[3][3], const float B[3][3])
Definition: math_matrix.c:388
MINLINE void madd_v3_v3fl(float r[3], const float a[3], float f)
MINLINE float normalize_v3(float r[3])
MINLINE void sub_v3_v3(float r[3], const float a[3])
MINLINE void sub_v3_v3v3(float r[3], const float a[3], const float b[3])
MINLINE void mul_v3_fl(float r[3], float f)
MINLINE void copy_v3_v3(float r[3], const float a[3])
MINLINE void negate_v3_v3(float r[3], const float a[3])
MINLINE void copy_v3_fl3(float v[3], float x, float y, float z)
MINLINE float dot_v3v3(const float a[3], const float b[3]) ATTR_WARN_UNUSED_RESULT
MINLINE void add_v3_v3v3(float r[3], const float a[3], const float b[3])
MINLINE void cross_v3_v3v3(float r[3], const float a[3], const float b[3])
MINLINE float normalize_v3_v3(float r[3], const float a[3])
MINLINE void madd_v3_v3v3fl(float r[3], const float a[3], const float b[3], float f)
MINLINE void zero_v3(float r[3])
MINLINE void mul_v3_v3fl(float r[3], const float a[3], float f)
MINLINE void add_v3_v3(float r[3], const float a[3])
MINLINE float len_v3(const float a[3]) ATTR_WARN_UNUSED_RESULT
#define UNUSED(x)
#define MAX2(a, b)
Object is a sort of wrapper for general info.
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble z
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble w _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat w _GL_VOID_RET _GL_VOID GLint GLint GLint w _GL_VOID_RET _GL_VOID GLshort GLshort GLshort w _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble y2 _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat y2 _GL_VOID_RET _GL_VOID GLint GLint GLint y2 _GL_VOID_RET _GL_VOID GLshort GLshort GLshort y2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLuint *buffer _GL_VOID_RET _GL_VOID GLdouble t _GL_VOID_RET _GL_VOID GLfloat t _GL_VOID_RET _GL_VOID GLint t _GL_VOID_RET _GL_VOID GLshort t _GL_VOID_RET _GL_VOID GLdouble GLdouble r _GL_VOID_RET _GL_VOID GLfloat GLfloat r _GL_VOID_RET _GL_VOID GLint GLint r _GL_VOID_RET _GL_VOID GLshort GLshort r _GL_VOID_RET _GL_VOID GLdouble GLdouble r
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint i1
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble u2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLdouble GLdouble v2 _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLdouble GLdouble nz _GL_VOID_RET _GL_VOID GLfloat GLfloat nz _GL_VOID_RET _GL_VOID GLint GLint nz _GL_VOID_RET _GL_VOID GLshort GLshort nz _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const GLfloat *values _GL_VOID_RET _GL_VOID GLsizei const GLushort *values _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID const GLuint const GLclampf *priorities _GL_VOID_RET _GL_VOID GLdouble y _GL_VOID_RET _GL_VOID GLfloat y _GL_VOID_RET _GL_VOID GLint y _GL_VOID_RET _GL_VOID GLshort y _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLfloat GLfloat z _GL_VOID_RET _GL_VOID GLint GLint z _GL_VOID_RET _GL_VOID GLshort GLshort z _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble w _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat w _GL_VOID_RET _GL_VOID GLint GLint GLint w _GL_VOID_RET _GL_VOID GLshort GLshort GLshort w _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble y2 _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat y2 _GL_VOID_RET _GL_VOID GLint GLint GLint y2 _GL_VOID_RET _GL_VOID GLshort GLshort GLshort y2 _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLdouble GLdouble z _GL_VOID_RET _GL_VOID GLuint *buffer _GL_VOID_RET _GL_VOID GLdouble t _GL_VOID_RET _GL_VOID GLfloat t _GL_VOID_RET _GL_VOID GLint t _GL_VOID_RET _GL_VOID GLshort t _GL_VOID_RET _GL_VOID GLdouble t
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble GLdouble GLdouble zFar _GL_VOID_RET _GL_UINT GLdouble *equation _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLenum GLfloat *v _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLfloat *values _GL_VOID_RET _GL_VOID GLushort *values _GL_VOID_RET _GL_VOID GLenum GLfloat *params _GL_VOID_RET _GL_VOID GLenum GLdouble *params _GL_VOID_RET _GL_VOID GLenum GLint *params _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_VOID GLsizei const void *pointer _GL_VOID_RET _GL_BOOL GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLenum GLfloat param _GL_VOID_RET _GL_VOID GLenum GLint param _GL_VOID_RET _GL_VOID GLushort pattern _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint const GLdouble *points _GL_VOID_RET _GL_VOID GLdouble GLdouble GLint GLint GLdouble v1
Read Guarded memory(de)allocation.
Platform independent time functions.
@ SIM_SOLVER_SUCCESS
@ SIM_SOLVER_NO_CONVERGENCE
ATTR_WARN_UNUSED_RESULT const BMVert * v2
ATTR_WARN_UNUSED_RESULT const BMLoop * l
ATTR_WARN_UNUSED_RESULT const BMVert * v
static DBVT_INLINE btScalar size(const btDbvtVolume &a)
Definition: btDbvt.cpp:52
SIMD_FORCE_INLINE const btScalar & w() const
Return the w value.
Definition: btQuadWord.h:119
SIMD_FORCE_INLINE btScalar angle(const btVector3 &v) const
Return the angle between this and another vector.
Definition: btVector3.h:356
StackEntry * from
SyclQueue void void size_t num_bytes void
int len
Definition: draw_manager.c:108
BLI_INLINE void print_lvector(const lVector3f &v)
Definition: eigen_utils.h:190
static float verts[][3]
uint nor
BLI_INLINE void implicit_print_matrix_elem(float v)
Definition: implicit.h:46
void SIM_mass_spring_force_gravity(Implicit_Data *data, int index, float mass, const float g[3])
void SIM_mass_spring_set_vertex_mass(Implicit_Data *data, int index, float mass)
void SIM_mass_spring_clear_forces(Implicit_Data *data)
BLI_INLINE void apply_spring(Implicit_Data *data, int i, int j, const float f[3], const float dfdx[3][3], const float dfdv[3][3])
DO_INLINE void add_lfvector_lfvector(float(*to)[3], float(*fLongVectorA)[3], float(*fLongVectorB)[3], unsigned int verts)
void SIM_mass_spring_clear_constraints(Implicit_Data *data)
BLI_INLINE float fbstar_jacobi(float length, float L, float kb, float cb)
void SIM_mass_spring_force_pressure(Implicit_Data *data, int v1, int v2, int v3, float common_pressure, const float *vertex_pressure, const float weights[3])
DO_INLINE void add_lfvector_lfvectorS(float(*to)[3], float(*fLongVectorA)[3], float(*fLongVectorB)[3], float bS, unsigned int verts)
DO_INLINE void add_fmatrix_fmatrix(float to[3][3], const float matrixA[3][3], const float matrixB[3][3])
void SIM_mass_spring_add_constraint_ndof0(Implicit_Data *data, int index, const float dV[3])
static void edge_wind_vertex(const float dir[3], float length, float radius, const float wind[3], float f[3], float UNUSED(dfdx[3][3]), float UNUSED(dfdv[3][3]))
DO_INLINE void mul_fvectorT_fvector(float to[3][3], const float vectorA[3], const float vectorB[3])
void SIM_mass_spring_get_motion_state(struct Implicit_Data *data, int index, float x[3], float v[3])
DO_INLINE void del_lfvector(float(*fLongVector)[3])
BLI_INLINE void cross_m3_v3m3(float r[3][3], const float v[3], const float m[3][3])
bool SIM_mass_spring_force_spring_bending_hair(Implicit_Data *data, int i, int j, int k, const float target[3], float stiffness, float damping)
BLI_INLINE float fb(float length, float L)
struct fmatrix3x3 fmatrix3x3
void SIM_mass_spring_set_velocity(Implicit_Data *data, int index, const float v[3])
DO_INLINE lfVector * create_lfvector(unsigned int verts)
DO_INLINE void mul_fvectorT_fvectorS(float to[3][3], float vectorA[3], float vectorB[3], float aS)
DO_INLINE void mul_fmatrix_S(float matrix[3][3], float scalar)
bool SIM_mass_spring_solve_positions(Implicit_Data *data, float dt)
DO_INLINE void initdiag_fmatrixS(float to[3][3], float aS)
Implicit_Data * SIM_mass_spring_solver_create(int numverts, int numsprings)
void SIM_mass_spring_set_new_velocity(struct Implicit_Data *data, int index, const float v[3])
BLI_INLINE float fbderiv(float length, float L)
BLI_INLINE void madd_m3_m3fl(float r[3][3], const float m[3][3], float f)
DO_INLINE void mul_lfvectorS(float(*to)[3], float(*fLongVector)[3], float scalar, unsigned int verts)
DO_INLINE void filter(lfVector *V, fmatrix3x3 *S)
DO_INLINE void muladd_fmatrix_fvector(float to[3], const float matrix[3][3], const float from[3])
DO_INLINE void zero_lfvector(float(*to)[3], unsigned int verts)
DO_INLINE void mul_bfmatrix_lfvector(float(*to)[3], fmatrix3x3 *from, lfVector *fLongVector)
BLI_INLINE void spring_hairbend_forces(Implicit_Data *data, int i, int j, int k, const float goal[3], float stiffness, float damping, int q, const float dx[3], const float dv[3], float r_f[3])
float SIM_tri_area(struct Implicit_Data *data, int v1, int v2, int v3)
bool SIM_mass_spring_force_spring_goal(Implicit_Data *data, int i, const float goal_x[3], const float goal_v[3], float stiffness, float damping)
void SIM_mass_spring_force_face_extern(Implicit_Data *data, int v1, int v2, int v3, const float(*forcevec)[3])
DO_INLINE void muladd_fmatrixT_fvector(float to[3], const float matrix[3][3], const float from[3])
bool SIM_mass_spring_solve_velocities(Implicit_Data *data, float dt, ImplicitSolverResult *result)
DO_INLINE void mul_fvector_fmatrix(float *to, const float *from, const float matrix[3][3])
void SIM_mass_spring_get_velocity(struct Implicit_Data *data, int index, float v[3])
BLI_INLINE void init_fmatrix(fmatrix3x3 *matrix, int r, int c)
void SIM_mass_spring_add_constraint_ndof2(Implicit_Data *data, int index, const float c1[3], const float dV[3])
BLI_INLINE bool spring_length(Implicit_Data *data, int i, int j, float r_extent[3], float r_dir[3], float *r_length, float r_vel[3])
void SIM_mass_spring_get_new_position(struct Implicit_Data *data, int index, float x[3])
static int SIM_mass_spring_add_block(Implicit_Data *data, int v1, int v2)
BLI_INLINE void edge_avg(lfVector *data, int i, int j, float r_avg[3])
void SIM_mass_spring_force_drag(Implicit_Data *data, float drag)
DO_INLINE fmatrix3x3 * create_bfmatrix(unsigned int verts, unsigned int springs)
DO_INLINE void subadd_fmatrixS_fmatrixS(float to[3][3], const float matrixA[3][3], float aS, const float matrixB[3][3], float bS)
BLI_INLINE void root_to_world_m3(Implicit_Data *data, int index, float r[3][3], const float m[3][3])
void SIM_mass_spring_solver_free(Implicit_Data *id)
BLI_INLINE void cross_v3_identity(float r[3][3], const float v[3])
DO_INLINE void init_lfvector(float(*fLongVector)[3], const float vector[3], unsigned int verts)
DO_INLINE void mul_fmatrix_fvector(float *to, const float matrix[3][3], const float from[3])
bool SIM_mass_spring_force_spring_linear(Implicit_Data *data, int i, int j, float restlen, float stiffness_tension, float damping_tension, float stiffness_compression, float damping_compression, bool resist_compress, bool new_compress, float clamp_force)
BLI_INLINE void outerproduct(float r[3][3], const float a[3], const float b[3])
bool SIM_mass_spring_force_spring_bending(Implicit_Data *data, int i, int j, float restlen, float kb, float cb)
DO_INLINE void mul_fvector_S(float to[3], const float from[3], float scalar)
DO_INLINE float dot_lfvector(float(*fLongVectorA)[3], float(*fLongVectorB)[3], unsigned int verts)
static float ZERO[3][3]
void SIM_mass_spring_set_motion_state(Implicit_Data *data, int index, const float x[3], const float v[3])
void SIM_mass_spring_set_position(Implicit_Data *data, int index, const float x[3])
DO_INLINE void sub_fmatrix_fmatrix(float to[3][3], const float matrixA[3][3], const float matrixB[3][3])
struct Implicit_Data Implicit_Data
DO_INLINE void add_lfvectorS_lfvectorS(float(*to)[3], float(*fLongVectorA)[3], float aS, float(*fLongVectorB)[3], float bS, unsigned int verts)
BLI_INLINE void dfdx_spring(float to[3][3], const float dir[3], float length, float L, float k)
void SIM_mass_spring_set_rest_transform(Implicit_Data *data, int index, float tfm[3][3])
void SIM_mass_spring_get_position(struct Implicit_Data *data, int index, float x[3])
static int cg_filtered(lfVector *ldV, fmatrix3x3 *lA, lfVector *lB, lfVector *z, fmatrix3x3 *S, ImplicitSolverResult *result)
DO_INLINE void cp_fmatrix(float to[3][3], const float from[3][3])
BLI_INLINE void poly_norm(lfVector *data, int i, int j, int *inds, int len, float r_dir[3])
void SIM_mass_spring_set_new_position(struct Implicit_Data *data, int index, const float x[3])
DO_INLINE void initdiag_bfmatrix(fmatrix3x3 *matrix, float m3[3][3])
void SIM_mass_spring_force_edge_wind(Implicit_Data *data, int v1, int v2, float radius1, float radius2, const float(*winvec)[3])
DO_INLINE void cp_bfmatrix(fmatrix3x3 *to, fmatrix3x3 *from)
BLI_INLINE void spring_grad_dir(Implicit_Data *data, int i, int j, float edge[3], float dir[3], float grad_dir[3][3])
void SIM_mass_spring_force_face_wind(Implicit_Data *data, int v1, int v2, int v3, const float(*winvec)[3])
static float calc_nor_area_tri(float nor[3], const float v1[3], const float v2[3], const float v3[3])
BLI_INLINE void root_to_world_v3(Implicit_Data *data, int index, float r[3], const float v[3])
BLI_INLINE void poly_avg(lfVector *data, const int *inds, int len, float r_avg[3])
void SIM_mass_spring_get_new_velocity(struct Implicit_Data *data, int index, float v[3])
DO_INLINE void subadd_bfmatrixS_bfmatrixS(fmatrix3x3 *to, fmatrix3x3 *from, float aS, fmatrix3x3 *matrix, float bS)
void SIM_mass_spring_add_constraint_ndof1(Implicit_Data *data, int index, const float c1[3], const float c2[3], const float dV[3])
BLI_INLINE void world_to_root_m3(Implicit_Data *data, int index, float r[3][3], const float m[3][3])
static float I[3][3]
DO_INLINE void sub_lfvector_lfvector(float(*to)[3], float(*fLongVectorA)[3], float(*fLongVectorB)[3], unsigned int verts)
DO_INLINE void submul_lfvectorS(float(*to)[3], float(*fLongVector)[3], float scalar, unsigned int verts)
BLI_INLINE void edge_norm(lfVector *data, int i, int j, float r_dir[3])
DO_INLINE void sub_lfvector_lfvectorS(float(*to)[3], float(*fLongVectorA)[3], float(*fLongVectorB)[3], float bS, unsigned int verts)
BLI_INLINE void spring_hairbend_estimate_dfdv(Implicit_Data *data, int i, int j, int k, const float goal[3], float stiffness, float damping, int q, float dfdv[3][3])
BLI_INLINE float fbstar(float length, float L, float kb, float cb)
DO_INLINE void del_bfmatrix(fmatrix3x3 *matrix)
void SIM_mass_spring_force_reference_frame(Implicit_Data *data, int index, const float acceleration[3], const float omega[3], const float domega_dt[3], float mass)
float SIM_tri_tetra_volume_signed_6x(Implicit_Data *data, int v1, int v2, int v3)
bool SIM_mass_spring_force_spring_angular(Implicit_Data *data, int i, int j, int *i_a, int *i_b, int len_a, int len_b, float restang, float stiffness, float damping)
BLI_INLINE void spring_hairbend_estimate_dfdx(Implicit_Data *data, int i, int j, int k, const float goal[3], float stiffness, float damping, int q, float dfdx[3][3])
BLI_INLINE void world_to_root_v3(Implicit_Data *data, int index, float r[3], const float v[3])
BLI_INLINE void dfdv_damp(float to[3][3], const float dir[3], float damping)
void SIM_mass_spring_force_extern(struct Implicit_Data *data, int i, const float f[3], float dfdx[3][3], float dfdv[3][3])
float lfVector[3]
void SIM_mass_spring_apply_result(Implicit_Data *data)
BLI_INLINE float bend_angle(const float dir_a[3], const float dir_b[3], const float dir_e[3])
BLI_INLINE void spring_angle(Implicit_Data *data, int i, int j, int *i_a, int *i_b, int len_a, int len_b, float r_dir_a[3], float r_dir_b[3], float *r_angle, float r_vel_a[3], float r_vel_b[3])
DO_INLINE void init_bfmatrix(fmatrix3x3 *matrix, float m3[3][3])
DO_INLINE void cp_lfvector(float(*to)[3], float(*from)[3], unsigned int verts)
void SIM_mass_spring_force_vertex_wind(Implicit_Data *data, int v, float UNUSED(radius), const float(*winvec)[3])
void(* MEM_freeN)(void *vmemh)
Definition: mallocn.c:27
void *(* MEM_callocN)(size_t len, const char *str)
Definition: mallocn.c:31
static float P(float k)
Definition: math_interp.c:25
#define M
#define L
#define atan2f(x, y)
Definition: metal/compat.h:227
#define sqrtf(x)
Definition: metal/compat.h:243
static unsigned c
Definition: RandGen.cpp:83
static unsigned a[3]
Definition: RandGen.cpp:78
INLINE Rall1d< T, V, S > cos(const Rall1d< T, V, S > &arg)
Definition: rall1d.h:319
INLINE Rall1d< T, V, S > sin(const Rall1d< T, V, S > &arg)
Definition: rall1d.h:311
static void area(int d1, int d2, int e1, int e2, float weights[2])
T length(const vec_base< T, Size > &a)
static const pxr::TfToken b("b", pxr::TfToken::Immortal)
static const pxr::TfToken g("g", pxr::TfToken::Immortal)
static const pxr::TfToken density("density", pxr::TfToken::Immortal)
fmatrix3x3 * M
fmatrix3x3 * bigI
fmatrix3x3 * A
fmatrix3x3 * Pinv
fmatrix3x3 * S
fmatrix3x3 * dFdX
fmatrix3x3 * dFdV
fmatrix3x3 * P
fmatrix3x3 * tfm
unsigned int r
float m[3][3]
unsigned int c
unsigned int scount
unsigned int vcount
double PIL_check_seconds_timer(void)
Definition: time.c:64
CCL_NAMESPACE_BEGIN struct Window V