* 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
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
* 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,
* 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,
}
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
}
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);
}
<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>
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_*
--- /dev/null
+//---------------------------- 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);
+ }
+}
+
+
+
--- /dev/null
+
+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
// 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
+ }
}
// $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
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;
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()
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
// $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
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(),
+DEAL::24 cells are distorted.
DEAL::OK
// $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
{
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
DEAL::Reading grid_in/2d.xda
+DEAL::25 cells are distorted.
DEAL:: 29037
DEAL:: 29037
DEAL:: 58529