From: Wolfgang Bangerth Date: Thu, 12 Mar 1998 16:59:42 +0000 (+0000) Subject: Add grid read functionality. X-Git-Tag: v8.0.0~23189 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c4cdfd17bf11ce76a000492459b2186257f13ad0;p=dealii.git Add grid read functionality. git-svn-id: https://svn.dealii.org/trunk@62 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index c3385844d5..684787ce82 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -627,6 +627,46 @@ class TriaDimensionInfo<2> { \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::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::read_*# functions. The appropriate function to be + called is #Triangulation::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 @@ -803,17 +843,35 @@ class Triangulation : public TriaDimensionInfo { * is called the last time. */ void set_boundary (const Boundary *boundary_object); + + /** + * Create a triangulation from a list + * of vertices and a list of cells, each of + * the latter being a list of #1< > &vertices, + const vector > &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.); @@ -1366,7 +1424,18 @@ class Triangulation : public TriaDimensionInfo { 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: /** diff --git a/deal.II/deal.II/include/numerics/data_io.h b/deal.II/deal.II/include/numerics/data_io.h index 7751111e5d..2789980bdc 100644 --- a/deal.II/deal.II/include/numerics/data_io.h +++ b/deal.II/deal.II/include/numerics/data_io.h @@ -16,11 +16,75 @@ #endif -template -class DoFHandler; +template class Triangulation; +template class DoFHandler; class dVector; +template class map; +template 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 +class DataIn { + public: + /** + * Constructor. + */ + DataIn (); + + /** + * Attach this triangulation + * to be fed with the grid data. + */ + void attach_triangulation (Triangulation *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 *tria; +}; + + + + + /** diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 1b4aa3cccf..45ad30ae2e 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -34,9 +34,114 @@ void Triangulation::set_boundary (const Boundary *boundary_object) { -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 > &v, + const vector > &cells) { + const unsigned int dim=1; + + Assert (vertices.size() == 0, ExcTriangulationNotEmpty()); + Assert (n_lines() == 0, ExcTriangulationNotEmpty()); + + // copy vertices + vertices = v; + vertices_used = vector (v.size(), true); + + // store the indices of the lines which + // are adjacent to a given vertex + vector > 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; cellused()) + ++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; ivertex_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 > &v, + const vector > &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); @@ -53,12 +158,17 @@ void Triangulation<1>::create_unit_hypercube () { -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 diff --git a/deal.II/deal.II/source/numerics/data_io.cc b/deal.II/deal.II/source/numerics/data_io.cc index b7aa92856a..edf2e2617c 100644 --- a/deal.II/deal.II/source/numerics/data_io.cc +++ b/deal.II/deal.II/source/numerics/data_io.cc @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -14,6 +15,116 @@ extern TriaActiveIterator<1,DoFCellAccessor<1> > __dummy566; // wait for gcc2.8 extern TriaActiveIterator<2,DoFCellAccessor<2> > __dummy567; + + + + +template +DataIn::DataIn () : + tria(0) {}; + + + +template +void DataIn::attach_triangulation (Triangulation *t) { + tria = t; +}; + + + +template +void DataIn::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 > vertices (n_vertices); + // set up mapping between numbering + // in ucd-file (key) and in the + // vertices vector + map > vertex_indices; + + for (unsigned int vertex=0; vertex> vertex_number + >> x[0] >> x[1] >> x[2]; + + // store vertex + for (unsigned int d=0; d > cells; + + for (unsigned int cell=0; cell> 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 (1<> cells.back()[i]; + + // transform from ucd to + // consecutive numbering + for (unsigned int i=0; i<(1<create_triangulation (vertices, cells); +}; + + + + + + + template DataOut::DataOut () : dofs(0) {}; @@ -227,5 +338,7 @@ void DataOut::write_gnuplot (ostream &out) const { //explicite instantiations +template class DataIn<1>; +template class DataIn<2>; template class DataOut<1>; template class DataOut<2>;