Skip to content

Commit c4f5188

Browse files
larsblatnylarsblatny
committed
added support for MUSL and other minor changes
1 parent 9f10811 commit c4f5188

File tree

9 files changed

+53
-157
lines changed

9 files changed

+53
-157
lines changed

‎src/sampling/sampling_particles.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -76,8 +76,8 @@
7676
void sampleParticles(S& sim, T kRadius, T ppc = 6, unsigned int crop_to_shape = 0, std::uint32_t attempts = 200, std::uint32_t seed = 42){
7777
const T Lx = sim.Lx;
7878
const T Ly = sim.Ly;
79-
std::uint32_t kAttempts = 200;
80-
std::uint32_t kSeed = 42;
79+
std::uint32_t kAttempts = attempts;
80+
std::uint32_t kSeed = seed;
8181
std::array<T, 2> kXMin = std::array<T, 2>{{0, 0}};
8282
std::array<T, 2> kXMax = std::array<T, 2>{{Lx, Ly}};
8383

‎src/simulation/deformation_update.cpp‎

Lines changed: 19 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ void Simulation::deformationUpdate(){
1616
#pragma omp parallel for reduction(+:plastic_count) num_threads(n_threads)
1717
for(int p=0; p<Np; p++){
1818

19-
TM sum = TM::Zero();
19+
TM v_grad = TM::Zero();
2020
TV xp = particles.x[p];
2121
unsigned int i_base = std::floor((xp(0)-grid.xc)*one_over_dx) - 1; // the subtraction of one is valid for both quadratic and cubic splines
2222
unsigned int j_base = std::floor((xp(1)-grid.yc)*one_over_dx) - 1;
@@ -26,21 +26,36 @@ void Simulation::deformationUpdate(){
2626

2727
for(int i = i_base; i < i_base+4; i++){
2828
T xi = grid.x[i];
29+
T wi = N((xp(0)-xi)*one_over_dx);
30+
T wi_grad = dNdu((xp(0) - xi) * one_over_dx) * one_over_dx;
2931
for(int j = j_base; j < j_base+4; j++){
3032
T yi = grid.y[j];
33+
T wj = N((xp(1) - yi)*one_over_dx);
34+
T wj_grad = dNdu((xp(1) - yi) * one_over_dx) * one_over_dx;
3135
#ifdef THREEDIM
3236
for(int k = k_base; k < k_base+4; k++){
3337
T zi = grid.z[k];
34-
sum += grid.v[ind(i,j,k)] * grad_wip(xp(0), xp(1), xp(2), xi, yi, zi, one_over_dx).transpose();
38+
T wk = N((xp(2) - zi)*one_over_dx);
39+
T wk_grad = dNdu((xp(2) - zi) * one_over_dx) * one_over_dx;
40+
T weight = wi * wj * wk;
41+
TV weight_grad;
42+
weight_grad << wi_grad*wj*wk,
43+
wi*wj_grad*wk,
44+
wi*wj*wk_grad;
45+
v_grad += grid.v[ind(i,j,k)] * weight_grad.transpose();
3546
} // end loop k
3647
#else
37-
sum += grid.v[ind(i,j)] * grad_wip(xp(0), xp(1), xi, yi, one_over_dx).transpose();
48+
T weight = wi * wj;
49+
TV weight_grad;
50+
weight_grad << wi_grad*wj,
51+
wi*wj_grad;
52+
v_grad += grid.v[ind(i,j)] * weight_grad.transpose();
3853
#endif
3954
} // end loop i
4055
} // end loop j
4156

4257
TM Fe_trial = particles.F[p];
43-
Fe_trial = Fe_trial + dt * sum * Fe_trial;
58+
Fe_trial = Fe_trial + dt * v_grad * Fe_trial;
4459
particles.F[p] = Fe_trial;
4560

4661
plasticity(p, plastic_count, Fe_trial);

‎src/simulation/explicit_euler_update.cpp‎

Lines changed: 0 additions & 60 deletions
Original file line numberDiff line numberDiff line change
@@ -10,66 +10,6 @@ void Simulation::explicitEulerUpdate(){
1010
debug("explicitEulerUpdate");
1111
#endif
1212

13-
// std::vector<TV> grid_force(grid_nodes, TV::Zero());
14-
15-
// #pragma omp parallel num_threads(n_threads)
16-
// {
17-
// std::vector<TV> grid_force_local(grid_nodes, TV::Zero());
18-
19-
// #pragma omp for nowait
20-
// for(int p = 0; p < Np; p++){
21-
22-
// TM Fe = particles.F[p];
23-
24-
// TM dPsidF;
25-
// if (elastic_model == ElasticModel::NeoHookean){
26-
// dPsidF = NeoHookeanPiola(Fe);
27-
// }
28-
// else if (elastic_model == ElasticModel::Hencky){ // St Venant Kirchhoff with Hencky strain
29-
// dPsidF = HenckyPiola(Fe);
30-
// }
31-
// else{
32-
// debug("You specified an unvalid ELASTIC model!");
33-
// }
34-
35-
// TM tau = dPsidF * Fe.transpose();
36-
37-
// TV xp = particles.x[p];
38-
// unsigned int i_base = std::max(0, int(std::floor((xp(0)-grid.xc)*one_over_dx)) - 1); // i_base = std::min(i_base, Nx-4); // the subtraction of one is valid for both quadratic and cubic splines
39-
// unsigned int j_base = std::max(0, int(std::floor((xp(1)-grid.yc)*one_over_dx)) - 1); // j_base = std::min(j_base, Ny-4);
40-
// #ifdef THREEDIM
41-
// unsigned int k_base = std::max(0, int(std::floor((xp(2)-grid.zc)*one_over_dx)) - 1); // k_base = std::min(k_base, Nz-4);
42-
// #endif
43-
44-
// for(int i = i_base; i < i_base+4; i++){
45-
// T xi = grid.x[i];
46-
// for(int j = j_base; j < j_base+4; j++){
47-
// T yi = grid.y[j];
48-
// #ifdef THREEDIM
49-
// for(int k = k_base; k < k_base+4; k++){
50-
// T zi = grid.z[k];
51-
// if ( grid.mass[ind(i,j,k)] > 0){
52-
// grid_force_local[ind(i,j,k)] += tau * grad_wip(xp(0), xp(1), xp(2), xi, yi, zi, one_over_dx);
53-
// } // end if non-zero grid mass
54-
// } // end for k
55-
// #else
56-
// if ( grid.mass[ind(i,j)] > 0){
57-
// grid_force_local[ind(i,j)] += tau * grad_wip(xp(0), xp(1), xi, yi, one_over_dx);
58-
// } // end if non-zero grid mass
59-
// #endif
60-
// } // end for j
61-
// } // end for i
62-
// } // end for particles
63-
64-
// #pragma omp critical
65-
// {
66-
// for (int l = 0; l<grid_nodes; l++){
67-
// grid_force[l] += grid_force_local[l];
68-
// } // end for l
69-
// } // end omp critical
70-
71-
// } // end omp parallel
72-
7313
//////////// if external grid gravity: //////////////////
7414
// std::pair<TMX, TMX> external_gravity_pair = createExternalGridGravity();
7515
////////////////////////////////////////////////////////

‎src/simulation/g2p.cpp‎

Lines changed: 11 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@ void Simulation::G2P(){
4646
T wk = N((xp(2) - zi)*one_over_dx);
4747
T wk_grad = dNdu((xp(2) - zi) * one_over_dx) * one_over_dx;
4848

49-
T weight = wi * wj * wk; //wip(xp(0), xp(1), xp(2), xi, yi, zi, one_over_dx);
49+
T weight = wi * wj * wk;
5050
TV weight_grad;
5151
weight_grad << wi_grad*wj*wk,
5252
wi*wj_grad*wk,
@@ -66,7 +66,7 @@ void Simulation::G2P(){
6666
}
6767
} // end loop k
6868
#else
69-
T weight = wi * wj; //wip(xp(0), xp(1), xi, yi, one_over_dx);
69+
T weight = wi * wj;
7070
TV weight_grad;
7171
weight_grad << wi_grad*wj,
7272
wi*wj_grad;
@@ -93,16 +93,19 @@ void Simulation::G2P(){
9393
particles.flip[p] = flipp;
9494
}
9595

96-
// Update material
97-
TM Fe_trial = particles.F[p];
98-
Fe_trial = (TM::Identity() + dt * v_grad )* Fe_trial;
99-
particles.F[p] = Fe_trial;
100-
101-
plasticity(p, plastic_count, Fe_trial);
96+
if (!use_musl){
97+
TM Fe_trial = particles.F[p];
98+
Fe_trial = Fe_trial + dt * v_grad * Fe_trial;
99+
particles.F[p] = Fe_trial;
102100

101+
plasticity(p, plastic_count, Fe_trial);
102+
}
103103

104104
} // end loop p
105105

106106
} // end omp paralell
107107

108+
if (!reduce_verbose && !use_musl)
109+
debug(" Proj: ", plastic_count, " / ", Np);
110+
108111
} // end G2P

‎src/simulation/musl.cpp‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -42,12 +42,15 @@ void Simulation::MUSL(){
4242

4343
for(int i = i_base; i < i_base+4; i++){
4444
T xi = grid.x[i];
45+
T wi = N((xp(0)-xi)*one_over_dx);
4546
for(int j = j_base; j < j_base+4; j++){
4647
T yi = grid.y[j];
48+
T wj = N((xp(1) - yi)*one_over_dx);
4749
#ifdef THREEDIM
4850
for(int k = k_base; k < k_base+4; k++){
4951
T zi = grid.z[k];
50-
T weight = wip(xp(0), xp(1), xp(2), xi, yi, zi, one_over_dx);
52+
T wk = N((xp(2) - zi)*one_over_dx);
53+
T weight = wi * wj * wk;
5154
if (weight > 1e-25){
5255
grid_v_local[ind(i,j,k)] += particles.v[p] * weight;
5356
if (flip_ratio < 0){ // APIC
@@ -60,9 +63,9 @@ void Simulation::MUSL(){
6063
}
6164
} // end for k
6265
#else
63-
T weight = wip(xp(0), xp(1), xi, yi, one_over_dx);
66+
T weight = wi * wj;
6467
if (weight > 1e-25){
65-
grid_v_local[ind(i,j)] += particles.v[p] * weight;
68+
grid_v_local[ind(i,j)] += particles.v[p] * weight;
6669
if (flip_ratio < 0){ // APIC
6770
TV posdiffvec = TV::Zero();
6871
posdiffvec(0) = xi-xp(0);

‎src/simulation/p2g.cpp‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -25,9 +25,8 @@ void Simulation::P2G(){
2525
unsigned int k_base = std::max(0, int(std::floor((xp(2)-grid.zc)*one_over_dx)) - 1); // k_base = std::min(k_base, Nz-4);
2626
#endif
2727

28-
// Stress
28+
// ALT 1: Compute stress
2929
// TM Fe = particles.F[p];
30-
3130
// TM dPsidF;
3231
// if (elastic_model == ElasticModel::NeoHookean){
3332
// dPsidF = NeoHookeanPiola(Fe);
@@ -38,9 +37,9 @@ void Simulation::P2G(){
3837
// else{
3938
// debug("You specified an unvalid ELASTIC model!");
4039
// }
41-
4240
// TM tau = dPsidF * Fe.transpose();
4341

42+
// ALT 2: Use the saved stress
4443
TM tau = particles.tau[p];
4544

4645
for(int i = i_base; i < i_base+4; i++){
@@ -57,7 +56,7 @@ void Simulation::P2G(){
5756
T wk = N((xp(2) - zi)*one_over_dx);
5857
T wk_grad = dNdu((xp(2) - zi) * one_over_dx) * one_over_dx;
5958

60-
T weight = wi * wj * wk; //wip(xp(0), xp(1), xp(2), xi, yi, zi, one_over_dx);
59+
T weight = wi * wj * wk;
6160
TV weight_grad;
6261
weight_grad << wi_grad*wj*wk,
6362
wi*wj_grad*wk,
@@ -79,7 +78,7 @@ void Simulation::P2G(){
7978
}
8079
} // end for k
8180
#else
82-
T weight = wi * wj; //wip(xp(0), xp(1), xi, yi, one_over_dx);
81+
T weight = wi * wj;
8382
TV weight_grad;
8483
weight_grad << wi_grad*wj,
8584
wi*wj_grad;

‎src/simulation/save_data.cpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,12 +15,14 @@ void Simulation::saveParticleData(std::string extra){
1515

1616
TM Fe = particles.F[p];
1717

18+
// ALT 1: Compute stress
1819
// TM tau;
1920
// if (elastic_model == ElasticModel::NeoHookean)
2021
// tau = NeoHookeanPiola(Fe) * Fe.transpose();
2122
// else if (elastic_model == ElasticModel::Hencky)
2223
// tau = HenckyPiola(Fe) * Fe.transpose();
2324

25+
// ALT 2: Use the saved stress
2426
TM tau = particles.tau[p];
2527

2628
T Je = Fe.determinant();
@@ -268,12 +270,14 @@ void Simulation::computeAvgData(TM& volavg_cauchy, TM& volavg_kirchh, T& Javg){
268270

269271
TM Fe = particles.F[p];
270272

273+
// ALT 1: Compute stress
271274
// TM tau;
272275
// if (elastic_model == ElasticModel::NeoHookean)
273276
// tau = NeoHookeanPiola(Fe) * Fe.transpose();
274277
// else if (elastic_model == ElasticModel::Hencky)
275278
// tau = HenckyPiola(Fe) * Fe.transpose();
276279

280+
// ALT 2: Use the saved stress
277281
TM tau = particles.tau[p];
278282

279283
T Je = Fe.determinant();

‎src/simulation/simulation.cpp‎

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -199,7 +199,8 @@ void Simulation::simulate(){
199199
debug("Runtime P2G = ", runtime_p2g * 1000.0, " milliseconds");
200200
debug("Runtime G2P = ", runtime_g2p * 1000.0, " milliseconds");
201201
debug("Runtime Euler = ", runtime_euler * 1000.0, " milliseconds");
202-
debug("Runtime DefGrad = ", runtime_defgrad * 1000.0, " milliseconds");
202+
if (use_musl)
203+
debug("Runtime DefGrad = ", runtime_defgrad * 1000.0, " milliseconds");
203204

204205
if (save_sim)
205206
saveTiming();
@@ -252,12 +253,13 @@ void Simulation::advanceStep(){
252253
G2P();
253254
t_g2p.stop(); runtime_g2p += t_g2p.get_timing();
254255

255-
if (use_musl == true)
256+
if (use_musl){
256257
MUSL();
257258

258-
timer t_defgrad; t_defgrad.start();
259-
// deformationUpdate();
260-
t_defgrad.stop(); runtime_defgrad += t_defgrad.get_timing();
259+
timer t_defgrad; t_defgrad.start();
260+
deformationUpdate();
261+
t_defgrad.stop(); runtime_defgrad += t_defgrad.get_timing();
262+
}
261263

262264
positionUpdate();
263265

‎src/tools.hpp‎

Lines changed: 0 additions & 70 deletions
Original file line numberDiff line numberDiff line change
@@ -216,76 +216,6 @@ inline T d2Ndu2(T u){
216216
#error Unsupported spline degree
217217
#endif
218218

219-
220-
#ifdef THREEDIM
221-
222-
inline T wip(T xp, T yp, T zp, T xi, T yi, T zi, T one_over_h){
223-
return N( (xp - xi) * one_over_h ) * N( (yp - yi) * one_over_h ) * N( (zp - zi) * one_over_h );
224-
}
225-
226-
inline TV grad_wip(T xp, T yp, T zp, T xi, T yi, T zi, T one_over_h){
227-
TV out;
228-
out << dNdu((xp - xi) * one_over_h) * N((yp - yi) * one_over_h) * N((zp - zi) * one_over_h) * one_over_h,
229-
dNdu((yp - yi) * one_over_h) * N((xp - xi) * one_over_h) * N((zp - zi) * one_over_h) * one_over_h,
230-
dNdu((zp - zi) * one_over_h) * N((xp - xi) * one_over_h) * N((yp - yi) * one_over_h) * one_over_h;
231-
return out;
232-
}
233-
234-
inline T laplace_wip(T xp, T yp, T zp, T xi, T yi, T zi, T one_over_h, T one_over_h_square){
235-
T term1 = d2Ndu2((xp - xi) * one_over_h) * N((yp - yi) * one_over_h) * N((zp - zi) * one_over_h);
236-
T term2 = d2Ndu2((yp - yi) * one_over_h) * N((xp - xi) * one_over_h) * N((zp - zi) * one_over_h);
237-
T term3 = d2Ndu2((zp - zi) * one_over_h) * N((xp - xi) * one_over_h) * N((yp - yi) * one_over_h);
238-
return ( term1 + term2 + term3 ) * one_over_h_square;
239-
}
240-
241-
#else // TWODIM
242-
243-
inline T wip(T xp, T yp, T xi, T yi, T one_over_h){
244-
return N( (xp - xi) * one_over_h ) * N( (yp - yi) * one_over_h );
245-
}
246-
247-
inline TV grad_wip(T xp, T yp, T xi, T yi, T one_over_h){
248-
TV out;
249-
out << dNdu((xp - xi) * one_over_h) * N((yp - yi) * one_over_h) * one_over_h,
250-
dNdu((yp - yi) * one_over_h) * N((xp - xi) * one_over_h) * one_over_h;
251-
return out;
252-
}
253-
254-
// inline T wip(T xp, T yp, T xi, T yi, T one_over_h){
255-
// T Lx = 0.1;
256-
// return N( fmod(xp - xi, Lx) * one_over_h ) * N( (yp - yi) * one_over_h );
257-
// }
258-
// inline TV grad_wip(T xp, T yp, T xi, T yi, T one_over_h){
259-
// T Lx = 0.1;
260-
// T xpminxi = xp-xi;
261-
// if (xpminxi <= -Lx || xpminxi >= Lx)
262-
// xpminxi = -fmod(xpminxi, Lx);
263-
// TV out;
264-
// out << dNdu( xpminxi * one_over_h) * N((yp - yi) * one_over_h) * one_over_h,
265-
// dNdu((yp - yi) * one_over_h) * N( xpminxi * one_over_h) * one_over_h;
266-
// return out;
267-
// }
268-
269-
inline T laplace_wip(T xp, T yp, T xi, T yi, T one_over_h, T one_over_h_square){
270-
T term1 = d2Ndu2((xp - xi) * one_over_h) * N((yp - yi) * one_over_h);
271-
T term2 = d2Ndu2((yp - yi) * one_over_h) * N((xp - xi) * one_over_h);
272-
return ( term1 + term2 ) * one_over_h_square;
273-
}
274-
275-
#endif
276-
277-
// Taken from: https://gist.github.com/lorenzoriano/5414671
278-
// linspace(a,b,N) return std::vector of length N in the closed interval [a,b], i.e., including b
279-
inline std::vector<T> linspace(T a, T b, size_t N) {
280-
T h = (b - a) / static_cast<T>(N-1);
281-
std::vector<T> xs(N);
282-
typename std::vector<T>::iterator x;
283-
T val;
284-
for (x = xs.begin(), val = a; x != xs.end(); ++x, val += h)
285-
*x = val;
286-
return xs;
287-
}
288-
289219
// Taken from: https://stackoverflow.com/questions/21216909/these-python-functions-in-c
290220
// Works like numpy.arange, does NOT include stop value
291221
inline std::vector<T> arange(T start, T stop, T step) {

0 commit comments

Comments
 (0)