]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix a bogus assertion.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Nov 2004 15:29:53 +0000 (15:29 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Nov 2004 15:29:53 +0000 (15:29 +0000)
git-svn-id: https://svn.dealii.org/trunk@9759 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/grid/tria_accessor.cc
deal.II/doc/news/c-5.0.html
tests/bits/mesh_3d_17.cc
tests/results/i686-pc-linux-gnu+gcc3.2/bits/mesh_3d_17.output

index d4bf89608942c03e5ca047e528620c68ea5123cd..caed8e0faf3d23a036930a008888549ab123ba06 100644 (file)
@@ -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;
 }
 
index cb1c4832ee2095c9a74088f537c27e9d8e3d6c3d..61638c1d971112f4ae81703cf91de35ac3bfcb70 100644 (file)
@@ -389,6 +389,17 @@ inconvenience this causes.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       Fixed: In rare cases of 3d meshes with a certain topology, we triggered
+       an assertion in
+       <code>TriaAccessor::neighbor_child_on_subface</code>. 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.
+       <br>
+       (WB 2004/11/12)
+       </p>
+
   <li> <p>
        Fixed: The <code>GridGenerator::cylinder</code> function in 3d
        assigned the wrong boundary value to the top and bottom part of
index 80fdbe6945651d7809efc25a18e9900dfdd7f42d..6df508ddff4a696dd430a8c2aa346a5944e4e906 100644 (file)
 #include <grid/grid_reordering.h>
 #include <grid/grid_in.h>
 #include <grid/grid_out.h>
-#include <dofs/dof_handler.h>
-#include <dofs/dof_accessor.h>
-#include <fe/fe_q.h>
 
 #include <fstream>
 #include <iostream>
 
-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<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
+  Triangulation<dim>::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<GeometryInfo<dim>::faces_per_cell; ++face_no)
+    for (unsigned int face_no=0;
+         face_no<GeometryInfo<dim>::faces_per_cell;
+         ++face_no)
       {
-       const DoFHandler<dim>::face_iterator
+       const Triangulation<dim>::face_iterator
          face=cell->face(face_no);
-           
-       if (!face->at_boundary())
-         {
-           Assert (cell->neighbor(face_no).state() == IteratorState::valid,
-                   ExcInternalError());
-           DoFHandler<dim>::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<GeometryInfo<dim>::subfaces_per_face;
-                    ++subface_no)
-                 {
-                   deallog << "subface_no=" << subface_no << endl;
-                   DoFHandler<dim>::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<dim>::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<dim>::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<GeometryInfo<dim>::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();
   
index 45152045f862aa30b6514335e4748233a3fd624b..aa1a5e014b445b5c2fb9f10a173ac6d4abac53ba 100644 (file)
@@ -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

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.