Wildmeshing Toolkit
CollapseEnergyBeforeInvariant.cpp
Go to the documentation of this file.
1 #include "CollapseEnergyBeforeInvariant.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/utils/orient.hpp>
8 
9 #include <wmtk/utils/Logger.hpp>
10 
11 namespace wmtk::invariants {
12 
13 CollapseEnergyBeforeInvariant::CollapseEnergyBeforeInvariant(
14  const Mesh& m,
15  const attribute::TypedAttributeHandle<Rational>& coordinate,
16  const attribute::TypedAttributeHandle<double>& energy,
17  int64_t collapse_type)
18  : Invariant(m, true, false, false)
19  , m_coordinate_handle(coordinate)
20  , m_energy_handle(energy)
21  , m_collapse_type(collapse_type)
22 {}
23 
24 bool CollapseEnergyBeforeInvariant::before(const simplex::Simplex& s) const
25 {
26  constexpr static PrimitiveType PV = PrimitiveType::Vertex;
27  constexpr static PrimitiveType PE = PrimitiveType::Edge;
28  constexpr static PrimitiveType PF = PrimitiveType::Triangle;
29  constexpr static PrimitiveType PT = PrimitiveType::Tetrahedron;
30 
31  if (mesh().top_simplex_type() == PT) {
32  auto position_accessor = mesh().create_const_accessor(m_coordinate_handle);
33 
34  auto amips_accessor = mesh().create_const_accessor(m_energy_handle);
35 
36  double old_energy_max = std::numeric_limits<double>::lowest();
37  double new_energy_max = std::numeric_limits<double>::lowest();
38 
39  switch (m_collapse_type) {
40  case 0: {
41  const Tuple v0 = s.tuple();
42  const Tuple v1 = mesh().switch_tuple(v0, PV);
43  const simplex::Simplex v0_s(mesh(), PV, v0);
44  const simplex::Simplex v1_s(mesh(), PV, v1);
45  const auto v1_pos = position_accessor.const_vector_attribute(v1);
46 
47  // TODO: check sort and clean
48  const auto& v0_incident_tets =
49  simplex::top_dimension_cofaces(mesh(), simplex::Simplex::vertex(mesh(), v0), false);
50  const auto& v1_incident_tets =
51  simplex::top_dimension_cofaces(mesh(), simplex::Simplex::vertex(mesh(), v1), false);
52 
53  // compute energy for v1 incident tets
54  for (const auto& t : v1_incident_tets) {
55  auto lv = mesh().orient_vertices(t.tuple()); // get orient tet vertices
56  assert(lv.size() == 4);
57 
58  // compute old max
59  old_energy_max =
60  std::max(old_energy_max, amips_accessor.const_scalar_attribute(t.tuple()));
61 
62  // skip if incident both vertices
63  bool incident_both = false;
64  for (int i = 0; i < 4; ++i) {
65  if (simplex::utils::SimplexComparisons::equal(
66  mesh(),
67  simplex::Simplex::vertex(mesh(), lv[i]),
68  v0_s)) {
69  incident_both = true;
70  break;
71  }
72  }
73  if (incident_both) continue;
74 
75  // position after the collapse
76  std::array<Eigen::Vector<Rational, 3>, 4> lv_pos = {
77  {position_accessor.const_vector_attribute(lv[0]),
78  position_accessor.const_vector_attribute(lv[1]),
79  position_accessor.const_vector_attribute(lv[2]),
80  position_accessor.const_vector_attribute(lv[3])}};
81 
82  for (int i = 0; i < 4; ++i) {
83  if (simplex::utils::SimplexComparisons::equal(
84  mesh(),
85  simplex::Simplex::vertex(mesh(), lv[i]),
86  v1_s)) {
87  // change the v1 entry to v0
88  lv_pos[i] = position_accessor.const_vector_attribute(v0);
89  break;
90  }
91  assert(i != 3); // should find one and break;
92  }
93 
94  // check inversion, any new tet inverted return false
95  if (wmtk::utils::wmtk_orient3d(lv_pos[0], lv_pos[1], lv_pos[2], lv_pos[3]) <= 0) {
96  return false;
97  }
98 
99  double new_energy = wmtk::function::utils::Tet_AMIPS_energy(
100  {{lv_pos[0][0].to_double(),
101  lv_pos[0][1].to_double(),
102  lv_pos[0][2].to_double(),
103  lv_pos[1][0].to_double(),
104  lv_pos[1][1].to_double(),
105  lv_pos[1][2].to_double(),
106  lv_pos[2][0].to_double(),
107  lv_pos[2][1].to_double(),
108  lv_pos[2][2].to_double(),
109  lv_pos[3][0].to_double(),
110  lv_pos[3][1].to_double(),
111  lv_pos[3][2].to_double()}});
112 
113  new_energy_max = std::max(new_energy_max, new_energy);
114  }
115 
116  // passed all inversion check, return true if v1 is not rounded
117  for (int64_t i = 0; i < 3; ++i) {
118  // if the collapse from vertex is not rounded, return true anyway
119  if (!v1_pos[i].is_rounded()) return true;
120  }
121 
122  // compute energy for v0 incident, old and new are the same because v0 is not moved
123  for (const auto& t : v0_incident_tets) {
124  auto lv = mesh().orient_vertices(t.tuple()); // get orient tet vertices
125 
126  // skip if incident both vertices
127  bool incident_both = false;
128  for (int i = 0; i < 4; ++i) {
129  if (simplex::utils::SimplexComparisons::equal(
130  mesh(),
131  simplex::Simplex::vertex(mesh(), lv[i]),
132  v1_s)) {
133  incident_both = true;
134  break;
135  }
136  }
137  if (incident_both) continue;
138 
139  // compute old max
140  const double energy = amips_accessor.const_scalar_attribute(t.tuple());
141  old_energy_max = std::max(old_energy_max, energy);
142 
143  // compute new energy
144  new_energy_max = std::max(new_energy_max, energy);
145  }
146 
147  break;
148  }
149  case 1: {
150  const Tuple v0 = s.tuple();
151  const Tuple v1 = mesh().switch_tuple(v0, PV);
152  const simplex::Simplex v0_s(mesh(), PV, v0);
153  const simplex::Simplex v1_s(mesh(), PV, v1);
154  const auto v0_pos = position_accessor.const_vector_attribute(v0);
155 
156  // TODO: check sort and clean
157  const auto& v0_incident_tets =
158  simplex::top_dimension_cofaces(mesh(), simplex::Simplex::vertex(mesh(), v0), false);
159  const auto& v1_incident_tets =
160  simplex::top_dimension_cofaces(mesh(), simplex::Simplex::vertex(mesh(), v1), false);
161 
162  // compute energy for v0 incident tets
163  for (const auto& t : v0_incident_tets) {
164  auto lv = mesh().orient_vertices(t.tuple()); // get orient tet vertices
165  assert(lv.size() == 4);
166 
167  // compute old max
168  old_energy_max =
169  std::max(old_energy_max, amips_accessor.const_scalar_attribute(t.tuple()));
170 
171  // compute new energy
172  // skip if incident both vertices
173  bool incident_both = false;
174  for (int i = 0; i < 4; ++i) {
175  if (simplex::utils::SimplexComparisons::equal(
176  mesh(),
177  simplex::Simplex::vertex(mesh(), lv[i]),
178  v1_s)) {
179  incident_both = true;
180  break;
181  }
182  }
183  if (incident_both) continue;
184 
185  // position after the collapse
186  std::array<Eigen::Vector<Rational, 3>, 4> lv_pos = {
187  {position_accessor.const_vector_attribute(lv[0]),
188  position_accessor.const_vector_attribute(lv[1]),
189  position_accessor.const_vector_attribute(lv[2]),
190  position_accessor.const_vector_attribute(lv[3])}};
191 
192  for (int i = 0; i < 4; ++i) {
193  if (simplex::utils::SimplexComparisons::equal(
194  mesh(),
195  simplex::Simplex::vertex(mesh(), lv[i]),
196  v0_s)) {
197  // change the v0 entry to v1
198  lv_pos[i] = position_accessor.const_vector_attribute(v1);
199  break;
200  }
201  assert(i != 3); // should find one and break;
202  }
203 
204  // check inversion, any new tet inverted return false
205  if (wmtk::utils::wmtk_orient3d(lv_pos[0], lv_pos[1], lv_pos[2], lv_pos[3]) <= 0) {
206  return false;
207  }
208 
209  double new_energy = wmtk::function::utils::Tet_AMIPS_energy(
210  {{lv_pos[0][0].to_double(),
211  lv_pos[0][1].to_double(),
212  lv_pos[0][2].to_double(),
213  lv_pos[1][0].to_double(),
214  lv_pos[1][1].to_double(),
215  lv_pos[1][2].to_double(),
216  lv_pos[2][0].to_double(),
217  lv_pos[2][1].to_double(),
218  lv_pos[2][2].to_double(),
219  lv_pos[3][0].to_double(),
220  lv_pos[3][1].to_double(),
221  lv_pos[3][2].to_double()}});
222  new_energy_max = std::max(new_energy_max, new_energy);
223  }
224 
225  // passed all inversion check, return true if v0 is not rounded
226  for (int64_t i = 0; i < 3; ++i) {
227  // if the collapse from vertex is not rounded, return true anyway
228  if (!v0_pos[i].is_rounded()) return true;
229  }
230 
231  // compute energy for v1 incident, old and new are the same because v0 is not moved
232  for (const auto& t : v1_incident_tets) {
233  auto lv = mesh().orient_vertices(t.tuple()); // get orient tet vertices
234 
235  // skip if incident both vertices
236  bool incident_both = false;
237  for (int i = 0; i < 4; ++i) {
238  if (simplex::utils::SimplexComparisons::equal(
239  mesh(),
240  simplex::Simplex::vertex(mesh(), lv[i]),
241  v0_s)) {
242  incident_both = true;
243  break;
244  }
245  }
246  if (incident_both) continue;
247 
248  // compute old max
249  const double energy = amips_accessor.const_scalar_attribute(t.tuple());
250  old_energy_max = std::max(old_energy_max, energy);
251 
252  // compute new energy
253  new_energy_max = std::max(new_energy_max, energy);
254  }
255 
256  break;
257  }
258 
259  case 2: {
260  throw std::runtime_error("not implemented");
261  }
262 
263  default: throw std::runtime_error("Invalid collapse type");
264  }
265 
266 
267  return new_energy_max <= old_energy_max;
268 
269  } else if (mesh().top_simplex_type() == PF) {
270  // TODO: not implement yet
271  return true;
272  }
273 
274  return true;
275 }
276 
277 } // namespace wmtk::invariants
wmtk::Mesh::orient_vertices
virtual std::vector< Tuple > orient_vertices(const Tuple &t) const =0
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::attribute::TypedAttributeHandle
Handle that represents attributes for some mesh.
Definition: TypedAttributeHandle.hpp:28
wmtk::invariants::CollapseEnergyBeforeInvariant::CollapseEnergyBeforeInvariant
CollapseEnergyBeforeInvariant(const Mesh &m, const attribute::TypedAttributeHandle< Rational > &coordinate, const attribute::TypedAttributeHandle< double > &energy, int64_t collapse_type=0)
Definition: CollapseEnergyBeforeInvariant.cpp:13
wmtk::invariants::CollapseEnergyBeforeInvariant::m_coordinate_handle
const attribute::TypedAttributeHandle< Rational > m_coordinate_handle
Definition: CollapseEnergyBeforeInvariant.hpp:20
wmtk::invariants::CollapseEnergyBeforeInvariant::m_energy_handle
const attribute::TypedAttributeHandle< double > m_energy_handle
Definition: CollapseEnergyBeforeInvariant.hpp:21
wmtk::invariants::CollapseEnergyBeforeInvariant::before
bool before(const simplex::Simplex &s) const override
Definition: CollapseEnergyBeforeInvariant.cpp:24
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