From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Thu, 14 Dec 2017 06:46:19 +0000 (+0100)
Subject: Fix corner case in comparison. Add test.
X-Git-Tag: v9.0.0-rc1~659^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2696c47d839e9cdc28b7df63f8e6be9b8efbca36;p=dealii.git

Fix corner case in comparison. Add test.
---

diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc
index 727a99eb4c..17d51850a3 100644
--- a/source/dofs/dof_tools.cc
+++ b/source/dofs/dof_tools.cc
@@ -60,9 +60,13 @@ namespace DoFTools
      * std::map.
      *
      * Comparison is done through an artificial downstream direction that
-     * tells directions apart through a factor of 1e-5. The final comparison
-     * is done without an epsilon, so points need to have identical floating
-     * point components to be considered equal.
+     * tells directions apart through a factor of 1e-5. Once we got the
+     * direction, we check for its value. In case the distance is exactly zero
+     * (without an epsilon), we might still have the case that two points
+     * combine in a particular way, e.g. the points (1.0, 1.0) and (1.0+1e-5,
+     * 0.0). Thus, compare the points component by component in the second
+     * step. Thus, points need to have identical floating point components to
+     * be considered equal.
      */
     template <int dim,typename Number=double>
     struct ComparisonHelper
@@ -79,10 +83,23 @@ namespace DoFTools
         double weight = 1.;
         for (unsigned int d=0; d<dim; ++d)
           {
-            downstream_size += (lhs[d] - rhs[d]) * weight;
+            downstream_size += (rhs[d] - lhs[d]) * weight;
             weight *= 1e-5;
           }
-        return downstream_size < 0;
+        if (downstream_size < 0)
+          return false;
+        else if (downstream_size > 0)
+          return true;
+        else
+          {
+            for (unsigned int d=0; d<dim; ++d)
+              {
+                if (lhs[d] == rhs[d])
+                  continue;
+                return lhs[d]<rhs[d];
+              }
+            return false;
+          }
       }
 
     };
diff --git a/tests/dofs/dof_tools_write_gnuplot_dof_support_point_info_04.cc b/tests/dofs/dof_tools_write_gnuplot_dof_support_point_info_04.cc
new file mode 100644
index 0000000000..6250b02629
--- /dev/null
+++ b/tests/dofs/dof_tools_write_gnuplot_dof_support_point_info_04.cc
@@ -0,0 +1,46 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Check DoFTools::write_gnuplot_dof_support_point_info for a particular set of points
+
+#include "../tests.h"
+
+#include <deal.II/dofs/dof_tools.h>
+
+void
+test (const double epsilon)
+{
+  constexpr unsigned int dim = 2;
+  Point<dim> p1(1.0, 1.0);
+  Point<dim> p2(1.0+epsilon, 0.0);
+
+  std::map<types::global_dof_index, Point<dim> > support_points;
+  support_points[0] = p1;
+  support_points[1] = p2;
+
+  DoFTools::write_gnuplot_dof_support_point_info(deallog.get_file_stream(),
+                                                 support_points);
+  deallog << std::endl;
+}
+
+
+
+int main()
+{
+  initlog();
+  test(1e-4);
+  test(1e-5);
+  test(1e-6);
+}
diff --git a/tests/dofs/dof_tools_write_gnuplot_dof_support_point_info_04.output b/tests/dofs/dof_tools_write_gnuplot_dof_support_point_info_04.output
new file mode 100644
index 0000000000..c591160267
--- /dev/null
+++ b/tests/dofs/dof_tools_write_gnuplot_dof_support_point_info_04.output
@@ -0,0 +1,10 @@
+
+1.00000 1.00000 "0"
+1.00010 0.00000 "1"
+DEAL::
+1.00000 1.00000 "0"
+1.00001 0.00000 "1"
+DEAL::
+1.00000 0.00000 "1"
+1.00000 1.00000 "0"
+DEAL::