Make sure that compute refinement thresholds are no larger than the max indicator.

This could happen because we search the threshold in an interval larger than the
actual indicators. If we end up with a threshold that is larger than the largest
indicator, we need to cut things back.

In the process of debugging, I also found a few oddities that have to do
with the fact that we do a lot of arithmetic in double precision, but often
pass in criteria as vectors of float. Fix some of those as well.

diff --git a/source/distributed/grid_refinement.cc b/source/distributed/grid_refinement.cc
index 714f25596276e483a3293b9278f1fb89da546db1..f4ab17334534a1c433f9e98d4ff99652aa000b9a 100644 (file)
--- a/source/distributed/grid_refinement.cc
+++ b/source/distributed/grid_refinement.cc
@@ -34,6 +34,7 @@
 #include <numeric>
 #include <algorithm>
 #include <limits>
+#include <iomanip>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -76,7 +77,7 @@ namespace
    * others get a pair of zeros.
    */
   template <typename number>
-  std::pair<double,double>
+  std::pair<number,number>
   compute_global_min_and_max_at_root (const Vector<number> &criteria,
                                       MPI_Comm              mpi_communicator)
   {
@@ -153,7 +154,7 @@ namespace
   void
   get_locally_owned_indicators (const parallel::distributed::Triangulation<dim,spacedim> &tria,
                                 const Vector &criteria,
-                                dealii::Vector<float> &locally_owned_indicators)
+                                dealii::Vector<typename Vector::value_type> &locally_owned_indicators)
   {
     Assert (locally_owned_indicators.size() == tria.n_locally_owned_active_cells(),
             ExcInternalError());
@@ -412,7 +413,22 @@ namespace
                      0, mpi_communicator);
 
           if (interesting_range[0] == interesting_range[1])
-            return interesting_range[0];
+            {
+              // so we have found our threshold. since we adjust
+              // the range at the top of the function to be slightly
+              // larger than the actual extremes of the refinement
+              // criteria values, we can end up in a situation where
+              // the threshold is in fact larger than the maximal
+              // refinement indicator. in such cases, we get no
+              // refinement at all. thus, cap the threshold by the
+              // actual largest value
+              double final_threshold =  std::min (interesting_range[0],
+                                                  global_min_and_max.second);
+              MPI_Bcast (&final_threshold, 1, MPI_DOUBLE,
+                         0, mpi_communicator);
+
+              return final_threshold;
+            }
 
           const double test_threshold
             = (interesting_range[0] > 0
@@ -500,7 +516,15 @@ namespace
                      0, mpi_communicator);
 
           if (interesting_range[0] == interesting_range[1])
-            return interesting_range[0];
+            {
+              // we got the final value. the master adjusts the
+              // value one more time, so receive it from there
+              double final_threshold = -1;
+              MPI_Bcast (&final_threshold, 1, MPI_DOUBLE,
+                         0, mpi_communicator);
+
+              return final_threshold;
+            }
 
           // count how many elements
           // there are that are bigger
@@ -569,7 +593,7 @@ namespace parallel
         // vector of indicators the
         // ones that correspond to
         // cells that we locally own
-        dealii::Vector<float>
+        dealii::Vector<typename Vector::value_type>
         locally_owned_indicators (tria.n_locally_owned_active_cells());
         get_locally_owned_indicators (tria,
                                       criteria,
@@ -583,7 +607,7 @@ namespace parallel
         // need it here, but it's a
         // collective communication
         // call
-        const std::pair<double,double> global_min_and_max
+        const std::pair<typename Vector::value_type,typename Vector::value_type> global_min_and_max
           = compute_global_min_and_max_at_root (locally_owned_indicators,
                                                 mpi_communicator);
 
@@ -681,7 +705,7 @@ namespace parallel
         // vector of indicators the
         // ones that correspond to
         // cells that we locally own
-        dealii::Vector<float>
+        dealii::Vector<typename Vector::value_type>
         locally_owned_indicators (tria.n_locally_owned_active_cells());
         get_locally_owned_indicators (tria,
                                       criteria,
@@ -718,7 +742,6 @@ namespace parallel
                                           top_fraction_of_error *
                                           total_error,
                                           mpi_communicator);
-
             // compute bottom
             // threshold only if
             // necessary. otherwise

diff --git a/tests/mpi/refine_and_coarsen_fixed_fraction_07.cc b/tests/mpi/refine_and_coarsen_fixed_fraction_07.cc
new file mode 100644 (file)
index 0000000..5612c11
--- /dev/null
+++ b/tests/mpi/refine_and_coarsen_fixed_fraction_07.cc
@@ -0,0 +1,221 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// A testcase by Andrea Bonito (in modified form): Given a particular
+// set of refinement indicators, we get no refinement at all because
+// the threshold is computed to be *larger* than the largest
+// refinement indicator
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/base/utilities.h>
+
+
+#include <fstream>
+
+
+void test()
+{
+  const unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+
+  // create a mesh with 120 cells (because that's what Andrea's
+  // original testcase had)
+  parallel::distributed::Triangulation<2> triangulation(MPI_COMM_WORLD);
+  std::vector<unsigned int> subdivisions(2);
+  subdivisions[0] = 120;
+  subdivisions[1] = 1;
+  GridGenerator::subdivided_hyper_rectangle (triangulation,
+                                            subdivisions,
+                                            Point<2>(),
+                                            Point<2>(1,1));
+  // initialize the refinement indicators with a set of particular values from the original
+  // testcase
+  const double values[] =
+    {
+         1.48589e-08,
+         3.31859e-06,
+         3.31859e-06,
+         0.0025918,
+         3.31859e-06,
+         0.0025918,
+         0.0025918,
+         3.31859e-06,
+         1.48589e-08,
+         0.0025918,
+         3.31859e-06,
+         0.0025918,
+         3.31859e-06,
+         0.0025918,
+         3.31859e-06,
+         0.0025918,
+         1.48589e-08,
+         3.31859e-06,
+         0.0025918,
+         3.31859e-06,
+         0.0025918,
+         0.0025918,
+         3.31859e-06,
+         3.31859e-06,
+         1.48589e-08,
+         0.0025918,
+         0.0025918,
+         3.31859e-06,
+         3.31859e-06,
+         0.0025918,
+         0.0025918,
+         1.48589e-08,
+         3.31859e-06,
+         3.31859e-06,
+         0.0025918,
+         0.0025918,
+         3.31859e-06,
+         0.0025918,
+         3.31859e-06,
+         1.48589e-08,
+         0.0025918,
+         3.31859e-06,
+         0.0025918,
+         3.31859e-06,
+         0.0025918,
+         3.31859e-06,
+         0.0025918,
+         1.48589e-08,
+         3.31859e-06,
+         0.0025918,
+         0.0025918,
+         3.31859e-06,
+         0.0025918,
+         3.31859e-06,
+         3.31859e-06,
+         1.48589e-08,
+         0.0120662,
+         0.0446999,
+         0.0446999,
+         0.0644361,
+         0.0446999,
+         0.0644361,
+         0.0644361,
+         3.14051,
+         0.0446999,
+         0.0120662,
+         0.0644361,
+         0.0446999,
+         0.0644361,
+         0.0446999,
+         3.14051,
+         0.0644361,
+         0.0446999,
+         0.0644361,
+         0.0120662,
+         0.0446999,
+         0.0644361,
+         3.14051,
+         0.0446999,
+         0.0644361,
+         0.0644361,
+         0.0446999,
+         0.0446999,
+         0.0120662,
+         3.14051,
+         0.0644361,
+         0.0644361,
+         0.0446999,
+         0.0446999,
+         0.0644361,
+         0.0644361,
+         3.14051,
+         0.0120662,
+         0.0446999,
+         0.0446999,
+         0.0644361,
+         0.0644361,
+         0.0446999,
+         3.14051,
+         0.0644361,
+         0.0446999,
+         0.0120662,
+         0.0644361,
+         0.0446999,
+         0.0644361,
+         3.14051,
+         0.0446999,
+         0.0644361,
+         0.0446999,
+         0.0644361,
+         0.0120662,
+         0.0446999,
+         3.14051,
+         0.0644361,
+         0.0644361,
+         0.0446999,
+         0.0644361,
+         0.0446999,
+         0.0446999,
+         0.0120662
+    };
+  
+  Vector<double> estimated_error_per_cell(&values[0],
+                                         &values[0] + 120);
+  
+  deallog << "Number of cells before refinement: " << triangulation.n_active_cells()
+         << std::endl;
+            
+  parallel::distributed::GridRefinement::
+    refine_and_coarsen_fixed_fraction (triangulation,
+                                      estimated_error_per_cell,
+                                      0.3, 0.0);
+  triangulation.execute_coarsening_and_refinement ();
+  
+  deallog << "Number of cells after refinement: " << triangulation.n_active_cells()
+         << std::endl;
+
+}
+
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+
+
+  deallog.push(Utilities::int_to_string(myid));
+
+  if (myid == 0)
+    {
+      std::ofstream logfile("output");
+      deallog.attach(logfile);
+      deallog.depth_console(0);
+      deallog.threshold_double(1.e-10);
+
+      test();
+    }
+  else
+    test();
+
+}

diff --git a/tests/mpi/refine_and_coarsen_fixed_fraction_07.mpirun=1.output b/tests/mpi/refine_and_coarsen_fixed_fraction_07.mpirun=1.output
new file mode 100644 (file)
index 0000000..5bf1b8a
--- /dev/null
+++ b/tests/mpi/refine_and_coarsen_fixed_fraction_07.mpirun=1.output
@@ -0,0 +1,3 @@
+
+DEAL:0::Number of cells before refinement: 120
+DEAL:0::Number of cells after refinement: 144