diff --git a/Polyhedron/demo/Polyhedron/Plugins/Operations_on_polyhedra/Clip_polyhedron_plugin.ui b/Polyhedron/demo/Polyhedron/Plugins/Operations_on_polyhedra/Clip_polyhedron_plugin.ui
index 3f001f77f78..7c5473bb039 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/Operations_on_polyhedra/Clip_polyhedron_plugin.ui
+++ b/Polyhedron/demo/Polyhedron/Plugins/Operations_on_polyhedra/Clip_polyhedron_plugin.ui
@@ -11,7 +11,7 @@
    </rect>
   </property>
   <property name="windowTitle">
-   <string>Clip polyhedra</string>
+   <string>Clip Polyhedra</string>
   </property>
   <widget class="QWidget" name="dockWidgetContents">
    <layout class="QVBoxLayout" name="verticalLayout_2">
@@ -20,7 +20,7 @@
       <item>
        <widget class="QLabel" name="label">
         <property name="text">
-         <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;This function allows to clip all the selected polyhedra by an half-space. The blue side of the clipping plane will be clipped out, and the yellow side will be kept.&lt;br/&gt;&lt;br/&gt;If the option &lt;span style=&quot; font-style:italic;&quot;&gt;keep closed&lt;/span&gt; is checked, the clipped part of each polyhedron will be closed. Otherwise, it will be left open.&lt;/p&gt;&lt;p&gt;&lt;br/&gt;&lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
+         <string>&lt;html&gt;&lt;head/&gt;&lt;body&gt;&lt;p&gt;This function allows to clip all the selected polyhedra against a halfspace. What is on the blue side of the clipping plane will be clipped, and what is on the yellow side will be kept.&lt;br/&gt;&lt;br/&gt;If the option &lt;span style=&quot; font-style:italic;&quot;&gt;keep closed&lt;/span&gt; is checked, the clipped part of each polyhedron will be closed, if it has a closed contour on the clipping plane. Otherwise, it will be left open.&lt;/p&gt;&lt;p&gt;&lt;br/&gt;&lt;/p&gt;&lt;/body&gt;&lt;/html&gt;</string>
         </property>
         <property name="wordWrap">
          <bool>true</bool>
diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_dialog.ui b/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_dialog.ui
index 51a81169f3d..3ddde77db51 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_dialog.ui
+++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_dialog.ui
@@ -162,6 +162,23 @@
         </property>
        </widget>
       </item>
+      <item row="4" column="1">
+       <widget class="QCheckBox" name="preserveDuplicates_checkbox">
+        <property name="text">
+         <string/>
+        </property>
+       </widget>
+      </item>
+      <item row="4" column="0">
+       <widget class="QLabel" name="preserveDuplicates_label">
+        <property name="text">
+         <string>Preserve duplicated edges</string>
+        </property>
+        <property name="alignment">
+         <set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
+        </property>
+       </widget>
+      </item>
      </layout>
     </widget>
    </item>
diff --git a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp
index e2549ee3fcf..8861170d0d4 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/PMP/Isotropic_remeshing_plugin.cpp
@@ -14,8 +14,12 @@
 #include <CGAL/iterator.h>
 #include <CGAL/Polygon_mesh_processing/remesh.h>
 #include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
+#include <CGAL/boost/graph/properties_Polyhedron_3.h>
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+#include <CGAL/utility.h>
 
 #include <boost/graph/graph_traits.hpp>
+#include <boost/unordered_set.hpp>
 #include <CGAL/property_map.h>
 
 #include <QTime>
@@ -39,6 +43,118 @@
 
 #include "ui_Isotropic_remeshing_dialog.h"
 
+// give a halfedge and a target edge length, put in `out` points
+// which the edge equally spaced such that splitting the edge
+// using the sequence of points make the edges shorter than
+// `target_length`
+template <class TriangleMesh, class PointPMap, class PointOutputIterator>
+PointOutputIterator
+sample_edge(
+  typename boost::graph_traits<TriangleMesh>::halfedge_descriptor hd,
+  TriangleMesh& triangle_mesh,
+  double target_length,
+  const PointPMap& pmap,
+  PointOutputIterator out)
+{
+  typedef typename boost::property_traits<PointPMap>::value_type Point_3;
+  typedef typename CGAL::Kernel_traits<Point_3>::Kernel::Vector_3 Vector_3;
+  typename boost::property_traits<PointPMap>::reference src=get(pmap, source(hd,triangle_mesh) );
+  typename boost::property_traits<PointPMap>::reference tgt=get(pmap, target(hd,triangle_mesh) );
+
+  double length = std::sqrt( CGAL::squared_distance(src, tgt) );
+  if ( length <= target_length ) return out;
+
+  double nb_points = std::floor( length / target_length );
+  Vector_3 unit = (tgt-src) / (nb_points+1);
+
+  for(double i=0; i<nb_points; ++i)
+    *out++=src+unit*(i+1);
+
+  return out;
+}
+
+// given a set of points that are expected to be on an edge, split
+// that edge and retriangulate the face incident to the edge
+// Points are sorted so that they are sorted from the source to the target
+// of the edge (the sequence does not contains edge endpoints)
+template <class TriangleMesh, class PointPMap, class PointRange, class EdgeOutputIterator>
+EdgeOutputIterator
+split_identical_edges(
+  typename boost::graph_traits<TriangleMesh>::halfedge_descriptor hd,
+  TriangleMesh& tm,
+  const PointPMap& pmap,
+  const PointRange& points,
+  EdgeOutputIterator out)
+{
+  typedef typename PointRange::value_type Point_3;
+  typedef boost::graph_traits<TriangleMesh> GT;
+  typedef typename GT::halfedge_descriptor halfedge_descriptor;
+
+  BOOST_FOREACH(const Point_3& p, points)
+  {
+    // split the edge
+    halfedge_descriptor new_hd=CGAL::Euler::split_edge(hd,tm);
+    // set the vertex point
+    put(pmap, target(new_hd, tm), p);
+    *out++=edge(new_hd, tm);
+  }
+  *out++=edge(hd, tm);
+  return out;
+}
+
+// HedgeRange is expected to be a range with value type being
+// std::pair<halfedge_descriptor, TriangleMesh*>
+// Given a set of halfedges representing different edges
+// but with identical endpoints, and a target edge length
+// we split all edges identically so that subedges are
+// or length <= length
+template <class HedgeRange, class Edges_to_protect>
+void split_long_duplicated_edge(const HedgeRange& hedge_range,
+                                double target_length,
+                                Edges_to_protect& edges_to_protect)
+{
+  typedef typename HedgeRange::value_type Pair;
+  typedef typename Pair::first_type halfedge_descriptor;
+  typedef typename boost::remove_pointer<
+    typename Pair::second_type>::type TriangleMesh;
+  typedef typename boost::property_map<TriangleMesh,
+    CGAL::vertex_point_t>::type PointPMap;
+  typedef typename boost::property_traits<PointPMap>::value_type Point_3;
+
+  if (hedge_range.empty()) return;
+
+  const Pair& p = *hedge_range.begin();
+  PointPMap pmap = get(boost::vertex_point, *p.second);
+
+  std::vector<Point_3> points;
+  halfedge_descriptor hd = p.first;
+
+  // collect points to be add inside the edges
+  sample_edge(hd, *p.second, target_length, pmap, std::back_inserter(points) );
+
+  CGAL_assertion_code(Point_3 src = get(pmap, source(hd, *p.second));)
+  CGAL_assertion_code(Point_3 tgt = get(pmap, target(hd, *p.second));)
+
+  // split the edges and collect faces to triangulate
+  BOOST_FOREACH(const Pair& h_and_p, hedge_range)
+  {
+    halfedge_descriptor hc=h_and_p.first;
+    TriangleMesh* polyc = h_and_p.second;
+
+    PointPMap pmap_2 = get(boost::vertex_point, *polyc);
+    //make sure halfedge are consistently oriented
+    CGAL_assertion( get(pmap_2, source(hc, *polyc)) == src );
+    CGAL_assertion( get(pmap_2, target(hc, *polyc)) == tgt );
+
+    typedef typename Edges_to_protect::value_type::second_type Edge_set;
+    Edge_set& edge_set = edges_to_protect[polyc];
+
+    // now split the halfedge and incident faces
+    split_identical_edges(hc,*polyc,pmap_2, points,
+                          std::inserter( edge_set, edge_set.begin()));
+  }
+}
+
 using namespace CGAL::Three;
 class Polyhedron_demo_isotropic_remeshing_plugin :
   public QObject,
@@ -48,6 +164,13 @@ class Polyhedron_demo_isotropic_remeshing_plugin :
   Q_INTERFACES(CGAL::Three::Polyhedron_demo_plugin_interface)
   Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0")
 
+  typedef boost::graph_traits<Polyhedron>::edge_descriptor edge_descriptor;
+  typedef boost::graph_traits<Polyhedron>::halfedge_descriptor halfedge_descriptor;
+  typedef boost::graph_traits<Polyhedron>::face_descriptor face_descriptor;
+
+  typedef boost::unordered_set<edge_descriptor, CGAL::Handle_hash_function>    Edge_set;
+  typedef Scene_polyhedron_selection_item::Is_constrained_map<Edge_set> Edge_constrained_pmap;
+
 public:
   void init(QMainWindow* mainWindow, Scene_interface* scene_interface)
   {
@@ -83,6 +206,58 @@ public:
     return false;
   }
 
+  void detect_and_split_duplicates(std::vector<Scene_polyhedron_item*>& selection,
+                                   std::map<Polyhedron*,Edge_set>& edges_to_protect,
+                                   double target_length)
+  {
+    typedef Polyhedron::Point_3 Point_3;
+    typedef std::pair<Point_3,Point_3> Segment_3;
+
+    typedef std::map< Segment_3,
+                      std::vector< std::pair<halfedge_descriptor, Polyhedron*> > > MapType;
+    MapType duplicated_edges;
+
+
+    BOOST_FOREACH(Scene_polyhedron_item* poly_item, selection){
+      Polyhedron& pmesh = *poly_item->polyhedron();
+      BOOST_FOREACH(edge_descriptor ed, edges(pmesh)){
+        halfedge_descriptor hd = halfedge(ed,pmesh);
+        Point_3 p = source(hd,pmesh)->point(), q = target(hd,pmesh)->point();
+        Segment_3 s = CGAL::make_sorted_pair(p,q);
+        if (s.first==q) hd=opposite(hd,pmesh); // make sure the halfedges are consistently oriented
+
+        duplicated_edges[s].push_back( std::make_pair(hd,&pmesh) );
+      }
+    }
+
+    // consistently split duplicate edges and triangulate incident faces
+    typedef std::pair<face_descriptor, Polyhedron*> Face_and_poly;
+    std::set< Face_and_poly > faces_to_triangulate;
+    BOOST_FOREACH(const MapType::value_type& p, duplicated_edges)
+      if (p.second.size()>1){
+        //collect faces to retriangulate
+        typedef std::pair<halfedge_descriptor, Polyhedron*> Pair_type;
+        BOOST_FOREACH(const Pair_type& h_and_p, p.second)
+        {
+          halfedge_descriptor hc=h_and_p.first;
+          Polyhedron* polyc = h_and_p.second;
+
+          if ( !is_border(hc, *polyc) )
+            faces_to_triangulate.insert( Face_and_poly(face(hc,*polyc), polyc) );
+          if ( !is_border(opposite(hc, *polyc), *polyc) )
+            faces_to_triangulate.insert(
+              Face_and_poly(face(opposite(hc, *polyc),*polyc), polyc) );
+        }
+        // split the edges
+        split_long_duplicated_edge(p.second, target_length, edges_to_protect);
+      }
+    // now retriangulate
+    namespace PMP=CGAL::Polygon_mesh_processing;
+    BOOST_FOREACH(Face_and_poly f_and_p, faces_to_triangulate)
+      PMP::triangulate_face(f_and_p.first, *f_and_p.second);
+  }
+
+
 public Q_SLOTS:
   void isotropic_remeshing()
   {
@@ -114,6 +289,7 @@ public Q_SLOTS:
         return;
       }
       bool edges_only = ui.splitEdgesOnly_checkbox->isChecked();
+      bool preserve_duplicates = ui.preserveDuplicates_checkbox->isChecked();
       double target_length = ui.edgeLength_dspinbox->value();
       unsigned int nb_iter = ui.nbIterations_spinbox->value();
       bool protect = ui.protect_checkbox->isChecked();
@@ -132,6 +308,15 @@ public Q_SLOTS:
       const Polyhedron& pmesh = (poly_item != NULL)
         ? *poly_item->polyhedron()
         : *selection_item->polyhedron();
+
+      // tricks to use the function detect_and_split_duplicates
+      // that uses several poly items
+      std::map<Polyhedron*,Edge_set > edges_to_protect_map;
+      std::vector<Scene_polyhedron_item*> poly_items(1,poly_item);
+      Edge_set& edges_to_protect=edges_to_protect_map[poly_item->polyhedron()];
+      if(preserve_duplicates)
+        detect_and_split_duplicates(poly_items, edges_to_protect_map, target_length);
+
       boost::property_map<Polyhedron, CGAL::face_index_t>::type fim
         = get(CGAL::face_index, pmesh);
       unsigned int id = 0;
@@ -174,13 +359,16 @@ public Q_SLOTS:
         }
         else
         {
+     Edge_constrained_pmap ecm(edges_to_protect);
+     CGAL::OR_property_map<Edge_constrained_pmap, Edge_constrained_pmap> etp(ecm, selection_item->constrained_edges_pmap());
+
          CGAL::Polygon_mesh_processing::isotropic_remeshing(
            selection_item->selected_facets
          , target_length
          , *selection_item->polyhedron()
          , CGAL::Polygon_mesh_processing::parameters::number_of_iterations(nb_iter)
          .protect_constraints(protect)
-         .edge_is_constrained_map(selection_item->constrained_edges_pmap())
+         .edge_is_constrained_map(etp)
          .smooth_along_features(smooth_features)
          .vertex_is_constrained_map(selection_item->constrained_vertices_pmap()));
         }
@@ -212,12 +400,14 @@ public Q_SLOTS:
         }
         else
         {
+        Scene_polyhedron_selection_item::Is_constrained_map<Edge_set > ecm(edges_to_protect);
         CGAL::Polygon_mesh_processing::isotropic_remeshing(
            faces(*poly_item->polyhedron())
          , target_length
          , *poly_item->polyhedron()
          , CGAL::Polygon_mesh_processing::parameters::number_of_iterations(nb_iter)
          .protect_constraints(protect)
+          .edge_is_constrained_map(ecm)
          .smooth_along_features(smooth_features));
         }
         poly_item->invalidateOpenGLBuffers();
@@ -236,7 +426,7 @@ public Q_SLOTS:
   void isotropic_remeshing_of_several_polyhedra()
   {
     // Remeshing parameters
-    bool edges_only = false;
+    bool edges_only = false, preserve_duplicates = false;
     double target_length = 0.;
     unsigned int nb_iter = 1;
     bool protect = false;
@@ -272,6 +462,7 @@ public Q_SLOTS:
         }
 
         edges_only = ui.splitEdgesOnly_checkbox->isChecked();
+        preserve_duplicates = ui.preserveDuplicates_checkbox->isChecked();
         target_length = ui.edgeLength_dspinbox->value();
         nb_iter = ui.nbIterations_spinbox->value();
         protect = ui.protect_checkbox->isChecked();
@@ -290,6 +481,13 @@ public Q_SLOTS:
     QApplication::setOverrideCursor(Qt::WaitCursor);
     int total_time = 0;
 
+
+    //     typedef boost::graph_traits<Polyhedron>::edge_descriptor edge_descriptor;
+    std::map<Polyhedron*,Edge_set > edges_to_protect;
+
+    if(preserve_duplicates)
+      detect_and_split_duplicates(selection, edges_to_protect, target_length);
+
 #ifdef CGAL_LINKED_WITH_TBB
     QTime time;
     time.start();
@@ -297,10 +495,12 @@ public Q_SLOTS:
       tbb::parallel_for(
         tbb::blocked_range<std::size_t>(0, selection.size()),
         Remesh_polyhedron_item_for_parallel_for<Remesh_polyhedron_item>(
-          selection, edges_only, target_length, nb_iter, protect, smooth_features));
+                                                                        selection, edges_to_protect, edges_only, target_length, nb_iter, protect, smooth_features));
 
     total_time = time.elapsed();
+
 #else
+
     Remesh_polyhedron_item remesher(edges_only,
       target_length, nb_iter, protect, smooth_features);
     BOOST_FOREACH(Scene_polyhedron_item* poly_item, selection)
@@ -308,7 +508,7 @@ public Q_SLOTS:
       QTime time;
       time.start();
 
-        remesher(poly_item);
+      remesher(poly_item, edges_to_protect[poly_item->polyhedron()]);
 
       total_time += time.elapsed();
       std::cout << "Remeshing of " << poly_item->name().data()
@@ -323,7 +523,7 @@ public Q_SLOTS:
       poly_item->invalidateOpenGLBuffers();
       Q_EMIT poly_item->itemChanged();
     }
-    
+
     // default cursor
     QApplication::restoreOverrideCursor();
   }
@@ -342,7 +542,8 @@ private:
     bool smooth_features_;
 
   protected:
-    void remesh(Scene_polyhedron_item* poly_item) const
+    void remesh(Scene_polyhedron_item* poly_item,
+                Edge_set& edges_to_protect) const
     {
       //fill face_index property map
       boost::property_map<Polyhedron, CGAL::face_index_t>::type fim
@@ -369,12 +570,14 @@ private:
       }
       else
       {
+        Scene_polyhedron_selection_item::Is_constrained_map<Edge_set > ecm(edges_to_protect);
         CGAL::Polygon_mesh_processing::isotropic_remeshing(
             faces(*poly_item->polyhedron())
           , target_length_
           , *poly_item->polyhedron()
           , CGAL::Polygon_mesh_processing::parameters::number_of_iterations(nb_iter_)
           .protect_constraints(protect_)
+          .edge_is_constrained_map(ecm)
           .smooth_along_features(smooth_features_));
       }
     }
@@ -401,9 +604,10 @@ private:
       , smooth_features_(remesh.smooth_features_)
     {}
 
-    void operator()(Scene_polyhedron_item* poly_item) const
+    void operator()(Scene_polyhedron_item* poly_item,
+                    Edge_set& edges_to_protect) const
     {
-      remesh(poly_item);
+      remesh(poly_item, edges_to_protect);
     }
   };
 
@@ -413,34 +617,35 @@ private:
     : RemeshFunctor
   {
     const std::vector<Scene_polyhedron_item*>& selection_;
+    std::map<Polyhedron*,Edge_set >& edges_to_protect_;
 
   public:
     // Constructor
     Remesh_polyhedron_item_for_parallel_for(
       const std::vector<Scene_polyhedron_item*>& selection,
+      std::map<Polyhedron*,Edge_set >& edges_to_protect,
       const bool edges_only,
       const double target_length,
       const unsigned int nb_iter,
       const bool protect,
       const bool smooth_features)
       : RemeshFunctor(edges_only, target_length, nb_iter, protect, smooth_features)
-      , selection_(selection)
-    {
-      ;
-    }
+      , selection_(selection), edges_to_protect_(edges_to_protect)
+    {}
 
     // Constructor
     Remesh_polyhedron_item_for_parallel_for(
       const Remesh_polyhedron_item_for_parallel_for &remesh)
       : RemeshFunctor(remesh)
       , selection_(remesh.selection_)
+      , edges_to_protect_(remesh.edges_to_protect_)
     {}
 
     // operator()
     void operator()(const tbb::blocked_range<size_t>& r) const
     {
       for (size_t i = r.begin(); i != r.end(); ++i)
-        RemeshFunctor::remesh(selection_[i]);
+        RemeshFunctor::remesh(selection_[i], edges_to_protect_[selection_[i]->polyhedron()]);
     }
   };
 #endif
@@ -464,6 +669,10 @@ private:
             ui.smooth1D_checkbox, SLOT(setDisabled(bool)));
     connect(ui.splitEdgesOnly_checkbox, SIGNAL(toggled(bool)),
             ui.smooth1D_checkbox, SLOT(setDisabled(bool)));
+    connect(ui.preserveDuplicates_checkbox, SIGNAL(toggled(bool)),
+            ui.protect_checkbox, SLOT(setChecked(bool)));
+    connect(ui.preserveDuplicates_checkbox, SIGNAL(toggled(bool)),
+            ui.protect_checkbox, SLOT(setDisabled(bool)));
 
     //Set default parameters
     Scene_interface::Bbox bbox = poly_item != NULL ? poly_item->bbox()
diff --git a/Property_map/include/CGAL/property_map.h b/Property_map/include/CGAL/property_map.h
index fa41272e98a..a4a8f871169 100644
--- a/Property_map/include/CGAL/property_map.h
+++ b/Property_map/include/CGAL/property_map.h
@@ -38,6 +38,40 @@ namespace CGAL {
 
 
 /// \cond SKIP_DOXYGEN
+
+
+template <typename PM1, typename PM2>
+class OR_property_map {
+  typedef typename PM1::key_type key_type;
+  typedef typename PM1::value_type value_type;
+  typedef typename PM1::reference reference;
+  typedef boost::read_write_property_map_tag category;
+
+  PM1 pm1;
+  PM2 pm2;
+
+ public:
+  OR_property_map(PM1 pm1, PM2 pm2)
+    : pm1(pm1),pm2(pm2)
+  {}
+
+  inline friend
+    value_type
+    get(const OR_property_map& pm, const key_type& k)
+  {
+    return get(pm.pm1,k) || get(pm.pm2,k);
+  }
+
+  inline friend
+    void
+  put(OR_property_map& pm, const key_type& k, const value_type& v)
+  {
+    put(pm.pm1,k, v);
+    put(pm.pm2,k, v);
+  }
+
+};
+
 /// Property map that accesses a value from an iterator
 ///
 /// \cgalModels `ReadablePropertyMap`