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