From: Peter Munch Date: Thu, 12 Nov 2020 17:24:28 +0000 (+0100) Subject: Simplex mesh (2D): implement refinement X-Git-Tag: v9.3.0-rc1~800^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11178%2Fhead;p=dealii.git Simplex mesh (2D): implement refinement --- diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index b2166d1424..0af2ce4e39 100644 --- a/include/deal.II/grid/tria_accessor.h +++ b/include/deal.II/grid/tria_accessor.h @@ -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; }; diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h index 07624ba4a4..79779bf023 100644 --- a/include/deal.II/grid/tria_accessor.templates.h +++ b/include/deal.II/grid/tria_accessor.templates.h @@ -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 */ diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 68778d1a82..d9d2b7737a 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -423,7 +423,7 @@ namespace std::vector line_cell_count(triangulation.n_raw_lines(), 0); for (const auto &cell : triangulation.cell_iterators()) - for (unsigned int l = 0; l < GeometryInfo::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 bool has_distorted_children( - const typename Triangulation::cell_iterator &cell, - std::integral_constant, - std::integral_constant) + const typename Triangulation::cell_iterator &cell) { Assert(cell->has_children(), ExcInternalError()); @@ -694,9 +692,7 @@ namespace template bool has_distorted_children( - const typename Triangulation::cell_iterator &, - std::integral_constant, - std::integral_constant) + const typename Triangulation::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::face_indices()) + for (auto f : cell->face_indices()) { const typename Triangulation::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::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 + static typename Triangulation::DistortedCellList + execute_refinement_isotropic(Triangulation &triangulation, + const bool check_for_distorted_cells) + { + AssertDimension(dim, 2); + + { + typename Triangulation::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::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()); + triangulation.vertices_used.resize(needed_vertices, false); + } + + unsigned int next_unused_vertex = 0; + + { + typename Triangulation::active_line_iterator + line = triangulation.begin_active_line(), + endl = triangulation.end_line(); + typename Triangulation::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::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::DistortedCellList + cells_with_distorted_children; + + typename Triangulation::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 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::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(new_vertices[vertex_no]) || + cell->line(face_no)->child(0)->vertex_index(1) == + static_cast(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::raw_cell_iterator + subcells[GeometryInfo::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(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(triangulation.levels.size()) - 1; + ++level) + { + typename Triangulation::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(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::raw_cell_iterator + cell = triangulation.begin(), + endc = triangulation.end(); + for (; cell != endc; ++cell) + if (cell->used()) + if (cell->refine_flag_set() == RefinementCase::cut_x || + cell->refine_flag_set() == RefinementCase::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(), - std::integral_constant())) + has_distorted_children(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(), - std::integral_constant())) + has_distorted_children(hex)) cells_with_distorted_children.distorted_cells.push_back( hex); @@ -9379,11 +9933,8 @@ namespace internal execute_refinement(Triangulation &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 @@ -9391,7 +9942,7 @@ namespace internal prevent_distorted_boundary_cells( Triangulation &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 &triangulation) { - AssertThrow(false, ExcNotImplemented()); - (void)triangulation; + Implementation::prepare_refinement_dim_dependent(triangulation); } template @@ -13995,7 +14545,7 @@ Triangulation::prepare_coarsening_and_refinement() if (cell->refine_flag_set()) { // loop over neighbors of cell - for (const unsigned int i : GeometryInfo::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 diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 29463ef5bd..6c9aceaaac 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1496,7 +1496,7 @@ TriaAccessor::set_bounding_object_indices( const ArrayView 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::set_bounding_object_indices( const ArrayView 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 index 0000000000..5346fe9a19 --- /dev/null +++ b/tests/simplex/refinement_01.cc @@ -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 +#include + +#include + +#include "../tests.h" + +#include "./simplex_grids.h" + +void +test(const unsigned int v) +{ + const unsigned int dim = 2; + + Triangulation 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 index 0000000000..46e83ffb56 --- /dev/null +++ b/tests/simplex/refinement_01.with_simplex_support=on.output @@ -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 + +