PolyFEM
Loading...
Searching...
No Matches
AMIPSForm.cpp
Go to the documentation of this file.
1#include "AMIPSForm.hpp"
2
3#include <igl/predicates/predicates.h>
4
5namespace polyfem
6{
7 namespace solver
8 {
9 AMIPSForm::AMIPSForm(
10 const Eigen::MatrixXd X_rest,
11 const Eigen::MatrixXd X)
12 : X_rest_(X_rest), X_(X)
13 {
14 }
15
16 bool AMIPSForm::is_step_valid(const Eigen::VectorXd &, const Eigen::VectorXd &x1) const
17 {
18 igl::predicates::exactinit();
19
20 // Use igl for checking orientation
21 igl::predicates::Orientation res = igl::predicates::orient2d(
22 Eigen::Vector2d(x1.head<2>()),
23 Eigen::Vector2d(X_rest_.row(1).head<2>()),
24 Eigen::Vector2d(X_rest_.row(2).head<2>()));
25
26 // The element is inverted if it not positive (i.e. it is negative or it is degenerate)
27 return res == igl::predicates::Orientation::POSITIVE;
28 }
29
30 double AMIPSForm::energy(
31 const Eigen::Vector2d &x0_rest,
32 const Eigen::Vector2d &x1_rest,
33 const Eigen::Vector2d &x2_rest,
34 const Eigen::Vector2d &x0,
35 const Eigen::Vector2d &x1,
36 const Eigen::Vector2d &x2)
37 {
38 return autogen::AMIPS2D_energy(
39 x0_rest.x(), x0_rest.y(),
40 x1_rest.x(), x1_rest.y(),
41 x2_rest.x(), x2_rest.y(),
42 x0.x(), x0.y(),
43 x1.x(), x1.y(),
44 x2.x(), x2.y());
45 }
46
47 double AMIPSForm::value_unweighted(const Eigen::VectorXd &x) const
48 {
49 if (x.size() == 2)
50 {
51 return energy(
52 x.head<2>(),
53 X_rest_.row(1).head<2>(),
54 X_rest_.row(2).head<2>(),
55 X_.row(0).head<2>(),
56 X_.row(1).head<2>(),
57 X_.row(2).head<2>());
58 }
59 else
60 {
61 assert(x.size() == 3);
62 throw std::runtime_error("AMIPSForm::first_derivative_unweighted not implemented for 3D");
63 }
64 }
65
66 void AMIPSForm::first_derivative_unweighted(
67 const Eigen::VectorXd &x, Eigen::VectorXd &grad) const
68 {
69 grad.resize(x.size());
70
71 if (x.size() == 2)
72 {
73 autogen::AMIPS2D_gradient(
74 x[0], x[1],
75 X_rest_(1, 0), X_rest_(1, 1),
76 X_rest_(2, 0), X_rest_(2, 1),
77 X_(0, 0), X_(0, 1),
78 X_(1, 0), X_(1, 1),
79 X_(2, 0), X_(2, 1),
80 grad.data());
81 }
82 else
83 {
84 assert(x.size() == 3);
85 throw std::runtime_error("AMIPSForm::first_derivative_unweighted not implemented for 3D");
86 }
87 }
88
89 void AMIPSForm::second_derivative_unweighted(
90 const Eigen::VectorXd &x, StiffnessMatrix &hessian) const
91 {
92 Eigen::MatrixXd H(x.size(), x.size());
93 if (x.size() == 2)
94 {
95 // NOTE: it doesnt matter if H is column major or row major because the hessian is symmetric
96 autogen::AMIPS2D_hessian(
97 x[0], x[1],
98 X_rest_(1, 0), X_rest_(1, 1),
99 X_rest_(2, 0), X_rest_(2, 1),
100 X_(0, 0), X_(0, 1),
101 X_(1, 0), X_(1, 1),
102 X_(2, 0), X_(2, 1),
103 H.data());
104 }
105 else
106 {
107 assert(x.size() == 3);
108 throw std::runtime_error("AMIPSForm::second_derivative_unweighted not implemented for 3D");
109 }
110 hessian = H.sparseView();
111 }
112 } // namespace solver
113
114 namespace autogen
115 {
116 double AMIPS2D_energy(double x0_rest, double y0_rest, double x1_rest, double y1_rest, double x2_rest, double y2_rest, double x0, double y0, double x1, double y1, double x2, double y2)
117 {
118 const auto t0 = x0 * y1 - x0 * y2 - x1 * y0 + x1 * y2 + x2 * y0 - x2 * y1;
119 const auto t1 = 1.0 / (x0_rest * y1_rest - x0_rest * y2_rest - x1_rest * y0_rest + x1_rest * y2_rest + x2_rest * y0_rest - x2_rest * y1_rest);
120 const auto t2 = x0 - x1;
121 const auto t3 = x0_rest - x2_rest;
122 const auto t4 = x0 - x2;
123 const auto t5 = x0_rest - x1_rest;
124 const auto t6 = y0_rest - y2_rest;
125 const auto t7 = y0_rest - y1_rest;
126 const auto t8 = y0 - y2;
127 const auto t9 = y0 - y1;
128 return ((t0 * t1 > 0) ? (
129 t1 * (std::pow(t2 * t3 - t4 * t5, 2) + std::pow(t2 * t6 - t4 * t7, 2) + std::pow(-t3 * t9 + t5 * t8, 2) + std::pow(t6 * t9 - t7 * t8, 2)) / t0)
130 : (
131 std::numeric_limits<double>::infinity()));
132 }
133
134 void AMIPS2D_gradient(double x0_rest, double y0_rest, double x1_rest, double y1_rest, double x2_rest, double y2_rest, double x0, double y0, double x1, double y1, double x2, double y2, double g[2])
135 {
136 const auto t0 = x0 - x1;
137 const auto t1 = -y2_rest;
138 const auto t2 = t1 + y0_rest;
139 const auto t3 = -x2;
140 const auto t4 = t3 + x0;
141 const auto t5 = y0_rest - y1_rest;
142 const auto t6 = t4 * t5;
143 const auto t7 = t0 * t2 - t6;
144 const auto t8 = std::pow(t7, 2);
145 const auto t9 = 1.0 / (x0_rest * y1_rest - x0_rest * y2_rest - x1_rest * y0_rest + x1_rest * y2_rest + x2_rest * y0_rest - x2_rest * y1_rest);
146 const auto t10 = t1 + y1_rest;
147 const auto t11 = t10 * t9;
148 const auto t12 = 2 * t11;
149 const auto t13 = y0 - y1;
150 const auto t14 = -y2;
151 const auto t15 = t14 + y0;
152 const auto t16 = t15 * t5;
153 const auto t17 = t13 * t2 - t16;
154 const auto t18 = std::pow(t17, 2);
155 const auto t19 = -x2_rest;
156 const auto t20 = t19 + x0_rest;
157 const auto t21 = t0 * t20;
158 const auto t22 = x0_rest - x1_rest;
159 const auto t23 = t22 * t4;
160 const auto t24 = t21 - t23;
161 const auto t25 = t3 + x1;
162 const auto t26 = t13 * t20;
163 const auto t27 = t15 * t22 - t26;
164 const auto t28 = t14 + y1;
165 const auto t29 = std::pow(t27, 2);
166 const auto t30 = t18 + t29 + t8;
167 const auto t31 = x0 * y1 - x0 * y2 - x1 * y0 + x1 * y2 + x2 * y0 - x2 * y1;
168 const auto t32 = t9 / t31;
169 const auto t33 = t31 * t9 > 0;
170 const auto t34 = std::pow(t24, 2);
171 const auto t35 = t19 + x1_rest;
172 const auto t36 = t35 * t9;
173 const auto t37 = 2 * t36;
174 if (t33)
175 {
176 g[0] = t32 * (t10 * t9 * (std::pow(t24, 2) + t30) - t12 * t18 - t12 * t8 - 2 * t24 * (t11 * t21 - t11 * t23 + t25) - 2 * t27 * (t10 * t15 * t22 * t9 - t11 * t26 - t28));
177 }
178 else
179 {
180 g[0] = 0;
181 }
182 if (t33)
183 {
184 g[1] = t32 * (2 * t17 * (t13 * t2 * t35 * t9 - t16 * t36 - t28) + t29 * t37 + t34 * t37 - t36 * (t30 + t34) + 2 * t7 * (t0 * t2 * t35 * t9 - t25 - t36 * t6));
185 }
186 else
187 {
188 g[1] = 0;
189 }
190 }
191
192 void AMIPS2D_hessian(double x0_rest, double y0_rest, double x1_rest, double y1_rest, double x2_rest, double y2_rest, double x0, double y0, double x1, double y1, double x2, double y2, double H[4])
193 {
194 const auto t0 = x0 - x1;
195 const auto t1 = -y2_rest;
196 const auto t2 = t1 + y0_rest;
197 const auto t3 = -x2;
198 const auto t4 = t3 + x0;
199 const auto t5 = y0_rest - y1_rest;
200 const auto t6 = t4 * t5;
201 const auto t7 = t0 * t2 - t6;
202 const auto t8 = std::pow(t7, 2);
203 const auto t9 = t1 + y1_rest;
204 const auto t10 = x0_rest * y1_rest - x0_rest * y2_rest - x1_rest * y0_rest + x1_rest * y2_rest + x2_rest * y0_rest - x2_rest * y1_rest;
205 const auto t11 = std::pow(t10, -2);
206 const auto t12 = 3 * t11;
207 const auto t13 = t12 * std::pow(t9, 2);
208 const auto t14 = y0 - y1;
209 const auto t15 = -y2;
210 const auto t16 = t15 + y0;
211 const auto t17 = t16 * t5;
212 const auto t18 = t14 * t2 - t17;
213 const auto t19 = std::pow(t18, 2);
214 const auto t20 = -x2_rest;
215 const auto t21 = t20 + x0_rest;
216 const auto t22 = t0 * t21;
217 const auto t23 = 1.0 / t10;
218 const auto t24 = t23 * t9;
219 const auto t25 = t22 * t24;
220 const auto t26 = x0_rest - x1_rest;
221 const auto t27 = t26 * t4;
222 const auto t28 = t24 * t27;
223 const auto t29 = t3 + x1;
224 const auto t30 = t25 - t28 + t29;
225 const auto t31 = t14 * t21;
226 const auto t32 = t24 * t31;
227 const auto t33 = t15 + y1;
228 const auto t34 = t16 * t23 * t26 * t9 - t32 - t33;
229 const auto t35 = t22 - t27;
230 const auto t36 = t30 * t35;
231 const auto t37 = 2 * t24;
232 const auto t38 = t16 * t26 - t31;
233 const auto t39 = t34 * t38;
234 const auto t40 = t19 * t24 + t24 * t8 + t36 + t39;
235 const auto t41 = x0 * y1 - x0 * y2 - x1 * y0 + x1 * y2 + x2 * y0 - x2 * y1;
236 const auto t42 = 2 / t41;
237 const auto t43 = t23 * t42;
238 const auto t44 = t23 * t41 > 0;
239 const auto t45 = t20 + x1_rest;
240 const auto t46 = t23 * t45;
241 const auto t47 = t46 * t6;
242 const auto t48 = t0 * t2 * t23 * t45 - t29 - t47;
243 const auto t49 = t48 * t7;
244 const auto t50 = t17 * t46;
245 const auto t51 = t14 * t2 * t23 * t45 - t33 - t50;
246 const auto t52 = t18 * t51;
247 const auto t53 = std::pow(t38, 2);
248 const auto t54 = t46 * t53 + t49 + t52;
249 const auto t55 = ((t44) ? (
250 t11 * t42 * (-t18 * t9 * (2 * t14 * t2 * t23 * t45 - t33 - 2 * t50) - t35 * t45 * (2 * t25 - 2 * t28 + t29) - t36 * t45 - t38 * t45 * (2 * t16 * t23 * t26 * t9 - 2 * t32 - t33) - t39 * t45 + t40 * t45 - t49 * t9 - t52 * t9 - t7 * t9 * (2 * t0 * t2 * t23 * t45 - t29 - 2 * t47) + t9 * (std::pow(t35, 2) * t46 + t54)))
251 : (
252 0));
253 const auto t56 = std::pow(t35, 2);
254 const auto t57 = t12 * std::pow(t45, 2);
255 const auto t58 = 2 * t46;
256 if (t44)
257 {
258 H[0] = t43 * (t13 * t19 + t13 * t8 + std::pow(t30, 2) + std::pow(t34, 2) + t36 * t37 + t37 * t39 - t37 * t40);
259 }
260 else
261 {
262 H[0] = 0;
263 }
264 H[1] = t55;
265 H[2] = t55;
266 if (t44)
267 {
268 H[3] = t43 * (std::pow(t48, 2) + t49 * t58 + std::pow(t51, 2) + t52 * t58 + t53 * t57 + t56 * t57 - t58 * (t46 * t56 + t54));
269 }
270 else
271 {
272 H[3] = 0;
273 }
274 }
275 } // namespace autogen
276
277} // namespace polyfem
x
int x
Definition SplineBasis3d.cpp:55
polyfem::solver::AMIPSForm::value_unweighted
double value_unweighted(const Eigen::VectorXd &x) const override
Compute the value of the form.
Definition AMIPSForm.cpp:47
polyfem::solver::AMIPSForm::energy
static double energy(const Eigen::Vector2d &x0_rest, const Eigen::Vector2d &x1_rest, const Eigen::Vector2d &x2_rest, const Eigen::Vector2d &x0, const Eigen::Vector2d &x1, const Eigen::Vector2d &x2)
Definition AMIPSForm.cpp:30
polyfem::solver::AMIPSForm::is_step_valid
bool is_step_valid(const Eigen::VectorXd &, const Eigen::VectorXd &x1) const override
Determine if a step from solution x0 to solution x1 is allowed.
Definition AMIPSForm.cpp:16
polyfem::solver::AMIPSForm::second_derivative_unweighted
void second_derivative_unweighted(const Eigen::VectorXd &x, StiffnessMatrix &hessian) const override
Compute the second derivative of the value wrt x.
Definition AMIPSForm.cpp:89
polyfem::solver::AMIPSForm::X_rest_
const Eigen::MatrixXd X_rest_
Definition AMIPSForm.hpp:44
polyfem::solver::AMIPSForm::X_
const Eigen::MatrixXd X_
Definition AMIPSForm.hpp:45
polyfem::solver::AMIPSForm::first_derivative_unweighted
void first_derivative_unweighted(const Eigen::VectorXd &x, Eigen::VectorXd &gradv) const override
Compute the first derivative of the value wrt x.
Definition AMIPSForm.cpp:66
polyfem::solver::AMIPSForm::AMIPSForm
AMIPSForm(const Eigen::MatrixXd X_rest, const Eigen::MatrixXd X)
Assumes the variable to optimize over is X_rest.row(0)
Definition AMIPSForm.cpp:9
polyfem::autogen::AMIPS2D_hessian
void AMIPS2D_hessian(double x0_rest, double y0_rest, double x1_rest, double y1_rest, double x2_rest, double y2_rest, double x0, double y0, double x1, double y1, double x2, double y2, double H[4])
Definition AMIPSForm.cpp:192
polyfem::autogen::AMIPS2D_gradient
void AMIPS2D_gradient(double x0_rest, double y0_rest, double x1_rest, double y1_rest, double x2_rest, double y2_rest, double x0, double y0, double x1, double y1, double x2, double y2, double g[2])
Definition AMIPSForm.cpp:134
polyfem::autogen::AMIPS2D_energy
double AMIPS2D_energy(double x0_rest, double y0_rest, double x1_rest, double y1_rest, double x2_rest, double y2_rest, double x0, double y0, double x1, double y1, double x2, double y2)
Definition AMIPSForm.cpp:116
polyfem::StiffnessMatrix
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:20
AMIPSForm.hpp