Wildmeshing Toolkit
Loading...
Searching...
No Matches
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
11namespace wmtk::invariants {
12
13CollapseEnergyBeforeInvariant::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
24bool 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::create_const_accessor
const attribute::Accessor< T, Mesh, attribute::CachingAttribute< T >, 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::orient_vertices
virtual std::vector< Tuple > orient_vertices(const Tuple &t) const =0
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::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::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