mirror of https://github.com/CGAL/cgal
wip: testing perfomance
This commit is contained in:
Konstantinos Katrioplas
parent
1eb7b22387
commit
802102fa2e
|
|
@ -33,6 +33,72 @@ namespace internal {
|
|||
}
|
||||
};
|
||||
|
||||
|
||||
template<typename PolygonMesh, typename VertexPointMap,
|
||||
typename CotangentValue = CGAL::internal::Cotangent_value_Meyer<PolygonMesh, VertexPointMap>>
|
||||
struct Edge_cotangent_weight : CotangentValue
|
||||
{
|
||||
Edge_cotangent_weight(PolygonMesh& pmesh_, VertexPointMap vpmap_)
|
||||
: CotangentValue(pmesh_, vpmap_)
|
||||
{}
|
||||
|
||||
PolygonMesh& pmesh()
|
||||
{
|
||||
return CotangentValue::pmesh();
|
||||
}
|
||||
|
||||
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
|
||||
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
|
||||
|
||||
double operator()(halfedge_descriptor he)
|
||||
{
|
||||
if(is_border_edge(he, pmesh()))
|
||||
{
|
||||
halfedge_descriptor h1 = next(he, pmesh());
|
||||
vertex_descriptor vs = source(he, pmesh());
|
||||
vertex_descriptor vt = target(he, pmesh());
|
||||
vertex_descriptor v1 = target(h1, pmesh());
|
||||
return (CotangentValue::operator ()(vs, v1, vt));
|
||||
}
|
||||
else
|
||||
{
|
||||
halfedge_descriptor h1 = next(he, pmesh());
|
||||
halfedge_descriptor h2 = prev(opposite(he, pmesh()), pmesh());
|
||||
vertex_descriptor vs = source(he, pmesh());
|
||||
vertex_descriptor vt = target(he, pmesh());
|
||||
vertex_descriptor v1 = target(h1, pmesh());
|
||||
vertex_descriptor v2 = source(h2, pmesh());
|
||||
return ( CotangentValue::operator()(vs, v1, vt) + CotangentValue::operator()(vs, v2, vt) ) / 2.0;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
template<typename PolygonMesh>
|
||||
struct Incident_area
|
||||
{
|
||||
Incident_area(PolygonMesh& mesh) : pmesh(mesh){}
|
||||
|
||||
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
|
||||
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
|
||||
|
||||
double operator()(halfedge_descriptor he)
|
||||
{
|
||||
halfedge_descriptor hopp = opposite(he, pmesh);
|
||||
face_descriptor f1 = face(he, pmesh);
|
||||
face_descriptor f2 = face(hopp, pmesh);
|
||||
|
||||
double A1 = f1 == boost::graph_traits<PolygonMesh>::null_face() ? 0 : face_area(f1, pmesh);
|
||||
double A2 = f2 == boost::graph_traits<PolygonMesh>::null_face() ? 0 : face_area(f2, pmesh);
|
||||
return A1 + A2;
|
||||
}
|
||||
PolygonMesh& pmesh;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -38,70 +38,12 @@
|
|||
#include <unordered_map>
|
||||
|
||||
|
||||
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/constraints_map.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Polygon_mesh_processing {
|
||||
namespace internal {
|
||||
/*
|
||||
template<typename PolygonMesh, typename VertexPointMap,
|
||||
typename CotangentValue = CGAL::internal::Cotangent_value_Meyer<PolygonMesh, VertexPointMap>>
|
||||
struct Edge_cotangent_weight : CotangentValue
|
||||
{
|
||||
Edge_cotangent_weight(PolygonMesh& pmesh_, VertexPointMap vpmap_)
|
||||
: CotangentValue(pmesh_, vpmap_)
|
||||
{}
|
||||
|
||||
PolygonMesh& pmesh()
|
||||
{
|
||||
return CotangentValue::pmesh();
|
||||
}
|
||||
|
||||
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
|
||||
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
|
||||
|
||||
double operator()(halfedge_descriptor he)
|
||||
{
|
||||
if(is_border_edge(he, pmesh()))
|
||||
{
|
||||
halfedge_descriptor h1 = next(he, pmesh());
|
||||
vertex_descriptor vs = source(he, pmesh());
|
||||
vertex_descriptor vt = target(he, pmesh());
|
||||
vertex_descriptor v1 = target(h1, pmesh());
|
||||
return (CotangentValue::operator ()(vs, v1, vt));
|
||||
}
|
||||
else
|
||||
{
|
||||
halfedge_descriptor h1 = next(he, pmesh());
|
||||
halfedge_descriptor h2 = prev(opposite(he, pmesh()), pmesh());
|
||||
vertex_descriptor vs = source(he, pmesh());
|
||||
vertex_descriptor vt = target(he, pmesh());
|
||||
vertex_descriptor v1 = target(h1, pmesh());
|
||||
vertex_descriptor v2 = source(h2, pmesh());
|
||||
return ( CotangentValue::operator()(vs, v1, vt) + CotangentValue::operator()(vs, v2, vt) ) / 2.0;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
template<typename PolygonMesh>
|
||||
struct Incident_area
|
||||
{
|
||||
Incident_area(PolygonMesh& mesh) : pmesh(mesh){}
|
||||
|
||||
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
|
||||
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
|
||||
|
||||
double operator()(halfedge_descriptor he)
|
||||
{
|
||||
halfedge_descriptor hopp = opposite(he, pmesh);
|
||||
face_descriptor f1 = face(he, pmesh);
|
||||
face_descriptor f2 = face(hopp, pmesh);
|
||||
|
||||
double A1 = f1 == boost::graph_traits<PolygonMesh>::null_face() ? 0 : face_area(f1, pmesh);
|
||||
double A2 = f2 == boost::graph_traits<PolygonMesh>::null_face() ? 0 : face_area(f2, pmesh);
|
||||
return A1 + A2;
|
||||
}
|
||||
PolygonMesh& pmesh;
|
||||
};
|
||||
*/
|
||||
template<typename PolygonMesh,
|
||||
typename VertexPointMap,
|
||||
typename VertexConstraintMap,
|
||||
|
|
@ -150,51 +92,42 @@ public:
|
|||
{
|
||||
calc_mass_matrix(D);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cerr << "compute coefficient matrix...\n";
|
||||
#endif
|
||||
|
||||
compute_coeff_matrix(A, L, D, time);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cerr << "conpute rhs...\n";
|
||||
#endif
|
||||
|
||||
compute_rhs(bx, by, bz, D);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cout<<"Done with setting up the system.\n";
|
||||
#endif
|
||||
}
|
||||
|
||||
void solve_system(const Eigen_matrix& A,
|
||||
Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz,
|
||||
const Eigen_vector& bx, const Eigen_vector& by, const Eigen_vector& bz)
|
||||
{
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cerr << "Preparing linear solver ...\n";
|
||||
#endif
|
||||
|
||||
Eigen_solver solver;
|
||||
solver.compute(A);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cerr << "solving...";
|
||||
#endif
|
||||
|
||||
Xx = solver.solve(bx);
|
||||
Xy = solver.solve(by);
|
||||
Xz = solver.solve(bz);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cerr << "ok." << std::endl;
|
||||
#endif
|
||||
|
||||
|
||||
if(solver.info() != Eigen::Success)
|
||||
{
|
||||
std::cerr << "Not Solved!" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
std::cout << "Xx old:\n";
|
||||
for(int i = 0; i < Xx.rows(); ++i)
|
||||
{
|
||||
std::cout << Xx[i] << std::endl;
|
||||
}
|
||||
std::cout <<"\n";
|
||||
*/
|
||||
}
|
||||
|
||||
void calc_stiff_matrix(Eigen_matrix& mat)
|
||||
|
|
|
|||
|
|
@ -39,71 +39,12 @@
|
|||
#include <unordered_set>
|
||||
#include <unordered_map>
|
||||
|
||||
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/constraints_map.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Polygon_mesh_processing {
|
||||
namespace internal {
|
||||
|
||||
template<typename PolygonMesh, typename VertexPointMap,
|
||||
typename CotangentValue = CGAL::internal::Cotangent_value_Meyer<PolygonMesh, VertexPointMap>>
|
||||
struct Edge_cotangent_weight : CotangentValue
|
||||
{
|
||||
Edge_cotangent_weight(PolygonMesh& pmesh_, VertexPointMap vpmap_)
|
||||
: CotangentValue(pmesh_, vpmap_)
|
||||
{}
|
||||
|
||||
PolygonMesh& pmesh()
|
||||
{
|
||||
return CotangentValue::pmesh();
|
||||
}
|
||||
|
||||
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
|
||||
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
|
||||
|
||||
double operator()(halfedge_descriptor he)
|
||||
{
|
||||
if(is_border_edge(he, pmesh()))
|
||||
{
|
||||
halfedge_descriptor h1 = next(he, pmesh());
|
||||
vertex_descriptor vs = source(he, pmesh());
|
||||
vertex_descriptor vt = target(he, pmesh());
|
||||
vertex_descriptor v1 = target(h1, pmesh());
|
||||
return (CotangentValue::operator ()(vs, v1, vt));
|
||||
}
|
||||
else
|
||||
{
|
||||
halfedge_descriptor h1 = next(he, pmesh());
|
||||
halfedge_descriptor h2 = prev(opposite(he, pmesh()), pmesh());
|
||||
vertex_descriptor vs = source(he, pmesh());
|
||||
vertex_descriptor vt = target(he, pmesh());
|
||||
vertex_descriptor v1 = target(h1, pmesh());
|
||||
vertex_descriptor v2 = source(h2, pmesh());
|
||||
return ( CotangentValue::operator()(vs, v1, vt) + CotangentValue::operator()(vs, v2, vt) ) / 2.0;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
template<typename PolygonMesh>
|
||||
struct Incident_area
|
||||
{
|
||||
Incident_area(PolygonMesh& mesh) : pmesh(mesh){}
|
||||
|
||||
typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
|
||||
typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
|
||||
|
||||
double operator()(halfedge_descriptor he)
|
||||
{
|
||||
halfedge_descriptor hopp = opposite(he, pmesh);
|
||||
face_descriptor f1 = face(he, pmesh);
|
||||
face_descriptor f2 = face(hopp, pmesh);
|
||||
|
||||
double A1 = f1 == boost::graph_traits<PolygonMesh>::null_face() ? 0 : face_area(f1, pmesh);
|
||||
double A2 = f2 == boost::graph_traits<PolygonMesh>::null_face() ? 0 : face_area(f2, pmesh);
|
||||
return A1 + A2;
|
||||
}
|
||||
PolygonMesh& pmesh;
|
||||
};
|
||||
|
||||
template<typename PolygonMesh,
|
||||
typename VertexPointMap,
|
||||
typename VertexConstraintMap,
|
||||
|
|
@ -147,62 +88,35 @@ public:
|
|||
check_face_range(face_range);
|
||||
}
|
||||
|
||||
void setup_system(Eigen_matrix& A, const Eigen_matrix& L, Eigen_matrix& D,
|
||||
void setup_system(Eigen_matrix& A,
|
||||
Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz,
|
||||
const double& time)
|
||||
{
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cerr << "compute coefficient matrix...\n";
|
||||
#endif
|
||||
|
||||
compute_coeff_matrix(A, L, D, time);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cerr << "conpute rhs...\n";
|
||||
#endif
|
||||
|
||||
compute_rhs(bx, by, bz, D);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cout<<"Done with setting up the system.\n";
|
||||
#endif
|
||||
compute_coefficient_matrix(A, time);
|
||||
compute_rhs(bx, by, bz);
|
||||
}
|
||||
|
||||
void solve_system(const Eigen_matrix& A,
|
||||
Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz,
|
||||
const Eigen_vector& bx, const Eigen_vector& by, const Eigen_vector& bz)
|
||||
{
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cerr << "Preparing linear solver ...\n";
|
||||
#endif
|
||||
SparseLinearSolver solver;
|
||||
NT Dx, Dy, Dz;
|
||||
|
||||
Eigen_solver solver;
|
||||
solver.compute(A);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cerr << "solving...";
|
||||
#endif
|
||||
|
||||
Xx = solver.solve(bx);
|
||||
Xy = solver.solve(by);
|
||||
Xz = solver.solve(bz);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_DEBUG
|
||||
std::cerr << "ok." << std::endl;
|
||||
#endif
|
||||
|
||||
if(solver.info() != Eigen::Success)
|
||||
if(!solver.linear_solver(A, bx, Xx, Dx) ||
|
||||
!solver.linear_solver(A, by, Xy, Dy) ||
|
||||
!solver.linear_solver(A, bz, Xz, Dz) )
|
||||
{
|
||||
std::cerr << "Not Solved!" << std::endl;
|
||||
std::cerr << "Could not solve linear system." << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void calculate_stiffness_matrix_elements()
|
||||
{
|
||||
CGAL_assertion(tripletList_.empty());
|
||||
tripletList.reserve(8 * nb_vert_); //estimation
|
||||
CGAL_assertion(stiffness_elements_.empty());
|
||||
stiffness_elements_.reserve(8 * nb_vert_); //estimation
|
||||
|
||||
BOOST_FOREACH(face_descriptor f, frange_)
|
||||
{
|
||||
|
|
@ -213,13 +127,13 @@ public:
|
|||
|
||||
if(!is_constrained(v_source) && !is_constrained(v_target))
|
||||
{
|
||||
auto i_source = vimap_[v_source];
|
||||
auto i_target = vimap_[v_target];
|
||||
int i_source = vimap_[v_source];
|
||||
int i_target = vimap_[v_target];
|
||||
NT Lij = weight_calculator_(hi);
|
||||
tripletList.push_back(Triplet(i_source, i_target, Lij));
|
||||
tripletList.push_back(Triplet(i_target, i_source, Lij));
|
||||
tripletList.push_back(Triplet(i_source, i_source, -Lij));
|
||||
tripletList.push_back(Triplet(i_target, i_target, -Lij));
|
||||
stiffness_elements_.push_back(Triplet(i_source, i_target, Lij));
|
||||
stiffness_elements_.push_back(Triplet(i_target, i_source, Lij));
|
||||
stiffness_elements_.push_back(Triplet(i_source, i_source, -Lij));
|
||||
stiffness_elements_.push_back(Triplet(i_target, i_target, -Lij));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -231,62 +145,105 @@ public:
|
|||
}
|
||||
|
||||
private:
|
||||
void calc_mass_matrix(Eigen_matrix& D)
|
||||
|
||||
void compute_coefficient_matrix(Eigen_matrix& A, const double& time)
|
||||
{
|
||||
CGAL_assertion(A.row_dimension() != 0);
|
||||
CGAL_assertion(A.column_dimension() != 0);
|
||||
CGAL_assertion(A.row_dimension() == nb_vert_);
|
||||
CGAL_assertion(A.column_dimension() == nb_vert_);
|
||||
|
||||
calculate_D_diagonal(diagonal_);
|
||||
CGAL_assertion(diagonal_.size() == nb_vert_);
|
||||
|
||||
|
||||
// fill A = D - time * L
|
||||
|
||||
|
||||
for(int idx = 0; idx < diagonal_.size(); ++idx)
|
||||
{
|
||||
A.set_coef(idx, idx, diagonal_[idx], true); // A is empty
|
||||
}
|
||||
|
||||
for(Triplet t : stiffness_elements_)
|
||||
{
|
||||
A.add_coef(t.get<0>(), t.get<1>(), -time * t.get<2>());
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
std::vector<bool> diag_set(diagonal_.size(), false);
|
||||
|
||||
for(Triplet t : stiffness_elements_)
|
||||
{
|
||||
if (t.get<0>() != t.get<1>())
|
||||
A.set_coef(t.get<0>(), t.get<1>(), -time * t.get<2>(), true);
|
||||
else
|
||||
{
|
||||
A.set_coef(t.get<0>(), t.get<1>(), diagonal_[t.get<0>()] -time * t.get<2>(), true);
|
||||
diag_set[t.get<0>()]=true;
|
||||
}
|
||||
}
|
||||
|
||||
for(int idx = 0; idx < diagonal_.size(); ++idx)
|
||||
{
|
||||
if (!diag_set[idx])
|
||||
{
|
||||
std::cout << "diag_NOT_set\n" ;
|
||||
A.set_coef(idx, idx, diagonal_[idx], true); // A is empty
|
||||
}
|
||||
}
|
||||
*/
|
||||
|
||||
|
||||
A.assemble_matrix(); // does setFromTriplets and some compression
|
||||
|
||||
}
|
||||
|
||||
void calculate_D_diagonal(std::vector<double>& diagonal)
|
||||
{
|
||||
diagonal.clear();
|
||||
//diagonal.resize(nb_vert_);
|
||||
diagonal.assign(nb_vert_, 0);
|
||||
for(face_descriptor f : frange_)
|
||||
{
|
||||
double area = face_area(f, mesh_);
|
||||
for(vertex_descriptor v : vertices_around_face(halfedge(f, mesh_), mesh_))
|
||||
{
|
||||
auto idx = vimap_[v];
|
||||
int idx = vimap_[v];
|
||||
if(!is_constrained(v))
|
||||
D.coeffRef(idx, idx) += 2.0 * area;
|
||||
{
|
||||
diagonal[idx] += 2.0 * area;
|
||||
}
|
||||
else
|
||||
D.coeffRef(idx, idx) = 1.0;
|
||||
{
|
||||
diagonal[idx] = 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
D /= 12.0;
|
||||
|
||||
for(int i = 0; i < diagonal.size(); ++i)
|
||||
{
|
||||
diagonal[i] /= 12.0; // constraints (= 1) as well ?
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void compute_coefficient_matrix(Eigen_matrix& A, const double& time)
|
||||
void compute_rhs(Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz)
|
||||
{
|
||||
assert(A.rows() != 0);
|
||||
assert(A.cols() != 0);
|
||||
assert(A.rows() == L.rows());
|
||||
assert(A.cols() == L.cols());
|
||||
CGAL_assertion(diagonal_.size() == nb_vert_);
|
||||
CGAL_assertion(bx.size() == nb_vert_);
|
||||
CGAL_assertion(by.size() == nb_vert_);
|
||||
CGAL_assertion(bz.size() == nb_vert_);
|
||||
|
||||
std::vector<Triple> diagonal;
|
||||
calculate_D_diagonal(diagonal);
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
//A = D - time * L;
|
||||
}
|
||||
|
||||
|
||||
void calculate_D_diagonal(std::vector<Triple>& diagonal)
|
||||
{
|
||||
// to fill
|
||||
}
|
||||
|
||||
|
||||
void compute_rhs(Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz,
|
||||
Eigen_matrix& D)
|
||||
{
|
||||
for(vertex_descriptor vi : vrange_)
|
||||
{
|
||||
int index = vimap_[vi];
|
||||
Point p = get(vpmap_, vi);
|
||||
bx.coeffRef(index) = p.x();
|
||||
by.coeffRef(index) = p.y();
|
||||
bz.coeffRef(index) = p.z();
|
||||
bx.set(index, diagonal_[index] * p.x());
|
||||
by.set(index, diagonal_[index] * p.y());
|
||||
bz.set(index, diagonal_[index] * p.z());
|
||||
}
|
||||
bx = D * bx;
|
||||
by = D * by;
|
||||
bz = D * bz;
|
||||
}
|
||||
|
||||
// update mesh
|
||||
|
|
@ -300,24 +257,14 @@ private:
|
|||
return Triangle(get(vpmap_, v1), get(vpmap_, v2), get(vpmap_, v3));
|
||||
}
|
||||
|
||||
/*
|
||||
void normalize_area(Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz)
|
||||
{
|
||||
NT surface_area = area(faces(mesh_), mesh_);
|
||||
Xx /= CGAL::sqrt(surface_area);
|
||||
Xy /= CGAL::sqrt(surface_area);
|
||||
Xz /= CGAL::sqrt(surface_area);
|
||||
}
|
||||
*/
|
||||
|
||||
void update_map(Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz)
|
||||
{
|
||||
for (vertex_descriptor v : vertices(mesh_))
|
||||
{
|
||||
int index = get(vimap_, v);
|
||||
NT x_new = Xx.coeffRef(index);
|
||||
NT y_new = Xy.coeffRef(index);
|
||||
NT z_new = Xz.coeffRef(index);
|
||||
NT x_new = Xx[index];
|
||||
NT y_new = Xy[index];
|
||||
NT z_new = Xz[index];
|
||||
put(vpmap_, v, Point(x_new, y_new, z_new));
|
||||
}
|
||||
}
|
||||
|
|
@ -346,7 +293,11 @@ private:
|
|||
PolygonMesh& mesh_;
|
||||
VertexPointMap& vpmap_;
|
||||
std::size_t nb_vert_;
|
||||
std::vector<Triplet> tripletList_;
|
||||
|
||||
// linear system data
|
||||
std::vector<Triplet> stiffness_elements_;
|
||||
std::vector<double> diagonal_; // index of vector -> index of vimap_
|
||||
|
||||
std::set<face_descriptor> frange_;
|
||||
std::set<vertex_descriptor> vrange_;
|
||||
IndexMap vimap_ = get(boost::vertex_index, mesh_);
|
||||
|
|
|
|||
|
|
@ -33,6 +33,9 @@
|
|||
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_explicit_impl.h>
|
||||
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/constraints_map.h>
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
#include <CGAL/Timer.h>
|
||||
#endif
|
||||
// new cgal traits
|
||||
#include <CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_new_impl.h>
|
||||
|
||||
|
|
@ -221,7 +224,7 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con
|
|||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
t.stop();
|
||||
std::cout << " done ("<< t.time() <<" sec)." << std::endl;
|
||||
std::cout << "Initializing..." << std::endl;
|
||||
std::cout << "Initializing...";
|
||||
t.reset(); t.start();
|
||||
#endif
|
||||
|
||||
|
|
@ -231,31 +234,61 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con
|
|||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
t.stop();
|
||||
std::cout << " done ("<< t.time() <<" sec)." << std::endl;
|
||||
std::cout << "#iter = " << nb_iterations << std::endl;
|
||||
std::cout << "Solving linear system..." << std::endl;
|
||||
std::cout << "Setup system...";
|
||||
t.reset(); t.start();
|
||||
#endif
|
||||
|
||||
for(std::size_t iter = 0; iter < nb_iterations; ++iter)
|
||||
{
|
||||
smoother.setup_system(A, stiffness_matrix, mass_matrix, bx, by, bz, time);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
t.stop();
|
||||
std::cout << " done ("<< t.time() <<" sec)." << std::endl;
|
||||
std::cout << "Solve system...";
|
||||
t.reset(); t.start();
|
||||
#endif
|
||||
|
||||
smoother.solve_system(A, Xx, Xy, Xz, bx, by, bz);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
t.stop();
|
||||
std::cout << " done ("<< t.time() <<" sec)." << std::endl;
|
||||
std::cout << "update_mesh...";
|
||||
t.reset(); t.start();
|
||||
#endif
|
||||
|
||||
|
||||
smoother.update_mesh(Xx, Xy, Xz);
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
t.stop();
|
||||
std::cout << "Shape smoothing done in ";
|
||||
std::cout << " done in ";
|
||||
std::cout << t.time() << " sec." << std::endl;
|
||||
std::cout<<std::endl;
|
||||
#endif
|
||||
}
|
||||
|
||||
template<typename PolygonMesh, typename NamedParameters>
|
||||
void smooth_along_curvature_flow(PolygonMesh& pmesh, const double& time,
|
||||
const NamedParameters& np)
|
||||
{
|
||||
smooth_along_curvature_flow(faces(pmesh), pmesh, time, np);
|
||||
}
|
||||
|
||||
template<typename PolygonMesh>
|
||||
void smooth_along_curvature_flow(PolygonMesh& pmesh, const double& time)
|
||||
{
|
||||
smooth_along_curvature_flow(faces(pmesh), pmesh, time,
|
||||
parameters::all_default());
|
||||
}
|
||||
|
||||
|
||||
// new with cgal solver traits
|
||||
template<typename PolygonMesh, typename FaceRange, typename NamedParameters>
|
||||
void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, const double& time,
|
||||
void smooth_along_curvature_flow_new(const FaceRange& faces, PolygonMesh& pmesh, const double& time,
|
||||
const NamedParameters& np)
|
||||
{
|
||||
using boost::choose_param;
|
||||
|
|
@ -263,7 +296,7 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con
|
|||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
CGAL::Timer t;
|
||||
std::cout << "Smoothing parameters...";
|
||||
std::cout << "Initializing...";
|
||||
std::cout.flush();
|
||||
t.start();
|
||||
#endif
|
||||
|
|
@ -296,7 +329,9 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con
|
|||
// implicit scheme
|
||||
#if defined(CGAL_EIGEN3_ENABLED)
|
||||
#if EIGEN_VERSION_AT_LEAST(3,2,0)
|
||||
typedef CGAL::Eigen_solver_traits<> Default_solver;
|
||||
typedef typename Eigen::SparseMatrix<double> Eigen_sparse_matrix;
|
||||
typedef typename Eigen::BiCGSTAB<Eigen_sparse_matrix, Eigen::IncompleteLUT<double> > Eigen_solver;
|
||||
typedef CGAL::Eigen_solver_traits<Eigen_solver> Default_solver;
|
||||
#else
|
||||
typedef bool Default_solver;//compilation should crash
|
||||
//if no solver is provided and Eigen version < 3.2
|
||||
|
|
@ -320,63 +355,68 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con
|
|||
typedef typename Default_solver::Vector Eigen_vector;
|
||||
|
||||
std::size_t n = static_cast<int>(vertices(pmesh).size());
|
||||
Eigen_matrix A(n, n), stiffness_matrix(n, n), mass_matrix(n, n);
|
||||
Eigen_matrix A(n, n);
|
||||
Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n);
|
||||
|
||||
internal::Shape_smoother_new<PolygonMesh, VertexPointMap, VCMap, Default_solver, GeomTraits>
|
||||
smoother(pmesh, vpmap, vcmap);
|
||||
|
||||
smoother.init_smoothing(faces);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
t.stop();
|
||||
std::cout << " done ("<< t.time() <<" sec)." << std::endl;
|
||||
std::cout << "Initializing..." << std::endl;
|
||||
std::cout << "calculate_stiffness_matrix_elements... ";
|
||||
t.reset(); t.start();
|
||||
#endif
|
||||
|
||||
smoother.init_smoothing(faces);
|
||||
smoother.calculate_stiffness_matrix_elements();
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
t.stop();
|
||||
std::cout << " done ("<< t.time() <<" sec)." << std::endl;
|
||||
std::cout << "#iter = " << nb_iterations << std::endl;
|
||||
std::cout << "Solving linear system..." << std::endl;
|
||||
t.reset(); t.start();
|
||||
#endif
|
||||
|
||||
for(std::size_t iter = 0; iter < nb_iterations; ++iter)
|
||||
{
|
||||
smoother.setup_system(A, stiffness_matrix, mass_matrix, bx, by, bz, time);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
std::cout << "setup_system... ";
|
||||
t.reset(); t.start();
|
||||
#endif
|
||||
|
||||
smoother.setup_system(A, bx, by, bz, time);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
t.stop();
|
||||
std::cout << " done ("<< t.time() <<" sec)." << std::endl;
|
||||
std::cout << "solve_system... ";
|
||||
t.reset(); t.start();
|
||||
#endif
|
||||
|
||||
smoother.solve_system(A, Xx, Xy, Xz, bx, by, bz);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
t.stop();
|
||||
std::cout << " done ("<< t.time() <<" sec)." << std::endl;
|
||||
std::cout << "update_mesh... ";
|
||||
t.reset(); t.start();
|
||||
#endif
|
||||
|
||||
smoother.update_mesh(Xx, Xy, Xz);
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
t.stop();
|
||||
std::cout << " done ("<< t.time() <<" sec)." << std::endl;
|
||||
#endif
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
t.stop();
|
||||
std::cout << "Shape smoothing done in ";
|
||||
std::cout << t.time() << " sec." << std::endl;
|
||||
std::cout<<std::endl;
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename PolygonMesh, typename NamedParameters>
|
||||
void smooth_along_curvature_flow(PolygonMesh& pmesh, const double& time,
|
||||
const NamedParameters& np)
|
||||
{
|
||||
smooth_along_curvature_flow(faces(pmesh), pmesh, time, np);
|
||||
}
|
||||
|
||||
template<typename PolygonMesh>
|
||||
void smooth_along_curvature_flow(PolygonMesh& pmesh, const double& time)
|
||||
{
|
||||
smooth_along_curvature_flow(faces(pmesh), pmesh, time,
|
||||
parameters::all_default());
|
||||
}
|
||||
|
||||
// demo helpers, undocumented
|
||||
template<typename PolygonMesh, typename FaceRange, typename NamedParameters>
|
||||
void setup_mcf_system(const FaceRange& faces, PolygonMesh& mesh,
|
||||
|
|
|
|||
|
|
@ -1,10 +1,11 @@
|
|||
#define CGAL_PMP_SMOOTHING_VERBOSE 1
|
||||
#define CGAL_PMP_SMOOTHING_DEBUG 1
|
||||
#include <fstream>
|
||||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
#include <CGAL/Polygon_mesh_processing/shape_smoothing.h>
|
||||
#include <boost/graph/graph_traits.hpp>
|
||||
|
||||
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
|
||||
typedef Kernel::Point_3 Point;
|
||||
typedef CGAL::Surface_mesh<Point> Mesh;
|
||||
|
|
@ -140,6 +141,29 @@ void test_curvature_flow(const char* filename)
|
|||
#endif
|
||||
}
|
||||
|
||||
void test_curvature_flow_new(const char* filename)
|
||||
{
|
||||
std::ifstream input(filename);
|
||||
Mesh mesh;
|
||||
input >> mesh;
|
||||
input.close();
|
||||
|
||||
boost::property_map<Mesh, CGAL::vertex_point_t>::type vpmap =
|
||||
get(CGAL::vertex_point, mesh);
|
||||
|
||||
const double time_step = 1.0;
|
||||
CGAL::Polygon_mesh_processing::smooth_along_curvature_flow_new(faces(mesh), mesh, time_step,
|
||||
CGAL::Polygon_mesh_processing::parameters::all_default());
|
||||
|
||||
#ifdef CGAL_PMP_SMOOTHING_VERBOSE
|
||||
std::ofstream out("data/output_precision_pyramid_new.off");
|
||||
out << mesh;
|
||||
out.close();
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
|
||||
void test_demo_helpers(const char* filename)
|
||||
{
|
||||
std::ifstream input(filename);
|
||||
|
|
@ -168,11 +192,12 @@ int main(int argc, char* argv[])
|
|||
const char* filename_devil = "data/mannequin-devil.off";
|
||||
const char* filename_pyramid = "data/simple_pyramid.off";
|
||||
|
||||
test_demo_helpers(filename_devil);
|
||||
test_curvature_flow_time_step(filename_devil);
|
||||
test_curvature_flow(filename_pyramid);
|
||||
test_implicit_constrained(filename_devil);
|
||||
test_explicit_scheme(filename_devil);
|
||||
//test_demo_helpers(filename_devil);
|
||||
//test_curvature_flow_time_step(filename_devil);
|
||||
test_curvature_flow(filename_devil);
|
||||
test_curvature_flow_new(filename_devil);
|
||||
//test_implicit_constrained(filename_devil);
|
||||
//test_explicit_scheme(filename_devil);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue