From: Wolfgang Bangerth Date: Tue, 24 Mar 1998 17:11:10 +0000 (+0000) Subject: Add support for material and boundary indicators. X-Git-Tag: v8.0.0~23157 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5d3e59b311733843becb45157dd363b91d698fbb;p=dealii.git Add support for material and boundary indicators. git-svn-id: https://svn.dealii.org/trunk@94 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 5bd827fa3a..a7fd173119 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -33,6 +33,9 @@ template class TriaActiveIterator; template class DoFHandler; +template struct CellData; +struct SubCellData; + class istream; class ostream; @@ -57,7 +60,16 @@ class TriangulationLevel; In #TriangulationLevel<0># all data is stored which is not dependant on the dimension, e.g. a field to store the refinement flag for the cells (what a cell actually is - is declared elsewhere), etc. + is declared elsewhere), etc. Actually, it is only cell-based + data, like neighborship info or refinement flags. There is another + field, which may fit in here, namely the material data (for cells) + or the boundary indicators (for faces), but since we need for a line + or quad either boundary information or material data, we store them + with the lines and quads rather than with the common data. We may, + however lose some memory in three dimensions, when we need the + material data for cell, boundary data for the quads, but nothing + for the lines. Since we only store one byte per line, quad or hex, + this is a minor loss and we can live with that. @memo Information belonging to one level of the multilevel hierarchy. */ @@ -125,7 +137,7 @@ class TriangulationLevel<0> { * to the mother cell of this cell). */ vector > neighbors; - + /** * Reserve enough space to accomodate * #total_cells# cells on this level. @@ -251,6 +263,29 @@ class TriangulationLevel<1> : public TriangulationLevel<0> { * #Triangulation<>::clear_user_flags()#. */ vector user_flags; + + /** + * Store boundary and material data. In + * one dimension, this field stores + * the material id of a line, which is + * number between 0 and 254. In more + * than one dimension, lines have no + * material id, but they may be at the + * boundary; then, we store the + * boundary indicator in this field, + * which denotes to which part of the + * boundary this line belongs and which + * boundary conditions hold on this + * part. The boundary indicator also + * is a number between zero and 254; + * the id 255 is reserved for lines + * in the interior and may be used + * to check whether a line is at the + * boundary or not, which otherwise + * is not possible if you don't know + * which cell it belongs to. + */ + vector material_id; }; public: @@ -337,6 +372,29 @@ class TriangulationLevel<2> : public TriangulationLevel<1> * Same as for #LineData::used#. */ vector user_flags; + + /** + * Store boundary and material data. In + * two dimension, this field stores + * the material id of a quad, which is + * number between 0 and 254. In more + * than two dimensions, quads have no + * material id, but they may be at the + * boundary; then, we store the + * boundary indicator in this field, + * which denotes to which part of the + * boundary this line belongs and which + * boundary conditions hold on this + * part. The boundary indicator also + * is a number between zero and 254; + * the id 255 is reserved for quads + * in the interior and may be used + * to check whether a quad is at the + * boundary or not, which otherwise + * is not possible if you don't know + * which cell it belongs to. + */ + vector material_id; }; public: @@ -704,6 +762,76 @@ class TriaDimensionInfo<2> { #Triangulation::create_triangulation (2)# function. \end{itemize} + The material id for each cell must be specified upon construction of + a triangulation. (There is a special section on material ids and + boundary indicators. See there for more information.) + The standard region functions (for hypercube, hyperball, + etc.) denote all cells the material id zero. You may change that afterwards, + but you should not use the material id 255. When reading a triangulation, + the material id must be specified in the input file (UCD format) or is + otherwise set to zero. When creating explicitely, the material id must + be given to the creation function. + + Regarding the boundary indicator for lines in two dimensions and quads + in three (subsummed by the word "faces"), all interior faces are denoted + the value 255. Trying to give an interior face another value results in + an error if in debug mode. Faces at the boundary of the domain are preset + with the boundary indicator zero, but you can give a list of faces with + different boundary indicators to the triangulation creation function. + The standard domain functions assume all faces to have boundary indicator + zero, which you may change manually afterwards. When reading from a file, + you have to give boundary indicators other than zero explicitely, e.g. in + UCD format by giving a list of lines with material id in the input file. + + Lines in two dimensions and quads in three dimensions inherit their + boundary indicator to their children upon refinement. You should therefore + make sure that if you have different boundary parts, the different parts + are separated by a vertex (in 2D) or a line (in 3D) such that each boundary + line or quad has a unique boundary indicator. + + Likewise, material data is inherited from mother to child cells. Place your + coarse level cells so, that the interface between cells is also the + interface between regions of different materials. + + + + {\bf Material and boundary information} + + Each line, quad, etc stores one byte of information denoting the material + a cell is made of (used in the computation of stiffness matrices) and to + which part of the boundary it belongs. Obviously, the material id is what + is needed for a cell, while for all structures with a dimension less than + the dimension of the domain (i.e. lines in 2D, lines and quads in 3D), the + boundary information is what is needed. Since either material or boundary + information is needed, but never both, only one field is used to store this + data, namely the #TriangulationLevel<1>::LinesData.material_id# and + #TriangulationLevel<2>::QuadsData.material_id# vectors. The data can be + read and written using line, quad and cell iterators. + + Material and boundary indicators are stored as one byte (an + #unsigned char#). They can therefore be in the range zero to 255, but + only zero to 254 is allowed. The value 255 is reserved to denote + interior lines (in 2D) and interior lines and quads (in 3D), which need + not have a boundary or material indicator. However, using this value, it + is possible to say whether a line in 2D is interior or not, which would + otherwise be impossible because the hierarchical structure of a + triangulation stores neighborship information and the like only with + cells. Finding out whether a line is an interior one would then only be + possible by looking at the cell it belongs to. There would be no way to + loop over all lines and for example do a contour integral, since there + would be no way to find out which of the lines we loop over actually are + on the contour. + + Since in one dimension, no substructures of lower dimension exist to + cells (of course apart from vertices, but these are handled + in another way than the structures and substructures with dimension one + or greater), there is no way to denote boundary indicators to boundary + vertices (the endpoints). This is not a big thing, however, since you + will normally not want to do a loop over all vertices, but rather work + on the cells, which in turn have a possibility to find out whether they + are at one of the endpoints. Only the handling of boundary values gets + a bit more complicated, but this seems to be the price to be paid for + the different handling of vertices from lines and quads. {\bf History of a triangulation} @@ -907,13 +1035,19 @@ class Triangulation : public TriaDimensionInfo { * the cell list should be useful (connected * domain, etc.). * + * Material data for the cells is given + * within the #cells# array, while boundary + * information is given in the + * #subcelldata# field. + * * The numbering of vertices within the * #cells# array is subject to some * constraints; see the general class * documentation for this. */ - void create_triangulation (const vector > &vertices, - const vector > &cells); + void create_triangulation (const vector > &vertices, + const vector > &cells, + const SubCellData &subcelldata); /** * Initialize the triangulation with a @@ -1603,6 +1737,18 @@ class Triangulation : public TriaDimensionInfo { /** * Exception */ + DeclException2 (ExcLineInexistant, + int, int, + << "When trying to give a boundary indicator to a line: " + << "the line with end vertices " << arg1 << " and " + << arg2 << " does not exist."); + /** + * Exception + */ + DeclException0 (ExcInteriorLineCantBeBoundary); + /** + * Exception + */ DeclException0 (ExcInternalError); //@} protected: diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h index 8478adb8df..bff8f3d57f 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.h +++ b/deal.II/deal.II/include/grid/tria_accessor.h @@ -206,6 +206,10 @@ class TriaAccessor { * Exception */ DeclException0 (ExcDereferenceInvalidObject); + /** + * Exception + */ + DeclException0 (ExcNotUsefulForThisDimension); /*@}*/ protected: @@ -365,6 +369,51 @@ class LineAccessor : public TriaAccessor { */ bool has_children () const; + /** + * Return the boundary indicator of this + * line. Since boundary data is only useful + * for structures with a dimension less + * than the dimension of a cell, this + * function issues an error if #dim<2#. + * + * If the return value is 255, then this + * line is in the interior of the domain. + */ + unsigned char boundary_indicator () const; + + /** + * Set the boundary indicator of this line. + * The same applies as for the + * #boundary_indicator()# function. + * + * Should be careful with this function + * and especially never try to set the + * boundary indicator to 255, unless + * you exactly know what you are doing, + * since this value is reserved for + * another purpose and algorithms may + * not work if boundary cells have have + * this boundary indicator or if interior + * cells have boundary indicators other + * than 255. + */ + void set_boundary_indicator (unsigned char) const; + + /** + * Return whether this line is at the + * boundary. This is checked via the + * the boundary indicator field, which + * is always 255 if the line is in the + * interior of the domain. Obviously, + * this is only possible for #dim>1#; + * however, for #dim==1#, a line is + * a cell and the #CellAccessor# class + * offers another possibility to + * determine whether a cell is at the + * boundary or not. + */ + bool at_boundary () const; + private: /** * Copy operator. This is normally @@ -544,6 +593,51 @@ class QuadAccessor : public TriaAccessor { */ bool has_children () const; + /** + * Return the boundary indicator of this + * line. Since boundary data is only useful + * for structures with a dimension less + * than the dimension of a cell, this + * function issues an error if #dim<3#. + * + * If the return value is 255, then this + * line is in the interior of the domain. + */ + unsigned char boundary_indicator () const; + + /** + * Set the boundary indicator of this line. + * The same applies as for the + * #boundary_indicator()# function. + * + * Should be careful with this function + * and especially never try to set the + * boundary indicator to 255, unless + * you exactly know what you are doing, + * since this value is reserved for + * another purpose and algorithms may + * not work if boundary cells have have + * this boundary indicator or if interior + * cells have boundary indicators other + * than 255. + */ + void set_boundary_indicator (unsigned char) const; + + /** + * Return whether this line is at the + * boundary. This is checked via the + * the boundary indicator field, which + * is always 255 if the line is in the + * interior of the domain. Obviously, + * this is only possible for #dim>2#; + * however, for #dim==2#, a quad is + * a cell and the #CellAccessor# class + * offers another possibility to + * determine whether a cell is at the + * boundary or not. + */ + bool at_boundary () const; + private: /** * Copy operator. This is normally @@ -668,7 +762,7 @@ class TriaSubstructAccessor<2> : public QuadAccessor<2> { /** - This class allows access to a cell, i.e. a line in one dimension, a quad + This class allows access to a cell: a line in one dimension, a quad in two dimension, etc. The following refers to any space dimension: @@ -722,7 +816,8 @@ class CellAccessor : public TriaSubstructAccessor { const TriaIterator > &pointer) const; /** - * Return whether the #i#th vertex is + * Return whether the #i#th vertex or + * face (depending on the dimension) is * part of the boundary. This is true, if * the #i#th neighbor does not exist. */ @@ -758,6 +853,16 @@ class CellAccessor : public TriaSubstructAccessor { */ TriaIterator > child (const unsigned int i) const; + /** + * Return the material id of this cell. + */ + unsigned char material_id () const; + + /** + * Set the material id of this cell. + */ + void set_material_id (const unsigned char) const; + /** * Test whether the cell has children * (this is the criterion for activity diff --git a/deal.II/deal.II/include/numerics/data_io.h b/deal.II/deal.II/include/numerics/data_io.h index b80676e563..6d35cb4956 100644 --- a/deal.II/deal.II/include/numerics/data_io.h +++ b/deal.II/deal.II/include/numerics/data_io.h @@ -25,6 +25,82 @@ template class map; template struct less; + + + +/** + Structure which is passed to the #Triangulation::create_triangulation# + function. It contains all data needed to construct a cell, namely the + indices of the vertices and the material indicator. +*/ +template +struct CellData { + int vertices[2<::create_triangulation# + function to describe boundary information. + + This structure is the same for all dimensions, since we use an input + function which is the same for all dimensions. The content of objects + of this structure varies with the dimensions, however. + + Since in one space dimension, there is no boundary information apart + from the two end points of the interval, this structure does not contain + anything and exists only for consistency, to allow a common interface + for all space dimensions. All fields should always be empty. + + Boundary data in 2D consists + of a list of lines which belong to a given boundary component. A + boundary component is a list of lines which are given a common + number describing the boundary condition to hold on this part of the + boundary. The triangulation creation function gives lines not in this + list either the boundary indicator zero (if on the boundary) or 255 + (if in the interior). Explicitely giving a line the indicator 255 + will result in an error, as well as giving an interior line a boundary + indicator. +*/ +struct SubCellData { + /** + * Each record of this vector describes + * a line on the boundary and its boundary + * indicator. + */ + vector > boundary_lines; + + /** + * Each record of this vector describes + * a quad on the boundary and its boundary + * indicator. + */ + vector > boundary_quads; + + /** + * This function checks whether the vectors + * which may not be used in a given + * dimension are really empty. I.e., + * whether the #boundary_*# arrays are + * empty when in one space dimension + * and whether the #boundary_quads# + * array is empty when in two dimensions. + * + * Since this structure is the same for all + * dimensions, the actual dimension has + * to be given as a parameter. + */ + bool check_consistency (const unsigned int dim) const; +}; + + + + + /** This class implements an input mechanism for grid data. It allows to read a grid structure into a triangulation object. Future versions @@ -38,12 +114,17 @@ template struct less; 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 + and line and quads in two dimensions are accepted. All other cell types + (e.g. triangles in two dimensions, quads and hexes in 3d) are rejected. + 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. + Material indicators are accepted to denote the material id of cells and + to denote boundary part indication for lines in 2D. Read the according + sections in the documentation of the \Ref{Triangulation} class for + further details. + {\bf Structure of input grid data} @@ -157,7 +238,12 @@ class DataIn { DeclException2 (ExcInvalidVertexIndex, int, int, << "Trying to access invalid vertex index " << arg2 - << " while creating cell " << arg1); + << " while creating cell " << arg1); + /** + * Exception + */ + DeclException0 (ExcInternalError); + private: /** * Store address of the triangulation to diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 2d0b3ba87d..fe955a0547 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -40,7 +41,12 @@ void Triangulation::set_boundary (const Boundary *boundary_object) { void Triangulation<1>::create_triangulation (const vector > &v, - const vector > &cells) { + const vector > &cells, + const SubCellData &) { + // note: since no boundary information + // can be given in one dimension, the + // #subcelldata# field is ignored. + const unsigned int dim=1; Assert (vertices.size() == 0, ExcTriangulationNotEmpty()); @@ -65,13 +71,14 @@ void Triangulation<1>::create_triangulation (const vector > &v, while (next_free_line->used()) ++next_free_line; - next_free_line->set (Line (cells[cell][0], cells[cell][1])); + next_free_line->set (Line (cells[cell].vertices[0], cells[cell].vertices[1])); next_free_line->set_used_flag (); + next_free_line->set_material_id (cells[cell].material_id); // 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); + lines_at_vertex[cells[cell].vertices[0]].push_back (cell); + lines_at_vertex[cells[cell].vertices[1]].push_back (cell); }; @@ -132,7 +139,8 @@ void Triangulation<1>::create_triangulation (const vector > &v, void Triangulation<2>::create_triangulation (const vector > &v, - const vector > &c) { + const vector > &c, + const SubCellData &subcelldata) { const unsigned int dim=2; Assert (vertices.size() == 0, ExcTriangulationNotEmpty()); @@ -145,7 +153,7 @@ void Triangulation<2>::create_triangulation (const vector > &v, // copy cells. This is needed since we // may need to change entries - vector > cells(c); + vector > cells(c); // make up a list of the needed lines @@ -165,16 +173,16 @@ void Triangulation<2>::create_triangulation (const vector > &v, // 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())); + Assert ((0<=cells[cell].vertices[vertex]) && + (cells[cell].vertices[vertex]<(signed int)vertices.size()), + ExcInvalidVertexIndex (cell, cells[cell].vertices[vertex], vertices.size())); #endif pair 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])}; + make_pair (cells[cell].vertices[0], cells[cell].vertices[1]), + make_pair (cells[cell].vertices[1], cells[cell].vertices[2]), + make_pair (cells[cell].vertices[0], cells[cell].vertices[3]), + make_pair (cells[cell].vertices[3], cells[cell].vertices[2])}; // note the following: if the sense // of the vertices of a cell is correct, @@ -193,17 +201,18 @@ void Triangulation<2>::create_triangulation (const vector > &v, // 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)) + line_vertices[line].first)) != needed_lines.end()) { // rotate vertex numbers - rotate (cells[cell].begin(), cells[cell].begin()+2, cells[cell].end()); + swap (cells[cell].vertices[0], cells[cell].vertices[2]); + swap (cells[cell].vertices[1], cells[cell].vertices[3]); // 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]); + line_vertices[0] = make_pair (cells[cell].vertices[0], cells[cell].vertices[1]); + line_vertices[1] = make_pair (cells[cell].vertices[1], cells[cell].vertices[2]); + line_vertices[2] = make_pair (cells[cell].vertices[0], cells[cell].vertices[3]); + line_vertices[3] = make_pair (cells[cell].vertices[3], cells[cell].vertices[2]); // allow for only one such // rotation break; @@ -290,7 +299,7 @@ void Triangulation<2>::create_triangulation (const vector > &v, // store for each line index // the adjacent cells - map, less > adjacent_cells; + map,less > adjacent_cells; // finally make up cells if (true) @@ -300,10 +309,10 @@ void Triangulation<2>::create_triangulation (const vector > &v, { // 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])]}; + needed_lines[make_pair(cells[c].vertices[0], cells[c].vertices[1])], + needed_lines[make_pair(cells[c].vertices[1], cells[c].vertices[2])], + needed_lines[make_pair(cells[c].vertices[3], cells[c].vertices[2])], + needed_lines[make_pair(cells[c].vertices[0], cells[c].vertices[3])]}; cell->set (Quad(lines[0]->index(), lines[1]->index(), @@ -311,7 +320,7 @@ void Triangulation<2>::create_triangulation (const vector > &v, lines[3]->index())); cell->set_used_flag (); - + cell->set_material_id (cells[c].material_id); // note that this cell is adjacent // to the four lines for (unsigned int line=0; line<4; ++line) @@ -331,18 +340,67 @@ void Triangulation<2>::create_triangulation (const vector > &v, }; -#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 + for (line_iterator line=begin_line(); line!=end_line(); ++line) + { + const unsigned int n_adj_cells = adjacent_cells[line->index()].size(); + // assert that every line has one or + // two adjacent cells + Assert ((n_adj_cells >= 1) && + (n_adj_cells <= 2), + ExcInternalError()); + + // if only one cell: line is at + // boundary -> give it the boundary + // indicator zero by default + if (n_adj_cells == 1) + line->set_boundary_indicator (0); + else + // interior line -> 255 + line->set_boundary_indicator (255); + }; + + // set boundary indicators where given + vector >::const_iterator boundary_line + = subcelldata.boundary_lines.begin(); + vector >::const_iterator end_boundary_line + = subcelldata.boundary_lines.end(); + for (; boundary_line!=end_boundary_line; ++boundary_line) + { + line_iterator line; + pair line_vertices(make_pair((*boundary_line).vertices[0], + (*boundary_line).vertices[1])); + if (needed_lines.find(line_vertices) != needed_lines.end()) + // line found in this direction + line = needed_lines[line_vertices]; + else + { + // look whether it exists in + // reverse direction + swap (line_vertices.first, line_vertices.second); + if (needed_lines.find(line_vertices) != needed_lines.end()) + line = needed_lines[line_vertices]; + else + { + // line does not exist + Assert (false, ExcLineInexistant(line_vertices.first, + line_vertices.second)); + line = end_line(); + }; + }; + + // Assert that only exterior lines + // are sgiven a boundary indicator + Assert (line->boundary_indicator() == 0, + ExcInteriorLineCantBeBoundary()); + + line->set_boundary_indicator ((*boundary_line).material_id); + }; + - // update neighborship info + + + // finally 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) @@ -367,21 +425,16 @@ 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.push_back (true); - vertices_used.push_back (true); + const Point<1> vertices[2] = { Point<1>(left), Point<1>(right) }; + vector > cells (1, CellData<1>()); + cells[0].vertices[0] = 0; + cells[0].vertices[1] = 1; + cells[0].material_id = 0; - levels[0]->lines.lines.push_back (Line(0,1)); - levels[0]->lines.children.push_back (-1); - levels[0]->lines.used.push_back (true); - levels[0]->lines.user_flags.push_back (false); - levels[0]->neighbors.push_back (make_pair(-1,-1)); - levels[0]->neighbors.push_back (make_pair(-1,-1)); - - levels[0]->refine_flags.push_back (false); + create_triangulation (vector >(&vertices[0], &vertices[2]), + cells, + SubCellData()); // no boundary information }; @@ -398,13 +451,14 @@ void Triangulation<2>::create_hypercube (const double left, Point<2>(right,right), Point<2>(left,right) }; const int cell_vertices[1][4] = { { 0,1,2,3 } }; - vector > cells (1, vector()); - cells[0].reserve (4); + vector > cells (1, CellData<2>()); for (unsigned int j=0; j<4; ++j) - cells[0].push_back (cell_vertices[0][j]); + cells[0].vertices[j] = cell_vertices[0][j]; + cells[0].material_id = 0; create_triangulation (vector >(&vertices[0], &vertices[4]), - cells); + cells, + SubCellData()); // no boundary information }; @@ -429,21 +483,18 @@ void Triangulation<2>::create_hyper_L (const double a, const double b) { {1, 2, 5, 4}, {3, 4, 7, 6}}; - // with gcc2.9, uncomment this line - // and delete following workaraound -// vector > cells (3, vector(4,-1)); - vector tmp; - tmp.reserve (4); - for (unsigned int i=0; i<4; ++i) - tmp.push_back (-1); - vector > cells (3, tmp); + vector > cells (3, CellData<2>()); - for (unsigned int i=0; i<3; ++i) - for (unsigned int j=0; j<4; ++j) - cells[i][j] = cell_vertices[i][j]; + for (unsigned int i=0; i<3; ++i) + { + for (unsigned int j=0; j<4; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + }; create_triangulation (vector >(&vertices[0], &vertices[8]), - cells); + cells, + SubCellData()); // no boundary information }; @@ -471,21 +522,18 @@ void Triangulation<2>::create_hyper_ball (const Point<2> p, const double radius) {1, 7, 5, 3}, {6, 4, 5, 7}}; - // with gcc2.9, uncomment this line - // and delete following workaraound -// vector > cells (5, vector(4,-1)); - vector tmp; - tmp.reserve (4); - for (unsigned int i=0; i<4; ++i) - tmp.push_back (-1); - vector > cells (5, tmp); - - for (unsigned int i=0; i<5; ++i) - for (unsigned int j=0; j<4; ++j) - cells[i][j] = cell_vertices[i][j]; + vector > cells (5, CellData<2>()); + + for (unsigned int i=0; i<5; ++i) + { + for (unsigned int j=0; j<4; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + }; create_triangulation (vector >(&vertices[0], &vertices[8]), - cells); + cells, + SubCellData()); // no boundary information }; @@ -1614,6 +1662,8 @@ void Triangulation<1>::execute_refinement () { cell->set_children (first_child->index()); first_child->clear_children (); first_child->set (Line (cell->vertex_index(0), next_unused_vertex)); + first_child->set_material_id (cell->material_id()); + // reset neighborship info // (refer to \Ref{TriangulationLevel<0>} // for details) @@ -1649,6 +1699,7 @@ void Triangulation<1>::execute_refinement () { second_child->clear_children (); second_child->set (Line (next_unused_vertex, cell->vertex_index(1))); second_child->set_neighbor (0, first_child); + second_child->set_material_id (cell->material_id()); if (cell->neighbor(1).state() != valid) second_child->set_neighbor (1, cell->neighbor(1)); else @@ -2143,6 +2194,24 @@ void Triangulation<2>::execute_refinement () { new_lines[11]->set_used_flag (); new_lines[11]->clear_children (); + // set the boundary indicators of + // the outer cells. + new_lines[0]->set_boundary_indicator (cell->line(0)->boundary_indicator()); + new_lines[1]->set_boundary_indicator (cell->line(0)->boundary_indicator()); + new_lines[2]->set_boundary_indicator (cell->line(1)->boundary_indicator()); + new_lines[3]->set_boundary_indicator (cell->line(1)->boundary_indicator()); + new_lines[4]->set_boundary_indicator (cell->line(2)->boundary_indicator()); + new_lines[5]->set_boundary_indicator (cell->line(2)->boundary_indicator()); + new_lines[6]->set_boundary_indicator (cell->line(3)->boundary_indicator()); + new_lines[7]->set_boundary_indicator (cell->line(3)->boundary_indicator()); + // inner cells have boundary + // indicator 255 + new_lines[8]->set_boundary_indicator (255); + new_lines[9]->set_boundary_indicator (255); + new_lines[10]->set_boundary_indicator (255); + new_lines[11]->set_boundary_indicator (255); + + // finally add the four new cells! // search for next unused cell @@ -2221,6 +2290,11 @@ void Triangulation<2>::execute_refinement () { subcells[3]->set_neighbor (1, subcells[2]); subcells[3]->set_neighbor (2, neighbors[5]); subcells[3]->set_neighbor (3, neighbors[6]); + + subcells[0]->set_material_id (cell->material_id()); + subcells[1]->set_material_id (cell->material_id()); + subcells[2]->set_material_id (cell->material_id()); + subcells[3]->set_material_id (cell->material_id()); }; }; #ifdef DEBUG @@ -2319,11 +2393,18 @@ void TriangulationLevel<1>::reserve_space (const unsigned int new_lines) { lines.used.insert (lines.used.end(), new_size-lines.used.size(), false); lines.user_flags.reserve (new_size); - lines.user_flags.insert (lines.user_flags.end(), new_size-lines.user_flags.size(), false); + lines.user_flags.insert (lines.user_flags.end(), + new_size-lines.user_flags.size(), false); lines.children.reserve (new_size); lines.children.insert (lines.children.end(), new_size-lines.children.size(), -1); + + lines.material_id.reserve (new_size); + lines.material_id.insert (lines.material_id.end(), + new_size-lines.material_id.size(), + 255); + }; @@ -2352,6 +2433,8 @@ void TriangulationLevel<1>::monitor_memory (const unsigned int true_dimension) c ExcMemoryInexact (lines.lines.size(), lines.user_flags.size())); Assert (lines.lines.size() == lines.children.size(), ExcMemoryInexact (lines.lines.size(), lines.children.size())); + Assert (lines.lines.size() == lines.material_id.size(), + ExcMemoryInexact (lines.lines.size(), lines.material_id.size())); Assert (lines.used[lines.used.size()-1]==true , ExcUnusedMemoryAtEnd()); @@ -2377,11 +2460,17 @@ void TriangulationLevel<2>::reserve_space (const unsigned int new_quads) { quads.used.insert (quads.used.end(), new_size-quads.used.size(), false); quads.user_flags.reserve (new_size); - quads.user_flags.insert (quads.user_flags.end(), new_size-quads.user_flags.size(), false); + quads.user_flags.insert (quads.user_flags.end(), + new_size-quads.user_flags.size(), false); quads.children.reserve (new_size); quads.children.insert (quads.children.end(), new_size-quads.children.size(), -1); + + quads.material_id.reserve (new_size); + quads.material_id.insert (quads.material_id.end(), + new_size-quads.material_id.size(), + 255); }; @@ -2409,6 +2498,8 @@ void TriangulationLevel<2>::monitor_memory (const unsigned int true_dimension) c ExcMemoryInexact (quads.quads.size(), quads.user_flags.size())); Assert (quads.quads.size() == quads.children.size(), ExcMemoryInexact (quads.quads.size(), quads.children.size())); + Assert (quads.quads.size() == quads.material_id.size(), + ExcMemoryInexact (quads.quads.size(), quads.material_id.size())); Assert (quads.used[quads.used.size()-1]==true , ExcUnusedMemoryAtEnd()); diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc index 4460a2816b..c9df0fa519 100644 --- a/deal.II/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/deal.II/source/grid/tria_accessor.cc @@ -241,6 +241,35 @@ void LineAccessor::operator -- () { +template +unsigned char LineAccessor::boundary_indicator () const { + Assert (dim>=2, ExcNotUsefulForThisDimension()); + Assert (used(), ExcRefineCellNotUsed()); + + return tria->levels[present_level]->lines.material_id[present_index]; +}; + + + +template +void LineAccessor::set_boundary_indicator (unsigned char boundary_ind) const { + Assert (dim>=2, ExcNotUsefulForThisDimension()); + Assert (used(), ExcRefineCellNotUsed()); + + tria->levels[present_level]->lines.material_id[present_index] = boundary_ind; +}; + + + +template +bool LineAccessor::at_boundary () const { + // error checking is done + // in boundary_indicator() + return (boundary_indicator() != 255); +}; + + + template class LineAccessor<1>; template class LineAccessor<2>; @@ -440,6 +469,33 @@ void QuadAccessor::operator -- () { +template +unsigned char QuadAccessor::boundary_indicator () const { + Assert (dim>=3, ExcNotUsefulForThisDimension()); + Assert (used(), ExcRefineCellNotUsed()); + + return tria->levels[present_level]->quads.material_id[present_index]; +}; + + + +template +void QuadAccessor::set_boundary_indicator (unsigned char boundary_ind) const { + Assert (dim>=3, ExcNotUsefulForThisDimension()); + Assert (used(), ExcRefineCellNotUsed()); + + tria->levels[present_level]->quads.material_id[present_index] = boundary_ind; +}; + + + +template +bool QuadAccessor::at_boundary () const { + // error checking is done + // in boundary_indicator() + return (boundary_indicator() != 255); +}; + template class QuadAccessor<2>; @@ -456,6 +512,22 @@ bool CellAccessor<1>::at_boundary () const { +unsigned char CellAccessor<1>::material_id () const { + Assert (used(), + typename TriaSubstructAccessor<1>::ExcRefineCellNotUsed()); + return tria->levels[present_level]->lines.material_id[present_index]; +}; + + + +void CellAccessor<1>::set_material_id (const unsigned char mat_id) const { + Assert (used(), + typename TriaSubstructAccessor<1>::ExcRefineCellNotUsed()); + tria->levels[present_level]->lines.material_id[present_index] + = mat_id; +}; + + /*------------------------ Functions: CellAccessor<2> -----------------------*/ @@ -466,6 +538,22 @@ bool CellAccessor<2>::at_boundary () const { +unsigned char CellAccessor<2>::material_id () const { + Assert (used(), + typename TriaSubstructAccessor<2>::ExcRefineCellNotUsed()); + return tria->levels[present_level]->quads.material_id[present_index]; +}; + + + +void CellAccessor<2>::set_material_id (const unsigned char mat_id) const { + Assert (used(), + typename TriaSubstructAccessor<2>::ExcRefineCellNotUsed()); + tria->levels[present_level]->quads.material_id[present_index] + = mat_id; +}; + + /*------------------------ Functions: CellAccessor -----------------------*/ @@ -505,7 +593,7 @@ bool CellAccessor::at_boundary (const unsigned int i) const { Assert (i<2*dim, typename TriaSubstructAccessor::ExcInvalidIndex (i,0,2*dim-1)); - return (neighbor(i).state() != valid); + return (neighbor_index(i) == -1); }; diff --git a/deal.II/deal.II/source/numerics/data_io.cc b/deal.II/deal.II/source/numerics/data_io.cc index 63d2417ca7..393f8af930 100644 --- a/deal.II/deal.II/source/numerics/data_io.cc +++ b/deal.II/deal.II/source/numerics/data_io.cc @@ -13,6 +13,21 @@ +bool SubCellData::check_consistency (const unsigned int dim) const { + switch (dim) + { + case 1: + return ((boundary_lines.size() == 0) && + (boundary_quads.size() == 0)); + case 2: + return (boundary_quads.size() == 0); + }; + return true; +}; + + + + template DataIn::DataIn () : @@ -83,43 +98,76 @@ void DataIn::read_ucd (istream &in) { }; // set up array of cells - vector > cells; + vector > cells; + SubCellData subcelldata; for (unsigned int cell=0; cell> dummy // cell number - >> dummy; // material id + >> 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 + if (((cell_type == "line") && (dim == 1)) || + ((cell_type == "quad") && (dim == 2))) + // found a cell { // allocate and read indices - cells.push_back (vector (1<()); for (unsigned int i=0; i<(1<> cells.back()[i]; + in >> cells.back().vertices[i]; + cells.back().material_id = material_id; // transform from ucd to // consecutive numbering for (unsigned int i=0; i<(1<()); + in >> subcelldata.boundary_lines.back().vertices[0] + >> subcelldata.boundary_lines.back().vertices[1]; + subcelldata.boundary_lines.back().material_id = material_id; + + // transform from ucd to + // consecutive numbering + for (unsigned int i=0; i<2; ++i) + if (vertex_indices.find (subcelldata.boundary_lines.back().vertices[i]) != + vertex_indices.end()) + // vertex with this index exists + subcelldata.boundary_lines.back().vertices[i] + = vertex_indices[subcelldata.boundary_lines.back().vertices[i]]; + else + { + // no such vertex index + Assert (false, + ExcInvalidVertexIndex(cell, + subcelldata.boundary_lines.back().vertices[i])); + subcelldata.boundary_lines.back().vertices[i] = -1; + }; + } + else + // cannot read this + Assert (false, ExcNotImplemented()); }; - tria->create_triangulation (vertices, cells); + // check that no forbidden arrays are used + Assert (subcelldata.check_consistency(dim), ExcInternalError()); + + tria->create_triangulation (vertices, cells, subcelldata); };