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_initialized() || 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,
237 const std::vector<int> &bounday_nodes,
238 const int resolution,
239 Eigen::MatrixXd &rhs) const
240 {
241 const int n_el = int(bases_.size());
242
243 Eigen::MatrixXd uv, samples, gtmp, rhs_fun;
244 Eigen::VectorXi global_primitive_ids;
245
246 const int actual_dim = problem_.is_scalar() ? 1 : mesh_.dimension();
247
248 Eigen::Matrix<bool, Eigen::Dynamic, 1> is_boundary(n_basis_);
249 is_boundary.setConstant(false);
250 int skipped_count = 0;
251 for (int b : bounday_nodes)
252 {
253 int bindex = b / actual_dim;
254
255 if (bindex < is_boundary.size())
256 is_boundary[bindex] = true;
257 else
258 skipped_count++;
259 }
260 assert(skipped_count <= 1);
261
262 for (int d = 0; d < size_; ++d)
263 {
264 int index = 0;
265 std::vector<int> indices;
266 indices.reserve(n_el * 10);
267 std::vector<int> tags;
268 tags.reserve(n_el * 10);
269
270 long total_size = 0;
271
272 Eigen::VectorXi global_index_to_col(n_basis_);
273 global_index_to_col.setConstant(-1);
274
275 std::vector<AssemblyValues> tmp_val;
276
277 for (const auto &lb : local_boundary)
278 {
279 const int e = lb.element_id();
280 bool has_samples = utils::BoundarySampler::sample_boundary(lb, resolution, mesh_, true, uv, samples, global_primitive_ids);
281
282 if (!has_samples)
283 continue;
284
285 const basis::ElementBases &bs = bases_[e];
286 bs.evaluate_bases(samples, tmp_val);
287 const int n_local_bases = int(bs.bases.size());
288 assert(global_primitive_ids.size() == samples.rows());
289
290 for (int s = 0; s < samples.rows(); ++s)
291 {
292 const int tag = mesh_.get_boundary_id(global_primitive_ids(s));
293 if (!problem_.all_dimensions_dirichlet() && !problem_.is_dimension_dirichet(tag, d))
294 continue;
295
296 total_size++;
297
298 for (int j = 0; j < n_local_bases; ++j)
299 {
300 const basis::Basis &b = bs.bases[j];
301 const double tmp = tmp_val[j].val(s);
302
303 if (fabs(tmp) < 1e-10)
304 continue;
305
306 for (std::size_t ii = 0; ii < b.global().size(); ++ii)
307 {
308 // pt found
309 if (is_boundary[b.global()[ii].index])
310 {
311 if (global_index_to_col(b.global()[ii].index) == -1)
312 {
313 global_index_to_col(b.global()[ii].index) = index++;
314 indices.push_back(b.global()[ii].index);
315 tags.push_back(tag);
316 assert(indices.size() == size_t(index));
317 }
318 }
319 }
320 }
321 }
322 }
323
324 Eigen::MatrixXd global_rhs = Eigen::MatrixXd::Zero(total_size, 1);
325
326 const long buffer_size = total_size * long(indices.size());
327 std::vector<Eigen::Triplet<double>> entries, entries_t;
328
329 index = 0;
330
331 int global_counter = 0;
332 Eigen::MatrixXd mapped;
333
334 for (const auto &lb : local_boundary)
335 {
336 const int e = lb.element_id();
337 bool has_samples = utils::BoundarySampler::sample_boundary(lb, resolution, mesh_, false, uv, samples, global_primitive_ids);
338
339 if (!has_samples)
340 continue;
341
342 const basis::ElementBases &bs = bases_[e];
343 const basis::ElementBases &gbs = gbases_[e];
344 const int n_local_bases = int(bs.bases.size());
345
346 gbs.eval_geom_mapping(samples, mapped);
347
348 bs.evaluate_bases(samples, tmp_val);
349 df(global_primitive_ids, uv, mapped, rhs_fun);
350
351 for (int s = 0; s < samples.rows(); ++s)
352 {
353 const int tag = mesh_.get_boundary_id(global_primitive_ids(s));
354 if (!problem_.all_dimensions_dirichlet() && !problem_.is_dimension_dirichet(tag, d))
355 continue;
356
357 for (int j = 0; j < n_local_bases; ++j)
358 {
359 const basis::Basis &b = bs.bases[j];
360 const double tmp = tmp_val[j].val(s);
361
362 for (std::size_t ii = 0; ii < b.global().size(); ++ii)
363 {
364 auto item = global_index_to_col(b.global()[ii].index);
365 if (item != -1)
366 {
367 entries.push_back(Eigen::Triplet<double>(global_counter, item, tmp * b.global()[ii].val));
368 entries_t.push_back(Eigen::Triplet<double>(item, global_counter, tmp * b.global()[ii].val));
369 }
370 }
371 }
372
373 global_rhs(global_counter) = rhs_fun(s, d);
374 global_counter++;
375 }
376 }
377
378 assert(global_counter == total_size);
379
380 if (total_size > 0)
381 {
382 const double mmin = global_rhs.minCoeff();
383 const double mmax = global_rhs.maxCoeff();
384
385 if (fabs(mmin) < 1e-8 && fabs(mmax) < 1e-8)
386 {
387 for (size_t i = 0; i < indices.size(); ++i)
388 {
389 const int tag = tags[i];
390 if (problem_.all_dimensions_dirichlet() || problem_.is_dimension_dirichet(tag, d))
391 rhs(indices[i] * size_ + d) = 0;
392 }
393 }
394 else
395 {
396 StiffnessMatrix mat(int(total_size), int(indices.size()));
397 mat.setFromTriplets(entries.begin(), entries.end());
398
399 StiffnessMatrix mat_t(int(indices.size()), int(total_size));
400 mat_t.setFromTriplets(entries_t.begin(), entries_t.end());
401
402 StiffnessMatrix A = mat_t * mat;
403 Eigen::VectorXd b = mat_t * global_rhs;
404
405 Eigen::VectorXd coeffs(b.rows(), 1);
406 auto solver = linear::Solver::create(solver_params_, logger());
407 logger().info("Solve RHS using {} linear solver", solver->name());
408 solver->analyze_pattern(A, A.rows());
409 solver->factorize(A);
410 coeffs.setZero();
411 solver->solve(b, coeffs);
412
413 logger().trace("RHS solve error {}", (A * coeffs - b).norm());
414
415 for (long i = 0; i < coeffs.rows(); ++i)
416 {
417 const int tag = tags[i];
418 if (problem_.all_dimensions_dirichlet() || problem_.is_dimension_dirichet(tag, d))
419 rhs(indices[i] * size_ + d) = coeffs(i);
420 }
421 }
422 }
423 }
424 }
425
426 void RhsAssembler::sample_bc(const std::function<void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &df,
427 const std::vector<LocalBoundary> &local_boundary, const std::vector<int> &bounday_nodes, Eigen::MatrixXd &rhs) const
428 {
429 const int n_el = int(bases_.size());
430
431 Eigen::MatrixXd rhs_fun;
432 Eigen::VectorXi global_primitive_ids(1);
433 Eigen::MatrixXd nans(1, 1);
434 nans(0) = std::nan("");
435
436#ifndef NDEBUG
437 Eigen::Matrix<bool, Eigen::Dynamic, 1> is_boundary(n_basis_);
438 is_boundary.setConstant(false);
439
440 const int actual_dim = problem_.is_scalar() ? 1 : mesh_.dimension();
441
442 int skipped_count = 0;
443 for (int b : bounday_nodes)
444 {
445 int bindex = b / actual_dim;
446
447 if (bindex < is_boundary.size())
448 is_boundary[bindex] = true;
449 else
450 skipped_count++;
451 }
452 assert(skipped_count <= 1);
453#endif
454
455 for (const auto &lb : local_boundary)
456 {
457 const int e = lb.element_id();
458 const basis::ElementBases &bs = bases_[e];
459
460 for (int i = 0; i < lb.size(); ++i)
461 {
462 global_primitive_ids(0) = lb.global_primitive_id(i);
463 const auto nodes = bs.local_nodes_for_primitive(global_primitive_ids(0), mesh_);
464 assert(global_primitive_ids.size() == 1);
465 const int tag = mesh_.get_boundary_id(global_primitive_ids(0));
466
467 for (long n = 0; n < nodes.size(); ++n)
468 {
469 const auto &b = bs.bases[nodes(n)];
470 const auto &glob = b.global();
471
472 for (size_t ii = 0; ii < glob.size(); ++ii)
473 {
474 assert(is_boundary[glob[ii].index]);
475
476 // TODO, missing UV!!!!
477 df(global_primitive_ids, nans, glob[ii].node, rhs_fun);
478
479 for (int d = 0; d < size_; ++d)
480 {
481 if (problem_.all_dimensions_dirichlet() || problem_.is_dimension_dirichet(tag, d))
482 {
483 assert(problem_.all_dimensions_dirichlet() || std::find(bounday_nodes.begin(), bounday_nodes.end(), glob[ii].index * size_ + d) != bounday_nodes.end());
484 rhs(glob[ii].index * size_ + d) = rhs_fun(0, d);
485 }
486 }
487 }
488 }
489 }
490 }
491 }
492
493 void RhsAssembler::set_bc(
494 const std::function<void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &df,
495 const std::function<void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &nf,
496 const std::vector<LocalBoundary> &local_boundary,
497 const std::vector<int> &bounday_nodes,
498 const QuadratureOrders &resolution,
499 const std::vector<LocalBoundary> &local_neumann_boundary,
500 const Eigen::MatrixXd &displacement,
501 const double t,
502 Eigen::MatrixXd &rhs) const
503 {
504 if (bc_method_ == "sample")
505 sample_bc(df, local_boundary, bounday_nodes, rhs);
506 else
507 lsq_bc(df, local_boundary, bounday_nodes, resolution[0], rhs);
508
509 if (bounday_nodes.size() > 0)
510 {
511 Eigen::MatrixXd tmp_val;
512 for (int n = 0; n < dirichlet_nodes_.size(); ++n)
513 {
514 const auto &n_id = dirichlet_nodes_[n];
515 const auto &pt = dirichlet_nodes_position_[n];
516
517 const int tag = mesh_.get_node_id(n_id);
518 problem_.dirichlet_nodal_value(mesh_, n_id, pt, t, tmp_val);
519 assert(tmp_val.size() == size_);
520
521 for (int d = 0; d < size_; ++d)
522 {
523 if (!problem_.is_nodal_dimension_dirichlet(n_id, tag, d))
524 continue;
525 const int g_index = n_id * size_ + d;
526 rhs(g_index) = tmp_val(d);
527 }
528 }
529 }
530
531 // Neumann
532 Eigen::MatrixXd uv, samples, gtmp, rhs_fun, deform_mat, trafo;
533 Eigen::VectorXi global_primitive_ids;
534 Eigen::MatrixXd points, normals;
535 Eigen::VectorXd weights;
536
537 ElementAssemblyValues vals;
538
539 for (const auto &lb : local_neumann_boundary)
540 {
541 const int e = lb.element_id();
542 const basis::ElementBases &gbs = gbases_[e];
543 const basis::ElementBases &bs = bases_[e];
544
545 for (int i = 0; i < lb.size(); ++i)
546 {
547 const int primitive_global_id = lb.global_primitive_id(i);
548 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
549 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
550
551 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
552
553 for (int n = 0; n < vals.jac_it.size(); ++n)
554 {
555 trafo = vals.jac_it[n].inverse();
556
557 if (displacement.size() > 0)
558 {
559 assert(size_ == 2 || size_ == 3);
560 deform_mat.resize(size_, size_);
561 deform_mat.setZero();
562 for (const auto &b : vals.basis_values)
563 {
564 for (const auto &g : b.global)
565 {
566 for (int d = 0; d < size_; ++d)
567 {
568 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
569 }
570 }
571 }
572
573 trafo += deform_mat;
574 }
575
576 normals.row(n) = normals.row(n) * trafo.inverse();
577 normals.row(n).normalize();
578 }
579
580 // problem_.neumann_bc(mesh_, global_primitive_ids, vals.val, t, rhs_fun);
581 nf(global_primitive_ids, uv, vals.val, normals, rhs_fun);
582
583 // UIState::ui_state().debug_data().add_points(vals.val, Eigen::RowVector3d(0,1,0));
584
585 for (int d = 0; d < size_; ++d)
586 rhs_fun.col(d) = rhs_fun.col(d).array() * weights.array();
587
588 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
589
590 for (long n = 0; n < nodes.size(); ++n)
591 {
592 // const auto &b = bs.bases[nodes(n)];
593 const AssemblyValues &v = vals.basis_values[nodes(n)];
594 for (int d = 0; d < size_; ++d)
595 {
596 const double rhs_value = (rhs_fun.col(d).array() * v.val.array()).sum();
597
598 for (size_t g = 0; g < v.global.size(); ++g)
599 {
600 const int g_index = v.global[g].index * size_ + d;
601 const bool is_neumann = std::find(bounday_nodes.begin(), bounday_nodes.end(), g_index) == bounday_nodes.end();
602
603 if (is_neumann)
604 {
605 rhs(g_index) += rhs_value * v.global[g].val;
606 }
607 }
608 }
609 }
610 }
611 }
612
613 // TODO add nodal neumann
614 }
615
616 void RhsAssembler::set_bc(const std::vector<LocalBoundary> &local_boundary,
617 const std::vector<int> &bounday_nodes,
618 const QuadratureOrders &resolution,
619 const std::vector<LocalBoundary> &local_neumann_boundary,
620 Eigen::MatrixXd &rhs,
621 const Eigen::MatrixXd &displacement,
622 const double t) const
623 {
624 set_bc(
625 [&](const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) {
626 problem_.dirichlet_bc(mesh_, global_ids, uv, pts, t, val);
627 },
628 [&](const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &normals, Eigen::MatrixXd &val) {
629 problem_.neumann_bc(mesh_, global_ids, uv, pts, normals, t, val);
630 },
631 local_boundary, bounday_nodes, resolution, local_neumann_boundary, displacement, t, rhs);
632
633 obstacle_.update_displacement(t, rhs);
634 }
635
636 void RhsAssembler::compute_energy_grad(const std::vector<LocalBoundary> &local_boundary,
637 const std::vector<int> &bounday_nodes,
638 const Density &density,
639 const QuadratureOrders &resolution,
640 const std::vector<LocalBoundary> &local_neumann_boundary,
641 const Eigen::MatrixXd &final_rhs,
642 const double t,
643 Eigen::MatrixXd &rhs) const
644 {
645 if (problem_.is_constant_in_time())
646 {
647 rhs = final_rhs;
648 }
649 else
650 {
651 assemble(density, rhs, t);
652 rhs *= -1;
653
654 if (rhs.size() != final_rhs.size())
655 {
656 const int prev_size = rhs.size();
657 rhs.conservativeResize(final_rhs.size(), rhs.cols());
658 // Zero initial pressure
659 rhs.block(prev_size, 0, final_rhs.size() - prev_size, rhs.cols()).setZero();
660 rhs(rhs.size() - 1) = 0;
661 }
662
663 assert(rhs.size() == final_rhs.size());
664 }
665 }
666
667 double RhsAssembler::compute_energy(const Eigen::MatrixXd &displacement,
668 const Eigen::MatrixXd &displacement_prev,
669 const std::vector<LocalBoundary> &local_neumann_boundary,
670 const Density &density,
671 const QuadratureOrders &resolution,
672 const double t) const
673 {
674
675 double res = 0;
676
677 if (!problem_.is_rhs_zero())
678 {
679 auto storage = create_thread_storage(LocalThreadScalarStorage());
680 const int n_bases = int(bases_.size());
681
682 maybe_parallel_for(n_bases, [&](int start, int end, int thread_id) {
683 LocalThreadScalarStorage &local_storage = get_local_thread_storage(storage, thread_id);
684 VectorNd local_displacement(size_);
685 Eigen::MatrixXd forces;
686
687 for (int e = start; e < end; ++e)
688 {
689 ElementAssemblyValues &vals = local_storage.vals;
690 // vals.compute(e, mesh_.is_volume(), bases_[e], gbases_[e]);
691 ass_vals_cache_.compute(e, mesh_.is_volume(), bases_[e], gbases_[e], vals);
692
693 const Quadrature &quadrature = vals.quadrature;
694 const Eigen::VectorXd da = vals.det.array() * quadrature.weights.array();
695
696 problem_.rhs(assembler_, vals.val, t, forces);
697 assert(forces.rows() == da.size());
698 assert(forces.cols() == size_);
699
700 for (long p = 0; p < da.size(); ++p)
701 {
702 local_displacement.setZero();
703
704 for (size_t i = 0; i < vals.basis_values.size(); ++i)
705 {
706 const auto &bs = vals.basis_values[i];
707 assert(bs.val.size() == da.size());
708 const double b_val = bs.val(p);
709
710 for (int d = 0; d < size_; ++d)
711 {
712 for (std::size_t ii = 0; ii < bs.global.size(); ++ii)
713 {
714 local_displacement(d) += (bs.global[ii].val * b_val) * displacement(bs.global[ii].index * size_ + d);
715 }
716 }
717 }
718 // const double rho = problem_.is_time_dependent() ? density(vals.quadrature.points.row(p), vals.val.row(p), vals.element_id) : 1;
719 const double rho = density(vals.quadrature.points.row(p), vals.val.row(p), t, vals.element_id);
720
721 for (int d = 0; d < size_; ++d)
722 {
723 local_storage.val += forces(p, d) * local_displacement(d) * da(p) * rho;
724 // res += forces(p, d) * local_displacement(d) * da(p);
725 }
726 }
727 }
728 });
729
730 // Serially merge local storages
731 for (const LocalThreadScalarStorage &local_storage : storage)
732 res += local_storage.val;
733 }
734
735 VectorNd local_displacement(size_);
736 Eigen::MatrixXd forces;
737
738 ElementAssemblyValues vals;
739 // Neumann
740 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
741 Eigen::VectorXd weights;
742 Eigen::VectorXi global_primitive_ids;
743 for (const auto &lb : local_neumann_boundary)
744 {
745 const int e = lb.element_id();
746 const basis::ElementBases &gbs = gbases_[e];
747 const basis::ElementBases &bs = bases_[e];
748
749 for (int i = 0; i < lb.size(); ++i)
750 {
751 const int primitive_global_id = lb.global_primitive_id(i);
752
753 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
754 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
755
756 if (!has_samples)
757 continue;
758
759 const basis::ElementBases &gbs = gbases_[e];
760 const basis::ElementBases &bs = bases_[e];
761
762 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
763
764 for (int n = 0; n < vals.jac_it.size(); ++n)
765 {
766 trafo = vals.jac_it[n].inverse();
767
768 if (displacement_prev.size() > 0)
769 {
770 assert(size_ == 2 || size_ == 3);
771 deform_mat.resize(size_, size_);
772 deform_mat.setZero();
773 for (const auto &b : vals.basis_values)
774 {
775 for (const auto &g : b.global)
776 {
777 for (int d = 0; d < size_; ++d)
778 {
779 deform_mat.row(d) += displacement_prev(g.index * size_ + d) * b.grad.row(n);
780 }
781 }
782 }
783
784 trafo += deform_mat;
785 }
786
787 normals.row(n) = normals.row(n) * trafo.inverse();
788 normals.row(n).normalize();
789 }
790 problem_.neumann_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, forces);
791
792 // UIState::ui_state().debug_data().add_points(vals.val, Eigen::RowVector3d(1,0,0));
793
794 for (long p = 0; p < weights.size(); ++p)
795 {
796 local_displacement.setZero();
797
798 for (size_t i = 0; i < vals.basis_values.size(); ++i)
799 {
800 const auto &vv = vals.basis_values[i];
801 assert(vv.val.size() == weights.size());
802 const double b_val = vv.val(p);
803
804 for (int d = 0; d < size_; ++d)
805 {
806 for (std::size_t ii = 0; ii < vv.global.size(); ++ii)
807 {
808 local_displacement(d) += (vv.global[ii].val * b_val) * displacement(vv.global[ii].index * size_ + d);
809 }
810 }
811 }
812
813 for (int d = 0; d < size_; ++d)
814 res -= forces(p, d) * local_displacement(d) * weights(p);
815 }
816 }
817 }
818
819 return res;
820 }
821
822 void RhsAssembler::compute_energy_hess(
823 const std::vector<int> &bounday_nodes,
824 const QuadratureOrders &resolution,
825 const std::vector<mesh::LocalBoundary> &local_neumann_boundary,
826 const Eigen::MatrixXd &displacement,
827 const double t,
828 const bool project_to_psd,
829 StiffnessMatrix &hess) const
830 {
831 hess.resize(n_basis_ * size_, n_basis_ * size_);
832 if (displacement.size() == 0)
833 return;
834
835 std::vector<Eigen::Triplet<double>> entries, entries_t;
836
837 ElementAssemblyValues vals;
838 Eigen::MatrixXd uv, samples, gtmp, rhs_fun, deform_mat, jac_mat, trafo;
839 Eigen::VectorXi global_primitive_ids;
840 Eigen::MatrixXd points, normals;
841 Eigen::VectorXd weights;
842 Eigen::MatrixXd local_hessian;
843
844 for (const auto &lb : local_neumann_boundary)
845 {
846 const int e = lb.element_id();
847 const basis::ElementBases &gbs = gbases_[e];
848 const basis::ElementBases &bs = bases_[e];
849
850 for (int i = 0; i < lb.size(); ++i)
851 {
852 const int primitive_global_id = lb.global_primitive_id(i);
853 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
854 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
855
856 if (!has_samples)
857 continue;
858
859 Eigen::MatrixXd reference_normals = normals;
860
861 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
862
863 std::vector<std::vector<Eigen::MatrixXd>> grad_normal;
864 for (int n = 0; n < vals.jac_it.size(); ++n)
865 {
866 trafo = vals.jac_it[n].inverse();
867
868 assert(size_ == 2 || size_ == 3);
869 deform_mat.resize(size_, size_);
870 deform_mat.setZero();
871 jac_mat.resize(size_, vals.basis_values.size());
872 int b_idx = 0;
873 for (const auto &b : vals.basis_values)
874 {
875 jac_mat.col(b_idx++) = b.grad.row(n);
876
877 for (const auto &g : b.global)
878 for (int d = 0; d < size_; ++d)
879 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
880 }
881
882 trafo += deform_mat;
883 trafo = trafo.inverse();
884
885 Eigen::VectorXd displaced_normal = normals.row(n) * trafo;
886 normals.row(n) = displaced_normal / displaced_normal.norm();
887
888 std::vector<Eigen::MatrixXd> grad;
889 {
890 Eigen::MatrixXd vec = -(jac_mat.transpose() * trafo * reference_normals.row(n).transpose());
891 // Gradient of the displaced normal computation
892 for (int k = 0; k < size_; ++k)
893 {
894 Eigen::MatrixXd grad_i(jac_mat.rows(), jac_mat.cols());
895 grad_i.setZero();
896 for (int m = 0; m < jac_mat.rows(); ++m)
897 for (int l = 0; l < jac_mat.cols(); ++l)
898 grad_i(m, l) = -(reference_normals.row(n) * trafo)(m) * (jac_mat.transpose() * trafo)(l, k);
899 grad.push_back(grad_i);
900 }
901 }
902
903 {
904 Eigen::MatrixXd normalization_chain_rule = (normals.row(n).transpose() * normals.row(n));
905 normalization_chain_rule = Eigen::MatrixXd::Identity(size_, size_) - normalization_chain_rule;
906 normalization_chain_rule /= displaced_normal.norm();
907
908 Eigen::VectorXd vec(size_);
909 b_idx = 0;
910 for (const auto &b : vals.basis_values)
911 {
912 for (int d = 0; d < size_; ++d)
913 {
914 for (int k = 0; k < size_; ++k)
915 vec(k) = grad[k](d, b_idx);
916 vec = normalization_chain_rule * vec;
917 for (int k = 0; k < size_; ++k)
918 grad[k](d, b_idx) = vec(k);
919 }
920 b_idx++;
921 }
922 }
923
924 grad_normal.push_back(grad);
925 }
926 Eigen::MatrixXd rhs_fun;
927 problem_.neumann_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, rhs_fun);
928
929 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
930
931 const bool is_pressure = problem_.is_boundary_pressure(mesh_.get_boundary_id(primitive_global_id));
932 if (!is_pressure)
933 continue;
934
935 local_hessian.setZero(vals.basis_values.size() * size_, vals.basis_values.size() * size_);
936
937 for (long n = 0; n < nodes.size(); ++n)
938 {
939 // const auto &b = bs.bases[nodes(n)];
940 const AssemblyValues &v = vals.basis_values[nodes(n)];
941 for (int d = 0; d < size_; ++d)
942 {
943 for (size_t g = 0; g < v.global.size(); ++g)
944 {
945 const int g_index = v.global[g].index * size_ + d;
946 const bool is_neumann = std::find(bounday_nodes.begin(), bounday_nodes.end(), g_index) == bounday_nodes.end();
947
948 if (is_neumann)
949 {
950 for (long ni = 0; ni < nodes.size(); ++ni)
951 {
952 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
953 for (int di = 0; di < size_; ++di)
954 {
955 for (size_t gi = 0; gi < vi.global.size(); ++gi)
956 {
957 const int gi_index = vi.global[gi].index * size_ + di;
958 double value = 0;
959
960 for (int q = 0; q < vals.jac_it.size(); ++q)
961 {
962 double pressure_val = rhs_fun.row(q).dot(normals.row(q));
963
964 // value += grad_normal[ni](d, nodes(ni) * size_ + di) * pressure_val * weights(q) * vi.val(q);
965 value += grad_normal[q][d](di, nodes(ni)) * pressure_val * weights(q) * vi.val(q);
966 }
967
968 value *= v.global[g].val;
969
970 const bool is_neumann_i = std::find(bounday_nodes.begin(), bounday_nodes.end(), gi_index) == bounday_nodes.end();
971
972 if (is_neumann_i)
973 {
974 local_hessian(nodes(n) * size_ + d, nodes(ni) * size_ + di) = value;
975 }
976 }
977 }
978 }
979 }
980 }
981 }
982 }
983
984 if (project_to_psd)
985 local_hessian = ipc::project_to_psd(local_hessian);
986
987 for (long n = 0; n < nodes.size(); ++n)
988 {
989 const AssemblyValues &v = vals.basis_values[nodes(n)];
990 for (int d = 0; d < size_; ++d)
991 {
992 for (size_t g = 0; g < v.global.size(); ++g)
993 {
994 const int g_index = v.global[g].index * size_ + d;
995
996 for (long ni = 0; ni < nodes.size(); ++ni)
997 {
998 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
999 for (int di = 0; di < size_; ++di)
1000 {
1001 for (size_t gi = 0; gi < vi.global.size(); ++gi)
1002 {
1003 const int gi_index = vi.global[gi].index * size_ + di;
1004 entries.push_back(Eigen::Triplet<double>(g_index, gi_index, local_hessian(nodes(n) * size_ + d, nodes(ni) * size_ + di)));
1005 }
1006 }
1007 }
1008 }
1009 }
1010 }
1011 }
1012 }
1013
1014 hess.setFromTriplets(entries.begin(), entries.end());
1015 }
1016 } // namespace assembler
1017} // 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
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:45
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::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::compute_energy_hess
void compute_energy_hess(const std::vector< int > &bounday_nodes, const QuadratureOrders &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:822
polyfem::assembler::RhsAssembler::set_bc
void set_bc(const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &bounday_nodes, const QuadratureOrders &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::mesh
const mesh::Mesh & mesh() const
Definition RhsAssembler.hpp:90
polyfem::assembler::RhsAssembler::gbases_
const std::vector< basis::ElementBases > & gbases_
basis functions associated with geometric mapping
Definition RhsAssembler.hpp:135
polyfem::assembler::RhsAssembler::size_
const int size_
dimension of problem
Definition RhsAssembler.hpp:133
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 QuadratureOrders &resolution, const double t) const
Definition RhsAssembler.cpp:667
polyfem::assembler::RhsAssembler::dirichlet_nodes_position_
const std::vector< RowVectorNd > & dirichlet_nodes_position_
Definition RhsAssembler.hpp:141
polyfem::assembler::RhsAssembler::ass_vals_cache_
const AssemblyValsCache & ass_vals_cache_
Definition RhsAssembler.hpp:136
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:131
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:134
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 QuadratureOrders &resolution, const std::vector< mesh::LocalBoundary > &local_neumann_boundary, const Eigen::MatrixXd &final_rhs, const double t, Eigen::MatrixXd &rhs) const
Definition RhsAssembler.cpp:636
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:140
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:426
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:40
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:481
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:152
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:490
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:715
polyfem::utils::BoundarySampler::boundary_quadrature
static bool boundary_quadrature(const mesh::LocalBoundary &local_boundary, const QuadratureOrders &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:644
generate_rotation_mtx.e
list e
Definition generate_rotation_mtx.py:49
generate_rotation_mtx.g
list g
Definition generate_rotation_mtx.py:56
p_bases.points
points
Definition p_bases.py:251
p_bases.nodes
str nodes
Definition p_bases.py:398
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:44
polyfem::QuadratureOrders
std::array< int, 2 > QuadratureOrders
Definition Types.hpp:19
polyfem::VectorNd
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 3, 1 > VectorNd
Definition Types.hpp:11
polyfem::json
nlohmann::json json
Definition Common.hpp:9
polyfem::StiffnessMatrix
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:24