Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FixedCorotational.cpp
Go to the documentation of this file.
2
5
6namespace polyfem::assembler
7{
8 namespace
9 {
10 bool delta(int i, int j)
11 {
12 return (i == j) ? true : false;
13 }
14
15 template <int dim>
16 Eigen::Matrix<double, dim, dim> hat(const Eigen::Matrix<double, dim, 1> &x)
17 {
18
19 Eigen::Matrix<double, dim, dim> prod;
20 prod.setZero();
21
22 prod(0, 1) = -x(2);
23 prod(0, 2) = x(1);
24 prod(1, 0) = x(2);
25 prod(1, 2) = -x(0);
26 prod(2, 0) = -x(1);
27 prod(2, 1) = x(0);
28
29 return prod;
30 }
31
32 template <int dim>
33 Eigen::Matrix<double, dim, 1> cross(const Eigen::Matrix<double, dim, 1> &x, const Eigen::Matrix<double, dim, 1> &y)
34 {
35
36 Eigen::Matrix<double, dim, 1> z;
37 z.setZero();
38
39 z(0) = x(1) * y(2) - x(2) * y(1);
40 z(1) = x(2) * y(0) - x(0) * y(2);
41 z(2) = x(0) * y(1) - x(1) * y(0);
42
43 return z;
44 }
45
46 template <int dim>
47 Eigen::Matrix<double, dim, dim> computeCofactorMtr(Eigen::Matrix<double, dim, dim>& F)
48 {
49 Eigen::Matrix<double, dim, dim> A;
50 A.setZero();
51 if constexpr (dim == 2)
52 {
53 A(0, 0) = F(1, 1);
54 A(0, 1) = -F(1, 0);
55 A(1, 0) = -F(0, 1);
56 A(1, 1) = F(0, 0);
57 }
58 else if constexpr (dim == 3)
59 {
60 A(0, 0) = F(1, 1) * F(2, 2) - F(1, 2) * F(2, 1);
61 A(0, 1) = F(1, 2) * F(2, 0) - F(1, 0) * F(2, 2);
62 A(0, 2) = F(1, 0) * F(2, 1) - F(1, 1) * F(2, 0);
63 A(1, 0) = F(0, 2) * F(2, 1) - F(0, 1) * F(2, 2);
64 A(1, 1) = F(0, 0) * F(2, 2) - F(0, 2) * F(2, 0);
65 A(1, 2) = F(0, 1) * F(2, 0) - F(0, 0) * F(2, 1);
66 A(2, 0) = F(0, 1) * F(1, 2) - F(0, 2) * F(1, 1);
67 A(2, 1) = F(0, 2) * F(1, 0) - F(0, 0) * F(1, 2);
68 A(2, 2) = F(0, 0) * F(1, 1) - F(0, 1) * F(1, 0);
69 }
70 return A;
71 }
72
73 } // namespace
74
78
79 void FixedCorotational::add_multimaterial(const int index, const json &params, const Units &units)
80 {
81 assert(size() == 2 || size() == 3);
82
83 params_.add_multimaterial(index, params, size() == 3, units.stress());
84 }
85
86 Eigen::VectorXd
88 {
89 Eigen::Matrix<double, Eigen::Dynamic, 1> gradient;
90
91 if (size() == 2)
92 {
93 switch (data.vals.basis_values.size())
94 {
95 case 3:
96 {
97 gradient.resize(6);
98 compute_energy_aux_gradient_fast<3, 2>(data, gradient);
99 break;
100 }
101 case 6:
102 {
103 gradient.resize(12);
104 compute_energy_aux_gradient_fast<6, 2>(data, gradient);
105 break;
106 }
107 case 10:
108 {
109 gradient.resize(20);
110 compute_energy_aux_gradient_fast<10, 2>(data, gradient);
111 break;
112 }
113 default:
114 {
115 gradient.resize(data.vals.basis_values.size() * 2);
116 compute_energy_aux_gradient_fast<Eigen::Dynamic, 2>(data, gradient);
117 break;
118 }
119 }
120 }
121 else // if (size() == 3)
122 {
123 assert(size() == 3);
124 switch (data.vals.basis_values.size())
125 {
126 case 4:
127 {
128 gradient.resize(12);
129 compute_energy_aux_gradient_fast<4, 3>(data, gradient);
130 break;
131 }
132 case 10:
133 {
134 gradient.resize(30);
135 compute_energy_aux_gradient_fast<10, 3>(data, gradient);
136 break;
137 }
138 case 20:
139 {
140 gradient.resize(60);
141 compute_energy_aux_gradient_fast<20, 3>(data, gradient);
142 break;
143 }
144 default:
145 {
146 gradient.resize(data.vals.basis_values.size() * 3);
147 compute_energy_aux_gradient_fast<Eigen::Dynamic, 3>(data, gradient);
148 break;
149 }
150 }
151 }
152
153 return gradient;
154 }
155
156 Eigen::MatrixXd
158 {
159 Eigen::MatrixXd hessian;
160
161 if (size() == 2)
162 {
163 switch (data.vals.basis_values.size())
164 {
165 case 3:
166 {
167 hessian.resize(6, 6);
168 hessian.setZero();
169 compute_energy_hessian_aux_fast<3, 2>(data, hessian);
170 break;
171 }
172 case 6:
173 {
174 hessian.resize(12, 12);
175 hessian.setZero();
176 compute_energy_hessian_aux_fast<6, 2>(data, hessian);
177 break;
178 }
179 case 10:
180 {
181 hessian.resize(20, 20);
182 hessian.setZero();
183 compute_energy_hessian_aux_fast<10, 2>(data, hessian);
184 break;
185 }
186 default:
187 {
188 hessian.resize(data.vals.basis_values.size() * 2, data.vals.basis_values.size() * 2);
189 hessian.setZero();
190 compute_energy_hessian_aux_fast<Eigen::Dynamic, 2>(data, hessian);
191 break;
192 }
193 }
194 }
195 else // if (size() == 3)
196 {
197 assert(size() == 3);
198 switch (data.vals.basis_values.size())
199 {
200 case 4:
201 {
202 hessian.resize(12, 12);
203 hessian.setZero();
204 compute_energy_hessian_aux_fast<4, 3>(data, hessian);
205 break;
206 }
207 case 10:
208 {
209 hessian.resize(30, 30);
210 hessian.setZero();
211 compute_energy_hessian_aux_fast<10, 3>(data, hessian);
212 break;
213 }
214 case 20:
215 {
216 hessian.resize(60, 60);
217 hessian.setZero();
218 compute_energy_hessian_aux_fast<20, 3>(data, hessian);
219 break;
220 }
221 default:
222 {
223 hessian.resize(data.vals.basis_values.size() * 3, data.vals.basis_values.size() * 3);
224 hessian.setZero();
225 compute_energy_hessian_aux_fast<Eigen::Dynamic, 3>(data, hessian);
226 break;
227 }
228 }
229 }
230
231 return hessian;
232 }
233
235 const int all_size,
236 const ElasticityTensorType &type,
237 Eigen::MatrixXd &all,
238 const std::function<Eigen::MatrixXd(const Eigen::MatrixXd &)> &fun) const
239 {
240 const auto &displacement = data.fun;
241 const auto &local_pts = data.local_pts;
242 const auto &bs = data.bs;
243 const auto &gbs = data.gbs;
244 const auto el_id = data.el_id;
245 const auto t = data.t;
246
247 Eigen::MatrixXd displacement_grad(size(), size());
248
249 assert(displacement.cols() == 1);
250
251 all.resize(local_pts.rows(), all_size);
252
254 vals.compute(el_id, size() == 3, local_pts, bs, gbs);
255 const auto I = Eigen::MatrixXd::Identity(size(), size());
256
257 for (long p = 0; p < local_pts.rows(); ++p)
258 {
259 compute_diplacement_grad(size(), bs, vals, local_pts, p, displacement, displacement_grad);
260
261 const Eigen::MatrixXd def_grad = I + displacement_grad;
262 if (type == ElasticityTensorType::F)
263 {
264 all.row(p) = fun(def_grad);
265 continue;
266 }
267
268 double lambda, mu;
269 params_.lambda_mu(local_pts.row(p), vals.val.row(p), t, vals.element_id, lambda, mu);
270
271 Eigen::MatrixXd stress_tensor;
272 if (size() == 2)
273 stress_tensor = compute_stress_from_def_grad<2>(def_grad, lambda, mu) * def_grad.transpose() / def_grad.determinant();
274 else
275 stress_tensor = compute_stress_from_def_grad<3>(def_grad, lambda, mu) * def_grad.transpose() / def_grad.determinant();
276 if (type == ElasticityTensorType::PK1)
277 stress_tensor = pk1_from_cauchy(stress_tensor, def_grad);
278 else if (type == ElasticityTensorType::PK2)
279 stress_tensor = pk2_from_cauchy(stress_tensor, def_grad);
280
281 all.row(p) = fun(stress_tensor);
282 }
283 }
284
286 {
287 if (size() == 2)
288 return compute_energy_aux<2>(data);
289 else
290 return compute_energy_aux<3>(data);
291 }
292
293 template <int dim>
295 {
296 Eigen::VectorXd local_disp;
297 get_local_disp(data, dim, local_disp);
298
299 Eigen::Matrix<double, dim, dim> def_grad;
300
301 double energy = 0.0;
302
303 const int n_pts = data.da.size();
304 for (long p = 0; p < n_pts; ++p)
305 {
306 compute_disp_grad_at_quad(data, local_disp, p, dim, def_grad);
307
308 // Id + grad d
309 def_grad.diagonal().array() += 1.0;
310
311 double lambda, mu;
312 params_.lambda_mu(data.vals.quadrature.points.row(p), data.vals.val.row(p), data.t, data.vals.element_id, lambda, mu);
313
314 const double val = compute_energy_from_def_grad(def_grad, lambda, mu);
315
316 energy += val * data.da(p);
317 }
318 return energy;
319 }
320
321 template <int n_basis, int dim>
322 void FixedCorotational::compute_energy_aux_gradient_fast(const NonLinearAssemblerData &data, Eigen::Matrix<double, Eigen::Dynamic, 1> &G_flattened) const
323 {
324 assert(data.x.cols() == 1);
325
326 const int n_pts = data.da.size();
327
328 Eigen::Matrix<double, n_basis, dim> local_disp(data.vals.basis_values.size(), size());
329 local_disp.setZero();
330 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
331 {
332 const auto &bs = data.vals.basis_values[i];
333 for (size_t ii = 0; ii < bs.global.size(); ++ii)
334 {
335 for (int d = 0; d < size(); ++d)
336 {
337 local_disp(i, d) += bs.global[ii].val * data.x(bs.global[ii].index * size() + d);
338 }
339 }
340 }
341
342 Eigen::Matrix<double, dim, dim> def_grad(size(), size());
343
344 Eigen::Matrix<double, n_basis, dim> G(data.vals.basis_values.size(), size());
345 G.setZero();
346
347 for (long p = 0; p < n_pts; ++p)
348 {
349 Eigen::Matrix<double, n_basis, dim> grad(data.vals.basis_values.size(), size());
350
351 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
352 {
353 grad.row(i) = data.vals.basis_values[i].grad.row(p);
354 }
355
356 const Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
357 const Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
358
359 // Id + grad d
360 def_grad = local_disp.transpose() * delF_delU + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
361
362 double lambda, mu;
363 params_.lambda_mu(data.vals.quadrature.points.row(p), data.vals.val.row(p), data.t, data.vals.element_id, lambda, mu);
364
365 Eigen::Matrix<double, dim, dim> gradient_temp = compute_stress_from_def_grad(def_grad, lambda, mu);
366
367 Eigen::Matrix<double, n_basis, dim> gradient = delF_delU * gradient_temp.transpose();
368
369 G.noalias() += gradient * data.da(p);
370 }
371
372 Eigen::Matrix<double, dim, n_basis> G_T = G.transpose();
373
374 constexpr int N = (n_basis == Eigen::Dynamic) ? Eigen::Dynamic : n_basis * dim;
375 Eigen::Matrix<double, N, 1> temp(Eigen::Map<Eigen::Matrix<double, N, 1>>(G_T.data(), G_T.size()));
376 G_flattened = temp;
377 }
378
379 template <int n_basis, int dim>
381 {
382 assert(data.x.cols() == 1);
383
384 constexpr int N = (n_basis == Eigen::Dynamic) ? Eigen::Dynamic : n_basis * dim;
385 const int n_pts = data.da.size();
386
387 Eigen::Matrix<double, n_basis, dim> local_disp(data.vals.basis_values.size(), size());
388 local_disp.setZero();
389 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
390 {
391 const auto &bs = data.vals.basis_values[i];
392 for (size_t ii = 0; ii < bs.global.size(); ++ii)
393 {
394 for (int d = 0; d < size(); ++d)
395 {
396 local_disp(i, d) += bs.global[ii].val * data.x(bs.global[ii].index * size() + d);
397 }
398 }
399 }
400
401 Eigen::Matrix<double, dim, dim> def_grad(size(), size());
402
403 for (long p = 0; p < n_pts; ++p)
404 {
405 Eigen::Matrix<double, n_basis, dim> grad(data.vals.basis_values.size(), size());
406
407 for (size_t i = 0; i < data.vals.basis_values.size(); ++i)
408 {
409 grad.row(i) = data.vals.basis_values[i].grad.row(p);
410 }
411
412 Eigen::Matrix<double, dim, dim> jac_it = data.vals.jac_it[p];
413 const Eigen::Matrix<double, n_basis, dim> delF_delU = grad * jac_it;
414
415 // Id + grad d
416 def_grad = local_disp.transpose() * delF_delU + Eigen::Matrix<double, dim, dim>::Identity(size(), size());
417
418 double lambda, mu;
419 params_.lambda_mu(data.vals.quadrature.points.row(p), data.vals.val.row(p), data.t, data.vals.element_id, lambda, mu);
420
421 Eigen::Matrix<double, dim * dim, dim * dim> hessian_temp = compute_stiffness_from_def_grad(def_grad, lambda, mu);
422
423 Eigen::Matrix<double, dim * dim, N> delF_delU_tensor = Eigen::Matrix<double, dim * dim, N>::Zero(jac_it.size(), grad.size());
424
425 for (size_t i = 0; i < local_disp.rows(); ++i)
426 {
427 for (size_t j = 0; j < dim; ++j)
428 for (size_t k = 0; k < dim; ++k)
429 delF_delU_tensor(dim * k + j, i * dim + j) = delF_delU(i, k);
430 }
431
432 Eigen::Matrix<double, N, N> hessian = delF_delU_tensor.transpose() * hessian_temp * delF_delU_tensor;
433
434 H += hessian * data.da(p);
435 }
436 }
437
439 const OptAssemblerData &data,
440 const Eigen::MatrixXd &mat,
441 Eigen::MatrixXd &stress,
442 Eigen::MatrixXd &result) const
443 {
444 const double t = data.t;
445 const int el_id = data.el_id;
446 const Eigen::MatrixXd &local_pts = data.local_pts;
447 const Eigen::MatrixXd &global_pts = data.global_pts;
448 const Eigen::MatrixXd &grad_u_i = data.grad_u_i;
449
450 double lambda, mu;
451 params_.lambda_mu(local_pts, global_pts, t, el_id, lambda, mu);
452
453 const Eigen::MatrixXd def_grad = Eigen::MatrixXd::Identity(grad_u_i.rows(), grad_u_i.cols()) + grad_u_i;
454 if (size() == 2)
455 {
456 stress = compute_stress_from_def_grad<2>(def_grad, lambda, mu);
457 result = (compute_stiffness_from_def_grad<2>(def_grad, lambda, mu) * mat.reshaped(size() * size(), 1)).reshaped(size(), size());
458 }
459 else
460 {
461 stress = compute_stress_from_def_grad<3>(def_grad, lambda, mu);
462 result = (compute_stiffness_from_def_grad<3>(def_grad, lambda, mu) * mat.reshaped(size() * size(), 1)).reshaped(size(), size());
463 }
464 }
465
467 const OptAssemblerData &data,
468 Eigen::MatrixXd &stress,
469 Eigen::MatrixXd &result) const
470 {
471 const double t = data.t;
472 const int el_id = data.el_id;
473 const Eigen::MatrixXd &local_pts = data.local_pts;
474 const Eigen::MatrixXd &global_pts = data.global_pts;
475 const Eigen::MatrixXd &grad_u_i = data.grad_u_i;
476
477 double lambda, mu;
478 params_.lambda_mu(local_pts, global_pts, t, el_id, lambda, mu);
479
480 Eigen::MatrixXd def_grad = Eigen::MatrixXd::Identity(grad_u_i.rows(), grad_u_i.cols()) + grad_u_i;
481 if (size() == 2)
482 {
483 stress = compute_stress_from_def_grad<2>(def_grad, lambda, mu);
484 result = (compute_stiffness_from_def_grad<2>(def_grad, lambda, mu) * stress.reshaped(size() * size(), 1)).reshaped(size(), size());
485 }
486 else
487 {
488 stress = compute_stress_from_def_grad<3>(def_grad, lambda, mu);
489 result = (compute_stiffness_from_def_grad<3>(def_grad, lambda, mu) * stress.reshaped(size() * size(), 1)).reshaped(size(), size());
490 }
491 }
492
494 const OptAssemblerData &data,
495 Eigen::MatrixXd &dstress_dmu,
496 Eigen::MatrixXd &dstress_dlambda) const
497 {
498 const double t = data.t;
499 const int el_id = data.el_id;
500 const Eigen::MatrixXd &local_pts = data.local_pts;
501 const Eigen::MatrixXd &global_pts = data.global_pts;
502 const Eigen::MatrixXd &grad_u_i = data.grad_u_i;
503
504 const Eigen::MatrixXd def_grad = Eigen::MatrixXd::Identity(size(), size()) + grad_u_i;
505
506 Eigen::MatrixXd UVT;
507 Eigen::VectorXd sigmas;
508 {
509 if (size() == 2)
510 {
511 utils::AutoFlipSVD<Eigen::Matrix<double, 2, 2>> svd(def_grad, Eigen::ComputeFullU | Eigen::ComputeFullV);
512 sigmas = svd.singularValues();
513 UVT = svd.matrixU() * svd.matrixV().transpose();
514 }
515 else
516 {
517 utils::AutoFlipSVD<Eigen::Matrix<double, 3, 3>> svd(def_grad, Eigen::ComputeFullU | Eigen::ComputeFullV);
518 sigmas = svd.singularValues();
519 UVT = svd.matrixU() * svd.matrixV().transpose();
520 }
521 }
522
523 Eigen::MatrixXd delJ_delF(size(), size());
524 delJ_delF.setZero();
525
526 if (size() == 2)
527 {
528 delJ_delF(0, 0) = def_grad(1, 1);
529 delJ_delF(0, 1) = -def_grad(1, 0);
530 delJ_delF(1, 0) = -def_grad(0, 1);
531 delJ_delF(1, 1) = def_grad(0, 0);
532 }
533 else if (size() == 3)
534 {
535 const Eigen::Matrix<double, 3, 1> u = def_grad.col(0);
536 const Eigen::Matrix<double, 3, 1> v = def_grad.col(1);
537 const Eigen::Matrix<double, 3, 1> w = def_grad.col(2);
538
539 delJ_delF.col(0) = cross<3>(v, w);
540 delJ_delF.col(1) = cross<3>(w, u);
541 delJ_delF.col(2) = cross<3>(u, v);
542 }
543
544 dstress_dmu = 2 * (def_grad - UVT);
545 dstress_dlambda = (sigmas.prod() - 1) * delJ_delF;
546 }
547
548 std::map<std::string, Assembler::ParamFunc> FixedCorotational::parameters() const
549 {
550 std::map<std::string, ParamFunc> res;
551 const auto &params = params_;
552 const int size = this->size();
553
554 res["lambda"] = [&params](const RowVectorNd &uv, const RowVectorNd &p, double t, int e) {
555 double lambda, mu;
556
557 params.lambda_mu(uv, p, t, e, lambda, mu);
558 return lambda;
559 };
560
561 res["mu"] = [&params](const RowVectorNd &uv, const RowVectorNd &p, double t, int e) {
562 double lambda, mu;
563
564 params.lambda_mu(uv, p, t, e, lambda, mu);
565 return mu;
566 };
567
568 res["E"] = [&params, size](const RowVectorNd &uv, const RowVectorNd &p, double t, int e) {
569 double lambda, mu;
570 params.lambda_mu(uv, p, t, e, lambda, mu);
571
572 if (size == 3)
573 return mu * (3.0 * lambda + 2.0 * mu) / (lambda + mu);
574 else
575 return 2 * mu * (2.0 * lambda + 2.0 * mu) / (lambda + 2.0 * mu);
576 };
577
578 res["nu"] = [&params, size](const RowVectorNd &uv, const RowVectorNd &p, double t, int e) {
579 double lambda, mu;
580
581 params.lambda_mu(uv, p, t, e, lambda, mu);
582
583 if (size == 3)
584 return lambda / (2.0 * (lambda + mu));
585 else
586 return lambda / (lambda + 2.0 * mu);
587 };
588
589 return res;
590 }
591
592 template <int dim>
593 double FixedCorotational::compute_energy_from_singular_values(const Eigen::Vector<double, dim> &sigmas, const double lambda, const double mu)
594 {
595 const double sigmam12Sum = (sigmas - Eigen::Matrix<double, dim, 1>::Ones()).squaredNorm();
596 const double sigmaProdm1 = sigmas.prod() - 1.0;
597
598 return mu * sigmam12Sum + lambda / 2.0 * sigmaProdm1 * sigmaProdm1;
599 }
600
601 template <int dim>
602 Eigen::Vector<double, dim> FixedCorotational::compute_stress_from_singular_values(const Eigen::Vector<double, dim> &sigmas, const double lambda, const double mu)
603 {
604 const double sigmaProdm1lambda = lambda * (sigmas.prod() - 1.0);
605 Eigen::Matrix<double, dim, 1> sigmaProd_noI;
606 if constexpr (dim == 2) {
607 sigmaProd_noI[0] = sigmas[1];
608 sigmaProd_noI[1] = sigmas[0];
609 }
610 else {
611 sigmaProd_noI[0] = sigmas[1] * sigmas[2];
612 sigmaProd_noI[1] = sigmas[2] * sigmas[0];
613 sigmaProd_noI[2] = sigmas[0] * sigmas[1];
614 }
615
616 const double _2mu = mu * 2;
617 Eigen::Vector<double, dim> dE_div_dsigma;
618 dE_div_dsigma.setZero();
619 dE_div_dsigma[0] = (_2mu * (sigmas[0] - 1.0) + sigmaProd_noI[0] * sigmaProdm1lambda);
620 dE_div_dsigma[1] = (_2mu * (sigmas[1] - 1.0) + sigmaProd_noI[1] * sigmaProdm1lambda);
621 if constexpr (dim == 3) {
622 dE_div_dsigma[2] = (_2mu * (sigmas[2] - 1.0) + sigmaProd_noI[2] * sigmaProdm1lambda);
623 }
624 return dE_div_dsigma;
625 }
626
627 template <int dim>
628 Eigen::Matrix<double, dim, dim> FixedCorotational::compute_stiffness_from_singular_values(const Eigen::Vector<double, dim> &sigmas, const double lambda, const double mu)
629 {
630 const double sigmaProd = sigmas.prod();
631 Eigen::Matrix<double, dim, 1> sigmaProd_noI;
632 if constexpr (dim == 2) {
633 sigmaProd_noI[0] = sigmas[1];
634 sigmaProd_noI[1] = sigmas[0];
635 }
636 else {
637 sigmaProd_noI[0] = sigmas[1] * sigmas[2];
638 sigmaProd_noI[1] = sigmas[2] * sigmas[0];
639 sigmaProd_noI[2] = sigmas[0] * sigmas[1];
640 }
641
642 double _2mu = mu * 2;
643 Eigen::Matrix<double, dim, dim> d2E_div_dsigma2;
644 d2E_div_dsigma2.setZero();
645 d2E_div_dsigma2(0, 0) = _2mu + lambda * sigmaProd_noI[0] * sigmaProd_noI[0];
646 d2E_div_dsigma2(1, 1) = _2mu + lambda * sigmaProd_noI[1] * sigmaProd_noI[1];
647 if constexpr (dim == 3) {
648 d2E_div_dsigma2(2, 2) = _2mu + lambda * sigmaProd_noI[2] * sigmaProd_noI[2];
649 }
650
651 if constexpr (dim == 2) {
652 d2E_div_dsigma2(0, 1) = d2E_div_dsigma2(1, 0) = lambda * ((sigmaProd - 1.0) + sigmaProd_noI[0] * sigmaProd_noI[1]);
653 }
654 else {
655 d2E_div_dsigma2(0, 1) = d2E_div_dsigma2(1, 0) = lambda * (sigmas[2] * (sigmaProd - 1.0) + sigmaProd_noI[0] * sigmaProd_noI[1]);
656 d2E_div_dsigma2(0, 2) = d2E_div_dsigma2(2, 0) = lambda * (sigmas[1] * (sigmaProd - 1.0) + sigmaProd_noI[0] * sigmaProd_noI[2]);
657 d2E_div_dsigma2(2, 1) = d2E_div_dsigma2(1, 2) = lambda * (sigmas[0] * (sigmaProd - 1.0) + sigmaProd_noI[2] * sigmaProd_noI[1]);
658 }
659
660 return d2E_div_dsigma2;
661 }
662
663 template <int dim>
664 double FixedCorotational::compute_energy_from_def_grad(const Eigen::Matrix<double, dim, dim> &F, const double lambda, const double mu)
665 {
666 const Eigen::Vector<double, dim> sigmas = utils::singular_values<dim>(F);
667 return compute_energy_from_singular_values(sigmas, lambda, mu);
668 }
669
670 template <int dim>
671 Eigen::Matrix<double, dim, dim> FixedCorotational::compute_stress_from_def_grad(const Eigen::Matrix<double, dim, dim> &F, const double lambda, const double mu)
672 {
673 utils::AutoFlipSVD<Eigen::Matrix<double, dim, dim>> svd(F, Eigen::ComputeFullU | Eigen::ComputeFullV);
674 const Eigen::Vector<double, dim> sigmas = svd.singularValues();
675
676 Eigen::Matrix<double, dim, dim> delJ_delF(dim, dim);
677 delJ_delF.setZero();
678
679 if (dim == 2)
680 {
681 delJ_delF(0, 0) = F(1, 1);
682 delJ_delF(0, 1) = -F(1, 0);
683 delJ_delF(1, 0) = -F(0, 1);
684 delJ_delF(1, 1) = F(0, 0);
685 }
686 else if (dim == 3)
687 {
688 const Eigen::Matrix<double, dim, 1> u = F.col(0);
689 const Eigen::Matrix<double, dim, 1> v = F.col(1);
690 const Eigen::Matrix<double, dim, 1> w = F.col(2);
691
692 delJ_delF.col(0) = cross<dim>(v, w);
693 delJ_delF.col(1) = cross<dim>(w, u);
694 delJ_delF.col(2) = cross<dim>(u, v);
695 }
696
697 return lambda * (sigmas.prod() - 1) * delJ_delF + mu * 2 * (F - svd.matrixU() * svd.matrixV().transpose());
698 }
699
700 template <int dim>
701 Eigen::Matrix<double, dim*dim, dim*dim> FixedCorotational::compute_stiffness_from_def_grad(const Eigen::Matrix<double, dim, dim> &F, const double lambda, const double mu)
702 {
703 utils::AutoFlipSVD<Eigen::Matrix<double, dim, dim>> svd(F, Eigen::ComputeFullU | Eigen::ComputeFullV);
704 const Eigen::Vector<double, dim> sigmas = svd.singularValues();
705 Eigen::Matrix<double, dim, 1> dE_div_dsigma = compute_stress_from_singular_values(sigmas, lambda, mu);
706 Eigen::Matrix<double, dim, dim> d2E_div_dsigma2 = compute_stiffness_from_singular_values(sigmas, lambda, mu);
707
708 constexpr int Cdim2 = dim * (dim - 1) / 2;
709 Eigen::Matrix<double, Cdim2, 1> BLeftCoef;
710 BLeftCoef.setZero();
711 {
712 const double tmp = lambda / 2.0 * (sigmas.prod() - 1);
713 if constexpr (dim == 2) {
714 BLeftCoef[0] = mu - tmp;
715 }
716 else {
717 BLeftCoef[0] = mu - tmp * sigmas[2];
718 BLeftCoef[1] = mu - tmp * sigmas[0];
719 BLeftCoef[2] = mu - tmp * sigmas[1];
720 }
721 }
722 Eigen::Matrix2d B[Cdim2];
723 for (int cI = 0; cI < Cdim2; cI++) {
724 B[cI].setZero();
725 int cI_post = (cI + 1) % dim;
726
727 double rightCoef = dE_div_dsigma[cI] + dE_div_dsigma[cI_post];
728 double sum_sigma = sigmas[cI] + sigmas[cI_post];
729 rightCoef /= 2.0 * std::max(sum_sigma, 1.0e-12);
730
731 const double& leftCoef = BLeftCoef[cI];
732 B[cI](0, 0) = B[cI](1, 1) = leftCoef + rightCoef;
733 B[cI](0, 1) = B[cI](1, 0) = leftCoef - rightCoef;
734 }
735
736 // compute M using A(d2E_div_dsigma2) and B
737 Eigen::Matrix<double, dim * dim, dim * dim> M;
738 M.setZero();
739 if constexpr (dim == 2) {
740 M(0, 0) = d2E_div_dsigma2(0, 0);
741 M(0, 3) = d2E_div_dsigma2(0, 1);
742 M.block(1, 1, 2, 2) = B[0];
743 M(3, 0) = d2E_div_dsigma2(1, 0);
744 M(3, 3) = d2E_div_dsigma2(1, 1);
745 }
746 else {
747 // A
748 M(0, 0) = d2E_div_dsigma2(0, 0);
749 M(0, 4) = d2E_div_dsigma2(0, 1);
750 M(0, 8) = d2E_div_dsigma2(0, 2);
751 M(4, 0) = d2E_div_dsigma2(1, 0);
752 M(4, 4) = d2E_div_dsigma2(1, 1);
753 M(4, 8) = d2E_div_dsigma2(1, 2);
754 M(8, 0) = d2E_div_dsigma2(2, 0);
755 M(8, 4) = d2E_div_dsigma2(2, 1);
756 M(8, 8) = d2E_div_dsigma2(2, 2);
757 // B01
758 M(1, 1) = B[0](0, 0);
759 M(1, 3) = B[0](0, 1);
760 M(3, 1) = B[0](1, 0);
761 M(3, 3) = B[0](1, 1);
762 // B12
763 M(5, 5) = B[1](0, 0);
764 M(5, 7) = B[1](0, 1);
765 M(7, 5) = B[1](1, 0);
766 M(7, 7) = B[1](1, 1);
767 // B20
768 M(2, 2) = B[2](1, 1);
769 M(2, 6) = B[2](1, 0);
770 M(6, 2) = B[2](0, 1);
771 M(6, 6) = B[2](0, 0);
772 }
773
774 // compute hessian
775 Eigen::Matrix<double, dim*dim, dim*dim> hessian;
776 hessian.setZero();
777 const Eigen::Matrix<double, dim, dim>& U = svd.matrixU();
778 const Eigen::Matrix<double, dim, dim>& V = svd.matrixV();
779 for (int i = 0; i < dim; i++) {
780 // int _dim_i = i * dim;
781 for (int j = 0; j < dim; j++) {
782 int ij = i + j * dim;
783 for (int r = 0; r < dim; r++) {
784 // int _dim_r = r * dim;
785 for (int s = 0; s < dim; s++) {
786 int rs = r + s * dim;
787 if (ij > rs) {
788 // bottom left, same as upper right
789 continue;
790 }
791
792 if constexpr (dim == 2) {
793 hessian(ij, rs) = M(0, 0) * U(i, 0) * V(j, 0) * U(r, 0) * V(s, 0) + M(0, 3) * U(i, 0) * V(j, 0) * U(r, 1) * V(s, 1) + M(1, 1) * U(i, 0) * V(j, 1) * U(r, 0) * V(s, 1) + M(1, 2) * U(i, 0) * V(j, 1) * U(r, 1) * V(s, 0) + M(2, 1) * U(i, 1) * V(j, 0) * U(r, 0) * V(s, 1) + M(2, 2) * U(i, 1) * V(j, 0) * U(r, 1) * V(s, 0) + M(3, 0) * U(i, 1) * V(j, 1) * U(r, 0) * V(s, 0) + M(3, 3) * U(i, 1) * V(j, 1) * U(r, 1) * V(s, 1);
794 }
795 else {
796 hessian(ij, rs) = M(0, 0) * U(i, 0) * V(j, 0) * U(r, 0) * V(s, 0) + M(0, 4) * U(i, 0) * V(j, 0) * U(r, 1) * V(s, 1) + M(0, 8) * U(i, 0) * V(j, 0) * U(r, 2) * V(s, 2) + M(4, 0) * U(i, 1) * V(j, 1) * U(r, 0) * V(s, 0) + M(4, 4) * U(i, 1) * V(j, 1) * U(r, 1) * V(s, 1) + M(4, 8) * U(i, 1) * V(j, 1) * U(r, 2) * V(s, 2) + M(8, 0) * U(i, 2) * V(j, 2) * U(r, 0) * V(s, 0) + M(8, 4) * U(i, 2) * V(j, 2) * U(r, 1) * V(s, 1) + M(8, 8) * U(i, 2) * V(j, 2) * U(r, 2) * V(s, 2) + M(1, 1) * U(i, 0) * V(j, 1) * U(r, 0) * V(s, 1) + M(1, 3) * U(i, 0) * V(j, 1) * U(r, 1) * V(s, 0) + M(3, 1) * U(i, 1) * V(j, 0) * U(r, 0) * V(s, 1) + M(3, 3) * U(i, 1) * V(j, 0) * U(r, 1) * V(s, 0) + M(5, 5) * U(i, 1) * V(j, 2) * U(r, 1) * V(s, 2) + M(5, 7) * U(i, 1) * V(j, 2) * U(r, 2) * V(s, 1) + M(7, 5) * U(i, 2) * V(j, 1) * U(r, 1) * V(s, 2) + M(7, 7) * U(i, 2) * V(j, 1) * U(r, 2) * V(s, 1) + M(2, 2) * U(i, 0) * V(j, 2) * U(r, 0) * V(s, 2) + M(2, 6) * U(i, 0) * V(j, 2) * U(r, 2) * V(s, 0) + M(6, 2) * U(i, 2) * V(j, 0) * U(r, 0) * V(s, 2) + M(6, 6) * U(i, 2) * V(j, 0) * U(r, 2) * V(s, 0);
797 }
798
799 if (ij < rs)
800 hessian(rs, ij) = hessian(ij, rs);
801 }
802 }
803 }
804 }
805
806 return hessian;
807 }
808} // namespace polyfem::assembler
int V
double val
Definition Assembler.cpp:86
ElementAssemblyValues vals
Definition Assembler.cpp:22
int y
int z
int x
std::string stress() const
Definition Units.hpp:27
stores per element basis values at given quadrature points and geometric mapping
void compute(const int el_index, const bool is_volume, const Eigen::MatrixXd &pts, const basis::ElementBases &basis, const basis::ElementBases &gbasis)
computes the per element values at the local (ref el) points (pts) sets basis_values,...
std::vector< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > > jac_it
static double compute_energy_from_def_grad(const Eigen::Matrix< double, dim, dim > &F, const double lambda, const double mu)
static Eigen::Matrix< double, dim, dim > compute_stress_from_def_grad(const Eigen::Matrix< double, dim, dim > &F, const double lambda, const double mu)
void compute_energy_aux_gradient_fast(const NonLinearAssemblerData &data, Eigen::VectorXd &G_flattened) const
Eigen::MatrixXd assemble_hessian(const NonLinearAssemblerData &data) const override
static Eigen::Vector< double, dim > compute_stress_from_singular_values(const Eigen::Vector< double, dim > &sigmas, const double lambda, const double mu)
static double compute_energy_from_singular_values(const Eigen::Vector< double, dim > &sigmas, const double lambda, const double mu)
void add_multimaterial(const int index, const json &params, const Units &units) override
static Eigen::Matrix< double, dim *dim, dim *dim > compute_stiffness_from_def_grad(const Eigen::Matrix< double, dim, dim > &F, const double lambda, const double mu)
double compute_energy_aux(const NonLinearAssemblerData &data) const
double compute_energy(const NonLinearAssemblerData &data) const override
static Eigen::Matrix< double, dim, dim > compute_stiffness_from_singular_values(const Eigen::Vector< double, dim > &sigmas, const double lambda, const double mu)
std::map< std::string, ParamFunc > parameters() const override
void assign_stress_tensor(const OutputData &data, const int all_size, const ElasticityTensorType &type, Eigen::MatrixXd &all, const std::function< Eigen::MatrixXd(const Eigen::MatrixXd &)> &fun) const override
void compute_stress_grad_multiply_mat(const OptAssemblerData &data, const Eigen::MatrixXd &mat, Eigen::MatrixXd &stress, Eigen::MatrixXd &result) const override
void compute_energy_hessian_aux_fast(const NonLinearAssemblerData &data, Eigen::MatrixXd &H) const
Eigen::VectorXd assemble_gradient(const NonLinearAssemblerData &data) const override
void compute_stress_grad_multiply_stress(const OptAssemblerData &data, Eigen::MatrixXd &stress, Eigen::MatrixXd &result) const override
void compute_dstress_dmu_dlambda(const OptAssemblerData &data, Eigen::MatrixXd &dstress_dmu, Eigen::MatrixXd &dstress_dlambda) const override
void lambda_mu(double px, double py, double pz, double x, double y, double z, double t, int el_id, double &lambda, double &mu) const
void add_multimaterial(const int index, const json &params, const bool is_volume, const std::string &stress_unit)
const basis::ElementBases & bs
const Eigen::MatrixXd & fun
const Eigen::MatrixXd & local_pts
const basis::ElementBases & gbs
const Eigen::JacobiSVD< MatrixType >::SingularValuesType & singularValues(void) const
Definition svd.hpp:317
const MatrixType & matrixV(void) const
Definition svd.hpp:335
const MatrixType & matrixU(void) const
Definition svd.hpp:326
Used for test only.
Eigen::Matrix< double, dim, dim > hat(const Eigen::Matrix< double, dim, 1 > &x)
Eigen::Matrix< double, dim, 1 > cross(const Eigen::Matrix< double, dim, 1 > &x, const Eigen::Matrix< double, dim, 1 > &y)
Eigen::MatrixXd pk2_from_cauchy(const Eigen::MatrixXd &stress, const Eigen::MatrixXd &F)
void compute_disp_grad_at_quad(const assembler::NonLinearAssemblerData &data, const AutoDiffVect &local_disp, const int p, const int size, AutoDiffGradMat &def_grad)
nlohmann::json json
Definition Common.hpp:9
void compute_diplacement_grad(const int size, const ElementAssemblyValues &vals, const Eigen::MatrixXd &local_pts, const int p, const Eigen::MatrixXd &displacement, Eigen::MatrixXd &displacement_grad)
void get_local_disp(const assembler::NonLinearAssemblerData &data, const int size, AutoDiffVect &local_disp)
Eigen::MatrixXd pk1_from_cauchy(const Eigen::MatrixXd &stress, const Eigen::MatrixXd &F)
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13