From a2d7de113d0f57398e167878892090d00db884fe Mon Sep 17 00:00:00 2001
From: David Wells <wellsd2@rpi.edu>
Date: Sun, 29 Apr 2018 14:08:50 -0400
Subject: [PATCH] Avoid writing colinear points in GridOut::write_gnuplot.

GridOut::write_gnuplot has an option where, for plotting curved grids, it can
plot multiple line segments per face to give the illusion of a cell with curved
faces. This procedure creates massive output files since every face of every
corresponds to many individual line segments.

This commit implements a compression scheme where we do not output intermediate
points that are colinear. This lowers the number of points in a polar grid
output file by about 50% (i.e., the radial cell faces are not curved and can be
plotted with two points instead of perhaps five or six; same for the interior
cell faces).
---
 source/grid/grid_out.cc | 67 ++++++++++++++++++++++++++++++++++++-----
 1 file changed, 59 insertions(+), 8 deletions(-)

diff --git a/source/grid/grid_out.cc b/source/grid/grid_out.cc
index 42b8ad118a..3f09caab8d 100644
--- a/source/grid/grid_out.cc
+++ b/source/grid/grid_out.cc
@@ -3157,6 +3157,34 @@ namespace internal
 {
   namespace
   {
+    /**
+     * GNUPlot output can, optionally, output multiple line segments on each
+     * grid line to make curved lines look curved. However, this is very
+     * wasteful when the line itself is straight since we do not need multiple
+     * line segments to draw a straight line. This function tries to identify
+     * whether or not the collection of points corresponds to a straight line:
+     * if it does then all but the first and last points can be removed.
+     */
+    template <int spacedim>
+    void
+    remove_colinear_points(std::vector<Point<spacedim>> &points)
+    {
+      while (points.size() > 2)
+        {
+          Tensor<1, spacedim> first_difference = points[1] - points[0];
+          first_difference /= first_difference.norm();
+          Tensor<1, spacedim> second_difference = points[2] - points[1];
+          second_difference /= second_difference.norm();
+          // If the three points are colinear then remove the middle one.
+          if ((first_difference - second_difference).norm() < 1e-10)
+            points.erase(points.begin() + 1);
+          else
+            break;
+        }
+    }
+
+
+
     template <int spacedim>
     void write_gnuplot (const dealii::Triangulation<1,spacedim> &tria,
                         std::ostream             &out,
@@ -3263,12 +3291,24 @@ namespace internal
                   face = cell->face(face_no);
                   if (face->at_boundary() || gnuplot_flags.curved_inner_cells)
                     {
+                      // Save the points on each face to a vector and then try
+                      // to remove colinear points that won't show up in the
+                      // generated plot.
+                      std::vector<Point<spacedim>> line_points;
                       // compute offset of quadrature points within set of
                       // projected points
                       const unsigned int offset=face_no*n_points;
-                      for (unsigned int i=0; i<n_points; ++i)
-                        out << (mapping->transform_unit_to_real_cell
-                                (cell, q_projector->point(offset+i)))
+                      // we don't need to transform the vertices: they are
+                      // already in physical space.
+                      line_points.push_back(face->vertex(0));
+                      for (unsigned int i=1; i<n_points - 1; ++i)
+                        line_points.push_back(mapping->transform_unit_to_real_cell
+                                              (cell, q_projector->point(offset+i)));
+                      line_points.push_back(face->vertex(1));
+                      internal::remove_colinear_points(line_points);
+
+                      for (const Point<spacedim> &point : line_points)
+                        out << point
                             << ' ' << cell->level()
                             << ' ' << static_cast<unsigned int>(cell->material_id())
                             << '\n';
@@ -3462,19 +3502,30 @@ namespace internal
                           const typename dealii::Triangulation<dim,spacedim>::line_iterator
                           line=face->line(l);
 
-                          const Point<spacedim> &v0=line->vertex(0),
-                                                 &v1=line->vertex(1);
+                          const Point<spacedim> &v0=line->vertex(0), &v1=line->vertex(1);
                           if (line->at_boundary() || gnuplot_flags.curved_inner_cells)
                             {
+                              // Save the points on each face to a vector and
+                              // then try to remove colinear points that won't
+                              // show up in the generated plot.
+                              std::vector<Point<spacedim>> line_points;
                               // transform_real_to_unit_cell could be replaced
                               // by using QProjector<dim>::project_to_line
                               // which is not yet implemented
                               const Point<spacedim> u0=mapping->transform_real_to_unit_cell(cell, v0),
                                                     u1=mapping->transform_real_to_unit_cell(cell, v1);
 
-                              for (unsigned int i=0; i<n_points; ++i)
-                                out << (mapping->transform_unit_to_real_cell
-                                        (cell, (1-boundary_points[i][0])*u0+boundary_points[i][0]*u1))
+                              const Point<spacedim> center;
+                              // we don't need to transform the vertices: they
+                              // are already in physical space.
+                              line_points.push_back(v0);
+                              for (unsigned int i=1; i<n_points - 1; ++i)
+                                line_points.push_back(mapping->transform_unit_to_real_cell
+                                                      (cell, (1-boundary_points[i][0])*u0 + boundary_points[i][0]*u1));
+                              line_points.push_back(v1);
+                              internal::remove_colinear_points(line_points);
+                              for (const Point<spacedim> &point : line_points)
+                                out << point
                                     << ' ' << cell->level()
                                     << ' ' << static_cast<unsigned int>(cell->material_id()) << '\n';
                             }
-- 
2.39.5