testing all cases

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/constraints_map.h

@@ -73,28 +73,6 @@ namespace internal {
       }
   };
 
-  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;
-  };
-
-
 
 
 

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_impl.h

@@ -18,28 +18,27 @@
 //
 // Author(s)     : Konstantinos Katrioplas (konst.katrioplas@gmail.com)
 
-#ifndef CGAL_POLYGON_MESH_PROCESSING_CURVATURE_FLOW_IMPL_H
-#define CGAL_POLYGON_MESH_PROCESSING_CURVATURE_FLOW_IMPL_H
 
+#ifndef CURVATURE_FLOW_NEW_IMPL_H
+#define CURVATURE_FLOW_NEW_IMPL_H
+
+#include <CGAL/Polygon_mesh_processing/measure.h>
+#include <CGAL/Polygon_mesh_processing/Weights.h>
 #include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
 #include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
-#include <CGAL/Polygon_mesh_processing/Weights.h>
+#include <CGAL/Polygon_mesh_processing/internal/Smoothing/constraints_map.h>
-#include <CGAL/Polygon_mesh_processing/measure.h>
 #include <boost/graph/graph_traits.hpp>
 #include <boost/property_map/property_map.hpp>
-
-#include <CGAL/Eigen_solver_traits.h>
-#include <CGAL/Eigen_matrix.h>
 #include <CGAL/barycenter.h>
-#include <Eigen/Sparse>
+#include <CGAL/utility.h>
-
+#if defined(CGAL_EIGEN3_ENABLED)
+#include <CGAL/Eigen_solver_traits.h>
+#endif
 #include <fstream>
 #include <unordered_set>
 #include <unordered_map>
 
 
-#include <CGAL/Polygon_mesh_processing/internal/Smoothing/constraints_map.h>
-
 namespace CGAL {
 namespace Polygon_mesh_processing {
 namespace internal {
@@ -47,6 +46,7 @@ namespace internal {
 template<typename PolygonMesh,
          typename VertexPointMap,
          typename VertexConstraintMap,
+         typename SparseLinearSolver,
          typename GeomTraits>
 class Shape_smoother{
 
@@ -54,30 +54,26 @@ private:
 
   typedef typename GeomTraits::FT NT;
   typedef typename GeomTraits::Point_3 Point;
-  typedef typename GeomTraits::Triangle_3 Triangle;
+  typedef CGAL::Triple<int, int, double> Triplet;
 
   typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
   typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
   typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
   typedef typename boost::graph_traits<PolygonMesh>::edge_descriptor edge_descriptor;
-
   typedef typename boost::property_map<PolygonMesh, boost::vertex_index_t>::type IndexMap;
 
-  typedef typename Eigen::SparseMatrix<double> Eigen_matrix;
+  // linear system
-  typedef typename Eigen::VectorXd Eigen_vector;
+  typedef typename SparseLinearSolver::Matrix Eigen_matrix;
-
+  typedef typename SparseLinearSolver::Vector Eigen_vector;
-  typedef typename Eigen::BiCGSTAB<Eigen_matrix, Eigen::IncompleteLUT<double> > Eigen_solver;
-  // typedef typename Eigen::SimplicialLDLT<Eigen_matrix> Eigen_solver;
 
 public:
   Shape_smoother(PolygonMesh& mesh,
-                 VertexPointMap& vpmap,
+                     VertexPointMap& vpmap,
-                 VertexConstraintMap& vcmap) :
+                     VertexConstraintMap& vcmap) :
       mesh_(mesh),
       vpmap_(vpmap),
       vcmap_(vcmap),
       weight_calculator_(mesh, vpmap),
-      inc_areas_calculator_(mesh),
       nb_vert_(static_cast<int>(vertices(mesh).size())) {}
 
   template<typename FaceRange>
@@ -86,216 +82,168 @@ 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,
+                    std::vector<Triplet>& stiffness_elements,
                     const double& time)
   {
-    calc_mass_matrix(D);
+    compute_coefficient_matrix(A, stiffness_elements, time);
-
+    compute_rhs(bx, by, bz);
-    compute_coeff_matrix(A, L, D, time);
-
-    compute_rhs(bx, by, bz, D);
-
   }
 
   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)
   {
+    SparseLinearSolver solver;
+    NT D;
 
-    Eigen_solver solver;
+    // calls compute once to factorize with the preconditioner
-    solver.compute(A);
+    if(!solver.factor(A, D))
-
-
-    Xx = solver.solve(bx);
-    Xy = solver.solve(by);
-    Xz = solver.solve(bz);
-
-
-
-    if(solver.info() != Eigen::Success)
     {
-      std::cerr << "Not Solved!" << std::endl;
+      std::cerr << "Could not factorize linear system with preconditioner." << std::endl;
       return;
     }
 
-
+    if(!solver.linear_solver(bx, Xx) ||
-    /*
+       !solver.linear_solver(by, Xy) ||
-    std::cout << "Xx old:\n";
+       !solver.linear_solver(bz, Xz) )
-    for(int i = 0; i < Xx.rows(); ++i)
     {
-      std::cout << Xx[i] << std::endl;
+      std::cerr << "Could not solve linear system." << std::endl;
+      return;
     }
-    std::cout <<"\n";
-    */
   }
 
-  void calc_stiff_matrix(Eigen_matrix& mat)
+  void calculate_stiffness_matrix_elements(std::vector<Triplet>& stiffness_elements)
   {
-    typedef Eigen::Triplet<double> Triplet;
+    CGAL_assertion(stiffness_elements.empty());
-    std::vector<Triplet> tripletList;
+    stiffness_elements.reserve(8 * nb_vert_); //estimation
-    tripletList.reserve(8 * nb_vert_);
-    // todo: calculate exactly how many non zero entries there will be.
 
-    for(face_descriptor f : frange_)
+    boost::unordered_map<std::size_t, NT> diag_coeff;
+    BOOST_FOREACH(face_descriptor f, frange_)
     {
-      for(halfedge_descriptor hi : halfedges_around_face(halfedge(f, mesh_), mesh_))
+      BOOST_FOREACH(halfedge_descriptor hi, halfedges_around_face(halfedge(f, mesh_), mesh_))
       {
+        halfedge_descriptor hi_opp = opposite(hi, mesh_);
+        if(!is_border(hi_opp, mesh_) && hi > hi_opp) continue;
         vertex_descriptor v_source = source(hi, mesh_);
         vertex_descriptor v_target = target(hi, mesh_);
 
         if(!is_constrained(v_source) && !is_constrained(v_target))
         {
-          auto i_source = vimap_[v_source];
+          int i_source = vimap_[v_source];
-          auto i_target = vimap_[v_target];
+          int i_target = vimap_[v_target];
           NT Lij = weight_calculator_(hi);
-          tripletList.push_back(Triplet(i_source, i_target, Lij));
+          if (!is_border(hi_opp, mesh_))
-          tripletList.push_back(Triplet(i_target, i_source, Lij));
+            Lij+= weight_calculator_(hi_opp);
-          tripletList.push_back(Triplet(i_source, i_source, -Lij));
+          stiffness_elements.push_back(Triplet(i_source, i_target, Lij));
-          tripletList.push_back(Triplet(i_target, i_target, -Lij));
+          stiffness_elements.push_back(Triplet(i_target, i_source, Lij));
+          diag_coeff.insert(std::make_pair(i_source,0)).first->second -= Lij;
+          diag_coeff.insert(std::make_pair(i_target,0)).first->second -= Lij;
         }
       }
     }
-    mat.setFromTriplets(tripletList.begin(), tripletList.end());
-
-
-
 
+    for(typename boost::unordered_map<std::size_t, NT>::iterator p = diag_coeff.begin();
+        p != diag_coeff.end(); ++p)
+    {
+      stiffness_elements.push_back(Triplet(p->first, p->first, p->second));
+    }
   }
 
   void update_mesh(Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz)
   {
-    update_map(Xx, Xy, Xz);
+    for (vertex_descriptor v : vertices(mesh_))
-
+    {
-    // normalize area
+      int index = get(vimap_, v);
-    //NT surface_area = area(faces(mesh_), mesh_);
+      NT x_new = Xx[index];
-    //std::cout << "area= " << surface_area << nl;
+      NT y_new = Xy[index];
-    //std::cout << "area sqrt= " << CGAL::sqrt(surface_area) << nl;
+      NT z_new = Xz[index];
-    //normalize_area(Xx, Xy, Xz);
+      put(vpmap_, v, Point(x_new, y_new, z_new));
-    //update_map(Xx, Xy, Xz);
+    }
-    //surface_area = area(faces(mesh_), mesh_);
-    //std::cout << "surface_area normalized= " << surface_area << nl;
   }
 
 private:
-  void calc_mass_matrix(Eigen_matrix& D)
+
+  // compute linear system
+  // -----------------------------------
+  void compute_coefficient_matrix(Eigen_matrix& A,
+                                  std::vector<Triplet>& stiffness_elements,
+                                  const double& time)
   {
-    for(face_descriptor f : frange_)
+    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
+    BOOST_FOREACH(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);
+      }
+    }
+    A.assemble_matrix(); // does setFromTriplets and some compression
+  }
+
+  void calculate_D_diagonal(std::vector<double>& diagonal)
+  {
+    diagonal.clear();
+    diagonal.assign(nb_vert_, 0);
+    std::vector<bool> constraints_flags(diagonal.size(), false);
+    BOOST_FOREACH(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;
+          constraints_flags[idx] = true;
+        }
       }
     }
-    D /= 12.0;
 
-
+    for(int i = 0; i < diagonal.size(); ++i)
-    /*
-    std::ofstream out_d("data/diag_eigen.cat");
-    for(int i = 0; i < D.rows(); ++i)
     {
-      out_d << D.coeffRef(i, i) << std::endl;
+      //if(!constraints_flags[i] == true)
+        diagonal[i] /= 12.0;
     }
-    */
-
-
-
-
   }
 
-  void compute_coeff_matrix(Eigen_matrix& A, const Eigen_matrix& L, const Eigen_matrix& D, const double& time)
+  void compute_rhs(Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz)
   {
-    assert(A.rows() != 0);
+    CGAL_assertion(diagonal_.size() == nb_vert_);
-    assert(A.cols() != 0);
+    CGAL_assertion(bx.size() == nb_vert_);
-    assert(A.rows() == L.rows());
+    CGAL_assertion(by.size() == nb_vert_);
-    assert(A.cols() == L.cols());
+    CGAL_assertion(bz.size() == nb_vert_);
-    assert(A.rows() == D.rows());
-    assert(A.cols() == D.cols());
-    A = D - time * L;
 
-    /*
+    BOOST_FOREACH(vertex_descriptor vi, vrange_)
-    std::ofstream out("data/eigen_A.dat");
-    for(auto j = 0 ; j < A.cols(); ++j)
-    {
-      for(auto i = 0; i < A.rows(); ++i)
-      {
-
-        out << A.coeffRef(i, j) << " ";
-      }
-      out << std::endl;
-    }
-    */
-
-
-  }
-
-  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();
+      bx.set(index, p.x());
-      by.coeffRef(index) = p.y();
+      by.set(index, p.y());
-      bz.coeffRef(index) = p.z();
+      bz.set(index, p.z());
-    }
-    bx = D * bx;
-    by = D * by;
-    bz = D * bz;
-
-
-    std::ofstream outbx("data/outb-eigen-x");
-    std::ofstream outby("data/outb-eigen-y");
-    std::ofstream outbz("data/outb-eigen-z");
-    for(int i = 0 ; i < bx.rows(); ++i)
-    {
-      outbx << bx[i] << std::endl;
-      outby << by[i] << std::endl;
-      outbz << bz[i] << std::endl;
-
     }
 
-
+    // multiply D * b
-
+    for(int i = 0; i < nb_vert_; ++i)
-  }
-
-  // update mesh
-  // -----------------------------------
-  Triangle triangle(face_descriptor f) const
-  {
-    halfedge_descriptor h = halfedge(f, mesh_);
-    vertex_descriptor v1  = target(h, mesh_);
-    vertex_descriptor v2  = target(next(h, mesh_), mesh_);
-    vertex_descriptor v3  = target(next(next(h, mesh_), mesh_), mesh_);
-    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);
+      bx[i] *= diagonal_[i];
-      NT x_new = Xx.coeffRef(index);
+      by[i] *= diagonal_[i];
-      NT y_new = Xy.coeffRef(index);
+      bz[i] *= diagonal_[i];
-      NT z_new = Xz.coeffRef(index);
-      put(vpmap_, v, Point(x_new, y_new, z_new));
     }
   }
 
@@ -320,15 +268,19 @@ private:
   private:
   // data members
   // ------------
+  std::size_t nb_vert_;
   PolygonMesh& mesh_;
   VertexPointMap& vpmap_;
-  std::size_t nb_vert_;
+  VertexConstraintMap& vcmap_;
+  IndexMap vimap_ = get(boost::vertex_index, mesh_);
+
+  // 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_);
-  VertexConstraintMap vcmap_;
   Edge_cotangent_weight<PolygonMesh, VertexPointMap> weight_calculator_;
-  Incident_area<PolygonMesh> inc_areas_calculator_;
 
 };
 
@@ -338,4 +290,17 @@ private:
 } // CGAL
 
 
-#endif // CGAL_POLYGON_MESH_PROCESSING_CURVATURE_FLOW_IMPL_H
+
+
+
+
+
+
+
+
+
+
+
+
+
+#endif // CURVATURE_FLOW_NEW_IMPL_H

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/curvature_flow_new_impl.h

@@ -1,305 +0,0 @@
-// Copyright (c) 2017 GeometryFactory (France).
-// All rights reserved.
-//
-// This file is part of CGAL (www.cgal.org).
-// You can redistribute it and/or modify it under the terms of the GNU
-// General Public License as published by the Free Software Foundation,
-// either version 3 of the License, or (at your option) any later version.
-//
-// Licensees holding a valid commercial license may use this file in
-// accordance with the commercial license agreement provided with the software.
-//
-// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
-// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
-//
-// $URL$
-// $Id$
-//
-//
-// Author(s)     : Konstantinos Katrioplas (konst.katrioplas@gmail.com)
-
-
-#ifndef CURVATURE_FLOW_NEW_IMPL_H
-#define CURVATURE_FLOW_NEW_IMPL_H
-
-#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
-#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
-#include <CGAL/Polygon_mesh_processing/internal/Smoothing/constraints_map.h>
-#include <CGAL/Polygon_mesh_processing/Weights.h>
-#include <CGAL/Polygon_mesh_processing/measure.h>
-#include <boost/graph/graph_traits.hpp>
-#include <boost/property_map/property_map.hpp>
-#include <CGAL/barycenter.h>
-#include <CGAL/utility.h>
-#if defined(CGAL_EIGEN3_ENABLED)
-#include <CGAL/Eigen_solver_traits.h>
-#endif
-#include <fstream>
-#include <unordered_set>
-#include <unordered_map>
-
-
-namespace CGAL {
-namespace Polygon_mesh_processing {
-namespace internal {
-
-template<typename PolygonMesh,
-         typename VertexPointMap,
-         typename VertexConstraintMap,
-         typename SparseLinearSolver,
-         typename GeomTraits>
-class Shape_smoother_new{
-
-private:
-
-  typedef typename GeomTraits::FT NT;
-  typedef typename GeomTraits::Point_3 Point;
-  typedef CGAL::Triple<int, int, double> Triplet;
-
-  typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
-  typedef typename boost::graph_traits<PolygonMesh>::halfedge_descriptor halfedge_descriptor;
-  typedef typename boost::graph_traits<PolygonMesh>::face_descriptor face_descriptor;
-  typedef typename boost::graph_traits<PolygonMesh>::edge_descriptor edge_descriptor;
-  typedef typename boost::property_map<PolygonMesh, boost::vertex_index_t>::type IndexMap;
-
-  // linear system
-  typedef typename SparseLinearSolver::Matrix Eigen_matrix;
-  typedef typename SparseLinearSolver::Vector Eigen_vector;
-
-public:
-  Shape_smoother_new(PolygonMesh& mesh,
-                     VertexPointMap& vpmap,
-                     VertexConstraintMap& vcmap) :
-      mesh_(mesh),
-      vpmap_(vpmap),
-      vcmap_(vcmap),
-      weight_calculator_(mesh, vpmap),
-      nb_vert_(static_cast<int>(vertices(mesh).size())) {}
-
-  template<typename FaceRange>
-  void init_smoothing(const FaceRange& face_range)
-  {
-    check_face_range(face_range);
-  }
-
-  void setup_system(Eigen_matrix& A,
-                    Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz,
-                    std::vector<Triplet>& stiffness_elements,
-                    const double& time)
-  {
-    compute_coefficient_matrix(A, stiffness_elements, 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)
-  {
-    SparseLinearSolver solver;
-    NT D;
-
-    // calls compute once to factorize with the preconditioner
-    if(!solver.factor(A, D))
-    {
-      std::cerr << "Could not factorize lienar system with preconditioner." << std::endl;
-    }
-
-    if(!solver.linear_solver(bx, Xx) ||
-       !solver.linear_solver(by, Xy) ||
-       !solver.linear_solver(bz, Xz) )
-    {
-      std::cerr << "Could not solve linear system." << std::endl;
-      return;
-    }
-  }
-
-  void calculate_stiffness_matrix_elements(std::vector<Triplet>& stiffness_elements)
-  {
-    CGAL_assertion(stiffness_elements.empty());
-    stiffness_elements.reserve(8 * nb_vert_); //estimation
-
-    boost::unordered_map<std::size_t, NT> diag_coeff;
-    BOOST_FOREACH(face_descriptor f, frange_)
-    {
-      BOOST_FOREACH(halfedge_descriptor hi, halfedges_around_face(halfedge(f, mesh_), mesh_))
-      {
-        halfedge_descriptor hi_opp = opposite(hi, mesh_);
-        if(!is_border(hi_opp, mesh_) && hi > hi_opp) continue;
-        vertex_descriptor v_source = source(hi, mesh_);
-        vertex_descriptor v_target = target(hi, mesh_);
-
-        if(!is_constrained(v_source) && !is_constrained(v_target))
-        {
-          int i_source = vimap_[v_source];
-          int i_target = vimap_[v_target];
-          NT Lij = weight_calculator_(hi);
-          if (!is_border(hi_opp, mesh_))
-            Lij+= weight_calculator_(hi_opp);
-          stiffness_elements.push_back(Triplet(i_source, i_target, Lij));
-          stiffness_elements.push_back(Triplet(i_target, i_source, Lij));
-          diag_coeff.insert(std::make_pair(i_source,0)).first->second -= Lij;
-          diag_coeff.insert(std::make_pair(i_target,0)).first->second -= Lij;
-        }
-      }
-    }
-
-    for(typename boost::unordered_map<std::size_t, NT>::iterator p = diag_coeff.begin();
-        p != diag_coeff.end(); ++p)
-    {
-      stiffness_elements.push_back(Triplet(p->first, p->first, p->second));
-    }
-  }
-
-  void update_mesh(Eigen_vector& Xx, Eigen_vector& Xy, Eigen_vector& Xz)
-  {
-    for (vertex_descriptor v : vertices(mesh_))
-    {
-      int index = get(vimap_, v);
-      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));
-    }
-  }
-
-private:
-
-  // compute linear system
-  // -----------------------------------
-  void compute_coefficient_matrix(Eigen_matrix& A,
-                                  std::vector<Triplet>& stiffness_elements,
-                                  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
-    BOOST_FOREACH(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);
-      }
-    }
-    A.assemble_matrix(); // does setFromTriplets and some compression
-  }
-
-  void calculate_D_diagonal(std::vector<double>& diagonal)
-  {
-    diagonal.clear();
-    diagonal.assign(nb_vert_, 0);
-    std::vector<bool> constraints_flags(diagonal.size(), false);
-    BOOST_FOREACH(face_descriptor f, frange_)
-    {
-      double area = face_area(f, mesh_);
-      for(vertex_descriptor v : vertices_around_face(halfedge(f, mesh_), mesh_))
-      {
-        int idx = vimap_[v];
-        if(!is_constrained(v))
-        {
-          diagonal[idx] += 2.0 * area;
-        }
-        else
-        {
-          diagonal[idx] = 1.0;
-          constraints_flags[idx] = true;
-        }
-      }
-    }
-
-    for(int i = 0; i < diagonal.size(); ++i)
-    {
-      if(!constraints_flags[i] == true)
-        diagonal[i] /= 12.0;
-    }
-  }
-
-  void compute_rhs(Eigen_vector& bx, Eigen_vector& by, Eigen_vector& bz)
-  {
-    CGAL_assertion(diagonal_.size() == nb_vert_);
-    CGAL_assertion(bx.size() == nb_vert_);
-    CGAL_assertion(by.size() == nb_vert_);
-    CGAL_assertion(bz.size() == nb_vert_);
-
-    BOOST_FOREACH(vertex_descriptor vi, vrange_)
-    {
-      int index = vimap_[vi];
-      Point p = get(vpmap_, vi);
-      bx.set(index, p.x());
-      by.set(index, p.y());
-      bz.set(index, p.z());
-    }
-
-    // multiply D * b
-    for(int i = 0; i < nb_vert_; ++i)
-    {
-      bx[i] *= diagonal_[i];
-      by[i] *= diagonal_[i];
-      bz[i] *= diagonal_[i];
-    }
-  }
-
-  // handling constrains
-  // -----------------------------------
-  bool is_constrained(const vertex_descriptor& v)
-  {
-    return get(vcmap_, v);
-  }
-
-  template<typename FaceRange>
-  void check_face_range(const FaceRange& face_range)
-  {
-    BOOST_FOREACH(face_descriptor f, face_range)
-    {
-      frange_.insert(f);
-      BOOST_FOREACH(vertex_descriptor v, vertices_around_face(halfedge(f, mesh_), mesh_))
-        vrange_.insert(v);
-    }
-  }
-
-  private:
-  // data members
-  // ------------
-  std::size_t nb_vert_;
-  PolygonMesh& mesh_;
-  VertexPointMap& vpmap_;
-  VertexConstraintMap& vcmap_;
-  IndexMap vimap_ = get(boost::vertex_index, mesh_);
-
-  // 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_;
-  Edge_cotangent_weight<PolygonMesh, VertexPointMap> weight_calculator_;
-
-};
-
-
-} // internal
-} // PMP
-} // CGAL
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-#endif // CURVATURE_FLOW_NEW_IMPL_H

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/shape_smoothing.h

@@ -25,7 +25,6 @@
 #include <boost/property_map/property_map.hpp>
 
 #if defined(CGAL_EIGEN3_ENABLED)
-#include <Eigen/Sparse>
 #include <CGAL/Eigen_solver_traits.h>
 #endif
 
@@ -37,8 +36,6 @@
 #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>
 
 namespace CGAL {
 namespace Polygon_mesh_processing {
@@ -179,122 +176,6 @@ void smooth_along_curvature_flow(const FaceRange& faces, PolygonMesh& pmesh, con
   using boost::choose_param;
   using boost::get_param;
 
-  #ifdef CGAL_PMP_SMOOTHING_VERBOSE
-  CGAL::Timer t;
-  std::cout << "Smoothing parameters...";
-  std::cout.flush();
-  t.start();
-  #endif
-
-  typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GeomTraits;
-
-  typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::type VertexPointMap;
-  VertexPointMap vpmap = choose_param(get_param(np, internal_np::vertex_point),
-                               get_property_map(CGAL::vertex_point, pmesh));
-
-  typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
-  typedef typename boost::lookup_named_param_def <
-      internal_np::vertex_is_constrained_t,
-      NamedParameters,
-      internal::Constrained_vertices_map<vertex_descriptor>
-    > ::type VCMap;
-  VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained),
-                             internal::Constrained_vertices_map<vertex_descriptor>());
-
-  std::size_t nb_iterations = choose_param(get_param(np, internal_np::number_of_iterations), 1);
-
-  bool use_explicit_scheme = choose_param(get_param(np, internal_np::use_explicit_scheme), false);
-
-  if(use_explicit_scheme)
-  {
-    smooth_curvature_flow_explicit(faces, pmesh, parameters::number_of_iterations(nb_iterations));
-  }
-  else
-  {
-    // implicit scheme
-    typedef typename Eigen::VectorXd Eigen_vector;
-    typedef typename Eigen::SparseMatrix<double> Eigen_matrix;
-
-    std::size_t n = static_cast<int>(vertices(pmesh).size());
-    Eigen_matrix A(n, n), stiffness_matrix(n, n), mass_matrix(n, n);
-    Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n);
-
-    internal::Shape_smoother<PolygonMesh, VertexPointMap, VCMap, GeomTraits>
-        smoother(pmesh, vpmap, vcmap);
-
-    #ifdef CGAL_PMP_SMOOTHING_VERBOSE
-    t.stop();
-    std::cout << " done ("<< t.time() <<" sec)." << std::endl;
-    std::cout << "Initializing...";
-    t.reset(); t.start();
-    #endif
-
-    smoother.init_smoothing(faces);
-    smoother.calc_stiff_matrix(stiffness_matrix);
-
-    #ifdef CGAL_PMP_SMOOTHING_VERBOSE
-    t.stop();
-    std::cout << " done ("<< t.time() <<" sec)." << 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 << " 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_new(const FaceRange& faces, PolygonMesh& pmesh, const double& time,
-                                 const NamedParameters& np)
-{
-  using boost::choose_param;
-  using boost::get_param;
-
   #ifdef CGAL_PMP_SMOOTHING_VERBOSE
   CGAL::Timer t;
   std::cout << "Initializing...";
@@ -360,7 +241,7 @@ void smooth_along_curvature_flow_new(const FaceRange& faces, PolygonMesh& pmesh,
     Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n);
     std::vector<CGAL::Triple<int, int, double> > stiffness;
 
-    internal::Shape_smoother_new<PolygonMesh, VertexPointMap, VCMap, Default_solver, GeomTraits>
+    internal::Shape_smoother<PolygonMesh, VertexPointMap, VCMap, Default_solver, GeomTraits>
         smoother(pmesh, vpmap, vcmap);
 
     smoother.init_smoothing(faces);
@@ -416,94 +297,20 @@ void smooth_along_curvature_flow_new(const FaceRange& faces, PolygonMesh& pmesh,
 }
 
 template<typename PolygonMesh, typename NamedParameters>
-void smooth_along_curvature_flow_new(PolygonMesh& pmesh, const double& time,
+void smooth_along_curvature_flow(PolygonMesh& pmesh, const double& time,
-                                     const NamedParameters& np)
+                                 const NamedParameters& np)
 {
-  smooth_along_curvature_flow_new(faces(pmesh), pmesh, time, np);
+  smooth_along_curvature_flow(faces(pmesh), pmesh, time, np);
 }
 
 template<typename PolygonMesh>
-void smooth_along_curvature_flow_new(PolygonMesh& pmesh, const double& time)
+void smooth_along_curvature_flow(PolygonMesh& pmesh, const double& time)
 {
-  smooth_along_curvature_flow_new(faces(pmesh), pmesh, time,
+  smooth_along_curvature_flow(faces(pmesh), pmesh, time,
-                                  parameters::all_default());
+                              parameters::all_default());
 }
 
 // API for Polyhedron demo plugin
-template<typename PolygonMesh, typename FaceRange, typename NamedParameters>
-void setup_mcf_system(const FaceRange& faces, PolygonMesh& mesh,
-                      Eigen::SparseMatrix<double>& stiffness_matrix, const NamedParameters& np)
-{
-  using boost::choose_param;
-  using boost::get_param;
-
-  typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GeomTraits;
-
-  typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::type VertexPointMap;
-  VertexPointMap vpmap = choose_param(get_param(np, internal_np::vertex_point),
-                               get_property_map(CGAL::vertex_point, mesh));
-
-  typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
-  typedef typename boost::lookup_named_param_def <
-      internal_np::vertex_is_constrained_t,
-      NamedParameters,
-      internal::Constrained_vertices_map<vertex_descriptor>
-    > ::type VCMap;
-  VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained),
-                             internal::Constrained_vertices_map<vertex_descriptor>());
-
-  std::size_t n = static_cast<int>(vertices(mesh).size());
-  stiffness_matrix.resize(n, n);
-
-  internal::Shape_smoother<PolygonMesh, VertexPointMap, VCMap, GeomTraits>
-      smoother(mesh, vpmap, vcmap);
-  smoother.init_smoothing(faces);
-  smoother.calc_stiff_matrix(stiffness_matrix);
-}
-
-template<typename PolygonMesh, typename FaceRange, typename NamedParameters>
-void solve_mcf_system(const FaceRange& faces, PolygonMesh& mesh, const double& time,
-                      Eigen::SparseMatrix<double>& stiffness_matrix, const NamedParameters& np)
-{
-
-  using boost::choose_param;
-  using boost::get_param;
-
-  //typedef typename GetGeomTraits<PolygonMesh>::type GeomTraits;
-  //typedef typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type VertexPointMap;
-  //VertexPointMap vpmap = get(CGAL::vertex_point, mesh);
-
-  typedef typename GetGeomTraits<PolygonMesh, NamedParameters>::type GeomTraits;
-
-  typedef typename GetVertexPointMap<PolygonMesh, NamedParameters>::type VertexPointMap;
-  VertexPointMap vpmap = choose_param(get_param(np, internal_np::vertex_point),
-                               get_property_map(CGAL::vertex_point, mesh));
-
-  typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
-  typedef typename boost::lookup_named_param_def <
-      internal_np::vertex_is_constrained_t,
-      NamedParameters,
-      internal::Constrained_vertices_map<vertex_descriptor>
-    > ::type VCMap;
-  VCMap vcmap = choose_param(get_param(np, internal_np::vertex_is_constrained),
-                             internal::Constrained_vertices_map<vertex_descriptor>());
-
-  typedef typename Eigen::VectorXd Eigen_vector;
-  typedef typename Eigen::SparseMatrix<double> Eigen_matrix;
-
-  std::size_t n = static_cast<int>(vertices(mesh).size());
-
-  Eigen_matrix A(n, n), mass_matrix(n, n);
-  Eigen_vector bx(n), by(n), bz(n), Xx(n), Xy(n), Xz(n);
-
-  internal::Shape_smoother<PolygonMesh, VertexPointMap, VCMap, GeomTraits>
-      smoother(mesh, vpmap, vcmap);
-  smoother.init_smoothing(faces);
-  smoother.setup_system(A, stiffness_matrix, mass_matrix, bx, by, bz, time);
-  smoother.solve_system(A, Xx, Xy, Xz, bx, by, bz);
-  smoother.update_mesh(Xx, Xy, Xz);
-}
-
 template<typename PolygonMesh, typename FaceRange, typename NamedParameters>
 void solve_mcf(const FaceRange& faces, PolygonMesh& mesh, const double& time,
                std::vector<CGAL::Triple<int, int, double> >& stiffness, bool compute_stiffness,
@@ -565,14 +372,14 @@ void solve_mcf(const FaceRange& faces, PolygonMesh& mesh, const double& time,
 
   if(compute_stiffness)
   {
-    internal::Shape_smoother_new<PolygonMesh, VertexPointMap, VCMap, Default_solver, GeomTraits>
+    internal::Shape_smoother<PolygonMesh, VertexPointMap, VCMap, Default_solver, GeomTraits>
         smoother(mesh, vpmap, vcmap);
     smoother.init_smoothing(faces);
     smoother.calculate_stiffness_matrix_elements(stiffness);
   }
   else
   {
-    internal::Shape_smoother_new<PolygonMesh, VertexPointMap, VCMap, Default_solver, GeomTraits>
+    internal::Shape_smoother<PolygonMesh, VertexPointMap, VCMap, Default_solver, GeomTraits>
         smoother(mesh, vpmap, vcmap);
     smoother.init_smoothing(faces);
     smoother.setup_system(A, bx, by, bz, stiffness, time);

Polygon_mesh_processing/test/Polygon_mesh_processing/test_shape_smoothing.cpp

@@ -3,6 +3,7 @@
 #include <fstream>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Surface_mesh.h>
+#include <CGAL/utility.h>
 #include <CGAL/Polygon_mesh_processing/shape_smoothing.h>
 #include <boost/graph/graph_traits.hpp>
 
@@ -47,6 +48,10 @@ public:
 
 void test_implicit_constrained(const char* filename)
 {
+  #ifdef CGAL_PMP_SMOOTHING_VERBOSE
+  std::cout << "-- test_implicit_constrained --" << std::endl;
+  #endif
+
   std::ifstream input(filename);
   Mesh mesh;
   input >> mesh;
@@ -70,7 +75,7 @@ void test_implicit_constrained(const char* filename)
   const double time_step = 1.0;
   CGAL::Polygon_mesh_processing::smooth_along_curvature_flow(mesh, time_step,
                             CGAL::Polygon_mesh_processing::parameters::vertex_is_constrained_map(vcmap).
-                                                                       number_of_iterations(5));
+                                                                       number_of_iterations(1));
   #ifdef CGAL_PMP_SMOOTHING_VERBOSE
   std::ofstream out("data/output_implicit_constrained.off");
   out << mesh;
@@ -80,6 +85,10 @@ void test_implicit_constrained(const char* filename)
 
 void test_explicit_scheme(const char* filename)
 {
+  #ifdef CGAL_PMP_SMOOTHING_VERBOSE
+  std::cout << "-- test_explicit_scheme --" << std::endl;
+  #endif
+
   std::ifstream input(filename);
   Mesh mesh;
   input >> mesh;
@@ -103,6 +112,10 @@ void test_explicit_scheme(const char* filename)
 
 void test_curvature_flow_time_step(const char* filename)
 {
+  #ifdef CGAL_PMP_SMOOTHING_VERBOSE
+  std::cout << "-- test_curvature_flow_time_step --" << std::endl;
+  #endif
+
   std::ifstream input(filename);
   Mesh mesh;
   input >> mesh;
@@ -123,6 +136,10 @@ void test_curvature_flow_time_step(const char* filename)
 
 void test_curvature_flow(const char* filename)
 {
+  #ifdef CGAL_PMP_SMOOTHING_VERBOSE
+  std::cout << "-- test_curvature_flow --" << std::endl;
+  #endif
+
   std::ifstream input(filename);
   Mesh mesh;
   input >> mesh;
@@ -141,31 +158,12 @@ 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)
 {
+  #ifdef CGAL_PMP_SMOOTHING_VERBOSE
+  std::cout << "-- test_demo_helpers --" << std::endl;
+  #endif
+
   std::ifstream input(filename);
   Mesh mesh;
   input >> mesh;
@@ -175,9 +173,16 @@ void test_demo_helpers(const char* filename)
           get(CGAL::vertex_point, mesh);
 
   const double time_step = 1e-2;
-  Eigen::SparseMatrix<double> stiff_matrix;
+  std::vector<CGAL::Triple<int, int, double> > stiffness;
-  CGAL::Polygon_mesh_processing::setup_mcf_system(faces(mesh), mesh, stiff_matrix, CGAL::Polygon_mesh_processing::parameters::all_default());
+
-  CGAL::Polygon_mesh_processing::solve_mcf_system(faces(mesh), mesh, time_step, stiff_matrix, CGAL::Polygon_mesh_processing::parameters::all_default());
+  bool compute_stiffness = true;
+  CGAL::Polygon_mesh_processing::solve_mcf(faces(mesh), mesh, time_step,
+                                           stiffness, compute_stiffness,
+                                           CGAL::Polygon_mesh_processing::parameters::all_default());
+  compute_stiffness = false;
+  CGAL::Polygon_mesh_processing::solve_mcf(faces(mesh), mesh, time_step,
+                                           stiffness, compute_stiffness,
+                                           CGAL::Polygon_mesh_processing::parameters::all_default());
 
   #ifdef CGAL_PMP_SMOOTHING_VERBOSE
   std::ofstream out("data/output_devil_demo_helpers.off");
@@ -194,9 +199,8 @@ int main(int argc, char* argv[])
 
   //test_demo_helpers(filename_devil);
   //test_curvature_flow_time_step(filename_devil);
-  test_curvature_flow(filename_devil);
+  //test_curvature_flow(filename_pyramid);
-  test_curvature_flow_new(filename_devil);
+  test_implicit_constrained(filename_devil);
-  //test_implicit_constrained(filename_devil);
   //test_explicit_scheme(filename_devil);
 
   return 0;