Wildmeshing Toolkit
AMIPS.cpp
Go to the documentation of this file.
1 #include "AMIPS.hpp"
2 
3 #include <wmtk/Mesh.hpp>
4 #include <wmtk/function/utils/SimplexGetter.hpp>
5 #include <wmtk/simplex/Simplex.hpp>
6 #include <wmtk/utils/Rational.hpp>
7 #include <wmtk/utils/orient.hpp>
8 
9 
10 namespace wmtk::function {
11 
12 double Tet_AMIPS_energy_aux(const std::array<double, 12>& T)
13 {
14  double helper_0[12];
15  helper_0[0] = T[0];
16  helper_0[1] = T[1];
17  helper_0[2] = T[2];
18  helper_0[3] = T[3];
19  helper_0[4] = T[4];
20  helper_0[5] = T[5];
21  helper_0[6] = T[6];
22  helper_0[7] = T[7];
23  helper_0[8] = T[8];
24  helper_0[9] = T[9];
25  helper_0[10] = T[10];
26  helper_0[11] = T[11];
27  double helper_1 = helper_0[2];
28  double helper_2 = helper_0[11];
29  double helper_3 = helper_0[0];
30  double helper_4 = helper_0[3];
31  double helper_5 = helper_0[9];
32  double helper_6 =
33  0.577350269189626 * helper_3 - 1.15470053837925 * helper_4 + 0.577350269189626 * helper_5;
34  double helper_7 = helper_0[1];
35  double helper_8 = helper_0[4];
36  double helper_9 = helper_0[7];
37  double helper_10 = helper_0[10];
38  double helper_11 = 0.408248290463863 * helper_10 + 0.408248290463863 * helper_7 +
39  0.408248290463863 * helper_8 - 1.22474487139159 * helper_9;
40  double helper_12 =
41  0.577350269189626 * helper_10 + 0.577350269189626 * helper_7 - 1.15470053837925 * helper_8;
42  double helper_13 = helper_0[6];
43  double helper_14 = -1.22474487139159 * helper_13 + 0.408248290463863 * helper_3 +
44  0.408248290463863 * helper_4 + 0.408248290463863 * helper_5;
45  double helper_15 = helper_0[5];
46  double helper_16 = helper_0[8];
47  double helper_17 = 0.408248290463863 * helper_1 + 0.408248290463863 * helper_15 -
48  1.22474487139159 * helper_16 + 0.408248290463863 * helper_2;
49  double helper_18 =
50  0.577350269189626 * helper_1 - 1.15470053837925 * helper_15 + 0.577350269189626 * helper_2;
51  double helper_19 = 0.5 * helper_13 + 0.5 * helper_4;
52  double helper_20 = 0.5 * helper_8 + 0.5 * helper_9;
53  double helper_21 = 0.5 * helper_15 + 0.5 * helper_16;
54  double helper_22 = (helper_1 - helper_2) * (helper_11 * helper_6 - helper_12 * helper_14) -
55  (-helper_10 + helper_7) * (-helper_14 * helper_18 + helper_17 * helper_6) +
56  (helper_3 - helper_5) * (-helper_11 * helper_18 + helper_12 * helper_17);
57  double res =
58  -(helper_1 * (-1.5 * helper_1 + 0.5 * helper_2 + helper_21) +
59  helper_10 * (-1.5 * helper_10 + helper_20 + 0.5 * helper_7) +
60  helper_13 * (-1.5 * helper_13 + 0.5 * helper_3 + 0.5 * helper_4 + 0.5 * helper_5) +
61  helper_15 * (0.5 * helper_1 - 1.5 * helper_15 + 0.5 * helper_16 + 0.5 * helper_2) +
62  helper_16 * (0.5 * helper_1 + 0.5 * helper_15 - 1.5 * helper_16 + 0.5 * helper_2) +
63  helper_2 * (0.5 * helper_1 - 1.5 * helper_2 + helper_21) +
64  helper_3 * (helper_19 - 1.5 * helper_3 + 0.5 * helper_5) +
65  helper_4 * (0.5 * helper_13 + 0.5 * helper_3 - 1.5 * helper_4 + 0.5 * helper_5) +
66  helper_5 * (helper_19 + 0.5 * helper_3 - 1.5 * helper_5) +
67  helper_7 * (0.5 * helper_10 + helper_20 - 1.5 * helper_7) +
68  helper_8 * (0.5 * helper_10 + 0.5 * helper_7 - 1.5 * helper_8 + 0.5 * helper_9) +
69  helper_9 * (0.5 * helper_10 + 0.5 * helper_7 + 0.5 * helper_8 - 1.5 * helper_9)) /
70  std::cbrt(helper_22 * helper_22);
71 
72  return res;
73 }
74 
75 void Tet_AMIPS_jacobian(const std::array<double, 12>& T, Eigen::Vector3d& result_0)
76 {
77  double helper_0[12];
78  helper_0[0] = T[0];
79  helper_0[1] = T[1];
80  helper_0[2] = T[2];
81  helper_0[3] = T[3];
82  helper_0[4] = T[4];
83  helper_0[5] = T[5];
84  helper_0[6] = T[6];
85  helper_0[7] = T[7];
86  helper_0[8] = T[8];
87  helper_0[9] = T[9];
88  helper_0[10] = T[10];
89  helper_0[11] = T[11];
90  double helper_1 = helper_0[1];
91  double helper_2 = helper_0[10];
92  double helper_3 = helper_1 - helper_2;
93  double helper_4 = helper_0[0];
94  double helper_5 = helper_0[3];
95  double helper_6 = helper_0[9];
96  double helper_7 =
97  0.577350269189626 * helper_4 - 1.15470053837925 * helper_5 + 0.577350269189626 * helper_6;
98  double helper_8 = helper_0[2];
99  double helper_9 = 0.408248290463863 * helper_8;
100  double helper_10 = helper_0[5];
101  double helper_11 = 0.408248290463863 * helper_10;
102  double helper_12 = helper_0[8];
103  double helper_13 = 1.22474487139159 * helper_12;
104  double helper_14 = helper_0[11];
105  double helper_15 = 0.408248290463863 * helper_14;
106  double helper_16 = helper_11 - helper_13 + helper_15 + helper_9;
107  double helper_17 = 0.577350269189626 * helper_8;
108  double helper_18 = 1.15470053837925 * helper_10;
109  double helper_19 = 0.577350269189626 * helper_14;
110  double helper_20 = helper_17 - helper_18 + helper_19;
111  double helper_21 = helper_0[6];
112  double helper_22 = -1.22474487139159 * helper_21 + 0.408248290463863 * helper_4 +
113  0.408248290463863 * helper_5 + 0.408248290463863 * helper_6;
114  double helper_23 = helper_16 * helper_7 - helper_20 * helper_22;
115  double helper_24 = -helper_14 + helper_8;
116  double helper_25 = 0.408248290463863 * helper_1;
117  double helper_26 = helper_0[4];
118  double helper_27 = 0.408248290463863 * helper_26;
119  double helper_28 = helper_0[7];
120  double helper_29 = 1.22474487139159 * helper_28;
121  double helper_30 = 0.408248290463863 * helper_2;
122  double helper_31 = helper_25 + helper_27 - helper_29 + helper_30;
123  double helper_32 = helper_31 * helper_7;
124  double helper_33 = 0.577350269189626 * helper_1;
125  double helper_34 = 1.15470053837925 * helper_26;
126  double helper_35 = 0.577350269189626 * helper_2;
127  double helper_36 = helper_33 - helper_34 + helper_35;
128  double helper_37 = helper_22 * helper_36;
129  double helper_38 = helper_4 - helper_6;
130  double helper_39 = helper_23 * helper_3 - helper_24 * (helper_32 - helper_37) -
131  helper_38 * (helper_16 * helper_36 - helper_20 * helper_31);
132  double helper_40 = pow(pow(helper_39, 2), -0.333333333333333);
133  double helper_41 = 0.707106781186548 * helper_10 - 0.707106781186548 * helper_12;
134  double helper_42 = 0.707106781186548 * helper_26 - 0.707106781186548 * helper_28;
135  double helper_43 = 0.5 * helper_21 + 0.5 * helper_5;
136  double helper_44 = 0.5 * helper_26 + 0.5 * helper_28;
137  double helper_45 = 0.5 * helper_10 + 0.5 * helper_12;
138  double helper_46 =
139  0.666666666666667 *
140  (helper_1 * (-1.5 * helper_1 + 0.5 * helper_2 + helper_44) +
141  helper_10 * (-1.5 * helper_10 + 0.5 * helper_12 + 0.5 * helper_14 + 0.5 * helper_8) +
142  helper_12 * (0.5 * helper_10 - 1.5 * helper_12 + 0.5 * helper_14 + 0.5 * helper_8) +
143  helper_14 * (-1.5 * helper_14 + helper_45 + 0.5 * helper_8) +
144  helper_2 * (0.5 * helper_1 - 1.5 * helper_2 + helper_44) +
145  helper_21 * (-1.5 * helper_21 + 0.5 * helper_4 + 0.5 * helper_5 + 0.5 * helper_6) +
146  helper_26 * (0.5 * helper_1 + 0.5 * helper_2 - 1.5 * helper_26 + 0.5 * helper_28) +
147  helper_28 * (0.5 * helper_1 + 0.5 * helper_2 + 0.5 * helper_26 - 1.5 * helper_28) +
148  helper_4 * (-1.5 * helper_4 + helper_43 + 0.5 * helper_6) +
149  helper_5 * (0.5 * helper_21 + 0.5 * helper_4 - 1.5 * helper_5 + 0.5 * helper_6) +
150  helper_6 * (0.5 * helper_4 + helper_43 - 1.5 * helper_6) +
151  helper_8 * (0.5 * helper_14 + helper_45 - 1.5 * helper_8)) /
152  helper_39;
153  double helper_47 = -0.707106781186548 * helper_21 + 0.707106781186548 * helper_5;
154  result_0[0] = -helper_40 * (1.0 * helper_21 - 3.0 * helper_4 +
155  helper_46 * (helper_41 * (-helper_1 + helper_2) -
156  helper_42 * (helper_14 - helper_8) -
157  (-helper_17 + helper_18 - helper_19) *
158  (-helper_25 - helper_27 + helper_29 - helper_30) +
159  (-helper_33 + helper_34 - helper_35) *
160  (-helper_11 + helper_13 - helper_15 - helper_9)) +
161  1.0 * helper_5 + 1.0 * helper_6);
162  result_0[1] =
163  helper_40 * (3.0 * helper_1 - 1.0 * helper_2 - 1.0 * helper_26 - 1.0 * helper_28 +
164  helper_46 * (helper_23 + helper_24 * helper_47 - helper_38 * helper_41));
165  result_0[2] =
166  helper_40 *
167  (-1.0 * helper_10 - 1.0 * helper_12 - 1.0 * helper_14 +
168  helper_46 * (-helper_3 * helper_47 - helper_32 + helper_37 + helper_38 * helper_42) +
169  3.0 * helper_8);
170 }
171 
172 void Tet_AMIPS_hessian(const std::array<double, 12>& T, Eigen::Matrix3d& result_0)
173 {
174  double helper_0[12];
175  helper_0[0] = T[0];
176  helper_0[1] = T[1];
177  helper_0[2] = T[2];
178  helper_0[3] = T[3];
179  helper_0[4] = T[4];
180  helper_0[5] = T[5];
181  helper_0[6] = T[6];
182  helper_0[7] = T[7];
183  helper_0[8] = T[8];
184  helper_0[9] = T[9];
185  helper_0[10] = T[10];
186  helper_0[11] = T[11];
187  double helper_1 = helper_0[2];
188  double helper_2 = helper_0[11];
189  double helper_3 = helper_1 - helper_2;
190  double helper_4 = helper_0[0];
191  double helper_5 = 0.577350269189626 * helper_4;
192  double helper_6 = helper_0[3];
193  double helper_7 = 1.15470053837925 * helper_6;
194  double helper_8 = helper_0[9];
195  double helper_9 = 0.577350269189626 * helper_8;
196  double helper_10 = helper_5 - helper_7 + helper_9;
197  double helper_11 = helper_0[1];
198  double helper_12 = 0.408248290463863 * helper_11;
199  double helper_13 = helper_0[4];
200  double helper_14 = 0.408248290463863 * helper_13;
201  double helper_15 = helper_0[7];
202  double helper_16 = 1.22474487139159 * helper_15;
203  double helper_17 = helper_0[10];
204  double helper_18 = 0.408248290463863 * helper_17;
205  double helper_19 = helper_12 + helper_14 - helper_16 + helper_18;
206  double helper_20 = helper_10 * helper_19;
207  double helper_21 = 0.577350269189626 * helper_11;
208  double helper_22 = 1.15470053837925 * helper_13;
209  double helper_23 = 0.577350269189626 * helper_17;
210  double helper_24 = helper_21 - helper_22 + helper_23;
211  double helper_25 = 0.408248290463863 * helper_4;
212  double helper_26 = 0.408248290463863 * helper_6;
213  double helper_27 = helper_0[6];
214  double helper_28 = 1.22474487139159 * helper_27;
215  double helper_29 = 0.408248290463863 * helper_8;
216  double helper_30 = helper_25 + helper_26 - helper_28 + helper_29;
217  double helper_31 = helper_24 * helper_30;
218  double helper_32 = helper_3 * (helper_20 - helper_31);
219  double helper_33 = helper_4 - helper_8;
220  double helper_34 = 0.408248290463863 * helper_1;
221  double helper_35 = helper_0[5];
222  double helper_36 = 0.408248290463863 * helper_35;
223  double helper_37 = helper_0[8];
224  double helper_38 = 1.22474487139159 * helper_37;
225  double helper_39 = 0.408248290463863 * helper_2;
226  double helper_40 = helper_34 + helper_36 - helper_38 + helper_39;
227  double helper_41 = helper_24 * helper_40;
228  double helper_42 = 0.577350269189626 * helper_1;
229  double helper_43 = 1.15470053837925 * helper_35;
230  double helper_44 = 0.577350269189626 * helper_2;
231  double helper_45 = helper_42 - helper_43 + helper_44;
232  double helper_46 = helper_19 * helper_45;
233  double helper_47 = helper_41 - helper_46;
234  double helper_48 = helper_33 * helper_47;
235  double helper_49 = helper_11 - helper_17;
236  double helper_50 = helper_10 * helper_40;
237  double helper_51 = helper_30 * helper_45;
238  double helper_52 = helper_50 - helper_51;
239  double helper_53 = helper_49 * helper_52;
240  double helper_54 = helper_32 + helper_48 - helper_53;
241  double helper_55 = pow(helper_54, 2);
242  double helper_56 = pow(helper_55, -0.333333333333333);
243  double helper_57 = 1.0 * helper_27 - 3.0 * helper_4 + 1.0 * helper_6 + 1.0 * helper_8;
244  double helper_58 = 0.707106781186548 * helper_13;
245  double helper_59 = 0.707106781186548 * helper_15;
246  double helper_60 = helper_58 - helper_59;
247  double helper_61 = helper_3 * helper_60;
248  double helper_62 = 0.707106781186548 * helper_35 - 0.707106781186548 * helper_37;
249  double helper_63 = helper_49 * helper_62;
250  double helper_64 = helper_47 + helper_61 - helper_63;
251  double helper_65 = 1.33333333333333 / helper_54;
252  double helper_66 = 1.0 / helper_55;
253  double helper_67 = 0.5 * helper_27 + 0.5 * helper_6;
254  double helper_68 = -1.5 * helper_4 + helper_67 + 0.5 * helper_8;
255  double helper_69 = 0.5 * helper_4 + helper_67 - 1.5 * helper_8;
256  double helper_70 = -1.5 * helper_27 + 0.5 * helper_4 + 0.5 * helper_6 + 0.5 * helper_8;
257  double helper_71 = 0.5 * helper_27 + 0.5 * helper_4 - 1.5 * helper_6 + 0.5 * helper_8;
258  double helper_72 = 0.5 * helper_13 + 0.5 * helper_15;
259  double helper_73 = -1.5 * helper_11 + 0.5 * helper_17 + helper_72;
260  double helper_74 = 0.5 * helper_11 - 1.5 * helper_17 + helper_72;
261  double helper_75 = 0.5 * helper_11 + 0.5 * helper_13 - 1.5 * helper_15 + 0.5 * helper_17;
262  double helper_76 = 0.5 * helper_11 - 1.5 * helper_13 + 0.5 * helper_15 + 0.5 * helper_17;
263  double helper_77 = 0.5 * helper_35 + 0.5 * helper_37;
264  double helper_78 = -1.5 * helper_1 + 0.5 * helper_2 + helper_77;
265  double helper_79 = 0.5 * helper_1 - 1.5 * helper_2 + helper_77;
266  double helper_80 = 0.5 * helper_1 + 0.5 * helper_2 + 0.5 * helper_35 - 1.5 * helper_37;
267  double helper_81 = 0.5 * helper_1 + 0.5 * helper_2 - 1.5 * helper_35 + 0.5 * helper_37;
268  double helper_82 = helper_1 * helper_78 + helper_11 * helper_73 + helper_13 * helper_76 +
269  helper_15 * helper_75 + helper_17 * helper_74 + helper_2 * helper_79 +
270  helper_27 * helper_70 + helper_35 * helper_81 + helper_37 * helper_80 +
271  helper_4 * helper_68 + helper_6 * helper_71 + helper_69 * helper_8;
272  double helper_83 = 0.444444444444444 * helper_66 * helper_82;
273  double helper_84 = helper_66 * helper_82;
274  double helper_85 = -helper_32 - helper_48 + helper_53;
275  double helper_86 = 1.0 / helper_85;
276  double helper_87 = helper_86 * pow(pow(helper_85, 2), -0.333333333333333);
277  double helper_88 = 0.707106781186548 * helper_6;
278  double helper_89 = 0.707106781186548 * helper_27;
279  double helper_90 = helper_88 - helper_89;
280  double helper_91 =
281  0.666666666666667 * helper_10 * helper_40 + 0.666666666666667 * helper_3 * helper_90 -
282  0.666666666666667 * helper_30 * helper_45 - 0.666666666666667 * helper_33 * helper_62;
283  double helper_92 = -3.0 * helper_11 + 1.0 * helper_13 + 1.0 * helper_15 + 1.0 * helper_17;
284  double helper_93 = -helper_11 + helper_17;
285  double helper_94 = -helper_1 + helper_2;
286  double helper_95 = -helper_21 + helper_22 - helper_23;
287  double helper_96 = -helper_34 - helper_36 + helper_38 - helper_39;
288  double helper_97 = -helper_42 + helper_43 - helper_44;
289  double helper_98 = -helper_12 - helper_14 + helper_16 - helper_18;
290  double helper_99 =
291  -0.666666666666667 * helper_60 * helper_94 + 0.666666666666667 * helper_62 * helper_93 +
292  0.666666666666667 * helper_95 * helper_96 - 0.666666666666667 * helper_97 * helper_98;
293  double helper_100 = helper_3 * helper_90;
294  double helper_101 = helper_33 * helper_62;
295  double helper_102 = helper_100 - helper_101 + helper_52;
296  double helper_103 = -helper_60 * helper_94 + helper_62 * helper_93 + helper_95 * helper_96 -
297  helper_97 * helper_98;
298  double helper_104 = 0.444444444444444 * helper_102 * helper_103 * helper_82 * helper_86 +
299  helper_57 * helper_91 - helper_92 * helper_99;
300  double helper_105 =
301  1.85037170770859e-17 * helper_1 * helper_78 + 1.85037170770859e-17 * helper_11 * helper_73 +
302  1.85037170770859e-17 * helper_13 * helper_76 +
303  1.85037170770859e-17 * helper_15 * helper_75 +
304  1.85037170770859e-17 * helper_17 * helper_74 + 1.85037170770859e-17 * helper_2 * helper_79 +
305  1.85037170770859e-17 * helper_27 * helper_70 +
306  1.85037170770859e-17 * helper_35 * helper_81 +
307  1.85037170770859e-17 * helper_37 * helper_80 + 1.85037170770859e-17 * helper_4 * helper_68 +
308  1.85037170770859e-17 * helper_6 * helper_71 + 1.85037170770859e-17 * helper_69 * helper_8;
309  double helper_106 = helper_64 * helper_82 * helper_86;
310  double helper_107 =
311  -0.666666666666667 * helper_10 * helper_19 + 0.666666666666667 * helper_24 * helper_30 +
312  0.666666666666667 * helper_33 * helper_60 - 0.666666666666667 * helper_49 * helper_90;
313  double helper_108 = -3.0 * helper_1 + 1.0 * helper_2 + 1.0 * helper_35 + 1.0 * helper_37;
314  double helper_109 = -helper_20 + helper_31 + helper_33 * helper_60 - helper_49 * helper_90;
315  double helper_110 = 0.444444444444444 * helper_109 * helper_82 * helper_86;
316  double helper_111 = helper_103 * helper_110 + helper_107 * helper_57 - helper_108 * helper_99;
317  double helper_112 = -helper_4 + helper_8;
318  double helper_113 = -helper_88 + helper_89;
319  double helper_114 = -helper_5 + helper_7 - helper_9;
320  double helper_115 = -helper_25 - helper_26 + helper_28 - helper_29;
321  double helper_116 = helper_82 * helper_86 *
322  (helper_112 * helper_62 + helper_113 * helper_94 + helper_114 * helper_96 -
323  helper_115 * helper_97);
324  double helper_117 = -helper_100 + helper_101 - helper_50 + helper_51;
325  double helper_118 = -helper_102 * helper_110 + helper_107 * helper_92 + helper_108 * helper_91;
326  double helper_119 = helper_82 * helper_86 *
327  (helper_112 * (-helper_58 + helper_59) - helper_113 * helper_93 -
328  helper_114 * helper_98 + helper_115 * helper_95);
329  result_0(0, 0) =
330  helper_56 * (helper_57 * helper_64 * helper_65 - pow(helper_64, 2) * helper_83 +
331  0.666666666666667 * helper_64 * helper_84 *
332  (-helper_41 + helper_46 - helper_61 + helper_63) +
333  3.0);
334  result_0(0, 1) = helper_87 * (helper_104 - helper_105 * helper_35 + helper_106 * helper_91);
335  result_0(0, 2) = helper_87 * (helper_106 * helper_107 + helper_111);
336  result_0(1, 0) = helper_87 * (helper_104 + helper_116 * helper_99);
337  result_0(1, 1) =
338  helper_56 * (-pow(helper_117, 2) * helper_83 + helper_117 * helper_65 * helper_92 +
339  helper_117 * helper_84 * helper_91 + 3.0);
340  result_0(1, 2) = helper_87 * (-helper_105 * helper_6 - helper_107 * helper_116 + helper_118);
341  result_0(2, 0) = helper_87 * (-helper_105 * helper_13 + helper_111 + helper_119 * helper_99);
342  result_0(2, 1) = helper_87 * (helper_118 - helper_119 * helper_91);
343  result_0(2, 2) = helper_56 * (-helper_108 * helper_109 * helper_65 -
344  1.11111111111111 * pow(helper_109, 2) * helper_84 + 3.0);
345 }
346 
347 /*bool tet_is_energy_unstable(const std::array<double, 12>& T, double res)
348 {
349  static const std::vector<std::array<int, 4>> combs = {
350  {{0, 1, 3, 2}}, {{0, 2, 1, 3}}, {{0, 2, 3, 1}}, {{0, 3, 1, 2}}, {{0, 3, 2, 1}},
351  {{1, 0, 2, 3}}, {{1, 0, 3, 2}}, {{1, 2, 0, 3}}, {{1, 2, 3, 0}}, {{1, 3, 0, 2}},
352  {{1, 3, 2, 0}}, {{2, 0, 1, 3}}, {{2, 0, 3, 1}}, {{2, 1, 0, 3}}, {{2, 1, 3, 0}},
353  {{2, 3, 0, 1}}, {{2, 3, 1, 0}}, {{3, 0, 1, 2}}, {{3, 0, 2, 1}}, {{3, 1, 0, 2}},
354  {{3, 1, 2, 0}}, {{3, 2, 0, 1}}, {{3, 2, 1, 0}}};
355 
356  double res0 = -1;
357  if (std::isinf(res)) return true;
358 
359  for (int i = 0; i < combs.size(); i++) {
360  std::array<double, 12> tmp_T;
361  for (int j = 0; j < 4; j++) {
362  for (int k = 0; k < 3; k++) tmp_T[j * 3 + k] = T[combs[i][j] * 3 + k];
363  }
364  double res1 = Tet_AMIPS_energy_aux(tmp_T);
365 
366  if (std::isinf(res1)) continue;
367 
368  if (res0 == -1) res0 = res1;
369 
370  if (std::abs(res1 - res0) / res0 > 0.01) return true;
371  }
372  return false;
373 }*/
374 
375 double Tet_AMIPS_energy(const std::array<double, 12>& T)
376 {
377  double res = Tet_AMIPS_energy_aux(T);
378 
379  // Maybe use tet_is_energy_unstable
380  if (res > 1e8) {
381  std::array<Rational, 12> r_T;
382  for (int j = 0; j < 12; j++) {
383  r_T[j] = T[j];
384  }
385 
386  const Rational twothird = Rational(2) / Rational(3);
387 
388  Rational tmp =
389  ((-r_T[1 + 2] + r_T[1 + 5]) * r_T[1 + 1] + r_T[1 + 2] * r_T[1 + 7] +
390  (r_T[1 + -1] - r_T[1 + 5]) * r_T[1 + 4] - r_T[1 + -1] * r_T[1 + 7]) *
391  r_T[1 + 9] +
392  ((r_T[1 + 2] - r_T[1 + 5]) * r_T[1 + 0] - r_T[1 + 2] * r_T[1 + 6] +
393  (-r_T[1 + -1] + r_T[1 + 5]) * r_T[1 + 3] + r_T[1 + -1] * r_T[1 + 6]) *
394  r_T[1 + 10] +
395  (-r_T[1 + 2] * r_T[1 + 7] + (-r_T[1 + 8] + r_T[1 + 5]) * r_T[1 + 4] +
396  r_T[1 + 8] * r_T[1 + 7]) *
397  r_T[1 + 0] +
398  (r_T[1 + 2] * r_T[1 + 6] + (r_T[1 + 8] - r_T[1 + 5]) * r_T[1 + 3] -
399  r_T[1 + 8] * r_T[1 + 6]) *
400  r_T[1 + 1] +
401  (r_T[1 + 3] * r_T[1 + 7] - r_T[1 + 4] * r_T[1 + 6]) * (r_T[1 + -1] - r_T[1 + 8]);
402 
403  if (tmp == 0) return std::numeric_limits<double>::infinity();
404 
405  auto res_r =
406  Rational(27) / 16 * pow(tmp, -2) *
407  pow(r_T[1 + 9] * r_T[1 + 9] +
408  (-twothird * r_T[1 + 0] - twothird * r_T[1 + 3] - twothird * r_T[1 + 6]) *
409  r_T[1 + 9] +
410  r_T[1 + 10] * r_T[1 + 10] +
411  (-twothird * r_T[1 + 1] - twothird * r_T[1 + 4] - twothird * r_T[1 + 7]) *
412  r_T[1 + 10] +
413  r_T[1 + 0] * r_T[1 + 0] +
414  (-twothird * r_T[1 + 3] - twothird * r_T[1 + 6]) * r_T[1 + 0] +
415  r_T[1 + 1] * r_T[1 + 1] +
416  (-twothird * r_T[1 + 4] - twothird * r_T[1 + 7]) * r_T[1 + 1] +
417  r_T[1 + 2] * r_T[1 + 2] +
418  (-twothird * r_T[1 + -1] - twothird * r_T[1 + 8] - twothird * r_T[1 + 5]) *
419  r_T[1 + 2] +
420  r_T[1 + 3] * r_T[1 + 3] - twothird * r_T[1 + 3] * r_T[1 + 6] +
421  r_T[1 + 4] * r_T[1 + 4] - twothird * r_T[1 + 4] * r_T[1 + 7] +
422  r_T[1 + 5] * r_T[1 + 5] +
423  (-twothird * r_T[1 + -1] - twothird * r_T[1 + 8]) * r_T[1 + 5] -
424  twothird * r_T[1 + -1] * r_T[1 + 8] + r_T[1 + -1] * r_T[1 + -1] +
425  r_T[1 + 8] * r_T[1 + 8] + r_T[1 + 6] * r_T[1 + 6] + r_T[1 + 7] * r_T[1 + 7],
426  3);
427 
428  return std::cbrt(res_r.to_double());
429  } else {
430  return res;
431  }
432 
433 } // namespace
434 
435 double Tri_AMIPS_energy_aux(const std::array<double, 6>& T)
436 {
437  double helper_0[6];
438  helper_0[0] = T[0];
439  helper_0[1] = T[1];
440  helper_0[2] = T[2];
441  helper_0[3] = T[3];
442  helper_0[4] = T[4];
443  helper_0[5] = T[5];
444  double helper_1 = helper_0[0];
445  double helper_2 = helper_0[2];
446  double helper_3 = helper_0[1];
447  double helper_4 = helper_0[3];
448  double helper_5 = helper_0[5];
449  double helper_6 = helper_0[4];
450  double helper_7 = 0.666666666666667 * helper_6;
451  double helper_8 = 0.666666666666667 * helper_5;
452  double denom =
453  ((helper_1 - helper_2) * (0.577350269189626 * helper_3 + 0.577350269189626 * helper_4 -
454  1.15470053837925 * helper_5) -
455  (helper_3 - helper_4) * (0.577350269189626 * helper_1 + 0.577350269189626 * helper_2 -
456  1.15470053837925 * helper_6));
457 
458  denom = std::abs(denom) * wmtk::utils::wmtk_orient2d(T[0], T[1], T[2], T[3], T[4], T[5]);
459 
460  // if (std::abs(denom) < 1e-12) return std::numeric_limits<double>::infinity();
461 
462  return -(helper_1 * (-1.33333333333333 * helper_1 + 0.666666666666667 * helper_2 + helper_7) +
463  helper_2 * (0.666666666666667 * helper_1 - 1.33333333333333 * helper_2 + helper_7) +
464  helper_3 * (-1.33333333333333 * helper_3 + 0.666666666666667 * helper_4 + helper_8) +
465  helper_4 * (0.666666666666667 * helper_3 - 1.33333333333333 * helper_4 + helper_8) +
466  helper_5 * (0.666666666666667 * helper_3 + 0.666666666666667 * helper_4 -
467  1.33333333333333 * helper_5) +
468  helper_6 * (0.666666666666667 * helper_1 + 0.666666666666667 * helper_2 -
469  1.33333333333333 * helper_6)) /
470  denom;
471 }
472 
473 double Tri_AMIPS_energy(const std::array<double, 6>& T)
474 {
475  static const double sqrt3 = sqrt(3);
476  static const Rational two_third = Rational(2) / Rational(3);
477 
478  double res = Tri_AMIPS_energy_aux(T);
479 
480  // Maybe use tet_is_energy_unstable
481  if (res > 1e8 || res <= 0) {
482  std::array<Rational, 6> r_T;
483  for (int j = 0; j < 6; j++) {
484  r_T[j] = T[j];
485  }
486 
487  const auto tmp = r_T[0] * r_T[3] - r_T[0] * r_T[5] - r_T[1] * r_T[2] + r_T[1] * r_T[4] +
488  r_T[2] * r_T[5] - r_T[3] * r_T[4];
489 
490  if (tmp <= 0) return std::numeric_limits<double>::infinity();
491 
492  const auto res_r = two_third * sqrt3 *
493  (r_T[0] * r_T[0] - r_T[0] * r_T[2] - r_T[0] * r_T[4] + r_T[1] * r_T[1] -
494  r_T[1] * r_T[3] - r_T[1] * r_T[5] + r_T[2] * r_T[2] - r_T[2] * r_T[4] +
495  r_T[3] * r_T[3] - r_T[3] * r_T[5] + r_T[4] * r_T[4] + r_T[5] * r_T[5]) /
496  tmp;
497  return res_r.to_double();
498  } else {
499  return res;
500  }
501 }
502 
503 void Tri_AMIPS_jacobian(const std::array<double, 6>& T, Eigen::Vector2d& result_0)
504 {
505  double helper_0[6];
506  helper_0[0] = T[0];
507  helper_0[1] = T[1];
508  helper_0[2] = T[2];
509  helper_0[3] = T[3];
510  helper_0[4] = T[4];
511  helper_0[5] = T[5];
512  double helper_1 = helper_0[0];
513  double helper_2 = helper_0[2];
514  double helper_3 = helper_0[1];
515  double helper_4 = helper_0[3];
516  double helper_5 = helper_0[5];
517  double helper_6 = helper_0[4];
518  double helper_7 =
519  1.0 /
520  ((helper_1 - helper_2) * (0.577350269189626 * helper_3 + 0.577350269189626 * helper_4 -
521  1.15470053837925 * helper_5) -
522  (helper_3 - helper_4) * (0.577350269189626 * helper_1 + 0.577350269189626 * helper_2 -
523  1.15470053837925 * helper_6));
524  double helper_8 = -1.33333333333333 * helper_6;
525  double helper_9 = 0.666666666666667 * helper_6;
526  double helper_10 = 0.666666666666667 * helper_5;
527  double helper_11 = 1.33333333333333 * helper_5;
528  double helper_12 =
529  1.15470053837925 * helper_7 *
530  (helper_1 * (-1.33333333333333 * helper_1 + 0.666666666666667 * helper_2 + helper_9) +
531  helper_2 * (0.666666666666667 * helper_1 - 1.33333333333333 * helper_2 + helper_9) +
532  helper_3 * (helper_10 - 1.33333333333333 * helper_3 + 0.666666666666667 * helper_4) +
533  helper_4 * (helper_10 + 0.666666666666667 * helper_3 - 1.33333333333333 * helper_4) +
534  helper_5 * (-helper_11 + 0.666666666666667 * helper_3 + 0.666666666666667 * helper_4) +
535  helper_6 * (0.666666666666667 * helper_1 + 0.666666666666667 * helper_2 + helper_8));
536  result_0(0) = helper_7 * (2.66666666666667 * helper_1 + helper_12 * (helper_4 - helper_5) -
537  1.33333333333333 * helper_2 + helper_8);
538  result_0(1) = -helper_7 * (helper_11 + helper_12 * (helper_2 - helper_6) -
539  2.66666666666667 * helper_3 + 1.33333333333333 * helper_4);
540 }
541 
542 void Tri_AMIPS_hessian(const std::array<double, 6>& T, Eigen::Matrix2d& result_0)
543 {
544  double helper_0[6];
545  helper_0[0] = T[0];
546  helper_0[1] = T[1];
547  helper_0[2] = T[2];
548  helper_0[3] = T[3];
549  helper_0[4] = T[4];
550  helper_0[5] = T[5];
551  double helper_1 = helper_0[0];
552  double helper_2 = helper_0[2];
553  double helper_3 = helper_0[1];
554  double helper_4 = helper_0[3];
555  double helper_5 = helper_0[5];
556  double helper_6 = helper_0[4];
557  double helper_7 =
558  (helper_1 - helper_2) * (0.577350269189626 * helper_3 + 0.577350269189626 * helper_4 -
559  1.15470053837925 * helper_5) -
560  (helper_3 - helper_4) * (0.577350269189626 * helper_1 + 0.577350269189626 * helper_2 -
561  1.15470053837925 * helper_6);
562  double helper_8 = 1.0 / helper_7;
563  double helper_9 = helper_4 - helper_5;
564  double helper_10 = 1.33333333333333 * helper_6;
565  double helper_11 = -2.66666666666667 * helper_1 + helper_10 + 1.33333333333333 * helper_2;
566  double helper_12 = 2.3094010767585 * helper_8;
567  double helper_13 = pow(helper_7, -2);
568  double helper_14 = 0.666666666666667 * helper_6;
569  double helper_15 = 0.666666666666667 * helper_5;
570  double helper_16 = 1.33333333333333 * helper_5;
571  double helper_17 =
572  helper_1 * (-1.33333333333333 * helper_1 + helper_14 + 0.666666666666667 * helper_2) +
573  helper_2 * (0.666666666666667 * helper_1 + helper_14 - 1.33333333333333 * helper_2) +
574  helper_3 * (helper_15 - 1.33333333333333 * helper_3 + 0.666666666666667 * helper_4) +
575  helper_4 * (helper_15 + 0.666666666666667 * helper_3 - 1.33333333333333 * helper_4) +
576  helper_5 * (-helper_16 + 0.666666666666667 * helper_3 + 0.666666666666667 * helper_4) +
577  helper_6 * (0.666666666666667 * helper_1 - helper_10 + 0.666666666666667 * helper_2);
578  double helper_18 = 2.66666666666667 * helper_13 * helper_17;
579  double helper_19 = helper_2 - helper_6;
580  double helper_20 = helper_16 - 2.66666666666667 * helper_3 + 1.33333333333333 * helper_4;
581  double helper_21 = helper_13 * (-1.15470053837925 * helper_11 * helper_19 +
582  2.66666666666667 * helper_17 * helper_19 * helper_8 * helper_9 +
583  1.15470053837925 * helper_20 * helper_9);
584  result_0(0) = helper_8 * (helper_11 * helper_12 * helper_9 - helper_18 * pow(helper_9, 2) +
585  2.66666666666667);
586  result_0(1) = helper_21;
587  result_0(2) = helper_21;
588  result_0(3) = helper_8 * (-helper_12 * helper_19 * helper_20 - helper_18 * pow(helper_19, 2) +
589  2.66666666666667);
590 }
591 
592 template <int64_t NV, int64_t DIM>
593 std::array<double, NV * DIM> unbox(
594  const std::vector<std::decay_t<typename attribute::ConstMapResult<double>>>& data,
595  const int64_t index)
596 {
597  std::array<double, NV * DIM> res;
598  assert(data.size() == NV);
599 
600  const size_t start = index < 0 ? 0 : index;
601 
602  for (size_t i = 0; i < NV; ++i) {
603  const size_t ii = (i + start) % NV;
604  assert(data[ii].size() == DIM);
605 
606  for (size_t j = 0; j < DIM; ++j) {
607  res[DIM * i + j] = data[ii][j];
608  }
609  }
610 
611  return res;
612 }
613 
614 template <int64_t NV, int64_t DIM>
615 std::array<Rational, NV * DIM> unbox(
616  const std::vector<std::decay_t<typename attribute::ConstMapResult<Rational>>>& data,
617  const int64_t index)
618 {
619  std::array<Rational, NV * DIM> res;
620  assert(data.size() == NV);
621 
622  const size_t start = index < 0 ? 0 : index;
623 
624  for (size_t i = 0; i < NV; ++i) {
625  const size_t ii = (i + start) % NV;
626  assert(data[ii].size() == DIM);
627 
628  for (size_t j = 0; j < DIM; ++j) {
629  res[DIM * i + j] = data[ii][j];
630  }
631  }
632 
633  return res;
634 }
635 
636 
637 template <int64_t NV, int64_t DIM>
638 std::array<double, NV * DIM> AMIPS::get_raw_coordinates(
639  const simplex::Simplex& domain_simplex,
640  const std::optional<simplex::Simplex>& variable_simplex) const
641 {
642  if (embedded_dimension() != DIM) throw std::runtime_error("AMIPS wrong dimension");
643 
644  if (attribute_handle().holds<double>()) {
645  const attribute::Accessor<double> accessor =
646  mesh().create_const_accessor(attribute_handle().as<double>());
647 
648  auto [attrs, index] = utils::get_simplex_attributes(
649  mesh(),
650  accessor,
651  m_primitive_type,
652  domain_simplex,
653  variable_simplex.has_value() ? variable_simplex->tuple() : std::optional<Tuple>());
654 
655  return unbox<NV, DIM>(attrs, index);
656  } else if (attribute_handle().holds<Rational>()) {
657  const attribute::Accessor<Rational> accessor =
658  mesh().create_const_accessor(attribute_handle().as<Rational>());
659 
660  auto [attrs, index] = utils::get_simplex_attributes(
661  mesh(),
662  accessor,
663  m_primitive_type,
664  domain_simplex,
665  variable_simplex.has_value() ? variable_simplex->tuple() : std::optional<Tuple>());
666 
667  const auto coord_rational = unbox<NV, DIM>(attrs, index);
668  std::array<double, coord_rational.size()> coord;
669 
670  for (int i = 0; i < coord_rational.size(); ++i) {
671  coord[i] = coord_rational[i].to_double();
672  }
673 
674  return coord;
675  }
676 
677  assert(false);
678  throw std::runtime_error("AMIP coordinate handle with invalid type");
679 }
680 
681 
682 AMIPS::AMIPS(const Mesh& mesh, const attribute::MeshAttributeHandle& attribute_handle)
683  : PerSimplexFunction(mesh, PrimitiveType::Vertex, attribute_handle)
684 
685 {}
686 
687 double AMIPS::get_value(const simplex::Simplex& domain_simplex) const
688 {
689  static constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
690 
691  double res = 0;
692  if (domain_simplex.primitive_type() == PrimitiveType::Tetrahedron)
693  res = Tet_AMIPS_energy(get_raw_coordinates<4, 3>(domain_simplex));
694  else if (domain_simplex.primitive_type() == PrimitiveType::Triangle)
695  res = Tri_AMIPS_energy(get_raw_coordinates<3, 2>(domain_simplex));
696  else
697  throw std::runtime_error("AMIPS wrong simplex type");
698 
699  assert(res >= 0);
700  return res;
701 }
702 
703 Eigen::VectorXd AMIPS::get_gradient(
704  const simplex::Simplex& domain_simplex,
705  const simplex::Simplex& variable_simplex) const
706 {
707  if (domain_simplex.primitive_type() == PrimitiveType::Tetrahedron) {
708  Eigen::Vector3d res;
709  Tet_AMIPS_jacobian(get_raw_coordinates<4, 3>(domain_simplex, variable_simplex), res);
710  return res;
711  } else if (domain_simplex.primitive_type() == PrimitiveType::Triangle) {
712  Eigen::Vector2d res;
713  Tri_AMIPS_jacobian(get_raw_coordinates<3, 2>(domain_simplex, variable_simplex), res);
714  return res;
715  } else
716  throw std::runtime_error("AMIPS wrong simplex type");
717 }
718 
719 Eigen::MatrixXd AMIPS::get_hessian(
720  const simplex::Simplex& domain_simplex,
721  const simplex::Simplex& variable_simplex) const
722 {
723  if (domain_simplex.primitive_type() == PrimitiveType::Tetrahedron) {
724  Eigen::Matrix3d res;
725  Tet_AMIPS_hessian(get_raw_coordinates<4, 3>(domain_simplex, variable_simplex), res);
726  return res;
727  } else if (domain_simplex.primitive_type() == PrimitiveType::Triangle) {
728  Eigen::Matrix2d res;
729  Tri_AMIPS_hessian(get_raw_coordinates<3, 2>(domain_simplex, variable_simplex), res);
730  return res;
731  } else
732  throw std::runtime_error("AMIPS wrong simplex type");
733 }
734 
735 } // namespace wmtk::function
wmtk::Mesh::create_const_accessor
const attribute::Accessor< T, Mesh, D > create_const_accessor(const attribute::MeshAttributeHandle &handle) const
wmtk::Rational::to_double
double to_double() const
Definition: Rational.cpp:317
wmtk::attribute::Accessor
A CachingAccessor that uses tuples for accessing attributes instead of indices.
Definition: Accessor.hpp:25
wmtk::function::AMIPS::get_gradient
Eigen::VectorXd get_gradient(const simplex::Simplex &domain_simplex, const simplex::Simplex &variable_simplex) const override
Definition: AMIPS.cpp:703
wmtk::function::AMIPS::get_raw_coordinates
std::array< double, NV *DIM > get_raw_coordinates(const simplex::Simplex &domain_simplex, const std::optional< simplex::Simplex > &variable_simplex={}) const
Definition: AMIPS.cpp:638
wmtk::function::AMIPS::AMIPS
AMIPS(const Mesh &mesh, const attribute::MeshAttributeHandle &attribute_handle)
Construct a new AMIPS function.
Definition: AMIPS.cpp:682
wmtk::function::AMIPS::get_hessian
Eigen::MatrixXd get_hessian(const simplex::Simplex &domain_simplex, const simplex::Simplex &variable_simplex) const override
Definition: AMIPS.cpp:719
wmtk::function::AMIPS::get_value
double get_value(const simplex::Simplex &domain_simplex) const override
This function is defined over a simplex (normally a triangle or tetrahedron).
Definition: AMIPS.cpp:687
wmtk::function::PerSimplexFunction::attribute_handle
const attribute::MeshAttributeHandle & attribute_handle() const
Definition: PerSimplexFunction.hpp:42
wmtk::simplex::Simplex::tuple
const Tuple & tuple() const
Definition: Simplex.hpp:53
wmtk::simplex::Simplex::primitive_type
PrimitiveType primitive_type() const
Definition: Simplex.hpp:51
wmtk::attribute::ConstMapResult
typename internal::VectorResult< T, R >::ConstMapType ConstMapResult
Definition: MapTypes.hpp:12
wmtk::function::utils::get_simplex_attributes
std::tuple< std::vector< std::decay_t< typename attribute::ConstMapResult< T > > >, int64_t > get_simplex_attributes(const Mesh &mesh, const wmtk::attribute::Accessor< T > &accessor, const PrimitiveType primitive_type, const simplex::Simplex &simplex_in, const std::optional< wmtk::Tuple > &vertex_marker)
get attributes from a simplex
Definition: SimplexGetter.cpp:9
wmtk::function::Tri_AMIPS_energy_aux
double Tri_AMIPS_energy_aux(const std::array< double, 6 > &T)
Definition: AMIPS.cpp:435
wmtk::function::Tet_AMIPS_jacobian
void Tet_AMIPS_jacobian(const std::array< double, 12 > &T, Eigen::Vector3d &result_0)
Definition: AMIPS.cpp:75
wmtk::function::Tet_AMIPS_energy_aux
double Tet_AMIPS_energy_aux(const std::array< double, 12 > &T)
Definition: AMIPS.cpp:12
wmtk::function::Tri_AMIPS_hessian
void Tri_AMIPS_hessian(const std::array< double, 6 > &T, Eigen::Matrix2d &result_0)
Definition: AMIPS.cpp:542
wmtk::function::Tri_AMIPS_jacobian
void Tri_AMIPS_jacobian(const std::array< double, 6 > &T, Eigen::Vector2d &result_0)
Definition: AMIPS.cpp:503
wmtk::function::Tet_AMIPS_energy
double Tet_AMIPS_energy(const std::array< double, 12 > &T)
Definition: AMIPS.cpp:375
wmtk::function::Tet_AMIPS_hessian
void Tet_AMIPS_hessian(const std::array< double, 12 > &T, Eigen::Matrix3d &result_0)
Definition: AMIPS.cpp:172
wmtk::function::Tri_AMIPS_energy
double Tri_AMIPS_energy(const std::array< double, 6 > &T)
Definition: AMIPS.cpp:473
wmtk::function::unbox
std::array< double, NV *DIM > unbox(const std::vector< std::decay_t< typename attribute::ConstMapResult< double >>> &data, const int64_t index)
Definition: AMIPS.cpp:593
wmtk::utils::wmtk_orient2d
int wmtk_orient2d(double p0x, double p0y, double p1x, double p1y, double p2x, double p2y)
Definition: orient.cpp:178
wmtk::abs
Rational abs(const Rational &r0)
Definition: Rational.cpp:328
wmtk::Vector3d
Vector< double, 3 > Vector3d
Definition: Types.hpp:39
wmtk::pow
Rational pow(const Rational &x, int p)
Definition: Rational.cpp:166