]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add in/out functionality.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 16 Mar 1998 08:35:35 +0000 (08:35 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 16 Mar 1998 08:35:35 +0000 (08:35 +0000)
git-svn-id: https://svn.dealii.org/trunk@68 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 45ad30ae2e3f1b60475cd3b69eb3adeced9ff5db..48417272815d2294ba5d97bfb5f1fbf9a5a7c79f 100644 (file)
@@ -1,12 +1,13 @@
 /* $Id$ */
 
+#include <grid/tria_boundary.h>
 #include <grid/tria_iterator.h>
 #include <grid/tria.h>
 #include <basic/magic_numbers.h>
 #include <iostream.h>
 #include <algo.h>
-
-
+#include <map.h>
+#include <cmath>
 
 template <int dim>
 Triangulation<dim>::Triangulation () {
@@ -128,9 +129,233 @@ void Triangulation<1>::create_triangulation (const vector<Point<1> >    &v,
 
 
 void Triangulation<2>::create_triangulation (const vector<Point<2> >    &v,
-                                            const vector<vector<int> > &cells) {
-  const unsigned int dim=1;
-  Assert (false, ExcInternalError());
+                                            const vector<vector<int> > &c) {
+  const unsigned int dim=2;
+
+  Assert (vertices.size() == 0, ExcTriangulationNotEmpty());
+  Assert (n_lines() == 0, ExcTriangulationNotEmpty());
+  Assert (n_quads() == 0, ExcTriangulationNotEmpty());
+
+                                  // copy vertices
+  vertices = v;
+  vertices_used = vector<bool> (v.size(), true);
+
+                                  // copy cells. This is needed since we
+                                  // may need to change entries
+  vector<vector<int> > cells(c);
+
+  
+                                  // make up a list of the needed lines
+                                  // each line is a pair of vertices. The list
+                                  // is kept sorted and it is guaranteed
+                                  // that each line is inserted only once.
+                                  // While the key of such an entry is the
+                                  // pair of vertices, the thing it points
+                                  // to is an iterator pointing to the line
+                                  // object itself. In the first run, these
+                                  // iterators are all invalid ones, but they
+                                  // are filled afterwards
+  map<pair<int,int>,line_iterator,less<pair<int,int> > > needed_lines;
+  for (unsigned int cell=0; cell<cells.size(); ++cell)
+    {
+#ifdef DEBUG
+                                      // in debug mode: check whether vertex
+                                      // indices are valid ones
+      for (unsigned int vertex=0; vertex<4; ++vertex)
+       Assert ((0<=cells[cell][vertex]) &&
+               (cells[cell][vertex]<(signed int)vertices.size()),
+               ExcInvalidVertexIndex (cell, cells[cell][vertex], vertices.size()));
+#endif
+      
+      pair<int,int> line_vertices[4] = {   // note the order of the vertices
+           make_pair (cells[cell][0], cells[cell][1]),
+           make_pair (cells[cell][1], cells[cell][2]),
+           make_pair (cells[cell][0], cells[cell][3]),
+           make_pair (cells[cell][3], cells[cell][2])};
+
+                                      // note the following: if the sense
+                                      // of the vertices of a cell is correct,
+                                      // but the vertices are given in an
+                                      // order which makes the sense of one line
+                                      // ambiguous when viewed from the two
+                                      // adjacent cells, we can heal this by
+                                      // shifting the vertex indices of one
+                                      // cell by two (diagonally exchanging
+                                      // the two vertices from which the
+                                      // four lines originate and to which
+                                      // they converge).
+                                      // If two lines are wrong, we could heal
+                                      // this by rotating by one or three
+                                      // vertices, but deciding this is
+                                      // difficult and not implemented.
+      for (unsigned int line=0; line<4; ++line)
+       if (needed_lines.find(make_pair(line_vertices[line].second,
+                                             line_vertices[line].first))
+           !=
+           needed_lines.end())
+         {
+                                            // rotate vertex numbers
+           rotate (cells[cell].begin(), cells[cell].begin()+2, cells[cell].end());
+                                            // remake lines
+           line_vertices[0] = make_pair (cells[cell][0], cells[cell][1]);
+           line_vertices[1] = make_pair (cells[cell][1], cells[cell][2]);
+           line_vertices[2] = make_pair (cells[cell][0], cells[cell][3]);
+           line_vertices[3] = make_pair (cells[cell][3], cells[cell][2]);
+                                            // allow for only one such
+                                            // rotation
+           break;
+         };
+      
+      
+      for (unsigned int line=0; line<4; ++line)
+       {
+                                          // assert that the line was not
+                                          // already inserted in reverse
+                                          // order. This happens in spite of
+                                          // the vertex rotation above, if the
+                                          // sense of the cell was incorrect.
+                                          //
+                                          // Here is what usually happened when
+                                          // this exception is thrown:
+                                          // consider these two cells
+                                          // and the vertices
+                                          //  3---4---5
+                                          //  |   |   |
+                                          //  0---1---2
+                                          // If in the input vector the
+                                          // two cells are given with
+                                          // vertices <0 1 4 3> and
+                                          // <4 1 2 5>, in the first cell
+                                          // the middle line would have
+                                          // direction 1->4, while in
+                                          // the second it would be 4->1.
+                                          // This will cause the exception.
+         Assert (needed_lines.find(make_pair(line_vertices[line].second,
+                                             line_vertices[line].first))
+                 ==
+                 needed_lines.end(),
+                 ExcGridHasInvalidCell(cell));
+                 
+                                          // insert line, with invalid iterator
+                                          // if line already exists, then
+                                          // nothing bad happens here
+         needed_lines[line_vertices[line]] = end_line();
+       };
+    };
+  
+
+#ifdef DEBUG
+                                  // in debug mode: check the every vertex has
+                                  // at least two adjacent lines
+  if (true) 
+    {
+      vector<unsigned short int> vertex_touch_count (v.size(), 0);
+      map<pair<int,int>,line_iterator,less<pair<int,int> > >::iterator i;
+      for (i=needed_lines.begin(); i!=needed_lines.end(); i++) 
+       {
+                                          // touch the vertices of this line
+         ++vertex_touch_count[(*i).first.first];
+         ++vertex_touch_count[(*i).first.second];
+       };
+
+                                      // assert minimum touch count is at
+                                      // least two
+      Assert (* (min_element(vertex_touch_count.begin(),
+                            vertex_touch_count.end())) >= 2,
+             ExcGridHasInvalidVertices());
+    };
+#endif
+       
+                                  // reserve enough space
+  levels[0]->TriangulationLevel<0>::reserve_space (cells.size(), dim);
+  levels[0]->TriangulationLevel<1>::reserve_space (needed_lines.size());
+  levels[0]->TriangulationLevel<2>::reserve_space (cells.size());
+
+                                  // make up lines
+  if (true) 
+    {
+      raw_line_iterator line = begin_raw_line();
+      map<pair<int,int>,line_iterator,less<pair<int,int> > >::iterator i;
+      for (i = needed_lines.begin(); line!=end_line(); ++line, ++i) 
+       {
+         line->set (Line((*i).first.first, (*i).first.second));
+         line->set_used_flag ();
+         (*i).second = line;
+       };
+    };
+
+
+                                  // store for each line index
+                                  // the adjacent cells
+  map<int,vector<cell_iterator>, less<int> > adjacent_cells;
+
+                                  // finally make up cells
+  if (true) 
+    {
+      raw_cell_iterator cell = begin_raw_quad();
+      for (unsigned int c=0; c<cells.size(); ++c, ++cell)
+       {
+                                          // list of iterators of lines
+         const line_iterator lines[4] = {
+               needed_lines[make_pair(cells[c][0], cells[c][1])],
+               needed_lines[make_pair(cells[c][1], cells[c][2])],
+               needed_lines[make_pair(cells[c][3], cells[c][2])],
+               needed_lines[make_pair(cells[c][0], cells[c][3])]};
+         
+         cell->set (Quad(lines[0]->index(),
+                         lines[1]->index(),
+                         lines[2]->index(),
+                         lines[3]->index()));
+         
+         cell->set_used_flag ();
+         
+                                          // note that this cell is adjacent
+                                          // to the four lines
+         for (unsigned int line=0; line<4; ++line)
+           adjacent_cells[lines[line]->index()].push_back (cell);
+         
+                                          // make some checks on the vertices
+                                          // and their ordering
+         Assert (lines[0]->vertex_index(0) == lines[3]->vertex_index(0),
+                 ExcInternalErrorOnCell(c));
+         Assert (lines[0]->vertex_index(1) == lines[1]->vertex_index(0),
+                 ExcInternalErrorOnCell(c));
+         Assert (lines[1]->vertex_index(1) == lines[2]->vertex_index(1),
+                 ExcInternalErrorOnCell(c));
+         Assert (lines[2]->vertex_index(0) == lines[3]->vertex_index(1),
+                 ExcInternalErrorOnCell(c));
+       };
+    };
+  
+
+#ifdef DEBUG
+  map<int,vector<cell_iterator>, less<int> >::const_iterator adj;
+  for (adj=adjacent_cells.begin(); adj!=adjacent_cells.end(); ++adj)
+                                    // assert that every line has one or
+                                    // two adjacent cells
+    Assert (((*adj).second.size() >= 1) &&
+           ((*adj).second.size() <= 2),
+           ExcInternalError());
+#endif
+
+
+                                  // update neighborship info
+  for (cell_iterator cell=begin(); cell!=end(); ++cell)
+    for (unsigned int side=0; side<4; ++side)
+      if (adjacent_cells[cell->line(side)->index()][0] == cell)
+                                        // first adjacent cell is this one
+       {
+         if (adjacent_cells[cell->line(side)->index()].size() == 2)
+                                            // there is another adjacent cell
+           cell->set_neighbor (side,
+                               adjacent_cells[cell->line(side)->index()][1]);
+       }
+                                  // first adjacent cell is not this one,
+                                  // -> it must be the neighbor we are
+                                  // looking for
+      else
+       cell->set_neighbor (side,
+                           adjacent_cells[cell->line(side)->index()][0]);
 };
 
 
@@ -208,6 +433,69 @@ void Triangulation<2>::create_hypercube (const double left,
 
 
 
+void Triangulation<1>::create_hyper_L (const double, const double) {
+  Assert (false, ExcInternalError());
+};
+
+
+
+void Triangulation<2>::create_hyper_L (const double a, const double b) {
+  const unsigned int dim=2;
+  const Point<dim> vertices[8] = { Point<dim> (a,a),
+                                  Point<dim> ((a+b)/2,a),
+                                  Point<dim> (b,a),
+                                  Point<dim> (a,(a+b)/2),
+                                  Point<dim> ((a+b)/2,(a+b)/2),
+                                  Point<dim> (b,(a+b)/2),
+                                  Point<dim> (a,b),
+                                  Point<dim> ((a+b)/2,b)  };
+  const int cell_vertices[3][4] = {{0, 1, 4, 3},
+                                  {1, 2, 5, 4},
+                                  {3, 4, 7, 6}};
+  vector<vector<int> > cells (3, vector<int>(4,-1));
+  for (unsigned int i=0; i<3; ++i)
+    for (unsigned int j=0; j<4; ++j)
+      cells[i][j] = cell_vertices[i][j];
+  
+  create_triangulation (vector<Point<dim> >(&vertices[0], &vertices[8]),
+                       cells);
+};
+
+
+
+void Triangulation<1>::create_hyper_ball (const Point<1>, const double) {
+  Assert (false, ExcInternalError());
+};
+
+
+
+void Triangulation<2>::create_hyper_ball (const Point<2> p, const double radius) {
+  const unsigned int dim=2;
+  const double a = 0.453;         // equilibrate cell sizes
+  const Point<dim> vertices[8] = { p+Point<dim>(-1,-1)*radius/sqrt(2),
+                                  p+Point<dim>(+1,-1)*radius/sqrt(2),
+                                  p+Point<dim>(-1,-1)*radius/sqrt(2)*a,
+                                  p+Point<dim>(+1,-1)*radius/sqrt(2)*a,
+                                  p+Point<dim>(-1,+1)*radius/sqrt(2)*a,
+                                  p+Point<dim>(+1,+1)*radius/sqrt(2)*a,
+                                  p+Point<dim>(-1,+1)*radius/sqrt(2),
+                                  p+Point<dim>(+1,+1)*radius/sqrt(2) };
+  const int cell_vertices[5][4] = {{0, 1, 3, 2},
+                                  {0, 2, 4, 6},
+                                  {2, 3, 5, 4},
+                                  {1, 7, 5, 3},
+                                  {6, 4, 5, 7}};
+  vector<vector<int> > cells (5, vector<int>(4,-1));
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<4; ++j)
+      cells[i][j] = cell_vertices[i][j];
+  
+  create_triangulation (vector<Point<dim> >(&vertices[0], &vertices[8]),
+                       cells);
+};
+
+
+
 
 template <int dim>
 void Triangulation<dim>::set_all_refine_flags () {
index edf2e2617c066bea7fbdd6c16bf6d8f3211ac421..d610124da038b35e24674115fd106b8113cfe1e4 100644 (file)
@@ -112,7 +112,15 @@ void DataIn<dim>::read_ucd (istream &in) {
                                           // transform from ucd to
                                           // consecutive numbering
          for (unsigned int i=0; i<(1<<dim); ++i)
-           cells.back()[i] = vertex_indices[cells.back()[i]];
+           if (vertex_indices.find (cells.back()[i]) != vertex_indices.end())
+                                              // vertex with this index exists
+             cells.back()[i] = vertex_indices[cells.back()[i]];
+           else 
+             {
+                                                // no such vertex index
+               Assert (false, ExcInvalidVertexIndex(cell, cells.back()[i]));
+               cells.back()[i] = -1;
+             };
        };
     };
 
@@ -139,10 +147,10 @@ void DataOut<dim>::attach_dof_handler (DoFHandler<dim> &d) {
 
 
 template <int dim>
-void DataOut<dim>::add_data_vector (dVector &vec,
-                                   String  &name,
-                                   String  &units) {
-  DataEntry new_entry = {&vec, name, units};
+void DataOut<dim>::add_data_vector (const dVector &vec,
+                                   const String  &name,
+                                   const String  &units) {
+  DataEntry new_entry (&vec, name, units);
   data.push_back (new_entry);
 };
 
@@ -313,7 +321,8 @@ void DataOut<dim>::write_gnuplot (ostream &out) const {
                    out << cell->vertex(vertex)
                        << "   ";
                    for (unsigned int i=0; i!=data.size(); ++i)
-                     out << (*data[i].data)(cell->vertex_dof_index(vertex,0)) << ' ';
+                     out << (*data[i].data)(cell->vertex_dof_index(vertex,0))
+                         << ' ';
                    out << endl;
                  };
                
@@ -321,10 +330,26 @@ void DataOut<dim>::write_gnuplot (ostream &out) const {
 
          case 2:
                                                 // two dimension: output
-                                                // grid
-               for (unsigned int vertex=0; vertex<4; ++vertex)
-                 out << cell->vertex(vertex) << endl;
-               out << cell->vertex(0) << endl << endl;
+                                                // grid and data as a sequence
+                                                // of lines in 3d
+               for (unsigned int vertex=0; vertex<4; ++vertex) 
+                 {
+                   out << cell->vertex(vertex) << "   ";
+                   for (unsigned int i=0; i!=data.size(); ++i)
+                     out << (*data[i].data)(cell->vertex_dof_index(vertex,0))
+                         << ' ';
+                   out << endl;
+                 };
+                                                // first vertex again
+               out << cell->vertex(0) << "   ";
+               for (unsigned int i=0; i!=data.size(); ++i)
+                 out << (*data[i].data)(cell->vertex_dof_index(0,0))
+                     << ' ';
+               out << endl
+                   << endl
+                   << endl;      // end of cell; two newlines, since this
+                                                // stops continuous drawing
+                                                // of lines
 
                break;
 

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.