Wildmeshing Toolkit
Swap32EnergyBeforeInvariantDouble.cpp
Go to the documentation of this file.
1 #include "Swap32EnergyBeforeInvariantDouble.hpp"
2 #include <wmtk/Mesh.hpp>
3 #include <wmtk/function/utils/amips.hpp>
4 #include <wmtk/utils/orient.hpp>
5 
6 namespace wmtk {
7 Swap32EnergyBeforeInvariantDouble::Swap32EnergyBeforeInvariantDouble(
8  const Mesh& m,
9  const attribute::TypedAttributeHandle<double>& coordinate,
10  double eps)
11  : Invariant(m, true, false, false)
12  , m_coordinate_handle(coordinate)
13  , m_eps(eps)
14 {}
15 
16 bool Swap32EnergyBeforeInvariantDouble::before(const simplex::Simplex& t) const
17 {
18  constexpr static PrimitiveType PV = PrimitiveType::Vertex;
19  constexpr static PrimitiveType PE = PrimitiveType::Edge;
20  constexpr static PrimitiveType PF = PrimitiveType::Triangle;
21  constexpr static PrimitiveType PT = PrimitiveType::Tetrahedron;
22 
23  auto accessor = mesh().create_const_accessor(m_coordinate_handle);
24 
25  // get the coords of the vertices
26  // input edge end points
27  const Tuple v0 = t.tuple();
28  const Tuple v1 = mesh().switch_tuple(v0, PV);
29  // other 3 vertices
30  const Tuple v2 = mesh().switch_tuples(v0, {PE, PV});
31  const Tuple v3 = mesh().switch_tuples(v0, {PF, PE, PV});
32  const Tuple v4 = mesh().switch_tuples(v0, {PT, PF, PE, PV});
33 
34 
35  std::array<Eigen::Vector3<double>, 5> positions = {
36  {accessor.const_vector_attribute(v0),
37  accessor.const_vector_attribute(v1),
38  accessor.const_vector_attribute(v2),
39  accessor.const_vector_attribute(v3),
40  accessor.const_vector_attribute(v4)}};
41  std::array<Eigen::Vector3d, 5> positions_double = {
42  {accessor.const_vector_attribute(v0),
43  accessor.const_vector_attribute(v1),
44  accessor.const_vector_attribute(v2),
45  accessor.const_vector_attribute(v3),
46  accessor.const_vector_attribute(v4)}};
47 
48  std::array<std::array<int, 4>, 3> old_tets = {{{{0, 1, 2, 3}}, {{0, 1, 3, 4}}, {{0, 1, 4, 2}}}};
49  std::array<std::array<int, 4>, 2> new_tets = {{{{2, 3, 4, 0}}, {{2, 4, 3, 1}}}};
50 
51 
52  double old_energy_max = std::numeric_limits<double>::lowest();
53  double new_energy_max = std::numeric_limits<double>::lowest();
54 
55  for (int i = 0; i < 3; ++i) {
56  if (utils::wmtk_orient3d(
57  positions[old_tets[i][0]],
58  positions[old_tets[i][1]],
59  positions[old_tets[i][2]],
60  positions[old_tets[i][3]]) > 0) {
61  auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
62  positions_double[old_tets[i][0]][0],
63  positions_double[old_tets[i][0]][1],
64  positions_double[old_tets[i][0]][2],
65  positions_double[old_tets[i][1]][0],
66  positions_double[old_tets[i][1]][1],
67  positions_double[old_tets[i][1]][2],
68  positions_double[old_tets[i][2]][0],
69  positions_double[old_tets[i][2]][1],
70  positions_double[old_tets[i][2]][2],
71  positions_double[old_tets[i][3]][0],
72  positions_double[old_tets[i][3]][1],
73  positions_double[old_tets[i][3]][2],
74  }});
75 
76 
77  if (energy > old_energy_max) old_energy_max = energy;
78  } else {
79  auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
80  positions_double[old_tets[i][1]][0],
81  positions_double[old_tets[i][1]][1],
82  positions_double[old_tets[i][1]][2],
83  positions_double[old_tets[i][0]][0],
84  positions_double[old_tets[i][0]][1],
85  positions_double[old_tets[i][0]][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  }
97  }
98 
99  for (int i = 0; i < 2; ++i) {
100  if (utils::wmtk_orient3d(
101  positions[new_tets[i][0]],
102  positions[new_tets[i][1]],
103  positions[new_tets[i][2]],
104  positions[new_tets[i][3]]) > 0) {
105  auto energy = wmtk::function::utils::Tet_AMIPS_energy({{
106  positions_double[new_tets[i][0]][0],
107  positions_double[new_tets[i][0]][1],
108  positions_double[new_tets[i][0]][2],
109  positions_double[new_tets[i][1]][0],
110  positions_double[new_tets[i][1]][1],
111  positions_double[new_tets[i][1]][2],
112  positions_double[new_tets[i][2]][0],
113  positions_double[new_tets[i][2]][1],
114  positions_double[new_tets[i][2]][2],
115  positions_double[new_tets[i][3]][0],
116  positions_double[new_tets[i][3]][1],
117  positions_double[new_tets[i][3]][2],
118  }});
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 
139  if (energy > new_energy_max) new_energy_max = energy;
140  }
141  }
142 
143  return old_energy_max > new_energy_max * m_eps;
144 }
145 
146 } // 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:967
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::Swap32EnergyBeforeInvariantDouble::before
bool before(const simplex::Simplex &t) const override
Definition: Swap32EnergyBeforeInvariantDouble.cpp:16
wmtk::Swap32EnergyBeforeInvariantDouble::m_coordinate_handle
const attribute::TypedAttributeHandle< double > m_coordinate_handle
Definition: Swap32EnergyBeforeInvariantDouble.hpp:19
wmtk::Swap32EnergyBeforeInvariantDouble::Swap32EnergyBeforeInvariantDouble
Swap32EnergyBeforeInvariantDouble(const Mesh &m, const attribute::TypedAttributeHandle< double > &coordinate, double eps=1.0001)
Definition: Swap32EnergyBeforeInvariantDouble.cpp:7
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::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
Definition: Accessor.hpp:6
wmtk::PV
constexpr PrimitiveType PV
Definition: EnvelopeInvariant.cpp:24
wmtk::PE
constexpr PrimitiveType PE
Definition: EnvelopeInvariant.cpp:25