template <int dim> class DoFHandler;
+template <int dim> struct CellData;
+struct SubCellData;
+
class istream;
class ostream;
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.
*/
* to the mother cell of this cell).
*/
vector<pair<int,int> > neighbors;
-
+
/**
* Reserve enough space to accomodate
* #total_cells# cells on this level.
* #Triangulation<>::clear_user_flags()#.
*/
vector<bool> 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<unsigned char> material_id;
};
public:
* Same as for #LineData::used#.
*/
vector<bool> 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<unsigned char> material_id;
};
public:
#Triangulation<dim>::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}
* 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<Point<dim> > &vertices,
- const vector<vector<int> > &cells);
+ void create_triangulation (const vector<Point<dim> > &vertices,
+ const vector<CellData<dim> > &cells,
+ const SubCellData &subcelldata);
/**
* Initialize the triangulation with a
/**
* 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:
* Exception
*/
DeclException0 (ExcDereferenceInvalidObject);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNotUsefulForThisDimension);
/*@}*/
protected:
*/
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
*/
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
/**
- 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:
const TriaIterator<dim,CellAccessor<dim> > &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.
*/
*/
TriaIterator<dim,CellAccessor<dim> > 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
template <class T> struct less;
+
+
+
+/**
+ Structure which is passed to the #Triangulation<dim>::create_triangulation#
+ function. It contains all data needed to construct a cell, namely the
+ indices of the vertices and the material indicator.
+*/
+template <int dim>
+struct CellData {
+ int vertices[2<<dim];
+ unsigned char material_id;
+};
+
+
+
+
+
+/**
+ Structure to be passed to the #Triangulation<dim>::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<CellData<1> > boundary_lines;
+
+ /**
+ * Each record of this vector describes
+ * a quad on the boundary and its boundary
+ * indicator.
+ */
+ vector<CellData<2> > 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
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}
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
#include <grid/tria_iterator.h>
#include <grid/tria.h>
#include <basic/magic_numbers.h>
+#include <basic/data_io.h>
#include <iostream.h>
#include <algo.h>
#include <map.h>
void Triangulation<1>::create_triangulation (const vector<Point<1> > &v,
- const vector<vector<int> > &cells) {
+ const vector<CellData<1> > &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());
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);
};
void Triangulation<2>::create_triangulation (const vector<Point<2> > &v,
- const vector<vector<int> > &c) {
+ const vector<CellData<2> > &c,
+ const SubCellData &subcelldata) {
const unsigned int dim=2;
Assert (vertices.size() == 0, ExcTriangulationNotEmpty());
// copy cells. This is needed since we
// may need to change entries
- vector<vector<int> > cells(c);
+ vector<CellData<2> > cells(c);
// make up a list of the needed lines
// 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<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])};
+ 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,
// 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;
// store for each line index
// the adjacent cells
- map<int,vector<cell_iterator>, less<int> > adjacent_cells;
+ map<int,vector<cell_iterator>,less<int> > adjacent_cells;
// finally make up cells
if (true)
{
// 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(),
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)
};
-#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
+ 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<CellData<1> >::const_iterator boundary_line
+ = subcelldata.boundary_lines.begin();
+ vector<CellData<1> >::const_iterator end_boundary_line
+ = subcelldata.boundary_lines.end();
+ for (; boundary_line!=end_boundary_line; ++boundary_line)
+ {
+ line_iterator line;
+ pair<int,int> 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)
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<CellData<1> > 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<Point<1> >(&vertices[0], &vertices[2]),
+ cells,
+ SubCellData()); // no boundary information
};
Point<2>(right,right),
Point<2>(left,right) };
const int cell_vertices[1][4] = { { 0,1,2,3 } };
- vector<vector<int> > cells (1, vector<int>());
- cells[0].reserve (4);
+ vector<CellData<2> > 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<Point<2> >(&vertices[0], &vertices[4]),
- cells);
+ cells,
+ SubCellData()); // no boundary information
};
{1, 2, 5, 4},
{3, 4, 7, 6}};
- // with gcc2.9, uncomment this line
- // and delete following workaraound
-// vector<vector<int> > cells (3, vector<int>(4,-1));
- vector<int> tmp;
- tmp.reserve (4);
- for (unsigned int i=0; i<4; ++i)
- tmp.push_back (-1);
- vector<vector<int> > cells (3, tmp);
+ vector<CellData<2> > 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<Point<dim> >(&vertices[0], &vertices[8]),
- cells);
+ cells,
+ SubCellData()); // no boundary information
};
{1, 7, 5, 3},
{6, 4, 5, 7}};
- // with gcc2.9, uncomment this line
- // and delete following workaraound
-// vector<vector<int> > cells (5, vector<int>(4,-1));
- vector<int> tmp;
- tmp.reserve (4);
- for (unsigned int i=0; i<4; ++i)
- tmp.push_back (-1);
- vector<vector<int> > 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<CellData<2> > 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<Point<dim> >(&vertices[0], &vertices[8]),
- cells);
+ cells,
+ SubCellData()); // no boundary information
};
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)
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
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
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
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);
+
};
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());
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);
};
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());
+template <int dim>
+unsigned char LineAccessor<dim>::boundary_indicator () const {
+ Assert (dim>=2, ExcNotUsefulForThisDimension());
+ Assert (used(), ExcRefineCellNotUsed());
+
+ return tria->levels[present_level]->lines.material_id[present_index];
+};
+
+
+
+template <int dim>
+void LineAccessor<dim>::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 <int dim>
+bool LineAccessor<dim>::at_boundary () const {
+ // error checking is done
+ // in boundary_indicator()
+ return (boundary_indicator() != 255);
+};
+
+
+
template class LineAccessor<1>;
template class LineAccessor<2>;
+template <int dim>
+unsigned char QuadAccessor<dim>::boundary_indicator () const {
+ Assert (dim>=3, ExcNotUsefulForThisDimension());
+ Assert (used(), ExcRefineCellNotUsed());
+
+ return tria->levels[present_level]->quads.material_id[present_index];
+};
+
+
+
+template <int dim>
+void QuadAccessor<dim>::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 <int dim>
+bool QuadAccessor<dim>::at_boundary () const {
+ // error checking is done
+ // in boundary_indicator()
+ return (boundary_indicator() != 255);
+};
+
template class QuadAccessor<2>;
+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> -----------------------*/
+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<dim> -----------------------*/
Assert (i<2*dim,
typename TriaSubstructAccessor<dim>::ExcInvalidIndex (i,0,2*dim-1));
- return (neighbor(i).state() != valid);
+ return (neighbor_index(i) == -1);
};
+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 <int dim>
DataIn<dim>::DataIn () :
};
// set up array of cells
- vector<vector<int> > cells;
+ vector<CellData<dim> > cells;
+ SubCellData subcelldata;
for (unsigned int cell=0; cell<n_cells; ++cell)
{
String cell_type;
-
+ int material_id;
+
in >> 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<int> (1<<dim));
+ cells.push_back (CellData<dim>());
for (unsigned int i=0; i<(1<<dim); ++i)
- in >> 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<<dim); ++i)
- if (vertex_indices.find (cells.back()[i]) != vertex_indices.end())
+ if (vertex_indices.find (cells.back().vertices[i]) != vertex_indices.end())
// vertex with this index exists
- cells.back()[i] = vertex_indices[cells.back()[i]];
+ cells.back().vertices[i] = vertex_indices[cells.back().vertices[i]];
else
{
// no such vertex index
- Assert (false, ExcInvalidVertexIndex(cell, cells.back()[i]));
- cells.back()[i] = -1;
+ Assert (false, ExcInvalidVertexIndex(cell, cells.back().vertices[i]));
+ cells.back().vertices[i] = -1;
};
- };
+ }
+ else
+ if ((cell_type == "line") && (dim == 2))
+ // boundary info
+ {
+ subcelldata.boundary_lines.push_back (CellData<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);
};