const hp::FECollection<DoFHandlerType::dimension,DoFHandlerType::space_dimension> fe_collection (dof.get_fe());
- // first, for each finite element, build a mask for each dof, not like
- // the one given which represents components. make sure we do the right
- // thing also with respect to non-primitive shape functions, which
- // takes some additional thought
- std::vector<Table<2,bool> > dof_mask(fe_collection.size());
-
- // check whether the table of couplings contains only true arguments,
- // i.e., we do not exclude any index. that is the easy case, since we
- // don't have to set up the tables
- bool need_dof_mask = false;
- for (unsigned int i=0; i<couplings.n_rows(); ++i)
- for (unsigned int j=0; j<couplings.n_cols(); ++j)
- if (couplings(i,j) == none)
- need_dof_mask = true;
-
- if (need_dof_mask == true)
- for (unsigned int f=0; f<fe_collection.size(); ++f)
- {
- const unsigned int dofs_per_cell = fe_collection[f].dofs_per_cell;
-
- dof_mask[f].reinit (dofs_per_cell, dofs_per_cell);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (fe_collection[f].is_primitive(i) &&
- fe_collection[f].is_primitive(j))
- dof_mask[f](i,j)
- = (couplings(fe_collection[f].system_to_component_index(i).first,
- fe_collection[f].system_to_component_index(j).first) != none);
- else
- {
- const unsigned int first_nonzero_comp_i
- = fe_collection[f].get_nonzero_components(i).first_selected_component();
- const unsigned int first_nonzero_comp_j
- = fe_collection[f].get_nonzero_components(j).first_selected_component();
- Assert (first_nonzero_comp_i < fe_collection[f].n_components(),
- ExcInternalError());
- Assert (first_nonzero_comp_j < fe_collection[f].n_components(),
- ExcInternalError());
-
- dof_mask[f](i,j)
- = (couplings(first_nonzero_comp_i,first_nonzero_comp_j) != none);
- }
- }
-
+ const std::vector<Table<2,Coupling> > dof_mask //(fe_collection.size())
+ = dof_couplings_from_component_couplings (fe_collection,
+ couplings);
+
+ // Convert the dof_mask to bool_dof_mask so we can pass it
+ // to constraints.add_entries_local_to_global()
+ std::vector<Table<2,bool> > bool_dof_mask (fe_collection.size());
+ for (unsigned int f=0; f<fe_collection.size(); ++f)
+ {
+ bool_dof_mask[f].reinit (TableIndices<2>(fe_collection[f].dofs_per_cell,fe_collection[f].dofs_per_cell));
+ bool_dof_mask[f].fill (false);
+ for (unsigned int i=0; i<fe_collection[f].dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe_collection[f].dofs_per_cell; ++j)
+ if (dof_mask[f](i,j) != none)
+ bool_dof_mask[f](i,j) = true;
+ }
std::vector<types::global_dof_index> dofs_on_this_cell(fe_collection.max_dofs_per_cell());
typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(),
constraints.add_entries_local_to_global (dofs_on_this_cell,
sparsity,
keep_constrained_dofs,
- dof_mask[fe_index]);
+ bool_dof_mask[fe_index]);
}
}
// non-hp DoFHandlers
template <typename DoFHandlerType, typename SparsityPatternType>
void
- make_flux_sparsity_pattern (const DoFHandlerType &dof,
- SparsityPatternType &sparsity,
- const Table<2,Coupling> &int_mask,
- const Table<2,Coupling> &flux_mask)
+ make_flux_sparsity_pattern (const DoFHandlerType &dof,
+ SparsityPatternType &sparsity,
+ const ConstraintMatrix &constraints,
+ const bool keep_constrained_dofs,
+ const Table<2,Coupling> &int_mask,
+ const Table<2,Coupling> &flux_mask,
+ const types::subdomain_id subdomain_id)
{
const FiniteElement<DoFHandlerType::dimension,DoFHandlerType::space_dimension>
&fe = dof.get_fe();
for (unsigned int f=0; f<GeometryInfo<DoFHandlerType::dimension>::faces_per_cell; ++f)
support_on_face(i,f) = fe.has_support_on_face(i,f);
+ // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
+ // to constraints.add_entries_local_to_global()
+ Table<2,bool> bool_int_dof_mask (fe.dofs_per_cell,fe.dofs_per_cell);
+ bool_int_dof_mask.fill (false);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ if (int_dof_mask(i,j) != none)
+ bool_int_dof_mask(i,j) = true;
+
typename DoFHandlerType::active_cell_iterator cell = dof.begin_active (),
endc = dof.end ();
for (; cell!=endc; ++cell)
- if (cell->is_locally_owned ())
+ if (((subdomain_id == numbers::invalid_subdomain_id)
+ ||
+ (subdomain_id == cell->subdomain_id()))
+ &&
+ cell->is_locally_owned())
{
cell->get_dof_indices (dofs_on_this_cell);
// make sparsity pattern for this cell
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- if (int_dof_mask (i,j) != none)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
-
+ constraints.add_entries_local_to_global (dofs_on_this_cell,
+ sparsity,
+ keep_constrained_dofs,
+ bool_int_dof_mask);
// Loop over all interior neighbors
for (unsigned int face_n = 0;
face_n < GeometryInfo<DoFHandlerType::dimension>::faces_per_cell;
void
make_flux_sparsity_pattern (const dealii::hp::DoFHandler<dim,spacedim> &dof,
SparsityPatternType &sparsity,
+ const ConstraintMatrix &constraints,
+ const bool keep_constrained_dofs,
const Table<2,Coupling> &int_mask,
- const Table<2,Coupling> &flux_mask)
+ const Table<2,Coupling> &flux_mask,
+ const types::subdomain_id subdomain_id)
{
// while the implementation above is quite optimized and caches a
// lot of data (see e.g. the int/flux_dof_mask tables), this is no
std::vector<types::global_dof_index> dofs_on_this_cell(DoFTools::max_dofs_per_cell(dof));
std::vector<types::global_dof_index> dofs_on_other_cell(DoFTools::max_dofs_per_cell(dof));
- const std::vector<Table<2,Coupling> >
- int_dof_mask
- = dof_couplings_from_component_couplings(fe, int_mask);
+ std::vector<Table<2,Coupling> > int_dof_mask (fe.size());
+
+ int_dof_mask = dof_couplings_from_component_couplings (fe, int_mask);
+
+ // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
+ // to constraints.add_entries_local_to_global()
+ std::vector<Table<2,bool> > bool_int_dof_mask (fe.size());
+ for (unsigned int f=0; f<fe.size(); ++f)
+ {
+ bool_int_dof_mask[f].reinit (TableIndices<2>(fe[f].dofs_per_cell,fe[f].dofs_per_cell));
+ bool_int_dof_mask[f].fill (false);
+ for (unsigned int i=0; i<fe[f].dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe[f].dofs_per_cell; ++j)
+ if (int_dof_mask[f](i,j) != none)
+ bool_int_dof_mask[f](i,j) = true;
+ }
+
typename dealii::hp::DoFHandler<dim,spacedim>::active_cell_iterator
cell = dof.begin_active(),
endc = dof.end();
for (; cell!=endc; ++cell)
- {
- dofs_on_this_cell.resize (cell->get_fe().dofs_per_cell);
- cell->get_dof_indices (dofs_on_this_cell);
+ if (((subdomain_id == numbers::invalid_subdomain_id)
+ ||
+ (subdomain_id == cell->subdomain_id()))
+ &&
+ cell->is_locally_owned())
+ {
+ dofs_on_this_cell.resize (cell->get_fe().dofs_per_cell);
+ cell->get_dof_indices (dofs_on_this_cell);
- // make sparsity pattern for this cell
- for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- for (unsigned int j=0; j<cell->get_fe().dofs_per_cell; ++j)
- if (int_dof_mask[cell->active_fe_index()](i,j) != none)
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
-
- // Loop over all interior neighbors
- for (unsigned int face = 0;
- face < GeometryInfo<dim>::faces_per_cell;
- ++face)
- {
- const typename dealii::hp::DoFHandler<dim,spacedim>::face_iterator
- cell_face = cell->face(face);
- if (cell_face->user_flag_set ())
- continue;
-
- const bool periodic_neighbor = cell->has_periodic_neighbor (face);
-
- if (cell->at_boundary (face) && (!periodic_neighbor))
- {
- for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- for (unsigned int j=0; j<cell->get_fe().dofs_per_cell; ++j)
- if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
- cell->get_fe().system_to_component_index(j).first)
- == always)
- ||
- (flux_mask(cell->get_fe().system_to_component_index(i).first,
- cell->get_fe().system_to_component_index(j).first)
- == nonzero))
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- }
- else
- {
- typename dealii::hp::DoFHandler<dim,spacedim>::level_cell_iterator
- neighbor = cell->neighbor_or_periodic_neighbor(face);
-
- // Refinement edges are taken care of by coarser cells
- if ((!periodic_neighbor && cell->neighbor_is_coarser(face)) ||
- (periodic_neighbor && cell->periodic_neighbor_is_coarser(face)))
- continue;
-
- const unsigned int
- neighbor_face = periodic_neighbor?
- cell->periodic_neighbor_of_periodic_neighbor(face):
- cell->neighbor_of_neighbor(face);
-
- if (cell_face->has_children())
- {
- for (unsigned int sub_nr = 0;
- sub_nr != cell_face->n_children();
- ++sub_nr)
- {
- const typename dealii::hp::DoFHandler<dim,spacedim>::level_cell_iterator
- sub_neighbor
- = periodic_neighbor?
- cell->periodic_neighbor_child_on_subface (face, sub_nr):
- cell->neighbor_child_on_subface (face, sub_nr);
-
- dofs_on_other_cell.resize (sub_neighbor->get_fe().dofs_per_cell);
- sub_neighbor->get_dof_indices (dofs_on_other_cell);
- for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<sub_neighbor->get_fe().dofs_per_cell;
- ++j)
- {
- if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
- sub_neighbor->get_fe().system_to_component_index(j).first)
- == always)
- ||
- (flux_mask(cell->get_fe().system_to_component_index(i).first,
- sub_neighbor->get_fe().system_to_component_index(j).first)
- == nonzero))
- {
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_other_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_other_cell[j]);
- }
-
- if ((flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first,
- cell->get_fe().system_to_component_index(i).first)
- == always)
- ||
- (flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first,
- cell->get_fe().system_to_component_index(i).first)
- == nonzero))
- {
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_other_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_other_cell[i]);
- }
- }
- }
- sub_neighbor->face(neighbor_face)->set_user_flag ();
- }
- }
- else
- {
- dofs_on_other_cell.resize (neighbor->get_fe().dofs_per_cell);
- neighbor->get_dof_indices (dofs_on_other_cell);
- for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<neighbor->get_fe().dofs_per_cell; ++j)
- {
- if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
- neighbor->get_fe().system_to_component_index(j).first)
- == always)
- ||
- (flux_mask(cell->get_fe().system_to_component_index(i).first,
- neighbor->get_fe().system_to_component_index(j).first)
- == nonzero))
- {
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_other_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_this_cell[i],
- dofs_on_this_cell[j]);
- sparsity.add (dofs_on_other_cell[i],
- dofs_on_other_cell[j]);
- }
-
- if ((flux_mask(neighbor->get_fe().system_to_component_index(j).first,
- cell->get_fe().system_to_component_index(i).first)
- == always)
- ||
- (flux_mask(neighbor->get_fe().system_to_component_index(j).first,
- cell->get_fe().system_to_component_index(i).first)
- == nonzero))
- {
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_other_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_this_cell[j],
- dofs_on_this_cell[i]);
- sparsity.add (dofs_on_other_cell[j],
- dofs_on_other_cell[i]);
- }
- }
- }
- neighbor->face(neighbor_face)->set_user_flag ();
- }
- }
- }
- }
+ // make sparsity pattern for this cell
+ constraints.add_entries_local_to_global (dofs_on_this_cell,
+ sparsity,
+ keep_constrained_dofs,
+ bool_int_dof_mask[cell->active_fe_index()]);
+ // Loop over all interior neighbors
+ for (unsigned int face = 0;
+ face < GeometryInfo<dim>::faces_per_cell;
+ ++face)
+ {
+ const typename dealii::hp::DoFHandler<dim,spacedim>::face_iterator
+ cell_face = cell->face(face);
+ if (cell_face->user_flag_set ())
+ continue;
+
+ const bool periodic_neighbor = cell->has_periodic_neighbor (face);
+
+ if (cell->at_boundary (face) && (!periodic_neighbor))
+ {
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+ for (unsigned int j=0; j<cell->get_fe().dofs_per_cell; ++j)
+ if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
+ cell->get_fe().system_to_component_index(j).first)
+ == always)
+ ||
+ (flux_mask(cell->get_fe().system_to_component_index(i).first,
+ cell->get_fe().system_to_component_index(j).first)
+ == nonzero))
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ }
+ else
+ {
+ typename dealii::hp::DoFHandler<dim,spacedim>::level_cell_iterator
+ neighbor = cell->neighbor_or_periodic_neighbor(face);
+
+ // Refinement edges are taken care of by coarser cells
+ if ((!periodic_neighbor && cell->neighbor_is_coarser(face)) ||
+ (periodic_neighbor && cell->periodic_neighbor_is_coarser(face)))
+ continue;
+
+ const unsigned int
+ neighbor_face = periodic_neighbor?
+ cell->periodic_neighbor_of_periodic_neighbor(face):
+ cell->neighbor_of_neighbor(face);
+
+ if (cell_face->has_children())
+ {
+ for (unsigned int sub_nr = 0;
+ sub_nr != cell_face->n_children();
+ ++sub_nr)
+ {
+ const typename dealii::hp::DoFHandler<dim,spacedim>::level_cell_iterator
+ sub_neighbor
+ = periodic_neighbor?
+ cell->periodic_neighbor_child_on_subface (face, sub_nr):
+ cell->neighbor_child_on_subface (face, sub_nr);
+
+ dofs_on_other_cell.resize (sub_neighbor->get_fe().dofs_per_cell);
+ sub_neighbor->get_dof_indices (dofs_on_other_cell);
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<sub_neighbor->get_fe().dofs_per_cell;
+ ++j)
+ {
+ if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
+ sub_neighbor->get_fe().system_to_component_index(j).first)
+ == always)
+ ||
+ (flux_mask(cell->get_fe().system_to_component_index(i).first,
+ sub_neighbor->get_fe().system_to_component_index(j).first)
+ == nonzero))
+ {
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_other_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_other_cell[j]);
+ }
+
+ if ((flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first,
+ cell->get_fe().system_to_component_index(i).first)
+ == always)
+ ||
+ (flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first,
+ cell->get_fe().system_to_component_index(i).first)
+ == nonzero))
+ {
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_other_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_other_cell[i]);
+ }
+ }
+ }
+ sub_neighbor->face(neighbor_face)->set_user_flag ();
+ }
+ }
+ else
+ {
+ dofs_on_other_cell.resize (neighbor->get_fe().dofs_per_cell);
+ neighbor->get_dof_indices (dofs_on_other_cell);
+ for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<neighbor->get_fe().dofs_per_cell; ++j)
+ {
+ if ((flux_mask(cell->get_fe().system_to_component_index(i).first,
+ neighbor->get_fe().system_to_component_index(j).first)
+ == always)
+ ||
+ (flux_mask(cell->get_fe().system_to_component_index(i).first,
+ neighbor->get_fe().system_to_component_index(j).first)
+ == nonzero))
+ {
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_other_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_this_cell[i],
+ dofs_on_this_cell[j]);
+ sparsity.add (dofs_on_other_cell[i],
+ dofs_on_other_cell[j]);
+ }
+
+ if ((flux_mask(neighbor->get_fe().system_to_component_index(j).first,
+ cell->get_fe().system_to_component_index(i).first)
+ == always)
+ ||
+ (flux_mask(neighbor->get_fe().system_to_component_index(j).first,
+ cell->get_fe().system_to_component_index(i).first)
+ == nonzero))
+ {
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_other_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_this_cell[j],
+ dofs_on_this_cell[i]);
+ sparsity.add (dofs_on_other_cell[j],
+ dofs_on_other_cell[i]);
+ }
+ }
+ }
+ neighbor->face(neighbor_face)->set_user_flag ();
+ }
+ }
+ }
+ }
}
}
template <typename DoFHandlerType, typename SparsityPatternType>
void
- make_flux_sparsity_pattern (const DoFHandlerType &dof,
- SparsityPatternType &sparsity,
- const Table<2,Coupling> &int_mask,
- const Table<2,Coupling> &flux_mask)
+ make_flux_sparsity_pattern (const DoFHandlerType &dof,
+ SparsityPatternType &sparsity,
+ const Table<2,Coupling> &int_mask,
+ const Table<2,Coupling> &flux_mask,
+ const types::subdomain_id subdomain_id)
+ {
+ ConstraintMatrix constraints;
+ const bool keep_constrained_dofs = true;
+ make_flux_sparsity_pattern (dof, sparsity,
+ constraints, keep_constrained_dofs,
+ int_mask, flux_mask,
+ subdomain_id);
+ }
+
+ template <typename DoFHandlerType, typename SparsityPatternType>
+ void
+ make_flux_sparsity_pattern (const DoFHandlerType &dof,
+ SparsityPatternType &sparsity,
+ const ConstraintMatrix &constraints,
+ const bool keep_constrained_dofs,
+ const Table<2,Coupling> &int_mask,
+ const Table<2,Coupling> &flux_mask,
+ const types::subdomain_id subdomain_id)
{
// do the error checking and frame code here, and then pass on to more
// specialized functions in the internal namespace
Assert (flux_mask.n_cols() == n_comp,
ExcDimensionMismatch (flux_mask.n_cols(), n_comp));
+ // If we have a distributed::Triangulation only allow locally_owned
+ // subdomain. Not setting a subdomain is also okay, because we skip
+ // ghost cells in the loop below.
+ Assert (
+ (dof.get_triangulation().locally_owned_subdomain() == numbers::invalid_subdomain_id)
+ ||
+ (subdomain_id == numbers::invalid_subdomain_id)
+ ||
+ (subdomain_id == dof.get_triangulation().locally_owned_subdomain()),
+ ExcMessage ("For parallel::distributed::Triangulation objects and "
+ "associated DoF handler objects, asking for any subdomain other "
+ "than the locally owned one does not make sense."));
+
// Clear user flags because we will need them. But first we save them
// and make sure that we restore them later such that at the end of
// this function the Triangulation will be in the same state as it was
(dof.get_triangulation()).clear_user_flags ();
internal::make_flux_sparsity_pattern (dof, sparsity,
- int_mask, flux_mask);
+ constraints, keep_constrained_dofs,
+ int_mask, flux_mask,
+ subdomain_id);
// finally restore the user flags
const_cast<Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &>
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+//
+// This test was written by Sam Cox.
+
+// Tests the behaviour of DoFTools::make_flux_sparsity_pattern (DoFHandler,
+// SparsityPattern, ConstraintMatrix, bool,
+// coupling, flux_coupling, subdomain_id)
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+#include <list>
+#include <set>
+#include <cstdio>
+
+
+void test ()
+{
+
+ unsigned int np = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+ unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+ const int dim = 2;
+ // Setup system
+ dealii::parallel::distributed::Triangulation<dim> triangulation(MPI_COMM_WORLD);
+
+ GridGenerator::hyper_rectangle (triangulation,
+ Point<dim>(0,0),
+ Point<dim>(1,1),
+ true);
+
+ triangulation.refine_global(1);
+ // Extra refinement to generate hanging nodes
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ cell != triangulation.end(); ++cell)
+ if (cell->center()(0) > 0.49)
+ cell->set_refine_flag ();
+
+ triangulation.prepare_coarsening_and_refinement();
+ triangulation.execute_coarsening_and_refinement ();
+ triangulation.repartition();
+
+ const FESystem<dim> fe_system (FE_Q<dim>(2),1,
+ FE_DGQ<dim>(2),1);
+
+ DoFHandler<dim> dh (triangulation);
+
+ dh.distribute_dofs (fe_system);
+
+ // Couple the internal DoFs of both finite elements.
+ // Only couple the face terms of the discontinuous element.
+ Table<2,DoFTools::Coupling> coupling (2, 2);
+ Table<2,DoFTools::Coupling> flux_coupling (2, 2);
+
+ for (unsigned int i=0; i<2; ++i)
+ for (unsigned int j=0; j<2; ++j)
+ {
+ coupling[i][j] = DoFTools::none;
+ flux_coupling[i][j] = DoFTools::none;
+ }
+
+ coupling[0][0] = DoFTools::always;
+ coupling[1][1] = DoFTools::always;
+ flux_coupling[1][1] = DoFTools::always;
+
+ IndexSet relevant_partitioning (dh.n_dofs());
+ DoFTools::extract_locally_relevant_dofs (dh, relevant_partitioning);
+
+ // Generate hanging node constraints
+ ConstraintMatrix constraints;
+ constraints.clear();
+ DoFTools::make_hanging_node_constraints (dh,
+ constraints);
+ constraints.close();
+
+ // Generate sparsity pattern
+ DynamicSparsityPattern sp (relevant_partitioning);
+ DoFTools::make_flux_sparsity_pattern (dh, sp,
+ constraints, false,
+ coupling,
+ flux_coupling,
+ Utilities::MPI::this_mpi_process(MPI_COMM_WORLD));
+ SparsityTools::distribute_sparsity_pattern (sp,
+ dh.n_locally_owned_dofs_per_processor(),
+ MPI_COMM_WORLD,
+ relevant_partitioning);
+
+ // Output
+ MPI_Barrier(MPI_COMM_WORLD);
+
+ deallog.push(Utilities::int_to_string(myid));
+
+ deallog << "**** proc " << myid << ": \n\n";
+ deallog << "Sparsity pattern:\n";
+ sp.print_gnuplot(deallog.get_file_stream());
+
+ MPI_Barrier(MPI_COMM_WORLD);
+
+}
+
+int main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads());
+ const unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+ MPILogInitAll log;
+
+ deallog.push(Utilities::int_to_string(myid));
+
+ test();
+ deallog.pop();
+}