also make make_hanging_node_constraints work in 3d with FE_RaviartThomas

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

diff --git a/deal.II/source/dofs/dof_tools_constraints.cc b/deal.II/source/dofs/dof_tools_constraints.cc
index 2f32ccfacf637adb3481395b19b4a6ecb5b81d73..57bfd9c0b21795f64ea6e9c491b7d4af61d9c7ca 100644 (file)
--- a/deal.II/source/dofs/dof_tools_constraints.cc
+++ b/deal.II/source/dofs/dof_tools_constraints.cc
@@ -972,7 +972,10 @@ namespace DoFTools
                 //TODO[TL]: think about this and the following in case of anisotropic refinement
 
                 dofs_on_mother.resize (n_dofs_on_mother);
-                dofs_on_children.resize (n_dofs_on_children);
+                // we might not use all of those in case of artificial cells,
+                // so do not resize(), but reserve() and use push_back later.
+                dofs_on_children.clear();
+                dofs_on_children.reserve (n_dofs_on_children);
 
                 Assert(n_dofs_on_mother == fe.constraints().n(),
                        ExcDimensionMismatch(n_dofs_on_mother,
@@ -998,9 +1001,7 @@ namespace DoFTools
                   dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index);
                 AssertDimension (next_index, dofs_on_mother.size());
 
-                next_index = 0;
-
-                // assert some consistency assumptions
+                //TODO: assert some consistency assumptions
 
                 //TODO[TL]: think about this in case of anisotropic
                 //refinement
@@ -1013,46 +1014,56 @@ namespace DoFTools
                          (this_face->child(0)->vertex_index(3) ==
                           this_face->child(3)->vertex_index(0))),
                         ExcInternalError());
+
                 for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
-                  dofs_on_children[next_index++]
-                    = this_face->child(0)->vertex_dof_index(3,dof);
+                  dofs_on_children.push_back(
+                    this_face->child(0)->vertex_dof_index(3,dof));
 
                 // dof numbers on the centers of the lines bounding this
                 // face
                 for (unsigned int line=0; line<4; ++line)
                   for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
-                    dofs_on_children[next_index++]
-                      = this_face->line(line)->child(0)->vertex_dof_index(1,dof, fe_index);
+                    dofs_on_children.push_back(
+                      this_face->line(line)->child(0)->vertex_dof_index(1,dof, fe_index));
 
                 // next the dofs on the lines interior to the face; the
                 // order of these lines is laid down in the FiniteElement
                 // class documentation
                 for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-                  dofs_on_children[next_index++]
-                    = this_face->child(0)->line(1)->dof_index(dof, fe_index);
+                  dofs_on_children.push_back(
+                    this_face->child(0)->line(1)->dof_index(dof, fe_index));
                 for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-                  dofs_on_children[next_index++]
-                    = this_face->child(2)->line(1)->dof_index(dof, fe_index);
+                  dofs_on_children.push_back(
+                    this_face->child(2)->line(1)->dof_index(dof, fe_index));
                 for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-                  dofs_on_children[next_index++]
-                    = this_face->child(0)->line(3)->dof_index(dof, fe_index);
+                  dofs_on_children.push_back(
+                    this_face->child(0)->line(3)->dof_index(dof, fe_index));
                 for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-                  dofs_on_children[next_index++]
-                    = this_face->child(1)->line(3)->dof_index(dof, fe_index);
+                  dofs_on_children.push_back(
+                    this_face->child(1)->line(3)->dof_index(dof, fe_index));
 
                 // dofs on the bordering lines
                 for (unsigned int line=0; line<4; ++line)
                   for (unsigned int child=0; child<2; ++child)
-                    for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
-                      dofs_on_children[next_index++]
-                        = this_face->line(line)->child(child)->dof_index(dof, fe_index);
+                    {
+                      for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+                        dofs_on_children.push_back(
+                          this_face->line(line)->child(child)->dof_index(dof, fe_index));
+                    }
 
                 // finally, for the dofs interior to the four child faces
                 for (unsigned int child=0; child<4; ++child)
-                  for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
-                    dofs_on_children[next_index++]
-                      = this_face->child(child)->dof_index(dof, fe_index);
-                AssertDimension (next_index, dofs_on_children.size());
+                  {
+                    // skip artificial cells
+                    if (cell->neighbor_child_on_subface (face, child)->is_artificial())
+                      continue;
+                    for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
+                        dofs_on_children.push_back(
+                          this_face->child(child)->dof_index(dof, fe_index));
+                  }
+
+                // note: can get fewer DoFs when we have artificial cells:
+                Assert(dofs_on_children.size() <= n_dofs_on_children, ExcInternalError());
 
                 // for each row in the constraint matrix for this line:
                 for (unsigned int row=0; row!=dofs_on_children.size(); ++row)

diff --git a/tests/mpi/crash_06.cc b/tests/mpi/crash_06.cc
new file mode 100644 (file)
index 0000000..f6fcc5d
--- /dev/null
+++ b/tests/mpi/crash_06.cc
@@ -0,0 +1,104 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2008 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test different elements for hanging_node constraints. This extends crash_05.cc
+// to make sure stuff is working.
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/intergrid_map.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <fstream>
+#include <cstdlib>
+#include <numeric>
+
+
+template<int dim>
+void test(FiniteElement<dim> &fe)
+{
+  deallog << "dim=" << dim << std::endl;
+  
+    parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+    GridGenerator::hyper_ball(triangulation);
+    triangulation.refine_global(1);
+    
+    if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+      {
+       typename Triangulation<dim>::active_cell_iterator
+         cell = triangulation.begin_active();
+       for (;cell!=triangulation.end();++cell)
+         if (Testing::rand()%2)
+           cell->set_refine_flag ();
+      }
+
+    triangulation.execute_coarsening_and_refinement();
+    deallog << "n_cells: " << triangulation.n_global_active_cells() << std::endl;
+    DoFHandler<dim> dof_handler(triangulation);
+
+    dof_handler.distribute_dofs(fe);
+       
+    ConstraintMatrix constraints;
+    DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+    
+    if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
+      deallog << "OK" << std::endl;
+}
+
+template <int dim>
+void testit()
+{
+  std::vector<std_cxx1x::shared_ptr<FiniteElement<dim> > > fes;
+  fes.push_back(std_cxx1x::shared_ptr<FiniteElement<dim> >(new FE_RaviartThomas<dim>(0)));
+  fes.push_back(std_cxx1x::shared_ptr<FiniteElement<dim> >(new FE_RaviartThomas<dim>(1)));
+  fes.push_back(std_cxx1x::shared_ptr<FiniteElement<dim> >(new FE_Nedelec<dim>(0)));
+  fes.push_back(std_cxx1x::shared_ptr<FiniteElement<dim> >(new FE_Nedelec<dim>(1)));
+  fes.push_back(std_cxx1x::shared_ptr<FiniteElement<dim> >(new FE_Q<dim>(3)));
+  fes.push_back(std_cxx1x::shared_ptr<FiniteElement<dim> >(new FE_DGQ<dim>(2)));
+  fes.push_back(std_cxx1x::shared_ptr<FiniteElement<dim> >(new FE_Q_DG0<dim>(2)));
+
+  for (unsigned int i=0;i<fes.size();++i)
+    {
+      deallog << fes[i]->get_name() << std::endl;
+      test<dim>(*fes[i]);
+    }
+  
+}
+
+
+int main(int argc, char *argv[])
+{  
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll log;
+
+  testit<2>();
+  testit<3>();
+}

diff --git a/tests/mpi/crash_06.mpirun=2.output b/tests/mpi/crash_06.mpirun=2.output
new file mode 100644 (file)
index 0000000..172d74a
--- /dev/null
+++ b/tests/mpi/crash_06.mpirun=2.output
@@ -0,0 +1,101 @@
+
+DEAL:0::FE_RaviartThomas<2>(0)
+DEAL:0::dim=2
+DEAL:0::n_cells: 41
+DEAL:0::OK
+DEAL:0::FE_RaviartThomas<2>(1)
+DEAL:0::dim=2
+DEAL:0::n_cells: 38
+DEAL:0::OK
+DEAL:0::FE_Nedelec<2>(0)
+DEAL:0::dim=2
+DEAL:0::n_cells: 41
+DEAL:0::OK
+DEAL:0::FE_Nedelec<2>(1)
+DEAL:0::dim=2
+DEAL:0::n_cells: 41
+DEAL:0::OK
+DEAL:0::FE_Q<2>(3)
+DEAL:0::dim=2
+DEAL:0::n_cells: 32
+DEAL:0::OK
+DEAL:0::FE_DGQ<2>(2)
+DEAL:0::dim=2
+DEAL:0::n_cells: 38
+DEAL:0::OK
+DEAL:0::FE_Q_DG0<2>(2)
+DEAL:0::dim=2
+DEAL:0::n_cells: 50
+DEAL:0::OK
+DEAL:0::FE_RaviartThomas<3>(0)
+DEAL:0::dim=3
+DEAL:0::n_cells: 189
+DEAL:0::OK
+DEAL:0::FE_RaviartThomas<3>(1)
+DEAL:0::dim=3
+DEAL:0::n_cells: 147
+DEAL:0::OK
+DEAL:0::FE_Nedelec<3>(0)
+DEAL:0::dim=3
+DEAL:0::n_cells: 140
+DEAL:0::OK
+DEAL:0::FE_Nedelec<3>(1)
+DEAL:0::dim=3
+DEAL:0::n_cells: 161
+DEAL:0::OK
+DEAL:0::FE_Q<3>(3)
+DEAL:0::dim=3
+DEAL:0::n_cells: 161
+DEAL:0::OK
+DEAL:0::FE_DGQ<3>(2)
+DEAL:0::dim=3
+DEAL:0::n_cells: 182
+DEAL:0::OK
+DEAL:0::FE_Q_DG0<3>(2)
+DEAL:0::dim=3
+DEAL:0::n_cells: 147
+DEAL:0::OK
+
+DEAL:1::FE_RaviartThomas<2>(0)
+DEAL:1::dim=2
+DEAL:1::n_cells: 41
+DEAL:1::FE_RaviartThomas<2>(1)
+DEAL:1::dim=2
+DEAL:1::n_cells: 38
+DEAL:1::FE_Nedelec<2>(0)
+DEAL:1::dim=2
+DEAL:1::n_cells: 41
+DEAL:1::FE_Nedelec<2>(1)
+DEAL:1::dim=2
+DEAL:1::n_cells: 41
+DEAL:1::FE_Q<2>(3)
+DEAL:1::dim=2
+DEAL:1::n_cells: 32
+DEAL:1::FE_DGQ<2>(2)
+DEAL:1::dim=2
+DEAL:1::n_cells: 38
+DEAL:1::FE_Q_DG0<2>(2)
+DEAL:1::dim=2
+DEAL:1::n_cells: 50
+DEAL:1::FE_RaviartThomas<3>(0)
+DEAL:1::dim=3
+DEAL:1::n_cells: 189
+DEAL:1::FE_RaviartThomas<3>(1)
+DEAL:1::dim=3
+DEAL:1::n_cells: 147
+DEAL:1::FE_Nedelec<3>(0)
+DEAL:1::dim=3
+DEAL:1::n_cells: 140
+DEAL:1::FE_Nedelec<3>(1)
+DEAL:1::dim=3
+DEAL:1::n_cells: 161
+DEAL:1::FE_Q<3>(3)
+DEAL:1::dim=3
+DEAL:1::n_cells: 161
+DEAL:1::FE_DGQ<3>(2)
+DEAL:1::dim=3
+DEAL:1::n_cells: 182
+DEAL:1::FE_Q_DG0<3>(2)
+DEAL:1::dim=3
+DEAL:1::n_cells: 147
+