From 39a1268c963a9178a0b197fe589c0555355c1faa Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 3 Jan 2011 18:10:50 +0000 Subject: [PATCH] Fix a problem where the update_neighbors function in the Triangulation class needed to learn about the newly introduced direction_flag used in codim-one meshes. git-svn-id: https://svn.dealii.org/trunk@23115 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/grid/tria.cc | 53 ++++++++++-- tests/codim_one/hanging_nodes_01.cc | 73 +++++++++++++++++ tests/codim_one/hanging_nodes_01/cmp/generic | 9 +++ tests/codim_one/hanging_nodes_02.cc | 85 ++++++++++++++++++++ tests/codim_one/hanging_nodes_02/cmp/generic | 65 +++++++++++++++ tests/codim_one/hanging_nodes_03.cc | 74 +++++++++++++++++ tests/codim_one/hanging_nodes_03/cmp/generic | 31 +++++++ 7 files changed, 382 insertions(+), 8 deletions(-) create mode 100644 tests/codim_one/hanging_nodes_01.cc create mode 100644 tests/codim_one/hanging_nodes_01/cmp/generic create mode 100644 tests/codim_one/hanging_nodes_02.cc create mode 100644 tests/codim_one/hanging_nodes_02/cmp/generic create mode 100644 tests/codim_one/hanging_nodes_03.cc create mode 100644 tests/codim_one/hanging_nodes_03/cmp/generic diff --git a/deal.II/source/grid/tria.cc b/deal.II/source/grid/tria.cc index 5a43056e79..177f8aa2cc 100644 --- a/deal.II/source/grid/tria.cc +++ b/deal.II/source/grid/tria.cc @@ -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::faces_per_cell; ++f) { - const typename Triangulation::face_iterator face=cell->face(f); - const unsigned int offset=left_right_offset[dim-2][f][cell->face_orientation(f)]; + const typename Triangulation::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::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 index 0000000000..af757103b0 --- /dev/null +++ b/tests/codim_one/hanging_nodes_01.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include + + +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 boundary_mesh; + + /*****************************************************************/ + // Create Surface Mesh: Boundary of hypercube without one face + { + Triangulation volume_mesh; + GridGenerator::hyper_cube(volume_mesh); + Triangulation::active_cell_iterator + cell = volume_mesh.begin_active(); + + cell->face(0)->set_all_boundary_indicators (1); + std::set 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 fe(1); + DoFHandler 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 index 0000000000..cec94304bd --- /dev/null +++ b/tests/codim_one/hanging_nodes_01/cmp/generic @@ -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 index 0000000000..f36d96e426 --- /dev/null +++ b/tests/codim_one/hanging_nodes_02.cc @@ -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 +#include +#include +#include +#include + + +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 boundary_mesh; + + /*****************************************************************/ + // Create Surface Mesh: Boundary of hypercube without one face + { + Triangulation volume_mesh; + GridGenerator::hyper_cube(volume_mesh); + Triangulation::active_cell_iterator + cell = volume_mesh.begin_active(); + + cell->face(0)->set_all_boundary_indicators (1); + std::set 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::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::faces_per_cell; ++face) + { + deallog << " face = " << face + << " (neighbor = " << cell->neighbor(face) << ")" + << std::endl; + + if (cell->face(face)->has_children()) + for (unsigned int c=0; cface(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 index 0000000000..c37aa459bc --- /dev/null +++ b/tests/codim_one/hanging_nodes_02/cmp/generic @@ -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 index 0000000000..8fa36d7da4 --- /dev/null +++ b/tests/codim_one/hanging_nodes_03.cc @@ -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 +#include +#include +#include +#include + + +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 boundary_mesh; + + /*****************************************************************/ + // Create Surface Mesh: Boundary of hypercube without one face + { + Triangulation volume_mesh; + GridGenerator::hyper_cube(volume_mesh); + Triangulation::active_cell_iterator + cell = volume_mesh.begin_active(); + + cell->face(0)->set_all_boundary_indicators (1); + std::set boundary_ids; + boundary_ids.insert(0); + GridTools::extract_boundary_mesh (volume_mesh, boundary_mesh, boundary_ids); + } + + Triangulation::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::faces_per_cell; ++face) + { + deallog << " face = " << face + << " (neighbor = " << cell->neighbor(face) << ")" + << std::endl; + + if (cell->face(face)->has_children()) + for (unsigned int c=0; cface(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 index 0000000000..fc1fbf7ac6 --- /dev/null +++ b/tests/codim_one/hanging_nodes_03/cmp/generic @@ -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) -- 2.39.5