Commit 4c8c687d authored by Pit Nestle's avatar Pit Nestle

some changes

parent 6ee346d3
Pipeline #7481 failed with stages
in 5 minutes and 8 seconds
......@@ -26,6 +26,7 @@
#include "NProblemInterface.hh"
#include <CoMISo/Config/CoMISoDefines.hh>
#include <CoMISo/Utils/StopWatch.hh>
#include "TapeIDSingleton.hh"
......@@ -122,6 +123,8 @@ public:
}
else
{
COMISO::StopWatch sw;
sw.start();
std::cerr << "re-tape..." << std::endl;
adouble y_d = 0.0;
......@@ -135,8 +138,12 @@ public:
// Call virtual function to compute
// functional value
double t1 = sw.restart();
y_d = eval_f_adouble(xa);
double t2 = sw.stop();
std::cout << "ADOLC: " << t1 / 1e3 << "s , eval_f_double: " << t2 / 1e3 << "s" << std::endl;
y_d >>= y;
trace_off();
......@@ -159,19 +166,28 @@ public:
if(!tape_available_ || always_retape_)
{
// Evaluate original functional to get new tape
COMISO::StopWatch sw;
sw.start();
eval_f(_x);
double t1 = sw.restart();
// evaluate gradient
int ec = gradient(tape_, this->n_unknowns(), _x, _g);
double t2 = sw.stop();
std::cout << "eval_gradient_!tape : " << t1 / 1e3 << "s, " << t2 / 1e3 << "s" << std::endl;
#ifdef ADOLC_RET_CODES
std::cout << "Info: gradient() returned code " << ec << std::endl;
#endif
}
else
{
COMISO::StopWatch sw;
sw.start();
// evaluate gradient
int ec = gradient(tape_, this->n_unknowns(), _x, _g);
double t1 = sw.stop();
std::cout << "eval_gradient_!tape : " << t1 / 1e3 << "s" << std::endl;
// check if retaping is required
if(ec < 0)
{
......@@ -187,9 +203,11 @@ public:
// automatic evaluation of hessian
virtual void eval_hessian(const double* _x, SMatrixNP& _H)
{
COMISO::StopWatch sw;
_H.resize(this->n_unknowns(), this->n_unknowns());
_H.setZero();
sw.start();
// tape update required?
if(!tape_available_ || always_retape_)
eval_f(_x);
......@@ -294,6 +312,7 @@ ec = sparse_hess( tape_,
}
_H.setFromTriplets(triplets.begin(), triplets.end());
}
std::cout << "ADOLC Hessian " << sw.stop() << std::endl;
}
// automatic evaluation of hessian
......
......@@ -98,6 +98,7 @@ void SLIMProblemInterface::init()
{
if (iteration == 0) // otherwise everthing should be set already
{
std::cout << "Init start" << std::endl;
size_t f = n_faces();
size_t v = n_vertices();
......@@ -115,6 +116,7 @@ void SLIMProblemInterface::init()
this->lambda.resize(f, 4);
this->area.resize(f);
this->sqrt_area.resize(f);
for (size_t i = 0; i < v; i++)
{
......@@ -171,6 +173,8 @@ void SLIMProblemInterface::init()
std::cerr << "Error: area is not positive: " << tmp_area << std::endl;
area(i) = tmp_area; //
sqrt_area(i) = std::sqrt(tmp_area); //
totalArea += tmp_area;
double p_0x = 0; //p_0
double p_0y = 0;
......@@ -191,8 +195,8 @@ void SLIMProblemInterface::init()
this->dx.setFromTriplets(_dx.begin(), _dx.end());
this->dy.setFromTriplets(_dy.begin(), _dy.end());
std::cout << "Init finished" << std::endl;
}// iteration == 0
iteration++;
}
......
......@@ -51,7 +51,7 @@ public:
/// Default constructor
SLIMProblemInterface() : iteration(0) {}
SLIMProblemInterface() : iteration(0), totalArea(0) {}
/// Destructor
virtual ~SLIMProblemInterface() {};
......@@ -105,7 +105,6 @@ private:
bool resultAvailable;
Eigen::MatrixXd res;
int iteration;
public:
//== COMPUTE ===================================================================
......@@ -123,6 +122,9 @@ public:
Eigen::MatrixXd sv; //Singular Values of the jacobians, #F x 2
Eigen::VectorXd area;
Eigen::VectorXd sqrt_area;
double totalArea;
int iteration;
//Relation between vertices
std::vector<int> transformation;
......
......@@ -16,44 +16,54 @@ 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;
calculateJacobians(problem, problem.par, j);
if (!checkFlipFree(problem, j))
return -1;
calculateJacobians(problem, problem.par, j);
COMISO::StopWatch sw_total;
COMISO::StopWatch sw_iteration;
sw_total.start();
for (int i = 0; i < iterations; i++)
{
std::cout << "solve():calculateJacobians()" << std::endl;
sw_iteration.start();
std::cout << "Iteration : " << problem.iteration;
//std::cout << "solve():calculateJacobians()" << std::endl;
calculateJacobians(problem, problem.par, j);
std::cout << "solve():calculateWeightAndLambda()" << std::endl;
// std::cout << "solve():calculateWeightAndLambda()" << std::endl;
calculateWeightAndLambda(problem, j, problem.lambda);
std::cout << "solve():constructLinearSystem()" << std::endl;
// std::cout << "solve():constructLinearSystem()" << std::endl;
constructLinearSystem(problem, ata, atb);
std::cout << "solve():solveLinearSystem()" << std::endl;
// 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;
// std::cout << "solve():calculateSearchDirection()" << std::endl;
calculateSearchDirection(problem, p, dir);
std::cout << "solve():lineSearch()" << std::endl;
// std::cout << "solve():lineSearch()" << std::endl;
lineSearch(problem, dir, 1);
std::cout << "Energy : ";
std::cout << problem.computeEnergySV(problem.sv, problem.area) << std::endl;
// std::cout << "Energy : ";
// std::cout << problem.computeEnergySV(problem.sv, problem.area) << std::endl;
std::cout << " in " << sw_iteration.stop() / 1e3 << "s" << std::endl;
problem.iteration++;
}
double totalTime = sw_total.stop();
std::cout << "Total time for " << iterations << " took " << totalTime / 1e3 << "s";
// calculateJacobians(problem, problem.par, j);
// calculateWeightAndLambda(problem, j, problem.lambda);
// constructLinearSystem(problem, ata, atb);
......@@ -187,10 +197,14 @@ void SLIMSolver::constructLinearSystem(SLIMProblemInterface &p, Eigen::SparseMat
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));
a_trip.push_back(Eigen::Triplet<double>(row+0, col+0, val*w11*p.sqrt_area(row)));
a_trip.push_back(Eigen::Triplet<double>(row+0, col+v, val*w12*p.sqrt_area(row)));
a_trip.push_back(Eigen::Triplet<double>(row+f, col+0, val*w21*p.sqrt_area(row)));
a_trip.push_back(Eigen::Triplet<double>(row+f, col+v, val*w22*p.sqrt_area(row)));
// 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++)
......@@ -206,35 +220,29 @@ void SLIMSolver::constructLinearSystem(SLIMProblemInterface &p, Eigen::SparseMat
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_trip.push_back(Eigen::Triplet<double>(row+2*f, col+0, val*w11*p.sqrt_area(row)));
a_trip.push_back(Eigen::Triplet<double>(row+2*f, col+v, val*w12*p.sqrt_area(row)));
a_trip.push_back(Eigen::Triplet<double>(row+3*f, col+0, val*w21*p.sqrt_area(row)));
a_trip.push_back(Eigen::Triplet<double>(row+3*f, col+v, val*w22*p.sqrt_area(row)));
// 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(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)*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);
b(0*f+i) = (p.w(i, 0)*p.lambda(i, 0)
+p.w(i, 1)*p.lambda(i, 1))*std::sqrt(p.area(i));
b(1*f+i) = (p.w(i, 2)*p.lambda(i, 0)
+p.w(i, 3)*p.lambda(i, 1))*std::sqrt(p.area(i));
b(2*f+i) = (p.w(i, 0)*p.lambda(i, 2)
+p.w(i, 1)*p.lambda(i, 3))*std::sqrt(p.area(i));
b(3*f+i) = (p.w(i, 2)*p.lambda(i, 2)
+p.w(i, 3)*p.lambda(i, 3))*std::sqrt(p.area(i));
}
......@@ -353,7 +361,7 @@ 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);
std::cout << " , Energy : " << curEnergy / (p.totalArea);
Eigen::MatrixXd& curParam = p.par;
Eigen::MatrixXd newParam;
......@@ -364,11 +372,18 @@ 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);
// std::printf("lineSearch(), Iteration %d, Energy %f \n", i, newEnergy);
alphaMax *= 0.5;
} while (newEnergy > curEnergy && ++i < 200);
if (newEnergy < curEnergy)
{
p.par = newParam;
std::cout << " -> " << newEnergy / (p.totalArea) << " in " << i << " it.";
}
else
{
std::cout << " -> sab";
}
}
void SLIMSolver::calculateSearchDirection(SLIMProblemInterface &p, Eigen::VectorXd &proxy, Eigen::MatrixXd &dir)
......
......@@ -25,22 +25,46 @@ public:
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);
Eigen::Matrix2d j;// for (int i = 0; i < v; i++)
// {
// res += par.col(i, 0) * par(i, 0);
// res += par(i, 1) * par(i, 1);
// }
double j00 = _j(i, 0);
double j01 = _j(i, 1);
double j10 = _j(i, 2);
double j11 = _j(i, 3);
// 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 norm_sqrd = j00*j00
+j01*j01
+j10*j10
+j11*j11;
j00 /= det;
j01 /= -det;
j10 /= -det;
j11 /= det;
double norm_inv_sqrd = j00*j00
+j01*j01
+j10*j10
+j11*j11;
res += _a(i) * (norm_sqrd + norm_inv_sqrd - 4);
double j_n = j.norm();
double ji_n= j.inverse().norm();
// 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));
// 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;
......
......@@ -44,12 +44,15 @@ public:
}
}
inline
virtual adouble eval_f_adouble(const adouble* _x)
{
int v = n_vertices();
int f = n_faces();
adouble res = 0.0;
MatrixXad par;
par.resize(v, 2);
for (int i = 0; i < v; i++)
......@@ -60,34 +63,63 @@ public:
MatrixXad jacobian;
jacobian.resize(f, 4);
jacobian.col(0) = this->dx * par.col(0);
jacobian.col(1) = this->dx * par.col(0);
jacobian.col(1) = this->dy * par.col(0);
jacobian.col(2) = this->dx * par.col(1);
jacobian.col(3) = this->dx * par.col(1);
jacobian.col(3) = this->dy * 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);
adouble j00 = jacobian(i, 0);
adouble j01 = jacobian(i, 1);
adouble j10 = jacobian(i, 2);
adouble j11 = jacobian(i, 3);
//std::cout << "\n ########## Jacobian ##########" << std::endl;
adouble j_n = j00*j00
+j01*j01
+j10*j10
+j11*j11;
if (ji.determinant() <= 0)
adouble det = j00*j11 - j01*j10; // a*d - b*c
if (det.value() <= 0)
{
return std::numeric_limits<adouble>::infinity();
}
adouble j_n = ji.norm();
adouble ji_n= ji.inverse().norm();
j00 /= det;
j01 /= det;
j10 /= det;
j11 /= det;
adouble ji_n = j00*j00
+j01*j01
+j10*j10
+j11*j11;
//res += _a(i) * (j_n + ji_n);
res += area(i) * ((j_n*j_n + ji_n*ji_n - 4));
res += area(i) * ((j_n+ji_n) - 4);
//std::printf("i=%d, j_norm=%f, j⁻1_norm=%f, res=%f\n", i, j_n.value(), ji_n.value(), res.value());
}
return 0;
double lambda = 1e-4;
for (int i = 0; i < v; i++)
{
res += ((_x[i]-this->par(i, 0))*(_x[i]-this->par(i, 0)))*lambda;
res += ((_x[i+v]-this->par(i, 1))*(_x[i+v]-this->par(i, 1)))*lambda;
}
for (int i = 0; i < v; i++)
{
this->par(i, 0) = _x[i+0];
this->par(i, 1) = _x[i+v];
}
return res;
}
......@@ -97,8 +129,10 @@ public:
int v = n_vertices();
res.resize(v, 2);
// std::cout << "\n ########## RES ##########\n" << std::endl;
for (int i = 0; i < v; i++)
{
// printf(" %f \t %f \n", _x[i+0], _x[i+v]);
res(i, 0) = _x[i+0];
res(i, 1) = _x[i+v];
}
......@@ -107,7 +141,7 @@ public:
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;}
virtual bool sparse_hessian () { return true;}
void init()
{
......@@ -131,23 +165,32 @@ public:
this->area.resize(f);
// std::cout << "V : " << v << " | F : " << f << std::endl;
// std::cout << "\n ########## Initial Param ##########" << std::endl;
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;
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;
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;
// std::cout << par(i, 0).value() << " | " << par(i, 1).value() << std::endl;
}
//Local Transformation and FE Gradient Matrices
std::vector<Eigen::Triplet<adouble>> _dx;
std::vector<Eigen::Triplet<adouble>> _dy;
// std::cout << "\n ########## Local Calculation ##########" << std::endl;
for (size_t i = 0; i < f; i++)
{
int a, b, c;
......@@ -159,26 +202,28 @@ public:
//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();
Eigen::Vector3d ab = (this->ref.row(b) - this->ref.row(a)).transpose();
Eigen::Vector3d ac = (this->ref.row(c) - this->ref.row(a)).transpose();
//Length of edges
adouble ab_n = ab.norm();
adouble ac_n = ac.norm();
double ab_n = ab.norm();
double ac_n = ac.norm();
//Angle between edges
adouble angle = std::acos(((ab.dot(ac)) / (ab_n * ac_n)).value());
double angle = std::acos(((ab.dot(ac)) / (ab_n * ac_n)));
//local vector ab
adouble ab_x = ab_n;
adouble ab_y = 0;
double ab_x = ab_n;
double ab_y = 0;
//local vector ac
adouble ac_x = std::cos(angle.value()) * ac_n;
adouble ac_y = std::sin(angle.value()) * ac_n;
double ac_x = std::cos(angle) * ac_n;
double ac_y = std::sin(angle) * ac_n;
adouble tmp_area = 0.5 * ab_x * ac_y;
double tmp_area = 0.5 * ab_x * ac_y;
// assert(tmp_area > 0.0);
......@@ -187,12 +232,12 @@ public:
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;
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<adouble>(i, a, (p_1y - p_2y)*(1/(2*area(i)))));
......@@ -202,11 +247,22 @@ public:
_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)))));
// std::cout << "\n Triangle Edge Vectors" << std::endl;
// std::cout << "ab : ( " << ab(0) << " | " << ab(1) << " ), ac : ( " << ac(0) << " | " << ac(1) << " )" << std::endl;
// std::cout << "\n Triangle Edge Length" << std::endl;
// std::cout << "ab : " << ab_n << ", ac : " << ac_n << std::endl;
}
this->dx.setFromTriplets(_dx.begin(), _dx.end());
this->dy.setFromTriplets(_dy.begin(), _dy.end());
}// iteration == 0
}
iteration++;
}
......@@ -281,6 +337,5 @@ private:
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