]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add grid read functionality.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 12 Mar 1998 16:59:42 +0000 (16:59 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 12 Mar 1998 16:59:42 +0000 (16:59 +0000)
git-svn-id: https://svn.dealii.org/trunk@62 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/grid/tria.h
deal.II/deal.II/include/numerics/data_io.h
deal.II/deal.II/source/grid/tria.cc
deal.II/deal.II/source/numerics/data_io.cc

index c3385844d54aaae817324684812585df11ab7b87..684787ce82fcf5fe2fe258c67b80710c0c0d58b0 100644 (file)
@@ -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<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
@@ -803,17 +843,35 @@ class Triangulation : public TriaDimensionInfo<dim> {
                                      *  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.);
 
 
 
@@ -1366,7 +1424,18 @@ class Triangulation : public TriaDimensionInfo<dim> {
                    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:
                                     /**
index 7751111e5d537fc543d026770c8ec29ecd24dfcf..2789980bdcb8a932a7d67a10aa1a2bc811d58642 100644 (file)
 #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;
+};
+
+
+
+
+
 
 
 /**
index 1b4aa3cccffb2827f8b111a0ef3784fe8a1d90d3..45ad30ae2e3f1b60475cd3b69eb3adeced9ff5db 100644 (file)
@@ -34,9 +34,114 @@ void Triangulation<dim>::set_boundary (const Boundary<dim> *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<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);
   
@@ -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
index b7aa92856af32d420cce0ff7297ede3b8023fdb0..edf2e2617c066bea7fbdd6c16bf6d8f3211ac421 100644 (file)
@@ -5,6 +5,7 @@
 #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>
@@ -14,6 +15,116 @@ extern TriaActiveIterator<1,DoFCellAccessor<1> > __dummy566; // wait for gcc2.8
 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) {};
@@ -227,5 +338,7 @@ void DataOut<dim>::write_gnuplot (ostream &out) const {
 
 
 //explicite instantiations
+template class DataIn<1>;
+template class DataIn<2>;
 template class DataOut<1>;
 template class DataOut<2>;

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.