]> https://gitweb.dealii.org/ - dealii.git/commitdiff
no normal flux on levels
authorJiaqi Zhang <jiaqi2@clemson.edu>
Thu, 28 Jul 2022 14:57:28 +0000 (10:57 -0400)
committerJiaqi Zhang <jiaqi2@clemson.edu>
Wed, 17 Aug 2022 16:09:20 +0000 (12:09 -0400)
include/deal.II/numerics/vector_tools_constraints.h
include/deal.II/numerics/vector_tools_constraints.templates.h
source/numerics/vector_tools_constraints.inst.in
tests/numerics/no_flux_16.cc [new file with mode: 0644]
tests/numerics/no_flux_16.output [new file with mode: 0644]

index 5b3aa88d08b8cdc15ca539c4e7c3f37f5a584f7d..aa36668448bd610f12481b5cc5e600d54def29c1 100644 (file)
@@ -21,6 +21,8 @@
 
 #include <deal.II/dofs/dof_handler.h>
 
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
 #include <map>
 #include <set>
 
@@ -316,6 +318,33 @@ namespace VectorTools
 #endif
          ));
 
+  /**
+   * This function does the same as the
+   * compute_no_normal_flux_constraints(), but for the case of level meshes
+   * in the multigrid method.
+   * @ingroup constraints
+   *
+   * @see
+   * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
+   */
+  template <int dim, int spacedim>
+  void
+  compute_no_normal_flux_constraints_on_level(
+    const DoFHandler<dim, spacedim> &   dof_handler,
+    const MGConstrainedDoFs &           mg_constrained_dofs,
+    const unsigned int                  level,
+    const unsigned int                  first_vector_component,
+    const std::set<types::boundary_id> &boundary_ids,
+    AffineConstraints<double> &         constraints,
+    const Mapping<dim> &                mapping =
+      (ReferenceCells::get_hypercube<dim>()
+#ifndef _MSC_VER
+         .template get_default_linear_mapping<dim, spacedim>()
+#else
+         .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
+#endif
+         ));
+
   /**
    * Compute the constraints that correspond to boundary conditions of the
    * form $\vec u \times \vec n=\vec u_\Gamma \times \vec n$, i.e., tangential
index a556af03fe82ef68a1fbc4da687a5ddac17d82dc..08c2a44a70ebf9a9e0c43976d4ad75336edc7d7b 100644 (file)
@@ -1293,6 +1293,433 @@ namespace VectorTools
                                             mapping);
   }
 
+  template <int dim, int spacedim>
+  void
+  compute_no_normal_flux_constraints_on_level(
+    const DoFHandler<dim, spacedim> &   dof_handler,
+    const MGConstrainedDoFs &           mg_constrained_dofs,
+    const unsigned int                  level,
+    const unsigned int                  first_vector_component,
+    const std::set<types::boundary_id> &boundary_ids,
+    AffineConstraints<double> &         constraints,
+    const Mapping<dim> &                mapping)
+  {
+    // Copied from compute_nonzero_normal_flux_constraints()
+    // Only changed active_cell_iterator to cell_iterator in DoFToNormalsMap,
+    // CellToNormalsMap, CellContributions, and set inhomogeneity to 0.
+    const IndexSet &refinement_edge_indices =
+      mg_constrained_dofs.get_refinement_edge_indices(level);
+
+    const auto &                       fe = dof_handler.get_fe();
+    const std::vector<Point<dim - 1>> &unit_support_points =
+      fe.get_unit_face_support_points();
+    const Quadrature<dim - 1>            quadrature(unit_support_points);
+    const unsigned int                   dofs_per_face = fe.dofs_per_face;
+    std::vector<types::global_dof_index> face_dofs(dofs_per_face);
+
+
+    FEFaceValues<dim, spacedim> fe_face_values(mapping,
+                                               fe,
+                                               quadrature,
+                                               update_quadrature_points |
+                                                 update_normal_vectors);
+
+    std::set<types::boundary_id>::iterator b_id;
+    using DoFToNormalsMap = std::multimap<
+      internal::VectorDoFTuple<dim>,
+      std::pair<Tensor<1, dim>,
+                typename DoFHandler<dim, spacedim>::cell_iterator>>;
+    DoFToNormalsMap dof_to_normals_map;
+    for (const auto &cell : dof_handler.cell_iterators_on_level(level))
+      if (cell->level_subdomain_id() != numbers::artificial_subdomain_id &&
+          cell->level_subdomain_id() != numbers::invalid_subdomain_id)
+        for (const unsigned int face_no : cell->face_indices())
+          if ((b_id = boundary_ids.find(cell->face(face_no)->boundary_id())) !=
+              boundary_ids.end())
+            {
+              typename DoFHandler<dim, spacedim>::level_face_iterator face =
+                cell->face(face_no);
+              face->get_mg_dof_indices(level, face_dofs);
+              fe_face_values.reinit(cell, face_no);
+
+              for (unsigned int i = 0; i < face_dofs.size(); ++i)
+                if (fe.face_system_to_component_index(i).first ==
+                    first_vector_component)
+                  // Refinement edge indices are going to be constrained to 0
+                  // during a multigrid cycle and do not need no-normal-flux
+                  // constraints, so skip them:
+                  if (!refinement_edge_indices.is_element(face_dofs[i]))
+                    {
+                      const Point<dim> position =
+                        fe_face_values.quadrature_point(i);
+                      internal::VectorDoFTuple<dim> vector_dofs;
+                      vector_dofs.dof_indices[0] = face_dofs[i];
+                      for (unsigned int k = 0; k < dofs_per_face; ++k)
+                        if ((k != i) &&
+                            (quadrature.point(k) == quadrature.point(i)) &&
+                            (fe.face_system_to_component_index(k).first >=
+                             first_vector_component) &&
+                            (fe.face_system_to_component_index(k).first <
+                             first_vector_component + dim))
+                          vector_dofs.dof_indices
+                            [fe.face_system_to_component_index(k).first -
+                             first_vector_component] = face_dofs[k];
+
+                      Tensor<1, dim> normal_vector =
+                        cell->face(face_no)->get_manifold().normal_vector(
+                          cell->face(face_no), position);
+
+                      // make sure the normal vector is pointing to the right
+                      // direction, more dietails can be found in
+                      // compute_nonzero_normal_flux_constraints()
+                      if (normal_vector * fe_face_values.normal_vector(i) < 0)
+                        normal_vector *= -1;
+                      Assert(std::fabs(normal_vector.norm() - 1) < 1e-14,
+                             ExcInternalError());
+
+                      // remove small entries:
+                      for (unsigned int d = 0; d < dim; ++d)
+                        if (std::fabs(normal_vector[d]) < 1e-13)
+                          normal_vector[d] = 0;
+                      normal_vector /= normal_vector.norm();
+
+                      dof_to_normals_map.insert(
+                        std::make_pair(vector_dofs,
+                                       std::make_pair(normal_vector, cell)));
+#ifdef DEBUG_NO_NORMAL_FLUX
+                      std::cout << "Adding normal vector:" << std::endl
+                                << "   dofs=" << vector_dofs << std::endl
+                                << "   cell=" << cell << " at "
+                                << cell->center() << std::endl
+                                << "   normal=" << normal_vector << std::endl;
+#endif
+                    }
+            }
+    // Now do something with the collected information. To this end, loop
+    // through all sets of pairs (dofs,normal_vector) and identify which
+    // entries belong to the same set of dofs and then do as described in the
+    // documentation, i.e. either average the normal vector or don't for this
+    // particular set of dofs
+    typename DoFToNormalsMap::const_iterator p = dof_to_normals_map.begin();
+
+    while (p != dof_to_normals_map.end())
+      {
+        // first find the range of entries in the multimap that corresponds to
+        // the same vector-dof tuple. as usual, we define the range
+        // half-open. the first entry of course is 'p'
+        typename DoFToNormalsMap::const_iterator same_dof_range[2] = {p};
+        for (++p; p != dof_to_normals_map.end(); ++p)
+          if (p->first != same_dof_range[0]->first)
+            {
+              same_dof_range[1] = p;
+              break;
+            }
+        if (p == dof_to_normals_map.end())
+          same_dof_range[1] = dof_to_normals_map.end();
+
+#ifdef DEBUG_NO_NORMAL_FLUX
+        std::cout << "For dof indices <" << p->first
+                  << ">, found the following normals" << std::endl;
+        for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
+             q != same_dof_range[1];
+             ++q)
+          std::cout << "   " << q->second.first << " from cell "
+                    << q->second.second << std::endl;
+#endif
+
+
+        // now compute the reverse mapping: for each of the cells that
+        // contributed to the current set of vector dofs, add up the normal
+        // vectors. the values of the map are pairs of normal vectors and
+        // number of cells that have contributed
+        using CellToNormalsMap =
+          std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
+                   std::pair<Tensor<1, dim>, unsigned int>>;
+
+        CellToNormalsMap cell_to_normals_map;
+        for (typename DoFToNormalsMap::const_iterator q = same_dof_range[0];
+             q != same_dof_range[1];
+             ++q)
+          if (cell_to_normals_map.find(q->second.second) ==
+              cell_to_normals_map.end())
+            cell_to_normals_map[q->second.second] =
+              std::make_pair(q->second.first, 1U);
+          else
+            {
+              const Tensor<1, dim> old_normal =
+                cell_to_normals_map[q->second.second].first;
+              const unsigned int old_count =
+                cell_to_normals_map[q->second.second].second;
+
+              Assert(old_count > 0, ExcInternalError());
+
+              // in the same entry, store again the now averaged normal vector
+              // and the new count
+              cell_to_normals_map[q->second.second] =
+                std::make_pair((old_normal * old_count + q->second.first) /
+                                 (old_count + 1),
+                               old_count + 1);
+            }
+        Assert(cell_to_normals_map.size() >= 1, ExcInternalError());
+
+#ifdef DEBUG_NO_NORMAL_FLUX
+        std::cout << "   cell_to_normals_map:" << std::endl;
+        for (typename CellToNormalsMap::const_iterator x =
+               cell_to_normals_map.begin();
+             x != cell_to_normals_map.end();
+             ++x)
+          std::cout << "      " << x->first << " -> (" << x->second.first << ','
+                    << x->second.second << ')' << std::endl;
+#endif
+
+        // count the maximum number of contributions from each cell
+        unsigned int max_n_contributions_per_cell = 1;
+        for (typename CellToNormalsMap::const_iterator x =
+               cell_to_normals_map.begin();
+             x != cell_to_normals_map.end();
+             ++x)
+          max_n_contributions_per_cell =
+            std::max(max_n_contributions_per_cell, x->second.second);
+
+        // verify that each cell can have only contributed at most dim times,
+        // since that is the maximum number of faces that come together at a
+        // single place
+        Assert(max_n_contributions_per_cell <= dim, ExcInternalError());
+
+        switch (max_n_contributions_per_cell)
+          {
+            case 1:
+              {
+                // compute the average normal vector from all the ones that
+                // have the same set of dofs. we could add them up and divide
+                // them by the number of additions, or simply normalize them
+                // right away since we want them to have unit length anyway
+                Tensor<1, dim> normal;
+                for (typename CellToNormalsMap::const_iterator x =
+                       cell_to_normals_map.begin();
+                     x != cell_to_normals_map.end();
+                     ++x)
+                  normal += x->second.first;
+                normal /= normal.norm();
+
+                // normalize again
+                for (unsigned int d = 0; d < dim; ++d)
+                  if (std::fabs(normal[d]) < 1e-13)
+                    normal[d] = 0;
+                normal /= normal.norm();
+
+                const auto &dof_indices = same_dof_range[0]->first;
+                internal::add_constraint<dim>(dof_indices,
+                                              normal,
+                                              constraints,
+                                              0.0);
+
+                break;
+              }
+
+            case dim:
+              {
+                // assert that indeed only a single cell has contributed
+                Assert(cell_to_normals_map.size() == 1, ExcInternalError());
+
+                // check linear independence by computing the determinant of
+                // the matrix created from all the normal vectors. if they are
+                // linearly independent, then the determinant is nonzero. if
+                // they are orthogonal, then the matrix is in fact equal to 1
+                // (since they are all unit vectors); make sure the
+                // determinant is larger than 1e-3 to avoid cases where cells
+                // are degenerate
+                {
+                  Tensor<2, dim> t;
+
+                  typename DoFToNormalsMap::const_iterator x =
+                    same_dof_range[0];
+                  for (unsigned int i = 0; i < dim; ++i, ++x)
+                    for (unsigned int j = 0; j < dim; ++j)
+                      t[i][j] = x->second.first[j];
+
+                  Assert(
+                    std::fabs(determinant(t)) > 1e-3,
+                    ExcMessage(
+                      "Found a set of normal vectors that are nearly collinear."));
+                }
+
+                // so all components of this vector dof are constrained. enter
+                // this into the AffineConstraints object
+                //
+                // ignore dofs already constrained
+                const auto &dof_indices = same_dof_range[0]->first;
+
+                for (unsigned int i = 0; i < dim; ++i)
+                  if (!constraints.is_constrained(
+                        same_dof_range[0]->first.dof_indices[i]) &&
+                      constraints.can_store_line(
+                        same_dof_range[0]->first.dof_indices[i]))
+                    {
+                      const types::global_dof_index line =
+                        dof_indices.dof_indices[i];
+                      constraints.add_line(line);
+                    }
+
+                break;
+              }
+
+            // this is the case of an edge contribution in 3d, i.e. the vector
+            // is constrained in two directions but not the third.
+            default:
+              {
+                Assert(dim >= 3, ExcNotImplemented());
+                Assert(max_n_contributions_per_cell == 2, ExcInternalError());
+
+                // as described in the documentation, let us first collect
+                // what each of the cells contributed at the current point. we
+                // use a std::list instead of a std::set (which would be more
+                // natural) because std::set requires that the stored elements
+                // are comparable with operator<
+                using CellContributions =
+                  std::map<typename DoFHandler<dim, spacedim>::cell_iterator,
+                           std::list<Tensor<1, dim>>>;
+                CellContributions cell_contributions;
+
+                for (typename DoFToNormalsMap::const_iterator q =
+                       same_dof_range[0];
+                     q != same_dof_range[1];
+                     ++q)
+                  cell_contributions[q->second.second].push_back(
+                    q->second.first);
+                Assert(cell_contributions.size() >= 1, ExcInternalError());
+
+                // now for each cell that has contributed determine the number
+                // of normal vectors it has contributed. we currently only
+                // implement if this is dim-1 for all cells (if a single cell
+                // has contributed dim, or if all adjacent cells have
+                // contributed 1 normal vector, this is already handled
+                // above).
+                //
+                // we only implement the case that all cells contribute
+                // dim-1 because we assume that we are following an edge
+                // of the domain (think: we are looking at a vertex
+                // located on one of the edges of a refined cube where the
+                // boundary indicators of the two adjacent faces of the
+                // cube are both listed in the set of boundary indicators
+                // passed to this function). in that case, all cells along
+                // that edge of the domain are assumed to have contributed
+                // dim-1 normal vectors. however, there are cases where
+                // this assumption is not justified (see the lengthy
+                // explanation in test no_flux_12.cc) and in those cases
+                // we simply ignore the cell that contributes only
+                // once. this is also discussed at length in the
+                // documentation of this function.
+                //
+                // for each contributing cell compute the tangential vector
+                // that remains unconstrained
+                std::list<Tensor<1, dim>> tangential_vectors;
+                for (typename CellContributions::const_iterator contribution =
+                       cell_contributions.begin();
+                     contribution != cell_contributions.end();
+                     ++contribution)
+                  {
+#ifdef DEBUG_NO_NORMAL_FLUX
+                    std::cout
+                      << "   Treating edge case with dim-1 contributions."
+                      << std::endl
+                      << "   Looking at cell " << contribution->first
+                      << " which has contributed these normal vectors:"
+                      << std::endl;
+                    for (typename std::list<Tensor<1, dim>>::const_iterator t =
+                           contribution->second.begin();
+                         t != contribution->second.end();
+                         ++t)
+                      std::cout << "      " << *t << std::endl;
+#endif
+
+                    // as mentioned above, simply ignore cells that only
+                    // contribute once
+                    if (contribution->second.size() < dim - 1)
+                      continue;
+
+                    Tensor<1, dim> normals[dim - 1];
+                    {
+                      unsigned int index = 0;
+                      for (typename std::list<Tensor<1, dim>>::const_iterator
+                             t = contribution->second.begin();
+                           t != contribution->second.end();
+                           ++t, ++index)
+                        normals[index] = *t;
+                      Assert(index == dim - 1, ExcInternalError());
+                    }
+
+                    // calculate the tangent as the outer product of the
+                    // normal vectors. since these vectors do not need to be
+                    // orthogonal (think, for example, the case of the
+                    // deal.II/no_flux_07 test: a sheared cube in 3d, with Q2
+                    // elements, where we have constraints from the two normal
+                    // vectors of two faces of the sheared cube that are not
+                    // perpendicular to each other), we have to normalize the
+                    // outer product
+                    Tensor<1, dim> tangent;
+                    switch (dim)
+                      {
+                        case 3:
+                          // take cross product between normals[0] and
+                          // normals[1]. write it in the current form (with
+                          // [dim-2]) to make sure that compilers don't warn
+                          // about out-of-bounds accesses -- the warnings are
+                          // bogus since we get here only for dim==3, but at
+                          // least one isn't quite smart enough to notice this
+                          // and warns when compiling the function in 2d
+                          tangent =
+                            cross_product_3d(normals[0], normals[dim - 2]);
+                          break;
+                        default:
+                          Assert(false, ExcNotImplemented());
+                      }
+
+                    Assert(
+                      std::fabs(tangent.norm()) > 1e-12,
+                      ExcMessage(
+                        "Two normal vectors from adjacent faces are almost "
+                        "parallel."));
+                    tangent /= tangent.norm();
+
+                    tangential_vectors.push_back(tangent);
+                  }
+
+                // go through the list of tangents and make sure that they all
+                // roughly point in the same direction as the first one (i.e.
+                // have an angle less than 90 degrees); if they don't then
+                // flip their sign
+                {
+                  const Tensor<1, dim> first_tangent =
+                    tangential_vectors.front();
+                  typename std::list<Tensor<1, dim>>::iterator t =
+                    tangential_vectors.begin();
+                  ++t;
+                  for (; t != tangential_vectors.end(); ++t)
+                    if (*t * first_tangent < 0)
+                      *t *= -1;
+                }
+
+                // now compute the average tangent and normalize it
+                Tensor<1, dim> average_tangent;
+                for (typename std::list<Tensor<1, dim>>::const_iterator t =
+                       tangential_vectors.begin();
+                     t != tangential_vectors.end();
+                     ++t)
+                  average_tangent += *t;
+                average_tangent /= average_tangent.norm();
+
+                const auto &dof_indices = same_dof_range[0]->first;
+                internal::add_tangentiality_constraints(dof_indices,
+                                                        average_tangent,
+                                                        constraints);
+              }
+          }
+      }
+  }
+
+
+
   template <int dim, int spacedim>
   void
   compute_normal_flux_constraints(
index d23919f0ddca532cf83b88015968c78e4073e652..4c112dd8c08819afe3a18ef2b59f0dd10442fa64 100644 (file)
@@ -34,6 +34,16 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
         AffineConstraints<double> &       constraints,
         const Mapping<deal_II_dimension> &mapping);
 
+      template void
+      compute_no_normal_flux_constraints_on_level(
+        const DoFHandler<deal_II_dimension> &dof_handler,
+        const MGConstrainedDoFs &            mg_constrained_dofs,
+        const unsigned int                   level,
+        const unsigned int                   first_vector_component,
+        const std::set<types::boundary_id> & boundary_ids,
+        AffineConstraints<double> &          constraints,
+        const Mapping<deal_II_dimension> &   mapping);
+
       template void
       compute_nonzero_tangential_flux_constraints(
         const DoFHandler<deal_II_dimension> &dof_handler,
diff --git a/tests/numerics/no_flux_16.cc b/tests/numerics/no_flux_16.cc
new file mode 100644 (file)
index 0000000..f9acc1f
--- /dev/null
@@ -0,0 +1,105 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// add tests for compute_no_normal_flux_constraints_on_level()
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/vector_tools_constraints.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+run(const Triangulation<dim> &          triangulation,
+    const std::set<types::boundary_id> &no_flux_boundary)
+{
+  FESystem<dim>   fe(FE_Q<dim>(1), dim);
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(fe);
+  dof_handler.distribute_mg_dofs();
+
+  MGConstrainedDoFs mg_constrained_dofs;
+  mg_constrained_dofs.initialize(dof_handler);
+
+  const unsigned int n_levels = triangulation.n_global_levels();
+
+  MappingQ<dim> mapping(4);
+
+  for (unsigned int level = 0; level < n_levels; ++level)
+    {
+      AffineConstraints<double> user_level_constraints;
+
+      VectorTools::compute_no_normal_flux_constraints_on_level(
+        dof_handler,
+        mg_constrained_dofs,
+        level,
+        0,
+        no_flux_boundary,
+        user_level_constraints,
+        mapping);
+
+      user_level_constraints.print(deallog.get_file_stream());
+
+      deallog.get_file_stream() << std::flush;
+      user_level_constraints.close();
+
+      deallog << "Level " << level << " OK" << std::endl;
+    }
+}
+
+
+int
+main()
+{
+  initlog();
+  deallog.get_file_stream().precision(7);
+  deallog.get_file_stream().setf(std::ios::fixed);
+
+  {
+    const unsigned int dim = 2;
+    Triangulation<dim> triangulation(
+      Triangulation<dim>::limit_level_difference_at_vertices);
+    GridGenerator::quarter_hyper_ball(triangulation);
+    triangulation.refine_global(1);
+    std::set<types::boundary_id> no_flux_boundary{0, 1};
+    run<dim>(triangulation, no_flux_boundary);
+  }
+  {
+    const unsigned int dim = 3;
+    Triangulation<dim> triangulation(
+      Triangulation<dim>::limit_level_difference_at_vertices);
+    GridGenerator::quarter_hyper_ball(triangulation);
+    triangulation.refine_global(1);
+    std::set<types::boundary_id> no_flux_boundary{0, 1};
+    run<dim>(triangulation, no_flux_boundary);
+  }
+}
diff --git a/tests/numerics/no_flux_16.output b/tests/numerics/no_flux_16.output
new file mode 100644 (file)
index 0000000..5dee504
--- /dev/null
@@ -0,0 +1,141 @@
+
+    0 = 0
+    1 = 0
+    3 = 0
+    4 = 0
+    8 = 0
+    9 = 0
+    11 10:  -1.0000000
+    12 = 0
+    13 = 0
+DEAL::Level 0 OK
+    0 = 0
+    1 = 0
+    3 = 0
+    4 = 0
+    9 = 0
+    12 = 0
+    18 = 0
+    19 = 0
+    20 21:  -0.4142136
+    23 = 0
+    27 26:  -1.0000000
+    30 = 0
+    31 = 0
+    32 = 0
+    35 34:  -0.4142136
+DEAL::Level 1 OK
+    0 = 0
+    1 = 0
+    2 = 0
+    4 = 0
+    5 = 0
+    6 = 0
+    8 = 0
+    11 = 0
+    12 = 0
+    13 = 0
+    16 = 0
+    18 = 0
+    24 = 0
+    25 = 0
+    26 = 0
+    28 27:  -1.0000000
+    29 = 0
+    31 = 0
+    32 30:  -1.0000000
+    35 33:  -1.0000000
+    35 34:  -1.0000000
+    36 = 0
+    37 = 0
+    38 = 0
+    39 = 0
+    41 40:  -1.0000000
+    42 = 0
+    43 = 0
+    44 = 0
+DEAL::Level 0 OK
+    0 = 0
+    1 = 0
+    2 = 0
+    4 = 0
+    5 = 0
+    6 = 0
+    8 = 0
+    11 = 0
+    12 = 0
+    13 = 0
+    16 = 0
+    18 = 0
+    25 = 0
+    26 = 0
+    29 = 0
+    31 = 0
+    36 = 0
+    38 = 0
+    41 = 0
+    42 = 0
+    50 = 0
+    54 = 0
+    55 = 0
+    58 = 0
+    60 = 0
+    67 = 0
+    72 = 0
+    81 = 0
+    82 = 0
+    83 = 0
+    84 85:  -0.4142136
+    86 = 0
+    88 = 0
+    89 = 0
+    92 = 0
+    93 95:  -0.4142136
+    94 = 0
+    96 97:  -0.4237876
+    96 98:  -0.4237876
+    100 = 0
+    106 105:  -1.0000000
+    107 = 0
+    110 = 0
+    112 111:  -1.0000000
+    112 113:  -0.4494897
+    118 = 0
+    119 117:  -1.0000000
+    122 120:  -1.0000000
+    122 121:  -0.4494897
+    124 = 0
+    131 129:  -1.0000000
+    131 130:  -1.0000000
+    135 = 0
+    136 = 0
+    137 = 0
+    138 = 0
+    140 = 0
+    142 141:  -0.4142136
+    143 = 0
+    146 = 0
+    147 = 0
+    148 149:  -0.4142136
+    150 = 0
+    154 153:  -0.4237876
+    154 155:  -0.4237876
+    159 = 0
+    161 160:  -1.0000000
+    162 = 0
+    167 165:  -0.4494897
+    167 166:  -1.0000000
+    171 = 0
+    172 = 0
+    175 = 0
+    177 = 0
+    183 = 0
+    184 = 0
+    185 = 0
+    187 = 0
+    188 186:  -0.4142136
+    189 = 0
+    191 190:  -0.4142136
+    194 192:  -0.4237876
+    194 193:  -0.4237876
+DEAL::Level 1 OK

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.