From d8e175d65804e29c364bba9382828b0c05bf444e Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 24 Nov 2010 20:25:38 +0000 Subject: [PATCH] Trying to interpolate boundary values on a closed loop 1d mesh in 2d resulted in an infinite loop. This should now be fixed. git-svn-id: https://svn.dealii.org/trunk@22856 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/numerics/vectors.templates.h | 66 ++++++++--- ...erpolate_boundary_values_1d_closed_ring.cc | 108 ++++++++++++++++++ .../cmp/generic | 4 + 3 files changed, 164 insertions(+), 14 deletions(-) create mode 100644 tests/codim_one/interpolate_boundary_values_1d_closed_ring.cc create mode 100644 tests/codim_one/interpolate_boundary_values_1d_closed_ring/cmp/generic diff --git a/deal.II/include/deal.II/numerics/vectors.templates.h b/deal.II/include/deal.II/numerics/vectors.templates.h index 8bd2b0d612..ec4ac460eb 100644 --- a/deal.II/include/deal.II/numerics/vectors.templates.h +++ b/deal.II/include/deal.II/numerics/vectors.templates.h @@ -1277,27 +1277,65 @@ namespace internal Assert (function_map.find(255) == function_map.end(), dealii::VectorTools::ExcInvalidBoundaryIndicator()); + // if this 1d mesh is embedded + // in a higher dimensional + // space, it may be that it is + // a circle; in that case, + // there are no boundary nodes + // and we want to catch this + // now before we get into what + // turns out to be an endless + // loop further down below + // where we try to find the + // cell at the left or right + // edge of the domain + if (DH::space_dimension > 1) + { + typename DH::cell_iterator cell = dof.begin(0); + unsigned int count = 0; + while (cell->neighbor(1).state() == IteratorState::valid) + { + cell = cell->neighbor(1); + ++count; + + if (count > dof.get_tria().n_cells(0)) + // we looped back + // around. there is + // nothing for us to + // do here + return; + } + } + + // now do the actual work for (typename FunctionMap::type::const_iterator i=function_map.begin(); i!=function_map.end(); ++i) { - // check whether boundary values at - // the left or right boundary of - // the line are - // requested. direction denotes - // the neighboring direction in - // which we seek the boundary, - // i.e. 0 is left boundary and 1 is - // right. + // check whether boundary values at + // the left or right boundary of + // the line are + // requested. direction denotes + // the neighboring direction in + // which we seek the boundary, + // i.e. 0 is left boundary and 1 is + // right. const unsigned int direction = i->first; Assert (direction < 2, dealii::VectorTools::ExcInvalidBoundaryIndicator()); const Function & boundary_function = *(i->second); - // first find the outermost active - // cell by first traversing the coarse - // grid to its end and then going - // to the children + // first find the outermost active + // cell by first traversing the coarse + // grid to its end and then going + // to the children. + // + // note that we have + // already ruled out the + // case that the mesh is a + // topological circle + // above, so we don't have + // to check here any more typename DH::cell_iterator outermost_cell = dof.begin(0); while (outermost_cell->neighbor(direction).state() == IteratorState::valid) outermost_cell = outermost_cell->neighbor(direction); @@ -2886,7 +2924,7 @@ namespace internals { { 2, 0 }, { 0, 1 }, { 0, 1 } }; - + // The projection is // divided into two steps. // In the first step we @@ -3547,7 +3585,7 @@ project_boundary_values_curl_conforming (const hp::DoFHandler& dof_handler, constraints.add_line (dof); constraints.set_inhomogeneity (dof, computed_constraints[dof]); } - + break; } diff --git a/tests/codim_one/interpolate_boundary_values_1d_closed_ring.cc b/tests/codim_one/interpolate_boundary_values_1d_closed_ring.cc new file mode 100644 index 0000000000..ba5d43a9e9 --- /dev/null +++ b/tests/codim_one/interpolate_boundary_values_1d_closed_ring.cc @@ -0,0 +1,108 @@ + +//---------------------------- template.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2008, 2010 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. +// +//---------------------------- template.cc --------------------------- + + +// test VectorTools::interpolate_boundary_values for codim=1. like _01 +// but for 1d triangulations +// +// the function ran into an endless loop whenever the 1d codim=1 mesh +// was topologically a circle, since it tries to find the left or +// right most vertex by starting at an arbitrary cell and always going +// to the left/right. If the mesh is closed, this of course doesn't +// work so well. + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace dealii; + + +std::ofstream logfile("interpolate_boundary_values_1d_closed_ring/output"); + +void test() { + const unsigned int dim = 1; + const unsigned int spacedim=2; + + Triangulation tria; + std::map::cell_iterator,Triangulation::face_iterator > + surface_to_volume_mapping; + Triangulation volume_mesh; + GridGenerator::hyper_cube(volume_mesh); + + GridTools::extract_boundary_mesh (volume_mesh, tria, + surface_to_volume_mapping); + + FE_Q fe(2); + DoFHandler dof_handler (tria); + dof_handler.distribute_dofs (fe); + + deallog << dof_handler.n_dofs() << " degrees of freedom" << std::endl; + + // test left and right boundary + // separatel + for (unsigned int boundary_id=0; boundary_id<2; ++boundary_id) + { + std::map bv; + VectorTools::interpolate_boundary_values (dof_handler, + boundary_id, + Functions::SquareFunction(), + bv); + deallog << bv.size() << " boundary degrees of freedom" << std::endl; + + for (std::map::const_iterator i = bv.begin(); + i != bv.end(); ++i) + deallog << i->first << ' ' << i->second << std::endl; + + for (DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->at_boundary(f) && + (cell->face(f)->boundary_indicator() == boundary_id)) + for (unsigned int v=0; v::vertices_per_face; ++v) + for (unsigned int i=0; iface(f)->vertex_dof_index(v,i)) + != bv.end(), + ExcInternalError()); + Assert (bv[cell->face(f)->vertex_dof_index(v,i)] + == + Functions::SquareFunction() + .value(cell->face(f)->vertex(v),i), + ExcInternalError()); + } + } +} + + + +int main () +{ + deallog.attach(logfile); + deallog.depth_console(0); + + test(); + + return 0; +} + diff --git a/tests/codim_one/interpolate_boundary_values_1d_closed_ring/cmp/generic b/tests/codim_one/interpolate_boundary_values_1d_closed_ring/cmp/generic new file mode 100644 index 0000000000..858bbf138a --- /dev/null +++ b/tests/codim_one/interpolate_boundary_values_1d_closed_ring/cmp/generic @@ -0,0 +1,4 @@ + +DEAL::8 degrees of freedom +DEAL::0 boundary degrees of freedom +DEAL::0 boundary degrees of freedom -- 2.39.5