Commit ab635677 authored by Pit Nestle's avatar Pit Nestle

slim solver and problem interface implemented, symmetric dirichlet example added

parent c9cbb879
Pipeline #7474 failed with stages
in 7 minutes and 56 seconds
......@@ -535,6 +535,9 @@ if (COMISO_BUILD_EXAMPLES )
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/distortion_adolc/CMakeLists.txt" )
add_subdirectory (Examples/distortion_adolc)
endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/slim_symmetric_dirichlet/CMakeLists.txt" )
add_subdirectory (Examples/slim_symmetric_dirichlet)
endif()
endif (COMISO_BUILD_EXAMPLES )
......
include (CoMISoExample)
acg_add_executable (slim_symmetric_dirichlet ${sources} ${headers} )
# enable rpath linking
set_target_properties(slim_symmetric_dirichlet PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
target_link_libraries (slim_symmetric_dirichlet
CoMISo
${COMISO_LINK_LIBRARIES}
)
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#include <iostream>
#include <cmath>
#if (COMISO_EIGEN3_AVAILABLE)
#include <Eigen/Eigen>
#include <CoMISo/Utils/StopWatch.hh>
#include <vector>
#include "CoMISo/NSolver/SLIMProblemInterface.hh"
#include "CoMISo/NSolver/SLIMSolver.hh"
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 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::Matrix2Xd& _sv, Eigen::Matrix2Xd& _g)
{
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) {}
};
// Example main
int main(void)
{
SLIM_Symmetric_Dirichlet p;
p.addReferenceVertex(0, 0);
p.addReferenceVertex(1, 0);
p.addReferenceVertex(0, 1);
p.addReferenceVertex(1, 0, 1);
// p.addParamVertex(-0.25, -0.25);
// p.addParamVertex(0.75, -0.25);
// p.addParamVertex(-0.25, 0.75);
p.addParamVertex(0, 0);
p.addParamVertex(2, 0);
p.addParamVertex(0, 2);
p.addParamVertex(2, 2);
p.addFace(0, 1, 2);
p.addFace(3, 1, 2);
COMISO::SLIMSolver ss;
ss.solve(p);
return 0;
}
#else
int main(void)
{
std::cerr << "Warning: Example cannot be executed since EIGEN3 is not available..." << std::endl;
return 0;
}
#endif
//=============================================================================
//
// CLASS SLIMProblemInterface - IMPLEMENTATION
//
//=============================================================================
//== INCLUDES =================================================================
#include "SLIMProblemInterface.hh"
//== NAMESPACES ===============================================================
namespace COMISO {
//== IMPLEMENTATION ==========================================================
int SLIMProblemInterface::n_faces()
{
return this->i_faces.size();
}
int SLIMProblemInterface::n_vertices()
{
return this->i_reference_mesh_vertices.size();
}
void SLIMProblemInterface::addReferenceVertex(double _x, double _y, double _z)
{
this->i_reference_mesh_vertices.push_back({_x, _y, _z});
}
void SLIMProblemInterface::addReferenceVertex(double _x, double _y)
{
this->i_reference_mesh_vertices.push_back({_x, _y, 0});
}
void SLIMProblemInterface::addParamVertex(double _x, double _y)
{
this->i_param_mesh_vertices.push_back({_x, _y});
this->param_vertex_count++;
}
void SLIMProblemInterface::addFace(int _a, int _b, int _c)
{
this->i_faces.push_back({_a, _b, _c});
}
void SLIMProblemInterface::getReferenceVertex(size_t i, double &_x, double &_y, double &_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);
}
void SLIMProblemInterface::getParamVertex(size_t i, double &_x, double &_y)
{
_x = this->i_param_mesh_vertices[i](0);
_y = this->i_param_mesh_vertices[i](1);
}
void SLIMProblemInterface::getFace(size_t i, int &_a, int &_b, int &_c)
{
_a = this->i_faces[i](0);
_b = this->i_faces[i](1);
_c = this->i_faces[i](2);
}
void SLIMProblemInterface::storeResult(Eigen::MatrixXd &res)
{
this->res = res;
this->resultAvailable = true;
}
bool SLIMProblemInterface::isResultAvailable()
{
return this->resultAvailable;
}
void SLIMProblemInterface::getResVertex(size_t &i, double &x, double &y)
{
x = this->res.col(i)[0];
y = this->res.col(i)[1];
}
void SLIMProblemInterface::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(4*f, 1);
this->area.resize(f);
for (size_t i = 0; i < v; i++)
{
double 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;
double 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<double>> _dx;
std::vector<Eigen::Triplet<double>> _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
double ab_n = ab.norm();
double ac_n = ac.norm();
//Angle between edges
double angle = std::acos((ab.dot(ac)) / (ab_n * ac_n));
//local vector ab
double ab_x = ab_n;
double ab_y = 0;
//local vector ac
double ac_x = std::cos(angle) * ac_n;
double ac_y = std::sin(angle) * ac_n;
area(i) = 0.5 * ab_x * ac_y; //
double p_0x = 0; //p_0
double p_0y = 0;
double p_1x = ab_x; //p_1
double p_1y = ab_y;
double p_2x = ac_x; //p_2
double p_2y = ac_y;
_dx.push_back(Eigen::Triplet<double>(i, a, (p_1y - p_2y)*(1/(2*area(i)))));
_dx.push_back(Eigen::Triplet<double>(i, b, (p_2y - p_0y)*(1/(2*area(i)))));
_dx.push_back(Eigen::Triplet<double>(i, c, (p_0y - p_1y)*(1/(2*area(i)))));
_dy.push_back(Eigen::Triplet<double>(i, a, (p_1x - p_2x)*(1/(2*area(i)))));
_dy.push_back(Eigen::Triplet<double>(i, b, (p_2x - p_0x)*(1/(2*area(i)))));
_dy.push_back(Eigen::Triplet<double>(i, c, (p_0x - p_1x)*(1/(2*area(i)))));
}
this->dx.setFromTriplets(_dx.begin(), _dx.end());
this->dy.setFromTriplets(_dy.begin(), _dy.end());
}// iteration == 0
iteration++;
}
//=============================================================================
} // namespace COMISO
//=============================================================================
......@@ -39,7 +39,7 @@ namespace COMISO {
/** \class NProblemInterface NProblemGmmInterface.hh <COMISO/NSolver/NPRoblemInterface.hh>
/** \class SLIMProblemInterface SLIMProblemInterface.hh <COMISO/NSolver/SLIMProblemInterface.hh>
Brief Description.
......@@ -50,27 +50,78 @@ class COMISODLLEXPORT SLIMProblemInterface
public:
/// Default constructor
SLIMProblemInterface();
SLIMProblemInterface() : iteration(0) {}
/// Destructor
virtual ~SLIMProblemInterface();
virtual ~SLIMProblemInterface() {};
int n_triangles();
int n_faces();
int n_vertices();
void addReferenceVertex(double _x, double _y, double _z);
void addReferenceVertex(double _x, double _y);
void addParamVertex(double _x, double _y);
void addTriangle(int _a, int _b, int _c);
void addFace(int _a, int _b, int _c);
void getReferenceVertex(size_t i, double& _x, double& _y, double& _z);
void getParamVertex(size_t i, double& _x, double& _y);
void getFace(size_t i, int& _a, int& _b, int& _c);
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;
virtual void computeEnergyGradientSV(Eigen::Matrix2Xd& _sv, Eigen::Matrix2Xd& _g) = 0;
virtual double computeEnergyGradientSV_i(double _sv) = 0;
virtual double s_lambda_i(double sigma_i) = 0;
virtual void s_lambda(double sigma_1, double sigma_2, Eigen::Matrix2d& s_lambda) = 0;
virtual bool s_lambda_constant() { return false; }
virtual void storeResult(Eigen::MatrixXd& res);
bool isResultAvailable();
void getResVertex(size_t& i, double& x, double& y);
void init(); //Called by the solver to format the data into Eigen::Matrix etc.
int getIteration();
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;
Eigen::MatrixXd res;
int iteration;
public:
//== COMPUTE ===================================================================
Eigen::MatrixXd ref; //Reference Mesh, #V x 2
Eigen::MatrixXd par; //Parametrization, #V x 2
Eigen::MatrixXi fac; //Faces, #F x 3
Eigen::MatrixXd w; //Weight Matrix, colums represent diagonal of W_ii, (w11 | w12 | w21 | w22) ,#F x 4
Eigen::SparseMatrix<double> dx; //FE Gradient Matrices, #F x #V
Eigen::SparseMatrix<double> dy;
Eigen::MatrixXd lambda; //Rotation depending on energy, #F*4 x 1
void getJacobian(int _f, Eigen::Matrix2d& j);
void getRotation(int _f, Eigen::Matrix2d& r);
void getWeightMat(int _f, Eigen::Matrix2d& w);
void getSingularValues(int _f, double& s1, double& s2);
Eigen::MatrixXd sv; //Singular Values of the jacobians, #F x 2
virtual double evaluateEnergy() = 0;
virtual double evaluateGradient() = 0;
Eigen::VectorXd area;
};
......
//=============================================================================
//
// CLASS SLIMSolver - IMPLEMENTATION
//
//=============================================================================
//== INCLUDES =================================================================
#include "SLIMSolver.hh"
//== NAMESPACES ===============================================================
namespace COMISO {
//== IMPLEMENTATION ==========================================================
int SLIMSolver::solve(SLIMProblemInterface &problem, int iterations)
{
problem.init(); //Some precomputations if this is the first iteration
Eigen::MatrixXd j;
Eigen::SparseMatrix<double> ata;
Eigen::VectorXd atb;
Eigen::VectorXd p;
Eigen::MatrixXd dir;
for (int i = 0; i < 6; i++)
{
calculateJacobians(problem, problem.par, j);
calculateWeightAndLambda(problem, j, problem.lambda);
constructLinearSystem(problem, ata, atb);
solveLinearSystem(ata, atb, p);
calculateSearchDirection(problem, p, dir);
lineSearch(problem, dir, 1);
std::cout << "Energy : ";
std::cout << problem.computeEnergySV(problem.sv, problem.area) << std::endl;
}
// calculateJacobians(problem, problem.par, j);
// calculateWeightAndLambda(problem, j, problem.lambda);
// constructLinearSystem(problem, ata, atb);
// solveLinearSystem(ata, atb, p);
// calculateSearchDirection(problem, p, dir);
// lineSearch(problem, dir, 1);
// std::cout << "### Jacobian ###" << std::endl;
// std::cout << j << std::endl;
// std::cout << "### SV ###" << std::endl;
// std::cout << problem.sv << std::endl;
// std::cout << "### W ###" << std::endl;
// std::cout << problem.w << std::endl;
// std::cout << "### LAMBDA ###" << std::endl;
// std::cout << problem.lambda << std::endl;
// std::cout << "### Energy ###" << std::endl;
// std::cout << problem.computeEnergySV(problem.sv, problem.area) << std::endl;
// std::cout << "### DIR ###" << std::endl;
// std::cout << p << std::endl;
std::cout << "### RES ###" << std::endl;
std::cout << problem.par << std::endl;
return -1;
}
int SLIMSolver::solve(SLIMProblemInterface &problem)
{
return this->solve(problem, this->max_iters_);
}
void SLIMSolver::calculateJacobians(COMISO::SLIMProblemInterface &p, Eigen::MatrixXd& param, Eigen::MatrixXd &j)
{
size_t f = p.n_faces();
j.resize(f, 4);
j.col(0) = p.dx * param.col(0); // j_00
j.col(1) = p.dy * param.col(0); // j_01
j.col(2) = p.dx * param.col(1); // j_10
j.col(3) = p.dy * param.col(1); // j_11
}
void SLIMSolver::calculateWeightAndLambda(SLIMProblemInterface &p, Eigen::MatrixXd &j, Eigen::MatrixXd& lambda)
{
size_t f = p.n_faces();
lambda.resize(4*f, 1);
for (size_t i = 0; i < f; i++)
{
int a = p.fac(i, 0);
int b = p.fac(i, 1);
int c = p.fac(i, 2);
Eigen::Matrix2d j_i(2, 2);
j_i << j(i, 0), j(i, 1), j(i, 2), j(i, 3);
Eigen::JacobiSVD<Eigen::Matrix2d> svd(j_i, Eigen::ComputeFullU | Eigen::ComputeFullV);
Eigen::Matrix2d u_i = svd.matrixU();
Eigen::Matrix2d v_i = svd.matrixV(); //Possible Error
Eigen::Vector2d sv;
Eigen::Vector2d sl;
Eigen::Vector2d sw;
sv = svd.singularValues();
sl << p.s_lambda_i(sv(0)),
p.s_lambda_i(sv(1));
if (sv(0) == sl(0))
sw(0) = 1;
else
sw(0) = std::sqrt(p.computeEnergyGradientSV_i(sv(0)) / (sv(0) - sl(0)));
if (sv(1) == sl(1))
sw(1) = 1;
else
sw(1) = std::sqrt(p.computeEnergyGradientSV_i(sv(1)) / (sv(1) - sl(1)));
// sw << std::sqrt(p.computeEnergyGradientSV_i(sv(0)) / (sv(0) - sl(0))),
// std::sqrt(p.computeEnergyGradientSV_i(sv(1)) / (sv(1) - sl(1)));
Eigen::Matrix2d w_i = u_i * sw.asDiagonal() * u_i.transpose(); // U*Sw*U^t
Eigen::Matrix2d lambda_i = u_i * sl.asDiagonal() * v_i.transpose(); // U*Sl*V^t
p.w.row(i) << w_i(0,0), w_i(0,1), w_i(1,0), w_i(1,1); //Fill Weight Matrix in large matrix
p.sv.row(i) << sv(0), sv(1); //Save the singular values for Energy computation
lambda(i*4 + 0) = lambda_i(0, 0); //small lambda in big lambda
lambda(i*4 + 1) = lambda_i(0, 1);
lambda(i*4 + 2) = lambda_i(1, 0);
lambda(i*4 + 3) = lambda_i(1, 1);
//std::cout << lambda_i << std::endl;
}
}
void SLIMSolver::constructLinearSystem(SLIMProblemInterface &p, Eigen::SparseMatrix<double> &ata, Eigen::VectorXd &atb)
{
int f = p.n_faces();
int v = p.n_vertices();
Eigen::SparseMatrix<double> a;
a.resize(4*f, 2*v);
Eigen::VectorXd b = p.lambda;
ata.resize(2*v, 2*v);
atb.resize(2*v);
//Building of through triplets, since i have no idea how else without mega mat. multiplication
std::vector<Eigen::Triplet<double>> a_trip;
a_trip.reserve(p.dx.size()*4 + p.dy.size()*4); //The D Matrices are only weighted, therefore no chnage in size.
for (int k = 0; k < p.dx.outerSize(); k++)
{
for (Eigen::SparseMatrix<double>::InnerIterator it(p.dx, k); it; ++it)
{
int col = it.col();
int row = it.row();
double w11 = p.w(row, 0);
double w12 = p.w(row, 1);
double w21 = p.w(row, 2);
double w22 = p.w(row, 3);
double val = it.value();
a_trip.push_back(Eigen::Triplet<double>(row+0, col+0, val*w11));
a_trip.push_back(Eigen::Triplet<double>(row+0, col+v, val*w12));
a_trip.push_back(Eigen::Triplet<double>(row+f, col+0, val*w21));
a_trip.push_back(Eigen::Triplet<double>(row+f, col+v, val*w22));
}
}
for (int k = 0; k < p.dy.outerSize(); k++)
{
for (Eigen::SparseMatrix<double>::InnerIterator it(p.dy, k); it; ++it)
{
int col = it.col();
int row = it.row();
double w11 = p.w(row, 0);
double w12 = p.w(row, 1);
double w21 = p.w(row, 2);
double w22 = p.w(row, 3);
double val = it.value();
a_trip.push_back(Eigen::Triplet<double>(row+2*f, col+0, val*w11));
a_trip.push_back(Eigen::Triplet<double>(row+2*f, col+v, val*w12));
a_trip.push_back(Eigen::Triplet<double>(row+3*f, col+0, val*w21));
a_trip.push_back(Eigen::Triplet<double>(row+3*f, col+v, val*w22));
}
}
a.setFromTriplets(a_trip.begin(), a_trip.end());
#warning possible error, W*Lmabda for the right side?
//Adjusting weights of b/r/lambda
for (int i = 0; i < f; i++)
{
b(i*4+0) *= p.w(i, 0);
b(i*4+1) *= p.w(i, 1);
b(i*4+2) *= p.w(i, 2);
b(i*4+3) *= p.w(i, 3);
}
#warning proximal term and soft constraint to be added
ata = a.transpose()*a;
//ata.makeCompressed();
atb = a.transpose()*b;
}
void SLIMSolver::solveLinearSystem(Eigen::SparseMatrix<double> &ata, Eigen::VectorXd &atb, Eigen::VectorXd &res)
{
Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> solver;
solver.compute(ata);
if (solver.info() != Eigen::Success)
std::cerr << "SLIMSolver::solveLinearSystem(...), AtA decomposition failed" << std::endl;
res = solver.solve(atb);
}
void SLIMSolver::lineSearch(SLIMProblemInterface &p, Eigen::MatrixXd &dir, double alphaMax)
{
double curEnergy = p.computeEnergySV(p.sv, p.area); //Reference Energy
double newEnergy;
Eigen::MatrixXd& curParam = p.par;
Eigen::MatrixXd newParam;
int i = 0;
do
{
newParam = curParam + (alphaMax * dir);
Eigen::MatrixXd j;
calculateJacobians(p, newParam, j);
newEnergy = p.computeEnergy(j, p.area);
alphaMax *= 0.5;
} while (newEnergy > curEnergy && ++i < 100);
p.par = newParam;
}
void SLIMSolver::calculateSearchDirection(SLIMProblemInterface &p, Eigen::VectorXd &proxy, Eigen::MatrixXd &dir)
{
int v = p.n_vertices();
dir.resize(v, 2);
dir.col(0) = proxy.block(0,0,v,1);
dir.col(1) = proxy.block(v,0,v,1);
dir = dir - p.par;
}
double SLIMSolver::findAlphaMax(SLIMProblemInterface &p, Eigen::VectorXd& dir)
{
#warning not implemented
return 1;
}
//=============================================================================
} // namespace COMISO
//=============================================================================
......@@ -52,15 +52,46 @@ public:
/// Default constructor
SLIMSolver(const double _lambda = 1e-4, const int _max_iters = 100)
: lambda_(_lambda), max_iters_(_max_iters)
: iteration(0),
lambda_(_lambda),