]> https://gitweb.dealii.org/ - dealii.git/commitdiff
take over #2273 from tjhei/split_instantiations 2277/head
authorMatthias Maier <tamiko@43-1.org>
Mon, 29 Feb 2016 19:41:43 +0000 (13:41 -0600)
committerTimo Heister <timo.heister@gmail.com>
Mon, 29 Feb 2016 21:06:37 +0000 (16:06 -0500)
Split instantiations of matrix_tools.cc
Conflicts:
source/numerics/matrix_tools.cc
source/numerics/matrix_tools.inst.in

include/deal.II/numerics/matrix_creator.templates.h [new file with mode: 0644]
source/numerics/CMakeLists.txt
source/numerics/matrix_creator.cc [new file with mode: 0644]
source/numerics/matrix_creator.inst.in [new file with mode: 0644]
source/numerics/matrix_creator_inst2.cc [new file with mode: 0644]
source/numerics/matrix_creator_inst3.cc [new file with mode: 0644]
source/numerics/matrix_tools.cc
source/numerics/matrix_tools.inst.in
source/numerics/matrix_tools_once.cc [new file with mode: 0644]

diff --git a/include/deal.II/numerics/matrix_creator.templates.h b/include/deal.II/numerics/matrix_creator.templates.h
new file mode 100644 (file)
index 0000000..30953c1
--- /dev/null
@@ -0,0 +1,1894 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__matrix_creator_templates_h
+#define dealii__matrix_creator_templates_h
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+
+#ifdef DEAL_II_WITH_PETSC
+#  include <deal.II/lac/petsc_parallel_sparse_matrix.h>
+#  include <deal.II/lac/petsc_sparse_matrix.h>
+#  include <deal.II/lac/petsc_parallel_vector.h>
+#  include <deal.II/lac/petsc_vector.h>
+#  include <deal.II/lac/petsc_parallel_block_sparse_matrix.h>
+#endif
+
+#ifdef DEAL_II_WITH_TRILINOS
+#  include <deal.II/lac/trilinos_sparse_matrix.h>
+#  include <deal.II/lac/trilinos_vector.h>
+#  include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#  include <deal.II/lac/trilinos_block_vector.h>
+#endif
+
+#include <algorithm>
+
+
+#include <algorithm>
+#include <set>
+#include <cmath>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace MatrixCreator
+{
+  namespace internal
+  {
+    /**
+     * Convenience abbreviation for
+     * pairs of DoF handler cell
+     * iterators. This type works
+     * just like a
+     * <tt>std::pair<iterator,iterator></tt>
+     * but is templatized on the
+     * dof handler that should be used.
+     */
+    template <typename DoFHandlerType>
+    struct IteratorRange
+    {
+      /**
+       * Typedef for the iterator type.
+       */
+      typedef typename DoFHandlerType::active_cell_iterator active_cell_iterator;
+
+      /**
+       * Abbreviation for a pair of
+       * iterators.
+       */
+      typedef std::pair<active_cell_iterator,active_cell_iterator> iterator_pair;
+
+      /**
+       * Constructor. Initialize
+       * the two values by the
+       * given values.
+       */
+      IteratorRange (const active_cell_iterator &first,
+                     const active_cell_iterator &second);
+
+      /**
+       * Constructor taking a pair
+       * of values for
+       * initialization.
+       */
+      IteratorRange (const iterator_pair &ip);
+
+      /**
+       * Pair of iterators denoting
+       * a half-open range.
+       */
+      active_cell_iterator first, second;
+    };
+
+
+
+
+    template <typename DoFHandlerType>
+    inline
+    IteratorRange<DoFHandlerType>::
+    IteratorRange (const active_cell_iterator &first,
+                   const active_cell_iterator &second)
+      :
+      first (first),
+      second (second)
+    {}
+
+
+
+    template <typename DoFHandlerType>
+    inline
+    IteratorRange<DoFHandlerType>::IteratorRange (const iterator_pair &ip)
+      :
+      first (ip.first),
+      second (ip.second)
+    {}
+
+
+
+    namespace AssemblerData
+    {
+      template <int dim,
+                int spacedim>
+      struct Scratch
+      {
+        Scratch (const ::dealii::hp::FECollection<dim,spacedim> &fe,
+                 const UpdateFlags         update_flags,
+                 const Function<spacedim> *coefficient,
+                 const Function<spacedim> *rhs_function,
+                 const ::dealii::hp::QCollection<dim> &quadrature,
+                 const ::dealii::hp::MappingCollection<dim,spacedim> &mapping)
+          :
+          fe_collection (fe),
+          quadrature_collection (quadrature),
+          mapping_collection (mapping),
+          x_fe_values (mapping_collection,
+                       fe_collection,
+                       quadrature_collection,
+                       update_flags),
+          coefficient_values(quadrature_collection.max_n_quadrature_points()),
+          coefficient_vector_values (quadrature_collection.max_n_quadrature_points(),
+                                     dealii::Vector<double> (fe_collection.n_components())),
+          rhs_values(quadrature_collection.max_n_quadrature_points()),
+          rhs_vector_values(quadrature_collection.max_n_quadrature_points(),
+                            dealii::Vector<double> (fe_collection.n_components())),
+          coefficient (coefficient),
+          rhs_function (rhs_function),
+          update_flags (update_flags)
+        {}
+
+        Scratch (const Scratch &data)
+          :
+          fe_collection (data.fe_collection),
+          quadrature_collection (data.quadrature_collection),
+          mapping_collection (data.mapping_collection),
+          x_fe_values (mapping_collection,
+                       fe_collection,
+                       quadrature_collection,
+                       data.update_flags),
+          coefficient_values (data.coefficient_values),
+          coefficient_vector_values (data.coefficient_vector_values),
+          rhs_values (data.rhs_values),
+          rhs_vector_values (data.rhs_vector_values),
+          coefficient (data.coefficient),
+          rhs_function (data.rhs_function),
+          update_flags (data.update_flags)
+        {}
+
+        const ::dealii::hp::FECollection<dim,spacedim>      &fe_collection;
+        const ::dealii::hp::QCollection<dim>                &quadrature_collection;
+        const ::dealii::hp::MappingCollection<dim,spacedim> &mapping_collection;
+
+        ::dealii::hp::FEValues<dim,spacedim> x_fe_values;
+
+        std::vector<double>                  coefficient_values;
+        std::vector<dealii::Vector<double> > coefficient_vector_values;
+        std::vector<double>                  rhs_values;
+        std::vector<dealii::Vector<double> > rhs_vector_values;
+
+        const Function<spacedim>   *coefficient;
+        const Function<spacedim>   *rhs_function;
+
+        const UpdateFlags update_flags;
+      };
+
+
+      template <typename number>
+      struct CopyData
+      {
+        std::vector<types::global_dof_index> dof_indices;
+        FullMatrix<number>        cell_matrix;
+        dealii::Vector<number>    cell_rhs;
+        const ConstraintMatrix   *constraints;
+      };
+    }
+
+
+    template <int dim,
+              int spacedim,
+              typename CellIterator,
+              typename number>
+    void mass_assembler (const CellIterator &cell,
+                         MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim> &data,
+                         MatrixCreator::internal::AssemblerData::CopyData<number>      &copy_data)
+    {
+      data.x_fe_values.reinit (cell);
+      const FEValues<dim,spacedim> &fe_values = data.x_fe_values.get_present_fe_values ();
+
+      const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                         n_q_points    = fe_values.n_quadrature_points;
+      const FiniteElement<dim,spacedim> &fe  = fe_values.get_fe();
+      const unsigned int n_components  = fe.n_components();
+
+      Assert(data.rhs_function == 0 ||
+             data.rhs_function->n_components==1 ||
+             data.rhs_function->n_components==n_components,
+             ::dealii::MatrixCreator::ExcComponentMismatch());
+      Assert(data.coefficient == 0 ||
+             data.coefficient->n_components==1 ||
+             data.coefficient->n_components==n_components,
+             ::dealii::MatrixCreator::ExcComponentMismatch());
+
+      copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+      copy_data.cell_rhs.reinit (dofs_per_cell);
+
+      copy_data.dof_indices.resize (dofs_per_cell);
+      cell->get_dof_indices (copy_data.dof_indices);
+
+      const bool use_rhs_function = data.rhs_function != 0;
+      if (use_rhs_function)
+        {
+          if (data.rhs_function->n_components==1)
+            {
+              data.rhs_values.resize (n_q_points);
+              data.rhs_function->value_list (fe_values.get_quadrature_points(),
+                                             data.rhs_values);
+            }
+          else
+            {
+              data.rhs_vector_values.resize (n_q_points,
+                                             dealii::Vector<double>(n_components));
+              data.rhs_function->vector_value_list (fe_values.get_quadrature_points(),
+                                                    data.rhs_vector_values);
+            }
+        }
+
+      const bool use_coefficient = data.coefficient != 0;
+      if (use_coefficient)
+        {
+          if (data.coefficient->n_components==1)
+            {
+              data.coefficient_values.resize (n_q_points);
+              data.coefficient->value_list (fe_values.get_quadrature_points(),
+                                            data.coefficient_values);
+            }
+          else
+            {
+              data.coefficient_vector_values.resize (n_q_points,
+                                                     dealii::Vector<double>(n_components));
+              data.coefficient->vector_value_list (fe_values.get_quadrature_points(),
+                                                   data.coefficient_vector_values);
+            }
+        }
+
+
+      double add_data;
+      const std::vector<double> &JxW = fe_values.get_JxW_values();
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        if (fe.is_primitive ())
+          {
+            const unsigned int component_i =
+              fe.system_to_component_index(i).first;
+            const double *phi_i = &fe_values.shape_value(i,0);
+            add_data = 0;
+
+            // use symmetry in the mass matrix here:
+            // just need to calculate the diagonal
+            // and half of the elements above the
+            // diagonal
+            for (unsigned int j=i; j<dofs_per_cell; ++j)
+              if ((n_components==1) ||
+                  (fe.system_to_component_index(j).first ==
+                   component_i))
+                {
+                  const double *phi_j = &fe_values.shape_value(j,0);
+                  add_data = 0;
+                  if (use_coefficient)
+                    {
+                      if (data.coefficient->n_components==1)
+                        for (unsigned int point=0; point<n_q_points; ++point)
+                          add_data += (phi_i[point] * phi_j[point] * JxW[point] *
+                                       data.coefficient_values[point]);
+                      else
+                        for (unsigned int point=0; point<n_q_points; ++point)
+                          add_data += (phi_i[point] * phi_j[point] * JxW[point] *
+                                       data.coefficient_vector_values[point](component_i));
+                    }
+                  else
+                    for (unsigned int point=0; point<n_q_points; ++point)
+                      add_data += phi_i[point] * phi_j[point] * JxW[point];
+
+                  // this is even ok for i==j, since then
+                  // we just write the same value twice.
+                  copy_data.cell_matrix(i,j) = add_data;
+                  copy_data.cell_matrix(j,i) = add_data;
+                }
+
+            if (use_rhs_function)
+              {
+                add_data = 0;
+                if (data.rhs_function->n_components==1)
+                  for (unsigned int point=0; point<n_q_points; ++point)
+                    add_data += phi_i[point] * JxW[point] *
+                                data.rhs_values[point];
+                else
+                  for (unsigned int point=0; point<n_q_points; ++point)
+                    add_data += phi_i[point] * JxW[point] *
+                                data.rhs_vector_values[point](component_i);
+                copy_data.cell_rhs(i) = add_data;
+              }
+          }
+        else
+          {
+            // non-primitive vector-valued FE, using
+            // symmetry again
+            for (unsigned int j=i; j<dofs_per_cell; ++j)
+              {
+                add_data = 0;
+                for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                  if (fe.get_nonzero_components(i)[comp_i] &&
+                      fe.get_nonzero_components(j)[comp_i])
+                    {
+                      if (use_coefficient)
+                        {
+                          if (data.coefficient->n_components==1)
+                            for (unsigned int point=0; point<n_q_points; ++point)
+                              add_data += (fe_values.shape_value_component(i,point,comp_i) *
+                                           fe_values.shape_value_component(j,point,comp_i) *
+                                           JxW[point] *
+                                           data.coefficient_values[point]);
+                          else
+                            for (unsigned int point=0; point<n_q_points; ++point)
+                              add_data += (fe_values.shape_value_component(i,point,comp_i) *
+                                           fe_values.shape_value_component(j,point,comp_i) *
+                                           JxW[point] *
+                                           data.coefficient_vector_values[point](comp_i));
+                        }
+                      else
+                        for (unsigned int point=0; point<n_q_points; ++point)
+                          add_data += fe_values.shape_value_component(i,point,comp_i) *
+                                      fe_values.shape_value_component(j,point,comp_i) * JxW[point];
+                    }
+
+                copy_data.cell_matrix(i,j) = add_data;
+                copy_data.cell_matrix(j,i) = add_data;
+              }
+
+            if (use_rhs_function)
+              {
+                add_data = 0;
+                for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                  if (fe.get_nonzero_components(i)[comp_i])
+                    {
+                      if (data.rhs_function->n_components==1)
+                        for (unsigned int point=0; point<n_q_points; ++point)
+                          add_data += fe_values.shape_value_component(i,point,comp_i) *
+                                      JxW[point] * data.rhs_values[point];
+                      else
+                        for (unsigned int point=0; point<n_q_points; ++point)
+                          add_data += fe_values.shape_value_component(i,point,comp_i) *
+                                      JxW[point] * data.rhs_vector_values[point](comp_i);
+                    }
+                copy_data.cell_rhs(i) = add_data;
+              }
+          }
+    }
+
+
+
+    template <int dim,
+              int spacedim,
+              typename CellIterator>
+    void laplace_assembler (const CellIterator &cell,
+                            MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim> &data,
+                            MatrixCreator::internal::AssemblerData::CopyData<double>      &copy_data)
+    {
+      data.x_fe_values.reinit (cell);
+      const FEValues<dim,spacedim> &fe_values = data.x_fe_values.get_present_fe_values ();
+
+      const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                         n_q_points    = fe_values.n_quadrature_points;
+      const FiniteElement<dim,spacedim>    &fe  = fe_values.get_fe();
+      const unsigned int n_components  = fe.n_components();
+
+      Assert(data.rhs_function == 0 ||
+             data.rhs_function->n_components==1 ||
+             data.rhs_function->n_components==n_components,
+             ::dealii::MatrixCreator::ExcComponentMismatch());
+      Assert(data.coefficient == 0 ||
+             data.coefficient->n_components==1 ||
+             data.coefficient->n_components==n_components,
+             ::dealii::MatrixCreator::ExcComponentMismatch());
+
+      copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+      copy_data.cell_rhs.reinit (dofs_per_cell);
+      copy_data.dof_indices.resize (dofs_per_cell);
+      cell->get_dof_indices (copy_data.dof_indices);
+
+
+      const bool use_rhs_function = data.rhs_function != 0;
+      if (use_rhs_function)
+        {
+          if (data.rhs_function->n_components==1)
+            {
+              data.rhs_values.resize (n_q_points);
+              data.rhs_function->value_list (fe_values.get_quadrature_points(),
+                                             data.rhs_values);
+            }
+          else
+            {
+              data.rhs_vector_values.resize (n_q_points,
+                                             dealii::Vector<double>(n_components));
+              data.rhs_function->vector_value_list (fe_values.get_quadrature_points(),
+                                                    data.rhs_vector_values);
+            }
+        }
+
+      const bool use_coefficient = data.coefficient != 0;
+      if (use_coefficient)
+        {
+          if (data.coefficient->n_components==1)
+            {
+              data.coefficient_values.resize (n_q_points);
+              data.coefficient->value_list (fe_values.get_quadrature_points(),
+                                            data.coefficient_values);
+            }
+          else
+            {
+              data.coefficient_vector_values.resize (n_q_points,
+                                                     dealii::Vector<double>(n_components));
+              data.coefficient->vector_value_list (fe_values.get_quadrature_points(),
+                                                   data.coefficient_vector_values);
+            }
+        }
+
+
+      const std::vector<double> &JxW = fe_values.get_JxW_values();
+      double add_data;
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+        if (fe.is_primitive ())
+          {
+            const unsigned int component_i =
+              fe.system_to_component_index(i).first;
+            const Tensor<1,spacedim> *grad_phi_i =
+              &fe_values.shape_grad(i,0);
+
+            // can use symmetry
+            for (unsigned int j=i; j<dofs_per_cell; ++j)
+              if ((n_components==1) ||
+                  (fe.system_to_component_index(j).first ==
+                   component_i))
+                {
+                  const Tensor<1,spacedim> *grad_phi_j =
+                    & fe_values.shape_grad(j,0);
+                  add_data = 0;
+                  if (use_coefficient)
+                    {
+                      if (data.coefficient->n_components==1)
+                        for (unsigned int point=0; point<n_q_points; ++point)
+                          add_data += ((grad_phi_i[point]*grad_phi_j[point]) *
+                                       JxW[point] *
+                                       data.coefficient_values[point]);
+                      else
+                        for (unsigned int point=0; point<n_q_points; ++point)
+                          add_data += ((grad_phi_i[point]*grad_phi_j[point]) *
+                                       JxW[point] *
+                                       data.coefficient_vector_values[point](component_i));
+                    }
+                  else
+                    for (unsigned int point=0; point<n_q_points; ++point)
+                      add_data += (grad_phi_i[point]*grad_phi_j[point]) *
+                                  JxW[point];
+
+                  copy_data.cell_matrix(i,j) = add_data;
+                  copy_data.cell_matrix(j,i) = add_data;
+                }
+
+            if (use_rhs_function)
+              {
+                const double *phi_i = &fe_values.shape_value(i,0);
+                add_data = 0;
+                if (data.rhs_function->n_components==1)
+                  for (unsigned int point=0; point<n_q_points; ++point)
+                    add_data += phi_i[point] * JxW[point] *
+                                data.rhs_values[point];
+                else
+                  for (unsigned int point=0; point<n_q_points; ++point)
+                    add_data += phi_i[point] * JxW[point] *
+                                data.rhs_vector_values[point](component_i);
+                copy_data.cell_rhs(i) = add_data;
+              }
+          }
+        else
+          {
+            // non-primitive vector-valued FE
+            for (unsigned int j=i; j<dofs_per_cell; ++j)
+              {
+                add_data = 0;
+                for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                  if (fe.get_nonzero_components(i)[comp_i] &&
+                      fe.get_nonzero_components(j)[comp_i])
+                    {
+                      if (use_coefficient)
+                        {
+                          if (data.coefficient->n_components==1)
+                            for (unsigned int point=0; point<n_q_points; ++point)
+                              add_data += ((fe_values.shape_grad_component(i,point,comp_i) *
+                                            fe_values.shape_grad_component(j,point,comp_i)) *
+                                           JxW[point] *
+                                           data.coefficient_values[point]);
+                          else
+                            for (unsigned int point=0; point<n_q_points; ++point)
+                              add_data += ((fe_values.shape_grad_component(i,point,comp_i) *
+                                            fe_values.shape_grad_component(j,point,comp_i)) *
+                                           JxW[point] *
+                                           data.coefficient_vector_values[point](comp_i));
+                        }
+                      else
+                        for (unsigned int point=0; point<n_q_points; ++point)
+                          add_data += (fe_values.shape_grad_component(i,point,comp_i) *
+                                       fe_values.shape_grad_component(j,point,comp_i)) *
+                                      JxW[point];
+                    }
+
+                copy_data.cell_matrix(i,j) = add_data;
+                copy_data.cell_matrix(j,i) = add_data;
+              }
+
+            if (use_rhs_function)
+              {
+                add_data = 0;
+                for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                  if (fe.get_nonzero_components(i)[comp_i])
+                    {
+                      if (data.rhs_function->n_components==1)
+                        for (unsigned int point=0; point<n_q_points; ++point)
+                          add_data += fe_values.shape_value_component(i,point,comp_i) *
+                                      JxW[point] * data.rhs_values[point];
+                      else
+                        for (unsigned int point=0; point<n_q_points; ++point)
+                          add_data += fe_values.shape_value_component(i,point,comp_i) *
+                                      JxW[point] * data.rhs_vector_values[point](comp_i);
+                    }
+                copy_data.cell_rhs(i) = add_data;
+              }
+          }
+    }
+
+
+
+    template <typename number,
+              typename MatrixType,
+              typename VectorType>
+    void copy_local_to_global (const AssemblerData::CopyData<number> &data,
+                               MatrixType *matrix,
+                               VectorType *right_hand_side)
+    {
+      const unsigned int dofs_per_cell = data.dof_indices.size();
+      (void)dofs_per_cell;
+
+      Assert (data.cell_matrix.m() == dofs_per_cell,
+              ExcInternalError());
+      Assert (data.cell_matrix.n() == dofs_per_cell,
+              ExcInternalError());
+      Assert ((right_hand_side == 0)
+              ||
+              (data.cell_rhs.size() == dofs_per_cell),
+              ExcInternalError());
+
+      if (right_hand_side != 0)
+        data.constraints->distribute_local_to_global(data.cell_matrix,
+                                                     data.cell_rhs,
+                                                     data.dof_indices,
+                                                     *matrix, *right_hand_side);
+      else
+        data.constraints->distribute_local_to_global(data.cell_matrix,
+                                                     data.dof_indices,
+                                                     *matrix);
+    }
+
+
+
+    namespace AssemblerBoundary
+    {
+      struct Scratch
+      {
+        Scratch() {}
+      };
+
+      template <typename DoFHandlerType>
+      struct CopyData
+      {
+        CopyData() {};
+
+        CopyData(CopyData const &data);
+
+        unsigned int dofs_per_cell;
+        std::vector<types::global_dof_index> dofs;
+        std::vector<std::vector<bool> > dof_is_on_face;
+        typename DoFHandlerType::active_cell_iterator cell;
+        std::vector<FullMatrix<double> > cell_matrix;
+        std::vector<Vector<double> > cell_vector;
+      };
+
+      template <typename DoFHandlerType>
+      CopyData<DoFHandlerType>::CopyData(CopyData const &data) :
+        dofs_per_cell(data.dofs_per_cell),
+        dofs(data.dofs),
+        dof_is_on_face(data.dof_is_on_face),
+        cell(data.cell),
+        cell_matrix(data.cell_matrix),
+        cell_vector(data.cell_vector)
+      {}
+    }
+  }
+}
+
+
+namespace MatrixCreator
+{
+
+  template <int dim, typename number, int spacedim>
+  void create_mass_matrix (const Mapping<dim,spacedim>       &mapping,
+                           const DoFHandler<dim,spacedim>    &dof,
+                           const Quadrature<dim>    &q,
+                           SparseMatrix<number>     &matrix,
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
+  {
+    Assert (matrix.m() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+    Assert (matrix.n() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
+
+    hp::FECollection<dim,spacedim>      fe_collection (dof.get_fe());
+    hp::QCollection<dim>                q_collection (q);
+    hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
+    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
+    assembler_data (fe_collection,
+                    update_values | update_JxW_values |
+                    (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)),
+                    coefficient, /*rhs_function=*/0,
+                    q_collection, mapping_collection);
+
+    MatrixCreator::internal::AssemblerData::CopyData<number> copy_data;
+    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
+                                  assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.constraints = &constraints;
+
+    WorkStream::run (dof.begin_active(),
+                     static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
+                     &MatrixCreator::internal::mass_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator,number>,
+                     std_cxx11::bind (&MatrixCreator::internal::
+                                      copy_local_to_global<number,SparseMatrix<number>, Vector<number> >,
+                                      std_cxx11::_1, &matrix, (Vector<number> *)0),
+                     assembler_data,
+                     copy_data);
+  }
+
+
+
+  template <int dim, typename number, int spacedim>
+  void create_mass_matrix (const DoFHandler<dim,spacedim>    &dof,
+                           const Quadrature<dim>    &q,
+                           SparseMatrix<number>     &matrix,
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
+  {
+    create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof,
+                       q, matrix, coefficient, constraints);
+  }
+
+
+
+  template <int dim, typename number, int spacedim>
+  void create_mass_matrix (const Mapping<dim,spacedim>       &mapping,
+                           const DoFHandler<dim,spacedim>    &dof,
+                           const Quadrature<dim>    &q,
+                           SparseMatrix<number>     &matrix,
+                           const Function<spacedim>      &rhs,
+                           Vector<number>           &rhs_vector,
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
+  {
+    Assert (matrix.m() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+    Assert (matrix.n() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
+
+    hp::FECollection<dim,spacedim>      fe_collection (dof.get_fe());
+    hp::QCollection<dim>                q_collection (q);
+    hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
+    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
+    assembler_data (fe_collection,
+                    update_values |
+                    update_JxW_values | update_quadrature_points,
+                    coefficient, &rhs,
+                    q_collection, mapping_collection);
+    MatrixCreator::internal::AssemblerData::CopyData<number> copy_data;
+    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
+                                  assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.constraints = &constraints;
+
+    WorkStream::run (dof.begin_active(),
+                     static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
+                     &MatrixCreator::internal::mass_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator,number>,
+                     std_cxx11::bind(&MatrixCreator::internal::
+                                     copy_local_to_global<number,SparseMatrix<number>, Vector<number> >,
+                                     std_cxx11::_1, &matrix, &rhs_vector),
+                     assembler_data,
+                     copy_data);
+  }
+
+
+
+  template <int dim, typename number, int spacedim>
+  void create_mass_matrix (const DoFHandler<dim,spacedim>    &dof,
+                           const Quadrature<dim>    &q,
+                           SparseMatrix<number>     &matrix,
+                           const Function<spacedim>      &rhs,
+                           Vector<number>           &rhs_vector,
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
+  {
+    create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping,
+                       dof, q, matrix, rhs, rhs_vector, coefficient,
+                       constraints);
+  }
+
+
+
+  template <int dim, typename number, int spacedim>
+  void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
+                           const hp::DoFHandler<dim,spacedim>    &dof,
+                           const hp::QCollection<dim>    &q,
+                           SparseMatrix<number>     &matrix,
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
+  {
+    Assert (matrix.m() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+    Assert (matrix.n() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
+
+    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
+    assembler_data (dof.get_fe(),
+                    update_values | update_JxW_values |
+                    (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)),
+                    coefficient, /*rhs_function=*/0,
+                    q, mapping);
+    MatrixCreator::internal::AssemblerData::CopyData<number> copy_data;
+    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
+                                  assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.constraints = &constraints;
+
+    WorkStream::run (dof.begin_active(),
+                     static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
+                     &MatrixCreator::internal::mass_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,number>,
+                     std_cxx11::bind (&MatrixCreator::internal::
+                                      copy_local_to_global<number,SparseMatrix<number>, Vector<number> >,
+                                      std_cxx11::_1, &matrix, (Vector<number> *)0),
+                     assembler_data,
+                     copy_data);
+  }
+
+
+
+  template <int dim, typename number, int spacedim>
+  void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
+                           const hp::QCollection<dim> &q,
+                           SparseMatrix<number>     &matrix,
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
+  {
+    create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
+                       dof, q, matrix, coefficient, constraints);
+  }
+
+
+
+  template <int dim, typename number, int spacedim>
+  void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
+                           const hp::DoFHandler<dim,spacedim> &dof,
+                           const hp::QCollection<dim>    &q,
+                           SparseMatrix<number>     &matrix,
+                           const Function<spacedim>      &rhs,
+                           Vector<number>           &rhs_vector,
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
+  {
+    Assert (matrix.m() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+    Assert (matrix.n() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
+
+    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
+    assembler_data (dof.get_fe(),
+                    update_values |
+                    update_JxW_values | update_quadrature_points,
+                    coefficient, &rhs,
+                    q, mapping);
+    MatrixCreator::internal::AssemblerData::CopyData<number> copy_data;
+    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
+                                  assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.constraints = &constraints;
+
+    WorkStream::run (dof.begin_active(),
+                     static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
+                     &MatrixCreator::internal::mass_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,number>,
+                     std_cxx11::bind (&MatrixCreator::internal::
+                                      copy_local_to_global<number,SparseMatrix<number>, Vector<number> >,
+                                      std_cxx11::_1, &matrix, &rhs_vector),
+                     assembler_data,
+                     copy_data);
+  }
+
+
+
+  template <int dim, typename number, int spacedim>
+  void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
+                           const hp::QCollection<dim> &q,
+                           SparseMatrix<number>     &matrix,
+                           const Function<spacedim> &rhs,
+                           Vector<number>           &rhs_vector,
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
+  {
+    create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
+                       matrix, rhs, rhs_vector, coefficient, constraints);
+  }
+
+
+
+  namespace internal
+  {
+    template <int dim, int spacedim>
+    void
+    static inline
+    create_boundary_mass_matrix_1 (typename DoFHandler<dim,spacedim>::active_cell_iterator const &cell,
+                                   MatrixCreator::internal::AssemblerBoundary::Scratch const &,
+                                   MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,
+                                   spacedim> > &copy_data,
+                                   Mapping<dim, spacedim> const &mapping,
+                                   FiniteElement<dim,spacedim> const &fe,
+                                   Quadrature<dim-1> const &q,
+                                   typename FunctionMap<spacedim>::type const &boundary_functions,
+                                   Function<spacedim> const *const coefficient,
+                                   std::vector<unsigned int> const &component_mapping)
+
+    {
+      // All assertions for this function
+      // are in the calling function
+      // before creating threads.
+      const unsigned int n_components = fe.n_components();
+      const unsigned int n_function_components = boundary_functions.begin()->second->n_components;
+      const bool         fe_is_system = (n_components != 1);
+      const bool         fe_is_primitive = fe.is_primitive();
+
+      const unsigned int dofs_per_face = fe.dofs_per_face;
+
+      copy_data.cell = cell;
+      copy_data.dofs_per_cell = fe.dofs_per_cell;
+
+      UpdateFlags update_flags = UpdateFlags (update_values     |
+                                              update_JxW_values |
+                                              update_normal_vectors |
+                                              update_quadrature_points);
+      FEFaceValues<dim,spacedim> fe_values (mapping, fe, q, update_flags);
+
+      // two variables for the coefficient,
+      // one for the two cases indicated in
+      // the name
+      std::vector<double>          coefficient_values (fe_values.n_quadrature_points, 1.);
+      std::vector<Vector<double> > coefficient_vector_values (fe_values.n_quadrature_points,
+                                                              Vector<double>(n_components));
+      const bool coefficient_is_vector = (coefficient != 0 && coefficient->n_components != 1)
+                                         ? true : false;
+
+      std::vector<double>          rhs_values_scalar (fe_values.n_quadrature_points);
+      std::vector<Vector<double> > rhs_values_system (fe_values.n_quadrature_points,
+                                                      Vector<double>(n_function_components));
+
+      copy_data.dofs.resize(copy_data.dofs_per_cell);
+      cell->get_dof_indices (copy_data.dofs);
+
+      std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face);
+
+      // Because CopyData objects are reused and that push_back is used,
+      // dof_is_on_face, cell_matrix, and cell_vector must be cleared before
+      // they are reused
+      copy_data.dof_is_on_face.clear();
+      copy_data.cell_matrix.clear();
+      copy_data.cell_vector.clear();
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        // check if this face is on that part of
+        // the boundary we are interested in
+        if (boundary_functions.find(cell->face(face)->boundary_id()) !=
+            boundary_functions.end())
+          {
+            copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell,
+                                                                copy_data.dofs_per_cell));
+            copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell));
+            fe_values.reinit (cell, face);
+
+            if (fe_is_system)
+              // FE has several components
+              {
+                boundary_functions.find(cell->face(face)->boundary_id())
+                ->second->vector_value_list (fe_values.get_quadrature_points(),
+                                             rhs_values_system);
+
+                if (coefficient_is_vector)
+                  // If coefficient is
+                  // vector valued, fill
+                  // all components
+                  coefficient->vector_value_list (fe_values.get_quadrature_points(),
+                                                  coefficient_vector_values);
+                else
+                  {
+                    // If a scalar
+                    // function is
+                    // given, update
+                    // the values, if
+                    // not, use the
+                    // default one set
+                    // in the
+                    // constructor above
+                    if (coefficient != 0)
+                      coefficient->value_list (fe_values.get_quadrature_points(),
+                                               coefficient_values);
+                    // Copy scalar
+                    // values into vector
+                    for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                      coefficient_vector_values[point] = coefficient_values[point];
+                  }
+
+                // Special treatment
+                // for Hdiv and Hcurl
+                // elements, where only
+                // the normal or
+                // tangential component
+                // should be projected.
+                std::vector<std::vector<double> > normal_adjustment(fe_values.n_quadrature_points,
+                                                                    std::vector<double>(n_components, 1.));
+
+                for (unsigned int comp = 0; comp<n_components; ++comp)
+                  {
+                    const FiniteElement<dim,spacedim> &base = fe.base_element(fe.component_to_base_index(comp).first);
+                    const unsigned int bcomp = fe.component_to_base_index(comp).second;
+
+                    if (!base.conforms(FiniteElementData<dim>::H1) &&
+                        base.conforms(FiniteElementData<dim>::Hdiv))
+                      for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                        normal_adjustment[point][comp] = fe_values.normal_vector(point)[bcomp]
+                                                         * fe_values.normal_vector(point)[bcomp];
+                  }
+
+                for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                  {
+                    const double weight = fe_values.JxW(point);
+                    for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                      if (fe_is_primitive)
+                        {
+                          for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                            {
+                              if (fe.system_to_component_index(j).first
+                                  == fe.system_to_component_index(i).first)
+                                {
+                                  copy_data.cell_matrix.back()(i,j)
+                                  += weight
+                                     * fe_values.shape_value(j,point)
+                                     * fe_values.shape_value(i,point)
+                                     * coefficient_vector_values[point](fe.system_to_component_index(i).first);
+                                }
+                            }
+                          copy_data.cell_vector.back()(i) += fe_values.shape_value(i,point)
+                                                             * rhs_values_system[point](component_mapping[fe.system_to_component_index(i).first])
+                                                             * weight;
+                        }
+                      else
+                        {
+                          for (unsigned int comp=0; comp<n_components; ++comp)
+                            {
+                              for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                                copy_data.cell_matrix.back()(i,j)
+                                += fe_values.shape_value_component(j,point,comp)
+                                   * fe_values.shape_value_component(i,point,comp)
+                                   * normal_adjustment[point][comp]
+                                   * weight * coefficient_vector_values[point](comp);
+                              copy_data.cell_vector.back()(i) += fe_values.shape_value_component(i,point,comp) *
+                                                                 rhs_values_system[point](component_mapping[comp])
+                                                                 * normal_adjustment[point][comp]
+                                                                 * weight;
+                            }
+                        }
+                  }
+              }
+            else
+              // FE is a scalar one
+              {
+                boundary_functions.find(cell->face(face)->boundary_id())
+                ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
+
+                if (coefficient != 0)
+                  coefficient->value_list (fe_values.get_quadrature_points(),
+                                           coefficient_values);
+                for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                  {
+                    const double weight = fe_values.JxW(point);
+                    for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                      {
+                        const double v = fe_values.shape_value(i,point);
+                        for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                          {
+                            const double u = fe_values.shape_value(j,point);
+                            copy_data.cell_matrix.back()(i,j) += (u*v*weight*coefficient_values[point]);
+                          }
+                        copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight;
+                      }
+                  }
+              }
+
+
+            cell->face(face)->get_dof_indices (dofs_on_face_vector);
+            // for each dof on the cell, have a
+            // flag whether it is on the face
+            copy_data.dof_is_on_face.push_back(std::vector<bool> (copy_data.dofs_per_cell));
+            // check for each of the dofs on this cell
+            // whether it is on the face
+            for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+              copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(),
+                                                              dofs_on_face_vector.end(),
+                                                              copy_data.dofs[i])
+                                                    !=
+                                                    dofs_on_face_vector.end());
+          }
+    }
+
+    template <int dim,int spacedim>
+    void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,
+                                     spacedim> > const &copy_data,
+                                     typename FunctionMap<spacedim>::type const &boundary_functions,
+                                     std::vector<types::global_dof_index> const &dof_to_boundary_mapping,
+                                     SparseMatrix<double> &matrix,
+                                     Vector<double> &rhs_vector)
+    {
+      // now transfer cell matrix and vector to the whole boundary matrix
+      //
+      // in the following: dof[i] holds the global index of the i-th degree of
+      // freedom on the present cell. If it is also a dof on the boundary, it
+      // must be a nonzero entry in the dof_to_boundary_mapping and then
+      // the boundary index of this dof is dof_to_boundary_mapping[dof[i]].
+      //
+      // if dof[i] is not on the boundary, it should be zero on the boundary
+      // therefore on all quadrature points and finally all of its
+      // entries in the cell matrix and vector should be zero. If not, we
+      // throw an error (note: because of the evaluation of the shape
+      // functions only up to machine precision, the term "must be zero"
+      // really should mean: "should be very small". since this is only an
+      // assertion and not part of the code, we may choose "very small"
+      // quite arbitrarily)
+      //
+      // the main problem here is that the matrix or vector entry should also
+      // be zero if the degree of freedom dof[i] is on the boundary, but not
+      // on the present face, i.e. on another face of the same cell also
+      // on the boundary. We can therefore not rely on the
+      // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to
+      // determine whether dof[i] is a dof on the present face. We do so
+      // by getting the dofs on the face into @p{dofs_on_face_vector},
+      // a vector as always. Usually, searching in a vector is
+      // inefficient, so we copy the dofs into a set, which enables binary
+      // searches.
+      unsigned int pos(0);
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        {
+          // check if this face is on that part of
+          // the boundary we are interested in
+          if (boundary_functions.find(copy_data.cell->face(face)->boundary_id()) !=
+              boundary_functions.end())
+            {
+              for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+                {
+                  if (copy_data.dof_is_on_face[pos][i] &&
+                      dof_to_boundary_mapping[copy_data.dofs[i]] != numbers::invalid_dof_index)
+                    {
+                      for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+                        if (copy_data.dof_is_on_face[pos][j] &&
+                            dof_to_boundary_mapping[copy_data.dofs[j]] != numbers::invalid_dof_index)
+                          {
+                            AssertIsFinite(copy_data.cell_matrix[pos](i,j));
+                            matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]],
+                                       dof_to_boundary_mapping[copy_data.dofs[j]],
+                                       copy_data.cell_matrix[pos](i,j));
+                          }
+                      AssertIsFinite(copy_data.cell_vector[pos](i));
+                      rhs_vector(dof_to_boundary_mapping[copy_data.dofs[i]]) += copy_data.cell_vector[pos](i);
+                    }
+                }
+              ++pos;
+            }
+        }
+    }
+
+
+    template <>
+    void
+    inline
+    create_boundary_mass_matrix_1<1,3> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/,
+                                        MatrixCreator::internal::AssemblerBoundary::Scratch const &,
+                                        MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<1,
+                                        3> > &/*copy_data*/,
+                                        Mapping<1,3> const &,
+                                        FiniteElement<1,3> const &,
+                                        Quadrature<0> const &,
+                                        FunctionMap<3>::type const &/*boundary_functions*/,
+                                        Function<3> const *const /*coefficient*/,
+                                        std::vector<unsigned int> const &/*component_mapping*/)
+    {
+      Assert(false,ExcNotImplemented());
+    }
+  }
+
+
+
+
+  template <int dim, int spacedim>
+  void
+  create_boundary_mass_matrix (const Mapping<dim, spacedim>  &mapping,
+                               const DoFHandler<dim,spacedim> &dof,
+                               const Quadrature<dim-1>  &q,
+                               SparseMatrix<double>  &matrix,
+                               const typename FunctionMap<spacedim>::type  &boundary_functions,
+                               Vector<double>            &rhs_vector,
+                               std::vector<types::global_dof_index> &dof_to_boundary_mapping,
+                               const Function<spacedim> *const coefficient,
+                               std::vector<unsigned int> component_mapping)
+  {
+    // what would that be in 1d? the
+    // identity matrix on the boundary
+    // dofs?
+    if (dim == 1)
+      {
+        Assert (false, ExcNotImplemented());
+        return;
+      }
+
+    const FiniteElement<dim,spacedim> &fe = dof.get_fe();
+    const unsigned int n_components  = fe.n_components();
+
+    Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions),
+            ExcInternalError());
+    Assert (matrix.n() == matrix.m(), ExcInternalError());
+    Assert (matrix.n() == rhs_vector.size(), ExcInternalError());
+    Assert (boundary_functions.size() != 0, ExcInternalError());
+    Assert (dof_to_boundary_mapping.size() == dof.n_dofs(),
+            ExcInternalError());
+    Assert (coefficient ==0 ||
+            coefficient->n_components==1 ||
+            coefficient->n_components==n_components, ExcComponentMismatch());
+
+    if (component_mapping.size() == 0)
+      {
+        AssertDimension (n_components, boundary_functions.begin()->second->n_components);
+        for (unsigned int i=0; i<n_components; ++i)
+          component_mapping.push_back(i);
+      }
+    else
+      AssertDimension (n_components, component_mapping.size());
+
+    MatrixCreator::internal::AssemblerBoundary::Scratch scratch;
+    MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > copy_data;
+
+    WorkStream::run(dof.begin_active(),dof.end(),
+                    static_cast<std_cxx11::function<void (typename DoFHandler<dim,spacedim>::active_cell_iterator
+                                                          const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
+                                                          MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > &)> >
+                    (std_cxx11::bind(&internal::create_boundary_mass_matrix_1<dim,spacedim>,std_cxx11::_1,std_cxx11::_2,
+                                     std_cxx11::_3,
+                                     std_cxx11::cref(mapping),std_cxx11::cref(fe),std_cxx11::cref(q),
+                                     std_cxx11::cref(boundary_functions),coefficient,
+                                     std_cxx11::cref(component_mapping))),
+                    static_cast<std_cxx11::function<void (MatrixCreator::internal::AssemblerBoundary
+                                                          ::CopyData<DoFHandler<dim,spacedim> > const &)> > (std_cxx11::bind(
+                                                                &internal::copy_boundary_mass_matrix_1<dim,spacedim>,
+                                                                std_cxx11::_1,
+                                                                std_cxx11::cref(boundary_functions),
+                                                                std_cxx11::cref(dof_to_boundary_mapping),
+                                                                std_cxx11::ref(matrix),
+                                                                std_cxx11::ref(rhs_vector))),
+                    scratch,
+                    copy_data);
+  }
+
+
+
+  namespace
+  {
+
+    template <int dim, int spacedim>
+    void
+    create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator const
+                                      &cell,
+                                      MatrixCreator::internal::AssemblerBoundary::Scratch const &,
+                                      MatrixCreator::internal::AssemblerBoundary
+                                      ::CopyData<hp::DoFHandler<dim,spacedim> > &copy_data,
+                                      hp::MappingCollection<dim,spacedim> const &mapping,
+                                      hp::FECollection<dim,spacedim> const &fe_collection,
+                                      hp::QCollection<dim-1> const &q,
+                                      const typename FunctionMap<spacedim>::type &boundary_functions,
+                                      Function<spacedim> const *const coefficient,
+                                      std::vector<unsigned int> const &component_mapping)
+    {
+      const unsigned int n_components  = fe_collection.n_components();
+      const unsigned int n_function_components = boundary_functions.begin()->second->n_components;
+      const bool         fe_is_system  = (n_components != 1);
+      const FiniteElement<dim,spacedim> &fe = cell->get_fe();
+      const unsigned int dofs_per_face = fe.dofs_per_face;
+
+      copy_data.cell = cell;
+      copy_data.dofs_per_cell = fe.dofs_per_cell;
+      copy_data.dofs.resize(copy_data.dofs_per_cell);
+      cell->get_dof_indices (copy_data.dofs);
+
+
+      UpdateFlags update_flags = UpdateFlags (update_values     |
+                                              update_JxW_values |
+                                              update_quadrature_points);
+      hp::FEFaceValues<dim,spacedim> x_fe_values (mapping, fe_collection, q, update_flags);
+
+      // two variables for the coefficient,
+      // one for the two cases indicated in
+      // the name
+      std::vector<double>          coefficient_values;
+      std::vector<Vector<double> > coefficient_vector_values;
+
+      std::vector<double>          rhs_values_scalar;
+      std::vector<Vector<double> > rhs_values_system;
+
+      std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face);
+
+      copy_data.dofs.resize(copy_data.dofs_per_cell);
+      cell->get_dof_indices (copy_data.dofs);
+
+      // Because CopyData objects are reused and that push_back is used,
+      // dof_is_on_face, cell_matrix, and cell_vector must be cleared before
+      // they are reused
+      copy_data.dof_is_on_face.clear();
+      copy_data.cell_matrix.clear();
+      copy_data.cell_vector.clear();
+
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        // check if this face is on that part of
+        // the boundary we are interested in
+        if (boundary_functions.find(cell->face(face)->boundary_id()) !=
+            boundary_functions.end())
+          {
+            x_fe_values.reinit (cell, face);
+
+            const FEFaceValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values ();
+
+            copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell,
+                                                                copy_data.dofs_per_cell));
+            copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell));
+
+            if (fe_is_system)
+              // FE has several components
+              {
+                rhs_values_system.resize (fe_values.n_quadrature_points,
+                                          Vector<double>(n_function_components));
+                boundary_functions.find(cell->face(face)->boundary_id())
+                ->second->vector_value_list (fe_values.get_quadrature_points(),
+                                             rhs_values_system);
+
+                if (coefficient != 0)
+                  {
+                    if (coefficient->n_components==1)
+                      {
+                        coefficient_values.resize (fe_values.n_quadrature_points);
+                        coefficient->value_list (fe_values.get_quadrature_points(),
+                                                 coefficient_values);
+                        for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                          {
+                            const double weight = fe_values.JxW(point);
+                            for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                              {
+                                const double v = fe_values.shape_value(i,point);
+                                for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                                  if (fe.system_to_component_index(i).first ==
+                                      fe.system_to_component_index(j).first)
+                                    {
+                                      const double u = fe_values.shape_value(j,point);
+                                      copy_data.cell_matrix.back()(i,j)
+                                      += (u * v * weight * coefficient_values[point]);
+                                    }
+
+                                copy_data.cell_vector.back()(i) += v *
+                                                                   rhs_values_system[point](
+                                                                     component_mapping[fe.system_to_component_index(i).first]) * weight;
+                              }
+                          }
+                      }
+                    else
+                      {
+                        coefficient_vector_values.resize (fe_values.n_quadrature_points,
+                                                          Vector<double>(n_components));
+                        coefficient->vector_value_list (fe_values.get_quadrature_points(),
+                                                        coefficient_vector_values);
+                        for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                          {
+                            const double weight = fe_values.JxW(point);
+                            for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                              {
+                                const double v = fe_values.shape_value(i,point);
+                                const unsigned int component_i=
+                                  fe.system_to_component_index(i).first;
+                                for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                                  if (fe.system_to_component_index(j).first ==
+                                      component_i)
+                                    {
+                                      const double u = fe_values.shape_value(j,point);
+                                      copy_data.cell_matrix.back()(i,j) +=
+                                        (u * v * weight * coefficient_vector_values[point](component_i));
+                                    }
+                                copy_data.cell_vector.back()(i) += v *
+                                                                   rhs_values_system[point](component_mapping[component_i]) * weight;
+                              }
+                          }
+                      }
+                  }
+                else  //      if (coefficient == 0)
+                  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                    {
+                      const double weight = fe_values.JxW(point);
+                      for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                        {
+                          const double v = fe_values.shape_value(i,point);
+                          for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                            if (fe.system_to_component_index(i).first ==
+                                fe.system_to_component_index(j).first)
+                              {
+                                const double u = fe_values.shape_value(j,point);
+                                copy_data.cell_matrix.back()(i,j) += (u * v * weight);
+                              }
+                          copy_data.cell_vector.back()(i) += v *
+                                                             rhs_values_system[point](
+                                                               fe.system_to_component_index(i).first) *
+                                                             weight;
+                        }
+                    }
+              }
+            else
+              // FE is a scalar one
+              {
+                rhs_values_scalar.resize (fe_values.n_quadrature_points);
+                boundary_functions.find(cell->face(face)->boundary_id())
+                ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
+
+                if (coefficient != 0)
+                  {
+                    coefficient_values.resize (fe_values.n_quadrature_points);
+                    coefficient->value_list (fe_values.get_quadrature_points(),
+                                             coefficient_values);
+                    for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                      {
+                        const double weight = fe_values.JxW(point);
+                        for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                          {
+                            const double v = fe_values.shape_value(i,point);
+                            for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                              {
+                                const double u = fe_values.shape_value(j,point);
+                                copy_data.cell_matrix.back()(i,j) += (u * v * weight *
+                                                                      coefficient_values[point]);
+                              }
+                            copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight;
+                          }
+                      }
+                  }
+                else
+                  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+                    {
+                      const double weight = fe_values.JxW(point);
+                      for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
+                        {
+                          const double v = fe_values.shape_value(i,point);
+                          for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
+                            {
+                              const double u = fe_values.shape_value(j,point);
+                              copy_data.cell_matrix.back()(i,j) += (u * v * weight);
+                            }
+                          copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] * weight;
+                        }
+                    }
+              }
+
+            cell->face(face)->get_dof_indices (dofs_on_face_vector,
+                                               cell->active_fe_index());
+            // for each dof on the cell, have a
+            // flag whether it is on the face
+            copy_data.dof_is_on_face.push_back(std::vector<bool> (copy_data.dofs_per_cell));
+            // check for each of the dofs on this cell
+            // whether it is on the face
+            for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+              copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(),
+                                                              dofs_on_face_vector.end(),
+                                                              copy_data.dofs[i])
+                                                    !=
+                                                    dofs_on_face_vector.end());
+          }
+    }
+
+
+
+    template <int dim,int spacedim>
+    void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary
+                                        ::CopyData<hp::DoFHandler<dim,spacedim> > const &copy_data,
+                                        typename FunctionMap<spacedim>::type const &boundary_functions,
+                                        std::vector<types::global_dof_index> const &dof_to_boundary_mapping,
+                                        SparseMatrix<double> &matrix,
+                                        Vector<double> &rhs_vector)
+    {
+      // now transfer cell matrix and vector to the whole boundary matrix
+      //
+      // in the following: dof[i] holds the  global index of the i-th degree of
+      // freedom on the present cell. If it is also a dof on the boundary, it
+      // must be a nonzero entry in the dof_to_boundary_mapping and then
+      // the boundary index of this dof is dof_to_boundary_mapping[dof[i]].
+      //
+      // if dof[i] is not on the boundary, it should be zero on the boundary
+      // therefore on all quadrature points and finally all of its
+      // entries in the cell matrix and vector should be zero. If not, we
+      // throw an error (note: because of the evaluation of the shape
+      // functions only up to machine precision, the term "must be zero"
+      // really should mean: "should be very small". since this is only an
+      // assertion and not part of the code, we may choose "very small"
+      // quite arbitrarily)
+      //
+      // the main problem here is that the matrix or vector entry should also
+      // be zero if the degree of freedom dof[i] is on the boundary, but not
+      // on the present face, i.e. on another face of the same cell also
+      // on the boundary. We can therefore not rely on the
+      // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to
+      // determine whether dof[i] is a dof on the present face. We do so
+      // by getting the dofs on the face into @p{dofs_on_face_vector},
+      // a vector as always. Usually, searching in a vector is
+      // inefficient, so we copy the dofs into a set, which enables binary
+      // searches.
+      unsigned int pos(0);
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        {
+          // check if this face is on that part of
+          // the boundary we are interested in
+          if (boundary_functions.find(copy_data.cell->face(face)->boundary_id()) !=
+              boundary_functions.end())
+            {
+#ifdef DEBUG
+              // in debug mode: compute an element in the matrix which is
+              // guaranteed to belong to a boundary dof. We do this to check that the
+              // entries in the cell matrix are guaranteed to be zero if the
+              // respective dof is not on the boundary. Since because of
+              // round-off, the actual value of the matrix entry may be
+              // only close to zero, we assert that it is small relative to an element
+              // which is guaranteed to be nonzero. (absolute smallness does not
+              // suffice since the size of the domain scales in here)
+              //
+              // for this purpose we seek the diagonal of the matrix, where there
+              // must be an element belonging to the boundary. we take the maximum
+              // diagonal entry.
+              types::global_dof_index max_element = static_cast<types::global_dof_index>(0);
+              for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
+                   i!=dof_to_boundary_mapping.end(); ++i)
+                if ((*i != hp::DoFHandler<dim,spacedim>::invalid_dof_index) &&
+                    (*i > max_element))
+                  max_element = *i;
+              Assert (max_element  == matrix.n()-1, ExcInternalError());
+
+              double max_diag_entry = 0;
+              for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+                if (std::fabs(copy_data.cell_matrix[pos](i,i)) > max_diag_entry)
+                  max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i));
+#endif
+
+              for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+                for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+                  {
+                    if (copy_data.dof_is_on_face[pos][i] && copy_data.dof_is_on_face[pos][j])
+                      matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]],
+                                 dof_to_boundary_mapping[copy_data.dofs[j]],
+                                 copy_data.cell_matrix[pos](i,j));
+                    else
+                      {
+                        // assume that all shape functions that are nonzero on the boundary
+                        // are also listed in the @p{dof_to_boundary} mapping. if that
+                        // is not the case, then the boundary mass matrix does not
+                        // make that much sense anyway, as it only contains entries for
+                        // parts of the functions living on the boundary
+                        //
+                        // these, we may compare here for relative smallness of all
+                        // entries in the local matrix which are not taken over to
+                        // the global one
+                        Assert (std::fabs(copy_data.cell_matrix[pos](i,j)) <= 1e-10 * max_diag_entry,
+                                ExcInternalError ());
+                      }
+                  }
+
+              for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+                if (copy_data.dof_is_on_face[pos][j])
+                  rhs_vector(dof_to_boundary_mapping[copy_data.dofs[j]]) += copy_data.cell_vector[pos](j);
+                else
+                  {
+                    // compare here for relative
+                    // smallness
+                    Assert (std::fabs(copy_data.cell_vector[pos](j)) <= 1e-10 * max_diag_entry,
+                            ExcInternalError());
+                  }
+              ++pos;
+            }
+        }
+    }
+  }
+
+
+
+  template <int dim, int spacedim>
+  void create_boundary_mass_matrix (const DoFHandler<dim,spacedim>     &dof,
+                                    const Quadrature<dim-1>   &q,
+                                    SparseMatrix<double>      &matrix,
+                                    const typename FunctionMap<spacedim>::type &rhs,
+                                    Vector<double>            &rhs_vector,
+                                    std::vector<types::global_dof_index> &dof_to_boundary_mapping,
+                                    const Function<spacedim> *const a,
+                                    std::vector<unsigned int> component_mapping)
+  {
+    create_boundary_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q,
+                                matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void
+  create_boundary_mass_matrix (const hp::MappingCollection<dim,spacedim>        &mapping,
+                               const hp::DoFHandler<dim,spacedim>     &dof,
+                               const hp::QCollection<dim-1>   &q,
+                               SparseMatrix<double>      &matrix,
+                               const typename FunctionMap<spacedim>::type         &boundary_functions,
+                               Vector<double>            &rhs_vector,
+                               std::vector<types::global_dof_index> &dof_to_boundary_mapping,
+                               const Function<spacedim> *const coefficient,
+                               std::vector<unsigned int> component_mapping)
+  {
+    // what would that be in 1d? the
+    // identity matrix on the boundary
+    // dofs?
+    if (dim == 1)
+      {
+        Assert (false, ExcNotImplemented());
+        return;
+      }
+
+    const hp::FECollection<dim,spacedim> &fe_collection = dof.get_fe();
+    const unsigned int n_components  = fe_collection.n_components();
+
+    Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions),
+            ExcInternalError());
+    Assert (matrix.n() == matrix.m(), ExcInternalError());
+    Assert (matrix.n() == rhs_vector.size(), ExcInternalError());
+    Assert (boundary_functions.size() != 0, ExcInternalError());
+    Assert (dof_to_boundary_mapping.size() == dof.n_dofs(),
+            ExcInternalError());
+    Assert (coefficient ==0 ||
+            coefficient->n_components==1 ||
+            coefficient->n_components==n_components, ExcComponentMismatch());
+
+    if (component_mapping.size() == 0)
+      {
+        AssertDimension (n_components, boundary_functions.begin()->second->n_components);
+        for (unsigned int i=0; i<n_components; ++i)
+          component_mapping.push_back(i);
+      }
+    else
+      AssertDimension (n_components, component_mapping.size());
+
+    MatrixCreator::internal::AssemblerBoundary::Scratch scratch;
+    MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > copy_data;
+
+    WorkStream::run(dof.begin_active(),dof.end(),
+                    static_cast<std_cxx11::function<void (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
+                                                          const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
+                                                          MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > &)> >
+                    (std_cxx11::bind( &create_hp_boundary_mass_matrix_1<dim,spacedim>,std_cxx11::_1,std_cxx11::_2,
+                                      std_cxx11::_3,
+                                      std_cxx11::cref(mapping),std_cxx11::cref(fe_collection),std_cxx11::cref(q),
+                                      std_cxx11::cref(boundary_functions),coefficient,
+                                      std_cxx11::cref(component_mapping))),
+                    static_cast<std_cxx11::function<void (MatrixCreator::internal::AssemblerBoundary
+                                                          ::CopyData<hp::DoFHandler<dim,spacedim> > const &)> > (
+                      std_cxx11::bind( &copy_hp_boundary_mass_matrix_1<dim,spacedim>,
+                                       std_cxx11::_1,
+                                       std_cxx11::cref(boundary_functions),
+                                       std_cxx11::cref(dof_to_boundary_mapping),
+                                       std_cxx11::ref(matrix),
+                                       std_cxx11::ref(rhs_vector))),
+                    scratch,
+                    copy_data);
+  }
+
+
+
+
+  template <int dim, int spacedim>
+  void create_boundary_mass_matrix (const hp::DoFHandler<dim,spacedim>     &dof,
+                                    const hp::QCollection<dim-1>   &q,
+                                    SparseMatrix<double>      &matrix,
+                                    const typename FunctionMap<spacedim>::type &rhs,
+                                    Vector<double>            &rhs_vector,
+                                    std::vector<types::global_dof_index> &dof_to_boundary_mapping,
+                                    const Function<spacedim> *const a,
+                                    std::vector<unsigned int> component_mapping)
+  {
+    create_boundary_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
+                                matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void create_laplace_matrix (const Mapping<dim, spacedim>       &mapping,
+                              const DoFHandler<dim,spacedim>    &dof,
+                              const Quadrature<dim>    &q,
+                              SparseMatrix<double>     &matrix,
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
+  {
+    Assert (matrix.m() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+    Assert (matrix.n() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
+
+    hp::FECollection<dim,spacedim>      fe_collection (dof.get_fe());
+    hp::QCollection<dim>                q_collection (q);
+    hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
+    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
+    assembler_data (fe_collection,
+                    update_gradients  | update_JxW_values |
+                    (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)),
+                    coefficient, /*rhs_function=*/0,
+                    q_collection, mapping_collection);
+    MatrixCreator::internal::AssemblerData::CopyData<double> copy_data;
+    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
+                                  assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.constraints = &constraints;
+
+    WorkStream::run (dof.begin_active(),
+                     static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
+                     &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator>,
+                     std_cxx11::bind (&MatrixCreator::internal::
+                                      copy_local_to_global<double,SparseMatrix<double>, Vector<double> >,
+                                      std_cxx11::_1,
+                                      &matrix,
+                                      (Vector<double> *)NULL),
+                     assembler_data,
+                     copy_data);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void create_laplace_matrix (const DoFHandler<dim,spacedim>    &dof,
+                              const Quadrature<dim>    &q,
+                              SparseMatrix<double>     &matrix,
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix &constraints)
+  {
+    create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping,
+                          dof, q, matrix, coefficient, constraints);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void create_laplace_matrix (const Mapping<dim, spacedim>       &mapping,
+                              const DoFHandler<dim,spacedim>    &dof,
+                              const Quadrature<dim>    &q,
+                              SparseMatrix<double>     &matrix,
+                              const Function<spacedim>      &rhs,
+                              Vector<double>           &rhs_vector,
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
+  {
+    Assert (matrix.m() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+    Assert (matrix.n() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
+
+    hp::FECollection<dim,spacedim>      fe_collection (dof.get_fe());
+    hp::QCollection<dim>                q_collection (q);
+    hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
+    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
+    assembler_data (fe_collection,
+                    update_gradients  | update_values |
+                    update_JxW_values | update_quadrature_points,
+                    coefficient, &rhs,
+                    q_collection, mapping_collection);
+    MatrixCreator::internal::AssemblerData::CopyData<double> copy_data;
+    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
+                                  assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.constraints = &constraints;
+
+    WorkStream::run (dof.begin_active(),
+                     static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
+                     &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator>,
+                     std_cxx11::bind (&MatrixCreator::internal::
+                                      copy_local_to_global<double,SparseMatrix<double>, Vector<double> >,
+                                      std_cxx11::_1,
+                                      &matrix,
+                                      &rhs_vector),
+                     assembler_data,
+                     copy_data);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void create_laplace_matrix (const DoFHandler<dim,spacedim>    &dof,
+                              const Quadrature<dim>    &q,
+                              SparseMatrix<double>     &matrix,
+                              const Function<spacedim>      &rhs,
+                              Vector<double>           &rhs_vector,
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix &constraints)
+  {
+    create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q,
+                          matrix, rhs, rhs_vector, coefficient, constraints);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void create_laplace_matrix (const hp::MappingCollection<dim,spacedim>       &mapping,
+                              const hp::DoFHandler<dim,spacedim>    &dof,
+                              const hp::QCollection<dim>    &q,
+                              SparseMatrix<double>     &matrix,
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
+  {
+    Assert (matrix.m() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+    Assert (matrix.n() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
+
+    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
+    assembler_data (dof.get_fe(),
+                    update_gradients  | update_JxW_values |
+                    (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)),
+                    coefficient, /*rhs_function=*/0,
+                    q, mapping);
+    MatrixCreator::internal::AssemblerData::CopyData<double> copy_data;
+    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
+                                  assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.constraints = &constraints;
+
+    WorkStream::run (dof.begin_active(),
+                     static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
+                     &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>,
+                     std_cxx11::bind (&MatrixCreator::internal::
+                                      copy_local_to_global<double,SparseMatrix<double>, Vector<double> >,
+                                      std_cxx11::_1,
+                                      &matrix,
+                                      (Vector<double> *)0),
+                     assembler_data,
+                     copy_data);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void create_laplace_matrix (const hp::DoFHandler<dim,spacedim>    &dof,
+                              const hp::QCollection<dim>    &q,
+                              SparseMatrix<double>     &matrix,
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix &constraints)
+  {
+    create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
+                          matrix, coefficient, constraints);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void create_laplace_matrix (const hp::MappingCollection<dim,spacedim>       &mapping,
+                              const hp::DoFHandler<dim,spacedim>    &dof,
+                              const hp::QCollection<dim>    &q,
+                              SparseMatrix<double>     &matrix,
+                              const Function<spacedim>      &rhs,
+                              Vector<double>           &rhs_vector,
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
+  {
+    Assert (matrix.m() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
+    Assert (matrix.n() == dof.n_dofs(),
+            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
+
+    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
+    assembler_data (dof.get_fe(),
+                    update_gradients  | update_values |
+                    update_JxW_values | update_quadrature_points,
+                    coefficient, &rhs,
+                    q, mapping);
+    MatrixCreator::internal::AssemblerData::CopyData<double> copy_data;
+    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
+                                  assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
+    copy_data.constraints = &constraints;
+
+    WorkStream::run (dof.begin_active(),
+                     static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
+                     &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>,
+                     std_cxx11::bind (&MatrixCreator::internal::
+                                      copy_local_to_global<double,SparseMatrix<double>, Vector<double> >,
+                                      std_cxx11::_1,
+                                      &matrix,
+                                      &rhs_vector),
+                     assembler_data,
+                     copy_data);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void create_laplace_matrix (const hp::DoFHandler<dim,spacedim>    &dof,
+                              const hp::QCollection<dim>    &q,
+                              SparseMatrix<double>     &matrix,
+                              const Function<spacedim>      &rhs,
+                              Vector<double>           &rhs_vector,
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix &constraints)
+  {
+    create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
+                          matrix, rhs, rhs_vector, coefficient, constraints);
+  }
+
+}  // namespace MatrixCreator
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 086368040c8ab7cb027056b467d8ba6c73a5d3f1..c8adfa339bfbe8687a5e975b3fb4b83bb716e773 100644 (file)
@@ -29,6 +29,10 @@ SET(_src
   error_estimator_inst2.cc
   fe_field_function.cc
   histogram.cc
+  matrix_creator.cc
+  matrix_creator_inst2.cc
+  matrix_creator_inst3.cc
+  matrix_tools_once.cc
   matrix_tools.cc
   point_value_history.cc
   solution_transfer.cc
@@ -59,6 +63,7 @@ SET(_inst
   error_estimator_1d.inst.in
   error_estimator.inst.in
   fe_field_function.inst.in
+  matrix_creator.inst.in
   matrix_tools.inst.in
   point_value_history.inst.in
   solution_transfer.inst.in
diff --git a/source/numerics/matrix_creator.cc b/source/numerics/matrix_creator.cc
new file mode 100644 (file)
index 0000000..1bb6205
--- /dev/null
@@ -0,0 +1,28 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/numerics/matrix_creator.templates.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+// explicit instantiations
+#define SPLIT_INSTANTIATIONS_COUNT 3
+#ifndef SPLIT_INSTANTIATIONS_INDEX
+#define SPLIT_INSTANTIATIONS_INDEX 0
+#endif
+#include "matrix_creator.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/numerics/matrix_creator.inst.in b/source/numerics/matrix_creator.inst.in
new file mode 100644 (file)
index 0000000..a84cafc
--- /dev/null
@@ -0,0 +1,315 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2010 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+  {
+#if deal_II_dimension <= deal_II_space_dimension
+
+// non-hp version of create_mass_matrix
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
+      (const Mapping<deal_II_dimension,deal_II_space_dimension>       &mapping,
+       const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
+      (const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
+      (const Mapping<deal_II_dimension,deal_II_space_dimension>       &mapping,
+       const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_space_dimension>      &rhs,
+       Vector<double>           &rhs_vector,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
+      (const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_space_dimension>      &rhs,
+       Vector<double>           &rhs_vector,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+
+    template
+      void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension>
+      (const DoFHandler<deal_II_dimension,deal_II_space_dimension>     &dof,
+       const Quadrature<deal_II_dimension-1>   &q,
+       SparseMatrix<double>      &matrix,
+       const FunctionMap<deal_II_space_dimension>::type &rhs,
+       Vector<double>            &rhs_vector,
+       std::vector<types::global_dof_index> &dof_to_boundary_mapping,
+       const Function<deal_II_space_dimension> * const a,
+       std::vector<unsigned int>);
+
+    template
+      void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension>
+      (const Mapping<deal_II_dimension,deal_II_space_dimension> &,
+       const DoFHandler<deal_II_dimension,deal_II_space_dimension>     &dof,
+       const Quadrature<deal_II_dimension-1>   &q,
+       SparseMatrix<double>      &matrix,
+       const FunctionMap<deal_II_space_dimension>::type &rhs,
+       Vector<double>            &rhs_vector,
+       std::vector<types::global_dof_index> &dof_to_boundary_mapping,
+       const Function<deal_II_space_dimension> * const a,
+       std::vector<unsigned int>);
+
+    template
+      void
+      MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension>
+      (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>&,
+       const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
+       const hp::QCollection<deal_II_dimension-1>&,
+       SparseMatrix<double>&,
+       const FunctionMap<deal_II_space_dimension>::type&,
+       Vector<double>&,
+       std::vector<types::global_dof_index>&,
+       const Function<deal_II_space_dimension> * const,
+       std::vector<unsigned int>);
+
+    template
+      void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension>
+      (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
+       const hp::QCollection<deal_II_dimension-1>&,
+       SparseMatrix<double>&,
+       const FunctionMap<deal_II_space_dimension>::type&,
+       Vector<double>&,
+       std::vector<types::global_dof_index>&,
+       const Function<deal_II_space_dimension> * const,
+       std::vector<unsigned int>);
+
+// same for float
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
+      (const Mapping<deal_II_dimension,deal_II_space_dimension>       &mapping,
+       const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<float>     &matrix,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
+      (const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<float>     &matrix,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
+      (const Mapping<deal_II_dimension,deal_II_space_dimension>       &mapping,
+       const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<float>     &matrix,
+       const Function<deal_II_space_dimension>      &rhs,
+       Vector<float>           &rhs_vector,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
+      (const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<float>     &matrix,
+       const Function<deal_II_space_dimension>      &rhs,
+       Vector<float>           &rhs_vector,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+#endif
+  }
+
+
+//TODO[SP]: replace <deal_II_dimension> by <deal_II_dimension, deal_II_space_dimension>
+// where applicable and move to codimension cases above also when applicable
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+  {
+// hp versions of functions
+#if deal_II_dimension <= deal_II_space_dimension
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
+      (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
+       const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
+      (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
+       const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_space_dimension>      &rhs,
+       Vector<double>           &rhs_vector,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+
+#endif
+
+#if deal_II_dimension == deal_II_space_dimension
+
+       
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension>
+      (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension>
+      (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_space_dimension>      &rhs,
+       Vector<double>           &rhs_vector,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+// non-hp versions of create_laplace_matrix
+    template
+      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
+      (const DoFHandler<deal_II_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
+      (const Mapping<deal_II_dimension>       &mapping,
+       const DoFHandler<deal_II_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
+      (const Mapping<deal_II_dimension>       &mapping,
+       const DoFHandler<deal_II_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_dimension>      &rhs,
+       Vector<double>           &rhs_vector,
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
+      (const DoFHandler<deal_II_dimension>    &dof,
+       const Quadrature<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_dimension>      &rhs,
+       Vector<double>           &rhs_vector,
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+// hp versions of create_laplace_matrix
+    template
+      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
+      (const hp::DoFHandler<deal_II_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
+      (const hp::MappingCollection<deal_II_dimension>       &mapping,
+       const hp::DoFHandler<deal_II_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
+      (const hp::MappingCollection<deal_II_dimension>       &mapping,
+       const hp::DoFHandler<deal_II_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_dimension>      &rhs,
+       Vector<double>           &rhs_vector,
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+    template
+      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
+      (const hp::DoFHandler<deal_II_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<double>     &matrix,
+       const Function<deal_II_dimension>      &rhs,
+       Vector<double>           &rhs_vector,
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+#endif
+
+// same for float
+#if deal_II_dimension <= deal_II_space_dimension
+
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
+      (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
+       const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<float>     &matrix,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
+      (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
+       const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<float>     &matrix,
+       const Function<deal_II_space_dimension>      &rhs,
+       Vector<float>           &rhs_vector,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+
+#endif
+
+#if deal_II_dimension == deal_II_space_dimension
+
+       
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension>
+      (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<float>     &matrix,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+    template
+      void MatrixCreator::create_mass_matrix<deal_II_dimension>
+      (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
+       const hp::QCollection<deal_II_dimension>    &q,
+       SparseMatrix<float>     &matrix,
+       const Function<deal_II_space_dimension>      &rhs,
+       Vector<float>           &rhs_vector,
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
+
+#endif
diff --git a/source/numerics/matrix_creator_inst2.cc b/source/numerics/matrix_creator_inst2.cc
new file mode 100644 (file)
index 0000000..73c34c5
--- /dev/null
@@ -0,0 +1,20 @@
+// ---------------------------------------------------------------------
+//
+// 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 file compiles a part of the instantiations from matrix_creator.cc
+// to reduce the compilation unit (and memory consumption)
+
+#define SPLIT_INSTANTIATIONS_INDEX 1
+#include "matrix_creator.cc"
diff --git a/source/numerics/matrix_creator_inst3.cc b/source/numerics/matrix_creator_inst3.cc
new file mode 100644 (file)
index 0000000..2cc0df2
--- /dev/null
@@ -0,0 +1,20 @@
+// ---------------------------------------------------------------------
+//
+// 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 file compiles a part of the instantiations from matrix_creator.cc
+// to reduce the compilation unit (and memory consumption)
+
+#define SPLIT_INSTANTIATIONS_INDEX 2
+#include "matrix_creator.cc"
index ddd53584222366612c4e9751779270d681456d97..28b46756c126a6937ad9645648d320f454582c64 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
-
-namespace MatrixCreator
-{
-  namespace internal
-  {
-    /**
-     * Convenience abbreviation for
-     * pairs of DoF handler cell
-     * iterators. This type works
-     * just like a
-     * <tt>std::pair<iterator,iterator></tt>
-     * but is templatized on the
-     * dof handler that should be used.
-     */
-    template <typename DoFHandlerType>
-    struct IteratorRange
-    {
-      /**
-       * Typedef for the iterator type.
-       */
-      typedef typename DoFHandlerType::active_cell_iterator active_cell_iterator;
-
-      /**
-       * Abbreviation for a pair of
-       * iterators.
-       */
-      typedef std::pair<active_cell_iterator,active_cell_iterator> iterator_pair;
-
-      /**
-       * Constructor. Initialize
-       * the two values by the
-       * given values.
-       */
-      IteratorRange (const active_cell_iterator &first,
-                     const active_cell_iterator &second);
-
-      /**
-       * Constructor taking a pair
-       * of values for
-       * initialization.
-       */
-      IteratorRange (const iterator_pair &ip);
-
-      /**
-       * Pair of iterators denoting
-       * a half-open range.
-       */
-      active_cell_iterator first, second;
-    };
-
-
-
-
-    template <typename DoFHandlerType>
-    inline
-    IteratorRange<DoFHandlerType>::
-    IteratorRange (const active_cell_iterator &first,
-                   const active_cell_iterator &second)
-      :
-      first (first),
-      second (second)
-    {}
-
-
-
-    template <typename DoFHandlerType>
-    inline
-    IteratorRange<DoFHandlerType>::IteratorRange (const iterator_pair &ip)
-      :
-      first (ip.first),
-      second (ip.second)
-    {}
-
-
-
-    namespace AssemblerData
-    {
-      template <int dim,
-                int spacedim>
-      struct Scratch
-      {
-        Scratch (const ::dealii::hp::FECollection<dim,spacedim> &fe,
-                 const UpdateFlags         update_flags,
-                 const Function<spacedim> *coefficient,
-                 const Function<spacedim> *rhs_function,
-                 const ::dealii::hp::QCollection<dim> &quadrature,
-                 const ::dealii::hp::MappingCollection<dim,spacedim> &mapping)
-          :
-          fe_collection (fe),
-          quadrature_collection (quadrature),
-          mapping_collection (mapping),
-          x_fe_values (mapping_collection,
-                       fe_collection,
-                       quadrature_collection,
-                       update_flags),
-          coefficient_values(quadrature_collection.max_n_quadrature_points()),
-          coefficient_vector_values (quadrature_collection.max_n_quadrature_points(),
-                                     dealii::Vector<double> (fe_collection.n_components())),
-          rhs_values(quadrature_collection.max_n_quadrature_points()),
-          rhs_vector_values(quadrature_collection.max_n_quadrature_points(),
-                            dealii::Vector<double> (fe_collection.n_components())),
-          coefficient (coefficient),
-          rhs_function (rhs_function),
-          update_flags (update_flags)
-        {}
-
-        Scratch (const Scratch &data)
-          :
-          fe_collection (data.fe_collection),
-          quadrature_collection (data.quadrature_collection),
-          mapping_collection (data.mapping_collection),
-          x_fe_values (mapping_collection,
-                       fe_collection,
-                       quadrature_collection,
-                       data.update_flags),
-          coefficient_values (data.coefficient_values),
-          coefficient_vector_values (data.coefficient_vector_values),
-          rhs_values (data.rhs_values),
-          rhs_vector_values (data.rhs_vector_values),
-          coefficient (data.coefficient),
-          rhs_function (data.rhs_function),
-          update_flags (data.update_flags)
-        {}
-
-        const ::dealii::hp::FECollection<dim,spacedim>      &fe_collection;
-        const ::dealii::hp::QCollection<dim>                &quadrature_collection;
-        const ::dealii::hp::MappingCollection<dim,spacedim> &mapping_collection;
-
-        ::dealii::hp::FEValues<dim,spacedim> x_fe_values;
-
-        std::vector<double>                  coefficient_values;
-        std::vector<dealii::Vector<double> > coefficient_vector_values;
-        std::vector<double>                  rhs_values;
-        std::vector<dealii::Vector<double> > rhs_vector_values;
-
-        const Function<spacedim>   *coefficient;
-        const Function<spacedim>   *rhs_function;
-
-        const UpdateFlags update_flags;
-      };
-
-
-      template <typename number>
-      struct CopyData
-      {
-        std::vector<types::global_dof_index> dof_indices;
-        FullMatrix<number>        cell_matrix;
-        dealii::Vector<number>    cell_rhs;
-        const ConstraintMatrix   *constraints;
-      };
-    }
-
-
-    template <int dim,
-              int spacedim,
-              typename CellIterator,
-              typename number>
-    void mass_assembler (const CellIterator &cell,
-                         MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim> &data,
-                         MatrixCreator::internal::AssemblerData::CopyData<number>      &copy_data)
-    {
-      data.x_fe_values.reinit (cell);
-      const FEValues<dim,spacedim> &fe_values = data.x_fe_values.get_present_fe_values ();
-
-      const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                         n_q_points    = fe_values.n_quadrature_points;
-      const FiniteElement<dim,spacedim> &fe  = fe_values.get_fe();
-      const unsigned int n_components  = fe.n_components();
-
-      Assert(data.rhs_function == 0 ||
-             data.rhs_function->n_components==1 ||
-             data.rhs_function->n_components==n_components,
-             ::dealii::MatrixCreator::ExcComponentMismatch());
-      Assert(data.coefficient == 0 ||
-             data.coefficient->n_components==1 ||
-             data.coefficient->n_components==n_components,
-             ::dealii::MatrixCreator::ExcComponentMismatch());
-
-      copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
-      copy_data.cell_rhs.reinit (dofs_per_cell);
-
-      copy_data.dof_indices.resize (dofs_per_cell);
-      cell->get_dof_indices (copy_data.dof_indices);
-
-      const bool use_rhs_function = data.rhs_function != 0;
-      if (use_rhs_function)
-        {
-          if (data.rhs_function->n_components==1)
-            {
-              data.rhs_values.resize (n_q_points);
-              data.rhs_function->value_list (fe_values.get_quadrature_points(),
-                                             data.rhs_values);
-            }
-          else
-            {
-              data.rhs_vector_values.resize (n_q_points,
-                                             dealii::Vector<double>(n_components));
-              data.rhs_function->vector_value_list (fe_values.get_quadrature_points(),
-                                                    data.rhs_vector_values);
-            }
-        }
-
-      const bool use_coefficient = data.coefficient != 0;
-      if (use_coefficient)
-        {
-          if (data.coefficient->n_components==1)
-            {
-              data.coefficient_values.resize (n_q_points);
-              data.coefficient->value_list (fe_values.get_quadrature_points(),
-                                            data.coefficient_values);
-            }
-          else
-            {
-              data.coefficient_vector_values.resize (n_q_points,
-                                                     dealii::Vector<double>(n_components));
-              data.coefficient->vector_value_list (fe_values.get_quadrature_points(),
-                                                   data.coefficient_vector_values);
-            }
-        }
-
-
-      double add_data;
-      const std::vector<double> &JxW = fe_values.get_JxW_values();
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        if (fe.is_primitive ())
-          {
-            const unsigned int component_i =
-              fe.system_to_component_index(i).first;
-            const double *phi_i = &fe_values.shape_value(i,0);
-            add_data = 0;
-
-            // use symmetry in the mass matrix here:
-            // just need to calculate the diagonal
-            // and half of the elements above the
-            // diagonal
-            for (unsigned int j=i; j<dofs_per_cell; ++j)
-              if ((n_components==1) ||
-                  (fe.system_to_component_index(j).first ==
-                   component_i))
-                {
-                  const double *phi_j = &fe_values.shape_value(j,0);
-                  add_data = 0;
-                  if (use_coefficient)
-                    {
-                      if (data.coefficient->n_components==1)
-                        for (unsigned int point=0; point<n_q_points; ++point)
-                          add_data += (phi_i[point] * phi_j[point] * JxW[point] *
-                                       data.coefficient_values[point]);
-                      else
-                        for (unsigned int point=0; point<n_q_points; ++point)
-                          add_data += (phi_i[point] * phi_j[point] * JxW[point] *
-                                       data.coefficient_vector_values[point](component_i));
-                    }
-                  else
-                    for (unsigned int point=0; point<n_q_points; ++point)
-                      add_data += phi_i[point] * phi_j[point] * JxW[point];
-
-                  // this is even ok for i==j, since then
-                  // we just write the same value twice.
-                  copy_data.cell_matrix(i,j) = add_data;
-                  copy_data.cell_matrix(j,i) = add_data;
-                }
-
-            if (use_rhs_function)
-              {
-                add_data = 0;
-                if (data.rhs_function->n_components==1)
-                  for (unsigned int point=0; point<n_q_points; ++point)
-                    add_data += phi_i[point] * JxW[point] *
-                                data.rhs_values[point];
-                else
-                  for (unsigned int point=0; point<n_q_points; ++point)
-                    add_data += phi_i[point] * JxW[point] *
-                                data.rhs_vector_values[point](component_i);
-                copy_data.cell_rhs(i) = add_data;
-              }
-          }
-        else
-          {
-            // non-primitive vector-valued FE, using
-            // symmetry again
-            for (unsigned int j=i; j<dofs_per_cell; ++j)
-              {
-                add_data = 0;
-                for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-                  if (fe.get_nonzero_components(i)[comp_i] &&
-                      fe.get_nonzero_components(j)[comp_i])
-                    {
-                      if (use_coefficient)
-                        {
-                          if (data.coefficient->n_components==1)
-                            for (unsigned int point=0; point<n_q_points; ++point)
-                              add_data += (fe_values.shape_value_component(i,point,comp_i) *
-                                           fe_values.shape_value_component(j,point,comp_i) *
-                                           JxW[point] *
-                                           data.coefficient_values[point]);
-                          else
-                            for (unsigned int point=0; point<n_q_points; ++point)
-                              add_data += (fe_values.shape_value_component(i,point,comp_i) *
-                                           fe_values.shape_value_component(j,point,comp_i) *
-                                           JxW[point] *
-                                           data.coefficient_vector_values[point](comp_i));
-                        }
-                      else
-                        for (unsigned int point=0; point<n_q_points; ++point)
-                          add_data += fe_values.shape_value_component(i,point,comp_i) *
-                                      fe_values.shape_value_component(j,point,comp_i) * JxW[point];
-                    }
-
-                copy_data.cell_matrix(i,j) = add_data;
-                copy_data.cell_matrix(j,i) = add_data;
-              }
-
-            if (use_rhs_function)
-              {
-                add_data = 0;
-                for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-                  if (fe.get_nonzero_components(i)[comp_i])
-                    {
-                      if (data.rhs_function->n_components==1)
-                        for (unsigned int point=0; point<n_q_points; ++point)
-                          add_data += fe_values.shape_value_component(i,point,comp_i) *
-                                      JxW[point] * data.rhs_values[point];
-                      else
-                        for (unsigned int point=0; point<n_q_points; ++point)
-                          add_data += fe_values.shape_value_component(i,point,comp_i) *
-                                      JxW[point] * data.rhs_vector_values[point](comp_i);
-                    }
-                copy_data.cell_rhs(i) = add_data;
-              }
-          }
-    }
-
-
-
-    template <int dim,
-              int spacedim,
-              typename CellIterator>
-    void laplace_assembler (const CellIterator &cell,
-                            MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim> &data,
-                            MatrixCreator::internal::AssemblerData::CopyData<double>      &copy_data)
-    {
-      data.x_fe_values.reinit (cell);
-      const FEValues<dim,spacedim> &fe_values = data.x_fe_values.get_present_fe_values ();
-
-      const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                         n_q_points    = fe_values.n_quadrature_points;
-      const FiniteElement<dim,spacedim>    &fe  = fe_values.get_fe();
-      const unsigned int n_components  = fe.n_components();
-
-      Assert(data.rhs_function == 0 ||
-             data.rhs_function->n_components==1 ||
-             data.rhs_function->n_components==n_components,
-             ::dealii::MatrixCreator::ExcComponentMismatch());
-      Assert(data.coefficient == 0 ||
-             data.coefficient->n_components==1 ||
-             data.coefficient->n_components==n_components,
-             ::dealii::MatrixCreator::ExcComponentMismatch());
-
-      copy_data.cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
-      copy_data.cell_rhs.reinit (dofs_per_cell);
-      copy_data.dof_indices.resize (dofs_per_cell);
-      cell->get_dof_indices (copy_data.dof_indices);
-
-
-      const bool use_rhs_function = data.rhs_function != 0;
-      if (use_rhs_function)
-        {
-          if (data.rhs_function->n_components==1)
-            {
-              data.rhs_values.resize (n_q_points);
-              data.rhs_function->value_list (fe_values.get_quadrature_points(),
-                                             data.rhs_values);
-            }
-          else
-            {
-              data.rhs_vector_values.resize (n_q_points,
-                                             dealii::Vector<double>(n_components));
-              data.rhs_function->vector_value_list (fe_values.get_quadrature_points(),
-                                                    data.rhs_vector_values);
-            }
-        }
-
-      const bool use_coefficient = data.coefficient != 0;
-      if (use_coefficient)
-        {
-          if (data.coefficient->n_components==1)
-            {
-              data.coefficient_values.resize (n_q_points);
-              data.coefficient->value_list (fe_values.get_quadrature_points(),
-                                            data.coefficient_values);
-            }
-          else
-            {
-              data.coefficient_vector_values.resize (n_q_points,
-                                                     dealii::Vector<double>(n_components));
-              data.coefficient->vector_value_list (fe_values.get_quadrature_points(),
-                                                   data.coefficient_vector_values);
-            }
-        }
-
-
-      const std::vector<double> &JxW = fe_values.get_JxW_values();
-      double add_data;
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        if (fe.is_primitive ())
-          {
-            const unsigned int component_i =
-              fe.system_to_component_index(i).first;
-            const Tensor<1,spacedim> *grad_phi_i =
-              &fe_values.shape_grad(i,0);
-
-            // can use symmetry
-            for (unsigned int j=i; j<dofs_per_cell; ++j)
-              if ((n_components==1) ||
-                  (fe.system_to_component_index(j).first ==
-                   component_i))
-                {
-                  const Tensor<1,spacedim> *grad_phi_j =
-                    & fe_values.shape_grad(j,0);
-                  add_data = 0;
-                  if (use_coefficient)
-                    {
-                      if (data.coefficient->n_components==1)
-                        for (unsigned int point=0; point<n_q_points; ++point)
-                          add_data += ((grad_phi_i[point]*grad_phi_j[point]) *
-                                       JxW[point] *
-                                       data.coefficient_values[point]);
-                      else
-                        for (unsigned int point=0; point<n_q_points; ++point)
-                          add_data += ((grad_phi_i[point]*grad_phi_j[point]) *
-                                       JxW[point] *
-                                       data.coefficient_vector_values[point](component_i));
-                    }
-                  else
-                    for (unsigned int point=0; point<n_q_points; ++point)
-                      add_data += (grad_phi_i[point]*grad_phi_j[point]) *
-                                  JxW[point];
-
-                  copy_data.cell_matrix(i,j) = add_data;
-                  copy_data.cell_matrix(j,i) = add_data;
-                }
-
-            if (use_rhs_function)
-              {
-                const double *phi_i = &fe_values.shape_value(i,0);
-                add_data = 0;
-                if (data.rhs_function->n_components==1)
-                  for (unsigned int point=0; point<n_q_points; ++point)
-                    add_data += phi_i[point] * JxW[point] *
-                                data.rhs_values[point];
-                else
-                  for (unsigned int point=0; point<n_q_points; ++point)
-                    add_data += phi_i[point] * JxW[point] *
-                                data.rhs_vector_values[point](component_i);
-                copy_data.cell_rhs(i) = add_data;
-              }
-          }
-        else
-          {
-            // non-primitive vector-valued FE
-            for (unsigned int j=i; j<dofs_per_cell; ++j)
-              {
-                add_data = 0;
-                for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-                  if (fe.get_nonzero_components(i)[comp_i] &&
-                      fe.get_nonzero_components(j)[comp_i])
-                    {
-                      if (use_coefficient)
-                        {
-                          if (data.coefficient->n_components==1)
-                            for (unsigned int point=0; point<n_q_points; ++point)
-                              add_data += ((fe_values.shape_grad_component(i,point,comp_i) *
-                                            fe_values.shape_grad_component(j,point,comp_i)) *
-                                           JxW[point] *
-                                           data.coefficient_values[point]);
-                          else
-                            for (unsigned int point=0; point<n_q_points; ++point)
-                              add_data += ((fe_values.shape_grad_component(i,point,comp_i) *
-                                            fe_values.shape_grad_component(j,point,comp_i)) *
-                                           JxW[point] *
-                                           data.coefficient_vector_values[point](comp_i));
-                        }
-                      else
-                        for (unsigned int point=0; point<n_q_points; ++point)
-                          add_data += (fe_values.shape_grad_component(i,point,comp_i) *
-                                       fe_values.shape_grad_component(j,point,comp_i)) *
-                                      JxW[point];
-                    }
-
-                copy_data.cell_matrix(i,j) = add_data;
-                copy_data.cell_matrix(j,i) = add_data;
-              }
-
-            if (use_rhs_function)
-              {
-                add_data = 0;
-                for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-                  if (fe.get_nonzero_components(i)[comp_i])
-                    {
-                      if (data.rhs_function->n_components==1)
-                        for (unsigned int point=0; point<n_q_points; ++point)
-                          add_data += fe_values.shape_value_component(i,point,comp_i) *
-                                      JxW[point] * data.rhs_values[point];
-                      else
-                        for (unsigned int point=0; point<n_q_points; ++point)
-                          add_data += fe_values.shape_value_component(i,point,comp_i) *
-                                      JxW[point] * data.rhs_vector_values[point](comp_i);
-                    }
-                copy_data.cell_rhs(i) = add_data;
-              }
-          }
-    }
-
-
-
-    template <typename number,
-              typename MatrixType,
-              typename VectorType>
-    void copy_local_to_global (const AssemblerData::CopyData<number> &data,
-                               MatrixType *matrix,
-                               VectorType *right_hand_side)
-    {
-      const unsigned int dofs_per_cell = data.dof_indices.size();
-      (void)dofs_per_cell;
-
-      Assert (data.cell_matrix.m() == dofs_per_cell,
-              ExcInternalError());
-      Assert (data.cell_matrix.n() == dofs_per_cell,
-              ExcInternalError());
-      Assert ((right_hand_side == 0)
-              ||
-              (data.cell_rhs.size() == dofs_per_cell),
-              ExcInternalError());
-
-      if (right_hand_side != 0)
-        data.constraints->distribute_local_to_global(data.cell_matrix,
-                                                     data.cell_rhs,
-                                                     data.dof_indices,
-                                                     *matrix, *right_hand_side);
-      else
-        data.constraints->distribute_local_to_global(data.cell_matrix,
-                                                     data.dof_indices,
-                                                     *matrix);
-    }
-
-
-
-    namespace AssemblerBoundary
-    {
-      struct Scratch
-      {
-        Scratch() {}
-      };
-
-      template <typename DoFHandlerType>
-      struct CopyData
-      {
-        CopyData() {};
-
-        CopyData(CopyData const &data);
-
-        unsigned int dofs_per_cell;
-        std::vector<types::global_dof_index> dofs;
-        std::vector<std::vector<bool> > dof_is_on_face;
-        typename DoFHandlerType::active_cell_iterator cell;
-        std::vector<FullMatrix<double> > cell_matrix;
-        std::vector<Vector<double> > cell_vector;
-      };
-
-      template <typename DoFHandlerType>
-      CopyData<DoFHandlerType>::CopyData(CopyData const &data) :
-        dofs_per_cell(data.dofs_per_cell),
-        dofs(data.dofs),
-        dof_is_on_face(data.dof_is_on_face),
-        cell(data.cell),
-        cell_matrix(data.cell_matrix),
-        cell_vector(data.cell_vector)
-      {}
-    }
-  }
-}
-
-
-namespace MatrixCreator
-{
-
-  template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const Mapping<dim,spacedim>       &mapping,
-                           const DoFHandler<dim,spacedim>    &dof,
-                           const Quadrature<dim>    &q,
-                           SparseMatrix<number>     &matrix,
-                           const Function<spacedim> *const coefficient,
-                           const ConstraintMatrix   &constraints)
-  {
-    Assert (matrix.m() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
-    Assert (matrix.n() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
-
-    hp::FECollection<dim,spacedim>      fe_collection (dof.get_fe());
-    hp::QCollection<dim>                q_collection (q);
-    hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
-    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
-    assembler_data (fe_collection,
-                    update_values | update_JxW_values |
-                    (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)),
-                    coefficient, /*rhs_function=*/0,
-                    q_collection, mapping_collection);
-
-    MatrixCreator::internal::AssemblerData::CopyData<number> copy_data;
-    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
-                                  assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.constraints = &constraints;
-
-    WorkStream::run (dof.begin_active(),
-                     static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
-                     &MatrixCreator::internal::mass_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator,number>,
-                     std_cxx11::bind (&MatrixCreator::internal::
-                                      copy_local_to_global<number,SparseMatrix<number>, Vector<number> >,
-                                      std_cxx11::_1, &matrix, (Vector<number> *)0),
-                     assembler_data,
-                     copy_data);
-  }
-
-
-
-  template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const DoFHandler<dim,spacedim>    &dof,
-                           const Quadrature<dim>    &q,
-                           SparseMatrix<number>     &matrix,
-                           const Function<spacedim> *const coefficient,
-                           const ConstraintMatrix   &constraints)
-  {
-    create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof,
-                       q, matrix, coefficient, constraints);
-  }
-
-
-
-  template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const Mapping<dim,spacedim>       &mapping,
-                           const DoFHandler<dim,spacedim>    &dof,
-                           const Quadrature<dim>    &q,
-                           SparseMatrix<number>     &matrix,
-                           const Function<spacedim>      &rhs,
-                           Vector<number>           &rhs_vector,
-                           const Function<spacedim> *const coefficient,
-                           const ConstraintMatrix   &constraints)
-  {
-    Assert (matrix.m() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
-    Assert (matrix.n() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
-
-    hp::FECollection<dim,spacedim>      fe_collection (dof.get_fe());
-    hp::QCollection<dim>                q_collection (q);
-    hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
-    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
-    assembler_data (fe_collection,
-                    update_values |
-                    update_JxW_values | update_quadrature_points,
-                    coefficient, &rhs,
-                    q_collection, mapping_collection);
-    MatrixCreator::internal::AssemblerData::CopyData<number> copy_data;
-    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
-                                  assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.constraints = &constraints;
-
-    WorkStream::run (dof.begin_active(),
-                     static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
-                     &MatrixCreator::internal::mass_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator,number>,
-                     std_cxx11::bind(&MatrixCreator::internal::
-                                     copy_local_to_global<number,SparseMatrix<number>, Vector<number> >,
-                                     std_cxx11::_1, &matrix, &rhs_vector),
-                     assembler_data,
-                     copy_data);
-  }
-
-
-
-  template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const DoFHandler<dim,spacedim>    &dof,
-                           const Quadrature<dim>    &q,
-                           SparseMatrix<number>     &matrix,
-                           const Function<spacedim>      &rhs,
-                           Vector<number>           &rhs_vector,
-                           const Function<spacedim> *const coefficient,
-                           const ConstraintMatrix   &constraints)
-  {
-    create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping,
-                       dof, q, matrix, rhs, rhs_vector, coefficient,
-                       constraints);
-  }
-
-
-
-  template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
-                           const hp::DoFHandler<dim,spacedim>    &dof,
-                           const hp::QCollection<dim>    &q,
-                           SparseMatrix<number>     &matrix,
-                           const Function<spacedim> *const coefficient,
-                           const ConstraintMatrix   &constraints)
-  {
-    Assert (matrix.m() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
-    Assert (matrix.n() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
-
-    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
-    assembler_data (dof.get_fe(),
-                    update_values | update_JxW_values |
-                    (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)),
-                    coefficient, /*rhs_function=*/0,
-                    q, mapping);
-    MatrixCreator::internal::AssemblerData::CopyData<number> copy_data;
-    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
-                                  assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.constraints = &constraints;
-
-    WorkStream::run (dof.begin_active(),
-                     static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
-                     &MatrixCreator::internal::mass_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,number>,
-                     std_cxx11::bind (&MatrixCreator::internal::
-                                      copy_local_to_global<number,SparseMatrix<number>, Vector<number> >,
-                                      std_cxx11::_1, &matrix, (Vector<number> *)0),
-                     assembler_data,
-                     copy_data);
-  }
-
-
-
-  template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
-                           const hp::QCollection<dim> &q,
-                           SparseMatrix<number>     &matrix,
-                           const Function<spacedim> *const coefficient,
-                           const ConstraintMatrix   &constraints)
-  {
-    create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
-                       dof, q, matrix, coefficient, constraints);
-  }
-
-
-
-  template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
-                           const hp::DoFHandler<dim,spacedim> &dof,
-                           const hp::QCollection<dim>    &q,
-                           SparseMatrix<number>     &matrix,
-                           const Function<spacedim>      &rhs,
-                           Vector<number>           &rhs_vector,
-                           const Function<spacedim> *const coefficient,
-                           const ConstraintMatrix   &constraints)
-  {
-    Assert (matrix.m() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
-    Assert (matrix.n() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
-
-    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
-    assembler_data (dof.get_fe(),
-                    update_values |
-                    update_JxW_values | update_quadrature_points,
-                    coefficient, &rhs,
-                    q, mapping);
-    MatrixCreator::internal::AssemblerData::CopyData<number> copy_data;
-    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
-                                  assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.constraints = &constraints;
-
-    WorkStream::run (dof.begin_active(),
-                     static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
-                     &MatrixCreator::internal::mass_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,number>,
-                     std_cxx11::bind (&MatrixCreator::internal::
-                                      copy_local_to_global<number,SparseMatrix<number>, Vector<number> >,
-                                      std_cxx11::_1, &matrix, &rhs_vector),
-                     assembler_data,
-                     copy_data);
-  }
-
-
-
-  template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
-                           const hp::QCollection<dim> &q,
-                           SparseMatrix<number>     &matrix,
-                           const Function<spacedim> &rhs,
-                           Vector<number>           &rhs_vector,
-                           const Function<spacedim> *const coefficient,
-                           const ConstraintMatrix   &constraints)
-  {
-    create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
-                       matrix, rhs, rhs_vector, coefficient, constraints);
-  }
-
-
-
-  namespace
-  {
-    template <int dim, int spacedim>
-    void
-    create_boundary_mass_matrix_1 (typename DoFHandler<dim,spacedim>::active_cell_iterator const &cell,
-                                   MatrixCreator::internal::AssemblerBoundary::Scratch const &,
-                                   MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,
-                                   spacedim> > &copy_data,
-                                   Mapping<dim, spacedim> const &mapping,
-                                   FiniteElement<dim,spacedim> const &fe,
-                                   Quadrature<dim-1> const &q,
-                                   typename FunctionMap<spacedim>::type const &boundary_functions,
-                                   Function<spacedim> const *const coefficient,
-                                   std::vector<unsigned int> const &component_mapping)
-
-    {
-      // All assertions for this function
-      // are in the calling function
-      // before creating threads.
-      const unsigned int n_components = fe.n_components();
-      const unsigned int n_function_components = boundary_functions.begin()->second->n_components;
-      const bool         fe_is_system = (n_components != 1);
-      const bool         fe_is_primitive = fe.is_primitive();
-
-      const unsigned int dofs_per_face = fe.dofs_per_face;
-
-      copy_data.cell = cell;
-      copy_data.dofs_per_cell = fe.dofs_per_cell;
-
-      UpdateFlags update_flags = UpdateFlags (update_values     |
-                                              update_JxW_values |
-                                              update_normal_vectors |
-                                              update_quadrature_points);
-      FEFaceValues<dim,spacedim> fe_values (mapping, fe, q, update_flags);
-
-      // two variables for the coefficient,
-      // one for the two cases indicated in
-      // the name
-      std::vector<double>          coefficient_values (fe_values.n_quadrature_points, 1.);
-      std::vector<Vector<double> > coefficient_vector_values (fe_values.n_quadrature_points,
-                                                              Vector<double>(n_components));
-      const bool coefficient_is_vector = (coefficient != 0 && coefficient->n_components != 1)
-                                         ? true : false;
-
-      std::vector<double>          rhs_values_scalar (fe_values.n_quadrature_points);
-      std::vector<Vector<double> > rhs_values_system (fe_values.n_quadrature_points,
-                                                      Vector<double>(n_function_components));
-
-      copy_data.dofs.resize(copy_data.dofs_per_cell);
-      cell->get_dof_indices (copy_data.dofs);
-
-      std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face);
-
-      // Because CopyData objects are reused and that push_back is used,
-      // dof_is_on_face, cell_matrix, and cell_vector must be cleared before
-      // they are reused
-      copy_data.dof_is_on_face.clear();
-      copy_data.cell_matrix.clear();
-      copy_data.cell_vector.clear();
-
-      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-        // check if this face is on that part of
-        // the boundary we are interested in
-        if (boundary_functions.find(cell->face(face)->boundary_id()) !=
-            boundary_functions.end())
-          {
-            copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell,
-                                                                copy_data.dofs_per_cell));
-            copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell));
-            fe_values.reinit (cell, face);
-
-            if (fe_is_system)
-              // FE has several components
-              {
-                boundary_functions.find(cell->face(face)->boundary_id())
-                ->second->vector_value_list (fe_values.get_quadrature_points(),
-                                             rhs_values_system);
-
-                if (coefficient_is_vector)
-                  // If coefficient is
-                  // vector valued, fill
-                  // all components
-                  coefficient->vector_value_list (fe_values.get_quadrature_points(),
-                                                  coefficient_vector_values);
-                else
-                  {
-                    // If a scalar
-                    // function is
-                    // given, update
-                    // the values, if
-                    // not, use the
-                    // default one set
-                    // in the
-                    // constructor above
-                    if (coefficient != 0)
-                      coefficient->value_list (fe_values.get_quadrature_points(),
-                                               coefficient_values);
-                    // Copy scalar
-                    // values into vector
-                    for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                      coefficient_vector_values[point] = coefficient_values[point];
-                  }
-
-                // Special treatment
-                // for Hdiv and Hcurl
-                // elements, where only
-                // the normal or
-                // tangential component
-                // should be projected.
-                std::vector<std::vector<double> > normal_adjustment(fe_values.n_quadrature_points,
-                                                                    std::vector<double>(n_components, 1.));
-
-                for (unsigned int comp = 0; comp<n_components; ++comp)
-                  {
-                    const FiniteElement<dim,spacedim> &base = fe.base_element(fe.component_to_base_index(comp).first);
-                    const unsigned int bcomp = fe.component_to_base_index(comp).second;
-
-                    if (!base.conforms(FiniteElementData<dim>::H1) &&
-                        base.conforms(FiniteElementData<dim>::Hdiv))
-                      for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                        normal_adjustment[point][comp] = fe_values.normal_vector(point)[bcomp]
-                                                         * fe_values.normal_vector(point)[bcomp];
-                  }
-
-                for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                  {
-                    const double weight = fe_values.JxW(point);
-                    for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
-                      if (fe_is_primitive)
-                        {
-                          for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                            {
-                              if (fe.system_to_component_index(j).first
-                                  == fe.system_to_component_index(i).first)
-                                {
-                                  copy_data.cell_matrix.back()(i,j)
-                                  += weight
-                                     * fe_values.shape_value(j,point)
-                                     * fe_values.shape_value(i,point)
-                                     * coefficient_vector_values[point](fe.system_to_component_index(i).first);
-                                }
-                            }
-                          copy_data.cell_vector.back()(i) += fe_values.shape_value(i,point)
-                                                             * rhs_values_system[point](component_mapping[fe.system_to_component_index(i).first])
-                                                             * weight;
-                        }
-                      else
-                        {
-                          for (unsigned int comp=0; comp<n_components; ++comp)
-                            {
-                              for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                                copy_data.cell_matrix.back()(i,j)
-                                += fe_values.shape_value_component(j,point,comp)
-                                   * fe_values.shape_value_component(i,point,comp)
-                                   * normal_adjustment[point][comp]
-                                   * weight * coefficient_vector_values[point](comp);
-                              copy_data.cell_vector.back()(i) += fe_values.shape_value_component(i,point,comp) *
-                                                                 rhs_values_system[point](component_mapping[comp])
-                                                                 * normal_adjustment[point][comp]
-                                                                 * weight;
-                            }
-                        }
-                  }
-              }
-            else
-              // FE is a scalar one
-              {
-                boundary_functions.find(cell->face(face)->boundary_id())
-                ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
-
-                if (coefficient != 0)
-                  coefficient->value_list (fe_values.get_quadrature_points(),
-                                           coefficient_values);
-                for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                  {
-                    const double weight = fe_values.JxW(point);
-                    for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
-                      {
-                        const double v = fe_values.shape_value(i,point);
-                        for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                          {
-                            const double u = fe_values.shape_value(j,point);
-                            copy_data.cell_matrix.back()(i,j) += (u*v*weight*coefficient_values[point]);
-                          }
-                        copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight;
-                      }
-                  }
-              }
-
-
-            cell->face(face)->get_dof_indices (dofs_on_face_vector);
-            // for each dof on the cell, have a
-            // flag whether it is on the face
-            copy_data.dof_is_on_face.push_back(std::vector<bool> (copy_data.dofs_per_cell));
-            // check for each of the dofs on this cell
-            // whether it is on the face
-            for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
-              copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(),
-                                                              dofs_on_face_vector.end(),
-                                                              copy_data.dofs[i])
-                                                    !=
-                                                    dofs_on_face_vector.end());
-          }
-    }
-
-    template <int dim,int spacedim>
-    void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,
-                                     spacedim> > const &copy_data,
-                                     typename FunctionMap<spacedim>::type const &boundary_functions,
-                                     std::vector<types::global_dof_index> const &dof_to_boundary_mapping,
-                                     SparseMatrix<double> &matrix,
-                                     Vector<double> &rhs_vector)
-    {
-      // now transfer cell matrix and vector to the whole boundary matrix
-      //
-      // in the following: dof[i] holds the global index of the i-th degree of
-      // freedom on the present cell. If it is also a dof on the boundary, it
-      // must be a nonzero entry in the dof_to_boundary_mapping and then
-      // the boundary index of this dof is dof_to_boundary_mapping[dof[i]].
-      //
-      // if dof[i] is not on the boundary, it should be zero on the boundary
-      // therefore on all quadrature points and finally all of its
-      // entries in the cell matrix and vector should be zero. If not, we
-      // throw an error (note: because of the evaluation of the shape
-      // functions only up to machine precision, the term "must be zero"
-      // really should mean: "should be very small". since this is only an
-      // assertion and not part of the code, we may choose "very small"
-      // quite arbitrarily)
-      //
-      // the main problem here is that the matrix or vector entry should also
-      // be zero if the degree of freedom dof[i] is on the boundary, but not
-      // on the present face, i.e. on another face of the same cell also
-      // on the boundary. We can therefore not rely on the
-      // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to
-      // determine whether dof[i] is a dof on the present face. We do so
-      // by getting the dofs on the face into @p{dofs_on_face_vector},
-      // a vector as always. Usually, searching in a vector is
-      // inefficient, so we copy the dofs into a set, which enables binary
-      // searches.
-      unsigned int pos(0);
-      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-        {
-          // check if this face is on that part of
-          // the boundary we are interested in
-          if (boundary_functions.find(copy_data.cell->face(face)->boundary_id()) !=
-              boundary_functions.end())
-            {
-              for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
-                {
-                  if (copy_data.dof_is_on_face[pos][i] &&
-                      dof_to_boundary_mapping[copy_data.dofs[i]] != numbers::invalid_dof_index)
-                    {
-                      for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
-                        if (copy_data.dof_is_on_face[pos][j] &&
-                            dof_to_boundary_mapping[copy_data.dofs[j]] != numbers::invalid_dof_index)
-                          {
-                            AssertIsFinite(copy_data.cell_matrix[pos](i,j));
-                            matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]],
-                                       dof_to_boundary_mapping[copy_data.dofs[j]],
-                                       copy_data.cell_matrix[pos](i,j));
-                          }
-                      AssertIsFinite(copy_data.cell_vector[pos](i));
-                      rhs_vector(dof_to_boundary_mapping[copy_data.dofs[i]]) += copy_data.cell_vector[pos](i);
-                    }
-                }
-              ++pos;
-            }
-        }
-    }
-
-
-    template <>
-    void
-    create_boundary_mass_matrix_1<1,3> (DoFHandler<1,3>::active_cell_iterator const &/*cell*/,
-                                        MatrixCreator::internal::AssemblerBoundary::Scratch const &,
-                                        MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<1,
-                                        3> > &/*copy_data*/,
-                                        Mapping<1,3> const &,
-                                        FiniteElement<1,3> const &,
-                                        Quadrature<0> const &,
-                                        FunctionMap<3>::type const &/*boundary_functions*/,
-                                        Function<3> const *const /*coefficient*/,
-                                        std::vector<unsigned int> const &/*component_mapping*/)
-    {
-      Assert(false,ExcNotImplemented());
-    }
-  }
-
-
-
-
-  template <int dim, int spacedim>
-  void
-  create_boundary_mass_matrix (const Mapping<dim, spacedim>  &mapping,
-                               const DoFHandler<dim,spacedim> &dof,
-                               const Quadrature<dim-1>  &q,
-                               SparseMatrix<double>  &matrix,
-                               const typename FunctionMap<spacedim>::type  &boundary_functions,
-                               Vector<double>            &rhs_vector,
-                               std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                               const Function<spacedim> *const coefficient,
-                               std::vector<unsigned int> component_mapping)
-  {
-    // what would that be in 1d? the
-    // identity matrix on the boundary
-    // dofs?
-    if (dim == 1)
-      {
-        Assert (false, ExcNotImplemented());
-        return;
-      }
-
-    const FiniteElement<dim,spacedim> &fe = dof.get_fe();
-    const unsigned int n_components  = fe.n_components();
-
-    Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions),
-            ExcInternalError());
-    Assert (matrix.n() == matrix.m(), ExcInternalError());
-    Assert (matrix.n() == rhs_vector.size(), ExcInternalError());
-    Assert (boundary_functions.size() != 0, ExcInternalError());
-    Assert (dof_to_boundary_mapping.size() == dof.n_dofs(),
-            ExcInternalError());
-    Assert (coefficient ==0 ||
-            coefficient->n_components==1 ||
-            coefficient->n_components==n_components, ExcComponentMismatch());
-
-    if (component_mapping.size() == 0)
-      {
-        AssertDimension (n_components, boundary_functions.begin()->second->n_components);
-        for (unsigned int i=0; i<n_components; ++i)
-          component_mapping.push_back(i);
-      }
-    else
-      AssertDimension (n_components, component_mapping.size());
-
-    MatrixCreator::internal::AssemblerBoundary::Scratch scratch;
-    MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > copy_data;
-
-    WorkStream::run(dof.begin_active(),dof.end(),
-                    static_cast<std_cxx11::function<void (typename DoFHandler<dim,spacedim>::active_cell_iterator
-                                                          const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
-                                                          MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > &)> >
-                    (std_cxx11::bind(&create_boundary_mass_matrix_1<dim,spacedim>,std_cxx11::_1,std_cxx11::_2,
-                                     std_cxx11::_3,
-                                     std_cxx11::cref(mapping),std_cxx11::cref(fe),std_cxx11::cref(q),
-                                     std_cxx11::cref(boundary_functions),coefficient,
-                                     std_cxx11::cref(component_mapping))),
-                    static_cast<std_cxx11::function<void (MatrixCreator::internal::AssemblerBoundary
-                                                          ::CopyData<DoFHandler<dim,spacedim> > const &)> > (std_cxx11::bind(
-                                                                &copy_boundary_mass_matrix_1<dim,spacedim>,
-                                                                std_cxx11::_1,
-                                                                std_cxx11::cref(boundary_functions),
-                                                                std_cxx11::cref(dof_to_boundary_mapping),
-                                                                std_cxx11::ref(matrix),
-                                                                std_cxx11::ref(rhs_vector))),
-                    scratch,
-                    copy_data);
-  }
-
-
-
-  namespace
-  {
-
-    template <int dim, int spacedim>
-    void
-    create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator const
-                                      &cell,
-                                      MatrixCreator::internal::AssemblerBoundary::Scratch const &,
-                                      MatrixCreator::internal::AssemblerBoundary
-                                      ::CopyData<hp::DoFHandler<dim,spacedim> > &copy_data,
-                                      hp::MappingCollection<dim,spacedim> const &mapping,
-                                      hp::FECollection<dim,spacedim> const &fe_collection,
-                                      hp::QCollection<dim-1> const &q,
-                                      const typename FunctionMap<spacedim>::type &boundary_functions,
-                                      Function<spacedim> const *const coefficient,
-                                      std::vector<unsigned int> const &component_mapping)
-    {
-      const unsigned int n_components  = fe_collection.n_components();
-      const unsigned int n_function_components = boundary_functions.begin()->second->n_components;
-      const bool         fe_is_system  = (n_components != 1);
-      const FiniteElement<dim,spacedim> &fe = cell->get_fe();
-      const unsigned int dofs_per_face = fe.dofs_per_face;
-
-      copy_data.cell = cell;
-      copy_data.dofs_per_cell = fe.dofs_per_cell;
-      copy_data.dofs.resize(copy_data.dofs_per_cell);
-      cell->get_dof_indices (copy_data.dofs);
-
-
-      UpdateFlags update_flags = UpdateFlags (update_values     |
-                                              update_JxW_values |
-                                              update_quadrature_points);
-      hp::FEFaceValues<dim,spacedim> x_fe_values (mapping, fe_collection, q, update_flags);
-
-      // two variables for the coefficient,
-      // one for the two cases indicated in
-      // the name
-      std::vector<double>          coefficient_values;
-      std::vector<Vector<double> > coefficient_vector_values;
-
-      std::vector<double>          rhs_values_scalar;
-      std::vector<Vector<double> > rhs_values_system;
-
-      std::vector<types::global_dof_index> dofs_on_face_vector (dofs_per_face);
-
-      copy_data.dofs.resize(copy_data.dofs_per_cell);
-      cell->get_dof_indices (copy_data.dofs);
-
-      // Because CopyData objects are reused and that push_back is used,
-      // dof_is_on_face, cell_matrix, and cell_vector must be cleared before
-      // they are reused
-      copy_data.dof_is_on_face.clear();
-      copy_data.cell_matrix.clear();
-      copy_data.cell_vector.clear();
-
-
-      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-        // check if this face is on that part of
-        // the boundary we are interested in
-        if (boundary_functions.find(cell->face(face)->boundary_id()) !=
-            boundary_functions.end())
-          {
-            x_fe_values.reinit (cell, face);
-
-            const FEFaceValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values ();
-
-            copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell,
-                                                                copy_data.dofs_per_cell));
-            copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell));
-
-            if (fe_is_system)
-              // FE has several components
-              {
-                rhs_values_system.resize (fe_values.n_quadrature_points,
-                                          Vector<double>(n_function_components));
-                boundary_functions.find(cell->face(face)->boundary_id())
-                ->second->vector_value_list (fe_values.get_quadrature_points(),
-                                             rhs_values_system);
-
-                if (coefficient != 0)
-                  {
-                    if (coefficient->n_components==1)
-                      {
-                        coefficient_values.resize (fe_values.n_quadrature_points);
-                        coefficient->value_list (fe_values.get_quadrature_points(),
-                                                 coefficient_values);
-                        for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                          {
-                            const double weight = fe_values.JxW(point);
-                            for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
-                              {
-                                const double v = fe_values.shape_value(i,point);
-                                for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                                  if (fe.system_to_component_index(i).first ==
-                                      fe.system_to_component_index(j).first)
-                                    {
-                                      const double u = fe_values.shape_value(j,point);
-                                      copy_data.cell_matrix.back()(i,j)
-                                      += (u * v * weight * coefficient_values[point]);
-                                    }
-
-                                copy_data.cell_vector.back()(i) += v *
-                                                                   rhs_values_system[point](
-                                                                     component_mapping[fe.system_to_component_index(i).first]) * weight;
-                              }
-                          }
-                      }
-                    else
-                      {
-                        coefficient_vector_values.resize (fe_values.n_quadrature_points,
-                                                          Vector<double>(n_components));
-                        coefficient->vector_value_list (fe_values.get_quadrature_points(),
-                                                        coefficient_vector_values);
-                        for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                          {
-                            const double weight = fe_values.JxW(point);
-                            for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
-                              {
-                                const double v = fe_values.shape_value(i,point);
-                                const unsigned int component_i=
-                                  fe.system_to_component_index(i).first;
-                                for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                                  if (fe.system_to_component_index(j).first ==
-                                      component_i)
-                                    {
-                                      const double u = fe_values.shape_value(j,point);
-                                      copy_data.cell_matrix.back()(i,j) +=
-                                        (u * v * weight * coefficient_vector_values[point](component_i));
-                                    }
-                                copy_data.cell_vector.back()(i) += v *
-                                                                   rhs_values_system[point](component_mapping[component_i]) * weight;
-                              }
-                          }
-                      }
-                  }
-                else  //      if (coefficient == 0)
-                  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                    {
-                      const double weight = fe_values.JxW(point);
-                      for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
-                        {
-                          const double v = fe_values.shape_value(i,point);
-                          for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                            if (fe.system_to_component_index(i).first ==
-                                fe.system_to_component_index(j).first)
-                              {
-                                const double u = fe_values.shape_value(j,point);
-                                copy_data.cell_matrix.back()(i,j) += (u * v * weight);
-                              }
-                          copy_data.cell_vector.back()(i) += v *
-                                                             rhs_values_system[point](
-                                                               fe.system_to_component_index(i).first) *
-                                                             weight;
-                        }
-                    }
-              }
-            else
-              // FE is a scalar one
-              {
-                rhs_values_scalar.resize (fe_values.n_quadrature_points);
-                boundary_functions.find(cell->face(face)->boundary_id())
-                ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar);
-
-                if (coefficient != 0)
-                  {
-                    coefficient_values.resize (fe_values.n_quadrature_points);
-                    coefficient->value_list (fe_values.get_quadrature_points(),
-                                             coefficient_values);
-                    for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                      {
-                        const double weight = fe_values.JxW(point);
-                        for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
-                          {
-                            const double v = fe_values.shape_value(i,point);
-                            for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                              {
-                                const double u = fe_values.shape_value(j,point);
-                                copy_data.cell_matrix.back()(i,j) += (u * v * weight *
-                                                                      coefficient_values[point]);
-                              }
-                            copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight;
-                          }
-                      }
-                  }
-                else
-                  for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
-                    {
-                      const double weight = fe_values.JxW(point);
-                      for (unsigned int i=0; i<fe_values.dofs_per_cell; ++i)
-                        {
-                          const double v = fe_values.shape_value(i,point);
-                          for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
-                            {
-                              const double u = fe_values.shape_value(j,point);
-                              copy_data.cell_matrix.back()(i,j) += (u * v * weight);
-                            }
-                          copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] * weight;
-                        }
-                    }
-              }
-
-            cell->face(face)->get_dof_indices (dofs_on_face_vector,
-                                               cell->active_fe_index());
-            // for each dof on the cell, have a
-            // flag whether it is on the face
-            copy_data.dof_is_on_face.push_back(std::vector<bool> (copy_data.dofs_per_cell));
-            // check for each of the dofs on this cell
-            // whether it is on the face
-            for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
-              copy_data.dof_is_on_face.back()[i] = (std::find(dofs_on_face_vector.begin(),
-                                                              dofs_on_face_vector.end(),
-                                                              copy_data.dofs[i])
-                                                    !=
-                                                    dofs_on_face_vector.end());
-          }
-    }
-
-
-
-    template <int dim,int spacedim>
-    void copy_hp_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary
-                                        ::CopyData<hp::DoFHandler<dim,spacedim> > const &copy_data,
-                                        typename FunctionMap<spacedim>::type const &boundary_functions,
-                                        std::vector<types::global_dof_index> const &dof_to_boundary_mapping,
-                                        SparseMatrix<double> &matrix,
-                                        Vector<double> &rhs_vector)
-    {
-      // now transfer cell matrix and vector to the whole boundary matrix
-      //
-      // in the following: dof[i] holds the  global index of the i-th degree of
-      // freedom on the present cell. If it is also a dof on the boundary, it
-      // must be a nonzero entry in the dof_to_boundary_mapping and then
-      // the boundary index of this dof is dof_to_boundary_mapping[dof[i]].
-      //
-      // if dof[i] is not on the boundary, it should be zero on the boundary
-      // therefore on all quadrature points and finally all of its
-      // entries in the cell matrix and vector should be zero. If not, we
-      // throw an error (note: because of the evaluation of the shape
-      // functions only up to machine precision, the term "must be zero"
-      // really should mean: "should be very small". since this is only an
-      // assertion and not part of the code, we may choose "very small"
-      // quite arbitrarily)
-      //
-      // the main problem here is that the matrix or vector entry should also
-      // be zero if the degree of freedom dof[i] is on the boundary, but not
-      // on the present face, i.e. on another face of the same cell also
-      // on the boundary. We can therefore not rely on the
-      // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to
-      // determine whether dof[i] is a dof on the present face. We do so
-      // by getting the dofs on the face into @p{dofs_on_face_vector},
-      // a vector as always. Usually, searching in a vector is
-      // inefficient, so we copy the dofs into a set, which enables binary
-      // searches.
-      unsigned int pos(0);
-      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-        {
-          // check if this face is on that part of
-          // the boundary we are interested in
-          if (boundary_functions.find(copy_data.cell->face(face)->boundary_id()) !=
-              boundary_functions.end())
-            {
-#ifdef DEBUG
-              // in debug mode: compute an element in the matrix which is
-              // guaranteed to belong to a boundary dof. We do this to check that the
-              // entries in the cell matrix are guaranteed to be zero if the
-              // respective dof is not on the boundary. Since because of
-              // round-off, the actual value of the matrix entry may be
-              // only close to zero, we assert that it is small relative to an element
-              // which is guaranteed to be nonzero. (absolute smallness does not
-              // suffice since the size of the domain scales in here)
-              //
-              // for this purpose we seek the diagonal of the matrix, where there
-              // must be an element belonging to the boundary. we take the maximum
-              // diagonal entry.
-              types::global_dof_index max_element = static_cast<types::global_dof_index>(0);
-              for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
-                   i!=dof_to_boundary_mapping.end(); ++i)
-                if ((*i != hp::DoFHandler<dim,spacedim>::invalid_dof_index) &&
-                    (*i > max_element))
-                  max_element = *i;
-              Assert (max_element  == matrix.n()-1, ExcInternalError());
-
-              double max_diag_entry = 0;
-              for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
-                if (std::fabs(copy_data.cell_matrix[pos](i,i)) > max_diag_entry)
-                  max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i));
-#endif
-
-              for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
-                for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
-                  {
-                    if (copy_data.dof_is_on_face[pos][i] && copy_data.dof_is_on_face[pos][j])
-                      matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]],
-                                 dof_to_boundary_mapping[copy_data.dofs[j]],
-                                 copy_data.cell_matrix[pos](i,j));
-                    else
-                      {
-                        // assume that all shape functions that are nonzero on the boundary
-                        // are also listed in the @p{dof_to_boundary} mapping. if that
-                        // is not the case, then the boundary mass matrix does not
-                        // make that much sense anyway, as it only contains entries for
-                        // parts of the functions living on the boundary
-                        //
-                        // these, we may compare here for relative smallness of all
-                        // entries in the local matrix which are not taken over to
-                        // the global one
-                        Assert (std::fabs(copy_data.cell_matrix[pos](i,j)) <= 1e-10 * max_diag_entry,
-                                ExcInternalError ());
-                      }
-                  }
-
-              for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
-                if (copy_data.dof_is_on_face[pos][j])
-                  rhs_vector(dof_to_boundary_mapping[copy_data.dofs[j]]) += copy_data.cell_vector[pos](j);
-                else
-                  {
-                    // compare here for relative
-                    // smallness
-                    Assert (std::fabs(copy_data.cell_vector[pos](j)) <= 1e-10 * max_diag_entry,
-                            ExcInternalError());
-                  }
-              ++pos;
-            }
-        }
-    }
-  }
-
-
-
-  template <int dim, int spacedim>
-  void create_boundary_mass_matrix (const DoFHandler<dim,spacedim>     &dof,
-                                    const Quadrature<dim-1>   &q,
-                                    SparseMatrix<double>      &matrix,
-                                    const typename FunctionMap<spacedim>::type &rhs,
-                                    Vector<double>            &rhs_vector,
-                                    std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                                    const Function<spacedim> *const a,
-                                    std::vector<unsigned int> component_mapping)
-  {
-    create_boundary_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q,
-                                matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping);
-  }
-
-
-
-  template <int dim, int spacedim>
-  void
-  create_boundary_mass_matrix (const hp::MappingCollection<dim,spacedim>        &mapping,
-                               const hp::DoFHandler<dim,spacedim>     &dof,
-                               const hp::QCollection<dim-1>   &q,
-                               SparseMatrix<double>      &matrix,
-                               const typename FunctionMap<spacedim>::type         &boundary_functions,
-                               Vector<double>            &rhs_vector,
-                               std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                               const Function<spacedim> *const coefficient,
-                               std::vector<unsigned int> component_mapping)
-  {
-    // what would that be in 1d? the
-    // identity matrix on the boundary
-    // dofs?
-    if (dim == 1)
-      {
-        Assert (false, ExcNotImplemented());
-        return;
-      }
-
-    const hp::FECollection<dim,spacedim> &fe_collection = dof.get_fe();
-    const unsigned int n_components  = fe_collection.n_components();
-
-    Assert (matrix.n() == dof.n_boundary_dofs(boundary_functions),
-            ExcInternalError());
-    Assert (matrix.n() == matrix.m(), ExcInternalError());
-    Assert (matrix.n() == rhs_vector.size(), ExcInternalError());
-    Assert (boundary_functions.size() != 0, ExcInternalError());
-    Assert (dof_to_boundary_mapping.size() == dof.n_dofs(),
-            ExcInternalError());
-    Assert (coefficient ==0 ||
-            coefficient->n_components==1 ||
-            coefficient->n_components==n_components, ExcComponentMismatch());
-
-    if (component_mapping.size() == 0)
-      {
-        AssertDimension (n_components, boundary_functions.begin()->second->n_components);
-        for (unsigned int i=0; i<n_components; ++i)
-          component_mapping.push_back(i);
-      }
-    else
-      AssertDimension (n_components, component_mapping.size());
-
-    MatrixCreator::internal::AssemblerBoundary::Scratch scratch;
-    MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > copy_data;
-
-    WorkStream::run(dof.begin_active(),dof.end(),
-                    static_cast<std_cxx11::function<void (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
-                                                          const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
-                                                          MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > &)> >
-                    (std_cxx11::bind( &create_hp_boundary_mass_matrix_1<dim,spacedim>,std_cxx11::_1,std_cxx11::_2,
-                                      std_cxx11::_3,
-                                      std_cxx11::cref(mapping),std_cxx11::cref(fe_collection),std_cxx11::cref(q),
-                                      std_cxx11::cref(boundary_functions),coefficient,
-                                      std_cxx11::cref(component_mapping))),
-                    static_cast<std_cxx11::function<void (MatrixCreator::internal::AssemblerBoundary
-                                                          ::CopyData<hp::DoFHandler<dim,spacedim> > const &)> > (
-                      std_cxx11::bind( &copy_hp_boundary_mass_matrix_1<dim,spacedim>,
-                                       std_cxx11::_1,
-                                       std_cxx11::cref(boundary_functions),
-                                       std_cxx11::cref(dof_to_boundary_mapping),
-                                       std_cxx11::ref(matrix),
-                                       std_cxx11::ref(rhs_vector))),
-                    scratch,
-                    copy_data);
-  }
-
-
-
-
-  template <int dim, int spacedim>
-  void create_boundary_mass_matrix (const hp::DoFHandler<dim,spacedim>     &dof,
-                                    const hp::QCollection<dim-1>   &q,
-                                    SparseMatrix<double>      &matrix,
-                                    const typename FunctionMap<spacedim>::type &rhs,
-                                    Vector<double>            &rhs_vector,
-                                    std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-                                    const Function<spacedim> *const a,
-                                    std::vector<unsigned int> component_mapping)
-  {
-    create_boundary_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
-                                matrix,rhs, rhs_vector, dof_to_boundary_mapping, a, component_mapping);
-  }
-
-
-
-  template <int dim, int spacedim>
-  void create_laplace_matrix (const Mapping<dim, spacedim>       &mapping,
-                              const DoFHandler<dim,spacedim>    &dof,
-                              const Quadrature<dim>    &q,
-                              SparseMatrix<double>     &matrix,
-                              const Function<spacedim> *const coefficient,
-                              const ConstraintMatrix   &constraints)
-  {
-    Assert (matrix.m() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
-    Assert (matrix.n() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
-
-    hp::FECollection<dim,spacedim>      fe_collection (dof.get_fe());
-    hp::QCollection<dim>                q_collection (q);
-    hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
-    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
-    assembler_data (fe_collection,
-                    update_gradients  | update_JxW_values |
-                    (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)),
-                    coefficient, /*rhs_function=*/0,
-                    q_collection, mapping_collection);
-    MatrixCreator::internal::AssemblerData::CopyData<double> copy_data;
-    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
-                                  assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.constraints = &constraints;
-
-    WorkStream::run (dof.begin_active(),
-                     static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
-                     &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator>,
-                     std_cxx11::bind (&MatrixCreator::internal::
-                                      copy_local_to_global<double,SparseMatrix<double>, Vector<double> >,
-                                      std_cxx11::_1,
-                                      &matrix,
-                                      (Vector<double> *)NULL),
-                     assembler_data,
-                     copy_data);
-  }
-
-
-
-  template <int dim, int spacedim>
-  void create_laplace_matrix (const DoFHandler<dim,spacedim>    &dof,
-                              const Quadrature<dim>    &q,
-                              SparseMatrix<double>     &matrix,
-                              const Function<spacedim> *const coefficient,
-                              const ConstraintMatrix &constraints)
-  {
-    create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping,
-                          dof, q, matrix, coefficient, constraints);
-  }
-
-
-
-  template <int dim, int spacedim>
-  void create_laplace_matrix (const Mapping<dim, spacedim>       &mapping,
-                              const DoFHandler<dim,spacedim>    &dof,
-                              const Quadrature<dim>    &q,
-                              SparseMatrix<double>     &matrix,
-                              const Function<spacedim>      &rhs,
-                              Vector<double>           &rhs_vector,
-                              const Function<spacedim> *const coefficient,
-                              const ConstraintMatrix   &constraints)
-  {
-    Assert (matrix.m() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
-    Assert (matrix.n() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
-
-    hp::FECollection<dim,spacedim>      fe_collection (dof.get_fe());
-    hp::QCollection<dim>                q_collection (q);
-    hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
-    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
-    assembler_data (fe_collection,
-                    update_gradients  | update_values |
-                    update_JxW_values | update_quadrature_points,
-                    coefficient, &rhs,
-                    q_collection, mapping_collection);
-    MatrixCreator::internal::AssemblerData::CopyData<double> copy_data;
-    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
-                                  assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.constraints = &constraints;
-
-    WorkStream::run (dof.begin_active(),
-                     static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
-                     &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator>,
-                     std_cxx11::bind (&MatrixCreator::internal::
-                                      copy_local_to_global<double,SparseMatrix<double>, Vector<double> >,
-                                      std_cxx11::_1,
-                                      &matrix,
-                                      &rhs_vector),
-                     assembler_data,
-                     copy_data);
-  }
-
-
-
-  template <int dim, int spacedim>
-  void create_laplace_matrix (const DoFHandler<dim,spacedim>    &dof,
-                              const Quadrature<dim>    &q,
-                              SparseMatrix<double>     &matrix,
-                              const Function<spacedim>      &rhs,
-                              Vector<double>           &rhs_vector,
-                              const Function<spacedim> *const coefficient,
-                              const ConstraintMatrix &constraints)
-  {
-    create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q,
-                          matrix, rhs, rhs_vector, coefficient, constraints);
-  }
-
-
-
-  template <int dim, int spacedim>
-  void create_laplace_matrix (const hp::MappingCollection<dim,spacedim>       &mapping,
-                              const hp::DoFHandler<dim,spacedim>    &dof,
-                              const hp::QCollection<dim>    &q,
-                              SparseMatrix<double>     &matrix,
-                              const Function<spacedim> *const coefficient,
-                              const ConstraintMatrix   &constraints)
-  {
-    Assert (matrix.m() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
-    Assert (matrix.n() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
-
-    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
-    assembler_data (dof.get_fe(),
-                    update_gradients  | update_JxW_values |
-                    (coefficient != 0 ? update_quadrature_points : UpdateFlags(0)),
-                    coefficient, /*rhs_function=*/0,
-                    q, mapping);
-    MatrixCreator::internal::AssemblerData::CopyData<double> copy_data;
-    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
-                                  assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.constraints = &constraints;
-
-    WorkStream::run (dof.begin_active(),
-                     static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
-                     &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>,
-                     std_cxx11::bind (&MatrixCreator::internal::
-                                      copy_local_to_global<double,SparseMatrix<double>, Vector<double> >,
-                                      std_cxx11::_1,
-                                      &matrix,
-                                      (Vector<double> *)0),
-                     assembler_data,
-                     copy_data);
-  }
-
-
-
-  template <int dim, int spacedim>
-  void create_laplace_matrix (const hp::DoFHandler<dim,spacedim>    &dof,
-                              const hp::QCollection<dim>    &q,
-                              SparseMatrix<double>     &matrix,
-                              const Function<spacedim> *const coefficient,
-                              const ConstraintMatrix &constraints)
-  {
-    create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
-                          matrix, coefficient, constraints);
-  }
-
-
-
-  template <int dim, int spacedim>
-  void create_laplace_matrix (const hp::MappingCollection<dim,spacedim>       &mapping,
-                              const hp::DoFHandler<dim,spacedim>    &dof,
-                              const hp::QCollection<dim>    &q,
-                              SparseMatrix<double>     &matrix,
-                              const Function<spacedim>      &rhs,
-                              Vector<double>           &rhs_vector,
-                              const Function<spacedim> *const coefficient,
-                              const ConstraintMatrix   &constraints)
-  {
-    Assert (matrix.m() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
-    Assert (matrix.n() == dof.n_dofs(),
-            ExcDimensionMismatch (matrix.n(), dof.n_dofs()));
-
-    MatrixCreator::internal::AssemblerData::Scratch<dim, spacedim>
-    assembler_data (dof.get_fe(),
-                    update_gradients  | update_values |
-                    update_JxW_values | update_quadrature_points,
-                    coefficient, &rhs,
-                    q, mapping);
-    MatrixCreator::internal::AssemblerData::CopyData<double> copy_data;
-    copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(),
-                                  assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell());
-    copy_data.constraints = &constraints;
-
-    WorkStream::run (dof.begin_active(),
-                     static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
-                     &MatrixCreator::internal::laplace_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>,
-                     std_cxx11::bind (&MatrixCreator::internal::
-                                      copy_local_to_global<double,SparseMatrix<double>, Vector<double> >,
-                                      std_cxx11::_1,
-                                      &matrix,
-                                      &rhs_vector),
-                     assembler_data,
-                     copy_data);
-  }
-
-
-
-  template <int dim, int spacedim>
-  void create_laplace_matrix (const hp::DoFHandler<dim,spacedim>    &dof,
-                              const hp::QCollection<dim>    &q,
-                              SparseMatrix<double>     &matrix,
-                              const Function<spacedim>      &rhs,
-                              Vector<double>           &rhs_vector,
-                              const Function<spacedim> *const coefficient,
-                              const ConstraintMatrix &constraints)
-  {
-    create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
-                          matrix, rhs, rhs_vector, coefficient, constraints);
-  }
-
-}  // namespace MatrixCreator
-
-
 namespace MatrixTools
 {
   namespace
@@ -2316,462 +489,6 @@ namespace MatrixTools
 
 
 
-#ifdef DEAL_II_WITH_PETSC
-
-  namespace internal
-  {
-    namespace PETScWrappers
-    {
-      template <typename PETScMatrix, typename PETScVector>
-      void
-      apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
-                             PETScMatrix      &matrix,
-                             PETScVector      &solution,
-                             PETScVector      &right_hand_side,
-                             const bool        eliminate_columns)
-      {
-        (void)eliminate_columns;
-        Assert (eliminate_columns == false, ExcNotImplemented());
-
-        Assert (matrix.n() == right_hand_side.size(),
-                ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
-        Assert (matrix.n() == solution.size(),
-                ExcDimensionMismatch(matrix.n(), solution.size()));
-
-        // if no boundary values are to be applied, then
-        // jump straight to the compress() calls that we still have
-        // to perform because they are collective operations
-        if (boundary_values.size() > 0)
-          {
-            const std::pair<types::global_dof_index, types::global_dof_index> local_range
-              = matrix.local_range();
-            Assert (local_range == right_hand_side.local_range(),
-                    ExcInternalError());
-            Assert (local_range == solution.local_range(),
-                    ExcInternalError());
-
-            // determine the first nonzero diagonal
-            // entry from within the part of the
-            // matrix that we can see. if we can't
-            // find such an entry, take one
-            PetscScalar average_nonzero_diagonal_entry = 1;
-            for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
-              if (matrix.diag_element(i) != PetscScalar ())
-                {
-                  average_nonzero_diagonal_entry = std::abs(matrix.diag_element(i));
-                  break;
-                }
-
-            // figure out which rows of the matrix we
-            // have to eliminate on this processor
-            std::vector<types::global_dof_index> constrained_rows;
-            for (std::map<types::global_dof_index,double>::const_iterator
-                 dof  = boundary_values.begin();
-                 dof != boundary_values.end();
-                 ++dof)
-              if ((dof->first >= local_range.first) &&
-                  (dof->first < local_range.second))
-                constrained_rows.push_back (dof->first);
-
-            // then eliminate these rows and set
-            // their diagonal entry to what we have
-            // determined above. note that for petsc
-            // matrices interleaving read with write
-            // operations is very expensive. thus, we
-            // here always replace the diagonal
-            // element, rather than first checking
-            // whether it is nonzero and in that case
-            // preserving it. this is different from
-            // the case of deal.II sparse matrices
-            // treated in the other functions.
-            matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
-
-            std::vector<types::global_dof_index> indices;
-            std::vector<PetscScalar>  solution_values;
-            for (std::map<types::global_dof_index,double>::const_iterator
-                 dof  = boundary_values.begin();
-                 dof != boundary_values.end();
-                 ++dof)
-              if ((dof->first >= local_range.first) &&
-                  (dof->first < local_range.second))
-                {
-                  indices.push_back (dof->first);
-                  solution_values.push_back (dof->second);
-                }
-            solution.set (indices, solution_values);
-
-            // now also set appropriate values for
-            // the rhs
-            for (unsigned int i=0; i<solution_values.size(); ++i)
-              solution_values[i] *= average_nonzero_diagonal_entry;
-
-            right_hand_side.set (indices, solution_values);
-          }
-        else
-          {
-            // clear_rows() is a collective operation so we still have to call
-            // it:
-            std::vector<types::global_dof_index> constrained_rows;
-            matrix.clear_rows (constrained_rows, 1.);
-          }
-
-        // clean up
-        solution.compress (VectorOperation::insert);
-        right_hand_side.compress (VectorOperation::insert);
-      }
-    }
-  }
-
-
-
-  void
-  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
-                         PETScWrappers::SparseMatrix   &matrix,
-                         PETScWrappers::Vector   &solution,
-                         PETScWrappers::Vector   &right_hand_side,
-                         const bool        eliminate_columns)
-  {
-    // simply redirect to the generic function
-    // used for both petsc matrix types
-    internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
-                                                    right_hand_side, eliminate_columns);
-  }
-
-
-
-  void
-  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
-                         PETScWrappers::MPI::SparseMatrix   &matrix,
-                         PETScWrappers::MPI::Vector   &solution,
-                         PETScWrappers::MPI::Vector   &right_hand_side,
-                         const bool        eliminate_columns)
-  {
-    // simply redirect to the generic function
-    // used for both petsc matrix types
-    internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
-                                                    right_hand_side, eliminate_columns);
-  }
-
-
-  void
-  apply_boundary_values (const std::map<types::global_dof_index,double>  &boundary_values,
-                         PETScWrappers::MPI::BlockSparseMatrix &matrix,
-                         PETScWrappers::MPI::BlockVector        &solution,
-                         PETScWrappers::MPI::BlockVector        &right_hand_side,
-                         const bool                            eliminate_columns)
-  {
-    Assert (matrix.n() == right_hand_side.size(),
-            ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
-    Assert (matrix.n() == solution.size(),
-            ExcDimensionMismatch(matrix.n(), solution.size()));
-    Assert (matrix.n_block_rows() == matrix.n_block_cols(),
-            ExcNotQuadratic());
-
-    const unsigned int n_blocks = matrix.n_block_rows();
-
-    // We need to find the subdivision
-    // into blocks for the boundary values.
-    // To this end, generate a vector of
-    // maps with the respective indices.
-    std::vector<std::map<dealii::types::global_dof_index,double> > block_boundary_values(n_blocks);
-    {
-      int block = 0;
-      dealii::types::global_dof_index offset = 0;
-      for (std::map<types::global_dof_index,double>::const_iterator
-           dof  = boundary_values.begin();
-           dof != boundary_values.end();
-           ++dof)
-        {
-          if (dof->first >= matrix.block(block,0).m() + offset)
-            {
-              offset += matrix.block(block,0).m();
-              block++;
-            }
-          const types::global_dof_index index = dof->first - offset;
-          block_boundary_values[block].insert(std::pair<types::global_dof_index, double> (index,dof->second));
-        }
-    }
-
-    // Now call the non-block variants on
-    // the diagonal subblocks and the
-    // solution/rhs.
-    for (unsigned int block=0; block<n_blocks; ++block)
-      internal::PETScWrappers::apply_boundary_values(block_boundary_values[block],
-                                                     matrix.block(block,block),
-                                                     solution.block(block),
-                                                     right_hand_side.block(block),
-                                                     eliminate_columns);
-
-    // Finally, we need to do something
-    // about the off-diagonal matrices. This
-    // is luckily not difficult. Just clear
-    // the whole row.
-    for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
-      {
-        const std::pair<types::global_dof_index, types::global_dof_index> local_range
-          = matrix.block(block_m,0).local_range();
-
-        std::vector<types::global_dof_index> constrained_rows;
-        for (std::map<types::global_dof_index,double>::const_iterator
-             dof  = block_boundary_values[block_m].begin();
-             dof != block_boundary_values[block_m].end();
-             ++dof)
-          if ((dof->first >= local_range.first) &&
-              (dof->first < local_range.second))
-            constrained_rows.push_back (dof->first);
-
-        for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
-          if (block_m != block_n)
-            matrix.block(block_m,block_n).clear_rows(constrained_rows);
-      }
-  }
-
-#endif
-
-
-
-#ifdef DEAL_II_WITH_TRILINOS
-
-  namespace internal
-  {
-    namespace TrilinosWrappers
-    {
-      template <typename TrilinosMatrix, typename TrilinosVector>
-      void
-      apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
-                             TrilinosMatrix      &matrix,
-                             TrilinosVector      &solution,
-                             TrilinosVector      &right_hand_side,
-                             const bool           eliminate_columns)
-      {
-        Assert (eliminate_columns == false, ExcNotImplemented());
-        (void)eliminate_columns;
-
-        Assert (matrix.n() == right_hand_side.size(),
-                ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
-        Assert (matrix.n() == solution.size(),
-                ExcDimensionMismatch(matrix.m(), solution.size()));
-
-        // if no boundary values are to be applied, then
-        // jump straight to the compress() calls that we still have
-        // to perform because they are collective operations
-        if (boundary_values.size() > 0)
-          {
-            const std::pair<types::global_dof_index, types::global_dof_index> local_range
-              = matrix.local_range();
-            Assert (local_range == right_hand_side.local_range(),
-                    ExcInternalError());
-            Assert (local_range == solution.local_range(),
-                    ExcInternalError());
-
-            // determine the first nonzero diagonal
-            // entry from within the part of the
-            // matrix that we can see. if we can't
-            // find such an entry, take one
-            TrilinosScalar average_nonzero_diagonal_entry = 1;
-            for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
-              if (matrix.diag_element(i) != 0)
-                {
-                  average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
-                  break;
-                }
-
-            // figure out which rows of the matrix we
-            // have to eliminate on this processor
-            std::vector<types::global_dof_index> constrained_rows;
-            for (std::map<types::global_dof_index,double>::const_iterator
-                 dof  = boundary_values.begin();
-                 dof != boundary_values.end();
-                 ++dof)
-              if ((dof->first >= local_range.first) &&
-                  (dof->first < local_range.second))
-                constrained_rows.push_back (dof->first);
-
-            // then eliminate these rows and
-            // set their diagonal entry to
-            // what we have determined
-            // above. if the value already is
-            // nonzero, it will be preserved,
-            // in accordance with the basic
-            // matrix classes in deal.II.
-            matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
-
-            std::vector<types::global_dof_index> indices;
-            std::vector<TrilinosScalar>  solution_values;
-            for (std::map<types::global_dof_index,double>::const_iterator
-                 dof  = boundary_values.begin();
-                 dof != boundary_values.end();
-                 ++dof)
-              if ((dof->first >= local_range.first) &&
-                  (dof->first < local_range.second))
-                {
-                  indices.push_back (dof->first);
-                  solution_values.push_back (dof->second);
-                }
-            solution.set (indices, solution_values);
-
-            // now also set appropriate
-            // values for the rhs
-            for (unsigned int i=0; i<solution_values.size(); ++i)
-              solution_values[i] *= matrix.diag_element(indices[i]);
-
-            right_hand_side.set (indices, solution_values);
-          }
-        else
-          {
-            // clear_rows() is a collective operation so we still have to call
-            // it:
-            std::vector<types::global_dof_index> constrained_rows;
-            matrix.clear_rows (constrained_rows, 1.);
-          }
-
-        // clean up
-        matrix.compress (VectorOperation::insert);
-        solution.compress (VectorOperation::insert);
-        right_hand_side.compress (VectorOperation::insert);
-      }
-
-
-
-      template <typename TrilinosMatrix, typename TrilinosBlockVector>
-      void
-      apply_block_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
-                                   TrilinosMatrix      &matrix,
-                                   TrilinosBlockVector &solution,
-                                   TrilinosBlockVector &right_hand_side,
-                                   const bool          eliminate_columns)
-      {
-        Assert (eliminate_columns == false, ExcNotImplemented());
-
-        Assert (matrix.n() == right_hand_side.size(),
-                ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
-        Assert (matrix.n() == solution.size(),
-                ExcDimensionMismatch(matrix.n(), solution.size()));
-        Assert (matrix.n_block_rows() == matrix.n_block_cols(),
-                ExcNotQuadratic());
-
-        const unsigned int n_blocks = matrix.n_block_rows();
-
-        // We need to find the subdivision
-        // into blocks for the boundary values.
-        // To this end, generate a vector of
-        // maps with the respective indices.
-        std::vector<std::map<types::global_dof_index,double> > block_boundary_values(n_blocks);
-        {
-          int block=0;
-          types::global_dof_index offset = 0;
-          for (std::map<types::global_dof_index,double>::const_iterator
-               dof  = boundary_values.begin();
-               dof != boundary_values.end();
-               ++dof)
-            {
-              if (dof->first >= matrix.block(block,0).m() + offset)
-                {
-                  offset += matrix.block(block,0).m();
-                  block++;
-                }
-              const types::global_dof_index index = dof->first - offset;
-              block_boundary_values[block].insert(
-                std::pair<types::global_dof_index, double> (index,dof->second));
-            }
-        }
-
-        // Now call the non-block variants on
-        // the diagonal subblocks and the
-        // solution/rhs.
-        for (unsigned int block=0; block<n_blocks; ++block)
-          TrilinosWrappers::apply_boundary_values(block_boundary_values[block],
-                                                  matrix.block(block,block),
-                                                  solution.block(block),
-                                                  right_hand_side.block(block),
-                                                  eliminate_columns);
-
-        // Finally, we need to do something
-        // about the off-diagonal matrices. This
-        // is luckily not difficult. Just clear
-        // the whole row.
-        for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
-          {
-            const std::pair<types::global_dof_index, types::global_dof_index> local_range
-              = matrix.block(block_m,0).local_range();
-
-            std::vector<types::global_dof_index> constrained_rows;
-            for (std::map<types::global_dof_index,double>::const_iterator
-                 dof  = block_boundary_values[block_m].begin();
-                 dof != block_boundary_values[block_m].end();
-                 ++dof)
-              if ((dof->first >= local_range.first) &&
-                  (dof->first < local_range.second))
-                constrained_rows.push_back (dof->first);
-
-            for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
-              if (block_m != block_n)
-                matrix.block(block_m,block_n).clear_rows(constrained_rows);
-          }
-      }
-    }
-  }
-
-
-
-
-  void
-  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
-                         TrilinosWrappers::SparseMatrix   &matrix,
-                         TrilinosWrappers::Vector         &solution,
-                         TrilinosWrappers::Vector         &right_hand_side,
-                         const bool        eliminate_columns)
-  {
-    // simply redirect to the generic function
-    // used for both trilinos matrix types
-    internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
-                                                       right_hand_side, eliminate_columns);
-  }
-
-
-
-  void
-  apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
-                         TrilinosWrappers::SparseMatrix   &matrix,
-                         TrilinosWrappers::MPI::Vector    &solution,
-                         TrilinosWrappers::MPI::Vector    &right_hand_side,
-                         const bool        eliminate_columns)
-  {
-    // simply redirect to the generic function
-    // used for both trilinos matrix types
-    internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
-                                                       right_hand_side, eliminate_columns);
-  }
-
-
-
-  void
-  apply_boundary_values (const std::map<types::global_dof_index,double>  &boundary_values,
-                         TrilinosWrappers::BlockSparseMatrix &matrix,
-                         TrilinosWrappers::BlockVector        &solution,
-                         TrilinosWrappers::BlockVector        &right_hand_side,
-                         const bool                            eliminate_columns)
-  {
-    internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
-                                                             solution, right_hand_side,
-                                                             eliminate_columns);
-  }
-
-
-
-  void
-  apply_boundary_values (const std::map<types::global_dof_index,double>  &boundary_values,
-                         TrilinosWrappers::BlockSparseMatrix &matrix,
-                         TrilinosWrappers::MPI::BlockVector   &solution,
-                         TrilinosWrappers::MPI::BlockVector   &right_hand_side,
-                         const bool                            eliminate_columns)
-  {
-    internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
-                                                             solution, right_hand_side,
-                                                             eliminate_columns);
-  }
-
-#endif
 
 
 
@@ -2894,37 +611,5 @@ namespace MatrixTools
 // explicit instantiations
 #include "matrix_tools.inst"
 
-namespace MatrixTools
-{
-  template
-  void
-  apply_boundary_values<double> (const std::map<types::global_dof_index,double> &boundary_values,
-                                 SparseMatrix<double>  &matrix,
-                                 Vector<double>   &solution,
-                                 Vector<double>   &right_hand_side,
-                                 const bool        eliminate_columns);
-  template
-  void
-  apply_boundary_values<float> (const std::map<types::global_dof_index,double> &boundary_values,
-                                SparseMatrix<float>  &matrix,
-                                Vector<float>   &solution,
-                                Vector<float>   &right_hand_side,
-                                const bool        eliminate_columns);
-
-  template
-  void
-  apply_boundary_values<double> (const std::map<types::global_dof_index,double> &boundary_values,
-                                 BlockSparseMatrix<double>  &matrix,
-                                 BlockVector<double>   &solution,
-                                 BlockVector<double>   &right_hand_side,
-                                 const bool        eliminate_columns);
-  template
-  void
-  apply_boundary_values<float> (const std::map<types::global_dof_index,double> &boundary_values,
-                                BlockSparseMatrix<float>  &matrix,
-                                BlockVector<float>   &solution,
-                                BlockVector<float>   &right_hand_side,
-                                const bool        eliminate_columns);
-}
 
 DEAL_II_NAMESPACE_CLOSE
index a2e9346dc5bed08b1300eb1bfadd57afb3d9c9c2..3516f400f36e962ec91a8f77567368ea09e15c07 100644 (file)
 //
 // ---------------------------------------------------------------------
 
-for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
-  {
-#if deal_II_dimension <= deal_II_space_dimension
-
-// non-hp version of create_mass_matrix
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
-      (const Mapping<deal_II_dimension,deal_II_space_dimension>       &mapping,
-       const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
-      (const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
-      (const Mapping<deal_II_dimension,deal_II_space_dimension>       &mapping,
-       const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_space_dimension>      &rhs,
-       Vector<double>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
-      (const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_space_dimension>      &rhs,
-       Vector<double>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
-
-    template
-      void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension>
-      (const DoFHandler<deal_II_dimension,deal_II_space_dimension>     &dof,
-       const Quadrature<deal_II_dimension-1>   &q,
-       SparseMatrix<double>      &matrix,
-       const FunctionMap<deal_II_space_dimension>::type &rhs,
-       Vector<double>            &rhs_vector,
-       std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-       const Function<deal_II_space_dimension> * const a,
-       std::vector<unsigned int>);
-
-    template
-      void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension>
-      (const Mapping<deal_II_dimension,deal_II_space_dimension> &,
-       const DoFHandler<deal_II_dimension,deal_II_space_dimension>     &dof,
-       const Quadrature<deal_II_dimension-1>   &q,
-       SparseMatrix<double>      &matrix,
-       const FunctionMap<deal_II_space_dimension>::type &rhs,
-       Vector<double>            &rhs_vector,
-       std::vector<types::global_dof_index> &dof_to_boundary_mapping,
-       const Function<deal_II_space_dimension> * const a,
-       std::vector<unsigned int>);
-
-    template
-      void
-      MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension>
-      (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>&,
-       const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
-       const hp::QCollection<deal_II_dimension-1>&,
-       SparseMatrix<double>&,
-       const FunctionMap<deal_II_space_dimension>::type&,
-       Vector<double>&,
-       std::vector<types::global_dof_index>&,
-       const Function<deal_II_space_dimension> * const,
-       std::vector<unsigned int>);
-
-    template
-      void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension>
-      (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
-       const hp::QCollection<deal_II_dimension-1>&,
-       SparseMatrix<double>&,
-       const FunctionMap<deal_II_space_dimension>::type&,
-       Vector<double>&,
-       std::vector<types::global_dof_index>&,
-       const Function<deal_II_space_dimension> * const,
-       std::vector<unsigned int>);
-
-// same for float
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
-      (const Mapping<deal_II_dimension,deal_II_space_dimension>       &mapping,
-       const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<float>     &matrix,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
-      (const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<float>     &matrix,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
-      (const Mapping<deal_II_dimension,deal_II_space_dimension>       &mapping,
-       const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<float>     &matrix,
-       const Function<deal_II_space_dimension>      &rhs,
-       Vector<float>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
-      (const DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<float>     &matrix,
-       const Function<deal_II_space_dimension>      &rhs,
-       Vector<float>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
-#endif
-  }
-
 
-//TODO[SP]: replace <deal_II_dimension> by <deal_II_dimension, deal_II_space_dimension>
-// where applicable and move to codimension cases above also when applicable
-for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+for (number: REAL_SCALARS)
   {
-// hp versions of functions
-#if deal_II_dimension <= deal_II_space_dimension
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
-      (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
-       const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
     template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
-      (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
-       const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_space_dimension>      &rhs,
-       Vector<double>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
-
-#endif
-
-#if deal_II_dimension == deal_II_space_dimension
+      void MatrixTools::apply_boundary_values
+      (const std::map<types::global_dof_index,double> &boundary_values,
+       SparseMatrix<number>  &matrix,
+       Vector<number>   &solution,
+       Vector<number>   &right_hand_side,
+       const bool        eliminate_columns);
 
-       
     template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension>
-      (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension>
-      (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_space_dimension>      &rhs,
-       Vector<double>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
-// non-hp versions of create_laplace_matrix
-    template
-      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
-      (const DoFHandler<deal_II_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
-      (const Mapping<deal_II_dimension>       &mapping,
-       const DoFHandler<deal_II_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
-      (const Mapping<deal_II_dimension>       &mapping,
-       const DoFHandler<deal_II_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_dimension>      &rhs,
-       Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
-      (const DoFHandler<deal_II_dimension>    &dof,
-       const Quadrature<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_dimension>      &rhs,
-       Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
-// hp versions of create_laplace_matrix
-    template
-      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
-      (const hp::DoFHandler<deal_II_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
-      (const hp::MappingCollection<deal_II_dimension>       &mapping,
-       const hp::DoFHandler<deal_II_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
-      (const hp::MappingCollection<deal_II_dimension>       &mapping,
-       const hp::DoFHandler<deal_II_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_dimension>      &rhs,
-       Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-    template
-      void MatrixCreator::create_laplace_matrix<deal_II_dimension>
-      (const hp::DoFHandler<deal_II_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<double>     &matrix,
-       const Function<deal_II_dimension>      &rhs,
-       Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
-#endif
-
-// same for float
-#if deal_II_dimension <= deal_II_space_dimension
-
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
-      (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
-       const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<float>     &matrix,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
-      (const hp::MappingCollection<deal_II_dimension,deal_II_space_dimension>       &mapping,
-       const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<float>     &matrix,
-       const Function<deal_II_space_dimension>      &rhs,
-       Vector<float>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
-
-#endif
-
-#if deal_II_dimension == deal_II_space_dimension
-
-       
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension>
-      (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<float>     &matrix,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
-    template
-      void MatrixCreator::create_mass_matrix<deal_II_dimension>
-      (const hp::DoFHandler<deal_II_dimension,deal_II_space_dimension>    &dof,
-       const hp::QCollection<deal_II_dimension>    &q,
-       SparseMatrix<float>     &matrix,
-       const Function<deal_II_space_dimension>      &rhs,
-       Vector<float>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient,
-       const ConstraintMatrix   &constraints);
-
-#endif
+      void MatrixTools::apply_boundary_values 
+      (const std::map<types::global_dof_index,double> &boundary_values,
+       BlockSparseMatrix<number>  &matrix,
+       BlockVector<number>   &solution,
+       BlockVector<number>   &right_hand_side,
+       const bool        eliminate_columns);
 
   }
-
diff --git a/source/numerics/matrix_tools_once.cc b/source/numerics/matrix_tools_once.cc
new file mode 100644 (file)
index 0000000..d9bd2d1
--- /dev/null
@@ -0,0 +1,520 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+
+#ifdef DEAL_II_WITH_PETSC
+#  include <deal.II/lac/petsc_parallel_sparse_matrix.h>
+#  include <deal.II/lac/petsc_sparse_matrix.h>
+#  include <deal.II/lac/petsc_parallel_vector.h>
+#  include <deal.II/lac/petsc_vector.h>
+#  include <deal.II/lac/petsc_parallel_block_sparse_matrix.h>
+#endif
+
+#ifdef DEAL_II_WITH_TRILINOS
+#  include <deal.II/lac/trilinos_sparse_matrix.h>
+#  include <deal.II/lac/trilinos_vector.h>
+#  include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#  include <deal.II/lac/trilinos_block_vector.h>
+#endif
+
+#include <algorithm>
+
+
+#include <algorithm>
+#include <set>
+#include <cmath>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+#ifdef DEAL_II_WITH_PETSC
+
+namespace internal
+{
+  namespace PETScWrappers
+  {
+    template <typename PETScMatrix, typename PETScVector>
+    void
+    apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
+                           PETScMatrix      &matrix,
+                           PETScVector      &solution,
+                           PETScVector      &right_hand_side,
+                           const bool        eliminate_columns)
+    {
+      (void)eliminate_columns;
+      Assert (eliminate_columns == false, ExcNotImplemented());
+
+      Assert (matrix.n() == right_hand_side.size(),
+              ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+      Assert (matrix.n() == solution.size(),
+              ExcDimensionMismatch(matrix.n(), solution.size()));
+
+      // if no boundary values are to be applied, then
+      // jump straight to the compress() calls that we still have
+      // to perform because they are collective operations
+      if (boundary_values.size() > 0)
+        {
+          const std::pair<types::global_dof_index, types::global_dof_index> local_range
+            = matrix.local_range();
+          Assert (local_range == right_hand_side.local_range(),
+                  ExcInternalError());
+          Assert (local_range == solution.local_range(),
+                  ExcInternalError());
+
+          // determine the first nonzero diagonal
+          // entry from within the part of the
+          // matrix that we can see. if we can't
+          // find such an entry, take one
+          PetscScalar average_nonzero_diagonal_entry = 1;
+          for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
+            if (matrix.diag_element(i) != PetscScalar ())
+              {
+                average_nonzero_diagonal_entry = std::abs(matrix.diag_element(i));
+                break;
+              }
+
+          // figure out which rows of the matrix we
+          // have to eliminate on this processor
+          std::vector<types::global_dof_index> constrained_rows;
+          for (std::map<types::global_dof_index,double>::const_iterator
+               dof  = boundary_values.begin();
+               dof != boundary_values.end();
+               ++dof)
+            if ((dof->first >= local_range.first) &&
+                (dof->first < local_range.second))
+              constrained_rows.push_back (dof->first);
+
+          // then eliminate these rows and set
+          // their diagonal entry to what we have
+          // determined above. note that for petsc
+          // matrices interleaving read with write
+          // operations is very expensive. thus, we
+          // here always replace the diagonal
+          // element, rather than first checking
+          // whether it is nonzero and in that case
+          // preserving it. this is different from
+          // the case of deal.II sparse matrices
+          // treated in the other functions.
+          matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
+
+          std::vector<types::global_dof_index> indices;
+          std::vector<PetscScalar>  solution_values;
+          for (std::map<types::global_dof_index,double>::const_iterator
+               dof  = boundary_values.begin();
+               dof != boundary_values.end();
+               ++dof)
+            if ((dof->first >= local_range.first) &&
+                (dof->first < local_range.second))
+              {
+                indices.push_back (dof->first);
+                solution_values.push_back (dof->second);
+              }
+          solution.set (indices, solution_values);
+
+          // now also set appropriate values for
+          // the rhs
+          for (unsigned int i=0; i<solution_values.size(); ++i)
+            solution_values[i] *= average_nonzero_diagonal_entry;
+
+          right_hand_side.set (indices, solution_values);
+        }
+      else
+        {
+          // clear_rows() is a collective operation so we still have to call
+          // it:
+          std::vector<types::global_dof_index> constrained_rows;
+          matrix.clear_rows (constrained_rows, 1.);
+        }
+
+      // clean up
+      solution.compress (VectorOperation::insert);
+      right_hand_side.compress (VectorOperation::insert);
+    }
+  }
+}
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
+                       PETScWrappers::SparseMatrix   &matrix,
+                       PETScWrappers::Vector   &solution,
+                       PETScWrappers::Vector   &right_hand_side,
+                       const bool        eliminate_columns)
+{
+  // simply redirect to the generic function
+  // used for both petsc matrix types
+  internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
+                                                  right_hand_side, eliminate_columns);
+}
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
+                       PETScWrappers::MPI::SparseMatrix   &matrix,
+                       PETScWrappers::MPI::Vector   &solution,
+                       PETScWrappers::MPI::Vector   &right_hand_side,
+                       const bool        eliminate_columns)
+{
+  // simply redirect to the generic function
+  // used for both petsc matrix types
+  internal::PETScWrappers::apply_boundary_values (boundary_values, matrix, solution,
+                                                  right_hand_side, eliminate_columns);
+}
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,double>  &boundary_values,
+                       PETScWrappers::MPI::BlockSparseMatrix &matrix,
+                       PETScWrappers::MPI::BlockVector        &solution,
+                       PETScWrappers::MPI::BlockVector        &right_hand_side,
+                       const bool                            eliminate_columns)
+{
+  Assert (matrix.n() == right_hand_side.size(),
+          ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+  Assert (matrix.n() == solution.size(),
+          ExcDimensionMismatch(matrix.n(), solution.size()));
+  Assert (matrix.n_block_rows() == matrix.n_block_cols(),
+          ExcNotQuadratic());
+
+  const unsigned int n_blocks = matrix.n_block_rows();
+
+  // We need to find the subdivision
+  // into blocks for the boundary values.
+  // To this end, generate a vector of
+  // maps with the respective indices.
+  std::vector<std::map<dealii::types::global_dof_index,double> > block_boundary_values(n_blocks);
+  {
+    int block = 0;
+    dealii::types::global_dof_index offset = 0;
+    for (std::map<types::global_dof_index,double>::const_iterator
+         dof  = boundary_values.begin();
+         dof != boundary_values.end();
+         ++dof)
+      {
+        if (dof->first >= matrix.block(block,0).m() + offset)
+          {
+            offset += matrix.block(block,0).m();
+            block++;
+          }
+        const types::global_dof_index index = dof->first - offset;
+        block_boundary_values[block].insert(std::pair<types::global_dof_index, double> (index,dof->second));
+      }
+  }
+
+  // Now call the non-block variants on
+  // the diagonal subblocks and the
+  // solution/rhs.
+  for (unsigned int block=0; block<n_blocks; ++block)
+    internal::PETScWrappers::apply_boundary_values(block_boundary_values[block],
+                                                   matrix.block(block,block),
+                                                   solution.block(block),
+                                                   right_hand_side.block(block),
+                                                   eliminate_columns);
+
+  // Finally, we need to do something
+  // about the off-diagonal matrices. This
+  // is luckily not difficult. Just clear
+  // the whole row.
+  for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
+    {
+      const std::pair<types::global_dof_index, types::global_dof_index> local_range
+        = matrix.block(block_m,0).local_range();
+
+      std::vector<types::global_dof_index> constrained_rows;
+      for (std::map<types::global_dof_index,double>::const_iterator
+           dof  = block_boundary_values[block_m].begin();
+           dof != block_boundary_values[block_m].end();
+           ++dof)
+        if ((dof->first >= local_range.first) &&
+            (dof->first < local_range.second))
+          constrained_rows.push_back (dof->first);
+
+      for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
+        if (block_m != block_n)
+          matrix.block(block_m,block_n).clear_rows(constrained_rows);
+    }
+}
+
+#endif
+
+
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+namespace internal
+{
+  namespace TrilinosWrappers
+  {
+    template <typename TrilinosMatrix, typename TrilinosVector>
+    void
+    apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
+                           TrilinosMatrix      &matrix,
+                           TrilinosVector      &solution,
+                           TrilinosVector      &right_hand_side,
+                           const bool           eliminate_columns)
+    {
+      Assert (eliminate_columns == false, ExcNotImplemented());
+      (void)eliminate_columns;
+
+      Assert (matrix.n() == right_hand_side.size(),
+              ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+      Assert (matrix.n() == solution.size(),
+              ExcDimensionMismatch(matrix.m(), solution.size()));
+
+      // if no boundary values are to be applied, then
+      // jump straight to the compress() calls that we still have
+      // to perform because they are collective operations
+      if (boundary_values.size() > 0)
+        {
+          const std::pair<types::global_dof_index, types::global_dof_index> local_range
+            = matrix.local_range();
+          Assert (local_range == right_hand_side.local_range(),
+                  ExcInternalError());
+          Assert (local_range == solution.local_range(),
+                  ExcInternalError());
+
+          // determine the first nonzero diagonal
+          // entry from within the part of the
+          // matrix that we can see. if we can't
+          // find such an entry, take one
+          TrilinosScalar average_nonzero_diagonal_entry = 1;
+          for (types::global_dof_index i=local_range.first; i<local_range.second; ++i)
+            if (matrix.diag_element(i) != 0)
+              {
+                average_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
+                break;
+              }
+
+          // figure out which rows of the matrix we
+          // have to eliminate on this processor
+          std::vector<types::global_dof_index> constrained_rows;
+          for (std::map<types::global_dof_index,double>::const_iterator
+               dof  = boundary_values.begin();
+               dof != boundary_values.end();
+               ++dof)
+            if ((dof->first >= local_range.first) &&
+                (dof->first < local_range.second))
+              constrained_rows.push_back (dof->first);
+
+          // then eliminate these rows and
+          // set their diagonal entry to
+          // what we have determined
+          // above. if the value already is
+          // nonzero, it will be preserved,
+          // in accordance with the basic
+          // matrix classes in deal.II.
+          matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
+
+          std::vector<types::global_dof_index> indices;
+          std::vector<TrilinosScalar>  solution_values;
+          for (std::map<types::global_dof_index,double>::const_iterator
+               dof  = boundary_values.begin();
+               dof != boundary_values.end();
+               ++dof)
+            if ((dof->first >= local_range.first) &&
+                (dof->first < local_range.second))
+              {
+                indices.push_back (dof->first);
+                solution_values.push_back (dof->second);
+              }
+          solution.set (indices, solution_values);
+
+          // now also set appropriate
+          // values for the rhs
+          for (unsigned int i=0; i<solution_values.size(); ++i)
+            solution_values[i] *= matrix.diag_element(indices[i]);
+
+          right_hand_side.set (indices, solution_values);
+        }
+      else
+        {
+          // clear_rows() is a collective operation so we still have to call
+          // it:
+          std::vector<types::global_dof_index> constrained_rows;
+          matrix.clear_rows (constrained_rows, 1.);
+        }
+
+      // clean up
+      matrix.compress (VectorOperation::insert);
+      solution.compress (VectorOperation::insert);
+      right_hand_side.compress (VectorOperation::insert);
+    }
+
+
+
+    template <typename TrilinosMatrix, typename TrilinosBlockVector>
+    void
+    apply_block_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
+                                 TrilinosMatrix      &matrix,
+                                 TrilinosBlockVector &solution,
+                                 TrilinosBlockVector &right_hand_side,
+                                 const bool          eliminate_columns)
+    {
+      Assert (eliminate_columns == false, ExcNotImplemented());
+
+      Assert (matrix.n() == right_hand_side.size(),
+              ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+      Assert (matrix.n() == solution.size(),
+              ExcDimensionMismatch(matrix.n(), solution.size()));
+      Assert (matrix.n_block_rows() == matrix.n_block_cols(),
+              ExcNotQuadratic());
+
+      const unsigned int n_blocks = matrix.n_block_rows();
+
+      // We need to find the subdivision
+      // into blocks for the boundary values.
+      // To this end, generate a vector of
+      // maps with the respective indices.
+      std::vector<std::map<types::global_dof_index,double> > block_boundary_values(n_blocks);
+      {
+        int block=0;
+        types::global_dof_index offset = 0;
+        for (std::map<types::global_dof_index,double>::const_iterator
+             dof  = boundary_values.begin();
+             dof != boundary_values.end();
+             ++dof)
+          {
+            if (dof->first >= matrix.block(block,0).m() + offset)
+              {
+                offset += matrix.block(block,0).m();
+                block++;
+              }
+            const types::global_dof_index index = dof->first - offset;
+            block_boundary_values[block].insert(
+              std::pair<types::global_dof_index, double> (index,dof->second));
+          }
+      }
+
+      // Now call the non-block variants on
+      // the diagonal subblocks and the
+      // solution/rhs.
+      for (unsigned int block=0; block<n_blocks; ++block)
+        TrilinosWrappers::apply_boundary_values(block_boundary_values[block],
+                                                matrix.block(block,block),
+                                                solution.block(block),
+                                                right_hand_side.block(block),
+                                                eliminate_columns);
+
+      // Finally, we need to do something
+      // about the off-diagonal matrices. This
+      // is luckily not difficult. Just clear
+      // the whole row.
+      for (unsigned int block_m=0; block_m<n_blocks; ++block_m)
+        {
+          const std::pair<types::global_dof_index, types::global_dof_index> local_range
+            = matrix.block(block_m,0).local_range();
+
+          std::vector<types::global_dof_index> constrained_rows;
+          for (std::map<types::global_dof_index,double>::const_iterator
+               dof  = block_boundary_values[block_m].begin();
+               dof != block_boundary_values[block_m].end();
+               ++dof)
+            if ((dof->first >= local_range.first) &&
+                (dof->first < local_range.second))
+              constrained_rows.push_back (dof->first);
+
+          for (unsigned int block_n=0; block_n<n_blocks; ++block_n)
+            if (block_m != block_n)
+              matrix.block(block_m,block_n).clear_rows(constrained_rows);
+        }
+    }
+  }
+}
+
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
+                       TrilinosWrappers::SparseMatrix   &matrix,
+                       TrilinosWrappers::Vector         &solution,
+                       TrilinosWrappers::Vector         &right_hand_side,
+                       const bool        eliminate_columns)
+{
+  // simply redirect to the generic function
+  // used for both trilinos matrix types
+  internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
+                                                     right_hand_side, eliminate_columns);
+}
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,double> &boundary_values,
+                       TrilinosWrappers::SparseMatrix   &matrix,
+                       TrilinosWrappers::MPI::Vector    &solution,
+                       TrilinosWrappers::MPI::Vector    &right_hand_side,
+                       const bool        eliminate_columns)
+{
+  // simply redirect to the generic function
+  // used for both trilinos matrix types
+  internal::TrilinosWrappers::apply_boundary_values (boundary_values, matrix, solution,
+                                                     right_hand_side, eliminate_columns);
+}
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,double>  &boundary_values,
+                       TrilinosWrappers::BlockSparseMatrix &matrix,
+                       TrilinosWrappers::BlockVector        &solution,
+                       TrilinosWrappers::BlockVector        &right_hand_side,
+                       const bool                            eliminate_columns)
+{
+  internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
+                                                           solution, right_hand_side,
+                                                           eliminate_columns);
+}
+
+
+
+void
+apply_boundary_values (const std::map<types::global_dof_index,double>  &boundary_values,
+                       TrilinosWrappers::BlockSparseMatrix &matrix,
+                       TrilinosWrappers::MPI::BlockVector   &solution,
+                       TrilinosWrappers::MPI::BlockVector   &right_hand_side,
+                       const bool                            eliminate_columns)
+{
+  internal::TrilinosWrappers::apply_block_boundary_values (boundary_values, matrix,
+                                                           solution, right_hand_side,
+                                                           eliminate_columns);
+}
+
+#endif
+
+
+DEAL_II_NAMESPACE_CLOSE

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.