PolyFEM
Loading...
Searching...
No Matches
RhsAssembler.cpp
Go to the documentation of this file.
1#include "RhsAssembler.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
10namespace polyfem
11{
12 using namespace polysolve;
13 using namespace mesh;
14 using namespace quadrature;
15 using namespace utils;
16
17 namespace assembler
18 {
19 namespace
20 {
21 class LocalThreadScalarStorage
22 {
23 public:
24 double val;
25 ElementAssemblyValues vals;
26
27 LocalThreadScalarStorage()
28 {
29 val = 0;
30 }
31 };
32 } // namespace
33
34 RhsAssembler::RhsAssembler(const Assembler &assembler, const Mesh &mesh, const Obstacle &obstacle,
35 const std::vector<int> &dirichlet_nodes, const std::vector<int> &neumann_nodes,
36 const std::vector<RowVectorNd> &dirichlet_nodes_position, const std::vector<RowVectorNd> &neumann_nodes_position,
37 const int n_basis, const int size,
38 const std::vector<basis::ElementBases> &bases, const std::vector<basis::ElementBases> &gbases, const AssemblyValsCache &ass_vals_cache,
39 const Problem &problem,
40 const std::string bc_method,
41 const json &solver_params)
42 : assembler_(assembler),
43 mesh_(mesh),
44 obstacle_(obstacle),
45 n_basis_(n_basis),
46 size_(size),
47 bases_(bases),
48 gbases_(gbases),
49 ass_vals_cache_(ass_vals_cache),
50 problem_(problem),
51 bc_method_(bc_method),
52 solver_params_(solver_params),
53 dirichlet_nodes_(dirichlet_nodes),
54 dirichlet_nodes_position_(dirichlet_nodes_position),
55 neumann_nodes_(neumann_nodes),
56 neumann_nodes_position_(neumann_nodes_position)
57 {
58 assert(ass_vals_cache_.is_mass());
59 }
60
61 void RhsAssembler::assemble(const Density &density, Eigen::MatrixXd &rhs, const double t) const
62 {
63 // set size of rhs to the number of basis functions * the dimension of the problem
64 rhs = Eigen::MatrixXd::Zero(n_basis_ * size_, 1);
65 if (!problem_.is_rhs_zero())
66 {
67 Eigen::MatrixXd rhs_fun;
68
69 const int n_elements = int(bases_.size());
70 ElementAssemblyValues vals;
71 for (int e = 0; e < n_elements; ++e)
72 {
73 // vals.compute(e, mesh_.is_volume(), bases_[e], gbases_[e]);
74
75 // compute geometric mapping
76 // evaluate and store basis functions/their gradients at quadrature points
77 ass_vals_cache_.compute(e, mesh_.is_volume(), bases_[e], gbases_[e], vals);
78
79 const Quadrature &quadrature = vals.quadrature;
80
81 // compute rhs values in physical space
82 problem_.rhs(assembler_, vals.val, t, rhs_fun);
83
84 for (int d = 0; d < size_; ++d)
85 {
86 // rhs_fun.col(d) = rhs_fun.col(d).array() * vals.det.array() * quadrature.weights.array();
87 for (int q = 0; q < quadrature.weights.size(); ++q)
88 {
89 // const double rho = problem_.is_time_dependent() ? density(vals.quadrature.points.row(q), vals.val.row(q), vals.element_id) : 1;
90 const double rho = density(vals.quadrature.points.row(q), vals.val.row(q), t, vals.element_id);
91 // prepare for integration by weighing rhs by determinant and quadrature weights
92 rhs_fun(q, d) *= vals.det(q) * quadrature.weights(q) * rho;
93 }
94 }
95
96 const int n_loc_bases_ = int(vals.basis_values.size());
97 for (int i = 0; i < n_loc_bases_; ++i)
98 {
99 const AssemblyValues &v = vals.basis_values[i];
100
101 for (int d = 0; d < size_; ++d)
102 {
103 // integrate rhs function times the given local basis
104 const double rhs_value = (rhs_fun.col(d).array() * v.val.array()).sum();
105 for (std::size_t ii = 0; ii < v.global.size(); ++ii)
106 // add local contribution to the global rhs vector (with some weight for non-conforming bases)
107 rhs(v.global[ii].index * size_ + d) += rhs_value * v.global[ii].val;
108 }
109 }
110 }
111 }
112 }
113
114 void RhsAssembler::initial_solution(Eigen::MatrixXd &sol) const
115 {
116 time_bc([&](const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) {
117 problem_.initial_solution(mesh, global_ids, pts, val);
118 },
119 sol);
120 }
121
122 void RhsAssembler::initial_velocity(Eigen::MatrixXd &sol) const
123 {
124 time_bc([&](const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) {
125 problem_.initial_velocity(mesh, global_ids, pts, val);
126 },
127 sol);
128 }
129
130 void RhsAssembler::initial_acceleration(Eigen::MatrixXd &sol) const
131 {
132 time_bc([&](const Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) {
133 problem_.initial_acceleration(mesh, global_ids, pts, val);
134 },
135 sol);
136 }
137
138 void RhsAssembler::time_bc(const std::function<void(const Mesh &, const Eigen::MatrixXi &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &fun, Eigen::MatrixXd &sol) const
139 {
140 sol = Eigen::MatrixXd::Zero(n_basis_ * size_, 1);
141 Eigen::MatrixXd loc_sol;
142
143 const int n_elements = int(bases_.size());
144 ElementAssemblyValues vals;
145 Eigen::MatrixXi ids;
146
147 if (bc_method_ == "sample")
148 {
149 for (int e = 0; e < n_elements; ++e)
150 {
151 const basis::ElementBases &bs = bases_[e];
152 // vals.compute(e, mesh_.is_volume(), bases_[e], gbases_[e]);
153 ass_vals_cache_.compute(e, mesh_.is_volume(), bases_[e], gbases_[e], vals);
154 ids.resize(1, 1);
155 ids.setConstant(e);
156
157 for (long i = 0; i < bs.bases.size(); ++i)
158 {
159 const auto &b = bs.bases[i];
160 const auto &glob = b.global();
161 // assert(glob.size() == 1);
162 for (size_t ii = 0; ii < glob.size(); ++ii)
163 {
164 fun(mesh_, ids, glob[ii].node, loc_sol);
165
166 for (int d = 0; d < size_; ++d)
167 {
168 sol(glob[ii].index * size_ + d) = loc_sol(d) * glob[ii].val;
169 }
170 }
171 }
172 }
173 }
174 else
175 {
176
177 for (int e = 0; e < n_elements; ++e)
178 {
179 // vals.compute(e, mesh_.is_volume(), bases_[e], gbases_[e]);
180 ass_vals_cache_.compute(e, mesh_.is_volume(), bases_[e], gbases_[e], vals);
181 ids.resize(vals.val.rows(), 1);
182 ids.setConstant(e);
183
184 const Quadrature &quadrature = vals.quadrature;
185 // problem_.initial_solution(vals.val, loc_sol);
186 fun(mesh_, ids, vals.val, loc_sol);
187
188 for (int d = 0; d < size_; ++d)
189 loc_sol.col(d) = loc_sol.col(d).array() * vals.det.array() * quadrature.weights.array();
190
191 const int n_loc_bases_ = int(vals.basis_values.size());
192 for (int i = 0; i < n_loc_bases_; ++i)
193 {
194 const AssemblyValues &v = vals.basis_values[i];
195
196 for (int d = 0; d < size_; ++d)
197 {
198 const double sol_value = (loc_sol.col(d).array() * v.val.array()).sum();
199 for (std::size_t ii = 0; ii < v.global.size(); ++ii)
200 sol(v.global[ii].index * size_ + d) += sol_value * v.global[ii].val;
201 }
202 }
203 }
204
205 Eigen::MatrixXd b = sol;
206 sol.setZero();
207
208 const double mmin = b.minCoeff();
209 const double mmax = b.maxCoeff();
210
211 if (fabs(mmin) > 1e-8 || fabs(mmax) > 1e-8)
212 {
213 assembler::Mass mass_mat_assembler;
214 mass_mat_assembler.set_size(assembler_.size());
215 mass_mat_assembler.add_multimaterial(0, json({}), Units());
216 StiffnessMatrix mass;
217 const int n_fe_basis = n_basis_ - obstacle_.n_vertices();
218 mass_mat_assembler.assemble(size_ == 3, n_fe_basis, bases_, gbases_, ass_vals_cache_, 0, mass, true);
219 assert(mass.rows() == n_basis_ * size_ - obstacle_.ndof() && mass.cols() == n_basis_ * size_ - obstacle_.ndof());
220
221 auto solver = linear::Solver::create(solver_params_, logger());
222 logger().info("Solve RHS using {} linear solver", solver->name());
223 solver->analyze_pattern(mass, mass.rows());
224 solver->factorize(mass);
225
226 for (long i = 0; i < b.cols(); ++i)
227 {
228 solver->solve(b.block(0, i, mass.rows(), 1), sol.block(0, i, mass.rows(), 1));
229 }
230 logger().trace("mass matrix error {}", (mass * sol - b).norm());
231 }
232 }
233 }
234
235 void RhsAssembler::lsq_bc(const std::function<void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &df,
236 const std::vector<LocalBoundary> &local_boundary, const std::vector<int> &bounday_nodes, const int resolution, Eigen::MatrixXd &rhs) const
237 {
238 const int n_el = int(bases_.size());
239
240 Eigen::MatrixXd uv, samples, gtmp, rhs_fun;
241 Eigen::VectorXi global_primitive_ids;
242
243 const int actual_dim = problem_.is_scalar() ? 1 : mesh_.dimension();
244
245 Eigen::Matrix<bool, Eigen::Dynamic, 1> is_boundary(n_basis_);
246 is_boundary.setConstant(false);
247 int skipped_count = 0;
248 for (int b : bounday_nodes)
249 {
250 int bindex = b / actual_dim;
251
252 if (bindex < is_boundary.size())
253 is_boundary[bindex] = true;
254 else
255 skipped_count++;
256 }
257 assert(skipped_count <= 1);
258
259 for (int d = 0; d < size_; ++d)
260 {
261 int index = 0;
262 std::vector<int> indices;
263 indices.reserve(n_el * 10);
264 std::vector<int> tags;
265 tags.reserve(n_el * 10);
266
267 long total_size = 0;
268
269 Eigen::VectorXi global_index_to_col(n_basis_);
270 global_index_to_col.setConstant(-1);
271
272 std::vector<AssemblyValues> tmp_val;
273
274 for (const auto &lb : local_boundary)
275 {
276 const int e = lb.element_id();
277 bool has_samples = utils::BoundarySampler::sample_boundary(lb, resolution, mesh_, true, uv, samples, global_primitive_ids);
278
279 if (!has_samples)
280 continue;
281
282 const basis::ElementBases &bs = bases_[e];
283 bs.evaluate_bases(samples, tmp_val);
284 const int n_local_bases = int(bs.bases.size());
285 assert(global_primitive_ids.size() == samples.rows());
286
287 for (int s = 0; s < samples.rows(); ++s)
288 {
289 const int tag = mesh_.get_boundary_id(global_primitive_ids(s));
290 if (!problem_.all_dimensions_dirichlet() && !problem_.is_dimension_dirichet(tag, d))
291 continue;
292
293 total_size++;
294
295 for (int j = 0; j < n_local_bases; ++j)
296 {
297 const basis::Basis &b = bs.bases[j];
298 const double tmp = tmp_val[j].val(s);
299
300 if (fabs(tmp) < 1e-10)
301 continue;
302
303 for (std::size_t ii = 0; ii < b.global().size(); ++ii)
304 {
305 // pt found
306 if (is_boundary[b.global()[ii].index])
307 {
308 if (global_index_to_col(b.global()[ii].index) == -1)
309 {
310 global_index_to_col(b.global()[ii].index) = index++;
311 indices.push_back(b.global()[ii].index);
312 tags.push_back(tag);
313 assert(indices.size() == size_t(index));
314 }
315 }
316 }
317 }
318 }
319 }
320
321 Eigen::MatrixXd global_rhs = Eigen::MatrixXd::Zero(total_size, 1);
322
323 const long buffer_size = total_size * long(indices.size());
324 std::vector<Eigen::Triplet<double>> entries, entries_t;
325
326 index = 0;
327
328 int global_counter = 0;
329 Eigen::MatrixXd mapped;
330
331 for (const auto &lb : local_boundary)
332 {
333 const int e = lb.element_id();
334 bool has_samples = utils::BoundarySampler::sample_boundary(lb, resolution, mesh_, false, uv, samples, global_primitive_ids);
335
336 if (!has_samples)
337 continue;
338
339 const basis::ElementBases &bs = bases_[e];
340 const basis::ElementBases &gbs = gbases_[e];
341 const int n_local_bases = int(bs.bases.size());
342
343 gbs.eval_geom_mapping(samples, mapped);
344
345 bs.evaluate_bases(samples, tmp_val);
346 df(global_primitive_ids, uv, mapped, rhs_fun);
347
348 for (int s = 0; s < samples.rows(); ++s)
349 {
350 const int tag = mesh_.get_boundary_id(global_primitive_ids(s));
351 if (!problem_.all_dimensions_dirichlet() && !problem_.is_dimension_dirichet(tag, d))
352 continue;
353
354 for (int j = 0; j < n_local_bases; ++j)
355 {
356 const basis::Basis &b = bs.bases[j];
357 const double tmp = tmp_val[j].val(s);
358
359 for (std::size_t ii = 0; ii < b.global().size(); ++ii)
360 {
361 auto item = global_index_to_col(b.global()[ii].index);
362 if (item != -1)
363 {
364 entries.push_back(Eigen::Triplet<double>(global_counter, item, tmp * b.global()[ii].val));
365 entries_t.push_back(Eigen::Triplet<double>(item, global_counter, tmp * b.global()[ii].val));
366 }
367 }
368 }
369
370 global_rhs(global_counter) = rhs_fun(s, d);
371 global_counter++;
372 }
373 }
374
375 assert(global_counter == total_size);
376
377 if (total_size > 0)
378 {
379 const double mmin = global_rhs.minCoeff();
380 const double mmax = global_rhs.maxCoeff();
381
382 if (fabs(mmin) < 1e-8 && fabs(mmax) < 1e-8)
383 {
384 for (size_t i = 0; i < indices.size(); ++i)
385 {
386 const int tag = tags[i];
387 if (problem_.all_dimensions_dirichlet() || problem_.is_dimension_dirichet(tag, d))
388 rhs(indices[i] * size_ + d) = 0;
389 }
390 }
391 else
392 {
393 StiffnessMatrix mat(int(total_size), int(indices.size()));
394 mat.setFromTriplets(entries.begin(), entries.end());
395
396 StiffnessMatrix mat_t(int(indices.size()), int(total_size));
397 mat_t.setFromTriplets(entries_t.begin(), entries_t.end());
398
399 StiffnessMatrix A = mat_t * mat;
400 Eigen::VectorXd b = mat_t * global_rhs;
401
402 Eigen::VectorXd coeffs(b.rows(), 1);
403 auto solver = linear::Solver::create(solver_params_, logger());
404 logger().info("Solve RHS using {} linear solver", solver->name());
405 solver->analyze_pattern(A, A.rows());
406 solver->factorize(A);
407 coeffs.setZero();
408 solver->solve(b, coeffs);
409
410 logger().trace("RHS solve error {}", (A * coeffs - b).norm());
411
412 for (long i = 0; i < coeffs.rows(); ++i)
413 {
414 const int tag = tags[i];
415 if (problem_.all_dimensions_dirichlet() || problem_.is_dimension_dirichet(tag, d))
416 rhs(indices[i] * size_ + d) = coeffs(i);
417 }
418 }
419 }
420 }
421 }
422
423 void RhsAssembler::sample_bc(const std::function<void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &df,
424 const std::vector<LocalBoundary> &local_boundary, const std::vector<int> &bounday_nodes, Eigen::MatrixXd &rhs) const
425 {
426 const int n_el = int(bases_.size());
427
428 Eigen::MatrixXd rhs_fun;
429 Eigen::VectorXi global_primitive_ids(1);
430 Eigen::MatrixXd nans(1, 1);
431 nans(0) = std::nan("");
432
433#ifndef NDEBUG
434 Eigen::Matrix<bool, Eigen::Dynamic, 1> is_boundary(n_basis_);
435 is_boundary.setConstant(false);
436
437 const int actual_dim = problem_.is_scalar() ? 1 : mesh_.dimension();
438
439 int skipped_count = 0;
440 for (int b : bounday_nodes)
441 {
442 int bindex = b / actual_dim;
443
444 if (bindex < is_boundary.size())
445 is_boundary[bindex] = true;
446 else
447 skipped_count++;
448 }
449 assert(skipped_count <= 1);
450#endif
451
452 for (const auto &lb : local_boundary)
453 {
454 const int e = lb.element_id();
455 const basis::ElementBases &bs = bases_[e];
456
457 for (int i = 0; i < lb.size(); ++i)
458 {
459 global_primitive_ids(0) = lb.global_primitive_id(i);
460 const auto nodes = bs.local_nodes_for_primitive(global_primitive_ids(0), mesh_);
461 assert(global_primitive_ids.size() == 1);
462 const int tag = mesh_.get_boundary_id(global_primitive_ids(0));
463
464 for (long n = 0; n < nodes.size(); ++n)
465 {
466 const auto &b = bs.bases[nodes(n)];
467 const auto &glob = b.global();
468
469 for (size_t ii = 0; ii < glob.size(); ++ii)
470 {
471 assert(is_boundary[glob[ii].index]);
472
473 // TODO, missing UV!!!!
474 df(global_primitive_ids, nans, glob[ii].node, rhs_fun);
475
476 for (int d = 0; d < size_; ++d)
477 {
478 if (problem_.all_dimensions_dirichlet() || problem_.is_dimension_dirichet(tag, d))
479 {
480 assert(problem_.all_dimensions_dirichlet() || std::find(bounday_nodes.begin(), bounday_nodes.end(), glob[ii].index * size_ + d) != bounday_nodes.end());
481 rhs(glob[ii].index * size_ + d) = rhs_fun(0, d);
482 }
483 }
484 }
485 }
486 }
487 }
488 }
489
490 void RhsAssembler::integrate_bc(const std::function<void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &df,
491 const std::vector<LocalBoundary> &local_boundary, const std::vector<int> &bounday_nodes, const int resolution, Eigen::MatrixXd &rhs) const
492 {
493 assert(false);
494 Eigen::MatrixXd uv, samples, rhs_fun, normals, mapped;
495 Eigen::VectorXd weights;
496
497 Eigen::VectorXi global_primitive_ids;
498 std::vector<AssemblyValues> tmp_val;
499
500 Eigen::Matrix<bool, Eigen::Dynamic, 1> is_boundary(n_basis_);
501 is_boundary.setConstant(false);
502
503 Eigen::MatrixXd areas(rhs.rows(), 1);
504 areas.setZero();
505
506 const int actual_dim = problem_.is_scalar() ? 1 : mesh_.dimension();
507
508 int skipped_count = 0;
509 for (int b : bounday_nodes)
510 {
511 rhs(b) = 0;
512 int bindex = b / actual_dim;
513
514 if (bindex < is_boundary.size())
515 is_boundary[bindex] = true;
516 else
517 skipped_count++;
518 }
519 assert(skipped_count <= 1);
520 ElementAssemblyValues vals;
521
522 for (const auto &lb : local_boundary)
523 {
524 const int e = lb.element_id();
525 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, false, uv, samples, normals, weights, global_primitive_ids);
526
527 if (!has_samples)
528 continue;
529
530 const basis::ElementBases &bs = bases_[e];
531 const basis::ElementBases &gbs = gbases_[e];
532
533 vals.compute(e, mesh_.is_volume(), samples, bs, gbs);
534
535 df(global_primitive_ids, uv, vals.val, rhs_fun);
536
537 for (int d = 0; d < size_; ++d)
538 rhs_fun.col(d) = rhs_fun.col(d).array() * weights.array();
539
540 for (int i = 0; i < lb.size(); ++i)
541 {
542 const int primitive_global_id = lb.global_primitive_id(i);
543 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
544
545 for (long n = 0; n < nodes.size(); ++n)
546 {
547 // const auto &b = bs.bases[nodes(n)];
548 const AssemblyValues &v = vals.basis_values[nodes(n)];
549 const double area = (weights.array() * v.val.array()).sum();
550 for (int d = 0; d < size_; ++d)
551 {
552 const double rhs_value = (rhs_fun.col(d).array() * v.val.array()).sum();
553
554 for (size_t g = 0; g < v.global.size(); ++g)
555 {
556 const int g_index = v.global[g].index * size_ + d;
557 if (problem_.all_dimensions_dirichlet() || std::find(bounday_nodes.begin(), bounday_nodes.end(), g_index) != bounday_nodes.end())
558 {
559 rhs(g_index) += rhs_value * v.global[g].val;
560 areas(g_index) += area * v.global[g].val;
561 }
562 }
563 }
564 }
565 }
566 }
567
568 for (int b : bounday_nodes)
569 {
570 assert(areas(b) != 0);
571 rhs(b) /= areas(b);
572 }
573 }
574
575 void RhsAssembler::set_bc(
576 const std::function<void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &df,
577 const std::function<void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &nf,
578 const std::vector<LocalBoundary> &local_boundary, const std::vector<int> &bounday_nodes,
579 const int resolution, const std::vector<LocalBoundary> &local_neumann_boundary,
580 const Eigen::MatrixXd &displacement, const double t,
581 Eigen::MatrixXd &rhs) const
582 {
583 if (bc_method_ == "sample")
584 sample_bc(df, local_boundary, bounday_nodes, rhs);
585 else if (bc_method_ == "integrate")
586 integrate_bc(df, local_boundary, bounday_nodes, resolution, rhs);
587 else
588 lsq_bc(df, local_boundary, bounday_nodes, resolution, rhs);
589
590 if (bounday_nodes.size() > 0)
591 {
592 Eigen::MatrixXd tmp_val;
593 for (int n = 0; n < dirichlet_nodes_.size(); ++n)
594 {
595 const auto &n_id = dirichlet_nodes_[n];
596 const auto &pt = dirichlet_nodes_position_[n];
597
598 const int tag = mesh_.get_node_id(n_id);
599 problem_.dirichlet_nodal_value(mesh_, n_id, pt, t, tmp_val);
600 assert(tmp_val.size() == size_);
601
602 for (int d = 0; d < size_; ++d)
603 {
604 if (!problem_.is_nodal_dimension_dirichlet(n_id, tag, d))
605 continue;
606 const int g_index = n_id * size_ + d;
607 rhs(g_index) = tmp_val(d);
608 }
609 }
610 }
611
612 // Neumann
613 Eigen::MatrixXd uv, samples, gtmp, rhs_fun, deform_mat, trafo;
614 Eigen::VectorXi global_primitive_ids;
615 Eigen::MatrixXd points, normals;
616 Eigen::VectorXd weights;
617
618 ElementAssemblyValues vals;
619
620 for (const auto &lb : local_neumann_boundary)
621 {
622 const int e = lb.element_id();
623 const basis::ElementBases &gbs = gbases_[e];
624 const basis::ElementBases &bs = bases_[e];
625
626 for (int i = 0; i < lb.size(); ++i)
627 {
628 const int primitive_global_id = lb.global_primitive_id(i);
629
630 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
631 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
632
633 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
634
635 for (int n = 0; n < vals.jac_it.size(); ++n)
636 {
637 trafo = vals.jac_it[n].inverse();
638
639 if (displacement.size() > 0)
640 {
641 assert(size_ == 2 || size_ == 3);
642 deform_mat.resize(size_, size_);
643 deform_mat.setZero();
644 for (const auto &b : vals.basis_values)
645 {
646 for (const auto &g : b.global)
647 {
648 for (int d = 0; d < size_; ++d)
649 {
650 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
651 }
652 }
653 }
654
655 trafo += deform_mat;
656 }
657
658 normals.row(n) = normals.row(n) * trafo.inverse();
659 normals.row(n).normalize();
660 }
661
662 // problem_.neumann_bc(mesh_, global_primitive_ids, vals.val, t, rhs_fun);
663 nf(global_primitive_ids, uv, vals.val, normals, rhs_fun);
664
665 // UIState::ui_state().debug_data().add_points(vals.val, Eigen::RowVector3d(0,1,0));
666
667 for (int d = 0; d < size_; ++d)
668 rhs_fun.col(d) = rhs_fun.col(d).array() * weights.array();
669
670 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
671
672 for (long n = 0; n < nodes.size(); ++n)
673 {
674 // const auto &b = bs.bases[nodes(n)];
675 const AssemblyValues &v = vals.basis_values[nodes(n)];
676 for (int d = 0; d < size_; ++d)
677 {
678 const double rhs_value = (rhs_fun.col(d).array() * v.val.array()).sum();
679
680 for (size_t g = 0; g < v.global.size(); ++g)
681 {
682 const int g_index = v.global[g].index * size_ + d;
683 const bool is_neumann = std::find(bounday_nodes.begin(), bounday_nodes.end(), g_index) == bounday_nodes.end();
684
685 if (is_neumann)
686 {
687 rhs(g_index) += rhs_value * v.global[g].val;
688 }
689 }
690 }
691 }
692 }
693 }
694
695 // TODO add nodal neumann
696 }
697
698 void RhsAssembler::set_bc(const std::vector<LocalBoundary> &local_boundary, const std::vector<int> &bounday_nodes, const int resolution, const std::vector<LocalBoundary> &local_neumann_boundary, Eigen::MatrixXd &rhs, const Eigen::MatrixXd &displacement, const double t) const
699 {
700 set_bc(
701 [&](const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) {
702 problem_.dirichlet_bc(mesh_, global_ids, uv, pts, t, val);
703 },
704 [&](const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &normals, Eigen::MatrixXd &val) {
705 problem_.neumann_bc(mesh_, global_ids, uv, pts, normals, t, val);
706 },
707 local_boundary, bounday_nodes, resolution, local_neumann_boundary, displacement, t, rhs);
708
709 obstacle_.update_displacement(t, rhs);
710 }
711
712 void RhsAssembler::compute_energy_grad(const std::vector<LocalBoundary> &local_boundary, const std::vector<int> &bounday_nodes, const Density &density, const int resolution, const std::vector<LocalBoundary> &local_neumann_boundary, const Eigen::MatrixXd &final_rhs, const double t, Eigen::MatrixXd &rhs) const
713 {
714 if (problem_.is_constant_in_time())
715 {
716 rhs = final_rhs;
717 }
718 else
719 {
720 assemble(density, rhs, t);
721 rhs *= -1;
722
723 if (rhs.size() != final_rhs.size())
724 {
725 const int prev_size = rhs.size();
726 rhs.conservativeResize(final_rhs.size(), rhs.cols());
727 // Zero initial pressure
728 rhs.block(prev_size, 0, final_rhs.size() - prev_size, rhs.cols()).setZero();
729 rhs(rhs.size() - 1) = 0;
730 }
731
732 assert(rhs.size() == final_rhs.size());
733 }
734 }
735
736 double RhsAssembler::compute_energy(const Eigen::MatrixXd &displacement,
737 const Eigen::MatrixXd &displacement_prev,
738 const std::vector<LocalBoundary> &local_neumann_boundary,
739 const Density &density,
740 const int resolution,
741 const double t) const
742 {
743
744 double res = 0;
745
746 if (!problem_.is_rhs_zero())
747 {
748 auto storage = create_thread_storage(LocalThreadScalarStorage());
749 const int n_bases = int(bases_.size());
750
751 maybe_parallel_for(n_bases, [&](int start, int end, int thread_id) {
752 LocalThreadScalarStorage &local_storage = get_local_thread_storage(storage, thread_id);
753 VectorNd local_displacement(size_);
754 Eigen::MatrixXd forces;
755
756 for (int e = start; e < end; ++e)
757 {
758 ElementAssemblyValues &vals = local_storage.vals;
759 // vals.compute(e, mesh_.is_volume(), bases_[e], gbases_[e]);
760 ass_vals_cache_.compute(e, mesh_.is_volume(), bases_[e], gbases_[e], vals);
761
762 const Quadrature &quadrature = vals.quadrature;
763 const Eigen::VectorXd da = vals.det.array() * quadrature.weights.array();
764
765 problem_.rhs(assembler_, vals.val, t, forces);
766 assert(forces.rows() == da.size());
767 assert(forces.cols() == size_);
768
769 for (long p = 0; p < da.size(); ++p)
770 {
771 local_displacement.setZero();
772
773 for (size_t i = 0; i < vals.basis_values.size(); ++i)
774 {
775 const auto &bs = vals.basis_values[i];
776 assert(bs.val.size() == da.size());
777 const double b_val = bs.val(p);
778
779 for (int d = 0; d < size_; ++d)
780 {
781 for (std::size_t ii = 0; ii < bs.global.size(); ++ii)
782 {
783 local_displacement(d) += (bs.global[ii].val * b_val) * displacement(bs.global[ii].index * size_ + d);
784 }
785 }
786 }
787 // const double rho = problem_.is_time_dependent() ? density(vals.quadrature.points.row(p), vals.val.row(p), vals.element_id) : 1;
788 const double rho = density(vals.quadrature.points.row(p), vals.val.row(p), t, vals.element_id);
789
790 for (int d = 0; d < size_; ++d)
791 {
792 local_storage.val += forces(p, d) * local_displacement(d) * da(p) * rho;
793 // res += forces(p, d) * local_displacement(d) * da(p);
794 }
795 }
796 }
797 });
798
799 // Serially merge local storages
800 for (const LocalThreadScalarStorage &local_storage : storage)
801 res += local_storage.val;
802 }
803
804 VectorNd local_displacement(size_);
805 Eigen::MatrixXd forces;
806
807 ElementAssemblyValues vals;
808 // Neumann
809 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
810 Eigen::VectorXd weights;
811 Eigen::VectorXi global_primitive_ids;
812 for (const auto &lb : local_neumann_boundary)
813 {
814 const int e = lb.element_id();
815 const basis::ElementBases &gbs = gbases_[e];
816 const basis::ElementBases &bs = bases_[e];
817
818 for (int i = 0; i < lb.size(); ++i)
819 {
820 const int primitive_global_id = lb.global_primitive_id(i);
821
822 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
823 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
824
825 if (!has_samples)
826 continue;
827
828 const basis::ElementBases &gbs = gbases_[e];
829 const basis::ElementBases &bs = bases_[e];
830
831 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
832
833 for (int n = 0; n < vals.jac_it.size(); ++n)
834 {
835 trafo = vals.jac_it[n].inverse();
836
837 if (displacement_prev.size() > 0)
838 {
839 assert(size_ == 2 || size_ == 3);
840 deform_mat.resize(size_, size_);
841 deform_mat.setZero();
842 for (const auto &b : vals.basis_values)
843 {
844 for (const auto &g : b.global)
845 {
846 for (int d = 0; d < size_; ++d)
847 {
848 deform_mat.row(d) += displacement_prev(g.index * size_ + d) * b.grad.row(n);
849 }
850 }
851 }
852
853 trafo += deform_mat;
854 }
855
856 normals.row(n) = normals.row(n) * trafo.inverse();
857 normals.row(n).normalize();
858 }
859 problem_.neumann_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, forces);
860
861 // UIState::ui_state().debug_data().add_points(vals.val, Eigen::RowVector3d(1,0,0));
862
863 for (long p = 0; p < weights.size(); ++p)
864 {
865 local_displacement.setZero();
866
867 for (size_t i = 0; i < vals.basis_values.size(); ++i)
868 {
869 const auto &vv = vals.basis_values[i];
870 assert(vv.val.size() == weights.size());
871 const double b_val = vv.val(p);
872
873 for (int d = 0; d < size_; ++d)
874 {
875 for (std::size_t ii = 0; ii < vv.global.size(); ++ii)
876 {
877 local_displacement(d) += (vv.global[ii].val * b_val) * displacement(vv.global[ii].index * size_ + d);
878 }
879 }
880 }
881
882 for (int d = 0; d < size_; ++d)
883 res -= forces(p, d) * local_displacement(d) * weights(p);
884 }
885 }
886 }
887
888 return res;
889 }
890
891 void RhsAssembler::compute_energy_hess(
892 const std::vector<int> &bounday_nodes,
893 const int resolution,
894 const std::vector<mesh::LocalBoundary> &local_neumann_boundary,
895 const Eigen::MatrixXd &displacement,
896 const double t,
897 const bool project_to_psd,
898 StiffnessMatrix &hess) const
899 {
900 hess.resize(n_basis_ * size_, n_basis_ * size_);
901 if (displacement.size() == 0)
902 return;
903
904 std::vector<Eigen::Triplet<double>> entries, entries_t;
905
906 ElementAssemblyValues vals;
907 Eigen::MatrixXd uv, samples, gtmp, rhs_fun, deform_mat, jac_mat, trafo;
908 Eigen::VectorXi global_primitive_ids;
909 Eigen::MatrixXd points, normals;
910 Eigen::VectorXd weights;
911 Eigen::MatrixXd local_hessian;
912
913 for (const auto &lb : local_neumann_boundary)
914 {
915 const int e = lb.element_id();
916 const basis::ElementBases &gbs = gbases_[e];
917 const basis::ElementBases &bs = bases_[e];
918
919 for (int i = 0; i < lb.size(); ++i)
920 {
921 const int primitive_global_id = lb.global_primitive_id(i);
922
923 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
924 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
925
926 if (!has_samples)
927 continue;
928
929 Eigen::MatrixXd reference_normals = normals;
930
931 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
932
933 std::vector<std::vector<Eigen::MatrixXd>> grad_normal;
934 for (int n = 0; n < vals.jac_it.size(); ++n)
935 {
936 trafo = vals.jac_it[n].inverse();
937
938 assert(size_ == 2 || size_ == 3);
939 deform_mat.resize(size_, size_);
940 deform_mat.setZero();
941 jac_mat.resize(size_, vals.basis_values.size());
942 int b_idx = 0;
943 for (const auto &b : vals.basis_values)
944 {
945 jac_mat.col(b_idx++) = b.grad.row(n);
946
947 for (const auto &g : b.global)
948 for (int d = 0; d < size_; ++d)
949 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
950 }
951
952 trafo += deform_mat;
953 trafo = trafo.inverse();
954
955 Eigen::VectorXd displaced_normal = normals.row(n) * trafo;
956 normals.row(n) = displaced_normal / displaced_normal.norm();
957
958 std::vector<Eigen::MatrixXd> grad;
959 {
960 Eigen::MatrixXd vec = -(jac_mat.transpose() * trafo * reference_normals.row(n).transpose());
961 // Gradient of the displaced normal computation
962 for (int k = 0; k < size_; ++k)
963 {
964 Eigen::MatrixXd grad_i(jac_mat.rows(), jac_mat.cols());
965 grad_i.setZero();
966 for (int m = 0; m < jac_mat.rows(); ++m)
967 for (int l = 0; l < jac_mat.cols(); ++l)
968 grad_i(m, l) = -(reference_normals.row(n) * trafo)(m) * (jac_mat.transpose() * trafo)(l, k);
969 grad.push_back(grad_i);
970 }
971 }
972
973 {
974 Eigen::MatrixXd normalization_chain_rule = (normals.row(n).transpose() * normals.row(n));
975 normalization_chain_rule = Eigen::MatrixXd::Identity(size_, size_) - normalization_chain_rule;
976 normalization_chain_rule /= displaced_normal.norm();
977
978 Eigen::VectorXd vec(size_);
979 b_idx = 0;
980 for (const auto &b : vals.basis_values)
981 {
982 for (int d = 0; d < size_; ++d)
983 {
984 for (int k = 0; k < size_; ++k)
985 vec(k) = grad[k](d, b_idx);
986 vec = normalization_chain_rule * vec;
987 for (int k = 0; k < size_; ++k)
988 grad[k](d, b_idx) = vec(k);
989 }
990 b_idx++;
991 }
992 }
993
994 grad_normal.push_back(grad);
995 }
996 Eigen::MatrixXd rhs_fun;
997 problem_.neumann_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, rhs_fun);
998
999 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
1000
1001 const bool is_pressure = problem_.is_boundary_pressure(mesh_.get_boundary_id(primitive_global_id));
1002 if (!is_pressure)
1003 continue;
1004
1005 local_hessian.setZero(vals.basis_values.size() * size_, vals.basis_values.size() * size_);
1006
1007 for (long n = 0; n < nodes.size(); ++n)
1008 {
1009 // const auto &b = bs.bases[nodes(n)];
1010 const AssemblyValues &v = vals.basis_values[nodes(n)];
1011 for (int d = 0; d < size_; ++d)
1012 {
1013 for (size_t g = 0; g < v.global.size(); ++g)
1014 {
1015 const int g_index = v.global[g].index * size_ + d;
1016 const bool is_neumann = std::find(bounday_nodes.begin(), bounday_nodes.end(), g_index) == bounday_nodes.end();
1017
1018 if (is_neumann)
1019 {
1020 for (long ni = 0; ni < nodes.size(); ++ni)
1021 {
1022 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
1023 for (int di = 0; di < size_; ++di)
1024 {
1025 for (size_t gi = 0; gi < vi.global.size(); ++gi)
1026 {
1027 const int gi_index = vi.global[gi].index * size_ + di;
1028 double value = 0;
1029
1030 for (int q = 0; q < vals.jac_it.size(); ++q)
1031 {
1032 double pressure_val = rhs_fun.row(q).dot(normals.row(q));
1033
1034 // value += grad_normal[ni](d, nodes(ni) * size_ + di) * pressure_val * weights(q) * vi.val(q);
1035 value += grad_normal[q][d](di, nodes(ni)) * pressure_val * weights(q) * vi.val(q);
1036 }
1037
1038 value *= v.global[g].val;
1039
1040 const bool is_neumann_i = std::find(bounday_nodes.begin(), bounday_nodes.end(), gi_index) == bounday_nodes.end();
1041
1042 if (is_neumann_i)
1043 {
1044 local_hessian(nodes(n) * size_ + d, nodes(ni) * size_ + di) = value;
1045 }
1046 }
1047 }
1048 }
1049 }
1050 }
1051 }
1052 }
1053
1054 if (project_to_psd)
1055 local_hessian = ipc::project_to_psd(local_hessian);
1056
1057 for (long n = 0; n < nodes.size(); ++n)
1058 {
1059 const AssemblyValues &v = vals.basis_values[nodes(n)];
1060 for (int d = 0; d < size_; ++d)
1061 {
1062 for (size_t g = 0; g < v.global.size(); ++g)
1063 {
1064 const int g_index = v.global[g].index * size_ + d;
1065
1066 for (long ni = 0; ni < nodes.size(); ++ni)
1067 {
1068 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
1069 for (int di = 0; di < size_; ++di)
1070 {
1071 for (size_t gi = 0; gi < vi.global.size(); ++gi)
1072 {
1073 const int gi_index = vi.global[gi].index * size_ + di;
1074 entries.push_back(Eigen::Triplet<double>(g_index, gi_index, local_hessian(nodes(n) * size_ + d, nodes(ni) * size_ + di)));
1075 }
1076 }
1077 }
1078 }
1079 }
1080 }
1081 }
1082 }
1083
1084 hess.setFromTriplets(entries.begin(), entries.end());
1085 }
1086 } // namespace assembler
1087} // namespace polyfem
vec
Eigen::MatrixXd vec
Definition Assembler.cpp:72
val
double val
Definition Assembler.cpp:86
da
QuadratureVector da
Definition Assembler.cpp:23
vals
ElementAssemblyValues vals
Definition Assembler.cpp:22
quadrature
Quadrature quadrature
Definition MassMatrixAssembler.cpp:30
entries
std::vector< Eigen::Triplet< double > > entries
Definition MassMatrixAssembler.cpp:25
mat
Eigen::MatrixXd mat
Definition TransientForm.cpp:24
polyfem::assembler::Assembler::set_size
virtual void set_size(const int size)
Definition Assembler.hpp:64
polyfem::assembler::AssemblyValsCache
Caches basis evaluation and geometric mapping at every element.
Definition AssemblyValsCache.hpp:11
polyfem::assembler::AssemblyValsCache::compute
void compute(const int el_index, const bool is_volume, const basis::ElementBases &basis, const basis::ElementBases &gbasis, ElementAssemblyValues &vals) const
retrieves cached basis evaluation and geometric for the given element if it doesn't exist,...
Definition AssemblyValsCache.cpp:33
polyfem::assembler::AssemblyValues
stores per local bases evaluations
Definition AssemblyValues.hpp:14
polyfem::assembler::AssemblyValues::global
std::vector< basis::Local2Global > global
Definition AssemblyValues.hpp:20
polyfem::assembler::ElementAssemblyValues
stores per element basis values at given quadrature points and geometric mapping
Definition ElementAssemblyValues.hpp:14
polyfem::assembler::Mass::add_multimaterial
void add_multimaterial(const int index, const json &params, const Units &units) override
inialize material parameter
Definition Mass.cpp:33
polyfem::assembler::Mass::assemble
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 9, 1 > assemble(const LinearAssemblerData &data) const override
computes and returns local stiffness matrix (1x1) for bases i,j (where i,j is passed in through data)...
Definition Mass.cpp:5
polyfem::assembler::Problem::neumann_bc
virtual void neumann_bc(const mesh::Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &normals, const double t, Eigen::MatrixXd &val) const
Definition Problem.hpp:31
polyfem::assembler::Problem::is_rhs_zero
virtual bool is_rhs_zero() const =0
polyfem::assembler::Problem::is_dimension_dirichet
virtual bool is_dimension_dirichet(const int tag, const int dim) const
Definition Problem.hpp:61
polyfem::assembler::Problem::all_dimensions_dirichlet
virtual bool all_dimensions_dirichlet() const
Definition Problem.hpp:64
polyfem::assembler::Problem::initial_acceleration
virtual void initial_acceleration(const mesh::Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) const
Definition Problem.hpp:55
polyfem::assembler::Problem::is_scalar
virtual bool is_scalar() const =0
polyfem::assembler::Problem::rhs
virtual void rhs(const assembler::Assembler &assembler, const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const =0
polyfem::assembler::Problem::dirichlet_nodal_value
virtual void dirichlet_nodal_value(const mesh::Mesh &mesh, const int node_id, const RowVectorNd &pt, const double t, Eigen::MatrixXd &val) const
Definition Problem.hpp:37
polyfem::assembler::Problem::initial_velocity
virtual void initial_velocity(const mesh::Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) const
Definition Problem.hpp:54
polyfem::assembler::Problem::is_constant_in_time
virtual bool is_constant_in_time() const
Definition Problem.hpp:51
polyfem::assembler::Problem::is_nodal_dimension_dirichlet
virtual bool is_nodal_dimension_dirichlet(const int n_id, const int tag, const int dim) const
Definition Problem.hpp:62
polyfem::assembler::Problem::is_boundary_pressure
virtual bool is_boundary_pressure(const int boundary_id) const
Definition Problem.hpp:35
polyfem::assembler::Problem::initial_solution
virtual void initial_solution(const mesh::Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) const
Definition Problem.hpp:53
polyfem::assembler::Problem::dirichlet_bc
virtual void dirichlet_bc(const mesh::Mesh &mesh, const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const double t, Eigen::MatrixXd &val) const =0
polyfem::assembler::RhsAssembler::set_bc
void set_bc(const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &bounday_nodes, const int resolution, const std::vector< mesh::LocalBoundary > &local_neumann_boundary, Eigen::MatrixXd &rhs, const Eigen::MatrixXd &displacement=Eigen::MatrixXd(), const double t=1) const
polyfem::assembler::RhsAssembler::time_bc
void time_bc(const std::function< void(const mesh::Mesh &, const Eigen::MatrixXi &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &fun, Eigen::MatrixXd &sol) const
Definition RhsAssembler.cpp:138
polyfem::assembler::RhsAssembler::mesh
const mesh::Mesh & mesh() const
Definition RhsAssembler.hpp:84
polyfem::assembler::RhsAssembler::gbases_
const std::vector< basis::ElementBases > & gbases_
basis functions associated with geometric mapping
Definition RhsAssembler.hpp:123
polyfem::assembler::RhsAssembler::size_
const int size_
dimension of problem
Definition RhsAssembler.hpp:121
polyfem::assembler::RhsAssembler::compute_energy_grad
void compute_energy_grad(const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &bounday_nodes, const Density &density, const int resolution, const std::vector< mesh::LocalBoundary > &local_neumann_boundary, const Eigen::MatrixXd &final_rhs, const double t, Eigen::MatrixXd &rhs) const
Definition RhsAssembler.cpp:712
polyfem::assembler::RhsAssembler::dirichlet_nodes_position_
const std::vector< RowVectorNd > & dirichlet_nodes_position_
Definition RhsAssembler.hpp:129
polyfem::assembler::RhsAssembler::compute_energy
double compute_energy(const Eigen::MatrixXd &displacement, const Eigen::MatrixXd &displacement_prev, const std::vector< mesh::LocalBoundary > &local_neumann_boundary, const Density &density, const int resolution, const double t) const
Definition RhsAssembler.cpp:736
polyfem::assembler::RhsAssembler::compute_energy_hess
void compute_energy_hess(const std::vector< int > &bounday_nodes, const int resolution, const std::vector< mesh::LocalBoundary > &local_neumann_boundary, const Eigen::MatrixXd &displacement, const double t, const bool project_to_psd, StiffnessMatrix &hess) const
Definition RhsAssembler.cpp:891
polyfem::assembler::RhsAssembler::ass_vals_cache_
const AssemblyValsCache & ass_vals_cache_
Definition RhsAssembler.hpp:124
polyfem::assembler::RhsAssembler::initial_solution
void initial_solution(Eigen::MatrixXd &sol) const
Definition RhsAssembler.cpp:114
polyfem::assembler::RhsAssembler::obstacle_
const mesh::Obstacle & obstacle_
Definition RhsAssembler.hpp:119
polyfem::assembler::RhsAssembler::lsq_bc
void lsq_bc(const std::function< void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &df, const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &bounday_nodes, const int resolution, Eigen::MatrixXd &rhs) const
Definition RhsAssembler.cpp:235
polyfem::assembler::RhsAssembler::bases_
const std::vector< basis::ElementBases > & bases_
basis functions associated with solution
Definition RhsAssembler.hpp:122
polyfem::assembler::RhsAssembler::initial_acceleration
void initial_acceleration(Eigen::MatrixXd &sol) const
Definition RhsAssembler.cpp:130
polyfem::assembler::RhsAssembler::initial_velocity
void initial_velocity(Eigen::MatrixXd &sol) const
Definition RhsAssembler.cpp:122
polyfem::assembler::RhsAssembler::assemble
void assemble(const Density &density, Eigen::MatrixXd &rhs, const double t=1) const
Definition RhsAssembler.cpp:61
polyfem::assembler::RhsAssembler::dirichlet_nodes_
const std::vector< int > & dirichlet_nodes_
Definition RhsAssembler.hpp:128
polyfem::assembler::RhsAssembler::integrate_bc
void integrate_bc(const std::function< void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &df, const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &bounday_nodes, const int resolution, Eigen::MatrixXd &rhs) const
Definition RhsAssembler.cpp:490
polyfem::assembler::RhsAssembler::RhsAssembler
RhsAssembler(const Assembler &assembler, const mesh::Mesh &mesh, const mesh::Obstacle &obstacle, const std::vector< int > &dirichlet_nodes, const std::vector< int > &neumann_nodes, const std::vector< RowVectorNd > &dirichlet_nodes_position, const std::vector< RowVectorNd > &neumann_nodes_position, const int n_basis, const int size, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, const AssemblyValsCache &ass_vals_cache, const Problem &problem, const std::string bc_method, const json &solver_params)
Definition RhsAssembler.cpp:34
polyfem::assembler::RhsAssembler::sample_bc
void sample_bc(const std::function< void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &df, const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &bounday_nodes, Eigen::MatrixXd &rhs) const
Definition RhsAssembler.cpp:423
polyfem::basis::Basis
Represents one basis function and its gradient.
Definition Basis.hpp:44
polyfem::basis::Basis::global
const std::vector< Local2Global > & global() const
Definition Basis.hpp:104
polyfem::basis::ElementBases
Stores the basis functions for a given element in a mesh (facet in 2d, cell in 3d).
Definition ElementBases.hpp:17
polyfem::basis::ElementBases::eval_geom_mapping
void eval_geom_mapping(const Eigen::MatrixXd &samples, Eigen::MatrixXd &mapped) const
Map the sample positions in the parametric domain to the object domain (if the element has no paramet...
Definition ElementBases.cpp:19
polyfem::basis::ElementBases::evaluate_bases
void evaluate_bases(const Eigen::MatrixXd &uv, std::vector< assembler::AssemblyValues > &basis_values) const
evaluate stored bases at given points on the reference element saves results to basis_values
Definition ElementBases.hpp:71
polyfem::basis::ElementBases::local_nodes_for_primitive
Eigen::VectorXi local_nodes_for_primitive(const int local_index, const mesh::Mesh &mesh) const
Definition ElementBases.hpp:35
polyfem::basis::ElementBases::bases
std::vector< Basis > bases
one basis function per node in the element
Definition ElementBases.hpp:27
polyfem::mesh::Mesh
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:39
polyfem::mesh::Mesh::get_boundary_id
virtual int get_boundary_id(const int primitive) const
Get the boundary selection of an element (face in 3d, edge in 2d)
Definition Mesh.hpp:478
polyfem::mesh::Mesh::is_volume
virtual bool is_volume() const =0
checks if mesh is volume
polyfem::mesh::Mesh::dimension
int dimension() const
utily for dimension
Definition Mesh.hpp:148
polyfem::mesh::Mesh::get_node_id
virtual int get_node_id(const int node_id) const
Get the boundary selection of a node.
Definition Mesh.hpp:487
polyfem::mesh::Obstacle::update_displacement
void update_displacement(const double t, Eigen::MatrixXd &sol) const
Definition Obstacle.cpp:212
polyfem::utils::BoundarySampler::sample_boundary
static bool sample_boundary(const mesh::LocalBoundary &local_boundary, const int n_samples, const mesh::Mesh &mesh, const bool skip_computation, Eigen::MatrixXd &uv, Eigen::MatrixXd &samples, Eigen::VectorXi &global_primitive_ids)
Definition BoundarySampler.cpp:530
polyfem::utils::BoundarySampler::boundary_quadrature
static bool boundary_quadrature(const mesh::LocalBoundary &local_boundary, const int order, const mesh::Mesh &mesh, const bool skip_computation, Eigen::MatrixXd &uv, Eigen::MatrixXd &points, Eigen::MatrixXd &normals, Eigen::VectorXd &weights, Eigen::VectorXi &global_primitive_ids)
Definition BoundarySampler.cpp:463
p_bases.points
points
Definition p_bases.py:225
p_bases.nodes
str nodes
Definition p_bases.py:372
polyfem::utils::get_local_thread_storage
auto & get_local_thread_storage(Storages &storage, int thread_id)
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::VectorNd
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 3, 1 > VectorNd
Definition Types.hpp:9
polyfem::json
nlohmann::json json
Definition Common.hpp:9
polyfem::StiffnessMatrix
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:20