Skip to content

Commit 747e6ce

Browse files
RoushelfyRoushelfy
committed
opt
1 parent 593bc99 commit 747e6ce

File tree

8 files changed

+487
-460
lines changed

8 files changed

+487
-460
lines changed

‎simulators/1_mass_spring/src/simulator.cu‎

Lines changed: 44 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -12,8 +12,9 @@ struct MassSpringSimulator<T, dim>::Impl
1212
int n_seg;
1313
T h, rho, side_len, initial_stretch, m, tol;
1414
int resolution = 900, scale = 200, offset = resolution / 2, radius = 5;
15-
std::vector<T> x, x_tilde, v, k, l2;
16-
std::vector<int> e;
15+
std::vector<T> host_x, host_k, host_l2;
16+
std::vector<int> host_e;
17+
DeviceBuffer<T> device_x, device_v;
1718
sf::RenderWindow window;
1819
InertialEnergy<T, dim> inertialenergy;
1920
MassSpringEnergy<T, dim> massspringenergy;
@@ -49,29 +50,29 @@ MassSpringSimulator<T, dim>::MassSpringSimulator(T rho, T side_len, T initial_st
4950
template <typename T, int dim>
5051
MassSpringSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch, T K, T h_, T tol_, int n_seg) : tol(tol_), h(h_), window(sf::VideoMode(resolution, resolution), "MassSpringSimulator")
5152
{
52-
generate(side_len, n_seg, x, e);
53-
v.resize(x.size(), 0);
54-
k.resize(e.size() / 2, K);
55-
l2.resize(e.size() / 2);
56-
for (int i = 0; i < e.size() / 2; i++)
53+
generate(side_len, n_seg, host_x, host_e);
54+
host_k.resize(host_e.size() / 2, K);
55+
host_l2.resize(host_e.size() / 2);
56+
device_x.resize(host_x.size(), 0);
57+
device_v.resize(host_x.size(), 0);
58+
for (int i = 0; i < host_e.size() / 2; i++)
5759
{
5860
T diff = 0;
59-
int idx1 = e[2 * i], idx2 = e[2 * i + 1];
61+
int idx1 = host_e[2 * i], idx2 = host_e[2 * i + 1];
6062
for (int d = 0; d < dim; d++)
6163
{
62-
diff += (x[idx1 * dim + d] - x[idx2 * dim + d]) * (x[idx1 * dim + d] - x[idx2 * dim + d]);
64+
diff += (host_x[idx1 * dim + d] - host_x[idx2 * dim + d]) * (host_x[idx1 * dim + d] - host_x[idx2 * dim + d]);
6365
}
64-
l2[i] = diff;
66+
host_l2[i] = diff;
6567
}
6668
m = rho * side_len * side_len / ((n_seg + 1) * (n_seg + 1));
6769
// initial stretch
68-
int N = x.size() / dim;
70+
int N = host_x.size() / dim;
6971
for (int i = 0; i < N; i++)
70-
x[i * dim + 0] *= initial_stretch;
72+
host_x[i * dim + 0] *= initial_stretch;
7173
inertialenergy = InertialEnergy<T, dim>(N, m);
72-
massspringenergy = MassSpringEnergy<T, dim>(x, e, l2, k);
73-
DeviceBuffer<T> x_device(x);
74-
update_x(x_device);
74+
massspringenergy = MassSpringEnergy<T, dim>(host_x, host_e, host_l2, host_k);
75+
device_x.copy_from(host_x);
7576
}
7677
template <typename T, int dim>
7778
void MassSpringSimulator<T, dim>::run()
@@ -105,36 +106,35 @@ void MassSpringSimulator<T, dim>::run()
105106
template <typename T, int dim>
106107
void MassSpringSimulator<T, dim>::Impl::step_forward()
107108
{
108-
DeviceBuffer<T> x_tilde(x.size()); // Predictive position
109-
update_x_tilde(add_vector<T>(x, v, 1, h));
110-
DeviceBuffer<T> x_n = x; // Copy current positions to x_n
109+
update_x_tilde(add_vector<T>(device_x, device_v, 1, h));
110+
DeviceBuffer<T> device_x_n = device_x; // Copy current positions to device_x_n
111111
int iter = 0;
112112
T E_last = IP_val();
113-
DeviceBuffer<T> p = search_direction();
114-
T residual = max_vector(p) / h;
113+
DeviceBuffer<T> device_p = search_direction();
114+
T residual = max_vector(device_p) / h;
115115
while (residual > tol)
116116
{
117117
std::cout << "Iteration " << iter << " residual " << residual << "E_last" << E_last << "\n";
118118
// Line search
119119
T alpha = 1;
120-
DeviceBuffer<T> x0 = x;
121-
update_x(add_vector<T>(x0, p, 1.0, alpha));
120+
DeviceBuffer<T> device_x0 = device_x;
121+
update_x(add_vector<T>(device_x0, device_p, 1.0, alpha));
122122
while (IP_val() > E_last)
123123
{
124124
alpha /= 2;
125-
update_x(add_vector<T>(x0, p, 1.0, alpha));
125+
update_x(add_vector<T>(device_x0, device_p, 1.0, alpha));
126126
}
127127
std::cout << "step size = " << alpha << "\n";
128128
E_last = IP_val();
129-
p = search_direction();
130-
residual = max_vector(p) / h;
129+
device_p = search_direction();
130+
residual = max_vector(device_p) / h;
131131
iter += 1;
132132
}
133-
update_v(add_vector<T>(x, x_n, 1 / h, -1 / h));
133+
update_v(add_vector<T>(device_x, device_x_n, 1 / h, -1 / h));
134134
}
135135
// ANCHOR_END: step_forward
136136

137-
.template <typename T, int dim>
137+
template <typename T, int dim>
138138
T MassSpringSimulator<T, dim>::Impl::screen_projection_x(T point)
139139
{
140140
return offset + scale * point;
@@ -144,57 +144,64 @@ T MassSpringSimulator<T, dim>::Impl::screen_projection_y(T point)
144144
{
145145
return resolution - (offset + scale * point);
146146
}
147+
148+
// ANCHOR: update_x
147149
template <typename T, int dim>
148150
void MassSpringSimulator<T, dim>::Impl::update_x(const DeviceBuffer<T> &new_x)
149151
{
150152
inertialenergy.update_x(new_x);
151153
massspringenergy.update_x(new_x);
152-
new_x.copy_to(x);
154+
device_x = new_x;
153155
}
156+
// ANCHOR_END: update_x
157+
154158
template <typename T, int dim>
155159
void MassSpringSimulator<T, dim>::Impl::update_x_tilde(const DeviceBuffer<T> &new_x_tilde)
156160
{
157161
inertialenergy.update_x_tilde(new_x_tilde);
158-
new_x_tilde.copy_to(x_tilde);
159162
}
160163
template <typename T, int dim>
161164
void MassSpringSimulator<T, dim>::Impl::update_v(const DeviceBuffer<T> &new_v)
162165
{
163-
new_v.copy_to(v);
166+
device_v = new_v;
164167
}
165168
template <typename T, int dim>
166169
void MassSpringSimulator<T, dim>::Impl::draw()
167170
{
171+
device_x.copy_to(host_x);
168172
window.clear(sf::Color::White); // Clear the previous frame
169173

170174
// Draw springs as lines
171-
for (int i = 0; i < e.size() / 2; ++i)
175+
for (int i = 0; i < host_e.size() / 2; ++i)
172176
{
173177
sf::Vertex line[] = {
174-
sf::Vertex(sf::Vector2f(screen_projection_x(x[e[i * 2] * dim]), screen_projection_y(x[e[i * 2] * dim + 1])), sf::Color::Blue),
175-
sf::Vertex(sf::Vector2f(screen_projection_x(x[e[i * 2 + 1] * dim]), screen_projection_y(x[e[i * 2 + 1] * dim + 1])), sf::Color::Blue)};
178+
sf::Vertex(sf::Vector2f(screen_projection_x(host_x[host_e[i * 2] * dim]), screen_projection_y(host_x[host_e[i * 2] * dim + 1])), sf::Color::Blue),
179+
sf::Vertex(sf::Vector2f(screen_projection_x(host_x[host_e[i * 2 + 1] * dim]), screen_projection_y(host_x[host_e[i * 2 + 1] * dim + 1])), sf::Color::Blue)};
176180
window.draw(line, 2, sf::Lines);
177181
}
178182

179183
// Draw masses as circles
180-
for (int i = 0; i < x.size() / dim; ++i)
184+
for (int i = 0; i < host_x.size() / dim; ++i)
181185
{
182186
sf::CircleShape circle(radius); // Set a fixed radius for each mass
183187
circle.setFillColor(sf::Color::Red);
184-
circle.setPosition(screen_projection_x(x[i * dim]) - radius, screen_projection_y(x[i * dim + 1]) - radius); // Center the circle on the mass
188+
circle.setPosition(screen_projection_x(host_x[i * dim]) - radius, screen_projection_y(host_x[i * dim + 1]) - radius); // Center the circle on the mass
185189
window.draw(circle);
186190
}
187191

188192
window.display(); // Display the rendered frame
189193
}
190194

195+
// ANCHOR: IP_val
191196
template <typename T, int dim>
192197
T MassSpringSimulator<T, dim>::Impl::IP_val()
193198
{
194199

195200
return inertialenergy.val() + massspringenergy.val() * h * h;
196201
}
202+
// ANCHOR_END: IP_val
197203

204+
// ANCHOR: IP_grad and IP_hess
198205
template <typename T, int dim>
199206
DeviceBuffer<T> MassSpringSimulator<T, dim>::Impl::IP_grad()
200207
{
@@ -209,11 +216,12 @@ DeviceTripletMatrix<T, 1> MassSpringSimulator<T, dim>::Impl::IP_hess()
209216
DeviceTripletMatrix<T, 1> hess = add_triplet<T>(inertial_hess, massspring_hess, 1.0, h * h);
210217
return hess;
211218
}
219+
// ANCHOR_END: IP_grad and IP_hess
212220
template <typename T, int dim>
213221
DeviceBuffer<T> MassSpringSimulator<T, dim>::Impl::search_direction()
214222
{
215223
DeviceBuffer<T> dir;
216-
dir.resize(x.size());
224+
dir.resize(host_x.size());
217225
search_dir(IP_grad(), IP_hess(), dir);
218226
return dir;
219227
}

‎simulators/2_dirichlet/src/simulator.cu‎

Lines changed: 36 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -13,8 +13,9 @@ struct DirichletSimulator<T, dim>::Impl
1313
int n_seg;
1414
T h, rho, side_len, initial_stretch, m, tol;
1515
int resolution = 900, scale = 200, offset = resolution / 2, radius = 5;
16-
std::vector<T> x, x_tilde, v, k, l2;
17-
std::vector<int> e;
16+
std::vector<T> host_x, host_k, host_l2;
17+
std::vector<int> host_e;
18+
DeviceBuffer<T> device_x, device_v;
1819
DeviceBuffer<int> device_DBC;
1920
sf::RenderWindow window;
2021
InertialEnergy<T, dim> inertialenergy;
@@ -52,34 +53,34 @@ DirichletSimulator<T, dim>::DirichletSimulator(T rho, T side_len, T initial_stre
5253
template <typename T, int dim>
5354
DirichletSimulator<T, dim>::Impl::Impl(T rho, T side_len, T initial_stretch, T K, T h_, T tol_, int n_seg) : tol(tol_), h(h_), window(sf::VideoMode(resolution, resolution), "DirichletSimulator")
5455
{
55-
generate(side_len, n_seg, x, e);
56-
std::vector<int> DBC(x.size() / dim, 0);
56+
generate(side_len, n_seg, host_x, host_e);
57+
std::vector<int> DBC(host_x.size() / dim, 0);
5758
DBC[n_seg] = 1;
5859
DBC[(n_seg + 1) * (n_seg + 1) - 1] = 1;
59-
v.resize(x.size(), 0);
60-
k.resize(e.size() / 2, K);
61-
l2.resize(e.size() / 2);
62-
for (int i = 0; i < e.size() / 2; i++)
60+
device_x.resize(host_x.size(), 0);
61+
device_v.resize(host_x.size(), 0);
62+
host_k.resize(host_e.size() / 2, K);
63+
host_l2.resize(host_e.size() / 2);
64+
for (int i = 0; i < host_e.size() / 2; i++)
6365
{
6466
T diff = 0;
65-
int idx1 = e[2 * i], idx2 = e[2 * i + 1];
67+
int idx1 = host_e[2 * i], idx2 = host_e[2 * i + 1];
6668
for (int d = 0; d < dim; d++)
6769
{
68-
diff += (x[idx1 * dim + d] - x[idx2 * dim + d]) * (x[idx1 * dim + d] - x[idx2 * dim + d]);
70+
diff += (host_x[idx1 * dim + d] - host_x[idx2 * dim + d]) * (host_x[idx1 * dim + d] - host_x[idx2 * dim + d]);
6971
}
70-
l2[i] = diff;
72+
host_l2[i] = diff;
7173
}
7274
m = rho * side_len * side_len / ((n_seg + 1) * (n_seg + 1));
7375
// initial stretch
74-
int N = x.size() / dim;
76+
int N = host_x.size() / dim;
7577
for (int i = 0; i < N; i++)
76-
x[i * dim + 0] *= initial_stretch;
78+
host_x[i * dim + 0] *= initial_stretch;
7779
inertialenergy = InertialEnergy<T, dim>(N, m);
78-
massspringenergy = MassSpringEnergy<T, dim>(x, e, l2, k);
80+
massspringenergy = MassSpringEnergy<T, dim>(host_x, host_e, host_l2, host_k);
7981
gravityenergy = GravityEnergy<T, dim>(N, m);
80-
DeviceBuffer<T> x_device(x);
81-
update_x(x_device);
8282
device_DBC = DeviceBuffer<int>(DBC);
83+
device_x.copy_from(host_x);
8384
}
8485
template <typename T, int dim>
8586

@@ -110,34 +111,33 @@ void DirichletSimulator<T, dim>::run()
110111
template <typename T, int dim>
111112
void DirichletSimulator<T, dim>::Impl::step_forward()
112113
{
113-
DeviceBuffer<T> x_tilde(x.size()); // Predictive position
114-
update_x_tilde(add_vector<T>(x, v, 1, h));
115-
DeviceBuffer<T> x_n = x; // Copy current positions to x_n
114+
update_x_tilde(add_vector<T>(device_x, device_v, 1, h));
115+
DeviceBuffer<T> device_x_n = device_x; // Copy current positions to device_x_n
116116
int iter = 0;
117117
T E_last = IP_val();
118-
DeviceBuffer<T> p = search_direction();
119-
T residual = max_vector(p) / h;
118+
DeviceBuffer<T> device_p = search_direction();
119+
T residual = max_vector(device_p) / h;
120120
// std::cout << "Initial residual " << residual << "\n";
121121
while (residual > tol)
122122
{
123123
std::cout << "Iteration " << iter << " residual " << residual << "E_last" << E_last << "\n";
124124
// Line search
125125
T alpha = 1;
126-
DeviceBuffer<T> x0 = x;
127-
update_x(add_vector<T>(x0, p, 1.0, alpha));
126+
DeviceBuffer<T> x0 = device_x;
127+
update_x(add_vector<T>(x0, device_p, 1.0, alpha));
128128
while (IP_val() > E_last)
129129
{
130130
alpha /= 2;
131-
update_x(add_vector<T>(x0, p, 1.0, alpha));
131+
update_x(add_vector<T>(x0, device_p, 1.0, alpha));
132132
}
133133
// std::cout << "step size = " << alpha << "\n";
134134
E_last = IP_val();
135135
// std::cout << "Iteration " << iter << " residual " << residual << "E_last" << E_last << "\n";
136-
p = search_direction();
137-
residual = max_vector(p) / h;
136+
device_p = search_direction();
137+
residual = max_vector(device_p) / h;
138138
iter += 1;
139139
}
140-
update_v(add_vector<T>(x, x_n, 1 / h, -1 / h));
140+
update_v(add_vector<T>(device_x, device_x_n, 1 / h, -1 / h));
141141
}
142142
template <typename T, int dim>
143143
T DirichletSimulator<T, dim>::Impl::screen_projection_x(T point)
@@ -155,39 +155,39 @@ void DirichletSimulator<T, dim>::Impl::update_x(const DeviceBuffer<T> &new_x)
155155
inertialenergy.update_x(new_x);
156156
massspringenergy.update_x(new_x);
157157
gravityenergy.update_x(new_x);
158-
new_x.copy_to(x);
158+
device_x = new_x;
159159
}
160160
template <typename T, int dim>
161161
void DirichletSimulator<T, dim>::Impl::update_x_tilde(const DeviceBuffer<T> &new_x_tilde)
162162
{
163163
inertialenergy.update_x_tilde(new_x_tilde);
164-
new_x_tilde.copy_to(x_tilde);
165164
}
166165
template <typename T, int dim>
167166
void DirichletSimulator<T, dim>::Impl::update_v(const DeviceBuffer<T> &new_v)
168167
{
169-
new_v.copy_to(v);
168+
device_v = new_v;
170169
}
171170
template <typename T, int dim>
172171
void DirichletSimulator<T, dim>::Impl::draw()
173172
{
173+
device_x.copy_to(host_x);
174174
window.clear(sf::Color::White); // Clear the previous frame
175175

176176
// Draw springs as lines
177-
for (int i = 0; i < e.size() / 2; ++i)
177+
for (int i = 0; i < host_e.size() / 2; ++i)
178178
{
179179
sf::Vertex line[] = {
180-
sf::Vertex(sf::Vector2f(screen_projection_x(x[e[i * 2] * dim]), screen_projection_y(x[e[i * 2] * dim + 1])), sf::Color::Blue),
181-
sf::Vertex(sf::Vector2f(screen_projection_x(x[e[i * 2 + 1] * dim]), screen_projection_y(x[e[i * 2 + 1] * dim + 1])), sf::Color::Blue)};
180+
sf::Vertex(sf::Vector2f(screen_projection_x(host_x[host_e[i * 2] * dim]), screen_projection_y(host_x[host_e[i * 2] * dim + 1])), sf::Color::Blue),
181+
sf::Vertex(sf::Vector2f(screen_projection_x(host_x[host_e[i * 2 + 1] * dim]), screen_projection_y(host_x[host_e[i * 2 + 1] * dim + 1])), sf::Color::Blue)};
182182
window.draw(line, 2, sf::Lines);
183183
}
184184

185185
// Draw masses as circles
186-
for (int i = 0; i < x.size() / dim; ++i)
186+
for (int i = 0; i < host_x.size() / dim; ++i)
187187
{
188188
sf::CircleShape circle(radius); // Set a fixed radius for each mass
189189
circle.setFillColor(sf::Color::Red);
190-
circle.setPosition(screen_projection_x(x[i * dim]) - radius, screen_projection_y(x[i * dim + 1]) - radius); // Center the circle on the mass
190+
circle.setPosition(screen_projection_x(host_x[i * dim]) - radius, screen_projection_y(host_x[i * dim + 1]) - radius); // Center the circle on the mass
191191
window.draw(circle);
192192
}
193193

@@ -219,7 +219,7 @@ template <typename T, int dim>
219219
DeviceBuffer<T> DirichletSimulator<T, dim>::Impl::search_direction()
220220
{
221221
DeviceBuffer<T> dir;
222-
dir.resize(x.size());
222+
dir.resize(device_x.size());
223223
DeviceBuffer<T> grad = IP_grad();
224224
DeviceTripletMatrix<T, 1> hess = IP_hess();
225225
search_dir<T, dim>(grad, hess, dir, device_DBC);

0 commit comments

Comments
 (0)