18#include <igl/per_vertex_normals.h>
25 using namespace assembler;
28 using namespace utils;
35 const int max_num_kernels = 300;
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)
43 std::vector<int> local_to_global;
45 for (
int lf = 0; lf < mesh.n_cell_faces(c); ++lf)
47 auto index = mesh.get_index_from_element(c, lf, 0);
48 const int c2 = mesh.switch_element(index).element;
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)
55 for (
const auto &
x : b.bases[other_local_basis_id].global())
57 const int global_node_id =
x.index;
58 local_to_global.push_back(global_node_id);
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));
67 return local_to_global;
73 void compute_canonical_pattern(
int n_samples_per_edge, Eigen::MatrixXd &
V, Eigen::MatrixXi &
F)
85 constexpr int lv0 = 3;
89 const Mesh3D &mesh,
int c, Eigen::MatrixXd &
V, Eigen::MatrixXi &
F)
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));
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)
99 face_g2l.emplace(mesh.get_index_from_element(c, lf, lv0).face, lf);
101 for (
int lf = 0; lf < mesh.n_cell_faces(c); ++lf)
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)
107 if (!vertex_g2l.count(index.vertex))
109 vertex_g2l.emplace(index.vertex, num_vertices++);
111 quads[lf][lv] = vertex_g2l.at(index.vertex);
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)
120 if (index3.edge == index2.edge)
122 std::get<1>(
adj[lf][lv]) = lv2;
123 if (index2.vertex != index3.vertex)
125 assert(mesh.switch_vertex(index3).vertex == index2.vertex);
126 std::get<2>(
adj[lf][lv]) =
true;
130 std::get<2>(
adj[lf][lv]) =
false;
133 index3 = mesh.next_around_face(index3);
136 index = mesh.next_around_face(index);
139 V.resize(num_vertices, 3);
140 for (
const auto &kv : vertex_g2l)
142 V.row(kv.second) = mesh.point(kv.first);
144 F.resize(quads.size(), 4);
148 F.row(f++) << q[0], q[1], q[2], q[3];
151 return [
adj](
int q,
int lv) {
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,
163 Eigen::MatrixXd PV, KN;
166 compute_canonical_pattern(n_kernels_per_edge, PV, PF);
171 if (
true || KV.rows() < max_num_kernels)
173 igl::per_vertex_normals(KV, KF, KN);
180 kernel_centers += eps * volume * KN;
187 assert(kernel_centers.cols() == 3);
189 std::vector<Eigen::RowVector3d> remap;
190 std::vector<Eigen::RowVector3d> rejected;
191 for (
int v = 0; v < kernel_centers.rows(); ++v)
193 if (std::sqrt(D(v)) > 0.8 * eps * volume)
195 remap.push_back(kernel_centers.row(v));
202 kernel_centers.resize(remap.size(), 3);
203 for (
int v = 0; v < kernel_centers.rows(); ++v)
205 kernel_centers.row(v) = remap[v];
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,
220 const std::map<int, InterfaceData> &poly_face_to_data,
221 const std::vector<ElementBases> &bases,
222 const std::vector<ElementBases> &gbases,
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,
233 Eigen::RowVector3d &translation)
236 local_to_global = compute_nonzero_bases_ids(mesh, element_index, bases, poly_face_to_data);
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);
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);
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);
268 Eigen::MatrixXd QV, KV;
269 Eigen::MatrixXi QF, KF;
270 auto getAdjLocalEdge = compute_quad_mesh_from_cell(mesh, element_index, QV, QF);
273 compute_offset_kernels(QV, QF, n_kernels_per_edge, eps, kernel_centers, KV, KF,
274 evalFuncGeom, getAdjLocalEdge);
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);
284 instantiate_pattern(QV, QF, PV, PF, collocation_points, CF,
nullptr, evalFuncGeom, getAdjLocalEdge);
286 reorder_mesh(collocation_points, CF, uv_sources, uv_ranges);
288 assert(uv_ranges.size() == mesh.n_cell_faces(element_index) + 1);
291 compute_canonical_pattern(n_quadrature_vertices_per_edge, PV, PF);
293 nullptr, evalFuncGeom, getAdjLocalEdge);
323 Eigen::MatrixXd samples;
324 std::vector<AssemblyValues> basis_val;
325 rhs.resize(UV.rows(), local_to_global.size());
327 for (
int lf = 0; lf < mesh.n_cell_faces(element_index); ++lf)
329 auto index = mesh.get_index_from_element(element_index, lf, 0);
330 const int c2 = mesh.switch_element(index).element;
333 const InterfaceData &bdata = poly_face_to_data.at(index.face);
334 const ElementBases &b = bases[
c2];
336 samples = UV.middleRows(uv_ranges(lf), uv_ranges(lf + 1) - uv_ranges(lf));
337 b.evaluate_bases(samples, basis_val);
340 for (
int other_local_basis_id : bdata.local_indices)
344 for (
const auto &
x : b.bases[other_local_basis_id].global())
346 const int global_node_id =
x.index;
347 const double weight =
x.val;
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;
357 Eigen::MatrixXd NV = triangulated_vertices;
361 translation.setZero();
364 quadrature_order, quadrature);
366 mass_quadrature_order, mass_quadrature);
373 triangulated_vertices = KV;
374 triangulated_faces = KF;
386 const std::vector<ElementBases> &bases,
387 const std::vector<ElementBases> &gbases,
388 Eigen::MatrixXd &basis_integrals)
392 basis_integrals.resize(n_bases, 9);
393 basis_integrals.setZero();
394 Eigen::MatrixXd rhs(n_bases, 9);
398 for (
int e = 0; e < n_elements; ++e)
410 const int n_local_bases = int(
vals.basis_values.size());
411 for (
int j = 0; j < n_local_bases; ++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();
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();
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();
426 const double area = (v.
val.array() *
vals.det.array() *
vals.quadrature.weights.array()).sum();
428 for (
size_t ii = 0; ii < v.
global.size(); ++ii)
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;
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;
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;
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;
449 basis_integrals -= rhs;
480 const int nn_samples_per_edge,
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)
492 if (poly_face_to_data.empty())
496 int n_kernels_per_edge = 4;
497 int n_samples_per_edge = 3 * n_kernels_per_edge;
500 Eigen::MatrixXd basis_integrals;
516 std::vector<int> local_to_global;
517 Eigen::MatrixXd collocation_points, kernel_centers, triangulated_vertices;
518 Eigen::MatrixXi triangulated_faces;
524 Quadrature tmp_quadrature, tmp_mass_quadrature;
526 Eigen::RowVector3d translation;
527 sample_polyhedra(e, 2, n_kernels_per_edge, n_samples_per_edge,
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);
571 Eigen::MatrixXd local_basis_integrals(rhs.cols(), basis_integrals.cols());
572 for (
long k = 0; k < rhs.cols(); ++k)
574 local_basis_integrals.row(k) = -basis_integrals.row(local_to_global[k]);
576 auto set_rbf = [&b](
auto rbf) {
577 b.
set_bases_func([rbf](
const Eigen::MatrixXd &uv, std::vector<AssemblyValues> &
val) {
579 rbf->bases_values(uv, tmp);
580 val.resize(tmp.cols());
581 assert(tmp.rows() == uv.rows());
583 for (
size_t i = 0; i < tmp.cols(); ++i)
585 val[i].val = tmp.col(i);
588 b.
set_grads_func([rbf](
const Eigen::MatrixXd &uv, std::vector<AssemblyValues> &
val) {
589 Eigen::MatrixXd tmpx, tmpy, tmpz;
591 rbf->bases_grads(0, uv, tmpx);
592 rbf->bases_grads(1, uv, tmpy);
593 rbf->bases_grads(2, uv, tmpz);
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)
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);
608 if (integral_constraints == 0)
610 set_rbf(std::make_shared<RBFWithLinear>(
611 kernel_centers, collocation_points, local_basis_integrals, tmp_quadrature, rhs,
false));
613 else if (integral_constraints == 1)
615 set_rbf(std::make_shared<RBFWithLinear>(
616 kernel_centers, collocation_points, local_basis_integrals, tmp_quadrature, rhs));
618 else if (integral_constraints == 2)
620 set_rbf(std::make_shared<RBFWithQuadratic>(
622 assembler, kernel_centers, collocation_points, local_basis_integrals, tmp_quadrature, rhs));
626 throw std::runtime_error(fmt::format(
"Unsupported constraint order: {:d}", integral_constraints));
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)
634 b.
bases[i].init(-2, local_to_global[i], i, Eigen::MatrixXd::Constant(1, 3, std::nan(
"")));
639 mapped_boundary[e].first = triangulated_vertices;
640 mapped_boundary[e].second = triangulated_faces;
vector< list< int > > adj
ElementAssemblyValues vals
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
bool has_parameterization
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
int n_elements() const
utitlity to return the number of elements, cells or faces in 3d and 2d
bool is_polytope(const int el_id) const
checks if element is polygon compatible
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.
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.