Wildmeshing Toolkit
Loading...
Searching...
No Matches
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
10namespace wmtk::invariants {
11
12CollapseSoftEnergyBeforeInvariant::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
29bool 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:953
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::Tuple
The Tuple is the basic navigation tool in our mesh data structure.
Definition Tuple.hpp:19
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::vertex
static Simplex vertex(const Mesh &m, const Tuple &t)
Definition Simplex.hpp:56
wmtk::simplex::Simplex::tuple
const Tuple & tuple() const
Definition Simplex.hpp:53
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:425
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::PE
constexpr PrimitiveType PE
Definition EnvelopeInvariant.cpp:25
wmtk::PV
constexpr PrimitiveType PV
Definition EnvelopeInvariant.cpp:24
top_dimension_cofaces.hpp