Loading [MathJax]/extensions/tex2jax.js
PolyFEM
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SplineBasis2d.cpp
Go to the documentation of this file.
1#include "SplineBasis2d.hpp"
2
3#include "LagrangeBasis2d.hpp"
4#include "function/QuadraticBSpline2d.hpp"
5
6#include <polyfem/quadrature/QuadQuadrature.hpp>
7#include <polyfem/mesh/MeshNodes.hpp>
8
9#include <polyfem/assembler/AssemblerUtils.hpp>
10
11#include <polysolve/linear/Solver.hpp>
12
13#include <polyfem/utils/Types.hpp>
14
15#include <polyfem/Common.hpp>
16#include <polyfem/autogen/auto_q_bases.hpp>
17
18#include <Eigen/Sparse>
19
20#include <cassert>
21#include <iostream>
22#include <vector>
23#include <array>
24#include <map>
25
26// TODO carefull with simplices
27
28namespace polyfem
29{
30 using namespace polysolve;
31 using namespace Eigen;
32 using namespace assembler;
33 using namespace mesh;
34 using namespace quadrature;
35
36 namespace basis
37 {
38 namespace
39 {
40 typedef Matrix<std::vector<int>, 3, 3> SpaceMatrix;
41 typedef Matrix<RowVectorNd, 3, 3> NodeMatrix;
42
43 void print_local_space(const SpaceMatrix &space)
44 {
45 std::stringstream ss;
46 ss << std::endl;
47 for (int j = 2; j >= 0; --j)
48 {
49 for (int i = 0; i < 3; ++i)
50 {
51 if (space(i, j).size() > 0)
52 {
53 for (std::size_t l = 0; l < space(i, j).size(); ++l)
54 ss << space(i, j)[l] << ",";
55
56 ss << "\t";
57 }
58 else
59 ss << "x\t";
60 }
61 ss << std::endl;
62 }
63
64 logger().trace("Local space:\n{}", ss.str());
65 }
66
67 int node_id_from_edge_index(const Mesh2D &mesh, MeshNodes &mesh_nodes, const Navigation::Index &index)
68 {
69 const int face_id = mesh.switch_face(index).face;
70 if (face_id >= 0 && mesh.is_cube(face_id))
71 return mesh_nodes.node_id_from_face(face_id);
72
73 return mesh_nodes.node_id_from_edge(index.edge);
74 }
75
76 void explore_direction(const Navigation::Index &index, const Mesh2D &mesh, MeshNodes &mesh_nodes, const int x, const int y, const bool is_x, const bool invert, LocalBoundary &lb, SpaceMatrix &space, NodeMatrix &node, std::map<int, InterfaceData> &poly_edge_to_data)
77 {
78 int node_id = node_id_from_edge_index(mesh, mesh_nodes, index);
79 // bool real_boundary = mesh.node_id_from_edge_index(index, node_id);
80
81 assert(std::find(space(x, y).begin(), space(x, y).end(), node_id) == space(x, y).end());
82 space(x, y).push_back(node_id);
83 node(x, y) = mesh_nodes.node_position(node_id);
84 assert(node(x, y).size() == 2);
85
86 const int x1 = is_x ? x : (invert ? 2 : 0);
87 const int y1 = !is_x ? y : (invert ? 2 : 0);
88
89 const int x2 = is_x ? x : (invert ? 0 : 2);
90 const int y2 = !is_x ? y : (invert ? 0 : 2);
91
92 const bool is_boundary = mesh_nodes.is_boundary(node_id);
93 const bool is_interface = mesh_nodes.is_interface(node_id);
94
95 if (is_boundary)
96 {
97 lb.add_boundary_primitive(index.edge, LagrangeBasis2d::quad_edge_local_nodes(2, mesh, index)[1] - 4);
98 // lb.add_boundary_primitive(index.edge, LagrangeBasis2d::quadr_quad_edge_local_nodes(mesh, index)[1]-4);
99 // bounday_nodes.push_back(node_id);
100 }
101 else if (is_interface)
102 {
103 InterfaceData &data = poly_edge_to_data[index.edge];
104 data.local_indices.push_back(y * 3 + x);
105 }
106 else
107 {
108 assert(!is_boundary && !is_interface);
109
110 Navigation::Index start_index = mesh.switch_face(index);
111 assert(start_index.vertex == index.vertex);
112 assert(start_index.face >= 0);
113
114 Navigation::Index edge1 = mesh.switch_edge(start_index);
115 node_id = node_id_from_edge_index(mesh, mesh_nodes, edge1);
116 // if(mesh_nodes.is_boundary(node_id))
117 // bounday_nodes.push_back(node_id);
118
119 if (std::find(space(x1, y1).begin(), space(x1, y1).end(), node_id) == space(x1, y1).end())
120 {
121 space(x1, y1).push_back(node_id);
122 node(x1, y1) = mesh_nodes.node_position(node_id);
123 }
124
125 Navigation::Index edge2 = mesh.switch_edge(mesh.switch_vertex(start_index));
126 node_id = node_id_from_edge_index(mesh, mesh_nodes, edge2);
127 // if(mesh_nodes.is_boundary(node_id))
128 // bounday_nodes.push_back(node_id);
129 if (std::find(space(x2, y2).begin(), space(x2, y2).end(), node_id) == space(x2, y2).end())
130 {
131 space(x2, y2).push_back(node_id);
132 node(x2, y2) = mesh_nodes.node_position(node_id);
133 // node(x2, y2).push_back(mesh.node_from_edge_index(edge2));
134 }
135 }
136 }
137
138 void add_id_for_poly(const Navigation::Index &index, const int x1, const int y1, const int x2, const int y2, const SpaceMatrix &space, std::map<int, InterfaceData> &poly_edge_to_data)
139 {
140 auto it = poly_edge_to_data.find(index.edge);
141 if (it != poly_edge_to_data.end())
142 {
143 InterfaceData &data = it->second;
144
145 assert(space(x1, y1).size() == 1);
146 data.local_indices.push_back(y1 * 3 + x1);
147
148 assert(space(x2, y2).size() == 1);
149 data.local_indices.push_back(y2 * 3 + x2);
150 }
151 }
152
153 void build_local_space(const Mesh2D &mesh, MeshNodes &mesh_nodes, const int el_index, SpaceMatrix &space, NodeMatrix &node, std::vector<LocalBoundary> &local_boundary, std::map<int, InterfaceData> &poly_edge_to_data)
154 {
155 assert(!mesh.is_volume());
156
157 Navigation::Index index;
158 // space.setConstant(-1);
159
160 const int el_node_id = mesh_nodes.node_id_from_face(el_index);
161 space(1, 1).push_back(el_node_id);
162 node(1, 1) = mesh_nodes.node_position(el_node_id);
163 // (mesh.node_from_face(el_index));
164
165 LocalBoundary lb(el_index, BoundaryType::QUAD_LINE);
166
168 index = mesh.get_index_from_face(el_index);
169 explore_direction(index, mesh, mesh_nodes, 1, 0, false, false, lb, space, node, poly_edge_to_data);
170
172 index = mesh.next_around_face(index);
173 explore_direction(index, mesh, mesh_nodes, 2, 1, true, false, lb, space, node, poly_edge_to_data);
174
176 index = mesh.next_around_face(index);
177 explore_direction(index, mesh, mesh_nodes, 1, 2, false, true, lb, space, node, poly_edge_to_data);
178
180 index = mesh.next_around_face(index);
181 explore_direction(index, mesh, mesh_nodes, 0, 1, true, true, lb, space, node, poly_edge_to_data);
182
184 if (mesh_nodes.is_boundary_or_interface(space(1, 2).front()) && mesh_nodes.is_boundary_or_interface(space(2, 1).front()))
185 {
186 assert(space(2, 2).empty());
187
188 Navigation::Index start_index = mesh.get_index_from_face(el_index);
189 start_index = mesh.next_around_face(start_index);
190 start_index = mesh.next_around_face(start_index);
191
192 const int node_id = mesh_nodes.node_id_from_vertex(start_index.vertex);
193 // mesh.vertex_node_id(start_index.vertex);
194 space(2, 2).push_back(node_id);
195 node(2, 2) = mesh_nodes.node_position(node_id);
196 // node(2,2).push_back(mesh.node_from_vertex(start_index.vertex));
197
198 // bounday_nodes.push_back(node_id);
199 }
200
201 if (mesh_nodes.is_boundary_or_interface(space(1, 0).front()) && mesh_nodes.is_boundary_or_interface(space(2, 1).front()))
202 {
203 assert(space(2, 0).empty());
204
205 Navigation::Index start_index = mesh.get_index_from_face(el_index);
206 start_index = mesh.next_around_face(start_index);
207
208 const int node_id = mesh_nodes.node_id_from_vertex(start_index.vertex);
209 // mesh.vertex_node_id(start_index.vertex);
210 space(2, 0).push_back(node_id);
211 node(2, 0) = mesh_nodes.node_position(node_id);
212 // .push_back(mesh.node_from_vertex(start_index.vertex));
213
214 // bounday_nodes.push_back(node_id);
215 }
216
217 if (mesh_nodes.is_boundary_or_interface(space(1, 2).front()) && mesh_nodes.is_boundary_or_interface(space(0, 1).front()))
218 {
219 assert(space(0, 2).empty());
220
221 Navigation::Index start_index = mesh.get_index_from_face(el_index);
222 start_index = mesh.next_around_face(start_index);
223 start_index = mesh.next_around_face(start_index);
224 start_index = mesh.next_around_face(start_index);
225
226 // const int node_id = mesh.vertex_node_id(start_index.vertex);
227 const int node_id = mesh_nodes.node_id_from_vertex(start_index.vertex);
228 space(0, 2).push_back(node_id);
229 node(0, 2) = mesh_nodes.node_position(node_id);
230 // .push_back(mesh.node_from_vertex(start_index.vertex));
231
232 // bounday_nodes.push_back(node_id);
233 }
234
235 if (mesh_nodes.is_boundary_or_interface(space(1, 0).front()) && mesh_nodes.is_boundary_or_interface(space(0, 1).front()))
236 {
237 Navigation::Index start_index = mesh.get_index_from_face(el_index);
238
239 // const int node_id = mesh.vertex_node_id(start_index.vertex);
240 const int node_id = mesh_nodes.node_id_from_vertex(start_index.vertex);
241 space(0, 0).push_back(node_id);
242 node(0, 0) = mesh_nodes.node_position(node_id);
243 //.push_back(mesh.node_from_vertex(start_index.vertex));
244
245 // bounday_nodes.push_back(node_id);
246 }
247
248 // print_local_space(space);
249
251 index = mesh.get_index_from_face(el_index);
252 add_id_for_poly(index, 0, 0, 2, 0, space, poly_edge_to_data);
253
254 index = mesh.next_around_face(index);
255 add_id_for_poly(index, 2, 0, 2, 2, space, poly_edge_to_data);
256
257 index = mesh.next_around_face(index);
258 add_id_for_poly(index, 2, 2, 0, 2, space, poly_edge_to_data);
259
260 index = mesh.next_around_face(index);
261 add_id_for_poly(index, 0, 2, 0, 0, space, poly_edge_to_data);
262
263 if (!lb.empty())
264 local_boundary.emplace_back(lb);
265 }
266
267 void setup_knots_vectors(MeshNodes &mesh_nodes, const SpaceMatrix &space, std::array<std::array<double, 4>, 3> &h_knots, std::array<std::array<double, 4>, 3> &v_knots)
268 {
269 // left and right neigh are absent
270 if (mesh_nodes.is_boundary_or_interface(space(0, 1).front()) && mesh_nodes.is_boundary_or_interface(space(2, 1).front()))
271 {
272 h_knots[0] = {{0, 0, 0, 1}};
273 h_knots[1] = {{0, 0, 1, 1}};
274 h_knots[2] = {{0, 1, 1, 1}};
275 }
276 // left neigh is absent
277 else if (mesh_nodes.is_boundary_or_interface(space(0, 1).front()))
278 {
279 h_knots[0] = {{0, 0, 0, 1}};
280 h_knots[1] = {{0, 0, 1, 2}};
281 h_knots[2] = {{0, 1, 2, 3}};
282 }
283 // right neigh is absent
284 else if (mesh_nodes.is_boundary_or_interface(space(2, 1).front()))
285 {
286 h_knots[0] = {{-2, -1, 0, 1}};
287 h_knots[1] = {{-1, 0, 1, 1}};
288 h_knots[2] = {{0, 1, 1, 1}};
289 }
290 else
291 {
292 h_knots[0] = {{-2, -1, 0, 1}};
293 h_knots[1] = {{-1, 0, 1, 2}};
294 h_knots[2] = {{0, 1, 2, 3}};
295 }
296
297 // top and bottom neigh are absent
298 if (mesh_nodes.is_boundary_or_interface(space(1, 0).front()) && mesh_nodes.is_boundary_or_interface(space(1, 2).front()))
299 {
300 v_knots[0] = {{0, 0, 0, 1}};
301 v_knots[1] = {{0, 0, 1, 1}};
302 v_knots[2] = {{0, 1, 1, 1}};
303 }
304 // bottom neigh is absent
305 else if (mesh_nodes.is_boundary_or_interface(space(1, 0).front()))
306 {
307 v_knots[0] = {{0, 0, 0, 1}};
308 v_knots[1] = {{0, 0, 1, 2}};
309 v_knots[2] = {{0, 1, 2, 3}};
310 }
311 // top neigh is absent
312 else if (mesh_nodes.is_boundary_or_interface(space(1, 2).front()))
313 {
314 v_knots[0] = {{-2, -1, 0, 1}};
315 v_knots[1] = {{-1, 0, 1, 1}};
316 v_knots[2] = {{0, 1, 1, 1}};
317 }
318 else
319 {
320 v_knots[0] = {{-2, -1, 0, 1}};
321 v_knots[1] = {{-1, 0, 1, 2}};
322 v_knots[2] = {{0, 1, 2, 3}};
323 }
324 }
325
326 void basis_for_regular_quad(const SpaceMatrix &space, const NodeMatrix &loc_nodes, const std::array<std::array<double, 4>, 3> &h_knots, const std::array<std::array<double, 4>, 3> &v_knots, ElementBases &b)
327 {
328 for (int y = 0; y < 3; ++y)
329 {
330 for (int x = 0; x < 3; ++x)
331 {
332 if (space(x, y).size() == 1)
333 {
334 const int global_index = space(x, y).front();
335 const Eigen::MatrixXd &node = loc_nodes(x, y);
336 assert(node.size() == 2);
337
338 const int local_index = y * 3 + x;
339 b.bases[local_index].init(2, global_index, local_index, node);
340
341 const QuadraticBSpline2d spline(h_knots[x], v_knots[y]);
342 b.bases[local_index].set_basis([spline](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { spline.interpolate(uv, val); });
343 b.bases[local_index].set_grad([spline](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { spline.derivative(uv, val); });
344 }
345 }
346 }
347 }
348
349 void basis_for_irregulard_quad(const int el_id, const Mesh2D &mesh, MeshNodes &mesh_nodes, const SpaceMatrix &space, const NodeMatrix &loc_nodes, const std::array<std::array<double, 4>, 3> &h_knots, const std::array<std::array<double, 4>, 3> &v_knots, ElementBases &b)
350 {
351 for (int y = 0; y < 3; ++y)
352 {
353 for (int x = 0; x < 3; ++x)
354 {
355 if (space(x, y).size() > 1)
356 {
357 const int mpx = 1;
358 const int mpy = y;
359
360 const int mmx = x;
361 const int mmy = 1;
362
363 std::vector<int> other_indices;
364 const auto &center = b.bases[1 * 3 + 1].global().front();
365
366 const auto &el1 = b.bases[mpy * 3 + mpx].global().front();
367 const auto &el2 = b.bases[mmy * 3 + mmx].global().front();
368
369 Navigation::Index start_index = mesh.get_index_from_face(el_id);
370 bool found = false;
371 for (int i = 0; i < 4; ++i)
372 {
373 other_indices.clear();
374 int n_neighs = 0;
375 Navigation::Index index = start_index;
376 do
377 {
378 const int f_index = mesh_nodes.node_id_from_face(index.face);
379 if (f_index != el1.index && f_index != el2.index && f_index != center.index)
380 other_indices.push_back(f_index);
381
382 ++n_neighs;
383 index = mesh.next_around_vertex(index);
384 } while (index.face != start_index.face);
385 if (n_neighs != 4)
386 {
387 found = true;
388 break;
389 }
390
391 start_index = mesh.next_around_face(start_index);
392 }
393 assert(found);
394
395 const int local_index = y * 3 + x;
396 auto &base = b.bases[local_index];
397
398 const int k = int(other_indices.size()) + 3;
399
400 base.global().resize(k);
401
402 base.global()[0].index = center.index;
403 base.global()[0].val = (4. - k) / k;
404 base.global()[0].node = center.node;
405
406 base.global()[1].index = el1.index;
407 base.global()[1].val = (4. - k) / k;
408 base.global()[1].node = el1.node;
409
410 base.global()[2].index = el2.index;
411 base.global()[2].val = (4. - k) / k;
412 base.global()[2].node = el2.node;
413
414 for (std::size_t n = 0; n < other_indices.size(); ++n)
415 {
416 base.global()[3 + n].index = other_indices[n];
417 base.global()[3 + n].val = 4. / k;
418 base.global()[3 + n].node = mesh_nodes.node_position(other_indices[n]);
419 }
420
421 const QuadraticBSpline2d spline(h_knots[x], v_knots[y]);
422 b.bases[local_index].set_basis([spline](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { spline.interpolate(uv, val); });
423 b.bases[local_index].set_grad([spline](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { spline.derivative(uv, val); });
424 }
425 }
426 }
427 }
428
429 void create_q2_nodes(const Mesh2D &mesh, const int el_index, std::set<int> &vertex_id, std::set<int> &edge_id, ElementBases &b, std::vector<LocalBoundary> &local_boundary, int &n_bases)
430 {
431 b.bases.resize(9);
432
433 LocalBoundary lb(el_index, BoundaryType::QUAD_LINE);
434
435 Navigation::Index index = mesh.get_index_from_face(el_index);
436 for (int j = 0; j < 4; ++j)
437 {
438 int current_vertex_node_id = -1;
439 int current_edge_node_id = -1;
440 Eigen::Matrix<double, 1, 2> current_edge_node;
441 Eigen::MatrixXd current_vertex_node;
442
443 // auto e2l = LagrangeBasis2d::quadr_quad_edge_local_nodes(mesh, index);
444 auto e2l = LagrangeBasis2d::quad_edge_local_nodes(2, mesh, index);
445
446 int vertex_basis_id = e2l[0];
447 int edge_basis_id = e2l[1];
448
449 const int opposite_face = mesh.switch_face(index).face;
450
451 // if the edge/vertex is boundary the it is a Q2 edge
452 bool is_vertex_q2 = true;
453 bool is_vertex_boundary = false;
454
455 Navigation::Index vindex = index;
456
457 do
458 {
459 if (vindex.face < 0)
460 {
461 is_vertex_boundary = true;
462 break;
463 }
464 if (mesh.is_spline_compatible(vindex.face))
465 {
466 is_vertex_q2 = false;
467 break;
468 }
469 vindex = mesh.next_around_vertex(vindex);
470 } while (vindex.edge != index.edge);
471
472 if (is_vertex_q2)
473 {
474 vindex = mesh.switch_face(index);
475 do
476 {
477 if (vindex.face < 0)
478 {
479 is_vertex_boundary = true;
480 break;
481 }
482
483 if (mesh.is_spline_compatible(vindex.face))
484 {
485 is_vertex_q2 = false;
486 break;
487 }
488 vindex = mesh.next_around_vertex(vindex);
489 } while (vindex.edge != index.edge);
490 }
491
492 const bool is_edge_q2 = opposite_face < 0 || !mesh.is_spline_compatible(opposite_face);
493
494 if (is_edge_q2)
495 {
496 const bool is_new_edge = edge_id.insert(index.edge).second;
497
498 if (is_new_edge)
499 {
500 current_edge_node_id = n_bases++;
501 current_edge_node = mesh.edge_barycenter(index.edge);
502
503 if (opposite_face < 0)
504 {
505 // bounday_nodes.push_back(current_edge_node_id);
506 lb.add_boundary_primitive(index.edge, edge_basis_id - 4);
507 }
508 }
509 }
510
511 if (is_vertex_q2)
512 {
513 assert(is_edge_q2);
514 const bool is_new_vertex = vertex_id.insert(index.vertex).second;
515
516 if (is_new_vertex)
517 {
518 current_vertex_node_id = n_bases++;
519 current_vertex_node = mesh.point(index.vertex);
520
521 // if(is_vertex_boundary)//mesh.is_vertex_boundary(index.vertex))
522 // bounday_nodes.push_back(current_vertex_node_id);
523 }
524 }
525
526 // init new Q2 nodes
527 if (current_vertex_node_id >= 0)
528 b.bases[vertex_basis_id].init(2, current_vertex_node_id, vertex_basis_id, current_vertex_node);
529
530 if (current_edge_node_id >= 0)
531 b.bases[edge_basis_id].init(2, current_edge_node_id, edge_basis_id, current_edge_node);
532
533 // set the basis functions
534 b.bases[vertex_basis_id].set_basis([vertex_basis_id](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { polyfem::autogen::q_basis_value_2d(2, vertex_basis_id, uv, val); });
535 b.bases[vertex_basis_id].set_grad([vertex_basis_id](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { polyfem::autogen::q_grad_basis_value_2d(2, vertex_basis_id, uv, val); });
536
537 b.bases[edge_basis_id].set_basis([edge_basis_id](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { polyfem::autogen::q_basis_value_2d(2, edge_basis_id, uv, val); });
538 b.bases[edge_basis_id].set_grad([edge_basis_id](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { polyfem::autogen::q_grad_basis_value_2d(2, edge_basis_id, uv, val); });
539
540 index = mesh.next_around_face(index);
541 }
542
543 // central node always present
544 const int face_basis_id = 8;
545 b.bases[face_basis_id].init(2, n_bases++, face_basis_id, mesh.face_barycenter(el_index));
546 b.bases[face_basis_id].set_basis([face_basis_id](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { polyfem::autogen::q_basis_value_2d(2, face_basis_id, uv, val); });
547 b.bases[face_basis_id].set_grad([face_basis_id](const Eigen::MatrixXd &uv, Eigen::MatrixXd &val) { polyfem::autogen::q_grad_basis_value_2d(2, face_basis_id, uv, val); });
548
549 if (!lb.empty())
550 local_boundary.emplace_back(lb);
551 }
552
553 void insert_into_global(const Local2Global &data, std::vector<Local2Global> &vec)
554 {
555 // ignore small weights
556 if (fabs(data.val) < 1e-10)
557 return;
558
559 bool found = false;
560
561 for (std::size_t i = 0; i < vec.size(); ++i)
562 {
563 if (vec[i].index == data.index)
564 {
565 // logger().trace("{} {} {}", vec[i].val, data.val, fabs(vec[i].val - data.val));
566 assert(fabs(vec[i].val - data.val) < 1e-10);
567 assert((vec[i].node - data.node).norm() < 1e-10);
568 found = true;
569 break;
570 }
571 }
572
573 if (!found)
574 vec.push_back(data);
575 }
576
577 void assign_q2_weights(const Mesh2D &mesh, const int el_index, std::vector<ElementBases> &bases)
578 {
579 // Eigen::MatrixXd eval_p;
580 std::vector<AssemblyValues> eval_p;
581 Navigation::Index index = mesh.get_index_from_face(el_index);
582 ElementBases &b = bases[el_index];
583
584 for (int j = 0; j < 4; ++j)
585 {
586 const int opposite_face = mesh.switch_face(index).face;
587
588 if (opposite_face < 0 || !mesh.is_cube(opposite_face))
589 {
590 index = mesh.next_around_face(index);
591 continue;
592 }
593
594 // const auto param_p = LagrangeBasis2d::quadr_quad_edge_local_nodes_coordinates(mesh, mesh.switch_face(index));
595 // const auto indices = LagrangeBasis2d::quadr_quad_edge_local_nodes(mesh, index);
596
597 const auto indices = LagrangeBasis2d::quad_edge_local_nodes(2, mesh, index);
598 Eigen::Matrix<double, 3, 2> param_p;
599
600 {
601 Eigen::MatrixXd quad_loc_nodes;
602 polyfem::autogen::q_nodes_2d(2, quad_loc_nodes);
603 const auto opposite_indices = LagrangeBasis2d::quad_edge_local_nodes(2, mesh, mesh.switch_face(index));
604 for (int k = 0; k < 3; ++k)
605 param_p.row(k) = quad_loc_nodes.row(opposite_indices[k]);
606 }
607
608 const int i0 = indices[0];
609 const int i1 = indices[1];
610 const int i2 = indices[2];
611
612 const auto &other_bases = bases[opposite_face];
613 other_bases.evaluate_bases(param_p, eval_p);
614
615 for (std::size_t i = 0; i < other_bases.bases.size(); ++i)
616 {
617 const auto &other_b = other_bases.bases[i];
618
619 if (other_b.global().empty())
620 continue;
621
622 assert(eval_p[i].val.size() == 3);
623
624 // basis i of element opposite face is zero on this elements
625 if (eval_p[i].val.cwiseAbs().maxCoeff() <= 1e-10)
626 continue;
627
628 for (std::size_t k = 0; k < other_b.global().size(); ++k)
629 {
630 // auto glob0 = other_b.global()[k]; glob0.val *= eval_p(0,i);
631 // auto glob1 = other_b.global()[k]; glob1.val *= eval_p(1,i);
632 // auto glob2 = other_b.global()[k]; glob2.val *= eval_p(2,i);
633
634 auto glob0 = other_b.global()[k];
635 glob0.val *= eval_p[i].val(0);
636 auto glob1 = other_b.global()[k];
637 glob1.val *= eval_p[i].val(1);
638 auto glob2 = other_b.global()[k];
639 glob2.val *= eval_p[i].val(2);
640
641 insert_into_global(glob0, b.bases[i0].global());
642 insert_into_global(glob1, b.bases[i1].global());
643 insert_into_global(glob2, b.bases[i2].global());
644 }
645 }
646
647 index = mesh.next_around_face(index);
648 }
649 }
650
651 void setup_data_for_polygons(const Mesh2D &mesh, const int el_index, const ElementBases &b, std::map<int, InterfaceData> &poly_edge_to_data)
652 {
653 Navigation::Index index = mesh.get_index_from_face(el_index);
654 for (int j = 0; j < 4; ++j)
655 {
656 const int opposite_face = mesh.switch_face(index).face;
657 const bool is_neigh_poly = opposite_face >= 0 && mesh.is_polytope(opposite_face);
658
659 if (is_neigh_poly)
660 {
661 // auto e2l = LagrangeBasis2d::quadr_quad_edge_local_nodes(mesh, index);
662 auto e2l = LagrangeBasis2d::quad_edge_local_nodes(2, mesh, index);
663 const int vertex_basis_id = e2l[0];
664 const int edge_basis_id = e2l[1];
665 const int vertex_basis_id2 = e2l[2];
666
667 InterfaceData &data = poly_edge_to_data[index.edge];
668
669 data.local_indices.push_back(edge_basis_id);
670 data.local_indices.push_back(vertex_basis_id);
671 data.local_indices.push_back(vertex_basis_id2);
672 }
673
674 index = mesh.next_around_face(index);
675 }
676 }
677 } // namespace
678
679 int SplineBasis2d::build_bases(const Mesh2D &mesh,
680 const std::string &assembler,
681 const int quadrature_order, const int mass_quadrature_order, std::vector<ElementBases> &bases, std::vector<LocalBoundary> &local_boundary, std::map<int, InterfaceData> &poly_edge_to_data)
682 {
683 using std::max;
684 assert(!mesh.is_volume());
685
686 MeshNodes mesh_nodes(mesh, true, true, 1, 1);
687
688 const int n_els = mesh.n_elements();
689 bases.resize(n_els);
690
691 local_boundary.clear();
692
693 // QuadQuadrature quad_quadrature;
694
695 for (int e = 0; e < n_els; ++e)
696 {
697 if (!mesh.is_spline_compatible(e))
698 continue;
699
700 SpaceMatrix space;
701 NodeMatrix loc_nodes;
702
703 // const int max_local_base =
704 build_local_space(mesh, mesh_nodes, e, space, loc_nodes, local_boundary, poly_edge_to_data);
705 // n_bases = max(n_bases, max_local_base);
706
707 ElementBases &b = bases[e];
708 // quad_quadrature.get_quadrature(quadrature_order, b.quadrature);
709 const int real_order = quadrature_order > 0 ? quadrature_order : AssemblerUtils::quadrature_order(assembler, 2, AssemblerUtils::BasisType::SPLINE, 2);
710 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order : AssemblerUtils::quadrature_order("Mass", 2, AssemblerUtils::BasisType::SPLINE, 2);
711
712 b.set_quadrature([real_order](Quadrature &quad) {
713 QuadQuadrature quad_quadrature;
714 quad_quadrature.get_quadrature(real_order, quad);
715 });
716 b.set_mass_quadrature([real_mass_order](Quadrature &quad) {
717 QuadQuadrature quad_quadrature;
718 quad_quadrature.get_quadrature(real_mass_order, quad);
719 });
720 b.bases.resize(9);
721
722 b.set_local_node_from_primitive_func([e](const int primitive_id, const Mesh &mesh) {
723 Eigen::VectorXi res(3);
724 const auto &mesh2d = dynamic_cast<const Mesh2D &>(mesh);
725 auto index = mesh2d.get_index_from_face(e);
726 int le;
727 for (le = 0; le < mesh2d.n_face_vertices(e); ++le)
728 {
729 if (index.edge == primitive_id)
730 break;
731 index = mesh2d.next_around_face(index);
732 }
733 assert(index.edge == primitive_id);
734
735 switch (le)
736 {
737 case 3:
738 res << (3 * 0 + 0), (3 * 1 + 0), (3 * 2 + 0);
739 break;
740 case 0:
741 res << (3 * 0 + 0), (3 * 0 + 1), (3 * 0 + 2);
742 break;
743 case 1:
744 res << (3 * 0 + 2), (3 * 1 + 2), (3 * 2 + 2);
745 break;
746 case 2:
747 res << (3 * 2 + 0), (3 * 2 + 1), (3 * 2 + 2);
748 break;
749 default:
750 assert(false);
751 }
752
753 return res;
754 });
755
756 std::array<std::array<double, 4>, 3> h_knots;
757 std::array<std::array<double, 4>, 3> v_knots;
758
759 setup_knots_vectors(mesh_nodes, space, h_knots, v_knots);
760
761 // print_local_space(space);
762
763 basis_for_regular_quad(space, loc_nodes, h_knots, v_knots, b);
764 basis_for_irregulard_quad(e, mesh, mesh_nodes, space, loc_nodes, h_knots, v_knots, b);
765 }
766
767 std::set<int> edge_id;
768 std::set<int> vertex_id;
769
770 int n_bases = mesh_nodes.n_nodes();
771
772 for (int e = 0; e < n_els; ++e)
773 {
774 if (mesh.is_polytope(e) || mesh.is_spline_compatible(e))
775 continue;
776
777 ElementBases &b = bases[e];
778 // quad_quadrature.get_quadrature(quadrature_order, b.quadrature);
779 const int real_order = quadrature_order > 0 ? quadrature_order : AssemblerUtils::quadrature_order(assembler, 2, AssemblerUtils::BasisType::CUBE_LAGRANGE, 2);
780 const int real_mass_order = mass_quadrature_order > 0 ? mass_quadrature_order : AssemblerUtils::quadrature_order("Mass", 2, AssemblerUtils::BasisType::CUBE_LAGRANGE, 2);
781
782 b.set_quadrature([real_order](Quadrature &quad) {
783 QuadQuadrature quad_quadrature;
784 quad_quadrature.get_quadrature(real_order, quad);
785 });
786 b.set_mass_quadrature([real_mass_order](Quadrature &quad) {
787 QuadQuadrature quad_quadrature;
788 quad_quadrature.get_quadrature(real_mass_order, quad);
789 });
790
791 b.set_local_node_from_primitive_func([e](const int primitive_id, const Mesh &mesh) {
792 const auto &mesh2d = dynamic_cast<const Mesh2D &>(mesh);
793 auto index = mesh2d.get_index_from_face(e);
794
795 for (int le = 0; le < mesh2d.n_face_vertices(e); ++le)
796 {
797 if (index.edge == primitive_id)
798 break;
799 index = mesh2d.next_around_face(index);
800 }
801 assert(index.edge == primitive_id);
802
803 // const auto indices = LagrangeBasis2d::quadr_quad_edge_local_nodes(mesh2d, index);
804 const auto indices = LagrangeBasis2d::quad_edge_local_nodes(2, mesh2d, index);
805 Eigen::VectorXi res(indices.size());
806
807 for (size_t i = 0; i < indices.size(); ++i)
808 res(i) = indices[i];
809
810 return res;
811 });
812
813 create_q2_nodes(mesh, e, vertex_id, edge_id, b, local_boundary, n_bases);
814 }
815
816 bool missing_bases = false;
817 do
818 {
819 missing_bases = false;
820 for (int e = 0; e < n_els; ++e)
821 {
822 if (mesh.is_polytope(e) || mesh.is_spline_compatible(e))
823 continue;
824
825 auto &b = bases[e];
826 if (b.is_complete())
827 continue;
828
829 assign_q2_weights(mesh, e, bases);
830
831 missing_bases = missing_bases || b.is_complete();
832 }
833 } while (missing_bases);
834
835 for (int e = 0; e < n_els; ++e)
836 {
837 if (mesh.is_polytope(e) || mesh.is_spline_compatible(e))
838 continue;
839 const ElementBases &b = bases[e];
840 setup_data_for_polygons(mesh, e, b, poly_edge_to_data);
841 }
842
843 return n_bases;
844 }
845
846 void SplineBasis2d::fit_nodes(const Mesh2D &mesh, const int n_bases, std::vector<ElementBases> &gbases)
847 {
848 assert(false);
849 // const int dim = 2;
850 // const int n_constraints = 9;
851 // const int n_elements = mesh.n_elements();
852
853 // std::vector< Eigen::Triplet<double> > entries, entries_t;
854
855 // MeshNodes nodes(mesh, 1, 1, 0);
856 // // Eigen::MatrixXd tmp;
857 // std::vector<AssemblyValues> tmp_val;
858
859 // Eigen::MatrixXd node_rhs(n_constraints*n_elements, dim);
860 // Eigen::MatrixXd samples(n_constraints, dim);
861 // polyfem::autogen::q_nodes_2d(2, samples);
862 // // for(int i = 0; i < n_constraints; ++i)
863 // // samples.row(i) = LagrangeBasis2d::quadr_quad_local_node_coordinates(i);
864
865 // for(int i = 0; i < n_elements; ++i)
866 // {
867 // auto &base = gbases[i];
868
869 // if(!mesh.is_cube(i))
870 // continue;
871
872 // auto global_ids = LagrangeBasis2d::quadr_quad_local_to_global(mesh, i);
873 // assert(global_ids.size() == n_constraints);
874
875 // for(int j = 0; j < n_constraints; ++j)
876 // {
877 // auto n_id = nodes.node_id_from_primitive(global_ids[j]);
878 // auto n = nodes.node_position(n_id);
879 // for(int d = 0; d < dim; ++d)
880 // node_rhs(n_constraints*i + j, d) = n(d);
881 // }
882
883 // base.evaluate_bases(samples, tmp_val);
884 // const auto &lbs = base.bases;
885
886 // const int n_local_bases = int(lbs.size());
887 // for(int j = 0; j < n_local_bases; ++j)
888 // {
889 // const Basis &b = lbs[j];
890 // const auto &tmp = tmp_val[j].val;
891
892 // for(std::size_t ii = 0; ii < b.global().size(); ++ii)
893 // {
894 // for (long k = 0; k < tmp.size(); ++k)
895 // {
896 // entries.emplace_back(n_constraints*i + k, b.global()[ii].index, tmp(k)*b.global()[ii].val);
897 // entries_t.emplace_back(b.global()[ii].index, n_constraints*i + k, tmp(k)*b.global()[ii].val);
898 // }
899 // }
900 // }
901 // }
902
903 // Eigen::MatrixXd new_nodes(n_bases, dim);
904
905 // {
906 // StiffnessMatrix mat(n_constraints*n_elements, n_bases);
907 // StiffnessMatrix mat_t(n_bases, n_constraints*n_elements);
908
909 // mat.setFromTriplets(entries.begin(), entries.end());
910 // mat_t.setFromTriplets(entries_t.begin(), entries_t.end());
911
912 // StiffnessMatrix A = mat_t * mat;
913 // Eigen::MatrixXd b = mat_t * node_rhs;
914
915 // json params = {
916 // {"mtype", -2}, // matrix type for Pardiso (2 = SPD)
917 // // {"max_iter", 0}, // for iterative solvers
918 // // {"tolerance", 1e-9}, // for iterative solvers
919 // };
920 // auto solver = LinearSolver::create("", "");
921 // solver->setParameters(params);
922 // solver->analyzePattern(A);
923 // solver->factorize(A);
924
925 // for(int d = 0; d < dim; ++d)
926 // solver->solve(b.col(d), new_nodes.col(d));
927 // }
928
929 // for(int i = 0; i < n_elements; ++i)
930 // {
931 // auto &base = gbases[i];
932
933 // if(!mesh.is_cube(i))
934 // continue;
935
936 // auto &lbs = base.bases;
937 // const int n_local_bases = int(lbs.size());
938 // for(int j = 0; j < n_local_bases; ++j)
939 // {
940 // Basis &b = lbs[j];
941
942 // for(std::size_t ii = 0; ii < b.global().size(); ++ii)
943 // {
944 // // if(nodes.is_primitive_boundary(b.global()[ii].index))
945 // // continue;
946
947 // for(int d = 0; d < dim; ++d)
948 // {
949 // b.global()[ii].node(d) = new_nodes(b.global()[ii].index, d);
950 // }
951 // }
952 // }
953 // }
954 }
955 } // namespace basis
956} // namespace polyfem
vec
Eigen::MatrixXd vec
Definition Assembler.cpp:72
val
double val
Definition Assembler.cpp:86
quadrature
Quadrature quadrature
Definition MassMatrixAssembler.cpp:30
y
int y
Definition SplineBasis3d.cpp:55
edge_id
int edge_id
Definition SplineBasis3d.cpp:56
x
int x
Definition SplineBasis3d.cpp:55
auto_q_bases.hpp
polyfem::assembler::AssemblerUtils::quadrature_order
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...
Definition AssemblerUtils.cpp:189
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::set_quadrature
void set_quadrature(const QuadratureFunction &fun)
Definition ElementBases.hpp:66
polyfem::basis::ElementBases::set_local_node_from_primitive_func
void set_local_node_from_primitive_func(LocalNodeFromPrimitiveFunc fun)
sets mapping from local nodes to global nodes
Definition ElementBases.hpp:100
polyfem::basis::ElementBases::bases
std::vector< Basis > bases
one basis function per node in the element
Definition ElementBases.hpp:27
polyfem::basis::ElementBases::set_mass_quadrature
void set_mass_quadrature(const QuadratureFunction &fun)
Definition ElementBases.hpp:67
polyfem::basis::LagrangeBasis2d::quad_edge_local_nodes
static Eigen::VectorXi quad_edge_local_nodes(const int q, const mesh::Mesh2D &mesh, mesh::Navigation::Index index)
Definition LagrangeBasis2d.cpp:642
polyfem::basis::SplineBasis2d::fit_nodes
static void fit_nodes(const mesh::Mesh2D &mesh, const int n_bases, std::vector< ElementBases > &gbases)
Definition SplineBasis2d.cpp:846
polyfem::basis::SplineBasis2d::build_bases
static int build_bases(const mesh::Mesh2D &mesh, const std::string &assembler, const int quadrature_order, const int mass_quadrature_order, std::vector< ElementBases > &bases, std::vector< mesh::LocalBoundary > &local_boundary, std::map< int, InterfaceData > &poly_edge_to_data)
Definition SplineBasis2d.cpp:679
polyfem::mesh::Mesh2D::is_volume
bool is_volume() const override
checks if mesh is volume
Definition Mesh2D.hpp:25
polyfem::mesh::Mesh2D::get_index_from_face
virtual Navigation::Index get_index_from_face(int f, int lv=0) const =0
polyfem::mesh::Mesh
Abstract mesh class to capture 2d/3d conforming and non-conforming meshes.
Definition Mesh.hpp:39
polyfem::mesh::Mesh::n_elements
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
Definition Mesh.hpp:161
polyfem::mesh::Mesh::is_polytope
bool is_polytope(const int el_id) const
checks if element is polygon compatible
Definition Mesh.cpp:365
polyfem::mesh::Mesh::is_spline_compatible
bool is_spline_compatible(const int el_id) const
checks if element is spline compatible
Definition Mesh.cpp:332
polyfem::quadrature::QuadQuadrature::get_quadrature
void get_quadrature(const int order, Quadrature &quad)
Definition QuadQuadrature.cpp:16
p_bases.base
str base
Definition p_bases.py:428
p_bases.indices
list indices
Definition p_bases.py:258
polyfem::autogen::q_nodes_2d
void q_nodes_2d(const int q, Eigen::MatrixXd &val)
Definition auto_q_bases_2d_nodes.cpp:67
polyfem::autogen::q_grad_basis_value_2d
void q_grad_basis_value_2d(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_q_bases_2d_grad.cpp:182
polyfem::autogen::q_basis_value_2d
void q_basis_value_2d(const int q, const int local_index, const Eigen::MatrixXd &uv, Eigen::MatrixXd &val)
Definition auto_q_bases_2d_val.cpp:88
polyfem::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:44