Wildmeshing Toolkit
Loading...
Searching...
No Matches
Swap56EnergyBeforeInvariant.cpp
Go to the documentation of this file.
1#include "Swap56EnergyBeforeInvariant.hpp"
2#include <wmtk/Mesh.hpp>
3#include <wmtk/function/utils/amips.hpp>
4#include <wmtk/simplex/top_dimension_cofaces.hpp>
5#include <wmtk/utils/orient.hpp>
6
7namespace wmtk {
8Swap56EnergyBeforeInvariant::Swap56EnergyBeforeInvariant(
9 const Mesh& m,
10 const attribute::TypedAttributeHandle<Rational>& coordinate,
11 int64_t collapse_index,
12 double eps)
13 : Invariant(m, true, false, false)
14 , m_coordinate_handle(coordinate)
15 , m_collapse_index(collapse_index)
16 , m_eps(eps)
17{}
18
19bool Swap56EnergyBeforeInvariant::before(const simplex::Simplex& t) const
20{
21 assert(simplex::top_dimension_cofaces(mesh(), t).size() == 5);
22 constexpr static PrimitiveType PV = PrimitiveType::Vertex;
23 constexpr static PrimitiveType PE = PrimitiveType::Edge;
24 constexpr static PrimitiveType PF = PrimitiveType::Triangle;
25 constexpr static PrimitiveType PT = PrimitiveType::Tetrahedron;
26
27 auto accessor = mesh().create_const_accessor(m_coordinate_handle);
28
29 const Tuple e0 = t.tuple();
30 const Tuple e1 = mesh().switch_tuple(e0, PV);
31
32 std::array<Tuple, 5> v;
33 auto iter_tuple = e0;
34 for (int64_t i = 0; i < 5; ++i) {
35 v[i] = mesh().switch_tuples(iter_tuple, {PE, PV});
36 iter_tuple = mesh().switch_tuples(iter_tuple, {PF, PT});
37 }
38 assert(iter_tuple == e0);
39
40 // five iterable vertices remap to 0-4 by m_collapse_index, 0: m_collapse_index, 5: e0, 6: e1
41 std::array<Eigen::Vector3<Rational>, 7> positions = {
42 {accessor.const_vector_attribute(v[(m_collapse_index + 0) % 5]),
43 accessor.const_vector_attribute(v[(m_collapse_index + 1) % 5]),
44 accessor.const_vector_attribute(v[(m_collapse_index + 2) % 5]),
45 accessor.const_vector_attribute(v[(m_collapse_index + 3) % 5]),
46 accessor.const_vector_attribute(v[(m_collapse_index + 4) % 5]),
47 accessor.const_vector_attribute(e0),
48 accessor.const_vector_attribute(e1)}};
49
50 std::array<Eigen::Vector3d, 7> positions_double = {
51 {positions[0].cast<double>(),
52 positions[1].cast<double>(),
53 positions[2].cast<double>(),
54 positions[3].cast<double>(),
55 positions[4].cast<double>(),
56 positions[5].cast<double>(),
57 positions[6].cast<double>()}};
58
59 std::array<std::array<int, 4>, 5> old_tets = {
60 {{{0, 1, 5, 6}}, {{1, 2, 5, 6}}, {{2, 3, 5, 6}}, {{3, 4, 5, 6}}, {{4, 0, 5, 6}}}};
61
62 std::array<std::array<int, 4>, 6> new_tets = {
63 {{{0, 1, 2, 5}},
64 {{0, 2, 3, 5}},
65 {{0, 3, 4, 5}},
66 {{0, 1, 2, 6}},
67 {{0, 2, 3, 6}},
68 {{0, 3, 4, 6}}}};
69
70 double old_energy_max = std::numeric_limits<double>::lowest();
71 double new_energy_max = std::numeric_limits<double>::lowest();
72
73 for (int i = 0; i < 5; ++i) {
74 if (utils::wmtk_orient3d(
75 positions[old_tets[i][0]],
76 positions[old_tets[i][1]],
77 positions[old_tets[i][2]],
78 positions[old_tets[i][3]]) > 0) {
79 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
80 positions_double[old_tets[i][0]][0],
81 positions_double[old_tets[i][0]][1],
82 positions_double[old_tets[i][0]][2],
83 positions_double[old_tets[i][1]][0],
84 positions_double[old_tets[i][1]][1],
85 positions_double[old_tets[i][1]][2],
86 positions_double[old_tets[i][2]][0],
87 positions_double[old_tets[i][2]][1],
88 positions_double[old_tets[i][2]][2],
89 positions_double[old_tets[i][3]][0],
90 positions_double[old_tets[i][3]][1],
91 positions_double[old_tets[i][3]][2],
92 }});
93
94
95 if (energy > old_energy_max) old_energy_max = energy;
96 } else {
97 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
98 positions_double[old_tets[i][1]][0],
99 positions_double[old_tets[i][1]][1],
100 positions_double[old_tets[i][1]][2],
101 positions_double[old_tets[i][0]][0],
102 positions_double[old_tets[i][0]][1],
103 positions_double[old_tets[i][0]][2],
104 positions_double[old_tets[i][2]][0],
105 positions_double[old_tets[i][2]][1],
106 positions_double[old_tets[i][2]][2],
107 positions_double[old_tets[i][3]][0],
108 positions_double[old_tets[i][3]][1],
109 positions_double[old_tets[i][3]][2],
110 }});
111
112
113 if (energy > old_energy_max) old_energy_max = energy;
114 }
115 }
116
117 for (int i = 0; i < 6; ++i) {
118 if (utils::wmtk_orient3d(
119 positions[new_tets[i][0]],
120 positions[new_tets[i][1]],
121 positions[new_tets[i][2]],
122 positions[new_tets[i][3]]) > 0) {
123 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
124 positions_double[new_tets[i][0]][0],
125 positions_double[new_tets[i][0]][1],
126 positions_double[new_tets[i][0]][2],
127 positions_double[new_tets[i][1]][0],
128 positions_double[new_tets[i][1]][1],
129 positions_double[new_tets[i][1]][2],
130 positions_double[new_tets[i][2]][0],
131 positions_double[new_tets[i][2]][1],
132 positions_double[new_tets[i][2]][2],
133 positions_double[new_tets[i][3]][0],
134 positions_double[new_tets[i][3]][1],
135 positions_double[new_tets[i][3]][2],
136 }});
137
138
139 if (energy > new_energy_max) new_energy_max = energy;
140 } else {
141 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
142 positions_double[new_tets[i][1]][0],
143 positions_double[new_tets[i][1]][1],
144 positions_double[new_tets[i][1]][2],
145 positions_double[new_tets[i][0]][0],
146 positions_double[new_tets[i][0]][1],
147 positions_double[new_tets[i][0]][2],
148 positions_double[new_tets[i][2]][0],
149 positions_double[new_tets[i][2]][1],
150 positions_double[new_tets[i][2]][2],
151 positions_double[new_tets[i][3]][0],
152 positions_double[new_tets[i][3]][1],
153 positions_double[new_tets[i][3]][2],
154 }});
155
156
157 if (energy > new_energy_max) new_energy_max = energy;
158 }
159 }
160
161 return old_energy_max > new_energy_max * m_eps;
162}
163} // namespace wmtk
wmtk::Mesh::create_const_accessor
const attribute::Accessor< T, Mesh, D > create_const_accessor(const attribute::MeshAttributeHandle &handle) const
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::Swap56EnergyBeforeInvariant::before
bool before(const simplex::Simplex &t) const override
Definition Swap56EnergyBeforeInvariant.cpp:19
wmtk::Swap56EnergyBeforeInvariant::m_coordinate_handle
const attribute::TypedAttributeHandle< Rational > m_coordinate_handle
Definition Swap56EnergyBeforeInvariant.hpp:20
wmtk::Swap56EnergyBeforeInvariant::Swap56EnergyBeforeInvariant
Swap56EnergyBeforeInvariant(const Mesh &m, const attribute::TypedAttributeHandle< Rational > &coordinate, int64_t collapse_index, double eps=1.0001)
Definition Swap56EnergyBeforeInvariant.cpp:8
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::Invariant::mesh
const Mesh & mesh() const
Definition Invariant.cpp:35
wmtk::simplex::Simplex::tuple
const Tuple & tuple() const
Definition Simplex.hpp:53
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