diff --git a/BGL/examples/BGL_surface_mesh/surface_mesh_partition.cpp b/BGL/examples/BGL_surface_mesh/surface_mesh_partition.cpp
index 3fa6921f638..629fddfb640 100644
--- a/BGL/examples/BGL_surface_mesh/surface_mesh_partition.cpp
+++ b/BGL/examples/BGL_surface_mesh/surface_mesh_partition.cpp
@@ -44,9 +44,7 @@ int main(int argc, char** argv)
   CGAL::copy_face_graph(filtered_sm, part_sm);
 
   // Output the mesh extracted from subpart n°0
-  std::ofstream out("sm_part_0.off");
-  out.precision(17);
-  CGAL::write_OFF(out, part_sm);
+  CGAL::write_polygon_mesh("sm_part_0.off", part_sm, CGAL::parameters::stream_precision(17));
 
   // Output all the vertices that are in the part n°0
   std::ofstream outxyz("out.xyz");
diff --git a/Point_set_3/examples/Point_set_3/point_set_read_ply.cpp b/Point_set_3/examples/Point_set_3/point_set_read_ply.cpp
index 9b60f42c2fa..5a58de114cb 100644
--- a/Point_set_3/examples/Point_set_3/point_set_read_ply.cpp
+++ b/Point_set_3/examples/Point_set_3/point_set_read_ply.cpp
@@ -14,12 +14,12 @@ typedef CGAL::Point_set_3<Point> Point_set;
 
 int main (int argc, char** argv)
 {
-  std::ifstream f (argc > 1 ? argv[1] : "data/example.ply",
-                   std::ios_base::binary); // Mandatory on Windows if input is binary PLY
+  std::ifstream f(argc > 1 ? argv[1] : "data/example.ply",
+                  std::ios_base::binary); // Mandatory on Windows if input is binary PLY
 
   Point_set point_set;
 
-  if(!CGAL::read_PLY (f, point_set)) // same as `f >> point_set`
+  if(!CGAL::read_PLY(f, point_set)) // same as `f >> point_set`
   {
     std::cerr << "Can't read input file " << std::endl;
     return EXIT_FAILURE;
@@ -28,33 +28,29 @@ int main (int argc, char** argv)
   // Shows which properties are defined
   std::vector<std::string> properties = point_set.properties();
   std::cerr << "Properties:" << std::endl;
-  for (std::size_t i = 0; i < properties.size(); ++ i)
+  for(std::size_t i = 0; i < properties.size(); ++ i)
     std::cerr << " * " << properties[i] << std::endl;
 
   // Recover "label" property of type int
   Point_set::Property_map<boost::int32_t> label_prop;
   bool found = false;
-  boost::tie (label_prop, found)  = point_set.property_map<boost::int32_t> ("label");
+  boost::tie(label_prop, found)  = point_set.property_map<boost::int32_t> ("label");
 
-  if (found)
-    {
-      std::cerr << "Point set has an integer \"label\" property with values:" << std::endl;
-      for (Point_set::iterator it = point_set.begin(); it != point_set.end(); ++ it)
-        std::cerr << " * " << label_prop[*it] << std::endl;
-    }
-
-  if (argc > 2 && strcmp (argv[2], "-b") == 0) // Optional binary output
+  if(found)
   {
-    std::ofstream out ("out.ply",
-                       std::ios_base::binary); // Mandatory on Windows
-    CGAL::set_binary_mode (out); // Select binary mode (ASCII is default)
-    out << point_set; // same as `CGAL::write_PLY(out, point_set)`
+    std::cerr << "Point set has an integer \"label\" property with values:" << std::endl;
+    for (Point_set::iterator it = point_set.begin(); it != point_set.end(); ++ it)
+      std::cerr << " * " << label_prop[*it] << std::endl;
+  }
+
+  if(argc > 2 && strcmp (argv[2], "-b") == 0) // Optional binary output
+  {
+    CGAL::write_PLY("out.ply", point_set, CGAL::parameters::stream_precision(17));
   }
   else // ASCII output
   {
-    std::ofstream out ("out.ply");
-    out.precision(17); // Use sufficient precision in ASCII
-    CGAL::write_PLY (out, point_set); // same as `out << point_set`
+    CGAL::write_PLY("out.ply", point_set, CGAL::parameters::stream_precision(17)
+                                                           .use_binary_mode(false));
   }
 
   return 0;
diff --git a/Point_set_processing_3/examples/Point_set_processing_3/bilateral_smooth_point_set_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/bilateral_smooth_point_set_example.cpp
index 74e24fd4269..dc37c7a7eb8 100644
--- a/Point_set_processing_3/examples/Point_set_processing_3/bilateral_smooth_point_set_example.cpp
+++ b/Point_set_processing_3/examples/Point_set_processing_3/bilateral_smooth_point_set_example.cpp
@@ -50,7 +50,7 @@ int main(int argc, char*argv[])
       k,
       CGAL::parameters::point_map(CGAL::First_of_pair_property_map<PointVectorPair>())
                        .normal_map(CGAL::Second_of_pair_property_map<PointVectorPair>())
-                       .sharpness_angle (sharpness_angle));
+                       .sharpness_angle(sharpness_angle));
   }
 
   //// Save point set.
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_SM.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_SM.cpp
index f9de85eef70..e4660e6fdef 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_SM.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_SM.cpp
@@ -34,12 +34,8 @@ int main(int argc, char* argv[])
             << num_vertices(mesh1) << " and "
             << num_vertices(mesh2) << "\n";
 
-  std::ofstream output("mesh1_refined.off");
-  output.precision(17);
-  output << mesh1;
-  output.close();
-  output.open("mesh2_refined.off");
-  output << mesh2;
+  CGAL::write_polygon_mesh("mesh1_refined.off", mesh1, CGAL::parameters::stream_precision(17));
+  CGAL::write_polygon_mesh("mesh2_refined.off", mesh2, CGAL::parameters::stream_precision(17));
 
   return 0;
 }
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_consecutive_bool_op.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_consecutive_bool_op.cpp
index 0338230f3fa..de555e8e25a 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_consecutive_bool_op.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_consecutive_bool_op.cpp
@@ -102,9 +102,7 @@ int main(int argc, char* argv[])
                                          params::vertex_point_map(mesh2_vpm) ) )
     {
       std::cout << "Intersection and union were successfully computed\n";
-      std::ofstream output("inter_union.off");
-      output.precision(17);
-      output << mesh2;
+      CGAL::write_polygon_mesh("inter_union.off", mesh2, CGAL::parameters::stream_precision(17));
       return 0;
     }
 
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_difference_remeshed.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_difference_remeshed.cpp
index 9b67d691b4f..54df02b8707 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_difference_remeshed.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_difference_remeshed.cpp
@@ -63,11 +63,10 @@ int main(int argc, char* argv[])
   if (valid_difference)
   {
     std::cout << "Difference was successfully computed\n";
-    std::ofstream output("difference.off");
-    output.precision(17);
-    output << mesh1;
+    CGAL::write_polygon_mesh("difference.off", mesh1, CGAL::parameters::stream_precision(17));
   }
-  else{
+  else
+  {
     std::cout << "Difference could not be computed\n";
     return 1;
   }
@@ -79,8 +78,7 @@ int main(int argc, char* argv[])
     if (is_constrained_map[e])
     {
       // insert all faces incident to the target vertex
-      for(halfedge_descriptor h :
-                    halfedges_around_target(halfedge(e,mesh1),mesh1))
+      for(halfedge_descriptor h : halfedges_around_target(halfedge(e,mesh1),mesh1))
       {
         if (!is_border(h, mesh1) )
         {
@@ -105,9 +103,7 @@ int main(int argc, char* argv[])
   PMP::isotropic_remeshing(selected_faces, 0.02, mesh1,
                            params::edge_is_constrained_map(is_constrained_map));
 
-  std::ofstream output("difference_remeshed.off");
-  output.precision(17);
-  output << mesh1;
+  CGAL::write_polygon_mesh("difference_remeshed.off", mesh1, CGAL::parameters::stream_precision(17));
 
   return 0;
 }
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union.cpp
index d4f92394a72..be53bb63b82 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union.cpp
@@ -29,9 +29,7 @@ int main(int argc, char* argv[])
   if(valid_union)
   {
     std::cout << "Union was successfully computed\n";
-    std::ofstream output("union.off");
-    output.precision(17);
-    output << out;
+    CGAL::write_polygon_mesh("union.off", out, CGAL::parameters::stream_precision(17));
     return 0;
   }
 
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union_and_intersection.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union_and_intersection.cpp
index 52fe89de94b..2ea59a2ac90 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union_and_intersection.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union_and_intersection.cpp
@@ -48,9 +48,7 @@ int main(int argc, char* argv[])
   if (res[PMP::Corefinement::UNION])
   {
     std::cout << "Union was successfully computed\n";
-    std::ofstream output("union.off");
-    output.precision(17);
-    output << out_union;
+    CGAL::write_polygon_mesh("union.off", out_union, CGAL::parameters::stream_precision(17));
   }
   else
     std::cout << "Union could not be computed\n";
@@ -58,9 +56,7 @@ int main(int argc, char* argv[])
   if (res[PMP::Corefinement::INTERSECTION])
   {
     std::cout << "Intersection was successfully computed\n";
-    std::ofstream output("intersection.off");
-    output.precision(17);
-    output << out_intersection;
+    CGAL::write_polygon_mesh("intersection.off", out_intersection, CGAL::parameters::stream_precision(17));
   }
   else
     std::cout << "Intersection could not be computed\n";
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union_with_attributes.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union_with_attributes.cpp
index b616dfd3948..787844be103 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union_with_attributes.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/corefinement_mesh_union_with_attributes.cpp
@@ -88,10 +88,7 @@ int main(int argc, char* argv[])
 
   if (valid_union)
   {
-    std::cout << "Union was successfully computed\n";
-    std::ofstream output("union.off");
-    output.precision(17);
-    output << out;
+    CGAL::write_polygon_mesh("union.off", out, CGAL::parameters::stream_precision(17));
     return 0;
   }
 
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hole_filling_example_SM.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hole_filling_example_SM.cpp
index e79778f984d..40689ff2054 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/hole_filling_example_SM.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/hole_filling_example_SM.cpp
@@ -90,11 +90,8 @@ int main(int argc, char* argv[])
   std::cout << std::endl;
   std::cout << nb_holes << " holes have been filled" << std::endl;
 
-  std::string outfile = "filled_SM.off";
-  std::ofstream out(outfile.c_str());
-  std::cout << "Mesh written to: " << outfile << std::endl;
-  out.precision(17);
-  out << mesh << std::endl;
+  CGAL::write_polygon_mesh("filled_SM.off", mesh, CGAL::parameters::stream_precision(17));
+  std::cout << "Mesh written to: filled_SM.off" << std::endl;
 
   return 0;
 }
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_of_patch_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_of_patch_example.cpp
index c61fbfa8904..b3a1395b2f7 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_of_patch_example.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/isotropic_remeshing_of_patch_example.cpp
@@ -75,9 +75,7 @@ int main(int argc, char* argv[])
                            .protect_constraints(true)//i.e. protect border, here
                            );
 
-  std::ofstream out("out.off");
-  out.precision(17);
-  out << mesh;
+  CGAL::write_polygon_mesh("out.off", mesh, CGAL::parameters::stream_precision(17));
   std::cout << "Remeshing done." << std::endl;
 
   return 0;
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_smoothing_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_smoothing_example.cpp
index 21046c672a7..afb47b3a465 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_smoothing_example.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/mesh_smoothing_example.cpp
@@ -47,9 +47,7 @@ int main(int argc, char** argv)
                                          .use_safety_constraints(false) // authorize all moves
                                          .edge_is_constrained_map(eif));
 
-  std::ofstream output("mesh_smoothed.off");
-  output.precision(17);
-  output << mesh;
+  CGAL::write_polygon_mesh("mesh_smoothed.off", mesh, CGAL::parameters::stream_precision(17));
 
   std::cout << "Done!" << std::endl;
 
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/random_perturbation_SM_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/random_perturbation_SM_example.cpp
index 70004cc1d79..d68467b62bf 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/random_perturbation_SM_example.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/random_perturbation_SM_example.cpp
@@ -31,10 +31,7 @@ int main(int argc, char* argv[])
 
   PMP::random_perturbation(mesh, max_size, CGAL::parameters::vertex_point_map(mesh.points()).geom_traits(K()));
 
-  std::ofstream out("data/eight_perturbed.off");
-  out.precision(17);
-  out << mesh;
-  out.close();
+  CGAL::write_polygon_mesh("data/eight_perturbed.off", mesh, CGAL::parameters::stream_precision(17));
 
   return 0;
 }
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/refine_fair_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/refine_fair_example.cpp
index 188c0968fa9..b850c28ae6e 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/refine_fair_example.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/refine_fair_example.cpp
@@ -1,5 +1,4 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-
 #include <CGAL/Polyhedron_3.h>
 
 #include <CGAL/Polygon_mesh_processing/refine.h>
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/shape_smoothing_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/shape_smoothing_example.cpp
index 850b6b1dd00..ae7ba6f4926 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/shape_smoothing_example.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/shape_smoothing_example.cpp
@@ -40,9 +40,7 @@ int main(int argc, char* argv[])
   PMP::smooth_shape(mesh, time, PMP::parameters::number_of_iterations(nb_iterations)
                                                 .vertex_is_constrained_map(vcmap));
 
-  std::ofstream output("mesh_shape_smoothed.off");
-  output.precision(17);
-  output << mesh;
+  CGAL::write_polygon_mesh("mesh_shape_smoothed.off", mesh, CGAL::parameters::stream_precision(17));
 
   std::cout << "Done!" << std::endl;
   return EXIT_SUCCESS;
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example.cpp
index 1106c1ed5a2..149bf6a5fee 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/stitch_borders_example.cpp
@@ -36,9 +36,7 @@ int main(int argc, char* argv[])
   std::cout << "\t Number of halfedges :\t" << mesh.size_of_halfedges() << std::endl;
   std::cout << "\t Number of facets    :\t" << mesh.size_of_facets() << std::endl;
 
-  std::ofstream output("mesh_stitched.off");
-  output.precision(17);
-  output << std::setprecision(17) << mesh;
+  CGAL::write_polygon_mesh("mesh_stitched.off", mesh, CGAL::parameters::stream_precision(17));
 
   return 0;
 }
diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_example.cpp
index c73bb39c787..70c6db3c96f 100644
--- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_example.cpp
+++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/triangulate_faces_example.cpp
@@ -34,9 +34,7 @@ int main(int argc, char* argv[])
     if(!CGAL::is_triangle(halfedge(f, mesh), mesh))
       std::cerr << "Error: non-triangular face left in mesh." << std::endl;
 
-  std::ofstream cube_off(outfilename);
-  cube_off.precision(17);
-  cube_off << mesh;
+  CGAL::write_polygon_mesh(outfilename, mesh, CGAL::parameters::stream_precision(17));
 
   return 0;
 }
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_degeneracies.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_degeneracies.h
index 62958206773..99fa07969c9 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_degeneracies.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_degeneracies.h
@@ -419,10 +419,7 @@ bool remove_almost_degenerate_faces(const FaceRange& face_range,
     std::cout << edges_to_collapse.size() << " needles and " << edges_to_flip.size() << " caps" << std::endl;
     std::ostringstream oss;
     oss << "degen_cleaning_iter_" << iter++ << ".off";
-    std::ofstream out(oss.str().c_str());
-    out << std::setprecision(17);
-    out << tmesh;
-    out.close();
+    CGAL::write_polygon_mesh(oss.str(), tmesh, CGAL::parameters::stream_precision(17));
 #endif
 
     if(edges_to_collapse.empty() && edges_to_flip.empty())
@@ -1638,9 +1635,7 @@ bool remove_degenerate_faces(const FaceRange& face_range,
 #ifdef CGAL_PMP_REMOVE_DEGENERATE_FACES_DEBUG
   {
     std::cout <<"Done with null edges.\n";
-    std::ofstream output("/tmp/no_null_edges.off");
-    output << std::setprecision(17) << tmesh << "\n";
-    output.close();
+    CGAL::write_polygon_mesh("/tmp/no_null_edges.off", tmesh, CGAL::parameters::stream_precision(17));
   }
 #endif
 
diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h
index def4359d132..5dbfa947364 100644
--- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h
+++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/repair_self_intersections.h
@@ -240,8 +240,7 @@ FaceOutputIterator replace_faces_with_patch(const std::vector<typename boost::gr
     remove_face(f, pmesh);
 
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT
-  std::ofstream res_out("results/last_patch_replacement.off");
-  res_out << std::setprecision(17) << pmesh;
+  CGAL::write_polygon_mesh("results/last_patch_replacement.off", pmesh, CGAL::parameters::stream_precision(17));
 #endif
 
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_DEBUG
@@ -463,9 +462,7 @@ bool remove_self_intersections_with_smoothing(std::set<typename boost::graph_tra
   CGAL::copy_face_graph(ffg, local_mesh, CP::vertex_point_map(vpm));
 
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT
-  std::ofstream out_p("results/local_mesh.off");
-  out_p << std::setprecision(17) << local_mesh;
-  out_p.close();
+  CGAL::write_polygon_mesh("results/local_mesh.off", local_mesh, CGAL::parameters::stream_precision(17));
 #endif
 
   // Constrain sharp and border edges
@@ -484,9 +481,7 @@ bool remove_self_intersections_with_smoothing(std::set<typename boost::graph_tra
                                                                          .use_safety_constraints(false));
 
 #ifdef CGAL_PMP_REMOVE_SELF_INTERSECTION_OUTPUT
-  std::ofstream out("results/post_smoothing_local_mesh.off");
-  out << std::setprecision(17) << local_mesh;
-  out.close();
+  CGAL::write_polygon_mesh("results/post_smoothing_local_mesh.off", local_mesh, CGAL::parameters::stream_precision(17));
 #endif
 
   if(does_self_intersect(local_mesh))
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/autorefinement_sm.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/autorefinement_sm.cpp
index e1112893d98..97f80e530ee 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/autorefinement_sm.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/autorefinement_sm.cpp
@@ -3,6 +3,7 @@
 
 #include <CGAL/Polygon_mesh_processing/intersection.h>
 #include <CGAL/Polygon_mesh_processing/corefinement.h>
+#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
 
 #include <fstream>
 
@@ -14,15 +15,13 @@ namespace PMP = CGAL::Polygon_mesh_processing;
 int main(int argc, char* argv[])
 {
   const char* filename = (argc > 1) ? argv[1] : "data/blobby.off";
-  std::ifstream input(filename);
 
   Mesh mesh;
-  if (!input || !(input >> mesh))
+  if(!PMP::read_polygon_mesh(filename, mesh))
   {
     std::cerr << "Input mesh is not a valid off file." << std::endl;
     return 1;
   }
-  input.close();
 
   std::cout << "Test surface_self_intersection\n";
   std::vector< std::vector<K::Point_3> >polylines;
@@ -43,21 +42,17 @@ int main(int argc, char* argv[])
   PMP::experimental::autorefine(mesh);
   std::cout << "Number of vertices after autorefinement " << mesh.number_of_vertices() << "\n";
 
-  output.open("mesh_autorefined.off");
-  output << mesh;
-  output.close();
+  CGAL::write_polygon_mesh("mesh_autorefined.off", mesh, CGAL::parameters::stream_precision(17));
+
+  clear(mesh);
+  CGAL::Polygon_mesh_processing::read_polygon_mesh(filename, mesh);
 
-  input.open(filename);
-  mesh.clear();
-  input >> mesh;
   std::cout << "Number of vertices before self-intersection removal " << mesh.number_of_vertices() << "\n";
   if (!PMP::experimental::autorefine_and_remove_self_intersections(mesh))
     std::cout << "WARNING: Cannot remove all self-intersections\n";
   std::cout << "Number of vertices after self-intersection removal " << mesh.number_of_vertices() << "\n";
 
-  output.open("mesh_fixed.off");
-  output << std::setprecision(17) << mesh;
-  output.close();
+  CGAL::write_polygon_mesh("mesh_fixed.off", mesh, CGAL::parameters::stream_precision(17));
 
-  return 0;
+  return EXIT_SUCCESS;
 }
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_corefinement_bool_op.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_corefinement_bool_op.cpp
index 646aedc9d4d..4cc5233a1ce 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_corefinement_bool_op.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_corefinement_bool_op.cpp
@@ -64,12 +64,11 @@ void run_boolean_operations(
   if ( res[CFR::UNION] ){
    std::cout << "  Union is a valid operation\n";
    assert(union_.is_valid());
-   #ifdef CGAL_COREFINEMENT_DEBUG
+#ifdef CGAL_COREFINEMENT_DEBUG
    std::stringstream fname;
    fname << scenario << "_tm1_union_tm2.off";
-   std::ofstream output(fname.str().c_str());
-   output << std::setprecision(17) << union_;
-   #endif
+   CGAL::write_polygon_mesh(fname.str(), union_, CGAL::parameters::stream_precision(17));
+#endif
   }
   else
     std::cout << "  Union is invalid\n";
@@ -77,12 +76,11 @@ void run_boolean_operations(
   if ( res[CFR::INTERSECTION] ){
    std::cout << "  Intersection is a valid operation\n";
    assert(inter.is_valid());
-   #ifdef CGAL_COREFINEMENT_DEBUG
+#ifdef CGAL_COREFINEMENT_DEBUG
    std::stringstream fname;
    fname << scenario << "_tm1_inter_tm2.off";
-   std::ofstream output(fname.str().c_str());
-   output << std::setprecision(17) << inter;
-   #endif
+   CGAL::write_polygon_mesh(fname.str(), inter, CGAL::parameters::stream_precision(17));
+#endif
   }
   else
     std::cout << "  Intersection is invalid\n";
@@ -90,12 +88,11 @@ void run_boolean_operations(
   if ( res[CFR::TM1_MINUS_TM2] ){
    std::cout << "  tm1-tm2 is a valid operation\n";
    assert(tm1_minus_tm2.is_valid());
-   #ifdef CGAL_COREFINEMENT_DEBUG
+#ifdef CGAL_COREFINEMENT_DEBUG
    std::stringstream fname;
    fname << scenario << "_tm1_minus_tm2.off";
-   std::ofstream output(fname.str().c_str());
-   output << std::setprecision(17) << tm1_minus_tm2;
-   #endif
+   CGAL::write_polygon_mesh(fname.str(), tm1_minus_tm2, CGAL::parameters::stream_precision(17));
+#endif
   }
   else
     std::cout << "  tm1-tm2 is invalid\n";
@@ -103,12 +100,11 @@ void run_boolean_operations(
   if ( res[CFR::TM2_MINUS_TM1] ){
    std::cout << "  tm2-tm1 is a valid operation\n";
    assert(tm2_minus_tm1.is_valid());
-   #ifdef CGAL_COREFINEMENT_DEBUG
+#ifdef CGAL_COREFINEMENT_DEBUG
    std::stringstream fname;
    fname << scenario << "_tm2_minus_tm1.off";
-   std::ofstream output(fname.str().c_str());
-   output << std::setprecision(17) << tm2_minus_tm1;
-   #endif
+   CGAL::write_polygon_mesh(fname.str(), tm2_minus_tm1, CGAL::parameters::stream_precision(17));
+#endif
   }
   else
     std::cout << "  tm2-tm1 is invalid\n";
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_merging_border_vertices.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_merging_border_vertices.cpp
index f1127d50665..a178117ab66 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_merging_border_vertices.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_merging_border_vertices.cpp
@@ -32,9 +32,7 @@ void test_merge_duplicated_vertices_in_boundary_cycles(const char* fname,
   if (expected_nb_vertices==0)
   {
     std::cout << "writing output to out1.off\n";
-    std::ofstream output("out1.off");
-    output << std::setprecision(17);
-    output << mesh;
+    CGAL::write_polygon_mesh("out1.off", mesh, CGAL::parameters::stream_precision(17));
   }
 }
 
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_remove_caps_needles.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_remove_caps_needles.cpp
index 397f4ffd7e4..6e28af873ee 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_remove_caps_needles.cpp
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_remove_caps_needles.cpp
@@ -36,8 +36,9 @@ int main(int argc, char** argv)
                                                     std::cos(160. / 180 * CGAL_PI),
                                                     4,
                                                     0.14);
-  std::ofstream out("cleaned_mesh.off");
-  out << std::setprecision(17) << mesh;
+
+
+  CGAL::write_polygon_mesh("cleaned_mesh.off", mesh, CGAL::parameters::stream_precision(17));
 
   std::cout << "Output mesh has " << edges(mesh).size() << " edges\n";
   if (PMP::does_self_intersect(mesh))
diff --git a/Surface_mesh/examples/Surface_mesh/sm_join.cpp b/Surface_mesh/examples/Surface_mesh/sm_join.cpp
index 62fdf2f369d..187c816e142 100644
--- a/Surface_mesh/examples/Surface_mesh/sm_join.cpp
+++ b/Surface_mesh/examples/Surface_mesh/sm_join.cpp
@@ -33,5 +33,5 @@ int main(int argc, char* argv[])
   for(vertex_descriptor vd : vertices(sm1))
     std::cerr << vd << " " <<  name1[vd] <<std::endl;
 
-  std::cout << std::setprecision(17) << sm1 << std::endl;
+  CGAL::write_OFF(std::cout, sm1, CGAL::parameters::stream_precision(17));
 }
diff --git a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_bounded_normal_change.cpp b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_bounded_normal_change.cpp
index 51734949d85..6d6346e380a 100644
--- a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_bounded_normal_change.cpp
+++ b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_bounded_normal_change.cpp
@@ -1,5 +1,3 @@
-#include <iostream>
-#include <fstream>
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Surface_mesh.h>
 
@@ -9,6 +7,9 @@
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_placement.h>
 #include <CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Bounded_normal_change_placement.h>
 
+#include <iostream>
+#include <fstream>
+
 typedef CGAL::Simple_cartesian<double> Kernel;
 typedef CGAL::Surface_mesh<Kernel::Point_3> Surface_mesh;
 
@@ -48,9 +49,7 @@ int main(int argc, char** argv)
                      CGAL::parameters::get_cost(SMS::LindstromTurk_cost<Surface_mesh>())
                                       .get_placement(Placement()));
 
-  std::ofstream os((argc > 3) ? argv[3] : "out.off");
-  os.precision(17);
-  os << surface_mesh;
+  CGAL::write_polygon_mesh((argc > 3) ? argv[3] : "out.off", surface_mesh, CGAL::parameters::stream_precision(17));
 
   return EXIT_SUCCESS;
 }
diff --git a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_constrained_border_surface_mesh.cpp b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_constrained_border_surface_mesh.cpp
index 92e8184e96c..eaba77051d5 100644
--- a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_constrained_border_surface_mesh.cpp
+++ b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_constrained_border_surface_mesh.cpp
@@ -90,9 +90,7 @@ int main(int argc, char** argv)
   std::cout << "\nFinished!\n" << r << " edges removed.\n"
             << surface_mesh.number_of_edges() << " final edges.\n";
 
-  std::ofstream os(argc > 2 ? argv[2] : "out.off");
-  os.precision(17);
-  os << surface_mesh;
+  CGAL::write_polygon_mesh((argc > 2) ? argv[2] : "out.off", surface_mesh, CGAL::parameters::stream_precision(17));
 
   // now check!
   for(halfedge_descriptor hd : halfedges(surface_mesh))
diff --git a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_garland_heckbert.cpp b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_garland_heckbert.cpp
index c8c5390664f..338a17aed37 100644
--- a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_garland_heckbert.cpp
+++ b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_garland_heckbert.cpp
@@ -68,9 +68,7 @@ int main(int argc, char** argv)
 
   std::cout << "\nFinished!\n" << r << " edges removed.\n" << surface_mesh.number_of_edges() << " final edges.\n";
 
-  std::ofstream os(argc > 3 ? argv[3] : "out.off");
-  os.precision(17);
-  os << surface_mesh;
+  CGAL::write_polygon_mesh((argc > 3) ? argv[3] : "out.off", surface_mesh, CGAL::parameters::stream_precision(17));
 
   return EXIT_SUCCESS;
 }
diff --git a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_surface_mesh.cpp b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_surface_mesh.cpp
index 3ffede53c30..063e6fd63c1 100644
--- a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_surface_mesh.cpp
+++ b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_surface_mesh.cpp
@@ -45,9 +45,7 @@ int main(int argc, char** argv)
   std::cout << "\nFinished!\n" << r << " edges removed.\n" << surface_mesh.number_of_edges() << " final edges.\n";
   std::cout << "Time elapsed: " << std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time).count() << "ms" << std::endl;
 
-  std::ofstream os(argc > 3 ? argv[3] : "out.off");
-  os.precision(17);
-  os << surface_mesh;
+  CGAL::write_polygon_mesh((argc > 3) ? argv[3] : "out.off", surface_mesh, CGAL::parameters::stream_precision(17));
 
   return EXIT_SUCCESS;
 }
diff --git a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_visitor_surface_mesh.cpp b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_visitor_surface_mesh.cpp
index a2f5eb0f98d..ffbe9574fe2 100644
--- a/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_visitor_surface_mesh.cpp
+++ b/Surface_mesh_simplification/examples/Surface_mesh_simplification/edge_collapse_visitor_surface_mesh.cpp
@@ -136,9 +136,7 @@ int main(int argc, char** argv)
   std::cout << "\nFinished!\n" << r << " edges removed.\n"
             << surface_mesh.number_of_edges() << " final edges.\n";
 
-  std::ofstream os(argc > 3 ? argv[3] : "out.off");
-  os.precision(17);
-  os << surface_mesh;
+  CGAL::write_polygon_mesh((argc > 3) ? argv[3] : "out.off", surface_mesh, CGAL::parameters::stream_precision(17));
 
   return EXIT_SUCCESS;
 }