//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,
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
(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)
--- /dev/null
+// ---------------------------------------------------------------------
+// $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>();
+}