Merge branch 'symmetric_dirichlet_problem' into 'master'

Symmetric dirichlet problem

See merge request !43

CMakeLists.txt

@@ -571,6 +571,9 @@ if (COMISO_BUILD_EXAMPLES )
   if( EXISTS "${CMAKE_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" )
+  add_subdirectory (Examples/small_symmetric_dirichlet)
+  endif()
   
 endif (COMISO_BUILD_EXAMPLES )
 

Examples/small_symmetric_dirichlet/CMakeLists.txt

+include (CoMISoExample)
+
+acg_add_executable (small_symmetric_dirichlet ${sources} ${headers} )
+
+# enable rpath linking
+set_target_properties(small_symmetric_dirichlet PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
+
+target_link_libraries (small_symmetric_dirichlet
+  CoMISo
+  ${COMISO_LINK_LIBRARIES}
+)

Examples/small_symmetric_dirichlet/main.cc

+/*===========================================================================*\
+ *                                                                           *
+ *                               CoMISo                                      *
+ *      Copyright (C) 2008-2019 by Computer Graphics Group, RWTH Aachen      *
+ *                           www.rwth-graphics.de                            *
+ *                                                                           *
+ *---------------------------------------------------------------------------*
+ *  This file is part of CoMISo.                                             *
+ *                                                                           *
+ *  CoMISo is free software: you can redistribute it and/or modify           *
+ *  it under the terms of the GNU General Public License as published by     *
+ *  the Free Software Foundation, either version 3 of the License, or        *
+ *  (at your option) any later version.                                      *
+ *                                                                           *
+ *  CoMISo is distributed in the hope that it will be useful,                *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of           *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            *
+ *  GNU General Public License for more details.                             *
+ *                                                                           *
+ *  You should have received a copy of the GNU General Public License        *
+ *  along with CoMISo.  If not, see <http://www.gnu.org/licenses/>.          *
+ *                                                                           *
+\*===========================================================================*/
+
+#include <CoMISo/Config/config.hh>
+#include <iostream>
+
+#if (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
+
+#include <CoMISo/Utils/StopWatch.hh>
+#include <CoMISo/NSolver/NewtonSolver.hh>
+#include <CoMISo/NSolver/SymmetricDirichletProblem.hh>
+#include <vector>
+
+//------------------------------------------------------------------------------------------------------
+
+// Example main
+int main(void)
+{
+  std::cout << "---------- 1) Get an instance of a SymmetricDirichletProblem..." << std::endl;
+
+  // then create finite element problem and add sets
+  unsigned int n_vertices = 4;
+  COMISO::SymmetricDirichletProblem sd_problem(n_vertices);
+  COMISO::SymmetricDirichletProblem::IndexVector indices(0,1,2);
+  COMISO::SymmetricDirichletProblem::ReferencePositionVector2D positions;
+  positions << 0, 0, // first point
+               1, 0, // second point
+               0, 1; // third point
+  sd_problem.add_triangle(indices, positions);
+  COMISO::SymmetricDirichletProblem::IndexVector indices2(3,2,1);
+  sd_problem.add_triangle(indices2, positions); // same reference positions can be used because we want both triangles to be isosceles
+
+  std::vector<double> initial_solution{0.0,0.0,
+                                       2.0,0.0,
+                                       2.0,2.0,
+                                       3.0,4.0};
+  sd_problem.x() = initial_solution;
+
+  std::cout << "---------- 2) Set up constraints..." << std::endl;
+  // fix first vertex to origin to fix translational degree of freedom
+  sd_problem.add_fix_point_constraint(0, 0.0, 0.0);
+  // fix v coordinate of second vertex to 0 to fix rotational degree of freedom
+  sd_problem.add_fix_coordinate_constraint(1, 1, 0.0);
+
+  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
+  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;
+
+  return 0;
+}
+
+#else // (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
+
+int main(void)
+{
+  std::cerr << "Warning: Example cannot be executed since either EIGEN3 or ADOLC is not available..." << std::endl;
+  return 0;
+}
+
+
+#endif // (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)

NSolver/FiniteElementProblem.hh

@@ -182,6 +182,7 @@ public:
       for(unsigned int j=0; j<NV; ++j)
         x_[j] = _x[instances_.index(i,j)];
 
+      triplets_.clear();
       element_.eval_hessian(x_, instances_.c(i), triplets_);
 
       for(unsigned int j=0; j<triplets_.size(); ++j)

NSolver/NewtonSolver.cc

@@ -160,10 +160,11 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
     double kkt_res2(0.0);
     double constraint_res2(0.0);
     int    reg_iters(0);
+    bool fact_ok = true;
     do
     {
       // get Newton search direction by solving LSE
-      bool fact_ok = factorize(_problem, _A, _b, x, regularize_hessian, regularize_constraints, first_factorization);
+      fact_ok = factorize(_problem, _A, _b, x, regularize_hessian, regularize_constraints, first_factorization);
       first_factorization = false;
 
       if(fact_ok)
@@ -204,7 +205,7 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
       }
       ++reg_iters;
     }
-    while( (kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2) && reg_iters < max_KKT_regularization_iters);
+    while( (!fact_ok || kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2) && reg_iters < max_KKT_regularization_iters);
 
     // no valid step could be found?
     if(kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2 || reg_iters >= max_KKT_regularization_iters)

NSolver/SymmetricDirichletProblem.cc

+/*===========================================================================*\
+ *                                                                           *
+ *                               CoMISo                                      *
+ *      Copyright (C) 2008-2019 by Computer Graphics Group, RWTH Aachen      *
+ *                           www.rwth-graphics.de                            *
+ *                                                                           *
+ *---------------------------------------------------------------------------*
+ *  This file is part of CoMISo.                                             *
+ *                                                                           *
+ *  CoMISo is free software: you can redistribute it and/or modify           *
+ *  it under the terms of the GNU General Public License as published by     *
+ *  the Free Software Foundation, either version 3 of the License, or        *
+ *  (at your option) any later version.                                      *
+ *                                                                           *
+ *  CoMISo is distributed in the hope that it will be useful,                *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of           *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            *
+ *  GNU General Public License for more details.                             *
+ *                                                                           *
+ *  You should have received a copy of the GNU General Public License        *
+ *  along with CoMISo.  If not, see <http://www.gnu.org/licenses/>.          *
+ *                                                                           *
+\*===========================================================================*/
+
+#include <CoMISo/Config/config.hh>
+
+#if (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
+
+#include "SymmetricDirichletProblem.hh"
+
+
+namespace COMISO { 
+
+double SymmetricDirichletElement::eval_f(const VecV& _x, const VecC& _c)
+{
+  double y = 0.0;
+
+  Vector12 x;
+  x << _x[0], _x[1], _x[2], _x[3], _x[4], _x[5],
+       _c[0], _c[1], _c[2], _c[3], _c[4], _c[5];
+
+  if(tape_available_)
+  {
+    int ec = function(tape_, 1, 12, const_cast<double*>(x.data()), &y);
+
+#ifdef ADOLC_RET_CODES
+    std::cout << "Info: function() returned code " << ec << std::endl;
+#endif
+
+    // tape not valid anymore? retape and evaluate again
+    if(ec < 0)
+    {
+      retape();
+      function(tape_, 1, 12, const_cast<double*>(x.data()), &y);
+    }
+  }
+  else
+  {
+    retape();
+    function(tape_, 1, 12, const_cast<double*>(x.data()), &y);
+  }
+
+  return y;
+}
+
+void SymmetricDirichletElement::eval_gradient(const VecV& _x, const VecC& _c, VecV& _g)
+{
+  Vector12 x;
+  x << _x[0], _x[1], _x[2], _x[3], _x[4], _x[5],
+       _c[0], _c[1], _c[2], _c[3], _c[4], _c[5];
+
+  Vector12 grad;
+  grad.setZero();
+
+  // evaluate gradient
+  int ec = gradient(tape_, 12, x.data(), grad.data());
+
+  // check if retaping is required
+  if(ec < 0)
+  {
+#ifdef ADOLC_RET_CODES
+  std::cout << __FUNCTION__ << " invokes retaping of function due to discontinuity! Return code: " << ec << std::endl;
+#endif
+    retape();
+    gradient(tape_, 12, x.data(), grad.data());
+  }
+
+  _g = grad.block(0,0,6,1);
+}
+
+void SymmetricDirichletElement::eval_hessian(const VecV& _x, const VecC& _c, std::vector<Triplet>& _triplets)
+{
+  Vector12 x;
+  x << _x[0], _x[1], _x[2], _x[3], _x[4], _x[5],
+       _c[0], _c[1], _c[2], _c[3], _c[4], _c[5];
+
+  // dense hessian
+  {
+    auto dense_hessian = new double*[12];
+    for(int i = 0; i < 12; ++i)
+      dense_hessian[i] = new double[i+1];
+
+    int ec = hessian(tape_, 12, const_cast<double*>(x.data()), dense_hessian);
+
+    if(ec < 0) {
+#ifdef ADOLC_RET_CODES
+      std::cout << __FUNCTION__ << " invokes retaping of function due to discontinuity! Return code: " << ec << std::endl;
+#endif
+      // Retape function if return code indicates discontinuity
+      retape();
+      ec = hessian(tape_, 12, const_cast<double*>(x.data()), dense_hessian);
+    }
+
+#ifdef ADOLC_RET_CODES
+    std::cout << "Info: hessian() returned code " << ec << std::endl;
+#endif
+
+    Eigen::MatrixXd H(6,6);
+    for (int i = 0; i < 6; ++i)
+    {
+      H(i,i) = dense_hessian[i][i];
+      for (int j = 0; j < i; ++j)
+      {
+        H(i,j) = dense_hessian[i][j];
+        H(j,i) = dense_hessian[i][j];
+      }
+    }
+
+    Eigen::MatrixXd Hspd(6,6);
+    project_hessian(H, Hspd, 1e-6);
+//    Hspd = H;
+
+    _triplets.reserve(6*6);
+    for (int i = 0; i < 6; ++i)
+      for (int j = 0; j < 6; ++j)
+        _triplets.push_back(Triplet(i,j,Hspd(i,j)));
+
+    for (int i = 0; i < 6; ++i)
+      delete dense_hessian[i];
+    delete[] dense_hessian;
+  }
+}
+
+void SymmetricDirichletElement::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 SymmetricDirichletElement::max_feasible_step(const VecV& _x, const VecV& _v, const VecC&)
+{
+  // get quadratic coefficients (ax^2 + b^x + c)
+  auto U11 = _x[0];
+  auto U12 = _x[1];
+  auto U21 = _x[2];
+  auto U22 = _x[3];
+  auto U31 = _x[4];
+  auto U32 = _x[5];
+
+  auto V11 = _v[0];
+  auto V12 = _v[1];
+  auto V21 = _v[2];
+  auto V22 = _v[3];
+  auto V31 = _v[4];
+  auto V32 = _v[5];
+
+  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();
+  }
+}
+
+adouble SymmetricDirichletElement::f_adouble(const adouble* _x)
+{
+  Matrix2x2ad B;
+  B(0,0) = _x[2]-_x[0];
+  B(0,1) = _x[4]-_x[0];
+  B(1,0) = _x[3]-_x[1];
+  B(1,1) = _x[5]-_x[1];
+  Matrix2x2ad Bin = B.inverse();
+
+  Matrix2x2ad R;
+  R(0,0) = _x[6+2]-_x[6+0];
+  R(0,1) = _x[6+4]-_x[6+0];
+  R(1,0) = _x[6+3]-_x[6+1];
+  R(1,1) = _x[6+5]-_x[6+1];
+  Matrix2x2ad Rin = R.inverse();
+
+  adouble area = 0.5 * R.determinant();
+  if (B.determinant() * area <= 0)
+  {
+    adouble res = std::numeric_limits<double>::max();
+    return res;
+  }
+
+  Matrix2x2ad J   =  B * Rin;
+  Matrix2x2ad Jin =  R * Bin;
+
+  adouble res = J.squaredNorm() + Jin.squaredNorm();
+
+  return area * (res - 4);
+}
+
+void SymmetricDirichletElement::retape()
+{
+  std::cerr << "re-tape..." << std::endl;
+  adouble y_d = 0.0;
+  double y = 0.0;
+
+
+  trace_on(tape_); // Start taping
+  adouble* xa = new adouble[12];
+
+  // Fill data vectors
+  xa[0] <<= 0.0;
+  xa[1] <<= 0.0;
+  xa[2] <<= 1.0;
+  xa[3] <<= 0.0;
+  xa[4] <<= 0.0;
+  xa[5] <<= 1.0;
+
+  xa[6+0] <<= 0.0;
+  xa[6+1] <<= 0.0;
+  xa[6+2] <<= 1.0;
+  xa[6+3] <<= 0.0;
+  xa[6+4] <<= 0.0;
+  xa[6+5] <<= 1.0;
+
+  y_d = f_adouble(xa);
+
+  y_d >>= y;
+
+  trace_off();
+
+#ifdef ADOLC_STATS
+  print_stats();
+#endif
+
+  delete[] xa;
+}
+
+void SymmetricDirichletElement::init_tape()
+{
+  static bool tape_initialized = false;
+  static short int tape;
+
+  if (!tape_initialized)
+  {
+    tape = static_cast<short int>(TapeIDSingleton::Instance()->requestId());
+    tape_ = tape;
+    retape();
+  }
+
+  tape_initialized = true;
+  tape_available_ = true;
+  tape_ = tape;
+}
+
+/// Default constructor
+SymmetricDirichletProblem::SymmetricDirichletProblem(const unsigned int _n_vertices)
+  :
+    FiniteElementProblem(2*_n_vertices)
+{
+  FiniteElementProblem::add_set(&element_set);
+}
+
+void SymmetricDirichletProblem::add_triangle(const IndexVector& _vertex_indices, const ReferencePositionVector2D& _reference_positions)
+{
+  SymmetricDirichletElement::VecI indices;
+  indices << 2*_vertex_indices[0], 2*_vertex_indices[0]+1,
+             2*_vertex_indices[1], 2*_vertex_indices[1]+1,
+             2*_vertex_indices[2], 2*_vertex_indices[2]+1;
+  SymmetricDirichletElement::VecC constants;
+  constants << _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);
+}
+
+void SymmetricDirichletProblem::add_fix_point_constraint(int _vertex_index, double _fix_u, double _fix_v)
+{
+  add_fix_coordinate_constraint(_vertex_index, 0, _fix_u);
+  add_fix_coordinate_constraint(_vertex_index, 1, _fix_v);
+}
+
+void SymmetricDirichletProblem::add_fix_coordinate_constraint(int _vertex_index, int _coordinate, double _fix_coordinate)
+{
+  fix_points.push_back(std::make_pair(2*_vertex_index + _coordinate, _fix_coordinate));
+}
+
+void SymmetricDirichletProblem::get_constraints(SMatrixD& _A, VectorD& _b)
+{
+  _A.resize(fix_points.size(), n_unknowns());
+  _b.resize(fix_points.size());
+  std::vector<Eigen::Triplet<double>> triplets;
+  for (size_t i = 0; i < fix_points.size(); ++i)
+  {
+    _b[i] = fix_points[i].second;
+    triplets.push_back(Eigen::Triplet<double>(i, fix_points[i].first, 1));
+  }
+  _A.setFromTriplets(triplets.begin(), triplets.end());
+}
+
+} // namespace COMISO
+
+#endif //(COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
+

NSolver/SymmetricDirichletProblem.hh

+/*===========================================================================*\
+ *                                                                           *
+ *                               CoMISo                                      *
+ *      Copyright (C) 2008-2019 by Computer Graphics Group, RWTH Aachen      *
+ *                           www.rwth-graphics.de                            *
+ *                                                                           *
+ *---------------------------------------------------------------------------*
+ *  This file is part of CoMISo.                                             *
+ *                                                                           *
+ *  CoMISo is free software: you can redistribute it and/or modify           *
+ *  it under the terms of the GNU General Public License as published by     *
+ *  the Free Software Foundation, either version 3 of the License, or        *
+ *  (at your option) any later version.                                      *
+ *                                                                           *
+ *  CoMISo is distributed in the hope that it will be useful,                *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of           *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            *
+ *  GNU General Public License for more details.                             *
+ *                                                                           *
+ *  You should have received a copy of the GNU General Public License        *
+ *  along with CoMISo.  If not, see <http://www.gnu.org/licenses/>.          *
+ *                                                                           *
+\*===========================================================================*/
+//=============================================================================
+//
+//  CLASS SymmetricDirichletProblem
+//
+//=============================================================================
+
+
+#ifndef COMISO_SYMMETRICDIRICHLETTPROBLEM_HH
+#define COMISO_SYMMETRICDIRICHLETTPROBLEM_HH
+
+
+//== INCLUDES =================================================================
+
+#include <CoMISo/Config/config.hh>
+
+#if (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
+
+#include <CoMISo/Config/CoMISoDefines.hh>
+#include "FiniteElementProblem.hh"
+
+#include <adolc/adolc.h>
+#include <adolc/adouble.h>
+#include <adolc/drivers/drivers.h>
+#include <adolc/sparse/sparsedrivers.h>
+#include <adolc/taping.h>
+#include <adolc/drivers/taylor.h>
+#include <CoMISo/NSolver/TapeIDSingleton.hh>
+
+
+//== FORWARDDECLARATIONS ======================================================
+
+//== NAMESPACES ===============================================================
+
+namespace COMISO { 
+
+//== CLASS DEFINITION =========================================================
+
+class SymmetricDirichletElement
+{
+public:
+
+  // define dimensions
+  const static int NV = 6; // the three u and v coordinates. u1, v1, u2, v2, u3, v3
+  const static int NC = 6; // 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,6,6> Matrix6x6;
+
+  SymmetricDirichletElement()
+    :
+      tape_(-1),
+      tape_available_(false)
+  {
+    init_tape();
+  }
+
+  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*/);
+
+  adouble f_adouble( const adouble* _x);
+
+  void retape();
+  void init_tape();
+
+private:
+
+  // index of tape
+  short int tape_;
+
+  // taping already done?
+  bool tape_available_;
+};
+
+
+/** \class SymmetricDirichletProblem
+
+    A problem that allows you to add triangles with reference positions for which
+    the symmetric dirichlet energy should be minimized
+*/
+
+class COMISODLLEXPORT SymmetricDirichletProblem : public FiniteElementProblem
+{
+public:
+
+  typedef FiniteElementSet<SymmetricDirichletElement> SymmetricDirichletElementSet;
+
+  typedef Eigen::Matrix<size_t,3,1> IndexVector;
+  typedef Eigen::Matrix<double,3,2> ReferencePositionVector2D;
+  typedef Eigen::Matrix<double,3,3> ReferencePositionVector3D;
+
+  typedef Eigen::VectorXd             VectorD;
+  typedef Eigen::SparseMatrix<double> SMatrixD;
+
+
+  /// Default constructor
+  SymmetricDirichletProblem(const unsigned int _n_vertices);
+
+  void add_triangle(const IndexVector& _vertex_indices, const ReferencePositionVector2D& _reference_positions);
+
+  /// add fix point constraint. Note that the final constraints still need to be passed
+  /// to the solver by you after retrieving them via get_constraints
+  void add_fix_point_constraint(int _vertex_index, double _fix_u, double _fix_v);
+
+  /// add fix point constraint for only one of the two coordinates of a vertex. Note that the
+  /// final constraints still need to be passed to the solver by you after retrieving them via get_constraints
+