PolyFEM
Loading...
Searching...
No Matches
PressureAssembler.cpp
Go to the documentation of this file.
1#include "PressureAssembler.hpp"
2
3#include <polyfem/assembler/Mass.hpp>
4#include <polyfem/utils/BoundarySampler.hpp>
5#include <polyfem/utils/Logger.hpp>
6#include <polyfem/utils/MaybeParallelFor.hpp>
7#include <ipc/utils/eigen_ext.hpp>
8#include <polysolve/linear/Solver.hpp>
9
10#include <polyfem/autogen/auto_p_bases.hpp>
11#include <polyfem/io/Evaluator.hpp>
12
13#include <igl/boundary_loop.h>
14
15namespace polyfem
16{
17 using namespace polysolve;
18 using namespace mesh;
19 using namespace quadrature;
20 using namespace utils;
21
22 namespace assembler
23 {
24 namespace
25 {
26 class LocalThreadScalarStorage
27 {
28 public:
29 double val;
30
31 LocalThreadScalarStorage()
32 {
33 val = 0;
34 }
35 };
36
37 class LocalThreadVecStorage
38 {
39 public:
40 Eigen::MatrixXd vec;
41
42 LocalThreadVecStorage(const int size)
43 {
44 vec.resize(size, 1);
45 vec.setZero();
46 }
47 };
48
49 class LocalThreadPrimitiveStorage
50 {
51 public:
52 std::vector<Eigen::VectorXi> faces;
53
54 LocalThreadPrimitiveStorage() {}
55 };
56
57 class LocalThreadMatStorage
58 {
59 public:
60 std::vector<Eigen::Triplet<double>> entries;
61 };
62
63 Eigen::Matrix3d wedge_product(const Eigen::Vector3d &a, const Eigen::Vector3d &b)
64 {
65 return a * b.transpose() - b * a.transpose();
66 }
67
68 void g_dual(const Eigen::Vector3d &g1, const Eigen::Vector3d &g2, Eigen::Vector3d &g1_up, Eigen::Vector3d &g2_up)
69 {
70 g1_up = (wedge_product(g1, g2) * g2) / (g1.dot(wedge_product(g1, g2) * g2));
71 g2_up = (wedge_product(g2, g1) * g1) / (g2.dot(wedge_product(g2, g1) * g1));
72 }
73
74 void boundary_loop_2d(const Eigen::MatrixXi &F, std::vector<int> &L)
75 {
76 std::set<int> vertices_set;
77 for (int i = 0; i < F.rows(); ++i)
78 for (int j = 0; j < 2; ++j)
79 {
80 if (vertices_set.find(F(i, j)) != vertices_set.end())
81 vertices_set.erase(F(i, j));
82 else
83 vertices_set.insert(F(i, j));
84 }
85
86 L = std::vector<int>(vertices_set.begin(), vertices_set.end());
87 }
88 } // namespace
89
90 double PressureAssembler::compute_volume(
91 const Eigen::MatrixXd &displacement,
92 const std::vector<mesh::LocalBoundary> &local_boundary,
93 const int resolution,
94 const double t,
95 const bool multiply_pressure) const
96 {
97 double res = 0;
98
99 auto storage = utils::create_thread_storage(LocalThreadScalarStorage());
100
101 utils::maybe_parallel_for(local_boundary.size(), [&](int start, int end, int thread_id) {
102 LocalThreadScalarStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
103
104 Eigen::MatrixXd pressure_vals, g_3;
105 ElementAssemblyValues vals;
106 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
107 Eigen::VectorXd weights;
108 Eigen::VectorXi global_primitive_ids;
109 for (int lb_id = start; lb_id < end; ++lb_id)
110 {
111 const auto &lb = local_boundary[lb_id];
112 const int e = lb.element_id();
113 const basis::ElementBases &gbs = gbases_[e];
114 const basis::ElementBases &bs = bases_[e];
115
116 for (int i = 0; i < lb.size(); ++i)
117 {
118 const int primitive_global_id = lb.global_primitive_id(i);
119 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
120
121 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
122 if (mesh_.is_volume())
123 weights /= 2 * mesh_.tri_area(primitive_global_id);
124 else
125 weights /= mesh_.edge_length(primitive_global_id);
126 g_3.setZero(normals.rows(), normals.cols());
127
128 if (!has_samples)
129 continue;
130
131 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
132
133 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
134
135 for (int n = 0; n < vals.jac_it.size(); ++n)
136 {
137 trafo = vals.jac_it[n].inverse();
138
139 if (displacement.size() > 0)
140 {
141 assert(size_ == 2 || size_ == 3);
142 deform_mat.resize(size_, size_);
143 deform_mat.setZero();
144 for (const auto &b : vals.basis_values)
145 {
146 for (const auto &g : b.global)
147 {
148 for (int d = 0; d < size_; ++d)
149 {
150 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
151 }
152 }
153 }
154
155 trafo += deform_mat;
156 }
157
158 normals.row(n) = normals.row(n) * trafo.inverse();
159 normals.row(n).normalize();
160
161 if (mesh_.is_volume())
162 {
163 Eigen::Vector3d g1, g2, g3;
164 auto endpoints = utils::BoundarySampler::tet_local_node_coordinates_from_face(lb[i]);
165 g1 = trafo * (endpoints.row(0) - endpoints.row(1)).transpose();
166 g2 = trafo * (endpoints.row(0) - endpoints.row(2)).transpose();
167 if (lb[i] == 0)
168 g1 *= -1;
169 g3 = g1.cross(g2);
170 g_3.row(n) = g3.transpose();
171 }
172 else
173 {
174 Eigen::Vector2d g1, g3;
175 auto endpoints = utils::BoundarySampler::tri_local_node_coordinates_from_edge(lb[i]);
176 g1 = trafo * (endpoints.row(0) - endpoints.row(1)).transpose();
177 g3(0) = -g1(1);
178 g3(1) = g1(0);
179 g_3.row(n) = g3.transpose();
180 }
181 }
182
183 if (multiply_pressure)
184 problem_.pressure_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, pressure_vals);
185 else
186 pressure_vals = Eigen::MatrixXd::Ones(weights.size(), 1);
187
188 Eigen::MatrixXd u, grad_u;
189 if (displacement.size() > 0)
190 io::Evaluator::interpolate_at_local_vals(mesh_, problem_.is_scalar(), bases_, gbases_, e, points, displacement, u, grad_u);
191 else
192 u.setZero(weights.size(), size_);
193 u += vals.val;
194
195 for (long p = 0; p < weights.size(); ++p)
196 for (int d = 0; d < size_; ++d)
197 local_storage.val += pressure_vals(p) * g_3(p, d) * u(p, d) * weights(p);
198 }
199 }
200 });
201
202 for (const LocalThreadScalarStorage &local_storage : storage)
203 res += local_storage.val;
204
205 res *= 1. / size_;
206
207 return res;
208 }
209
210 void PressureAssembler::compute_grad_volume(
211 const Eigen::MatrixXd &displacement,
212 const std::vector<mesh::LocalBoundary> &local_boundary,
213 const std::vector<int> &dirichlet_nodes,
214 const int resolution,
215 Eigen::VectorXd &grad,
216 const double t,
217 const bool multiply_pressure) const
218 {
219 grad.setZero(n_basis_ * size_);
220
221 auto storage = utils::create_thread_storage(LocalThreadVecStorage(grad.size()));
222
223 utils::maybe_parallel_for(local_boundary.size(), [&](int start, int end, int thread_id) {
224 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
225
226 Eigen::MatrixXd pressure_vals, g_3;
227 ElementAssemblyValues vals;
228 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
229 Eigen::VectorXd weights;
230 Eigen::VectorXi global_primitive_ids;
231 for (int lb_id = start; lb_id < end; ++lb_id)
232 {
233 const auto &lb = local_boundary[lb_id];
234 const int e = lb.element_id();
235 const basis::ElementBases &gbs = gbases_[e];
236 const basis::ElementBases &bs = bases_[e];
237
238 for (int i = 0; i < lb.size(); ++i)
239 {
240 const int primitive_global_id = lb.global_primitive_id(i);
241 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
242
243 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
244 if (mesh_.is_volume())
245 weights /= 2 * mesh_.tri_area(primitive_global_id);
246 else
247 weights /= mesh_.edge_length(primitive_global_id);
248 g_3.setZero(normals.rows(), normals.cols());
249
250 if (!has_samples)
251 continue;
252
253 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
254
255 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
256 for (int n = 0; n < vals.jac_it.size(); ++n)
257 {
258 trafo = vals.jac_it[n].inverse();
259
260 if (displacement.size() > 0)
261 {
262 assert(size_ == 2 || size_ == 3);
263 deform_mat.resize(size_, size_);
264 deform_mat.setZero();
265 for (const auto &b : vals.basis_values)
266 {
267 for (const auto &g : b.global)
268 {
269 for (int d = 0; d < size_; ++d)
270 {
271 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
272 }
273 }
274 }
275
276 trafo += deform_mat;
277 }
278
279 normals.row(n) = normals.row(n) * trafo.inverse();
280 normals.row(n).normalize();
281
282 if (mesh_.is_volume())
283 {
284 Eigen::Vector3d g1, g2, g3;
285 auto endpoints = utils::BoundarySampler::tet_local_node_coordinates_from_face(lb[i]);
286 g1 = trafo * (endpoints.row(0) - endpoints.row(1)).transpose();
287 g2 = trafo * (endpoints.row(0) - endpoints.row(2)).transpose();
288 if (lb[i] == 0)
289 g1 *= -1;
290 g3 = g1.cross(g2);
291 g_3.row(n) = g3.transpose();
292 }
293 else
294 {
295 Eigen::Vector2d g1, g3;
296 auto endpoints = utils::BoundarySampler::tri_local_node_coordinates_from_edge(lb[i]);
297 g1 = trafo * (endpoints.row(0) - endpoints.row(1)).transpose();
298 g3(0) = -g1(1);
299 g3(1) = g1(0);
300 g_3.row(n) = g3.transpose();
301 }
302 }
303
304 if (multiply_pressure)
305 problem_.pressure_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, pressure_vals);
306 else
307 pressure_vals = Eigen::MatrixXd::Ones(weights.size(), 1);
308
309 for (long n = 0; n < nodes.size(); ++n)
310 {
311 const AssemblyValues &v = vals.basis_values[nodes(n)];
312 for (int d = 0; d < size_; ++d)
313 {
314 for (size_t g = 0; g < v.global.size(); ++g)
315 {
316 const int g_index = v.global[g].index * size_ + d;
317 const bool is_dof_dirichlet = std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), g_index) != dirichlet_nodes.end();
318 if (is_dof_dirichlet)
319 continue;
320
321 for (long p = 0; p < weights.size(); ++p)
322 {
323 local_storage.vec(g_index) += pressure_vals(p) * g_3(p, d) * v.val(p) * weights(p);
324 }
325 }
326 }
327 }
328 }
329 }
330 });
331 for (const LocalThreadVecStorage &local_storage : storage)
332 grad += local_storage.vec;
333 }
334
335 void PressureAssembler::compute_grad_volume_id(
336 const Eigen::MatrixXd &displacement,
337 const int boundary_id,
338 const std::vector<mesh::LocalBoundary> &local_boundary,
339 const std::vector<int> &dirichlet_nodes,
340 const int resolution,
341 Eigen::VectorXd &grad,
342 const double t,
343 const bool multiply_pressure) const
344 {
345 grad.setZero(n_basis_ * size_);
346
347 auto storage = utils::create_thread_storage(LocalThreadVecStorage(grad.size()));
348
349 utils::maybe_parallel_for(local_boundary.size(), [&](int start, int end, int thread_id) {
350 LocalThreadVecStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
351
352 Eigen::MatrixXd pressure_vals, g_3;
353 ElementAssemblyValues vals;
354 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
355 Eigen::VectorXd weights;
356 Eigen::VectorXi global_primitive_ids;
357 for (int lb_id = start; lb_id < end; ++lb_id)
358 {
359 const auto &lb = local_boundary[lb_id];
360 const int e = lb.element_id();
361 const basis::ElementBases &gbs = gbases_[e];
362 const basis::ElementBases &bs = bases_[e];
363
364 for (int i = 0; i < lb.size(); ++i)
365 {
366 const int primitive_global_id = lb.global_primitive_id(i);
367 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
368 const int curr_boundary_id = mesh_.get_boundary_id(primitive_global_id);
369
370 if (curr_boundary_id != boundary_id)
371 continue;
372
373 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
374 if (mesh_.is_volume())
375 weights /= 2 * mesh_.tri_area(primitive_global_id);
376 else
377 weights /= mesh_.edge_length(primitive_global_id);
378 g_3.setZero(normals.rows(), normals.cols());
379
380 if (!has_samples)
381 continue;
382
383 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
384
385 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
386 for (int n = 0; n < vals.jac_it.size(); ++n)
387 {
388 trafo = vals.jac_it[n].inverse();
389
390 if (displacement.size() > 0)
391 {
392 assert(size_ == 2 || size_ == 3);
393 deform_mat.resize(size_, size_);
394 deform_mat.setZero();
395 for (const auto &b : vals.basis_values)
396 {
397 for (const auto &g : b.global)
398 {
399 for (int d = 0; d < size_; ++d)
400 {
401 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
402 }
403 }
404 }
405
406 trafo += deform_mat;
407 }
408
409 normals.row(n) = normals.row(n) * trafo.inverse();
410 normals.row(n).normalize();
411
412 if (mesh_.is_volume())
413 {
414 Eigen::Vector3d g1, g2, g3;
415 auto endpoints = utils::BoundarySampler::tet_local_node_coordinates_from_face(lb[i]);
416 g1 = trafo * (endpoints.row(0) - endpoints.row(1)).transpose();
417 g2 = trafo * (endpoints.row(0) - endpoints.row(2)).transpose();
418 if (lb[i] == 0)
419 g1 *= -1;
420 g3 = g1.cross(g2);
421 g_3.row(n) = g3.transpose();
422 }
423 else
424 {
425 Eigen::Vector2d g1, g3;
426 auto endpoints = utils::BoundarySampler::tri_local_node_coordinates_from_edge(lb[i]);
427 g1 = trafo * (endpoints.row(0) - endpoints.row(1)).transpose();
428 g3(0) = -g1(1);
429 g3(1) = g1(0);
430 g_3.row(n) = g3.transpose();
431 }
432 }
433
434 if (multiply_pressure)
435 problem_.pressure_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, pressure_vals);
436 else
437 pressure_vals = Eigen::MatrixXd::Ones(weights.size(), 1);
438
439 for (long n = 0; n < nodes.size(); ++n)
440 {
441 const AssemblyValues &v = vals.basis_values[nodes(n)];
442 for (int d = 0; d < size_; ++d)
443 {
444 for (size_t g = 0; g < v.global.size(); ++g)
445 {
446 const int g_index = v.global[g].index * size_ + d;
447 const bool is_dof_dirichlet = std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), g_index) != dirichlet_nodes.end();
448 if (is_dof_dirichlet)
449 continue;
450
451 for (long p = 0; p < weights.size(); ++p)
452 {
453 local_storage.vec(g_index) += pressure_vals(p) * g_3(p, d) * v.val(p) * weights(p);
454 }
455 }
456 }
457 }
458 }
459 }
460 });
461 for (const LocalThreadVecStorage &local_storage : storage)
462 grad += local_storage.vec;
463 }
464
465 void PressureAssembler::compute_hess_volume_3d(
466 const Eigen::MatrixXd &displacement,
467 const std::vector<mesh::LocalBoundary> &local_boundary,
468 const std::vector<int> &dirichlet_nodes,
469 const int resolution,
470 StiffnessMatrix &hess,
471 const double t,
472 const bool multiply_pressure) const
473 {
474 hess.resize(n_basis_ * size_, n_basis_ * size_);
475
476 auto storage = create_thread_storage(LocalThreadMatStorage());
477
478 utils::maybe_parallel_for(local_boundary.size(), [&](int start, int end, int thread_id) {
479 LocalThreadMatStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
480
481 Eigen::MatrixXd pressure_vals, g_1, g_2, g_3, local_hessian;
482 ElementAssemblyValues vals;
483 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
484 Eigen::VectorXd weights;
485 Eigen::VectorXi global_primitive_ids;
486 for (int lb_id = start; lb_id < end; ++lb_id)
487 {
488 const auto &lb = local_boundary[lb_id];
489 const int e = lb.element_id();
490 const basis::ElementBases &gbs = gbases_[e];
491 const basis::ElementBases &bs = bases_[e];
492
493 for (int i = 0; i < lb.size(); ++i)
494 {
495 const int primitive_global_id = lb.global_primitive_id(i);
496 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
497
498 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
499 std::vector<Eigen::VectorXd> param_chain_rule;
500 if (mesh_.is_volume())
501 {
502 weights /= 2 * mesh_.tri_area(primitive_global_id);
503 g_1.setZero(normals.rows(), normals.cols());
504 g_2.setZero(normals.rows(), normals.cols());
505 g_3.setZero(normals.rows(), normals.cols());
506
507 auto endpoints = utils::BoundarySampler::tet_local_node_coordinates_from_face(lb[i]);
508 param_chain_rule = {(endpoints.row(0) - endpoints.row(1)).transpose(), (endpoints.row(0) - endpoints.row(2)).transpose()};
509 if (lb[i] == 0)
510 param_chain_rule[0] *= -1;
511 }
512 else
513 assert(false);
514
515 if (!has_samples)
516 continue;
517
518 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
519
520 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
521 for (int n = 0; n < vals.jac_it.size(); ++n)
522 {
523 trafo = vals.jac_it[n].inverse();
524
525 if (displacement.size() > 0)
526 {
527 assert(size_ == 3);
528 deform_mat.resize(size_, size_);
529 deform_mat.setZero();
530 for (const auto &b : vals.basis_values)
531 {
532 for (const auto &g : b.global)
533 {
534 for (int d = 0; d < size_; ++d)
535 {
536 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
537 }
538 }
539 }
540
541 trafo += deform_mat;
542 }
543
544 normals.row(n) = normals.row(n) * trafo.inverse();
545 normals.row(n).normalize();
546
547 if (mesh_.is_volume())
548 {
549 Eigen::Vector3d g1, g2, g3;
550 auto endpoints = utils::BoundarySampler::tet_local_node_coordinates_from_face(lb[i]);
551 g1 = trafo * (endpoints.row(0) - endpoints.row(1)).transpose();
552 g2 = trafo * (endpoints.row(0) - endpoints.row(2)).transpose();
553 if (lb[i] == 0)
554 g1 *= -1;
555 g3 = g1.cross(g2);
556
557 g_1.row(n) = g1.transpose();
558 g_2.row(n) = g2.transpose();
559 g_3.row(n) = g3.transpose();
560 }
561 else
562 assert(false);
563 }
564
565 if (multiply_pressure)
566 problem_.pressure_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, pressure_vals);
567 else
568 pressure_vals = Eigen::MatrixXd::Ones(weights.size(), 1);
569
570 local_hessian.setZero(vals.basis_values.size() * size_, vals.basis_values.size() * size_);
571
572 for (long p = 0; p < weights.size(); ++p)
573 {
574 Eigen::Vector3d g_up_1, g_up_2;
575 g_dual(g_1.row(p).transpose(), g_2.row(p).transpose(), g_up_1, g_up_2);
576 std::vector<Eigen::MatrixXd> g_3_wedge_g_up = {wedge_product(g_3.row(p), g_up_1), wedge_product(g_3.row(p), g_up_2)};
577
578 for (long ni = 0; ni < nodes.size(); ++ni)
579 {
580 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
581 for (int di = 0; di < size_; ++di)
582 {
583 const int gi_index = vi.global[0].index * size_ + di;
584 const bool is_dof_i_dirichlet = std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), gi_index) != dirichlet_nodes.end();
585 if (is_dof_i_dirichlet)
586 continue;
587
588 Eigen::MatrixXd grad_phi_i;
589 grad_phi_i.setZero(3, 3);
590 grad_phi_i.row(di) = vi.grad.row(p);
591
592 for (long nj = 0; nj < nodes.size(); ++nj)
593 {
594 const AssemblyValues &vj = vals.basis_values[nodes(nj)];
595 for (int dj = 0; dj < size_; ++dj)
596 {
597 const int gj_index = vj.global[0].index * size_ + dj;
598 const bool is_dof_j_dirichlet = std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), gj_index) != dirichlet_nodes.end();
599 if (is_dof_j_dirichlet)
600 continue;
601
602 Eigen::MatrixXd grad_phi_j;
603 grad_phi_j.setZero(3, 3);
604 grad_phi_j.row(dj) = vj.grad.row(p);
605
606 double value = 0;
607 for (int alpha = 0; alpha < size_ - 1; ++alpha)
608 {
609 auto a = vj.val(p) * g_3_wedge_g_up[alpha].row(di) * (grad_phi_i * param_chain_rule[alpha]);
610 auto b = (g_3_wedge_g_up[alpha] * (grad_phi_j * param_chain_rule[alpha])).row(di) * vi.val(p);
611 value += pressure_vals(p) * (a(0) + b(0)) * weights(p);
612 }
613 local_hessian(nodes(ni) * size_ + di, nodes(nj) * size_ + dj) += value;
614 }
615 }
616 }
617 }
618 }
619
620 for (long ni = 0; ni < nodes.size(); ++ni)
621 {
622 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
623 for (int di = 0; di < size_; ++di)
624 {
625 const int gi_index = vi.global[0].index * size_ + di;
626 for (long nj = 0; nj < nodes.size(); ++nj)
627 {
628 const AssemblyValues &vj = vals.basis_values[nodes(nj)];
629 for (int dj = 0; dj < size_; ++dj)
630 {
631 const int gj_index = vj.global[0].index * size_ + dj;
632
633 local_storage.entries.push_back(Eigen::Triplet<double>(gi_index, gj_index, local_hessian(nodes(ni) * size_ + di, nodes(nj) * size_ + dj)));
634 }
635 }
636 }
637 }
638 }
639 }
640 });
641
642 std::vector<Eigen::Triplet<double>> all_entries;
643 // Serially merge local storages
644 for (LocalThreadMatStorage &local_storage : storage)
645 {
646 all_entries.insert(all_entries.end(), local_storage.entries.begin(), local_storage.entries.end());
647 }
648
649 hess.setFromTriplets(all_entries.begin(), all_entries.end());
650 }
651
652 void PressureAssembler::compute_hess_volume_2d(
653 const Eigen::MatrixXd &displacement,
654 const std::vector<mesh::LocalBoundary> &local_boundary,
655 const std::vector<int> &dirichlet_nodes,
656 const int resolution,
657 StiffnessMatrix &hess,
658 const double t,
659 const bool multiply_pressure) const
660 {
661 hess.resize(n_basis_ * size_, n_basis_ * size_);
662
663 auto storage = create_thread_storage(LocalThreadMatStorage());
664
665 utils::maybe_parallel_for(local_boundary.size(), [&](int start, int end, int thread_id) {
666 LocalThreadMatStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
667
668 Eigen::MatrixXd pressure_vals, g_3_grad, local_hessian;
669 ElementAssemblyValues vals;
670 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
671 Eigen::VectorXd weights;
672 Eigen::VectorXi global_primitive_ids;
673 for (int lb_id = start; lb_id < end; ++lb_id)
674 {
675 const auto &lb = local_boundary[lb_id];
676 const int e = lb.element_id();
677 const basis::ElementBases &gbs = gbases_[e];
678 const basis::ElementBases &bs = bases_[e];
679
680 for (int i = 0; i < lb.size(); ++i)
681 {
682 const int primitive_global_id = lb.global_primitive_id(i);
683 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
684
685 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
686 std::vector<Eigen::VectorXd> param_chain_rule;
687
688 if (!has_samples)
689 continue;
690
691 if (mesh_.is_volume())
692 assert(false);
693 else
694 weights /= mesh_.edge_length(primitive_global_id);
695
696 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
697
698 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
699 g_3_grad.setZero(normals.rows(), size_ * vals.basis_values.size() * size_);
700
701 for (int n = 0; n < vals.jac_it.size(); ++n)
702 {
703 trafo = vals.jac_it[n].inverse();
704
705 if (displacement.size() > 0)
706 {
707 assert(size_ == 2);
708 deform_mat.resize(size_, size_);
709 deform_mat.setZero();
710 for (const auto &b : vals.basis_values)
711 {
712 for (const auto &g : b.global)
713 {
714 for (int d = 0; d < size_; ++d)
715 {
716 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
717 }
718 }
719 }
720
721 trafo += deform_mat;
722 }
723
724 normals.row(n) = normals.row(n) * trafo.inverse();
725 normals.row(n).normalize();
726
727 if (mesh_.is_volume())
728 assert(false);
729 else
730 {
731 Eigen::MatrixXd g_3_grad_local(size_, vals.basis_values.size() * size_);
732 g_3_grad_local.setZero();
733 auto endpoints = utils::BoundarySampler::tri_local_node_coordinates_from_edge(lb[i]);
734 Eigen::VectorXd reference_normal = (endpoints.row(0) - endpoints.row(1)).transpose();
735 int b_count = 0;
736 for (const auto &b : vals.basis_values)
737 {
738 for (const auto &g : b.global)
739 for (int d = 0; d < size_; ++d)
740 for (int k = 0; k < size_; ++k)
741 g_3_grad_local(1 - d, b_count * size_ + d) += b.grad(n, k) * reference_normal(k);
742 ++b_count;
743 }
744 g_3_grad_local.row(0) *= -1;
745
746 for (int k = 0; k < size_; ++k)
747 for (int l = 0; l < vals.basis_values.size() * size_; ++l)
748 g_3_grad(n, k * (vals.basis_values.size() * size_) + l) = g_3_grad_local(k, l);
749 }
750 }
751
752 if (multiply_pressure)
753 problem_.pressure_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, pressure_vals);
754 else
755 pressure_vals = Eigen::MatrixXd::Ones(weights.size(), 1);
756
757 local_hessian.setZero(vals.basis_values.size() * size_, vals.basis_values.size() * size_);
758
759 for (long p = 0; p < weights.size(); ++p)
760 {
761 for (long ni = 0; ni < nodes.size(); ++ni)
762 {
763 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
764 for (int di = 0; di < size_; ++di)
765 {
766 const int gi_index = vi.global[0].index * size_ + di;
767 const bool is_dof_i_dirichlet = std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), gi_index) != dirichlet_nodes.end();
768 if (is_dof_i_dirichlet)
769 continue;
770
771 for (long nj = 0; nj < nodes.size(); ++nj)
772 {
773 const AssemblyValues &vj = vals.basis_values[nodes(nj)];
774 for (int dj = 0; dj < size_; ++dj)
775 {
776 const int gj_index = vj.global[0].index * size_ + dj;
777 const bool is_dof_j_dirichlet = std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), gj_index) != dirichlet_nodes.end();
778 if (is_dof_j_dirichlet)
779 continue;
780
781 double value = pressure_vals(p) * g_3_grad(p, (di * vals.basis_values.size() * size_) + nodes(nj) * size_ + dj) * vi.val(p) * weights(p);
782 local_hessian(nodes(ni) * size_ + di, nodes(nj) * size_ + dj) += value;
783 }
784 }
785 }
786 }
787 }
788
789 for (long ni = 0; ni < nodes.size(); ++ni)
790 {
791 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
792 for (int di = 0; di < size_; ++di)
793 {
794 const int gi_index = vi.global[0].index * size_ + di;
795 for (long nj = 0; nj < nodes.size(); ++nj)
796 {
797 const AssemblyValues &vj = vals.basis_values[nodes(nj)];
798 for (int dj = 0; dj < size_; ++dj)
799 {
800 const int gj_index = vj.global[0].index * size_ + dj;
801
802 local_storage.entries.push_back(Eigen::Triplet<double>(gi_index, gj_index, local_hessian(nodes(ni) * size_ + di, nodes(nj) * size_ + dj)));
803 }
804 }
805 }
806 }
807 }
808 }
809 });
810
811 std::vector<Eigen::Triplet<double>> all_entries;
812 // Serially merge local storages
813 for (LocalThreadMatStorage &local_storage : storage)
814 {
815 all_entries.insert(all_entries.end(), local_storage.entries.begin(), local_storage.entries.end());
816 }
817
818 hess.setFromTriplets(all_entries.begin(), all_entries.end());
819 }
820
821 bool PressureAssembler::is_closed_or_boundary_fixed(
822 const std::vector<mesh::LocalBoundary> &local_boundary,
823 const std::vector<int> &dirichlet_nodes) const
824 {
825 if (local_boundary.size() == 0)
826 return true;
827
828 auto storage = utils::create_thread_storage(LocalThreadPrimitiveStorage());
829
830 utils::maybe_parallel_for(local_boundary.size(), [&](int start, int end, int thread_id) {
831 LocalThreadPrimitiveStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
832
833 ElementAssemblyValues vals;
834 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
835 Eigen::VectorXd weights;
836 Eigen::VectorXi global_primitive_ids;
837 for (int lb_id = start; lb_id < end; ++lb_id)
838 {
839 const auto &lb = local_boundary[lb_id];
840 const int e = lb.element_id();
841 const basis::ElementBases &bs = bases_[e];
842 const basis::ElementBases &gbs = gbases_[e];
843
844 for (int i = 0; i < lb.size(); ++i)
845 {
846 const int primitive_global_id = lb.global_primitive_id(i);
847 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
848
849 points.resize(0, size_);
850 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
851
852 if (!((lb.type() == BoundaryType::TRI_LINE) || (lb.type() == BoundaryType::TRI)))
853 {
854 logger().warn("Detected element not triangle or tet, cannot check if pressure boundary is closed");
855 return;
856 }
857
858 Eigen::VectorXi face(size_);
859 // Assuming the geometric nodes are listed in the beginning
860 for (int n = 0; n < size_; ++n)
861 {
862 const AssemblyValues &v = vals.basis_values[nodes(n)];
863 assert(v.global.size() == 1);
864 face(n) = v.global[0].index;
865 }
866 local_storage.faces.push_back(face);
867 }
868 }
869 });
870
871 Eigen::MatrixXi F(0, size_);
872 int count = 0;
873 for (const LocalThreadPrimitiveStorage &local_storage : storage)
874 {
875 for (int i = 0; i < local_storage.faces.size(); ++i)
876 {
877 F.conservativeResize(count + 1, size_);
878 F.row(count) = local_storage.faces[i];
879 ++count;
880 }
881 }
882
883 std::vector<int> L;
884 if (size_ == 2)
885 boundary_loop_2d(F, L);
886 else
887 igl::boundary_loop(F, L);
888
889 for (const int &l : L)
890 if (std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), l * size_) == dirichlet_nodes.end())
891 return false;
892
893 return true;
894 }
895
896 PressureAssembler::PressureAssembler(const Assembler &assembler, const Mesh &mesh, const Obstacle &obstacle,
897 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
898 const std::unordered_map<int, std::vector<mesh::LocalBoundary>> &local_pressure_cavity,
899 const std::vector<int> &dirichlet_nodes,
900 const std::vector<int> &primitive_to_nodes,
901 const std::vector<int> &node_to_primitives,
902 const int n_basis, const int size,
903 const std::vector<basis::ElementBases> &bases, const std::vector<basis::ElementBases> &gbases,
904 const Problem &problem)
905 : assembler_(assembler),
906 mesh_(mesh),
907 obstacle_(obstacle),
908 n_basis_(n_basis),
909 size_(size),
910 bases_(bases),
911 gbases_(gbases),
912 problem_(problem),
913 primitive_to_nodes_(primitive_to_nodes),
914 node_to_primitives_(node_to_primitives)
915 {
916 for (const auto &v : local_pressure_cavity)
917 starting_volumes_[v.first] = compute_volume(Eigen::VectorXd(), v.second, 5, 0, false);
918
919 if (!is_closed_or_boundary_fixed(local_pressure_boundary, dirichlet_nodes))
920 logger().error("Pressure boundary condition must be applied to a closed volume or have dirichlet fixed boundary.");
921
922 for (const auto &b : local_pressure_cavity)
923 if (!is_closed_or_boundary_fixed(b.second, dirichlet_nodes))
924 logger().error("Pressure cavity boundary condition must be applied to a closed volume or have dirichlet fixed boundary.");
925
926 cavity_thermodynamics_ = std::make_unique<AdiabaticProcess>();
927 // cavity_thermodynamics_ = std::make_unique<IsothermalProcess>();
928 }
929
930 double PressureAssembler::compute_cavity_energy(
931 const Eigen::MatrixXd &displacement,
932 const std::unordered_map<int, std::vector<mesh::LocalBoundary>> &local_pressure_cavity,
933 const int resolution,
934 const double t) const
935 {
936 double energy_ = 0;
937 for (const auto &v : local_pressure_cavity)
938 {
939 const double start_pressure = problem_.pressure_cavity_bc(v.first, t);
940 const double start_volume = starting_volumes_.at(v.first);
941 const double curr_volume = compute_volume(displacement, v.second, resolution, t, false);
942
943 energy_ += cavity_thermodynamics_->energy(-start_pressure, -start_volume, -curr_volume);
944 }
945 return energy_;
946 }
947
948 void PressureAssembler::compute_cavity_energy_grad(
949 const Eigen::MatrixXd &displacement,
950 const std::unordered_map<int, std::vector<mesh::LocalBoundary>> &local_pressure_cavity,
951 const std::vector<int> &dirichlet_nodes,
952 const int resolution,
953 const double t,
954 Eigen::VectorXd &grad) const
955 {
956 grad.setZero(displacement.size());
957 for (const auto &v : local_pressure_cavity)
958 {
959 double start_pressure = problem_.pressure_cavity_bc(v.first, t);
960 double start_volume = starting_volumes_.at(v.first);
961 double curr_volume = compute_volume(displacement, v.second, resolution, t, false);
962 Eigen::VectorXd g;
963
964 compute_grad_volume(displacement, v.second, dirichlet_nodes, resolution, g, t, false);
965
966 double p = -cavity_thermodynamics_->pressure(
967 -start_pressure, -start_volume, -curr_volume);
968
969 grad += p * g;
970 }
971 }
972
973 void PressureAssembler::compute_cavity_energy_hess(
974 const Eigen::MatrixXd &displacement,
975 const std::unordered_map<int, std::vector<mesh::LocalBoundary>> &local_pressure_cavity,
976 const std::vector<int> &dirichlet_nodes,
977 const int resolution,
978 const double t,
979 const bool project_to_psd,
980 StiffnessMatrix &hess) const
981 {
982 if (project_to_psd && (local_pressure_cavity.size() > 0))
983 {
984 log_and_throw_error("Cannot project caivity pressure to PSD!");
985 }
986
987 hess.resize(displacement.size(), displacement.size());
988 for (const auto &v : local_pressure_cavity)
989 {
990 double start_pressure = problem_.pressure_cavity_bc(v.first, t);
991 double start_volume = starting_volumes_.at(v.first);
992 double curr_volume = compute_volume(displacement, v.second, resolution, t, false);
993 Eigen::VectorXd g;
994 StiffnessMatrix h;
995
996 compute_grad_volume(displacement, v.second, dirichlet_nodes, resolution, g, t, false);
997 if (size_ == 2)
998 compute_hess_volume_2d(displacement, v.second, dirichlet_nodes, resolution, h, t, false);
999 else if (size_ == 3)
1000 compute_hess_volume_3d(displacement, v.second, dirichlet_nodes, resolution, h, t, false);
1001
1002 double p = -cavity_thermodynamics_->pressure(
1003 -start_pressure, -start_volume, -curr_volume);
1004 double dp_dv = cavity_thermodynamics_->first_derivative(
1005 -start_pressure, -start_volume, -curr_volume);
1006
1007 hess += p * h;
1008 hess += dp_dv * (g * g.transpose()).sparseView();
1009 }
1010 }
1011
1012 double PressureAssembler::compute_energy(
1013 const Eigen::MatrixXd &displacement,
1014 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
1015 const int resolution,
1016 const double t) const
1017 {
1018 return compute_volume(displacement, local_pressure_boundary, resolution, t, true);
1019 }
1020
1021 void PressureAssembler::compute_energy_grad(
1022 const Eigen::MatrixXd &displacement,
1023 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
1024 const std::vector<int> &dirichlet_nodes,
1025 const int resolution,
1026 const double t,
1027 Eigen::VectorXd &grad) const
1028 {
1029 compute_grad_volume(displacement, local_pressure_boundary, dirichlet_nodes, resolution, grad, t, true);
1030 }
1031
1032 void PressureAssembler::compute_energy_hess(
1033 const Eigen::MatrixXd &displacement,
1034 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
1035 const std::vector<int> &dirichlet_nodes,
1036 const int resolution,
1037 const double t,
1038 const bool project_to_psd,
1039 StiffnessMatrix &hess) const
1040 {
1041 if (size_ == 2)
1042 compute_hess_volume_2d(displacement, local_pressure_boundary, dirichlet_nodes, resolution, hess, t, true);
1043 else if (size_ == 3)
1044 compute_hess_volume_3d(displacement, local_pressure_boundary, dirichlet_nodes, resolution, hess, t, true);
1045 }
1046
1047 void PressureAssembler::compute_force_jacobian(
1048 const Eigen::MatrixXd &displacement,
1049 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
1050 const std::vector<int> &dirichlet_nodes,
1051 const int resolution,
1052 const double t,
1053 const int n_vertices,
1054 StiffnessMatrix &hess) const
1055 {
1056 compute_energy_hess(displacement, local_pressure_boundary, dirichlet_nodes, resolution, t, false, hess);
1057 }
1058
1059 } // namespace assembler
1060} // namespace polyfem
vec
Eigen::MatrixXd vec
Definition Assembler.cpp:72
val
double val
Definition Assembler.cpp:86
quadrature
Quadrature quadrature
Definition MassMatrixAssembler.cpp:30
entries
std::vector< Eigen::Triplet< double > > entries
Definition MassMatrixAssembler.cpp:25
faces
std::vector< Eigen::VectorXi > faces
Definition PressureAssembler.cpp:52
auto_p_bases.hpp
polyfem::assembler::PressureAssembler::compute_energy
double compute_energy(const Eigen::MatrixXd &displacement, const std::vector< mesh::LocalBoundary > &local_pressure_boundary, const int resolution, const double t) const
Definition PressureAssembler.cpp:1012
polyfem::assembler::PressureAssembler::is_closed_or_boundary_fixed
bool is_closed_or_boundary_fixed(const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &dirichlet_nodes) const
Definition PressureAssembler.cpp:821
polyfem::assembler::PressureAssembler::cavity_thermodynamics_
std::unique_ptr< ThermodynamicProcess > cavity_thermodynamics_
Definition PressureAssembler.hpp:200
polyfem::assembler::PressureAssembler::compute_grad_volume_id
void compute_grad_volume_id(const Eigen::MatrixXd &displacement, const int boundary_id, const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &dirichlet_nodes, const int resolution, Eigen::VectorXd &grad, const double t=0, const bool multiply_pressure=false) const
Definition PressureAssembler.cpp:335
polyfem::assembler::PressureAssembler::compute_hess_volume_3d
void compute_hess_volume_3d(const Eigen::MatrixXd &displacement, const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &dirichlet_nodes, const int resolution, StiffnessMatrix &hess, const double t=0, const bool multiply_pressure=false) const
Definition PressureAssembler.cpp:465
polyfem::assembler::PressureAssembler::compute_hess_volume_2d
void compute_hess_volume_2d(const Eigen::MatrixXd &displacement, const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &dirichlet_nodes, const int resolution, StiffnessMatrix &hess, const double t=0, const bool multiply_pressure=false) const
Definition PressureAssembler.cpp:652
polyfem::assembler::PressureAssembler::compute_energy_hess
void compute_energy_hess(const Eigen::MatrixXd &displacement, const std::vector< mesh::LocalBoundary > &local_pressure_boundary, const std::vector< int > &dirichlet_nodes, const int resolution, const double t, const bool project_to_psd, StiffnessMatrix &hess) const
Definition PressureAssembler.cpp:1032
polyfem::assembler::PressureAssembler::compute_volume
double compute_volume(const Eigen::MatrixXd &displacement, const std::vector< mesh::LocalBoundary > &local_boundary, const int resolution, const double t=0, const bool multiply_pressure=false) const
Definition PressureAssembler.cpp:90
polyfem::assembler::PressureAssembler::compute_grad_volume
void compute_grad_volume(const Eigen::MatrixXd &displacement, const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &dirichlet_nodes, const int resolution, Eigen::VectorXd &grad, const double t=0, const bool multiply_pressure=false) const
Definition PressureAssembler.cpp:210
polyfem::assembler::PressureAssembler::compute_energy_grad
void compute_energy_grad(const Eigen::MatrixXd &displacement, const std::vector< mesh::LocalBoundary > &local_pressure_boundary, const std::vector< int > &dirichlet_nodes, const int resolution, const double t, Eigen::VectorXd &grad) const
Definition PressureAssembler.cpp:1021
polyfem::assembler::PressureAssembler::compute_force_jacobian
void compute_force_jacobian(const Eigen::MatrixXd &displacement, const std::vector< mesh::LocalBoundary > &local_pressure_boundary, const std::vector< int > &dirichlet_nodes, const int resolution, const double t, const int n_vertices, StiffnessMatrix &hess) const
Definition PressureAssembler.cpp:1047
polyfem::assembler::PressureAssembler::PressureAssembler
PressureAssembler(const Assembler &assembler, const mesh::Mesh &mesh, const mesh::Obstacle &obstacle, const std::vector< mesh::LocalBoundary > &local_pressure_boundary, const std::unordered_map< int, std::vector< mesh::LocalBoundary > > &local_pressure_cavity, const std::vector< int > &dirichlet_nodes, const std::vector< int > &primitive_to_nodes, const std::vector< int > &node_to_primitives, const int n_basis, const int size, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const Problem &problem)
Definition PressureAssembler.cpp:896
polyfem::assembler::PressureAssembler::starting_volumes_
std::unordered_map< int, double > starting_volumes_
Definition PressureAssembler.hpp:199
polyfem::assembler::PressureAssembler::compute_cavity_energy_grad
void compute_cavity_energy_grad(const Eigen::MatrixXd &displacement, const std::unordered_map< int, std::vector< mesh::LocalBoundary > > &local_pressure_cavity, const std::vector< int > &dirichlet_nodes, const int resolution, const double t, Eigen::VectorXd &grad) const
Definition PressureAssembler.cpp:948
polyfem::assembler::PressureAssembler::compute_cavity_energy_hess
void compute_cavity_energy_hess(const Eigen::MatrixXd &displacement, const std::unordered_map< int, std::vector< mesh::LocalBoundary > > &local_pressure_cavity, const std::vector< int > &dirichlet_nodes, const int resolution, const double t, const bool project_to_psd, StiffnessMatrix &hess) const
Definition PressureAssembler.cpp:973
polyfem::assembler::PressureAssembler::compute_cavity_energy
double compute_cavity_energy(const Eigen::MatrixXd &displacement, const std::unordered_map< int, std::vector< mesh::LocalBoundary > > &local_pressure_cavity, const int resolution, const double t) const
Definition PressureAssembler.cpp:930
polyfem::assembler::Problem::pressure_cavity_bc
virtual double pressure_cavity_bc(const int boundary_id, const double t) const
Definition Problem.hpp:33
polyfem::mesh::Mesh
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:39
polyfem::mesh::find
int find(const Eigen::VectorXi &vec, int x)
Definition NCMesh2D.cpp:289
polyfem::utils::create_thread_storage
auto create_thread_storage(const LocalStorage &initial_local_storage)
polyfem::utils::maybe_parallel_for
void maybe_parallel_for(int size, const std::function< void(int, int, int)> &partial_for)
polyfem::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:42
polyfem::log_and_throw_error
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:71
polyfem::StiffnessMatrix
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22