static const unsigned int hexes_per_cell = (2*GeometryInfo<dim-1>::hexes_per_cell
+ GeometryInfo<dim-1>::quads_per_cell);
+ /**
+ * List of numbers which is
+ * denote which face is opposite
+ * to a given face. In 1d, this
+ * list is #{1,0}#, in 2d #{2, 3, 0, 1}#,
+ * in 3d #{1, 0, 4, 5, 2, 3}#.
+ */
+ static const unsigned int opposite_face[faces_per_cell];
+
/**
* This field store which child cells
* are adjacent to a certain face of
* Exception
*/
DeclException0 (ExcCantCompareIterators);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNeighborIsCoarser);
/*@}*/
protected:
* the possibility to check whether they are at the boundary etc. This class
* offers access to all this data.
*
- * @author Wolfgang Bangerth, 1998
+ * @author Wolfgang Bangerth, 1998, 1999, 2000
*/
template <int dim>
-class CellAccessor : public TriaObjectAccessor<dim,dim> {
+class CellAccessor : public TriaObjectAccessor<dim,dim>
+{
public:
/**
* Propagate the AccessorData type
void set_neighbor (const unsigned int i,
const TriaIterator<dim,CellAccessor<dim> > &pointer) const;
+ /**
+ * Return the how-many'th
+ * neighbor this cell is of
+ * #cell->neighbor(neighbor)#,
+ * i.e. return the number #n#
+ * such that
+ * #cell->neighbor(neighbor)->neighbor(n)==cell#. This
+ * function is the right one if
+ * you want to know how to get
+ * back from a neighbor to the
+ * present cell.
+ *
+ * Note that this operation is
+ * only useful if the neighbor is
+ * not on a coarser level than
+ * the present cell
+ * (i.e. #cell->neighbor(neighbor)->level()#
+ * needs to be equal to
+ * #cell->level()#, since
+ * otherwise the neighbors of the
+ * neighbor cell are on a coarser
+ * level than the present one and
+ * you can't get back from there
+ * to this cell.
+ */
+ unsigned int neighbor_of_neighbor (const unsigned int neighbor) const;
+
/**
* Return whether the #i#th vertex or
* face (depending on the dimension) is
//const unsigned int GeometryInfo<deal_II_dimension>::children_per_cell;
+template <>
+const unsigned int GeometryInfo<1>::opposite_face[GeometryInfo<1>::faces_per_cell]
+= { 0, 1 };
+
+
+template <>
+const unsigned int GeometryInfo<2>::opposite_face[GeometryInfo<2>::faces_per_cell]
+= { 2, 3, 0, 1 };
+
+
+template <>
+const unsigned int GeometryInfo<3>::opposite_face[GeometryInfo<3>::faces_per_cell]
+= { 1, 0, 4, 5, 2, 3 };
// two children which
// we can use.
cell_iterator neighbor = cell->neighbor(nb);
- for (unsigned int nb_nb=0; nb_nb<4; ++nb_nb)
- if (neighbor->neighbor(nb_nb)==cell)
// this cell is the nb_nb-th
// neighbor or neighbor(nb)
- {
- neighbors_neighbor[2*nb] = neighbors_neighbor[2*nb+1] = nb_nb;
- // vertex 1 of child 0
- // is always the interior
- // one
- new_vertices[2*nb+1] = neighbor->line(nb_nb)
- ->child(0)->vertex_index(1);
-
- if (nb < 2)
- {
- new_lines[2*nb] = neighbor->line(nb_nb)->child(0);
- new_lines[2*nb+1]= neighbor->line(nb_nb)->child(1);
- } else {
- // lines 2 and 3 have
- // opposite sense
- new_lines[2*nb] = neighbor->line(nb_nb)->child(1);
- new_lines[2*nb+1]= neighbor->line(nb_nb)->child(0);
- };
-
- // finally find out which
- // are the two neighbor
- // subcells, adjacent to
- // the two sublines
- const int child_mapping[4][2] = {{0,1},{1,2},{3,2},{0,3}};
- if (nb < 2)
- {
- neighbors[2*nb] = neighbor->child(child_mapping[nb_nb][0]);
- neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][1]);
- } else {
- neighbors[2*nb] = neighbor->child(child_mapping[nb_nb][1]);
- neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][0]);
- };
- };
+ const unsigned int nb_nb = cell->neighbor_of_neighbor (nb);
+
+ neighbors_neighbor[2*nb] = neighbors_neighbor[2*nb+1] = nb_nb;
+ // vertex 1 of child 0
+ // is always the interior
+ // one
+ new_vertices[2*nb+1] = neighbor->line(nb_nb)
+ ->child(0)->vertex_index(1);
+
+ if (nb < 2)
+ {
+ new_lines[2*nb] = neighbor->line(nb_nb)->child(0);
+ new_lines[2*nb+1]= neighbor->line(nb_nb)->child(1);
+ } else {
+ // lines 2 and 3 have
+ // opposite sense
+ new_lines[2*nb] = neighbor->line(nb_nb)->child(1);
+ new_lines[2*nb+1]= neighbor->line(nb_nb)->child(0);
+ };
+
+ // finally find out which
+ // are the two neighbor
+ // subcells, adjacent to
+ // the two sublines
+ static const unsigned int child_mapping[4][2] = {{0,1},{1,2},{3,2},{0,3}};
+ if (nb < 2)
+ {
+ neighbors[2*nb] = neighbor->child(child_mapping[nb_nb][0]);
+ neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][1]);
+ } else {
+ neighbors[2*nb] = neighbor->child(child_mapping[nb_nb][1]);
+ neighbors[2*nb+1]= neighbor->child(child_mapping[nb_nb][0]);
+ };
}
else
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell;
++face)
{
- cell_iterator neighbor = hex->neighbor(face);
+ const cell_iterator neighbor = hex->neighbor(face);
// if no neighbor
if (neighbor.state() != valid)
// of the neighbor
// adjacent to which
// the present cell is
- unsigned int nb_nb;
- for (nb_nb=0; nb_nb<GeometryInfo<dim>::faces_per_cell;
- ++nb_nb)
- if (neighbor->neighbor(nb_nb) == hex)
- break;
+ const unsigned int nb_nb = hex->neighbor_of_neighbor(face);
Assert (nb_nb<GeometryInfo<dim>::faces_per_cell,
ExcInternalError());
+template <int dim>
+unsigned int CellAccessor<dim>::neighbor_of_neighbor (const unsigned int neighbor) const
+{
+ // make sure that the neighbor is
+ // not on a coarser level
+ Assert (neighbor_level(neighbor) == present_level,
+ ExcNeighborIsCoarser());
+ Assert (neighbor < GeometryInfo<dim>::faces_per_cell,
+ ExcInvalidNeighbor(neighbor));
+
+ const TriaIterator<dim,CellAccessor<dim> > neighbor_cell = this->neighbor(neighbor);
+
+ // usually, on regular patches of
+ // the grid, this cell is just on
+ // the opposite side of the
+ // neighbor that the neighbor is of
+ // this cell. for example in 2d, if
+ // we want to know the
+ // neighbor_of_neighbor if
+ // neighbor==1 (the right
+ // neighbor), then we will get 3
+ // (the left neighbor) in most
+ // cases. look up this relationship
+ // in the table provided by
+ // GeometryInfo and try it
+ const unsigned int neighbor_guess
+ = GeometryInfo<dim>::opposite_face[neighbor];
+
+ if ((neighbor_cell->neighbor_index (neighbor_guess) == present_index) &&
+ (neighbor_cell->neighbor_level (neighbor_guess) == present_level))
+ return neighbor_guess;
+ else
+ // if the guess was false, then
+ // we need to loop over all
+ // neighbors and find the number
+ // the hard way
+ {
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if ((neighbor_cell->neighbor_index (face) == present_index) &&
+ (neighbor_cell->neighbor_level (face) == present_level))
+ return face;
+
+ // we should never get here,
+ // since then we did not find
+ // our way back...
+ Assert (false, ExcInternalError());
+ return static_cast<unsigned int>(-1);
+ };
+};
+
+
+
template <int dim>
bool CellAccessor<dim>::at_boundary (const unsigned int i) const {
Assert (used(), ExcCellNotUsed());
// of gradient across this face
{
Assert (cell->neighbor(face_no).state() == valid,
- ExcInternalError());
- unsigned int neighbor_neighbor;
- DoFHandler<dim>::active_cell_iterator neighbor = cell->neighbor(face_no);
+ ExcInternalError());
+
+ const DoFHandler<dim>::active_cell_iterator neighbor = cell->neighbor(face_no);
// find which number the current
// face has relative to the neighboring
// cell
- for (neighbor_neighbor=0; neighbor_neighbor<GeometryInfo<dim>::faces_per_cell;
- ++neighbor_neighbor)
- if (neighbor->neighbor(neighbor_neighbor) == cell)
- break;
-
+ const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell, ExcInternalError());
// get restriction of finite element
FESubfaceValues<dim> &fe_subface_values)
{
const DoFHandler<dim>::cell_iterator neighbor = cell->neighbor(face_no);
+
Assert (neighbor.state() == valid, ExcInternalError());
Assert (neighbor->has_children(), ExcInternalError());
// set up a vector of the gradients
// store which number #cell# has in the
// list of neighbors of #neighbor#
- unsigned int neighbor_neighbor;
- for (neighbor_neighbor=0; neighbor_neighbor<GeometryInfo<dim>::faces_per_cell;
- ++neighbor_neighbor)
- if (neighbor->neighbor(neighbor_neighbor) == cell)
- break;
+ const unsigned int neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
Assert (neighbor_neighbor<GeometryInfo<dim>::faces_per_cell, ExcInternalError());
// loop over all subfaces