...
 
Commits (10)
......@@ -2,3 +2,4 @@ Config/config.hh
*.orig
*.rej
.project
CMakeLists.txt.user
......@@ -530,9 +530,15 @@ if (COMISO_BUILD_EXAMPLES )
add_subdirectory (Examples/small_mosek_native)
endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_mosek_fusion_sdp/CMakeLists.txt" )
add_subdirectory (Examples/small_mosek_fusion_sdp)
add_subdirectory (Examples/small_mosek_fusion_sdp)
endif()
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 )
# Only create install target, when we are building CoMISo standalone
......
include (CoMISoExample)
acg_add_executable (distortion_adolc ${sources} ${headers} )
# enable rpath linking
set_target_properties(distortion_adolc PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
target_link_libraries (distortion_adolc
CoMISo
${COMISO_LINK_LIBRARIES}
)
//== 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 Problem_Distortion : public COMISO::NProblemInterfaceADOLC
{
public:
typedef Eigen::Matrix<adouble, 2, 2> Matrix2ad;
typedef Eigen::Matrix<adouble, 2, 1> Vector2ad;
//Base Mesh
std::vector<Eigen::Vector3d> mesh_vertices;
int mesh_vertex_count;
//Parametrization Mesh
std::vector<Eigen::Vector2d> param_vertices;
int param_vertex_count;
//Triangles
struct Triangle { int a,b,c; };
std::vector<Triangle> triangles;
//Precalculated local Triangle Vectors and Area
std::vector<Eigen::Matrix2d> mesh_local_triangle_vectors;
std::vector<double> mesh_triangle_area;
Problem_Distortion() : COMISO::NProblemInterfaceADOLC(),
mesh_vertex_count(0),
param_vertex_count(0)
{}
virtual int n_unknowns()
{
return param_vertex_count*2;
}
virtual void initial_x(double* _x)
{
if (this->param_vertex_count != this->mesh_vertex_count)
{
std::cerr << "Problem_Distorion::initial_x() : different vertex counts" << std::endl;
return;
}
for (int i = 0; i < this->param_vertex_count; i++)
{
_x[i*2] = this->param_vertices[i](0);
_x[i*2 + 1] = this->param_vertices[i](1);
std::cout << "Init \n" << mesh_vertices[i] << std::endl;
}
//local coordinate calculation for each triangle of the base mesh
for (Triangle t : triangles)
{
Eigen::Vector3d ab = mesh_vertices[t.b] - mesh_vertices[t.a]; // Vector 0 -> 1
Eigen::Vector3d ac = mesh_vertices[t.c] - mesh_vertices[t.a]; // Vector 0 -> 2
Eigen::Matrix2d local;
//Length of each vector
double ab_n = ab.norm();
double ac_n = ac.norm();
//Angle between ab and ac
double angle = std::acos((ab.dot(ac)) / (ab_n * ac_n));
//local ab vector
double ab_l_x = ab_n;
double ab_l_y = 0;
//local ac vector
double ac_l_x = std::cos(angle) * ac_n;
double ac_l_y = std::sin(angle) * ac_n;
//Triangle Area
double area = ab_l_x * ac_l_y * 0.5;
local << ab_l_x, ac_l_x, ab_l_y, ac_l_y; // x1, x2, y1, y2
this->mesh_local_triangle_vectors.push_back(local);
this->mesh_triangle_area.push_back(area);
}
}
virtual adouble eval_f_adouble(const adouble* _x)
{
adouble result = 0;
for (int i = 0; i < triangles.size(); i++)
{
Triangle t = triangles[i];
Vector2ad a;
Vector2ad b;
Vector2ad c;
a << _x[t.a*2],
_x[t.a*2 + 1];
b << _x[t.b*2],
_x[t.b*2 + 1];
c << _x[t.c*2],
_x[t.c*2 + 1];
Vector2ad ab = b - a;
Vector2ad ac = c - a;
Matrix2ad r;
r << ab, ac;
Matrix2ad q;
q(0, 0) = mesh_local_triangle_vectors[i](0, 0);
q(0, 1) = mesh_local_triangle_vectors[i](0, 1);
q(1, 0) = mesh_local_triangle_vectors[i](1, 0);
q(1, 1) = mesh_local_triangle_vectors[i](1, 1);
Matrix2ad j = q*r.inverse();
adouble j_norm = j.norm();
adouble j_inverse_norm = j.inverse().norm();
//Symmetric Dirichlet
adouble energy = j_norm*j_norm + j_inverse_norm*j_inverse_norm;
//
energy = (energy - 4) * mesh_triangle_area[i];
result += energy;
}
std::cout << "Energy : " << result.value() << std::endl;
return result;
}
virtual void store_result(const double* _x)
{
for (int i = 0; i < this->n_unknowns(); i++)
{
std::cout << _x[i] << std::endl;
}
}
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 true;}
void addMeshVertex(double x, double y, double z)
{
this->mesh_vertices.push_back(Eigen::Vector3d(x,y,z));
this->mesh_vertex_count++;
}
void addParamVertex(double x, double y)
{
this->param_vertices.push_back(Eigen::Vector2d(x,y));
this->param_vertex_count++;
}
bool addTriangle(int v0, int v1, int v2)
{
if (mesh_vertex_count < v0 || mesh_vertex_count < v1 || mesh_vertex_count < v2)
{
std::cerr << "Problem_Distorion::addTriangle() : Vertex index out of range." << std::endl;
return false;
}
this->triangles.push_back({v0, v1, v2});
return true;
}
};
// Example main
int main(void)
{
Problem_Distortion pd;
pd.addMeshVertex(0, 0, 0);
pd.addMeshVertex(4, 0, 0);
pd.addMeshVertex(0, 3, 0);
pd.addMeshVertex(3, 4, 0);
pd.addParamVertex(1, 1);
pd.addParamVertex(2, 1);
pd.addParamVertex(1, 2);
pd.addParamVertex(2, 2);
pd.addTriangle(0, 1, 2);
pd.addTriangle(1, 2, 3);
COMISO::NewtonSolver ns;
ns.solve(&pd);
return 0;
}
#else
int main(void)
{
std::cerr << "Warning: Example cannot be executed since either EIGEN3 or ADOLC is not available..." << std::endl;
return 0;
}
#endif
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 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)
{
_g.resize(sv.rows(), 2);
for (int i = 0; i < _sv.rows(); i++)
{
double s1 = _sv(i, 0);
double s2 = _sv(i, 1);
_g.row(i) << 2 * (s1 - std::pow(s1, -3)),
2 * (s2 - std::pow(s2, -3));
}
}
double computeEnergyGradientSV_i(double _sv)
{
return 2 * (_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
......@@ -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
......
//=============================================================================
//
// 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::fixVertex(int i)
{
this->softconstraintIndices.push_back(i);
}
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(tr);
this->translation.push_back({tu, tv});
}
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
{
std::cout << "Init start" << std::endl;
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);
this->sqrt_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;
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; //
sqrt_area(i) = std::sqrt(tmp_area); //
totalArea += tmp_area;
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_2x - p_1x)*(1/(2*area(i)))));
_dy.push_back(Eigen::Triplet<double>(i, b, (p_0x - p_2x)*(1/(2*area(i)))));
_dy.push_back(Eigen::Triplet<double>(i, c, (p_1x - p_0x)*(1/(2*area(i)))));
}
this->dx.setFromTriplets(_dx.begin(), _dx.end());
this->dy.setFromTriplets(_dy.begin(), _dy.end());
std::cout << "Init finished" << std::endl;
}// iteration == 0
}
//=============================================================================
} // namespace COMISO
//=============================================================================
//=============================================================================
//
// CLASS NProblemInterface
//
//=============================================================================
#ifndef COMISO_SLIMPROBLEMINTERFACE_HH
#define COMISO_SLIMPROBLEMINTERFACE_HH
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#if COMISO_EIGEN3_AVAILABLE
//== INCLUDES =================================================================
#include <Base/Code/Quality.hh>
#include <iostream>
#include <cfloat>
LOW_CODE_QUALITY_SECTION_BEGIN
#include <Eigen/Eigen>
#if !(EIGEN_VERSION_AT_LEAST(3,1,0))
#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
#endif
#include <Eigen/Sparse>
LOW_CODE_QUALITY_SECTION_END
#include <CoMISo/Config/CoMISoDefines.hh>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class SLIMProblemInterface SLIMProblemInterface.hh <COMISO/NSolver/SLIMProblemInterface.hh>
Brief Description.
A more elaborate description follows.
*/
class COMISODLLEXPORT SLIMProblemInterface
{
public:
/// Default constructor
SLIMProblemInterface() : iteration(0), totalArea(0) {}
/// Destructor
virtual ~SLIMProblemInterface() {};
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 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);
void fixVertex(int i);
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;
virtual void computeEnergyGradientSV(Eigen::MatrixXd& _sv, Eigen::MatrixXd& _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;
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 x 4
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;
std::vector<Eigen::Vector2d> translation;
std::vector<Eigen::Vector2i> relationIndices;
//Fixation/Softconstraints
std::vector<int> softconstraintIndices;
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_EIGEN3_AVAILABLE
//=============================================================================
#endif // COMISO_SLIMPROBLEMINTERFACE_HH defined
//=============================================================================
This diff is collapsed.
//=============================================================================
//
// CLASS SLIMSolver
//
//=============================================================================
#ifndef COMISO_SLIMSOLVER_HH
#define COMISO_SLIMSOLVER_HH
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include <CoMISo/Utils/StopWatch.hh>
#include "SLIMProblemInterface.hh"
//#include <Base/Debug/DebTime.hh>
#if COMISO_SUITESPARSE_AVAILABLE
#include <Eigen/UmfPackSupport>
#include <Eigen/CholmodSupport>
#endif
// ToDo: why is Metis not working yet?
//#if COMISO_METIS_AVAILABLE
// #include <Eigen/MetisSupport>
//#endif
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class NewtonSolver NewtonSolver.hh <ACG/.../NewtonSolver.hh>
Brief Description.
A more elaborate description follows.
*/
class COMISODLLEXPORT SLIMSolver
{
public:
/// Default constructor
SLIMSolver(const double _lambda = 1e-4, const int _max_iters = 20)
: iteration(0),
lambda_(_lambda),
max_iters_(_max_iters)
{
//Empty
}
int solve(SLIMProblemInterface& problem, int iterations);
int solve(SLIMProblemInterface& problem);
private:
void precompute();
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 addSoftConstraint(COMISO::SLIMProblemInterface& p, Eigen::SparseMatrix<double>& ata, Eigen::VectorXd& atb);
void addRelationConstraint(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);
void calculateSearchDirection(COMISO::SLIMProblemInterface& p, Eigen::VectorXd& proxy, Eigen::MatrixXd& dir);
bool checkFlipFree(COMISO::SLIMProblemInterface& p, Eigen::MatrixXd& j);
private:
Eigen::Matrix3Xd vertices;
Eigen::Matrix2Xd param;
Eigen::Matrix3Xd faces;
std::vector<Eigen::Matrix2d> localFaceVectors;
std::vector<double> faceArea;
//Iterative Data
int iteration;
//Parameters
double lambda_;
int max_iters_;
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // ACG_NEWTONSOLVER_HH defined
//=============================================================================
//== 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;// 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();
// 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();
}
}
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++)
{
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->dy * par.col(0);
jacobian.col(2) = this->dx * par.col(1);
jacobian.col(3) = this->dy * par.col(1);
for(int i = 0; i < f; i++)
{
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;
adouble det = j00*j11 - j01*j10; // a*d - b*c
if (det.value() <= 0)
{
return std::numeric_limits<adouble>::infinity();
}
j00 /= det;
j01 /= det;
j10 /= det;
j11 /= det;
adouble ji_n = j00*j00
+j01*j01
+j10*j10
+j11*j11;
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());
}
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;
}
virtual void store_result(const double* _x)
{
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];
}
}
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 true;}
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);
// std::cout << "V : " << v << " | F : " << f << std::endl;
// std::cout << "\n ########## Initial Param ##########" << std::endl;
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;
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;
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
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
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;
double 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;
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)))));
_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)))));
// 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++;
}
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