From: wolf Date: Mon, 29 Sep 2003 15:18:11 +0000 (+0000) Subject: Fix a number of problems in code that didn't recognize face orientation issues yet. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f343fdc7646998802a25ebaef7f54b74df5f57ad;p=dealii-svn.git Fix a number of problems in code that didn't recognize face orientation issues yet. git-svn-id: https://svn.dealii.org/trunk@8053 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/tria_accessor.templates.h b/deal.II/deal.II/include/grid/tria_accessor.templates.h index 1a436935e9..98e5833172 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.templates.h +++ b/deal.II/deal.II/include/grid/tria_accessor.templates.h @@ -507,33 +507,36 @@ TriaObjectAccessor<3,dim>::line (const unsigned int i) const Assert (this->used(), typename TriaAccessor::ExcCellNotUsed()); Assert (i<12, ExcIndexRange (i,0,12)); - // egcs 1.1.2 gets into trouble if we - // omit the this-> here, if one tries to - // inline this function into another - // function where the variable name - // quad is also used.. It then complains - // that the name look-up for - // for-loop-variables has changed. using - // this-> as here works around the problem - if (i<4) - return this->quad(0)->line(i); - else - if (i<8) - return this->quad(1)->line(i-4); - else - switch (i) - { - case 8: - return this->quad(2)->line(3); - case 9: - return this->quad(2)->line(1); - case 10: - return this->quad(4)->line(1); - case 11: - return this->quad(4)->line(3); - }; - Assert (false, ExcIndexRange(i,0,12)); - return TriaIterator >(this->tria, -1, -1, 0); + // get the line index by asking the + // quads. make sure we handle + // reverted faces correctly + // + // so set up a table that for each + // line describes a) from which + // quad to take it, b) which line + // therein it is if the face is + // oriented correctly, and c) if in + // the opposite direction + static const unsigned int lookup_table[12][3] = + { { 0, 0, 3 }, // take first four lines from front face + { 0, 1, 2 }, + { 0, 2, 1 }, + { 0, 3, 0 }, + + { 1, 0, 3 }, // second four lines from back face + { 1, 1, 2 }, + { 1, 2, 1 }, + { 1, 3, 0 }, + + { 2, 3, 0 }, // the rest randomly + { 2, 1, 2 }, + { 4, 1, 2 }, + { 4, 3, 0 }}; + + return (this->quad(lookup_table[i][0]) + ->line(get_face_orientation(lookup_table[i][0]) ? + lookup_table[i][1] : + lookup_table[i][2])); } @@ -562,25 +565,36 @@ TriaObjectAccessor<3,dim>::line_index (const unsigned int i) const { Assert (i<12, ExcIndexRange(i,0,12)); - if (i<4) - return quad(0)->line_index(i); - else - if (i<8) - return quad(1)->line_index(i-4); - else - switch (i) - { - case 8: - return quad(2)->line_index(3); - case 9: - return quad(2)->line_index(1); - case 10: - return quad(4)->line_index(1); - case 11: - return quad(4)->line_index(3); - }; - Assert (false, ExcIndexRange(i,0,12)); - return 0; + // get the line index by asking the + // quads. make sure we handle + // reverted faces correctly + // + // so set up a table that for each + // line describes a) from which + // quad to take it, b) which line + // therein it is if the face is + // oriented correctly, and c) if in + // the opposite direction + static const unsigned int lookup_table[12][3] = + { { 0, 0, 3 }, // take first four lines from front face + { 0, 1, 2 }, + { 0, 2, 1 }, + { 0, 3, 0 }, + + { 1, 0, 3 }, // second four lines from back face + { 1, 1, 2 }, + { 1, 2, 1 }, + { 1, 3, 0 }, + + { 2, 3, 0 }, // the rest randomly + { 2, 1, 2 }, + { 4, 1, 2 }, + { 4, 3, 0 }}; + + return (this->quad(lookup_table[i][0]) + ->line_index(get_face_orientation(lookup_table[i][0]) ? + lookup_table[i][1] : + lookup_table[i][2])); } diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 53aad9c6f3..eaa9e49dfa 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -5679,7 +5679,7 @@ Triangulation<3>::execute_refinement () // it might be that the // quad itself is not // at the boundary, but - // that one of its line + // that one of its lines // actually is. in this // case, the newly // created vertices at @@ -5958,7 +5958,11 @@ Triangulation<3>::execute_refinement () vertices[next_unused_vertex] += hex->line(line)->child(0)->vertex(1) * 7./192.; // finally add centers of - // faces + // faces. note that + // vertex 2 is an + // invariant with respect + // to the face + // orientation for (unsigned int face=0; face::faces_per_cell; ++face) vertices[next_unused_vertex] += hex->face(face)->child(0)->vertex(2) * @@ -5967,7 +5971,7 @@ Triangulation<3>::execute_refinement () // now that we created // the right point, make // up the six lines - // interior to the quad + // interior to the hex // (++ takes care of the // end of the vector) raw_line_iterator new_lines[6]; @@ -5979,7 +5983,8 @@ Triangulation<3>::execute_refinement () new_lines[i] = next_unused_line; ++next_unused_line; - Assert (new_lines[i]->used() == false, ExcCellShouldBeUnused()); + Assert (new_lines[i]->used() == false, + ExcCellShouldBeUnused()); }; // set the data of the @@ -6004,6 +6009,14 @@ Triangulation<3>::execute_refinement () // *3-*-1* *3-*-1* // / 0 / | 4 | // *--*--* *--*--* + // + // note that both asking + // for child 0 and for + // vertex 2 within that + // is invariant with + // respect to the face + // orientation, so we do + // not have to ask here const unsigned int vertex_indices[7] = { hex->face(0)->child(0)->vertex_index(2), hex->face(3)->child(0)->vertex_index(2), @@ -6102,46 +6115,85 @@ Triangulation<3>::execute_refinement () // | * | * // |/ |/ // * * - + // + // this time we have to + // take into account + // whether the different + // faces are oriented + // correctly or in the + // opposite direction, so + // store that up front + const bool face_orientation[6] + = { hex->get_face_orientation (0), + hex->get_face_orientation (1), + hex->get_face_orientation (2), + hex->get_face_orientation (3), + hex->get_face_orientation (4), + hex->get_face_orientation (5) }; + const unsigned int line_indices[30] = { - hex->face(0)->child(0)->line_index(1), //0 - hex->face(0)->child(1)->line_index(2), //1 - hex->face(0)->child(2)->line_index(3), //2 - hex->face(0)->child(3)->line_index(0), //3 - - hex->face(2)->child(0)->line_index(2), //4 - hex->face(2)->child(1)->line_index(2), //5 - hex->face(3)->child(0)->line_index(1), //6 - hex->face(3)->child(3)->line_index(1), //7 - - hex->face(4)->child(1)->line_index(2), //8 - hex->face(4)->child(0)->line_index(2), //9 - hex->face(5)->child(3)->line_index(1), //10 - hex->face(5)->child(0)->line_index(1), //11 + hex->face(0)->child(0 ) + ->line_index(face_orientation[0] ? 1 : 2), //0 + hex->face(0)->child(face_orientation[0] ? 1 : 3) + ->line_index(face_orientation[0] ? 2 : 1), //1 + hex->face(0)->child(2 ) + ->line_index(face_orientation[0] ? 3 : 0), //2 + hex->face(0)->child(face_orientation[0] ? 3 : 1) + ->line_index(face_orientation[0] ? 0 : 3), //3 + + hex->face(2)->child(0 ) + ->line_index(face_orientation[2] ? 2 : 1), //4 + hex->face(2)->child(face_orientation[2] ? 1 : 3) + ->line_index(face_orientation[2] ? 2 : 1), //5 + hex->face(3)->child(0 ) + ->line_index(face_orientation[3] ? 1 : 2), //6 + hex->face(3)->child(face_orientation[3] ? 3 : 1) + ->line_index(face_orientation[3] ? 1 : 2), //7 + + hex->face(4)->child(face_orientation[4] ? 1 : 3) + ->line_index(face_orientation[4] ? 2 : 1), //8 + hex->face(4)->child(0 ) + ->line_index(face_orientation[4] ? 2 : 1), //9 + hex->face(5)->child(face_orientation[5] ? 3 : 1) + ->line_index(face_orientation[5] ? 1 : 2), //10 + hex->face(5)->child(0 ) + ->line_index(face_orientation[5] ? 1 : 2), //11 new_lines[4]->index(), //12 new_lines[1]->index(), //13 new_lines[5]->index(), //14 new_lines[3]->index(), //15 - hex->face(1)->child(0)->line_index(1), //16 - hex->face(1)->child(1)->line_index(2), //17 - hex->face(1)->child(2)->line_index(3), //18 - hex->face(1)->child(3)->line_index(0), //19 - - hex->face(5)->child(0)->line_index(2), //20 - hex->face(5)->child(1)->line_index(2), //21 - hex->face(4)->child(0)->line_index(1), //22 - hex->face(4)->child(3)->line_index(1), //23 + hex->face(1)->child(0 ) + ->line_index(face_orientation[1] ? 1 : 2), //16 + hex->face(1)->child(face_orientation[1] ? 1 : 3) + ->line_index(face_orientation[1] ? 2 : 1), //17 + hex->face(1)->child(2 ) + ->line_index(face_orientation[1] ? 3 : 0), //18 + hex->face(1)->child(face_orientation[1] ? 3 : 1) + ->line_index(face_orientation[1] ? 0 : 3), //19 + + hex->face(5)->child(0 ) + ->line_index(face_orientation[5] ? 2 : 1), //20 + hex->face(5)->child(face_orientation[5] ? 1 : 3) + ->line_index(face_orientation[5] ? 2 : 1), //21 + hex->face(4)->child(0 ) + ->line_index(face_orientation[4] ? 1 : 2), //22 + hex->face(4)->child(face_orientation[4] ? 3 : 1) + ->line_index(face_orientation[4] ? 1 : 2), //23 new_lines[0]->index(), //24 new_lines[2]->index(), //25 - hex->face(2)->child(0)->line_index(1), //26 - hex->face(2)->child(3)->line_index(1), //27 - - hex->face(3)->child(0)->line_index(2), //28 - hex->face(3)->child(1)->line_index(2) //29 + hex->face(2)->child(0 ) + ->line_index(face_orientation[2] ? 1 : 2), //26 + hex->face(2)->child(face_orientation[2] ? 3 : 1) + ->line_index(face_orientation[2] ? 1 : 2), //27 + + hex->face(3)->child(0 ) + ->line_index(face_orientation[3] ? 2 : 1), //28 + hex->face(3)->child(face_orientation[3] ? 1 : 3) + ->line_index(face_orientation[3] ? 2 : 1) //29 }; // find some space for @@ -6156,7 +6208,8 @@ Triangulation<3>::execute_refinement () new_quads[i] = next_unused_quad; ++next_unused_quad; - Assert (new_quads[i]->used() == false, ExcCellShouldBeUnused()); + Assert (new_quads[i]->used() == false, + ExcCellShouldBeUnused()); }; // set up the 12 quads, @@ -6309,6 +6362,17 @@ Triangulation<3>::execute_refinement () // 32/ / | |24 // |/20 21/ |12 13|/ // *-------* *-------* + // + // note that we have to + // take care of the + // orientation of + // faces. as an + // optimization: asking + // for child 0 or 2 of a + // face is invariant + // under the orientation, + // so we don't have to + // ask for it then const unsigned int quad_indices[36] = { new_quads[0]->index(), //0 @@ -6325,34 +6389,34 @@ Triangulation<3>::execute_refinement () new_quads[11]->index(), //11 hex->face(0)->child_index(0), //12 - hex->face(0)->child_index(1), + hex->face(0)->child_index(face_orientation[0] ? 1 : 3), hex->face(0)->child_index(2), - hex->face(0)->child_index(3), + hex->face(0)->child_index(face_orientation[0] ? 3 : 1), hex->face(1)->child_index(0), //16 - hex->face(1)->child_index(1), + hex->face(1)->child_index(face_orientation[1] ? 1 : 3), hex->face(1)->child_index(2), - hex->face(1)->child_index(3), + hex->face(1)->child_index(face_orientation[1] ? 3 : 1), hex->face(2)->child_index(0), //20 - hex->face(2)->child_index(1), + hex->face(2)->child_index(face_orientation[2] ? 1 : 3), hex->face(2)->child_index(2), - hex->face(2)->child_index(3), + hex->face(2)->child_index(face_orientation[2] ? 3 : 1), hex->face(3)->child_index(0), //24 - hex->face(3)->child_index(1), + hex->face(3)->child_index(face_orientation[3] ? 1 : 3), hex->face(3)->child_index(2), - hex->face(3)->child_index(3), + hex->face(3)->child_index(face_orientation[3] ? 3 : 1), hex->face(4)->child_index(0), //28 - hex->face(4)->child_index(1), + hex->face(4)->child_index(face_orientation[4] ? 1 : 3), hex->face(4)->child_index(2), - hex->face(4)->child_index(3), + hex->face(4)->child_index(face_orientation[4] ? 3 : 1), hex->face(5)->child_index(0), //32 - hex->face(5)->child_index(1), + hex->face(5)->child_index(face_orientation[5] ? 1 : 3), hex->face(5)->child_index(2), - hex->face(5)->child_index(3) + hex->face(5)->child_index(face_orientation[5] ? 3 : 1) }; @@ -6371,7 +6435,8 @@ Triangulation<3>::execute_refinement () for (unsigned int i=0; i<8; ++i) { new_hexes[i] = next_unused_hex; - Assert (new_hexes[i]->used() == false, ExcCellShouldBeUnused()); + Assert (new_hexes[i]->used() == false, + ExcCellShouldBeUnused()); ++next_unused_hex; };