Wildmeshing Toolkit
Loading...
Searching...
No Matches
EnvelopeInvariant.cpp
Go to the documentation of this file.
1#if defined(__GNUG__) && !defined(__clang__)
2#pragma GCC diagnostic push
3#pragma GCC diagnostic ignored "-Wnull-dereference"
4#endif
5#include <Eigen/Dense>
6#if defined(__GNUG__) && !defined(__clang__)
7#pragma GCC diagnostic pop
8#endif
9#include "EnvelopeInvariant.hpp"
10
11#include <wmtk/Mesh.hpp>
12#include <wmtk/simplex/SimplexCollection.hpp>
13#include <wmtk/simplex/faces_single_dimension.hpp>
14#include <wmtk/simplex/utils/tuple_vector_to_homogeneous_simplex_vector.hpp>
15#include <wmtk/utils/Logger.hpp>
16
17
18#include <fastenvelope/FastEnvelope.h>
19#include <SimpleBVH/BVH.hpp>
20
21
22namespace wmtk {
23
24constexpr PrimitiveType PV = PrimitiveType::Vertex;
25constexpr PrimitiveType PE = PrimitiveType::Edge;
26} // namespace wmtk
27namespace wmtk::invariants {
28
29EnvelopeInvariant::EnvelopeInvariant(
30 const attribute::MeshAttributeHandle& envelope_mesh_coordinate,
31 double envelope_size,
32 const attribute::MeshAttributeHandle& coordinate)
33 : Invariant(coordinate.mesh())
34 , m_coordinate_handle(coordinate)
35 , m_envelope_size(envelope_size)
36{
37 const auto& envelope_mesh = envelope_mesh_coordinate.mesh();
38
39 assert(envelope_mesh_coordinate.holds<Rational>() || envelope_mesh_coordinate.holds<double>());
40
41 if (envelope_mesh_coordinate.holds<Rational>()) {
42 // for rational
43 const attribute::Accessor<Rational> accessor =
44 envelope_mesh.create_const_accessor(envelope_mesh_coordinate.as<Rational>());
45
46
47 if (envelope_mesh.top_simplex_type() == PrimitiveType::Triangle) {
48 std::vector<Eigen::Vector3d> vertices;
49 std::vector<Eigen::Vector3i> faces;
50
51 int count = 0;
52 assert(accessor.dimension() == 3);
53
54 const std::vector<Tuple>& facest = envelope_mesh.get_all(wmtk::PrimitiveType::Triangle);
55 for (const auto& f : facest) {
56 Eigen::Vector3d p0 = accessor.const_vector_attribute(f).cast<double>();
57 Eigen::Vector3d p1 =
58 accessor.const_vector_attribute(envelope_mesh.switch_tuple(f, PV))
59 .cast<double>();
60 Eigen::Vector3d p2 =
61 accessor.const_vector_attribute(envelope_mesh.switch_tuples(f, {PE, PV}))
62 .cast<double>();
63
64 faces.emplace_back(count, count + 1, count + 2);
65 vertices.push_back(p0);
66 vertices.push_back(p1);
67 vertices.push_back(p2);
68
69 count += 3;
70 }
71
72 m_envelope =
73 std::make_shared<fastEnvelope::FastEnvelope>(vertices, faces, envelope_size);
74
75 } else if (envelope_mesh.top_simplex_type() == PrimitiveType::Edge) {
76 logger().warn("Envelope for edge mesh is using sampling");
77
78 int64_t count = 0;
79 int64_t index = 0;
80
81 const std::vector<Tuple>& edgest = envelope_mesh.get_all(wmtk::PrimitiveType::Edge);
82
83 Eigen::MatrixXd vertices(2 * edgest.size(), accessor.dimension());
84 Eigen::MatrixXi edges(edgest.size(), 2);
85
86 for (const auto& e : edgest) {
87 auto p0 = accessor.const_vector_attribute(e).cast<double>();
88 auto p1 = accessor.const_vector_attribute(envelope_mesh.switch_tuple(e, PV))
89 .cast<double>();
90
91 edges.row(index) << count, count + 1;
92 vertices.row(2 * index) = p0;
93 vertices.row(2 * index + 1) = p1;
94
95 count += 2;
96 ++index;
97 }
98
99 m_bvh = std::make_shared<SimpleBVH::BVH>();
100 m_bvh->init(vertices, edges, 1e-10);
101 } else {
102 throw std::runtime_error("Envelope works only for tri/edges meshes");
103 }
104 } else if (envelope_mesh_coordinate.holds<double>()) {
105 // for double
106 const attribute::Accessor<double> accessor =
107 envelope_mesh.create_const_accessor(envelope_mesh_coordinate.as<double>());
108
109
110 if (envelope_mesh.top_simplex_type() == PrimitiveType::Triangle) {
111 std::vector<Eigen::Vector3d> vertices;
112 std::vector<Eigen::Vector3i> faces;
113
114 int count = 0;
115 assert(accessor.dimension() == 3);
116
117 const std::vector<Tuple>& facest = envelope_mesh.get_all(wmtk::PrimitiveType::Triangle);
118 for (const auto& f : facest) {
119 Eigen::Vector3d p0 = accessor.const_vector_attribute(f);
120 Eigen::Vector3d p1 =
121 accessor.const_vector_attribute(envelope_mesh.switch_tuple(f, PV));
122 Eigen::Vector3d p2 =
123 accessor.const_vector_attribute(envelope_mesh.switch_tuples(f, {PE, PV}));
124
125 faces.emplace_back(count, count + 1, count + 2);
126 vertices.push_back(p0);
127 vertices.push_back(p1);
128 vertices.push_back(p2);
129
130 count += 3;
131 }
132
133 m_envelope =
134 std::make_shared<fastEnvelope::FastEnvelope>(vertices, faces, envelope_size);
135
136 } else if (envelope_mesh.top_simplex_type() == PrimitiveType::Edge) {
137 logger().warn("Envelope for edge mesh is using sampling");
138
139 int64_t count = 0;
140 int64_t index = 0;
141
142 const std::vector<Tuple>& edgest = envelope_mesh.get_all(wmtk::PrimitiveType::Edge);
143
144 Eigen::MatrixXd vertices(2 * edgest.size(), accessor.dimension());
145 Eigen::MatrixXi edges(edgest.size(), 2);
146
147 for (const auto& e : edgest) {
148 auto p0 = accessor.const_vector_attribute(e);
149 auto p1 = accessor.const_vector_attribute(envelope_mesh.switch_tuple(e, PV));
150
151 edges.row(index) << count, count + 1;
152 vertices.row(2 * index) = p0;
153 vertices.row(2 * index + 1) = p1;
154
155 count += 2;
156 ++index;
157 }
158
159 m_bvh = std::make_shared<SimpleBVH::BVH>();
160 m_bvh->init(vertices, edges, 1e-10);
161 } else {
162 throw std::runtime_error("Envelope works only for tri/edges meshes");
163 }
164 } else {
165 throw std::runtime_error("Envelope mesh handle type invlid");
166 }
167}
168
169bool EnvelopeInvariant::after(
170 const std::vector<Tuple>& top_dimension_tuples_before,
171 const std::vector<Tuple>& top_dimension_tuples_after) const
172{
173 if (top_dimension_tuples_after.empty()) return true;
174
175 assert(m_coordinate_handle.holds<Rational>() || m_coordinate_handle.holds<double>());
176
177 if (m_coordinate_handle.holds<Rational>()) {
178 const attribute::Accessor<Rational> accessor =
179 mesh().create_const_accessor(m_coordinate_handle.as<Rational>());
180 const PrimitiveType type = mesh().top_simplex_type();
181
182 if (m_envelope) {
183 assert(accessor.dimension() == 3);
184
185 std::vector<Tuple> faces;
186
187 if (type == PrimitiveType::Triangle) {
188 std::array<Eigen::Vector3d, 3> triangle;
189
190 for (const Tuple& tuple : top_dimension_tuples_after) {
191 faces = faces_single_dimension_tuples(
192 mesh(),
193 simplex::Simplex(mesh(), type, tuple),
194 PrimitiveType::Vertex);
195
196 triangle[0] = accessor.const_vector_attribute(faces[0]).cast<double>();
197 triangle[1] = accessor.const_vector_attribute(faces[1]).cast<double>();
198 triangle[2] = accessor.const_vector_attribute(faces[2]).cast<double>();
199
200 if (m_envelope->is_outside(triangle)) {
201 wmtk::logger().debug("fail envelope check 1");
202 return false;
203 }
204 }
205
206 return true;
207 } else if (type == PrimitiveType::Edge) {
208 for (const Tuple& tuple : top_dimension_tuples_after) {
209 faces = faces_single_dimension_tuples(
210 mesh(),
211 simplex::Simplex(mesh(), type, tuple),
212 PrimitiveType::Vertex);
213
214 Eigen::Vector3d p0 = accessor.const_vector_attribute(faces[0]).cast<double>();
215 Eigen::Vector3d p1 = accessor.const_vector_attribute(faces[1]).cast<double>();
216
217 if (m_envelope->is_outside(p0, p1)) {
218 wmtk::logger().debug("fail envelope check 2");
219 return false;
220 }
221 }
222
223 return true;
224 } else if (type == PrimitiveType::Vertex) {
225 for (const Tuple& tuple : top_dimension_tuples_after) {
226 Eigen::Vector3d p = accessor.const_vector_attribute(tuple).cast<double>();
227
228 if (m_envelope->is_outside(p)) {
229 wmtk::logger().debug("fail envelope check 3");
230 return false;
231 }
232 }
233
234 return true;
235 } else {
236 throw std::runtime_error("Invalid mesh type");
237 }
238 return true;
239 } else {
240 assert(m_bvh);
241
242 SimpleBVH::VectorMax3d nearest_point;
243 double sq_dist;
244
245 const double d = m_envelope_size;
246 const double real_envelope = m_envelope_size - d / sqrt(accessor.dimension());
247 const double real_envelope_2 = real_envelope * real_envelope;
248
249 if (type == PrimitiveType::Edge) {
250 std::vector<SimpleBVH::VectorMax3d> pts;
251
252 for (const Tuple& tuple : top_dimension_tuples_after) {
253 const auto p0 = accessor.const_vector_attribute(tuple).cast<double>().eval();
254 const auto p1 = accessor.const_vector_attribute(mesh().switch_tuple(tuple, PV))
255 .cast<double>()
256 .eval();
257
258 const int64_t N = (p0 - p1).norm() / d + 1;
259 pts.reserve(pts.size() + N);
260
261 for (int64_t n = 0; n <= N; n++) {
262 auto tmp = p0 * (double(n) / N) + p1 * (N - double(n)) / N;
263 pts.push_back(tmp);
264 }
265 }
266
267 auto current_point = pts[0];
268
269 int prev_facet = m_bvh->nearest_facet(current_point, nearest_point, sq_dist);
270 if (sq_dist > real_envelope_2) {
271 wmtk::logger().debug("fail envelope check 4");
272 return false;
273 }
274
275 for (const auto& v : pts) {
276 sq_dist = (v - nearest_point).squaredNorm();
277 m_bvh->nearest_facet_with_hint(v, prev_facet, nearest_point, sq_dist);
278 if (sq_dist > real_envelope_2) {
279 wmtk::logger().debug("fail envelope check 5");
280 return false;
281 }
282 }
283
284 return true;
285 } else if (type == PrimitiveType::Vertex) {
286 for (const Tuple& tuple : top_dimension_tuples_after) {
287 Eigen::Vector3d p = accessor.const_vector_attribute(tuple).cast<double>();
288 m_bvh->nearest_facet(p, nearest_point, sq_dist);
289 if (sq_dist > m_envelope_size * m_envelope_size) {
290 wmtk::logger().debug("fail envelope check 6");
291 return false;
292 }
293 }
294
295 return true;
296 } else {
297 throw std::runtime_error("Invalid mesh type");
298 }
299 return true;
300 }
301 } else if (m_coordinate_handle.holds<double>()) {
302 const attribute::Accessor<double> accessor =
303 mesh().create_const_accessor(m_coordinate_handle.as<double>());
304 const auto type = mesh().top_simplex_type();
305
306 if (m_envelope) {
307 assert(accessor.dimension() == 3);
308
309 std::vector<Tuple> faces;
310
311 if (type == PrimitiveType::Triangle) {
312 std::array<Eigen::Vector3d, 3> triangle;
313
314 for (const Tuple& tuple : top_dimension_tuples_after) {
315 faces = faces_single_dimension_tuples(
316 mesh(),
317 simplex::Simplex(mesh(), type, tuple),
318 PrimitiveType::Vertex);
319
320 triangle[0] = accessor.const_vector_attribute(faces[0]);
321 triangle[1] = accessor.const_vector_attribute(faces[1]);
322 triangle[2] = accessor.const_vector_attribute(faces[2]);
323
324 if (m_envelope->is_outside(triangle)) {
325 wmtk::logger().debug("fail envelope check 7");
326 return false;
327 }
328 }
329
330 return true;
331 } else if (type == PrimitiveType::Edge) {
332 for (const Tuple& tuple : top_dimension_tuples_after) {
333 faces = faces_single_dimension_tuples(
334 mesh(),
335 simplex::Simplex(mesh(), type, tuple),
336 PrimitiveType::Vertex);
337
338 Eigen::Vector3d p0 = accessor.const_vector_attribute(faces[0]);
339 Eigen::Vector3d p1 = accessor.const_vector_attribute(faces[1]);
340
341 if (m_envelope->is_outside(p0, p1)) {
342 wmtk::logger().debug("fail envelope check 8");
343 return false;
344 }
345 }
346
347 return true;
348 } else if (type == PrimitiveType::Vertex) {
349 for (const Tuple& tuple : top_dimension_tuples_after) {
350 Eigen::Vector3d p = accessor.const_vector_attribute(tuple);
351
352 if (m_envelope->is_outside(p)) {
353 wmtk::logger().debug("fail envelope check 9");
354 return false;
355 }
356 }
357
358 return true;
359 } else {
360 throw std::runtime_error("Invalid mesh type");
361 }
362 return true;
363 } else {
364 assert(m_bvh);
365
366 SimpleBVH::VectorMax3d nearest_point;
367 double sq_dist;
368
369 const double d = m_envelope_size;
370 const double real_envelope = m_envelope_size - d / sqrt(accessor.dimension());
371 const double real_envelope_2 = real_envelope * real_envelope;
372
373 if (type == PrimitiveType::Edge) {
374 std::vector<SimpleBVH::VectorMax3d> pts;
375
376 for (const Tuple& tuple : top_dimension_tuples_after) {
377 SimpleBVH::VectorMax3d p0 =
378 accessor.const_vector_attribute(tuple).cast<double>();
379 SimpleBVH::VectorMax3d p1 =
380 accessor.const_vector_attribute(mesh().switch_tuple(tuple, PV))
381 .cast<double>();
382
383 const int64_t N = (p0 - p1).norm() / d + 1;
384 pts.reserve(pts.size() + N);
385
386 for (int64_t n = 0; n <= N; n++) {
387 auto tmp = p0 * (double(n) / N) + p1 * (N - double(n)) / N;
388 pts.push_back(tmp);
389 }
390 }
391
392 auto current_point = pts[0];
393
394 int prev_facet = m_bvh->nearest_facet(current_point, nearest_point, sq_dist);
395 if (sq_dist > real_envelope_2) {
396 wmtk::logger().debug("fail envelope check 10");
397 return false;
398 }
399
400 for (const auto& v : pts) {
401 sq_dist = (v - nearest_point).squaredNorm();
402 m_bvh->nearest_facet_with_hint(v, prev_facet, nearest_point, sq_dist);
403 if (sq_dist > real_envelope_2) {
404 wmtk::logger().debug("fail envelope check 11");
405 return false;
406 }
407 }
408
409 return true;
410 } else if (type == PrimitiveType::Vertex) {
411 for (const Tuple& tuple : top_dimension_tuples_after) {
412 Eigen::Vector3d p = accessor.const_vector_attribute(tuple).cast<double>();
413 m_bvh->nearest_facet(p, nearest_point, sq_dist);
414 if (sq_dist > m_envelope_size * m_envelope_size) {
415 wmtk::logger().debug("fail envelope check 12");
416 return false;
417 }
418 }
419
420 return true;
421 } else {
422 throw std::runtime_error("Invalid mesh type");
423 }
424 return true;
425 }
426 } else {
427 throw std::runtime_error("Envelope mesh handle type invlid");
428 }
429}
430
431
432} // namespace wmtk::invariants
wmtk::Tuple
The Tuple is the basic navigation tool in our mesh data structure.
Definition Tuple.hpp:19
wmtk::attribute::Accessor
A CachingAccessor that uses tuples for accessing attributes instead of indices.
Definition Accessor.hpp:28
wmtk::attribute::Accessor::const_vector_attribute
ConstMapResult< D > const_vector_attribute(const ArgType &t) const
wmtk::attribute::MeshAttributeHandle::as
auto as() const -> const held_handle_type< held_type_from_primitive< T >()> &
Definition MeshAttributeHandle.hpp:161
wmtk::invariants::EnvelopeInvariant::m_envelope
std::shared_ptr< fastEnvelope::FastEnvelope > m_envelope
Definition EnvelopeInvariant.hpp:32
wmtk::invariants::EnvelopeInvariant::EnvelopeInvariant
EnvelopeInvariant(const attribute::MeshAttributeHandle &envelope_mesh_coordinate, double envelope_size, const attribute::MeshAttributeHandle &coordinate)
Definition EnvelopeInvariant.cpp:29
wmtk::invariants::EnvelopeInvariant::m_bvh
std::shared_ptr< SimpleBVH::BVH > m_bvh
Definition EnvelopeInvariant.hpp:33
faces_single_dimension.hpp
wmtk::PE
constexpr PrimitiveType PE
Definition EnvelopeInvariant.cpp:25
wmtk::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:58
wmtk::PV
constexpr PrimitiveType PV
Definition EnvelopeInvariant.cpp:24
tuple_vector_to_homogeneous_simplex_vector.hpp