From: Wolfgang Bangerth Date: Mon, 28 Nov 2016 00:15:13 +0000 (-0700) Subject: Find a better way to be dimension independent. X-Git-Tag: v8.5.0-rc1~358^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e2dae1060ddde4a4e709c0644b8757b193cc0556;p=dealii.git Find a better way to be dimension independent. This patch revises commit 9386b1e30c in #3639. There, I tried to avoid instantiating types that have a template signature of the kind <3,spacedim> if spacedim is less than three, since that leads to invalid types. I did so by using <3,max(3,spacedim)>, which works because we can only get to these places if dim==3, and consequently spacedim>=3, but the compiler may not always recognize that if spacedim<3, then we must be in dead code. This patch goes the other way around, which I find conceptually clearer: in code where we have if (dim == 3) { typename Triangulation<3,spacedim>::cell_iterator cell = ...; we clearly only got into the 'if' block because dim==3. Consequently, we might as well have written if (dim == 3) { typename Triangulation::cell_iterator cell = ...; to the same effect. On the plus side, however, whenever dim<3 (i.e. in cases where this is dead code anyway), we still do not instantiate invalid types. This patch therefore undoes the changes made in #3639 and replaces them with the second option above. There are a number of places where we then have to replace calls to GeometryInfo<3>::... by GeometryInfo::... as well as one place where we have to do a cast from RefinementCase<3> to RefinementCase or the other way around; the latter is the identity cast if we do get into this block. --- diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 7765c5d7c5..a55f4db938 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1773,7 +1773,7 @@ CellAccessor::neighbor_of_coarser_neighbor (const unsigned int ne case 2: { const int this_face_index=face_index(neighbor); - const TriaIterator > neighbor_cell = this->neighbor(neighbor); + const TriaIterator > neighbor_cell = this->neighbor(neighbor); // usually, on regular patches of // the grid, this cell is just on @@ -1791,7 +1791,7 @@ CellAccessor::neighbor_of_coarser_neighbor (const unsigned int ne const unsigned int face_no_guess = GeometryInfo<2>::opposite_face[neighbor]; - const TriaIterator > face_guess + const TriaIterator > face_guess =neighbor_cell->face(face_no_guess); if (face_guess->has_children()) @@ -1807,7 +1807,7 @@ CellAccessor::neighbor_of_coarser_neighbor (const unsigned int ne { if (face_no!=face_no_guess) { - const TriaIterator > face + const TriaIterator > face =neighbor_cell->face(face_no); if (face->has_children()) for (unsigned int subface_no=0; subface_non_children(); ++subface_no) @@ -1827,7 +1827,7 @@ CellAccessor::neighbor_of_coarser_neighbor (const unsigned int ne case 3: { const int this_face_index=face_index(neighbor); - const TriaIterator > + const TriaIterator > neighbor_cell = this->neighbor(neighbor); // usually, on regular patches of the grid, this cell is just on the @@ -1839,7 +1839,7 @@ CellAccessor::neighbor_of_coarser_neighbor (const unsigned int ne const unsigned int face_no_guess = GeometryInfo<3>::opposite_face[neighbor]; - const TriaIterator > face_guess + const TriaIterator > face_guess =neighbor_cell->face(face_no_guess); if (face_guess->has_children()) @@ -1867,7 +1867,7 @@ CellAccessor::neighbor_of_coarser_neighbor (const unsigned int ne if (face_no==face_no_guess) continue; - const TriaIterator > face + const TriaIterator > face =neighbor_cell->face(face_no); if (!face->has_children()) @@ -2347,24 +2347,21 @@ neighbor_child_on_subface (const unsigned int face, // | | 1 | | 0 | | | 0 | | 0 | 1 | // *---*---* *---*---* *-------* *---*---* - // in the following line, we know that if we got here that we - // have spacedim>=3. avoid instantiating invalid classes - // otherwise, even though we know that it is in dead code - const typename Triangulation<3,(spacedim<3 ? 3 : spacedim)>::face_iterator + const typename Triangulation::face_iterator mother_face = this->face(face); const unsigned int total_children=mother_face->number_of_children(); Assert (subface::max_children_per_face, ExcInternalError()); unsigned int neighbor_neighbor; - TriaIterator > neighbor_child; - const TriaIterator > neighbor + TriaIterator > neighbor_child; + const TriaIterator > neighbor = this->neighbor(face); - const RefinementCase<2> mother_face_ref_case + const RefinementCase mother_face_ref_case = mother_face->refinement_case(); - if (mother_face_ref_case==RefinementCase<2>::cut_xy) // total_children==4 + if (mother_face_ref_case==static_cast >(RefinementCase<2>::cut_xy)) // total_children==4 { // this case is quite easy. we are sure, // that the neighbor is not coarser. @@ -2377,11 +2374,11 @@ neighbor_child_on_subface (const unsigned int face, // now use the info provided by GeometryInfo // to extract the neighbors child number const unsigned int neighbor_child_index - = GeometryInfo<3>::child_cell_on_face(neighbor->refinement_case(), - neighbor_neighbor, subface, - neighbor->face_orientation(neighbor_neighbor), - neighbor->face_flip(neighbor_neighbor), - neighbor->face_rotation(neighbor_neighbor)); + = GeometryInfo::child_cell_on_face(neighbor->refinement_case(), + neighbor_neighbor, subface, + neighbor->face_orientation(neighbor_neighbor), + neighbor->face_flip(neighbor_neighbor), + neighbor->face_rotation(neighbor_neighbor)); neighbor_child = neighbor->child(neighbor_child_index); // make sure that the neighbor child cell we @@ -2467,24 +2464,24 @@ neighbor_child_on_subface (const unsigned int face, iso_subface=first_child_to_find + 2*indices.second; } neighbor_child_index - = GeometryInfo<3>::child_cell_on_face(neighbor->refinement_case(), - neighbor_neighbor, - iso_subface, - neighbor->face_orientation(neighbor_neighbor), - neighbor->face_flip(neighbor_neighbor), - neighbor->face_rotation(neighbor_neighbor)); + = GeometryInfo::child_cell_on_face(neighbor->refinement_case(), + neighbor_neighbor, + iso_subface, + neighbor->face_orientation(neighbor_neighbor), + neighbor->face_flip(neighbor_neighbor), + neighbor->face_rotation(neighbor_neighbor)); } else //neighbor is not coarser { neighbor_neighbor=neighbor_of_neighbor(face); neighbor_child_index - = GeometryInfo<3>::child_cell_on_face(neighbor->refinement_case(), - neighbor_neighbor, - first_child_to_find, - neighbor->face_orientation(neighbor_neighbor), - neighbor->face_flip(neighbor_neighbor), - neighbor->face_rotation(neighbor_neighbor), - mother_face_ref_case); + = GeometryInfo::child_cell_on_face(neighbor->refinement_case(), + neighbor_neighbor, + first_child_to_find, + neighbor->face_orientation(neighbor_neighbor), + neighbor->face_flip(neighbor_neighbor), + neighbor->face_rotation(neighbor_neighbor), + mother_face_ref_case); } neighbor_child=neighbor->child(neighbor_child_index); @@ -2493,11 +2490,11 @@ neighbor_child_on_subface (const unsigned int face, // along the given subface. go down that // list and deliver the last of those. while (neighbor_child->has_children() && - GeometryInfo<3>::face_refinement_case(neighbor_child->refinement_case(), - neighbor_neighbor) + GeometryInfo::face_refinement_case(neighbor_child->refinement_case(), + neighbor_neighbor) == RefinementCase<2>::no_refinement) neighbor_child = - neighbor_child->child(GeometryInfo<3>:: + neighbor_child->child(GeometryInfo:: child_cell_on_face(neighbor_child->refinement_case(), neighbor_neighbor, 0)); @@ -2514,7 +2511,7 @@ neighbor_child_on_subface (const unsigned int face, { if (subface<2) neighbor_child = - neighbor_child->child(GeometryInfo<3>:: + neighbor_child->child(GeometryInfo:: child_cell_on_face(neighbor_child->refinement_case(), neighbor_neighbor,subface, neighbor_child->face_orientation(neighbor_neighbor), @@ -2527,7 +2524,7 @@ neighbor_child_on_subface (const unsigned int face, Assert(mother_face->child(1)->has_children(), ExcInternalError()); if (subface>0) neighbor_child = - neighbor_child->child(GeometryInfo<3>:: + neighbor_child->child(GeometryInfo:: child_cell_on_face(neighbor_child->refinement_case(), neighbor_neighbor,subface-1, neighbor_child->face_orientation(neighbor_neighbor), @@ -2539,7 +2536,7 @@ neighbor_child_on_subface (const unsigned int face, else if (total_children==4) { neighbor_child = - neighbor_child->child(GeometryInfo<3>:: + neighbor_child->child(GeometryInfo:: child_cell_on_face(neighbor_child->refinement_case(), neighbor_neighbor,subface%2, neighbor_child->face_orientation(neighbor_neighbor), @@ -2555,7 +2552,7 @@ neighbor_child_on_subface (const unsigned int face, // deliver the last of those. while (neighbor_child->has_children()) neighbor_child - = neighbor_child->child(GeometryInfo<3>:: + = neighbor_child->child(GeometryInfo:: child_cell_on_face(neighbor_child->refinement_case(), neighbor_neighbor, 0)); @@ -2563,11 +2560,7 @@ neighbor_child_on_subface (const unsigned int face, #ifdef DEBUG // check, whether the face neighbor_child matches the requested // subface. - // - // in the following line, we know that if we got here that we - // have spacedim>=3. avoid instantiating invalid classes - // otherwise, even though we know that it is in dead code - typename Triangulation<3,(spacedim<3 ? 3 : spacedim)>::face_iterator requested; + typename Triangulation::face_iterator requested; switch (this->subface_case(face)) { case internal::SubfaceCase<3>::case_x: