PolyFEM
Loading...
Searching...
No Matches
OperatorSplittingSolver.cpp
Go to the documentation of this file.
1#include "OperatorSplittingSolver.hpp"
2#include <unsupported/Eigen/SparseExtra>
3
4#include <polysolve/linear/FEMSolver.hpp>
5
6#ifdef POLYFEM_WITH_OPENVDB
7#include <openvdb/openvdb.h>
8#endif
9
10namespace polyfem::solver
11{
12 using namespace assembler;
13
14 void OperatorSplittingSolver::initialize_grid(const mesh::Mesh &mesh,
15 const std::vector<basis::ElementBases> &gbases,
16 const std::vector<basis::ElementBases> &bases,
17 const double &density_dx)
18 {
19 resolution = density_dx;
20
21 grid_cell_num = RowVectorNd::Zero(dim);
22 for (int d = 0; d < dim; d++)
23 {
24 grid_cell_num(d) = ceil((max_domain(d) - min_domain(d)) / resolution);
25 }
26 if (dim == 2)
27 density = Eigen::VectorXd::Zero((grid_cell_num(0) + 1) * (grid_cell_num(1) + 1));
28 else
29 density = Eigen::VectorXd::Zero((grid_cell_num(0) + 1) * (grid_cell_num(1) + 1) * (grid_cell_num(2) + 1));
30 }
31
32 void OperatorSplittingSolver::initialize_mesh(const mesh::Mesh &mesh,
33 const int shape, const int n_el,
34 const std::vector<mesh::LocalBoundary> &local_boundary)
35 {
36 dim = mesh.dimension();
37 mesh.bounding_box(min_domain, max_domain);
38
39 const int size = local_boundary.size();
40 boundary_elem_id.reserve(size);
41 for (int e = 0; e < size; e++)
42 {
43 boundary_elem_id.push_back(local_boundary[e].element_id());
44 }
45
46 T.resize(n_el, shape);
47 for (int e = 0; e < n_el; e++)
48 {
49 for (int i = 0; i < shape; i++)
50 {
51 T(e, i) = mesh.cell_vertex(e, i);
52 }
53 }
54 V = Eigen::MatrixXd::Zero(mesh.n_vertices(), 3);
55 for (int i = 0; i < V.rows(); i++)
56 {
57 auto p = mesh.point(i);
58 for (int d = 0; d < dim; d++)
59 {
60 V(i, d) = p(d);
61 }
62 if (dim == 2)
63 V(i, 2) = 0;
64 }
65 }
66
67 void OperatorSplittingSolver::initialize_hashtable(const mesh::Mesh &mesh)
68 {
69 Eigen::MatrixXd p0, p1, p;
70 mesh.get_edges(p0, p1);
71 p = p0 - p1;
72 double avg_edge_length = p.rowwise().norm().mean();
73
74 long total_cell_num = 1;
75 hash_table_cell_num.setZero(dim);
76 for (int d = 0; d < dim; d++)
77 {
78 hash_table_cell_num(d) = (long)std::round((max_domain(d) - min_domain(d)) / avg_edge_length) * 4;
79 logger().debug("hash grid in {} dimension: {}", d, hash_table_cell_num(d));
80 total_cell_num *= hash_table_cell_num(d);
81 }
82 hash_table.resize(total_cell_num);
83 for (int e = 0; e < T.rows(); e++)
84 {
85 Eigen::VectorXd min_ = V.row(T(e, 0));
86 Eigen::VectorXd max_ = min_;
87
88 for (int i = 1; i < T.cols(); i++)
89 {
90 Eigen::VectorXd p = V.row(T(e, i));
91 min_ = min_.cwiseMin(p);
92 max_ = max_.cwiseMax(p);
93 }
94
95 Eigen::Matrix<long, Eigen::Dynamic, 1, Eigen::ColMajor, 3, 1> min_int(dim), max_int(dim);
96
97 for (int d = 0; d < dim; d++)
98 {
99 double temp = hash_table_cell_num(d) / (max_domain(d) - min_domain(d));
100 min_int(d) = floor((min_(d) * (1 - 1e-14) - min_domain(d)) * temp);
101 max_int(d) = ceil((max_(d) * (1 + 1e-14) - min_domain(d)) * temp);
102
103 min_int = min_int.cwiseMax(0);
104 max_int = max_int.cwiseMin(hash_table_cell_num);
105 }
106
107 for (long x = min_int(0); x < max_int(0); x++)
108 {
109 for (long y = min_int(1); y < max_int(1); y++)
110 {
111 if (dim == 2)
112 {
113 long idx = x + y * hash_table_cell_num(0);
114 hash_table[idx].push_back(e);
115 }
116 else
117 {
118 for (long z = min_int(2); z < max_int(2); z++)
119 {
120 long idx = x + (y + z * hash_table_cell_num(1)) * hash_table_cell_num(0);
121 hash_table[idx].push_back(e);
122 }
123 }
124 }
125 }
126 }
127
128 float average_intersection_num = 0;
129 int max_intersection_num = 0;
130 for (int i = 0; i < hash_table.size(); i++)
131 {
132 average_intersection_num += hash_table[i].size();
133 (hash_table[i].size() > max_intersection_num) ? (max_intersection_num = hash_table[i].size()) : 1;
134 }
135 average_intersection_num /= hash_table.size();
136 logger().debug("average intersection number for hash grid: {}", average_intersection_num);
137 logger().debug("max intersection number for hash grid: {}", max_intersection_num);
138 }
139
140 void OperatorSplittingSolver::initialize_solver(const mesh::Mesh &mesh,
141 const int shape_, const int n_el_,
142 const std::vector<mesh::LocalBoundary> &local_boundary,
143 const std::vector<int> &bnd_nodes)
144 {
145 shape = shape_;
146 n_el = n_el_;
147 boundary_nodes = bnd_nodes;
148
149 initialize_mesh(mesh, shape, n_el, local_boundary);
150 initialize_hashtable(mesh);
151 }
152
153 OperatorSplittingSolver::OperatorSplittingSolver(const mesh::Mesh &mesh,
154 const int shape, const int n_el,
155 const std::vector<mesh::LocalBoundary> &local_boundary,
156 const std::vector<int> &bnd_nodes)
157 {
158 initialize_solver(mesh, shape, n_el, local_boundary, bnd_nodes);
159 }
160
161 OperatorSplittingSolver::OperatorSplittingSolver(const mesh::Mesh &mesh,
162 const int shape, const int n_el,
163 const std::vector<mesh::LocalBoundary> &local_boundary,
164 const std::vector<int> &boundary_nodes_,
165 const std::vector<int> &pressure_boundary_nodes,
166 const std::vector<int> &bnd_nodes,
167 const StiffnessMatrix &mass,
168 const StiffnessMatrix &stiffness_viscosity,
169 const StiffnessMatrix &stiffness_velocity,
170 const StiffnessMatrix &mass_velocity,
171 const double &dt,
172 const double &viscosity_,
173 const json &params)
174 {
175 initialize_solver(mesh, shape, n_el, local_boundary, bnd_nodes);
176
177 logger().info("Prefactorization begins...");
178
179 solver_mass = polysolve::linear::Solver::create(params, logger());
180 // if (solver_type == "Pardiso" || solver_type == "Eigen::SimplicialLDLT" || solver_type == "Eigen::SparseLU")
181 {
182 StiffnessMatrix mat1 = mass_velocity;
183 prefactorize(*solver_mass, mat1, boundary_nodes_, mat1.rows(), "");
184 }
185
186 mat_diffusion = mass + viscosity_ * dt * stiffness_viscosity;
187
188 solver_diffusion = polysolve::linear::Solver::create(params, logger());
189 // if (solver_type == "Pardiso" || solver_type == "Eigen::SimplicialLDLT" || solver_type == "Eigen::SparseLU")
190 {
191 StiffnessMatrix mat1 = mat_diffusion;
192 prefactorize(*solver_diffusion, mat1, bnd_nodes, mat1.rows(), "");
193 }
194
195 if (pressure_boundary_nodes.size() == 0)
196 mat_projection.resize(stiffness_velocity.rows() + 1, stiffness_velocity.cols() + 1);
197 else
198 mat_projection.resize(stiffness_velocity.rows(), stiffness_velocity.cols());
199
200 std::vector<Eigen::Triplet<double>> coefficients;
201 coefficients.reserve(stiffness_velocity.nonZeros() + 2 * stiffness_velocity.rows());
202
203 for (int i = 0; i < stiffness_velocity.outerSize(); i++)
204 {
205 for (StiffnessMatrix::InnerIterator it(stiffness_velocity, i); it; ++it)
206 {
207 coefficients.emplace_back(it.row(), it.col(), it.value());
208 }
209 }
210
211 if (pressure_boundary_nodes.size() == 0)
212 {
213 const double val = 1. / (mat_projection.rows() - 1);
214 for (int i = 0; i < mat_projection.rows() - 1; i++)
215 {
216 coefficients.emplace_back(i, mat_projection.cols() - 1, val);
217 coefficients.emplace_back(mat_projection.rows() - 1, i, val);
218 }
219 coefficients.emplace_back(mat_projection.rows() - 1, mat_projection.cols() - 1, 2);
220 }
221
222 mat_projection.setFromTriplets(coefficients.begin(), coefficients.end());
223 solver_projection = polysolve::linear::Solver::create(params, logger());
224 logger().info("{}...", solver_projection->name());
225 // if (solver_type == "Pardiso" || solver_type == "Eigen::SimplicialLDLT" || solver_type == "Eigen::SparseLU")
226 {
227 StiffnessMatrix mat2 = mat_projection;
228 prefactorize(*solver_projection, mat2, pressure_boundary_nodes, mat2.rows(), "");
229 }
230 logger().info("Prefactorization ends!");
231 }
232
233 int OperatorSplittingSolver::handle_boundary_advection(RowVectorNd &pos)
234 {
235 double dist = 1e10;
236 int idx = -1, local_idx = -1;
237 const int size = boundary_elem_id.size();
238#ifdef POLYFEM_WITH_TBB
239 tbb::parallel_for(0, size, 1, [&](int e)
240#else
241 for (int e = 0; e < size; e++)
242#endif
243 {
244 int elem_idx = boundary_elem_id[e];
245
246 for (int i = 0; i < shape; i++)
247 {
248 double dist_ = 0;
249 for (int d = 0; d < dim; d++)
250 {
251 dist_ += pow(pos(d) - V(T(elem_idx, i), d), 2);
252 }
253 dist_ = sqrt(dist_);
254 if (dist_ < dist)
255 {
256 dist = dist_;
257 idx = elem_idx;
258 local_idx = i;
259 }
260 }
261 }
262#ifdef POLYFEM_WITH_TBB
263 );
264#endif
265 for (int d = 0; d < dim; d++)
266 pos(d) = V(T(idx, local_idx), d);
267 return idx;
268 }
269
270 int OperatorSplittingSolver::trace_back(const std::vector<basis::ElementBases> &gbases,
271 const std::vector<basis::ElementBases> &bases,
272 const RowVectorNd &pos_1,
273 const RowVectorNd &vel_1,
274 RowVectorNd &pos_2,
275 RowVectorNd &vel_2,
276 Eigen::MatrixXd &local_pos,
277 const Eigen::MatrixXd &sol,
278 const double dt)
279 {
280 pos_2 = pos_1 - vel_1 * dt;
281
282 return interpolator(gbases, bases, pos_2, vel_2, local_pos, sol);
283 }
284
285 int OperatorSplittingSolver::interpolator(const std::vector<basis::ElementBases> &gbases,
286 const std::vector<basis::ElementBases> &bases,
287 const RowVectorNd &pos,
288 RowVectorNd &vel,
289 Eigen::MatrixXd &local_pos,
290 const Eigen::MatrixXd &sol)
291 {
292 bool insideDomain = true;
293
294 int new_elem;
295 if ((new_elem = search_cell(gbases, pos, local_pos)) == -1)
296 {
297 insideDomain = false;
298 RowVectorNd pos_ = pos;
299 new_elem = handle_boundary_advection(pos_);
300 calculate_local_pts(gbases[new_elem], new_elem, pos_, local_pos);
301 }
302
303 // interpolation
304 vel = RowVectorNd::Zero(dim);
305 assembler::ElementAssemblyValues vals;
306 vals.compute(new_elem, dim == 3, local_pos, bases[new_elem], gbases[new_elem]);
307 for (int d = 0; d < dim; d++)
308 {
309 for (int i = 0; i < vals.basis_values.size(); i++)
310 {
311 vel(d) += vals.basis_values[i].val(0) * sol(bases[new_elem].bases[i].global()[0].index * dim + d);
312 }
313 }
314 if (insideDomain)
315 return new_elem;
316 else
317 return -1;
318 }
319
320 void OperatorSplittingSolver::interpolator(const RowVectorNd &pos, double &val)
321 {
322 val = 0;
323
324 Eigen::Matrix<long, Eigen::Dynamic, 1> int_pos(dim);
325 Eigen::MatrixXd weights(2, dim);
326 for (int d = 0; d < dim; d++)
327 {
328 int_pos(d) = floor((pos(d) - min_domain(d)) / resolution);
329 if (int_pos(d) < 0 || int_pos(d) >= grid_cell_num(d))
330 return;
331 weights(1, d) = (pos(d) - min_domain(d)) / resolution - int_pos(d);
332 weights(0, d) = 1 - weights(1, d);
333 }
334
335 for (int d1 = 0; d1 < 2; d1++)
336 for (int d2 = 0; d2 < 2; d2++)
337 {
338 if (dim == 2)
339 {
340 const long idx = (int_pos(0) + d1) + (int_pos(1) + d2) * (grid_cell_num(0) + 1);
341 val += density(idx) * weights(d1, 0) * weights(d2, 1);
342 }
343 else
344 {
345 for (int d3 = 0; d3 < 2; d3++)
346 {
347 const long idx = (int_pos(0) + d1) + (int_pos(1) + d2 + (int_pos(2) + d3) * (grid_cell_num(1) + 1)) * (grid_cell_num(0) + 1);
348 val += density(idx) * weights(d1, 0) * weights(d2, 1) * weights(d3, 2);
349 }
350 }
351 }
352 }
353
354 void OperatorSplittingSolver::advection(const mesh::Mesh &mesh,
355 const std::vector<basis::ElementBases> &gbases,
356 const std::vector<basis::ElementBases> &bases,
357 Eigen::MatrixXd &sol,
358 const double dt,
359 const Eigen::MatrixXd &local_pts,
360 const int order,
361 const int RK)
362 {
363 // to store new velocity
364 Eigen::MatrixXd new_sol = Eigen::MatrixXd::Zero(sol.size(), 1);
365 // number of FEM nodes
366 const int n_vert = sol.size() / dim;
367 Eigen::VectorXi traversed = Eigen::VectorXi::Zero(n_vert);
368
369#ifdef POLYFEM_WITH_TBB
370 tbb::parallel_for(0, n_el, 1, [&](int e)
371#else
372 for (int e = 0; e < n_el; ++e)
373#endif
374 {
375 // to compute global position with barycentric coordinate
376 Eigen::MatrixXd mapped;
377 gbases[e].eval_geom_mapping(local_pts, mapped);
378
379 for (int i = 0; i < local_pts.rows(); i++)
380 {
381 // global index of this FEM node
382 int global = bases[e].bases[i].global()[0].index;
383
384 if (traversed(global))
385 continue;
386 traversed(global) = 1;
387
388 // velocity of this FEM node
389 RowVectorNd vel_ = sol.block(global * dim, 0, dim, 1).transpose();
390
391 // global position of this FEM node
392 RowVectorNd pos_ = RowVectorNd::Zero(1, dim);
393 for (int d = 0; d < dim; d++)
394 pos_(d) = mapped(i, d) - vel_(d) * dt;
395
396 Eigen::MatrixXd local_pos;
397 interpolator(gbases, bases, pos_, vel_, local_pos, sol);
398
399 new_sol.block(global * dim, 0, dim, 1) = vel_.transpose();
400 }
401 }
402#ifdef POLYFEM_WITH_TBB
403 );
404#endif
405 sol.swap(new_sol);
406 }
407
408 void OperatorSplittingSolver::advect_density_exact(const std::vector<basis::ElementBases> &gbases,
409 const std::vector<basis::ElementBases> &bases,
410 const std::shared_ptr<assembler::Problem> problem,
411 const double t,
412 const double dt,
413 const int RK)
414 {
415 Eigen::VectorXd new_density = Eigen::VectorXd::Zero(density.size());
416 const int Nx = grid_cell_num(0);
417#ifdef POLYFEM_WITH_TBB
418 tbb::parallel_for(0, Nx + 1, 1, [&](int i)
419#else
420 for (int i = 0; i <= Nx; i++)
421#endif
422 {
423 for (int j = 0; j <= grid_cell_num(1); j++)
424 {
425 if (dim == 2)
426 {
427 Eigen::MatrixXd pos(1, dim);
428 pos(0) = i * resolution + min_domain(0);
429 pos(1) = j * resolution + min_domain(1);
430 const long idx = i + (long)j * (grid_cell_num(0) + 1);
431
432 Eigen::MatrixXd vel1, pos_;
433 problem->exact(pos, t, vel1);
434 if (RK > 1)
435 {
436 Eigen::MatrixXd vel2, vel3;
437 problem->exact(pos - 0.5 * dt * vel1, t, vel2);
438 problem->exact(pos - 0.75 * dt * vel2, t, vel3);
439 pos_ = pos - (2 * vel1 + 3 * vel2 + 4 * vel3) * dt / 9;
440 }
441 else
442 {
443 pos_ = pos - vel1 * dt;
444 }
445 interpolator(pos_, new_density[idx]);
446 }
447 else
448 {
449 for (int k = 0; k <= grid_cell_num(2); k++)
450 {
451 RowVectorNd pos(1, dim);
452 pos(0) = i * resolution + min_domain(0);
453 pos(1) = j * resolution + min_domain(1);
454 pos(2) = k * resolution + min_domain(2);
455 const long idx = i + (j + (long)k * (grid_cell_num(1) + 1)) * (grid_cell_num(0) + 1);
456
457 Eigen::MatrixXd vel1, pos_;
458 problem->exact(pos, t, vel1);
459 if (RK > 1)
460 {
461 Eigen::MatrixXd vel2, vel3;
462 problem->exact(pos - 0.5 * dt * vel1, t, vel2);
463 problem->exact(pos - 0.75 * dt * vel2, t, vel3);
464 pos_ = pos - (2 * vel1 + 3 * vel2 + 4 * vel3) * dt / 9;
465 }
466 else
467 {
468 pos_ = pos - vel1 * dt;
469 }
470 interpolator(pos_, new_density[idx]);
471 }
472 }
473 }
474 }
475#ifdef POLYFEM_WITH_TBB
476 );
477#endif
478 density.swap(new_density);
479 }
480
481 void OperatorSplittingSolver::advect_density(const std::vector<basis::ElementBases> &gbases,
482 const std::vector<basis::ElementBases> &bases,
483 const Eigen::MatrixXd &sol,
484 const double dt,
485 const int RK)
486 {
487 Eigen::VectorXd new_density = Eigen::VectorXd::Zero(density.size());
488 const int Nx = grid_cell_num(0);
489#ifdef POLYFEM_WITH_TBB
490 tbb::parallel_for(0, Nx + 1, 1, [&](int i)
491#else
492 for (int i = 0; i <= Nx; i++)
493#endif
494 {
495 for (int j = 0; j <= grid_cell_num(1); j++)
496 {
497 Eigen::MatrixXd local_pos;
498 if (dim == 2)
499 {
500 RowVectorNd pos(1, dim);
501 pos(0) = i * resolution + min_domain(0);
502 pos(1) = j * resolution + min_domain(1);
503 const long idx = i + (long)j * (grid_cell_num(0) + 1);
504
505 RowVectorNd vel1, pos_;
506 interpolator(gbases, bases, pos, vel1, local_pos, sol);
507 if (RK > 1)
508 {
509 RowVectorNd vel2, vel3;
510 interpolator(gbases, bases, pos - 0.5 * dt * vel1, vel2, local_pos, sol);
511 interpolator(gbases, bases, pos - 0.75 * dt * vel2, vel3, local_pos, sol);
512 pos_ = pos - (2 * vel1 + 3 * vel2 + 4 * vel3) * dt / 9;
513 }
514 else
515 {
516 pos_ = pos - vel1 * dt;
517 }
518 interpolator(pos_, new_density[idx]);
519 }
520 else
521 {
522 for (int k = 0; k <= grid_cell_num(2); k++)
523 {
524 RowVectorNd pos(1, dim);
525 pos(0) = i * resolution + min_domain(0);
526 pos(1) = j * resolution + min_domain(1);
527 pos(2) = k * resolution + min_domain(2);
528 const long idx = i + (j + (long)k * (grid_cell_num(1) + 1)) * (grid_cell_num(0) + 1);
529
530 RowVectorNd vel1, pos_;
531 interpolator(gbases, bases, pos, vel1, local_pos, sol);
532 if (RK > 1)
533 {
534 RowVectorNd vel2, vel3;
535 interpolator(gbases, bases, pos - 0.5 * dt * vel1, vel2, local_pos, sol);
536 interpolator(gbases, bases, pos - 0.75 * dt * vel2, vel3, local_pos, sol);
537 pos_ = pos - (2 * vel1 + 3 * vel2 + 4 * vel3) * dt / 9;
538 }
539 else
540 {
541 pos_ = pos - vel1 * dt;
542 }
543 interpolator(pos_, new_density[idx]);
544 }
545 }
546 }
547 }
548#ifdef POLYFEM_WITH_TBB
549 );
550#endif
551 density.swap(new_density);
552 }
553
554 void OperatorSplittingSolver::advection_FLIP(const mesh::Mesh &mesh, const std::vector<basis::ElementBases> &gbases, const std::vector<basis::ElementBases> &bases, Eigen::MatrixXd &sol, const double dt, const Eigen::MatrixXd &local_pts, const int order)
555 {
556 const int ppe = shape; // particle per element
557 const double FLIPRatio = 1;
558 // initialize or resample particles and update velocity via g2p
559 if (position_particle.empty())
560 {
561 // initialize particles
562 position_particle.resize(n_el * ppe);
563 velocity_particle.resize(n_el * ppe);
564 cellI_particle.resize(n_el * ppe);
565#ifdef POLYFEM_WITH_TBB
566 tbb::parallel_for(0, n_el, 1, [&](int e)
567#else
568 for (int e = 0; e < n_el; ++e)
569#endif
570 {
571 // sample particle in element e
572 Eigen::MatrixXd local_pts_particle;
573 local_pts_particle.setRandom(ppe, dim);
574 local_pts_particle.array() += 1;
575 local_pts_particle.array() /= 2;
576 for (int ppeI = 0; ppeI < ppe; ++ppeI)
577 {
578 if (shape == 3 && dim == 2 && local_pts_particle.row(ppeI).sum() > 1)
579 {
580 double x = 1 - local_pts_particle(ppeI, 1);
581 local_pts_particle(ppeI, 1) = 1 - local_pts_particle(ppeI, 0);
582 local_pts_particle(ppeI, 0) = x;
583 // TODO: dim == 3
584 }
585 }
586
587 for (int i = 0; i < shape; ++i)
588 cellI_particle[e * ppe + i] = e;
589
590 // compute global position and velocity of particles
591 // construct interpolant (linear for position)
592 Eigen::MatrixXd mapped;
593 gbases[e].eval_geom_mapping(local_pts_particle, mapped);
594 // construct interpolant (for velocity)
595 assembler::ElementAssemblyValues vals;
596 vals.compute(e, dim == 3, local_pts_particle, bases[e], gbases[e]); // possibly higher-order
597 for (int j = 0; j < ppe; ++j)
598 {
599 position_particle[ppe * e + j].setZero(1, dim);
600 for (int d = 0; d < dim; d++)
601 {
602 position_particle[ppe * e + j](d) = mapped(j, d);
603 }
604
605 velocity_particle[e * ppe + j].setZero(1, dim);
606 for (int i = 0; i < vals.basis_values.size(); ++i)
607 {
608 velocity_particle[e * ppe + j] += vals.basis_values[i].val(j) * sol.block(bases[e].bases[i].global()[0].index * dim, 0, dim, 1).transpose();
609 }
610 }
611 }
612#ifdef POLYFEM_WITH_TBB
613 );
614#endif
615 }
616 else
617 {
618 // resample particles
619 // count particle per cell
620 std::vector<int> counter(n_el, 0);
621 std::vector<int> redundantPI;
622 std::vector<bool> isRedundant(cellI_particle.size(), false);
623 for (int pI = 0; pI < cellI_particle.size(); ++pI)
624 {
625 if (cellI_particle[pI] < 0)
626 {
627 redundantPI.emplace_back(pI);
628 isRedundant[pI] = true;
629 }
630 else
631 {
632 ++counter[cellI_particle[pI]];
633 if (counter[cellI_particle[pI]] > ppe)
634 {
635 redundantPI.emplace_back(pI);
636 isRedundant[pI] = true;
637 }
638 }
639 }
640 // g2p -- update velocity
641#ifdef POLYFEM_WITH_TBB
642 tbb::parallel_for(0, (int)cellI_particle.size(), 1, [&](int pI)
643#else
644 for (int pI = 0; pI < cellI_particle.size(); ++pI)
645#endif
646 {
647 if (!isRedundant[pI])
648 {
649 int e = cellI_particle[pI];
650 Eigen::MatrixXd local_pts_particle;
651 calculate_local_pts(gbases[e], e, position_particle[pI], local_pts_particle);
652
653 assembler::ElementAssemblyValues vals;
654 vals.compute(e, dim == 3, local_pts_particle, bases[e], gbases[e]); // possibly higher-order
655 RowVectorNd FLIPdVel, PICVel;
656 FLIPdVel.setZero(1, dim);
657 PICVel.setZero(1, dim);
658 for (int i = 0; i < vals.basis_values.size(); ++i)
659 {
660 FLIPdVel += vals.basis_values[i].val(0) * (sol.block(bases[e].bases[i].global()[0].index * dim, 0, dim, 1) - new_sol.block(bases[e].bases[i].global()[0].index * dim, 0, dim, 1)).transpose();
661 PICVel += vals.basis_values[i].val(0) * sol.block(bases[e].bases[i].global()[0].index * dim, 0, dim, 1).transpose();
662 }
663 velocity_particle[pI] = (1.0 - FLIPRatio) * PICVel + FLIPRatio * (velocity_particle[pI] + FLIPdVel);
664 }
665 }
666#ifdef POLYFEM_WITH_TBB
667 );
668#endif
669 // resample
670 for (int e = 0; e < n_el; ++e)
671 {
672 if (counter[e] >= ppe)
673 {
674 continue;
675 }
676
677 while (counter[e] < ppe)
678 {
679 int pI = redundantPI.back();
680 redundantPI.pop_back();
681
682 cellI_particle[pI] = e;
683
684 // sample particle in element e
685 Eigen::MatrixXd local_pts_particle;
686 local_pts_particle.setRandom(1, dim);
687 local_pts_particle.array() += 1;
688 local_pts_particle.array() /= 2;
689 if (shape == 3 && dim == 2 && local_pts_particle.sum() > 1)
690 {
691 double x = 1 - local_pts_particle(0, 1);
692 local_pts_particle(0, 1) = 1 - local_pts_particle(0, 0);
693 local_pts_particle(0, 0) = x;
694 // TODO: dim == 3
695 }
696
697 // compute global position and velocity of particles
698 // construct interpolant (linear for position)
699 Eigen::MatrixXd mapped;
700 gbases[e].eval_geom_mapping(local_pts_particle, mapped);
701 for (int d = 0; d < dim; d++)
702 position_particle[pI](d) = mapped(0, d);
703
704 // construct interpolant (for velocity)
705 assembler::ElementAssemblyValues vals;
706 vals.compute(e, dim == 3, local_pts_particle, bases[e], gbases[e]); // possibly higher-order
707 velocity_particle[pI].setZero(1, dim);
708 for (int i = 0; i < vals.basis_values.size(); ++i)
709 {
710 velocity_particle[pI] += vals.basis_values[i].val(0) * sol.block(bases[e].bases[i].global()[0].index * dim, 0, dim, 1).transpose();
711 }
712
713 ++counter[e];
714 }
715 }
716 }
717
718 // advect
719 std::vector<assembler::ElementAssemblyValues> velocity_interpolator(ppe * n_el);
720#ifdef POLYFEM_WITH_TBB
721 tbb::parallel_for(0, (int)(ppe * n_el), 1, [&](int pI)
722#else
723 for (int pI = 0; pI < ppe * n_el; ++pI)
724#endif
725 {
726 // update particle position via advection
727 RowVectorNd newvel;
728 Eigen::MatrixXd local_pos;
729 cellI_particle[pI] = trace_back(gbases, bases, position_particle[pI], velocity_particle[pI],
730 position_particle[pI], newvel, local_pos, sol, -dt);
731
732 // RK3:
733 // RowVectorNd bypass, vel2, vel3;
734 // trace_back( gbases, bases, position_particle[pI], velocity_particle[pI],
735 // bypass, vel2, sol, -0.5 * dt);
736 // trace_back( gbases, bases, position_particle[pI], vel2,
737 // bypass, vel3, sol, -0.75 * dt);
738 // trace_back( gbases, bases, position_particle[pI],
739 // 2 * velocity_particle[pI] + 3 * vel2 + 4 * vel3,
740 // position_particle[pI], bypass, sol, -dt / 9);
741
742 // prepare P2G
743 if (cellI_particle[pI] >= 0)
744 {
745 // construct interpolator (always linear for P2G, can use gaussian or bspline later)
746 velocity_interpolator[pI].compute(cellI_particle[pI], dim == 3, local_pos,
747 gbases[cellI_particle[pI]], gbases[cellI_particle[pI]]);
748 }
749 }
750#ifdef POLYFEM_WITH_TBB
751 );
752#endif
753
754 // P2G
755 new_sol = Eigen::MatrixXd::Zero(sol.size(), 1);
756 new_sol_w = Eigen::MatrixXd::Zero(sol.size() / dim, 1);
757 new_sol_w.array() += 1e-13;
758 for (int pI = 0; pI < ppe * n_el; ++pI)
759 {
760 int cellI = cellI_particle[pI];
761 if (cellI == -1)
762 {
763 continue;
764 }
765 assembler::ElementAssemblyValues &vals = velocity_interpolator[pI];
766 for (int i = 0; i < vals.basis_values.size(); ++i)
767 {
768 new_sol.block(bases[cellI].bases[i].global()[0].index * dim, 0, dim, 1) +=
769 vals.basis_values[i].val(0) * velocity_particle[pI].transpose();
770 new_sol_w(bases[cellI].bases[i].global()[0].index) += vals.basis_values[i].val(0);
771 }
772 }
773 // TODO: need to add up boundary velocities and weights because of perodic BC
774
775#ifdef POLYFEM_WITH_TBB
776 tbb::parallel_for(0, (int)new_sol.rows() / dim, 1, [&](int i)
777#else
778 for (int i = 0; i < new_sol.rows() / dim; ++i)
779#endif
780 {
781 new_sol.block(i * dim, 0, dim, 1) /= new_sol_w(i, 0);
782 sol.block(i * dim, 0, dim, 1) = new_sol.block(i * dim, 0, dim, 1);
783 }
784#ifdef POLYFEM_WITH_TBB
785 );
786#endif
787 }
788
789 void OperatorSplittingSolver::advection_PIC(const mesh::Mesh &mesh, const std::vector<basis::ElementBases> &gbases, const std::vector<basis::ElementBases> &bases, Eigen::MatrixXd &sol, const double dt, const Eigen::MatrixXd &local_pts, const int order)
790 {
791 // to store new velocity and weights for particle grid transfer
792 Eigen::MatrixXd new_sol = Eigen::MatrixXd::Zero(sol.size(), 1);
793 Eigen::MatrixXd new_sol_w = Eigen::MatrixXd::Zero(sol.size() / dim, 1);
794 new_sol_w.array() += 1e-13;
795
796 const int ppe = shape; // particle per element
797 std::vector<assembler::ElementAssemblyValues> velocity_interpolator(ppe * n_el);
798 position_particle.resize(ppe * n_el);
799 velocity_particle.resize(ppe * n_el);
800 cellI_particle.resize(ppe * n_el);
801#ifdef POLYFEM_WITH_TBB
802 tbb::parallel_for(0, n_el, 1, [&](int e)
803#else
804 for (int e = 0; e < n_el; ++e)
805#endif
806 {
807 // resample particle in element e
808 Eigen::MatrixXd local_pts_particle;
809 local_pts_particle.setRandom(ppe, dim);
810 local_pts_particle.array() += 1;
811 local_pts_particle.array() /= 2;
812
813 // geometry vertices of element e
814 std::vector<RowVectorNd> vert(shape);
815 for (int i = 0; i < shape; ++i)
816 {
817 vert[i] = mesh.point(mesh.cell_vertex(e, i));
818 }
819
820 // construct interpolant (linear for position)
821 assembler::ElementAssemblyValues gvals;
822 gvals.compute(e, dim == 3, local_pts_particle, gbases[e], gbases[e]);
823
824 // compute global position of particles
825 for (int i = 0; i < ppe; ++i)
826 {
827 position_particle[ppe * e + i].setZero(1, dim);
828 for (int j = 0; j < shape; ++j)
829 {
830 position_particle[ppe * e + i] += gvals.basis_values[j].val(i) * vert[j];
831 }
832 }
833
834 // compute velocity
835 assembler::ElementAssemblyValues vals;
836 vals.compute(e, dim == 3, local_pts_particle, bases[e], gbases[e]); // possibly higher-order
837 for (int j = 0; j < ppe; ++j)
838 {
839 velocity_particle[e * ppe + j].setZero(1, dim);
840 for (int i = 0; i < vals.basis_values.size(); ++i)
841 {
842 velocity_particle[e * ppe + j] += vals.basis_values[i].val(j) * sol.block(bases[e].bases[i].global()[0].index * dim, 0, dim, 1).transpose();
843 }
844 }
845
846 // update particle position via advection
847 for (int j = 0; j < ppe; ++j)
848 {
849 RowVectorNd newvel;
850 Eigen::MatrixXd local_pos;
851 cellI_particle[ppe * e + j] = trace_back(gbases, bases, position_particle[ppe * e + j], velocity_particle[e * ppe + j],
852 position_particle[ppe * e + j], newvel, local_pos, sol, -dt);
853
854 // RK3:
855 // RowVectorNd bypass, vel2, vel3;
856 // trace_back( gbases, bases, position_particle[ppe * e + i], velocity_particle[e * ppe + i],
857 // bypass, vel2, sol, -0.5 * dt);
858 // trace_back( gbases, bases, position_particle[ppe * e + i], vel2,
859 // bypass, vel3, sol, -0.75 * dt);
860 // trace_back( gbases, bases, position_particle[ppe * e + i],
861 // 2 * velocity_particle[e * ppe + i] + 3 * vel2 + 4 * vel3,
862 // position_particle[ppe * e + i], bypass, sol, -dt / 9);
863
864 // prepare P2G
865 if (cellI_particle[ppe * e + j] >= 0)
866 {
867 // construct interpolator (always linear for P2G, can use gaussian or bspline later)
868 velocity_interpolator[ppe * e + j].compute(cellI_particle[ppe * e + j], dim == 3, local_pos,
869 gbases[cellI_particle[ppe * e + j]], gbases[cellI_particle[ppe * e + j]]);
870 }
871 }
872 }
873#ifdef POLYFEM_WITH_TBB
874 );
875#endif
876
877 // P2G
878 for (int e = 0; e < n_el; ++e)
879 {
880 for (int j = 0; j < ppe; ++j)
881 {
882 int cellI = cellI_particle[ppe * e + j];
883 if (cellI == -1)
884 {
885 continue;
886 }
887 assembler::ElementAssemblyValues &vals = velocity_interpolator[ppe * e + j];
888 for (int i = 0; i < vals.basis_values.size(); ++i)
889 {
890 new_sol.block(bases[cellI].bases[i].global()[0].index * dim, 0, dim, 1) +=
891 vals.basis_values[i].val(0) * velocity_particle[ppe * e + j].transpose();
892 new_sol_w(bases[cellI].bases[i].global()[0].index) += vals.basis_values[i].val(0);
893 }
894 }
895 }
896 // TODO: need to add up boundary velocities and weights because of perodic BC
897
898#ifdef POLYFEM_WITH_TBB
899 tbb::parallel_for(0, (int)new_sol.rows() / dim, 1, [&](int i)
900#else
901 for (int i = 0; i < new_sol.rows() / dim; ++i)
902#endif
903 {
904 sol.block(i * dim, 0, dim, 1) = new_sol.block(i * dim, 0, dim, 1) / new_sol_w(i, 0);
905 }
906#ifdef POLYFEM_WITH_TBB
907 );
908#endif
909 // TODO: need to think about what to do with negative quadratic weight
910 }
911
912 void OperatorSplittingSolver::solve_diffusion_1st(const StiffnessMatrix &mass, const std::vector<int> &bnd_nodes, Eigen::MatrixXd &sol)
913 {
914 Eigen::VectorXd rhs;
915
916 for (int d = 0; d < dim; d++)
917 {
918 Eigen::VectorXd x(sol.size() / dim);
919 for (int j = 0; j < x.size(); j++)
920 {
921 x(j) = sol(j * dim + d);
922 }
923 rhs = mass * x;
924
925 // keep dirichlet bc
926 for (int i = 0; i < bnd_nodes.size(); i++)
927 {
928 rhs(bnd_nodes[i]) = x(bnd_nodes[i]);
929 }
930
931 // if (solver_type == "Pardiso" || solver_type == "Eigen::SimplicialLDLT" || solver_type == "Eigen::SparseLU")
932 {
933 dirichlet_solve_prefactorized(*solver_diffusion, mat_diffusion, rhs, bnd_nodes, x);
934 }
935 // else
936 // {
937 // dirichlet_solve(*solver_diffusion, mat_diffusion, rhs, bnd_nodes, x, mat_diffusion.rows(), "", false, false, false);
938 // }
939
940 for (int j = 0; j < x.size(); j++)
941 {
942 sol(j * dim + d) = x(j);
943 }
944 }
945 }
946
947 void OperatorSplittingSolver::external_force(const mesh::Mesh &mesh,
948 const assembler::Assembler &assembler,
949 const std::vector<basis::ElementBases> &gbases,
950 const std::vector<basis::ElementBases> &bases,
951 const double dt,
952 Eigen::MatrixXd &sol,
953 const Eigen::MatrixXd &local_pts,
954 const std::shared_ptr<assembler::Problem> problem,
955 const double time)
956 {
957#ifdef POLYFEM_WITH_TBB
958 tbb::parallel_for(0, n_el, 1, [&](int e)
959#else
960 for (int e = 0; e < n_el; e++)
961#endif
962 {
963 Eigen::MatrixXd mapped;
964 gbases[e].eval_geom_mapping(local_pts, mapped);
965
966 for (int local_idx = 0; local_idx < bases[e].bases.size(); local_idx++)
967 {
968 int global_idx = bases[e].bases[local_idx].global()[0].index;
969
970 Eigen::MatrixXd pos = Eigen::MatrixXd::Zero(1, dim);
971 for (int d = 0; d < dim; d++)
972 pos(0, d) = mapped(local_idx, d);
973
974 Eigen::MatrixXd val;
975 problem->rhs(assembler, pos, time, val);
976
977 for (int d = 0; d < dim; d++)
978 {
979 sol(global_idx * dim + d) += val(d) * dt;
980 }
981 }
982 }
983#ifdef POLYFEM_WITH_TBB
984 );
985#endif
986 }
987
988 void OperatorSplittingSolver::solve_pressure(const StiffnessMatrix &mixed_stiffness, const std::vector<int> &pressure_boundary_nodes, Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
989 {
990 Eigen::VectorXd rhs;
991 if (pressure_boundary_nodes.size() == 0)
992 rhs = Eigen::VectorXd::Zero(mixed_stiffness.rows() + 1); // mixed_stiffness * sol;
993 else
994 rhs = Eigen::VectorXd::Zero(mixed_stiffness.rows()); // mixed_stiffness * sol;
995
996 Eigen::VectorXd temp = mixed_stiffness * sol;
997 for (int i = 0; i < temp.rows(); i++)
998 {
999 rhs(i) = temp(i);
1000 }
1001
1002 Eigen::VectorXd x = Eigen::VectorXd::Zero(rhs.size());
1003 for (int i = 0; i < pressure.size(); i++)
1004 {
1005 x(i) = pressure(i);
1006 }
1007
1008 for (int i = 0; i < pressure_boundary_nodes.size(); i++)
1009 {
1010 rhs(pressure_boundary_nodes[i]) = 0;
1011 x(pressure_boundary_nodes[i]) = 0;
1012 }
1013 // if (solver_type == "Pardiso" || solver_type == "Eigen::SimplicialLDLT" || solver_type == "Eigen::SparseLU")
1014 // {
1015 dirichlet_solve_prefactorized(*solver_projection, mat_projection, rhs, pressure_boundary_nodes, x);
1016 // }
1017 // else
1018 // {
1019 // dirichlet_solve(*solver_projection, mat_projection, rhs, std::vector<int>(), x, mat_projection.rows() - 1, "", false, false, false);
1020 // }
1021
1022 if (pressure_boundary_nodes.size() == 0)
1023 pressure = x.head(x.size() - 1);
1024 else
1025 pressure = x;
1026 }
1027
1028 void OperatorSplittingSolver::projection(const StiffnessMatrix &velocity_mass, const StiffnessMatrix &mixed_stiffness, const std::vector<int> &boundary_nodes_, Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure)
1029 {
1030 Eigen::VectorXd rhs = mixed_stiffness.transpose() * pressure;
1031 Eigen::VectorXd dx = Eigen::VectorXd::Zero(sol.size());
1032
1033 for (int i = 0; i < boundary_nodes_.size(); i++)
1034 {
1035 rhs(boundary_nodes_[i]) = 0;
1036 }
1037
1038 if (solver_diffusion->name() == "Pardiso" || solver_diffusion->name() == "Eigen::SimplicialLDLT" || solver_diffusion->name() == "Eigen::SparseLU")
1039 {
1040 dirichlet_solve_prefactorized(*solver_mass, velocity_mass, rhs, boundary_nodes_, dx);
1041 }
1042 else
1043 {
1044 auto mat = velocity_mass;
1045 dirichlet_solve(*solver_mass, mat, rhs, boundary_nodes_, dx, velocity_mass.rows(), "", false, false, false);
1046 }
1047
1048 sol -= dx;
1049 }
1050
1051 void OperatorSplittingSolver::projection(int n_bases,
1052 const std::vector<basis::ElementBases> &gbases,
1053 const std::vector<basis::ElementBases> &bases,
1054 const std::vector<basis::ElementBases> &pressure_bases,
1055 const Eigen::MatrixXd &local_pts,
1056 Eigen::MatrixXd &pressure,
1057 Eigen::MatrixXd &sol)
1058 {
1059 Eigen::VectorXd grad_pressure = Eigen::VectorXd::Zero(n_bases * dim);
1060 Eigen::VectorXi traversed = Eigen::VectorXi::Zero(n_bases);
1061
1062 assembler::ElementAssemblyValues vals;
1063 for (int e = 0; e < n_el; ++e)
1064 {
1065 vals.compute(e, dim == 3, local_pts, pressure_bases[e], gbases[e]);
1066 for (int j = 0; j < local_pts.rows(); j++)
1067 {
1068 int global_ = bases[e].bases[j].global()[0].index;
1069 for (int i = 0; i < vals.basis_values.size(); i++)
1070 {
1071 for (int d = 0; d < dim; d++)
1072 {
1073 assert(pressure_bases[e].bases[i].global().size() == 1);
1074 grad_pressure(global_ * dim + d) += vals.basis_values[i].grad_t_m(j, d) * pressure(pressure_bases[e].bases[i].global()[0].index);
1075 }
1076 }
1077 traversed(global_)++;
1078 }
1079 }
1080 for (int i = 0; i < traversed.size(); i++)
1081 {
1082 for (int d = 0; d < dim; d++)
1083 {
1084 sol(i * dim + d) -= grad_pressure(i * dim + d) / traversed(i);
1085 }
1086 }
1087 }
1088
1089 void OperatorSplittingSolver::initialize_density(const std::shared_ptr<assembler::Problem> &problem)
1090 {
1091 Eigen::MatrixXd pts(1, dim);
1092 Eigen::MatrixXd tmp;
1093 for (long i = 0; i <= grid_cell_num(0); i++)
1094 {
1095 pts(0, 0) = i * resolution + min_domain(0);
1096 for (long j = 0; j <= grid_cell_num(1); j++)
1097 {
1098 pts(0, 1) = j * resolution + min_domain(1);
1099 if (dim == 2)
1100 {
1101 const long idx = i + j * (grid_cell_num(0) + 1);
1102 problem->initial_density(pts, tmp);
1103 density(idx) = tmp(0);
1104 }
1105 else
1106 {
1107 for (long k = 0; k <= grid_cell_num(2); k++)
1108 {
1109 pts(0, 2) = k * resolution + min_domain(2);
1110 const long idx = i + (j + k * (grid_cell_num(1) + 1)) * (grid_cell_num(0) + 1);
1111 problem->initial_density(pts, tmp);
1112 density(idx) = tmp(0);
1113 }
1114 }
1115 }
1116 }
1117 }
1118
1119 long OperatorSplittingSolver::search_cell(const std::vector<basis::ElementBases> &gbases, const RowVectorNd &pos, Eigen::MatrixXd &local_pts)
1120 {
1121 Eigen::Matrix<long, Eigen::Dynamic, 1> pos_int(dim);
1122 for (int d = 0; d < dim; d++)
1123 {
1124 pos_int(d) = floor((pos(d) - min_domain(d)) / (max_domain(d) - min_domain(d)) * hash_table_cell_num(d));
1125 if (pos_int(d) < 0)
1126 return -1;
1127 else if (pos_int(d) >= hash_table_cell_num(d))
1128 return -1;
1129 }
1130
1131 long idx = 0, dim_num = 1;
1132 for (int d = 0; d < dim; d++)
1133 {
1134 idx += pos_int(d) * dim_num;
1135 dim_num *= hash_table_cell_num(d);
1136 }
1137
1138 const std::vector<long> &list = hash_table[idx];
1139 for (auto it = list.begin(); it != list.end(); it++)
1140 {
1141 calculate_local_pts(gbases[*it], *it, pos, local_pts);
1142
1143 if (shape == dim + 1)
1144 {
1145 if (local_pts.minCoeff() > -1e-13 && local_pts.sum() < 1 + 1e-13)
1146 return *it;
1147 }
1148 else
1149 {
1150 if (local_pts.minCoeff() > -1e-13 && local_pts.maxCoeff() < 1 + 1e-13)
1151 return *it;
1152 }
1153 }
1154 return -1; // not inside any elem
1155 }
1156
1157 bool OperatorSplittingSolver::outside_quad(const std::vector<RowVectorNd> &vert, const RowVectorNd &pos)
1158 {
1159 double a = (vert[1](0) - vert[0](0)) * (pos(1) - vert[0](1)) - (vert[1](1) - vert[0](1)) * (pos(0) - vert[0](0));
1160 double b = (vert[2](0) - vert[1](0)) * (pos(1) - vert[1](1)) - (vert[2](1) - vert[1](1)) * (pos(0) - vert[1](0));
1161 double c = (vert[3](0) - vert[2](0)) * (pos(1) - vert[2](1)) - (vert[3](1) - vert[2](1)) * (pos(0) - vert[2](0));
1162 double d = (vert[0](0) - vert[3](0)) * (pos(1) - vert[3](1)) - (vert[0](1) - vert[3](1)) * (pos(0) - vert[3](0));
1163
1164 if ((a > 0 && b > 0 && c > 0 && d > 0) || (a < 0 && b < 0 && c < 0 && d < 0))
1165 return false;
1166 return true;
1167 }
1168
1169 void OperatorSplittingSolver::compute_gbase_val(const int elem_idx, const Eigen::MatrixXd &local_pos, Eigen::MatrixXd &pos)
1170 {
1171 pos = Eigen::MatrixXd::Zero(1, dim);
1172 Eigen::VectorXd weights(shape);
1173 const double x = local_pos(0);
1174 const double y = local_pos(1);
1175 if (dim == 2)
1176 {
1177 if (shape == 3)
1178 {
1179 weights << 1 - x - y, x, y;
1180 }
1181 else
1182 {
1183 weights << (1 - x) * (1 - y), x * (1 - y),
1184 x * y, y * (1 - x);
1185 }
1186 }
1187 else
1188 {
1189 const double z = local_pos(2);
1190 if (shape == 4)
1191 {
1192 weights << 1 - x - y - z, x, y, z;
1193 }
1194 else
1195 {
1196 weights << (1 - x) * (1 - y) * (1 - z),
1197 x * (1 - y) * (1 - z),
1198 x * y * (1 - z),
1199 (1 - x) * y * (1 - z),
1200 (1 - x) * (1 - y) * z,
1201 x * (1 - y) * z,
1202 x * y * z,
1203 (1 - x) * y * z;
1204 }
1205 }
1206
1207 for (int d = 0; d < dim; d++)
1208 {
1209 for (int j = 0; j < shape; j++)
1210 {
1211 pos(d) += V(T(elem_idx, j), d) * weights(j);
1212 }
1213 }
1214 }
1215
1216 void OperatorSplittingSolver::compute_gbase_jacobi(const int elem_idx, const Eigen::MatrixXd &local_pos, Eigen::MatrixXd &jacobi)
1217 {
1218 jacobi = Eigen::MatrixXd::Zero(dim, dim);
1219 Eigen::MatrixXd weights(shape, dim);
1220 const double x = local_pos(0);
1221 const double y = local_pos(1);
1222 if (dim == 2)
1223 {
1224 if (shape == 3)
1225 {
1226 weights << -1, -1,
1227 1, 0,
1228 0, 1;
1229 }
1230 else
1231 {
1232 weights << y - 1, x - 1,
1233 1 - y, -x,
1234 y, x,
1235 -y, 1 - x;
1236 }
1237 }
1238 else
1239 {
1240 const double z = local_pos(2);
1241 if (shape == 4)
1242 {
1243 weights << -1, -1, -1,
1244 1, 0, 0,
1245 0, 1, 0,
1246 0, 0, 1;
1247 }
1248 else
1249 {
1250 weights << -(1 - y) * (1 - z), -(1 - x) * (1 - z), -(1 - x) * (1 - y),
1251 (1 - y) * (1 - z), -x * (1 - z), -x * (1 - y),
1252 y * (1 - z), x * (1 - z), -x * y,
1253 -y * (1 - z), (1 - x) * (1 - z), -(1 - x) * y,
1254 -(1 - y) * z, -(1 - x) * z, (1 - x) * (1 - y),
1255 (1 - y) * z, -x * z, x * (1 - y),
1256 y * z, x * z, x * y,
1257 -y * z, (1 - x) * z, (1 - x) * y;
1258 }
1259 }
1260
1261 for (int d1 = 0; d1 < dim; d1++)
1262 {
1263 for (int d2 = 0; d2 < dim; d2++)
1264 {
1265 for (int i = 0; i < shape; i++)
1266 {
1267 jacobi(d1, d2) += weights(i, d2) * V(T(elem_idx, i), d1);
1268 }
1269 }
1270 }
1271 }
1272
1273 void OperatorSplittingSolver::calculate_local_pts(const basis::ElementBases &gbase,
1274 const int elem_idx,
1275 const RowVectorNd &pos,
1276 Eigen::MatrixXd &local_pos)
1277 {
1278 local_pos = Eigen::MatrixXd::Zero(1, dim);
1279
1280 // if(shape == 4 && dim == 2 && outside_quad(vert, pos))
1281 // {
1282 // local_pos(0) = local_pos(1) = -1;
1283 // return;
1284 // }
1285 Eigen::MatrixXd res;
1286 int iter_times = 0;
1287 int max_iter = 20;
1288 do
1289 {
1290 res = -pos;
1291
1292 Eigen::MatrixXd mapped;
1293 gbase.eval_geom_mapping(local_pos, mapped);
1294 for (int d = 0; d < dim; d++)
1295 res(d) += mapped(0, d);
1296
1297 std::vector<Eigen::MatrixXd> grads;
1298 gbase.eval_geom_mapping_grads(local_pos, grads);
1299 Eigen::MatrixXd jacobi = grads[0].transpose();
1300
1301 Eigen::VectorXd delta = jacobi.colPivHouseholderQr().solve(res.transpose());
1302 for (int d = 0; d < dim; d++)
1303 {
1304 local_pos(d) -= delta(d);
1305 }
1306 iter_times++;
1307 if (shape == dim + 1)
1308 break;
1309 } while (res.norm() > 1e-12 && iter_times < max_iter);
1310
1311 if (iter_times >= max_iter)
1312 {
1313 for (int d = 0; d < dim; d++)
1314 local_pos(d) = -1;
1315 }
1316 }
1317
1318 void OperatorSplittingSolver::save_density()
1319 {
1320#ifdef POLYFEM_WITH_OPENVDB
1321 openvdb::initialize();
1322 openvdb::FloatGrid::Ptr grid = openvdb::FloatGrid::create();
1323 openvdb::FloatGrid::Accessor accessor = grid->getAccessor();
1324
1325 for (int i = 0; i <= grid_cell_num(0); i++)
1326 {
1327 for (int j = 0; j <= grid_cell_num(1); j++)
1328 {
1329 if (dim == 2)
1330 {
1331 const int idx = i + j * (grid_cell_num(0) + 1);
1332 openvdb::Coord xyz(i, j, 0);
1333 if (density(idx) > 1e-8)
1334 accessor.setValue(xyz, density(idx));
1335 }
1336 else
1337 {
1338 for (int k = 0; k <= grid_cell_num(2); k++)
1339 {
1340 const int idx = i + (j + k * (grid_cell_num(1) + 1)) * (grid_cell_num(0) + 1);
1341 openvdb::Coord xyz(i, j, k);
1342 if (density(idx) > 1e-8)
1343 accessor.setValue(xyz, density(idx));
1344 }
1345 }
1346 }
1347 }
1348 grid->setName("density_smoke");
1349 grid->setGridClass(openvdb::GRID_FOG_VOLUME);
1350
1351 static int num_frame = 0;
1352 const std::string filename = "density" + std::to_string(num_frame) + ".vdb";
1353 openvdb::io::File file(filename.c_str());
1354 num_frame++;
1355
1356 openvdb::GridPtrVec(grids);
1357 grids.push_back(grid);
1358 file.write(grids);
1359 file.close();
1360#else
1361 static int num_frame = 0;
1362 std::string name = "density" + std::to_string(num_frame) + ".txt";
1363 std::ofstream file(name.c_str());
1364 num_frame++;
1365 for (int i = 0; i <= grid_cell_num(0); i++)
1366 {
1367 for (int j = 0; j <= grid_cell_num(1); j++)
1368 {
1369 if (dim == 2)
1370 {
1371 const int idx = i + j * (grid_cell_num(0) + 1);
1372 if (density(idx) < 1e-10)
1373 continue;
1374 file << i << " " << j << " " << density(idx) << std::endl;
1375 }
1376 else
1377 {
1378 for (int k = 0; k <= grid_cell_num(2); k++)
1379 {
1380 const int idx = i + (j + k * (grid_cell_num(1) + 1)) * (grid_cell_num(0) + 1);
1381 if (density(idx) < 1e-10)
1382 continue;
1383 file << i << " " << j << " " << k << " " << density(idx) << std::endl;
1384 }
1385 }
1386 }
1387 }
1388 file.close();
1389#endif
1390 };
1391} // namespace polyfem::solver
val
double val
Definition Assembler.cpp:86
vals
ElementAssemblyValues vals
Definition Assembler.cpp:22
gvals
assembler::ElementAssemblyValues gvals
Definition BodyForm.cpp:25
y
int y
Definition SplineBasis3d.cpp:55
z
int z
Definition SplineBasis3d.cpp:55
x
int x
Definition SplineBasis3d.cpp:55
mat
Eigen::MatrixXd mat
Definition TransientForm.cpp:24
polyfem::assembler::ElementAssemblyValues
stores per element basis values at given quadrature points and geometric mapping
Definition ElementAssemblyValues.hpp:14
polyfem::assembler::ElementAssemblyValues::compute
void compute(const int el_index, const bool is_volume, const Eigen::MatrixXd &pts, const basis::ElementBases &basis, const basis::ElementBases &gbasis)
computes the per element values at the local (ref el) points (pts) sets basis_values,...
Definition ElementAssemblyValues.cpp:151
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::eval_geom_mapping_grads
void eval_geom_mapping_grads(const Eigen::MatrixXd &samples, std::vector< Eigen::MatrixXd > &grads) const
Evaluate the gradients of the geometric mapping defined above.
Definition ElementBases.cpp:82
polyfem::mesh::Mesh
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:39
polyfem::mesh::Mesh::n_vertices
virtual int n_vertices() const =0
number of vertices
polyfem::mesh::Mesh::get_edges
virtual void get_edges(Eigen::MatrixXd &p0, Eigen::MatrixXd &p1) const =0
Get all the edges.
polyfem::mesh::Mesh::bounding_box
virtual void bounding_box(RowVectorNd &min, RowVectorNd &max) const =0
computes the bbox of the mesh
polyfem::mesh::Mesh::point
virtual RowVectorNd point(const int global_index) const =0
point coordinates
polyfem::mesh::Mesh::dimension
int dimension() const
utily for dimension
Definition Mesh.hpp:148
polyfem::mesh::Mesh::cell_vertex
virtual int cell_vertex(const int f_id, const int lv_id) const =0
id of the vertex of a cell
polyfem::solver::OperatorSplittingSolver::initialize_mesh
void initialize_mesh(const mesh::Mesh &mesh, const int shape, const int n_el, const std::vector< mesh::LocalBoundary > &local_boundary)
Definition OperatorSplittingSolver.cpp:32
polyfem::solver::OperatorSplittingSolver::solve_pressure
void solve_pressure(const StiffnessMatrix &mixed_stiffness, const std::vector< int > &pressure_boundary_nodes, Eigen::MatrixXd &sol, Eigen::MatrixXd &pressure)
Definition OperatorSplittingSolver.cpp:988
polyfem::solver::OperatorSplittingSolver::projection
void projection(const StiffnessMatrix &velocity_mass, const StiffnessMatrix &mixed_stiffness, const std::vector< int > &boundary_nodes_, Eigen::MatrixXd &sol, const Eigen::MatrixXd &pressure)
Definition OperatorSplittingSolver.cpp:1028
polyfem::solver::OperatorSplittingSolver::interpolator
int interpolator(const std::vector< basis::ElementBases > &gbases, const std::vector< basis::ElementBases > &bases, const RowVectorNd &pos, RowVectorNd &vel, Eigen::MatrixXd &local_pos, const Eigen::MatrixXd &sol)
Definition OperatorSplittingSolver.cpp:285
polyfem::solver::OperatorSplittingSolver::compute_gbase_val
void compute_gbase_val(const int elem_idx, const Eigen::MatrixXd &local_pos, Eigen::MatrixXd &pos)
Definition OperatorSplittingSolver.cpp:1169
polyfem::solver::OperatorSplittingSolver::external_force
void external_force(const mesh::Mesh &mesh, const assembler::Assembler &assembler, const std::vector< basis::ElementBases > &gbases, const std::vector< basis::ElementBases > &bases, const double dt, Eigen::MatrixXd &sol, const Eigen::MatrixXd &local_pts, const std::shared_ptr< assembler::Problem > problem, const double time)
Definition OperatorSplittingSolver.cpp:947
polyfem::solver::OperatorSplittingSolver::solver_diffusion
std::unique_ptr< polysolve::linear::Solver > solver_diffusion
Definition OperatorSplittingSolver.hpp:172
polyfem::solver::OperatorSplittingSolver::advection
void advection(const mesh::Mesh &mesh, const std::vector< basis::ElementBases > &gbases, const std::vector< basis::ElementBases > &bases, Eigen::MatrixXd &sol, const double dt, const Eigen::MatrixXd &local_pts, const int order=1, const int RK=1)
Definition OperatorSplittingSolver.cpp:354
polyfem::solver::OperatorSplittingSolver::initialize_density
void initialize_density(const std::shared_ptr< assembler::Problem > &problem)
Definition OperatorSplittingSolver.cpp:1089
polyfem::solver::OperatorSplittingSolver::advection_PIC
void advection_PIC(const mesh::Mesh &mesh, const std::vector< basis::ElementBases > &gbases, const std::vector< basis::ElementBases > &bases, Eigen::MatrixXd &sol, const double dt, const Eigen::MatrixXd &local_pts, const int order=1)
Definition OperatorSplittingSolver.cpp:789
polyfem::solver::OperatorSplittingSolver::position_particle
std::vector< Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > > position_particle
Definition OperatorSplittingSolver.hpp:163
polyfem::solver::OperatorSplittingSolver::search_cell
long search_cell(const std::vector< basis::ElementBases > &gbases, const RowVectorNd &pos, Eigen::MatrixXd &local_pts)
Definition OperatorSplittingSolver.cpp:1119
polyfem::solver::OperatorSplittingSolver::hash_table
std::vector< std::vector< long > > hash_table
Definition OperatorSplittingSolver.hpp:160
polyfem::solver::OperatorSplittingSolver::trace_back
int trace_back(const std::vector< basis::ElementBases > &gbases, const std::vector< basis::ElementBases > &bases, const RowVectorNd &pos_1, const RowVectorNd &vel_1, RowVectorNd &pos_2, RowVectorNd &vel_2, Eigen::MatrixXd &local_pos, const Eigen::MatrixXd &sol, const double dt)
Definition OperatorSplittingSolver.cpp:270
polyfem::solver::OperatorSplittingSolver::advect_density
void advect_density(const std::vector< basis::ElementBases > &gbases, const std::vector< basis::ElementBases > &bases, const Eigen::MatrixXd &sol, const double dt, const int RK=3)
Definition OperatorSplittingSolver.cpp:481
polyfem::solver::OperatorSplittingSolver::advection_FLIP
void advection_FLIP(const mesh::Mesh &mesh, const std::vector< basis::ElementBases > &gbases, const std::vector< basis::ElementBases > &bases, Eigen::MatrixXd &sol, const double dt, const Eigen::MatrixXd &local_pts, const int order=1)
Definition OperatorSplittingSolver.cpp:554
polyfem::solver::OperatorSplittingSolver::outside_quad
bool outside_quad(const std::vector< RowVectorNd > &vert, const RowVectorNd &pos)
Definition OperatorSplittingSolver.cpp:1157
polyfem::solver::OperatorSplittingSolver::grid_cell_num
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > grid_cell_num
Definition OperatorSplittingSolver.hpp:182
polyfem::solver::OperatorSplittingSolver::hash_table_cell_num
Eigen::Matrix< long, Eigen::Dynamic, 1, Eigen::ColMajor, 3, 1 > hash_table_cell_num
Definition OperatorSplittingSolver.hpp:161
polyfem::solver::OperatorSplittingSolver::solver_projection
std::unique_ptr< polysolve::linear::Solver > solver_projection
Definition OperatorSplittingSolver.hpp:173
polyfem::solver::OperatorSplittingSolver::velocity_particle
std::vector< Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > > velocity_particle
Definition OperatorSplittingSolver.hpp:164
polyfem::solver::OperatorSplittingSolver::calculate_local_pts
void calculate_local_pts(const basis::ElementBases &gbase, const int elem_idx, const RowVectorNd &pos, Eigen::MatrixXd &local_pos)
Definition OperatorSplittingSolver.cpp:1273
polyfem::solver::OperatorSplittingSolver::solver_mass
std::unique_ptr< polysolve::linear::Solver > solver_mass
Definition OperatorSplittingSolver.hpp:174
polyfem::solver::OperatorSplittingSolver::compute_gbase_jacobi
void compute_gbase_jacobi(const int elem_idx, const Eigen::MatrixXd &local_pos, Eigen::MatrixXd &jacobi)
Definition OperatorSplittingSolver.cpp:1216
polyfem::solver::OperatorSplittingSolver::advect_density_exact
void advect_density_exact(const std::vector< basis::ElementBases > &gbases, const std::vector< basis::ElementBases > &bases, const std::shared_ptr< assembler::Problem > problem, const double t, const double dt, const int RK=3)
Definition OperatorSplittingSolver.cpp:408
polyfem::solver::OperatorSplittingSolver::initialize_solver
void initialize_solver(const mesh::Mesh &mesh, const int shape_, const int n_el_, const std::vector< mesh::LocalBoundary > &local_boundary, const std::vector< int > &bnd_nodes)
Definition OperatorSplittingSolver.cpp:140
polyfem::solver::OperatorSplittingSolver::solve_diffusion_1st
void solve_diffusion_1st(const StiffnessMatrix &mass, const std::vector< int > &bnd_nodes, Eigen::MatrixXd &sol)
Definition OperatorSplittingSolver.cpp:912
polyfem::solver::OperatorSplittingSolver::initialize_grid
void initialize_grid(const mesh::Mesh &mesh, const std::vector< basis::ElementBases > &gbases, const std::vector< basis::ElementBases > &bases, const double &density_dx)
Definition OperatorSplittingSolver.cpp:14
polyfem::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:42
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:11
polyfem::StiffnessMatrix
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:20