GridTools::find_active_cell_around_point() and find_all_active_cells_around_point...

diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 7045acce08ef1c6d85fe2b91b2647122cc37d2b6..e840d8955d461fb0556cf6d8824601a32375d0e4 100644 (file)
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -2998,8 +2998,8 @@ namespace GridTools
                   {
                     const Point<dim> p_unit =
                       mapping.transform_real_to_unit_cell(*cell, p);
-                    if (GeometryInfo<dim>::is_inside_unit_cell(p_unit,
-                                                               tolerance))
+                    if ((*cell)->reference_cell().contains_point(p_unit,
+                                                                 tolerance))
                       {
                         cell_and_position.first  = *cell;
                         cell_and_position.second = p_unit;
@@ -3011,8 +3011,8 @@ namespace GridTools
                     // outside it is and whether we want to use this cell as a
                     // backup if we can't find a cell within which the point
                     // lies.
-                    const double dist =
-                      GeometryInfo<dim>::distance_to_unit_cell(p_unit);
+                    const double dist = p_unit.distance(
+                      (*cell)->reference_cell().closest_point(p_unit));
                     if (dist < best_distance)
                       {
                         best_distance                   = dist;

diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc
index 9a93f8668f413362f0691d1030b7a9164f2f612e..e4bde46a7641dfa78f3f7da2dbcd75b65a98e874 100644 (file)
--- a/source/grid/grid_tools_dof_handlers.cc
+++ b/source/grid/grid_tools_dof_handlers.cc
@@ -629,97 +629,137 @@ namespace GridTools
     // insert the fist cell and point into the vector
     cells_and_points.push_back(first_cell);
 
-    // check if the given point is on the surface of the unit cell. If yes,
-    // need to find all neighbors
     const Point<dim> unit_point = cells_and_points.front().second;
     const auto       my_cell    = cells_and_points.front().first;
 
-    Tensor<1, dim> distance_to_center;
-    unsigned int   n_dirs_at_threshold     = 0;
-    unsigned int   last_point_at_threshold = numbers::invalid_unsigned_int;
-    for (unsigned int d = 0; d < dim; ++d)
+    std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
+      cells_to_add;
+
+    if (my_cell->reference_cell().is_hyper_cube())
       {
-        distance_to_center[d] = std::abs(unit_point[d] - 0.5);
-        if (distance_to_center[d] > 0.5 - tolerance)
+        // check if the given point is on the surface of the unit cell. If yes,
+        // need to find all neighbors
+
+        Tensor<1, dim> distance_to_center;
+        unsigned int   n_dirs_at_threshold     = 0;
+        unsigned int   last_point_at_threshold = numbers::invalid_unsigned_int;
+        for (unsigned int d = 0; d < dim; ++d)
           {
-            ++n_dirs_at_threshold;
-            last_point_at_threshold = d;
+            distance_to_center[d] = std::abs(unit_point[d] - 0.5);
+            if (distance_to_center[d] > 0.5 - tolerance)
+              {
+                ++n_dirs_at_threshold;
+                last_point_at_threshold = d;
+              }
           }
-      }
 
-    std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
-      cells_to_add;
-    // point is within face -> only need neighbor
-    if (n_dirs_at_threshold == 1)
-      {
-        unsigned int neighbor_index =
-          2 * last_point_at_threshold +
-          (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
-        if (!my_cell->at_boundary(neighbor_index))
+        // point is within face -> only need neighbor
+        if (n_dirs_at_threshold == 1)
           {
-            const auto neighbor_cell = my_cell->neighbor(neighbor_index);
+            unsigned int neighbor_index =
+              2 * last_point_at_threshold +
+              (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
+            if (!my_cell->at_boundary(neighbor_index))
+              {
+                const auto neighbor_cell = my_cell->neighbor(neighbor_index);
 
-            if (neighbor_cell->is_active())
-              cells_to_add.push_back(neighbor_cell);
-            else
-              for (const auto &child_cell : neighbor_cell->child_iterators())
-                {
-                  if (child_cell->is_active())
-                    cells_to_add.push_back(child_cell);
-                }
+                if (neighbor_cell->is_active())
+                  cells_to_add.push_back(neighbor_cell);
+                else
+                  for (const auto &child_cell :
+                       neighbor_cell->child_iterators())
+                    {
+                      if (child_cell->is_active())
+                        cells_to_add.push_back(child_cell);
+                    }
+              }
           }
-      }
-    // corner point -> use all neighbors
-    else if (n_dirs_at_threshold == dim)
-      {
-        unsigned int local_vertex_index = 0;
-        for (unsigned int d = 0; d < dim; ++d)
-          local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d;
+        // corner point -> use all neighbors
+        else if (n_dirs_at_threshold == dim)
+          {
+            unsigned int local_vertex_index = 0;
+            for (unsigned int d = 0; d < dim; ++d)
+              local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d;
 
-        const auto fu = [&](const auto &tentative_cells) {
-          for (const auto &cell : tentative_cells)
-            if (cell != my_cell)
-              cells_to_add.push_back(cell);
-        };
+            const auto fu = [&](const auto &tentative_cells) {
+              for (const auto &cell : tentative_cells)
+                if (cell != my_cell)
+                  cells_to_add.push_back(cell);
+            };
 
-        const auto vertex_index = my_cell->vertex_index(local_vertex_index);
+            const auto vertex_index = my_cell->vertex_index(local_vertex_index);
 
-        if (vertex_to_cells != nullptr)
-          fu((*vertex_to_cells)[vertex_index]);
-        else
-          fu(find_cells_adjacent_to_vertex(mesh, vertex_index));
-      }
-    // point on line in 3d: We cannot simply take the intersection between
-    // the two vertices of cells because of hanging nodes. So instead we
-    // list the vertices around both points and then select the
-    // appropriate cells according to the result of read_to_unit_cell
-    // below.
-    else if (n_dirs_at_threshold == 2)
-      {
-        std::pair<unsigned int, unsigned int> vertex_indices[3];
-        unsigned int                          count_vertex_indices = 0;
-        unsigned int free_direction = numbers::invalid_unsigned_int;
-        for (unsigned int d = 0; d < dim; ++d)
+            if (vertex_to_cells != nullptr)
+              fu((*vertex_to_cells)[vertex_index]);
+            else
+              fu(find_cells_adjacent_to_vertex(mesh, vertex_index));
+          }
+        // point on line in 3d: We cannot simply take the intersection between
+        // the two vertices of cells because of hanging nodes. So instead we
+        // list the vertices around both points and then select the
+        // appropriate cells according to the result of read_to_unit_cell
+        // below.
+        else if (n_dirs_at_threshold == 2)
           {
-            if (distance_to_center[d] > 0.5 - tolerance)
+            std::pair<unsigned int, unsigned int> vertex_indices[3];
+            unsigned int                          count_vertex_indices = 0;
+            unsigned int free_direction = numbers::invalid_unsigned_int;
+            for (unsigned int d = 0; d < dim; ++d)
               {
-                vertex_indices[count_vertex_indices].first = d;
-                vertex_indices[count_vertex_indices].second =
-                  unit_point[d] > 0.5 ? 1 : 0;
-                count_vertex_indices++;
+                if (distance_to_center[d] > 0.5 - tolerance)
+                  {
+                    vertex_indices[count_vertex_indices].first = d;
+                    vertex_indices[count_vertex_indices].second =
+                      unit_point[d] > 0.5 ? 1 : 0;
+                    count_vertex_indices++;
+                  }
+                else
+                  free_direction = d;
               }
-            else
-              free_direction = d;
-          }
 
-        AssertDimension(count_vertex_indices, 2);
-        Assert(free_direction != numbers::invalid_unsigned_int,
-               ExcInternalError());
+            AssertDimension(count_vertex_indices, 2);
+            Assert(free_direction != numbers::invalid_unsigned_int,
+                   ExcInternalError());
+
+            const unsigned int first_vertex =
+              (vertex_indices[0].second << vertex_indices[0].first) +
+              (vertex_indices[1].second << vertex_indices[1].first);
+            for (unsigned int d = 0; d < 2; ++d)
+              {
+                const auto fu = [&](const auto &tentative_cells) {
+                  for (const auto &cell : tentative_cells)
+                    {
+                      bool cell_not_yet_present = true;
+                      for (const auto &other_cell : cells_to_add)
+                        if (cell == other_cell)
+                          {
+                            cell_not_yet_present = false;
+                            break;
+                          }
+                      if (cell_not_yet_present)
+                        cells_to_add.push_back(cell);
+                    }
+                };
+
+                const auto vertex_index =
+                  my_cell->vertex_index(first_vertex + (d << free_direction));
+
+                if (vertex_to_cells != nullptr)
+                  fu((*vertex_to_cells)[vertex_index]);
+                else
+                  fu(find_cells_adjacent_to_vertex(mesh, vertex_index));
+              }
+          }
+      }
+    else
+      {
+        // Note: The non-hypercube path takes a very naive approach and
+        // checks all possible neighbors. This can be made faster by 1)
+        // checking if the point is in the inner cell and 2) identifying
+        // the right lines/vertices so that the number of potential
+        // neighbors is reduced.
 
-        const unsigned int first_vertex =
-          (vertex_indices[0].second << vertex_indices[0].first) +
-          (vertex_indices[1].second << vertex_indices[1].first);
-        for (unsigned int d = 0; d < 2; ++d)
+        for (const auto v : my_cell->vertex_indices())
           {
             const auto fu = [&](const auto &tentative_cells) {
               for (const auto &cell : tentative_cells)
@@ -736,8 +776,7 @@ namespace GridTools
                 }
             };
 
-            const auto vertex_index =
-              my_cell->vertex_index(first_vertex + (d << free_direction));
+            const auto vertex_index = my_cell->vertex_index(v);
 
             if (vertex_to_cells != nullptr)
               fu((*vertex_to_cells)[vertex_index]);
@@ -746,8 +785,6 @@ namespace GridTools
           }
       }
 
-    const double original_distance_to_unit_cell =
-      my_cell->reference_cell().closest_point(unit_point).distance(unit_point);
     for (const auto &cell : cells_to_add)
       {
         if (cell != my_cell)
@@ -755,8 +792,7 @@ namespace GridTools
             {
               const Point<dim> p_unit =
                 mapping.transform_real_to_unit_cell(cell, p);
-              if (cell->reference_cell().closest_point(p_unit).distance(
-                    p_unit) < original_distance_to_unit_cell + tolerance)
+              if (cell->reference_cell().contains_point(p_unit, tolerance))
                 cells_and_points.emplace_back(cell, p_unit);
             }
           catch (typename Mapping<dim>::ExcTransformationFailed &)

diff --git a/tests/remote_point_evaluation/find_all_active_cells_around_point_01.cc b/tests/remote_point_evaluation/find_all_active_cells_around_point_01.cc
new file mode 100644 (file)
index 0000000..1d6046a
--- /dev/null
+++ b/tests/remote_point_evaluation/find_all_active_cells_around_point_01.cc
@@ -0,0 +1,72 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test GridTools::find_active_cell_around_point() and
+// GridTools::find_all_active_cells_around_point for simplices. These
+// functions are used in find_all_locally_owned_active_cells_around_point(),
+// which is used by distributed_compute_point_locations().
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test()
+{
+  const auto  reference_cell = ReferenceCells::get_simplex<dim>();
+  const auto &mapping =
+    reference_cell.template get_default_linear_mapping<dim>();
+
+  Triangulation<dim> tria;
+  if (false)
+    GridGenerator::reference_cell(tria, reference_cell);
+  else
+    GridGenerator::subdivided_hyper_cube_with_simplices(tria, 2);
+
+  GridTools::Cache<dim> cache(tria, mapping);
+
+  const unsigned n_subdivisions = 8;
+
+  for (unsigned int i = 0; i <= n_subdivisions; ++i)
+    for (unsigned int j = 0; j <= n_subdivisions; ++j)
+      {
+        Point<dim> p(1.0 * i / n_subdivisions, 1.0 * j / n_subdivisions);
+
+        const auto first_point =
+          GridTools::find_active_cell_around_point(mapping, tria, p, {}, 1e-6);
+
+        const auto result = GridTools::find_all_active_cells_around_point(
+          mapping, tria, p, 1e-6, {first_point.first, first_point.second});
+
+        deallog << result.size() << std::endl; // should be 2
+        for (const auto &cell : result)
+          deallog << cell.first->id() << " " << cell.second << std::endl;
+        deallog << std::endl;
+      }
+}
+
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  deallog.precision(8);
+
+  test<2>();
+}

diff --git a/tests/remote_point_evaluation/find_all_active_cells_around_point_01.output b/tests/remote_point_evaluation/find_all_active_cells_around_point_01.output
new file mode 100644 (file)
index 0000000..70c3599
--- /dev/null
+++ b/tests/remote_point_evaluation/find_all_active_cells_around_point_01.output
@@ -0,0 +1,283 @@
+
+DEAL:0::1
+DEAL:0::0_0: 0.0000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::3
+DEAL:0::0_0: 0.0000000 1.0000000
+DEAL:0::1_0: 1.0000000 0.0000000
+DEAL:0::4_0: 0.0000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::4_0: 0.0000000 1.0000000
+DEAL:0::5_0: 1.0000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::0_0: 0.25000000 0.75000000
+DEAL:0::1_0: 0.75000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.75000000 0.0000000
+DEAL:0::4_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::4_0: 0.25000000 0.75000000
+DEAL:0::5_0: 0.75000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::0_0: 0.50000000 0.50000000
+DEAL:0::1_0: 0.50000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::1_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.50000000 0.0000000
+DEAL:0::4_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::4_0: 0.50000000 0.50000000
+DEAL:0::5_0: 0.50000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::0_0: 0.75000000 0.25000000
+DEAL:0::1_0: 0.25000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::1_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::1_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.25000000 0.0000000
+DEAL:0::4_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::4_0: 0.75000000 0.25000000
+DEAL:0::5_0: 0.25000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::3
+DEAL:0::0_0: 1.0000000 0.0000000
+DEAL:0::1_0: 0.0000000 1.0000000
+DEAL:0::2_0: 0.0000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.0000000 0.75000000
+DEAL:0::2_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.0000000 0.50000000
+DEAL:0::2_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.0000000 0.25000000
+DEAL:0::2_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::6
+DEAL:0::1_0: 0.0000000 0.0000000
+DEAL:0::2_0: 0.0000000 1.0000000
+DEAL:0::3_0: 1.0000000 0.0000000
+DEAL:0::4_0: 1.0000000 0.0000000
+DEAL:0::5_0: 0.0000000 1.0000000
+DEAL:0::6_0: 0.0000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::5_0: 0.0000000 0.75000000
+DEAL:0::6_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::5_0: 0.0000000 0.50000000
+DEAL:0::6_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::5_0: 0.0000000 0.25000000
+DEAL:0::6_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::3
+DEAL:0::5_0: 0.0000000 0.0000000
+DEAL:0::6_0: 0.0000000 1.0000000
+DEAL:0::7_0: 1.0000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::2_0: 0.25000000 0.75000000
+DEAL:0::3_0: 0.75000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::3_0: 0.75000000 0.0000000
+DEAL:0::6_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::6_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::6_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::6_0: 0.25000000 0.75000000
+DEAL:0::7_0: 0.75000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::2_0: 0.50000000 0.50000000
+DEAL:0::3_0: 0.50000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::3_0: 0.50000000 0.0000000
+DEAL:0::6_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::6_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::6_0: 0.50000000 0.50000000
+DEAL:0::7_0: 0.50000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::2_0: 0.75000000 0.25000000
+DEAL:0::3_0: 0.25000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::3_0: 0.25000000 0.0000000
+DEAL:0::6_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::6_0: 0.75000000 0.25000000
+DEAL:0::7_0: 0.25000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::2_0: 1.0000000 0.0000000
+DEAL:0::3_0: 0.0000000 1.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::3
+DEAL:0::3_0: 0.0000000 0.0000000
+DEAL:0::6_0: 1.0000000 0.0000000
+DEAL:0::7_0: 0.0000000 1.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.0000000 0.0000000
+DEAL:0::