Commit 3c9a64cf authored by Patric Schmitz's avatar Patric Schmitz

Merge remote-tracking branch 'origin/master' into cmake-overhaul

preparing merge into master
parents baa9ad84 16189854
Pipeline #12794 passed with stages
in 7 minutes and 43 seconds
Subproject commit 5c54ef0065f46b07dbea44644515a00665563235 Subproject commit 3623fbdfb4eba14a65926ff4034bc3916eab2cf0
...@@ -571,6 +571,9 @@ if (COMISO_BUILD_EXAMPLES ) ...@@ -571,6 +571,9 @@ if (COMISO_BUILD_EXAMPLES )
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_mosek_fusion_sdp/CMakeLists.txt" ) if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_mosek_fusion_sdp/CMakeLists.txt" )
add_subdirectory (Examples/small_mosek_fusion_sdp) add_subdirectory (Examples/small_mosek_fusion_sdp)
endif() endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_symmetric_dirichlet/CMakeLists.txt" )
add_subdirectory (Examples/small_symmetric_dirichlet)
endif()
endif (COMISO_BUILD_EXAMPLES ) endif (COMISO_BUILD_EXAMPLES )
......
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}
)
/*===========================================================================*\
* *
* 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)
...@@ -76,6 +76,8 @@ solve(NProblemInterface* _problem, ...@@ -76,6 +76,8 @@ solve(NProblemInterface* _problem,
// move to next constraint // move to next constraint
++n_constraints; ++n_constraints;
} else {
std::cerr << "Warning: COMISOSolver received a problem with non-equality constraints!!!" << std::endl;
} }
// resize matrix to final number of constraints // resize matrix to final number of constraints
......
...@@ -80,7 +80,7 @@ private: ...@@ -80,7 +80,7 @@ private:
} // namespace COMISO } // namespace COMISO
//============================================================================= //=============================================================================
// support std vectors // support std vectors
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::ConeConstraint); EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::ConeConstraint)
//============================================================================= //=============================================================================
#endif // ACG_CONECONSTRAINT_HH defined #endif // ACG_CONECONSTRAINT_HH defined
//============================================================================= //=============================================================================
......
...@@ -182,6 +182,7 @@ public: ...@@ -182,6 +182,7 @@ public:
for(unsigned int j=0; j<NV; ++j) for(unsigned int j=0; j<NV; ++j)
x_[j] = _x[instances_.index(i,j)]; x_[j] = _x[instances_.index(i,j)];
triplets_.clear();
element_.eval_hessian(x_, instances_.c(i), triplets_); element_.eval_hessian(x_, instances_.c(i), triplets_);
for(unsigned int j=0; j<triplets_.size(); ++j) for(unsigned int j=0; j<triplets_.size(); ++j)
......
...@@ -144,7 +144,7 @@ void IPOPTSolverLean::set_enable_all_lazy_contraints(const bool ...@@ -144,7 +144,7 @@ void IPOPTSolverLean::set_enable_all_lazy_contraints(const bool
static void throw_ipopt_solve_failure(Ipopt::ApplicationReturnStatus const status) static void throw_ipopt_solve_failure(Ipopt::ApplicationReturnStatus const status)
{ {
DEB_enter_func DEB_enter_func
DEB_warning(1, " IPOPT solve failure code is " << status) DEB_error(" IPOPT solve failure code is " << status)
// TODO: we could translate these return codes, but will not do it for now // TODO: we could translate these return codes, but will not do it for now
// enum ApplicationReturnStatus // enum ApplicationReturnStatus
// { // {
......
...@@ -88,7 +88,7 @@ private: ...@@ -88,7 +88,7 @@ private:
} // namespace COMISO } // namespace COMISO
//============================================================================= //=============================================================================
// support std vectors // support std vectors
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::LinearConstraint); EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(COMISO::LinearConstraint)
//============================================================================= //=============================================================================
#endif // ACG_LINEARCONSTRAINT_HH defined #endif // ACG_LINEARCONSTRAINT_HH defined
//============================================================================= //=============================================================================
......
...@@ -23,6 +23,7 @@ solve(NProblemGmmInterface* _problem) ...@@ -23,6 +23,7 @@ solve(NProblemGmmInterface* _problem)
{ {
DEB_enter_func; DEB_enter_func;
#if COMISO_SUITESPARSE_AVAILABLE #if COMISO_SUITESPARSE_AVAILABLE
converged_ = true;
// get problem size // get problem size
int n = _problem->n_unknowns(); int n = _problem->n_unknowns();
...@@ -98,16 +99,19 @@ solve(NProblemGmmInterface* _problem) ...@@ -98,16 +99,19 @@ solve(NProblemGmmInterface* _problem)
_problem->store_result(P(x)); _problem->store_result(P(x));
DEB_line(2, "Newton solver reached max regularization but did not " DEB_line(2, "Newton solver reached max regularization but did not "
"converge"); "converge");
converged_ = false;
return false; return false;
} }
} }
} }
_problem->store_result(P(x)); _problem->store_result(P(x));
DEB_line(2, "Newton Solver did not converge!!! after iterations."); DEB_line(2, "Newton Solver did not converge!!! after iterations.");
converged_ = false;
return false; return false;
#else #else
DEB_warning(1,"NewtonSolver requires not-available CholmodSolver"); DEB_warning(1,"NewtonSolver requires not-available CholmodSolver");
converged_ = false;
return false; return false;
#endif #endif
} }
...@@ -120,6 +124,7 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A, ...@@ -120,6 +124,7 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
const VectorD& _b) const VectorD& _b)
{ {
DEB_time_func_def; DEB_time_func_def;
converged_ = false;
const double KKT_res_eps = 1e-6; const double KKT_res_eps = 1e-6;
const int max_KKT_regularization_iters = 40; const int max_KKT_regularization_iters = 40;
...@@ -160,10 +165,11 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A, ...@@ -160,10 +165,11 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
double kkt_res2(0.0); double kkt_res2(0.0);
double constraint_res2(0.0); double constraint_res2(0.0);
int reg_iters(0); int reg_iters(0);
bool fact_ok = true;
do do
{ {
// get Newton search direction by solving LSE // 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; first_factorization = false;
if(fact_ok) if(fact_ok)
...@@ -204,7 +210,7 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A, ...@@ -204,7 +210,7 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
} }
++reg_iters; ++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? // no valid step could be found?
if(kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2 || reg_iters >= max_KKT_regularization_iters) if(kkt_res2 > KKT_res_eps || constraint_res2 > max_allowed_constraint_violation2 || reg_iters >= max_KKT_regularization_iters)
...@@ -250,8 +256,11 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A, ...@@ -250,8 +256,11 @@ int NewtonSolver::solve(NProblemInterface* _problem, const SMatrixD& _A,
<< ", KKT residual^2 = " << kkt_res2); << ", KKT residual^2 = " << kkt_res2);
// converged? // converged?
if(newton_decrement < eps_ || std::abs(t) < eps_ls_) if(newton_decrement < eps_ || std::abs(t) < eps_ls_ )
{
converged_ = true;
break; break;
}
++iter; ++iter;
} }
...@@ -397,11 +406,11 @@ bool NewtonSolver::numerical_factorization(SMatrixD& _KKT) ...@@ -397,11 +406,11 @@ bool NewtonSolver::numerical_factorization(SMatrixD& _KKT)
DEB_enter_func; DEB_enter_func;
switch(solver_type_) switch(solver_type_)
{ {
case LS_EigenLU: case LS_EigenLU:
lu_solver_.factorize(_KKT); lu_solver_.factorize(_KKT);
return (lu_solver_.info() == Eigen::Success); return (lu_solver_.info() == Eigen::Success);
#if COMISO_SUITESPARSE_AVAILABLE #if COMISO_SUITESPARSE_AVAILABLE
case LS_Umfpack: case LS_Umfpack:
umfpack_solver_.factorize(_KKT); umfpack_solver_.factorize(_KKT);
return (umfpack_solver_.info() == Eigen::Success); return (umfpack_solver_.info() == Eigen::Success);
#endif #endif
......
...@@ -84,6 +84,8 @@ public: ...@@ -84,6 +84,8 @@ public:
solver_type_ = _st; solver_type_ = _st;
} }
bool converged() { return converged_; }
protected: protected:
bool factorize(NProblemInterface* _problem, const SMatrixD& _A, bool factorize(NProblemInterface* _problem, const SMatrixD& _A,
...@@ -149,6 +151,8 @@ private: ...@@ -149,6 +151,8 @@ private:
// deprecated // deprecated
bool constant_hessian_structure_; bool constant_hessian_structure_;
bool converged_;
}; };
......
...@@ -84,6 +84,10 @@ public: ...@@ -84,6 +84,10 @@ public:
x_ = x; x_ = x;
} }
const SMatrixNP &A() const {return A_;}
const Eigen::VectorXd &b() const {return b_;}
double c() const {return c_;}
// get current solution // get current solution
Eigen::VectorXd& x() { return x_;} Eigen::VectorXd& x() { return x_;}
......
/*===========================================================================*\
* *
* 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