Commit ac3920ef authored by Max Lyon's avatar Max Lyon
Browse files

Merge branch 'master' into exact_constraint_satisfaction

parents bf86d830 37557048
......@@ -13,19 +13,6 @@ macos-c++11:
tags:
- Apple
CoMISo-VS2013-Qt-5.5.1-x64:
variables:
BUILD_PLATFORM: "VS2013"
ARCHITECTURE: "x64"
QT_VERSION: "Qt5.5.1"
GIT_SUBMODULE_STRATEGY: recursive
COMPILER: "VS2013"
script: "CI\\Windows.bat"
tags:
- VS2013
- Qt551
CoMISo-VS2017-Qt-5.10.1-x64:
variables:
BUILD_PLATFORM: "VS2017"
......
Subproject commit 3623fbdfb4eba14a65926ff4034bc3916eab2cf0
Subproject commit 5e38f629fc7c7ae12316db1f19e6f878c1387bb2
......@@ -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()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_exact_constraint_satifaction_example/CMakeLists.txt" )
......
......@@ -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;
}
......
......@@ -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)
{
......
......@@ -3,15 +3,18 @@
#include <CoMISo/Config/config.hh>
#include <CoMISo/Utils/StopWatch.hh>
#include <vector>
#include <CoMISo/NSolver/NProblemInterface.hh>
#include <Base/Code/Quality.hh>
#include <Base/Debug/DebOut.hh>
LOW_CODE_QUALITY_SECTION_BEGIN
#include <Eigen/Eigen>
#include <Eigen/Core>
#include <Eigen/Sparse>
LOW_CODE_QUALITY_SECTION_END
#include <vector>
//== NAMESPACES ===============================================================
......@@ -26,36 +29,31 @@ public:
// typedef Eigen::DynamicSparseMatrix<double,Eigen::ColMajor> SMatrixNP;
QuadraticProblem()
: A_(0,0), b_(Eigen::VectorXd::Index(0)), c_(0.0)
{
}
: A_(0, 0), b_(Eigen::VectorXd::Index(0)), c_(0.0)
{}
QuadraticProblem(SMatrixNP& _A, Eigen::VectorXd& _b, const double _c)
: A_(_A), c_(_c)
QuadraticProblem(SMatrixNP& _A, const Eigen::VectorXd& _b, const double _c)
{
if(A_.rows() != A_.cols())
std::cerr << "Warning: matrix not square in QuadraticProblem" << std::endl;
b_ = _b;
x_ = Eigen::VectorXd::Zero(A_.cols());
set_A(_A);
set_b(_b);
set_c(_c);
}
// number of unknowns
virtual int n_unknowns()
{
return A_.rows();
return A_.rows();
}
// initial value where the optimization should start from
virtual void initial_x(double* _x)
{
for( int i=0; i<this->n_unknowns(); ++i)
_x[i] = x_[i];
for (int i = 0; i < this->n_unknowns(); ++i)
_x[i] = x_[i];
}
// function evaluation at location _x
virtual double eval_f( const double* _x )
virtual double eval_f(const double* _x)
{
Eigen::Map<const Eigen::VectorXd> x(_x, this->n_unknowns());
......@@ -63,22 +61,22 @@ public:
}
// gradient evaluation at location _x
virtual void eval_gradient( const double* _x, double* _g)
virtual void eval_gradient(const double* _x, double* _g)
{
Eigen::Map<const Eigen::VectorXd> x(_x, this->n_unknowns());
Eigen::Map<Eigen::VectorXd> g(_g, this->n_unknowns());
g = A_*x - b_;
}
}
// hessian matrix evaluation at location _x
virtual void eval_hessian ( const double* _x, SMatrixNP& _H)
virtual void eval_hessian(const double* _x, SMatrixNP& _H)
{
_H = A_;
}
// print result
virtual void store_result ( const double* _x )
virtual void store_result(const double* _x)
{
Eigen::Map<const Eigen::VectorXd> x(_x, this->n_unknowns());
x_ = x;
......@@ -89,16 +87,15 @@ public:
double c() const {return c_;}
// get current solution
Eigen::VectorXd& x() { return x_;}
Eigen::VectorXd& x() { return x_; }
// advanced properties
virtual bool constant_hessian() const { return true; }
virtual bool constant_hessian() const { return true; }
void set_A(const SMatrixNP& _A)
{
DEB_error_if(_A.rows() != _A.cols(), "Square matrix expected");
A_ = _A;
if(A_.rows() != A_.cols())
std::cerr << "Warning: matrix not square in QuadraticProblem" << std::endl;
x_ = Eigen::VectorXd::Zero(A_.cols());
}
......@@ -107,19 +104,18 @@ public:
b_ = _b;
}
void set_c( const double _c)
void set_c(const double _c)
{
c_ = _c;
}
private:
// quadratic problem 0.5*x^T A x -x^t b + c
SMatrixNP A_;
Eigen::VectorXd b_;
double c_;
// current solution, which is also used as initial value
Eigen::VectorXd x_;
SMatrixNP A_;
Eigen::VectorXd b_;
double c_;
// current solution, which is also used as initial value
Eigen::VectorXd x_;
};
//=============================================================================
......
......@@ -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)