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