Wildmeshing Toolkit
amips.cpp
Go to the documentation of this file.
1 #include "amips.hpp"
2 #include <Eigen/Dense>
3 #include <array>
4 #include <wmtk/utils/triangle_areas.hpp>
5 namespace wmtk::function::utils {
6 namespace detail {
7 namespace {
8 auto make_amips_target_triangle()
9 {
10  const static std::array<std::array<double, 2>, 3> m_target_triangle{{
11  // comments to keep formatting
12  std::array<double, 2>{{0., 0.}}, //
13  std::array<double, 2>{{1., 0.}}, //
14  std::array<double, 2>{{0.5, sqrt(3) / 2.}},
15  //
16  }};
17  auto map = Eigen::Matrix<double, 2, 3>::ConstMapType(m_target_triangle[0].data());
18 
19 
20 #if !defined(NDEBUG)
21  auto x = map.col(0);
22  auto y = map.col(1);
23  auto z = map.col(2);
24  assert(wmtk::utils::triangle_signed_2d_area(x, y, z) > 0);
25 #endif
26  return map;
27 }
28 
29 } // namespace
30 const Eigen::Matrix<double, 2, 3> amips_target_triangle = make_amips_target_triangle();
31 
32 namespace {
33 Eigen::Matrix2d make_amips_reference_to_barycentric()
34 {
35  const auto& A = amips_target_triangle;
36  Eigen::Matrix2d Ds = (A.rightCols<2>().colwise() - A.col(0));
37 
38  return Ds.inverse();
39 }
40 
41 } // namespace
42 const Eigen::Matrix2d amips_reference_to_barycentric = make_amips_reference_to_barycentric();
43 } // namespace detail
44 
45 double Tet_AMIPS_energy_aux(const std::array<double, 12>& T)
46 {
47  double helper_0[12];
48  helper_0[0] = T[0];
49  helper_0[1] = T[1];
50  helper_0[2] = T[2];
51  helper_0[3] = T[3];
52  helper_0[4] = T[4];
53  helper_0[5] = T[5];
54  helper_0[6] = T[6];
55  helper_0[7] = T[7];
56  helper_0[8] = T[8];
57  helper_0[9] = T[9];
58  helper_0[10] = T[10];
59  helper_0[11] = T[11];
60  double helper_1 = helper_0[2];
61  double helper_2 = helper_0[11];
62  double helper_3 = helper_0[0];
63  double helper_4 = helper_0[3];
64  double helper_5 = helper_0[9];
65  double helper_6 =
66  0.577350269189626 * helper_3 - 1.15470053837925 * helper_4 + 0.577350269189626 * helper_5;
67  double helper_7 = helper_0[1];
68  double helper_8 = helper_0[4];
69  double helper_9 = helper_0[7];
70  double helper_10 = helper_0[10];
71  double helper_11 = 0.408248290463863 * helper_10 + 0.408248290463863 * helper_7 +
72  0.408248290463863 * helper_8 - 1.22474487139159 * helper_9;
73  double helper_12 =
74  0.577350269189626 * helper_10 + 0.577350269189626 * helper_7 - 1.15470053837925 * helper_8;
75  double helper_13 = helper_0[6];
76  double helper_14 = -1.22474487139159 * helper_13 + 0.408248290463863 * helper_3 +
77  0.408248290463863 * helper_4 + 0.408248290463863 * helper_5;
78  double helper_15 = helper_0[5];
79  double helper_16 = helper_0[8];
80  double helper_17 = 0.408248290463863 * helper_1 + 0.408248290463863 * helper_15 -
81  1.22474487139159 * helper_16 + 0.408248290463863 * helper_2;
82  double helper_18 =
83  0.577350269189626 * helper_1 - 1.15470053837925 * helper_15 + 0.577350269189626 * helper_2;
84  double helper_19 = 0.5 * helper_13 + 0.5 * helper_4;
85  double helper_20 = 0.5 * helper_8 + 0.5 * helper_9;
86  double helper_21 = 0.5 * helper_15 + 0.5 * helper_16;
87  double helper_22 = (helper_1 - helper_2) * (helper_11 * helper_6 - helper_12 * helper_14) -
88  (-helper_10 + helper_7) * (-helper_14 * helper_18 + helper_17 * helper_6) +
89  (helper_3 - helper_5) * (-helper_11 * helper_18 + helper_12 * helper_17);
90  double res =
91  -(helper_1 * (-1.5 * helper_1 + 0.5 * helper_2 + helper_21) +
92  helper_10 * (-1.5 * helper_10 + helper_20 + 0.5 * helper_7) +
93  helper_13 * (-1.5 * helper_13 + 0.5 * helper_3 + 0.5 * helper_4 + 0.5 * helper_5) +
94  helper_15 * (0.5 * helper_1 - 1.5 * helper_15 + 0.5 * helper_16 + 0.5 * helper_2) +
95  helper_16 * (0.5 * helper_1 + 0.5 * helper_15 - 1.5 * helper_16 + 0.5 * helper_2) +
96  helper_2 * (0.5 * helper_1 - 1.5 * helper_2 + helper_21) +
97  helper_3 * (helper_19 - 1.5 * helper_3 + 0.5 * helper_5) +
98  helper_4 * (0.5 * helper_13 + 0.5 * helper_3 - 1.5 * helper_4 + 0.5 * helper_5) +
99  helper_5 * (helper_19 + 0.5 * helper_3 - 1.5 * helper_5) +
100  helper_7 * (0.5 * helper_10 + helper_20 - 1.5 * helper_7) +
101  helper_8 * (0.5 * helper_10 + 0.5 * helper_7 - 1.5 * helper_8 + 0.5 * helper_9) +
102  helper_9 * (0.5 * helper_10 + 0.5 * helper_7 + 0.5 * helper_8 - 1.5 * helper_9)) /
103  std::cbrt(helper_22 * helper_22);
104 
105  return res;
106 }
107 
108 void Tet_AMIPS_jacobian(const std::array<double, 12>& T, Eigen::Vector3d& result_0)
109 {
110  double helper_0[12];
111  helper_0[0] = T[0];
112  helper_0[1] = T[1];
113  helper_0[2] = T[2];
114  helper_0[3] = T[3];
115  helper_0[4] = T[4];
116  helper_0[5] = T[5];
117  helper_0[6] = T[6];
118  helper_0[7] = T[7];
119  helper_0[8] = T[8];
120  helper_0[9] = T[9];
121  helper_0[10] = T[10];
122  helper_0[11] = T[11];
123  double helper_1 = helper_0[1];
124  double helper_2 = helper_0[10];
125  double helper_3 = helper_1 - helper_2;
126  double helper_4 = helper_0[0];
127  double helper_5 = helper_0[3];
128  double helper_6 = helper_0[9];
129  double helper_7 =
130  0.577350269189626 * helper_4 - 1.15470053837925 * helper_5 + 0.577350269189626 * helper_6;
131  double helper_8 = helper_0[2];
132  double helper_9 = 0.408248290463863 * helper_8;
133  double helper_10 = helper_0[5];
134  double helper_11 = 0.408248290463863 * helper_10;
135  double helper_12 = helper_0[8];
136  double helper_13 = 1.22474487139159 * helper_12;
137  double helper_14 = helper_0[11];
138  double helper_15 = 0.408248290463863 * helper_14;
139  double helper_16 = helper_11 - helper_13 + helper_15 + helper_9;
140  double helper_17 = 0.577350269189626 * helper_8;
141  double helper_18 = 1.15470053837925 * helper_10;
142  double helper_19 = 0.577350269189626 * helper_14;
143  double helper_20 = helper_17 - helper_18 + helper_19;
144  double helper_21 = helper_0[6];
145  double helper_22 = -1.22474487139159 * helper_21 + 0.408248290463863 * helper_4 +
146  0.408248290463863 * helper_5 + 0.408248290463863 * helper_6;
147  double helper_23 = helper_16 * helper_7 - helper_20 * helper_22;
148  double helper_24 = -helper_14 + helper_8;
149  double helper_25 = 0.408248290463863 * helper_1;
150  double helper_26 = helper_0[4];
151  double helper_27 = 0.408248290463863 * helper_26;
152  double helper_28 = helper_0[7];
153  double helper_29 = 1.22474487139159 * helper_28;
154  double helper_30 = 0.408248290463863 * helper_2;
155  double helper_31 = helper_25 + helper_27 - helper_29 + helper_30;
156  double helper_32 = helper_31 * helper_7;
157  double helper_33 = 0.577350269189626 * helper_1;
158  double helper_34 = 1.15470053837925 * helper_26;
159  double helper_35 = 0.577350269189626 * helper_2;
160  double helper_36 = helper_33 - helper_34 + helper_35;
161  double helper_37 = helper_22 * helper_36;
162  double helper_38 = helper_4 - helper_6;
163  double helper_39 = helper_23 * helper_3 - helper_24 * (helper_32 - helper_37) -
164  helper_38 * (helper_16 * helper_36 - helper_20 * helper_31);
165  double helper_40 = pow(pow(helper_39, 2), -0.333333333333333);
166  double helper_41 = 0.707106781186548 * helper_10 - 0.707106781186548 * helper_12;
167  double helper_42 = 0.707106781186548 * helper_26 - 0.707106781186548 * helper_28;
168  double helper_43 = 0.5 * helper_21 + 0.5 * helper_5;
169  double helper_44 = 0.5 * helper_26 + 0.5 * helper_28;
170  double helper_45 = 0.5 * helper_10 + 0.5 * helper_12;
171  double helper_46 =
172  0.666666666666667 *
173  (helper_1 * (-1.5 * helper_1 + 0.5 * helper_2 + helper_44) +
174  helper_10 * (-1.5 * helper_10 + 0.5 * helper_12 + 0.5 * helper_14 + 0.5 * helper_8) +
175  helper_12 * (0.5 * helper_10 - 1.5 * helper_12 + 0.5 * helper_14 + 0.5 * helper_8) +
176  helper_14 * (-1.5 * helper_14 + helper_45 + 0.5 * helper_8) +
177  helper_2 * (0.5 * helper_1 - 1.5 * helper_2 + helper_44) +
178  helper_21 * (-1.5 * helper_21 + 0.5 * helper_4 + 0.5 * helper_5 + 0.5 * helper_6) +
179  helper_26 * (0.5 * helper_1 + 0.5 * helper_2 - 1.5 * helper_26 + 0.5 * helper_28) +
180  helper_28 * (0.5 * helper_1 + 0.5 * helper_2 + 0.5 * helper_26 - 1.5 * helper_28) +
181  helper_4 * (-1.5 * helper_4 + helper_43 + 0.5 * helper_6) +
182  helper_5 * (0.5 * helper_21 + 0.5 * helper_4 - 1.5 * helper_5 + 0.5 * helper_6) +
183  helper_6 * (0.5 * helper_4 + helper_43 - 1.5 * helper_6) +
184  helper_8 * (0.5 * helper_14 + helper_45 - 1.5 * helper_8)) /
185  helper_39;
186  double helper_47 = -0.707106781186548 * helper_21 + 0.707106781186548 * helper_5;
187  result_0[0] = -helper_40 * (1.0 * helper_21 - 3.0 * helper_4 +
188  helper_46 * (helper_41 * (-helper_1 + helper_2) -
189  helper_42 * (helper_14 - helper_8) -
190  (-helper_17 + helper_18 - helper_19) *
191  (-helper_25 - helper_27 + helper_29 - helper_30) +
192  (-helper_33 + helper_34 - helper_35) *
193  (-helper_11 + helper_13 - helper_15 - helper_9)) +
194  1.0 * helper_5 + 1.0 * helper_6);
195  result_0[1] =
196  helper_40 * (3.0 * helper_1 - 1.0 * helper_2 - 1.0 * helper_26 - 1.0 * helper_28 +
197  helper_46 * (helper_23 + helper_24 * helper_47 - helper_38 * helper_41));
198  result_0[2] =
199  helper_40 *
200  (-1.0 * helper_10 - 1.0 * helper_12 - 1.0 * helper_14 +
201  helper_46 * (-helper_3 * helper_47 - helper_32 + helper_37 + helper_38 * helper_42) +
202  3.0 * helper_8);
203 }
204 
205 void Tet_AMIPS_hessian(const std::array<double, 12>& T, Eigen::Matrix3d& result_0)
206 {
207  double helper_0[12];
208  helper_0[0] = T[0];
209  helper_0[1] = T[1];
210  helper_0[2] = T[2];
211  helper_0[3] = T[3];
212  helper_0[4] = T[4];
213  helper_0[5] = T[5];
214  helper_0[6] = T[6];
215  helper_0[7] = T[7];
216  helper_0[8] = T[8];
217  helper_0[9] = T[9];
218  helper_0[10] = T[10];
219  helper_0[11] = T[11];
220  double helper_1 = helper_0[2];
221  double helper_2 = helper_0[11];
222  double helper_3 = helper_1 - helper_2;
223  double helper_4 = helper_0[0];
224  double helper_5 = 0.577350269189626 * helper_4;
225  double helper_6 = helper_0[3];
226  double helper_7 = 1.15470053837925 * helper_6;
227  double helper_8 = helper_0[9];
228  double helper_9 = 0.577350269189626 * helper_8;
229  double helper_10 = helper_5 - helper_7 + helper_9;
230  double helper_11 = helper_0[1];
231  double helper_12 = 0.408248290463863 * helper_11;
232  double helper_13 = helper_0[4];
233  double helper_14 = 0.408248290463863 * helper_13;
234  double helper_15 = helper_0[7];
235  double helper_16 = 1.22474487139159 * helper_15;
236  double helper_17 = helper_0[10];
237  double helper_18 = 0.408248290463863 * helper_17;
238  double helper_19 = helper_12 + helper_14 - helper_16 + helper_18;
239  double helper_20 = helper_10 * helper_19;
240  double helper_21 = 0.577350269189626 * helper_11;
241  double helper_22 = 1.15470053837925 * helper_13;
242  double helper_23 = 0.577350269189626 * helper_17;
243  double helper_24 = helper_21 - helper_22 + helper_23;
244  double helper_25 = 0.408248290463863 * helper_4;
245  double helper_26 = 0.408248290463863 * helper_6;
246  double helper_27 = helper_0[6];
247  double helper_28 = 1.22474487139159 * helper_27;
248  double helper_29 = 0.408248290463863 * helper_8;
249  double helper_30 = helper_25 + helper_26 - helper_28 + helper_29;
250  double helper_31 = helper_24 * helper_30;
251  double helper_32 = helper_3 * (helper_20 - helper_31);
252  double helper_33 = helper_4 - helper_8;
253  double helper_34 = 0.408248290463863 * helper_1;
254  double helper_35 = helper_0[5];
255  double helper_36 = 0.408248290463863 * helper_35;
256  double helper_37 = helper_0[8];
257  double helper_38 = 1.22474487139159 * helper_37;
258  double helper_39 = 0.408248290463863 * helper_2;
259  double helper_40 = helper_34 + helper_36 - helper_38 + helper_39;
260  double helper_41 = helper_24 * helper_40;
261  double helper_42 = 0.577350269189626 * helper_1;
262  double helper_43 = 1.15470053837925 * helper_35;
263  double helper_44 = 0.577350269189626 * helper_2;
264  double helper_45 = helper_42 - helper_43 + helper_44;
265  double helper_46 = helper_19 * helper_45;
266  double helper_47 = helper_41 - helper_46;
267  double helper_48 = helper_33 * helper_47;
268  double helper_49 = helper_11 - helper_17;
269  double helper_50 = helper_10 * helper_40;
270  double helper_51 = helper_30 * helper_45;
271  double helper_52 = helper_50 - helper_51;
272  double helper_53 = helper_49 * helper_52;
273  double helper_54 = helper_32 + helper_48 - helper_53;
274  double helper_55 = pow(helper_54, 2);
275  double helper_56 = pow(helper_55, -0.333333333333333);
276  double helper_57 = 1.0 * helper_27 - 3.0 * helper_4 + 1.0 * helper_6 + 1.0 * helper_8;
277  double helper_58 = 0.707106781186548 * helper_13;
278  double helper_59 = 0.707106781186548 * helper_15;
279  double helper_60 = helper_58 - helper_59;
280  double helper_61 = helper_3 * helper_60;
281  double helper_62 = 0.707106781186548 * helper_35 - 0.707106781186548 * helper_37;
282  double helper_63 = helper_49 * helper_62;
283  double helper_64 = helper_47 + helper_61 - helper_63;
284  double helper_65 = 1.33333333333333 / helper_54;
285  double helper_66 = 1.0 / helper_55;
286  double helper_67 = 0.5 * helper_27 + 0.5 * helper_6;
287  double helper_68 = -1.5 * helper_4 + helper_67 + 0.5 * helper_8;
288  double helper_69 = 0.5 * helper_4 + helper_67 - 1.5 * helper_8;
289  double helper_70 = -1.5 * helper_27 + 0.5 * helper_4 + 0.5 * helper_6 + 0.5 * helper_8;
290  double helper_71 = 0.5 * helper_27 + 0.5 * helper_4 - 1.5 * helper_6 + 0.5 * helper_8;
291  double helper_72 = 0.5 * helper_13 + 0.5 * helper_15;
292  double helper_73 = -1.5 * helper_11 + 0.5 * helper_17 + helper_72;
293  double helper_74 = 0.5 * helper_11 - 1.5 * helper_17 + helper_72;
294  double helper_75 = 0.5 * helper_11 + 0.5 * helper_13 - 1.5 * helper_15 + 0.5 * helper_17;
295  double helper_76 = 0.5 * helper_11 - 1.5 * helper_13 + 0.5 * helper_15 + 0.5 * helper_17;
296  double helper_77 = 0.5 * helper_35 + 0.5 * helper_37;
297  double helper_78 = -1.5 * helper_1 + 0.5 * helper_2 + helper_77;
298  double helper_79 = 0.5 * helper_1 - 1.5 * helper_2 + helper_77;
299  double helper_80 = 0.5 * helper_1 + 0.5 * helper_2 + 0.5 * helper_35 - 1.5 * helper_37;
300  double helper_81 = 0.5 * helper_1 + 0.5 * helper_2 - 1.5 * helper_35 + 0.5 * helper_37;
301  double helper_82 = helper_1 * helper_78 + helper_11 * helper_73 + helper_13 * helper_76 +
302  helper_15 * helper_75 + helper_17 * helper_74 + helper_2 * helper_79 +
303  helper_27 * helper_70 + helper_35 * helper_81 + helper_37 * helper_80 +
304  helper_4 * helper_68 + helper_6 * helper_71 + helper_69 * helper_8;
305  double helper_83 = 0.444444444444444 * helper_66 * helper_82;
306  double helper_84 = helper_66 * helper_82;
307  double helper_85 = -helper_32 - helper_48 + helper_53;
308  double helper_86 = 1.0 / helper_85;
309  double helper_87 = helper_86 * pow(pow(helper_85, 2), -0.333333333333333);
310  double helper_88 = 0.707106781186548 * helper_6;
311  double helper_89 = 0.707106781186548 * helper_27;
312  double helper_90 = helper_88 - helper_89;
313  double helper_91 =
314  0.666666666666667 * helper_10 * helper_40 + 0.666666666666667 * helper_3 * helper_90 -
315  0.666666666666667 * helper_30 * helper_45 - 0.666666666666667 * helper_33 * helper_62;
316  double helper_92 = -3.0 * helper_11 + 1.0 * helper_13 + 1.0 * helper_15 + 1.0 * helper_17;
317  double helper_93 = -helper_11 + helper_17;
318  double helper_94 = -helper_1 + helper_2;
319  double helper_95 = -helper_21 + helper_22 - helper_23;
320  double helper_96 = -helper_34 - helper_36 + helper_38 - helper_39;
321  double helper_97 = -helper_42 + helper_43 - helper_44;
322  double helper_98 = -helper_12 - helper_14 + helper_16 - helper_18;
323  double helper_99 =
324  -0.666666666666667 * helper_60 * helper_94 + 0.666666666666667 * helper_62 * helper_93 +
325  0.666666666666667 * helper_95 * helper_96 - 0.666666666666667 * helper_97 * helper_98;
326  double helper_100 = helper_3 * helper_90;
327  double helper_101 = helper_33 * helper_62;
328  double helper_102 = helper_100 - helper_101 + helper_52;
329  double helper_103 = -helper_60 * helper_94 + helper_62 * helper_93 + helper_95 * helper_96 -
330  helper_97 * helper_98;
331  double helper_104 = 0.444444444444444 * helper_102 * helper_103 * helper_82 * helper_86 +
332  helper_57 * helper_91 - helper_92 * helper_99;
333  double helper_105 =
334  1.85037170770859e-17 * helper_1 * helper_78 + 1.85037170770859e-17 * helper_11 * helper_73 +
335  1.85037170770859e-17 * helper_13 * helper_76 +
336  1.85037170770859e-17 * helper_15 * helper_75 +
337  1.85037170770859e-17 * helper_17 * helper_74 + 1.85037170770859e-17 * helper_2 * helper_79 +
338  1.85037170770859e-17 * helper_27 * helper_70 +
339  1.85037170770859e-17 * helper_35 * helper_81 +
340  1.85037170770859e-17 * helper_37 * helper_80 + 1.85037170770859e-17 * helper_4 * helper_68 +
341  1.85037170770859e-17 * helper_6 * helper_71 + 1.85037170770859e-17 * helper_69 * helper_8;
342  double helper_106 = helper_64 * helper_82 * helper_86;
343  double helper_107 =
344  -0.666666666666667 * helper_10 * helper_19 + 0.666666666666667 * helper_24 * helper_30 +
345  0.666666666666667 * helper_33 * helper_60 - 0.666666666666667 * helper_49 * helper_90;
346  double helper_108 = -3.0 * helper_1 + 1.0 * helper_2 + 1.0 * helper_35 + 1.0 * helper_37;
347  double helper_109 = -helper_20 + helper_31 + helper_33 * helper_60 - helper_49 * helper_90;
348  double helper_110 = 0.444444444444444 * helper_109 * helper_82 * helper_86;
349  double helper_111 = helper_103 * helper_110 + helper_107 * helper_57 - helper_108 * helper_99;
350  double helper_112 = -helper_4 + helper_8;
351  double helper_113 = -helper_88 + helper_89;
352  double helper_114 = -helper_5 + helper_7 - helper_9;
353  double helper_115 = -helper_25 - helper_26 + helper_28 - helper_29;
354  double helper_116 = helper_82 * helper_86 *
355  (helper_112 * helper_62 + helper_113 * helper_94 + helper_114 * helper_96 -
356  helper_115 * helper_97);
357  double helper_117 = -helper_100 + helper_101 - helper_50 + helper_51;
358  double helper_118 = -helper_102 * helper_110 + helper_107 * helper_92 + helper_108 * helper_91;
359  double helper_119 = helper_82 * helper_86 *
360  (helper_112 * (-helper_58 + helper_59) - helper_113 * helper_93 -
361  helper_114 * helper_98 + helper_115 * helper_95);
362  result_0(0, 0) =
363  helper_56 * (helper_57 * helper_64 * helper_65 - pow(helper_64, 2) * helper_83 +
364  0.666666666666667 * helper_64 * helper_84 *
365  (-helper_41 + helper_46 - helper_61 + helper_63) +
366  3.0);
367  result_0(0, 1) = helper_87 * (helper_104 - helper_105 * helper_35 + helper_106 * helper_91);
368  result_0(0, 2) = helper_87 * (helper_106 * helper_107 + helper_111);
369  result_0(1, 0) = helper_87 * (helper_104 + helper_116 * helper_99);
370  result_0(1, 1) =
371  helper_56 * (-pow(helper_117, 2) * helper_83 + helper_117 * helper_65 * helper_92 +
372  helper_117 * helper_84 * helper_91 + 3.0);
373  result_0(1, 2) = helper_87 * (-helper_105 * helper_6 - helper_107 * helper_116 + helper_118);
374  result_0(2, 0) = helper_87 * (-helper_105 * helper_13 + helper_111 + helper_119 * helper_99);
375  result_0(2, 1) = helper_87 * (helper_118 - helper_119 * helper_91);
376  result_0(2, 2) = helper_56 * (-helper_108 * helper_109 * helper_65 -
377  1.11111111111111 * pow(helper_109, 2) * helper_84 + 3.0);
378 }
379 
380 double Tet_AMIPS_energy(const std::array<double, 12>& T)
381 {
382  double res = Tet_AMIPS_energy_aux(T);
383 
384  // Maybe use tet_is_energy_unstable
385  if (res > 1e8) {
386  std::array<Rational, 12> r_T;
387  for (int j = 0; j < 12; j++) {
388  r_T[j] = T[j];
389  }
390 
391  const Rational twothird = Rational(2) / Rational(3);
392 
393  Rational tmp =
394  ((-r_T[1 + 2] + r_T[1 + 5]) * r_T[1 + 1] + r_T[1 + 2] * r_T[1 + 7] +
395  (r_T[1 + -1] - r_T[1 + 5]) * r_T[1 + 4] - r_T[1 + -1] * r_T[1 + 7]) *
396  r_T[1 + 9] +
397  ((r_T[1 + 2] - r_T[1 + 5]) * r_T[1 + 0] - r_T[1 + 2] * r_T[1 + 6] +
398  (-r_T[1 + -1] + r_T[1 + 5]) * r_T[1 + 3] + r_T[1 + -1] * r_T[1 + 6]) *
399  r_T[1 + 10] +
400  (-r_T[1 + 2] * r_T[1 + 7] + (-r_T[1 + 8] + r_T[1 + 5]) * r_T[1 + 4] +
401  r_T[1 + 8] * r_T[1 + 7]) *
402  r_T[1 + 0] +
403  (r_T[1 + 2] * r_T[1 + 6] + (r_T[1 + 8] - r_T[1 + 5]) * r_T[1 + 3] -
404  r_T[1 + 8] * r_T[1 + 6]) *
405  r_T[1 + 1] +
406  (r_T[1 + 3] * r_T[1 + 7] - r_T[1 + 4] * r_T[1 + 6]) * (r_T[1 + -1] - r_T[1 + 8]);
407 
408  if (tmp == 0) return std::numeric_limits<double>::infinity();
409 
410  auto res_r =
411  Rational(27) / 16 * pow(tmp, -2) *
412  pow(r_T[1 + 9] * r_T[1 + 9] +
413  (-twothird * r_T[1 + 0] - twothird * r_T[1 + 3] - twothird * r_T[1 + 6]) *
414  r_T[1 + 9] +
415  r_T[1 + 10] * r_T[1 + 10] +
416  (-twothird * r_T[1 + 1] - twothird * r_T[1 + 4] - twothird * r_T[1 + 7]) *
417  r_T[1 + 10] +
418  r_T[1 + 0] * r_T[1 + 0] +
419  (-twothird * r_T[1 + 3] - twothird * r_T[1 + 6]) * r_T[1 + 0] +
420  r_T[1 + 1] * r_T[1 + 1] +
421  (-twothird * r_T[1 + 4] - twothird * r_T[1 + 7]) * r_T[1 + 1] +
422  r_T[1 + 2] * r_T[1 + 2] +
423  (-twothird * r_T[1 + -1] - twothird * r_T[1 + 8] - twothird * r_T[1 + 5]) *
424  r_T[1 + 2] +
425  r_T[1 + 3] * r_T[1 + 3] - twothird * r_T[1 + 3] * r_T[1 + 6] +
426  r_T[1 + 4] * r_T[1 + 4] - twothird * r_T[1 + 4] * r_T[1 + 7] +
427  r_T[1 + 5] * r_T[1 + 5] +
428  (-twothird * r_T[1 + -1] - twothird * r_T[1 + 8]) * r_T[1 + 5] -
429  twothird * r_T[1 + -1] * r_T[1 + 8] + r_T[1 + -1] * r_T[1 + -1] +
430  r_T[1 + 8] * r_T[1 + 8] + r_T[1 + 6] * r_T[1 + 6] + r_T[1 + 7] * r_T[1 + 7],
431  3);
432 
433  return std::cbrt(res_r.to_double());
434  } else {
435  return res;
436  }
437 }
438 
439 } // namespace wmtk::function::utils
wmtk::function::utils::detail::amips_reference_to_barycentric
const Eigen::Matrix2d amips_reference_to_barycentric
Definition: amips.cpp:42
wmtk::function::utils::detail::amips_target_triangle
const Eigen::Matrix< double, 2, 3 > amips_target_triangle
Definition: amips.cpp:30
wmtk::function::utils::Tet_AMIPS_energy_aux
double Tet_AMIPS_energy_aux(const std::array< double, 12 > &T)
Definition: amips.cpp:45
wmtk::function::utils::Tet_AMIPS_jacobian
void Tet_AMIPS_jacobian(const std::array< double, 12 > &T, Eigen::Vector3d &result_0)
Definition: amips.cpp:108
wmtk::function::utils::Tet_AMIPS_energy
double Tet_AMIPS_energy(const std::array< double, 12 > &T)
Definition: amips.cpp:380
wmtk::function::utils::Tet_AMIPS_hessian
void Tet_AMIPS_hessian(const std::array< double, 12 > &T, Eigen::Matrix3d &result_0)
Definition: amips.cpp:205
wmtk::utils::triangle_signed_2d_area
auto triangle_signed_2d_area(const Eigen::MatrixBase< ADerived > &a, const Eigen::MatrixBase< BDerived > &b, const Eigen::MatrixBase< CDerived > &c) -> typename ADerived::Scalar
Definition: triangle_areas.hpp:22
wmtk::Vector3d
Vector< double, 3 > Vector3d
Definition: Types.hpp:39
wmtk::pow
Rational pow(const Rational &x, int p)
Definition: Rational.cpp:166
triangle_areas.hpp