From f654ef024dc029585c1e708cf7df9db722f153b0 Mon Sep 17 00:00:00 2001 From: David Wells Date: Wed, 25 May 2022 10:43:17 -0400 Subject: [PATCH] Refactor the way we apply grid fixup functions. This exposes a few issues with the ExodusII test grids. --- .../incompatibilities/20220325DavidWells | 4 + source/grid/grid_in.cc | 201 ++++-------------- source/grid/grid_tools.cc | 8 +- ...odusii.with_trilinos_with_seacas=on.output | 84 +++----- tests/simplex/negative_measure_01.cc | 162 ++++++++++++++ tests/simplex/negative_measure_01.output | 88 ++++++++ 6 files changed, 338 insertions(+), 209 deletions(-) create mode 100644 doc/news/changes/incompatibilities/20220325DavidWells create mode 100644 tests/simplex/negative_measure_01.cc create mode 100644 tests/simplex/negative_measure_01.output diff --git a/doc/news/changes/incompatibilities/20220325DavidWells b/doc/news/changes/incompatibilities/20220325DavidWells new file mode 100644 index 0000000000..acb0802f6c --- /dev/null +++ b/doc/news/changes/incompatibilities/20220325DavidWells @@ -0,0 +1,4 @@ +Changed: All GridIn functions now remove unused vertices and will attempt to fix +pyramids and wedges with negative volumes. +
+(David Wells, 2022/05/25) diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc index b9a8edae92..882d1b0bda 100644 --- a/source/grid/grid_in.cc +++ b/source/grid/grid_in.cc @@ -97,6 +97,35 @@ namespace // vertices except in 1d Assert(dim != 1, ExcInternalError()); } + + /** + * Apply each of the grid fixup routines in the correct sequence. + */ + template + void + apply_grid_fixup_functions(std::vector> &vertices, + std::vector> & cells, + SubCellData & subcelldata) + { + // check that no forbidden arrays are used + Assert(subcelldata.check_consistency(dim), ExcInternalError()); + const auto n_hypercube_vertices = + ReferenceCells::get_hypercube().n_vertices(); + bool is_only_hypercube = true; + for (const CellData &cell : cells) + if (cell.vertices.size() != n_hypercube_vertices) + { + is_only_hypercube = false; + break; + } + + GridTools::delete_unused_vertices(vertices, cells, subcelldata); + if (dim == spacedim) + GridTools::invert_cells_with_negative_measure(vertices, cells); + + if (is_only_hypercube) + GridTools::consistently_order_cells(cells); + } } // namespace template @@ -195,9 +224,6 @@ GridIn::read_vtk(std::istream &in) unsigned int n_geometric_objects = 0; unsigned int n_ints; - bool is_quad_or_hex_mesh = false; - bool is_tria_or_tet_mesh = false; - if (keyword == "CELLS") { // jump to the `CELL_TYPES` section and read in cell types @@ -235,11 +261,6 @@ GridIn::read_vtk(std::istream &in) // VTK_TETRA is 10, VTK_HEXAHEDRON is 12 if (cell_types[count] == 10 || cell_types[count] == 12) { - if (cell_types[count] == 10) - is_tria_or_tet_mesh = true; - if (cell_types[count] == 12) - is_quad_or_hex_mesh = true; - // we assume that the file contains first all cells, // and only then any faces or lines AssertThrow(subcelldata.boundary_quads.size() == 0 && @@ -267,11 +288,6 @@ GridIn::read_vtk(std::istream &in) // VTK_TRIANGLE is 5, VTK_QUAD is 9 else if (cell_types[count] == 5 || cell_types[count] == 9) { - if (cell_types[count] == 5) - is_tria_or_tet_mesh = true; - if (cell_types[count] == 9) - is_quad_or_hex_mesh = true; - // we assume that the file contains first all cells, // then all faces, and finally all lines AssertThrow(subcelldata.boundary_lines.size() == 0, @@ -319,11 +335,6 @@ GridIn::read_vtk(std::istream &in) AssertThrow(subcelldata.boundary_lines.size() == 0, ExcNotImplemented()); - if (cell_types[count] == 5) - is_tria_or_tet_mesh = true; - if (cell_types[count] == 9) - is_quad_or_hex_mesh = true; - cells.emplace_back(n_vertices); for (unsigned int j = 0; j < n_vertices; @@ -556,29 +567,8 @@ GridIn::read_vtk(std::istream &in) } } - Assert(subcelldata.check_consistency(dim), ExcInternalError()); - - - // TODO: the functions below (GridTools::delete_unused_vertices(), - // GridTools::invert_all_negative_measure_cells(), - // GridTools::consistently_order_cells()) need to be - // revisited for simplex/mixed meshes - - if (dim == 1 || (is_quad_or_hex_mesh && !is_tria_or_tet_mesh)) - { - GridTools::delete_unused_vertices(vertices, cells, subcelldata); - - if (dim == spacedim) - GridTools::invert_all_negative_measure_cells(vertices, cells); - - GridTools::consistently_order_cells(cells); - tria->create_triangulation(vertices, cells, subcelldata); - } - else - { - // simplex or mixed mesh - tria->create_triangulation(vertices, cells, subcelldata); - } + apply_grid_fixup_functions(vertices, cells, subcelldata); + tria->create_triangulation(vertices, cells, subcelldata); } else AssertThrow(false, @@ -879,15 +869,7 @@ GridIn::read_unv(std::istream &in) } } - Assert(subcelldata.check_consistency(dim), ExcInternalError()); - - GridTools::delete_unused_vertices(vertices, cells, subcelldata); - - if (dim == spacedim) - GridTools::invert_all_negative_measure_cells(vertices, cells); - - GridTools::consistently_order_cells(cells); - + apply_grid_fixup_functions(vertices, cells, subcelldata); tria->create_triangulation(vertices, cells, subcelldata); } @@ -1108,18 +1090,9 @@ GridIn::read_ucd(std::istream &in, AssertThrow(false, ExcUnknownIdentifier(cell_type)); } - - // check that no forbidden arrays are used - Assert(subcelldata.check_consistency(dim), ExcInternalError()); - AssertThrow(in.fail() == false, ExcIO()); - // do some clean-up on vertices... - GridTools::delete_unused_vertices(vertices, cells, subcelldata); - // ... and cells - if (dim == spacedim) - GridTools::invert_all_negative_measure_cells(vertices, cells); - GridTools::consistently_order_cells(cells); + apply_grid_fixup_functions(vertices, cells, subcelldata); tria->create_triangulation(vertices, cells, subcelldata); } @@ -1361,18 +1334,9 @@ GridIn::read_dbmesh(std::istream &in) ; // ok, so we are not at the end of // the file, that's it, mostly - - - // check that no forbidden arrays are used - Assert(subcelldata.check_consistency(dim), ExcInternalError()); - AssertThrow(in.fail() == false, ExcIO()); - // do some clean-up on vertices... - GridTools::delete_unused_vertices(vertices, cells, subcelldata); - // ...and cells - GridTools::invert_all_negative_measure_cells(vertices, cells); - GridTools::consistently_order_cells(cells); + apply_grid_fixup_functions(vertices, cells, subcelldata); tria->create_triangulation(vertices, cells, subcelldata); } @@ -1436,11 +1400,7 @@ GridIn::read_xda(std::istream &in) } AssertThrow(in.fail() == false, ExcIO()); - // do some clean-up on vertices... - GridTools::delete_unused_vertices(vertices, cells, subcelldata); - // ... and cells - GridTools::invert_all_negative_measure_cells(vertices, cells); - GridTools::consistently_order_cells(cells); + apply_grid_fixup_functions(vertices, cells, subcelldata); tria->create_triangulation(vertices, cells, subcelldata); } @@ -2313,8 +2273,6 @@ GridIn::read_msh(std::istream &in) std::vector> cells; SubCellData subcelldata; std::map boundary_ids_1d; - bool is_quad_or_hex_mesh = false; - bool is_tria_or_tet_mesh = false; { unsigned int global_cell = 0; @@ -2465,25 +2423,13 @@ GridIn::read_msh(std::istream &in) if (cell_type == 1) // line vertices_per_cell = 2; else if (cell_type == 2) // tri - { - vertices_per_cell = 3; - is_tria_or_tet_mesh = true; - } + vertices_per_cell = 3; else if (cell_type == 3) // quad - { - vertices_per_cell = 4; - is_quad_or_hex_mesh = true; - } + vertices_per_cell = 4; else if (cell_type == 4) // tet - { - vertices_per_cell = 4; - is_tria_or_tet_mesh = true; - } + vertices_per_cell = 4; else if (cell_type == 5) // hex - { - vertices_per_cell = 8; - is_quad_or_hex_mesh = true; - } + vertices_per_cell = 8; AssertThrow(nod_num == vertices_per_cell, ExcMessage( @@ -2582,15 +2528,9 @@ GridIn::read_msh(std::istream &in) unsigned int vertices_per_cell = 0; // check cell type if (cell_type == 2) // tri - { - vertices_per_cell = 3; - is_tria_or_tet_mesh = true; - } + vertices_per_cell = 3; else if (cell_type == 3) // quad - { - vertices_per_cell = 4; - is_quad_or_hex_mesh = true; - } + vertices_per_cell = 4; subcelldata.boundary_quads.emplace_back(); @@ -2666,10 +2606,6 @@ GridIn::read_msh(std::istream &in) static const std::string end_elements_marker[] = {"$ENDELM", "$EndElements"}; AssertThrow(line == end_elements_marker[gmsh_file_format == 10 ? 0 : 1], ExcInvalidGMSHInput(line)); - - // check that no forbidden arrays are used - Assert(subcelldata.check_consistency(dim), ExcInternalError()); - AssertThrow(in.fail() == false, ExcIO()); // check that we actually read some cells. @@ -2677,25 +2613,7 @@ GridIn::read_msh(std::istream &in) ExcGmshNoCellInformation(subcelldata.boundary_lines.size(), subcelldata.boundary_quads.size())); - // TODO: the functions below (GridTools::delete_unused_vertices(), - // GridTools::invert_all_negative_measure_cells(), - // GridTools::consistently_order_cells()) need to be revisited - // for simplex/mixed meshes - - if (dim == 1 || (is_quad_or_hex_mesh && !is_tria_or_tet_mesh)) - { - // do some clean-up on vertices... - GridTools::delete_unused_vertices(vertices, cells, subcelldata); - // ... and cells - if (dim == spacedim) - GridTools::invert_cells_with_negative_measure(vertices, cells); - GridTools::consistently_order_cells(cells); - } - else if (is_tria_or_tet_mesh) - { - if (dim == spacedim) - GridTools::invert_cells_with_negative_measure(vertices, cells); - } + apply_grid_fixup_functions(vertices, cells, subcelldata); tria->create_triangulation(vertices, cells, subcelldata); // in 1d, we also have to attach boundary ids to vertices, which does not @@ -2903,8 +2821,7 @@ GridIn::read_msh(const std::string &fname) } } - Assert(subcelldata.check_consistency(dim), ExcInternalError()); - + apply_grid_fixup_functions(vertices, cells, subcelldata); tria->create_triangulation(vertices, cells, subcelldata); // in 1d, we also have to attach boundary ids to vertices, which does not @@ -3375,19 +3292,10 @@ GridIn<2>::read_tecplot(std::istream &in) for (const unsigned int j : GeometryInfo::vertex_indices()) in >> cells[i].vertices[GeometryInfo::ucd_to_deal[j]]; } - // do some clean-up on vertices - GridTools::delete_unused_vertices(vertices, cells, subcelldata); } - - // check that no forbidden arrays are - // used. as we do not read in any - // subcelldata, nothing should happen here. - Assert(subcelldata.check_consistency(dim), ExcInternalError()); AssertThrow(in.fail() == false, ExcIO()); - // do some cleanup on cells - GridTools::invert_all_negative_measure_cells(vertices, cells); - GridTools::consistently_order_cells(cells); + apply_grid_fixup_functions(vertices, cells, subcelldata); tria->create_triangulation(vertices, cells, subcelldata); } @@ -3539,10 +3447,7 @@ GridIn::read_assimp(const std::string &filename, } } - GridTools::delete_unused_vertices(vertices, cells, subcelldata); - if (dim == spacedim) - GridTools::invert_all_negative_measure_cells(vertices, cells); - GridTools::consistently_order_cells(cells); + apply_grid_fixup_functions(vertices, cells, subcelldata); tria->create_triangulation(vertices, cells, subcelldata); #else @@ -3842,8 +3747,6 @@ GridIn::read_exodusii( ierr = ex_get_ids(ex_id, EX_ELEM_BLOCK, element_block_ids.data()); AssertThrowExodusII(ierr); - bool is_only_quad_or_hex = true; - std::vector> cells; // Elements are grouped together by same reference cell type in element // blocks. There may be multiple blocks for a single reference cell type, @@ -3871,9 +3774,6 @@ GridIn::read_exodusii( const ReferenceCell type = exodusii_name_to_type(string_temp.data(), n_nodes_per_element); - if (type.is_simplex()) - is_only_quad_or_hex = false; - // The number of nodes per element may be larger than what we want to // read - for example, if the Exodus file contains a QUAD9 element, we // only want to read the first four values and ignore the rest. @@ -3908,16 +3808,7 @@ GridIn::read_exodusii( ierr = ex_close(ex_id); AssertThrowExodusII(ierr); - if (is_only_quad_or_hex) - { - // do some clean-up on vertices... - GridTools::delete_unused_vertices(vertices, cells, pair.first); - // ... and cells - if (dim == spacedim) - GridTools::invert_all_negative_measure_cells(vertices, cells); - GridTools::consistently_order_cells(cells); - } - + apply_grid_fixup_functions(vertices, cells, pair.first); tria->create_triangulation(vertices, cells, pair.first); ExodusIIData out; out.id_to_sideset_ids = std::move(pair.second); diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index c9da816fc3..f6f642ab1f 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -885,15 +885,16 @@ namespace GridTools for (auto &cell : cells) { const ArrayView vertices(cell.vertices); - if (GridTools::cell_measure(all_vertices, vertices) < 0) + // Some pathologically twisted cells can have exactly zero measure but + // we can still fix them + if (GridTools::cell_measure(all_vertices, vertices) <= 0) { + ++n_negative_cells; const auto reference_cell = ReferenceCell::n_vertices_to_type(dim, vertices.size()); if (reference_cell.is_hyper_cube()) { - ++n_negative_cells; - if (dim == 2) { // flip the cell across the y = x line in 2D @@ -910,7 +911,6 @@ namespace GridTools } else if (reference_cell.is_simplex()) { - ++n_negative_cells; // By basic rules for computing determinants we can just swap // two vertices to fix a negative volume. Arbitrarily pick the // last two. diff --git a/tests/grid/grid_in_exodusii.with_trilinos_with_seacas=on.output b/tests/grid/grid_in_exodusii.with_trilinos_with_seacas=on.output index 14d6585b06..b1b43d02d4 100644 --- a/tests/grid/grid_in_exodusii.with_trilinos_with_seacas=on.output +++ b/tests/grid/grid_in_exodusii.with_trilinos_with_seacas=on.output @@ -123,18 +123,18 @@ DEAL::face center = 1.50000 3.00000 1.50000 face boundary id = 0 DEAL::face center = 1.33333 0.333333 3.66667 face boundary id = 0 DEAL::face center = 0.333333 1.33333 3.66667 face boundary id = 0 DEAL::face center = 2.33333 0.666667 4.33333 face boundary id = 0 -DEAL::face center = 5.50000 0.500000 4.00000 face boundary id = 0 -DEAL::face center = 5.50000 1.50000 3.00000 face boundary id = 0 DEAL::face center = 1.33333 2.33333 3.66667 face boundary id = 0 DEAL::face center = 2.33333 1.66667 4.33333 face boundary id = 0 -DEAL::face center = 5.50000 2.00000 4.00000 face boundary id = 0 +DEAL::face center = 5.50000 0.500000 4.00000 face boundary id = 0 +DEAL::face center = 5.50000 1.50000 3.00000 face boundary id = 0 DEAL::face center = 8.00000 1.33333 3.66667 face boundary id = 0 +DEAL::face center = 5.50000 2.00000 4.00000 face boundary id = 0 DEAL::Number of vertices: 13 DEAL::Number of cells: 4 DEAL::cell 0 type = Hex volume = 27.0000 DEAL::cell 1 type = Pyramid volume = 6.00000 DEAL::cell 2 type = Tet volume = 2.00000 -DEAL::cell 3 type = Wedge volume = -15.0000 +DEAL::cell 3 type = Wedge volume = 15.0000 # vtk DataFile Version 3.0 Triangulation generated with deal.II ASCII @@ -158,7 +158,7 @@ CELLS 4 27 8 0 1 2 3 4 5 6 7 5 4 5 6 7 8 4 5 9 6 8 -6 6 5 9 11 10 12 +6 11 10 12 6 5 9 CELL_TYPES 4 12 14 10 13 @@ -189,23 +189,23 @@ DEAL::face center = 1.50000 3.00000 1.50000 face boundary id = 0 DEAL::face center = 1.33333 0.333333 3.66667 face boundary id = 0 DEAL::face center = 0.333333 1.33333 3.66667 face boundary id = 0 DEAL::face center = 2.33333 0.666667 4.33333 face boundary id = 0 -DEAL::face center = 5.50000 0.500000 4.00000 face boundary id = 0 -DEAL::face center = 5.50000 1.50000 3.00000 face boundary id = 0 DEAL::face center = 1.33333 2.33333 3.66667 face boundary id = 0 DEAL::face center = 2.33333 1.66667 4.33333 face boundary id = 0 -DEAL::face center = 5.50000 2.00000 4.00000 face boundary id = 0 +DEAL::face center = 5.50000 0.500000 4.00000 face boundary id = 0 +DEAL::face center = 5.50000 1.50000 3.00000 face boundary id = 0 DEAL::face center = 8.00000 1.33333 3.66667 face boundary id = 0 -DEAL::Number of vertices: 29 +DEAL::face center = 5.50000 2.00000 4.00000 face boundary id = 0 +DEAL::Number of vertices: 13 DEAL::Number of cells: 4 DEAL::cell 0 type = Hex volume = 27.0000 DEAL::cell 1 type = Pyramid volume = 6.00000 DEAL::cell 2 type = Tet volume = 2.00000 -DEAL::cell 3 type = Wedge volume = -15.0000 +DEAL::cell 3 type = Wedge volume = 15.0000 # vtk DataFile Version 3.0 Triangulation generated with deal.II ASCII DATASET UNSTRUCTURED_GRID -POINTS 29 double +POINTS 13 double 0.0 0.0 0.0 3.0 0.0 0.0 3.0 3.0 0.0 @@ -219,28 +219,12 @@ POINTS 29 double 8.0 0.0 3.0 8.0 3.0 3.0 8.0 1.0 5.0 -3.0 1.5 0.0 -1.5 3.0 0.0 -0.0 1.5 0.0 -1.5 0.0 0.0 -3.0 0.0 1.5 -3.0 3.0 1.5 -0.0 3.0 1.5 -0.0 0.0 1.5 -3.0 1.5 3.0 -1.5 3.0 3.0 -0.0 1.5 3.0 -1.5 0.0 3.0 -0.50 0.50 4.0 -2.0 0.50 4.0 -2.0 2.0 4.0 -0.50 2.0 4.0 CELLS 4 27 8 0 1 2 3 4 5 6 7 5 4 5 6 7 8 4 5 9 6 8 -6 6 5 9 11 10 12 +6 11 10 12 6 5 9 CELL_TYPES 4 12 14 10 13 @@ -277,18 +261,18 @@ DEAL::face center = 1.50000 3.00000 1.50000 face boundary id = 4 DEAL::face center = 1.33333 0.333333 3.66667 face boundary id = 1 DEAL::face center = 0.333333 1.33333 3.66667 face boundary id = 5 DEAL::face center = 2.33333 0.666667 4.33333 face boundary id = 1 -DEAL::face center = 5.50000 0.500000 4.00000 face boundary id = 1 -DEAL::face center = 5.50000 1.50000 3.00000 face boundary id = 4 DEAL::face center = 1.33333 2.33333 3.66667 face boundary id = 4 DEAL::face center = 2.33333 1.66667 4.33333 face boundary id = 3 -DEAL::face center = 5.50000 2.00000 4.00000 face boundary id = 3 +DEAL::face center = 5.50000 0.500000 4.00000 face boundary id = 1 +DEAL::face center = 5.50000 1.50000 3.00000 face boundary id = 4 DEAL::face center = 8.00000 1.33333 3.66667 face boundary id = 6 +DEAL::face center = 5.50000 2.00000 4.00000 face boundary id = 3 DEAL::Number of vertices: 13 DEAL::Number of cells: 4 DEAL::cell 0 type = Hex volume = 27.0000 DEAL::cell 1 type = Pyramid volume = 6.00000 DEAL::cell 2 type = Tet volume = 2.00000 -DEAL::cell 3 type = Wedge volume = -15.0000 +DEAL::cell 3 type = Wedge volume = 15.0000 # vtk DataFile Version 3.0 Triangulation generated with deal.II ASCII @@ -312,7 +296,7 @@ CELLS 18 91 8 0 1 2 3 4 5 6 7 5 4 5 6 7 8 4 5 9 6 8 -6 6 5 9 11 10 12 +6 11 10 12 6 5 9 4 0 4 5 1 4 0 1 2 3 4 0 3 7 4 @@ -321,23 +305,23 @@ CELLS 18 91 3 5 4 8 3 4 7 8 3 9 5 8 -4 5 9 12 10 -4 6 5 10 11 3 7 6 8 3 6 9 8 -4 9 6 11 12 -3 11 10 12 +4 10 12 9 5 +4 11 10 5 6 +3 10 11 12 +4 12 11 6 9 CELL_TYPES 18 12 14 10 13 -9 9 9 9 9 5 5 5 9 9 5 5 9 5 +9 9 9 9 9 5 5 5 5 5 9 9 5 9 CELL_DATA 18 SCALARS MaterialID int 1 LOOKUP_TABLE default 1 10 100 1000 -1 2 5 3 4 1 5 1 1 4 4 3 3 6 +1 2 5 3 4 1 5 1 4 3 1 4 6 3 SCALARS ManifoldID int 1 @@ -363,18 +347,18 @@ DEAL::face center = 1.50000 3.00000 1.50000 face manifold id = 4 DEAL::face center = 1.33333 0.333333 3.66667 face manifold id = 1 DEAL::face center = 0.333333 1.33333 3.66667 face manifold id = 5 DEAL::face center = 2.33333 0.666667 4.33333 face manifold id = 1 -DEAL::face center = 5.50000 0.500000 4.00000 face manifold id = 1 -DEAL::face center = 5.50000 1.50000 3.00000 face manifold id = 4 DEAL::face center = 1.33333 2.33333 3.66667 face manifold id = 4 DEAL::face center = 2.33333 1.66667 4.33333 face manifold id = 3 -DEAL::face center = 5.50000 2.00000 4.00000 face manifold id = 3 +DEAL::face center = 5.50000 0.500000 4.00000 face manifold id = 1 +DEAL::face center = 5.50000 1.50000 3.00000 face manifold id = 4 DEAL::face center = 8.00000 1.33333 3.66667 face manifold id = 6 +DEAL::face center = 5.50000 2.00000 4.00000 face manifold id = 3 DEAL::Number of vertices: 13 DEAL::Number of cells: 4 DEAL::cell 0 type = Hex volume = 27.0000 DEAL::cell 1 type = Pyramid volume = 6.00000 DEAL::cell 2 type = Tet volume = 2.00000 -DEAL::cell 3 type = Wedge volume = -15.0000 +DEAL::cell 3 type = Wedge volume = 15.0000 # vtk DataFile Version 3.0 Triangulation generated with deal.II ASCII @@ -398,7 +382,7 @@ CELLS 18 91 8 0 1 2 3 4 5 6 7 5 4 5 6 7 8 4 5 9 6 8 -6 6 5 9 11 10 12 +6 11 10 12 6 5 9 4 0 4 5 1 4 0 1 2 3 4 0 3 7 4 @@ -407,16 +391,16 @@ CELLS 18 91 3 5 4 8 3 4 7 8 3 9 5 8 -4 5 9 12 10 -4 6 5 10 11 3 7 6 8 3 6 9 8 -4 9 6 11 12 -3 11 10 12 +4 10 12 9 5 +4 11 10 5 6 +3 10 11 12 +4 12 11 6 9 CELL_TYPES 18 12 14 10 13 -9 9 9 9 9 5 5 5 9 9 5 5 9 5 +9 9 9 9 9 5 5 5 5 5 9 9 5 9 CELL_DATA 18 @@ -429,6 +413,6 @@ LOOKUP_TABLE default SCALARS ManifoldID int 1 LOOKUP_TABLE default -1 -1 -1 -1 -1 2 5 3 4 1 5 1 1 4 4 3 3 6 +1 2 5 3 4 1 5 1 4 3 1 4 6 3 DEAL::OK diff --git a/tests/simplex/negative_measure_01.cc b/tests/simplex/negative_measure_01.cc new file mode 100644 index 0000000000..f3c40974e5 --- /dev/null +++ b/tests/simplex/negative_measure_01.cc @@ -0,0 +1,162 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2002 - 2022 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 GridTools::invert_cells_with_negative_measure() with a mixed mesh. + +#include +#include + +#include + +#include "../tests.h" + +void +test() +{ + deal_II_exceptions::disable_abort_on_exception(); + + // test a grid with a pyramid that cannot be fixed + const static int dim = 3; + try + { + std::vector> vertices(10); + // hypercube + vertices[0] = Point(0, 0, 0); + vertices[1] = Point(1, 0, 0); + vertices[2] = Point(0, 1, 0); + vertices[3] = Point(1, 1, 0); + vertices[4] = Point(0, 0, 1); + vertices[5] = Point(1, 0, 1); + vertices[6] = Point(0, 1, 1); + vertices[7] = Point(1, 1, 1); + // pyramid on top + vertices[8] = Point(0.75, 0.5, 1.5); + // pyramid on right - a mess of an element + vertices[9] = Point(1.0, 1.25, 0.5); + + + std::vector> cells(3); + cells[0].vertices = {0, 1, 2, 3, 4, 5, 6, 7}; + // this one is twisted + cells[1].vertices = {4, 5, 7, 6, 8}; + // this one is twisted in an unfixable way + cells[2].vertices = {1, 3, 5, 7, 9}; + + SubCellData subcelldata; + const std::size_t n_cells_inverted = + GridTools::invert_cells_with_negative_measure(vertices, cells); + + deallog << "We inverted " << n_cells_inverted << " cell(s)." << std::endl; + } + catch (const ExcGridHasInvalidCell &exception) + { + deallog << "Successfully failed to set up a mesh with a pyramid " + << "pointing into another element" << std::endl; + } + + // test with a fixable pyramid + { + std::vector> vertices(10); + // hypercube + vertices[0] = Point(0, 0, 0); + vertices[1] = Point(1, 0, 0); + vertices[2] = Point(0, 1, 0); + vertices[3] = Point(1, 1, 0); + vertices[4] = Point(0, 0, 1); + vertices[5] = Point(1, 0, 1); + vertices[6] = Point(0, 1, 1); + vertices[7] = Point(1, 1, 1); + // pyramid on top + vertices[8] = Point(0.75, 0.5, 1.5); + // pyramid on right + vertices[9] = Point(1.25, 0.75, 0.5); + + + std::vector> cells(3); + cells[0].vertices = {0, 1, 2, 3, 4, 5, 6, 7}; + // this one is twisted + cells[1].vertices = {4, 5, 7, 6, 8}; + // this one is not twisted + cells[2].vertices = {5, 1, 7, 3, 9}; + + SubCellData subcelldata; + const std::size_t n_cells_inverted = + GridTools::invert_cells_with_negative_measure(vertices, cells); + + deallog << "We inverted " << n_cells_inverted << " cell(s)." << std::endl; + + Triangulation tria; + tria.create_triangulation(vertices, cells, subcelldata); + + std::ostream &logfile = deallog.get_file_stream(); + logfile << "---------------------------------------------" << std::endl + << std::endl + << std::endl; + GridOut grid_out; + grid_out.write_vtk(tria, logfile); + } + + // test with a twisted and untwisted wedge + { + std::vector> vertices(12); + // hypercube + vertices[0] = Point(0, 0, 0); + vertices[1] = Point(1, 0, 0); + vertices[2] = Point(0, 1, 0); + vertices[3] = Point(1, 1, 0); + vertices[4] = Point(0, 0, 1); + vertices[5] = Point(1, 0, 1); + vertices[6] = Point(0, 1, 1); + vertices[7] = Point(1, 1, 1); + // wedge on top + vertices[8] = Point(0.0, 0.5, 1.5); + vertices[9] = Point(1.0, 0.6, 1.6); + // wedge on right + vertices[10] = Point(1.0, 0.5, 0.0); + vertices[11] = Point(1.0, 0.6, 0.9); + + + std::vector> cells(3); + cells[0].vertices = {0, 1, 2, 3, 4, 5, 6, 7}; + // this one is twisted + cells[1].vertices = {5, 7, 9, 4, 6, 8}; + // this one is not twisted + cells[2].vertices = {1, 10, 3, 5, 11, 7}; + + SubCellData subcelldata; + const std::size_t n_cells_inverted = + GridTools::invert_cells_with_negative_measure(vertices, cells); + + deallog << "We inverted " << n_cells_inverted << " cell(s)." << std::endl; + + Triangulation tria; + tria.create_triangulation(vertices, cells, subcelldata); + + std::ostream &logfile = deallog.get_file_stream(); + logfile << "---------------------------------------------" << std::endl + << std::endl + << std::endl; + GridOut grid_out; + grid_out.write_vtk(tria, logfile); + } +} + +int +main() +{ + initlog(); + test(); +} diff --git a/tests/simplex/negative_measure_01.output b/tests/simplex/negative_measure_01.output new file mode 100644 index 0000000000..20e21352e6 --- /dev/null +++ b/tests/simplex/negative_measure_01.output @@ -0,0 +1,88 @@ + +DEAL::Successfully failed to set up a mesh with a pyramid pointing into another element +DEAL::We inverted 1 cell(s). +--------------------------------------------- + + +# 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 +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.750000 0.500000 1.50000 +1.25000 0.750000 0.500000 + +CELLS 3 21 +8 0 1 3 2 4 5 7 6 +5 4 5 7 6 8 +5 5 1 3 7 9 + +CELL_TYPES 3 +12 14 14 + + + +CELL_DATA 3 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +0 0 0 + + + +SCALARS ManifoldID int 1 +LOOKUP_TABLE default +-1 -1 -1 + + +DEAL::We inverted 1 cell(s). +--------------------------------------------- + + +# vtk DataFile Version 3.0 +Triangulation generated with deal.II +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 12 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.00000 0.500000 1.50000 +1.00000 0.600000 1.60000 +1.00000 0.500000 0.00000 +1.00000 0.600000 0.900000 + +CELLS 3 23 +8 0 1 3 2 4 5 7 6 +6 4 6 8 5 7 9 +6 1 10 3 5 11 7 + +CELL_TYPES 3 +12 13 13 + + + +CELL_DATA 3 +SCALARS MaterialID int 1 +LOOKUP_TABLE default +0 0 0 + + + +SCALARS ManifoldID int 1 +LOOKUP_TABLE default +-1 -1 -1 + + -- 2.39.5