\end{verbatim}
+ {\bf Creating a triangulation}
+
+ There are several possibilities to create a triangulation:
+ \begin{itemize}
+ \item Hypercube triangulations: a hypercube triangulation is a
+ domain which is the tensor product of an interval $[a,b]$ in
+ the given number of spatial dimensions. If you want to create such
+ a domain, which is a common test case for model problems, call
+ #Triangulation<dim>::create_hypercube (a,b)#, which produces a
+ hypercube domain triangulated with exactly one element. You can
+ get tensor product meshes by successive refinement of this cell.
+
+ \item Reading in a triangulation: By using an object of the \Ref{#DataIn#}
+ class, you can read in fairly general triangulations. See there for
+ more information. The mentionned class uses the interface described
+ directly below to transfer the data into the triangulation.
+
+ \item Explicitely creating a triangulation: you can create a triangulation
+ by providing a list of vertices and a list of cells. Each such cell
+ consists of a vector storing the indices of the vertices of this cell
+ in the vertex list. To see how this works, you can take a look at the
+ #DataIn<dim>::read_*# functions. The appropriate function to be
+ called is #Triangulation<dim>::create_triangulation (2)#.
+
+ Creating the hierarchical information needed for this library from
+ cells storing only vertex information can be a quite complex task.
+ For example in 2d, we have to create lines between vertices (but only
+ once, though there are two cells which link these two vertices) and
+ we have to create neighborship information. Grids being read in
+ should therefore not be too large, reading refined grids would be
+ inefficient. Apart from the performance aspect, refined grids do not
+ lend too well to multigrid algorithms, since solving on the coarsest
+ level is expensive. It is wiser in any case to read in a grid as coarse
+ as possible and then do the needed refinement steps.
+
+ It is your duty to guarantee that cells have the correct orientation.
+ \end{itemize}
+
+
+
{\bf History of a triangulation}
It is possible to reconstruct a grid from its refinement history, which
* is called the last time.
*/
void set_boundary (const Boundary<dim> *boundary_object);
+
+ /**
+ * Create a triangulation from a list
+ * of vertices and a list of cells, each of
+ * the latter being a list of #1<<dim#
+ * vertex indices. The triangulation must
+ * be empty upon calling this function and
+ * the cell list should be useful (connected
+ * domain, etc.).
+ */
+ void create_triangulation (const vector<Point<dim> > &vertices,
+ const vector<vector<int> > &cells);
/**
* Initialize the triangulation with a
- * unit hypercube (unit line in 1D,
- * unit square in 2D, etc) consisting
- * of exactly one cell.
+ * hypercube (line in 1D, square in 2D, etc)
+ * consisting of exactly one cell. The
+ * hypercube volume is the tensor product
+ * of the intervall $[left,right]$ in the
+ * present number of dimensions, where
+ * the limits are given as arguments. They
+ * default to zero and unity, then producing
+ * the unit hypercube.
*
* The triangulation needs to be void
* upon calling this function.
*/
- void create_unit_hypercube ();
+ void create_hypercube (const double left = 0.,
+ const double right= 1.);
int, int,
<< "The present grid has " << arg1 << " active cells, "
<< "but the one in the file had " << arg2);
+ /**
+ * Exception
+ */
DeclException0 (ExcGridReadError);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcTriangulationNotEmpty);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcInternalError);
//@}
protected:
/**
#endif
-template <int dim>
-class DoFHandler;
+template <int dim> class Triangulation;
+template <int dim> class DoFHandler;
class dVector;
+template <class Key, class T, class Compare> class map;
+template <class T> struct less;
+
+
+/**
+ This class implements an input mechanism for grid data. It allows to
+ read a grid structure into a triangulation object. Future versions
+ will also allow to read data on this grid into vectors.
+
+ 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.
+
+ 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
+ which does not keep to the right format, the state of the triangulation
+ will be undefined afterwards. Upon input, only lines in one dimension
+ and quads in two dimensions are accepted. All other cell types (e.g. lines
+ or triangles in two dimensions, quads and hexes in 3d) are ignored. No
+ warning is issued. The vertex and cell numbering in the UCD file, which
+ need not be consecutively, is lost upon transfer to the triangulation
+ object, since this one needs consecutively numbered elements.
+ */
+template <int dim>
+class DataIn {
+ public:
+ /**
+ * Constructor.
+ */
+ DataIn ();
+
+ /**
+ * Attach this triangulation
+ * to be fed with the grid data.
+ */
+ void attach_triangulation (Triangulation<dim> *tria);
+
+ /**
+ * Read grid data from an ucd file.
+ * Numerical data is ignored.
+ */
+ void read_ucd (istream &);
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNotImplemented);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNoTriangulationSelected);
+
+ private:
+ /**
+ * Store address of the triangulation to
+ * be fed with the data read in.
+ */
+ Triangulation<dim> *tria;
+};
+
+
+
+
+
/**
-void Triangulation<1>::create_unit_hypercube () {
- vertices.push_back (Point<1> (0.0));
- vertices.push_back (Point<1> (1.0));
+
+void Triangulation<1>::create_triangulation (const vector<Point<1> > &v,
+ const vector<vector<int> > &cells) {
+ const unsigned int dim=1;
+
+ Assert (vertices.size() == 0, ExcTriangulationNotEmpty());
+ Assert (n_lines() == 0, ExcTriangulationNotEmpty());
+
+ // copy vertices
+ vertices = v;
+ vertices_used = vector<bool> (v.size(), true);
+
+ // store the indices of the lines which
+ // are adjacent to a given vertex
+ vector<vector<int> > lines_at_vertex (v.size());
+
+ // reserve enough space
+ levels[0]->TriangulationLevel<0>::reserve_space (cells.size(), dim);
+ levels[0]->TriangulationLevel<1>::reserve_space (cells.size());
+
+ // make up cells
+ raw_line_iterator next_free_line = begin_raw_line ();
+ for (unsigned int cell=0; cell<cells.size(); ++cell)
+ {
+ while (next_free_line->used())
+ ++next_free_line;
+
+ next_free_line->set (Line (cells[cell][0], cells[cell][1]));
+ next_free_line->set_used_flag ();
+
+ // note that this cell
+ // is adjacent to these vertices
+ lines_at_vertex[cells[cell][0]].push_back (cell);
+ lines_at_vertex[cells[cell][1]].push_back (cell);
+ };
+
+
+#ifdef DEBUG
+ // some security tests
+ unsigned int boundary_nodes = 0;
+ for (unsigned int i=0; i<lines_at_vertex.size(); ++i)
+ switch (lines_at_vertex[i].size())
+ {
+ case 1: // this vertex has only one adjacent line
+ ++boundary_nodes;
+ break;
+ case 2:
+ break;
+ default: // a node must have one or two adjacent lines
+ Assert (false, ExcInternalError());
+ };
+
+ // assert there are no more than two boundary
+ // nodes
+ Assert (boundary_nodes == 2, ExcInternalError());
+#endif
+
+
+ // update neighborship info
+ active_line_iterator line = begin_active_line ();
+ // for all lines
+ for (; line!=end(); ++line)
+ // for each of the two vertices
+ for (unsigned int vertex=0; vertex<(1<<dim); ++vertex)
+ // if first cell adjacent to this
+ // vertex is the present one, then
+ // the neighbor is the second adjacent
+ // cell and vice versa
+ if (lines_at_vertex[line->vertex_index(vertex)][0] == line->index())
+ if (lines_at_vertex[line->vertex_index(vertex)].size() == 2)
+ {
+ const cell_iterator neighbor ((Triangulation<1>*)this,
+ 0, // level
+ lines_at_vertex[line->vertex_index(vertex)][1]);
+ line->set_neighbor (vertex, neighbor);
+ }
+ else
+ // no second adjacent cell entered
+ // -> cell at boundary
+ line->set_neighbor (vertex, end());
+ else
+ // present line is not first adjacent
+ // one -> first adjacent one is neighbor
+ {
+ const cell_iterator neighbor ((Triangulation<1>*)this,
+ 0, // level
+ lines_at_vertex[line->vertex_index(vertex)][0]);
+ line->set_neighbor (vertex, neighbor);
+ };
+};
+
+
+
+void Triangulation<2>::create_triangulation (const vector<Point<2> > &v,
+ const vector<vector<int> > &cells) {
+ const unsigned int dim=1;
+ Assert (false, ExcInternalError());
+};
+
+
+
+void Triangulation<1>::create_hypercube (const double left,
+ const double right) {
+ Assert (vertices.size() == 0, ExcTriangulationNotEmpty());
+ Assert (n_lines() == 0, ExcTriangulationNotEmpty());
+
+ vertices.push_back (Point<1> (left));
+ vertices.push_back (Point<1> (right));
vertices_used.insert (vertices_used.end(), 2, true);
-void Triangulation<2>::create_unit_hypercube () {
+void Triangulation<2>::create_hypercube (const double left,
+ const double right) {
+ Assert (vertices.size() == 0, ExcTriangulationNotEmpty());
+ Assert (n_lines() == 0, ExcTriangulationNotEmpty());
+ Assert (n_quads() == 0, ExcTriangulationNotEmpty());
+
// create vertices
- vertices.push_back (Point<2> (0,0));
- vertices.push_back (Point<2> (1,0));
- vertices.push_back (Point<2> (1,1));
- vertices.push_back (Point<2> (0,1));
+ vertices.push_back (Point<2> (left,left));
+ vertices.push_back (Point<2> (right,left));
+ vertices.push_back (Point<2> (right,right));
+ vertices.push_back (Point<2> (left,right));
vertices_used.insert (vertices_used.end(), 4, true);
// create lines
#include <grid/dof_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/tria.h>
+#include <map.h>
#include <iostream.h>
#include <algo.h>
#include <time.h>
extern TriaActiveIterator<2,DoFCellAccessor<2> > __dummy567;
+
+
+
+
+template <int dim>
+DataIn<dim>::DataIn () :
+ tria(0) {};
+
+
+
+template <int dim>
+void DataIn<dim>::attach_triangulation (Triangulation<dim> *t) {
+ tria = t;
+};
+
+
+
+template <int dim>
+void DataIn<dim>::read_ucd (istream &in) {
+ Assert (tria != 0, ExcNoTriangulationSelected());
+ Assert ((1<=dim) && (dim<=2), ExcNotImplemented());
+
+
+ // skip comments at start of file
+ char c;
+ while (in.get(c), c=='#')
+ {
+ char line[256];
+ in.get (line, 255, '\n'); // ignore rest of line, at most 256 chars
+ in.get (c); // ignore '\n' at end of line.
+ };
+
+ // put back first character of
+ // first non-comment line
+ in.putback (c);
+
+
+ unsigned int n_vertices;
+ unsigned int n_cells;
+ int dummy;
+
+ in >> n_vertices
+ >> n_cells
+ >> dummy // number of data vectors
+ >> dummy // cell data
+ >> dummy; // model data
+
+ // set up array of vertices
+ vector<Point<dim> > vertices (n_vertices);
+ // set up mapping between numbering
+ // in ucd-file (key) and in the
+ // vertices vector
+ map<int,int,less<int> > vertex_indices;
+
+ for (unsigned int vertex=0; vertex<n_vertices; ++vertex)
+ {
+ int vertex_number;
+ double x[3];
+
+ // read vertex
+ in >> vertex_number
+ >> x[0] >> x[1] >> x[2];
+
+ // store vertex
+ for (unsigned int d=0; d<dim; ++d)
+ vertices[vertex](d) = x[d];
+ // store mapping; note that
+ // vertices_indices[i] is automatically
+ // created upon first usage
+ vertex_indices[vertex_number] = vertex;
+ };
+
+ // set up array of cells
+ vector<vector<int> > cells;
+
+ for (unsigned int cell=0; cell<n_cells; ++cell)
+ {
+ String cell_type;
+
+ in >> dummy // cell number
+ >> dummy; // material id
+ in >> cell_type;
+
+ if (((cell_type = "line") && (dim == 1)) ||
+ ((cell_type = "quad") && (dim == 2)))
+ // ignore lines in more than one
+ // dimension, quads in more than
+ // two, and triangles in any dimension
+ {
+ // allocate and read indices
+ cells.push_back (vector<int> (1<<dim));
+ for (unsigned int i=0; i<(1<<dim); ++i)
+ in >> cells.back()[i];
+
+ // transform from ucd to
+ // consecutive numbering
+ for (unsigned int i=0; i<(1<<dim); ++i)
+ cells.back()[i] = vertex_indices[cells.back()[i]];
+ };
+ };
+
+ tria->create_triangulation (vertices, cells);
+};
+
+
+
+
+
+
+
template <int dim>
DataOut<dim>::DataOut () :
dofs(0) {};
//explicite instantiations
+template class DataIn<1>;
+template class DataIn<2>;
template class DataOut<1>;
template class DataOut<2>;