Wildmeshing Toolkit
CollapseSoftEnergyBeforeInvariant.cpp
Go to the documentation of this file.
1 #include "CollapseSoftEnergyBeforeInvariant.hpp"
2 #include <wmtk/simplex/utils/SimplexComparisons.hpp>
3 
4 #include <wmtk/function/utils/amips.hpp>
5 #include <wmtk/simplex/Simplex.hpp>
6 #include <wmtk/simplex/top_dimension_cofaces.hpp>
7 #include <wmtk/simplex/utils/SimplexComparisons.hpp>
8 #include <wmtk/utils/orient.hpp>
9 
10 namespace wmtk::invariants {
11 
12 CollapseSoftEnergyBeforeInvariant::CollapseSoftEnergyBeforeInvariant(
13  const Mesh& m,
14  const attribute::TypedAttributeHandle<Rational>& coordinate,
15  const attribute::TypedAttributeHandle<double>& energy,
16  const attribute::TypedAttributeHandle<double>& edge_length,
17  const attribute::TypedAttributeHandle<double>& target_edge_length,
18  int64_t collapse_type,
19  double eps)
20  : Invariant(m, true, false, false)
21  , m_coordinate_handle(coordinate)
22  , m_energy_handle(energy)
23  , m_edge_length_handle(edge_length)
24  , m_target_edge_length_handle(target_edge_length)
25  , m_collapse_type(collapse_type)
26  , m_eps(eps)
27 {}
28 
29 bool CollapseSoftEnergyBeforeInvariant::before(const simplex::Simplex& s) const
30 {
31  constexpr static PrimitiveType PV = PrimitiveType::Vertex;
32  constexpr static PrimitiveType PE = PrimitiveType::Edge;
33  constexpr static PrimitiveType PF = PrimitiveType::Triangle;
34  constexpr static PrimitiveType PT = PrimitiveType::Tetrahedron;
35 
36  if (mesh().top_simplex_type() == PT) {
37  auto position_accessor = mesh().create_const_accessor(m_coordinate_handle);
38 
39  auto amips_accessor = mesh().create_const_accessor(m_energy_handle);
40  auto edge_length_accessor = mesh().create_const_accessor(m_edge_length_handle);
41  auto target_edge_length_accessor =
42  mesh().create_const_accessor(m_target_edge_length_handle);
43 
44  const double current_length = edge_length_accessor.const_scalar_attribute(s.tuple());
45  const double target_length = target_edge_length_accessor.const_scalar_attribute(s.tuple());
46  const bool mini_edge = current_length <= target_length * m_eps;
47 
48  double old_energy_max = std::numeric_limits<double>::lowest();
49  double new_energy_max = std::numeric_limits<double>::lowest();
50 
51  switch (m_collapse_type) {
52  case 0: {
53  const Tuple v0 = s.tuple();
54  const Tuple v1 = mesh().switch_tuple(v0, PV);
55  const auto v1_pos = position_accessor.const_vector_attribute(v1);
56 
57  // TODO: check sort and clean
58  const auto& v0_incident_tets =
59  simplex::top_dimension_cofaces(mesh(), simplex::Simplex::vertex(mesh(), v0), false);
60  const auto& v1_incident_tets =
61  simplex::top_dimension_cofaces(mesh(), simplex::Simplex::vertex(mesh(), v1), false);
62 
63  // compute energy for v1 incident tets
64  for (const auto& t : v1_incident_tets) {
65  const std::array<Tuple, 4> lv = {
66  {t.tuple(),
67  mesh().switch_tuple(t.tuple(), PV),
68  mesh().switch_tuples(t.tuple(), {PE, PV}),
69  mesh().switch_tuples(t.tuple(), {PF, PE, PV})}};
70 
71  // compute old max
72  old_energy_max =
73  std::max(old_energy_max, amips_accessor.const_scalar_attribute(t.tuple()));
74 
75  // compute new energy
76  const simplex::Simplex v0_s(mesh(), PV, v0);
77 
78  if (simplex::utils::SimplexComparisons::equal(
79  mesh(),
80  simplex::Simplex::vertex(mesh(), lv[1]),
81  v0_s) ||
82  simplex::utils::SimplexComparisons::equal(
83  mesh(),
84  simplex::Simplex::vertex(mesh(), lv[2]),
85  v0_s) ||
86  simplex::utils::SimplexComparisons::equal(
87  mesh(),
88  simplex::Simplex::vertex(mesh(), lv[3]),
89  v0_s)) {
90  // skip if incident both vertices
91  continue;
92  }
93 
94  // position after the collapse
95  const std::array<Eigen::Vector<Rational, 3>, 4> lv_pos = {
96  {position_accessor.const_vector_attribute(v0), // v1 collapsed to v0
97  position_accessor.const_vector_attribute(lv[1]),
98  position_accessor.const_vector_attribute(lv[2]),
99  position_accessor.const_vector_attribute(lv[3])}};
100 
101  // check inversion, any new tet inverted return false
102  bool is_ccw_tuple = mesh().is_ccw(lv[0]);
103  if (is_ccw_tuple) {
104  if (wmtk::utils::wmtk_orient3d(lv_pos[3], lv_pos[0], lv_pos[1], lv_pos[2]) <=
105  0) {
106  return false;
107  }
108  } else {
109  if (wmtk::utils::wmtk_orient3d(lv_pos[3], lv_pos[0], lv_pos[2], lv_pos[1]) <=
110  0) {
111  return false;
112  }
113  }
114 
115  double new_energy;
116  if (is_ccw_tuple) {
117  new_energy = wmtk::function::utils::Tet_AMIPS_energy(
118  {{lv_pos[3][0].to_double(),
119  lv_pos[3][1].to_double(),
120  lv_pos[3][2].to_double(),
121  lv_pos[0][0].to_double(),
122  lv_pos[0][1].to_double(),
123  lv_pos[0][2].to_double(),
124  lv_pos[1][0].to_double(),
125  lv_pos[1][1].to_double(),
126  lv_pos[1][2].to_double(),
127  lv_pos[2][0].to_double(),
128  lv_pos[2][1].to_double(),
129  lv_pos[2][2].to_double()}});
130  } else {
131  new_energy = wmtk::function::utils::Tet_AMIPS_energy(
132  {{lv_pos[3][0].to_double(),
133  lv_pos[3][1].to_double(),
134  lv_pos[3][2].to_double(),
135  lv_pos[0][0].to_double(),
136  lv_pos[0][1].to_double(),
137  lv_pos[0][2].to_double(),
138  lv_pos[2][0].to_double(),
139  lv_pos[2][1].to_double(),
140  lv_pos[2][2].to_double(),
141  lv_pos[1][0].to_double(),
142  lv_pos[1][1].to_double(),
143  lv_pos[1][2].to_double()}});
144  }
145  new_energy_max = std::max(new_energy_max, new_energy);
146  }
147 
148  // passed all inversion check, return true if v1 is not rounded
149  for (int64_t i = 0; i < 3; ++i) {
150  // if the collapse from vertex is not rounded, return true anyway
151  if (!v1_pos[i].is_rounded()) return true;
152  }
153 
154  // collapse a tiny edge anyway
155  if (mini_edge) {
156  return true;
157  }
158 
159  // compute energy for v0 incident, old and new are the same because v0 is not moved
160  for (const auto& t : v0_incident_tets) {
161  const std::array<Tuple, 4> lv = {
162  {t.tuple(),
163  mesh().switch_tuple(t.tuple(), PV),
164  mesh().switch_tuples(t.tuple(), {PE, PV}),
165  mesh().switch_tuples(t.tuple(), {PF, PE, PV})}};
166 
167  // skip tets incident the edge, old is already computed above
168  const simplex::Simplex v1_s(mesh(), PV, v1);
169  if (simplex::utils::SimplexComparisons::equal(
170  mesh(),
171  simplex::Simplex::vertex(mesh(), lv[1]),
172  v1_s) ||
173  simplex::utils::SimplexComparisons::equal(
174  mesh(),
175  simplex::Simplex::vertex(mesh(), lv[2]),
176  v1_s) ||
177  simplex::utils::SimplexComparisons::equal(
178  mesh(),
179  simplex::Simplex::vertex(mesh(), lv[3]),
180  v1_s)) {
181  // skip if incident both vertices
182  continue;
183  }
184 
185  // compute old max
186  const double energy = amips_accessor.const_scalar_attribute(t.tuple());
187  old_energy_max = std::max(old_energy_max, energy);
188 
189  // compute new energy
190  new_energy_max = std::max(new_energy_max, energy);
191  }
192 
193  break;
194  }
195  case 1: {
196  const Tuple v0 = s.tuple();
197  const Tuple v1 = mesh().switch_tuple(v0, PV);
198  const auto v0_pos = position_accessor.const_vector_attribute(v0);
199 
200  // TODO: check sort and clean
201  const auto& v0_incident_tets =
202  simplex::top_dimension_cofaces(mesh(), simplex::Simplex::vertex(mesh(), v0), false);
203  const auto& v1_incident_tets =
204  simplex::top_dimension_cofaces(mesh(), simplex::Simplex::vertex(mesh(), v1), false);
205 
206  // compute energy for v0 incident tets
207  for (const auto& t : v0_incident_tets) {
208  const std::array<Tuple, 4> lv = {
209  {t.tuple(),
210  mesh().switch_tuple(t.tuple(), PV),
211  mesh().switch_tuples(t.tuple(), {PE, PV}),
212  mesh().switch_tuples(t.tuple(), {PF, PE, PV})}};
213 
214  // compute old max
215  old_energy_max =
216  std::max(old_energy_max, amips_accessor.const_scalar_attribute(t.tuple()));
217 
218  // compute new energy
219  const simplex::Simplex v1_s(mesh(), PV, v1);
220  if (simplex::utils::SimplexComparisons::equal(
221  mesh(),
222  simplex::Simplex::vertex(mesh(), lv[1]),
223  v1_s) ||
224  simplex::utils::SimplexComparisons::equal(
225  mesh(),
226  simplex::Simplex::vertex(mesh(), lv[2]),
227  v1_s) ||
228  simplex::utils::SimplexComparisons::equal(
229  mesh(),
230  simplex::Simplex::vertex(mesh(), lv[3]),
231  v1_s)) {
232  // skip if incident both vertices
233  continue;
234  }
235 
236  // position after the collapse
237  const std::array<Eigen::Vector<Rational, 3>, 4> lv_pos = {
238  {position_accessor.const_vector_attribute(v1), // v0 collapsed to v1
239  position_accessor.const_vector_attribute(lv[1]),
240  position_accessor.const_vector_attribute(lv[2]),
241  position_accessor.const_vector_attribute(lv[3])}};
242 
243  // check inversion, any new tet inverted return false
244  bool is_ccw_tuple = mesh().is_ccw(lv[0]);
245  if (is_ccw_tuple) {
246  if (wmtk::utils::wmtk_orient3d(lv_pos[3], lv_pos[0], lv_pos[1], lv_pos[2]) <=
247  0) {
248  return false;
249  }
250  } else {
251  if (wmtk::utils::wmtk_orient3d(lv_pos[3], lv_pos[0], lv_pos[2], lv_pos[1]) <=
252  0) {
253  return false;
254  }
255  }
256 
257  double new_energy;
258  if (is_ccw_tuple) {
259  new_energy = wmtk::function::utils::Tet_AMIPS_energy(
260  {{lv_pos[3][0].to_double(),
261  lv_pos[3][1].to_double(),
262  lv_pos[3][2].to_double(),
263  lv_pos[0][0].to_double(),
264  lv_pos[0][1].to_double(),
265  lv_pos[0][2].to_double(),
266  lv_pos[1][0].to_double(),
267  lv_pos[1][1].to_double(),
268  lv_pos[1][2].to_double(),
269  lv_pos[2][0].to_double(),
270  lv_pos[2][1].to_double(),
271  lv_pos[2][2].to_double()}});
272  } else {
273  new_energy = wmtk::function::utils::Tet_AMIPS_energy(
274  {{lv_pos[3][0].to_double(),
275  lv_pos[3][1].to_double(),
276  lv_pos[3][2].to_double(),
277  lv_pos[0][0].to_double(),
278  lv_pos[0][1].to_double(),
279  lv_pos[0][2].to_double(),
280  lv_pos[2][0].to_double(),
281  lv_pos[2][1].to_double(),
282  lv_pos[2][2].to_double(),
283  lv_pos[1][0].to_double(),
284  lv_pos[1][1].to_double(),
285  lv_pos[1][2].to_double()}});
286  }
287  new_energy_max = std::max(new_energy_max, new_energy);
288  }
289 
290  // passed all inversion check, return true if v1 is not rounded
291  for (int64_t i = 0; i < 3; ++i) {
292  // if the collapse from vertex is not rounded, return true anyway
293  if (!v0_pos[i].is_rounded()) return true;
294  }
295 
296  // collapse a tiny edge anyway
297  if (mini_edge) {
298  return true;
299  }
300 
301  // compute energy for v1 incident, old and new are the same because v0 is not moved
302  for (const auto& t : v1_incident_tets) {
303  const std::array<Tuple, 4> lv = {
304  {t.tuple(),
305  mesh().switch_tuple(t.tuple(), PV),
306  mesh().switch_tuples(t.tuple(), {PE, PV}),
307  mesh().switch_tuples(t.tuple(), {PF, PE, PV})}};
308 
309  // skip tets incident the edge, old is already computed above
310  const simplex::Simplex v0_s(mesh(), PV, v0);
311  if (simplex::utils::SimplexComparisons::equal(
312  mesh(),
313  simplex::Simplex::vertex(mesh(), lv[1]),
314  v0_s) ||
315  simplex::utils::SimplexComparisons::equal(
316  mesh(),
317  simplex::Simplex::vertex(mesh(), lv[2]),
318  v0_s) ||
319  simplex::utils::SimplexComparisons::equal(
320  mesh(),
321  simplex::Simplex::vertex(mesh(), lv[3]),
322  v0_s)) {
323  // skip if incident both vertices
324  continue;
325  }
326 
327  // compute old max
328  const double energy = amips_accessor.const_scalar_attribute(t.tuple());
329  old_energy_max = std::max(old_energy_max, energy);
330 
331  // compute new energy
332  new_energy_max = std::max(new_energy_max, energy);
333  }
334 
335  break;
336  }
337 
338  case 2: {
339  throw std::runtime_error("not implemented");
340  }
341 
342  default: throw std::runtime_error("Invalid collapse type");
343  }
344 
345  return new_energy_max <= old_energy_max;
346 
347  } else if (mesh().top_simplex_type() == PF) {
348  // TODO: not implement yet
349  return true;
350  }
351 
352  return true;
353 }
354 
355 } // namespace wmtk::invariants
wmtk::Mesh::create_const_accessor
const attribute::Accessor< T, Mesh, D > create_const_accessor(const attribute::MeshAttributeHandle &handle) const
wmtk::Mesh::is_ccw
virtual bool is_ccw(const Tuple &tuple) const =0
TODO this needs dimension?
wmtk::Mesh::switch_tuples
Tuple switch_tuples(const Tuple &tuple, const ContainerType &op_sequence) const
Definition: Mesh.hpp:967
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::attribute::TypedAttributeHandle
Handle that represents attributes for some mesh.
Definition: TypedAttributeHandle.hpp:28
wmtk::invariants::CollapseSoftEnergyBeforeInvariant::m_coordinate_handle
const attribute::TypedAttributeHandle< Rational > m_coordinate_handle
Definition: CollapseSoftEnergyBeforeInvariant.hpp:23
wmtk::invariants::CollapseSoftEnergyBeforeInvariant::m_target_edge_length_handle
const attribute::TypedAttributeHandle< double > m_target_edge_length_handle
Definition: CollapseSoftEnergyBeforeInvariant.hpp:26
wmtk::invariants::CollapseSoftEnergyBeforeInvariant::m_energy_handle
const attribute::TypedAttributeHandle< double > m_energy_handle
Definition: CollapseSoftEnergyBeforeInvariant.hpp:24
wmtk::invariants::CollapseSoftEnergyBeforeInvariant::CollapseSoftEnergyBeforeInvariant
CollapseSoftEnergyBeforeInvariant(const Mesh &m, const attribute::TypedAttributeHandle< Rational > &coordinate, const attribute::TypedAttributeHandle< double > &energy, const attribute::TypedAttributeHandle< double > &edge_length, const attribute::TypedAttributeHandle< double > &target_edge_length, int64_t collapse_type=0, double eps=0.1)
Definition: CollapseSoftEnergyBeforeInvariant.cpp:12
wmtk::invariants::CollapseSoftEnergyBeforeInvariant::m_edge_length_handle
const attribute::TypedAttributeHandle< double > m_edge_length_handle
Definition: CollapseSoftEnergyBeforeInvariant.hpp:25
wmtk::invariants::Invariant::mesh
const Mesh & mesh() const
Definition: Invariant.cpp:35
wmtk::simplex::Simplex::tuple
const Tuple & tuple() const
Definition: Simplex.hpp:53
wmtk::simplex::Simplex::vertex
static Simplex vertex(const Mesh &m, const Tuple &t)
Definition: Simplex.hpp:56
wmtk::simplex::utils::SimplexComparisons::equal
static bool equal(const Mesh &m, const Simplex &s0, const Simplex &s1)
Definition: SimplexComparisons.hxx:8
PT
constexpr wmtk::PrimitiveType PT
Definition: faces_single_dimension.cpp:6
PF
constexpr wmtk::PrimitiveType PF
Definition: faces_single_dimension.cpp:5
wmtk::function::utils::Tet_AMIPS_energy
double Tet_AMIPS_energy(const std::array< double, 12 > &T)
Definition: amips.cpp:380
wmtk::simplex::top_dimension_cofaces
void top_dimension_cofaces(const Simplex &simplex, SimplexCollection &simplex_collection, const bool sort_and_clean)
Get all top dimension cofaces of the given simplex.
Definition: top_dimension_cofaces.cpp:428
wmtk::utils::wmtk_orient3d
int wmtk_orient3d(const Eigen::Ref< const Eigen::Vector3< Rational >> &p0, const Eigen::Ref< const Eigen::Vector3< Rational >> &p1, const Eigen::Ref< const Eigen::Vector3< Rational >> &p2, const Eigen::Ref< const Eigen::Vector3< Rational >> &p3)
Definition: orient.cpp:75
wmtk::PV
constexpr PrimitiveType PV
Definition: EnvelopeInvariant.cpp:24
wmtk::PE
constexpr PrimitiveType PE
Definition: EnvelopeInvariant.cpp:25
top_dimension_cofaces.hpp