22 using typename polysolve::nonlinear::Problem::Scalar;
23 using typename polysolve::nonlinear::Problem::THessian;
24 using typename polysolve::nonlinear::Problem::TVector;
30 double value(
const TVector& x)
override;
31 void gradient(
const TVector& x, TVector& gradv)
override;
34 throw std::runtime_error(
"Sparse functions do not exist, use dense solver");
40 bool is_step_valid(
const TVector& x0,
const TVector& x1)
override;
43 std::vector<std::array<double, S>>
m_cells;
145 Eigen::MatrixXd& hessian)
148 constexpr int64_t size = S == 6 ? 2 : 3;
149 hessian.resize(size, size);
151 Eigen::Matrix<double, size, size> tmp;
154 Eigen::Matrix<double, size, 1> tmpj(size);
157 for (
auto c : m_cells) {
158 if constexpr (S == 6) {
159 assert(x.size() == 2);
168 assert(x.size() == 3);
179 hessian += (p) * (p - 1) * std::pow(tmp2, p - 2) * tmpj * tmpj.transpose() +
180 p * std::pow(tmp2, p - 1) * tmp;
271 assert(child_simplices.size() > 0);
272 if (child_simplices.size() == 1) {
276 mesh().top_simplex_type());
281 std::vector<std::array<double, 12>> cells;
288 cell.primitive_type(),
293 const auto vertices =
297 assert(vertices.size() == 4);
300 vertices.simplex_vector()[0],
302 std::cout <<
"error here" << std::endl;
306 std::array<double, 12> single_cell;
307 std::vector<Vector3d> ps;
308 for (
size_t i = 0; i < 4; ++i) {
311 const auto p = accessor.const_vector_attribute(v);
313 single_cell[3 * i + 0] = p[0];
314 single_cell[3 * i + 1] = p[1];
315 single_cell[3 * i + 2] = p[2];
321 cells.emplace_back(single_cell);
328 Eigen::Vector3d x1(cells[0][0], cells[0][1], cells[0][2]);
330 std::cout <<
"step is not valid!!!!!!!!!!!!!!!!" << std::endl;
340 }
catch (
const std::exception&) {
359 accessor.vector_attribute(simplex.
tuple())[d] = (1 - alpha) * x0[d] + alpha * x[d];
366 }
else if (child_simplices.size() == 2) {
376 auto pos_0 = accessor.const_vector_attribute(child_simplices[0].tuple());
377 auto pos_1 = accessor.const_vector_attribute(child_simplices[1].tuple());
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])) {
388 if (
abs(dist[1]) >
abs(dist[2])) {
398 std::vector<std::array<double, 12>> cells;
405 cell.primitive_type(),
410 const auto vertices =
414 assert(vertices.size() == 4);
416 std::array<double, 12> single_cell;
417 std::vector<Vector3d> ps;
418 for (
size_t i = 0; i < 4; ++i) {
421 const auto p = accessor.const_vector_attribute(v);
423 single_cell[3 * i + 0] = p[0];
424 single_cell[3 * i + 1] = p[1];
425 single_cell[3 * i + 2] = p[2];
431 cells.emplace_back(single_cell);
441 cell.primitive_type(),
446 const auto vertices =
450 assert(vertices.size() == 4);
452 std::array<double, 12> single_cell;
453 std::vector<Vector3d> ps;
454 for (
size_t i = 0; i < 4; ++i) {
457 const auto p = accessor.const_vector_attribute(v);
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;
467 cells.emplace_back(single_cell);
476 Eigen::Vector3d x1(cells[0][0], cells[0][1], cells[0][2]);
478 std::cout <<
"step is not valid!!!!!!!!!!!!!!!!" << std::endl;
488 }
catch (
const std::exception&) {
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);
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;