Commit 6ee346d3 authored by Pit Nestle's avatar Pit Nestle

SLIM Solver softconstraints added

parent 112e49c3
Pipeline #7480 failed with stages
in 8 minutes and 17 seconds
......@@ -70,11 +70,11 @@ void SLIMProblemInterface::fixVertex(int i)
this->softconstraintIndices.push_back(i);
}
void SLIMProblemInterface::addRelation(int a, int b, Eigen::Matrix2d transformation, Eigen::Vector2d translation)
void SLIMProblemInterface::addRelation(int a, int b, int tr, double tu, double tv)
{
this->relationIndices.push_back(Eigen::Vector2i(a, b));
this->transformation.push_back(transformation);
this->translation.push_back(translation);
this->transformation.push_back(tr);
this->translation.push_back({tu, tv});
}
void SLIMProblemInterface::storeResult(Eigen::MatrixXd &res)
......
......@@ -69,7 +69,7 @@ public:
void getFace(size_t i, int& _a, int& _b, int& _c);
void fixVertex(int i);
void addRelation(int a, int b, Eigen::Matrix2d transformation, Eigen::Vector2d translation);
void addRelation(int a, int b, int tr, double tu, double tv);
virtual double computeEnergy(Eigen::MatrixXd& _j, Eigen::VectorXd& _a) = 0; //Evaluates the energy. _j #Faces x 4
virtual double computeEnergySV(Eigen::MatrixXd& _sv, Eigen::VectorXd& _a) = 0;
......@@ -125,7 +125,7 @@ public:
Eigen::VectorXd area;
//Relation between vertices
std::vector<Eigen::Matrix2d> transformation;
std::vector<int> transformation;
std::vector<Eigen::Vector2d> translation;
std::vector<Eigen::Vector2i> relationIndices;
......
......@@ -40,6 +40,8 @@ int SLIMSolver::solve(SLIMProblemInterface &problem, int iterations)
std::cout << "solve():solveLinearSystem()" << std::endl;
addProximalTerm(problem, ata, atb);
addSoftConstraint(problem, ata, atb);
addRelationConstraint(problem, ata, atb);
solveLinearSystem(ata, atb, p);
std::cout << "solve():calculateSearchDirection()" << std::endl;
......@@ -258,7 +260,7 @@ void SLIMSolver::addProximalTerm(SLIMProblemInterface &p, Eigen::SparseMatrix<do
void SLIMSolver::addSoftConstraint(SLIMProblemInterface &p, Eigen::SparseMatrix<double> &ata, Eigen::VectorXd &atb)
{
double penalty = 10e4;
double penalty = 10e16;
int v = p.n_vertices();
int n = p.softconstraintIndices.size();
......@@ -283,7 +285,7 @@ void SLIMSolver::addSoftConstraint(SLIMProblemInterface &p, Eigen::SparseMatrix<
void SLIMSolver::addRelationConstraint(SLIMProblemInterface &p, Eigen::SparseMatrix<double> &ata, Eigen::VectorXd &atb)
{
double penalty = 10e4;
double penalty = 10e16;
int v = p.n_vertices();
int n = p.relationIndices.size();
......@@ -294,6 +296,15 @@ void SLIMSolver::addRelationConstraint(SLIMProblemInterface &p, Eigen::SparseMat
Eigen::VectorXd t;
t.resize(2*n);
std::vector<Eigen::Matrix2d> rotation(4);
rotation[0] << 1, 0, 0, 1;
Eigen::Matrix2d rot90;
rot90 << 0, -1, 1, 0;
for (int i=1; i < 4; i++)
rotation[i] = rot90*rotation[i-1];
std::vector<Eigen::Triplet<double>> trip_sc;
for (int i = 0; i < n; i++)
......@@ -301,27 +312,27 @@ void SLIMSolver::addRelationConstraint(SLIMProblemInterface &p, Eigen::SparseMat
int a = p.relationIndices[i](0);
int b = p.relationIndices[i](1);
double m00 = p.transformation[i](0, 0);
double m01 = p.transformation[i](0, 1);
double m10 = p.transformation[i](1, 0);
double m11 = p.transformation[i](1, 1);
double m00 = rotation[p.transformation[i]](0, 0);
double m01 = rotation[p.transformation[i]](0, 1);
double m10 = rotation[p.transformation[i]](1, 0);
double m11 = rotation[p.transformation[i]](1, 1);
//x_a*m00 + y_a*m01 - x_b = t_x
//x_a*m00 + y_a*m01 - x_b = -t_x
trip_sc.push_back(Eigen::Triplet<double>(i, a+0, m00)); //x_a*m00
trip_sc.push_back(Eigen::Triplet<double>(i, a+v, m01)); //y_a*m01
trip_sc.push_back(Eigen::Triplet<double>(i, b+0, -1)); //x_b* -1
t(i) = p.translation[i](0);
t(i) = -p.translation[i](0);
//x_a*m10 + y_a*m11 - y_b = t_y
//x_a*m10 + y_a*m11 - y_b = -t_y
trip_sc.push_back(Eigen::Triplet<double>(i+n, a+0, m10)); //x_a*m10
trip_sc.push_back(Eigen::Triplet<double>(i+n, a+v, m11)); //y_a*m11
trip_sc.push_back(Eigen::Triplet<double>(i+n, b+v, -1)); //y_b* -1
t(i+n) = p.translation[i](1);
t(i+n) = -p.translation[i](1);
}
sc.setFromTriplets(trip_sc.begin(), trip_sc.end());
atb += sc.transpose() * t;
ata += sc.transpose() * sc;
atb += (sc.transpose() * t) * penalty;
ata += (sc.transpose() * sc) * penalty;
}
......
......@@ -51,7 +51,7 @@ class COMISODLLEXPORT SLIMSolver
public:
/// Default constructor
SLIMSolver(const double _lambda = 1e-4, const int _max_iters = 100)
SLIMSolver(const double _lambda = 1e-4, const int _max_iters = 20)
: iteration(0),
lambda_(_lambda),
max_iters_(_max_iters)
......
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#include <iostream>
#include <cmath>
#include <Eigen/Eigen>
#include <CoMISo/Utils/StopWatch.hh>
#include <vector>
#include "CoMISo/NSolver/SLIMProblemInterface.hh"
#include "CoMISo/NSolver/SLIMSolver.hh"
namespace COMISO {
class SLIM_Symmetric_Dirichlet : public COMISO::SLIMProblemInterface
{
public:
SLIM_Symmetric_Dirichlet() {}
~SLIM_Symmetric_Dirichlet() {}
// SLIMProblemInterface interface
public:
double computeEnergy(Eigen::MatrixXd& _j, Eigen::VectorXd& _a)
{
double res = 0;
for (int i = 0; i < _j.rows(); i++)
{
Eigen::Matrix2d j;
j << _j(i, 0), _j(i, 1),
_j(i, 2), _j(i, 3);
double det = _j(i, 0) * _j(i, 3) - _j(i, 1) *_j(i, 2);
if (det <= 0)
return std::numeric_limits<double>::infinity();
double j_n = j.norm();
double ji_n= j.inverse().norm();
//res += _a(i) * (j_n + ji_n);
res += _a(i) * ((std::pow(j_n, 2) + std::pow(ji_n, 2) - 4));
}
//std::printf("SLIM_Symmetric_Dirichlet::computeEnergy(...) = %f\n", res);
return res;
}
double computeEnergySV(Eigen::MatrixXd& _sv, Eigen::VectorXd& _a)
{
double res = 0;
for (int i = 0; i < _sv.rows(); i++)
{
double s0 = _sv(i, 0);
double s1 = _sv(i, 1);
res += _a(i) * (std::pow(s0, 2) +
std::pow(s0, -2) +
std::pow(s1, 2) +
std::pow(s1, -2) - 4);
}
return res;
}
void computeEnergyGradientSV(Eigen::MatrixXd& _sv, Eigen::MatrixXd& _g) override
{
for (int i = 0; i < _sv.cols(); i++)
{
double s1 = _sv(0, i);
double s2 = _sv(1, i);
_g.col(i) << 2 * (s1 - std::pow(s1, -3)),
2 * (s2 - std::pow(s2, -3));
}
}
double computeEnergyGradientSV_i(double _sv)
{
return (_sv - std::pow(_sv, -3));
}
double s_lambda_i(double sigma_i)
{
return 1;
}
void s_lambda(double sigma_1, double sigma_2, Eigen::Matrix2d &s_lambda)
{
s_lambda << 1, 0, 0, 1;
}
bool s_lambda_constant()
{
return true;
}
// SLIMProblemInterface interface
public:
void storeResult(Eigen::Matrix3Xd &res) {}
};
}
#ifndef SYMMETRIC_DIRICHLET_ADOLC
#define SYMMETRIC_DIRICHLET_ADOLC
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#include <iostream>
#include <cmath>
#if (COMISO_ADOLC_AVAILABLE && COMISO_EIGEN3_AVAILABLE)
#include <CoMISo/Utils/StopWatch.hh>
#include <vector>
#include <CoMISo/NSolver/LinearConstraint.hh>
#include <CoMISo/NSolver/NPDerivativeChecker.hh>
#include <CoMISo/NSolver/NProblemInterfaceADOLC.hh>
#include <CoMISo/NSolver/IPOPTSolver.hh>
#include <CoMISo/NSolver/NewtonSolver.hh>
class Symmetric_Dirichlet_ADOLC : public COMISO::NProblemInterfaceADOLC
{
public:
typedef Eigen::Matrix<adouble, Eigen::Dynamic, Eigen::Dynamic> MatrixXad;
Symmetric_Dirichlet_ADOLC() : COMISO::NProblemInterfaceADOLC(),
mesh_vertex_count(0),
param_vertex_count(0),
iteration(0)
{}
virtual int n_unknowns()
{
return n_vertices()*2;
}
virtual void initial_x(double* _x)
{
int v = n_vertices();
for (int i = 0; i < v; i++)
{
_x[i+0] = this->par(i, 0).value();
_x[i+v] = this->par(i, 1).value();
}
}
virtual adouble eval_f_adouble(const adouble* _x)
{
int v = n_vertices();
int f = n_faces();
MatrixXad par;
par.resize(v, 2);
for (int i = 0; i < v; i++)
{
par(i, 0) = _x[i+0];
par(i, 1) = _x[i+v];
}
MatrixXad jacobian;
jacobian.resize(f, 4);
jacobian.col(0) = this->dx * par.col(0);
jacobian.col(1) = this->dx * par.col(0);
jacobian.col(2) = this->dx * par.col(1);
jacobian.col(3) = this->dx * par.col(1);
adouble res = 0;
for(int i = 0; i < f; i++)
{
Eigen::Matrix<adouble, 2, 2> ji;
ji(0, 0) = jacobian(i, 0);
ji(0, 1) = jacobian(i, 1);
ji(1, 0) = jacobian(i, 2);
ji(1, 1) = jacobian(i, 3);
if (ji.determinant() <= 0)
{
return std::numeric_limits<adouble>::infinity();
}
adouble j_n = ji.norm();
adouble ji_n= ji.inverse().norm();
//res += _a(i) * (j_n + ji_n);
res += area(i) * ((j_n*j_n + ji_n*ji_n - 4));
}
return 0;
}
virtual void store_result(const double* _x)
{
int v = n_vertices();
res.resize(v, 2);
for (int i = 0; i < v; i++)
{
res(i, 0) = _x[i+0];
res(i, 1) = _x[i+v];
}
}
virtual bool constant_gradient() const { return false; }
virtual bool constant_hessian() const { return false; }
virtual bool sparse_gradient() { return false;}
virtual bool sparse_hessian () { return false;}
void init()
{
if (iteration == 0) // otherwise everthing should be set already
{
size_t f = n_faces();
size_t v = n_vertices();
this->ref.resize(v, 3);
this->par.resize(v, 2);
this->fac.resize(f, 3);
this->w.resize(f, 4);
this->dx.resize(f, v);
this->dy.resize(f, v);
this->sv.resize(f, 2);
this->lambda.resize(f, 4);
this->area.resize(f);
for (size_t i = 0; i < v; i++)
{
// adouble x, y, z;
// x = this->i_reference_mesh_vertices[i][0];
// y = this->i_reference_mesh_vertices[i][1];
// z = this->i_reference_mesh_vertices[i][2];
// this->ref.row(i) << x, y, z;
adouble u, v;
u = this->i_param_mesh_vertices[i][0];
v = this->i_param_mesh_vertices[i][1];
this->par.row(i) << u, v;
}
//Local Transformation and FE Gradient Matrices
std::vector<Eigen::Triplet<adouble>> _dx;
std::vector<Eigen::Triplet<adouble>> _dy;
for (size_t i = 0; i < f; i++)
{
int a, b, c;
a = this->i_faces[i][0];
b = this->i_faces[i][1];
c = this->i_faces[i][2];
this->fac.row(i) << a, b, c;
//FE Gradient Matrices
//Triangle Edge Vectors
auto ab = (this->ref.row(b) - this->ref.row(a)).transpose();
auto ac = (this->ref.row(c) - this->ref.row(a)).transpose();
//Length of edges
adouble ab_n = ab.norm();
adouble ac_n = ac.norm();
//Angle between edges
adouble angle = std::acos(((ab.dot(ac)) / (ab_n * ac_n)).value());
//local vector ab
adouble ab_x = ab_n;
adouble ab_y = 0;
//local vector ac
adouble ac_x = std::cos(angle.value()) * ac_n;
adouble ac_y = std::sin(angle.value()) * ac_n;
adouble tmp_area = 0.5 * ab_x * ac_y;
// assert(tmp_area > 0.0);
if (tmp_area <= 0.0)
std::cerr << "Error: area is not positive: " << tmp_area << std::endl;
area(i) = tmp_area; //
adouble p_0x = 0; //p_0
adouble p_0y = 0;
adouble p_1x = ab_x; //p_1
adouble p_1y = ab_y;
adouble p_2x = ac_x; //p_2
adouble p_2y = ac_y;
_dx.push_back(Eigen::Triplet<adouble>(i, a, (p_1y - p_2y)*(1/(2*area(i)))));
_dx.push_back(Eigen::Triplet<adouble>(i, b, (p_2y - p_0y)*(1/(2*area(i)))));
_dx.push_back(Eigen::Triplet<adouble>(i, c, (p_0y - p_1y)*(1/(2*area(i)))));
_dy.push_back(Eigen::Triplet<adouble>(i, a, (p_2x - p_1x)*(1/(2*area(i)))));
_dy.push_back(Eigen::Triplet<adouble>(i, b, (p_0x - p_2x)*(1/(2*area(i)))));
_dy.push_back(Eigen::Triplet<adouble>(i, c, (p_1x - p_0x)*(1/(2*area(i)))));
}
this->dx.setFromTriplets(_dx.begin(), _dx.end());
this->dy.setFromTriplets(_dy.begin(), _dy.end());
}// iteration == 0
iteration++;
}
int n_faces()
{
return this->i_faces.size();
}
int n_vertices()
{
return this->i_reference_mesh_vertices.size();
}
void addReferenceVertex(double _x, double _y, double _z)
{
this->i_reference_mesh_vertices.push_back({_x, _y, _z});
}
void addReferenceVertex(double _x, double _y)
{
this->i_reference_mesh_vertices.push_back({_x, _y, 0});
}
void addParamVertex(double _x, double _y)
{
this->i_param_mesh_vertices.push_back({_x, _y});
this->param_vertex_count++;
}
void addFace(int _a, int _b, int _c)
{
this->i_faces.push_back({_a, _b, _c});
}
private:
//== INPUT ====================================================================
//Base Mesh
std::vector<Eigen::Vector3d> i_reference_mesh_vertices;
int mesh_vertex_count;
//Parametrization Mesh
std::vector<Eigen::Vector2d> i_param_mesh_vertices;
int param_vertex_count;
//Triangles
std::vector<Eigen::Vector3i> i_faces;
//Result
bool resultAvailable;
int iteration;
public:
//== COMPUTE ===================================================================
Eigen::MatrixXd ref; //Reference Mesh, #V x 2
MatrixXad par; //Parametrization, #V x 2
Eigen::MatrixXi fac; //Faces, #F x 3
Eigen::MatrixXd res;
Eigen::MatrixXd w; //Weight Matrix, colums represent diagonal of W_ii, (w11 | w12 | w21 | w22) ,#F x 4
Eigen::SparseMatrix<adouble> dx; //FE Gradient Matrices, #F x #V
Eigen::SparseMatrix<adouble> dy;
Eigen::MatrixXd lambda; //Rotation depending on energy, #F x 4
Eigen::MatrixXd sv; //Singular Values of the jacobians, #F x 2
Eigen::Matrix<adouble, Eigen::Dynamic, 1> area;
};
#endif
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment