From dbde32d6b40877d03460b032a575873d6056195c Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 15 Nov 2004 15:29:53 +0000 Subject: [PATCH] Fix a bogus assertion. git-svn-id: https://svn.dealii.org/trunk@9759 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/grid/tria_accessor.cc | 44 +++++- deal.II/doc/news/c-5.0.html | 11 ++ tests/bits/mesh_3d_17.cc | 147 ++++++++++-------- .../bits/mesh_3d_17.output | 28 +++- 4 files changed, 156 insertions(+), 74 deletions(-) diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc index d4bf896089..caed8e0faf 100644 --- a/deal.II/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/deal.II/source/grid/tria_accessor.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -2200,12 +2200,46 @@ neighbor_child_on_subface (const unsigned int face, (face_orientations_match ? subface : subface_translation[subface]))); - TriaIterator<3,CellAccessor<3> > neighbor_child= + const TriaIterator<3,CellAccessor<3> > neighbor_child= this->neighbor(face)->child(neighbor_child_index); - Assert(this->face(face)->child(subface)== - neighbor_child->face(neighbor_neighbor), ExcInternalError()); - + // if the face on the side of the present + // cell is in the correct order, then make + // sure that the neighbor child cell we + // have found shares the desired subface. + // + // otherwise: if the face is turned + // inside out when viewed from the + // present cell, then the subface we are + // interested in is not + // this->face(face)->child(subface), but + // instead + // this->face(face)->child(subface'), + // where subface' is the translated + // subface number. this is so because we + // only store the face only once, so its + // children are defined in terms of its + // own (circular) orientation, not in + // terms of face_orientation as viewed + // from one of the adjacent cells. in + // that case, cell->face(f)->subface(sf) + // may yield unexpected results; in that + // case, a caller may need to adjust + // according to the face_orientation + // flag, though one in general only wants + // to loop over all subfaces, and not + // pick a particular one + Assert(((this->face_orientation(face) == true) + && + (this->face(face)->child(subface) == + neighbor_child->face(neighbor_neighbor))) + || + ((this->face_orientation(face) == false) + && + (this->face(face)->child(subface_translation[subface]) == + neighbor_child->face(neighbor_neighbor))), + ExcInternalError()); + return neighbor_child; } diff --git a/deal.II/doc/news/c-5.0.html b/deal.II/doc/news/c-5.0.html index cb1c4832ee..61638c1d97 100644 --- a/deal.II/doc/news/c-5.0.html +++ b/deal.II/doc/news/c-5.0.html @@ -389,6 +389,17 @@ inconvenience this causes.

deal.II

    +

  1. + Fixed: In rare cases of 3d meshes with a certain topology, we triggered + an assertion in + TriaAccessor::neighbor_child_on_subface. It turns out that + the code actually computes the correct answer, but the assertion had a + condition that doesn't always have to be satisfied. This bogus + assertion is now fixed. +
    + (WB 2004/11/12) +

    +

  2. Fixed: The GridGenerator::cylinder function in 3d assigned the wrong boundary value to the top and bottom part of diff --git a/tests/bits/mesh_3d_17.cc b/tests/bits/mesh_3d_17.cc index 80fdbe6945..6df508ddff 100644 --- a/tests/bits/mesh_3d_17.cc +++ b/tests/bits/mesh_3d_17.cc @@ -9,91 +9,116 @@ #include #include #include -#include -#include -#include #include #include -using namespace std; -void check_this (Triangulation<3> &tria) + +void check (Triangulation<3> &tria) { - FE_Q<3> fe(1); - - DoFHandler<3> dof_handler (tria); - dof_handler.distribute_dofs (fe); + deallog << "Coarse cell 0 vertices:" << std::endl; + for (unsigned int i=0; i<8; ++i) + deallog << ' ' << tria.begin_active()->vertex_index(i); + deallog << std::endl; + + deallog << "Coarse cell 1 vertices:" << std::endl; + for (unsigned int i=0; i<8; ++i) + deallog << ' ' << (++tria.begin_active())->vertex_index(i); + deallog << std::endl; + + tria.begin_active()->set_refine_flag (); + tria.execute_coarsening_and_refinement (); + const unsigned int dim=3; - DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + Triangulation::active_cell_iterator + cell = tria.begin_active(), + endc = tria.end(); + + // loop over all interior faces of active + // cells with children. there should be + // exactly one, namely the middle face of + // the big cell (the middle faces of the + // active cells on the refined side don't + // have any children) for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + for (unsigned int face_no=0; + face_no::faces_per_cell; + ++face_no) { - const DoFHandler::face_iterator + const Triangulation::face_iterator face=cell->face(face_no); - - if (!face->at_boundary()) - { - Assert (cell->neighbor(face_no).state() == IteratorState::valid, - ExcInternalError()); - DoFHandler::cell_iterator neighbor= - cell->neighbor(face_no); - - if (face->has_children()) - { - const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); - - deallog << "cell_no=" << cell_no << ", face_no=" << face_no << endl; - deallog << " cell->face_orientation(face_no)=" - << cell->face_orientation(face_no) << endl; - deallog << " neighbor2=" << neighbor2 << endl; - deallog << " neighbor->face_orientation(neighbor2)=" - << neighbor->face_orientation(neighbor2) << endl; - - - for (unsigned int subface_no=0; - subface_no::subfaces_per_face; - ++subface_no) - { - deallog << "subface_no=" << subface_no << endl; - DoFHandler::cell_iterator neighbor_child= - cell->neighbor_child_on_subface(face_no, subface_no); - } - } - } - } -} + if (!face->at_boundary() && face->has_children()) + { + deallog << "Interior face with children:" << std::endl; + deallog << " cell=" << cell << std::endl; + deallog << " face_no=" << face_no << std::endl; + deallog << " vertices=" + << face->vertex_index(0) << ' ' + << face->vertex_index(1) << ' ' + << face->vertex_index(2) << ' ' + << face->vertex_index(3) << std::endl; + deallog << " face orientation=" + << (cell->face_orientation(face_no) ? "true" : "false") + << std::endl; + + const Triangulation::cell_iterator neighbor + = cell->neighbor(face_no); + const unsigned int nb_of_nb + = cell->neighbor_of_neighbor(face_no); -void check (Triangulation<3> &tria) -{ - tria.begin_active()->set_refine_flag (); - tria.execute_coarsening_and_refinement (); - - if (false) - { - GridOut grid_out; - ofstream grid_stream("grid.gnuplot"); - grid_out.write_gnuplot(tria, grid_stream); - } - - check_this (tria); + // check that + // neighbor_of_neighbor works: + Assert (neighbor->neighbor(nb_of_nb) == cell, + ExcInternalError()); + + const Triangulation::face_iterator neighbor_face + = neighbor->face(nb_of_nb); + + deallog << " neighbor face vertices=" + << neighbor_face->vertex_index(0) << ' ' + << neighbor_face->vertex_index(1) << ' ' + << neighbor_face->vertex_index(2) << ' ' + << neighbor_face->vertex_index(3) + << std::endl; + deallog << " neighbor->face_orientation(nb_of_nb)=" + << (neighbor->face_orientation(nb_of_nb) ? + "true" : "false") + << std::endl; + + + deallog << " checking subfaces:" << std::endl; + for (unsigned int subface_no=0; + subface_no::subfaces_per_face; + ++subface_no) + { + deallog << " subface_no=" << subface_no << std::endl; + deallog << " center=" << face->child(subface_no)->center() + << std::endl; + + // we used to abort in the + // following call, due to a + // wrong Assert condition: + cell->neighbor_child_on_subface(face_no, subface_no); + } + } + } } + int main () { - ofstream logfile("mesh_3d_17.output"); + std::ofstream logfile("mesh_3d_17.output"); deallog.attach(logfile); deallog.depth_console(0); Triangulation<3> coarse_grid; GridIn<3> in; in.attach_triangulation(coarse_grid); - ifstream ucd_file("two_cubes.inp"); + std::ifstream ucd_file("two_cubes.inp"); in.read_ucd(ucd_file); ucd_file.close(); diff --git a/tests/results/i686-pc-linux-gnu+gcc3.2/bits/mesh_3d_17.output b/tests/results/i686-pc-linux-gnu+gcc3.2/bits/mesh_3d_17.output index 45152045f8..aa1a5e014b 100644 --- a/tests/results/i686-pc-linux-gnu+gcc3.2/bits/mesh_3d_17.output +++ b/tests/results/i686-pc-linux-gnu+gcc3.2/bits/mesh_3d_17.output @@ -1,9 +1,21 @@ -DEAL::cell_no=0, face_no=2 -DEAL:: cell->face_orientation(face_no)=0 -DEAL:: neighbor2=5 -DEAL:: neighbor->face_orientation(neighbor2)=1 -DEAL::subface_no=0 -DEAL::subface_no=1 -DEAL::subface_no=2 -DEAL::subface_no=3 +DEAL::Coarse cell 0 vertices: +DEAL:: 1 2 5 4 7 8 11 10 +DEAL::Coarse cell 1 vertices: +DEAL:: 1 4 3 0 7 10 9 6 +DEAL::Interior face with children: +DEAL:: cell=0.1 +DEAL:: face_no=2 +DEAL:: vertices=1 7 10 4 +DEAL:: face orientation=false +DEAL:: neighbor face vertices=1 7 10 4 +DEAL:: neighbor->face_orientation(nb_of_nb)=true +DEAL:: checking subfaces: +DEAL:: subface_no=0 +DEAL:: center=0.00000 0.250000 0.750000 +DEAL:: subface_no=1 +DEAL:: center=0.00000 0.250000 0.250000 +DEAL:: subface_no=2 +DEAL:: center=0.00000 0.750000 0.250000 +DEAL:: subface_no=3 +DEAL:: center=0.00000 0.750000 0.750000 -- 2.39.5