Merge branch 'Progressive_Embedding' into 'master'

Progressive embedding

See merge request !53

CMakeLists.txt

@@ -491,85 +491,89 @@ configure_file ("${CMAKE_CURRENT_SOURCE_DIR}/Config/config.hh.in"
 # of the library are already there
 #######################################################################
 
-set (COMISO_BUILD_EXAMPLES TRUE CACHE BOOL "Build CoMISo Examples")
+if(${PROJECT_NAME} MATCHES "CoMISo")
+  set (COMISO_BUILD_EXAMPLES TRUE CACHE BOOL "Build CoMISo Examples")
+else()
+  set (COMISO_BUILD_EXAMPLES FALSE CACHE BOOL "Build CoMISo Examples")
+endif()
 
 if (COMISO_BUILD_EXAMPLES )
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/factored_solver/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/factored_solver/CMakeLists.txt" )
     add_subdirectory (Examples/factored_solver)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/quadratic_solver/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/quadratic_solver/CMakeLists.txt" )
     add_subdirectory (Examples/quadratic_solver)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/test2/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/test2/CMakeLists.txt" )
     add_subdirectory (Examples/test2)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_quadratic_example/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_quadratic_example/CMakeLists.txt" )
     add_subdirectory (Examples/small_quadratic_example)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_factored_example/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_factored_example/CMakeLists.txt" )
     add_subdirectory (Examples/small_factored_example)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/super_sparse_matrix/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/super_sparse_matrix/CMakeLists.txt" )
     add_subdirectory (Examples/super_sparse_matrix)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/eigen_solver/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/eigen_solver/CMakeLists.txt" )
     add_subdirectory (Examples/eigen_solver)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_nsolver/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_nsolver/CMakeLists.txt" )
     add_subdirectory (Examples/small_nsolver)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_eigenproblem/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_eigenproblem/CMakeLists.txt" )
     add_subdirectory (Examples/small_eigenproblem)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_miqp/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_miqp/CMakeLists.txt" )
     add_subdirectory (Examples/small_miqp)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_nleast_squares/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_nleast_squares/CMakeLists.txt" )
     add_subdirectory (Examples/small_nleast_squares)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_sparseqr/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_sparseqr/CMakeLists.txt" )
     add_subdirectory (Examples/small_sparseqr)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_quadratic_resolve_example/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_quadratic_resolve_example/CMakeLists.txt" )
     add_subdirectory (Examples/small_quadratic_resolve_example)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_cplex_soc/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_cplex_soc/CMakeLists.txt" )
     add_subdirectory (Examples/small_cplex_soc)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_adolc/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_adolc/CMakeLists.txt" )
     add_subdirectory (Examples/small_adolc)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_dco/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_dco/CMakeLists.txt" )
     add_subdirectory (Examples/small_dco)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/vector1_adolc/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/vector1_adolc/CMakeLists.txt" )
     add_subdirectory (Examples/vector1_adolc)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_linear_problem/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_linear_problem/CMakeLists.txt" )
     add_subdirectory (Examples/small_linear_problem)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/crossfield3d/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/crossfield3d/CMakeLists.txt" )
     add_subdirectory (Examples/crossfield3d)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/crossfield3d/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/crossfield3d/CMakeLists.txt" )
     add_subdirectory (Examples/crossfield3d_dco)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_finite_element/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_finite_element/CMakeLists.txt" )
     add_subdirectory (Examples/small_finite_element)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_AQP/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_AQP/CMakeLists.txt" )
     add_subdirectory (Examples/small_AQP)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/finite_element_integrability_problem/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/finite_element_integrability_problem/CMakeLists.txt" )
     add_subdirectory (Examples/finite_element_integrability_problem)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_mosek_native/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_mosek_native/CMakeLists.txt" )
     add_subdirectory (Examples/small_mosek_native)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_mosek_fusion_sdp/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_mosek_fusion_sdp/CMakeLists.txt" )
     add_subdirectory (Examples/small_mosek_fusion_sdp)
   endif()
-  if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_symmetric_dirichlet/CMakeLists.txt" )
+  if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Examples/small_symmetric_dirichlet/CMakeLists.txt" )
   add_subdirectory (Examples/small_symmetric_dirichlet)
   endif()
 

Examples/small_symmetric_dirichlet/main.cc

@@ -34,9 +34,10 @@
 
 //------------------------------------------------------------------------------------------------------
 
-// Example main
-int main(void)
+
+void symmetric_dirichlet_problem_example()
 {
+  std::cout << "---------- Example of symmetric dirichlet problem:" << std::endl;
   std::cout << "---------- 1) Get an instance of a SymmetricDirichletProblem..." << std::endl;
 
   // then create finite element problem and add sets
@@ -72,7 +73,139 @@ int main(void)
 
   // print result
   for (unsigned int i = 0; i < n_vertices; ++i)
-    std::cerr << "p" << i << " = ( " << sd_problem.x()[2*i+0] << ", " << sd_problem.x()[2*i+1] << ")"  << std::endl;
+    std::cout << "p" << i << " = ( " << sd_problem.x()[2*i+0] << ", " << sd_problem.x()[2*i+1] << ")"  << std::endl;
+}
+
+
+void symmetric_dirichlet_problem_example_one_ring_with_constraints(bool verbose = true)
+{
+  if (verbose)
+  {
+    std::cout << "---------- Example of symmetric dirichlet problem of a one ring with fixed boundary:" << std::endl;
+    std::cout << "---------- 1) Get an instance of a SymmetricDirichletProblem..." << std::endl;
+  }
+
+  // then create finite element problem and add sets
+  unsigned int n_vertices = 7;
+  COMISO::SymmetricDirichletProblem sd_problem(n_vertices);
+  auto reference_positions = sd_problem.get_equilateral_refernce_positions();
+  for (int i = 0; i < 6; ++i)
+  {
+    COMISO::SymmetricDirichletProblem::IndexVector indices(0,i+1,(i+2-1)%6+1);
+    sd_problem.add_triangle(indices, reference_positions);
+  }
+
+  std::vector<double> initial_solution;
+  // center vertex position
+  initial_solution.push_back(-0.5);
+  initial_solution.push_back( 0.3);
+
+  Eigen::Rotation2D<double> rot(60.0/180.0*M_PI);
+  Eigen::Vector2d pos{1.0,0.0};
+  for (int i = 0; i < 6; ++i)
+  {
+    initial_solution.push_back(pos(0));
+    initial_solution.push_back(pos(1));
+    pos = rot * pos;
+  }
+
+  sd_problem.x() = initial_solution;
+
+  if (verbose)
+    std::cout << "---------- 2) Set up constraints..." << std::endl;
+  // fix all boundary positions
+  for (int i = 1; i < 7; ++i)
+    sd_problem.add_fix_point_constraint(i, initial_solution[2*i], initial_solution[2*i+1]);
+
+  if (verbose)
+    std::cout << "---------- 3) Solve with Newton Solver..." << std::endl;
+  COMISO::SymmetricDirichletProblem::VectorD b;
+  COMISO::SymmetricDirichletProblem::SMatrixD A;
+  sd_problem.get_constraints(A, b);
+  COMISO::NewtonSolver nsolver;
+  nsolver.solve(&sd_problem, A, b);
+
+  // print result
+  if (verbose)
+    for (unsigned int i = 0; i < n_vertices; ++i)
+      std::cout << "p" << i << " = ( " << sd_problem.x()[2*i+0] << ", " << sd_problem.x()[2*i+1] << ")"  << std::endl;
+}
+
+
+void symmetric_dirichlet_problem_example_one_ring_with_specialized_problem(bool verbose = true)
+{
+  if (verbose)
+  {
+    std::cout << "---------- Example of specialized symmetric dirichlet problem for one ring optimization:" << std::endl;
+    std::cout << "---------- 1) Get an instance of a SymmetricDirichletProblem..." << std::endl;
+  }
+
+  // then create finite element problem and add sets
+  unsigned int n_vertices = 7;
+  COMISO::SymmetricDirichletOneRingProblem sd_problem;
+  auto reference_positions = sd_problem.get_equilateral_refernce_positions();
+
+  std::vector<double> initial_solution;
+  // center vertex position
+  initial_solution.push_back(-0.5);
+  initial_solution.push_back( 0.3);
+
+  Eigen::Rotation2D<double> rot(60.0/180.0*M_PI);
+  Eigen::Vector2d pos{1.0,0.0};
+  for (int i = 0; i < 6; ++i)
+  {
+    initial_solution.push_back(pos(0));
+    initial_solution.push_back(pos(1));
+    pos = rot * pos;
+  }
+
+  for (int i = 1; i < 7; ++i)
+  {
+    COMISO::SymmetricDirichletOneRingProblem::InputPositionVector2D current_positions;
+    int id0 = 0;
+    int id1 = i;
+    int id2 = i%6+1;
+    current_positions << initial_solution[2*id0], initial_solution[2*id0+1],
+                         initial_solution[2*id1], initial_solution[2*id1+1],
+                         initial_solution[2*id2], initial_solution[2*id2+1];
+    sd_problem.add_triangle(current_positions, reference_positions);
+  }
+
+  // only first positions required because we only optimize center vertex
+  initial_solution.resize(2);
+  sd_problem.x() = initial_solution;
+
+  if (verbose)
+    std::cout << "---------- 3) Solve with Newton Solver..." << std::endl;
+  COMISO::NewtonSolver nsolver;
+  nsolver.solve(&sd_problem);
+
+  // print result
+
+  if (verbose)
+    std::cout << "p" << 0 << " = ( " << sd_problem.x()[2*0+0] << ", " << sd_problem.x()[2*0+1] << ")"  << std::endl;
+}
+
+// Example main
+int main(void)
+{
+  symmetric_dirichlet_problem_example();
+  symmetric_dirichlet_problem_example_one_ring_with_constraints();
+  symmetric_dirichlet_problem_example_one_ring_with_specialized_problem();
+
+  COMISO::StopWatch sw_constraints;
+  sw_constraints.start();
+  for (int i = 0; i < 1000; ++i)
+    symmetric_dirichlet_problem_example_one_ring_with_constraints(false);
+  double time_constraints = sw_constraints.stop();
+  std::cout << "Solve with constraints took " << time_constraints/1000.0 << " seconds" << std::endl;
+
+  COMISO::StopWatch sw_specialized_problem;
+  sw_specialized_problem.start();
+  for (int i = 0; i < 1000; ++i)
+    symmetric_dirichlet_problem_example_one_ring_with_specialized_problem(false);
+  double time_specialized_problem = sw_specialized_problem.stop();
+  std::cout << "Solve with specialized problem took " << time_specialized_problem/1000.0 << " seconds" << std::endl;
 
   return 0;
 }

NSolver/NewtonSolver.cc

@@ -189,7 +189,10 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
 
       if(!fact_ok || kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2)
       {
-        DEB_warning(2, "Numerical issues in KKT system"); 
+        DEB_warning(2, "Numerical issues in KKT system");
+        DEB_warning_if(!fact_ok, 2, "Factorization not ok");
+        DEB_warning_if(kkt_res2 > KKT_res_eps, 2, "KKT Residuum too high: " << kkt_res2);
+        DEB_warning_if(constraint_res2 > max_allowed_constraint_violation2, 2, "Constraint violation too high: " << constraint_res2);
         // alternate hessian and constraints regularization
         if(reg_iters % 2 == 0 || regularize_constraints >= regularize_constraints_limit)
         {

NSolver/SymmetricDirichletProblem.cc

@@ -334,6 +334,218 @@ void SymmetricDirichletProblem::get_constraints(SMatrixD& _A, VectorD& _b)
   _A.setFromTriplets(triplets.begin(), triplets.end());
 }
 
+SymmetricDirichletProblem::ReferencePositionVector2D SymmetricDirichletProblem::get_equilateral_refernce_positions(double _area)
+{
+  ReferencePositionVector2D equilateral_reference;
+  equilateral_reference << 0.0, 0.0,
+                           1.0, 0.0,
+                           0.5, 0.5*std::sqrt(3.0);
+  equilateral_reference *= _area / 0.5*0.5*std::sqrt(3.0);
+  return equilateral_reference;
+}
+
+
+
+double SymmetricDirichletOneVertexElement::eval_f(const VecV& _x, const VecC& _c)
+{
+  Matrix2x2d B;
+  B(0,0) = _c[0]-_x[0];
+  B(0,1) = _c[2]-_x[0];
+  B(1,0) = _c[1]-_x[1];
+  B(1,1) = _c[3]-_x[1];
+  Matrix2x2d Bin = B.inverse();
+
+  Matrix2x2d R;
+  R(0,0) = _c[4+2]-_c[4+0];
+  R(0,1) = _c[4+4]-_c[4+0];
+  R(1,0) = _c[4+3]-_c[4+1];
+  R(1,1) = _c[4+5]-_c[4+1];
+  Matrix2x2d Rin = R.inverse();
+
+  double area = 0.5 * R.determinant();
+  if (B.determinant() * area <= 0)
+  {
+    double res = std::numeric_limits<double>::max();
+    return res;
+  }
+
+  Matrix2x2d J   =  B * Rin;
+  Matrix2x2d Jin =  R * Bin;
+
+  double res = J.squaredNorm() + Jin.squaredNorm();
+
+  return area * (res - 4);
+}
+
+void SymmetricDirichletOneVertexElement::eval_gradient(const VecV& _x, const VecC& _c, VecV& _g)
+{
+  const double a_x = _x[0];
+  const double a_y = _x[1];
+  const double b_x = _c[0];
+  const double b_y = _c[1];
+  const double c_x = _c[2];
+  const double c_y = _c[3];
+  const double d_x = _c[4];
+  const double d_y = _c[5];
+  const double e_x = _c[6];
+  const double e_y = _c[7];
+  const double f_x = _c[8];
+  const double f_y = _c[9];
+
+  const double det = d_x*e_y - d_x*f_y - d_y*e_x + d_y*f_x + e_x*f_y - e_y*f_x;
+
+  _g[0] = (2*(-(d_y - f_y)*(b_x - c_x) + (e_y - f_y)*(a_x - c_x))*(e_y - f_y) + 2*(-(a_x - c_x)*(e_x - f_x) + (b_x - c_x)*(d_x - f_x))*(-e_x + f_x))/
+          std::pow(d_x*e_y - d_x*f_y - d_y*e_x + d_y*f_x + e_x*f_y - e_y*f_x,2)
+      -  2*(std::pow(-(a_y - c_y)*(e_x - f_x) + (b_y - c_y)*(d_x - f_x),2) +
+            std::pow(-(b_x - c_x)*(d_x - f_x) + (a_x - c_x)*(e_x - f_x),2) +
+            std::pow( (b_y - c_y)*(d_y - f_y) - (a_y - c_y)*(e_y - f_y),2) +
+            std::pow(-(d_y - f_y)*(b_x - c_x) + (e_y - f_y)*(a_x - c_x),2)) *
+            ( b_y - c_y)/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,3)
+      + (2*(-(b_x - c_x)*(d_x - f_x) + (a_x - c_x)*(e_x - f_x))*( e_x - f_x) +
+         2*(-(d_y - f_y)*(b_x - c_x) + (e_y - f_y)*(a_x - c_x))*( e_y - f_y))/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,2);
+
+  _g[1] = (2*( (a_y - c_y)*(e_y - f_y) - (b_y - c_y)*(d_y - f_y))*(e_y - f_y) + 2*(-(a_y - c_y)*(e_x - f_x) + (b_y - c_y)*(d_x - f_x))*(-e_x + f_x))/
+          std::pow(d_x*e_y - d_x*f_y - d_y*e_x + d_y*f_x + e_x*f_y - e_y*f_x,2)
+      - 2*(std::pow(-(a_y - c_y)*(e_x - f_x) + (b_y - c_y)*(d_x - f_x),2) +
+           std::pow(-(b_x - c_x)*(d_x - f_x) + (a_x - c_x)*(e_x - f_x),2) +
+           std::pow( (b_y - c_y)*(d_y - f_y) - (a_y - c_y)*(e_y - f_y),2) +
+           std::pow(-(d_y - f_y)*(b_x - c_x) + (e_y - f_y)*(a_x - c_x),2)) *
+           (-b_x + c_x)/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,3)
+      + (2*(-(a_y - c_y)*(e_x - f_x) + (b_y - c_y)*(d_x - f_x))*(-e_x + f_x) +
+         2*( (b_y - c_y)*(d_y - f_y) - (a_y - c_y)*(e_y - f_y))*(-e_y + f_y))/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,2);
+
+  // area weight
+  _g[0] *= 0.5*det;
+  _g[1] *= 0.5*det;
+}
+
+void SymmetricDirichletOneVertexElement::eval_hessian(const VecV& _x, const VecC& _c, std::vector<Triplet>& _triplets)
+{
+  const double a_x = _x[0];
+  const double a_y = _x[1];
+  const double b_x = _c[0];
+  const double b_y = _c[1];
+  const double c_x = _c[2];
+  const double c_y = _c[3];
+  const double d_x = _c[4];
+  const double d_y = _c[5];
+  const double e_x = _c[6];
+  const double e_y = _c[7];
+  const double f_x = _c[8];
+  const double f_y = _c[9];
+
+  const double det = d_x*e_y - d_x*f_y - d_y*e_x + d_y*f_x + e_x*f_y - e_y*f_x;
+
+  Eigen::MatrixXd H(2,2);
+  H(0,0) = (2*std::pow(e_y - f_y,2) + 2*std::pow(-e_x + f_x,2))/std::pow(d_x*e_y - d_x*f_y - d_y*e_x + d_y*f_x + e_x*f_y - e_y*f_x,2) + 6*(std::pow(-(a_y - c_y)*(e_x - f_x) + (b_y - c_y)*(d_x - f_x),2) + std::pow(-(b_x - c_x)*(d_x - f_x) + (a_x - c_x)*(e_x - f_x),2) + std::pow((b_y - c_y)*(d_y - f_y) - (a_y - c_y)*(e_y - f_y),2) + std::pow(-(d_y - f_y)*(b_x - c_x) + (e_y - f_y)*(a_x - c_x),2))*std::pow( b_y - c_y,2)/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,4) - 4*(2*(-(b_x - c_x)*(d_x - f_x) + (a_x - c_x)*(e_x - f_x))*( e_x - f_x) + 2*(-(d_y - f_y)*(b_x - c_x) + (e_y - f_y)*(a_x - c_x))*( e_y - f_y))*( b_y - c_y)/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,3) + (2*std::pow( e_x - f_x,2) + 2*std::pow( e_y - f_y,2))/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,2);
+  H(1,1) = (2*std::pow(e_y - f_y,2) + 2*std::pow(-e_x + f_x,2))/std::pow(d_x*e_y - d_x*f_y - d_y*e_x + d_y*f_x + e_x*f_y - e_y*f_x,2) + 6*(std::pow(-(a_y - c_y)*(e_x - f_x) + (b_y - c_y)*(d_x - f_x),2) + std::pow(-(b_x - c_x)*(d_x - f_x) + (a_x - c_x)*(e_x - f_x),2) + std::pow((b_y - c_y)*(d_y - f_y) - (a_y - c_y)*(e_y - f_y),2) + std::pow(-(d_y - f_y)*(b_x - c_x) + (e_y - f_y)*(a_x - c_x),2))*std::pow(-b_x + c_x,2)/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,4) - 4*(2*(-(a_y - c_y)*(e_x - f_x) + (b_y - c_y)*(d_x - f_x))*(-e_x + f_x) + 2*( (b_y - c_y)*(d_y - f_y) - (a_y - c_y)*(e_y - f_y))*(-e_y + f_y))*(-b_x + c_x)/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,3) + (2*std::pow(-e_x + f_x,2) + 2*std::pow(-e_y + f_y,2))/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,2);
+  H(1,0) = 6*(std::pow(-(a_y - c_y)*(e_x - f_x) + (b_y - c_y)*(d_x - f_x),2) + std::pow(-(b_x - c_x)*(d_x - f_x) + (a_x - c_x)*(e_x - f_x),2) + std::pow((b_y - c_y)*(d_y - f_y) - (a_y - c_y)*(e_y - f_y),2) + std::pow(-(d_y - f_y)*(b_x - c_x) + (e_y - f_y)*(a_x - c_x),2))*(b_y - c_y)*(-b_x + c_x)/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,4) - 2*(2*(-(a_y - c_y)*(e_x - f_x) + (b_y - c_y)*(d_x - f_x))*(-e_x + f_x) + 2*((b_y - c_y)*(d_y - f_y) - (a_y - c_y)*(e_y - f_y))*(-e_y + f_y))*(b_y - c_y)/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,3) - 2*(2*(-(b_x - c_x)*(d_x - f_x) + (a_x - c_x)*(e_x - f_x))*(e_x - f_x) + 2*(-(d_y - f_y)*(b_x - c_x) + (e_y - f_y)*(a_x - c_x))*(e_y - f_y))*(-b_x + c_x)/std::pow(a_x*b_y - a_x*c_y - a_y*b_x + a_y*c_x + b_x*c_y - b_y*c_x,3);
+  H(0,1) = H(1,0);
+  H *=  0.5 * det;
+
+  Eigen::MatrixXd Hspd(2,2);
+  project_hessian(H, Hspd, 1e-6);
+
+  _triplets.reserve(2*2);
+  for (int i = 0; i < 2; ++i)
+    for (int j = 0; j < 2; ++j)
+      _triplets.push_back(Triplet(i,j,Hspd(i,j)));
+
+}
+
+void SymmetricDirichletOneVertexElement::project_hessian(Eigen::MatrixXd& H_orig, Eigen::MatrixXd& H_spd, double eps)
+{
+  // Compute eigen-decomposition (of symmetric matrix)
+  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eig(H_orig);
+  Eigen::MatrixXd V = eig.eigenvectors();
+  Eigen::MatrixXd D = eig.eigenvalues().asDiagonal();
+
+  // Clamp all eigenvalues to eps
+  for (int i = 0; i < H_orig.rows(); ++i)
+      D(i, i) = std::max(eps, D(i, i));
+  H_spd = V * D * V.inverse();
+}
+
+double SymmetricDirichletOneVertexElement::max_feasible_step(const VecV& _x, const VecV& _v, const VecC& _c)
+{
+  // get quadratic coefficients (ax^2 + b^x + c)
+  auto U11 = _x[0];
+  auto U12 = _x[1];
+  auto U21 = _c[0];
+  auto U22 = _c[1];
+  auto U31 = _c[2];
+  auto U32 = _c[3];
+
+  auto V11 = _v[0];
+  auto V12 = _v[1];
+  const VecV::Scalar V21 = 0.0;
+  const VecV::Scalar V22 = 0.0;
+  const VecV::Scalar V31 = 0.0;
+  const VecV::Scalar V32 = 0.0;
+
+  double a = V11*V22 - V12*V21 - V11*V32 + V12*V31 + V21*V32 - V22*V31;
+  double b = U11*V22 - U12*V21 - U21*V12 + U22*V11 - U11*V32 + U12*V31 + U31*V12 - U32*V11 + U21*V32 - U22*V31 - U31*V22 + U32*V21;
+  double c = U11*U22 - U12*U21 - U11*U32 + U12*U31 + U21*U32 - U22*U31;
+
+  double delta_in = pow(b,2) - 4*a*c;
+  if (delta_in < 0) {
+    return std::numeric_limits<double>::max();
+  }
+  double delta = sqrt(delta_in);
+  double t1 = (-b + delta)/ (2*a);
+  double t2 = (-b - delta)/ (2*a);
+
+  double tmp_n = std::min(t1,t2);
+  t1 = std::max(t1,t2); t2 = tmp_n;
+  // return the smallest negative root if it exists, otherwise return infinity
+  if (t1 > 0)
+  {
+    if (t2 > 0)
+    {
+      return 0.999 * t2;
+    }
+    else
+    {
+      return 0.999 * t1;
+    }
+  }
+  else
+  {
+    return std::numeric_limits<double>::max();
+  }
+}
+
+SymmetricDirichletOneRingProblem::SymmetricDirichletOneRingProblem()
+  :
+    FiniteElementProblem(2)
+{
+  FiniteElementProblem::add_set(&element_set);
+}
+
+void SymmetricDirichletOneRingProblem::add_triangle(const InputPositionVector2D& _current_positions, const SymmetricDirichletOneRingProblem::ReferencePositionVector2D& _reference_positions)
+{
+
+  SymmetricDirichletOneVertexElement::VecI indices;
+  indices << 0,1;
+  SymmetricDirichletOneVertexElement::VecC constants;
+  constants << _current_positions  (1, 0),  _current_positions (1, 1),
+               _current_positions  (2, 0),  _current_positions (2, 1),
+               _reference_positions(0, 0), _reference_positions(0, 1),
+               _reference_positions(1, 0), _reference_positions(1, 1),
+               _reference_positions(2, 0), _reference_positions(2, 1);
+  element_set.instances().add_element(indices, constants);
+}
+
+ SymmetricDirichletOneRingProblem::ReferencePositionVector2D SymmetricDirichletOneRingProblem::get_equilateral_refernce_positions(double _area)
+{
+  ReferencePositionVector2D equilateral_reference;
+  equilateral_reference << 0.0, 0.0,
+                           1.0, 0.0,
+                           0.5, 0.5*std::sqrt(3.0);
+  equilateral_reference *= _area/(0.5*0.5*std::sqrt(3.0));
+  return equilateral_reference;
+}
+
 } // namespace COMISO
 
 #endif //(COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)

NSolver/SymmetricDirichletProblem.hh

@@ -146,6 +146,9 @@ public:
   /// remove all constraints
   void clear_constraints() { fix_points.clear();}
 
+  /// get reference positions for an equilateral triangle with the given area
+  ReferencePositionVector2D get_equilateral_refernce_positions(double _area = 1.0);
+
 private:
   SymmetricDirichletElementSet element_set;
 
@@ -153,6 +156,72 @@ private:
 };
 
 
+class SymmetricDirichletOneVertexElement
+{
+public:
+
+  // define dimensions
+  const static int NV = 2; // the one u and v coordinate of the single vertex that is optimized. u1, v1
+  const static int NC = 10; // the two other positions of the triangle u2, v2, u3, and v3, and the three reference positions of the triangle, x1, y1, x2, y2, x3, y3
+
+  typedef Eigen::Matrix<size_t,NV,1> VecI;
+  typedef Eigen::Matrix<double,NV,1> VecV;
+  typedef Eigen::Matrix<double,NC,1> VecC;
+  typedef Eigen::Triplet<double> Triplet;
+
+  typedef Eigen::Matrix<double,12,1> Vector12;
+  typedef Eigen::Matrix<adouble,2,2> Matrix2x2ad;
+  typedef Eigen::Matrix<double,2,2> Matrix2x2d;
+  typedef Eigen::Matrix<double,6,6> Matrix6x6;
+
+  SymmetricDirichletOneVertexElement(){}
+
+  double eval_f       (const VecV& _x, const VecC& _c);
+  void   eval_gradient(const VecV& _x, const VecC& _c, VecV& _g);
+  void   eval_hessian (const VecV& _x, const VecC& _c, std::vector<Triplet>& _triplets);
+
+  void project_hessian(Eigen::MatrixXd& H_orig, Eigen::MatrixXd& H_spd, double eps);
+
+  double max_feasible_step(const VecV& _x, const VecV& _v, const VecC& /*_c*/);
+};
+
+
+/** \class SymmetricDirichletOneRingProblem
+
+    A problem that allows you to add triangles with reference positions for which
+    the symmetric dirichlet energy should be minimized. Secial case of SymmetricDirichletProblem
+    that optimizes only a single vertex for which the user inputs all adjacent triangles.
+*/
+
+class COMISODLLEXPORT SymmetricDirichletOneRingProblem : public FiniteElementProblem
+{
+public:
+
+  typedef FiniteElementSet<SymmetricDirichletOneVertexElement> SymmetricDirichletOneVertexElementSet;
+
+  typedef Eigen::Matrix<size_t,3,1> IndexVector;
+  typedef Eigen::Matrix<double,3,2> InputPositionVector2D;