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<spacedim>::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<DH::space_dimension> & 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);
{ 2, 0 },
{ 0, 1 },
{ 0, 1 } };
-
+
// The projection is
// divided into two steps.
// In the first step we
constraints.add_line (dof);
constraints.set_inhomogeneity (dof, computed_constraints[dof]);
}
-
+
break;
}
--- /dev/null
+
+//---------------------------- 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 <fstream>
+#include <base/logstream.h>
+#include <base/function_lib.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_tools.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <numerics/vectors.h>
+
+#include <string>
+
+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<dim,spacedim> tria;
+ std::map<Triangulation<dim,spacedim>::cell_iterator,Triangulation<spacedim,spacedim>::face_iterator >
+ surface_to_volume_mapping;
+ Triangulation<spacedim> volume_mesh;
+ GridGenerator::hyper_cube(volume_mesh);
+
+ GridTools::extract_boundary_mesh (volume_mesh, tria,
+ surface_to_volume_mapping);
+
+ FE_Q<dim,spacedim> fe(2);
+ DoFHandler<dim,spacedim> 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<unsigned int, double> bv;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ boundary_id,
+ Functions::SquareFunction<spacedim>(),
+ bv);
+ deallog << bv.size() << " boundary degrees of freedom" << std::endl;
+
+ for (std::map<unsigned int, double>::const_iterator i = bv.begin();
+ i != bv.end(); ++i)
+ deallog << i->first << ' ' << i->second << std::endl;
+
+ for (DoFHandler<dim,spacedim>::active_cell_iterator
+ cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell->at_boundary(f) &&
+ (cell->face(f)->boundary_indicator() == boundary_id))
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+ for (unsigned int i=0; i<fe.dofs_per_vertex; ++i)
+ {
+ Assert (bv.find(cell->face(f)->vertex_dof_index(v,i))
+ != bv.end(),
+ ExcInternalError());
+ Assert (bv[cell->face(f)->vertex_dof_index(v,i)]
+ ==
+ Functions::SquareFunction<spacedim>()
+ .value(cell->face(f)->vertex(v),i),
+ ExcInternalError());
+ }
+ }
+}
+
+
+
+int main ()
+{
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ test();
+
+ return 0;
+}
+
projects / dealii-svn.git / commitdiff