]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Delete unused vertices. Debug output.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 1 Oct 2000 22:50:23 +0000 (22:50 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 1 Oct 2000 22:50:23 +0000 (22:50 +0000)
git-svn-id: https://svn.dealii.org/trunk@3364 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/grid_in.h
deal.II/deal.II/source/grid/grid_in.cc
deal.II/doc/news/2000/c-3-0.html

index 4088de7361e5ac4abc99d495fe7ff68a7eebe853..c83b75d2e9e83728f94e676ac6a34f96b5266ac5 100644 (file)
@@ -25,9 +25,13 @@ template <int dim> class Triangulation;
 /**
  * This class implements an input mechanism for grid data. It allows
  * to read a grid structure into a triangulation object. At present,
- * only UCD (unstructured cell data) is supported as input format for
- * grid data. Any numerical data after the block of topological
- * information is ignored.
+ * only UCD (unstructured cell data) and DB Mesh is supported as input
+ * format for grid data. Any numerical data after the block of
+ * topological information is ignored.
+ *
+ * Note: if you experience unexpected problems with the use of this
+ * class, be sure to read the documentation right until the end, and
+ * also read the documentation of the @ref{GridReordering} class.
  *
  * To read grid data, the triangulation to be fed with has to be empty.
  * When giving a file which does not contain the assumed information or
@@ -67,13 +71,14 @@ template <int dim> class Triangulation;
  * @end{itemize}
  *
  *
- * @sect3{Structure of input grid data}
+ * @sect3{Structure of input grid data. The GridReordering class}
  * 
- * It is your duty to use a correct numbering of vertices in the cell list,
- * i.e. for lines, you have to first give the vertex with the lower coordinate
- * value, then that with the higher coordinate value. For quadrilaterals in
- * two dimensions, the vertex indices in the @p{quad} list have to be such that
- * the vertices are numbered in counter-clockwise sense.
+ * It is your duty to use a correct numbering of vertices in the cell
+ * list, i.e. for lines in 1d, you have to first give the vertex with
+ * the lower coordinate value, then that with the higher coordinate
+ * value. For quadrilaterals in two dimensions, the vertex indices in
+ * the @p{quad} list have to be such that the vertices are numbered in
+ * counter-clockwise sense.
  *
  * In two dimensions, another difficulty occurs, which has to do with the sense
  * of a quadrilateral. A quad consists of four lines which have a direction,
@@ -92,60 +97,34 @@ template <int dim> class Triangulation;
  *   |   |   |
  *   0---1---2
  * @end{verbatim}
- * may be characterised by the vertex numbers (0 1 4 3) and (1 2 5 4), since
- * the middle line would get the direction @p{1->4} when viewed from both cells.
- * The numbering (0 1 4 3) and (5 4 1 2) would not be allowed, since the left
- * quad would give the common line the direction @p{1->4}, while the right one
- * would want to use @p{4->1}, leading to ambiguity. The @ref{Triangulation} object
- * is capable of detecting this special case, which can be eliminated by
- * rotating the indices of the right quad by two. However, it would not
- * know what to do if you gave the vertex indices (4 1 2 5), since then it
- * would have to rotate by one element or three, the decision which to take is
- * not yet implemented.
- *
- * There are more ambiguous cases, where the triangulation may not know what
- * to do at all without the use of very sophisticated algorithms. On such example
- * is the following:
- * @begin{verbatim}
- *   9---10-----11
- *   |   |    / |
- *   6---7---8  |
- *   |   |   |  |
- *   3---4---5  |
- *   |   |    \ |
- *   0---1------2
- * @end{verbatim}
- * Assume that you had numbered the vertices in the cells at the left boundary
- * in a way, that the following line directions are induced:
- * @begin{verbatim}
- *   9->-10-----11
- *   ^   ^    / |
- *   6->-7---8  |
- *   ^   ^   |  |
- *   3->-4---5  |
- *   ^   ^    \ |
- *   0->-1------2
- * @end{verbatim}
- * (This could for example be done by using the indices (0 1 4 3), (3 4 7 6),
- * (6 7 10 9) for the three cells). Now, you will not find a way of giving
- * indices for the right cells, without introducing either ambiguity for
- * one line or other, or without violating that within each cells, there must be
- * one vertex from which both lines are directed away and the opposite one to
- * which both adjacent lines point to.
+ * may be characterised by the vertex numbers @p{(0 1 4 3)} and
+ * @p{(1 2 5 4)}, since the middle line would get the direction @p{1->4}
+ * when viewed from both cells.  The numbering @p{(0 1 4 3)} and
+ * @p{(5 4 1 2)} would not be allowed, since the left quad would give the
+ * common line the direction @p{1->4}, while the right one would want
+ * to use @p{4->1}, leading to an ambiguity. The @ref{Triangulation}
+ * object is capable of detecting this special case, which can be
+ * eliminated by rotating the indices of the right quad by
+ * two. However, it would not know what to do if you gave the vertex
+ * indices @p{(4 1 2 5)}, since then it would have to rotate by one
+ * element or three, the decision which to take is not yet
+ * implemented.
  *
- * The solution in this case is to renumber one of the three left cells, e.g.
- * by reverting the sense of the line between vertices 7 and 10 by numbering
- * the top left cell by (9 6 7 10).
+ * There are more ambiguous cases, where the triangulation may not
+ * know what to do at all without the use of sophisticated
+ * algorithms. Furthermore, similar problems exist in three space
+ * dimensions, where faces and lines have orientations that need to be
+ * taken care of.
  *
- * But this is a thing that the triangulation
- * object can't do for you, since it would involve backtracking to cells
- * already created when we find that we can't number the indices of one of
- * the rightmost cells consistently. It is neither clear how to do this
- * backtracking nor whether it can be done with a stopping algorithm, if
- * possible within polynomial time. This kind of numbering must be made
- * upon construction of the coarse grid, unfortunately.
+ * For this reason, the @p{read_*} functions of this class that read
+ * in grids in various input formats call the @ref{GridReordering}
+ * class to bring the order of vertices that define the cells into an
+ * ordering that satisfies the requiremenets of the
+ * @ref{Triangulation} class. Be sure to read the documentation of
+ * that class if you experience unexpected problems when reading grids
+ * through this class.
  *
- * @author Wolfgang Bangerth, 1998
+ * @author Wolfgang Bangerth, 1998, 2000
  */
 template <int dim>
 class GridIn
@@ -249,6 +228,86 @@ class GridIn
                                      */
     static void skip_comment_lines (istream    &in,
                                    const char  comment_start);
+
+                                    /**
+                                     * Remove vertices that are not
+                                     * referenced by any of the
+                                     * cells. This function is called
+                                     * by all @p{read_*} functions to
+                                     * eliminate vertices that are
+                                     * listed in the input files but
+                                     * are not used by the cells in
+                                     * the input file. While these
+                                     * vertices should not be in the
+                                     * input from the beginning, they
+                                     * sometimes are, most often when
+                                     * some cells have been removed
+                                     * by hand without wanting to
+                                     * update the vertex lists, as
+                                     * they might be lengthy.
+                                     *
+                                     * This function is called by all
+                                     * @p{read_*} functions as the
+                                     * triangulation class requires
+                                     * them to be called with used
+                                     * vertices only. This is so,
+                                     * since the vertices are copied
+                                     * verbatim by that class, so we
+                                     * have to eliminate unused
+                                     * vertices beforehand.
+                                     */
+    static void delete_unused_vertices (vector<Point<dim> >    &vertices,
+                                       vector<CellData<dim> > &cells,
+                                       SubCellData            &subcelldata);
+
+                                    /**
+                                     * This function can write the
+                                     * raw cell data objects created
+                                     * by the @p{read_*} functions in
+                                     * Gnuplot format to a
+                                     * stream. This is sometimes
+                                     * handy if one would like to see
+                                     * what actually was created, if
+                                     * it is known that the data is
+                                     * not correct in some way, but
+                                     * the @red{Triangulation} class
+                                     * refuses to generate a
+                                     * triangulation because of these
+                                     * errors. In particular, the
+                                     * output of this class writes
+                                     * out the cell numbers along
+                                     * with the direction of the
+                                     * faces of each cell. In
+                                     * particular the latter
+                                     * information is needed to
+                                     * verify whether the cell data
+                                     * objects follow the
+                                     * requirements of the ordering
+                                     * of cells and their faces,
+                                     * i.e. that all faces need to
+                                     * have unique directions and
+                                     * specified orientations with
+                                     * respect to neighboring cells
+                                     * (see the documentations to
+                                     * this class and the
+                                     * @ref{GridReordering} class).
+                                     *
+                                     * The output of this function
+                                     * consists of vectors for each
+                                     * line bounding the cells
+                                     * indicating the direction it
+                                     * has with respect to the
+                                     * orientation of this cell, and
+                                     * the cell number. The whole
+                                     * output is in a form such that
+                                     * it can be read in by Gnuplot
+                                     * and generate the full plot
+                                     * without further ado by the
+                                     * user.
+                                     */
+    static void debug_output_grid (const vector<CellData<dim> > &cells,
+                                  const vector<Point<dim> >    &vertices,
+                                  ostream                      &out);
 };
 
 
index d38f6a1b250562ff247096ccbae3f04059f58486..513fa54729c47eb08f71057d88f919974305e109 100644 (file)
 #include <algorithm>
 
 
-template <int dim>
-void
-delete_unused_vertices (vector<Point<dim> >    &vertices,
-                       vector<CellData<dim> > &cells,
-                       SubCellData            &subcelldata)
-{
-                                  // first check which vertices are
-                                  // actually used
-  vector<bool> vertex_used (vertices.size(), false);
-  for (unsigned int c=0; c<cells.size(); ++c)
-    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-      vertex_used[cells[c].vertices[v]] = true;
-
-                                  // then renumber the vertices that
-                                  // are actually used in the same
-                                  // order as they were beforehand
-  const unsigned int invalid_vertex = static_cast<unsigned int>(-1);
-  vector<unsigned int> new_vertex_numbers (vertices.size(), invalid_vertex);
-  unsigned int next_free_number = 0;
-  for (unsigned int i=0; i<vertices.size(); ++i)
-    if (vertex_used[i] == true)
-      {
-       new_vertex_numbers[i] = next_free_number;
-       ++next_free_number;
-      };
-
-                                  // next replace old vertex numbers
-                                  // by the new ones
-  for (unsigned int c=0; c<cells.size(); ++c)
-    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
-      cells[c].vertices[v] = new_vertex_numbers[cells[c].vertices[v]];
-
-                                  // same for boundary data
-  for (unsigned int c=0; c<subcelldata.boundary_lines.size(); ++c)
-    for (unsigned int v=0; v<GeometryInfo<1>::vertices_per_cell; ++v)
-      subcelldata.boundary_lines[c].vertices[v]
-       = new_vertex_numbers[subcelldata.boundary_lines[c].vertices[v]];
-  for (unsigned int c=0; c<subcelldata.boundary_quads.size(); ++c)
-    for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
-      subcelldata.boundary_quads[c].vertices[v]
-       = new_vertex_numbers[subcelldata.boundary_quads[c].vertices[v]];
-
-                                  // finally copy over the vertices
-                                  // which we really need to a new
-                                  // array and replace the old one by
-                                  // the new one
-  vector<Point<dim> > tmp;
-  tmp.reserve (count(vertex_used.begin(), vertex_used.end(), true));
-  for (unsigned int v=0; v<vertices.size(); ++v)
-    if (vertex_used[v] == true)
-      tmp.push_back (vertices[v]);
-  swap (vertices, tmp);
-};
 
 
 
 #include <fstream>
-void
-debug_output (const vector<CellData<3> > &cells,
-             const vector<Point<3> >    &vertices,
-             const string               &filename);
-void
-debug_output (const vector<CellData<1> > &cells,
-             const vector<Point<1> >    &vertices,
-             const string               &filename);
-void
-debug_output (const vector<CellData<2> > &cells,
-             const vector<Point<2> >    &vertices,
-             const string               &filename)
-{
-  double min_x = vertices[cells[0].vertices[0]](0),
-        max_x = vertices[cells[0].vertices[0]](0),
-        min_y = vertices[cells[0].vertices[0]](1),
-        max_y = vertices[cells[0].vertices[0]](1);
-  
-  ofstream x(filename.c_str());
-  for (unsigned int i=0; i<cells.size(); ++i)
-    {
-      for (unsigned int v=0; v<4; ++v)
-       {
-         const Point<2> &p = vertices[cells[i].vertices[v]];
-         
-         if (p(0) < min_x)
-           min_x = p(0);
-         if (p(0) > max_x)
-           max_x = p(0);
-         if (p(1) < min_y)
-           min_y = p(1);
-         if (p(1) > max_y)
-           max_y = p(1);
-       };
-
-      x << "# cell " << i << endl;
-      Point<2> center;
-      for (unsigned int f=0; f<4; ++f)
-       center += vertices[cells[i].vertices[f]];
-      center /= 4;
-
-      x << "set label \"" << i << "\" at "
-       << center(0) << ',' << center(1)
-       << " center"
-       << endl;
-
-                                      // first two line right direction
-      for (unsigned int f=0; f<2; ++f)
-       x << "set arrow from "
-         << vertices[cells[i].vertices[f]](0) << ',' << vertices[cells[i].vertices[f]](1)
-         << " to "
-         << vertices[cells[i].vertices[(f+1)%4]](0) << ',' << vertices[cells[i].vertices[(f+1)%4]](1)
-         << endl;
-                                      // other two lines reverse direction
-      for (unsigned int f=2; f<4; ++f)
-       x << "set arrow from "
-         << vertices[cells[i].vertices[(f+1)%4]](0) << ',' << vertices[cells[i].vertices[(f+1)%4]](1)
-         << " to "
-         << vertices[cells[i].vertices[f]](0) << ',' << vertices[cells[i].vertices[f]](1)
-         << endl;
-      x << endl;
-    };
-  
-
-  x << endl
-    << "pl [" << min_x << ':' << max_x << "]["
-    << min_y << ':' << max_y <<  "] "
-    << min_y << endl
-    << "pause -1" << endl;
-};
-
-      
-
-
-
 
 
 template <int dim>
@@ -285,6 +157,8 @@ void GridIn<dim>::read_ucd (istream &in)
 
   AssertThrow (in, ExcIO());
 
+                                  // do some clean-up on vertices
+  delete_unused_vertices (vertices, cells, subcelldata);
   tria->create_triangulation (vertices, cells, subcelldata);
 };
 
@@ -447,10 +321,12 @@ void GridIn<dim>::read_dbmesh (istream &in)
 
                                   // do some clean-up on vertices
   delete_unused_vertices (vertices, cells, subcelldata);
-  
-  debug_output (cells, vertices, "x1");
+
+  ofstream x1("x1");
+  debug_output_grid (cells, vertices, x1);
   GridReordering<dim>::reorder_cells (cells);
-  debug_output (cells, vertices, "x2");
+  ofstream x2("x2");
+  debug_output_grid (cells, vertices, x2);
   tria->create_triangulation (vertices, cells, subcelldata);
 };
 
@@ -503,6 +379,143 @@ void GridIn<dim>::skip_comment_lines (istream    &in,
 
 
 
+template <int dim>
+void
+GridIn<dim>::delete_unused_vertices (vector<Point<dim> >    &vertices,
+                                    vector<CellData<dim> > &cells,
+                                    SubCellData            &subcelldata)
+{
+                                  // first check which vertices are
+                                  // actually used
+  vector<bool> vertex_used (vertices.size(), false);
+  for (unsigned int c=0; c<cells.size(); ++c)
+    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+      vertex_used[cells[c].vertices[v]] = true;
+
+                                  // then renumber the vertices that
+                                  // are actually used in the same
+                                  // order as they were beforehand
+  const unsigned int invalid_vertex = static_cast<unsigned int>(-1);
+  vector<unsigned int> new_vertex_numbers (vertices.size(), invalid_vertex);
+  unsigned int next_free_number = 0;
+  for (unsigned int i=0; i<vertices.size(); ++i)
+    if (vertex_used[i] == true)
+      {
+       new_vertex_numbers[i] = next_free_number;
+       ++next_free_number;
+      };
+
+                                  // next replace old vertex numbers
+                                  // by the new ones
+  for (unsigned int c=0; c<cells.size(); ++c)
+    for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+      cells[c].vertices[v] = new_vertex_numbers[cells[c].vertices[v]];
+
+                                  // same for boundary data
+  for (unsigned int c=0; c<subcelldata.boundary_lines.size(); ++c)
+    for (unsigned int v=0; v<GeometryInfo<1>::vertices_per_cell; ++v)
+      subcelldata.boundary_lines[c].vertices[v]
+       = new_vertex_numbers[subcelldata.boundary_lines[c].vertices[v]];
+  for (unsigned int c=0; c<subcelldata.boundary_quads.size(); ++c)
+    for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
+      subcelldata.boundary_quads[c].vertices[v]
+       = new_vertex_numbers[subcelldata.boundary_quads[c].vertices[v]];
+
+                                  // finally copy over the vertices
+                                  // which we really need to a new
+                                  // array and replace the old one by
+                                  // the new one
+  vector<Point<dim> > tmp;
+  tmp.reserve (count(vertex_used.begin(), vertex_used.end(), true));
+  for (unsigned int v=0; v<vertices.size(); ++v)
+    if (vertex_used[v] == true)
+      tmp.push_back (vertices[v]);
+  swap (vertices, tmp);
+};
+
+
+
+template <int dim>
+void GridIn<dim>::debug_output_grid (const vector<CellData<dim> > &/*cells*/,
+                                    const vector<Point<dim> >    &/*vertices*/,
+                                    ostream                      &/*out*/)
+{
+  Assert (false, ExcNotImplemented());
+};
+
+
+#if deal_II_dimension == 2
+
+template <>
+void
+GridIn<2>::debug_output_grid (const vector<CellData<2> > &cells,
+                             const vector<Point<2> >    &vertices,
+                             ostream                    &out)
+{
+  double min_x = vertices[cells[0].vertices[0]](0),
+        max_x = vertices[cells[0].vertices[0]](0),
+        min_y = vertices[cells[0].vertices[0]](1),
+        max_y = vertices[cells[0].vertices[0]](1);
+  
+  for (unsigned int i=0; i<cells.size(); ++i)
+    {
+      for (unsigned int v=0; v<4; ++v)
+       {
+         const Point<2> &p = vertices[cells[i].vertices[v]];
+         
+         if (p(0) < min_x)
+           min_x = p(0);
+         if (p(0) > max_x)
+           max_x = p(0);
+         if (p(1) < min_y)
+           min_y = p(1);
+         if (p(1) > max_y)
+           max_y = p(1);
+       };
+
+      out << "# cell " << i << endl;
+      Point<2> center;
+      for (unsigned int f=0; f<4; ++f)
+       center += vertices[cells[i].vertices[f]];
+      center /= 4;
+
+      out << "set label \"" << i << "\" at "
+         << center(0) << ',' << center(1)
+         << " center"
+         << endl;
+
+                                      // first two line right direction
+      for (unsigned int f=0; f<2; ++f)
+       out << "set arrow from "
+           << vertices[cells[i].vertices[f]](0) << ',' 
+           << vertices[cells[i].vertices[f]](1)
+           << " to "
+           << vertices[cells[i].vertices[(f+1)%4]](0) << ',' 
+           << vertices[cells[i].vertices[(f+1)%4]](1)
+           << endl;
+                                      // other two lines reverse direction
+      for (unsigned int f=2; f<4; ++f)
+       out << "set arrow from "
+           << vertices[cells[i].vertices[(f+1)%4]](0) << ',' 
+           << vertices[cells[i].vertices[(f+1)%4]](1)
+           << " to "
+           << vertices[cells[i].vertices[f]](0) << ',' 
+           << vertices[cells[i].vertices[f]](1)
+           << endl;
+      out << endl;
+    };
+  
+
+  out << endl
+      << "set nokey" << endl
+      << "pl [" << min_x << ':' << max_x << "]["
+      << min_y << ':' << max_y <<  "] "
+      << min_y << endl
+      << "pause -1" << endl;
+};
+
+#endif
+
 
 //explicit instantiations
 template class GridIn<deal_II_dimension>;
index 6cdb24b434d69bb8e156f931310ff329115efa6f..f50c45014db32017ee55118d71dc99d46bb43276 100644 (file)
@@ -469,6 +469,20 @@ documentation, etc</a>.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       Extended: The <code class="member">GridIn::delete_unused_vertices</code>
+       function now eliminates vertices from the input that are not
+       referenced by any of the cells in the input file. This makes is
+       simpler to delete some cells from the input file by hand,
+       without the need to update the vertex lists, which can be
+       tiring as several cells usually use each vertex. All functions
+       in the <code class="class">GridIn</code> reading grids in
+       several input files call this function before passing the data
+       to the triangulation object.
+       <br>
+       (WB 2000/09/26)
+       </p>       
+
   <li> <p>
        New: The <code class="class">GridIn</code>
        class can now read the basics of grids in DB Mesh format.

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.