]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix a problem in GridTools::find_cells_adjacent_to_vertex with anisotropic refinement.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 8 Oct 2012 17:31:55 +0000 (17:31 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 8 Oct 2012 17:31:55 +0000 (17:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@27010 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/source/grid/grid_tools.cc
tests/aniso/anisotropic_crash.cc [moved from tests/fail/anisotropic_crash.cc with 82% similarity]
tests/aniso/anisotropic_crash/cmp/generic [new file with mode: 0644]
tests/fail/anisotropic_crash/cmp/generic [deleted file]

index d98241be0ed88b5dcea8400a63ef98972667478e..414c76aae946840d642d4de5f1e8dde6be1e77a0 100644 (file)
@@ -87,6 +87,11 @@ DoFHandler, in particular removal of specializations.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: GridTools::find_cells_adjacent_to_vertex got into
+trouble with anisotropically refined meshes. This is now fixed.
+<br>
+(Abner Salgado, Tobias Leicht, Wolfgang Bangerth, 2012/10/08)
+
 <li> Fixed: FESystem can now deal with n_elements==0 for a block.
 <br>
 (Timo Heister, 2012/09/28)
index 750aad1dcd71e0feef5f4538e16a444712d2b1d8..668d238b18bb0abe87a8e6736d7ad923fc428230 100644 (file)
@@ -679,7 +679,7 @@ namespace GridTools
       cell = container.begin_active(),
       endc = container.end();
 
-    // go through all active cells and look. if the vertex is part of that cell
+    // go through all active cells and look if the vertex is part of that cell
     //
     // in 1d, this is all we need to care about. in 2d/3d we also need to worry
     // that the vertex might be a hanging node on a face or edge of a cell; in
@@ -701,6 +701,17 @@ namespace GridTools
     // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we
     // haven't yet found the vertex for the current cell we also need to
     // look at the mid-points of edges
+    //
+    // as a final note, deciding whether a neighbor is actually coarses is
+    // simple in the case of isotropic refinement (we just need to look at
+    // the level of the current and the neighboring cell). however, this
+    // isn't so simple if we have used anisotropic refinement since then
+    // the level of a cell is not indicative of whether it is coarser or
+    // not than the current cell. ultimately, we want to add all cells on
+    // which the vertex is, independent of whether they are coarser or
+    // finer and so in the 2d case below we simply add *any* *active* neighbor.
+    // in the worst case, we add cells multiple times to the adjacent_cells
+    // list, but std::set throws out those cells already entered
     for (; cell != endc; ++cell)
       {
         for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; v++)
@@ -712,9 +723,9 @@ namespace GridTools
               adjacent_cells.insert(cell);
 
               // as explained above, in 2+d we need to check whether
-              // this vertex is on a face behind which there is a coarser
-              // neighbor. if this is the case, then we need to also
-              // add this neighbor
+              // this vertex is on a face behind which there is a
+              // (possibly) coarser neighbor. if this is the case,
+              // then we need to also add this neighbor
               if (dim >= 2)
                 for (unsigned int vface = 0; vface < dim; vface++)
                   {
@@ -723,16 +734,17 @@ namespace GridTools
 
                     if (!cell->at_boundary(face)
                         &&
-                        (cell->neighbor(face)->level() < cell->level()))
+                        cell->neighbor(face)->active())
                       {
-                        // there is a coarser cell behind a face to which the
-                        // vertex belongs. the vertex we are looking at is
-                        // then either a vertex of that coarser neighbor,
-                        // or it is a hanging node on one of the faces of
-                        // that cell. in either case, it is adjacent
-                        // to the vertex, so add it to the list as well
-                        // (if the cell was already in the list then
-                        // the std::set makes sure that we get it only
+                        // there is a (possibly )coarser cell behind a
+                        // face to which the vertex belongs. the
+                        // vertex we are looking at is then either a
+                        // vertex of that coarser neighbor, or it is a
+                        // hanging node on one of the faces of that
+                        // cell. in either case, it is adjacent to the
+                        // vertex, so add it to the list as well (if
+                        // the cell was already in the list then the
+                        // std::set makes sure that we get it only
                         // once)
                         adjacent_cells.insert (cell->neighbor(face));
                         goto next_cell;
@@ -786,7 +798,7 @@ namespace GridTools
   {
        template <int dim, template<int, int> class Container, int spacedim>
        void find_active_cell_around_point_internal(const Container<dim,spacedim>& container,
-         std::set<typename Container<dim,spacedim>::active_cell_iterator>& searched_cells, 
+         std::set<typename Container<dim,spacedim>::active_cell_iterator>& searched_cells,
          std::set<typename Container<dim,spacedim>::active_cell_iterator>& adjacent_cells)
        {
          typedef typename Container<dim,spacedim>::active_cell_iterator cell_iterator;
@@ -794,10 +806,10 @@ namespace GridTools
                                        // update the searched cells
          searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end());
                                        // now we to collect all neighbors
-                                       // of the cells in adjacent_cells we 
+                                       // of the cells in adjacent_cells we
                                        // have not yet searched.
          std::set<cell_iterator> adjacent_cells_new;
-               
+
          typename std::set<cell_iterator>::const_iterator
          cell = adjacent_cells.begin(),
          endc = adjacent_cells.end();
@@ -826,7 +838,7 @@ namespace GridTools
                        adjacent_cells.insert(it);
                         break;
                  }
-       }         
+       }
        }
   }
 
@@ -862,7 +874,7 @@ namespace GridTools
 
     std::vector<cell_iterator> adjacent_cells_tmp =
       find_cells_adjacent_to_vertex(container, vertex);
-         
+
                                                                        // Make sure that we have found
                                                                        // at least one cell adjacent to vertex.
          Assert(adjacent_cells_tmp.size()>0, ExcInternalError());
@@ -875,7 +887,7 @@ namespace GridTools
                                                                        // in the grid.
                                                                        // As long as we have not found
                                                                        // the cell and have not searched
-                                                                       // every cell in the triangulation, 
+                                                                       // every cell in the triangulation,
                                                                        // we keep on looking.
          const unsigned int n_cells =get_tria(container).n_cells();
          bool found = false;
@@ -925,9 +937,9 @@ namespace GridTools
                  }
                                                        //udpate the number of cells searched
                  cells_searched += adjacent_cells.size();
-                                                       // if we have not found the cell in 
+                                                       // if we have not found the cell in
                                                        // question and have not yet searched every
-                                                       // cell, we expand our search to 
+                                                       // cell, we expand our search to
                                                        // all the not already searched neighbors of
                                                        // the cells in adjacent_cells. This is
                                                        // what find_active_cell_around_poin_internal
@@ -984,7 +996,7 @@ namespace GridTools
 
                std::vector<cell_iterator> adjacent_cells_tmp =
         find_cells_adjacent_to_vertex(container, vertex);
-               
+
                                                                        // Make sure that we have found
                                                                        // at least one cell adjacent to vertex.
          Assert(adjacent_cells_tmp.size()>0, ExcInternalError());
@@ -997,7 +1009,7 @@ namespace GridTools
                                                                        // in the grid.
                                                                        // As long as we have not found
                                                                        // the cell and have not searched
-                                                                       // every cell in the triangulation, 
+                                                                       // every cell in the triangulation,
                                                                        // we keep on looking.
          const unsigned int n_cells =get_tria(container).n_cells();
          bool found = false;
@@ -1013,7 +1025,7 @@ namespace GridTools
                          {
                                const Point<dim> p_cell = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p);
 
-                               
+
                                // calculate the infinity norm of
                                // the distance vector to the unit cell.
                                const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
@@ -1048,9 +1060,9 @@ namespace GridTools
                  }
                                                                                //udpate the number of cells searched
                  cells_searched += adjacent_cells.size();
-                                                       // if we have not found the cell in 
+                                                       // if we have not found the cell in
                                                        // question and have not yet searched every
-                                                       // cell, we expand our search to 
+                                                       // cell, we expand our search to
                                                        // all the not already searched neighbors of
                                                        // the cells in adjacent_cells.
                  if(!found && cells_searched < n_cells)
similarity index 82%
rename from tests/fail/anisotropic_crash.cc
rename to tests/aniso/anisotropic_crash.cc
index 70d332fbbed1d690963912acfe8296f8ad80ef76..adced6408947b675a5820602697df28cedb201be 100644 (file)
@@ -12,8 +12,8 @@
 //----------------------------  anisotropic_crash.cc  ---------------------------
 
 
-// Trying to catch a bug in the construction of patches when
-// there is anisotropic refinement
+// GridTools::find_cells_adjacent_to_vertex had a problem in that it
+// wasn't prepared to deal with anisotropic refinement
 
 
 #include <deal.II/grid/tria.h>
@@ -79,7 +79,15 @@ int main()
 
                                   /// For each vertex find the patch of cells
                                   /// that surrounds it
-  for( unsigned v=0; v<tri.n_used_vertices(); ++v )
+  for( unsigned v=0; v<tri.n_vertices(); ++v )
     if (tri.get_used_vertices()[v] == true)
-      GridTools::find_cells_adjacent_to_vertex( tri, v );
+      {
+       deallog << "Vertex=" << v << std::endl;
+
+       const std::vector<Triangulation<2>::active_cell_iterator>
+         tmp = GridTools::find_cells_adjacent_to_vertex( tri, v );
+
+       for (unsigned int i=0; i<tmp.size(); ++i)
+         deallog << "   " << tmp[i] << std::endl;
+      }
 }
diff --git a/tests/aniso/anisotropic_crash/cmp/generic b/tests/aniso/anisotropic_crash/cmp/generic
new file mode 100644 (file)
index 0000000..ba2c502
--- /dev/null
@@ -0,0 +1,209 @@
+
+DEAL::Vertex=0
+DEAL::   2.0
+DEAL::Vertex=1
+DEAL::   2.5
+DEAL::Vertex=2
+DEAL::   3.20
+DEAL::Vertex=3
+DEAL::   3.29
+DEAL::Vertex=4
+DEAL::   3.1
+DEAL::   3.6
+DEAL::Vertex=5
+DEAL::   3.5
+DEAL::   3.14
+DEAL::Vertex=6
+DEAL::   2.13
+DEAL::   3.13
+DEAL::Vertex=7
+DEAL::   2.11
+DEAL::   2.14
+DEAL::Vertex=8
+DEAL::   2.3
+DEAL::   3.10
+DEAL::   3.16
+DEAL::   3.22
+DEAL::Vertex=9
+DEAL::   2.0
+DEAL::   3.0
+DEAL::Vertex=10
+DEAL::   2.5
+DEAL::   3.6
+DEAL::Vertex=11
+DEAL::   2.0
+DEAL::   3.4
+DEAL::Vertex=12
+DEAL::   3.15
+DEAL::   3.18
+DEAL::Vertex=13
+DEAL::   2.5
+DEAL::   3.13
+DEAL::Vertex=14
+DEAL::   2.13
+DEAL::   3.27
+DEAL::Vertex=15
+DEAL::   2.11
+DEAL::   3.21
+DEAL::Vertex=16
+DEAL::   2.14
+DEAL::   3.28
+DEAL::Vertex=17
+DEAL::   2.3
+DEAL::   3.3
+DEAL::   3.7
+DEAL::   3.8
+DEAL::Vertex=18
+DEAL::   2.11
+DEAL::   2.14
+DEAL::   3.17
+DEAL::   3.24
+DEAL::Vertex=19
+DEAL::   2.3
+DEAL::   3.5
+DEAL::   3.14
+DEAL::   3.16
+DEAL::Vertex=20
+DEAL::   2.13
+DEAL::   3.11
+DEAL::   3.12
+DEAL::   3.23
+DEAL::Vertex=21
+DEAL::   2.0
+DEAL::   2.3
+DEAL::   3.2
+DEAL::   3.4
+DEAL::Vertex=22
+DEAL::   2.5
+DEAL::   3.7
+DEAL::   3.9
+DEAL::   3.12
+DEAL::Vertex=23
+DEAL::   2.11
+DEAL::   3.15
+DEAL::   3.17
+DEAL::   3.19
+DEAL::Vertex=24
+DEAL::   2.13
+DEAL::   2.14
+DEAL::   3.25
+DEAL::   3.26
+DEAL::Vertex=25
+DEAL::   3.0
+DEAL::   3.1
+DEAL::Vertex=26
+DEAL::   3.4
+DEAL::   3.5
+DEAL::Vertex=27
+DEAL::   3.14
+DEAL::   3.15
+DEAL::Vertex=28
+DEAL::   3.18
+DEAL::   3.20
+DEAL::Vertex=29
+DEAL::   3.27
+DEAL::   3.29
+DEAL::Vertex=30
+DEAL::   3.20
+DEAL::   3.21
+DEAL::Vertex=31
+DEAL::   3.28
+DEAL::   3.29
+DEAL::Vertex=32
+DEAL::   3.1
+DEAL::   3.3
+DEAL::   3.6
+DEAL::   3.7
+DEAL::Vertex=33
+DEAL::   2.3
+DEAL::   3.8
+DEAL::   3.10
+DEAL::Vertex=34
+DEAL::   3.16
+DEAL::   3.17
+DEAL::   3.22
+DEAL::   3.24
+DEAL::Vertex=35
+DEAL::   3.10
+DEAL::   3.11
+DEAL::   3.22
+DEAL::   3.23
+DEAL::Vertex=36
+DEAL::   3.12
+DEAL::   3.13
+DEAL::Vertex=37
+DEAL::   2.0
+DEAL::   3.0
+DEAL::   3.2
+DEAL::Vertex=38
+DEAL::   2.3
+DEAL::   3.4
+DEAL::   3.5
+DEAL::Vertex=39
+DEAL::   3.2
+DEAL::   3.3
+DEAL::Vertex=40
+DEAL::   2.5
+DEAL::   3.6
+DEAL::   3.7
+DEAL::Vertex=41
+DEAL::   3.9
+DEAL::   3.11
+DEAL::   3.12
+DEAL::Vertex=42
+DEAL::   3.8
+DEAL::   3.9
+DEAL::Vertex=43
+DEAL::   3.12
+DEAL::   3.13
+DEAL::Vertex=44
+DEAL::   3.14
+DEAL::   3.15
+DEAL::   3.16
+DEAL::   3.17
+DEAL::Vertex=45
+DEAL::   2.11
+DEAL::   3.19
+DEAL::   3.21
+DEAL::Vertex=46
+DEAL::   3.18
+DEAL::   3.19
+DEAL::Vertex=47
+DEAL::   2.13
+DEAL::   3.23
+DEAL::   3.25
+DEAL::Vertex=48
+DEAL::   2.14
+DEAL::   3.26
+DEAL::   3.28
+DEAL::Vertex=49
+DEAL::   3.24
+DEAL::   3.25
+DEAL::Vertex=50
+DEAL::   3.26
+DEAL::   3.27
+DEAL::Vertex=51
+DEAL::   3.0
+DEAL::   3.1
+DEAL::   3.2
+DEAL::   3.3
+DEAL::Vertex=52
+DEAL::   3.8
+DEAL::   3.9
+DEAL::   3.10
+DEAL::   3.11
+DEAL::Vertex=53
+DEAL::   3.18
+DEAL::   3.19
+DEAL::   3.20
+DEAL::   3.21
+DEAL::Vertex=54
+DEAL::   3.22
+DEAL::   3.23
+DEAL::   3.24
+DEAL::   3.25
+DEAL::Vertex=55
+DEAL::   3.26
+DEAL::   3.27
+DEAL::   3.28
+DEAL::   3.29
diff --git a/tests/fail/anisotropic_crash/cmp/generic b/tests/fail/anisotropic_crash/cmp/generic
deleted file mode 100644 (file)
index 9494292..0000000
+++ /dev/null
@@ -1,110 +0,0 @@
-
-DEAL::Checking Q1 in 1d:
-DEAL::Checking Q1 in 2d:
-DEAL::Checking Q1 in 3d:
-DEAL::Checking Q2 in 1d:
-DEAL::Checking Q2 in 2d:
-DEAL::Checking Q2 in 3d:
-DEAL::Checking Q3 in 1d:
-DEAL::Checking Q3 in 2d:
-DEAL::Checking Q3 in 3d:
-DEAL::Checking DGQ0 in 1d:
-DEAL::Checking DGQ0 in 2d:
-DEAL::Checking DGQ0 in 3d:
-DEAL::Checking DGQ1 in 1d:
-DEAL::Checking DGQ1 in 2d:
-DEAL::Checking DGQ1 in 3d:
-DEAL::Checking DGQ3 in 1d:
-DEAL::Checking DGQ3 in 2d:
-DEAL::Checking DGQ3 in 3d:
-DEAL::Checking DGP0 in 1d:
-DEAL::Checking DGP0 in 2d:
-DEAL::Checking DGP0 in 3d:
-DEAL::Checking DGP1 in 1d:
-DEAL::Checking DGP1 in 2d:
-DEAL::Checking DGP1 in 3d:
-DEAL::Checking DGP3 in 1d:
-DEAL::Checking DGP3 in 2d:
-DEAL::Checking DGP3 in 3d:
-DEAL::Checking Nedelec1 in 2d:
-DEAL::Checking Nedelec1 in 3d:
-DEAL::Checking RaviartThomas0 in 2d:
-DEAL::4
-DEAL::1
-DEAL:cg::Starting value 1.2
-DEAL:cg::Convergence step 10 value 0
-DEAL::L2_Error : 0
-DEAL::Checking RaviartThomas1 in 2d:
-DEAL::8
-DEAL::2
-DEAL:cg::Starting value 1.4
-DEAL:cg::Convergence step 12 value 0
-DEAL::L2_Error : 0
-DEAL::Checking RaviartThomas2 in 2d:
-DEAL::12
-DEAL::2
-DEAL:cg::Starting value 1.4
-DEAL:cg::Convergence step 14 value 0
-DEAL::L2_Error : 0
-DEAL::Checking RaviartThomasNodal0 in 2d:
-DEAL::4
-DEAL::1
-DEAL:cg::Starting value 0.49
-DEAL:cg::Convergence step 10 value 0
-DEAL::L2_Error : 0
-DEAL::Checking RaviartThomasNodal1 in 2d:
-DEAL::8
-DEAL::2
-DEAL:cg::Starting value 0.29
-DEAL:cg::Convergence step 12 value 0
-DEAL::L2_Error : 0
-DEAL::Checking RaviartThomasNodal2 in 2d:
-DEAL::12
-DEAL::3
-DEAL:cg::Starting value 0.23
-DEAL:cg::Convergence step 16 value 0
-DEAL::L2_Error : 0
-DEAL::Checking RaviartThomasNodal0 in 3d:
-DEAL::16
-DEAL::1
-DEAL:cg::Starting value 0.47
-DEAL:cg::Convergence step 9 value 0
-DEAL::L2_Error : 0
-DEAL::Checking RaviartThomasNodal1 in 3d:
-DEAL::64
-DEAL::4
-DEAL:cg::Starting value 0.20
-DEAL:cg::Convergence step 12 value 0
-DEAL::L2_Error : 0
-DEAL::Checking RaviartThomasNodal2 in 3d:
-DEAL::144
-DEAL::9
-DEAL:cg::Starting value 0.13
-DEAL:cg::Convergence step 15 value 0
-DEAL::L2_Error : 0
-DEAL::Checking FE_Q<1>(1)3 in 1d:
-DEAL::Checking FE_DGQ<1>(2)2 in 1d:
-DEAL::Checking FE_DGP<1>(3)1 in 1d:
-DEAL::Checking FE_Q<2>(1)3 in 2d:
-DEAL::Checking FE_DGQ<2>(2)2 in 2d:
-DEAL::Checking FE_DGP<2>(3)1 in 2d:
-DEAL::Checking FE_Q<3>(1)3 in 3d:
-DEAL::Checking FE_DGQ<3>(2)2 in 3d:
-DEAL::Checking FE_DGP<3>(3)1 in 3d:
-DEAL::Checking FE_Q<1>(1)3FE_DGQ<1>(2)2 in 1d:
-DEAL::Checking FE_DGQ<1>(2)2FE_DGP<1>(3)1 in 1d:
-DEAL::Checking FE_DGP<1>(3)1FE_DGQ<1>(2)2 in 1d:
-DEAL::Checking FE_Q<2>(1)3FE_DGQ<2>(2)2 in 2d:
-DEAL::Checking FE_DGQ<2>(2)2FE_DGP<2>(3)1 in 2d:
-DEAL::Checking FE_DGP<2>(3)1FE_DGQ<2>(2)2 in 2d:
-DEAL::Checking FE_Q<1>(1)3FE_DGP<1>(3)1FE_Q<1>(1)3 in 1d:
-DEAL::Checking FE_DGQ<1>(2)2FE_DGQ<1>(2)2FE_Q<1>(3)3 in 1d:
-DEAL::Checking FE_DGP<1>(3)1FE_DGP<1>(3)1FE_Q<1>(2)3 in 1d:
-DEAL::Checking FE_Q<2>(1)3FE_DGP<2>(3)1FE_Q<2>(1)3 in 2d:
-DEAL::Checking FE_DGQ<2>(2)2FE_DGQ<2>(2)2FE_Q<2>(3)3 in 2d:
-DEAL::Checking FE_DGP<2>(3)1FE_DGP<2>(3)1FE_Q<2>(2)3 in 2d:
-DEAL::Checking FE_DGQ<3>(1)3FE_DGP<3>(3)1FE_Q<3>(1)3 in 3d:
-DEAL::Checking (FESystem<2>(FE_Q<2>(1),3))3FE_DGQ<2>(0)1FE_Q<2>(1)3 in 2d:
-DEAL::Checking FE_DGQ<2>(3)1FESystem<2>(FE_DGQ<2>(0),3)1FESystem<2>(FE_Q<2>(2),1, FE_DGQ<2>(0),1)2 in 2d:
-DEAL::Checking FE_DGQ<2>(3)1FE_Nedelec<2>(1)2 in 2d:
-DEAL::Checking FE_Nedelec<2>(1)1FESystem<2>(FE_DGQ<2>(1),2)1FESystem<2>(FE_Q<2>(2),1, FE_Nedelec<2>(1),2)2 in 2d:

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.