]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Simplex mesh (2D): implement refinement 11178/head
authorPeter Munch <peterrmuench@gmail.com>
Thu, 12 Nov 2020 17:24:28 +0000 (18:24 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Fri, 4 Dec 2020 18:51:07 +0000 (19:51 +0100)
include/deal.II/grid/tria_accessor.h
include/deal.II/grid/tria_accessor.templates.h
source/grid/tria.cc
source/grid/tria_accessor.cc
tests/simplex/refinement_01.cc [new file with mode: 0644]
tests/simplex/refinement_01.with_simplex_support=on.output [new file with mode: 0644]

index b2166d14241fd511fa3e709076dd34ee41ff9307..0af2ce4e39ccef02a7c3f78c075408985397742f 100644 (file)
@@ -61,6 +61,7 @@ namespace internal
   {
     class TriaObjects;
     struct Implementation;
+    struct ImplementationMixedMesh;
   } // namespace TriangulationImplementation
 
   namespace TriaAccessorImplementation
@@ -1810,6 +1811,8 @@ private:
   friend class Triangulation;
 
   friend struct dealii::internal::TriangulationImplementation::Implementation;
+  friend struct dealii::internal::TriangulationImplementation::
+    ImplementationMixedMesh;
   friend struct dealii::internal::TriaAccessorImplementation::Implementation;
 };
 
@@ -3807,6 +3810,8 @@ private:
   friend class parallel::TriangulationBase;
 
   friend struct dealii::internal::TriangulationImplementation::Implementation;
+  friend struct dealii::internal::TriangulationImplementation::
+    ImplementationMixedMesh;
 };
 
 
index 07624ba4a4749bbf2e4a3b8647816b7f1b7c0bfc..79779bf02349f4c6f4c490c2682d63934adf5246 100644 (file)
@@ -584,7 +584,7 @@ namespace internal
                        const unsigned int)
       {
         /*
-         * Default implementation used in 1d and 2d
+         * Default implementation used in 1d
          *
          * In 1d, face_orientation is always true
          */
index 68778d1a82d196d634b755cf26417fabf72cf599..d9d2b7737ae7eeab3ec49fcc45df63fb639ca370 100644 (file)
@@ -423,7 +423,7 @@ namespace
         std::vector<unsigned int> line_cell_count(triangulation.n_raw_lines(),
                                                   0);
         for (const auto &cell : triangulation.cell_iterators())
-          for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+          for (unsigned int l = 0; l < cell->n_lines(); ++l)
             ++line_cell_count[cell->line_index(l)];
         return line_cell_count;
       }
@@ -659,9 +659,7 @@ namespace
   template <int dim>
   bool
   has_distorted_children(
-    const typename Triangulation<dim, dim>::cell_iterator &cell,
-    std::integral_constant<int, dim>,
-    std::integral_constant<int, dim>)
+    const typename Triangulation<dim, dim>::cell_iterator &cell)
   {
     Assert(cell->has_children(), ExcInternalError());
 
@@ -694,9 +692,7 @@ namespace
   template <int dim, int spacedim>
   bool
   has_distorted_children(
-    const typename Triangulation<dim, spacedim>::cell_iterator &,
-    std::integral_constant<int, dim>,
-    std::integral_constant<int, spacedim>)
+    const typename Triangulation<dim, spacedim>::cell_iterator &)
   {
     return false;
   }
@@ -791,7 +787,7 @@ namespace
       adjacent_cells(2 * triangulation.n_raw_faces(), dummy);
 
     for (const auto &cell : triangulation.cell_iterators())
-      for (auto f : GeometryInfo<dim>::face_indices())
+      for (auto f : cell->face_indices())
         {
           const typename Triangulation<dim, spacedim>::face_iterator face =
             cell->face(f);
@@ -871,7 +867,7 @@ namespace
     // left_right_offset in this case as we want
     // the offset of the neighbor, not our own.
     for (const auto &cell : triangulation.cell_iterators())
-      for (auto f : GeometryInfo<dim>::face_indices())
+      for (auto f : cell->face_indices())
         {
           const unsigned int offset =
             (cell->direction_flag() ?
@@ -1478,6 +1474,16 @@ namespace internal
                   tria_level.face_orientations.size(),
                 true);
             }
+          else if (tria_level.dim == 2) // TODO: not needed for QUADs!!!
+            {                           // we need to pass the information here
+              tria_level.face_orientations.reserve(
+                total_cells * GeometryInfo<2>::faces_per_cell);
+              tria_level.face_orientations.insert(
+                tria_level.face_orientations.end(),
+                total_cells * GeometryInfo<2>::faces_per_cell -
+                  tria_level.face_orientations.size(),
+                true);
+            }
 
           if (tria_level.dim == 2 || tria_level.dim == 3)
             {
@@ -4047,6 +4053,540 @@ namespace internal
 
 
 
+      template <int dim, int spacedim>
+      static typename Triangulation<dim, spacedim>::DistortedCellList
+      execute_refinement_isotropic(Triangulation<dim, spacedim> &triangulation,
+                                   const bool check_for_distorted_cells)
+      {
+        AssertDimension(dim, 2);
+
+        {
+          typename Triangulation<dim, spacedim>::raw_cell_iterator
+            cell = triangulation.begin_active(triangulation.levels.size() - 1),
+            endc = triangulation.end();
+          for (; cell != endc; ++cell)
+            if (cell->used())
+              if (cell->refine_flag_set())
+                {
+                  triangulation.levels.push_back(
+                    std::make_unique<
+                      internal::TriangulationImplementation::TriaLevel>(dim));
+                  break;
+                }
+        }
+
+        for (typename Triangulation<dim, spacedim>::line_iterator line =
+               triangulation.begin_line();
+             line != triangulation.end_line();
+             ++line)
+          {
+            line->clear_user_flag();
+            line->clear_user_data();
+          }
+
+        unsigned int n_single_lines   = 0;
+        unsigned int n_lines_in_pairs = 0;
+        unsigned int needed_vertices  = 0;
+
+        for (int level = triangulation.levels.size() - 2; level >= 0; --level)
+          {
+            // count number of flagged cells on this level and compute
+            // how many new vertices and new lines will be needed
+            unsigned int needed_cells = 0;
+
+            for (const auto &cell :
+                 triangulation.active_cell_iterators_on_level(level))
+              if (cell->refine_flag_set())
+                {
+                  if (cell->reference_cell_type() == ReferenceCell::Type::Tri)
+                    {
+                      needed_cells += 4;
+                      needed_vertices += 0;
+                      n_single_lines += 3;
+                    }
+                  else if (cell->reference_cell_type() ==
+                           ReferenceCell::Type::Quad)
+                    {
+                      needed_cells += 4;
+                      needed_vertices += 1;
+                      n_single_lines += 4;
+                    }
+                  else
+                    {
+                      AssertThrow(false, ExcNotImplemented());
+                    }
+
+                  for (const auto line_no : cell->face_indices())
+                    {
+                      auto line = cell->line(line_no);
+                      if (line->has_children() == false)
+                        line->set_user_flag();
+                    }
+                }
+
+
+            const unsigned int used_cells =
+              std::count(triangulation.levels[level + 1]->cells.used.begin(),
+                         triangulation.levels[level + 1]->cells.used.end(),
+                         true);
+
+
+            reserve_space(*triangulation.levels[level + 1],
+                          used_cells + needed_cells,
+                          2,
+                          spacedim);
+
+            reserve_space(triangulation.levels[level + 1]->cells,
+                          needed_cells,
+                          0);
+          }
+
+        for (auto line = triangulation.begin_line();
+             line != triangulation.end_line();
+             ++line)
+          if (line->user_flag_set())
+            {
+              Assert(line->has_children() == false, ExcInternalError());
+              n_lines_in_pairs += 2;
+              needed_vertices += 1;
+            }
+
+        reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
+
+        needed_vertices += std::count(triangulation.vertices_used.begin(),
+                                      triangulation.vertices_used.end(),
+                                      true);
+
+        if (needed_vertices > triangulation.vertices.size())
+          {
+            triangulation.vertices.resize(needed_vertices, Point<spacedim>());
+            triangulation.vertices_used.resize(needed_vertices, false);
+          }
+
+        unsigned int next_unused_vertex = 0;
+
+        {
+          typename Triangulation<dim, spacedim>::active_line_iterator
+            line = triangulation.begin_active_line(),
+            endl = triangulation.end_line();
+          typename Triangulation<dim, spacedim>::raw_line_iterator
+            next_unused_line = triangulation.begin_raw_line();
+
+          for (; line != endl; ++line)
+            if (line->user_flag_set())
+              {
+                // this line needs to be refined
+
+                // find the next unused vertex and set it
+                // appropriately
+                while (triangulation.vertices_used[next_unused_vertex] == true)
+                  ++next_unused_vertex;
+                Assert(
+                  next_unused_vertex < triangulation.vertices.size(),
+                  ExcMessage(
+                    "Internal error: During refinement, the triangulation wants to access an element of the 'vertices' array but it turns out that the array is not large enough."));
+                triangulation.vertices_used[next_unused_vertex] = true;
+
+                triangulation.vertices[next_unused_vertex] = line->center(true);
+
+                bool pair_found = false;
+                (void)pair_found;
+                for (; next_unused_line != endl; ++next_unused_line)
+                  if (!next_unused_line->used() &&
+                      !(++next_unused_line)->used())
+                    {
+                      --next_unused_line;
+                      pair_found = true;
+                      break;
+                    }
+                Assert(pair_found, ExcInternalError());
+
+                line->set_children(0, next_unused_line->index());
+
+                const typename Triangulation<dim, spacedim>::raw_line_iterator
+                  children[2] = {next_unused_line, ++next_unused_line};
+
+                Assert(
+                  children[0]->used() == false,
+                  ExcMessage(
+                    "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
+                Assert(
+                  children[1]->used() == false,
+                  ExcMessage(
+                    "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
+
+                children[0]->set_bounding_object_indices(
+                  {line->vertex_index(0), next_unused_vertex});
+                children[1]->set_bounding_object_indices(
+                  {next_unused_vertex, line->vertex_index(1)});
+
+                children[0]->set_used_flag();
+                children[1]->set_used_flag();
+                children[0]->clear_children();
+                children[1]->clear_children();
+                children[0]->clear_user_data();
+                children[1]->clear_user_data();
+                children[0]->clear_user_flag();
+                children[1]->clear_user_flag();
+
+
+                children[0]->set_boundary_id_internal(line->boundary_id());
+                children[1]->set_boundary_id_internal(line->boundary_id());
+
+                children[0]->set_manifold_id(line->manifold_id());
+                children[1]->set_manifold_id(line->manifold_id());
+
+                line->clear_user_flag();
+              }
+        }
+
+        reserve_space(triangulation.faces->lines, 0, n_single_lines);
+
+        typename Triangulation<dim, spacedim>::DistortedCellList
+          cells_with_distorted_children;
+
+        typename Triangulation<dim, spacedim>::raw_line_iterator
+          next_unused_line = triangulation.begin_raw_line();
+
+        const auto create_children = [](auto &        triangulation,
+                                        unsigned int &next_unused_vertex,
+                                        auto &        next_unused_line,
+                                        auto &        next_unused_cell,
+                                        const auto &  cell) {
+          const auto ref_case = cell->refine_flag_set();
+          cell->clear_refine_flag();
+
+          unsigned int n_new_vertices = 0;
+
+          if (cell->reference_cell_type() == ReferenceCell::Type::Tri)
+            n_new_vertices = 6;
+          else if (cell->reference_cell_type() == ReferenceCell::Type::Quad)
+            n_new_vertices = 9;
+          else
+            AssertThrow(false, ExcNotImplemented());
+
+          std::vector<int> new_vertices(n_new_vertices);
+          for (unsigned int vertex_no = 0; vertex_no < cell->n_vertices();
+               ++vertex_no)
+            new_vertices[vertex_no] = cell->vertex_index(vertex_no);
+          for (unsigned int line_no = 0; line_no < cell->n_lines(); ++line_no)
+            if (cell->line(line_no)->has_children())
+              new_vertices[cell->n_vertices() + line_no] =
+                cell->line(line_no)->child(0)->vertex_index(1);
+
+          if (cell->reference_cell_type() == ReferenceCell::Type::Quad)
+            {
+              while (triangulation.vertices_used[next_unused_vertex] == true)
+                ++next_unused_vertex;
+              Assert(
+                next_unused_vertex < triangulation.vertices.size(),
+                ExcMessage(
+                  "Internal error: During refinement, the triangulation wants to access an element of the 'vertices' array but it turns out that the array is not large enough."));
+              triangulation.vertices_used[next_unused_vertex] = true;
+
+              new_vertices[8] = next_unused_vertex;
+
+              if (dim == spacedim)
+                {
+                  triangulation.vertices[next_unused_vertex] =
+                    cell->center(true);
+
+                  if (cell->user_flag_set())
+                    {
+                      cell->clear_user_flag();
+                      triangulation.vertices[next_unused_vertex] =
+                        cell->center(true, true);
+                    }
+                }
+              else
+                {
+                  cell->clear_user_flag();
+
+                  triangulation.vertices[next_unused_vertex] =
+                    cell->center(true, true);
+                }
+            }
+
+          std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
+                     12>
+                       new_lines;
+          unsigned int lmin = 0;
+          unsigned int lmax = 0;
+
+          if (cell->reference_cell_type() == ReferenceCell::Type::Tri)
+            {
+              lmin = 6;
+              lmax = 9;
+            }
+          else if (cell->reference_cell_type() == ReferenceCell::Type::Quad)
+            {
+              lmin = 8;
+              lmax = 12;
+            }
+          else
+            {
+              AssertThrow(false, ExcNotImplemented());
+            }
+
+          for (unsigned int l = lmin; l < lmax; ++l)
+            {
+              while (next_unused_line->used() == true)
+                ++next_unused_line;
+              new_lines[l] = next_unused_line;
+              ++next_unused_line;
+
+              Assert(
+                new_lines[l]->used() == false,
+                ExcMessage(
+                  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
+            }
+
+          if (true)
+            {
+              if (cell->reference_cell_type() == ReferenceCell::Type::Tri)
+                {
+                  // add lines in the right order [TODO: clean up]
+                  const auto ref = [&](const unsigned int face_no,
+                                       const unsigned int vertex_no) {
+                    if (cell->line(face_no)->child(0)->vertex_index(0) ==
+                          static_cast<unsigned int>(new_vertices[vertex_no]) ||
+                        cell->line(face_no)->child(0)->vertex_index(1) ==
+                          static_cast<unsigned int>(new_vertices[vertex_no]))
+                      {
+                        new_lines[2 * face_no + 0] =
+                          cell->line(face_no)->child(0);
+                        new_lines[2 * face_no + 1] =
+                          cell->line(face_no)->child(1);
+                      }
+                    else
+                      {
+                        new_lines[2 * face_no + 0] =
+                          cell->line(face_no)->child(1);
+                        new_lines[2 * face_no + 1] =
+                          cell->line(face_no)->child(0);
+                      }
+                  };
+
+                  ref(0, 0);
+                  ref(1, 1);
+                  ref(2, 2);
+
+                  new_lines[6]->set_bounding_object_indices(
+                    {new_vertices[3], new_vertices[4]});
+                  new_lines[7]->set_bounding_object_indices(
+                    {new_vertices[4], new_vertices[5]});
+                  new_lines[8]->set_bounding_object_indices(
+                    {new_vertices[5], new_vertices[3]});
+                }
+              else if (cell->reference_cell_type() == ReferenceCell::Type::Quad)
+                {
+                  unsigned int l = 0;
+                  for (const unsigned int face_no : cell->face_indices())
+                    for (unsigned int c = 0; c < 2; ++c, ++l)
+                      new_lines[l] = cell->line(face_no)->child(c);
+
+                  new_lines[8]->set_bounding_object_indices(
+                    {new_vertices[6], new_vertices[8]});
+                  new_lines[9]->set_bounding_object_indices(
+                    {new_vertices[8], new_vertices[7]});
+                  new_lines[10]->set_bounding_object_indices(
+                    {new_vertices[4], new_vertices[8]});
+                  new_lines[11]->set_bounding_object_indices(
+                    {new_vertices[8], new_vertices[5]});
+                }
+              else
+                {
+                  AssertThrow(false, ExcNotImplemented());
+                }
+            }
+
+
+          for (unsigned int l = lmin; l < lmax; ++l)
+            {
+              new_lines[l]->set_used_flag();
+              new_lines[l]->clear_user_flag();
+              new_lines[l]->clear_user_data();
+              new_lines[l]->clear_children();
+              // interior line
+              new_lines[l]->set_boundary_id_internal(
+                numbers::internal_face_boundary_id);
+              new_lines[l]->set_manifold_id(cell->manifold_id());
+            }
+
+          typename Triangulation<dim, spacedim>::raw_cell_iterator
+            subcells[GeometryInfo<dim>::max_children_per_cell];
+          while (next_unused_cell->used() == true)
+            ++next_unused_cell;
+
+          unsigned int n_children = 0;
+
+          if (cell->reference_cell_type() == ReferenceCell::Type::Tri)
+            n_children = 4;
+          else if (cell->reference_cell_type() == ReferenceCell::Type::Quad)
+            n_children = 4;
+          else
+            AssertThrow(false, ExcNotImplemented());
+
+          for (unsigned int i = 0; i < n_children; ++i)
+            {
+              Assert(
+                next_unused_cell->used() == false,
+                ExcMessage(
+                  "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be."));
+              subcells[i] = next_unused_cell;
+              ++next_unused_cell;
+              if (i % 2 == 1 && i < n_children - 1)
+                while (next_unused_cell->used() == true)
+                  ++next_unused_cell;
+            }
+
+          if (cell->reference_cell_type() == ReferenceCell::Type::Tri)
+            {
+              subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
+                                                        new_lines[8]->index(),
+                                                        new_lines[5]->index()});
+              subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
+                                                        new_lines[2]->index(),
+                                                        new_lines[6]->index()});
+              subcells[2]->set_bounding_object_indices({new_lines[7]->index(),
+                                                        new_lines[3]->index(),
+                                                        new_lines[4]->index()});
+              subcells[3]->set_bounding_object_indices({new_lines[6]->index(),
+                                                        new_lines[7]->index(),
+                                                        new_lines[8]->index()});
+
+              // subcell 0
+
+              const auto ref = [&](const unsigned int line_no,
+                                   const unsigned int vertex_no,
+                                   const unsigned int subcell_no,
+                                   const unsigned int subcell_line_no) {
+                if (new_lines[line_no]->vertex_index(1) !=
+                    static_cast<unsigned int>(new_vertices[vertex_no]))
+                  triangulation.levels[subcells[subcell_no]->level()]
+                    ->face_orientations[subcells[subcell_no]->index() *
+                                          GeometryInfo<2>::faces_per_cell +
+                                        subcell_line_no] = 0;
+              };
+
+              ref(0, 3, 0, 0);
+              ref(8, 5, 0, 1);
+              ref(5, 0, 0, 2);
+
+              ref(1, 1, 1, 0);
+              ref(2, 4, 1, 1);
+              ref(6, 3, 1, 2);
+
+              ref(7, 4, 2, 0);
+              ref(3, 2, 2, 1);
+              ref(4, 5, 2, 2);
+
+              ref(6, 4, 3, 0);
+              ref(7, 5, 3, 1);
+              ref(8, 3, 3, 2);
+
+              // triangulation.levels[subcells[1]->level()]->face_orientations[subcells[1]->index()
+              // * GeometryInfo<2>::faces_per_cell + 2] = 0;
+              // triangulation.levels[subcells[2]->level()]->face_orientations[subcells[2]->index()
+              // * GeometryInfo<2>::faces_per_cell + 0] = 0;
+            }
+          else if (cell->reference_cell_type() == ReferenceCell::Type::Quad)
+            {
+              subcells[0]->set_bounding_object_indices(
+                {new_lines[0]->index(),
+                 new_lines[8]->index(),
+                 new_lines[4]->index(),
+                 new_lines[10]->index()});
+              subcells[1]->set_bounding_object_indices(
+                {new_lines[8]->index(),
+                 new_lines[2]->index(),
+                 new_lines[5]->index(),
+                 new_lines[11]->index()});
+              subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
+                                                        new_lines[9]->index(),
+                                                        new_lines[10]->index(),
+                                                        new_lines[6]->index()});
+              subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
+                                                        new_lines[3]->index(),
+                                                        new_lines[11]->index(),
+                                                        new_lines[7]->index()});
+            }
+          else
+            {
+              AssertThrow(false, ExcNotImplemented());
+            }
+
+          types::subdomain_id subdomainid = cell->subdomain_id();
+
+          for (unsigned int i = 0; i < n_children; ++i)
+            {
+              subcells[i]->set_used_flag();
+              subcells[i]->clear_refine_flag();
+              subcells[i]->clear_user_flag();
+              subcells[i]->clear_user_data();
+              subcells[i]->clear_children();
+              // inherit material
+              // properties
+              subcells[i]->set_material_id(cell->material_id());
+              subcells[i]->set_manifold_id(cell->manifold_id());
+              subcells[i]->set_subdomain_id(subdomainid);
+
+              // TODO: here we assume that all children have the same reference
+              // cell type as the parent! This is justified for 2D.
+              triangulation.levels[subcells[i]->level()]
+                ->reference_cell_type[subcells[i]->index()] =
+                cell->reference_cell_type();
+
+              if (i % 2 == 0)
+                subcells[i]->set_parent(cell->index());
+            }
+
+          for (unsigned int i = 0; i < n_children / 2; ++i)
+            cell->set_children(2 * i, subcells[2 * i]->index());
+
+          cell->set_refinement_case(ref_case);
+
+          if (dim < spacedim)
+            for (unsigned int c = 0; c < n_children; ++c)
+              cell->child(c)->set_direction_flag(cell->direction_flag());
+        };
+
+        for (int level = 0;
+             level < static_cast<int>(triangulation.levels.size()) - 1;
+             ++level)
+          {
+            typename Triangulation<dim, spacedim>::raw_cell_iterator
+              next_unused_cell = triangulation.begin_raw(level + 1);
+
+            for (const auto &cell :
+                 triangulation.active_cell_iterators_on_level(level))
+              if (cell->refine_flag_set())
+                {
+                  if (cell->at_boundary())
+                    cell->set_user_flag();
+
+                  create_children(triangulation,
+                                  next_unused_vertex,
+                                  next_unused_line,
+                                  next_unused_cell,
+                                  cell);
+
+                  if (cell->reference_cell_type() ==
+                        ReferenceCell::Type::Quad &&
+                      check_for_distorted_cells &&
+                      has_distorted_children<dim, spacedim>(cell))
+                    cells_with_distorted_children.distorted_cells.push_back(
+                      cell);
+
+                  triangulation.signals.post_refinement_on_cell(cell);
+                }
+          }
+
+        return cells_with_distorted_children;
+      }
+
+
+
       /**
        * A function that performs the
        * refinement of a triangulation in 1d.
@@ -4298,6 +4838,26 @@ namespace internal
       {
         const unsigned int dim = 2;
 
+
+        {
+          bool flag_isotropic_mesh = true;
+          typename Triangulation<dim, spacedim>::raw_cell_iterator
+            cell = triangulation.begin(),
+            endc = triangulation.end();
+          for (; cell != endc; ++cell)
+            if (cell->used())
+              if (cell->refine_flag_set() == RefinementCase<dim>::cut_x ||
+                  cell->refine_flag_set() == RefinementCase<dim>::cut_y)
+                {
+                  flag_isotropic_mesh = false;
+                  break;
+                }
+
+          if (flag_isotropic_mesh)
+            return execute_refinement_isotropic(triangulation,
+                                                check_for_distorted_cells);
+        }
+
         // check whether a new level is needed we have to check for
         // this on the highest level only (on this, all used cells are
         // also active, so we only have to check for this)
@@ -4591,10 +5151,7 @@ namespace internal
                                   cell);
 
                   if (check_for_distorted_cells &&
-                      has_distorted_children(
-                        cell,
-                        std::integral_constant<int, dim>(),
-                        std::integral_constant<int, spacedim>()))
+                      has_distorted_children<dim, spacedim>(cell))
                     cells_with_distorted_children.distorted_cells.push_back(
                       cell);
                   // inform all listeners that cell refinement is done
@@ -8958,10 +9515,7 @@ namespace internal
                   // now see if we have created cells that are
                   // distorted and if so add them to our list
                   if (check_for_distorted_cells &&
-                      has_distorted_children(
-                        hex,
-                        std::integral_constant<int, dim>(),
-                        std::integral_constant<int, spacedim>()))
+                      has_distorted_children<dim, spacedim>(hex))
                     cells_with_distorted_children.distorted_cells.push_back(
                       hex);
 
@@ -9379,11 +9933,8 @@ namespace internal
       execute_refinement(Triangulation<dim, spacedim> &triangulation,
                          const bool check_for_distorted_cells)
       {
-        AssertThrow(false, ExcNotImplemented());
-        (void)triangulation;
-        (void)check_for_distorted_cells;
-
-        return {};
+        return Implementation::execute_refinement_isotropic(
+          triangulation, check_for_distorted_cells);
       }
 
       template <int dim, int spacedim>
@@ -9391,7 +9942,7 @@ namespace internal
       prevent_distorted_boundary_cells(
         Triangulation<dim, spacedim> &triangulation)
       {
-        AssertThrow(false, ExcNotImplemented());
+        // nothing to do since anisotropy is not supported
         (void)triangulation;
       }
 
@@ -9400,8 +9951,7 @@ namespace internal
       prepare_refinement_dim_dependent(
         Triangulation<dim, spacedim> &triangulation)
       {
-        AssertThrow(false, ExcNotImplemented());
-        (void)triangulation;
+        Implementation::prepare_refinement_dim_dependent(triangulation);
       }
 
       template <int dim, int spacedim>
@@ -13995,7 +14545,7 @@ Triangulation<dim, spacedim>::prepare_coarsening_and_refinement()
             if (cell->refine_flag_set())
               {
                 // loop over neighbors of cell
-                for (const unsigned int i : GeometryInfo<dim>::face_indices())
+                for (const auto i : cell->face_indices())
                   {
                     // only do something if the face is not at the
                     // boundary and if the face will be refined with
index 29463ef5bdc23e3c162e11ba47fa4a5e59bad690..6c9aceaaac917861ffecb92d7505a60dc05a4f0b 100644 (file)
@@ -1496,7 +1496,7 @@ TriaAccessor<structdim, dim, spacedim>::set_bounding_object_indices(
   const ArrayView<int> bounding_object_index_ref =
     this->objects().get_bounding_object_indices(this->present_index);
 
-  AssertDimension(bounding_object_index_ref.size(), new_indices.size());
+  AssertIndexRange(new_indices.size(), bounding_object_index_ref.size() + 1);
 
   unsigned int i = 0;
   for (const auto &new_index : new_indices)
@@ -1516,7 +1516,7 @@ TriaAccessor<structdim, dim, spacedim>::set_bounding_object_indices(
   const ArrayView<int> bounding_object_index_ref =
     this->objects().get_bounding_object_indices(this->present_index);
 
-  AssertDimension(bounding_object_index_ref.size(), new_indices.size());
+  AssertIndexRange(new_indices.size(), bounding_object_index_ref.size() + 1);
 
   unsigned int i = 0;
   for (const auto &new_index : new_indices)
diff --git a/tests/simplex/refinement_01.cc b/tests/simplex/refinement_01.cc
new file mode 100644 (file)
index 0000000..5346fe9
--- /dev/null
@@ -0,0 +1,78 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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 refinement of simplex mesh.
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/simplex/grid_generator.h>
+
+#include "../tests.h"
+
+#include "./simplex_grids.h"
+
+void
+test(const unsigned int v)
+{
+  const unsigned int dim = 2;
+
+  Triangulation<dim> tria;
+
+  if (v < 4)
+    GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+  else
+    GridGenerator::subdivided_hyper_cube_with_simplices_mix(tria, 2);
+
+
+  if (v == 1 || v == 4)
+    {
+      tria.begin()->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+  else if (v == 2 || v == 5)
+    {
+      tria.refine_global(1);
+    }
+  else if (v == 3)
+    {
+      tria.begin()->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+      tria.begin_active(1)->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+
+  GridOut grid_out;
+#if false
+  std::ofstream out("mesh." + std::to_string(v) + ".vtk");
+  grid_out.write_vtk(tria, out);
+#else
+  grid_out.write_vtk(tria, deallog.get_file_stream());
+#endif
+}
+
+int
+main()
+{
+  initlog();
+  test(0); // no refinement
+  test(1); // refinement of a single cell
+  test(2); // global refinement
+  test(3); // refine left bottom corner twice
+  test(4); // mixed mesh: refinement of a single cell
+  test(5); // mixed mesh: global refinement
+}
diff --git a/tests/simplex/refinement_01.with_simplex_support=on.output b/tests/simplex/refinement_01.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..46e83ff
--- /dev/null
@@ -0,0 +1,281 @@
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 4 double
+0.00000 0.00000 0
+1.00000 0.00000 0
+0.00000 1.00000 0
+1.00000 1.00000 0
+
+CELLS 2 8
+3 0 1 2
+3 3 2 1
+
+CELL_TYPES 2
+5 5 
+
+
+
+CELL_DATA 2
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 
+
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 7 double
+0.00000 0.00000 0
+1.00000 0.00000 0
+0.00000 1.00000 0
+1.00000 1.00000 0
+0.500000 0.00000 0
+0.500000 0.500000 0
+0.00000 0.500000 0
+
+CELLS 5 20
+3 3 2 1
+3 0 4 6
+3 4 1 5
+3 6 5 2
+3 4 5 6
+
+CELL_TYPES 5
+5 5 5 5 5 
+
+
+
+CELL_DATA 5
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 0 0 0 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 -1 
+
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 9 double
+0.00000 0.00000 0
+1.00000 0.00000 0
+0.00000 1.00000 0
+1.00000 1.00000 0
+0.500000 0.00000 0
+0.500000 0.500000 0
+1.00000 0.500000 0
+0.00000 0.500000 0
+0.500000 1.00000 0
+
+CELLS 8 32
+3 0 4 7
+3 4 1 5
+3 7 5 2
+3 4 5 7
+3 3 8 6
+3 8 2 5
+3 6 5 1
+3 8 5 6
+
+CELL_TYPES 8
+5 5 5 5 5 5 5 5 
+
+
+
+CELL_DATA 8
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 0 0 0 0 0 0 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 -1 -1 -1 -1 
+
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 10 double
+0.00000 0.00000 0
+1.00000 0.00000 0
+0.00000 1.00000 0
+1.00000 1.00000 0
+0.500000 0.00000 0
+0.500000 0.500000 0
+0.00000 0.500000 0
+0.250000 0.00000 0
+0.00000 0.250000 0
+0.250000 0.250000 0
+
+CELLS 8 32
+3 3 2 1
+3 4 1 5
+3 6 5 2
+3 4 5 6
+3 0 7 8
+3 7 4 9
+3 8 9 6
+3 7 9 8
+
+CELL_TYPES 8
+5 5 5 5 5 5 5 5 
+
+
+
+CELL_DATA 8
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 0 0 0 0 0 0 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-1 -1 -1 -1 -1 -1 -1 -1 
+
+
+# vtk DataFile Version 3.0
+Triangulation generated with deal.II
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 14 double
+0.00000 0.00000 0
+0.500000 0.00000 0
+1.00000 0.00000 0
+0.00000 0.500000 0
+0.500000 0.500000 0
+1.00000 0.500000 0
+0.00000 1.00000 0
+0.500000 1.00000 0
+1.00000 1.00000 0
+0.250000 0.00000 0
+0.00000 0.250000 0
+0.500000 0.250000 0
+0.250000 0.500000 0
+0.250000 0.250000 0
+
+CELLS 10 44
+3 1 2 4
+3 5 4 2
+3 3 4 6
+3 7 6 4
+3 4 5 7
+3 8 7 5
+4 0 9 13 10
+4 9 1 11 13
+4 10 13 12 3
+4 13 11 4 12
+
+CELL_TYPES 10
+5 5 5 5 5 5 9 9 9 9 
+
+
+
+CELL_DATA 10
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 0 0 0 0 0 0 0 0 
+
+
+
+SCALARS ManifoldID int 1
+LOOKUP_TABLE default
+-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 25 double
+0.00000 0.00000 0
+0.500000 0.00000 0
+1.00000 0.00000 0
+0.00000 0.500000 0
+0.500000 0.500000 0
+1.00000 0.500000 0
+0.00000 1.00000 0
+0.500000 1.00000 0
+1.00000 1.00000 0
+0.250000 0.00000 0
+0.00000 0.250000 0
+0.750000 0.00000 0
+0.500000 0.250000 0
+0.750000 0.250000 0
+1.00000 0.250000 0
+0.250000 0.500000 0
+0.250000 0.750000 0
+0.500000 0.750000 0
+0.750000 0.500000 0
+0.750000 0.750000 0
+1.00000 0.750000 0
+0.00000 0.750000 0
+0.250000 1.00000 0
+0.750000 1.00000 0
+0.250000 0.250000 0
+
+CELLS 28 116
+4 0 9 24 10
+4 9 1 12 24
+4 10 24 15 3
+4 24 12 4 15
+3 1 11 12
+3 11 2 13
+3 12 13 4
+3 11 13 12
+3 5 18 14
+3 18 4 13
+3 14 13 2
+3 18 13 14
+3 3 15 21
+3 15 4 16
+3 21 16 6
+3 15 16 21
+3 7 22 17
+3 22 6 16
+3 17 16 4
+3 22 16 17
+3 4 18 17
+3 18 5 19
+3 17 19 7
+3 18 19 17
+3 8 23 20
+3 23 7 19
+3 20 19 5
+3 23 19 20
+
+CELL_TYPES 28
+9 9 9 9 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 
+
+
+
+CELL_DATA 28
+SCALARS MaterialID int 1
+LOOKUP_TABLE default
+0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+
+
+
+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 
+
+

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.