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_hess_volume_3d(
336 const Eigen::MatrixXd &displacement,
337 const std::vector<mesh::LocalBoundary> &local_boundary,
338 const std::vector<int> &dirichlet_nodes,
339 const int resolution,
340 StiffnessMatrix &hess,
341 const double t,
342 const bool multiply_pressure) const
343 {
344 hess.resize(n_basis_ * size_, n_basis_ * size_);
345
346 auto storage = create_thread_storage(LocalThreadMatStorage());
347
348 utils::maybe_parallel_for(local_boundary.size(), [&](int start, int end, int thread_id) {
349 LocalThreadMatStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
350
351 Eigen::MatrixXd pressure_vals, g_1, g_2, g_3, local_hessian;
352 ElementAssemblyValues vals;
353 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
354 Eigen::VectorXd weights;
355 Eigen::VectorXi global_primitive_ids;
356 for (int lb_id = start; lb_id < end; ++lb_id)
357 {
358 const auto &lb = local_boundary[lb_id];
359 const int e = lb.element_id();
360 const basis::ElementBases &gbs = gbases_[e];
361 const basis::ElementBases &bs = bases_[e];
362
363 for (int i = 0; i < lb.size(); ++i)
364 {
365 const int primitive_global_id = lb.global_primitive_id(i);
366 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
367
368 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
369 std::vector<Eigen::VectorXd> param_chain_rule;
370 if (mesh_.is_volume())
371 {
372 weights /= 2 * mesh_.tri_area(primitive_global_id);
373 g_1.setZero(normals.rows(), normals.cols());
374 g_2.setZero(normals.rows(), normals.cols());
375 g_3.setZero(normals.rows(), normals.cols());
376
377 auto endpoints = utils::BoundarySampler::tet_local_node_coordinates_from_face(lb[i]);
378 param_chain_rule = {(endpoints.row(0) - endpoints.row(1)).transpose(), (endpoints.row(0) - endpoints.row(2)).transpose()};
379 if (lb[i] == 0)
380 param_chain_rule[0] *= -1;
381 }
382 else
383 assert(false);
384
385 if (!has_samples)
386 continue;
387
388 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
389
390 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
391 for (int n = 0; n < vals.jac_it.size(); ++n)
392 {
393 trafo = vals.jac_it[n].inverse();
394
395 if (displacement.size() > 0)
396 {
397 assert(size_ == 3);
398 deform_mat.resize(size_, size_);
399 deform_mat.setZero();
400 for (const auto &b : vals.basis_values)
401 {
402 for (const auto &g : b.global)
403 {
404 for (int d = 0; d < size_; ++d)
405 {
406 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
407 }
408 }
409 }
410
411 trafo += deform_mat;
412 }
413
414 normals.row(n) = normals.row(n) * trafo.inverse();
415 normals.row(n).normalize();
416
417 if (mesh_.is_volume())
418 {
419 Eigen::Vector3d g1, g2, g3;
420 auto endpoints = utils::BoundarySampler::tet_local_node_coordinates_from_face(lb[i]);
421 g1 = trafo * (endpoints.row(0) - endpoints.row(1)).transpose();
422 g2 = trafo * (endpoints.row(0) - endpoints.row(2)).transpose();
423 if (lb[i] == 0)
424 g1 *= -1;
425 g3 = g1.cross(g2);
426
427 g_1.row(n) = g1.transpose();
428 g_2.row(n) = g2.transpose();
429 g_3.row(n) = g3.transpose();
430 }
431 else
432 assert(false);
433 }
434
435 if (multiply_pressure)
436 problem_.pressure_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, pressure_vals);
437 else
438 pressure_vals = Eigen::MatrixXd::Ones(weights.size(), 1);
439
440 local_hessian.setZero(vals.basis_values.size() * size_, vals.basis_values.size() * size_);
441
442 for (long p = 0; p < weights.size(); ++p)
443 {
444 Eigen::Vector3d g_up_1, g_up_2;
445 g_dual(g_1.row(p).transpose(), g_2.row(p).transpose(), g_up_1, g_up_2);
446 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)};
447
448 for (long ni = 0; ni < nodes.size(); ++ni)
449 {
450 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
451 for (int di = 0; di < size_; ++di)
452 {
453 const int gi_index = vi.global[0].index * size_ + di;
454 const bool is_dof_i_dirichlet = std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), gi_index) != dirichlet_nodes.end();
455 if (is_dof_i_dirichlet)
456 continue;
457
458 Eigen::MatrixXd grad_phi_i;
459 grad_phi_i.setZero(3, 3);
460 grad_phi_i.row(di) = vi.grad.row(p);
461
462 for (long nj = 0; nj < nodes.size(); ++nj)
463 {
464 const AssemblyValues &vj = vals.basis_values[nodes(nj)];
465 for (int dj = 0; dj < size_; ++dj)
466 {
467 const int gj_index = vj.global[0].index * size_ + dj;
468 const bool is_dof_j_dirichlet = std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), gj_index) != dirichlet_nodes.end();
469 if (is_dof_j_dirichlet)
470 continue;
471
472 Eigen::MatrixXd grad_phi_j;
473 grad_phi_j.setZero(3, 3);
474 grad_phi_j.row(dj) = vj.grad.row(p);
475
476 double value = 0;
477 for (int alpha = 0; alpha < size_ - 1; ++alpha)
478 {
479 auto a = vj.val(p) * g_3_wedge_g_up[alpha].row(di) * (grad_phi_i * param_chain_rule[alpha]);
480 auto b = (g_3_wedge_g_up[alpha] * (grad_phi_j * param_chain_rule[alpha])).row(di) * vi.val(p);
481 value += pressure_vals(p) * (a(0) + b(0)) * weights(p);
482 }
483 local_hessian(nodes(ni) * size_ + di, nodes(nj) * size_ + dj) += value;
484 }
485 }
486 }
487 }
488 }
489
490 for (long ni = 0; ni < nodes.size(); ++ni)
491 {
492 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
493 for (int di = 0; di < size_; ++di)
494 {
495 const int gi_index = vi.global[0].index * size_ + di;
496 for (long nj = 0; nj < nodes.size(); ++nj)
497 {
498 const AssemblyValues &vj = vals.basis_values[nodes(nj)];
499 for (int dj = 0; dj < size_; ++dj)
500 {
501 const int gj_index = vj.global[0].index * size_ + dj;
502
503 local_storage.entries.push_back(Eigen::Triplet<double>(gi_index, gj_index, local_hessian(nodes(ni) * size_ + di, nodes(nj) * size_ + dj)));
504 }
505 }
506 }
507 }
508 }
509 }
510 });
511
512 std::vector<Eigen::Triplet<double>> all_entries;
513 // Serially merge local storages
514 for (LocalThreadMatStorage &local_storage : storage)
515 {
516 all_entries.insert(all_entries.end(), local_storage.entries.begin(), local_storage.entries.end());
517 }
518
519 hess.setFromTriplets(all_entries.begin(), all_entries.end());
520 }
521
522 void PressureAssembler::compute_hess_volume_2d(
523 const Eigen::MatrixXd &displacement,
524 const std::vector<mesh::LocalBoundary> &local_boundary,
525 const std::vector<int> &dirichlet_nodes,
526 const int resolution,
527 StiffnessMatrix &hess,
528 const double t,
529 const bool multiply_pressure) const
530 {
531 hess.resize(n_basis_ * size_, n_basis_ * size_);
532
533 auto storage = create_thread_storage(LocalThreadMatStorage());
534
535 utils::maybe_parallel_for(local_boundary.size(), [&](int start, int end, int thread_id) {
536 LocalThreadMatStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
537
538 Eigen::MatrixXd pressure_vals, g_3_grad, local_hessian;
539 ElementAssemblyValues vals;
540 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
541 Eigen::VectorXd weights;
542 Eigen::VectorXi global_primitive_ids;
543 for (int lb_id = start; lb_id < end; ++lb_id)
544 {
545 const auto &lb = local_boundary[lb_id];
546 const int e = lb.element_id();
547 const basis::ElementBases &gbs = gbases_[e];
548 const basis::ElementBases &bs = bases_[e];
549
550 for (int i = 0; i < lb.size(); ++i)
551 {
552 const int primitive_global_id = lb.global_primitive_id(i);
553 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
554
555 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
556 std::vector<Eigen::VectorXd> param_chain_rule;
557
558 if (!has_samples)
559 continue;
560
561 if (mesh_.is_volume())
562 assert(false);
563 else
564 weights /= mesh_.edge_length(primitive_global_id);
565
566 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
567
568 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
569 g_3_grad.setZero(normals.rows(), size_ * vals.basis_values.size() * size_);
570
571 for (int n = 0; n < vals.jac_it.size(); ++n)
572 {
573 trafo = vals.jac_it[n].inverse();
574
575 if (displacement.size() > 0)
576 {
577 assert(size_ == 2);
578 deform_mat.resize(size_, size_);
579 deform_mat.setZero();
580 for (const auto &b : vals.basis_values)
581 {
582 for (const auto &g : b.global)
583 {
584 for (int d = 0; d < size_; ++d)
585 {
586 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
587 }
588 }
589 }
590
591 trafo += deform_mat;
592 }
593
594 normals.row(n) = normals.row(n) * trafo.inverse();
595 normals.row(n).normalize();
596
597 if (mesh_.is_volume())
598 assert(false);
599 else
600 {
601 Eigen::MatrixXd g_3_grad_local(size_, vals.basis_values.size() * size_);
602 g_3_grad_local.setZero();
603 auto endpoints = utils::BoundarySampler::tri_local_node_coordinates_from_edge(lb[i]);
604 Eigen::VectorXd reference_normal = (endpoints.row(0) - endpoints.row(1)).transpose();
605 int b_count = 0;
606 for (const auto &b : vals.basis_values)
607 {
608 for (const auto &g : b.global)
609 for (int d = 0; d < size_; ++d)
610 for (int k = 0; k < size_; ++k)
611 g_3_grad_local(1 - d, b_count * size_ + d) += b.grad(n, k) * reference_normal(k);
612 ++b_count;
613 }
614 g_3_grad_local.row(0) *= -1;
615
616 for (int k = 0; k < size_; ++k)
617 for (int l = 0; l < vals.basis_values.size() * size_; ++l)
618 g_3_grad(n, k * (vals.basis_values.size() * size_) + l) = g_3_grad_local(k, l);
619 }
620 }
621
622 if (multiply_pressure)
623 problem_.pressure_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, pressure_vals);
624 else
625 pressure_vals = Eigen::MatrixXd::Ones(weights.size(), 1);
626
627 local_hessian.setZero(vals.basis_values.size() * size_, vals.basis_values.size() * size_);
628
629 for (long p = 0; p < weights.size(); ++p)
630 {
631 for (long ni = 0; ni < nodes.size(); ++ni)
632 {
633 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
634 for (int di = 0; di < size_; ++di)
635 {
636 const int gi_index = vi.global[0].index * size_ + di;
637 const bool is_dof_i_dirichlet = std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), gi_index) != dirichlet_nodes.end();
638 if (is_dof_i_dirichlet)
639 continue;
640
641 for (long nj = 0; nj < nodes.size(); ++nj)
642 {
643 const AssemblyValues &vj = vals.basis_values[nodes(nj)];
644 for (int dj = 0; dj < size_; ++dj)
645 {
646 const int gj_index = vj.global[0].index * size_ + dj;
647 const bool is_dof_j_dirichlet = std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), gj_index) != dirichlet_nodes.end();
648 if (is_dof_j_dirichlet)
649 continue;
650
651 double value = pressure_vals(p) * g_3_grad(p, (di * vals.basis_values.size() * size_) + nodes(nj) * size_ + dj) * vi.val(p) * weights(p);
652 local_hessian(nodes(ni) * size_ + di, nodes(nj) * size_ + dj) += value;
653 }
654 }
655 }
656 }
657 }
658
659 for (long ni = 0; ni < nodes.size(); ++ni)
660 {
661 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
662 for (int di = 0; di < size_; ++di)
663 {
664 const int gi_index = vi.global[0].index * size_ + di;
665 for (long nj = 0; nj < nodes.size(); ++nj)
666 {
667 const AssemblyValues &vj = vals.basis_values[nodes(nj)];
668 for (int dj = 0; dj < size_; ++dj)
669 {
670 const int gj_index = vj.global[0].index * size_ + dj;
671
672 local_storage.entries.push_back(Eigen::Triplet<double>(gi_index, gj_index, local_hessian(nodes(ni) * size_ + di, nodes(nj) * size_ + dj)));
673 }
674 }
675 }
676 }
677 }
678 }
679 });
680
681 std::vector<Eigen::Triplet<double>> all_entries;
682 // Serially merge local storages
683 for (LocalThreadMatStorage &local_storage : storage)
684 {
685 all_entries.insert(all_entries.end(), local_storage.entries.begin(), local_storage.entries.end());
686 }
687
688 hess.setFromTriplets(all_entries.begin(), all_entries.end());
689 }
690
691 bool PressureAssembler::is_closed_or_boundary_fixed(
692 const std::vector<mesh::LocalBoundary> &local_boundary,
693 const std::vector<int> &dirichlet_nodes) const
694 {
695 if (local_boundary.size() == 0)
696 return true;
697
698 auto storage = utils::create_thread_storage(LocalThreadPrimitiveStorage());
699
700 utils::maybe_parallel_for(local_boundary.size(), [&](int start, int end, int thread_id) {
701 LocalThreadPrimitiveStorage &local_storage = utils::get_local_thread_storage(storage, thread_id);
702
703 ElementAssemblyValues vals;
704 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
705 Eigen::VectorXd weights;
706 Eigen::VectorXi global_primitive_ids;
707 for (int lb_id = start; lb_id < end; ++lb_id)
708 {
709 const auto &lb = local_boundary[lb_id];
710 const int e = lb.element_id();
711 const basis::ElementBases &bs = bases_[e];
712 const basis::ElementBases &gbs = gbases_[e];
713
714 for (int i = 0; i < lb.size(); ++i)
715 {
716 const int primitive_global_id = lb.global_primitive_id(i);
717 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
718
719 points.resize(0, size_);
720 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
721
722 if (!((lb.type() == BoundaryType::TRI_LINE) || (lb.type() == BoundaryType::TRI)))
723 {
724 logger().warn("Detected element not triangle or tet, cannot check if pressure boundary is closed");
725 return;
726 }
727
728 Eigen::VectorXi face(size_);
729 // Assuming the geometric nodes are listed in the beginning
730 for (int n = 0; n < size_; ++n)
731 {
732 const AssemblyValues &v = vals.basis_values[nodes(n)];
733 assert(v.global.size() == 1);
734 face(n) = v.global[0].index;
735 }
736 local_storage.faces.push_back(face);
737 }
738 }
739 });
740
741 Eigen::MatrixXi F(0, size_);
742 int count = 0;
743 for (const LocalThreadPrimitiveStorage &local_storage : storage)
744 {
745 for (int i = 0; i < local_storage.faces.size(); ++i)
746 {
747 F.conservativeResize(count + 1, size_);
748 F.row(count) = local_storage.faces[i];
749 ++count;
750 }
751 }
752
753 std::vector<int> L;
754 if (size_ == 2)
755 boundary_loop_2d(F, L);
756 else
757 igl::boundary_loop(F, L);
758
759 for (const int &l : L)
760 if (std::find(dirichlet_nodes.begin(), dirichlet_nodes.end(), l * size_) == dirichlet_nodes.end())
761 return false;
762
763 return true;
764 }
765
766 PressureAssembler::PressureAssembler(const Assembler &assembler, const Mesh &mesh, const Obstacle &obstacle,
767 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
768 const std::unordered_map<int, std::vector<mesh::LocalBoundary>> &local_pressure_cavity,
769 const std::vector<int> &dirichlet_nodes,
770 const std::vector<int> &primitive_to_nodes,
771 const std::vector<int> &node_to_primitives,
772 const int n_basis, const int size,
773 const std::vector<basis::ElementBases> &bases, const std::vector<basis::ElementBases> &gbases,
774 const Problem &problem)
775 : assembler_(assembler),
776 mesh_(mesh),
777 obstacle_(obstacle),
778 n_basis_(n_basis),
779 size_(size),
780 bases_(bases),
781 gbases_(gbases),
782 problem_(problem),
783 primitive_to_nodes_(primitive_to_nodes),
784 node_to_primitives_(node_to_primitives)
785 {
786 for (const auto &v : local_pressure_cavity)
787 starting_volumes_[v.first] = compute_volume(Eigen::VectorXd(), v.second, 5, 0, false);
788
789 if (!is_closed_or_boundary_fixed(local_pressure_boundary, dirichlet_nodes))
790 logger().error("Pressure boundary condition must be applied to a closed volume or have dirichlet fixed boundary.");
791
792 for (const auto &b : local_pressure_cavity)
793 if (!is_closed_or_boundary_fixed(b.second, dirichlet_nodes))
794 logger().error("Pressure cavity boundary condition must be applied to a closed volume or have dirichlet fixed boundary.");
795
796 cavity_thermodynamics_ = std::make_unique<AdiabaticProcess>();
797 // cavity_thermodynamics_ = std::make_unique<IsothermalProcess>();
798 }
799
800 double PressureAssembler::compute_cavity_energy(
801 const Eigen::MatrixXd &displacement,
802 const std::unordered_map<int, std::vector<mesh::LocalBoundary>> &local_pressure_cavity,
803 const int resolution,
804 const double t) const
805 {
806 double energy_ = 0;
807 for (const auto &v : local_pressure_cavity)
808 {
809 const double start_pressure = problem_.pressure_cavity_bc(v.first, t);
810 const double start_volume = starting_volumes_.at(v.first);
811 const double curr_volume = compute_volume(displacement, v.second, resolution, t, false);
812
813 energy_ += cavity_thermodynamics_->energy(-start_pressure, -start_volume, -curr_volume);
814 }
815 return energy_;
816 }
817
818 void PressureAssembler::compute_cavity_energy_grad(
819 const Eigen::MatrixXd &displacement,
820 const std::unordered_map<int, std::vector<mesh::LocalBoundary>> &local_pressure_cavity,
821 const std::vector<int> &dirichlet_nodes,
822 const int resolution,
823 const double t,
824 Eigen::VectorXd &grad) const
825 {
826 grad.setZero(displacement.size());
827 for (const auto &v : local_pressure_cavity)
828 {
829 double start_pressure = problem_.pressure_cavity_bc(v.first, t);
830 double start_volume = starting_volumes_.at(v.first);
831 double curr_volume = compute_volume(displacement, v.second, resolution, t, false);
832 Eigen::VectorXd g;
833
834 compute_grad_volume(displacement, v.second, dirichlet_nodes, resolution, g, t, false);
835
836 double p = -cavity_thermodynamics_->pressure(
837 -start_pressure, -start_volume, -curr_volume);
838
839 grad += p * g;
840 }
841 }
842
843 void PressureAssembler::compute_cavity_energy_hess(
844 const Eigen::MatrixXd &displacement,
845 const std::unordered_map<int, std::vector<mesh::LocalBoundary>> &local_pressure_cavity,
846 const std::vector<int> &dirichlet_nodes,
847 const int resolution,
848 const double t,
849 const bool project_to_psd,
850 StiffnessMatrix &hess) const
851 {
852 if (project_to_psd && local_pressure_cavity.size() > 0)
853 {
854 log_and_throw_error("Cannot project caivity pressure to PSD!");
855 }
856
857 hess.resize(displacement.size(), displacement.size());
858 for (const auto &v : local_pressure_cavity)
859 {
860 double start_pressure = problem_.pressure_cavity_bc(v.first, t);
861 double start_volume = starting_volumes_.at(v.first);
862 double curr_volume = compute_volume(displacement, v.second, resolution, t, false);
863 Eigen::VectorXd g;
864 StiffnessMatrix h;
865
866 compute_grad_volume(displacement, v.second, dirichlet_nodes, resolution, g, t, false);
867 if (size_ == 2)
868 compute_hess_volume_2d(displacement, v.second, dirichlet_nodes, resolution, h, t, false);
869 else if (size_ == 3)
870 compute_hess_volume_3d(displacement, v.second, dirichlet_nodes, resolution, h, t, false);
871
872 double p = -cavity_thermodynamics_->pressure(
873 -start_pressure, -start_volume, -curr_volume);
874 double dp_dv = cavity_thermodynamics_->first_derivative(
875 -start_pressure, -start_volume, -curr_volume);
876
877 hess += p * h;
878 hess += dp_dv * (g * g.transpose()).sparseView();
879 }
880 }
881
882 double PressureAssembler::compute_energy(
883 const Eigen::MatrixXd &displacement,
884 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
885 const int resolution,
886 const double t) const
887 {
888 double vol = compute_volume(displacement, local_pressure_boundary, resolution, t, false);
889 return compute_volume(displacement, local_pressure_boundary, resolution, t, true);
890 }
891
892 void PressureAssembler::compute_energy_grad(
893 const Eigen::MatrixXd &displacement,
894 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
895 const std::vector<int> &dirichlet_nodes,
896 const int resolution,
897 const double t,
898 Eigen::VectorXd &grad) const
899 {
900 compute_grad_volume(displacement, local_pressure_boundary, dirichlet_nodes, resolution, grad, t, true);
901 }
902
903 void PressureAssembler::compute_energy_hess(
904 const Eigen::MatrixXd &displacement,
905 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
906 const std::vector<int> &dirichlet_nodes,
907 const int resolution,
908 const double t,
909 const bool project_to_psd,
910 StiffnessMatrix &hess) const
911 {
912 if (project_to_psd && local_pressure_boundary.size() > 0)
913 {
914 log_and_throw_error("Cannot project pressure to PSD!");
915 }
916
917 if (size_ == 2)
918 compute_hess_volume_2d(displacement, local_pressure_boundary, dirichlet_nodes, resolution, hess, t, true);
919 else if (size_ == 3)
920 compute_hess_volume_3d(displacement, local_pressure_boundary, dirichlet_nodes, resolution, hess, t, true);
921 }
922
923 void PressureAssembler::compute_force_jacobian(
924 const Eigen::MatrixXd &displacement,
925 const std::vector<mesh::LocalBoundary> &local_pressure_boundary,
926 const std::vector<int> &dirichlet_nodes,
927 const int resolution,
928 const double t,
929 const int n_vertices,
930 StiffnessMatrix &hess) const
931 {
932 compute_energy_hess(displacement, local_pressure_boundary, dirichlet_nodes, resolution, t, false, hess);
933 }
934
935 } // namespace assembler
936} // 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:882
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:691
polyfem::assembler::PressureAssembler::cavity_thermodynamics_
std::unique_ptr< ThermodynamicProcess > cavity_thermodynamics_
Definition PressureAssembler.hpp:190
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:335
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:522
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:903
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:892
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:923
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:766
polyfem::assembler::PressureAssembler::starting_volumes_
std::unordered_map< int, double > starting_volumes_
Definition PressureAssembler.hpp:189
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:818
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:843
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:800
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:20