]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Catch the case where we have distorted cells and throw an appropriate exception that...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 29 Jun 2009 18:31:15 +0000 (18:31 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 29 Jun 2009 18:31:15 +0000 (18:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@19002 0785d39b-7218-0410-832d-ea1e28bc413d

14 files changed:
deal.II/deal.II/include/grid/grid_in.h
deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/source/grid/tria.cc
deal.II/doc/news/changes.h
tests/bits/Makefile
tests/bits/distorted_cells_01.cc [new file with mode: 0644]
tests/bits/distorted_cells_01/cmp/generic [new file with mode: 0644]
tests/bits/neighboring_cells_at_two_faces.cc
tests/deal.II/grid_in.cc
tests/deal.II/grid_in/cmp/generic
tests/deal.II/grid_in_02.cc
tests/deal.II/grid_in_02/cmp/generic
tests/deal.II/number_cache_02.cc
tests/deal.II/number_cache_02/cmp/generic

index b520fa3bd36a8d3ee7a9f6b7af9fbbc8f289755e..1ad1b5a97abc508977ae357d94d5f5ef6bab6af3 100644 (file)
@@ -217,6 +217,24 @@ class SubCellData;
  * that class if you experience unexpected problems when reading grids
  * through this class.
  *
+ *
+ * <h3>Dealing with distorted mesh cells</h3>
+ *
+ * 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
index e275f3db6f080058b869f7215cefdaf9336b073b..f09ca48a3d0c5c3fc083efbab7588e262241b384 100644 (file)
@@ -1360,6 +1360,44 @@ class Triangulation : public Subscriptor
         copy_notification (const Triangulation<dim, spacedim> &old_tria,
                           const Triangulation<dim, spacedim> &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
+                                     * <i>deformed</i> 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<typename Triangulation<dim,spacedim>::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<Point<spacedim> >    &vertices,
                                       const std::vector<CellData<dim> > &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<Point<spacedim> >    &vertices,
index 099f73c54ece20691910f5bbc6cdec06c86f6e75..f7704793a3d87d3e7c60edb7254898d1641a1013 100644 (file)
@@ -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 <int dim>
+  typename Triangulation<dim,dim>::DistortedCellList
+  collect_distorted_coarse_cells (const Triangulation<dim,dim> &triangulation)
+  {
+    typename Triangulation<dim,dim>::DistortedCellList distorted_cells;
+    for (typename Triangulation<dim,dim>::cell_iterator
+          cell = triangulation.begin(0); cell != triangulation.end(0); ++cell)
+      {
+       Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
+       for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+         vertices[i] = cell->vertex(i);
+
+       double determinants[GeometryInfo<dim>::vertices_per_cell];
+       GeometryInfo<dim>::jacobian_determinants_at_vertices (vertices,
+                                                             determinants);
+
+       for (unsigned int i=0; i<GeometryInfo<dim>::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 <int dim, int spacedim>
+  typename Triangulation<dim,spacedim>::DistortedCellList
+  collect_distorted_coarse_cells (const Triangulation<dim,spacedim> &)
+  {
+    return typename Triangulation<dim,spacedim>::DistortedCellList();
+  }
   
 }// end of anonymous namespace
 
@@ -9879,6 +9936,21 @@ create_triangulation (const std::vector<Point<spacedim> >    &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);
 }
 
 
index 06c8f932abdd97e1913aae9e084ed58dc2827882..eda28ef6ce75c84a253fa7b0a74dfad1a29b9d89 100644 (file)
@@ -31,7 +31,18 @@ inconvenience this causes.
 
 <ol>
   <li>
-  <p> 
+  <p>
+  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.
+  <br>
+  (WB 2009/06/29)
   </p>
   </li>
 </ol>
index 88b61d211e1a9e3d33a0705629cd2d08e7c0aa97..b8f63adaa304a280c5c75f303f1519789a341b4e 100644 (file)
@@ -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 (file)
index 0000000..cb53cd8
--- /dev/null
@@ -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 <base/logstream.h>
+#include <base/quadrature_lib.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_reordering.h>
+#include <grid/grid_generator.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+
+#include <fstream>
+
+
+// create a (i) pinched cell (where two vertices coincide), or (ii)
+// twisted cell (where two vertices are swapped)
+template <int dim>
+void check (const unsigned int testcase)
+{
+  std::vector<Point<dim> > vertices;
+  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+    vertices.push_back (GeometryInfo<dim>::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<CellData<dim> > cells;
+  {
+    CellData<dim> cell;
+    for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
+      cell.vertices[j]   = j;
+    cells.push_back (cell);
+  }
+
+  Triangulation<dim> coarse_grid;
+
+  bool flag = false;
+  try
+    {
+      coarse_grid.create_triangulation (vertices, cells, SubCellData());
+    }
+  catch (typename Triangulation<dim>::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 (file)
index 0000000..c986483
--- /dev/null
@@ -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
index b64e5023afd95db432ae2942fef8ee804c804cde..bdb472e0378819ecf9c164cbcfa208fcc45b5c9f 100644 (file)
@@ -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
+    }
 }
 
 
index a4519ea768bfc9765517abfcfc47bde7b919e1c0..316148b10317250e369e4552a36480f85b29ebfa 100644 (file)
@@ -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<dim> gi;
   gi.attach_triangulation (tria);
   std::ifstream in ("grid_in/2d.xda");
-  gi.read_xda (in);
+  try
+    {
+      gi.read_xda (in);
+    }
+  catch (typename Triangulation<dim>::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<dim> tria;
   GridIn<dim> gi;
   gi.attach_triangulation (tria);
-  gi.read(name, format);
+  try
+    {
+      gi.read(name, format);
+    }
+  catch (typename Triangulation<dim>::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()
index 906ae5960cabcae249c7069f9af3439ba93c9109..86db31fde6c12732f7838935d5ea7f31480cf4ee 100644 (file)
@@ -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
index 31c9530738bef4bddd11335b6b4aef16543e641e..3c7531d03abcb508f4b35dc91f50c28fbb42713a 100644 (file)
@@ -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<dim> 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<dim>::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(),
index 0fd8fc12f0b442283edd8868867114c242b04d11..d75467cd096fa38bc6991d1bcaeca74565285f9f 100644 (file)
@@ -1,2 +1,3 @@
 
+DEAL::24 cells are distorted.
 DEAL::OK
index 2cfcaa190404569486913e50ad9f1c7c25fde4bb..eac4f99f7a37eb65e9954279bee7e76c0b76399d 100644 (file)
@@ -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<dim>::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
index 0ce9e3a50ac24de1bff7a3a7c1fba0f83b853f12..ecfd3af8c20d4cac1d12a05144ded37f5f76e5eb 100644 (file)
@@ -1,5 +1,6 @@
 
 DEAL::Reading grid_in/2d.xda
+DEAL::25 cells are distorted.
 DEAL::  29037
 DEAL::  29037
 DEAL::  58529

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.