Blender  V3.3
implicit_eigen.cpp
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_EIGEN
11 
12 //#define USE_EIGEN_CORE
13 # define USE_EIGEN_CONSTRAINED_CG
14 
15 # ifdef __GNUC__
16 # pragma GCC diagnostic push
17 /* XXX suppress verbose warnings in eigen */
18 //# pragma GCC diagnostic ignored "-Wlogical-op"
19 # endif
20 
21 # ifndef IMPLICIT_ENABLE_EIGEN_DEBUG
22 # ifdef NDEBUG
23 # define IMPLICIT_NDEBUG
24 # endif
25 # define NDEBUG
26 # endif
27 
28 # include <Eigen/Sparse>
29 # include <Eigen/src/Core/util/DisableStupidWarnings.h>
30 
31 # ifdef USE_EIGEN_CONSTRAINED_CG
33 # endif
34 
35 # ifndef IMPLICIT_ENABLE_EIGEN_DEBUG
36 # ifndef IMPLICIT_NDEBUG
37 # undef NDEBUG
38 # else
39 # undef IMPLICIT_NDEBUG
40 # endif
41 # endif
42 
43 # ifdef __GNUC__
44 # pragma GCC diagnostic pop
45 # endif
46 
47 # include "MEM_guardedalloc.h"
48 
49 extern "C" {
50 # include "DNA_meshdata_types.h"
51 # include "DNA_object_force_types.h"
52 # include "DNA_object_types.h"
53 # include "DNA_scene_types.h"
54 # include "DNA_texture_types.h"
55 
56 # include "BLI_linklist.h"
57 # include "BLI_math.h"
58 # include "BLI_utildefines.h"
59 
60 # include "BKE_cloth.h"
61 # include "BKE_collision.h"
62 # include "BKE_effect.h"
63 # include "BKE_global.h"
64 
65 # include "SIM_mass_spring.h"
66 }
67 
68 typedef float Scalar;
69 
70 static float I[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
71 
72 /* slightly extended Eigen vector class
73  * with conversion to/from plain C float array
74  */
75 class fVector : public Eigen::Vector3f {
76  public:
77  typedef float *ctype;
78 
79  fVector()
80  {
81  }
82 
83  fVector(const ctype &v)
84  {
85  for (int k = 0; k < 3; k++) {
86  coeffRef(k) = v[k];
87  }
88  }
89 
90  fVector &operator=(const ctype &v)
91  {
92  for (int k = 0; k < 3; k++) {
93  coeffRef(k) = v[k];
94  }
95  return *this;
96  }
97 
98  operator ctype()
99  {
100  return data();
101  }
102 };
103 
104 /* slightly extended Eigen matrix class
105  * with conversion to/from plain C float array
106  */
107 class fMatrix : public Eigen::Matrix3f {
108  public:
109  typedef float (*ctype)[3];
110 
111  fMatrix()
112  {
113  }
114 
115  fMatrix(const ctype &v)
116  {
117  for (int k = 0; k < 3; k++) {
118  for (int l = 0; l < 3; l++) {
119  coeffRef(l, k) = v[k][l];
120  }
121  }
122  }
123 
124  fMatrix &operator=(const ctype &v)
125  {
126  for (int k = 0; k < 3; k++) {
127  for (int l = 0; l < 3; l++) {
128  coeffRef(l, k) = v[k][l];
129  }
130  }
131  return *this;
132  }
133 
134  operator ctype()
135  {
136  return (ctype)data();
137  }
138 };
139 
140 /* Extension of dense Eigen vectors,
141  * providing 3-float block access for blenlib math functions
142  */
143 class lVector : public Eigen::VectorXf {
144  public:
145  typedef Eigen::VectorXf base_t;
146 
147  lVector()
148  {
149  }
150 
151  template<typename T> lVector &operator=(T rhs)
152  {
153  base_t::operator=(rhs);
154  return *this;
155  }
156 
157  float *v3(int vertex)
158  {
159  return &coeffRef(3 * vertex);
160  }
161 
162  const float *v3(int vertex) const
163  {
164  return &coeffRef(3 * vertex);
165  }
166 };
167 
168 typedef Eigen::Triplet<Scalar> Triplet;
169 typedef std::vector<Triplet> TripletList;
170 
171 typedef Eigen::SparseMatrix<Scalar> lMatrix;
172 
173 /* Constructor type that provides more convenient handling of Eigen triplets
174  * for efficient construction of sparse 3x3 block matrices.
175  * This should be used for building lMatrix instead of writing to such lMatrix directly (which is
176  * very inefficient). After all elements have been defined using the set() method, the actual
177  * matrix can be filled using construct().
178  */
179 struct lMatrixCtor {
180  lMatrixCtor()
181  {
182  }
183 
184  void reset()
185  {
186  m_trips.clear();
187  }
188 
189  void reserve(int numverts)
190  {
191  /* reserve for diagonal entries */
192  m_trips.reserve(numverts * 9);
193  }
194 
195  void add(int i, int j, const fMatrix &m)
196  {
197  i *= 3;
198  j *= 3;
199  for (int k = 0; k < 3; k++) {
200  for (int l = 0; l < 3; l++) {
201  m_trips.push_back(Triplet(i + k, j + l, m.coeff(l, k)));
202  }
203  }
204  }
205 
206  void sub(int i, int j, const fMatrix &m)
207  {
208  i *= 3;
209  j *= 3;
210  for (int k = 0; k < 3; k++) {
211  for (int l = 0; l < 3; l++) {
212  m_trips.push_back(Triplet(i + k, j + l, -m.coeff(l, k)));
213  }
214  }
215  }
216 
217  inline void construct(lMatrix &m)
218  {
219  m.setFromTriplets(m_trips.begin(), m_trips.end());
220  m_trips.clear();
221  }
222 
223  private:
224  TripletList m_trips;
225 };
226 
227 # ifdef USE_EIGEN_CORE
228 typedef Eigen::ConjugateGradient<lMatrix, Eigen::Lower, Eigen::DiagonalPreconditioner<Scalar>>
230 # endif
231 # ifdef USE_EIGEN_CONSTRAINED_CG
233  Eigen::Lower,
234  lMatrix,
235  Eigen::DiagonalPreconditioner<Scalar>>
236  ConstraintConjGrad;
237 # endif
238 using Eigen::ComputationInfo;
239 
240 static void print_lvector(const lVector &v)
241 {
242  for (int i = 0; i < v.rows(); i++) {
243  if (i > 0 && i % 3 == 0) {
244  printf("\n");
245  }
246 
247  printf("%f,\n", v[i]);
248  }
249 }
250 
251 static void print_lmatrix(const lMatrix &m)
252 {
253  for (int j = 0; j < m.rows(); j++) {
254  if (j > 0 && j % 3 == 0) {
255  printf("\n");
256  }
257 
258  for (int i = 0; i < m.cols(); i++) {
259  if (i > 0 && i % 3 == 0) {
260  printf(" ");
261  }
262 
263  implicit_print_matrix_elem(m.coeff(j, i));
264  }
265  printf("\n");
266  }
267 }
268 
269 BLI_INLINE void lMatrix_reserve_elems(lMatrix &m, int num)
270 {
271  m.reserve(Eigen::VectorXi::Constant(m.cols(), num));
272 }
273 
274 BLI_INLINE float *lVector_v3(lVector &v, int vertex)
275 {
276  return v.data() + 3 * vertex;
277 }
278 
279 BLI_INLINE const float *lVector_v3(const lVector &v, int vertex)
280 {
281  return v.data() + 3 * vertex;
282 }
283 
284 # if 0
285 BLI_INLINE void triplets_m3(TripletList &tlist, float m[3][3], int i, int j)
286 {
287  i *= 3;
288  j *= 3;
289  for (int l = 0; l < 3; l++) {
290  for (int k = 0; k < 3; k++) {
291  tlist.push_back(Triplet(i + k, j + l, m[k][l]));
292  }
293  }
294 }
295 
296 BLI_INLINE void triplets_m3fl(TripletList &tlist, float m[3][3], int i, int j, float factor)
297 {
298  i *= 3;
299  j *= 3;
300  for (int l = 0; l < 3; l++) {
301  for (int k = 0; k < 3; k++) {
302  tlist.push_back(Triplet(i + k, j + l, m[k][l] * factor));
303  }
304  }
305 }
306 
307 BLI_INLINE void lMatrix_add_triplets(lMatrix &r, const TripletList &tlist)
308 {
309  lMatrix t(r.rows(), r.cols());
310  t.setFromTriplets(tlist.begin(), tlist.end());
311  r += t;
312 }
313 
314 BLI_INLINE void lMatrix_madd_triplets(lMatrix &r, const TripletList &tlist, float f)
315 {
316  lMatrix t(r.rows(), r.cols());
317  t.setFromTriplets(tlist.begin(), tlist.end());
318  r += f * t;
319 }
320 
321 BLI_INLINE void lMatrix_sub_triplets(lMatrix &r, const TripletList &tlist)
322 {
323  lMatrix t(r.rows(), r.cols());
324  t.setFromTriplets(tlist.begin(), tlist.end());
325  r -= t;
326 }
327 # endif
328 
329 BLI_INLINE void outerproduct(float r[3][3], const float a[3], const float b[3])
330 {
331  mul_v3_v3fl(r[0], a, b[0]);
332  mul_v3_v3fl(r[1], a, b[1]);
333  mul_v3_v3fl(r[2], a, b[2]);
334 }
335 
336 BLI_INLINE void cross_m3_v3m3(float r[3][3], const float v[3], float m[3][3])
337 {
338  cross_v3_v3v3(r[0], v, m[0]);
339  cross_v3_v3v3(r[1], v, m[1]);
340  cross_v3_v3v3(r[2], v, m[2]);
341 }
342 
343 BLI_INLINE void cross_v3_identity(float r[3][3], const float v[3])
344 {
345  r[0][0] = 0.0f;
346  r[1][0] = v[2];
347  r[2][0] = -v[1];
348  r[0][1] = -v[2];
349  r[1][1] = 0.0f;
350  r[2][1] = v[0];
351  r[0][2] = v[1];
352  r[1][2] = -v[0];
353  r[2][2] = 0.0f;
354 }
355 
356 BLI_INLINE void madd_m3_m3fl(float r[3][3], float m[3][3], float f)
357 {
358  r[0][0] += m[0][0] * f;
359  r[0][1] += m[0][1] * f;
360  r[0][2] += m[0][2] * f;
361  r[1][0] += m[1][0] * f;
362  r[1][1] += m[1][1] * f;
363  r[1][2] += m[1][2] * f;
364  r[2][0] += m[2][0] * f;
365  r[2][1] += m[2][1] * f;
366  r[2][2] += m[2][2] * f;
367 }
368 
369 BLI_INLINE void madd_m3_m3m3fl(float r[3][3], float a[3][3], float b[3][3], float f)
370 {
371  r[0][0] = a[0][0] + b[0][0] * f;
372  r[0][1] = a[0][1] + b[0][1] * f;
373  r[0][2] = a[0][2] + b[0][2] * f;
374  r[1][0] = a[1][0] + b[1][0] * f;
375  r[1][1] = a[1][1] + b[1][1] * f;
376  r[1][2] = a[1][2] + b[1][2] * f;
377  r[2][0] = a[2][0] + b[2][0] * f;
378  r[2][1] = a[2][1] + b[2][1] * f;
379  r[2][2] = a[2][2] + b[2][2] * f;
380 }
381 
382 struct Implicit_Data {
383  typedef std::vector<fMatrix> fMatrixVector;
384 
385  Implicit_Data(int numverts)
386  {
387  resize(numverts);
388  }
389 
390  void resize(int numverts)
391  {
392  this->numverts = numverts;
393  int tot = 3 * numverts;
394 
395  M.resize(tot, tot);
396  F.resize(tot);
397  dFdX.resize(tot, tot);
398  dFdV.resize(tot, tot);
399 
400  tfm.resize(numverts, I);
401 
402  X.resize(tot);
403  Xnew.resize(tot);
404  V.resize(tot);
405  Vnew.resize(tot);
406 
407  A.resize(tot, tot);
408  B.resize(tot);
409 
410  dV.resize(tot);
411  z.resize(tot);
412  S.resize(tot, tot);
413 
414  iM.reserve(numverts);
415  idFdX.reserve(numverts);
416  idFdV.reserve(numverts);
417  iS.reserve(numverts);
418  }
419 
420  int numverts;
421 
422  /* inputs */
423  lMatrix M; /* masses */
424  lVector F; /* forces */
425  lMatrix dFdX, dFdV; /* force jacobians */
426 
427  fMatrixVector tfm; /* local coordinate transform */
428 
429  /* motion state data */
430  lVector X, Xnew; /* positions */
431  lVector V, Vnew; /* velocities */
432 
433  /* internal solver data */
434  lVector B; /* B for A*dV = B */
435  lMatrix A; /* A for A*dV = B */
436 
437  lVector dV; /* velocity change (solution of A*dV = B) */
438  lVector z; /* target velocity in constrained directions */
439  lMatrix S; /* filtering matrix for constraints */
440 
441  /* temporary constructors */
442  lMatrixCtor iM; /* masses */
443  lMatrixCtor idFdX, idFdV; /* force jacobians */
444  lMatrixCtor iS; /* filtering matrix for constraints */
445 };
446 
447 Implicit_Data *SIM_mass_spring_solver_create(int numverts, int numsprings)
448 {
449  Implicit_Data *id = new Implicit_Data(numverts);
450  return id;
451 }
452 
454 {
455  if (id) {
456  delete id;
457  }
458 }
459 
461 {
462  if (id) {
463  return id->numverts;
464  }
465  else {
466  return 0;
467  }
468 }
469 
470 /* ==== Transformation from/to root reference frames ==== */
471 
472 BLI_INLINE void world_to_root_v3(Implicit_Data *data, int index, float r[3], const float v[3])
473 {
474  copy_v3_v3(r, v);
475  mul_transposed_m3_v3(data->tfm[index], r);
476 }
477 
478 BLI_INLINE void root_to_world_v3(Implicit_Data *data, int index, float r[3], const float v[3])
479 {
480  mul_v3_m3v3(r, data->tfm[index], v);
481 }
482 
483 BLI_INLINE void world_to_root_m3(Implicit_Data *data, int index, float r[3][3], float m[3][3])
484 {
485  float trot[3][3];
486  copy_m3_m3(trot, data->tfm[index]);
487  transpose_m3(trot);
488  mul_m3_m3m3(r, trot, m);
489 }
490 
491 BLI_INLINE void root_to_world_m3(Implicit_Data *data, int index, float r[3][3], float m[3][3])
492 {
493  mul_m3_m3m3(r, data->tfm[index], m);
494 }
495 
496 /* ================================ */
497 
499 {
500 # ifdef USE_EIGEN_CORE
501  typedef ConjugateGradient solver_t;
502 # endif
503 # ifdef USE_EIGEN_CONSTRAINED_CG
504  typedef ConstraintConjGrad solver_t;
505 # endif
506 
507  data->iM.construct(data->M);
508  data->idFdX.construct(data->dFdX);
509  data->idFdV.construct(data->dFdV);
510  data->iS.construct(data->S);
511 
512  solver_t cg;
513  cg.setMaxIterations(100);
514  cg.setTolerance(0.01f);
515 
516 # ifdef USE_EIGEN_CONSTRAINED_CG
517  cg.filter() = data->S;
518 # endif
519 
520  data->A = data->M - dt * data->dFdV - dt * dt * data->dFdX;
521  cg.compute(data->A);
522 
523  data->B = dt * data->F + dt * dt * data->dFdX * data->V;
524 
525 # ifdef IMPLICIT_PRINT_SOLVER_INPUT_OUTPUT
526  printf("==== A ====\n");
527  print_lmatrix(id->A);
528  printf("==== z ====\n");
529  print_lvector(id->z);
530  printf("==== B ====\n");
531  print_lvector(id->B);
532  printf("==== S ====\n");
533  print_lmatrix(id->S);
534 # endif
535 
536 # ifdef USE_EIGEN_CORE
537  data->dV = cg.solve(data->B);
538 # endif
539 # ifdef USE_EIGEN_CONSTRAINED_CG
540  data->dV = cg.solveWithGuess(data->B, data->z);
541 # endif
542 
543 # ifdef IMPLICIT_PRINT_SOLVER_INPUT_OUTPUT
544  printf("==== dV ====\n");
545  print_lvector(id->dV);
546  printf("========\n");
547 # endif
548 
549  data->Vnew = data->V + data->dV;
550 
551  switch (cg.info()) {
552  case Eigen::Success:
553  result->status = SIM_SOLVER_SUCCESS;
554  break;
555  case Eigen::NoConvergence:
557  break;
558  case Eigen::InvalidInput:
560  break;
561  case Eigen::NumericalIssue:
563  break;
564  }
565 
566  result->iterations = cg.iterations();
567  result->error = cg.error();
568 
569  return cg.info() == Eigen::Success;
570 }
571 
573 {
574  data->Xnew = data->X + data->Vnew * dt;
575  return true;
576 }
577 
578 /* ================================ */
579 
581 {
582  data->X = data->Xnew;
583  data->V = data->Vnew;
584 }
585 
586 void SIM_mass_spring_set_vertex_mass(Implicit_Data *data, int index, float mass)
587 {
588  float m[3][3];
589  copy_m3_m3(m, I);
590  mul_m3_fl(m, mass);
591  data->iM.add(index, index, m);
592 }
593 
594 void SIM_mass_spring_set_rest_transform(Implicit_Data *data, int index, float tfm[3][3])
595 {
596 # ifdef CLOTH_ROOT_FRAME
597  copy_m3_m3(data->tfm[index], tfm);
598 # else
599  unit_m3(data->tfm[index]);
600  (void)tfm;
601 # endif
602 }
603 
605  int index,
606  const float x[3],
607  const float v[3])
608 {
609  world_to_root_v3(data, index, data->X.v3(index), x);
610  world_to_root_v3(data, index, data->V.v3(index), v);
611 }
612 
613 void SIM_mass_spring_set_position(Implicit_Data *data, int index, const float x[3])
614 {
615  world_to_root_v3(data, index, data->X.v3(index), x);
616 }
617 
618 void SIM_mass_spring_set_velocity(Implicit_Data *data, int index, const float v[3])
619 {
620  world_to_root_v3(data, index, data->V.v3(index), v);
621 }
622 
624  int index,
625  float x[3],
626  float v[3])
627 {
628  if (x) {
629  root_to_world_v3(data, index, x, data->X.v3(index));
630  }
631  if (v) {
632  root_to_world_v3(data, index, v, data->V.v3(index));
633  }
634 }
635 
636 void SIM_mass_spring_get_position(struct Implicit_Data *data, int index, float x[3])
637 {
638  root_to_world_v3(data, index, x, data->X.v3(index));
639 }
640 
641 void SIM_mass_spring_get_new_velocity(Implicit_Data *data, int index, float v[3])
642 {
643  root_to_world_v3(data, index, v, data->V.v3(index));
644 }
645 
646 void SIM_mass_spring_set_new_velocity(Implicit_Data *data, int index, const float v[3])
647 {
648  world_to_root_v3(data, index, data->V.v3(index), v);
649 }
650 
652 {
653  int numverts = data->numverts;
654  for (int i = 0; i < numverts; i++) {
655  data->iS.add(i, i, I);
656  zero_v3(data->z.v3(i));
657  }
658 }
659 
660 void SIM_mass_spring_add_constraint_ndof0(Implicit_Data *data, int index, const float dV[3])
661 {
662  data->iS.sub(index, index, I);
663 
664  world_to_root_v3(data, index, data->z.v3(index), dV);
665 }
666 
668  Implicit_Data *data, int index, const float c1[3], const float c2[3], const float dV[3])
669 {
670  float m[3][3], p[3], q[3], u[3], cmat[3][3];
671 
672  world_to_root_v3(data, index, p, c1);
673  outerproduct(cmat, p, p);
674  copy_m3_m3(m, cmat);
675 
676  world_to_root_v3(data, index, q, c2);
677  outerproduct(cmat, q, q);
678  add_m3_m3m3(m, m, cmat);
679 
680  /* XXX not sure but multiplication should work here */
681  data->iS.sub(index, index, m);
682  // mul_m3_m3m3(data->S[index].m, data->S[index].m, m);
683 
684  world_to_root_v3(data, index, u, dV);
685  add_v3_v3(data->z.v3(index), u);
686 }
687 
689  int index,
690  const float c1[3],
691  const float dV[3])
692 {
693  float m[3][3], p[3], u[3], cmat[3][3];
694 
695  world_to_root_v3(data, index, p, c1);
696  outerproduct(cmat, p, p);
697  copy_m3_m3(m, cmat);
698 
699  data->iS.sub(index, index, m);
700  // mul_m3_m3m3(data->S[index].m, data->S[index].m, m);
701 
702  world_to_root_v3(data, index, u, dV);
703  add_v3_v3(data->z.v3(index), u);
704 }
705 
707 {
708  data->F.setZero();
709  data->dFdX.setZero();
710  data->dFdV.setZero();
711 }
712 
714  int index,
715  const float acceleration[3],
716  const float omega[3],
717  const float domega_dt[3],
718  float mass)
719 {
720 # ifdef CLOTH_ROOT_FRAME
721  float acc[3], w[3], dwdt[3];
722  float f[3], dfdx[3][3], dfdv[3][3];
723  float euler[3], coriolis[3], centrifugal[3], rotvel[3];
724  float deuler[3][3], dcoriolis[3][3], dcentrifugal[3][3], drotvel[3][3];
725 
726  world_to_root_v3(data, index, acc, acceleration);
727  world_to_root_v3(data, index, w, omega);
728  world_to_root_v3(data, index, dwdt, domega_dt);
729 
730  cross_v3_v3v3(euler, dwdt, data->X.v3(index));
731  cross_v3_v3v3(coriolis, w, data->V.v3(index));
732  mul_v3_fl(coriolis, 2.0f);
733  cross_v3_v3v3(rotvel, w, data->X.v3(index));
734  cross_v3_v3v3(centrifugal, w, rotvel);
735 
736  sub_v3_v3v3(f, acc, euler);
737  sub_v3_v3(f, coriolis);
738  sub_v3_v3(f, centrifugal);
739 
740  mul_v3_fl(f, mass); /* F = m * a */
741 
742  cross_v3_identity(deuler, dwdt);
743  cross_v3_identity(dcoriolis, w);
744  mul_m3_fl(dcoriolis, 2.0f);
745  cross_v3_identity(drotvel, w);
746  cross_m3_v3m3(dcentrifugal, w, drotvel);
747 
748  add_m3_m3m3(dfdx, deuler, dcentrifugal);
749  negate_m3(dfdx);
750  mul_m3_fl(dfdx, mass);
751 
752  copy_m3_m3(dfdv, dcoriolis);
753  negate_m3(dfdv);
754  mul_m3_fl(dfdv, mass);
755 
756  add_v3_v3(data->F.v3(index), f);
757  data->idFdX.add(index, index, dfdx);
758  data->idFdV.add(index, index, dfdv);
759 # else
760  (void)data;
761  (void)index;
762  (void)acceleration;
763  (void)omega;
764  (void)domega_dt;
765 # endif
766 }
767 
768 void SIM_mass_spring_force_gravity(Implicit_Data *data, int index, float mass, const float g[3])
769 {
770  /* force = mass * acceleration (in this case: gravity) */
771  float f[3];
772  world_to_root_v3(data, index, f, g);
773  mul_v3_fl(f, mass);
774 
775  add_v3_v3(data->F.v3(index), f);
776 }
777 
779 {
780  int numverts = data->numverts;
781  for (int i = 0; i < numverts; i++) {
782  float tmp[3][3];
783 
784  /* NOTE: Uses root space velocity, no need to transform. */
785  madd_v3_v3fl(data->F.v3(i), data->V.v3(i), -drag);
786 
787  copy_m3_m3(tmp, I);
788  mul_m3_fl(tmp, -drag);
789  data->idFdV.add(i, i, tmp);
790  }
791 }
792 
794  struct Implicit_Data *data, int i, const float f[3], float dfdx[3][3], float dfdv[3][3])
795 {
796  float tf[3], tdfdx[3][3], tdfdv[3][3];
797  world_to_root_v3(data, i, tf, f);
798  world_to_root_m3(data, i, tdfdx, dfdx);
799  world_to_root_m3(data, i, tdfdv, dfdv);
800 
801  add_v3_v3(data->F.v3(i), tf);
802  data->idFdX.add(i, i, tdfdx);
803  data->idFdV.add(i, i, tdfdv);
804 }
805 
806 static float calc_nor_area_tri(float nor[3],
807  const float v1[3],
808  const float v2[3],
809  const float v3[3])
810 {
811  float n1[3], n2[3];
812 
813  sub_v3_v3v3(n1, v1, v2);
814  sub_v3_v3v3(n2, v2, v3);
815 
816  cross_v3_v3v3(nor, n1, n2);
817  return normalize_v3(nor) / 2.0f;
818 }
819 
820 /* XXX does not support force jacobians yet,
821  * since the effector system does not provide them either. */
823  Implicit_Data *data, int v1, int v2, int v3, const float (*winvec)[3])
824 {
825  const float effector_scale = 0.02f;
826  float win[3], nor[3], area;
827  float factor;
828 
829  /* calculate face normal and area */
830  area = calc_nor_area_tri(nor, data->X.v3(v1), data->X.v3(v2), data->X.v3(v3));
831  factor = effector_scale * area / 3.0f;
832 
833  world_to_root_v3(data, v1, win, winvec[v1]);
834  madd_v3_v3fl(data->F.v3(v1), nor, factor * dot_v3v3(win, nor));
835 
836  world_to_root_v3(data, v2, win, winvec[v2]);
837  madd_v3_v3fl(data->F.v3(v2), nor, factor * dot_v3v3(win, nor));
838 
839  world_to_root_v3(data, v3, win, winvec[v3]);
840  madd_v3_v3fl(data->F.v3(v3), nor, factor * dot_v3v3(win, nor));
841 }
842 
843 void SIM_mass_spring_force_edge_wind(Implicit_Data *data, int v1, int v2, const float (*winvec)[3])
844 {
845  const float effector_scale = 0.01;
846  float win[3], dir[3], nor[3], length;
847 
848  sub_v3_v3v3(dir, data->X.v3(v1), data->X.v3(v2));
849  length = normalize_v3(dir);
850 
851  world_to_root_v3(data, v1, win, winvec[v1]);
852  madd_v3_v3v3fl(nor, win, dir, -dot_v3v3(win, dir));
853  madd_v3_v3fl(data->F.v3(v1), nor, effector_scale * length);
854 
855  world_to_root_v3(data, v2, win, winvec[v2]);
856  madd_v3_v3v3fl(nor, win, dir, -dot_v3v3(win, dir));
857  madd_v3_v3fl(data->F.v3(v2), nor, effector_scale * length);
858 }
859 
860 BLI_INLINE void dfdx_spring(float to[3][3], const float dir[3], float length, float L, float k)
861 {
862  /* dir is unit length direction, rest is spring's restlength, k is spring constant. */
863  // return ((I - outerprod(dir, dir)) * Min(1.0f, rest / length) - I) * -k;
864  outerproduct(to, dir, dir);
865  sub_m3_m3m3(to, I, to);
866 
867  mul_m3_fl(to, (L / length));
868  sub_m3_m3m3(to, to, I);
869  mul_m3_fl(to, k);
870 }
871 
872 /* unused */
873 # if 0
874 BLI_INLINE void dfdx_damp(float to[3][3],
875  const float dir[3],
876  float length,
877  const float vel[3],
878  float rest,
879  float damping)
880 {
881  /* Inner spring damping vel is the relative velocity of the endpoints. */
882  // return (I - outerprod(dir, dir)) * (-damping * -(dot(dir, vel) / Max(length, rest)));
883  mul_fvectorT_fvector(to, dir, dir);
884  sub_fmatrix_fmatrix(to, I, to);
885  mul_fmatrix_S(to, (-damping * -(dot_v3v3(dir, vel) / MAX2(length, rest))));
886 }
887 # endif
888 
889 BLI_INLINE void dfdv_damp(float to[3][3], const float dir[3], float damping)
890 {
891  /* Derivative of force with regards to velocity. */
892  outerproduct(to, dir, dir);
893  mul_m3_fl(to, -damping);
894 }
895 
896 BLI_INLINE float fb(float length, float L)
897 {
898  float x = length / L;
899  return (-11.541f * powf(x, 4) + 34.193f * powf(x, 3) - 39.083f * powf(x, 2) + 23.116f * x -
900  9.713f);
901 }
902 
903 BLI_INLINE float fbderiv(float length, float L)
904 {
905  float x = length / L;
906 
907  return (-46.164f * powf(x, 3) + 102.579f * powf(x, 2) - 78.166f * x + 23.116f);
908 }
909 
910 BLI_INLINE float fbstar(float length, float L, float kb, float cb)
911 {
912  float tempfb_fl = kb * fb(length, L);
913  float fbstar_fl = cb * (length - L);
914 
915  if (tempfb_fl < fbstar_fl) {
916  return fbstar_fl;
917  }
918  else {
919  return tempfb_fl;
920  }
921 }
922 
923 /* Function to calculate bending spring force (taken from Choi & Co). */
924 BLI_INLINE float fbstar_jacobi(float length, float L, float kb, float cb)
925 {
926  float tempfb_fl = kb * fb(length, L);
927  float fbstar_fl = cb * (length - L);
928 
929  if (tempfb_fl < fbstar_fl) {
930  return -cb;
931  }
932  else {
933  return -kb * fbderiv(length, L);
934  }
935 }
936 
937 /* calculate elongation */
939  int i,
940  int j,
941  float r_extent[3],
942  float r_dir[3],
943  float *r_length,
944  float r_vel[3])
945 {
946  sub_v3_v3v3(r_extent, data->X.v3(j), data->X.v3(i));
947  sub_v3_v3v3(r_vel, data->V.v3(j), data->V.v3(i));
948  *r_length = len_v3(r_extent);
949 
950  if (*r_length > ALMOST_ZERO) {
951 # if 0
952  if (length > L) {
953  if ((clmd->sim_parms->flags & CSIMSETT_FLAG_TEARING_ENABLED) &&
954  (((length - L) * 100.0f / L) > clmd->sim_parms->maxspringlen)) {
955  /* cut spring! */
956  s->flags |= CSPRING_FLAG_DEACTIVATE;
957  return false;
958  }
959  }
960 # endif
961  mul_v3_v3fl(r_dir, r_extent, 1.0f / (*r_length));
962  }
963  else {
964  zero_v3(r_dir);
965  }
966 
967  return true;
968 }
969 
971  Implicit_Data *data, int i, int j, const float f[3], float dfdx[3][3], float dfdv[3][3])
972 {
973  add_v3_v3(data->F.v3(i), f);
974  sub_v3_v3(data->F.v3(j), f);
975 
976  data->idFdX.add(i, i, dfdx);
977  data->idFdX.add(j, j, dfdx);
978  data->idFdX.sub(i, j, dfdx);
979  data->idFdX.sub(j, i, dfdx);
980 
981  data->idFdV.add(i, i, dfdv);
982  data->idFdV.add(j, j, dfdv);
983  data->idFdV.sub(i, j, dfdv);
984  data->idFdV.sub(j, i, dfdv);
985 }
986 
988  int i,
989  int j,
990  float restlen,
991  float stiffness,
992  float damping,
993  bool no_compress,
994  float clamp_force,
995  float r_f[3],
996  float r_dfdx[3][3],
997  float r_dfdv[3][3])
998 {
999  float extent[3], length, dir[3], vel[3];
1000 
1001  /* calculate elongation */
1002  spring_length(data, i, j, extent, dir, &length, vel);
1003 
1004  if (length > restlen || no_compress) {
1005  float stretch_force, f[3], dfdx[3][3], dfdv[3][3];
1006 
1007  stretch_force = stiffness * (length - restlen);
1008  if (clamp_force > 0.0f && stretch_force > clamp_force) {
1009  stretch_force = clamp_force;
1010  }
1011  mul_v3_v3fl(f, dir, stretch_force);
1012 
1013  /* Ascher & Boxman, p.21: Damping only during elongation
1014  * something wrong with it... */
1015  madd_v3_v3fl(f, dir, damping * dot_v3v3(vel, dir));
1016 
1017  dfdx_spring(dfdx, dir, length, restlen, stiffness);
1018  dfdv_damp(dfdv, dir, damping);
1019 
1020  apply_spring(data, i, j, f, dfdx, dfdv);
1021 
1022  if (r_f) {
1023  copy_v3_v3(r_f, f);
1024  }
1025  if (r_dfdx) {
1026  copy_m3_m3(r_dfdx, dfdx);
1027  }
1028  if (r_dfdv) {
1029  copy_m3_m3(r_dfdv, dfdv);
1030  }
1031 
1032  return true;
1033  }
1034  else {
1035  if (r_f) {
1036  zero_v3(r_f);
1037  }
1038  if (r_dfdx) {
1039  zero_m3(r_dfdx);
1040  }
1041  if (r_dfdv) {
1042  zero_m3(r_dfdv);
1043  }
1044 
1045  return false;
1046  }
1047 }
1048 
1049 /* See "Stable but Responsive Cloth" (Choi, Ko 2005) */
1051  int i,
1052  int j,
1053  float restlen,
1054  float kb,
1055  float cb,
1056  float r_f[3],
1057  float r_dfdx[3][3],
1058  float r_dfdv[3][3])
1059 {
1060  float extent[3], length, dir[3], vel[3];
1061 
1062  /* calculate elongation */
1063  spring_length(data, i, j, extent, dir, &length, vel);
1064 
1065  if (length < restlen) {
1066  float f[3], dfdx[3][3], dfdv[3][3];
1067 
1068  mul_v3_v3fl(f, dir, fbstar(length, restlen, kb, cb));
1069 
1070  outerproduct(dfdx, dir, dir);
1071  mul_m3_fl(dfdx, fbstar_jacobi(length, restlen, kb, cb));
1072 
1073  /* XXX damping not supported */
1074  zero_m3(dfdv);
1075 
1076  apply_spring(data, i, j, f, dfdx, dfdv);
1077 
1078  if (r_f) {
1079  copy_v3_v3(r_f, f);
1080  }
1081  if (r_dfdx) {
1082  copy_m3_m3(r_dfdx, dfdx);
1083  }
1084  if (r_dfdv) {
1085  copy_m3_m3(r_dfdv, dfdv);
1086  }
1087 
1088  return true;
1089  }
1090  else {
1091  if (r_f) {
1092  zero_v3(r_f);
1093  }
1094  if (r_dfdx) {
1095  zero_m3(r_dfdx);
1096  }
1097  if (r_dfdv) {
1098  zero_m3(r_dfdv);
1099  }
1100 
1101  return false;
1102  }
1103 }
1104 
1105 /* Jacobian of a direction vector.
1106  * Basically the part of the differential orthogonal to the direction,
1107  * inversely proportional to the length of the edge.
1108  *
1109  * dD_ij/dx_i = -dD_ij/dx_j = (D_ij * D_ij^T - I) / len_ij
1110  */
1112  Implicit_Data *data, int i, int j, float edge[3], float dir[3], float grad_dir[3][3])
1113 {
1114  float length;
1115 
1116  sub_v3_v3v3(edge, data->X.v3(j), data->X.v3(i));
1117  length = normalize_v3_v3(dir, edge);
1118 
1119  if (length > ALMOST_ZERO) {
1120  outerproduct(grad_dir, dir, dir);
1121  sub_m3_m3m3(grad_dir, I, grad_dir);
1122  mul_m3_fl(grad_dir, 1.0f / length);
1123  }
1124  else {
1125  zero_m3(grad_dir);
1126  }
1127 }
1128 
1129 BLI_INLINE void spring_angbend_forces(Implicit_Data *data,
1130  int i,
1131  int j,
1132  int k,
1133  const float goal[3],
1134  float stiffness,
1135  float damping,
1136  int q,
1137  const float dx[3],
1138  const float dv[3],
1139  float r_f[3])
1140 {
1141  float edge_ij[3], dir_ij[3];
1142  float edge_jk[3], dir_jk[3];
1143  float vel_ij[3], vel_jk[3], vel_ortho[3];
1144  float f_bend[3], f_damp[3];
1145  float fk[3];
1146  float dist[3];
1147 
1148  zero_v3(fk);
1149 
1150  sub_v3_v3v3(edge_ij, data->X.v3(j), data->X.v3(i));
1151  if (q == i) {
1152  sub_v3_v3(edge_ij, dx);
1153  }
1154  if (q == j) {
1155  add_v3_v3(edge_ij, dx);
1156  }
1157  normalize_v3_v3(dir_ij, edge_ij);
1158 
1159  sub_v3_v3v3(edge_jk, data->X.v3(k), data->X.v3(j));
1160  if (q == j) {
1161  sub_v3_v3(edge_jk, dx);
1162  }
1163  if (q == k) {
1164  add_v3_v3(edge_jk, dx);
1165  }
1166  normalize_v3_v3(dir_jk, edge_jk);
1167 
1168  sub_v3_v3v3(vel_ij, data->V.v3(j), data->V.v3(i));
1169  if (q == i) {
1170  sub_v3_v3(vel_ij, dv);
1171  }
1172  if (q == j) {
1173  add_v3_v3(vel_ij, dv);
1174  }
1175 
1176  sub_v3_v3v3(vel_jk, data->V.v3(k), data->V.v3(j));
1177  if (q == j) {
1178  sub_v3_v3(vel_jk, dv);
1179  }
1180  if (q == k) {
1181  add_v3_v3(vel_jk, dv);
1182  }
1183 
1184  /* bending force */
1185  sub_v3_v3v3(dist, goal, edge_jk);
1186  mul_v3_v3fl(f_bend, dist, stiffness);
1187 
1188  add_v3_v3(fk, f_bend);
1189 
1190  /* damping force */
1191  madd_v3_v3v3fl(vel_ortho, vel_jk, dir_jk, -dot_v3v3(vel_jk, dir_jk));
1192  mul_v3_v3fl(f_damp, vel_ortho, damping);
1193 
1194  sub_v3_v3(fk, f_damp);
1195 
1196  copy_v3_v3(r_f, fk);
1197 }
1198 
1199 /* Finite Differences method for estimating the jacobian of the force */
1200 BLI_INLINE void spring_angbend_estimate_dfdx(Implicit_Data *data,
1201  int i,
1202  int j,
1203  int k,
1204  const float goal[3],
1205  float stiffness,
1206  float damping,
1207  int q,
1208  float dfdx[3][3])
1209 {
1210  const float delta = 0.00001f; /* TODO: find a good heuristic for this. */
1211  float dvec_null[3][3], dvec_pos[3][3], dvec_neg[3][3];
1212  float f[3];
1213  int a, b;
1214 
1215  zero_m3(dvec_null);
1216  unit_m3(dvec_pos);
1217  mul_m3_fl(dvec_pos, delta * 0.5f);
1218  copy_m3_m3(dvec_neg, dvec_pos);
1219  negate_m3(dvec_neg);
1220 
1221  /* XXX TODO: offset targets to account for position dependency. */
1222 
1223  for (a = 0; a < 3; a++) {
1224  spring_angbend_forces(
1225  data, i, j, k, goal, stiffness, damping, q, dvec_pos[a], dvec_null[a], f);
1226  copy_v3_v3(dfdx[a], f);
1227 
1228  spring_angbend_forces(
1229  data, i, j, k, goal, stiffness, damping, q, dvec_neg[a], dvec_null[a], f);
1230  sub_v3_v3(dfdx[a], f);
1231 
1232  for (b = 0; b < 3; b++) {
1233  dfdx[a][b] /= delta;
1234  }
1235  }
1236 }
1237 
1238 /* Finite Differences method for estimating the jacobian of the force */
1239 BLI_INLINE void spring_angbend_estimate_dfdv(Implicit_Data *data,
1240  int i,
1241  int j,
1242  int k,
1243  const float goal[3],
1244  float stiffness,
1245  float damping,
1246  int q,
1247  float dfdv[3][3])
1248 {
1249  const float delta = 0.00001f; /* TODO: find a good heuristic for this. */
1250  float dvec_null[3][3], dvec_pos[3][3], dvec_neg[3][3];
1251  float f[3];
1252  int a, b;
1253 
1254  zero_m3(dvec_null);
1255  unit_m3(dvec_pos);
1256  mul_m3_fl(dvec_pos, delta * 0.5f);
1257  copy_m3_m3(dvec_neg, dvec_pos);
1258  negate_m3(dvec_neg);
1259 
1260  /* XXX TODO: offset targets to account for position dependency. */
1261 
1262  for (a = 0; a < 3; a++) {
1263  spring_angbend_forces(
1264  data, i, j, k, goal, stiffness, damping, q, dvec_null[a], dvec_pos[a], f);
1265  copy_v3_v3(dfdv[a], f);
1266 
1267  spring_angbend_forces(
1268  data, i, j, k, goal, stiffness, damping, q, dvec_null[a], dvec_neg[a], f);
1269  sub_v3_v3(dfdv[a], f);
1270 
1271  for (b = 0; b < 3; b++) {
1272  dfdv[a][b] /= delta;
1273  }
1274  }
1275 }
1276 
1277 /* Angular spring that pulls the vertex toward the local target
1278  * See "Artistic Simulation of Curly Hair" (Pixar technical memo #12-03a)
1279  */
1280 bool SIM_mass_spring_force_spring_bending_angular(Implicit_Data *data,
1281  int i,
1282  int j,
1283  int k,
1284  const float target[3],
1285  float stiffness,
1286  float damping)
1287 {
1288  float goal[3];
1289  float fj[3], fk[3];
1290  float dfj_dxi[3][3], dfj_dxj[3][3], dfk_dxi[3][3], dfk_dxj[3][3], dfk_dxk[3][3];
1291  float dfj_dvi[3][3], dfj_dvj[3][3], dfk_dvi[3][3], dfk_dvj[3][3], dfk_dvk[3][3];
1292 
1293  const float vecnull[3] = {0.0f, 0.0f, 0.0f};
1294 
1295  world_to_root_v3(data, j, goal, target);
1296 
1297  spring_angbend_forces(data, i, j, k, goal, stiffness, damping, k, vecnull, vecnull, fk);
1298  negate_v3_v3(fj, fk); /* counterforce */
1299 
1300  spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, i, dfk_dxi);
1301  spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, j, dfk_dxj);
1302  spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, k, dfk_dxk);
1303  copy_m3_m3(dfj_dxi, dfk_dxi);
1304  negate_m3(dfj_dxi);
1305  copy_m3_m3(dfj_dxj, dfk_dxj);
1306  negate_m3(dfj_dxj);
1307 
1308  spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, i, dfk_dvi);
1309  spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, j, dfk_dvj);
1310  spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, k, dfk_dvk);
1311  copy_m3_m3(dfj_dvi, dfk_dvi);
1312  negate_m3(dfj_dvi);
1313  copy_m3_m3(dfj_dvj, dfk_dvj);
1314  negate_m3(dfj_dvj);
1315 
1316  /* add forces and jacobians to the solver data */
1317 
1318  add_v3_v3(data->F.v3(j), fj);
1319  add_v3_v3(data->F.v3(k), fk);
1320 
1321  data->idFdX.add(j, j, dfj_dxj);
1322  data->idFdX.add(k, k, dfk_dxk);
1323 
1324  data->idFdX.add(i, j, dfj_dxi);
1325  data->idFdX.add(j, i, dfj_dxi);
1326  data->idFdX.add(j, k, dfk_dxj);
1327  data->idFdX.add(k, j, dfk_dxj);
1328  data->idFdX.add(i, k, dfk_dxi);
1329  data->idFdX.add(k, i, dfk_dxi);
1330 
1331  data->idFdV.add(j, j, dfj_dvj);
1332  data->idFdV.add(k, k, dfk_dvk);
1333 
1334  data->idFdV.add(i, j, dfj_dvi);
1335  data->idFdV.add(j, i, dfj_dvi);
1336  data->idFdV.add(j, k, dfk_dvj);
1337  data->idFdV.add(k, j, dfk_dvj);
1338  data->idFdV.add(i, k, dfk_dvi);
1339  data->idFdV.add(k, i, dfk_dvi);
1340 
1341  /* XXX analytical calculation of derivatives below is incorrect.
1342  * This proved to be difficult, but for now just using the finite difference method for
1343  * estimating the jacobians should be sufficient.
1344  */
1345 # if 0
1346  float edge_ij[3], dir_ij[3], grad_dir_ij[3][3];
1347  float edge_jk[3], dir_jk[3], grad_dir_jk[3][3];
1348  float dist[3], vel_jk[3], vel_jk_ortho[3], projvel[3];
1349  float target[3];
1350  float tmp[3][3];
1351  float fi[3], fj[3], fk[3];
1352  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];
1353  float dfdvi[3][3];
1354 
1355  /* TESTING */
1356  damping = 0.0f;
1357 
1358  zero_v3(fi);
1359  zero_v3(fj);
1360  zero_v3(fk);
1361  zero_m3(dfi_dxi);
1362  zero_m3(dfj_dxi);
1363  zero_m3(dfk_dxi);
1364  zero_m3(dfk_dxj);
1365  zero_m3(dfk_dxk);
1366 
1367  /* jacobian of direction vectors */
1368  spring_grad_dir(data, i, j, edge_ij, dir_ij, grad_dir_ij);
1369  spring_grad_dir(data, j, k, edge_jk, dir_jk, grad_dir_jk);
1370 
1371  sub_v3_v3v3(vel_jk, data->V[k], data->V[j]);
1372 
1373  /* bending force */
1374  mul_v3_v3fl(target, dir_ij, restlen);
1375  sub_v3_v3v3(dist, target, edge_jk);
1376  mul_v3_v3fl(fk, dist, stiffness);
1377 
1378  /* damping force */
1379  madd_v3_v3v3fl(vel_jk_ortho, vel_jk, dir_jk, -dot_v3v3(vel_jk, dir_jk));
1380  madd_v3_v3fl(fk, vel_jk_ortho, damping);
1381 
1382  /* XXX this only holds true as long as we assume straight rest shape!
1383  * eventually will become a bit more involved since the opposite segment
1384  * gets its own target, under condition of having equal torque on both sides.
1385  */
1386  copy_v3_v3(fi, fk);
1387 
1388  /* counterforce on the middle point */
1389  sub_v3_v3(fj, fi);
1390  sub_v3_v3(fj, fk);
1391 
1392  /* === derivatives === */
1393 
1394  madd_m3_m3fl(dfk_dxi, grad_dir_ij, stiffness * restlen);
1395 
1396  madd_m3_m3fl(dfk_dxj, grad_dir_ij, -stiffness * restlen);
1397  madd_m3_m3fl(dfk_dxj, I, stiffness);
1398 
1399  madd_m3_m3fl(dfk_dxk, I, -stiffness);
1400 
1401  copy_m3_m3(dfi_dxi, dfk_dxk);
1402  negate_m3(dfi_dxi);
1403 
1404  /* dfj_dfi == dfi_dfj due to symmetry,
1405  * dfi_dfj == dfk_dfj due to fi == fk
1406  * XXX see comment above on future bent rest shapes
1407  */
1408  copy_m3_m3(dfj_dxi, dfk_dxj);
1409 
1410  /* dfj_dxj == -(dfi_dxj + dfk_dxj) due to fj == -(fi + fk) */
1411  sub_m3_m3m3(dfj_dxj, dfj_dxj, dfj_dxi);
1412  sub_m3_m3m3(dfj_dxj, dfj_dxj, dfk_dxj);
1413 
1414  /* add forces and jacobians to the solver data */
1415  add_v3_v3(data->F[i], fi);
1416  add_v3_v3(data->F[j], fj);
1417  add_v3_v3(data->F[k], fk);
1418 
1419  add_m3_m3m3(data->dFdX[i].m, data->dFdX[i].m, dfi_dxi);
1420  add_m3_m3m3(data->dFdX[j].m, data->dFdX[j].m, dfj_dxj);
1421  add_m3_m3m3(data->dFdX[k].m, data->dFdX[k].m, dfk_dxk);
1422 
1423  add_m3_m3m3(data->dFdX[block_ij].m, data->dFdX[block_ij].m, dfj_dxi);
1424  add_m3_m3m3(data->dFdX[block_jk].m, data->dFdX[block_jk].m, dfk_dxj);
1425  add_m3_m3m3(data->dFdX[block_ik].m, data->dFdX[block_ik].m, dfk_dxi);
1426 # endif
1427 
1428  return true;
1429 }
1430 
1432  int i,
1433  const float goal_x[3],
1434  const float goal_v[3],
1435  float stiffness,
1436  float damping,
1437  float r_f[3],
1438  float r_dfdx[3][3],
1439  float r_dfdv[3][3])
1440 {
1441  float root_goal_x[3], root_goal_v[3], extent[3], length, dir[3], vel[3];
1442  float f[3], dfdx[3][3], dfdv[3][3];
1443 
1444  /* goal is in world space */
1445  world_to_root_v3(data, i, root_goal_x, goal_x);
1446  world_to_root_v3(data, i, root_goal_v, goal_v);
1447 
1448  sub_v3_v3v3(extent, root_goal_x, data->X.v3(i));
1449  sub_v3_v3v3(vel, root_goal_v, data->V.v3(i));
1450  length = normalize_v3_v3(dir, extent);
1451 
1452  if (length > ALMOST_ZERO) {
1453  mul_v3_v3fl(f, dir, stiffness * length);
1454 
1455  /* Ascher & Boxman, p.21: Damping only during elongation
1456  * something wrong with it... */
1457  madd_v3_v3fl(f, dir, damping * dot_v3v3(vel, dir));
1458 
1459  dfdx_spring(dfdx, dir, length, 0.0f, stiffness);
1460  dfdv_damp(dfdv, dir, damping);
1461 
1462  add_v3_v3(data->F.v3(i), f);
1463  data->idFdX.add(i, i, dfdx);
1464  data->idFdV.add(i, i, dfdv);
1465 
1466  if (r_f) {
1467  copy_v3_v3(r_f, f);
1468  }
1469  if (r_dfdx) {
1470  copy_m3_m3(r_dfdx, dfdx);
1471  }
1472  if (r_dfdv) {
1473  copy_m3_m3(r_dfdv, dfdv);
1474  }
1475 
1476  return true;
1477  }
1478  else {
1479  if (r_f) {
1480  zero_v3(r_f);
1481  }
1482  if (r_dfdx) {
1483  zero_m3(r_dfdx);
1484  }
1485  if (r_dfdv) {
1486  zero_m3(r_dfdv);
1487  }
1488 
1489  return false;
1490  }
1491 }
1492 
1493 #endif /* IMPLICIT_SOLVER_EIGEN */
typedef float(TangentPoint)[2]
#define ALMOST_ZERO
Definition: BKE_cloth.h:31
#define BLI_INLINE
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 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 madd_m3_m3m3fl(float R[3][3], const float A[3][3], const float B[3][3], float f)
Definition: math_matrix.c:1054
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 float dot_v3v3(const float a[3], const float b[3]) ATTR_WARN_UNUSED_RESULT
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 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 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 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.
int SIM_mass_spring_solver_numvert(struct Implicit_Data *id)
@ SIM_SOLVER_SUCCESS
@ SIM_SOLVER_INVALID_INPUT
@ SIM_SOLVER_NUMERICAL_ISSUE
@ 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
btGeneric6DofConstraint & operator=(btGeneric6DofConstraint &other)
SIMD_FORCE_INLINE const btScalar & w() const
Return the w value.
Definition: btQuadWord.h:119
A conjugate gradient solver for sparse self-adjoint problems with additional constraints.
BLI_INLINE void madd_m3_m3fl(float r[3][3], const float m[3][3], float f)
Definition: cloth.c:1249
#define powf(x, y)
Definition: cuda/compat.h:103
SyclQueue void void size_t num_bytes void
Eigen::ConjugateGradient< lMatrix, Eigen::Lower, Eigen::DiagonalPreconditioner< Scalar > > ConjugateGradient
Definition: eigen_utils.h:186
std::vector< Triplet > TripletList
Definition: eigen_utils.h:127
Eigen::VectorXf lVector
Definition: eigen_utils.h:96
BLI_INLINE void print_lmatrix(const lMatrix &m)
Definition: eigen_utils.h:201
Eigen::SparseMatrix< Scalar > lMatrix
Definition: eigen_utils.h:129
BLI_INLINE void print_lvector(const lVector3f &v)
Definition: eigen_utils.h:190
float Scalar
Definition: eigen_utils.h:26
Eigen::Triplet< Scalar > Triplet
Definition: eigen_utils.h:126
uint nor
void SIM_mass_spring_add_constraint_ndof1(struct Implicit_Data *data, int index, const float c1[3], const float c2[3], const float dV[3])
void SIM_mass_spring_get_motion_state(struct Implicit_Data *data, int index, float x[3], float v[3])
void SIM_mass_spring_add_constraint_ndof0(struct Implicit_Data *data, int index, const float dV[3])
void SIM_mass_spring_force_reference_frame(struct Implicit_Data *data, int index, const float acceleration[3], const float omega[3], const float domega_dt[3], float mass)
void SIM_mass_spring_force_edge_wind(struct Implicit_Data *data, int v1, int v2, float radius1, float radius2, const float(*winvec)[3])
void SIM_mass_spring_set_new_velocity(struct Implicit_Data *data, int index, const float v[3])
bool SIM_mass_spring_force_spring_goal(struct Implicit_Data *data, int i, const float goal_x[3], const float goal_v[3], float stiffness, float damping)
void SIM_mass_spring_set_vertex_mass(struct Implicit_Data *data, int index, float mass)
void SIM_mass_spring_set_motion_state(struct Implicit_Data *data, int index, const float x[3], const float v[3])
bool SIM_mass_spring_force_spring_linear(struct 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)
bool SIM_mass_spring_force_spring_bending(struct Implicit_Data *data, int i, int j, float restlen, float kb, float cb)
void SIM_mass_spring_apply_result(struct Implicit_Data *data)
void SIM_mass_spring_set_position(struct Implicit_Data *data, int index, const float x[3])
void SIM_mass_spring_add_constraint_ndof2(struct Implicit_Data *data, int index, const float c1[3], const float dV[3])
BLI_INLINE void implicit_print_matrix_elem(float v)
Definition: implicit.h:46
void SIM_mass_spring_force_face_wind(struct Implicit_Data *data, int v1, int v2, int v3, const float(*winvec)[3])
void SIM_mass_spring_get_position(struct Implicit_Data *data, int index, float x[3])
bool SIM_mass_spring_solve_velocities(struct Implicit_Data *data, float dt, struct ImplicitSolverResult *result)
void SIM_mass_spring_force_drag(struct Implicit_Data *data, float drag)
void SIM_mass_spring_clear_constraints(struct Implicit_Data *data)
bool SIM_mass_spring_solve_positions(struct Implicit_Data *data, float dt)
void SIM_mass_spring_get_new_velocity(struct Implicit_Data *data, int index, float v[3])
void SIM_mass_spring_set_velocity(struct Implicit_Data *data, int index, const float v[3])
void SIM_mass_spring_force_extern(struct Implicit_Data *data, int i, const float f[3], float dfdx[3][3], float dfdv[3][3])
void SIM_mass_spring_force_gravity(struct Implicit_Data *data, int index, float mass, const float g[3])
void SIM_mass_spring_set_rest_transform(struct Implicit_Data *data, int index, float tfm[3][3])
void SIM_mass_spring_clear_forces(struct 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])
BLI_INLINE float fbstar_jacobi(float length, float L, float kb, float cb)
DO_INLINE void mul_fvectorT_fvector(float to[3][3], const float vectorA[3], const float vectorB[3])
BLI_INLINE void cross_m3_v3m3(float r[3][3], const float v[3], const float m[3][3])
BLI_INLINE float fb(float length, float L)
DO_INLINE void mul_fmatrix_S(float matrix[3][3], float scalar)
Implicit_Data * SIM_mass_spring_solver_create(int numverts, int numsprings)
BLI_INLINE float fbderiv(float length, float L)
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])
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])
BLI_INLINE void outerproduct(float r[3][3], const float a[3], const float b[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
BLI_INLINE void dfdx_spring(float to[3][3], const float dir[3], float length, float L, float k)
BLI_INLINE void spring_grad_dir(Implicit_Data *data, int i, int j, float edge[3], float dir[3], float grad_dir[3][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 world_to_root_m3(Implicit_Data *data, int index, float r[3][3], const float m[3][3])
BLI_INLINE float fbstar(float length, float L, float kb, float cb)
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)
ccl_gpu_kernel_postfix ccl_global float int int int int float bool reset
clear internal cached data and reset random seed
#define T
#define L
static unsigned a[3]
Definition: RandGen.cpp:78
bool add(void *owner, const AttributeIDRef &attribute_id, eAttrDomain domain, eCustomDataType data_type, const AttributeInit &initializer)
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)
#define I
fmatrix3x3 * M
fmatrix3x3 * A
fmatrix3x3 * S
fmatrix3x3 * dFdX
fmatrix3x3 * dFdV
fmatrix3x3 * tfm