Wildmeshing Toolkit
Loading...
Searching...
No Matches
Swap44EnergyBeforeInvariant.cpp
Go to the documentation of this file.
1#include "Swap44EnergyBeforeInvariant.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 {
8Swap44EnergyBeforeInvariant::Swap44EnergyBeforeInvariant(
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 Swap44EnergyBeforeInvariant::before(const simplex::Simplex& t) const
20{
21 assert(simplex::top_dimension_cofaces(mesh(), t).size() == 4);
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 // get the coords of the vertices
30 // input edge end points
31 const Tuple e0 = t.tuple();
32 const Tuple e1 = mesh().switch_tuple(e0, PV);
33 // other four vertices
34 std::array<Tuple, 4> v;
35 auto iter_tuple = e0;
36 for (int64_t i = 0; i < 4; ++i) {
37 v[i] = mesh().switch_tuples(iter_tuple, {PE, PV});
38 iter_tuple = mesh().switch_tuples(iter_tuple, {PF, PT});
39 }
40 assert(iter_tuple == e0);
41
42 std::array<Eigen::Vector3<Rational>, 6> positions = {
43 {accessor.const_vector_attribute(v[(m_collapse_index + 0) % 4]),
44 accessor.const_vector_attribute(v[(m_collapse_index + 1) % 4]),
45 accessor.const_vector_attribute(v[(m_collapse_index + 2) % 4]),
46 accessor.const_vector_attribute(v[(m_collapse_index + 3) % 4]),
47 accessor.const_vector_attribute(e0),
48 accessor.const_vector_attribute(e1)}};
49 std::array<Eigen::Vector3d, 6> positions_double = {
50 {positions[0].cast<double>(),
51 positions[1].cast<double>(),
52 positions[2].cast<double>(),
53 positions[3].cast<double>(),
54 positions[4].cast<double>(),
55 positions[5].cast<double>()}};
56
57 std::array<std::array<int, 4>, 4> old_tets = {
58 {{{0, 1, 4, 5}}, {{1, 2, 4, 5}}, {{2, 3, 4, 5}}, {{3, 0, 4, 5}}}};
59 std::array<std::array<int, 4>, 4> new_tets = {
60 {{{0, 1, 2, 4}}, {{0, 2, 3, 4}}, {{0, 1, 2, 5}}, {{0, 2, 3, 5}}}};
61
62 double old_energy_max = std::numeric_limits<double>::lowest();
63 double new_energy_max = std::numeric_limits<double>::lowest();
64
65 for (int i = 0; i < 4; ++i) {
66 if (utils::wmtk_orient3d(
67 positions[old_tets[i][0]],
68 positions[old_tets[i][1]],
69 positions[old_tets[i][2]],
70 positions[old_tets[i][3]]) > 0) {
71 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
72 positions_double[old_tets[i][0]][0],
73 positions_double[old_tets[i][0]][1],
74 positions_double[old_tets[i][0]][2],
75 positions_double[old_tets[i][1]][0],
76 positions_double[old_tets[i][1]][1],
77 positions_double[old_tets[i][1]][2],
78 positions_double[old_tets[i][2]][0],
79 positions_double[old_tets[i][2]][1],
80 positions_double[old_tets[i][2]][2],
81 positions_double[old_tets[i][3]][0],
82 positions_double[old_tets[i][3]][1],
83 positions_double[old_tets[i][3]][2],
84 }});
85
86 old_energy_max = std::max(energy, old_energy_max);
87 } else {
88 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
89 positions_double[old_tets[i][1]][0],
90 positions_double[old_tets[i][1]][1],
91 positions_double[old_tets[i][1]][2],
92 positions_double[old_tets[i][0]][0],
93 positions_double[old_tets[i][0]][1],
94 positions_double[old_tets[i][0]][2],
95 positions_double[old_tets[i][2]][0],
96 positions_double[old_tets[i][2]][1],
97 positions_double[old_tets[i][2]][2],
98 positions_double[old_tets[i][3]][0],
99 positions_double[old_tets[i][3]][1],
100 positions_double[old_tets[i][3]][2],
101 }});
102
103 old_energy_max = std::max(energy, old_energy_max);
104 }
105
106 if (utils::wmtk_orient3d(
107 positions[new_tets[i][0]],
108 positions[new_tets[i][1]],
109 positions[new_tets[i][2]],
110 positions[new_tets[i][3]]) > 0) {
111 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
112 positions_double[new_tets[i][0]][0],
113 positions_double[new_tets[i][0]][1],
114 positions_double[new_tets[i][0]][2],
115 positions_double[new_tets[i][1]][0],
116 positions_double[new_tets[i][1]][1],
117 positions_double[new_tets[i][1]][2],
118 positions_double[new_tets[i][2]][0],
119 positions_double[new_tets[i][2]][1],
120 positions_double[new_tets[i][2]][2],
121 positions_double[new_tets[i][3]][0],
122 positions_double[new_tets[i][3]][1],
123 positions_double[new_tets[i][3]][2],
124 }});
125
126 new_energy_max = std::max(energy, new_energy_max);
127 } else {
128 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
129 positions_double[new_tets[i][1]][0],
130 positions_double[new_tets[i][1]][1],
131 positions_double[new_tets[i][1]][2],
132 positions_double[new_tets[i][0]][0],
133 positions_double[new_tets[i][0]][1],
134 positions_double[new_tets[i][0]][2],
135 positions_double[new_tets[i][2]][0],
136 positions_double[new_tets[i][2]][1],
137 positions_double[new_tets[i][2]][2],
138 positions_double[new_tets[i][3]][0],
139 positions_double[new_tets[i][3]][1],
140 positions_double[new_tets[i][3]][2],
141 }});
142
143 new_energy_max = std::max(energy, new_energy_max);
144 }
145 }
146
147 return old_energy_max > new_energy_max * m_eps;
148}
149
150} // 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::Swap44EnergyBeforeInvariant::m_coordinate_handle
const attribute::TypedAttributeHandle< Rational > m_coordinate_handle
Definition Swap44EnergyBeforeInvariant.hpp:20
wmtk::Swap44EnergyBeforeInvariant::Swap44EnergyBeforeInvariant
Swap44EnergyBeforeInvariant(const Mesh &m, const attribute::TypedAttributeHandle< Rational > &coordinate, int64_t collapse_index, double eps=1.0001)
Definition Swap44EnergyBeforeInvariant.cpp:8
wmtk::Swap44EnergyBeforeInvariant::before
bool before(const simplex::Simplex &t) const override
Definition Swap44EnergyBeforeInvariant.cpp:19
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