This exposes a few issues with the ExodusII test grids.
--- /dev/null
+Changed: All GridIn functions now remove unused vertices and will attempt to fix
+pyramids and wedges with negative volumes.
+<br>
+(David Wells, 2022/05/25)
// vertices except in 1d
Assert(dim != 1, ExcInternalError());
}
+
+ /**
+ * Apply each of the grid fixup routines in the correct sequence.
+ */
+ template <int dim, int spacedim>
+ void
+ apply_grid_fixup_functions(std::vector<Point<spacedim>> &vertices,
+ std::vector<CellData<dim>> & cells,
+ SubCellData & subcelldata)
+ {
+ // check that no forbidden arrays are used
+ Assert(subcelldata.check_consistency(dim), ExcInternalError());
+ const auto n_hypercube_vertices =
+ ReferenceCells::get_hypercube<dim>().n_vertices();
+ bool is_only_hypercube = true;
+ for (const CellData<dim> &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 <int dim, int spacedim>
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
// 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 &&
// 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,
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;
}
}
- 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,
}
}
- 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);
}
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);
}
;
// 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);
}
}
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);
}
std::vector<CellData<dim>> cells;
SubCellData subcelldata;
std::map<unsigned int, types::boundary_id> boundary_ids_1d;
- bool is_quad_or_hex_mesh = false;
- bool is_tria_or_tet_mesh = false;
{
unsigned int global_cell = 0;
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(
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();
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.
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
}
}
- 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
for (const unsigned int j : GeometryInfo<dim>::vertex_indices())
in >> cells[i].vertices[GeometryInfo<dim>::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);
}
}
}
- 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
ierr = ex_get_ids(ex_id, EX_ELEM_BLOCK, element_block_ids.data());
AssertThrowExodusII(ierr);
- bool is_only_quad_or_hex = true;
-
std::vector<CellData<dim>> cells;
// Elements are grouped together by same reference cell type in element
// blocks. There may be multiple blocks for a single reference cell type,
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.
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);
for (auto &cell : cells)
{
const ArrayView<const unsigned int> 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
}
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.
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
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
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
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
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
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
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
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
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
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 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <string>
+
+#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<Point<dim>> vertices(10);
+ // hypercube
+ vertices[0] = Point<dim>(0, 0, 0);
+ vertices[1] = Point<dim>(1, 0, 0);
+ vertices[2] = Point<dim>(0, 1, 0);
+ vertices[3] = Point<dim>(1, 1, 0);
+ vertices[4] = Point<dim>(0, 0, 1);
+ vertices[5] = Point<dim>(1, 0, 1);
+ vertices[6] = Point<dim>(0, 1, 1);
+ vertices[7] = Point<dim>(1, 1, 1);
+ // pyramid on top
+ vertices[8] = Point<dim>(0.75, 0.5, 1.5);
+ // pyramid on right - a mess of an element
+ vertices[9] = Point<dim>(1.0, 1.25, 0.5);
+
+
+ std::vector<CellData<dim>> 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<Point<dim>> vertices(10);
+ // hypercube
+ vertices[0] = Point<dim>(0, 0, 0);
+ vertices[1] = Point<dim>(1, 0, 0);
+ vertices[2] = Point<dim>(0, 1, 0);
+ vertices[3] = Point<dim>(1, 1, 0);
+ vertices[4] = Point<dim>(0, 0, 1);
+ vertices[5] = Point<dim>(1, 0, 1);
+ vertices[6] = Point<dim>(0, 1, 1);
+ vertices[7] = Point<dim>(1, 1, 1);
+ // pyramid on top
+ vertices[8] = Point<dim>(0.75, 0.5, 1.5);
+ // pyramid on right
+ vertices[9] = Point<dim>(1.25, 0.75, 0.5);
+
+
+ std::vector<CellData<dim>> 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<dim> 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<Point<dim>> vertices(12);
+ // hypercube
+ vertices[0] = Point<dim>(0, 0, 0);
+ vertices[1] = Point<dim>(1, 0, 0);
+ vertices[2] = Point<dim>(0, 1, 0);
+ vertices[3] = Point<dim>(1, 1, 0);
+ vertices[4] = Point<dim>(0, 0, 1);
+ vertices[5] = Point<dim>(1, 0, 1);
+ vertices[6] = Point<dim>(0, 1, 1);
+ vertices[7] = Point<dim>(1, 1, 1);
+ // wedge on top
+ vertices[8] = Point<dim>(0.0, 0.5, 1.5);
+ vertices[9] = Point<dim>(1.0, 0.6, 1.6);
+ // wedge on right
+ vertices[10] = Point<dim>(1.0, 0.5, 0.0);
+ vertices[11] = Point<dim>(1.0, 0.6, 0.9);
+
+
+ std::vector<CellData<dim>> 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<dim> 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();
+}
--- /dev/null
+
+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
+
+
projects / dealii.git / commitdiff