From 52cef3374f71be7ad3982e49e24c3162df45fc58 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 30 Jan 2021 17:47:26 +0100 Subject: [PATCH] Enabling global refinement for TET meshes --- source/grid/tria.cc | 1264 +++++++++++++++++++++++++++- tests/simplex/refinement_04.cc | 72 ++ tests/simplex/refinement_04.output | 357 ++++++++ 3 files changed, 1675 insertions(+), 18 deletions(-) create mode 100644 tests/simplex/refinement_04.cc create mode 100644 tests/simplex/refinement_04.output diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 0c18c0957d..fb50cdca35 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -2445,9 +2445,11 @@ namespace internal // set face orientation if needed if (orientation_needed) - level.face_orientations - [cell * GeometryInfo::faces_per_cell + j] = - connectivity.entity_orientations(dim - 1)[i]; + { + level.face_orientations + [cell * GeometryInfo::faces_per_cell + j] = + connectivity.entity_orientations(dim - 1)[i]; + } } } } @@ -5112,6 +5114,1218 @@ namespace internal } + template + static typename Triangulation<3, spacedim>::DistortedCellList + execute_refinement_isotropic(Triangulation<3, spacedim> &triangulation, + const bool check_for_distorted_cells) + { + static const int dim = 3; + static const unsigned int X = numbers::invalid_unsigned_int; + + Assert(spacedim == 3, ExcNotImplemented()); + + { + 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; + } + } + + + triangulation.faces->quads.clear_user_data(); + + for (typename Triangulation::line_iterator line = + triangulation.begin_line(); + line != triangulation.end_line(); + ++line) + line->clear_user_flag(); + + for (typename Triangulation::quad_iterator quad = + triangulation.begin_quad(); + quad != triangulation.end_quad(); + ++quad) + quad->clear_user_flag(); + + unsigned int needed_vertices = 0; + unsigned int needed_lines_single = 0; + unsigned int needed_quads_single = 0; + unsigned int needed_lines_pair = 0; + unsigned int needed_quads_pair = 0; + for (int level = triangulation.levels.size() - 2; level >= 0; --level) + { + unsigned int new_cells = 0; + + for (const auto &acell : + triangulation.active_cell_iterators_on_level(level)) + if (acell->refine_flag_set()) + { + Assert(acell->refine_flag_set() == + RefinementCase::cut_xyz, + ExcInternalError()); + + if (acell->reference_cell() == ReferenceCells::Hexahedron) + { + ++needed_vertices; + needed_lines_single += 6; + needed_quads_single += 12; + new_cells += 8; + } + else if (acell->reference_cell() == + ReferenceCells::Tetrahedron) + { + needed_lines_single += 1; + needed_quads_single += 8; + new_cells += 8; + } + else + { + Assert(false, ExcInternalError()); + } + + for (const auto face : acell->face_indices()) + if (acell->face(face)->number_of_children() < 4) + acell->face(face)->set_user_flag(); + + for (const auto line : acell->line_indices()) + if (acell->line(line)->has_children() == false) + acell->line(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 + new_cells, + 3, + spacedim); + + reserve_space(triangulation.levels[level + 1]->cells, new_cells); + } + + for (typename Triangulation::quad_iterator quad = + triangulation.begin_quad(); + quad != triangulation.end_quad(); + ++quad) + { + if (quad->user_flag_set() == false) + continue; + + if (quad->reference_cell() == ReferenceCells::Quadrilateral) + { + needed_quads_pair += 4; + needed_lines_pair += 4; + needed_vertices += 1; + } + else if (quad->reference_cell() == ReferenceCells::Triangle) + { + needed_quads_pair += 4; + needed_lines_single += 3; + } + else + { + Assert(false, ExcInternalError()); + } + } + + for (typename Triangulation::line_iterator line = + triangulation.begin_line(); + line != triangulation.end_line(); + ++line) + { + if (line->user_flag_set() == false) + continue; + + needed_lines_pair += 2; + needed_vertices += 1; + } + + reserve_space(triangulation.faces->lines, + needed_lines_pair, + needed_lines_single); + reserve_space(*triangulation.faces, + needed_quads_pair, + needed_quads_single); + reserve_space(triangulation.faces->quads, + needed_quads_pair, + needed_quads_single); + + + 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; + + // LINES + { + 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() == false) + continue; + + 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); + + next_unused_line = + triangulation.faces->lines.template next_free_pair_object<1>( + triangulation); + Assert(next_unused_line.state() == IteratorState::valid, + 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(); + } + } + + // QUADS + { + typename Triangulation::quad_iterator + quad = triangulation.begin_quad(), + endq = triangulation.end_quad(); + typename Triangulation::raw_line_iterator + next_unused_line = triangulation.begin_raw_line(); + typename Triangulation::raw_quad_iterator + next_unused_quad = triangulation.begin_raw_quad(); + + for (; quad != endq; ++quad) + { + if (quad->user_flag_set() == false) + continue; + + const auto &reference_cell_type = quad->reference_cell(); + + // 1) create new vertex (at the center of the face) + if (reference_cell_type == ReferenceCells::Quadrilateral) + { + 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[next_unused_vertex] = + quad->center(true, true); + triangulation.vertices_used[next_unused_vertex] = true; + } + + // 2) create new lines (property is set later) + boost::container::small_vector< + typename Triangulation::raw_line_iterator, + 4> + new_lines(quad->n_lines()); + { + for (unsigned int i = 0; i < new_lines.size(); ++i) + { + if (quad->n_lines() % 2 == 0) + { + Assert(reference_cell_type == + ReferenceCells::Quadrilateral, + ExcNotImplemented()); + if (i % 2 == 0) + next_unused_line = + triangulation.faces->lines + .template next_free_pair_object<1>(triangulation); + } + else + { + Assert(reference_cell_type == ReferenceCells::Triangle, + ExcNotImplemented()); + next_unused_line = + triangulation.faces->lines + .template next_free_single_object<1>(triangulation); + } + + new_lines[i] = next_unused_line; + ++next_unused_line; + + Assert( + new_lines[i]->used() == false, + ExcMessage( + "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be.")); + } + } + + // 3) create new quads (properties are set below) + std::array< + typename Triangulation::raw_quad_iterator, + 4> + new_quads; + { + next_unused_quad = + triangulation.faces->quads.template next_free_pair_object<2>( + triangulation); + + new_quads[0] = next_unused_quad; + Assert( + new_quads[0]->used() == false, + ExcMessage( + "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be.")); + + ++next_unused_quad; + new_quads[1] = next_unused_quad; + Assert( + new_quads[1]->used() == false, + ExcMessage( + "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be.")); + + next_unused_quad = + triangulation.faces->quads.template next_free_pair_object<2>( + triangulation); + new_quads[2] = next_unused_quad; + Assert( + new_quads[2]->used() == false, + ExcMessage( + "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be.")); + + ++next_unused_quad; + new_quads[3] = next_unused_quad; + Assert( + new_quads[3]->used() == false, + ExcMessage( + "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be.")); + + quad->set_children(0, new_quads[0]->index()); + quad->set_children(2, new_quads[2]->index()); + quad->set_refinement_case(RefinementCase<2>::cut_xy); + } + + std::array vertex_indices; + { + for (auto &i : vertex_indices) + i = 0.0; + + unsigned int k = 0; + for (const auto i : quad->vertex_indices()) + vertex_indices[k++] = quad->vertex_index(i); + + for (const auto i : quad->line_indices()) + vertex_indices[k++] = + quad->line(i)->child(0)->vertex_index(1); + + vertex_indices[k++] = next_unused_vertex; + } + + boost::container::small_vector< + typename Triangulation::raw_line_iterator, + 12> + lines(quad->n_lines() * 3); + + { + unsigned int k = 0; + + for (unsigned int l = 0; l < quad->n_lines(); ++l) + for (unsigned int c = 0; c < 2; ++c) + { + static constexpr std::array, + 2> + index = { + {{{1, 0}}, // child 0, line_orientation=false and true + {{0, + 1}}}}; // child 1, line_orientation=false and true + + lines[k++] = quad->line(l)->child( + index[c][quad->line_orientation(l)]); + } + + for (unsigned int l = 0; l < new_lines.size(); ++l) + lines[k++] = new_lines[l]; + } + + boost::container::small_vector line_indices( + lines.size()); + for (unsigned int i = 0; i < line_indices.size(); ++i) + line_indices[i] = lines[i]->index(); + + static constexpr std::array, 12> + line_vertices_quad{{{{0, 4}}, + {{4, 2}}, + {{1, 5}}, + {{5, 3}}, + {{0, 6}}, + {{6, 1}}, + {{2, 7}}, + {{7, 3}}, + {{6, 8}}, + {{8, 7}}, + {{4, 8}}, + {{8, 5}}}}; + + static constexpr std::array, 4> + quad_lines_quad{{{{0, 8, 4, 10}}, + {{8, 2, 5, 11}}, + {{1, 9, 10, 6}}, + {{9, 3, 11, 7}}}}; + + static constexpr std:: + array, 4>, 4> + quad_line_vertices_quad{ + {{{{{0, 4}}, {{6, 8}}, {{0, 6}}, {{4, 8}}}}, + {{{{6, 8}}, {{1, 5}}, {{6, 1}}, {{8, 5}}}}, + {{{{4, 2}}, {{8, 7}}, {{4, 8}}, {{2, 7}}}}, + {{{{8, 7}}, {{5, 3}}, {{8, 5}}, {{7, 3}}}}}}; + + static constexpr std::array, 12> + line_vertices_tri{{{{0, 3}}, + {{3, 1}}, + {{1, 4}}, + {{4, 2}}, + {{2, 5}}, + {{5, 0}}, + {{3, 4}}, + {{4, 5}}, + {{3, 5}}, + {{X, X}}, + {{X, X}}, + {{X, X}}}}; + + static constexpr std::array, 4> + quad_lines_tri{{{{0, 8, 5, X}}, + {{1, 2, 6, X}}, + {{7, 3, 4, X}}, + {{6, 7, 8, X}}}}; + + static constexpr std:: + array, 4>, 4> + quad_line_vertices_tri{ + {{{{{0, 3}}, {{3, 5}}, {{5, 0}}, {{X, X}}}}, + {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}}, + {{{{5, 4}}, {{4, 2}}, {{2, 5}}, {{X, X}}}}, + {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}}}}; + + const auto &line_vertices = + (reference_cell_type == ReferenceCells::Quadrilateral) ? + line_vertices_quad : + line_vertices_tri; + const auto &quad_lines = + (reference_cell_type == ReferenceCells::Quadrilateral) ? + quad_lines_quad : + quad_lines_tri; + const auto &quad_line_vertices = + (reference_cell_type == ReferenceCells::Quadrilateral) ? + quad_line_vertices_quad : + quad_line_vertices_tri; + + // 4) set properties of lines + for (unsigned int i = 0, j = lines.size() - new_lines.size(); + i < new_lines.size(); + ++i, ++j) + { + auto &new_line = new_lines[i]; + new_line->set_bounding_object_indices( + {vertex_indices[line_vertices[j][0]], + vertex_indices[line_vertices[j][1]]}); + new_line->set_used_flag(); + new_line->clear_user_flag(); + new_line->clear_user_data(); + new_line->clear_children(); + new_line->set_boundary_id_internal(quad->boundary_id()); + new_line->set_manifold_id(quad->manifold_id()); + } + + // 5) set properties of quads + for (unsigned int i = 0; i < new_quads.size(); ++i) + { + auto &new_quad = new_quads[i]; + + // TODO: we assume here that all children have the same type + // as the parent + triangulation.faces->quad_reference_cell[new_quad->index()] = + reference_cell_type; + + if (new_quad->n_lines() == 3) + new_quad->set_bounding_object_indices( + {line_indices[quad_lines[i][0]], + line_indices[quad_lines[i][1]], + line_indices[quad_lines[i][2]]}); + else if (new_quad->n_lines() == 4) + new_quad->set_bounding_object_indices( + {line_indices[quad_lines[i][0]], + line_indices[quad_lines[i][1]], + line_indices[quad_lines[i][2]], + line_indices[quad_lines[i][3]]}); + else + Assert(false, ExcNotImplemented()); + + new_quad->set_used_flag(); + new_quad->clear_user_flag(); + new_quad->clear_user_data(); + new_quad->clear_children(); + new_quad->set_boundary_id_internal(quad->boundary_id()); + new_quad->set_manifold_id(quad->manifold_id()); + +#ifdef DEBUG + std::set s; +#endif + + // ... and fix orientation of faces (lines) of quad + for (const auto f : new_quad->line_indices()) + { + std::array vertices_0, vertices_1; + + for (unsigned int v = 0; v < 2; ++v) + vertices_0[v] = + lines[quad_lines[i][f]]->vertex_index(v); + + for (unsigned int v = 0; v < 2; ++v) + vertices_1[v] = + vertex_indices[quad_line_vertices[i][f][v]]; + + const auto orientation = + ReferenceCells::Line.compute_orientation(vertices_0, + vertices_1); + +#ifdef DEBUG + for (const auto i : vertices_0) + s.insert(i); + for (const auto i : vertices_1) + s.insert(i); +#endif + + new_quad->set_line_orientation(f, orientation); + } +#ifdef DEBUG + AssertDimension( + s.size(), + (reference_cell_type == ReferenceCells::Quadrilateral ? 4 : + 3)); +#endif + } + + quad->clear_user_flag(); + } + } + + typename Triangulation<3, spacedim>::DistortedCellList + cells_with_distorted_children; + + for (unsigned int level = 0; level != triangulation.levels.size() - 1; + ++level) + { + typename Triangulation::active_hex_iterator + hex = triangulation.begin_active_hex(level), + endh = triangulation.begin_active_hex(level + 1); + typename Triangulation::raw_hex_iterator + next_unused_hex = triangulation.begin_raw_hex(level + 1); + + for (; hex != endh; ++hex) + { + if (hex->refine_flag_set() == + RefinementCase::no_refinement) + continue; + + const auto &reference_cell_type = hex->reference_cell(); + + const RefinementCase ref_case = hex->refine_flag_set(); + hex->clear_refine_flag(); + hex->set_refinement_case(ref_case); + + unsigned int n_new_lines = 0; + unsigned int n_new_quads = 0; + unsigned int n_new_hexes = 0; + + if (reference_cell_type == ReferenceCells::Hexahedron) + { + n_new_lines = 6; + n_new_quads = 12; + n_new_hexes = 8; + } + else if (reference_cell_type == ReferenceCells::Tetrahedron) + { + n_new_lines = 1; + n_new_quads = 8; + n_new_hexes = 8; + } + else + Assert(false, ExcNotImplemented()); + + if (reference_cell_type == ReferenceCells::Hexahedron) + { + 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] = + hex->center(true, true); + } + + std::vector< + typename Triangulation::raw_line_iterator> + new_lines(n_new_lines); + for (unsigned int i = 0; i < n_new_lines; ++i) + { + new_lines[i] = + triangulation.faces->lines + .template next_free_single_object<1>(triangulation); + + Assert( + new_lines[i]->used() == false, + ExcMessage( + "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be.")); + new_lines[i]->set_used_flag(); + new_lines[i]->clear_user_flag(); + new_lines[i]->clear_user_data(); + new_lines[i]->clear_children(); + new_lines[i]->set_boundary_id_internal( + numbers::internal_face_boundary_id); + new_lines[i]->set_manifold_id(hex->manifold_id()); + } + + std::vector< + typename Triangulation::raw_quad_iterator> + new_quads(n_new_quads); + for (unsigned int i = 0; i < n_new_quads; ++i) + { + new_quads[i] = + triangulation.faces->quads + .template next_free_single_object<2>(triangulation); + + auto &new_quad = new_quads[i]; + + // TODO: faces of children have the same type as the faces + // of the parent + triangulation.faces + ->quad_reference_cell[new_quad->index()] = + (reference_cell_type == ReferenceCells::Hexahedron) ? + ReferenceCells::Quadrilateral : + ReferenceCells::Triangle; + + Assert( + new_quad->used() == false, + ExcMessage( + "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be.")); + new_quad->set_used_flag(); + new_quad->clear_user_flag(); + new_quad->clear_user_data(); + new_quad->clear_children(); + new_quad->set_boundary_id_internal( + numbers::internal_face_boundary_id); + new_quad->set_manifold_id(hex->manifold_id()); + for (const auto j : new_quads[i]->line_indices()) + new_quad->set_line_orientation(j, true); + } + + std::vector< + typename Triangulation::raw_hex_iterator> + new_hexes(n_new_hexes); + { + for (unsigned int i = 0; i < n_new_hexes; ++i) + { + if (i % 2 == 0) + next_unused_hex = + triangulation.levels[level + 1]->cells.next_free_hex( + triangulation, level + 1); + else + ++next_unused_hex; + + new_hexes[i] = next_unused_hex; + + auto &new_hex = new_hexes[i]; + + // TODO: children have the same type as the parent + triangulation.levels[new_hex->level()] + ->reference_cell[new_hex->index()] = + reference_cell_type; + + Assert( + new_hex->used() == false, + ExcMessage( + "Internal error: We want to use a cell during refinement that should be unused, but turns out not to be.")); + new_hex->set_used_flag(); + new_hex->clear_user_flag(); + new_hex->clear_user_data(); + new_hex->clear_children(); + new_hex->set_material_id(hex->material_id()); + new_hex->set_manifold_id(hex->manifold_id()); + new_hex->set_subdomain_id(hex->subdomain_id()); + + if (i % 2) + new_hex->set_parent(hex->index()); + for (const auto f : new_hex->face_indices()) + { + new_hex->set_face_orientation(f, true); + new_hex->set_face_flip(f, false); + new_hex->set_face_rotation(f, false); + } + } + // note these hexes as children to the present cell + for (unsigned int i = 0; i < n_new_hexes / 2; ++i) + hex->set_children(2 * i, new_hexes[2 * i]->index()); + } + + { + // load vertex indices + std::array vertex_indices; + for (auto &i : vertex_indices) + i = 0; + + { + unsigned int k = 0; + + for (const unsigned int i : hex->vertex_indices()) + vertex_indices[k++] = hex->vertex_index(i); + + for (const unsigned int i : hex->line_indices()) + vertex_indices[k++] = + hex->line(i)->child(0)->vertex_index(1); + + if (reference_cell_type == ReferenceCells::Hexahedron) + { + for (const unsigned int i : hex->face_indices()) + vertex_indices[k++] = + middle_vertex_index(hex->face(i)); + + vertex_indices[k++] = next_unused_vertex; + } + } + + // set up new lines + { + static constexpr std::array, 6> + new_line_vertices_hex = {{{{22, 26}}, + {{26, 23}}, + {{20, 26}}, + {{26, 21}}, + {{24, 26}}, + {{26, 25}}}}; + + static constexpr std::array, 6> + new_line_vertices_tet = {{{{6, 8}}, + {{X, X}}, + {{X, X}}, + {{X, X}}, + {{X, X}}, + {{X, X}}}}; + + const auto &new_line_vertices = + (reference_cell_type == ReferenceCells::Hexahedron) ? + new_line_vertices_hex : + new_line_vertices_tet; + + for (unsigned int i = 0; i < new_lines.size(); ++i) + new_lines[i]->set_bounding_object_indices( + {vertex_indices[new_line_vertices[i][0]], + vertex_indices[new_line_vertices[i][1]]}); + } + + // set up new quads + { + boost::container::small_vector< + typename Triangulation::raw_line_iterator, + 30> + relevant_lines(0); + + if (reference_cell_type == ReferenceCells::Hexahedron) + { + relevant_lines.resize(30); + for (unsigned int f = 0, k = 0; f < 6; ++f) + for (unsigned int c = 0; c < 4; ++c, ++k) + { + static constexpr std:: + array, 4> + temp = { + {{{0, 1}}, {{3, 0}}, {{0, 3}}, {{3, 2}}}}; + + relevant_lines[k] = + hex->face(f) + ->isotropic_child( + GeometryInfo:: + standard_to_real_face_vertex( + temp[c][0], + hex->face_orientation(f), + hex->face_flip(f), + hex->face_rotation(f))) + ->line(GeometryInfo:: + standard_to_real_face_line( + temp[c][1], + hex->face_orientation(f), + hex->face_flip(f), + hex->face_rotation(f))); + } + + for (unsigned int i = 0, k = 24; i < 6; ++i, ++k) + relevant_lines[k] = new_lines[i]; + } + else if (reference_cell_type == ReferenceCells::Tetrahedron) + { + relevant_lines.resize(13); + + unsigned int k = 0; + for (unsigned int f = 0; f < 4; ++f) + for (unsigned int l = 0; l < 3; ++l, ++k) + { + // TODO: add comment + static const std:: + array, 6> + table = {{{{1, 0, 2}}, // 0 + {{0, 1, 2}}, + {{0, 2, 1}}, // 2 + {{1, 2, 0}}, + {{2, 1, 0}}, // 4 + {{2, 0, 1}}}}; + + relevant_lines[k] = + hex->face(f) + ->child(3 /*center triangle*/) + ->line( + table[triangulation.levels[hex->level()] + ->face_orientations + [hex->index() * + GeometryInfo< + dim>::faces_per_cell + + f]][l]); + } + + relevant_lines[k++] = new_lines[0]; + + AssertDimension(k, 13); + } + else + Assert(false, ExcNotImplemented()); + + boost::container::small_vector + relevant_line_indices(relevant_lines.size()); + for (unsigned int i = 0; i < relevant_line_indices.size(); + ++i) + relevant_line_indices[i] = relevant_lines[i]->index(); + + static constexpr std::array, 12> + new_quad_lines_hex = {{{{10, 28, 16, 24}}, + {{28, 14, 17, 25}}, + {{11, 29, 24, 20}}, + {{29, 15, 25, 21}}, + {{18, 26, 0, 28}}, + {{26, 22, 1, 29}}, + {{19, 27, 28, 4}}, + {{27, 23, 29, 5}}, + {{2, 24, 8, 26}}, + {{24, 6, 9, 27}}, + {{3, 25, 26, 12}}, + {{25, 7, 27, 13}}}}; + + static constexpr std::array, 12> + new_quad_lines_tet = {{{{2, 3, 8, X}}, + {{0, 9, 5, X}}, + {{1, 6, 11, X}}, + {{4, 10, 7, X}}, + {{2, 12, 5, X}}, + {{1, 9, 12, X}}, + {{4, 8, 12, X}}, + {{6, 12, 10, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}; + + static constexpr std:: + array, 4>, 12> + table_hex = { + {{{{{10, 22}}, {{24, 26}}, {{10, 24}}, {{22, 26}}}}, + {{{{24, 26}}, {{11, 23}}, {{24, 11}}, {{26, 23}}}}, + {{{{22, 14}}, {{26, 25}}, {{22, 26}}, {{14, 25}}}}, + {{{{26, 25}}, {{23, 15}}, {{26, 23}}, {{25, 15}}}}, + {{{{8, 24}}, {{20, 26}}, {{8, 20}}, {{24, 26}}}}, + {{{{20, 26}}, {{12, 25}}, {{20, 12}}, {{26, 25}}}}, + {{{{24, 9}}, {{26, 21}}, {{24, 26}}, {{9, 21}}}}, + {{{{26, 21}}, {{25, 13}}, {{26, 25}}, {{21, 13}}}}, + {{{{16, 20}}, {{22, 26}}, {{16, 22}}, {{20, 26}}}}, + {{{{22, 26}}, {{17, 21}}, {{22, 17}}, {{26, 21}}}}, + {{{{20, 18}}, {{26, 23}}, {{20, 26}}, {{18, 23}}}}, + {{{{26, 23}}, {{21, 19}}, {{26, 21}}, {{23, 19}}}}}}; + + static constexpr std:: + array, 4>, 12> + table_tet = { + {{{{{6, 4}}, {{4, 7}}, {{7, 6}}, {{X, X}}}}, + {{{{4, 5}}, {{5, 8}}, {{8, 4}}, {{X, X}}}}, + {{{{5, 6}}, {{6, 9}}, {{9, 5}}, {{X, X}}}}, + {{{{7, 8}}, {{8, 9}}, {{9, 7}}, {{X, X}}}}, + {{{{4, 6}}, {{6, 8}}, {{8, 4}}, {{X, X}}}}, + {{{{6, 5}}, {{5, 8}}, {{8, 6}}, {{X, X}}}}, + {{{{8, 7}}, {{7, 6}}, {{6, 8}}, {{X, X}}}}, + {{{{9, 6}}, {{6, 8}}, {{8, 9}}, {{X, X}}}}, + {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}, + {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}, + {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}, + {{{{X, X}}, {{X, X}}, {{X, X}}, {{X, X}}}}}}; + + const auto &new_quad_lines = + (reference_cell_type == ReferenceCells::Hexahedron) ? + new_quad_lines_hex : + new_quad_lines_tet; + + const auto &table = + (reference_cell_type == ReferenceCells::Hexahedron) ? + table_hex : + table_tet; + + for (unsigned int q = 0; q < new_quads.size(); ++q) + { + for (unsigned int l = 0; l < 3; ++l) + { + std::array vertices_0, vertices_1; + + for (unsigned int v = 0; v < 2; ++v) + vertices_0[v] = + relevant_lines[new_quad_lines[q][l]] + ->vertex_index(v); + + for (unsigned int v = 0; v < 2; ++v) + vertices_1[v] = vertex_indices[table[q][l][v]]; + } + } + + for (unsigned int q = 0; q < new_quads.size(); ++q) + { + auto &new_quad = new_quads[q]; + + if (new_quad->n_lines() == 3) + new_quad->set_bounding_object_indices( + {relevant_line_indices[new_quad_lines[q][0]], + relevant_line_indices[new_quad_lines[q][1]], + relevant_line_indices[new_quad_lines[q][2]]}); + else if (new_quad->n_lines() == 4) + new_quad->set_bounding_object_indices( + {relevant_line_indices[new_quad_lines[q][0]], + relevant_line_indices[new_quad_lines[q][1]], + relevant_line_indices[new_quad_lines[q][2]], + relevant_line_indices[new_quad_lines[q][3]]}); + else + Assert(false, ExcNotImplemented()); + + for (const auto l : new_quad->line_indices()) + { + std::array vertices_0, vertices_1; + + for (unsigned int v = 0; v < 2; ++v) + vertices_0[v] = + relevant_lines[new_quad_lines[q][l]] + ->vertex_index(v); + + for (unsigned int v = 0; v < 2; ++v) + vertices_1[v] = vertex_indices[table[q][l][v]]; + + const auto orientation = + ReferenceCells::Line.compute_orientation( + vertices_0, vertices_1); + + new_quad->set_line_orientation(l, orientation); + } + } + } + + // set up new hex + { + std::array quad_indices; + + if (reference_cell_type == ReferenceCells::Hexahedron) + { + for (unsigned int i = 0; i < new_quads.size(); ++i) + quad_indices[i] = new_quads[i]->index(); + + for (unsigned int f = 0, k = new_quads.size(); f < 6; + ++f) + for (unsigned int c = 0; c < 4; ++c, ++k) + quad_indices[k] = + hex->face(f)->isotropic_child_index( + GeometryInfo::standard_to_real_face_vertex( + c, + hex->face_orientation(f), + hex->face_flip(f), + hex->face_rotation(f))); + } + else if (reference_cell_type == ReferenceCells::Tetrahedron) + { + for (unsigned int i = 0; i < new_quads.size(); ++i) + quad_indices[i] = new_quads[i]->index(); + + for (unsigned int f = 0, k = new_quads.size(); f < 4; + ++f) + for (unsigned int c = 0; c < 4; ++c, ++k) + { + quad_indices[k] = hex->face(f)->child_index( + (c == 3) ? + 3 : + reference_cell_type + .standard_to_real_face_vertex( + c, + f, + triangulation.levels[hex->level()] + ->face_orientations + [hex->index() * + GeometryInfo::faces_per_cell + + f])); + } + } + else + { + Assert(false, ExcNotImplemented()); + } + + static constexpr std::array, 8> + cell_quads_hex = {{ + {{12, 0, 20, 4, 28, 8}}, // bottom children + {{0, 16, 22, 6, 29, 9}}, // + {{13, 1, 4, 24, 30, 10}}, // + {{1, 17, 6, 26, 31, 11}}, // + {{14, 2, 21, 5, 8, 32}}, // top children + {{2, 18, 23, 7, 9, 33}}, // + {{15, 3, 5, 25, 10, 34}}, // + {{3, 19, 7, 27, 11, 35}} // + }}; + + static constexpr std::array, 8> + cell_quads_tet{{{{8, 13, 16, 0, X, X}}, + {{9, 12, 1, 21, X, X}}, + {{10, 2, 17, 20, X, X}}, + {{3, 14, 18, 22, X, X}}, + {{11, 1, 4, 5, X, X}}, + {{15, 0, 4, 6, X, X}}, + {{19, 7, 6, 3, X, X}}, + {{23, 5, 2, 7, X, X}}}}; + + static constexpr std:: + array, 6>, 8> + cell_face_vertices_hex{{{{{{0, 8, 16, 20}}, + {{10, 24, 22, 26}}, + {{0, 16, 10, 22}}, + {{8, 20, 24, 26}}, + {{0, 10, 8, 24}}, + {{16, 22, 20, 26}}}}, + {{{{10, 24, 22, 26}}, + {{1, 9, 17, 21}}, + {{10, 22, 1, 17}}, + {{24, 26, 9, 21}}, + {{10, 1, 24, 9}}, + {{22, 17, 26, 21}}}}, + {{{{8, 2, 20, 18}}, + {{24, 11, 26, 23}}, + {{8, 20, 24, 26}}, + {{2, 18, 11, 23}}, + {{8, 24, 2, 11}}, + {{20, 26, 18, 23}}}}, + {{{{24, 11, 26, 23}}, + {{9, 3, 21, 19}}, + {{24, 26, 9, 21}}, + {{11, 23, 3, 19}}, + {{24, 9, 11, 3}}, + {{26, 21, 23, 19}}}}, + {{{{16, 20, 4, 12}}, + {{22, 26, 14, 25}}, + {{16, 4, 22, 14}}, + {{20, 12, 26, 25}}, + {{16, 22, 20, 26}}, + {{4, 14, 12, 25}}}}, + {{{{22, 26, 14, 25}}, + {{17, 21, 5, 13}}, + {{22, 14, 17, 5}}, + {{26, 25, 21, 13}}, + {{22, 17, 26, 21}}, + {{14, 5, 25, 13}}}}, + {{{{20, 18, 12, 6}}, + {{26, 23, 25, 15}}, + {{20, 12, 26, 25}}, + {{18, 6, 23, 15}}, + {{20, 26, 18, 23}}, + {{12, 25, 6, 15}}}}, + {{{{26, 23, 25, 15}}, + {{21, 19, 13, 7}}, + {{26, 25, 21, 13}}, + {{23, 15, 19, 7}}, + {{26, 21, 23, 19}}, + {{25, 13, 15, 7}}}}}}; + + static constexpr std:: + array, 6>, 8> + cell_face_vertices_tet{{{{{{0, 4, 6, X}}, + {{4, 0, 7, X}}, + {{0, 6, 7, X}}, + {{6, 4, 7, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{4, 1, 5, X}}, + {{1, 4, 8, X}}, + {{4, 5, 8, X}}, + {{5, 1, 8, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{6, 5, 2, X}}, + {{5, 6, 9, X}}, + {{6, 2, 9, X}}, + {{2, 5, 9, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{7, 8, 9, X}}, + {{8, 7, 3, X}}, + {{7, 9, 3, X}}, + {{9, 8, 3, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{4, 5, 6, X}}, + {{5, 4, 8, X}}, + {{4, 6, 8, X}}, + {{6, 5, 8, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{4, 7, 8, X}}, + {{7, 4, 6, X}}, + {{4, 8, 6, X}}, + {{8, 7, 6, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{6, 9, 7, X}}, + {{9, 6, 8, X}}, + {{6, 7, 8, X}}, + {{7, 9, 8, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}, + {{{{5, 8, 9, X}}, + {{8, 5, 6, X}}, + {{5, 9, 6, X}}, + {{9, 8, 6, X}}, + {{X, X, X, X}}, + {{X, X, X, X}}}}}}; + + const auto &cell_quads = + (reference_cell_type == ReferenceCells::Hexahedron) ? + cell_quads_hex : + cell_quads_tet; + + const auto &cell_face_vertices = + (reference_cell_type == ReferenceCells::Hexahedron) ? + cell_face_vertices_hex : + cell_face_vertices_tet; + + for (unsigned int c = 0; c < 8; ++c) + { + auto &new_hex = new_hexes[c]; + + if (new_hex->n_faces() == 4) + new_hex->set_bounding_object_indices( + {quad_indices[cell_quads[c][0]], + quad_indices[cell_quads[c][1]], + quad_indices[cell_quads[c][2]], + quad_indices[cell_quads[c][3]]}); + else if (new_hex->n_faces() == 6) + new_hex->set_bounding_object_indices( + {quad_indices[cell_quads[c][0]], + quad_indices[cell_quads[c][1]], + quad_indices[cell_quads[c][2]], + quad_indices[cell_quads[c][3]], + quad_indices[cell_quads[c][4]], + quad_indices[cell_quads[c][5]]}); + else + Assert(false, ExcNotImplemented()); + + for (const auto f : new_hex->face_indices()) + { + std::array vertices_0, vertices_1; + + const auto &face = new_hex->face(f); + + for (const auto i : face->vertex_indices()) + vertices_0[i] = face->vertex_index(i); + + for (const auto i : face->vertex_indices()) + vertices_1[i] = + vertex_indices[cell_face_vertices[c][f][i]]; + + const auto orientation = + face->reference_cell().compute_orientation( + vertices_1, vertices_0); + + new_hex->set_face_orientation( + f, Utilities::get_bit(orientation, 0)); + new_hex->set_face_flip( + f, Utilities::get_bit(orientation, 2)); + new_hex->set_face_rotation( + f, Utilities::get_bit(orientation, 1)); + } + } + } + } + + if (check_for_distorted_cells && + has_distorted_children(hex)) + cells_with_distorted_children.distorted_cells.push_back(hex); + + triangulation.signals.post_refinement_on_cell(hex); + } + } + + triangulation.faces->quads.clear_user_data(); + + return cells_with_distorted_children; + } + /** * A function that performs the refinement of a triangulation in * 3d. @@ -5123,6 +6337,30 @@ namespace internal { const unsigned int dim = 3; + { + bool flag_isotropic_mesh = true; + typename Triangulation::raw_cell_iterator + cell = triangulation.begin(), + endc = triangulation.end(); + for (; cell != endc; ++cell) + if (cell->used()) + if (triangulation.get_anisotropic_refinement_flag() || + cell->refine_flag_set() == RefinementCase::cut_x || + cell->refine_flag_set() == RefinementCase::cut_y || + cell->refine_flag_set() == RefinementCase::cut_z || + cell->refine_flag_set() == RefinementCase::cut_xy || + cell->refine_flag_set() == RefinementCase::cut_xz || + cell->refine_flag_set() == RefinementCase::cut_yz) + { + flag_isotropic_mesh = false; + break; + } + + if (flag_isotropic_mesh) + return execute_refinement_isotropic(triangulation, + check_for_distorted_cells); + } + // this function probably also works for spacedim>3 but it // isn't tested. it will probably be necessary to pull new // vertices onto the manifold just as we do for the other @@ -9635,9 +10873,7 @@ namespace internal else if (cell->has_children() && !cell->child(0)->coarsen_flag_set()) { - for (unsigned int line = 0; - line < GeometryInfo::lines_per_cell; - ++line) + for (unsigned int line = 0; line < cell->n_lines(); ++line) if (GeometryInfo::line_refinement_case( cell->refinement_case(), line) == RefinementCase<1>::cut_x) @@ -9695,9 +10931,7 @@ namespace internal else cell->set_refine_flag(); - for (unsigned int l = 0; - l < GeometryInfo::lines_per_cell; - ++l) + for (unsigned int l = 0; l < cell->n_lines(); ++l) if (GeometryInfo::line_refinement_case( cell->refine_flag_set(), line) == RefinementCase<1>::cut_x) @@ -9711,9 +10945,7 @@ namespace internal // iterations if we flag all lines of // this cell now (and not at the outset // of the next iteration) for refinement - for (unsigned int l = 0; - l < GeometryInfo::lines_per_cell; - ++l) + for (unsigned int l = 0; l < cell->n_lines(); ++l) if (!cell->line(l)->has_children() && (GeometryInfo::line_refinement_case( cell->refine_flag_set(), l) != @@ -9749,9 +10981,7 @@ namespace internal --cell) { if (cell->user_flag_set()) - for (unsigned int line = 0; - line < GeometryInfo::lines_per_cell; - ++line) + for (unsigned int line = 0; line < cell->n_lines(); ++line) if (cell->line(line)->has_children() && (cell->line(line)->child(0)->user_flag_set() || cell->line(line)->child(1)->user_flag_set())) @@ -9759,9 +10989,7 @@ namespace internal for (unsigned int c = 0; c < cell->n_children(); ++c) cell->child(c)->clear_coarsen_flag(); cell->clear_user_flag(); - for (unsigned int l = 0; - l < GeometryInfo::lines_per_cell; - ++l) + for (unsigned int l = 0; l < cell->n_lines(); ++l) if (GeometryInfo::line_refinement_case( cell->refinement_case(), l) == RefinementCase<1>::cut_x) diff --git a/tests/simplex/refinement_04.cc b/tests/simplex/refinement_04.cc new file mode 100644 index 0000000000..466a1789c5 --- /dev/null +++ b/tests/simplex/refinement_04.cc @@ -0,0 +1,72 @@ +// --------------------------------------------------------------------- +// +// 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 3D simplex mesh. + +#include +#include +#include + +#include "../tests.h" + +template +void +print(Triangulation &tria, const std::string &label) +{ + GridOut grid_out; + + unsigned int counter = 0; + + for (const auto &cell : tria.active_cell_iterators()) + cell->set_material_id(counter++); + +#if false + std::ofstream out(label); + grid_out.write_vtk(tria, out); +#else + (void)label; + grid_out.write_vtk(tria, deallog.get_file_stream()); +#endif +} + +void +test(const unsigned int version) +{ + Triangulation<3> tria; + + if (version == 0 || version == 1) + GridGenerator::reference_cell(tria, ReferenceCells::Tetrahedron); + else if (version == 2) + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); + + print(tria, "tria.0.vtk"); + if (version == 0 || version == 2) + tria.refine_global(1); + else if (version == 1) + tria.refine_global(2); + + print(tria, "tria.1.vtk"); +} + +int +main() +{ + initlog(); + test(0); // coarse grid: 1 tet, n_refinements: 1 + test(1); // coarse grid: 1 tet, n_refinements: 2 + test(2); // coarse grid: 5 tets (cube), n_refinements: 1 +} diff --git a/tests/simplex/refinement_04.output b/tests/simplex/refinement_04.output new file mode 100644 index 0000000000..f9bfb8d951 --- /dev/null +++ b/tests/simplex/refinement_04.output @@ -0,0 +1,357 @@ + +# vtk DataFile Version 3.0 +Triangulation generated with deal.II +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 4 double +0.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +0.00000 1.00000 0.00000 +0.00000 0.00000 1.00000 + +CELLS 1 5 +4 0 1 2 3 + +CELL_TYPES 1 +10 + + + +CELL_DATA 1 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +0 + + + +SCALARS ManifoldID int 1 +LOOKUP_TABLE default +-1 + + +# vtk DataFile Version 3.0 +Triangulation generated with deal.II +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 10 double +0.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +0.00000 1.00000 0.00000 +0.00000 0.00000 1.00000 +0.500000 0.00000 0.00000 +0.00000 0.00000 0.500000 +0.500000 0.500000 0.00000 +0.500000 0.00000 0.500000 +0.00000 0.500000 0.00000 +0.00000 0.500000 0.500000 + +CELLS 8 40 +4 0 4 8 5 +4 4 1 6 7 +4 8 6 2 9 +4 5 7 9 3 +4 4 6 8 7 +4 4 5 7 8 +4 8 9 5 7 +4 6 7 9 8 + +CELL_TYPES 8 +10 10 10 10 10 10 10 10 + + + +CELL_DATA 8 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +0 1 2 3 4 5 6 7 + + + +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 4 double +0.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +0.00000 1.00000 0.00000 +0.00000 0.00000 1.00000 + +CELLS 1 5 +4 0 1 2 3 + +CELL_TYPES 1 +10 + + + +CELL_DATA 1 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +0 + + + +SCALARS ManifoldID int 1 +LOOKUP_TABLE default +-1 + + +# vtk DataFile Version 3.0 +Triangulation generated with deal.II +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 35 double +0.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +0.00000 1.00000 0.00000 +0.00000 0.00000 1.00000 +0.500000 0.00000 0.00000 +0.00000 0.00000 0.500000 +0.500000 0.500000 0.00000 +0.500000 0.00000 0.500000 +0.00000 0.500000 0.00000 +0.00000 0.500000 0.500000 +0.250000 0.00000 0.00000 +0.750000 0.00000 0.00000 +0.00000 0.00000 0.250000 +0.00000 0.00000 0.750000 +0.750000 0.250000 0.00000 +0.250000 0.750000 0.00000 +0.750000 0.00000 0.250000 +0.250000 0.00000 0.750000 +0.00000 0.750000 0.00000 +0.00000 0.250000 0.00000 +0.00000 0.750000 0.250000 +0.00000 0.250000 0.750000 +0.250000 0.250000 0.250000 +0.00000 0.250000 0.250000 +0.00000 0.250000 0.500000 +0.00000 0.500000 0.250000 +0.250000 0.500000 0.250000 +0.250000 0.250000 0.500000 +0.500000 0.250000 0.250000 +0.250000 0.250000 0.00000 +0.250000 0.500000 0.00000 +0.500000 0.250000 0.00000 +0.500000 0.00000 0.250000 +0.250000 0.00000 0.500000 +0.250000 0.00000 0.250000 + +CELLS 64 320 +4 0 10 19 12 +4 10 4 29 34 +4 19 29 8 23 +4 12 34 23 5 +4 10 29 19 34 +4 10 12 34 19 +4 19 23 12 34 +4 29 34 23 19 +4 4 11 31 32 +4 11 1 14 16 +4 31 14 6 28 +4 32 16 28 7 +4 11 14 31 16 +4 11 32 16 31 +4 31 28 32 16 +4 14 16 28 31 +4 8 30 18 25 +4 30 6 15 26 +4 18 15 2 20 +4 25 26 20 9 +4 30 15 18 26 +4 30 25 26 18 +4 18 20 25 26 +4 15 26 20 18 +4 5 33 24 13 +4 33 7 27 17 +4 24 27 9 21 +4 13 17 21 3 +4 33 27 24 17 +4 33 13 17 24 +4 24 21 13 17 +4 27 17 21 24 +4 4 31 29 32 +4 31 6 30 28 +4 29 30 8 22 +4 32 28 22 7 +4 31 30 29 28 +4 31 32 28 29 +4 29 22 32 28 +4 30 28 22 29 +4 4 34 32 29 +4 34 5 33 23 +4 32 33 7 22 +4 29 23 22 8 +4 34 33 32 23 +4 34 29 23 32 +4 32 22 29 23 +4 33 23 22 32 +4 8 25 23 22 +4 25 9 24 27 +4 23 24 5 33 +4 22 27 33 7 +4 25 24 23 27 +4 25 22 27 23 +4 23 33 22 27 +4 24 27 33 23 +4 6 28 26 30 +4 28 7 27 22 +4 26 27 9 25 +4 30 22 25 8 +4 28 27 26 22 +4 28 30 22 26 +4 26 25 30 22 +4 27 22 25 26 + +CELL_TYPES 64 +10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 + + + +CELL_DATA 64 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 + + + +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 -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 -1 -1 -1 -1 + + +# vtk DataFile Version 3.0 +Triangulation generated with deal.II +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 8 double +0.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +0.00000 1.00000 0.00000 +1.00000 1.00000 0.00000 +0.00000 0.00000 1.00000 +1.00000 0.00000 1.00000 +0.00000 1.00000 1.00000 +1.00000 1.00000 1.00000 + +CELLS 5 25 +4 0 1 2 4 +4 2 1 3 7 +4 1 4 5 7 +4 2 4 7 6 +4 1 2 4 7 + +CELL_TYPES 5 +10 10 10 10 10 + + + +CELL_DATA 5 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +0 1 2 3 4 + + + +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 26 double +0.00000 0.00000 0.00000 +1.00000 0.00000 0.00000 +0.00000 1.00000 0.00000 +1.00000 1.00000 0.00000 +0.00000 0.00000 1.00000 +1.00000 0.00000 1.00000 +0.00000 1.00000 1.00000 +1.00000 1.00000 1.00000 +0.500000 0.00000 0.00000 +0.00000 0.00000 0.500000 +0.500000 0.500000 0.00000 +1.00000 0.500000 0.00000 +0.500000 0.00000 0.500000 +1.00000 0.500000 0.500000 +0.00000 0.500000 0.00000 +0.00000 0.500000 0.500000 +0.00000 1.00000 0.500000 +0.500000 1.00000 0.500000 +0.500000 1.00000 0.00000 +1.00000 1.00000 0.500000 +0.500000 0.00000 1.00000 +0.00000 0.500000 1.00000 +0.500000 0.500000 1.00000 +1.00000 0.00000 0.500000 +1.00000 0.500000 1.00000 +0.500000 1.00000 1.00000 + +CELLS 40 200 +4 0 8 14 9 +4 8 1 10 12 +4 14 10 2 15 +4 9 12 15 4 +4 8 10 14 12 +4 8 9 12 14 +4 14 15 9 12 +4 10 12 15 14 +4 2 10 18 17 +4 10 1 11 13 +4 18 11 3 19 +4 17 13 19 7 +4 10 11 18 13 +4 10 17 13 18 +4 18 19 17 13 +4 11 13 19 18 +4 1 12 23 13 +4 12 4 20 22 +4 23 20 5 24 +4 13 22 24 7 +4 12 20 23 22 +4 12 13 22 23 +4 23 24 13 22 +4 20 22 24 23 +4 2 15 17 16 +4 15 4 22 21 +4 17 22 7 25 +4 16 21 25 6 +4 15 22 17 21 +4 15 16 21 17 +4 17 25 16 21 +4 22 21 25 17 +4 1 10 12 13 +4 10 2 15 17 +4 12 15 4 22 +4 13 17 22 7 +4 10 15 12 17 +4 10 13 17 12 +4 12 22 13 17 +4 15 17 22 12 + +CELL_TYPES 40 +10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 + + + +CELL_DATA 40 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 + + + +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 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 + + -- 2.39.5