]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add GridGenerator::alfeld_split_of_simplex_mesh
authorTimo Heister <timo.heister@gmail.com>
Fri, 30 Aug 2024 19:27:40 +0000 (15:27 -0400)
committerTimo Heister <timo.heister@gmail.com>
Fri, 30 Aug 2024 19:29:49 +0000 (15:29 -0400)
include/deal.II/grid/grid_generator.h
source/grid/grid_generator.cc
source/grid/grid_generator.inst.in
tests/simplex/alfeld_split_01.cc [new file with mode: 0644]
tests/simplex/alfeld_split_01.output [new file with mode: 0644]

index 28d2c743627802baae2080e7001eba7df5d96d94..d2a5eec06c49d3022a2cd0930f10fd82e9fa44cd 100644 (file)
@@ -2424,6 +2424,25 @@ namespace GridGenerator
                                     Triangulation<1, spacedim>       &out_tria);
 #endif
 
+
+  /**
+   * Perform an Alfeld split (also called barycentric refinement) of a simplex
+   * mesh.
+   *
+   * This function takes a simplex mesh (@p in_tria) and subdivides each
+   * triangle into three triangles with a single new vertex (the barycenter). In
+   * the process, the simplex mesh is flattened (no hierarchy is kept).
+   *
+   * @note Currently only implemented for @p dim = 2.
+   *
+   * Also see
+   * @ref simplex "Simplex support".
+   */
+  template <int dim, int spacedim>
+  void
+  alfeld_split_of_simplex_mesh(const Triangulation<dim, spacedim> &in_tria,
+                               Triangulation<dim, spacedim>       &out_tria);
+
   /**
    * Namespace Airfoil contains classes and functions in order to create a
    * C-type mesh for the (flow) field around Joukowski or NACA airfoils.
index 76229e99a54ffdd89515797207e82f559727752a..2c483da6e423f0fb88d6159bd0a3588c4d70ea44 100644 (file)
@@ -8403,6 +8403,209 @@ namespace GridGenerator
 
 
 
+  template <int dim, int spacedim>
+  void
+  alfeld_split_of_simplex_mesh(const Triangulation<dim, spacedim> &in_tria,
+                               Triangulation<dim, spacedim>       &out_tria)
+  {
+    Assert(dim == 2, ExcNotImplemented());
+
+    Triangulation<dim, spacedim> temp_tria;
+    if (in_tria.n_global_levels() > 1)
+      {
+        AssertThrow(!in_tria.has_hanging_nodes(), ExcNotImplemented());
+        GridGenerator::flatten_triangulation(in_tria, temp_tria);
+      }
+    const Triangulation<dim, spacedim> &ref_tria =
+      in_tria.n_global_levels() > 1 ? temp_tria : in_tria;
+
+    // Three triangles connecting to barycenter with vertex index 3:
+    static const ndarray<unsigned int, 3, 3> table_2D_cell = {
+      {{{0, 1, 3}}, {{1, 2, 3}}, {{2, 0, 3}}}};
+
+    /* Boundary-faces 2d:
+     * After converting, each of the 4 quadrilateral faces is defined by faces
+     * of 2 different triangles, i.e., lines. Note that lines are defined by 2
+     * vertices.
+     */
+    static const ndarray<unsigned int, 4, 2, 2>
+      vertex_ids_for_boundary_faces_2d = {
+        {{{{{0, 1}}}}, {{{{1, 2}}}}, {{{{2, 0}}}}}};
+
+
+    std::vector<Point<spacedim>> vertices;
+    std::vector<CellData<dim>>   cells;
+    SubCellData                  subcell_data;
+
+    // for each vertex we store the assigned index so that we only
+    // assign them a value once
+    std::vector<unsigned int> old_to_new_vertex_indices(
+      ref_tria.n_vertices(), numbers::invalid_unsigned_int);
+
+    // We first have to create all of the new vertices. To do this, we loop over
+    // all cells and on each cell
+    // (i) copy the existing vertex locations (and record their new indices in
+    // the 'old_to_new_vertex_indices' vector),
+    // (ii) create new barycenter vertex location
+    for (const auto &cell : ref_tria.cell_iterators())
+      {
+        // temporary array storing the global indices of each cell entity in the
+        // sequence: vertices, edges/faces, cell
+        std::array<unsigned int, dim == 2 ? 4 : 5> local_vertex_indices;
+
+        // (i) copy the existing vertex locations
+        Tensor<1, spacedim> barycenter;
+        for (const auto v : cell->vertex_indices())
+          {
+            const auto v_global = cell->vertex_index(v);
+
+            if (old_to_new_vertex_indices[v_global] ==
+                numbers::invalid_unsigned_int)
+              {
+                old_to_new_vertex_indices[v_global] = vertices.size();
+                vertices.push_back(cell->vertex(v));
+              }
+
+            AssertIndexRange(v, local_vertex_indices.size());
+            local_vertex_indices[v] = old_to_new_vertex_indices[v_global];
+
+            barycenter += vertices[local_vertex_indices[v]] - Point<spacedim>();
+          }
+
+        // (ii) barycenter:
+        local_vertex_indices[3] = vertices.size();
+        vertices.push_back(Point<spacedim>() + barycenter / 3.);
+
+        // helper function for creating cells and subcells
+        const auto add_cell = [&](const unsigned int struct_dim,
+                                  const auto        &index_vertices,
+                                  const unsigned int material_or_boundary_id,
+                                  const unsigned int manifold_id = 0) {
+          // sub-cell data only has to be stored if the information differs
+          // from the default
+          if (struct_dim < dim &&
+              (material_or_boundary_id == numbers::internal_face_boundary_id &&
+               manifold_id == numbers::flat_manifold_id))
+            return;
+
+          if (struct_dim == dim) // cells
+            {
+              if (dim == 2)
+                {
+                  AssertDimension(index_vertices.size(), 3);
+                }
+              else if (dim == 3)
+                {
+                  AssertDimension(index_vertices.size(), 4);
+                }
+
+              CellData<dim> cell_data(index_vertices.size());
+              for (unsigned int i = 0; i < index_vertices.size(); ++i)
+                {
+                  AssertIndexRange(index_vertices[i],
+                                   local_vertex_indices.size());
+                  cell_data.vertices[i] =
+                    local_vertex_indices[index_vertices[i]];
+                  cell_data.material_id =
+                    material_or_boundary_id; // inherit material id
+                  cell_data.manifold_id =
+                    manifold_id; // inherit cell-manifold id
+                }
+              cells.push_back(cell_data);
+            }
+          else if (dim == 2 && struct_dim == 1) // an edge of a simplex
+            {
+              Assert(index_vertices.size() == 2, ExcInternalError());
+              CellData<1> boundary_line(2);
+              boundary_line.boundary_id = material_or_boundary_id;
+              boundary_line.manifold_id = manifold_id;
+              for (unsigned int i = 0; i < index_vertices.size(); ++i)
+                {
+                  AssertIndexRange(index_vertices[i],
+                                   local_vertex_indices.size());
+                  boundary_line.vertices[i] =
+                    local_vertex_indices[index_vertices[i]];
+                }
+              subcell_data.boundary_lines.push_back(boundary_line);
+            }
+          else if (dim == 3 && struct_dim == 2) // a face of a tetrahedron
+            {
+              Assert(index_vertices.size() == 3, ExcInternalError());
+              CellData<2> boundary_quad(3);
+              boundary_quad.material_id = material_or_boundary_id;
+              boundary_quad.manifold_id = manifold_id;
+              for (unsigned int i = 0; i < index_vertices.size(); ++i)
+                {
+                  AssertIndexRange(index_vertices[i],
+                                   local_vertex_indices.size());
+                  boundary_quad.vertices[i] =
+                    local_vertex_indices[index_vertices[i]];
+                }
+              subcell_data.boundary_quads.push_back(boundary_quad);
+            }
+          else if (dim == 3 && struct_dim == 1) // an edge of a tetrahedron
+            {
+              Assert(index_vertices.size() == 2, ExcInternalError());
+              CellData<1> boundary_line(2);
+              boundary_line.boundary_id = material_or_boundary_id;
+              boundary_line.manifold_id = manifold_id;
+              for (unsigned int i = 0; i < index_vertices.size(); ++i)
+                {
+                  AssertIndexRange(index_vertices[i],
+                                   local_vertex_indices.size());
+                  boundary_line.vertices[i] =
+                    local_vertex_indices[index_vertices[i]];
+                }
+              subcell_data.boundary_lines.push_back(boundary_line);
+            }
+          else
+            {
+              DEAL_II_NOT_IMPLEMENTED();
+            }
+        };
+
+        const auto material_id_cell = cell->material_id();
+
+        // create cells one by one
+        if (dim == 2)
+          {
+            // get cell-manifold id from current quad cell
+            const auto manifold_id_cell = cell->manifold_id();
+
+            for (const auto &cell_vertices : table_2D_cell)
+              add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
+          }
+        else
+          DEAL_II_NOT_IMPLEMENTED();
+
+        // Set up sub-cell data.
+        for (const auto f : cell->face_indices())
+          {
+            const auto bid = cell->face(f)->boundary_id();
+            const auto mid = cell->face(f)->manifold_id();
+
+            // process boundary-faces: set boundary and manifold ids
+            if (dim == 2) // 2d boundary-faces
+              {
+                for (const auto &face_vertices :
+                     vertex_ids_for_boundary_faces_2d[f])
+                  add_cell(1, face_vertices, bid, mid);
+              }
+            else
+              DEAL_II_NOT_IMPLEMENTED();
+          }
+      }
+
+    out_tria.clear();
+    out_tria.create_triangulation(vertices, cells, subcell_data);
+
+    for (const auto i : out_tria.get_manifold_ids())
+      if (i != numbers::flat_manifold_id)
+        out_tria.set_manifold(i, FlatManifold<dim, spacedim>());
+  }
+
+
+
   template <template <int, int> class MeshType, int dim, int spacedim>
   DEAL_II_CXX20_REQUIRES(
     (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
index 0e928a8173d691626ac78e437d14ba0c1f028c3d..75a7d06b57795da0ba738bdf7ca3e6cf16efd0f2 100644 (file)
@@ -280,6 +280,11 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
       convert_hypercube_to_simplex_mesh(
         const Triangulation<deal_II_dimension, deal_II_space_dimension> &,
         Triangulation<deal_II_dimension, deal_II_space_dimension> &);
+
+      template void
+      alfeld_split_of_simplex_mesh(
+        const Triangulation<deal_II_dimension, deal_II_space_dimension> &,
+        Triangulation<deal_II_dimension, deal_II_space_dimension> &);
 #endif
     \}
   }
diff --git a/tests/simplex/alfeld_split_01.cc b/tests/simplex/alfeld_split_01.cc
new file mode 100644 (file)
index 0000000..7537e27
--- /dev/null
@@ -0,0 +1,95 @@
+/* ------------------------------------------------------------------------
+ *
+ * SPDX-License-Identifier: LGPL-2.1-or-later
+ * Copyright (C) 2020 - 2024 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * Part of the source code is dual licensed under Apache-2.0 WITH
+ * LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+ * governing the source code and code contributions can be found in
+ * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+ *
+ * ------------------------------------------------------------------------
+ *
+ * Test function GridGenerator::alfeld_split_of_simplex_mesh() in 2D.
+ */
+
+#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.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+template <int dim, int spacedim>
+void
+create_triangulation(Triangulation<dim, spacedim> &triangulation)
+{
+  GridGenerator::subdivided_hyper_cube(triangulation, 4);
+}
+
+template <int dim>
+void
+create_triangulation(Triangulation<dim, dim> &triangulation)
+{
+  GridGenerator::quarter_hyper_ball(triangulation);
+}
+
+template <int dim, int spacedim>
+void
+check()
+{
+  Triangulation<dim, spacedim> in_tria, out_tria;
+  GridGenerator::subdivided_hyper_cube_with_simplices(
+    in_tria, 2, 0.0, 1.0, true);
+
+  // make each cell a different material id
+  unsigned int m_id = 0;
+  for (const auto &cell : in_tria)
+    {
+      cell.set_material_id(m_id++);
+    }
+
+  GridGenerator::alfeld_split_of_simplex_mesh(in_tria, out_tria);
+
+  // write 2 outputs (total mesh and only surface mesh)
+  const auto grid_out = [](const auto &tria, const bool surface_mesh_only) {
+    GridOutFlags::Vtk flags;
+
+    if (surface_mesh_only)
+      {
+        flags.output_cells         = false;
+        flags.output_faces         = true;
+        flags.output_edges         = false;
+        flags.output_only_relevant = false;
+      }
+
+    GridOut grid_out;
+    grid_out.set_flags(flags);
+
+    grid_out.write_vtk(tria, deallog.get_file_stream());
+  };
+
+  grid_out(out_tria, false); // total mesh
+  grid_out(out_tria, true);  // only surface mesh
+
+  deallog << "OK!" << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+
+  deallog.push("2d");
+  check<2, 2>();
+  deallog.pop();
+}
diff --git a/tests/simplex/alfeld_split_01.output b/tests/simplex/alfeld_split_01.output
new file mode 100644 (file)
index 0000000..b533eed
--- /dev/null
@@ -0,0 +1,120 @@
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 17 double
+0.00000 0.00000 0
+0.500000 0.00000 0
+0.00000 0.500000 0
+0.166667 0.166667 0
+0.500000 0.500000 0
+0.333333 0.333333 0
+1.00000 0.00000 0
+0.666667 0.166667 0
+1.00000 0.500000 0
+0.833333 0.333333 0
+0.00000 1.00000 0
+0.166667 0.666667 0
+0.500000 1.00000 0
+0.333333 0.833333 0
+0.666667 0.666667 0
+1.00000 1.00000 0
+0.833333 0.833333 0
+
+CELLS 30 114
+3 0 1 3
+3 1 2 3
+3 2 0 3
+3 4 2 5
+3 2 1 5
+3 1 4 5
+3 1 6 7
+3 6 4 7
+3 4 1 7
+3 8 4 9
+3 4 6 9
+3 6 8 9
+3 2 4 11
+3 4 10 11
+3 10 2 11
+3 12 10 13
+3 10 4 13
+3 4 12 13
+3 4 8 14
+3 8 12 14
+3 12 4 14
+3 15 12 16
+3 12 8 16
+3 8 15 16
+2 0 1
+2 1 6
+2 6 8
+2 8 15
+2 12 10
+2 15 12
+
+CELL_TYPES 30
+5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 
+3 3 3 3 3 3 
+
+
+CELL_DATA 30
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 
+2 2 1 1 3 3 
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 
+-1 -1 -1 -1 -1 -1 
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 17 double
+0.00000 0.00000 0
+0.500000 0.00000 0
+0.00000 0.500000 0
+0.166667 0.166667 0
+0.500000 0.500000 0
+0.333333 0.333333 0
+1.00000 0.00000 0
+0.666667 0.166667 0
+1.00000 0.500000 0
+0.833333 0.333333 0
+0.00000 1.00000 0
+0.166667 0.666667 0
+0.500000 1.00000 0
+0.333333 0.833333 0
+0.666667 0.666667 0
+1.00000 1.00000 0
+0.833333 0.833333 0
+
+CELLS 8 24
+2 0 1
+2 1 6
+2 2 0
+2 6 8
+2 8 15
+2 10 2
+2 12 10
+2 15 12
+
+CELL_TYPES 8
+3 3 3 3 3 3 3 3 
+
+
+CELL_DATA 8
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+2 2 0 1 1 0 3 3 
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 -1 -1 -1 -1 
+DEAL:2d::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.