Fix a bug in DoFTools::make_flux_sparsity_pattern originally reported by Baerbel.

git-svn-id: https://svn.dealii.org/trunk@27655 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index db2c26e5a215106e688e914c1e610afe17ad4dfb..3a5b562c217d73c0bb37f8e12ec565d6dd42d42f 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -154,6 +154,13 @@ DoFHandler, in particular removal of specializations.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: The DoFTools::make_flux_sparsity_pattern() function
+had a bug that triggered in 1d whenever there were neighboring
+cells that differed in refinement level by more than one. This
+is now fixed.
+<br>
+(Wolfgang Bangerth, 2012/11/20)
+
 <li> New: The TrilinosWrappers::PreconditionBase class now has
 a function TrilinosWrappers::PreconditionBase::Tvmult that
 allows applying the transpose preconditioner.

diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc
index 98e49222b91fdfdd8a4922fd995b7c3ee4b18db9..8a0860406513186b273c99ef726f45617cc3b9cb 100644 (file)
--- a/deal.II/source/dofs/dof_tools.cc
+++ b/deal.II/source/dofs/dof_tools.cc
@@ -581,6 +581,14 @@ namespace DoFTools
                 {
                   typename DH::cell_iterator neighbor = cell->neighbor(face);
 
+                  // in 1d, we do not need to worry whether the neighbor
+                  // might have children and then loop over those children.
+                  // rather, we may as well go straight to to cell behind
+                  // this particular cell's most terminal child
+                  if (DH::dimension==1)
+                    while (neighbor->has_children())
+                      neighbor = neighbor->child(face==0 ? 1 : 0);
+
                   if (neighbor->has_children())
                     {
                       for (unsigned int sub_nr = 0;

diff --git a/tests/bits/dof_tools_18a_1d.cc b/tests/bits/dof_tools_18a_1d.cc
new file mode 100644 (file)
index 0000000..813f258
--- /dev/null
+++ b/tests/bits/dof_tools_18a_1d.cc
@@ -0,0 +1,154 @@
+//----------------------------  dof_tools_1a.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2003, 2004, 2005, 2006, 2012 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.
+//
+//----------------------------  dof_tools_1a.cc  ---------------------------
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <fstream>
+#include <iomanip>
+#include <iomanip>
+#include <string>
+
+// check
+//   DoFTools::
+//   make_flux_sparsity_pattern (const DoFHandler<dim> &,
+//                              SparsityPattern       &,
+//                              ConstraintMatrix      &);
+// this used to fail at one point because we forgot that in 1d
+// neighboring cells can be more than one level apart. with
+// FE_DGQ(0) elements, we get the following mesh and DoF indices:
+//
+// *-*-*-*-*-------*
+//  1 2 3 4    0
+//
+// we need to make sure that in row 0 there is a coupling to
+// DoF 4. this was previously missing
+
+
+std::string output_file_name = "dof_tools_18a_1d/output";
+
+
+template <int dim>
+void
+check_this (const DoFHandler<dim> &dof_handler)
+{
+                                   // create sparsity pattern
+  SparsityPattern sp (dof_handler.n_dofs(),
+                      dof_handler.max_couplings_between_dofs()*2);
+  DoFTools::make_flux_sparsity_pattern (dof_handler, sp,
+                                       ConstraintMatrix());
+  sp.compress ();
+
+                                   // write out 20 lines of this
+                                   // pattern (if we write out the
+                                   // whole pattern, the output file
+                                   // would be in the range of 40 MB)
+  for (unsigned int l=0; l<sp.n_rows(); ++l)
+    {
+      for (unsigned int c=0; c<sp.row_length(l); ++c)
+        deallog << sp.column_number(l,c) << " ";
+      deallog << std::endl;
+    }
+
+                                   // write out some other indicators
+  deallog << sp.bandwidth () << std::endl
+          << sp.max_entries_per_row () << std::endl
+          << sp.n_nonzero_elements () << std::endl;
+
+  unsigned int hash = 0;
+  for (unsigned int l=0; l<sp.n_rows(); ++l)
+    hash += l*(sp.row_length(l) +
+               sp.get_rowstart_indices()[l] +
+               sp.get_column_numbers()[sp.get_rowstart_indices()[l]
+                                       +
+                                       (sp.row_length(l)>1 ? 1 : 0)]);
+  deallog << hash << std::endl;
+}
+
+
+void
+check_this ()
+{
+                                  // create a mesh where two cells at levels
+                                  // 1 and 3 are adjacent
+  const int dim = 1;
+  Triangulation<dim> tr;
+  GridGenerator::hyper_cube (tr);
+  tr.refine_global(1);
+  tr.begin_active()->set_refine_flag();
+  tr.execute_coarsening_and_refinement();
+  for (Triangulation<dim>::active_cell_iterator
+        c=tr.begin_active(2); c!=tr.end_active(2); ++c)
+    c->set_refine_flag();
+  tr.execute_coarsening_and_refinement();
+
+  FE_DGQ<dim> fe(0);
+  DoFHandler<dim> dof_handler(tr);
+  dof_handler.distribute_dofs (fe);
+
+  check_this (dof_handler);
+}
+
+
+int
+main()
+{
+  try
+    {
+      std::ofstream logfile(output_file_name.c_str());
+      logfile << std::setprecision (2);
+      deallog << std::setprecision (2);
+      deallog.attach(logfile);
+      deallog.depth_console(0);
+      deallog.threshold_double(1.e-10);
+
+      check_this ();
+      return 0;
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      deallog << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      deallog << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}

diff --git a/tests/bits/dof_tools_18a_1d/cmp/generic b/tests/bits/dof_tools_18a_1d/cmp/generic
new file mode 100644 (file)
index 0000000..3e12661
--- /dev/null
+++ b/tests/bits/dof_tools_18a_1d/cmp/generic
@@ -0,0 +1,10 @@
+
+DEAL::0 4 
+DEAL::1 2 
+DEAL::2 1 3 
+DEAL::3 2 4 
+DEAL::4 0 3 
+DEAL::4
+DEAL::3
+DEAL::13
+DEAL::110