Wildmeshing Toolkit
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 
22 namespace wmtk {
23 
24 constexpr PrimitiveType PV = PrimitiveType::Vertex;
25 constexpr PrimitiveType PE = PrimitiveType::Edge;
26 } // namespace wmtk
27 namespace wmtk::invariants {
28 
29 EnvelopeInvariant::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 
169 bool 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::Mesh::create_const_accessor
const attribute::Accessor< T, Mesh, D > create_const_accessor(const attribute::MeshAttributeHandle &handle) const
wmtk::Mesh::switch_tuple
virtual Tuple switch_tuple(const Tuple &tuple, PrimitiveType type) const =0
switch the orientation of the Tuple of the given dimension
wmtk::Mesh::top_simplex_type
PrimitiveType top_simplex_type() const
Definition: Mesh.hpp:996
wmtk::attribute::Accessor
A CachingAccessor that uses tuples for accessing attributes instead of indices.
Definition: Accessor.hpp:25
wmtk::attribute::Accessor::dimension
int64_t dimension() const
Definition: AccessorBase.hxx:57
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:160
wmtk::invariants::EnvelopeInvariant::m_envelope
std::shared_ptr< fastEnvelope::FastEnvelope > m_envelope
Definition: EnvelopeInvariant.hpp:32
wmtk::invariants::EnvelopeInvariant::after
bool after(const std::vector< Tuple > &top_dimension_tuples_before, const std::vector< Tuple > &top_dimension_tuples_after) const override
Definition: EnvelopeInvariant.cpp:169
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
wmtk::invariants::EnvelopeInvariant::m_coordinate_handle
const attribute::MeshAttributeHandle m_coordinate_handle
Definition: EnvelopeInvariant.hpp:34
wmtk::invariants::Invariant::mesh
const Mesh & mesh() const
Definition: Invariant.cpp:35
faces_single_dimension.hpp
wmtk::simplex::vertices
std::vector< Tuple > vertices(const Mesh &m, const Simplex &simplex)
Definition: faces_single_dimension.cpp:9
wmtk::simplex::edges
std::vector< Tuple > edges(const Mesh &m, const Simplex &simplex)
Definition: faces_single_dimension.cpp:77
wmtk::simplex::faces_single_dimension_tuples
std::vector< Tuple > faces_single_dimension_tuples(const Mesh &mesh, const Simplex &simplex, const PrimitiveType face_type)
Definition: faces_single_dimension.cpp:214
wmtk::simplex::faces
SimplexCollection faces(const Mesh &mesh, const Simplex &simplex, const bool sort_and_clean)
Returns all faces of a simplex.
Definition: faces.cpp:10
wmtk
Definition: Accessor.hpp:6
wmtk::PV
constexpr PrimitiveType PV
Definition: EnvelopeInvariant.cpp:24
wmtk::logger
spdlog::logger & logger()
Retrieves the current logger.
Definition: Logger.cpp:58
wmtk::Vector3d
Vector< double, 3 > Vector3d
Definition: Types.hpp:39
wmtk::PE
constexpr PrimitiveType PE
Definition: EnvelopeInvariant.cpp:25
tuple_vector_to_homogeneous_simplex_vector.hpp