/**
* 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
* @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,
* | | |
* 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
*/
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);
};
#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>
AssertThrow (in, ExcIO());
+ // do some clean-up on vertices
+ delete_unused_vertices (vertices, cells, subcelldata);
tria->create_triangulation (vertices, cells, subcelldata);
};
// 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);
};
+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>;
<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.
projects / dealii-svn.git / commitdiff