From 534fd9ed41c7c9a594579a51d8caaa1ef9ec9b81 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 29 Jun 2009 18:31:15 +0000 Subject: [PATCH] Catch the case where we have distorted cells and throw an appropriate exception that can be caught. git-svn-id: https://svn.dealii.org/trunk@19002 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/grid_in.h | 18 ++++ deal.II/deal.II/include/grid/tria.h | 90 +++++++++++++++- deal.II/deal.II/source/grid/tria.cc | 72 +++++++++++++ deal.II/doc/news/changes.h | 13 ++- tests/bits/Makefile | 3 +- tests/bits/distorted_cells_01.cc | 102 +++++++++++++++++++ tests/bits/distorted_cells_01/cmp/generic | 13 +++ tests/bits/neighboring_cells_at_two_faces.cc | 9 +- tests/deal.II/grid_in.cc | 28 ++++- tests/deal.II/grid_in/cmp/generic | 2 + tests/deal.II/grid_in_02.cc | 15 ++- tests/deal.II/grid_in_02/cmp/generic | 1 + tests/deal.II/number_cache_02.cc | 10 +- tests/deal.II/number_cache_02/cmp/generic | 1 + 14 files changed, 363 insertions(+), 14 deletions(-) create mode 100644 tests/bits/distorted_cells_01.cc create mode 100644 tests/bits/distorted_cells_01/cmp/generic diff --git a/deal.II/deal.II/include/grid/grid_in.h b/deal.II/deal.II/include/grid/grid_in.h index b520fa3bd3..1ad1b5a97a 100644 --- a/deal.II/deal.II/include/grid/grid_in.h +++ b/deal.II/deal.II/include/grid/grid_in.h @@ -217,6 +217,24 @@ class SubCellData; * that class if you experience unexpected problems when reading grids * through this class. * + * + *

Dealing with distorted mesh cells

+ * + * For each of the mesh reading functions, the last call is always to + * Triangulation::create_triangulation(). That function checks whether + * all the cells it creates as part of the coarse mesh are distorted or + * not (where distortion here means that the Jacobian of the mapping + * from the reference cell to the real cell has a non-positive determinant, + * i.e. the cell is pinched or twisted). If it finds any such cells, it + * throws an exception. This exception is not caught in the grid reader + * functions of the current class, and so will propagate through to the + * function that called it. There, you can catch and ignore the + * exception if you are certain that there is no harm in dealing with + * such cells. If you were not aware that your mesh had such cells, your + * results will likely be of dubious quality at best if you ignore the + * exception. + * + * * @ingroup grid * @ingroup input * @author Wolfgang Bangerth, 1998, 2000, Luca Heltai, 2004, 2007 diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index e275f3db6f..f09ca48a3d 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -1360,6 +1360,44 @@ class Triangulation : public Subscriptor copy_notification (const Triangulation &old_tria, const Triangulation &new_tria); }; + + /** + * A structure that is used as an + * exception object by the + * create_triangulation() + * function to indicate which + * cells among the coarse mesh + * cells are inverted or severely + * distorted. + * + * Objects of this kind are + * thrown by the + * create_triangulation() + * function, and they can be + * caught in user code if this + * condition is to be ignored. + * + * A cell is called + * deformed if the + * determinant of the Jacobian of + * the mapping from reference + * cell to real cell is negative + * at least at one vertex. This + * computation is done using the + * GeometryInfo::jacobian_determinants_at_vertices + * function. + */ + struct DistortedCellList : public dealii::ExceptionBase + { + /** + * A list of those cells + * among the coarse mesh + * cells that are deformed. + */ + std::list::active_cell_iterator> + distorted_cells; + }; + /** * Make the dimension available @@ -1599,17 +1637,53 @@ class Triangulation : public Subscriptor * see the general class * documentation for this. * - * This function is made - * @p virtual to allow derived - * classes to set up some data - * structures as well. - * * For conditions when this * function can generate a valid * triangulation, see the * documentation of this class, * and the GridIn and * GridReordering class. + * + * At the very end of its + * operation, this function walks + * over all cells and verifies + * that none of the cells is + * deformed, where we call a cell + * deformed if the determinant of + * the Jacobian of the mapping + * from reference cell to real + * cell is negative at least at + * one of the vertices (this + * computation is done using the + * GeometryInfo::jacobian_determinants_at_vertices + * function). If there are + * deformed cells, this function + * throws an exception of kind + * DistortedCellList. Since this + * happens after all data + * structures have been set up, + * you can catch and ignore this + * exception if you know what you + * do -- for example, it may be + * that the determinant is zero + * (indicating that you have + * collapsed edges in a cell) but + * that this is ok because you + * didn't intend to integrate on + * this cell anyway. On the other + * hand, deformed cells are often + * a sign of a mesh that is too + * coarse to resolve the geometry + * of the domain, and in this + * case ignoring the exception is + * probably unwise. + * + * @note: The check for distorted + * cells is only done if + * dim==spacedim, as otherwise + * cells can legitimately be + * twisted if the manifold they + * describe is twisted. */ virtual void create_triangulation (const std::vector > &vertices, const std::vector > &cells, @@ -1624,6 +1698,12 @@ class Triangulation : public Subscriptor * new (lexicographic) ordering * and calls * create_triangulation(). + * + * @note This function internally + * calls create_triangulation and + * therefore can throw the same + * exception as the other + * function. */ virtual void create_triangulation_compatibility ( const std::vector > &vertices, diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 099f73c54e..f7704793a3 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -1132,6 +1132,63 @@ namespace } return numbers::invalid_unsigned_int; } + + + /** + * Collect all coarse mesh cells + * with at least one vertex at + * which the determinant of the + * Jacobian is zero or + * negative. This is the function + * for the case dim==spacedim. + */ + template + typename Triangulation::DistortedCellList + collect_distorted_coarse_cells (const Triangulation &triangulation) + { + typename Triangulation::DistortedCellList distorted_cells; + for (typename Triangulation::cell_iterator + cell = triangulation.begin(0); cell != triangulation.end(0); ++cell) + { + Point vertices[GeometryInfo::vertices_per_cell]; + for (unsigned int i=0; i::vertices_per_cell; ++i) + vertices[i] = cell->vertex(i); + + double determinants[GeometryInfo::vertices_per_cell]; + GeometryInfo::jacobian_determinants_at_vertices (vertices, + determinants); + + for (unsigned int i=0; i::vertices_per_cell; ++i) + if (determinants[i] <= 1e-9 * std::pow (cell->diameter(), + 1.*dim)) + { + distorted_cells.distorted_cells.push_back (cell); + break; + } + } + + return distorted_cells; + } + + + /** + * Collect all coarse mesh cells + * with at least one vertex at + * which the determinant of the + * Jacobian is zero or + * negative. This is the function + * for the case dim!=spacedim, + * where we can not determine + * whether a cell is twisted as it + * may, for example, discretize a + * manifold with a twist. + */ + template + typename Triangulation::DistortedCellList + collect_distorted_coarse_cells (const Triangulation &) + { + return typename Triangulation::DistortedCellList(); + } }// end of anonymous namespace @@ -9879,6 +9936,21 @@ create_triangulation (const std::vector > &v, } compute_number_cache (*this, number_cache); + + // now verify that there are indeed + // no distorted cells. as per the + // documentation of this class, we + // first collect all distorted + // cells and then throw an + // exception if there are any + DistortedCellList distorted_cells = collect_distorted_coarse_cells (*this); + // throw the array (and fill the + // various location fields) if + // there are distorted + // cells. otherwise, just fall off + // the end of the function + AssertThrow (distorted_cells.distorted_cells.size() == 0, + distorted_cells); } diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 06c8f932ab..eda28ef6ce 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -31,7 +31,18 @@ inconvenience this causes.
  1. -

    +

    + Changed: Previously, the Triangulation::create_triangulation function silently + accepted input meshes with inverted cells (i.e. cells with a zero or negative + determinant of the Jacobian of the mapping from the reference cell to the + real cell). This has been changed now: The function checks whether cells + are distorted or inverted, and may throw a list of cells for which this + is the case as an exception. If you know that this is harmless, for example + if you have cells with collapsed vertices in your mesh but you do not intend + to integrate on them, then you can catch and ignore this message. In all + other cases, the output of your computations are likely to be wrong anyway. +
    + (WB 2009/06/29)

diff --git a/tests/bits/Makefile b/tests/bits/Makefile index 88b61d211e..b8f63adaa3 100644 --- a/tests/bits/Makefile +++ b/tests/bits/Makefile @@ -100,7 +100,8 @@ tests_x = grid_generator_?? \ neighboring_cells_at_two_faces \ solution_transfer \ make_boundary_constraints_* \ - count_dofs_per_block* + count_dofs_per_block* \ + distorted_cells_* # tests for the hp branch: # fe_collection_* diff --git a/tests/bits/distorted_cells_01.cc b/tests/bits/distorted_cells_01.cc new file mode 100644 index 0000000000..cb53cd8f2a --- /dev/null +++ b/tests/bits/distorted_cells_01.cc @@ -0,0 +1,102 @@ +//---------------------------- distorted_cells_01.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003, 2004, 2005, 2009 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- distorted_cells_01.cc --------------------------- + + +// check that indeed Triangulation::create_triangulation throws an +// exception if we have distorted cells + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +// create a (i) pinched cell (where two vertices coincide), or (ii) +// twisted cell (where two vertices are swapped) +template +void check (const unsigned int testcase) +{ + std::vector > vertices; + for (unsigned int v=0; v::vertices_per_cell; ++v) + vertices.push_back (GeometryInfo::unit_cell_vertex(v)); + + switch (testcase) + { + case 1: + deallog << "Pinched cell in " << dim << "d" << std::endl; + vertices[0] = vertices[1]; + break; + case 2: + deallog << "Twisted cell in " << dim << "d" << std::endl; + std::swap (vertices[0], vertices[1]); + break; + default: + Assert (false, ExcNotImplemented()); + } + + + std::vector > cells; + { + CellData cell; + for (unsigned int j=0; j::vertices_per_cell; ++j) + cell.vertices[j] = j; + cells.push_back (cell); + } + + Triangulation coarse_grid; + + bool flag = false; + try + { + coarse_grid.create_triangulation (vertices, cells, SubCellData()); + } + catch (typename Triangulation::DistortedCellList &dcv) + { + flag = true; + + deallog << dcv.distorted_cells.size() << " distorted cells" + << std::endl; + Assert (dcv.distorted_cells.front() == coarse_grid.begin(0), + ExcInternalError()); + } + + Assert (flag == true, ExcInternalError()); +} + + +int main () +{ + std::ofstream logfile("distorted_cells_01/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + for (unsigned int testcase=1; testcase<=2; ++testcase) + { + check<1> (testcase); + check<2> (testcase); + check<3> (testcase); + } +} + + + diff --git a/tests/bits/distorted_cells_01/cmp/generic b/tests/bits/distorted_cells_01/cmp/generic new file mode 100644 index 0000000000..c9864839e9 --- /dev/null +++ b/tests/bits/distorted_cells_01/cmp/generic @@ -0,0 +1,13 @@ + +DEAL::Pinched cell in 1d +DEAL::1 distorted cells +DEAL::Pinched cell in 2d +DEAL::1 distorted cells +DEAL::Pinched cell in 3d +DEAL::1 distorted cells +DEAL::Twisted cell in 1d +DEAL::1 distorted cells +DEAL::Twisted cell in 2d +DEAL::1 distorted cells +DEAL::Twisted cell in 3d +DEAL::1 distorted cells diff --git a/tests/bits/neighboring_cells_at_two_faces.cc b/tests/bits/neighboring_cells_at_two_faces.cc index b64e5023af..bdb472e037 100644 --- a/tests/bits/neighboring_cells_at_two_faces.cc +++ b/tests/bits/neighboring_cells_at_two_faces.cc @@ -67,7 +67,14 @@ void create_grid (Triangulation<2> &tria) // generate a triangulation // out of this GridReordering<2>::reorder_cells (cells); - tria.create_triangulation_compatibility (vertices, cells, SubCellData()); + try + { + tria.create_triangulation_compatibility (vertices, cells, SubCellData()); + } + catch (Triangulation<2>::DistortedCellList &dcv) + { + // ignore the exception + } } diff --git a/tests/deal.II/grid_in.cc b/tests/deal.II/grid_in.cc index a4519ea768..316148b103 100644 --- a/tests/deal.II/grid_in.cc +++ b/tests/deal.II/grid_in.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2002, 2003, 2004, 2005, 2007, 2008 by the deal.II authors +// Copyright (C) 2002, 2003, 2004, 2005, 2007, 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -69,7 +69,19 @@ void test2 () GridIn gi; gi.attach_triangulation (tria); std::ifstream in ("grid_in/2d.xda"); - gi.read_xda (in); + try + { + gi.read_xda (in); + } + catch (typename Triangulation::DistortedCellList &dcv) + { + // ignore the exception that we + // get because the mesh has + // distorted cells + deallog << dcv.distorted_cells.size() << " cells are distorted." + << std::endl; + } + int hash = 0; int index = 0; @@ -102,7 +114,17 @@ void check_file (const std::string name, Triangulation tria; GridIn gi; gi.attach_triangulation (tria); - gi.read(name, format); + try + { + gi.read(name, format); + } + catch (typename Triangulation::DistortedCellList &dcv) + { + // ignore the exception + deallog << dcv.distorted_cells.size() << " cells are distorted." + << std::endl; + } + deallog << name << '\t' << tria.n_vertices() << '\t' << tria.n_cells() diff --git a/tests/deal.II/grid_in/cmp/generic b/tests/deal.II/grid_in/cmp/generic index 906ae5960c..86db31fde6 100644 --- a/tests/deal.II/grid_in/cmp/generic +++ b/tests/deal.II/grid_in/cmp/generic @@ -32,6 +32,8 @@ 8 0 quad 6 14 9 2 9 0 quad 16 10 3 13 10 0 quad 15 16 13 4 +DEAL::25 cells are distorted. DEAL::1266405810 DEAL::grid_in/2d 16 10 +DEAL::25 cells are distorted. DEAL::grid_in/2d 29490 29037 diff --git a/tests/deal.II/grid_in_02.cc b/tests/deal.II/grid_in_02.cc index 31c9530738..3c7531d03a 100644 --- a/tests/deal.II/grid_in_02.cc +++ b/tests/deal.II/grid_in_02.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2002, 2003, 2004, 2005, 2007, 2008 by the deal.II authors +// Copyright (C) 2002, 2003, 2004, 2005, 2007, 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -53,7 +53,18 @@ void test2 () GridIn gi; gi.attach_triangulation (tria); std::ifstream in ("grid_in_02/2d.xda"); - gi.read_xda (in); + try + { + gi.read_xda (in); + } + catch (typename Triangulation::DistortedCellList &dcv) + { + // ignore the exception that we + // get because the mesh has + // distorted cells + deallog << dcv.distorted_cells.size() << " cells are distorted." + << std::endl; + } Triangulation<2>::active_cell_iterator cell = tria.begin_active(), diff --git a/tests/deal.II/grid_in_02/cmp/generic b/tests/deal.II/grid_in_02/cmp/generic index 0fd8fc12f0..d75467cd09 100644 --- a/tests/deal.II/grid_in_02/cmp/generic +++ b/tests/deal.II/grid_in_02/cmp/generic @@ -1,2 +1,3 @@ +DEAL::24 cells are distorted. DEAL::OK diff --git a/tests/deal.II/number_cache_02.cc b/tests/deal.II/number_cache_02.cc index 2cfcaa1904..eac4f99f7a 100644 --- a/tests/deal.II/number_cache_02.cc +++ b/tests/deal.II/number_cache_02.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2008 by the deal.II authors +// Copyright (C) 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -87,6 +87,14 @@ void test (const char *filename) { gi.read_xda (in); } + catch (typename Triangulation::DistortedCellList &dcv) + { + // ignore the exception that we + // get because the mesh has + // distorted cells + deallog << dcv.distorted_cells.size() << " cells are distorted." + << std::endl; + } catch (std::exception &exc) { deallog << " caught exception:" << std::endl diff --git a/tests/deal.II/number_cache_02/cmp/generic b/tests/deal.II/number_cache_02/cmp/generic index 0ce9e3a50a..ecfd3af8c2 100644 --- a/tests/deal.II/number_cache_02/cmp/generic +++ b/tests/deal.II/number_cache_02/cmp/generic @@ -1,5 +1,6 @@ DEAL::Reading grid_in/2d.xda +DEAL::25 cells are distorted. DEAL:: 29037 DEAL:: 29037 DEAL:: 58529 -- 2.39.5