Wildmeshing Toolkit
CollapseEnergyBeforeInvariantDouble.cpp
Go to the documentation of this file.
1 #include "CollapseEnergyBeforeInvariantDouble.hpp"
2 
3 #include <wmtk/function/utils/amips.hpp>
4 #include <wmtk/simplex/Simplex.hpp>
5 #include <wmtk/simplex/top_dimension_cofaces.hpp>
6 #include <wmtk/simplex/utils/SimplexComparisons.hpp>
7 #include <wmtk/utils/orient.hpp>
8 
9 #include <wmtk/utils/Logger.hpp>
10 
11 namespace wmtk::invariants {
12 
13 CollapseEnergyBeforeInvariantDouble::CollapseEnergyBeforeInvariantDouble(
14  const Mesh& m,
15  const attribute::TypedAttributeHandle<double>& 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 CollapseEnergyBeforeInvariantDouble::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<double, 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],
101  lv_pos[0][1],
102  lv_pos[0][2],
103  lv_pos[1][0],
104  lv_pos[1][1],
105  lv_pos[1][2],
106  lv_pos[2][0],
107  lv_pos[2][1],
108  lv_pos[2][2],
109  lv_pos[3][0],
110  lv_pos[3][1],
111  lv_pos[3][2]}});
112 
113  new_energy_max = std::max(new_energy_max, new_energy);
114  }
115 
116  // compute energy for v0 incident, old and new are the same because v0 is not moved
117  for (const auto& t : v0_incident_tets) {
118  auto lv = mesh().orient_vertices(t.tuple()); // get orient tet vertices
119 
120  // skip if incident both vertices
121  bool incident_both = false;
122  for (int i = 0; i < 4; ++i) {
123  if (simplex::utils::SimplexComparisons::equal(
124  mesh(),
125  simplex::Simplex::vertex(mesh(), lv[i]),
126  v1_s)) {
127  incident_both = true;
128  break;
129  }
130  }
131  if (incident_both) continue;
132 
133  // compute old max
134  const double energy = amips_accessor.const_scalar_attribute(t.tuple());
135  old_energy_max = std::max(old_energy_max, energy);
136 
137  // compute new energy
138  new_energy_max = std::max(new_energy_max, energy);
139  }
140 
141  break;
142  }
143  case 1: {
144  const Tuple v0 = s.tuple();
145  const Tuple v1 = mesh().switch_tuple(v0, PV);
146  const simplex::Simplex v0_s(mesh(), PV, v0);
147  const simplex::Simplex v1_s(mesh(), PV, v1);
148  // const auto v0_pos = position_accessor.const_vector_attribute(v0);
149 
150  // TODO: check sort and clean
151  const auto& v0_incident_tets =
152  simplex::top_dimension_cofaces(mesh(), simplex::Simplex::vertex(mesh(), v0), false);
153  const auto& v1_incident_tets =
154  simplex::top_dimension_cofaces(mesh(), simplex::Simplex::vertex(mesh(), v1), false);
155 
156  // compute energy for v0 incident tets
157  for (const auto& t : v0_incident_tets) {
158  auto lv = mesh().orient_vertices(t.tuple()); // get orient tet vertices
159  assert(lv.size() == 4);
160 
161  // compute old max
162  old_energy_max =
163  std::max(old_energy_max, amips_accessor.const_scalar_attribute(t.tuple()));
164 
165  // compute new energy
166  // skip if incident both vertices
167  bool incident_both = false;
168  for (int i = 0; i < 4; ++i) {
169  if (simplex::utils::SimplexComparisons::equal(
170  mesh(),
171  simplex::Simplex::vertex(mesh(), lv[i]),
172  v1_s)) {
173  incident_both = true;
174  break;
175  }
176  }
177  if (incident_both) continue;
178 
179  // position after the collapse
180  std::array<Eigen::Vector<double, 3>, 4> lv_pos = {
181  {position_accessor.const_vector_attribute(lv[0]),
182  position_accessor.const_vector_attribute(lv[1]),
183  position_accessor.const_vector_attribute(lv[2]),
184  position_accessor.const_vector_attribute(lv[3])}};
185 
186  for (int i = 0; i < 4; ++i) {
187  if (simplex::utils::SimplexComparisons::equal(
188  mesh(),
189  simplex::Simplex::vertex(mesh(), lv[i]),
190  v0_s)) {
191  // change the v0 entry to v1
192  lv_pos[i] = position_accessor.const_vector_attribute(v1);
193  break;
194  }
195  assert(i != 3); // should find one and break;
196  }
197 
198  // check inversion, any new tet inverted return false
199  if (wmtk::utils::wmtk_orient3d(lv_pos[0], lv_pos[1], lv_pos[2], lv_pos[3]) <= 0) {
200  return false;
201  }
202 
203  double new_energy = wmtk::function::utils::Tet_AMIPS_energy(
204  {{lv_pos[0][0],
205  lv_pos[0][1],
206  lv_pos[0][2],
207  lv_pos[1][0],
208  lv_pos[1][1],
209  lv_pos[1][2],
210  lv_pos[2][0],
211  lv_pos[2][1],
212  lv_pos[2][2],
213  lv_pos[3][0],
214  lv_pos[3][1],
215  lv_pos[3][2]}});
216  new_energy_max = std::max(new_energy_max, new_energy);
217  }
218 
219  // compute energy for v1 incident, old and new are the same because v0 is not moved
220  for (const auto& t : v1_incident_tets) {
221  auto lv = mesh().orient_vertices(t.tuple()); // get orient tet vertices
222 
223  // skip if incident both vertices
224  bool incident_both = false;
225  for (int i = 0; i < 4; ++i) {
226  if (simplex::utils::SimplexComparisons::equal(
227  mesh(),
228  simplex::Simplex::vertex(mesh(), lv[i]),
229  v0_s)) {
230  incident_both = true;
231  break;
232  }
233  }
234  if (incident_both) continue;
235 
236  // compute old max
237  const double energy = amips_accessor.const_scalar_attribute(t.tuple());
238  old_energy_max = std::max(old_energy_max, energy);
239 
240  // compute new energy
241  new_energy_max = std::max(new_energy_max, energy);
242  }
243 
244  break;
245  }
246 
247  case 2: {
248  throw std::runtime_error("not implemented");
249  }
250 
251  default: throw std::runtime_error("Invalid collapse type");
252  }
253 
254 
255  return new_energy_max <= old_energy_max;
256 
257  } else if (mesh().top_simplex_type() == PF) {
258  // TODO: not implement yet
259  return true;
260  }
261 
262  return true;
263 }
264 
265 } // 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::invariants::CollapseEnergyBeforeInvariantDouble::CollapseEnergyBeforeInvariantDouble
CollapseEnergyBeforeInvariantDouble(const Mesh &m, const attribute::TypedAttributeHandle< double > &coordinate, const attribute::TypedAttributeHandle< double > &energy, int64_t collapse_type=0)
Definition: CollapseEnergyBeforeInvariantDouble.cpp:13
wmtk::invariants::CollapseEnergyBeforeInvariantDouble::m_energy_handle
const attribute::TypedAttributeHandle< double > m_energy_handle
Definition: CollapseEnergyBeforeInvariantDouble.hpp:21
wmtk::invariants::CollapseEnergyBeforeInvariantDouble::m_coordinate_handle
const attribute::TypedAttributeHandle< double > m_coordinate_handle
Definition: CollapseEnergyBeforeInvariantDouble.hpp:20
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