Commit a77951fc authored by Pit Nestle's avatar Pit Nestle

Bugfixes

parent 0e80f283
Pipeline #7477 failed with stages
in 7 minutes and 49 seconds
......@@ -28,7 +28,7 @@ public:
j << _j(i, 0), _j(i, 1),
_j(i, 2), _j(i, 3);
double det = _j(i, 0) * _j(i, 0) - _j(i, 1) *_j(i, 2);
double det = _j(i, 0) * _j(i, 3) - _j(i, 1) *_j(i, 2);
if (det <= 0)
return std::numeric_limits<double>::infinity();
......
......@@ -100,7 +100,7 @@ void SLIMProblemInterface::init()
this->sv.resize(f, 2);
this->lambda.resize(4*f, 1);
this->lambda.resize(f, 4);
this->area.resize(f);
......@@ -151,7 +151,14 @@ void SLIMProblemInterface::init()
double ac_x = std::cos(angle) * ac_n;
double ac_y = std::sin(angle) * ac_n;
area(i) = 0.5 * ab_x * ac_y; //
auto 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; //
double p_0x = 0; //p_0
double p_0y = 0;
......
......@@ -117,7 +117,7 @@ public:
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
Eigen::MatrixXd lambda; //Rotation depending on energy, #F x 4
Eigen::MatrixXd sv; //Singular Values of the jacobians, #F x 2
......
......@@ -16,21 +16,33 @@ namespace COMISO {
int SLIMSolver::solve(SLIMProblemInterface &problem, int iterations)
{
std::cout << "Init start" << std::endl;
problem.init(); //Some precomputations if this is the first iteration
std::cout << "Init finished" << std::endl;
Eigen::MatrixXd j;
Eigen::SparseMatrix<double> ata;
Eigen::VectorXd atb;
Eigen::VectorXd p;
Eigen::MatrixXd dir;
for (int i = 0; i < 10; i++)
for (int i = 0; i < iterations; i++)
{
std::cout << "solve():calculateJacobians()" << std::endl;
calculateJacobians(problem, problem.par, j);
std::cout << "solve():calculateWeightAndLambda()" << std::endl;
calculateWeightAndLambda(problem, j, problem.lambda);
std::cout << "solve():constructLinearSystem()" << std::endl;
constructLinearSystem(problem, ata, atb);
std::cout << "solve():solveLinearSystem()" << std::endl;
// addProximalTerm(problem, ata, atb);
solveLinearSystem(ata, atb, p);
std::cout << "solve():calculateSearchDirection()" << std::endl;
calculateSearchDirection(problem, p, dir);
std::cout << "solve():lineSearch()" << std::endl;
lineSearch(problem, dir, 1);
std::cout << "Energy : ";
std::cout << problem.computeEnergySV(problem.sv, problem.area) << std::endl;
......@@ -129,10 +141,10 @@ void SLIMSolver::calculateWeightAndLambda(SLIMProblemInterface &p, Eigen::Matrix
p.sv.row(i) << sv(0), sv(1); //Save the singular values for Energy computation
lambda(0*f + i) = lambda_i(0, 0); //small lambda in big lambda
lambda(1*f + i) = lambda_i(0, 1);
lambda(2*f + i) = lambda_i(1, 0);
lambda(3*f + i) = lambda_i(1, 1);
lambda(i, 0) = lambda_i(0, 0); //small lambda in big lambda
lambda(i, 1) = lambda_i(1, 0);
lambda(i, 2) = lambda_i(0, 1);
lambda(i, 3) = lambda_i(1, 1);
//std::cout << lambda_i << std::endl;
}
......@@ -146,7 +158,8 @@ void SLIMSolver::constructLinearSystem(SLIMProblemInterface &p, Eigen::SparseMat
Eigen::SparseMatrix<double> a;
a.resize(4*f, 2*v);
Eigen::VectorXd b = p.lambda;
Eigen::VectorXd b;
b.resize(4*f);
ata.resize(2*v, 2*v);
atb.resize(2*v);
......@@ -198,22 +211,47 @@ void SLIMSolver::constructLinearSystem(SLIMProblemInterface &p, Eigen::SparseMat
#warning possible error, W*Lmabda for the right side?
//Adjusting weights of b/r/lambda
// for (int i = 0; i < f; i++)
// {
// b(0*f+i) *= p.w(i, 0);
// b(1*f+i) *= p.w(i, 1);
// b(2*f+i) *= p.w(i, 2);
// b(3*f+i) *= p.w(i, 3);
// }
for (int i = 0; i < f; i++)
{
b(0*f+i) *= p.w(i, 0);
b(1*f+i) *= p.w(i, 1);
b(2*f+i) *= p.w(i, 2);
b(3*f+i) *= p.w(i, 3);
b(0*f+i) = p.w(i, 0)*p.lambda(i, 0)
+p.w(i, 1)*p.lambda(i, 1);
b(1*f+i) = p.w(i, 2)*p.lambda(i, 0)
+p.w(i, 3)*p.lambda(i, 1);
b(2*f+i) = p.w(i, 0)*p.lambda(i, 2)
+p.w(i, 1)*p.lambda(i, 3);
b(3*f+i) = p.w(i, 2)*p.lambda(i, 2)
+p.w(i, 3)*p.lambda(i, 3);
}
#warning proximal term and soft constraint to be added
ata = a.transpose()*a;
//ata.makeCompressed();
ata.makeCompressed();
atb = a.transpose()*b;
}
void SLIMSolver::addProximalTerm(SLIMProblemInterface &p, Eigen::SparseMatrix<double> &ata, Eigen::VectorXd &atb)
{
int v = p.n_vertices();
double l = 10e-8;
Eigen::MatrixXd id = Eigen::MatrixXd::Identity(v, v) * l;
ata = id + ata;
atb.block(0,0,v,1) += p.par.col(0) * l;
atb.block(v,0,v,1) += p.par.col(1) * l;
}
void SLIMSolver::solveLinearSystem(Eigen::SparseMatrix<double> &ata, Eigen::VectorXd &atb, Eigen::VectorXd &res)
{
Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> solver;
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
solver.compute(ata);
if (solver.info() != Eigen::Success)
......@@ -227,6 +265,8 @@ void SLIMSolver::lineSearch(SLIMProblemInterface &p, Eigen::MatrixXd &dir, doubl
double curEnergy = p.computeEnergySV(p.sv, p.area); //Reference Energy
double newEnergy;
std::printf("lineSearch(), Iteration %d, Energy %f \n", -1, curEnergy);
Eigen::MatrixXd& curParam = p.par;
Eigen::MatrixXd newParam;
int i = 0;
......@@ -236,6 +276,7 @@ void SLIMSolver::lineSearch(SLIMProblemInterface &p, Eigen::MatrixXd &dir, doubl
Eigen::MatrixXd j;
calculateJacobians(p, newParam, j);
newEnergy = p.computeEnergy(j, p.area);
std::printf("lineSearch(), Iteration %d, Energy %f \n", i, newEnergy);
alphaMax *= 0.5;
} while (newEnergy > curEnergy && ++i < 100);
p.par = newParam;
......
......@@ -68,6 +68,7 @@ private:
void calculateJacobians(COMISO::SLIMProblemInterface& p, Eigen::MatrixXd& param, Eigen::MatrixXd& j);
void calculateWeightAndLambda(COMISO::SLIMProblemInterface& p, Eigen::MatrixXd& j, Eigen::MatrixXd &lambda);
void constructLinearSystem(COMISO::SLIMProblemInterface& p, Eigen::SparseMatrix<double>& ata, Eigen::VectorXd& atb);
void addProximalTerm(COMISO::SLIMProblemInterface& p, Eigen::SparseMatrix<double>& ata, Eigen::VectorXd& atb);
void solveLinearSystem(Eigen::SparseMatrix<double>& ata, Eigen::VectorXd& atb, Eigen::VectorXd& res);
double findAlphaMax(COMISO::SLIMProblemInterface& p, Eigen::VectorXd& dir);
void lineSearch(COMISO::SLIMProblemInterface&p, Eigen::MatrixXd& dir, double alphaMax);
......
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