1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
#include "vec3n.h"
//#include "console.h"
extern int numX;
//
// Cloth - Backward Integrated Spring Network
//
// (c) Stan Melax 2006
// http://www.melax.com/cloth
// freeware demo and source
// Although its free software, I'll gaurantee and support this software as much as is reasonable.
// However, if you choose to use any of this code, then you agree that
// I assume no financial liability should the software not meet your expectations.
// But do feel free to send any feedback.
//
// The core backward integration functionality has all been extracted into the SpringNetwork class.
// This makes it easy for you if you just want to look at or use the math and the algorithms.
// The remainder of the code builds a cloth system with basic render support, I/O, and manipulators,
// so its possible to make use of the technology within a 3D application.
// This code is separated from the SpringNetwork class in order to avoid pushing a particular style
// and prevent any dependancies of the algorithms onto unrelated systems.
// Feel free to adapt any of this into your own 3D engine/environment.
//
// Instead of having unique Hooke force and damping coefficients on each spring, the SpringNetwork
// code uses a spring 'type' that indexes a short list of shared named coefficients.
// This was just more practical for the typical application of this technology.
// Over-designed systems that are too general can be slower for
// the next guy to understand and more painful to use.
// Editing/creation is easier when only 1 number needs to be changed.
// Nonetheless, feel free to adapt to your own needs.
//
#include <stdio.h>
#include <float.h>
#include "vec3n.h"
//#include "console.h"
//#include "manipulatori.h"
//#include "object.h"
//#include "xmlparse.h"
static const float3x3 I(1,0,0,0,1,0,0,0,1);
inline float3x3 dfdx_spring(const float3 &dir,float length,float rest,float k)
{
// dir is unit length direction, rest is spring's restlength, k is spring constant.
return ( (I-outerprod(dir,dir))*Min(1.0f,rest/length) - I) * -k;
}
inline float3x3 dfdx_damp(const float3 &dir,float length,const float3& vel,float rest,float damping)
{
// inner spring damping vel is the relative velocity of the endpoints.
return (I-outerprod(dir,dir)) * (-damping * -(dot(dir,vel)/Max(length,rest)));
}
inline float3x3 dfdv_damp(const float3 &dir,float damping)
{
// derivative of force wrt velocity.
return outerprod(dir,dir) * damping;
}
#include "SpringNetwork.h"
SpringNetwork::SpringNetwork(int _n):X(_n),V(_n),F(_n),dV(_n),A(_n),dFdX(_n),dFdV(_n)
{
assert(SPRING_STRUCT==0);
assert(&spring_shear == &spring_struct +SPRING_SHEAR);
assert(&spring_bend == &spring_struct +SPRING_BEND);
assert(&spring_struct== &spring_k[SPRING_STRUCT]);
assert(&spring_shear == &spring_k[SPRING_SHEAR ]);
assert(&spring_bend == &spring_k[SPRING_BEND ]);
// spring_struct=1000000.0f;
// spring_shear=1000000.0f;
spring_struct=1000.0f;
spring_shear=100.0f;
spring_bend=25.0f;
spring_damp=5.0f;
spring_air=1.0f;
spring_air=1.0f;
cloth_step = 0.25f; // delta time for cloth
cloth_gravity=float3(0,-10,0);
cloth_sleepthreshold = 0.001f;
cloth_sleepcount = 100;
awake = cloth_sleepcount;
//fix/pin two points in worldspace
float3Nx3N::Block zero;
zero.m = float3x3(0,0,0,0,0,0,0,0,0);
zero.c = 0;
zero.r = 0;
S.blocks.Add(zero);
zero.r = numX-1;
S.blocks.Add(zero);
}
SpringNetwork::Spring &SpringNetwork::AddBlocks(Spring &s)
{
// Called during initial creation of springs in our spring network.
// Sets up the sparse matrices corresponding to connections.
// Note the indices (s.iab,s.iba) are also stored with spring to avoid looking them up each time a spring is applied
// All 3 matrices A,dFdX, and dFdV are contstructed identically so the block array layout will be the same for each.
s.iab = A.blocks.count; // added 'ab' blocks will have this index.
A.blocks.Add(float3Nx3N::Block(s.a,s.b));
dFdX.blocks.Add(float3Nx3N::Block(s.a,s.b));
dFdV.blocks.Add(float3Nx3N::Block(s.a,s.b));
s.iba = A.blocks.count; // added 'ba' blocks will have this index.
A.blocks.Add(float3Nx3N::Block(s.b,s.a));
dFdX.blocks.Add(float3Nx3N::Block(s.b,s.a));
dFdV.blocks.Add(float3Nx3N::Block(s.b,s.a));
return s;
}
void SpringNetwork::PreSolveSpring(const SpringNetwork::Spring &s)
{
// Adds this spring's contribution into force vector F and force derivitves dFdX and dFdV
// One optimization would be premultiply dfdx by dt*dt and F and dFdV by dt right here in this function.
// However, for educational purposes we wont do that now and intead just follow the paper directly.
//assert(dFdX.blocks[s.a].c==s.a); // delete this assert, no bugs here
//assert(dFdX.blocks[s.a].r==s.a);
float3 extent = X[s.b] - X[s.a];
float length = magnitude(extent);
float3 dir = (length==0)?float3(0,0,0): extent * 1.0f/length;
float3 vel = V[s.b] - V[s.a];
float k = spring_k[s.type];
float3 f = dir * ((k * (length-s.restlen) ) + spring_damp * dot(vel,dir)); // spring force + damping force
F[s.a] += f;
F[s.b] -= f;
float3x3 dfdx = dfdx_spring(dir,length,s.restlen,k) + dfdx_damp(dir,length,vel,s.restlen,spring_damp);
dFdX.blocks[s.a].m -= dfdx; // diagonal chunk dFdX[a,a]
dFdX.blocks[s.b].m -= dfdx; // diagonal chunk dFdX[b,b]
dFdX.blocks[s.iab].m += dfdx; // off-diag chunk dFdX[a,b]
dFdX.blocks[s.iba].m += dfdx; // off-diag chunk dFdX[b,a]
float3x3 dfdv = dfdv_damp(dir,spring_damp);
dFdV.blocks[s.a].m -= dfdv; // diagonal chunk dFdV[a,a]
dFdV.blocks[s.b].m -= dfdv; // diagonal chunk dFdV[b,b]
dFdV.blocks[s.iab].m += dfdv; // off-diag chunk dFdV[a,b]
dFdV.blocks[s.iba].m += dfdv; // off-diag chunk dFdV[b,a]
}
void SpringNetwork::CalcForces()
{
// Collect forces and derivatives: F,dFdX,dFdV
dFdX.Zero();
dFdV.InitDiagonal(-spring_air);
F.Init(cloth_gravity);
F.element[0]=float3(0,0,0);
F.element[numX-1]=float3(0,0,0);
F -= V * spring_air;
for(int i=0;i<springs.count;i++)
{
PreSolveSpring(springs[i]); // will add to F,dFdX,dFdV
}
}
void SpringNetwork::Simulate(float dt)
{
// Get ready for conjugate gradient iterative solver step.
// Initialize operands.
if(!awake) return;
CalcForces();
int n=X.count; // all our big vectors are of this size
float3N dFdXmV(n); // temp to store result of matrix multiply
float3N B(n);
dV.Zero();
A.Identity(); // build up the big matrix we feed to solver
A -= dFdV * dt + dFdX * (dt*dt) ;
dFdXmV = dFdX * V;
B = F * dt + dFdXmV * (dt*dt);
ConjGradientFiltered(dV,A,B,S);
V = V + dV;
// V.element[0] = float3(0,0,0);
// V.element[numX-1] = float3(0,0,0);
X = X + V*dt;
UpdateLimits();
awake = (dot(V,V)<cloth_sleepthreshold)?awake-1:awake=cloth_sleepcount;
}
|
