Wildmeshing Toolkit
Loading...
Searching...
No Matches
Swap44EnergyBeforeInvariantDouble.cpp
Go to the documentation of this file.
1#include "Swap44EnergyBeforeInvariantDouble.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 {
8Swap44EnergyBeforeInvariantDouble::Swap44EnergyBeforeInvariantDouble(
9 const Mesh& m,
10 const attribute::TypedAttributeHandle<double>& 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 Swap44EnergyBeforeInvariantDouble::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<double>, 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], positions[1], positions[2], positions[3], positions[4], positions[5]}};
51
52 std::array<std::array<int, 4>, 4> old_tets = {
53 {{{0, 1, 4, 5}}, {{1, 2, 4, 5}}, {{2, 3, 4, 5}}, {{3, 0, 4, 5}}}};
54 std::array<std::array<int, 4>, 4> new_tets = {
55 {{{0, 1, 2, 4}}, {{0, 2, 3, 4}}, {{0, 1, 2, 5}}, {{0, 2, 3, 5}}}};
56
57 double old_energy_max = std::numeric_limits<double>::lowest();
58 double new_energy_max = std::numeric_limits<double>::lowest();
59
60 for (int i = 0; i < 4; ++i) {
61 if (utils::wmtk_orient3d(
62 positions[old_tets[i][0]],
63 positions[old_tets[i][1]],
64 positions[old_tets[i][2]],
65 positions[old_tets[i][3]]) > 0) {
66 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
67 positions_double[old_tets[i][0]][0],
68 positions_double[old_tets[i][0]][1],
69 positions_double[old_tets[i][0]][2],
70 positions_double[old_tets[i][1]][0],
71 positions_double[old_tets[i][1]][1],
72 positions_double[old_tets[i][1]][2],
73 positions_double[old_tets[i][2]][0],
74 positions_double[old_tets[i][2]][1],
75 positions_double[old_tets[i][2]][2],
76 positions_double[old_tets[i][3]][0],
77 positions_double[old_tets[i][3]][1],
78 positions_double[old_tets[i][3]][2],
79 }});
80
81 if (energy > old_energy_max) old_energy_max = energy;
82 } else {
83 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
84 positions_double[old_tets[i][1]][0],
85 positions_double[old_tets[i][1]][1],
86 positions_double[old_tets[i][1]][2],
87 positions_double[old_tets[i][0]][0],
88 positions_double[old_tets[i][0]][1],
89 positions_double[old_tets[i][0]][2],
90 positions_double[old_tets[i][2]][0],
91 positions_double[old_tets[i][2]][1],
92 positions_double[old_tets[i][2]][2],
93 positions_double[old_tets[i][3]][0],
94 positions_double[old_tets[i][3]][1],
95 positions_double[old_tets[i][3]][2],
96 }});
97
98 if (energy > old_energy_max) old_energy_max = energy;
99 }
100
101 if (utils::wmtk_orient3d(
102 positions[new_tets[i][0]],
103 positions[new_tets[i][1]],
104 positions[new_tets[i][2]],
105 positions[new_tets[i][3]]) > 0) {
106 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
107 positions_double[new_tets[i][0]][0],
108 positions_double[new_tets[i][0]][1],
109 positions_double[new_tets[i][0]][2],
110 positions_double[new_tets[i][1]][0],
111 positions_double[new_tets[i][1]][1],
112 positions_double[new_tets[i][1]][2],
113 positions_double[new_tets[i][2]][0],
114 positions_double[new_tets[i][2]][1],
115 positions_double[new_tets[i][2]][2],
116 positions_double[new_tets[i][3]][0],
117 positions_double[new_tets[i][3]][1],
118 positions_double[new_tets[i][3]][2],
119 }});
120
121 if (energy > new_energy_max) new_energy_max = energy;
122 } else {
123 auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
124 positions_double[new_tets[i][1]][0],
125 positions_double[new_tets[i][1]][1],
126 positions_double[new_tets[i][1]][2],
127 positions_double[new_tets[i][0]][0],
128 positions_double[new_tets[i][0]][1],
129 positions_double[new_tets[i][0]][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 if (energy > new_energy_max) new_energy_max = energy;
139 }
140 }
141
142 return old_energy_max > new_energy_max * m_eps;
143}
144
145} // 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::Swap44EnergyBeforeInvariantDouble::before
bool before(const simplex::Simplex &t) const override
Definition Swap44EnergyBeforeInvariantDouble.cpp:19
wmtk::Swap44EnergyBeforeInvariantDouble::Swap44EnergyBeforeInvariantDouble
Swap44EnergyBeforeInvariantDouble(const Mesh &m, const attribute::TypedAttributeHandle< double > &coordinate, int64_t collapse_index, double eps=1.0001)
Definition Swap44EnergyBeforeInvariantDouble.cpp:8
wmtk::Swap44EnergyBeforeInvariantDouble::m_coordinate_handle
const attribute::TypedAttributeHandle< double > m_coordinate_handle
Definition Swap44EnergyBeforeInvariantDouble.hpp:20
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