]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implemented regularize corner cells 4403/head
authorLuca Heltai <luca.heltai@sissa.it>
Thu, 18 May 2017 13:38:27 +0000 (15:38 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Mon, 22 May 2017 10:03:34 +0000 (12:03 +0200)
doc/doxygen/images/regularize_mesh_01.png [new file with mode: 0644]
doc/doxygen/images/regularize_mesh_02.png [new file with mode: 0644]
doc/doxygen/images/regularize_mesh_03.png [new file with mode: 0644]
doc/news/changes/minor/20170518LucaHeltai [new file with mode: 0644]
include/deal.II/grid/grid_tools.h
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in
tests/grid/grid_tools_regularize_01.cc [new file with mode: 0644]
tests/grid/grid_tools_regularize_01.output [new file with mode: 0644]
tests/grid/grid_tools_regularize_02.cc [new file with mode: 0644]
tests/grid/grid_tools_regularize_02.output [new file with mode: 0644]

diff --git a/doc/doxygen/images/regularize_mesh_01.png b/doc/doxygen/images/regularize_mesh_01.png
new file mode 100644 (file)
index 0000000..3e681db
Binary files /dev/null and b/doc/doxygen/images/regularize_mesh_01.png differ
diff --git a/doc/doxygen/images/regularize_mesh_02.png b/doc/doxygen/images/regularize_mesh_02.png
new file mode 100644 (file)
index 0000000..9c5b6a5
Binary files /dev/null and b/doc/doxygen/images/regularize_mesh_02.png differ
diff --git a/doc/doxygen/images/regularize_mesh_03.png b/doc/doxygen/images/regularize_mesh_03.png
new file mode 100644 (file)
index 0000000..08969ee
Binary files /dev/null and b/doc/doxygen/images/regularize_mesh_03.png differ
diff --git a/doc/news/changes/minor/20170518LucaHeltai b/doc/news/changes/minor/20170518LucaHeltai
new file mode 100644 (file)
index 0000000..30f83a2
--- /dev/null
@@ -0,0 +1,5 @@
+New: Added a GridTools::regularize_corner_cells function that 
+detects if the boundary cells of a mesh at corner positions 
+(with dim adjacent faces on the boundary) need to be split into 
+cells with smaller angles.
+<br> (Luca Heltai, Martin Kronbichler, 2017/05/18)
index 7b90b8672e70bbc3a5c2454d200909ebde07be31..37198c9744287e2af9fd16a74d8ea9218d8d53e3 100644 (file)
@@ -461,6 +461,99 @@ namespace GridTools
                       const double max_ratio = 1.6180339887,
                       const unsigned int max_iterations = 5);
 
+  /**
+   * Analyze the boundary cells of a mesh, and if one cell is found at
+   * a corner position (with dim adjacent faces on the boundary), and its
+   * dim-dimensional angle fraction exceeds @p limit_angle_fraction,
+   * refine globally once, and replace the children of such cell
+   * with children where the corner is no longer offending the given angle
+   * fraction.
+   *
+   * If no boundary cells exist with two adjacent faces on the boundary, then
+   * the triangulation is left untouched. If instead we do have cells with dim
+   * adjacent faces on the boundary, then the fraction between the dim-dimensional
+   * solid angle and dim*pi/2 is checked against the parameter @p limit_angle_fraction.
+   * If it is higher, the grid is refined once, and the children of the
+   * offending cell are replaced with some cells that instead respect the limit. After
+   * this process the triangulation is flattened, and all Manifold objects are restored
+   * as they were in the original triangulation.
+   *
+   * An example is given by the following mesh, obtained by attaching a SphericalManifold
+   * to a mesh generated using GridGenerator::hyper_cube:
+   *
+   * @code
+   * const SphericalManifold<dim> m0;
+   * Triangulation<dim> tria;
+   * GridGenerator::hyper_cube(tria,-1,1);
+   * tria.set_all_manifold_ids_on_boundary(0);
+   * tria.set_manifold(0, m0);
+   * tria.refine_global(4);
+   * @endcode
+   *
+   * <p ALIGN="center">
+   * @image html regularize_mesh_01.png
+   * </p>
+   *
+   * The four cells that were originally the corners of a square will give you some troubles
+   * during computations, as the jacobian of the transformation from the reference cell to
+   * those cells will go to zero, affecting the error constants of the finite element estimates.
+   *
+   * Those cells have a corner with an angle that is very close to 180 degrees, i.e., an angle
+   * fraction very close to one.
+   *
+   * The same code, adding a call to regularize_corner_cells:
+   * @code
+   * const SphericalManifold<dim> m0;
+   * Triangulation<dim> tria;
+   * GridGenerator::hyper_cube(tria,-1,1);
+   * tria.set_all_manifold_ids_on_boundary(0);
+   * tria.set_manifold(0, m0);
+   * GridTools::regularize_corner_cells(tria);
+   * tria.refine_global(2);
+   * @endcode
+   * generates a mesh that has a much better behaviour w.r.t. the jacobian of the Mapping:
+   *
+   * <p ALIGN="center">
+   * @image html regularize_mesh_02.png
+   * </p>
+   *
+   * This mesh is very similar to the one obtained by GridGenerator::hyper_ball. However, using
+   * GridTools::regularize_corner_cells one has the freedom to choose when to apply the
+   * regularization, i.e., one could in principle first refine a few times, and then call the
+   * regularize_corner_cells function:
+   *
+   * @code
+   * const SphericalManifold<dim> m0;
+   * Triangulation<dim> tria;
+   * GridGenerator::hyper_cube(tria,-1,1);
+   * tria.set_all_manifold_ids_on_boundary(0);
+   * tria.set_manifold(0, m0);
+   * tria.refine_global(2);
+   * GridTools::regularize_corner_cells(tria);
+   * tria.refine_global(1);
+   * @endcode
+   *
+   * This generates the following mesh:
+   *
+   * <p ALIGN="center">
+   * @image html regularize_mesh_03.png
+   * </p>
+   *
+   * The function is currently implemented only for dim = 2 and
+   * will throw an exception if called with dim = 3.
+   *
+   * @param[in,out] tria Triangulation to regularize.
+   *
+   * @param[in] limit_angle_fraction Maximum ratio of angle or solid
+   * angle that is allowed for a corner element in the mesh.
+   *
+   * @author Luca Heltai, Martin Kronbichler, 2017
+   */
+  template <int dim, int spacedim>
+  void
+  regularize_corner_cells(Triangulation<dim,spacedim> &tria,
+                          const double limit_angle_fraction=.75);
+
   /*@}*/
   /**
    * @name Finding cells and vertices of a triangulation
index d4e54d2de6b792525aa26bfd12338769915a4fd6..b9616c22dea5bf121663cf708b0b527d55be2a5e 100644 (file)
@@ -35,6 +35,7 @@
 #include <deal.II/grid/tria_boundary.h>
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_reordering.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/dofs/dof_tools.h>
@@ -4451,6 +4452,312 @@ next_cell:
       }
   }
 
+
+  template<int dim, int spacedim>
+  void regularize_corner_cells (Triangulation<dim,spacedim> &tria,
+                                const double limit_angle_fraction)
+  {
+    if (dim == 1)
+      return; // Nothing to do
+
+    // Check that we don't have hanging nodes
+    AssertThrow(!tria.has_hanging_nodes(), ExcMessage("The input Triangulation cannot "
+                                                      "have hanging nodes."));
+
+
+    bool has_cells_with_more_than_dim_faces_on_boundary = true;
+    bool has_cells_with_dim_faces_on_boundary = false;
+
+    unsigned int refinement_cycles = 0;
+
+    while (has_cells_with_more_than_dim_faces_on_boundary)
+      {
+        has_cells_with_more_than_dim_faces_on_boundary = false;
+
+        for (auto cell: tria.active_cell_iterators())
+          {
+            unsigned int boundary_face_counter = 0;
+            for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+              if (cell->face(f)->at_boundary())
+                boundary_face_counter++;
+            if (boundary_face_counter > dim)
+              {
+                has_cells_with_more_than_dim_faces_on_boundary = true;
+                break;
+              }
+            else if (boundary_face_counter == dim)
+              has_cells_with_dim_faces_on_boundary = true;
+          }
+        if (has_cells_with_more_than_dim_faces_on_boundary)
+          {
+            tria.refine_global(1);
+            refinement_cycles++;
+          }
+      }
+
+    if (has_cells_with_dim_faces_on_boundary)
+      {
+        tria.refine_global(1);
+        refinement_cycles++;
+      }
+    else
+      {
+        while (refinement_cycles>0)
+          {
+            for (auto cell: tria.active_cell_iterators())
+              cell->set_coarsen_flag();
+            tria.execute_coarsening_and_refinement();
+            refinement_cycles--;
+          }
+        return;
+      }
+
+    std::vector<bool> cells_to_remove(tria.n_active_cells(), false);
+    std::vector<Point<spacedim> > vertices = tria.get_vertices();
+
+    std::vector<bool> faces_to_remove(tria.n_raw_faces(),false);
+
+    std::vector<CellData<dim> > cells_to_add;
+    SubCellData                 subcelldata_to_add;
+
+    // Trick compiler for dimension independent things
+    const unsigned int
+    v0 = 0, v1 = 1,
+    v2 = (dim > 1 ? 2:0), v3 = (dim > 1 ? 3:0),
+    v4 = (dim > 2 ? 4:0), v5 = (dim > 2 ? 5:0),
+    v6 = (dim > 2 ? 6:0), v7 = (dim > 2 ? 7:0);
+
+    for (auto cell : tria.active_cell_iterators())
+      {
+        double angle_fraction = 0;
+        unsigned int vertex_at_corner = numbers::invalid_unsigned_int;
+
+        if (dim == 2)
+          {
+            Tensor<1,spacedim> p0;
+            p0[spacedim > 1 ? 1 : 0] = 1;
+            Tensor<1,spacedim> p1;
+            p1[0] = 1;
+
+            if (cell->face(v0)->at_boundary() && cell->face(v3)->at_boundary())
+              {
+                p0 = cell->vertex(v0) -  cell->vertex(v2);
+                p1 = cell->vertex(v3) -  cell->vertex(v2);
+                vertex_at_corner = v2;
+              }
+            else  if (cell->face(v3)->at_boundary() && cell->face(v1)->at_boundary())
+              {
+                p0 = cell->vertex(v2) -  cell->vertex(v3);
+                p1 = cell->vertex(v1) -  cell->vertex(v3);
+                vertex_at_corner = v3;
+              }
+            else  if (cell->face(1)->at_boundary() && cell->face(2)->at_boundary())
+              {
+                p0 = cell->vertex(v0) -  cell->vertex(v1);
+                p1 = cell->vertex(v3) -  cell->vertex(v1);
+                vertex_at_corner = v1;
+              }
+            else  if (cell->face(2)->at_boundary() && cell->face(0)->at_boundary())
+              {
+                p0 = cell->vertex(v2) -  cell->vertex(v0);
+                p1 = cell->vertex(v1) -  cell->vertex(v0);
+                vertex_at_corner = v0;
+              }
+            p0 /= p0.norm();
+            p1 /= p1.norm();
+            angle_fraction = std::acos(p0*p1)/numbers::PI;
+
+          }
+        else
+          {
+            Assert(false, ExcNotImplemented());
+          }
+
+        if (angle_fraction > limit_angle_fraction)
+          {
+
+            auto flags_removal = [&](unsigned int f1, unsigned int f2,
+                                     unsigned int n1, unsigned int n2) -> void
+            {
+              cells_to_remove[cell->active_cell_index()] = true;
+              cells_to_remove[cell->neighbor(n1)->active_cell_index()] = true;
+              cells_to_remove[cell->neighbor(n2)->active_cell_index()] = true;
+
+              faces_to_remove[cell->face(f1)->index()] = true;
+              faces_to_remove[cell->face(f2)->index()] = true;
+
+              faces_to_remove[cell->neighbor(n1)->face(f1)->index()] = true;
+              faces_to_remove[cell->neighbor(n2)->face(f2)->index()] = true;
+            };
+
+            auto cell_creation = [&](
+                                   const unsigned int vv0,
+                                   const unsigned int vv1,
+                                   const unsigned int f0,
+                                   const unsigned int f1,
+
+                                   const unsigned int n0,
+                                   const unsigned int v0n0,
+                                   const unsigned int v1n0,
+
+                                   const unsigned int n1,
+                                   const unsigned int v0n1,
+                                   const unsigned int v1n1)
+            {
+              CellData<dim> c1, c2;
+              CellData<1> l1, l2;
+
+              c1.vertices[v0] = cell->vertex_index(vv0);
+              c1.vertices[v1] = cell->vertex_index(vv1);
+              c1.vertices[v2] = cell->neighbor(n0)->vertex_index(v0n0);
+              c1.vertices[v3] = cell->neighbor(n0)->vertex_index(v1n0);
+
+              c1.manifold_id = cell->manifold_id();
+              c1.material_id = cell->material_id();
+
+              c2.vertices[v0] = cell->vertex_index(vv0);
+              c2.vertices[v1] = cell->neighbor(n1)->vertex_index(v0n1);
+              c2.vertices[v2] = cell->vertex_index(vv1);
+              c2.vertices[v3] = cell->neighbor(n1)->vertex_index(v1n1);
+
+              c2.manifold_id = cell->manifold_id();
+              c2.material_id = cell->material_id();
+
+              l1.vertices[0] = cell->vertex_index(vv0);
+              l1.vertices[1] = cell->neighbor(n0)->vertex_index(v0n0);
+
+              l1.boundary_id = cell->line(f0)->boundary_id();
+              l1.manifold_id = cell->line(f0)->manifold_id();
+              subcelldata_to_add.boundary_lines.push_back(l1);
+
+              l2.vertices[0] = cell->vertex_index(vv0);
+              l2.vertices[1] = cell->neighbor(n1)->vertex_index(v0n1);
+
+              l2.boundary_id = cell->line(f1)->boundary_id();
+              l2.manifold_id = cell->line(f1)->manifold_id();
+              subcelldata_to_add.boundary_lines.push_back(l2);
+
+              cells_to_add.push_back(c1);
+              cells_to_add.push_back(c2);
+            };
+
+            if (dim == 2)
+              {
+                switch (vertex_at_corner)
+                  {
+                  case 0:
+                    flags_removal(0,2,3,1);
+                    cell_creation(0,3, 0,2, 3,2,3,  1,1,3);
+                    break;
+                  case 1:
+                    flags_removal(1,2,3,0);
+                    cell_creation(1,2, 2,1, 0,0,2, 3,3,2);
+                    break;
+                  case 2:
+                    flags_removal(3,0,1,2);
+                    cell_creation(2,1, 3,0, 1,3,1, 2,0,1);
+                    break;
+                  case 3:
+                    flags_removal(3,1,0,2);
+                    cell_creation(3,0, 1,3, 2,1,0, 0,2,0);
+                    break;
+                  }
+              }
+            else
+              {
+                Assert(false, ExcNotImplemented());
+              }
+          }
+      }
+
+    // if no cells need to be added, then no regularization is necessary. Restore things
+    // as they were before this function was called.
+    if (cells_to_add.size() == 0)
+      {
+        while (refinement_cycles>0)
+          {
+            for (auto cell: tria.active_cell_iterators())
+              cell->set_coarsen_flag();
+            tria.execute_coarsening_and_refinement();
+            refinement_cycles--;
+          }
+        return;
+      }
+
+    // add the cells that were not marked as skipped
+    for (auto cell : tria.active_cell_iterators())
+      {
+        if (cells_to_remove[cell->active_cell_index()] == false)
+          {
+            CellData<dim> c;
+            for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+              c.vertices[v] = cell->vertex_index(v);
+            c.manifold_id = cell->manifold_id();
+            c.material_id = cell->material_id();
+            cells_to_add.push_back(c);
+          }
+      }
+
+    // Face counter for both dim == 2 and dim == 3
+    typename Triangulation<dim,spacedim>::active_face_iterator
+    face = tria.begin_active_face(),
+    endf = tria.end_face();
+    for (; face != endf; ++face)
+      if ( (face->at_boundary() || face->manifold_id() != numbers::invalid_manifold_id)
+           && faces_to_remove[face->index()] == false)
+        {
+          for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
+            {
+              CellData<1> line;
+              if (dim == 2)
+                {
+                  for (unsigned int v=0; v<GeometryInfo<1>::vertices_per_cell; ++v)
+                    line.vertices[v] = face->vertex_index(v);
+                  line.boundary_id = face->boundary_id();
+                  line.manifold_id = face->manifold_id();
+                }
+              else
+                {
+                  for (unsigned int v=0; v<GeometryInfo<1>::vertices_per_cell; ++v)
+                    line.vertices[v] = face->line(l)->vertex_index(v);
+                  line.boundary_id = face->line(l)->boundary_id();
+                  line.manifold_id = face->line(l)->manifold_id();
+                }
+              subcelldata_to_add.boundary_lines.push_back(line);
+            }
+          if (dim == 3)
+            {
+              CellData<2> quad;
+              for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
+                quad.vertices[v] = face->vertex_index(v);
+              quad.boundary_id = face->boundary_id();
+              quad.manifold_id = face->manifold_id();
+              subcelldata_to_add.boundary_quads.push_back(quad);
+            }
+        }
+    GridTools::delete_unused_vertices(vertices, cells_to_add, subcelldata_to_add);
+    GridReordering<dim,spacedim>::reorder_cells(cells_to_add, true);
+
+    // Save manifolds
+    auto manifold_ids = tria.get_manifold_ids();
+    std::map<types::manifold_id, const Manifold<dim,spacedim>*> manifolds;
+    // Set manifolds in new Triangulation
+    for (auto manifold_id: manifold_ids)
+      if (manifold_id != numbers::invalid_manifold_id)
+        manifolds[manifold_id] = &tria.get_manifold(manifold_id);
+
+    tria.clear();
+
+    tria.create_triangulation(vertices, cells_to_add, subcelldata_to_add);
+
+    // Restore manifolds
+    for (auto manifold_id: manifold_ids)
+      if (manifold_id != numbers::invalid_manifold_id)
+        tria.set_manifold(manifold_id, *manifolds[manifold_id]);
+  }
+
+
+
 } /* namespace GridTools */
 
 
index 1b593ec313ff772422b68f754249a2c08cad3a43..f148328ef1421278ac2b80a11544b2b1757e44a8 100644 (file)
@@ -411,6 +411,9 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS
     template void copy_material_to_manifold_id<deal_II_dimension, deal_II_space_dimension>
     (Triangulation<deal_II_dimension, deal_II_space_dimension> &, const bool);
 
+    template
+    void regularize_corner_cells
+    (Triangulation<deal_II_dimension,deal_II_space_dimension> &, double);
     \}
 #endif
 }
diff --git a/tests/grid/grid_tools_regularize_01.cc b/tests/grid/grid_tools_regularize_01.cc
new file mode 100644 (file)
index 0000000..743077a
--- /dev/null
@@ -0,0 +1,50 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// GridTools::regularize_corner_cells
+
+#include "../tests.h"
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+
+template <int dim>
+void test ()
+{
+  const SphericalManifold<dim> m0;
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria,-1,1);
+  tria.set_all_manifold_ids_on_boundary(0);
+  tria.set_manifold(0, m0);
+
+  GridTools::regularize_corner_cells (tria);
+
+  GridOut grid_out;
+  grid_out.write_msh (tria, deallog.get_file_stream());
+}
+
+
+
+int main ()
+{
+  initlog();
+
+  test<2> ();
+
+  return 0;
+}
+
diff --git a/tests/grid/grid_tools_regularize_01.output b/tests/grid/grid_tools_regularize_01.output
new file mode 100644 (file)
index 0000000..214914d
--- /dev/null
@@ -0,0 +1,35 @@
+$NOD
+17
+1  -1.00000 -1.00000 0
+2  1.00000 -1.00000 0
+3  -1.00000 1.00000 0
+4  1.00000 1.00000 0
+5  0.00000 -1.41421 0
+6  -1.41421 0.00000 0
+7  1.41421 0.00000 0
+8  0.00000 1.41421 0
+9  0.00000 0.00000 0
+10  0.00000 -0.707107 0
+11  0.00000 0.707107 0
+12  -0.707107 0.00000 0
+13  0.707107 0.00000 0
+14  -0.621135 -0.621135 0
+15  0.621135 -0.621135 0
+16  -0.621135 0.621135 0
+17  0.621135 0.621135 0
+$ENDNOD
+$ELM
+12
+1 3 0 0 4 12 6 1 14 
+2 3 0 0 4 10 14 1 5 
+3 3 0 0 4 10 5 2 15 
+4 3 0 0 4 13 15 2 7 
+5 3 0 0 4 11 8 3 16 
+6 3 0 0 4 12 16 3 6 
+7 3 0 0 4 13 7 4 17 
+8 3 0 0 4 11 17 4 8 
+9 3 0 0 4 9 12 14 10 
+10 3 0 0 4 9 10 15 13 
+11 3 0 0 4 9 11 16 12 
+12 3 0 0 4 9 13 17 11 
+$ENDELM
diff --git a/tests/grid/grid_tools_regularize_02.cc b/tests/grid/grid_tools_regularize_02.cc
new file mode 100644 (file)
index 0000000..0398167
--- /dev/null
@@ -0,0 +1,55 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// GridTools::regularize_corner_cells on more complicated mesh
+
+#include "../tests.h"
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+
+int main()
+{
+  initlog();
+
+  Point<2> p0(0,0), p1(4,1);
+  Point<2> c0(.1, .5), c1(3.9,.5);
+
+  SphericalManifold<2> m0(c0);
+  SphericalManifold<2> m1(c1);
+
+  Triangulation<2> tria;
+  std::vector<unsigned int> subdivisions(2);
+  subdivisions[0] = 4;
+  subdivisions[1] = 1;
+  GridGenerator::subdivided_hyper_rectangle(tria,subdivisions, p0,p1,true);
+
+  GridTools::copy_boundary_to_manifold_id(tria);
+
+  tria.set_manifold(0, m0);
+  tria.set_manifold(1, m1);
+
+  GridTools::regularize_corner_cells (tria);
+  tria.refine_global(1);
+
+  GridOut grid_out;
+  grid_out.write_msh (tria, deallog.get_file_stream());
+
+  return 0;
+}
+
+
diff --git a/tests/grid/grid_tools_regularize_02.output b/tests/grid/grid_tools_regularize_02.output
new file mode 100644 (file)
index 0000000..5d2bb06
--- /dev/null
@@ -0,0 +1,520 @@
+
+$NOD
+273
+1  0.00000 0.00000 0
+2  1.00000 0.00000 0
+3  2.00000 0.00000 0
+4  3.00000 0.00000 0
+5  4.00000 0.00000 0
+6  0.00000 1.00000 0
+7  1.00000 1.00000 0
+8  2.00000 1.00000 0
+9  3.00000 1.00000 0
+10  4.00000 1.00000 0
+11  0.500000 0.00000 0
+12  -0.409902 0.500000 0
+13  1.50000 0.00000 0
+14  1.00000 0.500000 0
+15  2.50000 0.00000 0
+16  2.00000 0.500000 0
+17  3.50000 0.00000 0
+18  3.00000 0.500000 0
+19  4.40990 0.500000 0
+20  0.500000 1.00000 0
+21  1.50000 1.00000 0
+22  2.50000 1.00000 0
+23  3.50000 1.00000 0
+24  0.448762 0.500000 0
+25  1.50000 0.500000 0
+26  2.50000 0.500000 0
+27  3.55124 0.500000 0
+28  0.750000 0.00000 0
+29  1.25000 0.00000 0
+30  1.75000 0.00000 0
+31  1.00000 0.250000 0
+32  1.00000 0.750000 0
+33  2.25000 0.00000 0
+34  2.75000 0.00000 0
+35  2.00000 0.250000 0
+36  2.00000 0.750000 0
+37  3.25000 0.00000 0
+38  3.00000 0.250000 0
+39  3.00000 0.750000 0
+40  0.750000 1.00000 0
+41  1.25000 1.00000 0
+42  1.75000 1.00000 0
+43  2.25000 1.00000 0
+44  2.75000 1.00000 0
+45  3.25000 1.00000 0
+46  0.474381 0.250000 0
+47  0.474381 0.750000 0
+48  0.0194302 0.500000 0
+49  0.724381 0.500000 0
+50  1.50000 0.250000 0
+51  1.50000 0.750000 0
+52  1.25000 0.500000 0
+53  1.75000 0.500000 0
+54  2.50000 0.250000 0
+55  2.50000 0.750000 0
+56  2.25000 0.500000 0
+57  2.75000 0.500000 0
+58  3.52562 0.250000 0
+59  3.52562 0.750000 0
+60  3.27562 0.500000 0
+61  3.98057 0.500000 0
+62  0.123543 0.240841 0
+63  0.737191 0.250000 0
+64  0.123543 0.759159 0
+65  0.737191 0.750000 0
+66  1.25000 0.250000 0
+67  1.75000 0.250000 0
+68  1.25000 0.750000 0
+69  1.75000 0.750000 0
+70  2.25000 0.250000 0
+71  2.75000 0.250000 0
+72  2.25000 0.750000 0
+73  2.75000 0.750000 0
+74  3.26281 0.250000 0
+75  3.87646 0.240841 0
+76  3.26281 0.750000 0
+77  3.87646 0.759159 0
+78  0.875000 0.00000 0
+79  1.87500 0.00000 0
+80  2.87500 0.00000 0
+81  0.875000 1.00000 0
+82  1.87500 1.00000 0
+83  2.87500 1.00000 0
+84  0.250000 0.00000 0
+85  -0.294329 0.176728 0
+86  -0.294329 0.823272 0
+87  1.37500 0.00000 0
+88  1.00000 0.375000 0
+89  1.00000 0.625000 0
+90  0.862191 0.500000 0
+91  2.37500 0.00000 0
+92  2.00000 0.375000 0
+93  2.00000 0.625000 0
+94  1.87500 0.500000 0
+95  3.75000 0.00000 0
+96  3.37500 0.00000 0
+97  3.00000 0.375000 0
+98  3.00000 0.625000 0
+99  2.87500 0.500000 0
+100  4.29433 0.176728 0
+101  4.29433 0.823272 0
+102  0.250000 1.00000 0
+103  1.37500 1.00000 0
+104  2.37500 1.00000 0
+105  3.75000 1.00000 0
+106  3.37500 1.00000 0
+107  0.461572 0.375000 0
+108  0.461572 0.625000 0
+109  0.234096 0.500000 0
+110  1.50000 0.375000 0
+111  1.50000 0.625000 0
+112  1.37500 0.500000 0
+113  2.50000 0.375000 0
+114  2.50000 0.625000 0
+115  2.37500 0.500000 0
+116  3.53843 0.375000 0
+117  3.53843 0.625000 0
+118  3.41343 0.500000 0
+119  3.76590 0.500000 0
+120  0.625000 0.00000 0
+121  1.12500 0.00000 0
+122  1.62500 0.00000 0
+123  1.00000 0.125000 0
+124  0.868595 0.250000 0
+125  1.00000 0.875000 0
+126  0.868595 0.750000 0
+127  2.12500 0.00000 0
+128  2.62500 0.00000 0
+129  2.00000 0.125000 0
+130  1.87500 0.250000 0
+131  2.00000 0.875000 0
+132  1.87500 0.750000 0
+133  3.12500 0.00000 0
+134  3.00000 0.125000 0
+135  2.87500 0.250000 0
+136  3.00000 0.875000 0
+137  2.87500 0.750000 0
+138  0.625000 1.00000 0
+139  1.12500 1.00000 0
+140  1.62500 1.00000 0
+141  2.12500 1.00000 0
+142  2.62500 1.00000 0
+143  3.12500 1.00000 0
+144  0.487191 0.125000 0
+145  0.298962 0.245420 0
+146  0.487191 0.875000 0
+147  0.298962 0.754580 0
+148  -0.195236 0.500000 0
+149  0.0714865 0.370420 0
+150  0.0714865 0.629580 0
+151  0.586572 0.500000 0
+152  0.730786 0.375000 0
+153  0.730786 0.625000 0
+154  1.50000 0.125000 0
+155  1.37500 0.250000 0
+156  1.50000 0.875000 0
+157  1.37500 0.750000 0
+158  1.12500 0.500000 0
+159  1.25000 0.375000 0
+160  1.25000 0.625000 0
+161  1.62500 0.500000 0
+162  1.75000 0.375000 0
+163  1.75000 0.625000 0
+164  2.50000 0.125000 0
+165  2.37500 0.250000 0
+166  2.50000 0.875000 0
+167  2.37500 0.750000 0
+168  2.12500 0.500000 0
+169  2.25000 0.375000 0
+170  2.25000 0.625000 0
+171  2.62500 0.500000 0
+172  2.75000 0.375000 0
+173  2.75000 0.625000 0
+174  3.51281 0.125000 0
+175  3.39421 0.250000 0
+176  3.70104 0.245420 0
+177  3.51281 0.875000 0
+178  3.39421 0.750000 0
+179  3.70104 0.754580 0
+180  3.13781 0.500000 0
+181  3.26921 0.375000 0
+182  3.26921 0.625000 0
+183  4.19524 0.500000 0
+184  3.92851 0.370420 0
+185  3.92851 0.629580 0
+186  0.0617714 0.120420 0
+187  0.743595 0.125000 0
+188  0.605786 0.250000 0
+189  0.0617714 0.879580 0
+190  0.743595 0.875000 0
+191  0.605786 0.750000 0
+192  1.25000 0.125000 0
+193  1.12500 0.250000 0
+194  1.75000 0.125000 0
+195  1.62500 0.250000 0
+196  1.12500 0.750000 0
+197  1.25000 0.875000 0
+198  1.75000 0.875000 0
+199  1.62500 0.750000 0
+200  2.25000 0.125000 0
+201  2.12500 0.250000 0
+202  2.75000 0.125000 0
+203  2.62500 0.250000 0
+204  2.12500 0.750000 0
+205  2.25000 0.875000 0
+206  2.75000 0.875000 0
+207  2.62500 0.750000 0
+208  3.25640 0.125000 0
+209  3.13140 0.250000 0
+210  3.93823 0.120420 0
+211  3.13140 0.750000 0
+212  3.25640 0.875000 0
+213  3.93823 0.879580 0
+214  -0.111421 0.273574 0
+215  0.274481 0.122710 0
+216  0.274481 0.877290 0
+217  -0.111421 0.726426 0
+218  3.72552 0.122710 0
+219  4.11142 0.273574 0
+220  4.11142 0.726426 0
+221  3.72552 0.877290 0
+222  0.266529 0.372710 0
+223  0.615393 0.125000 0
+224  0.871798 0.125000 0
+225  0.596179 0.375000 0
+226  0.865393 0.375000 0
+227  0.266529 0.627290 0
+228  0.596179 0.625000 0
+229  0.865393 0.625000 0
+230  0.615393 0.875000 0
+231  0.871798 0.875000 0
+232  1.12500 0.125000 0
+233  1.37500 0.125000 0
+234  1.12500 0.375000 0
+235  1.37500 0.375000 0
+236  1.62500 0.125000 0
+237  1.87500 0.125000 0
+238  1.62500 0.375000 0
+239  1.87500 0.375000 0
+240  1.12500 0.625000 0
+241  1.37500 0.625000 0
+242  1.12500 0.875000 0
+243  1.37500 0.875000 0
+244  1.62500 0.625000 0
+245  1.87500 0.625000 0
+246  1.62500 0.875000 0
+247  1.87500 0.875000 0
+248  2.12500 0.125000 0
+249  2.37500 0.125000 0
+250  2.12500 0.375000 0
+251  2.37500 0.375000 0
+252  2.62500 0.125000 0
+253  2.87500 0.125000 0
+254  2.62500 0.375000 0
+255  2.87500 0.375000 0
+256  2.12500 0.625000 0
+257  2.37500 0.625000 0
+258  2.12500 0.875000 0
+259  2.37500 0.875000 0
+260  2.62500 0.625000 0
+261  2.87500 0.625000 0
+262  2.62500 0.875000 0
+263  2.87500 0.875000 0
+264  3.12820 0.125000 0
+265  3.38461 0.125000 0
+266  3.13461 0.375000 0
+267  3.40382 0.375000 0
+268  3.73347 0.372710 0
+269  3.13461 0.625000 0
+270  3.40382 0.625000 0
+271  3.12820 0.875000 0
+272  3.38461 0.875000 0
+273  3.73347 0.627290 0
+$ENDNOD
+$ELM
+240
+1 3 3 0 4 48 148 214 149 
+2 3 3 0 4 148 12 85 214 
+3 3 3 0 4 149 214 186 62 
+4 3 3 0 4 214 85 1 186 
+5 3 3 0 4 46 145 215 144 
+6 3 3 0 4 145 62 186 215 
+7 3 3 0 4 144 215 84 11 
+8 3 3 0 4 215 186 1 84 
+9 3 3 0 4 47 146 216 147 
+10 3 3 0 4 146 20 102 216 
+11 3 3 0 4 147 216 189 64 
+12 3 3 0 4 216 102 6 189 
+13 3 3 0 4 48 150 217 148 
+14 3 3 0 4 150 64 189 217 
+15 3 3 0 4 148 217 86 12 
+16 3 3 0 4 217 189 6 86 
+17 3 3 0 4 58 174 218 176 
+18 3 3 0 4 174 17 95 218 
+19 3 3 0 4 176 218 210 75 
+20 3 3 0 4 218 95 5 210 
+21 3 3 0 4 61 184 219 183 
+22 3 3 0 4 184 75 210 219 
+23 3 3 0 4 183 219 100 19 
+24 3 3 0 4 219 210 5 100 
+25 3 3 0 4 61 183 220 185 
+26 3 3 0 4 183 19 101 220 
+27 3 3 0 4 185 220 213 77 
+28 3 3 0 4 220 101 10 213 
+29 3 3 0 4 59 179 221 177 
+30 3 3 0 4 179 77 213 221 
+31 3 3 0 4 177 221 105 23 
+32 3 3 0 4 221 213 10 105 
+33 3 3 0 4 24 109 222 107 
+34 3 3 0 4 109 48 149 222 
+35 3 3 0 4 107 222 145 46 
+36 3 3 0 4 222 149 62 145 
+37 3 3 0 4 63 188 223 187 
+38 3 3 0 4 188 46 144 223 
+39 3 3 0 4 187 223 120 28 
+40 3 3 0 4 223 144 11 120 
+41 3 3 0 4 31 124 224 123 
+42 3 3 0 4 124 63 187 224 
+43 3 3 0 4 123 224 78 2 
+44 3 3 0 4 224 187 28 78 
+45 3 3 0 4 49 151 225 152 
+46 3 3 0 4 151 24 107 225 
+47 3 3 0 4 152 225 188 63 
+48 3 3 0 4 225 107 46 188 
+49 3 3 0 4 14 90 226 88 
+50 3 3 0 4 90 49 152 226 
+51 3 3 0 4 88 226 124 31 
+52 3 3 0 4 226 152 63 124 
+53 3 3 0 4 24 108 227 109 
+54 3 3 0 4 108 47 147 227 
+55 3 3 0 4 109 227 150 48 
+56 3 3 0 4 227 147 64 150 
+57 3 3 0 4 49 153 228 151 
+58 3 3 0 4 153 65 191 228 
+59 3 3 0 4 151 228 108 24 
+60 3 3 0 4 228 191 47 108 
+61 3 3 0 4 14 89 229 90 
+62 3 3 0 4 89 32 126 229 
+63 3 3 0 4 90 229 153 49 
+64 3 3 0 4 229 126 65 153 
+65 3 3 0 4 65 190 230 191 
+66 3 3 0 4 190 40 138 230 
+67 3 3 0 4 191 230 146 47 
+68 3 3 0 4 230 138 20 146 
+69 3 3 0 4 32 125 231 126 
+70 3 3 0 4 125 7 81 231 
+71 3 3 0 4 126 231 190 65 
+72 3 3 0 4 231 81 40 190 
+73 3 3 0 4 66 193 232 192 
+74 3 3 0 4 193 31 123 232 
+75 3 3 0 4 192 232 121 29 
+76 3 3 0 4 232 123 2 121 
+77 3 3 0 4 50 155 233 154 
+78 3 3 0 4 155 66 192 233 
+79 3 3 0 4 154 233 87 13 
+80 3 3 0 4 233 192 29 87 
+81 3 3 0 4 52 158 234 159 
+82 3 3 0 4 158 14 88 234 
+83 3 3 0 4 159 234 193 66 
+84 3 3 0 4 234 88 31 193 
+85 3 3 0 4 25 112 235 110 
+86 3 3 0 4 112 52 159 235 
+87 3 3 0 4 110 235 155 50 
+88 3 3 0 4 235 159 66 155 
+89 3 3 0 4 67 195 236 194 
+90 3 3 0 4 195 50 154 236 
+91 3 3 0 4 194 236 122 30 
+92 3 3 0 4 236 154 13 122 
+93 3 3 0 4 35 130 237 129 
+94 3 3 0 4 130 67 194 237 
+95 3 3 0 4 129 237 79 3 
+96 3 3 0 4 237 194 30 79 
+97 3 3 0 4 53 161 238 162 
+98 3 3 0 4 161 25 110 238 
+99 3 3 0 4 162 238 195 67 
+100 3 3 0 4 238 110 50 195 
+101 3 3 0 4 16 94 239 92 
+102 3 3 0 4 94 53 162 239 
+103 3 3 0 4 92 239 130 35 
+104 3 3 0 4 239 162 67 130 
+105 3 3 0 4 52 160 240 158 
+106 3 3 0 4 160 68 196 240 
+107 3 3 0 4 158 240 89 14 
+108 3 3 0 4 240 196 32 89 
+109 3 3 0 4 25 111 241 112 
+110 3 3 0 4 111 51 157 241 
+111 3 3 0 4 112 241 160 52 
+112 3 3 0 4 241 157 68 160 
+113 3 3 0 4 68 197 242 196 
+114 3 3 0 4 197 41 139 242 
+115 3 3 0 4 196 242 125 32 
+116 3 3 0 4 242 139 7 125 
+117 3 3 0 4 51 156 243 157 
+118 3 3 0 4 156 21 103 243 
+119 3 3 0 4 157 243 197 68 
+120 3 3 0 4 243 103 41 197 
+121 3 3 0 4 53 163 244 161 
+122 3 3 0 4 163 69 199 244 
+123 3 3 0 4 161 244 111 25 
+124 3 3 0 4 244 199 51 111 
+125 3 3 0 4 16 93 245 94 
+126 3 3 0 4 93 36 132 245 
+127 3 3 0 4 94 245 163 53 
+128 3 3 0 4 245 132 69 163 
+129 3 3 0 4 69 198 246 199 
+130 3 3 0 4 198 42 140 246 
+131 3 3 0 4 199 246 156 51 
+132 3 3 0 4 246 140 21 156 
+133 3 3 0 4 36 131 247 132 
+134 3 3 0 4 131 8 82 247 
+135 3 3 0 4 132 247 198 69 
+136 3 3 0 4 247 82 42 198 
+137 3 3 0 4 70 201 248 200 
+138 3 3 0 4 201 35 129 248 
+139 3 3 0 4 200 248 127 33 
+140 3 3 0 4 248 129 3 127 
+141 3 3 0 4 54 165 249 164 
+142 3 3 0 4 165 70 200 249 
+143 3 3 0 4 164 249 91 15 
+144 3 3 0 4 249 200 33 91 
+145 3 3 0 4 56 168 250 169 
+146 3 3 0 4 168 16 92 250 
+147 3 3 0 4 169 250 201 70 
+148 3 3 0 4 250 92 35 201 
+149 3 3 0 4 26 115 251 113 
+150 3 3 0 4 115 56 169 251 
+151 3 3 0 4 113 251 165 54 
+152 3 3 0 4 251 169 70 165 
+153 3 3 0 4 71 203 252 202 
+154 3 3 0 4 203 54 164 252 
+155 3 3 0 4 202 252 128 34 
+156 3 3 0 4 252 164 15 128 
+157 3 3 0 4 38 135 253 134 
+158 3 3 0 4 135 71 202 253 
+159 3 3 0 4 134 253 80 4 
+160 3 3 0 4 253 202 34 80 
+161 3 3 0 4 57 171 254 172 
+162 3 3 0 4 171 26 113 254 
+163 3 3 0 4 172 254 203 71 
+164 3 3 0 4 254 113 54 203 
+165 3 3 0 4 18 99 255 97 
+166 3 3 0 4 99 57 172 255 
+167 3 3 0 4 97 255 135 38 
+168 3 3 0 4 255 172 71 135 
+169 3 3 0 4 56 170 256 168 
+170 3 3 0 4 170 72 204 256 
+171 3 3 0 4 168 256 93 16 
+172 3 3 0 4 256 204 36 93 
+173 3 3 0 4 26 114 257 115 
+174 3 3 0 4 114 55 167 257 
+175 3 3 0 4 115 257 170 56 
+176 3 3 0 4 257 167 72 170 
+177 3 3 0 4 72 205 258 204 
+178 3 3 0 4 205 43 141 258 
+179 3 3 0 4 204 258 131 36 
+180 3 3 0 4 258 141 8 131 
+181 3 3 0 4 55 166 259 167 
+182 3 3 0 4 166 22 104 259 
+183 3 3 0 4 167 259 205 72 
+184 3 3 0 4 259 104 43 205 
+185 3 3 0 4 57 173 260 171 
+186 3 3 0 4 173 73 207 260 
+187 3 3 0 4 171 260 114 26 
+188 3 3 0 4 260 207 55 114 
+189 3 3 0 4 18 98 261 99 
+190 3 3 0 4 98 39 137 261 
+191 3 3 0 4 99 261 173 57 
+192 3 3 0 4 261 137 73 173 
+193 3 3 0 4 73 206 262 207 
+194 3 3 0 4 206 44 142 262 
+195 3 3 0 4 207 262 166 55 
+196 3 3 0 4 262 142 22 166 
+197 3 3 0 4 39 136 263 137 
+198 3 3 0 4 136 9 83 263 
+199 3 3 0 4 137 263 206 73 
+200 3 3 0 4 263 83 44 206 
+201 3 3 0 4 74 209 264 208 
+202 3 3 0 4 209 38 134 264 
+203 3 3 0 4 208 264 133 37 
+204 3 3 0 4 264 134 4 133 
+205 3 3 0 4 58 175 265 174 
+206 3 3 0 4 175 74 208 265 
+207 3 3 0 4 174 265 96 17 
+208 3 3 0 4 265 208 37 96 
+209 3 3 0 4 60 180 266 181 
+210 3 3 0 4 180 18 97 266 
+211 3 3 0 4 181 266 209 74 
+212 3 3 0 4 266 97 38 209 
+213 3 3 0 4 27 118 267 116 
+214 3 3 0 4 118 60 181 267 
+215 3 3 0 4 116 267 175 58 
+216 3 3 0 4 267 181 74 175 
+217 3 3 0 4 27 116 268 119 
+218 3 3 0 4 116 58 176 268 
+219 3 3 0 4 119 268 184 61 
+220 3 3 0 4 268 176 75 184 
+221 3 3 0 4 60 182 269 180 
+222 3 3 0 4 182 76 211 269 
+223 3 3 0 4 180 269 98 18 
+224 3 3 0 4 269 211 39 98 
+225 3 3 0 4 27 117 270 118 
+226 3 3 0 4 117 59 178 270 
+227 3 3 0 4 118 270 182 60 
+228 3 3 0 4 270 178 76 182 
+229 3 3 0 4 76 212 271 211 
+230 3 3 0 4 212 45 143 271 
+231 3 3 0 4 211 271 136 39 
+232 3 3 0 4 271 143 9 136 
+233 3 3 0 4 59 177 272 178 
+234 3 3 0 4 177 23 106 272 
+235 3 3 0 4 178 272 212 76 
+236 3 3 0 4 272 106 45 212 
+237 3 3 0 4 27 119 273 117 
+238 3 3 0 4 119 61 185 273 
+239 3 3 0 4 117 273 179 59 
+240 3 3 0 4 273 185 77 179 
+$ENDELM

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.