Wildmeshing Toolkit
Loading...
Searching...
No Matches
AMIPSOptimizationSmoothingPeriodic.cpp
Go to the documentation of this file.
2
3#include <polysolve/nonlinear/Problem.hpp>
4#include <wmtk/Mesh.hpp>
5#include <wmtk/Types.hpp>
10#include <wmtk/utils/orient.hpp>
11
12#include <polysolve/nonlinear/Solver.hpp>
13
14#include <Eigen/Core>
15
16namespace wmtk::operations {
17
18template <int S>
19class AMIPSOptimizationSmoothingPeriodic::WMTKAMIPSProblem : public polysolve::nonlinear::Problem
20{
21public:
22 using typename polysolve::nonlinear::Problem::Scalar;
23 using typename polysolve::nonlinear::Problem::THessian;
24 using typename polysolve::nonlinear::Problem::TVector;
25
26 WMTKAMIPSProblem(std::vector<std::array<double, S>>& cells);
27
28 TVector initial_value() const;
29
30 double value(const TVector& x) override;
31 void gradient(const TVector& x, TVector& gradv) override;
32 void hessian(const TVector& x, THessian& hessian) override
33 {
34 throw std::runtime_error("Sparse functions do not exist, use dense solver");
35 }
36 void hessian(const TVector& x, Eigen::MatrixXd& hessian) override;
37
38 void solution_changed(const TVector& new_x) override;
39
40 bool is_step_valid(const TVector& x0, const TVector& x1) override;
41
42private:
43 std::vector<std::array<double, S>> m_cells;
44// Enables the use of power of p to approximate minimizing max energy
45// #define WMTK_ENABLE_P
46#ifdef WMTK_ENABLE_P
47 const int p = 4;
48#endif
49};
50
51template <int S>
53 std::vector<std::array<double, S>>& cells)
54 : m_cells(cells)
55{}
56
57template <int S>
58Eigen::Matrix<double, -1, 1>
60{
61 constexpr int64_t size = S == 6 ? 2 : 3;
62 Eigen::Matrix<double, size, 1> tmp;
63
64 for (int64_t d = 0; d < size; ++d) {
65 tmp(d) = m_cells[0][d];
66 }
67
68 return tmp;
69}
70
71template <int S>
73{
74 double res = 0;
75 for (auto c : m_cells) {
76 if constexpr (S == 6) {
77 assert(x.size() == 2);
78 c[0] = x[0];
79 c[1] = x[1];
80#ifdef WMTK_ENABLE_P
81 res += std::pow(wmtk::function::Tri_AMIPS_energy(c), p);
82#else
84#endif
85 } else {
86 assert(x.size() == 3);
87 c[0] = x[0];
88 c[1] = x[1];
89 c[2] = x[2];
90#ifdef WMTK_ENABLE_P
91 res += std::pow(wmtk::function::Tet_AMIPS_energy(c), p);
92#else
94#endif
95 }
96 }
97
98 // std::cout << res << std::endl;
99 return res;
100}
101
102template <int S>
104 const TVector& x,
105 TVector& gradv)
106{
107 constexpr int64_t size = S == 6 ? 2 : 3;
108 gradv.resize(size);
109 gradv.setZero();
110 Eigen::Matrix<double, size, 1> tmp(size);
111#ifdef WMTK_ENABLE_P
112 double tmp2 = 0;
113#endif
114
115 for (auto c : m_cells) {
116 if constexpr (S == 6) {
117 assert(x.size() == 2);
118 c[0] = x[0];
119 c[1] = x[1];
121#ifdef WMTK_ENABLE_P
123#endif
124 } else {
125 assert(x.size() == 3);
126 c[0] = x[0];
127 c[1] = x[1];
128 c[2] = x[2];
130#ifdef WMTK_ENABLE_P
132#endif
133 }
134#ifdef WMTK_ENABLE_P
135 gradv += double(p) * pow(tmp2, p - 1) * tmp;
136#else
137 gradv += tmp;
138#endif
139 }
140}
141
142template <int S>
144 const TVector& x,
145 Eigen::MatrixXd& hessian)
146{
147 // std::cout << "error here !!! Hessian used" << std::endl;
148 constexpr int64_t size = S == 6 ? 2 : 3;
149 hessian.resize(size, size);
150 hessian.setZero();
151 Eigen::Matrix<double, size, size> tmp;
152#ifdef WMTK_ENABLE_P
153 double tmp2 = 0;
154 Eigen::Matrix<double, size, 1> tmpj(size);
155#endif
156
157 for (auto c : m_cells) {
158 if constexpr (S == 6) {
159 assert(x.size() == 2);
160 c[0] = x[0];
161 c[1] = x[1];
163#ifdef WMTK_ENABLE_P
166#endif
167 } else {
168 assert(x.size() == 3);
169 c[0] = x[0];
170 c[1] = x[1];
171 c[2] = x[2];
173#ifdef WMTK_ENABLE_P
176#endif
177 }
178#ifdef WMTK_ENABLE_P
179 hessian += (p) * (p - 1) * std::pow(tmp2, p - 2) * tmpj * tmpj.transpose() +
180 p * std::pow(tmp2, p - 1) * tmp;
181#else
182 hessian += tmp;
183#endif
184 // 12 f(x)^2 * f(x)'*f^T'(x) + 4 f(x)^3 * H
185 }
186}
187
188template <int S>
191
192template <int S>
194 const TVector& x0,
195 const TVector& x1)
196{
197 constexpr int64_t size = S == 6 ? 2 : 3;
198 Eigen::Matrix<double, size, 1> p0 = x1, p1, p2, p3;
199 if constexpr (S == 6) {
200 for (const auto& c : m_cells) {
201 p1 << c[2], c[3];
202 p2 << c[4], c[5];
203
204 if (wmtk::utils::wmtk_orient2d(p0, p1, p2) <= 0) return false;
205 }
206 } else {
207 for (const auto& c : m_cells) {
208 p1 << c[3], c[4], c[5];
209 p2 << c[6], c[7], c[8];
210 p3 << c[9], c[10], c[11];
211
212 // Eigen::Matrix<double, size, 1> x0m = x0;
213 // assert(wmtk::utils::wmtk_orient3d(x0m, p1, p2, p3) > 0);
214 // if (wmtk::utils::wmtk_orient3d(p0, p1, p2, p3) <= 0) {
215 // // std::cout << wmtk::utils::wmtk_orient3d(x0m, p1, p2, p3) << std::endl;
216 // return false;
217 // }
218
219 Eigen::Matrix<double, size, 1> x0m = x0;
220 assert(wmtk::utils::wmtk_orient3d(p3, x0m, p1, p2) > 0);
221 if (wmtk::utils::wmtk_orient3d(p3, p0, p1, p2) <= 0) {
222 // std::cout << wmtk::utils::wmtk_orient3d(p3, x0m, p1, p2) << std::endl;
223 return false;
224 }
225 }
226 }
227 return true;
228}
229
230
232 Mesh& periodic_mesh,
233 Mesh& position_mesh,
234 const attribute::MeshAttributeHandle& coords)
235 : AttributesUpdate(position_mesh)
236 , m_periodic_mesh(periodic_mesh)
237 , m_position_mesh(position_mesh)
238 , m_coordinate_handle(coords)
239 , m_amips(position_mesh, coords)
240{
241 // assert(m_coordinate_handle.holds<Rational>());
242
243 m_linear_solver_params = R"({"solver": "Eigen::LDLT"})"_json;
244 m_nonlinear_solver_params = R"({"solver": "DenseNewton", "max_iterations": 10})"_json;
245 // m_nonlinear_solver_params = R"({"solver": "L-BFGS", "max_iterations": 100,
246 // "advanced":{"apply_gradient_fd": "FullFiniteDiff"}})"_json;
247 // m_nonlinear_solver_params = R"({"solver": "GradientDescent", "max_iterations": 100})"_json;
248
249
251}
252
254{
255 m_solver = polysolve::nonlinear::Solver::create(
258 1,
259 opt_logger());
260}
261
262
264 const simplex::Simplex& simplex)
265{
266 auto accessor = m_position_mesh.create_accessor(m_coordinate_handle.as<double>());
267
268 auto parent_simplex = m_position_mesh.map_to_parent(simplex);
269 auto child_simplices = m_periodic_mesh.map_to_child(m_position_mesh, parent_simplex);
270
271 assert(child_simplices.size() > 0);
272 if (child_simplices.size() == 1) {
274 mesh(),
275 simplex,
276 mesh().top_simplex_type());
277
278
279 assert(mesh().top_simplex_type() == PrimitiveType::Tetrahedron);
280
281 std::vector<std::array<double, 12>> cells;
282
283 for (simplex::Simplex cell : neighs) {
284 if (!mesh().is_ccw(cell.tuple())) {
285 // switch any local id but NOT the vertex
286 cell = simplex::Simplex(
287 mesh(),
288 cell.primitive_type(),
289 mesh().switch_tuple(cell.tuple(), PrimitiveType::Edge));
290 }
291 assert(mesh().is_ccw(cell.tuple()));
292
293 const auto vertices =
295
296
297 assert(vertices.size() == 4);
299 mesh(),
300 vertices.simplex_vector()[0],
301 simplex)) {
302 std::cout << "error here" << std::endl;
303 }
304
305
306 std::array<double, 12> single_cell;
307 std::vector<Vector3d> ps;
308 for (size_t i = 0; i < 4; ++i) {
309 const simplex::Simplex& v = vertices.simplex_vector()[i];
310 // const auto p = accessor.const_vector_attribute(vertices[i]);
311 const auto p = accessor.const_vector_attribute(v);
312 ps.push_back(p);
313 single_cell[3 * i + 0] = p[0];
314 single_cell[3 * i + 1] = p[1];
315 single_cell[3 * i + 2] = p[2];
316 }
317 // if (wmtk::utils::wmtk_orient3d(ps[3], ps[0], ps[1], ps[2]) <= 0) {
318 // std::cout << "this is wrong" << std::endl;
319 // }
320
321 cells.emplace_back(single_cell);
322
323 // cells.emplace_back(m_amips.get_raw_coordinates<4, 3>(cell, simplex));
324 }
325 WMTKAMIPSProblem<12> problem(cells);
326
327 {
328 Eigen::Vector3d x1(cells[0][0], cells[0][1], cells[0][2]);
329 if (!problem.is_step_valid(x1, x1)) {
330 std::cout << "step is not valid!!!!!!!!!!!!!!!!" << std::endl;
331 }
332 }
333
334 auto x = problem.initial_value();
335 auto x0 = problem.initial_value();
336
337 try {
338 m_solver->minimize(problem, x);
339
340 } catch (const std::exception&) {
341 }
342
343 // Hack for surface only
344
345 auto child_meshes = mesh().get_child_meshes();
346
347
348 double alpha = 1.00;
349 // for (auto child_mesh : child_meshes) {
350 // if (!mesh().map_to_child(*child_mesh, simplex).empty()) {
351 // alpha = 0.01;
352 // break;
353 // }
354 // }
355
356
357 for (int64_t d = 0; d < m_coordinate_handle.dimension(); ++d) {
358 // accessor.vector_attribute(simplex.tuple())[d] = Rational(x[d], true);
359 accessor.vector_attribute(simplex.tuple())[d] = (1 - alpha) * x0[d] + alpha * x[d];
360 }
361
362
363 // assert(attribute_handle() == m_function.attribute_handle());
364
365 double res = 0;
366 } else if (child_simplices.size() == 2) {
369 child_simplices[0],
373 child_simplices[1],
375
376 auto pos_0 = accessor.const_vector_attribute(child_simplices[0].tuple());
377 auto pos_1 = accessor.const_vector_attribute(child_simplices[1].tuple());
378
379 double space_x = 0, space_y = 0, space_z = 0;
380 auto dist = pos_1 - pos_0;
381 if (abs(dist[0]) > abs(dist[1])) {
382 if (abs(dist[0]) > abs(dist[2])) {
383 space_x = dist[0];
384 } else {
385 space_z = dist[2];
386 }
387 } else {
388 if (abs(dist[1]) > abs(dist[2])) {
389 space_y = dist[1];
390 } else {
391 space_z = dist[2];
392 }
393 }
394
395
396 assert(mesh().top_simplex_type() == PrimitiveType::Tetrahedron);
397
398 std::vector<std::array<double, 12>> cells;
399
400 for (simplex::Simplex cell : neighs_0) {
401 if (!mesh().is_ccw(cell.tuple())) {
402 // switch any local id but NOT the vertex
403 cell = simplex::Simplex(
404 mesh(),
405 cell.primitive_type(),
406 mesh().switch_tuple(cell.tuple(), PrimitiveType::Edge));
407 }
408 assert(mesh().is_ccw(cell.tuple()));
409
410 const auto vertices =
412
413
414 assert(vertices.size() == 4);
415
416 std::array<double, 12> single_cell;
417 std::vector<Vector3d> ps;
418 for (size_t i = 0; i < 4; ++i) {
419 const simplex::Simplex& v = vertices.simplex_vector()[i];
420 // const auto p = accessor.const_vector_attribute(vertices[i]);
421 const auto p = accessor.const_vector_attribute(v);
422 ps.push_back(p);
423 single_cell[3 * i + 0] = p[0];
424 single_cell[3 * i + 1] = p[1];
425 single_cell[3 * i + 2] = p[2];
426 }
427 // if (wmtk::utils::wmtk_orient3d(ps[3], ps[0], ps[1], ps[2]) <= 0) {
428 // std::cout << "this is wrong" << std::endl;
429 // }
430
431 cells.emplace_back(single_cell);
432
433 // cells.emplace_back(m_amips.get_raw_coordinates<4, 3>(cell, simplex));
434 }
435
436 for (simplex::Simplex cell : neighs_1) {
437 if (!mesh().is_ccw(cell.tuple())) {
438 // switch any local id but NOT the vertex
439 cell = simplex::Simplex(
440 mesh(),
441 cell.primitive_type(),
442 mesh().switch_tuple(cell.tuple(), PrimitiveType::Edge));
443 }
444 assert(mesh().is_ccw(cell.tuple()));
445
446 const auto vertices =
448
449
450 assert(vertices.size() == 4);
451
452 std::array<double, 12> single_cell;
453 std::vector<Vector3d> ps;
454 for (size_t i = 0; i < 4; ++i) {
455 const simplex::Simplex& v = vertices.simplex_vector()[i];
456 // const auto p = accessor.const_vector_attribute(vertices[i]);
457 const auto p = accessor.const_vector_attribute(v);
458 ps.push_back(p);
459 single_cell[3 * i + 0] = p[0] - space_x;
460 single_cell[3 * i + 1] = p[1] - space_y;
461 single_cell[3 * i + 2] = p[2] - space_z;
462 }
463 // if (wmtk::utils::wmtk_orient3d(ps[3], ps[0], ps[1], ps[2]) <= 0) {
464 // std::cout << "this is wrong" << std::endl;
465 // }
466
467 cells.emplace_back(single_cell);
468
469 // cells.emplace_back(m_amips.get_raw_coordinates<4, 3>(cell, simplex));
470 }
471
472
473 WMTKAMIPSProblem<12> problem(cells);
474
475 {
476 Eigen::Vector3d x1(cells[0][0], cells[0][1], cells[0][2]);
477 if (!problem.is_step_valid(x1, x1)) {
478 std::cout << "step is not valid!!!!!!!!!!!!!!!!" << std::endl;
479 }
480 }
481
482 auto x = problem.initial_value();
483 auto x0 = problem.initial_value();
484
485 try {
486 m_solver->minimize(problem, x);
487
488 } catch (const std::exception&) {
489 }
490
491 // Hack for surface only
492
493
494 double alpha = 1.00;
495 // for (auto child_mesh : child_meshes) {
496 // if (!mesh().map_to_child(*child_mesh, simplex).empty()) {
497 // alpha = 0.01;
498 // break;
499 // }
500 // }
501
502 accessor.vector_attribute(child_simplices[0].tuple()) = x;
503 accessor.vector_attribute(child_simplices[1].tuple()) =
504 x + Vector3d(space_x, space_y, space_z);
505
506 // projection
507 if (space_x > 0) {
508 accessor.vector_attribute(child_simplices[0].tuple())[0] = pos_0[0];
509 accessor.vector_attribute(child_simplices[1].tuple())[0] = pos_0[0] + space_x;
510 } else if (space_y > 0) {
511 accessor.vector_attribute(child_simplices[0].tuple())[1] = pos_0[1];
512 accessor.vector_attribute(child_simplices[1].tuple())[1] = pos_0[1] + space_y;
513 } else if (space_z > 0) {
514 accessor.vector_attribute(child_simplices[0].tuple())[2] = pos_0[2];
515 accessor.vector_attribute(child_simplices[1].tuple())[2] = pos_0[2] + space_z;
516 }
517
518 // assert(attribute_handle() == m_function.attribute_handle());
519
520 double res = 0;
521 } else {
522 // throw std::runtime_error("not handled case on bbox edge and corner");
523 }
524
525
526 return AttributesUpdate::execute(simplex);
527}
528
529} // namespace wmtk::operations
virtual bool is_ccw(const Tuple &tuple) const =0
TODO this needs dimension?
attribute::Accessor< T, Mesh, D > create_accessor(const attribute::MeshAttributeHandle &handle)
simplex::Simplex map_to_parent(const simplex::Simplex &my_simplex) const
optimized map from a simplex from this mesh to its direct parent
std::vector< simplex::Simplex > map_to_child(const Mesh &child_mesh, const simplex::Simplex &my_simplex) const
optimized map fromsimplex from this mesh to one of its direct children
std::vector< std::shared_ptr< Mesh > > get_child_meshes() const
returns the direct multimesh child meshes for the current mesh
auto as() const -> const held_handle_type< held_type_from_primitive< T >()> &
AMIPSOptimizationSmoothingPeriodic(Mesh &periodic_mesh, Mesh &position_mesh, const attribute::MeshAttributeHandle &coords)
std::vector< simplex::Simplex > execute(const simplex::Simplex &simplex) override
returns an empty vector in case of failure
virtual std::vector< simplex::Simplex > execute(const simplex::Simplex &simplex) override
returns an empty vector in case of failure
const Mesh & mesh() const
Definition Operation.hpp:45
const Tuple & tuple() const
Definition Simplex.hpp:53
static bool equal(const Mesh &m, const Simplex &s0, const Simplex &s1)
bool is_ccw(PrimitiveType pt, const Tuple &t)
Definition is_ccw.cpp:8
void Tet_AMIPS_jacobian(const std::array< double, 12 > &T, Eigen::Vector3d &result_0)
Definition AMIPS.cpp:75
void Tri_AMIPS_hessian(const std::array< double, 6 > &T, Eigen::Matrix2d &result_0)
Definition AMIPS.cpp:542
void Tri_AMIPS_jacobian(const std::array< double, 6 > &T, Eigen::Vector2d &result_0)
Definition AMIPS.cpp:503
double Tet_AMIPS_energy(const std::array< double, 12 > &T)
Definition AMIPS.cpp:375
void Tet_AMIPS_hessian(const std::array< double, 12 > &T, Eigen::Matrix3d &result_0)
Definition AMIPS.cpp:172
double Tri_AMIPS_energy(const std::array< double, 6 > &T)
Definition AMIPS.cpp:473
std::vector< Simplex > cofaces_single_dimension_simplices(const Mesh &mesh, const Simplex &simplex, PrimitiveType cofaces_type)
SimplexCollection faces_single_dimension(const Mesh &mesh, const Simplex &simplex, const PrimitiveType face_type)
Returns a vector with all faces in the boundary of a simplex of the given dimension.
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
int wmtk_orient2d(double p0x, double p0y, double p1x, double p1y, double p2x, double p2y)
Definition orient.cpp:178
Rational abs(const Rational &r0)
Definition Rational.cpp:328
Vector< double, 3 > Vector3d
Definition Types.hpp:39
spdlog::logger & opt_logger()
Retrieves the logger for the optimization.
Definition Logger.cpp:75
Rational pow(const Rational &x, int p)
Definition Rational.cpp:166