Wildmeshing Toolkit
Loading...
Searching...
No Matches
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
11namespace wmtk::invariants {
12
13CollapseEnergyBeforeInvariantDouble::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
24bool 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::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::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::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_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::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