const std::vector<Point<spacedim>> &all_vertices,
std::vector<CellData<dim>> & cells);
+ /**
+ * Check the given cells and inverts any cell that is considered to have
+ * negative measure/volume in the orientation required by deal.II.
+ *
+ * This function is identical to invert_all_negative_measure_cells() except it
+ * does not throw an error if only some of the cells are inverted. Instead,
+ * this function returns how many cells were inverted. Additionally, it will
+ * always throw an exception outside of codimension 0.
+ */
+ template <int dim, int spacedim>
+ std::size_t
+ invert_cells_with_negative_measure(
+ const std::vector<Point<spacedim>> &all_vertices,
+ std::vector<CellData<dim>> & cells);
+
/**
* Given a vector of CellData objects describing a mesh, reorder their
* vertices so that all lines are consistently oriented.
GridTools::delete_unused_vertices(vertices, cells, subcelldata);
// ... and cells
if (dim == spacedim)
- GridTools::invert_all_negative_measure_cells(vertices, cells);
+ GridTools::invert_cells_with_negative_measure(vertices, cells);
GridTools::consistently_order_cells(cells);
}
tria->create_triangulation(vertices, cells, subcelldata);
template <int dim, int spacedim>
- void
- invert_all_negative_measure_cells(
+ std::size_t
+ invert_cells_with_negative_measure(
const std::vector<Point<spacedim>> &all_vertices,
std::vector<CellData<dim>> & cells)
{
if (dim == 1)
- return;
+ return 0;
if (dim == 2 && spacedim == 3)
Assert(false, ExcNotImplemented());
ExcInternalError());
}
}
+ return n_negative_cells;
+ }
+
+
+ template <int dim, int spacedim>
+ void
+ invert_all_negative_measure_cells(
+ const std::vector<Point<spacedim>> &all_vertices,
+ std::vector<CellData<dim>> & cells)
+ {
+ const std::size_t n_negative_cells =
+ invert_cells_with_negative_measure(all_vertices, cells);
// We assume that all cells of a grid have
// either positive or negative volumes but
const std::vector<Point<deal_II_space_dimension>> &,
std::vector<CellData<deal_II_dimension>> &);
+ template std::size_t
+ invert_cells_with_negative_measure(
+ const std::vector<Point<deal_II_space_dimension>> &,
+ std::vector<CellData<deal_II_dimension>> &);
+
# if deal_II_dimension == deal_II_space_dimension
template void
consistently_order_cells(std::vector<CellData<deal_II_dimension>> &);
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 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 loading a .msh file from gmsh with some cells having negative volume
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <fstream>
+#include <iostream>
+#include <map>
+
+using namespace dealii;
+
+
+template <int dim>
+void
+print_mesh_info(const Triangulation<dim> &triangulation,
+ const std::string & filename)
+{
+ deallog << "Mesh info:" << std::endl
+ << " dimension: " << dim << std::endl
+ << " no. of cells: " << triangulation.n_active_cells() << std::endl;
+
+ {
+ std::map<types::boundary_id, unsigned int> boundary_count;
+ for (const auto &face : triangulation.active_face_iterators())
+ if (face->at_boundary())
+ boundary_count[face->boundary_id()]++;
+
+ deallog << " boundary indicators: ";
+ for (const std::pair<const types::boundary_id, unsigned int> &pair :
+ boundary_count)
+ {
+ deallog << pair.first << "(" << pair.second << " times) ";
+ }
+ deallog << std::endl;
+ }
+
+ // Finally, produce a graphical representation of the mesh to an output
+ // file:
+ std::ofstream out(filename);
+ GridOut grid_out;
+ grid_out.write_vtu(triangulation, out);
+ deallog << " written to " << filename << std::endl << std::endl;
+}
+
+
+void
+grid()
+{
+ Triangulation<2> triangulation;
+
+ GridIn<2> gridin;
+ gridin.attach_triangulation(triangulation);
+ std::ifstream f(SOURCE_DIR "/grid_in_gmsh_02.msh");
+ gridin.read_msh(f);
+
+ print_mesh_info(triangulation, "grid.vtu");
+}
+
+
+int
+main()
+{
+ initlog();
+ grid();
+}
--- /dev/null
+meshsize = 2;
+
+Point(1) = {-1, 0, 0, meshsize};
+Point(2) = {1, 0, 0, meshsize};
+Point(3) = {1, 2, 0, meshsize};
+Point(4) = {-1, 2, 0, meshsize};
+Point(5) = {-1, 4, 0, meshsize};
+Point(6) = {1, 4, 0, meshsize};
+Line(1) = {5, 6};
+Line(2) = {4, 3};
+Line(3) = {1, 2};
+Line(4) = {5, 4};
+Line(5) = {4, 1};
+Line(6) = {6, 3};
+Line(7) = {3, 2};
+
+Line Loop(9) = {1, 6, -2, -4};
+Plane Surface(9) = {9};
+Line Loop(11) = {3, -7, -2, 5};
+Plane Surface(11) = {11};
+
+Physical Curve("1") = {1};
+Physical Curve("2") = {3};
+Physical Surface(1) = {9, 11};
+
+Mesh.SubdivisionAlgorithm = 1;
--- /dev/null
+$MeshFormat
+4.1 0 8
+$EndMeshFormat
+$PhysicalNames
+2
+1 1 "1"
+1 2 "2"
+$EndPhysicalNames
+$Entities
+6 7 2 0
+1 -1 0 0 0
+2 1 0 0 0
+3 1 2 0 0
+4 -1 2 0 0
+5 -1 4 0 0
+6 1 4 0 0
+1 -1 4 0 1 4 0 1 1 2 5 -6
+2 -1 2 0 1 2 0 0 2 4 -3
+3 -1 0 0 1 0 0 1 2 2 1 -2
+4 -1 2 0 -1 4 0 0 2 5 -4
+5 -1 0 0 -1 2 0 0 2 4 -1
+6 1 2 0 1 4 0 0 2 6 -3
+7 1 0 0 1 2 0 0 2 3 -2
+9 -1 2 0 1 4 0 1 1 4 1 6 -2 -4
+11 -1 0 0 1 2 0 1 1 4 3 -7 -2 5
+$EndEntities
+$Nodes
+15 31 1 31
+0 1 0 1
+1
+-1 0 0
+0 2 0 1
+2
+1 0 0
+0 3 0 1
+3
+1 2 0
+0 4 0 1
+4
+-1 2 0
+0 5 0 1
+5
+-1 4 0
+0 6 0 1
+6
+1 4 0
+1 1 0 1
+7
+-2.625677453238495e-12 4 0
+1 2 0 1
+8
+-2.625677453238495e-12 2 0
+1 3 0 1
+9
+-2.625677453238495e-12 0 0
+1 4 0 1
+10
+-1 3.000000000002758 0
+1 5 0 1
+11
+-1 1.000000000002661 0
+1 6 0 1
+12
+1 3.000000000002758 0
+1 7 0 1
+13
+1 1.000000000002661 0
+2 9 0 9
+14
+15
+16
+17
+18
+19
+20
+21
+22
+0 3 0
+-0.5 2.5 0
+0.5 2.5 0
+0.5 3.5 0
+-0.5 3.5 0
+-4.376082829438133e-13 2.333333333333333 0
+0.6666666666666667 3.00000000000046 0
+-0.6666666666666667 3.00000000000046 0
+-4.376082829438133e-13 3.666666666666667 0
+2 11 0 9
+23
+24
+25
+26
+27
+28
+29
+30
+31
+0 1 0
+0.5 0.5 0
+-0.5 0.5 0
+-0.5 1.5 0
+0.5 1.5 0
+-4.376082829438133e-13 0.3333333333333333 0
+-0.6666666666666667 1.000000000000444 0
+0.6666666666666667 1.000000000000443 0
+-4.376082829438133e-13 1.666666666666667 0
+$EndNodes
+$Elements
+4 28 1 28
+1 1 1 2
+1 5 7
+2 7 6
+1 3 1 2
+3 1 9
+4 9 2
+2 9 3 12
+5 3 8 19 16
+6 8 4 15 19
+7 16 19 15 14
+8 6 12 20 17
+9 12 3 16 20
+10 17 20 16 14
+11 4 10 21 15
+12 10 5 18 21
+13 15 21 18 14
+14 5 7 22 18
+15 7 6 17 22
+16 18 22 17 14
+2 11 3 12
+17 1 9 28 25
+18 9 2 24 28
+19 25 28 24 23
+20 4 11 29 26
+21 11 1 25 29
+22 26 29 25 23
+23 2 13 30 24
+24 13 3 27 30
+25 24 30 27 23
+26 3 8 31 27
+27 8 4 26 31
+28 27 31 26 23
+$EndElements
--- /dev/null
+
+DEAL::Mesh info:
+DEAL:: dimension: 2
+DEAL:: no. of cells: 24
+DEAL:: boundary indicators: 0(8 times) 1(2 times) 2(2 times)
+DEAL:: written to grid.vtu
+DEAL::
projects / dealii.git / commitdiff