]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix a problem where the update_neighbors function in the Triangulation class needed...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 3 Jan 2011 18:10:50 +0000 (18:10 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 3 Jan 2011 18:10:50 +0000 (18:10 +0000)
git-svn-id: https://svn.dealii.org/trunk@23115 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/grid/tria.cc
tests/codim_one/hanging_nodes_01.cc [new file with mode: 0644]
tests/codim_one/hanging_nodes_01/cmp/generic [new file with mode: 0644]
tests/codim_one/hanging_nodes_02.cc [new file with mode: 0644]
tests/codim_one/hanging_nodes_02/cmp/generic [new file with mode: 0644]
tests/codim_one/hanging_nodes_03.cc [new file with mode: 0644]
tests/codim_one/hanging_nodes_03/cmp/generic [new file with mode: 0644]

index 5a43056e79961f96a570d097d431a6785a702aeb..177f8aa2cc92cf4edc844fa205d583fa47c57e0d 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -1051,8 +1051,8 @@ namespace
                                     // intrinsic normal we define the left
                                     // neighbor as the one for which the face
                                     // normal points outward, and store that
-                                    // one first, with an offset of one we
-                                    // store the right neighbor for which the
+                                    // one first; the second one is then
+                                    // the right neighbor for which the
                                     // face normal points inward. This
                                     // information depends on the type of cell
                                     // and local number of face for the
@@ -1073,6 +1073,27 @@ namespace
                                     // line points to the right if the line
                                     // points upwards.
                                     //
+                                    // There is one more point to
+                                    // consider, however: if we have
+                                    // dim<spacedim, then we may have
+                                    // cases where cells are
+                                    // inverted. In effect, both
+                                    // cells think they are the left
+                                    // neighbor of an edge, for
+                                    // example, which leads us to
+                                    // forget neighborship
+                                    // information (a case that shows
+                                    // this is
+                                    // codim_one/hanging_nodes_02). We
+                                    // store whether a cell is
+                                    // inverted using the
+                                    // direction_flag, so if a cell
+                                    // has a false direction_flag,
+                                    // then we need to invert our
+                                    // selection whether we are a
+                                    // left or right neighbor in all
+                                    // following computations.
+                                    //
                                     // first index:  dimension (minus 2)
                                     // second index: local face index
                                     // third index:  face_orientation (false and true)
@@ -1111,9 +1132,18 @@ namespace
     for (; cell != endc; ++cell)
       for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
        {
-         const typename Triangulation<dim,spacedim>::face_iterator face=cell->face(f);
-         const unsigned int offset=left_right_offset[dim-2][f][cell->face_orientation(f)];
+         const typename Triangulation<dim,spacedim>::face_iterator
+           face=cell->face(f);
+
+         const unsigned int
+           offset = (cell->direction_flag()
+                     ?
+                     left_right_offset[dim-2][f][cell->face_orientation(f)]
+                     :
+                     1-left_right_offset[dim-2][f][cell->face_orientation(f)]);
+
          adjacent_cells[2*face->index() + offset] = cell;
+
                                           // if this cell is not refined, but the
                                           // face is, then we'll have to set our
                                           // cell as neighbor for the cild faces
@@ -1178,9 +1208,16 @@ namespace
                                     // the offset of the neighbor, not our own.
     for (cell=triangulation.begin(); cell != endc; ++cell)
       for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-       cell->set_neighbor(f,
-                          adjacent_cells[2*cell->face(f)->index() + 1
-                                         - left_right_offset[dim-2][f][cell->face_orientation(f)]]);
+       {
+         const unsigned int
+           offset = (cell->direction_flag()
+                     ?
+                     left_right_offset[dim-2][f][cell->face_orientation(f)]
+                     :
+                     1-left_right_offset[dim-2][f][cell->face_orientation(f)]);
+         cell->set_neighbor(f,
+                            adjacent_cells[2*cell->face(f)->index() + 1 - offset]);
+       }
   }
 
 }// end of anonymous namespace
diff --git a/tests/codim_one/hanging_nodes_01.cc b/tests/codim_one/hanging_nodes_01.cc
new file mode 100644 (file)
index 0000000..af75710
--- /dev/null
@@ -0,0 +1,73 @@
+//----------------------------  hanging_nodes_01.cc  ---------------------------
+//    $Id: bem.cc 22693 2010-11-11 20:11:47Z kanschat $
+//    Version: $Name$
+//
+//    Copyright (C) 2010, 2011 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  hanging_nodes_01.cc  ---------------------------
+
+
+// we crash when building hanging nodes. however, as elucidated by the
+// _02 testcase, the reason is actually somewhere entirely different:
+// upon refinement of one cell, we forget who the neighbors of another
+// cell are. from there, no recovery is possible
+
+#include "../tests.h"
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_out.h>
+#include <grid/grid_tools.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_q.h>
+#include <lac/constraint_matrix.h>
+
+
+int main ()
+{
+  std::ofstream logfile("hanging_nodes_01/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  const unsigned int spacedim = 3;
+  const unsigned int dim = spacedim-1;
+
+  Triangulation<dim,spacedim> boundary_mesh;
+
+                                  /*****************************************************************/
+                                  // Create Surface Mesh:  Boundary of hypercube without one face
+  {
+    Triangulation<spacedim> volume_mesh;
+    GridGenerator::hyper_cube(volume_mesh);
+    Triangulation<spacedim>::active_cell_iterator
+      cell = volume_mesh.begin_active();
+
+    cell->face(0)->set_all_boundary_indicators (1);
+    std::set<unsigned char> boundary_ids;
+    boundary_ids.insert(0);
+    GridTools::extract_boundary_mesh (volume_mesh, boundary_mesh, boundary_ids);
+  }
+  boundary_mesh.begin_active()->set_refine_flag ();
+  boundary_mesh.execute_coarsening_and_refinement ();
+
+                                  /*****************************************************************/
+  FE_Q<dim,spacedim> fe(1);
+  DoFHandler<dim,spacedim>  dh (boundary_mesh);
+  ConstraintMatrix     hanging_node_constraints;
+
+  dh.distribute_dofs(fe);
+  hanging_node_constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dh,hanging_node_constraints);
+  hanging_node_constraints.close ();
+
+  hanging_node_constraints.print (deallog.get_file_stream());
+
+  return 0;
+}
diff --git a/tests/codim_one/hanging_nodes_01/cmp/generic b/tests/codim_one/hanging_nodes_01/cmp/generic
new file mode 100644 (file)
index 0000000..cec9430
--- /dev/null
@@ -0,0 +1,9 @@
+
+    8 2:  0.500000
+    8 6:  0.500000
+    9 2:  0.500000
+    9 3:  0.500000
+    11 6:  0.500000
+    11 7:  0.500000
+    12 3:  0.500000
+    12 7:  0.500000
diff --git a/tests/codim_one/hanging_nodes_02.cc b/tests/codim_one/hanging_nodes_02.cc
new file mode 100644 (file)
index 0000000..f36d96e
--- /dev/null
@@ -0,0 +1,85 @@
+//----------------------------  hanging_nodes_02.cc  ---------------------------
+//    $Id: bem.cc 22693 2010-11-11 20:11:47Z kanschat $
+//    Version: $Name$
+//
+//    Copyright (C) 2010, 2011 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  hanging_nodes_02.cc  ---------------------------
+
+
+// an extract of the _01 test that shows the essence of what's going
+// wrong. compared to the _03 testcase, upon refinement of the middle
+// cell (0.0), cell 0.1 forgets who some of its neighbors are. this is
+// clearly not good
+//
+// the underlying reason was that when we update neighborship
+// information in the triangulation class, we have to take into
+// account the direction_flag of cells
+
+#include "../tests.h"
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_out.h>
+#include <grid/grid_tools.h>
+
+
+int main ()
+{
+  std::ofstream logfile("hanging_nodes_02/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  const unsigned int spacedim = 3;
+  const unsigned int dim = spacedim-1;
+
+  Triangulation<dim,spacedim> boundary_mesh;
+
+                                  /*****************************************************************/
+                                  // Create Surface Mesh:  Boundary of hypercube without one face
+  {
+    Triangulation<spacedim> volume_mesh;
+    GridGenerator::hyper_cube(volume_mesh);
+    Triangulation<spacedim>::active_cell_iterator
+      cell = volume_mesh.begin_active();
+
+    cell->face(0)->set_all_boundary_indicators (1);
+    std::set<unsigned char> boundary_ids;
+    boundary_ids.insert(0);
+    GridTools::extract_boundary_mesh (volume_mesh, boundary_mesh, boundary_ids);
+  }
+  boundary_mesh.begin_active()->set_refine_flag ();
+  boundary_mesh.execute_coarsening_and_refinement ();
+
+  Triangulation<dim,spacedim>::active_cell_iterator
+    cell = boundary_mesh.begin_active();
+  for (; cell!=boundary_mesh.end(); ++cell)
+    {
+      deallog << "Cell = " << cell << std::endl;
+      deallog << "  direction_flag = " << cell->direction_flag() << std::endl;
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+       {
+         deallog << "  face = " << face
+                 << "  (neighbor = " << cell->neighbor(face) << ")"
+                 << std::endl;
+
+         if (cell->face(face)->has_children())
+           for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
+             {
+               deallog << "    subface = " << c << std::endl;
+               deallog << "              "
+                       << cell->neighbor_child_on_subface(face, c)
+                       << std::endl;
+             }
+       }
+    }
+
+  return 0;
+}
diff --git a/tests/codim_one/hanging_nodes_02/cmp/generic b/tests/codim_one/hanging_nodes_02/cmp/generic
new file mode 100644 (file)
index 0000000..c37aa45
--- /dev/null
@@ -0,0 +1,65 @@
+
+DEAL::Cell = 0.1
+DEAL::  direction_flag = 0
+DEAL::  face = 0  (neighbor = 0.3)
+DEAL::  face = 1  (neighbor = 0.4)
+DEAL::  face = 2  (neighbor = -1.-1)
+DEAL::  face = 3  (neighbor = 0.0)
+DEAL::    subface = 0
+DEAL::              1.0
+DEAL::    subface = 1
+DEAL::              1.2
+DEAL::Cell = 0.2
+DEAL::  direction_flag = 1
+DEAL::  face = 0  (neighbor = 0.3)
+DEAL::  face = 1  (neighbor = 0.4)
+DEAL::  face = 2  (neighbor = -1.-1)
+DEAL::  face = 3  (neighbor = 0.0)
+DEAL::    subface = 0
+DEAL::              1.1
+DEAL::    subface = 1
+DEAL::              1.3
+DEAL::Cell = 0.3
+DEAL::  direction_flag = 0
+DEAL::  face = 0  (neighbor = -1.-1)
+DEAL::  face = 1  (neighbor = 0.0)
+DEAL::    subface = 0
+DEAL::              1.0
+DEAL::    subface = 1
+DEAL::              1.1
+DEAL::  face = 2  (neighbor = 0.1)
+DEAL::  face = 3  (neighbor = 0.2)
+DEAL::Cell = 0.4
+DEAL::  direction_flag = 1
+DEAL::  face = 0  (neighbor = -1.-1)
+DEAL::  face = 1  (neighbor = 0.0)
+DEAL::    subface = 0
+DEAL::              1.2
+DEAL::    subface = 1
+DEAL::              1.3
+DEAL::  face = 2  (neighbor = 0.1)
+DEAL::  face = 3  (neighbor = 0.2)
+DEAL::Cell = 1.0
+DEAL::  direction_flag = 1
+DEAL::  face = 0  (neighbor = 0.1)
+DEAL::  face = 1  (neighbor = 1.1)
+DEAL::  face = 2  (neighbor = 0.3)
+DEAL::  face = 3  (neighbor = 1.2)
+DEAL::Cell = 1.1
+DEAL::  direction_flag = 1
+DEAL::  face = 0  (neighbor = 1.0)
+DEAL::  face = 1  (neighbor = 0.2)
+DEAL::  face = 2  (neighbor = 0.3)
+DEAL::  face = 3  (neighbor = 1.3)
+DEAL::Cell = 1.2
+DEAL::  direction_flag = 1
+DEAL::  face = 0  (neighbor = 0.1)
+DEAL::  face = 1  (neighbor = 1.3)
+DEAL::  face = 2  (neighbor = 1.0)
+DEAL::  face = 3  (neighbor = 0.4)
+DEAL::Cell = 1.3
+DEAL::  direction_flag = 1
+DEAL::  face = 0  (neighbor = 1.2)
+DEAL::  face = 1  (neighbor = 0.2)
+DEAL::  face = 2  (neighbor = 1.1)
+DEAL::  face = 3  (neighbor = 0.4)
diff --git a/tests/codim_one/hanging_nodes_03.cc b/tests/codim_one/hanging_nodes_03.cc
new file mode 100644 (file)
index 0000000..8fa36d7
--- /dev/null
@@ -0,0 +1,74 @@
+//----------------------------  hanging_nodes_03.cc  ---------------------------
+//    $Id: bem.cc 22693 2010-11-11 20:11:47Z kanschat $
+//    Version: $Name$
+//
+//    Copyright (C) 2010, 2011 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  hanging_nodes_03.cc  ---------------------------
+
+
+// a further extract of the _02 test. the results here are correct and
+// show the relationship between the various cells of the mesh
+
+#include "../tests.h"
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_out.h>
+#include <grid/grid_tools.h>
+
+
+int main ()
+{
+  std::ofstream logfile("hanging_nodes_03/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  const unsigned int spacedim = 3;
+  const unsigned int dim = spacedim-1;
+
+  Triangulation<dim,spacedim> boundary_mesh;
+
+                                  /*****************************************************************/
+                                  // Create Surface Mesh:  Boundary of hypercube without one face
+  {
+    Triangulation<spacedim> volume_mesh;
+    GridGenerator::hyper_cube(volume_mesh);
+    Triangulation<spacedim>::active_cell_iterator
+      cell = volume_mesh.begin_active();
+
+    cell->face(0)->set_all_boundary_indicators (1);
+    std::set<unsigned char> boundary_ids;
+    boundary_ids.insert(0);
+    GridTools::extract_boundary_mesh (volume_mesh, boundary_mesh, boundary_ids);
+  }
+
+  Triangulation<dim,spacedim>::active_cell_iterator
+    cell = boundary_mesh.begin_active();
+  for (; cell!=boundary_mesh.end(); ++cell)
+    {
+      deallog << "Cell = " << cell << std::endl;
+      deallog << "  direction_flag = " << cell->direction_flag() << std::endl;
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+       {
+         deallog << "  face = " << face
+                 << "  (neighbor = " << cell->neighbor(face) << ")"
+                 << std::endl;
+
+         if (cell->face(face)->has_children())
+           for (unsigned int c=0; c<cell->face(face)->n_children(); ++c)
+             {
+               deallog << "    subface = " << c << std::endl;
+             }
+       }
+    }
+
+  return 0;
+}
diff --git a/tests/codim_one/hanging_nodes_03/cmp/generic b/tests/codim_one/hanging_nodes_03/cmp/generic
new file mode 100644 (file)
index 0000000..fc1fbf7
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::Cell = 0.0
+DEAL::  direction_flag = 1
+DEAL::  face = 0  (neighbor = 0.1)
+DEAL::  face = 1  (neighbor = 0.2)
+DEAL::  face = 2  (neighbor = 0.3)
+DEAL::  face = 3  (neighbor = 0.4)
+DEAL::Cell = 0.1
+DEAL::  direction_flag = 0
+DEAL::  face = 0  (neighbor = 0.3)
+DEAL::  face = 1  (neighbor = 0.4)
+DEAL::  face = 2  (neighbor = -1.-1)
+DEAL::  face = 3  (neighbor = 0.0)
+DEAL::Cell = 0.2
+DEAL::  direction_flag = 1
+DEAL::  face = 0  (neighbor = 0.3)
+DEAL::  face = 1  (neighbor = 0.4)
+DEAL::  face = 2  (neighbor = -1.-1)
+DEAL::  face = 3  (neighbor = 0.0)
+DEAL::Cell = 0.3
+DEAL::  direction_flag = 0
+DEAL::  face = 0  (neighbor = -1.-1)
+DEAL::  face = 1  (neighbor = 0.0)
+DEAL::  face = 2  (neighbor = 0.1)
+DEAL::  face = 3  (neighbor = 0.2)
+DEAL::Cell = 0.4
+DEAL::  direction_flag = 1
+DEAL::  face = 0  (neighbor = -1.-1)
+DEAL::  face = 1  (neighbor = 0.0)
+DEAL::  face = 2  (neighbor = 0.1)
+DEAL::  face = 3  (neighbor = 0.2)

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.