Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PolygonalBasis3d.cpp
Go to the documentation of this file.
1
3#include "LagrangeBasis3d.hpp"
4
15
17
18#include <igl/per_vertex_normals.h>
19#include <random>
20#include <memory>
22
23namespace polyfem
24{
25 using namespace assembler;
26 using namespace mesh;
27 using namespace quadrature;
28 using namespace utils;
29
30 namespace basis
31 {
32 namespace
33 {
34
35 const int max_num_kernels = 300;
36
37 // -----------------------------------------------------------------------------
38
39 std::vector<int> compute_nonzero_bases_ids(const Mesh3D &mesh, const int c,
40 const std::vector<ElementBases> &bases,
41 const std::map<int, InterfaceData> &poly_face_to_data)
42 {
43 std::vector<int> local_to_global;
44
45 for (int lf = 0; lf < mesh.n_cell_faces(c); ++lf)
46 {
47 auto index = mesh.get_index_from_element(c, lf, 0);
48 const int c2 = mesh.switch_element(index).element;
49 assert(c2 >= 0); // no boundary polytope
50 assert(poly_face_to_data.count(index.face) > 0);
51 const InterfaceData &bdata = poly_face_to_data.at(index.face);
52 const ElementBases &b = bases[c2];
53 for (int other_local_basis_id : bdata.local_indices)
54 {
55 for (const auto &x : b.bases[other_local_basis_id].global())
56 {
57 const int global_node_id = x.index;
58 local_to_global.push_back(global_node_id);
59 }
60 }
61 }
62
63 std::sort(local_to_global.begin(), local_to_global.end());
64 auto it = std::unique(local_to_global.begin(), local_to_global.end());
65 local_to_global.resize(std::distance(local_to_global.begin(), it));
66
67 return local_to_global;
68 }
69
70 // -----------------------------------------------------------------------------
71
72 // Canonical triangle mesh in parametric domain
73 void compute_canonical_pattern(int n_samples_per_edge, Eigen::MatrixXd &V, Eigen::MatrixXi &F)
74 {
75 regular_2d_grid(n_samples_per_edge, false, V, F);
76
77 // igl::opengl::glfw::Viewer viewer;
78 // viewer.data().set_mesh(V, F);
79 // viewer.launch();
80 }
81
82 // -----------------------------------------------------------------------------
83
84 // Needs to be consistent between `evalFunc` and `compute_quad_mesh_from_cell`
85 constexpr int lv0 = 3;
86
87 // Assemble the surface quad mesh (V, F) corresponding to the polyhedron c
88 GetAdjacentLocalEdge compute_quad_mesh_from_cell(
89 const Mesh3D &mesh, int c, Eigen::MatrixXd &V, Eigen::MatrixXi &F)
90 {
91 std::vector<std::array<int, 4>> quads(mesh.n_cell_faces(c));
92 typedef std::tuple<int, int, bool> QuadLocalEdge;
93 std::vector<std::array<QuadLocalEdge, 4>> adj(mesh.n_cell_faces(c));
94 int num_vertices = 0;
95 std::map<int, int> vertex_g2l;
96 std::map<int, int> face_g2l;
97 for (int lf = 0; lf < mesh.n_cell_faces(c); ++lf)
98 {
99 face_g2l.emplace(mesh.get_index_from_element(c, lf, lv0).face, lf);
100 }
101 for (int lf = 0; lf < mesh.n_cell_faces(c); ++lf)
102 {
103 auto index = mesh.get_index_from_element(c, lf, lv0);
104 assert(mesh.n_face_vertices(index.face) == 4);
105 for (int lv = 0; lv < 4; ++lv)
106 {
107 if (!vertex_g2l.count(index.vertex))
108 {
109 vertex_g2l.emplace(index.vertex, num_vertices++);
110 }
111 quads[lf][lv] = vertex_g2l.at(index.vertex);
112
113 // Set adjacency info
114 auto index2 = mesh.switch_face(index);
115 int lf2 = face_g2l.at(index2.face);
116 std::get<0>(adj[lf][lv]) = lf2;
117 auto index3 = mesh.get_index_from_element(c, lf2, lv0);
118 for (int lv2 = 0; lv2 < 4; ++lv2)
119 {
120 if (index3.edge == index2.edge)
121 {
122 std::get<1>(adj[lf][lv]) = lv2;
123 if (index2.vertex != index3.vertex)
124 {
125 assert(mesh.switch_vertex(index3).vertex == index2.vertex);
126 std::get<2>(adj[lf][lv]) = true;
127 }
128 else
129 {
130 std::get<2>(adj[lf][lv]) = false;
131 }
132 }
133 index3 = mesh.next_around_face(index3);
134 }
135
136 index = mesh.next_around_face(index);
137 }
138 }
139 V.resize(num_vertices, 3);
140 for (const auto &kv : vertex_g2l)
141 {
142 V.row(kv.second) = mesh.point(kv.first);
143 }
144 F.resize(quads.size(), 4);
145 int f = 0;
146 for (auto q : quads)
147 {
148 F.row(f++) << q[0], q[1], q[2], q[3];
149 }
150
151 return [adj](int q, int lv) {
152 return adj[q][lv];
153 };
154 }
155
156 // -----------------------------------------------------------------------------
157
158 void compute_offset_kernels(const Eigen::MatrixXd &QV, const Eigen::MatrixXi &QF,
159 int n_kernels_per_edge, double eps, Eigen::MatrixXd &kernel_centers,
160 Eigen::MatrixXd &KV, Eigen::MatrixXi &KF,
161 EvalParametersFunc evalFuncGeom, GetAdjacentLocalEdge getAdjLocalEdge)
162 {
163 Eigen::MatrixXd PV, KN;
164 Eigen::MatrixXi PF;
165 Eigen::VectorXd D;
166 compute_canonical_pattern(n_kernels_per_edge, PV, PF);
167 instantiate_pattern(QV, QF, PV, PF, KV, KF, nullptr, evalFuncGeom, getAdjLocalEdge);
168 orient_closed_surface(KV, KF);
169 double volume = std::pow(signed_volume(KV, KF), 1.0 / 3.0);
170
171 if (true || KV.rows() < max_num_kernels)
172 {
173 igl::per_vertex_normals(KV, KF, KN);
174 kernel_centers = KV;
175 }
176 else
177 {
178 polyfem::sample_surface(KV, KF, max_num_kernels, kernel_centers, &KN, 10, 10);
179 }
180 kernel_centers += eps * volume * KN;
181
182 // std::default_random_engine gen;
183 // std::uniform_real_distribution<double> dist(-1.0, 1.0);
184 // for (int v = 0; v < kernel_centers.rows(); ++v) {
185 // kernel_centers.row(v) = KV.row(v) + dist(gen) * KN.row(v);
186 // }
187 assert(kernel_centers.cols() == 3);
188 signed_squared_distances(KV, KF, kernel_centers, D);
189 std::vector<Eigen::RowVector3d> remap;
190 std::vector<Eigen::RowVector3d> rejected;
191 for (int v = 0; v < kernel_centers.rows(); ++v)
192 {
193 if (std::sqrt(D(v)) > 0.8 * eps * volume)
194 {
195 remap.push_back(kernel_centers.row(v));
196 }
197 else
198 {
199 // rejected.push_back(kernel_centers.row(v));
200 }
201 }
202 kernel_centers.resize(remap.size(), 3);
203 for (int v = 0; v < kernel_centers.rows(); ++v)
204 {
205 kernel_centers.row(v) = remap[v];
206 }
207 }
208
209 // -----------------------------------------------------------------------------
210
212 void sample_polyhedra(
213 const int element_index,
214 const int n_quadrature_vertices_per_edge,
215 const int n_kernels_per_edge,
216 int n_samples_per_edge,
217 const int quadrature_order,
218 const int mass_quadrature_order,
219 const Mesh3D &mesh,
220 const std::map<int, InterfaceData> &poly_face_to_data,
221 const std::vector<ElementBases> &bases,
222 const std::vector<ElementBases> &gbases,
223 const double eps,
224 std::vector<int> &local_to_global,
225 Eigen::MatrixXd &collocation_points,
226 Eigen::MatrixXd &kernel_centers,
227 Eigen::MatrixXd &rhs,
228 Eigen::MatrixXd &triangulated_vertices,
229 Eigen::MatrixXi &triangulated_faces,
230 Quadrature &quadrature,
231 Quadrature &mass_quadrature,
232 double &scaling,
233 Eigen::RowVector3d &translation)
234 {
235 // Local ids of nonzero bases over the polygon
236 local_to_global = compute_nonzero_bases_ids(mesh, element_index, bases, poly_face_to_data);
237
238 // Compute the image of the canonical pattern vertices through the geometric mapping
239 // of the given local face
240 auto evalFunc = [&](const Eigen::MatrixXd &uv, Eigen::MatrixXd &mapped, int lf) {
241 const auto &u = uv.col(0).array();
242 const auto &v = uv.col(1).array();
243 auto index = mesh.get_index_from_element(element_index, lf, lv0);
244 index = mesh.switch_element(index);
245 // Eigen::MatrixXd abcd = LagrangeBasis3d::linear_hex_face_local_nodes_coordinates(mesh, index);
246 const auto indices = LagrangeBasis3d::hex_face_local_nodes(false, 1, mesh, index);
247 assert(indices.size() == 4);
248 Eigen::MatrixXd abcd;
250 Eigen::RowVector3d a = abcd.row(indices(0));
251 Eigen::RowVector3d b = abcd.row(indices(1));
252 Eigen::RowVector3d c = abcd.row(indices(2));
253 Eigen::RowVector3d d = abcd.row(indices(3));
254 mapped = ((1 - u) * (1 - v)).matrix() * a + (u * (1 - v)).matrix() * b + (u * v).matrix() * c + ((1 - u) * v).matrix() * d;
255 mapped = mapped.array().max(0.0).min(1.0);
256 assert(mapped.maxCoeff() >= 0.0);
257 assert(mapped.maxCoeff() <= 1.0);
258 };
259 auto evalFuncGeom = [&](const Eigen::MatrixXd &uv, Eigen::MatrixXd &mapped, int lf) {
260 Eigen::MatrixXd samples;
261 evalFunc(uv, samples, lf);
262 auto index = mesh.get_index_from_element(element_index, lf, lv0);
263 index = mesh.switch_element(index);
264 const ElementBases &gb = gbases[index.element];
265 gb.eval_geom_mapping(samples, mapped);
266 };
267
268 Eigen::MatrixXd QV, KV;
269 Eigen::MatrixXi QF, KF;
270 auto getAdjLocalEdge = compute_quad_mesh_from_cell(mesh, element_index, QV, QF);
271
272 // Compute kernel centers
273 compute_offset_kernels(QV, QF, n_kernels_per_edge, eps, kernel_centers, KV, KF,
274 evalFuncGeom, getAdjLocalEdge);
275 // if (KV.rows() >= max_num_kernels) { n_samples_per_edge = 5; }
276
277 // Compute collocation points
278 Eigen::MatrixXd PV, UV;
279 Eigen::MatrixXi PF, CF, UF;
280 Eigen::VectorXi uv_sources, uv_ranges;
281 compute_canonical_pattern(n_samples_per_edge, PV, PF);
282 instantiate_pattern(QV, QF, PV, PF, UV, UF, &uv_sources, evalFunc, getAdjLocalEdge);
283 orient_closed_surface(UV, UF);
284 instantiate_pattern(QV, QF, PV, PF, collocation_points, CF, nullptr, evalFuncGeom, getAdjLocalEdge);
285 orient_closed_surface(collocation_points, CF);
286 reorder_mesh(collocation_points, CF, uv_sources, uv_ranges);
287 reorder_mesh(UV, UF, uv_sources, uv_ranges);
288 assert(uv_ranges.size() == mesh.n_cell_faces(element_index) + 1);
289
290 // Compute coarse surface surface for visualization
291 compute_canonical_pattern(n_quadrature_vertices_per_edge, PV, PF);
292 instantiate_pattern(QV, QF, PV, PF, triangulated_vertices, triangulated_faces,
293 nullptr, evalFuncGeom, getAdjLocalEdge);
294 orient_closed_surface(triangulated_vertices, triangulated_faces);
295
296 // for (int f = 0; f < KF.rows(); ++f) {
297 // triangulated_faces.row(f) = KF.row(f).reverse();
298 // }
299
300 // {
301 // Eigen::MatrixXd V;
302 // evalFuncGeom(PV, V, 0);
303 // igl::write_triangle_mesh("foo_dense.obj", collocation_points, CF);
304 // igl::write_triangle_mesh("foo_small.obj", triangulated_vertices, triangulated_faces);
305 // igl::opengl::glfw::Viewer viewer;
306 // viewer.data().set_points(kernel_centers, Eigen::RowVector3d(1,0,1));
307 // viewer.data().set_mesh(KV, KF);
308 // viewer.launch();
309 // }
310
311 // igl::opengl::glfw::Viewer viewer;
312 // viewer.data().set_mesh(collocation_points, CF);
313 // viewer.data().add_points(kernel_centers, Eigen::RowVector3d(0,1,1));
314 // for (int lf = 0; lf < mesh.n_cell_faces(element_index); ++lf) {
315 // Eigen::MatrixXd samples;
316 // samples = UV.middleRows(uv_ranges(lf), uv_ranges(lf+1) - uv_ranges(lf));
317 // Eigen::RowVector3d c = Eigen::RowVector3d::Random();
318 // viewer.data().add_points(samples, c);
319 // }
320 // viewer.launch();
321
322 // Compute right-hand side constraints for setting the harmonic kernels
323 Eigen::MatrixXd samples;
324 std::vector<AssemblyValues> basis_val;
325 rhs.resize(UV.rows(), local_to_global.size());
326 rhs.setZero();
327 for (int lf = 0; lf < mesh.n_cell_faces(element_index); ++lf)
328 {
329 auto index = mesh.get_index_from_element(element_index, lf, 0);
330 const int c2 = mesh.switch_element(index).element;
331 assert(c2 >= 0); // no boundary polytope
332
333 const InterfaceData &bdata = poly_face_to_data.at(index.face);
334 const ElementBases &b = bases[c2];
335
336 samples = UV.middleRows(uv_ranges(lf), uv_ranges(lf + 1) - uv_ranges(lf));
337 b.evaluate_bases(samples, basis_val);
338
339 // Evaluate field basis and set up the rhs
340 for (int other_local_basis_id : bdata.local_indices)
341 {
342 // b.bases[other_local_basis_id].basis(samples, basis_val);
343
344 for (const auto &x : b.bases[other_local_basis_id].global())
345 {
346 const int global_node_id = x.index;
347 const double weight = x.val;
348
349 const int poly_local_basis_id = std::distance(local_to_global.begin(),
350 std::find(local_to_global.begin(), local_to_global.end(), global_node_id));
351 rhs.block(uv_ranges(lf), poly_local_basis_id, basis_val[other_local_basis_id].val.size(), 1) += basis_val[other_local_basis_id].val * weight;
352 }
353 }
354 }
355
356 // Compute quadrature points + normalize kernels and collocation points
357 Eigen::MatrixXd NV = triangulated_vertices;
358 // scaling = (NV.colwise().maxCoeff() - NV.colwise().minCoeff()).maxCoeff();
359 // translation = NV.colwise().minCoeff();
360 scaling = 1.0;
361 translation.setZero();
362 // NV = (NV.rowwise() - translation) / scaling;
363 PolyhedronQuadrature::get_quadrature(NV, triangulated_faces, mesh.kernel(element_index),
364 quadrature_order, quadrature);
365 PolyhedronQuadrature::get_quadrature(NV, triangulated_faces, mesh.kernel(element_index),
366 mass_quadrature_order, mass_quadrature);
367
368 // Normalization
369 // collocation_points = (collocation_points.rowwise() - translation) / scaling;
370 // kernel_centers = (kernel_centers.rowwise() - translation) / scaling;
371 // KV = (KV.rowwise() - translation) / scaling;
372
373 triangulated_vertices = KV;
374 triangulated_faces = KF;
375 }
376
377 } // anonymous namespace
378
380
381 // Compute the integral constraints for each basis of the mesh
383 const Assembler &assembler,
384 const Mesh3D &mesh,
385 const int n_bases,
386 const std::vector<ElementBases> &bases,
387 const std::vector<ElementBases> &gbases,
388 Eigen::MatrixXd &basis_integrals)
389 {
390 assert(mesh.is_volume());
391
392 basis_integrals.resize(n_bases, 9);
393 basis_integrals.setZero();
394 Eigen::MatrixXd rhs(n_bases, 9);
395 rhs.setZero();
396
397 const int n_elements = mesh.n_elements();
398 for (int e = 0; e < n_elements; ++e)
399 {
400 if (mesh.is_polytope(e))
401 {
402 continue;
403 }
404 // ElementAssemblyValues vals = values[e];
405 // const ElementAssemblyValues &gvals = gvalues[e];
407 vals.compute(e, mesh.is_volume(), bases[e], gbases[e]);
408
409 // Computes the discretized integral of the PDE over the element
410 const int n_local_bases = int(vals.basis_values.size());
411 for (int j = 0; j < n_local_bases; ++j)
412 {
413 const AssemblyValues &v = vals.basis_values[j];
414 const double integral_100 = (v.grad_t_m.col(0).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
415 const double integral_010 = (v.grad_t_m.col(1).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
416 const double integral_001 = (v.grad_t_m.col(2).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
417
418 const double integral_110 = ((vals.val.col(1).array() * v.grad_t_m.col(0).array() + vals.val.col(0).array() * v.grad_t_m.col(1).array()) * vals.det.array() * vals.quadrature.weights.array()).sum();
419 const double integral_011 = ((vals.val.col(2).array() * v.grad_t_m.col(1).array() + vals.val.col(1).array() * v.grad_t_m.col(2).array()) * vals.det.array() * vals.quadrature.weights.array()).sum();
420 const double integral_101 = ((vals.val.col(0).array() * v.grad_t_m.col(2).array() + vals.val.col(2).array() * v.grad_t_m.col(0).array()) * vals.det.array() * vals.quadrature.weights.array()).sum();
421
422 const double integral_200 = 2 * (vals.val.col(0).array() * v.grad_t_m.col(0).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
423 const double integral_020 = 2 * (vals.val.col(1).array() * v.grad_t_m.col(1).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
424 const double integral_002 = 2 * (vals.val.col(2).array() * v.grad_t_m.col(2).array() * vals.det.array() * vals.quadrature.weights.array()).sum();
425
426 const double area = (v.val.array() * vals.det.array() * vals.quadrature.weights.array()).sum();
427
428 for (size_t ii = 0; ii < v.global.size(); ++ii)
429 {
430 basis_integrals(v.global[ii].index, 0) += integral_100 * v.global[ii].val;
431 basis_integrals(v.global[ii].index, 1) += integral_010 * v.global[ii].val;
432 basis_integrals(v.global[ii].index, 2) += integral_001 * v.global[ii].val;
433
434 basis_integrals(v.global[ii].index, 3) += integral_110 * v.global[ii].val;
435 basis_integrals(v.global[ii].index, 4) += integral_011 * v.global[ii].val;
436 basis_integrals(v.global[ii].index, 5) += integral_101 * v.global[ii].val;
437
438 basis_integrals(v.global[ii].index, 6) += integral_200 * v.global[ii].val;
439 basis_integrals(v.global[ii].index, 7) += integral_020 * v.global[ii].val;
440 basis_integrals(v.global[ii].index, 8) += integral_002 * v.global[ii].val;
441
442 rhs(v.global[ii].index, 6) += -2.0 * area * v.global[ii].val;
443 rhs(v.global[ii].index, 7) += -2.0 * area * v.global[ii].val;
444 rhs(v.global[ii].index, 8) += -2.0 * area * v.global[ii].val;
445 }
446 }
447 }
448
449 basis_integrals -= rhs;
450 }
451
452 // -----------------------------------------------------------------------------
453
454 // Distance from harmonic kernels to polygon boundary
455 double compute_epsilon(const Mesh3D &mesh, int e)
456 {
457 // double area = 0;
458 // const int n_edges = mesh.n_element_vertices(e);
459 // for (int i = 0; i < n_edges; ++i) {
460 // const int ip = (i + 1) % n_edges;
461
462 // Eigen::RowVector2d p0 = mesh.point(mesh.vertex_global_index(e, i));
463 // Eigen::RowVector2d p1 = mesh.point(mesh.vertex_global_index(e, ip));
464
465 // Eigen::Matrix2d det_mat;
466 // det_mat.row(0) = p0;
467 // det_mat.row(1) = p1;
468
469 // area += det_mat.determinant();
470 // }
471 // area = std::fabs(area);
472 // // const double eps = use_harmonic ? (0.08*area) : 0;
473 // const double eps = 0.08*area;
474
475 return 0.1; // will be relative to the volume of the poly
476 }
477
479 const LinearAssembler &assembler,
480 const int nn_samples_per_edge,
481 const Mesh3D &mesh,
482 const int n_bases,
483 const int quadrature_order,
484 const int mass_quadrature_order,
485 const int integral_constraints,
486 std::vector<ElementBases> &bases,
487 const std::vector<ElementBases> &gbases,
488 const std::map<int, InterfaceData> &poly_face_to_data,
489 std::map<int, std::pair<Eigen::MatrixXd, Eigen::MatrixXi>> &mapped_boundary)
490 {
491 assert(mesh.is_volume());
492 if (poly_face_to_data.empty())
493 {
494 return 0;
495 }
496 int n_kernels_per_edge = 4; //(int) std::round(n_samples_per_edge / 3.0);
497 int n_samples_per_edge = 3 * n_kernels_per_edge;
498
499 // Step 1: Compute integral constraints
500 Eigen::MatrixXd basis_integrals;
501 compute_integral_constraints(assembler, mesh, n_bases, bases, gbases, basis_integrals);
502
503 // Step 2: Compute the rest =)
504 for (int e = 0; e < mesh.n_elements(); ++e)
505 {
506 if (!mesh.is_polytope(e))
507 {
508 continue;
509 }
510 // No boundary polytope
511 // assert(element_type[e] != ElementType::BOUNDARY_POLYTOPE);
512
513 // Kernel distance to polygon boundary
514 const double eps = compute_epsilon(mesh, e);
515
516 std::vector<int> local_to_global; // map local basis id (the ones that are nonzero on the polygon boundary) to global basis id
517 Eigen::MatrixXd collocation_points, kernel_centers, triangulated_vertices;
518 Eigen::MatrixXi triangulated_faces;
519 Eigen::MatrixXd rhs; // 1 row per collocation point, 1 column per basis that is nonzero on the polygon boundary
520
521 ElementBases &b = bases[e];
522 b.has_parameterization = false;
523
524 Quadrature tmp_quadrature, tmp_mass_quadrature;
525 double scaling;
526 Eigen::RowVector3d translation;
527 sample_polyhedra(e, 2, n_kernels_per_edge, n_samples_per_edge,
528 quadrature_order > 0 ? quadrature_order : AssemblerUtils::quadrature_order(assembler.name(), 2, AssemblerUtils::BasisType::POLY, 3),
529 mass_quadrature_order > 0 ? mass_quadrature_order : AssemblerUtils::quadrature_order("Mass", 2, AssemblerUtils::BasisType::POLY, 3),
530 mesh, poly_face_to_data, bases, gbases, eps, local_to_global,
531 collocation_points, kernel_centers, rhs, triangulated_vertices,
532 triangulated_faces, tmp_quadrature, tmp_mass_quadrature, scaling, translation);
533
534 b.set_quadrature([tmp_quadrature](Quadrature &quad) { quad = tmp_quadrature; });
535 b.set_mass_quadrature([tmp_mass_quadrature](Quadrature &quad) { quad = tmp_mass_quadrature; });
536 // b.scaling_ = scaling;
537 // b.translation_ = translation;
538
539 // igl::opengl::glfw::Viewer & viewer = UIState::ui_state().viewer;
540 // viewer.data().clear();
541 // viewer.data().set_mesh(triangulated_vertices, triangulated_faces);
542 // viewer.data().add_points(kernel_centers, Eigen::Vector3d(0,1,1).transpose());
543 // add_spheres(viewer, kernel_centers, 0.005);
544
545 // Eigen::MatrixXd pts = triangulated_vertices, normals;
546 // Eigen::MatrixXi tris = triangulated_faces;
547 // igl::per_corner_normals(pts, tris, 20, normals);
548 // viewer.data().set_normals(normals);
549 // viewer.data().set_face_based(false);
550 // viewer.launch();
551
552 // for(int a = 0; rhs.cols();++a)
553 // {
554 // igl::opengl::glfw::Viewer viewer;
555 // Eigen::MatrixXd asd(collocation_points.rows(), 3);
556 // asd.col(0)=collocation_points.col(0);
557 // asd.col(1)=collocation_points.col(1);
558 // asd.col(2)=collocation_points.col(2);
559 // Eigen::VectorXd S = rhs.col(a);
560 // Eigen::MatrixXd C;
561 // igl::colormap(igl::COLOR_MAP_TYPE_VIRIDIS, S, true, C);
562 // viewer.data().add_points(asd, C);
563 // viewer.launch();
564 // }
565
566 // for(int asd = 0; asd < collocation_points.rows(); ++asd) {
567 // viewer.data().add_label(collocation_points.row(asd), std::to_string(asd));
568 // }
569
570 // Compute the weights of the RBF kernels
571 Eigen::MatrixXd local_basis_integrals(rhs.cols(), basis_integrals.cols());
572 for (long k = 0; k < rhs.cols(); ++k)
573 {
574 local_basis_integrals.row(k) = -basis_integrals.row(local_to_global[k]);
575 }
576 auto set_rbf = [&b](auto rbf) {
577 b.set_bases_func([rbf](const Eigen::MatrixXd &uv, std::vector<AssemblyValues> &val) {
578 Eigen::MatrixXd tmp;
579 rbf->bases_values(uv, tmp);
580 val.resize(tmp.cols());
581 assert(tmp.rows() == uv.rows());
582
583 for (size_t i = 0; i < tmp.cols(); ++i)
584 {
585 val[i].val = tmp.col(i);
586 }
587 });
588 b.set_grads_func([rbf](const Eigen::MatrixXd &uv, std::vector<AssemblyValues> &val) {
589 Eigen::MatrixXd tmpx, tmpy, tmpz;
590
591 rbf->bases_grads(0, uv, tmpx);
592 rbf->bases_grads(1, uv, tmpy);
593 rbf->bases_grads(2, uv, tmpz);
594
595 val.resize(tmpx.cols());
596 assert(tmpx.cols() == tmpy.cols());
597 assert(tmpx.cols() == tmpz.cols());
598 assert(tmpx.rows() == uv.rows());
599 for (size_t i = 0; i < tmpx.cols(); ++i)
600 {
601 val[i].grad.resize(uv.rows(), uv.cols());
602 val[i].grad.col(0) = tmpx.col(i);
603 val[i].grad.col(1) = tmpy.col(i);
604 val[i].grad.col(2) = tmpz.col(i);
605 }
606 });
607 };
608 if (integral_constraints == 0)
609 {
610 set_rbf(std::make_shared<RBFWithLinear>(
611 kernel_centers, collocation_points, local_basis_integrals, tmp_quadrature, rhs, false));
612 }
613 else if (integral_constraints == 1)
614 {
615 set_rbf(std::make_shared<RBFWithLinear>(
616 kernel_centers, collocation_points, local_basis_integrals, tmp_quadrature, rhs));
617 }
618 else if (integral_constraints == 2)
619 {
620 set_rbf(std::make_shared<RBFWithQuadratic>(
621 // set_rbf(std::make_shared<RBFWithQuadraticLagrange>(
622 assembler, kernel_centers, collocation_points, local_basis_integrals, tmp_quadrature, rhs));
623 }
624 else
625 {
626 throw std::runtime_error(fmt::format("Unsupported constraint order: {:d}", integral_constraints));
627 }
628
629 // Set the bases which are nonzero inside the polygon
630 const int n_poly_bases = int(local_to_global.size());
631 b.bases.resize(n_poly_bases);
632 for (int i = 0; i < n_poly_bases; ++i)
633 {
634 b.bases[i].init(-2, local_to_global[i], i, Eigen::MatrixXd::Constant(1, 3, std::nan("")));
635 }
636
637 // Polygon boundary after geometric mapping from neighboring elements
638 orient_closed_surface(triangulated_vertices, triangulated_faces, false); // stupid viewer is flipping all the faces
639 mapped_boundary[e].first = triangulated_vertices;
640 mapped_boundary[e].second = triangulated_faces;
641 }
642
643 return 0;
644 }
645 } // namespace basis
646} // namespace polyfem
vector< list< int > > adj
int V
double val
Definition Assembler.cpp:86
ElementAssemblyValues vals
Definition Assembler.cpp:22
Quadrature quadrature
int x
virtual std::string name() const =0
static int quadrature_order(const std::string &assembler, const int basis_degree, const BasisType &b_type, const int dim)
utility for retrieving the needed quadrature order to precisely integrate the given form on the given...
stores per local bases evaluations
std::vector< basis::Local2Global > global
stores per element basis values at given quadrature points and geometric mapping
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,...
assemble matrix based on the local assembler local assembler is eg Laplace, LinearElasticity etc
Stores the basis functions for a given element in a mesh (facet in 2d, cell in 3d).
void set_bases_func(EvalBasesFunc fun)
void set_quadrature(const QuadratureFunction &fun)
void set_grads_func(EvalBasesFunc fun)
std::vector< Basis > bases
one basis function per node in the element
void set_mass_quadrature(const QuadratureFunction &fun)
static Eigen::VectorXi hex_face_local_nodes(const bool serendipity, const int q, const mesh::Mesh3D &mesh, mesh::Navigation3D::Index index)
static int build_bases(const assembler::LinearAssembler &assembler, const int n_samples_per_edge, const mesh::Mesh3D &mesh, const int n_bases, const int quadrature_order, const int mass_quadrature_order, const int integral_constraints, std::vector< ElementBases > &bases, const std::vector< ElementBases > &gbases, const std::map< int, InterfaceData > &poly_face_to_data, std::map< int, std::pair< Eigen::MatrixXd, Eigen::MatrixXi > > &mapped_boundary)
Build bases over the remaining polygons of a mesh.
static void compute_integral_constraints(const assembler::Assembler &assembler, const mesh::Mesh3D &mesh, const int n_bases, const std::vector< ElementBases > &bases, const std::vector< ElementBases > &gbases, Eigen::MatrixXd &basis_integrals)
bool is_volume() const override
checks if mesh is volume
Definition Mesh3D.hpp:28
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
Definition Mesh.hpp:161
bool is_polytope(const int el_id) const
checks if element is polygon compatible
Definition Mesh.cpp:365
static void get_quadrature(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::RowVector3d &kernel, const int order, Quadrature &quadr)
Gets the quadrature points & weights for a polyhedron.
list indices
Definition p_bases.py:258
void q_nodes_3d(const int q, Eigen::MatrixXd &val)
double compute_epsilon(const Mesh2D &mesh, int e)
void sample_surface(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, int num_samples, Eigen::MatrixXd &P, Eigen::MatrixXd *N=nullptr, int num_lloyd=10, int num_newton=10)
Samples points on a surface.
void orient_closed_surface(const Eigen::MatrixXd &V, Eigen::MatrixXi &F, bool positive=true)
Orient a triangulated surface to have positive volume.
double signed_volume(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F)
Compute the signed volume of a surface mesh.
void reorder_mesh(Eigen::MatrixXd &V, Eigen::MatrixXi &F, const Eigen::VectorXi &C, Eigen::VectorXi &R)
Reorder vertices of a mesh using color tags, so that vertices are ordered by increasing colors.
void signed_squared_distances(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const Eigen::MatrixXd &P, Eigen::VectorXd &D)
Computes the signed squared distance from a list of points to a triangle mesh.
bool instantiate_pattern(const Eigen::MatrixXd &IV, const Eigen::MatrixXi &IF, const Eigen::MatrixXd &PV, const Eigen::MatrixXi &PF, Eigen::MatrixXd &OV, Eigen::MatrixXi &OF, Eigen::VectorXi *SF=nullptr, EvalParametersFunc evalFunc=nullptr, GetAdjacentLocalEdge getAdjLocalEdge=nullptr)
std::function< void(const Eigen::MatrixXd &, Eigen::MatrixXd &, int)> EvalParametersFunc
std::function< std::tuple< int, int, bool >(int, int)> GetAdjacentLocalEdge
void regular_2d_grid(const int n, bool tri, Eigen::MatrixXd &V, Eigen::MatrixXi &F)
Generate a canonical triangle/quad subdivided from a regular grid.