]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add TriangulationDescription::Utilities::create_description_from_triangulation()... 12312/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 29 May 2021 12:56:44 +0000 (14:56 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Sun, 30 May 2021 14:00:22 +0000 (16:00 +0200)
include/deal.II/grid/tria_description.h
source/grid/tria_description.cc
source/grid/tria_description.inst.in
tests/fullydistributed_grids/repartitioning_01.cc [new file with mode: 0644]
tests/fullydistributed_grids/repartitioning_01.mpirun=4.output [new file with mode: 0644]
tests/fullydistributed_grids/tests.h

index a2b38df3b6e7e7aeee88575c180a0a5318f855fd..2d1ad8f1f55a1a757aa7745078cdeeb2c9d3d374 100644 (file)
@@ -23,6 +23,8 @@
 #include <deal.II/grid/cell_id.h>
 #include <deal.II/grid/tria.h>
 
+#include <deal.II/lac/la_parallel_vector.h>
+
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -446,6 +448,24 @@ namespace TriangulationDescription
         TriangulationDescription::Settings::default_setting,
       const unsigned int my_rank_in = numbers::invalid_unsigned_int);
 
+    /**
+     * Similar to the above function but the owner of active cells are provided
+     * by a cell vector (see also
+     * parallel::TriangulationBase::global_active_cell_index_partitioner() and
+     * CellAccessor::global_active_cell_index()). This function allows to
+     * repartition distributed Triangulation objects.
+     *
+     * @note The communicator is extracted from the Triangulation @p tria.
+     *
+     * @note The multgrid levels are currently not constructed, since
+     *   @p partition only describes the partitioning of the active level.
+     */
+    template <int dim, int spacedim>
+    Description<dim, spacedim>
+    create_description_from_triangulation(
+      const Triangulation<dim, spacedim> &              tria,
+      const LinearAlgebra::distributed::Vector<double> &partition);
+
 
     /**
      * Construct a TriangulationDescription::Description. In contrast
index 5f6ee1e00dcdd8dec1cd60fa56a6823fdf0fae42..89d9dd6c94197f9cb533af8db5967fe657479129 100644 (file)
@@ -15,6 +15,7 @@
 
 #include <deal.II/base/geometry_info.h>
 #include <deal.II/base/mpi.h>
+#include <deal.II/base/mpi_consensus_algorithms.h>
 
 #include <deal.II/distributed/tria.h>
 
@@ -525,6 +526,386 @@ namespace TriangulationDescription
 #endif
     }
 
+    template <int dim, int spacedim>
+    struct DescriptionTemp
+    {
+      template <class Archive>
+      void
+      serialize(Archive &ar, const unsigned int /*version*/)
+      {
+        ar &coarse_cells;
+        ar &coarse_cell_vertices;
+        ar &coarse_cell_index_to_coarse_cell_id;
+        ar &cell_infos;
+      }
+
+      void
+      merge(const DescriptionTemp<dim, spacedim> &other)
+      {
+        this->cell_infos.resize(
+          std::max(other.cell_infos.size(), this->cell_infos.size()));
+
+        this->coarse_cells.insert(this->coarse_cells.end(),
+                                  other.coarse_cells.begin(),
+                                  other.coarse_cells.end());
+        this->coarse_cell_vertices.insert(this->coarse_cell_vertices.end(),
+                                          other.coarse_cell_vertices.begin(),
+                                          other.coarse_cell_vertices.end());
+        this->coarse_cell_index_to_coarse_cell_id.insert(
+          this->coarse_cell_index_to_coarse_cell_id.end(),
+          other.coarse_cell_index_to_coarse_cell_id.begin(),
+          other.coarse_cell_index_to_coarse_cell_id.end());
+
+        for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
+          this->cell_infos[i].insert(this->cell_infos[i].end(),
+                                     other.cell_infos[i].begin(),
+                                     other.cell_infos[i].end());
+      }
+
+      void
+      reduce()
+      {
+        // make coarse cells unique
+        {
+          std::vector<std::tuple<types::coarse_cell_id,
+                                 dealii::CellData<dim>,
+                                 unsigned int>>
+            temp;
+
+          for (unsigned int i = 0; i < this->coarse_cells.size(); ++i)
+            temp.emplace_back(this->coarse_cell_index_to_coarse_cell_id[i],
+                              this->coarse_cells[i],
+                              i);
+
+          std::sort(temp.begin(), temp.end(), [](const auto &a, const auto &b) {
+            return std::get<0>(a) < std::get<0>(b);
+          });
+          temp.erase(std::unique(temp.begin(),
+                                 temp.end(),
+                                 [](const auto &a, const auto &b) {
+                                   return std::get<0>(a) == std::get<0>(b);
+                                 }),
+                     temp.end());
+          std::sort(temp.begin(), temp.end(), [](const auto &a, const auto &b) {
+            return std::get<2>(a) < std::get<2>(b);
+          });
+
+          this->coarse_cell_index_to_coarse_cell_id.resize(temp.size());
+          this->coarse_cells.resize(temp.size());
+
+          for (unsigned int i = 0; i < temp.size(); ++i)
+            {
+              this->coarse_cell_index_to_coarse_cell_id[i] =
+                std::get<0>(temp[i]);
+              this->coarse_cells[i] = std::get<1>(temp[i]);
+            }
+        }
+
+        // make coarse cell vertices unique
+        {
+          std::sort(this->coarse_cell_vertices.begin(),
+                    this->coarse_cell_vertices.end(),
+                    [](const auto &a, const auto &b) {
+                      return a.first < b.first;
+                    });
+          this->coarse_cell_vertices.erase(
+            std::unique(this->coarse_cell_vertices.begin(),
+                        this->coarse_cell_vertices.end(),
+                        [](const auto &a, const auto &b) {
+                          return a.first == b.first;
+                        }),
+            this->coarse_cell_vertices.end());
+        }
+
+        // make cells unique
+        for (unsigned int i = 0; i < this->cell_infos.size(); ++i)
+          {
+            std::sort(this->cell_infos[i].begin(),
+                      this->cell_infos[i].end(),
+                      [](const auto &a, const auto &b) { return a.id < b.id; });
+            this->cell_infos[i].erase(std::unique(this->cell_infos[i].begin(),
+                                                  this->cell_infos[i].end(),
+                                                  [](const auto &a,
+                                                     const auto &b) {
+                                                    return a.id == b.id;
+                                                  }),
+                                      this->cell_infos[i].end());
+          }
+      }
+
+      Description<dim, spacedim>
+      convert(const MPI_Comm comm,
+              const typename Triangulation<dim, spacedim>::MeshSmoothing
+                mesh_smoothing)
+      {
+        Description<dim, spacedim> description;
+
+        // copy communicator
+        description.comm = comm;
+
+        // set settings (no multigrid levels are set for now)
+        description.settings =
+          TriangulationDescription::Settings::default_setting;
+
+        // use mesh smoothing from base triangulation
+        description.smoothing = mesh_smoothing;
+
+        std::map<unsigned int, unsigned int> map;
+
+        for (unsigned int i = 0; i < this->coarse_cell_vertices.size(); ++i)
+          {
+            description.coarse_cell_vertices.push_back(
+              this->coarse_cell_vertices[i].second);
+            map[this->coarse_cell_vertices[i].first] = i;
+          }
+
+        description.coarse_cells = this->coarse_cells;
+
+        for (auto &cell : description.coarse_cells)
+          for (unsigned int v = 0; v < cell.vertices.size(); ++v)
+            cell.vertices[v] = map[cell.vertices[v]];
+
+        description.coarse_cell_index_to_coarse_cell_id =
+          this->coarse_cell_index_to_coarse_cell_id;
+        description.cell_infos = this->cell_infos;
+
+        return description;
+      }
+
+      std::vector<dealii::CellData<dim>> coarse_cells;
+
+      std::vector<std::pair<unsigned int, Point<spacedim>>>
+        coarse_cell_vertices;
+
+      std::vector<types::coarse_cell_id> coarse_cell_index_to_coarse_cell_id;
+
+      std::vector<std::vector<CellData<dim>>> cell_infos;
+    };
+
+
+
+    template <int dim, int spacedim>
+    void
+    fill_cell_infos(const TriaIterator<CellAccessor<dim, spacedim>> &cell,
+                    std::vector<std::vector<CellData<dim>>> &        cell_infos,
+                    const LinearAlgebra::distributed::Vector<double> &partition)
+    {
+      if (cell->user_flag_set())
+        return; // cell has been already added -> nothing to do
+
+      cell->set_user_flag();
+
+      CellData<dim> cell_info;
+
+      // save coarse-cell id
+      cell_info.id = cell->id().template to_binary<dim>();
+
+      // save boundary_ids of each face of this cell
+      for (const auto f : cell->face_indices())
+        {
+          types::boundary_id boundary_ind = cell->face(f)->boundary_id();
+          if (boundary_ind != numbers::internal_face_boundary_id)
+            cell_info.boundary_ids.emplace_back(f, boundary_ind);
+        }
+
+      // save manifold id
+      {
+        // ... of cell
+        cell_info.manifold_id = cell->manifold_id();
+
+        // ... of lines
+        if (dim >= 2)
+          for (const auto line : cell->line_indices())
+            cell_info.manifold_line_ids[line] = cell->line(line)->manifold_id();
+
+        // ... of quads
+        if (dim == 3)
+          for (const auto f : cell->face_indices())
+            cell_info.manifold_quad_ids[f] = cell->quad(f)->manifold_id();
+      }
+
+      // subdomain and level subdomain id
+      if (cell->is_active())
+        cell_info.subdomain_id = static_cast<unsigned int>(
+          partition[cell->global_active_cell_index()]);
+      else
+        cell_info.subdomain_id = numbers::artificial_subdomain_id;
+
+      cell_info.level_subdomain_id = numbers::artificial_subdomain_id;
+
+      cell_infos[cell->level()].emplace_back(cell_info);
+
+      if (cell->level() != 0) // proceed with parent
+        fill_cell_infos(cell->parent(), cell_infos, partition);
+    }
+
+
+
+    template <int dim, int spacedim>
+    Description<dim, spacedim>
+    create_description_from_triangulation(
+      const Triangulation<dim, spacedim> &              tria,
+      const LinearAlgebra::distributed::Vector<double> &partition)
+    {
+      // 1) determine processes owning locally owned cells
+      const std::vector<unsigned int> relevant_processes = [&]() {
+        std::set<unsigned int> relevant_processes;
+
+        for (unsigned int i = 0; i < partition.local_size(); ++i)
+          relevant_processes.insert(
+            static_cast<unsigned int>(partition.local_element(i)));
+
+        return std::vector<unsigned int>(relevant_processes.begin(),
+                                         relevant_processes.end());
+      }();
+
+      // 2) determine coinciding vertices (important for periodicity)
+      std::map<unsigned int, std::vector<unsigned int>>
+                                           coinciding_vertex_groups;
+      std::map<unsigned int, unsigned int> vertex_to_coinciding_vertex_group;
+
+      GridTools::collect_coinciding_vertices(tria,
+                                             coinciding_vertex_groups,
+                                             vertex_to_coinciding_vertex_group);
+
+      const auto add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells =
+        [&coinciding_vertex_groups, &vertex_to_coinciding_vertex_group](
+          const auto &       cell,
+          std::vector<bool> &vertices_owned_by_locally_owned_cells) {
+          for (const auto v : cell->vertex_indices())
+            {
+              vertices_owned_by_locally_owned_cells[cell->vertex_index(v)] =
+                true;
+              const auto coinciding_vertex_group =
+                vertex_to_coinciding_vertex_group.find(cell->vertex_index(v));
+              if (coinciding_vertex_group !=
+                  vertex_to_coinciding_vertex_group.end())
+                for (const auto &co_vertex : coinciding_vertex_groups.at(
+                       coinciding_vertex_group->second))
+                  vertices_owned_by_locally_owned_cells[co_vertex] = true;
+            }
+        };
+
+      // 3) create a description (locally owned cell and a layer of ghost cells
+      // and all their parents)
+      std::vector<DescriptionTemp<dim, spacedim>> description_temp(
+        relevant_processes.size());
+
+      for (unsigned int i = 0; i < description_temp.size(); ++i)
+        {
+          const unsigned int proc               = relevant_processes[i];
+          auto &             description_temp_i = description_temp[i];
+          description_temp_i.cell_infos.resize(
+            tria.get_triangulation().n_global_levels());
+
+          // clear user_flags
+          std::vector<bool> old_user_flags;
+          tria.save_user_flags(old_user_flags);
+
+          const_cast<dealii::Triangulation<dim, spacedim> &>(tria)
+            .clear_user_flags();
+
+          // mark all vertices attached to locally owned cells
+          std::vector<bool> vertices_owned_by_locally_owned_cells(
+            tria.n_vertices());
+          for (const auto &cell : tria.active_cell_iterators())
+            if (cell->is_locally_owned() &&
+                static_cast<unsigned int>(
+                  partition[cell->global_active_cell_index()]) == proc)
+              add_vertices_of_cell_to_vertices_owned_by_locally_owned_cells(
+                cell, vertices_owned_by_locally_owned_cells);
+
+          // helper function if a cell is locally relevant (active and
+          // connected to cell via a vertex)
+          const auto is_locally_relevant_on_level = [&](const auto &cell) {
+            for (const auto v : cell->vertex_indices())
+              if (vertices_owned_by_locally_owned_cells[cell->vertex_index(v)])
+                return true;
+            return false;
+          };
+
+          // collect locally relevant cells (including their parents)
+          for (const auto &cell : tria.active_cell_iterators())
+            if (is_locally_relevant_on_level(cell))
+              fill_cell_infos(cell, description_temp_i.cell_infos, partition);
+
+          // collect coarse-grid cells
+          std::vector<bool> vertices_locally_relevant(tria.n_vertices(), false);
+
+          for (const auto &cell : tria.cell_iterators_on_level(0))
+            {
+              if (!cell->user_flag_set())
+                continue;
+
+              // extract cell definition
+              dealii::CellData<dim> cell_data(cell->n_vertices());
+              cell_data.material_id = cell->material_id();
+              cell_data.manifold_id = cell->manifold_id();
+              for (const auto v : cell->vertex_indices())
+                cell_data.vertices[v] = cell->vertex_index(v);
+              description_temp_i.coarse_cells.push_back(cell_data);
+
+              // mark cell vertices as relevant
+              for (const auto v : cell->vertex_indices())
+                vertices_locally_relevant[cell->vertex_index(v)] = true;
+
+              // save translation for corase grid: lid -> gid
+              description_temp_i.coarse_cell_index_to_coarse_cell_id.push_back(
+                cell->id().get_coarse_cell_id());
+            }
+
+          // collect coarse-grid vertices
+          for (unsigned int i = 0; i < vertices_locally_relevant.size(); ++i)
+            if (vertices_locally_relevant[i])
+              description_temp_i.coarse_cell_vertices.emplace_back(
+                i, tria.get_vertices()[i]);
+
+          // restore flags
+          const_cast<dealii::Triangulation<dim, spacedim> &>(tria)
+            .load_user_flags(old_user_flags);
+        }
+
+      // collect description from all processes that used to own locally-owned
+      // active cells of this process in a single description
+      DescriptionTemp<dim, spacedim> description_merged;
+
+      dealii::Utilities::MPI::ConsensusAlgorithms::AnonymousProcess<char, char>
+        process([&]() { return relevant_processes; },
+                [&](const unsigned int other_rank,
+                    std::vector<char> &send_buffer) {
+                  const auto ptr = std::find(relevant_processes.begin(),
+                                             relevant_processes.end(),
+                                             other_rank);
+
+                  Assert(ptr != relevant_processes.end(), ExcInternalError());
+
+                  const auto other_rank_index =
+                    std::distance(relevant_processes.begin(), ptr);
+
+                  send_buffer =
+                    dealii::Utilities::pack(description_temp[other_rank_index],
+                                            false);
+                },
+                [&](const unsigned int &,
+                    const std::vector<char> &recv_buffer,
+                    std::vector<char> &) {
+                  description_merged.merge(
+                    dealii::Utilities::unpack<DescriptionTemp<dim, spacedim>>(
+                      recv_buffer, false));
+                });
+
+      dealii::Utilities::MPI::ConsensusAlgorithms::Selector<char, char>(
+        process, tria.get_communicator())
+        .run();
+
+      // remove redundant entries
+      description_merged.reduce();
+
+      // convert to actual description
+      return description_merged.convert(tria.get_communicator(),
+                                        tria.get_mesh_smoothing());
+    }
+
   } // namespace Utilities
 } // namespace TriangulationDescription
 
index 9ca20f01570eb265897f6bc3315444d2c4f210a8..fca624f43ed1e2d1ba361b990b31ac61658f0f82 100644 (file)
@@ -45,6 +45,11 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS)
                                        deal_II_space_dimension>::MeshSmoothing
             smoothing,
           const TriangulationDescription::Settings);
+
+        template Description<deal_II_dimension, deal_II_space_dimension>
+        create_description_from_triangulation(
+          const Triangulation<deal_II_dimension, deal_II_space_dimension> &tria,
+          const LinearAlgebra::distributed::Vector<double> &partition);
 #endif
       \}
     \}
diff --git a/tests/fullydistributed_grids/repartitioning_01.cc b/tests/fullydistributed_grids/repartitioning_01.cc
new file mode 100644 (file)
index 0000000..f840ea9
--- /dev/null
@@ -0,0 +1,113 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 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
+// TriangulationDescription::Utilities::create_description_from_triangulation()
+// with repartitioning capabilities.
+
+#include <deal.II/base/mpi_consensus_algorithms.h>
+
+#include <deal.II/distributed/fully_distributed_tria.h>
+#include <deal.II/distributed/tria.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/grid/grid_tools.h>
+#include <deal.II/grid/tria_description.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/numerics/data_out.h>
+
+#include "./tests.h"
+
+using namespace dealii;
+
+
+template <int dim, int spacedim>
+LinearAlgebra::distributed::Vector<double>
+partition_distributed_triangulation(const Triangulation<dim, spacedim> &tria_in,
+                                    const unsigned int n_partitions)
+{
+  const auto tria =
+    dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria_in);
+
+  Assert(tria, ExcNotImplemented());
+
+  LinearAlgebra::distributed::Vector<double> partition(
+    tria->global_active_cell_index_partitioner().lock());
+
+  for (const auto &cell : tria_in.active_cell_iterators())
+    if (cell->is_locally_owned())
+      partition[cell->global_active_cell_index()] =
+        std::floor(cell->center()[0] * n_partitions);
+
+  partition.update_ghost_values();
+
+  return partition;
+}
+
+
+template <int dim>
+void
+test(const MPI_Comm comm, const unsigned int n_partitions)
+{
+  parallel::distributed::Triangulation<dim> tria(comm);
+  GridGenerator::subdivided_hyper_cube(tria, 4);
+  tria.refine_global(3);
+
+  const auto partition_new =
+    partition_distributed_triangulation(tria, n_partitions);
+
+  // repartition triangulation so that it has strided partitioning
+  const auto construction_data =
+    TriangulationDescription::Utilities::create_description_from_triangulation(
+      tria, partition_new);
+
+  parallel::fullydistributed::Triangulation<dim> tria_pft(comm);
+  tria_pft.create_triangulation(construction_data);
+
+  FE_Q<dim>       fe(2);
+  DoFHandler<dim> dof_handler(tria_pft);
+  dof_handler.distribute_dofs(fe);
+
+  // print statistics
+  print_statistics(tria_pft);
+  print_statistics(dof_handler);
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  MPI_Comm comm = MPI_COMM_WORLD;
+
+  deallog.push("all");
+  test<2>(comm, Utilities::MPI::n_mpi_processes(comm));
+  deallog.pop();
+
+  // test that we can eliminate processes
+  deallog.push("reduced");
+  test<2>(comm,
+          std::max<unsigned int>(1, Utilities::MPI::n_mpi_processes(comm) / 2));
+  deallog.pop();
+}
diff --git a/tests/fullydistributed_grids/repartitioning_01.mpirun=4.output b/tests/fullydistributed_grids/repartitioning_01.mpirun=4.output
new file mode 100644 (file)
index 0000000..9360fea
--- /dev/null
@@ -0,0 +1,63 @@
+
+DEAL:0:all::n_levels:                  4
+DEAL:0:all::n_cells:                   456
+DEAL:0:all::n_active_cells:            344
+DEAL:0:all::
+DEAL:0:all::n_dofs:                             4225
+DEAL:0:all::n_locally_owned_dofs:               1105
+DEAL:0:all::
+DEAL:0:reduced::n_levels:                  4
+DEAL:0:reduced::n_cells:                   796
+DEAL:0:reduced::n_active_cells:            600
+DEAL:0:reduced::
+DEAL:0:reduced::n_dofs:                             4225
+DEAL:0:reduced::n_locally_owned_dofs:               2145
+DEAL:0:reduced::
+
+DEAL:1:all::n_levels:                  4
+DEAL:1:all::n_cells:                   572
+DEAL:1:all::n_active_cells:            432
+DEAL:1:all::
+DEAL:1:all::n_dofs:                             4225
+DEAL:1:all::n_locally_owned_dofs:               1040
+DEAL:1:all::
+DEAL:1:reduced::n_levels:                  4
+DEAL:1:reduced::n_cells:                   796
+DEAL:1:reduced::n_active_cells:            600
+DEAL:1:reduced::
+DEAL:1:reduced::n_dofs:                             4225
+DEAL:1:reduced::n_locally_owned_dofs:               2080
+DEAL:1:reduced::
+
+
+DEAL:2:all::n_levels:                  4
+DEAL:2:all::n_cells:                   572
+DEAL:2:all::n_active_cells:            432
+DEAL:2:all::
+DEAL:2:all::n_dofs:                             4225
+DEAL:2:all::n_locally_owned_dofs:               1040
+DEAL:2:all::
+DEAL:2:reduced::n_levels:                  1
+DEAL:2:reduced::n_cells:                   1
+DEAL:2:reduced::n_active_cells:            1
+DEAL:2:reduced::
+DEAL:2:reduced::n_dofs:                             4225
+DEAL:2:reduced::n_locally_owned_dofs:               0
+DEAL:2:reduced::
+
+
+DEAL:3:all::n_levels:                  4
+DEAL:3:all::n_cells:                   456
+DEAL:3:all::n_active_cells:            344
+DEAL:3:all::
+DEAL:3:all::n_dofs:                             4225
+DEAL:3:all::n_locally_owned_dofs:               1040
+DEAL:3:all::
+DEAL:3:reduced::n_levels:                  1
+DEAL:3:reduced::n_cells:                   1
+DEAL:3:reduced::n_active_cells:            1
+DEAL:3:reduced::
+DEAL:3:reduced::n_dofs:                             4225
+DEAL:3:reduced::n_locally_owned_dofs:               0
+DEAL:3:reduced::
+
index d0ae2ca1319142df5c3b115f2f1b61b288171bfa..4029e7f7b5f223d4e424ee06f25799dcb60fe057 100644 (file)
@@ -21,9 +21,7 @@
 
 template <int dim, int spacedim>
 void
-print_statistics(
-  const parallel::fullydistributed::Triangulation<dim, spacedim> &tria,
-  bool                                                            do_mg = false)
+print_statistics(const Triangulation<dim, spacedim> &tria, bool do_mg = false)
 {
   deallog << "n_levels:                  " << tria.n_levels() << std::endl;
   deallog << "n_cells:                   " << tria.n_cells() << std::endl;

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.