Skip to content

Commit fd53a05

Browse files
RoushelfyRoushelfy
committed
5
1 parent 38a61b7 commit fd53a05

7 files changed

Lines changed: 158 additions & 49 deletions

File tree

‎simulators/4_friction/CMakeLists.txt‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,18 @@
1-
add_executable(4_friction_debug)
1+
add_executable(4_friction)
22

3-
target_compile_options(4_friction_debug PRIVATE -g)
4-
set_target_properties(4_friction_debug PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
3+
target_compile_options(4_friction PRIVATE -g)
4+
set_target_properties(4_friction PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
55

6-
target_link_libraries(4_friction_debug PRIVATE muda cusolver cublas cusparse )
6+
target_link_libraries(4_friction PRIVATE muda cusolver cublas cusparse )
77

88
option(BUILD_SHARED_LIBS "Build shared libraries" OFF)
99

1010
include_directories(include)
1111

12-
file(GLOB_RECURSE 4_friction_debug_CU_SOURCE CONFIGURE_DEPENDS "src/*.cu")
13-
target_sources(4_friction_debug PRIVATE ${4_friction_debug_CU_SOURCE})
12+
file(GLOB_RECURSE 4_friction_CU_SOURCE CONFIGURE_DEPENDS "src/*.cu")
13+
target_sources(4_friction PRIVATE ${4_friction_CU_SOURCE})
1414

15-
file(GLOB_RECURSE 4_friction_debug_CPP_SOURCE CONFIGURE_DEPENDS "src/*.cpp")
16-
target_sources(4_friction_debug PRIVATE ${4_friction_debug_CPP_SOURCE})
15+
file(GLOB_RECURSE 4_friction_CPP_SOURCE CONFIGURE_DEPENDS "src/*.cpp")
16+
target_sources(4_friction PRIVATE ${4_friction_CPP_SOURCE})
1717

18-
target_link_libraries(4_friction_debug PRIVATE sfml-graphics)
18+
target_link_libraries(4_friction PRIVATE sfml-graphics)

‎simulators/4_friction/src/main.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
int main()
44
{
5-
float rho = 1000, k = 4e4, initial_stretch = 1, n_seg = 8, h = 0.01, side_len = 1, tol = 0.01, mu = 0.11;
5+
float rho = 1000, k = 4e4, initial_stretch = 1, n_seg = 14, h = 0.01, side_len = 1, tol = 0.01, mu = 0.11;
66
// printf("Running mass-spring simulator with parameters: rho = %f, k = %f, initial_stretch = %f, n_seg = %d, h = %f, side_len = %f, tol = %f\n", rho, k, initial_stretch, n_seg, h, side_len, tol);
77
FrictionSimulator<float, 2> simulator(rho, side_len, initial_stretch, k, h, tol, mu, n_seg);
88
simulator.run();

‎simulators/5_mov_dirichlet/include/SpringEnergy.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ class SpringEnergy
1818
SpringEnergy &operator=(SpringEnergy &&rhs);
1919

2020
void update_x(const DeviceBuffer<T> &x);
21+
void update_DBC_target();
2122
T val(); // Calculate the value of the energy
2223
const DeviceBuffer<T> &grad(); // Calculate the gradient of the energy
2324
const DeviceTripletMatrix<T, 1> &hess(); // Calculate the Hessian matrix of the energy

‎simulators/5_mov_dirichlet/include/simulator.h‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -5,14 +5,14 @@
55
#include "square_mesh.h"
66
#include <iostream>
77
template <typename T, int dim>
8-
class FrictionSimulator
8+
class MovDirichletSimulator
99
{
1010
public:
11-
FrictionSimulator();
12-
~FrictionSimulator();
13-
FrictionSimulator(FrictionSimulator &&rhs);
14-
FrictionSimulator &operator=(FrictionSimulator &&rhs);
15-
FrictionSimulator(T rho, T side_len, T initial_stretch, T K, T h, T tol, T mu, int n_seg);
11+
MovDirichletSimulator();
12+
~MovDirichletSimulator();
13+
MovDirichletSimulator(MovDirichletSimulator &&rhs);
14+
MovDirichletSimulator &operator=(MovDirichletSimulator &&rhs);
15+
MovDirichletSimulator(T rho, T side_len, T initial_stretch, T K, T h, T tol, T mu, int n_seg);
1616
void run();
1717

1818
private:

‎simulators/5_mov_dirichlet/src/BarrierEnergy.cu‎

Lines changed: 90 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ template <typename T, int dim>
1111
struct BarrierEnergy<T, dim>::Impl
1212
{
1313
DeviceBuffer<T> device_x;
14-
DeviceBuffer<T> device_contact_area, device_n, device_o;
14+
DeviceBuffer<T> device_contact_area, device_n,device_n_ceil, device_o;
1515
int N;
1616
DeviceBuffer<T> device_grad;
1717
DeviceTripletMatrix<T, 1> device_hess;
@@ -38,9 +38,12 @@ BarrierEnergy<T, dim>::BarrierEnergy(const std::vector<T> &x, const std::vector<
3838
pimpl_->N = x.size() / dim;
3939
pimpl_->device_x.copy_from(x);
4040
pimpl_->device_contact_area.copy_from(contact_area);
41+
std::vector<T> n_ceil(dim);
42+
n_ceil[1]=-1;
43+
pimpl_->device_n_ceil.copy_from(n_ceil);
4144
pimpl_->device_n.copy_from(n);
4245
pimpl_->device_o.copy_from(o);
43-
pimpl_->device_hess.resize_triplets(pimpl_->N * dim * dim);
46+
pimpl_->device_hess.resize_triplets((pimpl_->N *2-1)* dim * dim);
4447
pimpl_->device_hess.reshape(x.size(), x.size());
4548
pimpl_->device_grad.resize(pimpl_->N * dim);
4649
}
@@ -57,20 +60,32 @@ T BarrierEnergy<T, dim>::val()
5760
auto &device_x = pimpl_->device_x;
5861
auto &device_contact_area = pimpl_->device_contact_area;
5962
auto &device_n = pimpl_->device_n;
63+
auto &device_n_ceil = pimpl_->device_n_ceil;
6064
auto &device_o = pimpl_->device_o;
6165
int N = device_x.size() / dim;
62-
DeviceBuffer<T> device_val(N);
63-
ParallelFor(256).apply(N, [device_val = device_val.viewer(), device_x = device_x.cviewer(), device_contact_area = device_contact_area.cviewer(), device_n = device_n.cviewer(), device_o = device_o.cviewer()] __device__(int i) mutable
66+
DeviceBuffer<T> device_val1(N);
67+
DeviceBuffer<T> device_val2(N);
68+
ParallelFor(256).apply(N, [device_val1 = device_val1.viewer(), device_x = device_x.cviewer(), device_contact_area = device_contact_area.cviewer(), device_n = device_n.cviewer(), device_o = device_o.cviewer()] __device__(int i) mutable
6469
{ T d = 0;
6570
for(int j=0;j<dim;j++){
6671
d += device_n(j)*(device_x(i*dim+j)-device_o(j));
6772
}
6873
if(d<dhat){
6974
T s = d / dhat;
70-
device_val(i)= kappa * device_contact_area(i) * dhat/2*(s-1)*log(s);
75+
device_val1(i)= kappa * device_contact_area(i) * dhat/2*(s-1)*log(s);
7176
} })
7277
.wait();
73-
return devicesum(device_val);
78+
ParallelFor(256).apply(N-1, [device_val2 = device_val2.viewer(), device_x = device_x.cviewer(), device_contact_area = device_contact_area.cviewer(), device_n_ceil = device_n_ceil.cviewer(), device_o = device_o.cviewer(),N] __device__(int i) mutable
79+
{ T d = 0;
80+
for(int j=0;j<dim;j++){
81+
d += device_n_ceil(j)*(device_x(i*dim+j)-device_x((N-1)*dim+j));
82+
}
83+
if(d<dhat){
84+
T s = d / dhat;
85+
device_val2(i)= kappa * device_contact_area(i) * dhat/2*(s-1)*log(s);
86+
} })
87+
.wait();
88+
return devicesum(device_val1)+devicesum(device_val2);
7489
} // Calculate the energy
7590

7691
template <typename T, int dim>
@@ -80,6 +95,7 @@ const DeviceBuffer<T> &BarrierEnergy<T, dim>::grad()
8095
auto &device_contact_area = pimpl_->device_contact_area;
8196
int N = device_x.size() / dim;
8297
auto &device_n = pimpl_->device_n;
98+
auto &device_n_ceil = pimpl_->device_n_ceil;
8399
auto &device_o = pimpl_->device_o;
84100
auto &device_grad = pimpl_->device_grad;
85101
device_grad.fill(0);
@@ -99,6 +115,23 @@ const DeviceBuffer<T> &BarrierEnergy<T, dim>::grad()
99115
}
100116
} })
101117
.wait();
118+
ParallelFor(256).apply(N-1, [device_x = device_x.cviewer(), device_contact_area = device_contact_area.cviewer(), device_grad = device_grad.viewer(), device_n_ceil = device_n_ceil.cviewer(), device_o = device_o.cviewer(),N] __device__(int i) mutable
119+
120+
{
121+
T d = 0;
122+
for(int j=0;j<dim;j++){
123+
d += device_n_ceil(j)*(device_x(i*dim+j)-device_x((N-1)*dim+j));
124+
}
125+
if (d < dhat)
126+
{
127+
T s = d / dhat;
128+
for (int j = 0; j < dim; j++)
129+
{
130+
T grad =device_contact_area(i) * dhat * (kappa / 2 * (log(s) / dhat + (s - 1) / d)) * device_n_ceil(j);
131+
device_grad(i * dim + j) += grad;
132+
device_grad((N-1) * dim + j) -= grad;
133+
}
134+
} }).wait();
102135
return device_grad;
103136
}
104137

@@ -108,6 +141,7 @@ const DeviceTripletMatrix<T, 1> &BarrierEnergy<T, dim>::hess()
108141
auto &device_x = pimpl_->device_x;
109142
auto &device_contact_area = pimpl_->device_contact_area;
110143
auto &device_n = pimpl_->device_n;
144+
auto &device_n_ceil = pimpl_->device_n_ceil;
111145
auto &device_o = pimpl_->device_o;
112146
auto &device_hess = pimpl_->device_hess;
113147
auto device_hess_row_idx = device_hess.row_indices();
@@ -144,6 +178,35 @@ const DeviceTripletMatrix<T, 1> &BarrierEnergy<T, dim>::hess()
144178
}
145179
} })
146180
.wait();
181+
ParallelFor(256).apply(N-1, [device_x = device_x.cviewer(), device_contact_area = device_contact_area.cviewer(), device_hess_val = device_hess_val.viewer(), device_hess_row_idx = device_hess_row_idx.viewer(), device_hess_col_idx = device_hess_col_idx.viewer(), N, device_n_ceil = device_n_ceil.cviewer(), device_o = device_o.cviewer()] __device__(int i) mutable
182+
{
183+
T d = 0;
184+
for (int j = 0; j < dim; j++)
185+
{
186+
d += device_n_ceil(j) * (device_x(i * dim + j) - device_x((N-1) * dim + j));
187+
}
188+
if (d < dhat)
189+
for (int j = 0; j < dim; j++)
190+
{
191+
for (int k = 0; k < dim; k++)
192+
{
193+
int idx =N*dim*dim+ i * dim * dim + j * dim + k;
194+
device_hess_row_idx(idx) = (N-1) * dim + j;
195+
device_hess_col_idx(idx) = (N-1) * dim + k;
196+
device_hess_val(idx) = device_contact_area(i) * dhat * kappa / (2 * d * d * dhat) * (d + dhat) * device_n_ceil(j) * device_n_ceil(k);
197+
}
198+
}
199+
else
200+
for (int j = 0; j < dim; j++)
201+
{
202+
for (int k = 0; k < dim; k++)
203+
{
204+
int idx = N*dim*dim+i * dim * dim + j * dim + k;
205+
device_hess_row_idx(idx) = (N-1) * dim + j;
206+
device_hess_col_idx(idx) = (N-1) * dim + k;
207+
device_hess_val(idx) = 0;
208+
}
209+
} }).wait();
147210
return device_hess;
148211

149212
} // Calculate the Hessian of the energy
@@ -153,6 +216,7 @@ T BarrierEnergy<T, dim>::init_step_size(const DeviceBuffer<T> &p)
153216
{
154217
auto &device_x = pimpl_->device_x;
155218
auto &device_n = pimpl_->device_n;
219+
auto &device_n_ceil = pimpl_->device_n_ceil;
156220
auto &device_o = pimpl_->device_o;
157221
int N = device_x.size() / dim;
158222
DeviceBuffer<T> device_alpha(N);
@@ -176,6 +240,26 @@ T BarrierEnergy<T, dim>::init_step_size(const DeviceBuffer<T> &p)
176240
device_alpha(i) = min(device_alpha(i), 0.9 * alpha / -p_n);
177241
} })
178242
.wait();
243+
244+
ParallelFor(256)
245+
.apply(N-1, [device_x = device_x.cviewer(), p = p.cviewer(), device_alpha = device_alpha.viewer(), device_n_ceil = device_n_ceil.cviewer(), device_o = device_o.cviewer(),N] __device__(int i) mutable
246+
247+
{
248+
T p_n = 0;
249+
for (int j = 0; j < dim; j++)
250+
{
251+
p_n += p(i * dim + j) * device_n_ceil(j);
252+
}
253+
if (p_n < 0)
254+
{
255+
T alpha = 0;
256+
for (int j = 0; j < dim; j++)
257+
{
258+
alpha += device_n_ceil(j) * (device_x(i * dim + j) - device_x((N-1) * dim + j));
259+
}
260+
device_alpha(i) = min(device_alpha(i), 0.9 * alpha / -p_n);
261+
} })
262+
.wait();
179263
return min_vector(device_alpha);
180264
}
181265
template class BarrierEnergy<float, 2>;

‎simulators/5_mov_dirichlet/src/SpringEnergy.cu‎

Lines changed: 20 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ struct SpringEnergy<T, dim>::Impl
1515
DeviceBuffer<Eigen::Matrix<T, dim, 1>> device_DBC_target;
1616
DeviceBuffer<T> device_grad;
1717
DeviceTripletMatrix<T, 1> device_hess;
18-
T k;
18+
T k,h;
1919
int N;
2020
};
2121

@@ -36,7 +36,7 @@ SpringEnergy<T, dim>::SpringEnergy(const SpringEnergy<T, dim> &rhs)
3636
: pimpl_{std::make_unique<Impl>(*rhs.pimpl_)} {}
3737

3838
template <typename T, int dim>
39-
SpringEnergy<T, dim>::SpringEnergy(const std::vector<T> &x, const std::vector<T> &m, const std::vector<int> &DBC, const std::vector<Eigen::Matrix<T, dim, 1>> &DBC_target, T k)
39+
SpringEnergy<T, dim>::SpringEnergy(const std::vector<T> &x, const std::vector<T> &m, const std::vector<int> &DBC, const std::vector<T> &DBC_target, T k,T h)
4040
: pimpl_{std::make_unique<Impl>()}
4141
{
4242
pimpl_->N = x.size() / dim;
@@ -45,6 +45,7 @@ SpringEnergy<T, dim>::SpringEnergy(const std::vector<T> &x, const std::vector<T>
4545
pimpl_->device_DBC.copy_from(DBC);
4646
pimpl_->device_DBC_target.copy_from(DBC_target);
4747
pimpl_->k = k;
48+
pimpl_->h = h;
4849
pimpl_->device_grad.resize(pimpl_->N * dim);
4950
pimpl_->device_hess.resize_triplets(pimpl_->N * dim * dim);
5051
pimpl_->device_hess.reshape(x.size(), x.size());
@@ -56,6 +57,21 @@ void SpringEnergy<T, dim>::update_x(const DeviceBuffer<T> &x)
5657
pimpl_->device_x.view().copy_from(x);
5758
}
5859

60+
template <typename T, int dim>
61+
void SpringEnergy<T, dim>::update_DBC_target()
62+
{
63+
// for i in range(0, len(DBC)):
64+
// if (DBC_limit[i] - x_n[DBC[i]]).dot(DBC_v[i]) > 0:
65+
// DBC_target.append(x_n[DBC[i]] + h * DBC_v[i])
66+
// else:
67+
// DBC_target.append(x_n[DBC[i]])
68+
auto &device_x = pimpl_->device_x;
69+
auto &device_DBC = pimpl_->device_DBC;
70+
auto &device_DBC_target = pimpl_->device_DBC_target;
71+
auto h = pimpl_->h;
72+
73+
}
74+
5975
template <typename T, int dim>
6076
T SpringEnergy<T, dim>::val()
6177
{
@@ -75,7 +91,7 @@ T SpringEnergy<T, dim>::val()
7591
Eigen::Matrix<T, dim, 1> diff;
7692
for (int j = 0; j < dim; ++j)
7793
{
78-
diff(j) = device_x(idx * dim + j) - device_DBC_target(i)(j);
94+
diff(j) = device_x(idx * dim + j) - device_DBC_target(i*dim + j);
7995
}
8096
device_val(i) = 0.5 * k * device_m(idx) * diff.dot(diff); })
8197
.wait();
@@ -102,7 +118,7 @@ const DeviceBuffer<T> &SpringEnergy<T, dim>::grad()
102118
Eigen::Matrix<T, dim, 1> grad;
103119
for (int j = 0; j < dim; ++j)
104120
{
105-
grad(j) = device_x(idx * dim + j) - device_DBC_target(i)(j);
121+
grad(j) = device_x(idx * dim + j) - device_DBC_target(i*dim + j);
106122
}
107123
grad *= k * device_m(idx);
108124
for (int j = 0; j < dim; ++j)

0 commit comments

Comments
 (0)