]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add support for material and boundary indicators.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 24 Mar 1998 17:11:10 +0000 (17:11 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 24 Mar 1998 17:11:10 +0000 (17:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@94 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 5bd827fa3a55455fc9b53b0e142c4367b9179ceb..a7fd1731199cc94ab965d777bfde57aef1491d63 100644 (file)
@@ -33,6 +33,9 @@ template <int dim, class Accessor> class TriaActiveIterator;
 
 template <int dim> class DoFHandler;
 
+template <int dim> 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<pair<int,int> > 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<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:
@@ -337,6 +372,29 @@ class TriangulationLevel<2> :  public TriangulationLevel<1>
                                          *  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:
@@ -704,6 +762,76 @@ class TriaDimensionInfo<2> {
         #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}
@@ -907,13 +1035,19 @@ class Triangulation : public TriaDimensionInfo<dim> {
                                      * 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
@@ -1603,6 +1737,18 @@ class Triangulation : public TriaDimensionInfo<dim> {
                                     /**
                                      * 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:
index 8478adb8df4e42285ee72135ec70111b4c029d1b..bff8f3d57f8e3896f3e4a3046a6f37a4c94bd2fb 100644 (file)
@@ -206,6 +206,10 @@ class TriaAccessor {
                                      *  Exception
                                      */
     DeclException0 (ExcDereferenceInvalidObject);
+                                    /**
+                                     *  Exception
+                                     */
+    DeclException0 (ExcNotUsefulForThisDimension);
                                     /*@}*/
        
   protected:
@@ -365,6 +369,51 @@ class LineAccessor :  public TriaAccessor<dim> {
                                      */
     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<dim> {
                                      */
     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<dim> {
                       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.
                                      */
@@ -758,6 +853,16 @@ class CellAccessor :  public TriaSubstructAccessor<dim> {
                                      */
     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
index b80676e5631ad92d3de1d169f5081c12550a348f..6d35cb4956114736b7c3bb45c4ae7d9d5215580e 100644 (file)
@@ -25,6 +25,82 @@ template <class Key, class T, class Compare, class Alloc = alloc> class map;
 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
@@ -38,12 +114,17 @@ template <class T> 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
index 2d0b3ba87de0a2a079e5f363474f2302687332b3..fe955a0547b639ed5019f614c18ec573ad39fecc 100644 (file)
@@ -4,6 +4,7 @@
 #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>
@@ -40,7 +41,12 @@ void Triangulation<dim>::set_boundary (const Boundary<dim> *boundary_object) {
 
 
 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());
@@ -65,13 +71,14 @@ void Triangulation<1>::create_triangulation (const vector<Point<1> >    &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<Point<1> >    &v,
 
 
 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());
@@ -145,7 +153,7 @@ void Triangulation<2>::create_triangulation (const vector<Point<2> >    &v,
 
                                   // 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
@@ -165,16 +173,16 @@ void Triangulation<2>::create_triangulation (const vector<Point<2> >    &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<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,
@@ -193,17 +201,18 @@ void Triangulation<2>::create_triangulation (const vector<Point<2> >    &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<Point<2> >    &v,
 
                                   // 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) 
@@ -300,10 +309,10 @@ void Triangulation<2>::create_triangulation (const vector<Point<2> >    &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<Point<2> >    &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<Point<2> >    &v,
     };
   
 
-#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)
@@ -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<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
 };
 
 
@@ -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<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
 };
 
 
@@ -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<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
 };
 
 
@@ -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<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
 };
 
 
@@ -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());
 
index 4460a2816bda33525c0a8a40f4e1d2d79b4d5496..c9df0fa51959909138c3a00f79c0973c5dcbe7ac 100644 (file)
@@ -241,6 +241,35 @@ void LineAccessor<dim>::operator -- () {
 
 
 
+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>;
@@ -440,6 +469,33 @@ void QuadAccessor<dim>::operator -- () {
 
 
 
+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>;
 
@@ -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<dim> -----------------------*/
 
 
@@ -505,7 +593,7 @@ bool CellAccessor<dim>::at_boundary (const unsigned int i) const {
   Assert (i<2*dim,
          typename TriaSubstructAccessor<dim>::ExcInvalidIndex (i,0,2*dim-1));
   
-  return (neighbor(i).state() != valid);
+  return (neighbor_index(i) == -1);
 };
 
 
index 63d2417ca70bda79c4a26aa8f1254d615907ab2f..393f8af9302fe3824cb1fe830dd560417e81b8fc 100644 (file)
 
 
 
+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 () :
@@ -83,43 +98,76 @@ void DataIn<dim>::read_ucd (istream &in) {
     };
 
                                   // 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);
 };
 
 

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.