]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add internal::MatrixFreeFunctions::DoFInfo::get_dof_indices_on_cell_batch() 7641/head
authorDenis Davydov <davydden@gmail.com>
Fri, 25 Jan 2019 21:18:17 +0000 (22:18 +0100)
committerDenis Davydov <davydden@gmail.com>
Wed, 30 Jan 2019 17:50:22 +0000 (18:50 +0100)
doc/doxygen/images/dofinfo_get_dof_indices.png [new file with mode: 0644]
doc/news/changes/minor/20190130DenisDavydov [new file with mode: 0644]
include/deal.II/matrix_free/dof_info.h
include/deal.II/matrix_free/dof_info.templates.h
tests/matrix_free/dof_info_01.cc [new file with mode: 0644]
tests/matrix_free/dof_info_01.with_mpi=true.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/matrix_free/dof_info_02.cc [new file with mode: 0644]
tests/matrix_free/dof_info_02.with_mpi=true.with_p4est=true.mpirun=1.output [new file with mode: 0644]

diff --git a/doc/doxygen/images/dofinfo_get_dof_indices.png b/doc/doxygen/images/dofinfo_get_dof_indices.png
new file mode 100644 (file)
index 0000000..8366cac
Binary files /dev/null and b/doc/doxygen/images/dofinfo_get_dof_indices.png differ
diff --git a/doc/news/changes/minor/20190130DenisDavydov b/doc/news/changes/minor/20190130DenisDavydov
new file mode 100644 (file)
index 0000000..4cf4ec9
--- /dev/null
@@ -0,0 +1,4 @@
+New: Add internal::MatrixFreeFunctions::DoFInfo::get_dof_indices_on_cell_batch() to
+return locally owned DoFs used by matrix-free framework on a given cell.
+<br>
+(Denis Davydov, 2019/01/30)
index 6839cd1d353da7ae290bd55930ed7cbe6a1830b1..8f80cf23342cbb9f894e53723aa31405d3cacc78 100644 (file)
@@ -102,6 +102,30 @@ namespace internal
       fe_index_from_degree(const unsigned int first_selected_component,
                            const unsigned int fe_degree) const;
 
+      /**
+       * Populate the vector @p locall_indices with locally owned degrees of freedom
+       * stored on the cell block @p cell.
+       * If @p with_constraints is `true`, then the returned vector will contain indices
+       * required to resolve constraints.
+       *
+       * The image below illustrates the output of this function for cell blocks
+       * zero and one with zero Dirichlet boundary conditions at the bottom of
+       * the domain. Note that due to the presence of constraints, the DoFs
+       * returned by this function for the case `with_constraints = true` are
+       * not a simple union
+       * of per cell DoFs on the cell block @p cell.
+       *
+       * @image html dofinfo_get_dof_indices.png
+       *
+       * @note The returned indices may contain duplicates. The unique set can be
+       * obtain using `std::sort()` followed by `std::unique()` and
+       * `std::vector::erase()`.
+       */
+      void
+      get_dof_indices_on_cell_batch(std::vector<unsigned int> &locall_indices,
+                                    const unsigned int         cell,
+                                    const bool with_constraints = true) const;
+
       /**
        * This internal method takes the local indices on a cell and fills them
        * into this class. It resolves the constraints and distributes the
index 153e6a8236fce0d83179452077505548bd2e8193..a2a7c6ae68f94f8eac4dc3cb4eb7e75b5ae29f48 100644 (file)
@@ -183,6 +183,73 @@ namespace internal
 
 
 
+    void
+    DoFInfo::get_dof_indices_on_cell_batch(std::vector<unsigned int> &my_rows,
+                                           const unsigned int         cell,
+                                           const bool apply_constraints) const
+    {
+      const unsigned int n_fe_components = start_components.back();
+      const unsigned int fe_index =
+        dofs_per_cell.size() == 1 ? 0 : cell_active_fe_index[cell];
+      const unsigned int dofs_this_cell = dofs_per_cell[fe_index];
+
+      const unsigned int n_vectorization  = vectorization_length;
+      constexpr auto     dof_access_index = dof_access_cell;
+      AssertIndexRange(cell,
+                       n_vectorization_lanes_filled[dof_access_index].size());
+      const unsigned int n_vectorization_actual =
+        n_vectorization_lanes_filled[dof_access_index][cell];
+
+      // we might have constraints, so the final number
+      // of indices is not known a priori.
+      // conservatively reserve the maximum without constraints
+      my_rows.reserve(n_vectorization * dofs_this_cell);
+      my_rows.resize(0);
+      unsigned int total_size = 0;
+      for (unsigned int v = 0; v < n_vectorization_actual; ++v)
+        {
+          const unsigned int ib =
+            (cell * n_vectorization + v) * n_fe_components;
+          const unsigned int ie =
+            (cell * n_vectorization + v + 1) * n_fe_components;
+
+          // figure out constraints by comparing constraint_indicator row
+          // shift for this cell within the block as compared to the next
+          // one
+          const bool has_constraints =
+            row_starts[ib].second != row_starts[ib + n_fe_components].second;
+
+          auto do_copy = [&](const unsigned int *begin,
+                             const unsigned int *end) {
+            const unsigned int shift = total_size;
+            total_size += (end - begin);
+            my_rows.resize(total_size);
+            std::copy(begin, end, my_rows.begin() + shift);
+          };
+
+          if (!has_constraints || apply_constraints)
+            {
+              const unsigned int *begin =
+                dof_indices.data() + row_starts[ib].first;
+              const unsigned int *end =
+                dof_indices.data() + row_starts[ie].first;
+              do_copy(begin, end);
+            }
+          else
+            {
+              Assert(row_starts_plain_indices[cell * n_vectorization + v] !=
+                       numbers::invalid_unsigned_int,
+                     ExcNotInitialized());
+              const unsigned int *begin =
+                plain_dof_indices.data() +
+                row_starts_plain_indices[cell * n_vectorization + v];
+              const unsigned int *end = begin + dofs_this_cell;
+              do_copy(begin, end);
+            }
+        }
+    }
+
+
     template <typename number>
     void
     DoFInfo ::read_dof_indices(
diff --git a/tests/matrix_free/dof_info_01.cc b/tests/matrix_free/dof_info_01.cc
new file mode 100644 (file)
index 0000000..462abca
--- /dev/null
@@ -0,0 +1,203 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests DoFInfo::get_dof_indices()
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim, int fe_degree, typename number = double>
+void
+test(const bool adaptive_ref = true)
+{
+  MPI_Comm           mpi_communicator(MPI_COMM_WORLD);
+  const unsigned int this_mpi_core =
+    dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
+
+  parallel::distributed::Triangulation<dim> tria(mpi_communicator);
+  GridGenerator::hyper_cube(tria, 0, 1, true);
+  tria.refine_global(1);
+  if (adaptive_ref)
+    {
+      for (auto cell : tria.active_cell_iterators())
+        if (cell->is_locally_owned())
+          if (cell->center().norm() < 0.5)
+            cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+      for (auto cell : tria.active_cell_iterators())
+        if (cell->is_locally_owned())
+          if (cell->center()[0] < 0.2)
+            cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+  else
+    {
+      tria.refine_global(1);
+    }
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+
+  IndexSet owned_set = dof.locally_owned_dofs();
+  IndexSet relevant_set;
+  DoFTools::extract_locally_relevant_dofs(dof, relevant_set);
+
+  AffineConstraints<double> constraints(relevant_set);
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  // constrain bottom part of the boundary (lower in y direction)
+  VectorTools::interpolate_boundary_values(dof,
+                                           2,
+                                           Functions::ZeroFunction<dim>(),
+                                           constraints);
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  std::shared_ptr<MatrixFree<dim, number>> mf_data(
+    new MatrixFree<dim, number>());
+  {
+    const QGauss<1>                                  quad(fe_degree + 1);
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+    data.tasks_block_size      = 7;
+    mf_data->reinit(dof, constraints, quad, data);
+  }
+
+  const unsigned int     n_cells  = mf_data->n_macro_cells();
+  const auto &           dof_info = mf_data->get_dof_info();
+  constexpr unsigned int n_vectorization =
+    VectorizedArray<number>::n_array_elements;
+
+  std::vector<unsigned int> my_rows;
+  my_rows.reserve(fe.dofs_per_cell * n_vectorization);
+
+  constexpr auto dof_access_index =
+    internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
+
+  bool checked_interleaved = false;
+  for (unsigned int cell = 0; cell < n_cells; ++cell)
+    {
+      deallog << "Cell: " << cell << std::endl;
+      checked_interleaved |=
+        (dof_info.index_storage_variants[dof_access_index][cell] ==
+         internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+           interleaved);
+      auto get_and_log = [&](const bool apply_constraints) {
+        dof_info.get_dof_indices_on_cell_batch(my_rows,
+                                               cell,
+                                               apply_constraints);
+        // sort and make unique to make visual inspection easier:
+        std::sort(my_rows.begin(), my_rows.end());
+        my_rows.erase(std::unique(my_rows.begin(), my_rows.end()),
+                      my_rows.end());
+        for (auto el : my_rows)
+          deallog << " " << el;
+        deallog << std::endl;
+      };
+      get_and_log(true);
+      get_and_log(false);
+    }
+
+  // if we don't have adaptive refinement, we should have
+  // some cells interleaved (at least that is the intention of this
+  // part of the test)
+  Assert(checked_interleaved || adaptive_ref, ExcInternalError());
+
+  // output in Gnuplot
+  if (dim == 2)
+    {
+      std::map<types::global_dof_index, Point<dim>> support_points;
+      MappingQ1<dim>                                mapping;
+      DoFTools::map_dofs_to_support_points(mapping, dof, support_points);
+
+      const std::string prefix =
+        std::is_same<float, number>::value ? "float_" : "double_";
+      const std::string href = (adaptive_ref ? "" : "global_");
+      const std::string base_filename =
+        prefix + href + "grid" + dealii::Utilities::int_to_string(dim) + "_" +
+        dealii::Utilities::int_to_string(fe_degree) + "_p" +
+        dealii::Utilities::int_to_string(this_mpi_core);
+
+      const std::string filename = base_filename + ".gp";
+      std::ofstream     f(filename.c_str());
+
+      f << "set terminal png size 400,410 enhanced font \"Helvetica,8\""
+        << std::endl
+        << "set output \"" << base_filename << ".png\"" << std::endl
+        << "set size square" << std::endl
+        << "set view equal xy" << std::endl
+        << "unset xtics" << std::endl
+        << "unset ytics" << std::endl
+        << "unset grid" << std::endl
+        << "unset border" << std::endl
+        << "plot '-' using 1:2 with lines notitle, '-' with labels tc rgb 'red' nopoint notitle, '-' with labels point pt 4 offset 1,1 notitle"
+        << std::endl;
+      GridOut().write_gnuplot(tria, f);
+      f << "e" << std::endl;
+
+      // output cell blocks:
+      for (unsigned int cell = 0; cell < n_cells; ++cell)
+        for (unsigned int c = 0; c < mf_data->n_components_filled(cell); ++c)
+          {
+            const auto dof_cell = mf_data->get_cell_iterator(cell, c);
+            f << dof_cell->center() << " \"" << cell << "\"\n";
+          }
+
+      f << std::flush;
+      f << "e" << std::endl << std::endl;
+
+      DoFTools::write_gnuplot_dof_support_point_info(f, support_points);
+
+      f << "e" << std::endl;
+    }
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+  MPILogInitAll mpi_init_log;
+
+  deallog.push("2d");
+  test<2, 1>();
+  test<2, 1, float>();
+  test<2, 1>(false);
+  deallog.pop();
+}
diff --git a/tests/matrix_free/dof_info_01.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/dof_info_01.with_mpi=true.with_p4est=true.mpirun=1.output
new file mode 100644 (file)
index 0000000..226a54f
--- /dev/null
@@ -0,0 +1,34 @@
+
+DEAL:0:2d::Testing FE_Q<2>(1)
+DEAL:0:2d::Cell: 0
+DEAL:0:2d:: 7 18 19 21 22
+DEAL:0:2d:: 6 7 16 17 18 19 20 21 22
+DEAL:0:2d::Cell: 1
+DEAL:0:2d:: 2 7 9 10 21 22 23 24
+DEAL:0:2d:: 0 6 7 8 10 21 22 23 24 25 26
+DEAL:0:2d::Cell: 2
+DEAL:0:2d:: 2 3 7 9 10 11 12 24
+DEAL:0:2d:: 0 1 2 3 7 8 9 10 11 12 24 25 26
+DEAL:0:2d::Cell: 3
+DEAL:0:2d:: 2 3 4 5 9 11 12 14 15
+DEAL:0:2d:: 2 3 4 5 9 11 12 13 14 15
+DEAL:0:2d::Testing FE_Q<2>(1)
+DEAL:0:2d::Cell: 0
+DEAL:0:2d:: 2 7 9 10 18 19 21 22 23 24
+DEAL:0:2d:: 0 6 7 8 10 16 17 18 19 20 21 22 23 24 25 26
+DEAL:0:2d::Cell: 1
+DEAL:0:2d:: 2 3 4 5 7 9 10 11 12 14 15 24
+DEAL:0:2d:: 0 1 2 3 4 5 7 8 9 10 11 12 13 14 15 24 25 26
+DEAL:0:2d::Testing FE_Q<2>(1)
+DEAL:0:2d::Cell: 0
+DEAL:0:2d:: 2 3 5 6 7 8
+DEAL:0:2d:: 0 1 2 3 4 5 6 7 8
+DEAL:0:2d::Cell: 1
+DEAL:0:2d:: 5 8 10 12 13 14
+DEAL:0:2d:: 4 5 8 9 10 11 12 13 14
+DEAL:0:2d::Cell: 2
+DEAL:0:2d:: 6 7 8 15 16 17 18 19 20
+DEAL:0:2d:: 6 7 8 15 16 17 18 19 20
+DEAL:0:2d::Cell: 3
+DEAL:0:2d:: 8 13 14 17 20 21 22 23 24
+DEAL:0:2d:: 8 13 14 17 20 21 22 23 24
diff --git a/tests/matrix_free/dof_info_02.cc b/tests/matrix_free/dof_info_02.cc
new file mode 100644 (file)
index 0000000..a980641
--- /dev/null
@@ -0,0 +1,205 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests DoFInfo::get_dof_indices() for multi-component system
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim, int fe_degree, typename number = double>
+void
+test(const bool adaptive_ref = true)
+{
+  MPI_Comm           mpi_communicator(MPI_COMM_WORLD);
+  const unsigned int this_mpi_core =
+    dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
+
+  parallel::distributed::Triangulation<dim> tria(mpi_communicator);
+  GridGenerator::hyper_cube(tria, 0, 1, true);
+  tria.refine_global(1);
+  if (adaptive_ref)
+    {
+      for (auto cell : tria.active_cell_iterators())
+        if (cell->is_locally_owned())
+          if (cell->center().norm() < 0.5)
+            cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+      for (auto cell : tria.active_cell_iterators())
+        if (cell->is_locally_owned())
+          if (cell->center()[0] < 0.2)
+            cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+  else
+    {
+      tria.refine_global(1);
+    }
+
+  FE_Q<dim>       feq(fe_degree);
+  FESystem<dim>   fe(feq, dim);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+
+  IndexSet owned_set = dof.locally_owned_dofs();
+  IndexSet relevant_set;
+  DoFTools::extract_locally_relevant_dofs(dof, relevant_set);
+
+  AffineConstraints<double> constraints(relevant_set);
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  // constrain bottom part of the boundary (lower in y direction)
+  VectorTools::interpolate_boundary_values(dof,
+                                           2,
+                                           Functions::ZeroFunction<dim>(2),
+                                           constraints);
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  std::shared_ptr<MatrixFree<dim, number>> mf_data(
+    new MatrixFree<dim, number>());
+  {
+    const QGauss<1>                                  quad(fe_degree + 1);
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+    data.tasks_block_size      = 7;
+    mf_data->reinit(dof, constraints, quad, data);
+  }
+
+  const unsigned int     n_cells  = mf_data->n_macro_cells();
+  const auto &           dof_info = mf_data->get_dof_info();
+  constexpr unsigned int n_vectorization =
+    VectorizedArray<number>::n_array_elements;
+
+  std::vector<unsigned int> my_rows;
+  my_rows.reserve(fe.dofs_per_cell * n_vectorization);
+
+  constexpr auto dof_access_index =
+    internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
+
+  bool checked_interleaved = false;
+  for (unsigned int cell = 0; cell < n_cells; ++cell)
+    {
+      deallog << "Cell: " << cell << std::endl;
+      checked_interleaved |=
+        (dof_info.index_storage_variants[dof_access_index][cell] ==
+         internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+           interleaved);
+      auto get_and_log = [&](const bool apply_constraints) {
+        dof_info.get_dof_indices_on_cell_batch(my_rows,
+                                               cell,
+                                               apply_constraints);
+        // sort and make unique to make visual inspection easier:
+        std::sort(my_rows.begin(), my_rows.end());
+        my_rows.erase(std::unique(my_rows.begin(), my_rows.end()),
+                      my_rows.end());
+        for (auto el : my_rows)
+          deallog << " " << el;
+        deallog << std::endl;
+      };
+      get_and_log(true);
+      get_and_log(false);
+    }
+
+  // if we don't have adaptive refinement, we should have
+  // some cells interleaved (at least that is the intention of this
+  // part of the test)
+  Assert(checked_interleaved || adaptive_ref, ExcInternalError());
+
+  // output in Gnuplot
+  if (dim == 2)
+    {
+      std::map<types::global_dof_index, Point<dim>> support_points;
+      MappingQ1<dim>                                mapping;
+      DoFTools::map_dofs_to_support_points(mapping, dof, support_points);
+
+      const std::string prefix =
+        std::is_same<float, number>::value ? "float_" : "double_";
+      const std::string href = (adaptive_ref ? "" : "global_");
+      const std::string base_filename =
+        prefix + href + "grid" + dealii::Utilities::int_to_string(dim) + "_" +
+        dealii::Utilities::int_to_string(fe_degree) + "_p" +
+        dealii::Utilities::int_to_string(this_mpi_core);
+
+      const std::string filename = base_filename + ".gp";
+      std::ofstream     f(filename.c_str());
+
+      f << "set terminal png size 400,410 enhanced font \"Helvetica,8\""
+        << std::endl
+        << "set output \"" << base_filename << ".png\"" << std::endl
+        << "set size square" << std::endl
+        << "set view equal xy" << std::endl
+        << "unset xtics" << std::endl
+        << "unset ytics" << std::endl
+        << "unset grid" << std::endl
+        << "unset border" << std::endl
+        << "plot '-' using 1:2 with lines notitle, '-' with labels tc rgb 'red' nopoint notitle, '-' with labels point pt 4 offset 1,1 notitle"
+        << std::endl;
+      GridOut().write_gnuplot(tria, f);
+      f << "e" << std::endl;
+
+      // output cell blocks:
+      for (unsigned int cell = 0; cell < n_cells; ++cell)
+        for (unsigned int c = 0; c < mf_data->n_components_filled(cell); ++c)
+          {
+            const auto dof_cell = mf_data->get_cell_iterator(cell, c);
+            f << dof_cell->center() << " \"" << cell << "\"\n";
+          }
+
+      f << std::flush;
+      f << "e" << std::endl << std::endl;
+
+      DoFTools::write_gnuplot_dof_support_point_info(f, support_points);
+
+      f << "e" << std::endl;
+    }
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+  MPILogInitAll mpi_init_log;
+
+  deallog.push("2d");
+  test<2, 1>();
+  test<2, 1, float>();
+  test<2, 1>(false);
+  deallog.pop();
+}
diff --git a/tests/matrix_free/dof_info_02.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/dof_info_02.with_mpi=true.with_p4est=true.mpirun=1.output
new file mode 100644 (file)
index 0000000..1da9dc6
--- /dev/null
@@ -0,0 +1,34 @@
+
+DEAL:0:2d::Testing FESystem<2>[FE_Q<2>(1)^2]
+DEAL:0:2d::Cell: 0
+DEAL:0:2d:: 14 15 36 37 38 39 42 43 44 45
+DEAL:0:2d:: 12 13 14 15 32 33 34 35 36 37 38 39 40 41 42 43 44 45
+DEAL:0:2d::Cell: 1
+DEAL:0:2d:: 4 5 14 15 18 19 20 21 42 43 44 45 46 47 48 49
+DEAL:0:2d:: 0 1 12 13 14 15 16 17 20 21 42 43 44 45 46 47 48 49 50 51 52 53
+DEAL:0:2d::Cell: 2
+DEAL:0:2d:: 4 5 6 7 14 15 18 19 20 21 22 23 24 25 48 49
+DEAL:0:2d:: 0 1 2 3 4 5 6 7 14 15 16 17 18 19 20 21 22 23 24 25 48 49 50 51 52 53
+DEAL:0:2d::Cell: 3
+DEAL:0:2d:: 4 5 6 7 8 9 10 11 18 19 22 23 24 25 28 29 30 31
+DEAL:0:2d:: 4 5 6 7 8 9 10 11 18 19 22 23 24 25 26 27 28 29 30 31
+DEAL:0:2d::Testing FESystem<2>[FE_Q<2>(1)^2]
+DEAL:0:2d::Cell: 0
+DEAL:0:2d:: 4 5 14 15 18 19 20 21 36 37 38 39 42 43 44 45 46 47 48 49
+DEAL:0:2d:: 0 1 12 13 14 15 16 17 20 21 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
+DEAL:0:2d::Cell: 1
+DEAL:0:2d:: 4 5 6 7 8 9 10 11 14 15 18 19 20 21 22 23 24 25 28 29 30 31 48 49
+DEAL:0:2d:: 0 1 2 3 4 5 6 7 8 9 10 11 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 48 49 50 51 52 53
+DEAL:0:2d::Testing FESystem<2>[FE_Q<2>(1)^2]
+DEAL:0:2d::Cell: 0
+DEAL:0:2d:: 4 5 6 7 10 11 12 13 14 15 16 17
+DEAL:0:2d:: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
+DEAL:0:2d::Cell: 1
+DEAL:0:2d:: 10 11 16 17 20 21 24 25 26 27 28 29
+DEAL:0:2d:: 8 9 10 11 16 17 18 19 20 21 22 23 24 25 26 27 28 29
+DEAL:0:2d::Cell: 2
+DEAL:0:2d:: 12 13 14 15 16 17 30 31 32 33 34 35 36 37 38 39 40 41
+DEAL:0:2d:: 12 13 14 15 16 17 30 31 32 33 34 35 36 37 38 39 40 41
+DEAL:0:2d::Cell: 3
+DEAL:0:2d:: 16 17 26 27 28 29 34 35 40 41 42 43 44 45 46 47 48 49
+DEAL:0:2d:: 16 17 26 27 28 29 34 35 40 41 42 43 44 45 46 47 48 49

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.