]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Trying to interpolate boundary values on a closed loop 1d mesh in 2d resulted in...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Nov 2010 20:25:38 +0000 (20:25 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Nov 2010 20:25:38 +0000 (20:25 +0000)
git-svn-id: https://svn.dealii.org/trunk@22856 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/numerics/vectors.templates.h
tests/codim_one/interpolate_boundary_values_1d_closed_ring.cc [new file with mode: 0644]
tests/codim_one/interpolate_boundary_values_1d_closed_ring/cmp/generic [new file with mode: 0644]

index 8bd2b0d6126b719c83960210d0b2fb2b9f5bda2f..ec4ac460eb73fdd61ed5ba7e3061ba9454f9d476 100644 (file)
@@ -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<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);
@@ -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<dim>& 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 (file)
index 0000000..ba5d43a
--- /dev/null
@@ -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 <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;
+}
+
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 (file)
index 0000000..858bbf1
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::8 degrees of freedom
+DEAL::0 boundary degrees of freedom
+DEAL::0 boundary degrees of freedom

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.