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::set_bc(
491 const std::function<void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &df,
492 const std::function<void(const Eigen::MatrixXi &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, const Eigen::MatrixXd &, Eigen::MatrixXd &)> &nf,
493 const std::vector<LocalBoundary> &local_boundary, const std::vector<int> &bounday_nodes,
494 const int resolution, const std::vector<LocalBoundary> &local_neumann_boundary,
495 const Eigen::MatrixXd &displacement, const double t,
496 Eigen::MatrixXd &rhs) const
497 {
498 if (bc_method_ == "sample")
499 sample_bc(df, local_boundary, bounday_nodes, rhs);
500 else
501 lsq_bc(df, local_boundary, bounday_nodes, resolution, rhs);
502
503 if (bounday_nodes.size() > 0)
504 {
505 Eigen::MatrixXd tmp_val;
506 for (int n = 0; n < dirichlet_nodes_.size(); ++n)
507 {
508 const auto &n_id = dirichlet_nodes_[n];
509 const auto &pt = dirichlet_nodes_position_[n];
510
511 const int tag = mesh_.get_node_id(n_id);
512 problem_.dirichlet_nodal_value(mesh_, n_id, pt, t, tmp_val);
513 assert(tmp_val.size() == size_);
514
515 for (int d = 0; d < size_; ++d)
516 {
517 if (!problem_.is_nodal_dimension_dirichlet(n_id, tag, d))
518 continue;
519 const int g_index = n_id * size_ + d;
520 rhs(g_index) = tmp_val(d);
521 }
522 }
523 }
524
525 // Neumann
526 Eigen::MatrixXd uv, samples, gtmp, rhs_fun, deform_mat, trafo;
527 Eigen::VectorXi global_primitive_ids;
528 Eigen::MatrixXd points, normals;
529 Eigen::VectorXd weights;
530
531 ElementAssemblyValues vals;
532
533 for (const auto &lb : local_neumann_boundary)
534 {
535 const int e = lb.element_id();
536 const basis::ElementBases &gbs = gbases_[e];
537 const basis::ElementBases &bs = bases_[e];
538
539 for (int i = 0; i < lb.size(); ++i)
540 {
541 const int primitive_global_id = lb.global_primitive_id(i);
542
543 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
544 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
545
546 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
547
548 for (int n = 0; n < vals.jac_it.size(); ++n)
549 {
550 trafo = vals.jac_it[n].inverse();
551
552 if (displacement.size() > 0)
553 {
554 assert(size_ == 2 || size_ == 3);
555 deform_mat.resize(size_, size_);
556 deform_mat.setZero();
557 for (const auto &b : vals.basis_values)
558 {
559 for (const auto &g : b.global)
560 {
561 for (int d = 0; d < size_; ++d)
562 {
563 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
564 }
565 }
566 }
567
568 trafo += deform_mat;
569 }
570
571 normals.row(n) = normals.row(n) * trafo.inverse();
572 normals.row(n).normalize();
573 }
574
575 // problem_.neumann_bc(mesh_, global_primitive_ids, vals.val, t, rhs_fun);
576 nf(global_primitive_ids, uv, vals.val, normals, rhs_fun);
577
578 // UIState::ui_state().debug_data().add_points(vals.val, Eigen::RowVector3d(0,1,0));
579
580 for (int d = 0; d < size_; ++d)
581 rhs_fun.col(d) = rhs_fun.col(d).array() * weights.array();
582
583 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
584
585 for (long n = 0; n < nodes.size(); ++n)
586 {
587 // const auto &b = bs.bases[nodes(n)];
588 const AssemblyValues &v = vals.basis_values[nodes(n)];
589 for (int d = 0; d < size_; ++d)
590 {
591 const double rhs_value = (rhs_fun.col(d).array() * v.val.array()).sum();
592
593 for (size_t g = 0; g < v.global.size(); ++g)
594 {
595 const int g_index = v.global[g].index * size_ + d;
596 const bool is_neumann = std::find(bounday_nodes.begin(), bounday_nodes.end(), g_index) == bounday_nodes.end();
597
598 if (is_neumann)
599 {
600 rhs(g_index) += rhs_value * v.global[g].val;
601 }
602 }
603 }
604 }
605 }
606 }
607
608 // TODO add nodal neumann
609 }
610
611 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
612 {
613 set_bc(
614 [&](const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, Eigen::MatrixXd &val) {
615 problem_.dirichlet_bc(mesh_, global_ids, uv, pts, t, val);
616 },
617 [&](const Eigen::MatrixXi &global_ids, const Eigen::MatrixXd &uv, const Eigen::MatrixXd &pts, const Eigen::MatrixXd &normals, Eigen::MatrixXd &val) {
618 problem_.neumann_bc(mesh_, global_ids, uv, pts, normals, t, val);
619 },
620 local_boundary, bounday_nodes, resolution, local_neumann_boundary, displacement, t, rhs);
621
622 obstacle_.update_displacement(t, rhs);
623 }
624
625 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
626 {
627 if (problem_.is_constant_in_time())
628 {
629 rhs = final_rhs;
630 }
631 else
632 {
633 assemble(density, rhs, t);
634 rhs *= -1;
635
636 if (rhs.size() != final_rhs.size())
637 {
638 const int prev_size = rhs.size();
639 rhs.conservativeResize(final_rhs.size(), rhs.cols());
640 // Zero initial pressure
641 rhs.block(prev_size, 0, final_rhs.size() - prev_size, rhs.cols()).setZero();
642 rhs(rhs.size() - 1) = 0;
643 }
644
645 assert(rhs.size() == final_rhs.size());
646 }
647 }
648
649 double RhsAssembler::compute_energy(const Eigen::MatrixXd &displacement,
650 const Eigen::MatrixXd &displacement_prev,
651 const std::vector<LocalBoundary> &local_neumann_boundary,
652 const Density &density,
653 const int resolution,
654 const double t) const
655 {
656
657 double res = 0;
658
659 if (!problem_.is_rhs_zero())
660 {
661 auto storage = create_thread_storage(LocalThreadScalarStorage());
662 const int n_bases = int(bases_.size());
663
664 maybe_parallel_for(n_bases, [&](int start, int end, int thread_id) {
665 LocalThreadScalarStorage &local_storage = get_local_thread_storage(storage, thread_id);
666 VectorNd local_displacement(size_);
667 Eigen::MatrixXd forces;
668
669 for (int e = start; e < end; ++e)
670 {
671 ElementAssemblyValues &vals = local_storage.vals;
672 // vals.compute(e, mesh_.is_volume(), bases_[e], gbases_[e]);
673 ass_vals_cache_.compute(e, mesh_.is_volume(), bases_[e], gbases_[e], vals);
674
675 const Quadrature &quadrature = vals.quadrature;
676 const Eigen::VectorXd da = vals.det.array() * quadrature.weights.array();
677
678 problem_.rhs(assembler_, vals.val, t, forces);
679 assert(forces.rows() == da.size());
680 assert(forces.cols() == size_);
681
682 for (long p = 0; p < da.size(); ++p)
683 {
684 local_displacement.setZero();
685
686 for (size_t i = 0; i < vals.basis_values.size(); ++i)
687 {
688 const auto &bs = vals.basis_values[i];
689 assert(bs.val.size() == da.size());
690 const double b_val = bs.val(p);
691
692 for (int d = 0; d < size_; ++d)
693 {
694 for (std::size_t ii = 0; ii < bs.global.size(); ++ii)
695 {
696 local_displacement(d) += (bs.global[ii].val * b_val) * displacement(bs.global[ii].index * size_ + d);
697 }
698 }
699 }
700 // const double rho = problem_.is_time_dependent() ? density(vals.quadrature.points.row(p), vals.val.row(p), vals.element_id) : 1;
701 const double rho = density(vals.quadrature.points.row(p), vals.val.row(p), t, vals.element_id);
702
703 for (int d = 0; d < size_; ++d)
704 {
705 local_storage.val += forces(p, d) * local_displacement(d) * da(p) * rho;
706 // res += forces(p, d) * local_displacement(d) * da(p);
707 }
708 }
709 }
710 });
711
712 // Serially merge local storages
713 for (const LocalThreadScalarStorage &local_storage : storage)
714 res += local_storage.val;
715 }
716
717 VectorNd local_displacement(size_);
718 Eigen::MatrixXd forces;
719
720 ElementAssemblyValues vals;
721 // Neumann
722 Eigen::MatrixXd points, uv, normals, deform_mat, trafo;
723 Eigen::VectorXd weights;
724 Eigen::VectorXi global_primitive_ids;
725 for (const auto &lb : local_neumann_boundary)
726 {
727 const int e = lb.element_id();
728 const basis::ElementBases &gbs = gbases_[e];
729 const basis::ElementBases &bs = bases_[e];
730
731 for (int i = 0; i < lb.size(); ++i)
732 {
733 const int primitive_global_id = lb.global_primitive_id(i);
734
735 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
736 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
737
738 if (!has_samples)
739 continue;
740
741 const basis::ElementBases &gbs = gbases_[e];
742 const basis::ElementBases &bs = bases_[e];
743
744 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
745
746 for (int n = 0; n < vals.jac_it.size(); ++n)
747 {
748 trafo = vals.jac_it[n].inverse();
749
750 if (displacement_prev.size() > 0)
751 {
752 assert(size_ == 2 || size_ == 3);
753 deform_mat.resize(size_, size_);
754 deform_mat.setZero();
755 for (const auto &b : vals.basis_values)
756 {
757 for (const auto &g : b.global)
758 {
759 for (int d = 0; d < size_; ++d)
760 {
761 deform_mat.row(d) += displacement_prev(g.index * size_ + d) * b.grad.row(n);
762 }
763 }
764 }
765
766 trafo += deform_mat;
767 }
768
769 normals.row(n) = normals.row(n) * trafo.inverse();
770 normals.row(n).normalize();
771 }
772 problem_.neumann_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, forces);
773
774 // UIState::ui_state().debug_data().add_points(vals.val, Eigen::RowVector3d(1,0,0));
775
776 for (long p = 0; p < weights.size(); ++p)
777 {
778 local_displacement.setZero();
779
780 for (size_t i = 0; i < vals.basis_values.size(); ++i)
781 {
782 const auto &vv = vals.basis_values[i];
783 assert(vv.val.size() == weights.size());
784 const double b_val = vv.val(p);
785
786 for (int d = 0; d < size_; ++d)
787 {
788 for (std::size_t ii = 0; ii < vv.global.size(); ++ii)
789 {
790 local_displacement(d) += (vv.global[ii].val * b_val) * displacement(vv.global[ii].index * size_ + d);
791 }
792 }
793 }
794
795 for (int d = 0; d < size_; ++d)
796 res -= forces(p, d) * local_displacement(d) * weights(p);
797 }
798 }
799 }
800
801 return res;
802 }
803
804 void RhsAssembler::compute_energy_hess(
805 const std::vector<int> &bounday_nodes,
806 const int resolution,
807 const std::vector<mesh::LocalBoundary> &local_neumann_boundary,
808 const Eigen::MatrixXd &displacement,
809 const double t,
810 const bool project_to_psd,
811 StiffnessMatrix &hess) const
812 {
813 hess.resize(n_basis_ * size_, n_basis_ * size_);
814 if (displacement.size() == 0)
815 return;
816
817 std::vector<Eigen::Triplet<double>> entries, entries_t;
818
819 ElementAssemblyValues vals;
820 Eigen::MatrixXd uv, samples, gtmp, rhs_fun, deform_mat, jac_mat, trafo;
821 Eigen::VectorXi global_primitive_ids;
822 Eigen::MatrixXd points, normals;
823 Eigen::VectorXd weights;
824 Eigen::MatrixXd local_hessian;
825
826 for (const auto &lb : local_neumann_boundary)
827 {
828 const int e = lb.element_id();
829 const basis::ElementBases &gbs = gbases_[e];
830 const basis::ElementBases &bs = bases_[e];
831
832 for (int i = 0; i < lb.size(); ++i)
833 {
834 const int primitive_global_id = lb.global_primitive_id(i);
835
836 bool has_samples = utils::BoundarySampler::boundary_quadrature(lb, resolution, mesh_, i, false, uv, points, normals, weights);
837 global_primitive_ids.setConstant(weights.size(), primitive_global_id);
838
839 if (!has_samples)
840 continue;
841
842 Eigen::MatrixXd reference_normals = normals;
843
844 vals.compute(e, mesh_.is_volume(), points, bs, gbs);
845
846 std::vector<std::vector<Eigen::MatrixXd>> grad_normal;
847 for (int n = 0; n < vals.jac_it.size(); ++n)
848 {
849 trafo = vals.jac_it[n].inverse();
850
851 assert(size_ == 2 || size_ == 3);
852 deform_mat.resize(size_, size_);
853 deform_mat.setZero();
854 jac_mat.resize(size_, vals.basis_values.size());
855 int b_idx = 0;
856 for (const auto &b : vals.basis_values)
857 {
858 jac_mat.col(b_idx++) = b.grad.row(n);
859
860 for (const auto &g : b.global)
861 for (int d = 0; d < size_; ++d)
862 deform_mat.row(d) += displacement(g.index * size_ + d) * b.grad.row(n);
863 }
864
865 trafo += deform_mat;
866 trafo = trafo.inverse();
867
868 Eigen::VectorXd displaced_normal = normals.row(n) * trafo;
869 normals.row(n) = displaced_normal / displaced_normal.norm();
870
871 std::vector<Eigen::MatrixXd> grad;
872 {
873 Eigen::MatrixXd vec = -(jac_mat.transpose() * trafo * reference_normals.row(n).transpose());
874 // Gradient of the displaced normal computation
875 for (int k = 0; k < size_; ++k)
876 {
877 Eigen::MatrixXd grad_i(jac_mat.rows(), jac_mat.cols());
878 grad_i.setZero();
879 for (int m = 0; m < jac_mat.rows(); ++m)
880 for (int l = 0; l < jac_mat.cols(); ++l)
881 grad_i(m, l) = -(reference_normals.row(n) * trafo)(m) * (jac_mat.transpose() * trafo)(l, k);
882 grad.push_back(grad_i);
883 }
884 }
885
886 {
887 Eigen::MatrixXd normalization_chain_rule = (normals.row(n).transpose() * normals.row(n));
888 normalization_chain_rule = Eigen::MatrixXd::Identity(size_, size_) - normalization_chain_rule;
889 normalization_chain_rule /= displaced_normal.norm();
890
891 Eigen::VectorXd vec(size_);
892 b_idx = 0;
893 for (const auto &b : vals.basis_values)
894 {
895 for (int d = 0; d < size_; ++d)
896 {
897 for (int k = 0; k < size_; ++k)
898 vec(k) = grad[k](d, b_idx);
899 vec = normalization_chain_rule * vec;
900 for (int k = 0; k < size_; ++k)
901 grad[k](d, b_idx) = vec(k);
902 }
903 b_idx++;
904 }
905 }
906
907 grad_normal.push_back(grad);
908 }
909 Eigen::MatrixXd rhs_fun;
910 problem_.neumann_bc(mesh_, global_primitive_ids, uv, vals.val, normals, t, rhs_fun);
911
912 const auto nodes = bs.local_nodes_for_primitive(primitive_global_id, mesh_);
913
914 const bool is_pressure = problem_.is_boundary_pressure(mesh_.get_boundary_id(primitive_global_id));
915 if (!is_pressure)
916 continue;
917
918 local_hessian.setZero(vals.basis_values.size() * size_, vals.basis_values.size() * size_);
919
920 for (long n = 0; n < nodes.size(); ++n)
921 {
922 // const auto &b = bs.bases[nodes(n)];
923 const AssemblyValues &v = vals.basis_values[nodes(n)];
924 for (int d = 0; d < size_; ++d)
925 {
926 for (size_t g = 0; g < v.global.size(); ++g)
927 {
928 const int g_index = v.global[g].index * size_ + d;
929 const bool is_neumann = std::find(bounday_nodes.begin(), bounday_nodes.end(), g_index) == bounday_nodes.end();
930
931 if (is_neumann)
932 {
933 for (long ni = 0; ni < nodes.size(); ++ni)
934 {
935 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
936 for (int di = 0; di < size_; ++di)
937 {
938 for (size_t gi = 0; gi < vi.global.size(); ++gi)
939 {
940 const int gi_index = vi.global[gi].index * size_ + di;
941 double value = 0;
942
943 for (int q = 0; q < vals.jac_it.size(); ++q)
944 {
945 double pressure_val = rhs_fun.row(q).dot(normals.row(q));
946
947 // value += grad_normal[ni](d, nodes(ni) * size_ + di) * pressure_val * weights(q) * vi.val(q);
948 value += grad_normal[q][d](di, nodes(ni)) * pressure_val * weights(q) * vi.val(q);
949 }
950
951 value *= v.global[g].val;
952
953 const bool is_neumann_i = std::find(bounday_nodes.begin(), bounday_nodes.end(), gi_index) == bounday_nodes.end();
954
955 if (is_neumann_i)
956 {
957 local_hessian(nodes(n) * size_ + d, nodes(ni) * size_ + di) = value;
958 }
959 }
960 }
961 }
962 }
963 }
964 }
965 }
966
967 if (project_to_psd)
968 local_hessian = ipc::project_to_psd(local_hessian);
969
970 for (long n = 0; n < nodes.size(); ++n)
971 {
972 const AssemblyValues &v = vals.basis_values[nodes(n)];
973 for (int d = 0; d < size_; ++d)
974 {
975 for (size_t g = 0; g < v.global.size(); ++g)
976 {
977 const int g_index = v.global[g].index * size_ + d;
978
979 for (long ni = 0; ni < nodes.size(); ++ni)
980 {
981 const AssemblyValues &vi = vals.basis_values[nodes(ni)];
982 for (int di = 0; di < size_; ++di)
983 {
984 for (size_t gi = 0; gi < vi.global.size(); ++gi)
985 {
986 const int gi_index = vi.global[gi].index * size_ + di;
987 entries.push_back(Eigen::Triplet<double>(g_index, gi_index, local_hessian(nodes(n) * size_ + d, nodes(ni) * size_ + di)));
988 }
989 }
990 }
991 }
992 }
993 }
994 }
995 }
996
997 hess.setFromTriplets(entries.begin(), entries.end());
998 }
999 } // namespace assembler
1000} // 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: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:119
polyfem::assembler::RhsAssembler::size_
const int size_
dimension of problem
Definition RhsAssembler.hpp:117
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:625
polyfem::assembler::RhsAssembler::dirichlet_nodes_position_
const std::vector< RowVectorNd > & dirichlet_nodes_position_
Definition RhsAssembler.hpp:125
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:649
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:804
polyfem::assembler::RhsAssembler::ass_vals_cache_
const AssemblyValsCache & ass_vals_cache_
Definition RhsAssembler.hpp:120
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:115
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:118
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:124
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:475
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:151
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:484
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
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: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:11
polyfem::json
nlohmann::json json
Definition Common.hpp:9
polyfem::StiffnessMatrix
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22