Loading [MathJax]/jax/output/HTML-CSS/config.js
PolyFEM
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ViscousDamping.cpp
Go to the documentation of this file.
1#include "ViscousDamping.hpp"
2
3namespace polyfem::assembler
4{
5 void ViscousDamping::compute_stress_aux(const Eigen::MatrixXd &F, const Eigen::MatrixXd &dFdt, Eigen::MatrixXd &dRdF, Eigen::MatrixXd &dRdFdot) const
6 {
7 const int size = F.rows();
8
9 Eigen::MatrixXd dEdt = 0.5 * (dFdt.transpose() * F + F.transpose() * dFdt);
10
11 Eigen::MatrixXd tmp = 2 * damping_params_[0] * dEdt + damping_params_[1] * dEdt.trace() * Eigen::MatrixXd::Identity(size, size);
12 dRdF = dFdt * tmp;
13 dRdFdot = F * tmp; // Fdot is dFdt
14 }
15
16 void ViscousDamping::compute_stress_grad_aux(const Eigen::MatrixXd &F, const Eigen::MatrixXd &dFdt, Eigen::MatrixXd &d2RdF2, Eigen::MatrixXd &d2RdFdFdot, Eigen::MatrixXd &d2RdFdot2) const
17 {
18 const int size = F.rows();
19
20 Eigen::MatrixXd dEdt = 0.5 * (dFdt.transpose() * F + F.transpose() * dFdt);
21
22 Eigen::MatrixXd dEdotdF, dEdotdFdot;
23 dEdotdF.setZero(size * size, size * size);
24 dEdotdFdot.setZero(size * size, size * size);
25 for (size_t i = 0; i < size; ++i)
26 {
27 for (size_t j = 0; j < size; ++j)
28 {
29 for (size_t p = 0; p < size; ++p)
30 {
31 dEdotdF(i * size + j, p * size + j) += dFdt(p, i) / 2.;
32 dEdotdFdot(i * size + j, p * size + j) += F(p, i) / 2.;
33 dEdotdF(i * size + j, p * size + i) += dFdt(p, j) / 2.;
34 dEdotdFdot(i * size + j, p * size + i) += F(p, j) / 2.;
35 }
36 }
37 }
38
39 d2RdF2.setZero(size * size, size * size);
40 d2RdFdFdot.setZero(size * size, size * size);
41 d2RdFdot2.setZero(size * size, size * size);
42 for (size_t i = 0; i < size; ++i)
43 {
44 for (size_t j = 0; j < size; ++j)
45 {
46 const int idx = i * size + j;
47 for (size_t k = 0; k < size; ++k)
48 {
49 d2RdF2.row(idx) += (2 * damping_params_[0]) * dFdt(i, k) * dEdotdF.row(k * size + j) + damping_params_[1] * dFdt(i, j) * dEdotdF.row(k * size + k);
50
51 d2RdFdot2.row(idx) += (2 * damping_params_[0]) * F(i, k) * dEdotdFdot.row(k * size + j) + damping_params_[1] * F(i, j) * dEdotdFdot.row(k * size + k);
52
53 d2RdFdFdot(idx, i * size + k) += 2 * damping_params_[0] * dEdt(k, j);
54 d2RdFdFdot.row(idx) += (2 * damping_params_[0]) * (dFdt(i, k) * dEdotdFdot.row(k * size + j)) + damping_params_[1] * (dEdotdFdot.row(k * size + k) * dFdt(i, j));
55 d2RdFdFdot(idx, idx) += damping_params_[1] * dEdt(k, k);
56 }
57 }
58 }
59 }
60
61 void ViscousDamping::add_multimaterial(const int index, const json &params, const Units &units)
62 {
63 // TODO add units
64 assert(size() == 2 || size() == 3);
65
66 if (params.contains("psi"))
67 damping_params_[0] = params["psi"];
68 if (params.contains("phi"))
69 damping_params_[1] = params["phi"];
70 }
71
72 // E := 0.5(F^T F - I), Compute \int F * (2\psi dE/dt + \phi Tr(dE/dt) I) : gradv du
73 Eigen::VectorXd
74 ViscousDampingPrev::assemble_gradient(const NonLinearAssemblerData &data) const
75 {
76 if (data.x_prev.size() != data.x.size())
77 return Eigen::VectorXd::Zero(data.vals.basis_values.size() * size());
78
79 Eigen::MatrixXd local_disp = Eigen::MatrixXd::Zero(data.vals.basis_values.size(), size());
80 Eigen::MatrixXd local_prev_disp = Eigen::MatrixXd::Zero(data.vals.basis_values.size(), size());
81 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
82 {
83 const auto &bs = data.vals.basis_values[i];
84 for (size_t ii = 0; ii < bs.global.size(); ++ii)
85 {
86 local_disp.row(i) += bs.global[ii].val * data.x.block(bs.global[ii].index * size(), 0, size(), 1).transpose();
87 local_prev_disp.row(i) += bs.global[ii].val * data.x_prev.block(bs.global[ii].index * size(), 0, size(), 1).transpose();
88 }
89 }
90
91 Eigen::MatrixXd local_vel = (local_disp - local_prev_disp) / data.dt;
92
93 Eigen::MatrixXd G;
94 G.setZero(data.vals.basis_values.size(), size());
95
96 const int n_pts = data.da.size();
97
98 Eigen::MatrixXd def_grad, vel_def_grad;
99 Eigen::MatrixXd delF_delU;
100 Eigen::MatrixXd dRdF, dRdFdot;
101 for (long p = 0; p < n_pts; ++p)
102 {
103 Eigen::MatrixXd grad(data.vals.basis_values.size(), size());
104
105 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
106 {
107 grad.row(i) = data.vals.basis_values[i].grad.row(p);
108 }
109
110 delF_delU = grad * data.vals.jac_it[p];
111 def_grad = local_disp.transpose() * delF_delU + Eigen::MatrixXd::Identity(size(), size());
112 vel_def_grad = local_vel.transpose() * delF_delU;
113
114 compute_stress_aux(def_grad, vel_def_grad, dRdF, dRdFdot);
115 G -= delF_delU * dRdFdot * (data.da(p) / data.dt);
116 }
117
118 G.transposeInPlace();
119
120 Eigen::VectorXd temp(Eigen::Map<Eigen::VectorXd>(G.data(), G.size()));
121
122 return temp;
123 }
124
125 // E := 0.5(F^T F - I), Compute Stress = \int dF/dt * (2\psi dE/dt + \phi Tr(dE/dt) I) : gradv du
126 Eigen::VectorXd
127 ViscousDamping::assemble_gradient(const NonLinearAssemblerData &data) const
128 {
129 if (data.x_prev.size() != data.x.size())
130 return Eigen::VectorXd::Zero(data.vals.basis_values.size() * size());
131
132 Eigen::MatrixXd local_disp = Eigen::MatrixXd::Zero(data.vals.basis_values.size(), size());
133 Eigen::MatrixXd local_prev_disp = Eigen::MatrixXd::Zero(data.vals.basis_values.size(), size());
134 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
135 {
136 const auto &bs = data.vals.basis_values[i];
137 for (size_t ii = 0; ii < bs.global.size(); ++ii)
138 {
139 local_disp.row(i) += bs.global[ii].val * data.x.block(bs.global[ii].index * size(), 0, size(), 1).transpose();
140 local_prev_disp.row(i) += bs.global[ii].val * data.x_prev.block(bs.global[ii].index * size(), 0, size(), 1).transpose();
141 }
142 }
143
144 Eigen::MatrixXd local_vel = (local_disp - local_prev_disp) / data.dt;
145
146 Eigen::MatrixXd G;
147 G.setZero(data.vals.basis_values.size(), size());
148
149 const int n_pts = data.da.size();
150
151 Eigen::MatrixXd def_grad, dFdt, dEdt;
152 for (long p = 0; p < n_pts; ++p)
153 {
154 Eigen::MatrixXd grad(data.vals.basis_values.size(), size());
155
156 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
157 {
158 grad.row(i) = data.vals.basis_values[i].grad.row(p);
159 }
160
161 Eigen::MatrixXd delF_delU = grad * data.vals.jac_it[p];
162
163 def_grad = local_disp.transpose() * delF_delU + Eigen::MatrixXd::Identity(size(), size());
164 dFdt = local_vel.transpose() * delF_delU;
165
166 dEdt = dFdt.transpose() * def_grad;
167 dEdt = (dEdt + dEdt.transpose()).eval() / 2.;
168
169 Eigen::MatrixXd tmp = (2 * damping_params_[0]) * dEdt + (damping_params_[1] * dEdt.trace()) * Eigen::MatrixXd::Identity(size(), size());
170 G += delF_delU * ((dFdt + def_grad / data.dt) * tmp).transpose() * data.da(p);
171 }
172
173 G.transposeInPlace();
174
175 Eigen::VectorXd temp(Eigen::Map<Eigen::VectorXd>(G.data(), G.size()));
176
177 return temp;
178 }
179
180 // Compute \int grad phi_i : dStress / dF^n : grad phi_j dx
181 Eigen::MatrixXd
182 ViscousDamping::assemble_hessian(const NonLinearAssemblerData &data) const
183 {
184 Eigen::MatrixXd hessian;
185 hessian.setZero(data.vals.basis_values.size() * size(), data.vals.basis_values.size() * size());
186 if (data.x_prev.size() != data.x.size())
187 return hessian;
188
189 Eigen::MatrixXd local_disp = Eigen::MatrixXd::Zero(data.vals.basis_values.size(), size());
190 Eigen::MatrixXd local_prev_disp = Eigen::MatrixXd::Zero(data.vals.basis_values.size(), size());
191 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
192 {
193 const auto &bs = data.vals.basis_values[i];
194 for (size_t ii = 0; ii < bs.global.size(); ++ii)
195 {
196 local_disp.row(i) += bs.global[ii].val * data.x.block(bs.global[ii].index * size(), 0, size(), 1).transpose();
197 local_prev_disp.row(i) += bs.global[ii].val * data.x_prev.block(bs.global[ii].index * size(), 0, size(), 1).transpose();
198 }
199 }
200
201 Eigen::MatrixXd local_vel = (local_disp - local_prev_disp) / data.dt;
202
203 const int n_pts = data.da.size();
204
205 Eigen::MatrixXd def_grad, dFdt;
206 Eigen::MatrixXd d2RdF2, d2RdFdFdot, d2RdFdot2;
207 Eigen::MatrixXd hessian_temp, delF_delU_tensor;
208 for (long p = 0; p < n_pts; ++p)
209 {
210 Eigen::MatrixXd grad(data.vals.basis_values.size(), size());
211
212 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
213 {
214 grad.row(i) = data.vals.basis_values[i].grad.row(p);
215 }
216
217 Eigen::MatrixXd delF_delU = grad * data.vals.jac_it[p];
218
219 def_grad = local_disp.transpose() * delF_delU + Eigen::MatrixXd::Identity(size(), size());
220 dFdt = local_vel.transpose() * delF_delU;
221
222 compute_stress_grad_aux(def_grad, dFdt, d2RdF2, d2RdFdFdot, d2RdFdot2);
223 hessian_temp = d2RdF2 + (1. / data.dt) * (d2RdFdFdot + d2RdFdFdot.transpose()) + (1. / data.dt / data.dt) * d2RdFdot2;
224
225 delF_delU_tensor = Eigen::MatrixXd::Zero(size() * size(), grad.size());
226 for (size_t i = 0; i < local_disp.rows(); ++i)
227 for (size_t j = 0; j < size(); ++j)
228 for (size_t d = 0; d < size(); d++)
229 delF_delU_tensor(size() * j + d, i * size() + j) = delF_delU(i, d);
230
231 hessian += delF_delU_tensor.transpose() * hessian_temp * delF_delU_tensor * data.da(p);
232 }
233
234 return hessian;
235 }
236
237 // Compute \int grad phi_i : dStress / dF^{n-1} : grad phi_j dx
238 Eigen::MatrixXd
239 ViscousDampingPrev::assemble_hessian(const NonLinearAssemblerData &data) const
240 {
241 Eigen::MatrixXd stress_grad_Ut;
242 stress_grad_Ut.setZero(data.vals.basis_values.size() * size(), data.vals.basis_values.size() * size());
243 if (data.x_prev.size() != data.x.size())
244 return stress_grad_Ut;
245
246 Eigen::MatrixXd local_disp = Eigen::MatrixXd::Zero(data.vals.basis_values.size(), size());
247 Eigen::MatrixXd local_prev_disp = Eigen::MatrixXd::Zero(data.vals.basis_values.size(), size());
248 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
249 {
250 const auto &bs = data.vals.basis_values[i];
251 for (size_t ii = 0; ii < bs.global.size(); ++ii)
252 {
253 local_disp.row(i) += bs.global[ii].val * data.x.block(bs.global[ii].index * size(), 0, size(), 1).transpose();
254 local_prev_disp.row(i) += bs.global[ii].val * data.x_prev.block(bs.global[ii].index * size(), 0, size(), 1).transpose();
255 }
256 }
257
258 Eigen::MatrixXd local_vel = (local_disp - local_prev_disp) / data.dt;
259
260 const int n_pts = data.da.size();
261
262 Eigen::MatrixXd def_grad, dFdt;
263 Eigen::MatrixXd d2RdF2, d2RdFdFdot, d2RdFdot2;
264 Eigen::MatrixXd stress_grad_Ut_temp, delF_delU_tensor;
265 for (long p = 0; p < n_pts; ++p)
266 {
267 Eigen::MatrixXd grad(data.vals.basis_values.size(), size());
268
269 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
270 {
271 grad.row(i) = data.vals.basis_values[i].grad.row(p);
272 }
273
274 Eigen::MatrixXd delF_delU = grad * data.vals.jac_it[p];
275
276 def_grad = local_disp.transpose() * delF_delU + Eigen::MatrixXd::Identity(size(), size());
277 dFdt = local_vel.transpose() * delF_delU;
278
279 compute_stress_grad_aux(def_grad, dFdt, d2RdF2, d2RdFdFdot, d2RdFdot2);
280 stress_grad_Ut_temp = -(1. / data.dt) * d2RdFdFdot - (1. / data.dt / data.dt) * d2RdFdot2;
281
282 delF_delU_tensor = Eigen::MatrixXd::Zero(size() * size(), grad.size());
283 for (size_t i = 0; i < local_disp.rows(); ++i)
284 for (size_t j = 0; j < size(); ++j)
285 for (size_t d = 0; d < size(); d++)
286 delF_delU_tensor(size() * j + d, i * size() + j) = delF_delU(i, d);
287
288 stress_grad_Ut += delF_delU_tensor.transpose() * stress_grad_Ut_temp * delF_delU_tensor * data.da(p);
289 }
290
291 return stress_grad_Ut;
292 }
293
294 // E := 0.5(F^T F - I), Compute Energy = \int \psi \| \frac{\partial E}{\partial t} \|^2 + 0.5 \phi (Tr(\frac{\partial E}{\partial t}))^2 du
295 double ViscousDamping::compute_energy(const NonLinearAssemblerData &data) const
296 {
297 if (data.x_prev.size() != data.x.size())
298 return 0;
299
300 Eigen::MatrixXd local_disp = Eigen::MatrixXd::Zero(data.vals.basis_values.size(), size());
301 Eigen::MatrixXd local_prev_disp = Eigen::MatrixXd::Zero(data.vals.basis_values.size(), size());
302 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
303 {
304 const auto &bs = data.vals.basis_values[i];
305 for (size_t ii = 0; ii < bs.global.size(); ++ii)
306 {
307 local_disp.row(i) += bs.global[ii].val * data.x.block(bs.global[ii].index * size(), 0, size(), 1).transpose();
308 local_prev_disp.row(i) += bs.global[ii].val * data.x_prev.block(bs.global[ii].index * size(), 0, size(), 1).transpose();
309 }
310 }
311
312 Eigen::MatrixXd local_vel = (local_disp - local_prev_disp) / data.dt;
313
314 double energy = 0;
315 const int n_pts = data.da.size();
316
317 Eigen::MatrixXd def_grad, dFdt, dEdt;
318 for (long p = 0; p < n_pts; ++p)
319 {
320 Eigen::MatrixXd grad(data.vals.basis_values.size(), size());
321
322 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
323 {
324 grad.row(i) = data.vals.basis_values[i].grad.row(p);
325 }
326
327 const Eigen::MatrixXd delF_delU = grad * data.vals.jac_it[p];
328 def_grad = local_disp.transpose() * delF_delU + Eigen::MatrixXd::Identity(size(), size());
329 dFdt = local_vel.transpose() * delF_delU;
330
331 dEdt = dFdt.transpose() * def_grad;
332 dEdt = (dEdt + dEdt.transpose()).eval() / 2.;
333
334 double val = damping_params_[0] * dEdt.squaredNorm() + 0.5 * damping_params_[1] * pow(dEdt.trace(), 2);
335 energy += val * data.da(p);
336 }
337
338 return energy;
339 }
340
341 void ViscousDamping::compute_stress_grad(
342 const OptAssemblerData &data,
343 const Eigen::MatrixXd &prev_grad_u_i,
344 Eigen::MatrixXd &stress,
345 Eigen::MatrixXd &result) const
346 {
347 const double dt = data.dt;
348 const Eigen::MatrixXd &grad_u_i = data.grad_u_i;
349
350 const Eigen::MatrixXd F = Eigen::MatrixXd::Identity(size(), size()) + grad_u_i;
351 const Eigen::MatrixXd Fdot = (grad_u_i - prev_grad_u_i) / dt;
352 Eigen::MatrixXd d2RdF2, d2RdFdFdot, d2RdFdot2;
353
354 Eigen::MatrixXd dEdt = Fdot.transpose() * F;
355 dEdt = (dEdt + dEdt.transpose()).eval() / 2.;
356
357 Eigen::MatrixXd tmp = (2 * damping_params_[0]) * dEdt + (damping_params_[1] * dEdt.trace()) * Eigen::MatrixXd::Identity(size(), size());
358 stress = (Fdot + (1. / dt) * F) * tmp;
359
360 compute_stress_grad_aux(F, Fdot, d2RdF2, d2RdFdFdot, d2RdFdot2);
361 result = d2RdF2 + (1. / dt) * (d2RdFdFdot + d2RdFdFdot.transpose()) + (1. / dt / dt) * d2RdFdot2;
362 }
363
364 void ViscousDamping::compute_stress_prev_grad(
365 const OptAssemblerData &data,
366 const Eigen::MatrixXd &prev_grad_u_i,
367 Eigen::MatrixXd &result) const
368 {
369 const double dt = data.dt;
370 const Eigen::MatrixXd &grad_u_i = data.grad_u_i;
371
372 const Eigen::MatrixXd def_grad = Eigen::MatrixXd::Identity(grad_u_i.rows(), grad_u_i.cols()) + grad_u_i;
373 const Eigen::MatrixXd def_grad_fd = (grad_u_i - prev_grad_u_i) / dt;
374 Eigen::MatrixXd d2RdF2, d2RdFdFdot, d2RdFdot2;
375
376 compute_stress_grad_aux(def_grad, def_grad_fd, d2RdF2, d2RdFdFdot, d2RdFdot2);
377 result = -(1. / dt) * d2RdFdFdot - (1. / dt / dt) * d2RdFdot2;
378 }
379
380 void ViscousDamping::compute_dstress_dpsi_dphi(
381 const OptAssemblerData &data,
382 const Eigen::MatrixXd &prev_grad_u_i,
383 Eigen::MatrixXd &dstress_dpsi,
384 Eigen::MatrixXd &dstress_dphi)
385 {
386 const double dt = data.dt;
387 const Eigen::MatrixXd &grad_u_i = data.grad_u_i;
388
389 const Eigen::MatrixXd F = Eigen::MatrixXd::Identity(grad_u_i.rows(), grad_u_i.cols()) + grad_u_i;
390 const Eigen::MatrixXd dFdt = (grad_u_i - prev_grad_u_i) / dt;
391 const Eigen::MatrixXd dEdt = 0.5 * (dFdt.transpose() * F + F.transpose() * dFdt);
392
393 dstress_dpsi = 2 * (dFdt + F / dt) * dEdt;
394 dstress_dphi = dEdt.trace() * (dFdt + F / dt);
395 }
396} // namespace polyfem::assembler
val
double val
Definition Assembler.cpp:86
polyfem::assembler::ElementAssemblyValues::jac_it
std::vector< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > > jac_it
Definition ElementAssemblyValues.hpp:22
polyfem::assembler::ViscousDamping::assemble_gradient
virtual Eigen::VectorXd assemble_gradient(const NonLinearAssemblerData &data) const override
Definition ViscousDamping.cpp:127
polyfem::assembler::ViscousDamping::compute_stress_aux
void compute_stress_aux(const Eigen::MatrixXd &F, const Eigen::MatrixXd &dFdt, Eigen::MatrixXd &dRdF, Eigen::MatrixXd &dRdFdot) const
Definition ViscousDamping.cpp:5
polyfem::assembler::ViscousDamping::compute_energy
double compute_energy(const NonLinearAssemblerData &data) const override
Definition ViscousDamping.cpp:295
polyfem::assembler::ViscousDamping::compute_dstress_dpsi_dphi
static void compute_dstress_dpsi_dphi(const OptAssemblerData &data, const Eigen::MatrixXd &prev_grad_u_i, Eigen::MatrixXd &dstress_dpsi, Eigen::MatrixXd &dstress_dphi)
Definition ViscousDamping.cpp:380
polyfem::assembler::ViscousDamping::compute_stress_grad_aux
void compute_stress_grad_aux(const Eigen::MatrixXd &F, const Eigen::MatrixXd &dFdt, Eigen::MatrixXd &d2RdF2, Eigen::MatrixXd &d2RdFdFdot, Eigen::MatrixXd &d2RdFdot2) const
Definition ViscousDamping.cpp:16
polyfem::assembler::ViscousDamping::assemble_hessian
virtual Eigen::MatrixXd assemble_hessian(const NonLinearAssemblerData &data) const override
Definition ViscousDamping.cpp:182
polyfem::assembler::ViscousDamping::add_multimaterial
void add_multimaterial(const int index, const json &params, const Units &units) override
Definition ViscousDamping.cpp:61
polyfem::assembler::ViscousDamping::compute_stress_grad
void compute_stress_grad(const OptAssemblerData &data, const Eigen::MatrixXd &prev_grad_u_i, Eigen::MatrixXd &stress, Eigen::MatrixXd &result) const override
Definition ViscousDamping.cpp:341
polyfem::assembler::ViscousDamping::compute_stress_prev_grad
void compute_stress_prev_grad(const OptAssemblerData &data, const Eigen::MatrixXd &prev_grad_u_i, Eigen::MatrixXd &result) const override
Definition ViscousDamping.cpp:364
polyfem::assembler::ViscousDampingPrev::assemble_gradient
Eigen::VectorXd assemble_gradient(const NonLinearAssemblerData &data) const override
Definition ViscousDamping.cpp:74
polyfem::assembler::ViscousDampingPrev::assemble_hessian
Eigen::MatrixXd assemble_hessian(const NonLinearAssemblerData &data) const override
Definition ViscousDamping.cpp:239
polyfem::assembler
Used for test only.
Definition AMIPSEnergy.cpp:6
polyfem::json
nlohmann::json json
Definition Common.hpp:9