// $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
// 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
// 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)
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
// 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
--- /dev/null
+//---------------------------- 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;
+}
--- /dev/null
+
+ 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
--- /dev/null
+//---------------------------- 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;
+}
--- /dev/null
+
+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)
--- /dev/null
+//---------------------------- 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;
+}
--- /dev/null
+
+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)
projects / dealii-svn.git / commitdiff