From 87df4cbcac81c6688339358a2d1cbcf071f99469 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 16 Mar 1998 08:35:35 +0000 Subject: [PATCH] Add in/out functionality. git-svn-id: https://svn.dealii.org/trunk@68 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/grid/tria.cc | 298 ++++++++++++++++++++- deal.II/deal.II/source/numerics/data_io.cc | 45 +++- 2 files changed, 328 insertions(+), 15 deletions(-) diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 45ad30ae2e..4841727281 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -1,12 +1,13 @@ /* $Id$ */ +#include #include #include #include #include #include - - +#include +#include template Triangulation::Triangulation () { @@ -128,9 +129,233 @@ void Triangulation<1>::create_triangulation (const vector > &v, void Triangulation<2>::create_triangulation (const vector > &v, - const vector > &cells) { - const unsigned int dim=1; - Assert (false, ExcInternalError()); + const vector > &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 (v.size(), true); + + // copy cells. This is needed since we + // may need to change entries + vector > 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,line_iterator,less > > needed_lines; + for (unsigned int cell=0; cell 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 vertex_touch_count (v.size(), 0); + map,line_iterator,less > >::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,line_iterator,less > >::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, less > adjacent_cells; + + // finally make up cells + if (true) + { + raw_cell_iterator cell = begin_raw_quad(); + for (unsigned int c=0; cset (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, less >::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 vertices[8] = { Point (a,a), + Point ((a+b)/2,a), + Point (b,a), + Point (a,(a+b)/2), + Point ((a+b)/2,(a+b)/2), + Point (b,(a+b)/2), + Point (a,b), + Point ((a+b)/2,b) }; + const int cell_vertices[3][4] = {{0, 1, 4, 3}, + {1, 2, 5, 4}, + {3, 4, 7, 6}}; + vector > cells (3, vector(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 >(&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 vertices[8] = { p+Point(-1,-1)*radius/sqrt(2), + p+Point(+1,-1)*radius/sqrt(2), + p+Point(-1,-1)*radius/sqrt(2)*a, + p+Point(+1,-1)*radius/sqrt(2)*a, + p+Point(-1,+1)*radius/sqrt(2)*a, + p+Point(+1,+1)*radius/sqrt(2)*a, + p+Point(-1,+1)*radius/sqrt(2), + p+Point(+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 > cells (5, vector(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 >(&vertices[0], &vertices[8]), + cells); +}; + + + template void Triangulation::set_all_refine_flags () { diff --git a/deal.II/deal.II/source/numerics/data_io.cc b/deal.II/deal.II/source/numerics/data_io.cc index edf2e2617c..d610124da0 100644 --- a/deal.II/deal.II/source/numerics/data_io.cc +++ b/deal.II/deal.II/source/numerics/data_io.cc @@ -112,7 +112,15 @@ void DataIn::read_ucd (istream &in) { // transform from ucd to // consecutive numbering for (unsigned int i=0; i<(1<::attach_dof_handler (DoFHandler &d) { template -void DataOut::add_data_vector (dVector &vec, - String &name, - String &units) { - DataEntry new_entry = {&vec, name, units}; +void DataOut::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::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::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; -- 2.39.5