]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Extend MatrixTools::create_mass_matrix and MatrixTools::create_laplace_matrix.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 9 May 2013 12:55:58 +0000 (12:55 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 9 May 2013 12:55:58 +0000 (12:55 +0000)
git-svn-id: https://svn.dealii.org/trunk@29486 0785d39b-7218-0410-832d-ea1e28bc413d

33 files changed:
deal.II/doc/news/changes.h
deal.II/include/deal.II/numerics/matrix_tools.h
deal.II/include/deal.II/numerics/vector_tools.templates.h
deal.II/source/numerics/matrix_tools.cc
deal.II/source/numerics/matrix_tools.inst.in
tests/deal.II/create_laplace_matrix_constraints_01.cc [new file with mode: 0644]
tests/deal.II/create_laplace_matrix_constraints_01/cmp/generic [new file with mode: 0644]
tests/deal.II/create_laplace_matrix_constraints_02.cc [new file with mode: 0644]
tests/deal.II/create_laplace_matrix_constraints_02/cmp/generic [new file with mode: 0644]
tests/deal.II/create_laplace_matrix_constraints_02b.cc [new file with mode: 0644]
tests/deal.II/create_laplace_matrix_constraints_02b/cmp/generic [new file with mode: 0644]
tests/deal.II/create_laplace_matrix_constraints_03.cc [new file with mode: 0644]
tests/deal.II/create_laplace_matrix_constraints_03/cmp/generic [new file with mode: 0644]
tests/deal.II/create_laplace_matrix_constraints_04.cc [new file with mode: 0644]
tests/deal.II/create_laplace_matrix_constraints_04/cmp/generic [new file with mode: 0644]
tests/deal.II/create_mass_matrix_constraints_01.cc [new file with mode: 0644]
tests/deal.II/create_mass_matrix_constraints_01/cmp/generic [new file with mode: 0644]
tests/deal.II/create_mass_matrix_constraints_02.cc [new file with mode: 0644]
tests/deal.II/create_mass_matrix_constraints_02/cmp/generic [new file with mode: 0644]
tests/deal.II/create_mass_matrix_constraints_02b.cc [new file with mode: 0644]
tests/deal.II/create_mass_matrix_constraints_02b/cmp/generic [new file with mode: 0644]
tests/deal.II/create_mass_matrix_constraints_03.cc [new file with mode: 0644]
tests/deal.II/create_mass_matrix_constraints_03/cmp/generic [new file with mode: 0644]
tests/deal.II/create_mass_matrix_constraints_04.cc [new file with mode: 0644]
tests/deal.II/create_mass_matrix_constraints_04/cmp/generic [new file with mode: 0644]
tests/hp/create_laplace_matrix_constraints_02.cc [new file with mode: 0644]
tests/hp/create_laplace_matrix_constraints_02/cmp/generic [new file with mode: 0644]
tests/hp/create_laplace_matrix_constraints_02b.cc [new file with mode: 0644]
tests/hp/create_laplace_matrix_constraints_02b/cmp/generic [new file with mode: 0644]
tests/hp/create_mass_matrix_constraints_02.cc [new file with mode: 0644]
tests/hp/create_mass_matrix_constraints_02/cmp/generic [new file with mode: 0644]
tests/hp/create_mass_matrix_constraints_02b.cc [new file with mode: 0644]
tests/hp/create_mass_matrix_constraints_02b/cmp/generic [new file with mode: 0644]

index ace3025a0996bed73679f690e9ce89c5190a85fc..26c06d08d1a3b7c809bb3087eb3cc681542f5d46 100644 (file)
@@ -117,6 +117,15 @@ this function.
 
 <ol>
 
+<li> Improved: The functions MatrixTools::create_mass_matrix and
+MatrixTools::create_laplace_matrix take now an optional ConstraintMatrix
+argument that allows to directly apply the constraints. This also helps
+VectorTools::project. Note that not providing constraints remains the default
+and recommended way to ensure consistency when several matrices are added.
+<br>
+(Martin Kronbichler, 2013/05/08)
+</li>
+
 <li> New: The classes TrilinosWrappers::SparseMatrix and
 TrilinosWrappers::BlockSparseMatrix now fully implement vmult and Tvmult with
 deal.II's own vector classes Vector<double> and
index 38be9336e27a9a4abd2f6b779c09ba9bf3637d22..abb91c8194bc2651fdd119b98caa74cd3f64d94d 100644 (file)
@@ -26,6 +26,8 @@ DEAL_II_NAMESPACE_OPEN
 // forward declarations
 template <int dim> class Quadrature;
 
+class ConstraintMatrix;
+
 template<typename number> class Vector;
 template<typename number> class FullMatrix;
 template<typename number> class SparseMatrix;
@@ -105,17 +107,27 @@ namespace TrilinosWrappers
  * the matrix wasn't empty before. Therefore you may want to clear the matrix
  * before assemblage.
  *
- * All created matrices are `raw': they are not condensed,
+ * By default, all created matrices are `raw': they are not condensed,
  * i.e. hanging nodes are not eliminated. The reason is that you may
  * want to add several matrices and could then condense afterwards
  * only once, instead of for every matrix. To actually do computations
  * with these matrices, you have to condense the matrix using the
  * ConstraintMatrix::condense function; you also have to
  * condense the right hand side accordingly and distribute the
- * solution afterwards.
+ * solution afterwards. Alternatively, you can give an optional argument
+ * ConstraintMatrix that writes cell matrix (and vector) entries with
+ * distribute_local_to_global into the global matrix and vector. This way,
+ * adding several matrices from different sources is more complicated and
+ * you should make sure that you do not mix different ways of applying
+ * constraints. Particular caution is necessary when the given
+ * constraint matrix contains inhomogeneous constraints: In that case, the
+ * matrix assembled this way must be the only matrix (or you need to
+ * assemble the <b>same</b> right hand side for <b>every</b> matrix
+ * you generate and add together).
  *
  * If you want to use boundary conditions with the matrices generated
- * by the functions of this class, you have to use a function like
+ * by the functions of this class in addition to the ones in a possible
+ * constraint matrix, you have to use a function like
  * <tt>ProblemBase<>::apply_dirichlet_bc</tt> to matrix and right hand
  * side.
  *
@@ -210,154 +222,141 @@ namespace TrilinosWrappers
 namespace MatrixCreator
 {
   /**
-   * Assemble the mass matrix. If no
-   * coefficient is given, it is assumed
-   * to be unity.
+   * Assemble the mass matrix. If no coefficient is given, it is assumed to be
+   * unity.
+   *
+   * If the library is configured to use multithreading, this function works
+   * in parallel.
    *
-   * If the library is configured
-   * to use multithreading, this
-   * function works in parallel.
+   * The optional argument @p constraints allows to apply constraints on the
+   * resulting matrix directly. Be careful when combining several matrices and
+   * using inhomogeneous constraints.
    *
-   * See the general doc of this class
-   * for more information.
+   * See the general doc of this class for more information.
    */
   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 a = 0);
+                           const Function<spacedim> *const a = 0,
+                           const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Calls the create_mass_matrix()
-   * function, see above, with
+   * Calls the create_mass_matrix() function, see above, with
    * <tt>mapping=MappingQ1@<dim@>()</tt>.
    */
   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 a = 0);
+                           const Function<spacedim> *const a = 0,
+                           const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Assemble the mass matrix and a
-   * right hand side vector. If no
-   * coefficient is given, it is
-   * assumed to be unity.
+   * Assemble the mass matrix and a right hand side vector. If no coefficient
+   * is given, it is assumed to be unity.
+   *
+   * If the library is configured to use multithreading, this function works
+   * in parallel.
    *
-   * If the library is configured
-   * to use multithreading, this
-   * function works in parallel.
+   * The optional argument @p constraints allows to apply constraints on the
+   * resulting matrix directly. Be careful when combining several matrices and
+   * using inhomogeneous constraints.
    *
-   * See the general doc of this
-   * class for more information.
+   * See the general doc of this class for more information.
    */
   template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const Mapping<dim, spacedim>       &mapping,
-                           const DoFHandler<dim,spacedim>    &dof,
+  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,
+                           const Function<spacedim> &rhs,
                            Vector<double>           &rhs_vector,
-                           const Function<spacedim> *const a = 0);
+                           const Function<spacedim> *const a = 0,
+                           const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Calls the create_mass_matrix()
-   * function, see above, with
+   * Calls the create_mass_matrix() function, see above, with
    * <tt>mapping=MappingQ1@<dim@>()</tt>.
    */
   template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const DoFHandler<dim,spacedim>    &dof,
+  void create_mass_matrix (const DoFHandler<dim,spacedim> &dof,
                            const Quadrature<dim>    &q,
                            SparseMatrix<number>     &matrix,
-                           const Function<spacedim>      &rhs,
+                           const Function<spacedim> &rhs,
                            Vector<double>           &rhs_vector,
-                           const Function<spacedim> *const a = 0);
+                           const Function<spacedim> *const a = 0,
+                           const ConstraintMatrix   &constraints = ConstraintMatrix());
+
   /**
-   * Same function as above, but for hp
-   * objects.
+   * Same function as above, but for hp objects.
    */
   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 a = 0);
+                           const Function<spacedim> *const a = 0,
+                           const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Same function as above, but for hp
-   * objects.
+   * Same function as above, but for hp objects.
    */
   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 a = 0);
+                           const Function<spacedim> *const a = 0,
+                           const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Same function as above, but for hp
-   * objects.
+   * Same function as above, but for hp objects.
    */
   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,
+  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,
+                           const Function<spacedim> &rhs,
                            Vector<double>           &rhs_vector,
-                           const Function<spacedim> *const a = 0);
+                           const Function<spacedim> *const a = 0,
+                           const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Same function as above, but for hp
-   * objects.
+   * Same function as above, but for hp objects.
    */
   template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const hp::DoFHandler<dim,spacedim>    &dof,
-                           const hp::QCollection<dim>    &q,
+  void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
+                           const hp::QCollection<dim> &q,
                            SparseMatrix<number>     &matrix,
-                           const Function<spacedim>      &rhs,
+                           const Function<spacedim> &rhs,
                            Vector<double>           &rhs_vector,
-                           const Function<spacedim> *const a = 0);
+                           const Function<spacedim> *const a = 0,
+                           const ConstraintMatrix   &constraints = ConstraintMatrix());
 
 
   /**
-   * Assemble the mass matrix and a
-   * right hand side vector along
-   * the boundary.
+   * Assemble the mass matrix and a right hand side vector along the boundary.
    *
-   * The matrix is assumed to
-   * already be initialized with a
-   * suiting sparsity pattern (the
-   * DoFHandler provides an
-   * appropriate function).
+   * The matrix is assumed to already be initialized with a suiting sparsity
+   * pattern (the DoFHandler provides an appropriate function).
    *
-   * If the library is configured
-   * to use multithreading, this
-   * function works in parallel.
+   * If the library is configured to use multithreading, this function works
+   * in parallel.
    *
-   * @arg @p weight: an optional
-   * weight for the computation of
-   * the mass matrix. If no weight
-   * is given, it is set to one.
+   * @arg @p weight: an optional weight for the computation of the mass
+   * matrix. If no weight is given, it is set to one.
    *
-   * @arg @p component_mapping: if
-   * the components in @p
-   * boundary_functions and @p dof
-   * do not coincide, this vector
-   * allows them to be remapped. If
-   * the vector is not empty, it
-   * has to have one entry for each
-   * component in @p dof. This
-   * entry is the component number
-   * in @p boundary_functions that
-   * should be used for this
-   * component in @p dof. By
-   * default, no remapping is
+   * @arg @p component_mapping: if the components in @p boundary_functions and
+   * @p dof do not coincide, this vector allows them to be remapped. If the
+   * vector is not empty, it has to have one entry for each component in @p
+   * dof. This entry is the component number in @p boundary_functions that
+   * should be used for this component in @p dof. By default, no remapping is
    * applied.
    *
-   * @todo This function does not
-   * work for finite elements with
-   * cell-dependent shape functions.
+   * @todo This function does not work for finite elements with cell-dependent
+   * shape functions.
    */
   template <int dim, int spacedim>
 
@@ -396,9 +395,7 @@ namespace MatrixCreator
 
 
   /**
-   * Calls the
-   * create_boundary_mass_matrix()
-   * function, see above, with
+   * Calls the create_boundary_mass_matrix() function, see above, with
    * <tt>mapping=MappingQ1@<dim@>()</tt>.
    */
   template <int dim, int spacedim>
@@ -413,8 +410,7 @@ namespace MatrixCreator
                                     std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
 
   /**
-   * Same function as above, but for hp
-   * objects.
+   * Same function as above, but for hp objects.
    */
   template <int dim, int spacedim>
 
@@ -429,8 +425,7 @@ namespace MatrixCreator
                                     std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
 
   /**
-   * Same function as above, but for hp
-   * objects.
+   * Same function as above, but for hp objects.
    */
 //
 //     void create_boundary_mass_matrix (const hp::MappingCollection<1,1>       &mapping,
@@ -443,8 +438,7 @@ namespace MatrixCreator
 //                                    const Function<1> * const a = 0);
 
   /**
-   * Same function as above, but for hp
-   * objects.
+   * Same function as above, but for hp objects.
    */
   template <int dim, int spacedim>
 
@@ -458,123 +452,122 @@ namespace MatrixCreator
                                     std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
 
   /**
-   * Assemble the Laplace
-   * matrix. If no coefficient is
-   * given, it is assumed to be
-   * constant one.
+   * Assemble the Laplace matrix. If no coefficient is given, it is assumed to
+   * be constant one.
+   *
+   * If the library is configured to use multithreading, this function works
+   * in parallel.
    *
-   * If the library is configured
-   * to use multithreading, this
-   * function works in parallel.
+   * The optional argument @p constraints allows to apply constraints on the
+   * resulting matrix directly. Be careful when combining several matrices and
+   * using inhomogeneous constraints.
    *
-   * See the general doc of this
-   * class for more information.
+   * See the general doc of this class for more information.
    */
   template <int dim, int spacedim>
-  void create_laplace_matrix (const Mapping<dim, spacedim>       &mapping,
-                              const DoFHandler<dim,spacedim>    &dof,
+  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 a = 0);
+                              const Function<spacedim> *const a = 0,
+                              const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Calls the
-   * create_laplace_matrix()
-   * function, see above, with
+   * Calls the create_laplace_matrix() function, see above, with
    * <tt>mapping=MappingQ1@<dim@>()</tt>.
    */
   template <int dim, int spacedim>
-  void create_laplace_matrix (const DoFHandler<dim,spacedim>    &dof,
+  void create_laplace_matrix (const DoFHandler<dim,spacedim> &dof,
                               const Quadrature<dim>    &q,
                               SparseMatrix<double>     &matrix,
-                              const Function<spacedim> *const a = 0);
+                              const Function<spacedim> *const a = 0,
+                              const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Assemble the Laplace matrix
-   * and a right hand side
-   * vector. If no coefficient is
-   * given, it is assumed to be
-   * constant one.
+   * Assemble the Laplace matrix and a right hand side vector. If no
+   * coefficient is given, it is assumed to be constant one.
+   *
+   * If the library is configured to use multithreading, this function works
+   * in parallel.
    *
-   * If the library is configured
-   * to use multithreading, this
-   * function works in parallel.
+   * The optional argument @p constraints allows to apply constraints on the
+   * resulting matrix directly. Be careful when combining several matrices and
+   * using inhomogeneous constraints.
    *
-   * See the general doc of this
-   * class for more information.
+   * See the general doc of this class for more information.
    */
   template <int dim, int spacedim>
-  void create_laplace_matrix (const Mapping<dim, spacedim>       &mapping,
-                              const DoFHandler<dim,spacedim>    &dof,
+  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,
+                              const Function<spacedim> &rhs,
                               Vector<double>           &rhs_vector,
-                              const Function<spacedim> *const a = 0);
+                              const Function<spacedim> *const a = 0,
+                              const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Calls the
-   * create_laplace_matrix()
-   * function, see above, with
+   * Calls the create_laplace_matrix() function, see above, with
    * <tt>mapping=MappingQ1@<dim@>()</tt>.
    */
   template <int dim, int spacedim>
-  void create_laplace_matrix (const DoFHandler<dim,spacedim>    &dof,
+  void create_laplace_matrix (const DoFHandler<dim,spacedim> &dof,
                               const Quadrature<dim>    &q,
                               SparseMatrix<double>     &matrix,
-                              const Function<spacedim>      &rhs,
+                              const Function<spacedim> &rhs,
                               Vector<double>           &rhs_vector,
-                              const Function<spacedim> *const a = 0);
+                              const Function<spacedim> *const a = 0,
+                              const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Like the functions above, but for hp
-   * dof handlers, mappings, and quadrature
-   * collections.
+   * Like the functions above, but for hp dof handlers, mappings, and
+   * quadrature collections.
    */
   template <int dim, int spacedim>
-  void create_laplace_matrix (const hp::MappingCollection<dim,spacedim>       &mapping,
-                              const hp::DoFHandler<dim,spacedim>    &dof,
+  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 a = 0);
+                              const Function<spacedim> *const a = 0,
+                              const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Like the functions above, but for hp
-   * dof handlers, mappings, and quadrature
-   * collections.
+   * Like the functions above, but for hp dof handlers, mappings, and
+   * quadrature collections.
    */
   template <int dim, int spacedim>
-  void create_laplace_matrix (const hp::DoFHandler<dim,spacedim>    &dof,
+  void create_laplace_matrix (const hp::DoFHandler<dim,spacedim> &dof,
                               const hp::QCollection<dim>    &q,
                               SparseMatrix<double>     &matrix,
-                              const Function<spacedim> *const a = 0);
+                              const Function<spacedim> *const a = 0,
+                              const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Like the functions above, but for hp
-   * dof handlers, mappings, and quadrature
-   * collections.
+   * Like the functions above, but for hp dof handlers, mappings, and
+   * quadrature collections.
    */
   template <int dim, int spacedim>
-  void create_laplace_matrix (const hp::MappingCollection<dim,spacedim>       &mapping,
-                              const hp::DoFHandler<dim,spacedim>    &dof,
+  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 a = 0);
+                              const Function<spacedim> *const a = 0,
+                              const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
-   * Like the functions above, but for hp
-   * dof handlers, mappings, and quadrature
-   * collections.
+   * Like the functions above, but for hp dof handlers, mappings, and
+   * quadrature collections.
    */
   template <int dim, int spacedim>
-  void create_laplace_matrix (const hp::DoFHandler<dim,spacedim>    &dof,
+  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 a = 0);
+                              const Function<spacedim> *const a = 0,
+                              const ConstraintMatrix   &constraints = ConstraintMatrix());
 
   /**
    * Exception
index 129cc5f779dacacc40e0018b34fac0722fc00ae3..3dec6628b36ed2d659121e02570d883ca24fdf0e 100644 (file)
@@ -557,11 +557,11 @@ namespace VectorTools
 
 
   template <int dim, class VECTOR, int spacedim>
-  void project (const Mapping<dim, spacedim>       &mapping,
-                const DoFHandler<dim,spacedim>    &dof,
+  void project (const Mapping<dim, spacedim>   &mapping,
+                const DoFHandler<dim,spacedim> &dof,
                 const ConstraintMatrix   &constraints,
                 const Quadrature<dim>    &quadrature,
-                const Function<spacedim>      &function,
+                const Function<spacedim> &function,
                 VECTOR                   &vec_result,
                 const bool                enforce_zero_boundary,
                 const Quadrature<dim-1>  &q_boundary,
@@ -607,11 +607,27 @@ namespace VectorTools
     Vector<double> vec (dof.n_dofs());
     SparsityPattern sparsity;
 
+    // check if constraints are compatible (see below)
+    bool constraints_are_compatible = true;
+    {
+      for (std::map<unsigned int,double>::iterator it=boundary_values.begin();
+           it != boundary_values.end(); ++it)
+        if (constraints.is_constrained(it->first))
+          if (!(constraints.get_constraint_entries(it->first)->size() > 0
+                ||
+                (constraints.get_inhomogeneity(it->first) == it->second)))
+            {
+              constraints_are_compatible = false;
+              break;
+            }
+    }
+
     // use csp to consume less memory and to
     // still be fast
     {
       CompressedSimpleSparsityPattern csp (dof.n_dofs(), dof.n_dofs());
-      DoFTools::make_sparsity_pattern (dof, csp, constraints);
+      DoFTools::make_sparsity_pattern (dof, csp, constraints,
+                                       !constraints_are_compatible);
 
       sparsity.copy_from (csp);
     }
@@ -619,52 +635,60 @@ namespace VectorTools
     SparseMatrix<double> mass_matrix (sparsity);
     Vector<double> tmp (mass_matrix.n());
 
-    // create mass matrix and rhs at once,
-    // which is faster.
-    MatrixCreator::create_mass_matrix (mapping, dof, quadrature, mass_matrix,
-                                       function, tmp);
-
-    constraints.condense (mass_matrix);
-    constraints.condense (tmp);
-    if (boundary_values.size() != 0)
-      MatrixTools::apply_boundary_values (boundary_values,
-                                          mass_matrix, vec, tmp,
-                                          true);
-    // Allow for a maximum of 5*n
-    // steps to reduce the residual by
-    // 10^-12. n steps may not be
-    // sufficient, since roundoff
-    // errors may accumulate for badly
-    // conditioned matrices
-    ReductionControl        control(5*tmp.size(), 0., 1e-12, false, false);
+    // If the constraint matrix does not conflict with the given boundary
+    // values (i.e., it either does not contain boundary values or it contains
+    // the same as boundary_values), we can let it call
+    // distribute_local_to_global straight away, otherwise we need to first
+    // interpolate the boundary values and then condense the matrix and vector
+    if (constraints_are_compatible)
+      {
+        const Function<spacedim>* dummy = 0;
+        MatrixCreator::create_mass_matrix (mapping, dof, quadrature,
+                                           mass_matrix, function, tmp,
+                                           dummy, constraints);
+        if (boundary_values.size() > 0)
+          MatrixTools::apply_boundary_values (boundary_values,
+                                              mass_matrix, vec, tmp,
+                                              true);
+      }
+    else
+      {
+        // create mass matrix and rhs at once, which is faster.
+        MatrixCreator::create_mass_matrix (mapping, dof, quadrature,
+                                           mass_matrix, function, tmp);
+        MatrixTools::apply_boundary_values (boundary_values,
+                                            mass_matrix, vec, tmp,
+                                            true);
+        constraints.condense(mass_matrix, tmp);
+      }
+
+    // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n
+    // steps may not be sufficient, since roundoff errors may accumulate for
+    // badly conditioned matrices
+    ReductionControl      control(5*tmp.size(), 0., 1e-12, false, false);
     GrowingVectorMemory<> memory;
-    SolverCG<>              cg(control,memory);
+    SolverCG<>            cg(control,memory);
 
     PreconditionSSOR<> prec;
     prec.initialize(mass_matrix, 1.2);
-    // solve
+
     cg.solve (mass_matrix, vec, tmp, prec);
 
-    // distribute solution
     constraints.distribute (vec);
 
-    // copy vec into vec_result. we
-    // can't use ve_result itself
-    // above, since it may be of
-    // another type than Vector<double>
-    // and that wouldn't necessarily go
-    // together with the matrix and
-    // other functions
+    // copy vec into vec_result. we can't use vec_result itself above, since
+    // it may be of another type than Vector<double> and that wouldn't
+    // necessarily go together with the matrix and other functions
     for (unsigned int i=0; i<vec.size(); ++i)
       vec_result(i) = vec(i);
   }
 
 
   template <int dim, class VECTOR, int spacedim>
-  void project (const DoFHandler<dim,spacedim>    &dof,
+  void project (const DoFHandler<dim,spacedim> &dof,
                 const ConstraintMatrix   &constraints,
                 const Quadrature<dim>    &quadrature,
-                const Function<spacedim>      &function,
+                const Function<spacedim> &function,
                 VECTOR                   &vec,
                 const bool                enforce_zero_boundary,
                 const Quadrature<dim-1>  &q_boundary,
@@ -3063,25 +3087,12 @@ namespace VectorTools
             { 0, 1 }
           };
 
-          // The projection is
-          // divided into two steps.
-          // In the first step we
-          // project the boundary
-          // function on the
-          // horizontal shape
-          // functions.  Then the
-          // bounary function is
-          // projected on the
-          // vertical shape
-          // functions.  We begin
-          // with the horizontal
-          // shape functions and
-          // set up a linear system
-          // of equations to get
-          // the values for degrees
-          // of freedom associated
-          // with the interior of
-          // the face.
+          // The projection is divided into two steps.  In the first step we
+          // project the boundary function on the horizontal shape functions.
+          // Then the boundary function is projected on the vertical shape
+          // functions.  We begin with the horizontal shape functions and set
+          // up a linear system of equations to get the values for degrees of
+          // freedom associated with the interior of the face.
           for (unsigned int q_point = 0;
                q_point < fe_values.n_quadrature_points; ++q_point)
             {
@@ -4269,19 +4280,16 @@ namespace VectorTools
 
     std::vector<unsigned int> face_dofs;
 
-    // create FE and mapping
-    // collections for all elements in
-    // use by this DoFHandler
+    // create FE and mapping collections for all elements in use by this
+    // DoFHandler
     hp::FECollection<dim,spacedim>      fe_collection (dof_handler.get_fe());
     hp::MappingCollection<dim,spacedim> mapping_collection;
     for (unsigned int i=0; i<fe_collection.size(); ++i)
       mapping_collection.push_back (mapping);
 
-    // now also create a quadrature
-    // collection for the faces of a
-    // cell. fill it with a quadrature
-    // formula with the support points
-    // on faces for each FE
+    // now also create a quadrature collection for the faces of a cell. fill
+    // it with a quadrature formula with the support points on faces for each
+    // FE
     hp::QCollection<dim-1> face_quadrature_collection;
     for (unsigned int i=0; i<fe_collection.size(); ++i)
       {
@@ -4295,38 +4303,22 @@ namespace VectorTools
         .push_back (Quadrature<dim-1> (unit_support_points));
       }
 
-    // now create the object with which
-    // we will generate the normal
-    // vectors
+    // now create the object with which we will generate the normal vectors
     hp::FEFaceValues<dim,spacedim> x_fe_face_values (mapping_collection,
                                                      fe_collection,
                                                      face_quadrature_collection,
                                                      update_q_points);
 
-    // have a map that stores normal
-    // vectors for each vector-dof
-    // tuple we want to
-    // constrain. since we can get at
-    // the same vector dof tuple more
-    // than once (for example if it is
-    // located at a vertex that we
-    // visit from all adjacent cells),
-    // we will want to average later on
-    // the normal vectors computed on
-    // different cells as described in
-    // the documentation of this
-    // function. however, we can only
-    // average if the contributions
-    // came from different cells,
-    // whereas we want to constrain
-    // twice or more in case the
-    // contributions came from
-    // different faces of the same cell
-    // (i.e. constrain not just the
-    // *average normal direction* but
-    // *all normal directions* we
-    // find). consequently, we also
-    // have to store which cell a
+    // have a map that stores normal vectors for each vector-dof tuple we want
+    // to constrain. since we can get at the same vector dof tuple more than
+    // once (for example if it is located at a vertex that we visit from all
+    // adjacent cells), we will want to average later on the normal vectors
+    // computed on different cells as described in the documentation of this
+    // function. however, we can only average if the contributions came from
+    // different cells, whereas we want to constrain twice or more in case the
+    // contributions came from different faces of the same cell
+    // (i.e. constrain not just the *average normal direction* but *all normal
+    // directions* we find). consequently, we also have to store which cell a
     // normal vector was computed on
     typedef
     std::multimap<internal::VectorDoFTuple<dim>,
@@ -4349,24 +4341,20 @@ namespace VectorTools
               const FiniteElement<dim> &fe = cell->get_fe ();
               typename DH<dim,spacedim>::face_iterator face = cell->face(face_no);
 
-              // get the indices of the
-              // dofs on this cell...
+              // get the indices of the dofs on this cell...
               face_dofs.resize (fe.dofs_per_face);
               face->get_dof_indices (face_dofs, cell->active_fe_index());
 
               x_fe_face_values.reinit (cell, face_no);
               const FEFaceValues<dim> &fe_values = x_fe_face_values.get_present_fe_values();
 
-              // then identify which of
-              // them correspond to the
-              // selected set of vector
-              // components
+              // then identify which of them correspond to the selected set of
+              // vector components
               for (unsigned int i=0; i<face_dofs.size(); ++i)
                 if (fe.face_system_to_component_index(i).first ==
                     first_vector_component)
                   {
-                    // find corresponding other
-                    // components of vector
+                    // find corresponding other components of vector
                     internal::VectorDoFTuple<dim> vector_dofs;
                     vector_dofs.dof_indices[0] = face_dofs[i];
 
@@ -4389,59 +4377,33 @@ namespace VectorTools
                       Assert (vector_dofs.dof_indices[d] < dof_handler.n_dofs(),
                               ExcInternalError());
 
-                    // we need the normal
-                    // vector on this face. we
-                    // know that it is a vector
-                    // of length 1 but at least
-                    // with higher order
-                    // mappings it isn't always
-                    // possible to guarantee
-                    // that each component is
-                    // exact up to zero
-                    // tolerance. in
-                    // particular, as shown in
-                    // the deal.II/no_flux_06
-                    // test, if we just take
-                    // the normal vector as
-                    // given by the fe_values
-                    // object, we can get
-                    // entries in the normal
-                    // vectors of the unit cube
-                    // that have entries up to
+                    // we need the normal vector on this face. we know that it
+                    // is a vector of length 1 but at least with higher order
+                    // mappings it isn't always possible to guarantee that
+                    // each component is exact up to zero tolerance. in
+                    // particular, as shown in the deal.II/no_flux_06 test, if
+                    // we just take the normal vector as given by the
+                    // fe_values object, we can get entries in the normal
+                    // vectors of the unit cube that have entries up to
                     // several times 1e-14.
                     //
-                    // the problem with this is
-                    // that this later yields
-                    // constraints that are
-                    // circular (e.g., in the
-                    // testcase, we get
-                    // constraints of the form
+                    // the problem with this is that this later yields
+                    // constraints that are circular (e.g., in the testcase,
+                    // we get constraints of the form
                     //
                     // x22 =  2.93099e-14*x21 + 2.93099e-14*x23
                     // x21 = -2.93099e-14*x22 + 2.93099e-14*x21
                     //
-                    // in both of these
-                    // constraints, the small
-                    // numbers should be zero
-                    // and the constraints
-                    // should simply be
-                    //    x22 = x21 = 0
+                    // in both of these constraints, the small numbers should
+                    // be zero and the constraints should simply be
+                    // x22 = x21 = 0
                     //
-                    // to achieve this, we
-                    // utilize that we know
-                    // that the normal vector
-                    // has (or should have)
-                    // length 1 and that we can
-                    // simply set small
-                    // elements to zero
-                    // (without having to check
-                    // that they are small
-                    // *relative to something
-                    // else*). we do this and
-                    // then normalize the
-                    // length of the vector
-                    // back to one, just to be
-                    // on the safe side
+                    // to achieve this, we utilize that we know that the
+                    // normal vector has (or should have) length 1 and that we
+                    // can simply set small elements to zero (without having
+                    // to check that they are small *relative to something
+                    // else*). we do this and then normalize the length of the
+                    // vector back to one, just to be on the safe side
                     Point<dim> normal_vector
                       = (cell->face(face_no)->get_boundary()
                          .normal_vector (cell->face(face_no),
@@ -4453,9 +4415,8 @@ namespace VectorTools
                         normal_vector[d] = 0;
                     normal_vector /= normal_vector.norm();
 
-                    // now enter the
-                    // (dofs,(normal_vector,cell))
-                    // entry into the map
+                    // now enter the (dofs,(normal_vector,cell)) entry into
+                    // the map
                     dof_to_normals_map
                     .insert (std::make_pair (vector_dofs,
                                              std::make_pair (normal_vector,
@@ -4463,27 +4424,19 @@ namespace VectorTools
                   }
             }
 
-    // Now do something with the
-    // collected information. To this
-    // end, loop through all sets of
-    // pairs (dofs,normal_vector) and
-    // identify which entries belong to
-    // the same set of dofs and then do
-    // as described in the
-    // documentation, i.e. either
-    // average the normal vector or
-    // don't for this particular set of
-    // dofs
+    // Now do something with the collected information. To this end, loop
+    // through all sets of pairs (dofs,normal_vector) and identify which
+    // entries belong to the same set of dofs and then do as described in the
+    // documentation, i.e. either average the normal vector or don't for this
+    // particular set of dofs
     typename DoFToNormalsMap::const_iterator
     p = dof_to_normals_map.begin();
 
     while (p != dof_to_normals_map.end())
       {
-        // first find the range of entries in
-        // the multimap that corresponds to the
-        // same vector-dof tuple. as usual, we
-        // define the range half-open. the
-        // first entry of course is 'p'
+        // first find the range of entries in the multimap that corresponds to
+        // the same vector-dof tuple. as usual, we define the range
+        // half-open. the first entry of course is 'p'
         typename DoFToNormalsMap::const_iterator same_dof_range[2]
           = { p };
         for (++p; p != dof_to_normals_map.end(); ++p)
@@ -4495,13 +4448,10 @@ namespace VectorTools
         if (p == dof_to_normals_map.end())
           same_dof_range[1] = dof_to_normals_map.end();
 
-        // now compute the reverse mapping: for
-        // each of the cells that contributed
-        // to the current set of vector dofs,
-        // add up the normal vectors. the
-        // values of the map are pairs of
-        // normal vectors and number of cells
-        // that have contributed
+        // now compute the reverse mapping: for each of the cells that
+        // contributed to the current set of vector dofs, add up the normal
+        // vectors. the values of the map are pairs of normal vectors and
+        // number of cells that have contributed
         typedef
         std::map
         <typename DH<dim,spacedim>::active_cell_iterator,
@@ -4525,9 +4475,7 @@ namespace VectorTools
 
               Assert (old_count > 0, ExcInternalError());
 
-              // in the same entry,
-              // store again the now
-              // averaged normal vector
+              // in the same entry, store again the now averaged normal vector
               // and the new count
               cell_to_normals_map[q->second.second]
                 = std::make_pair ((old_normal * old_count + q->second.first) / (old_count + 1),
@@ -4536,8 +4484,7 @@ namespace VectorTools
 
         Assert (cell_to_normals_map.size() >= 1, ExcInternalError());
 
-        // count the maximum number of
-        // contributions from each cell
+        // count the maximum number of contributions from each cell
         unsigned int max_n_contributions_per_cell = 1;
         for (typename CellToNormalsMap::const_iterator
              x = cell_to_normals_map.begin();
@@ -4546,40 +4493,26 @@ namespace VectorTools
             = std::max (max_n_contributions_per_cell,
                         x->second.second);
 
-        // verify that each cell can have only
-        // contributed at most dim times, since
-        // that is the maximum number of faces
-        // that come together at a single place
+        // verify that each cell can have only contributed at most dim times,
+        // since that is the maximum number of faces that come together at a
+        // single place
         Assert (max_n_contributions_per_cell <= dim, ExcInternalError());
 
         switch (max_n_contributions_per_cell)
           {
-            // first deal with the case that a
-            // number of cells all have
-            // registered that they have a
-            // normal vector defined at the
-            // location of a given vector dof,
-            // and that each of them have
-            // encountered this vector dof
-            // exactly once while looping over
-            // all their faces. as stated in
-            // the documentation, this is the
-            // case where we want to simply
-            // average over all normal vectors
+            // first deal with the case that a number of cells all have
+            // registered that they have a normal vector defined at the
+            // location of a given vector dof, and that each of them have
+            // encountered this vector dof exactly once while looping over all
+            // their faces. as stated in the documentation, this is the case
+            // where we want to simply average over all normal vectors
           case 1:
           {
 
-            // compute the average
-            // normal vector from all
-            // the ones that have the
-            // same set of dofs. we
-            // could add them up and
-            // divide them by the
-            // number of additions,
-            // or simply normalize
-            // them right away since
-            // we want them to have
-            // unit length anyway
+            // compute the average normal vector from all the ones that have
+            // the same set of dofs. we could add them up and divide them by
+            // the number of additions, or simply normalize them right away
+            // since we want them to have unit length anyway
             Tensor<1,dim> normal;
             for (typename CellToNormalsMap::const_iterator
                  x = cell_to_normals_map.begin();
@@ -4593,8 +4526,7 @@ namespace VectorTools
                 normal[d] = 0;
             normal /= normal.norm();
 
-            // then construct constraints
-            // from this:
+            // then construct constraints from this:
             const internal::VectorDoFTuple<dim> &
             dof_indices = same_dof_range[0]->first;
             internal::add_constraint (dof_indices, normal,
@@ -4604,56 +4536,27 @@ namespace VectorTools
           }
 
 
-          // this is the slightly
-          // more complicated case
-          // that a single cell has
-          // contributed with exactly
-          // DIM normal vectors to
-          // the same set of vector
-          // dofs. this is what
-          // happens in a corner in
-          // 2d and 3d (but not on an
-          // edge in 3d, where we
-          // have only 2, i.e. <DIM,
-          // contributions. Here we
-          // do not want to average
-          // the normal
-          // vectors. Since we have
-          // DIM contributions, let's
-          // assume (and verify) that
-          // they are in fact all
-          // linearly independent; in
-          // that case, all vector
-          // components are
-          // constrained and we need
-          // to set them to zero
+          // this is the slightly more complicated case that a single cell has
+          // contributed with exactly DIM normal vectors to the same set of
+          // vector dofs. this is what happens in a corner in 2d and 3d (but
+          // not on an edge in 3d, where we have only 2, i.e. <DIM,
+          // contributions. Here we do not want to average the normal
+          // vectors. Since we have DIM contributions, let's assume (and
+          // verify) that they are in fact all linearly independent; in that
+          // case, all vector components are constrained and we need to set
+          // them to zero
           case dim:
           {
-            // assert that indeed
-            // only a single cell has
-            // contributed
+            // assert that indeed only a single cell has contributed
             Assert (cell_to_normals_map.size() == 1,
                     ExcInternalError());
 
-            // check linear
-            // independence by
-            // computing the
-            // determinant of the
-            // matrix created from
-            // all the normal
-            // vectors. if they are
-            // linearly independent,
-            // then the determinant
-            // is nonzero. if they
-            // are orthogonal, then
-            // the matrix is in fact
-            // equal to 1 (since they
-            // are all unit vectors);
-            // make sure the
-            // determinant is larger
-            // than 1e-3 to avoid
-            // cases where cells are
-            // degenerate
+            // check linear independence by computing the determinant of the
+            // matrix created from all the normal vectors. if they are
+            // linearly independent, then the determinant is nonzero. if they
+            // are orthogonal, then the matrix is in fact equal to 1 (since
+            // they are all unit vectors); make sure the determinant is larger
+            // than 1e-3 to avoid cases where cells are degenerate
             {
               Tensor<2,dim> t;
 
@@ -4666,14 +4569,10 @@ namespace VectorTools
                       ExcMessage("Found a set of normal vectors that are nearly collinear."));
             }
 
-            // so all components of
-            // this vector dof are
-            // constrained. enter
-            // this into the
-            // constraint matrix
+            // so all components of this vector dof are constrained. enter
+            // this into the constraint matrix
             //
-            // ignore dofs already
-            // constrained
+            // ignore dofs already constrained
             for (unsigned int i=0; i<dim; ++i)
               if (!constraints.is_constrained (same_dof_range[0]
                                                ->first.dof_indices[i])
@@ -4689,31 +4588,18 @@ namespace VectorTools
           }
 
 
-          // this is the case of an
-          // edge contribution in 3d,
-          // i.e. the vector is
-          // constrained in two
-          // directions but not the
-          // third.
+          // this is the case of an edge contribution in 3d, i.e. the vector
+          // is constrained in two directions but not the third.
           default:
           {
             Assert (dim >= 3, ExcNotImplemented());
             Assert (max_n_contributions_per_cell == 2, ExcInternalError());
 
-            // as described in the
-            // documentation, let us
-            // first collect what
-            // each of the cells
-            // contributed at the
-            // current point. we use
-            // a std::list instead of
-            // a std::set (which
-            // would be more natural)
-            // because std::set
-            // requires that the
-            // stored elements are
-            // comparable with
-            // operator<
+            // as described in the documentation, let us first collect what
+            // each of the cells contributed at the current point. we use a
+            // std::list instead of a std::set (which would be more natural)
+            // because std::set requires that the stored elements are
+            // comparable with operator<
             typedef
             std::map<typename DH<dim,spacedim>::active_cell_iterator, std::list<Tensor<1,dim> > >
             CellContributions;
@@ -4725,24 +4611,13 @@ namespace VectorTools
               cell_contributions[q->second.second].push_back (q->second.first);
             Assert (cell_contributions.size() >= 1, ExcInternalError());
 
-            // now for each cell that
-            // has contributed
-            // determine the number
-            // of normal vectors it
-            // has contributed. we
-            // currently only
-            // implement if this is
-            // dim-1 for all cells
-            // (if a single cell has
-            // contributed dim, or if
-            // all adjacent cells
-            // have contributed 1
-            // normal vector, this is
-            // already handled above)
+            // now for each cell that has contributed determine the number of
+            // normal vectors it has contributed. we currently only implement
+            // if this is dim-1 for all cells (if a single cell has
+            // contributed dim, or if all adjacent cells have contributed 1
+            // normal vector, this is already handled above)
             //
-            // for each contributing
-            // cell compute the
-            // tangential vector that
+            // for each contributing cell compute the tangential vector that
             // remains unconstrained
             std::list<Tensor<1,dim> > tangential_vectors;
             for (typename CellContributions::const_iterator
@@ -4763,81 +4638,24 @@ namespace VectorTools
                   Assert (index == dim-1, ExcInternalError());
                 }
 
-                // calculate the
-                // tangent as the
-                // outer product of
-                // the normal
-                // vectors. since
-                // these vectors do
-                // not need to be
-                // orthogonal (think,
-                // for example, the
-                // case of the
-                // deal.II/no_flux_07
-                // test: a sheared
-                // cube in 3d, with
-                // Q2 elements, where
-                // we have
-                // constraints from
-                // the two normal
-                // vectors of two
-                // faces of the
-                // sheared cube that
-                // are not
-                // perpendicular to
-                // each other), we
-                // have to normalize
-                // the outer product
+                // calculate the tangent as the outer product of the normal
+                // vectors. since these vectors do not need to be orthogonal
+                // (think, for example, the case of the deal.II/no_flux_07
+                // test: a sheared cube in 3d, with Q2 elements, where we have
+                // constraints from the two normal vectors of two faces of the
+                // sheared cube that are not perpendicular to each other), we
+                // have to normalize the outer product
                 Tensor<1,dim> tangent;
                 switch (dim)
                   {
                   case 3:
-                    // take
-                    // cross
-                    // product
-                    // between
-                    // normals[0]
-                    // and
-                    // normals[1]. write
-                    // it in
-                    // the
-                    // current
-                    // form
-                    // (with
-                    // [dim-2])
-                    // to make
-                    // sure
-                    // that
-                    // compilers
-                    // don't
-                    // warn
-                    // about
-                    // out-of-bounds
-                    // accesses
-                    // -- the
-                    // warnings
-                    // are
-                    // bogus
-                    // since we
-                    // get here
-                    // only for
-                    // dim==3,
-                    // but at
-                    // least
-                    // one
-                    // isn't
-                    // quite
-                    // smart
-                    // enough
-                    // to
-                    // notice
-                    // this and
-                    // warns
-                    // when
-                    // compiling
-                    // the
-                    // function
-                    // in 2d
+                    // take cross product between normals[0] and
+                    // normals[1]. write it in the current form (with [dim-2])
+                    // to make sure that compilers don't warn about
+                    // out-of-bounds accesses -- the warnings are bogus since
+                    // we get here only for dim==3, but at least one isn't
+                    // quite smart enough to notice this and warns when
+                    // compiling the function in 2d
                     cross_product (tangent, normals[0], normals[dim-2]);
                     break;
                   default:
@@ -4852,15 +4670,9 @@ namespace VectorTools
                 tangential_vectors.push_back (tangent);
               }
 
-            // go through the list of
-            // tangents and make sure
-            // that they all roughly
-            // point in the same
-            // direction as the first
-            // one (i.e. have an
-            // angle less than 90
-            // degrees); if they
-            // don't then flip their
+            // go through the list of tangents and make sure that they all
+            // roughly point in the same direction as the first one (i.e. have
+            // an angle less than 90 degrees); if they don't then flip their
             // sign
             {
               const Tensor<1,dim> first_tangent = tangential_vectors.front();
@@ -4872,9 +4684,7 @@ namespace VectorTools
                   *t *= -1;
             }
 
-            // now compute the
-            // average tangent and
-            // normalize it
+            // now compute the average tangent and normalize it
             Tensor<1,dim> average_tangent;
             for (typename std::list<Tensor<1,dim> >::const_iterator
                  t = tangential_vectors.begin();
@@ -4883,11 +4693,8 @@ namespace VectorTools
               average_tangent += *t;
             average_tangent /= average_tangent.norm();
 
-            // now all that is left
-            // is that we add the
-            // constraints that the
-            // vector is parallel
-            // to the tangent
+            // now all that is left is that we add the constraints that the
+            // vector is parallel to the tangent
             const internal::VectorDoFTuple<dim> &
             dof_indices = same_dof_range[0]->first;
             internal::add_tangentiality_constraints (dof_indices,
index 776031f0f55f748f3022edb412ce5709a19c8a84..77fa2af43d032a0aeead33dcd217b1db0778eb0d 100644 (file)
@@ -209,6 +209,7 @@ namespace MatrixCreator
         std::vector<unsigned int> dof_indices;
         FullMatrix<double>        cell_matrix;
         dealii::Vector<double>    cell_rhs;
+        const ConstraintMatrix   *constraints;
       };
     }
 
@@ -592,10 +593,15 @@ namespace MatrixCreator
               (data.cell_rhs.size() == dofs_per_cell),
               ExcInternalError());
 
-      matrix->add(data.dof_indices, data.cell_matrix);
       if (right_hand_side != 0)
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
-          (*right_hand_side)(data.dof_indices[i]) += data.cell_rhs(i);
+        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);
     }
   }
 }
@@ -609,7 +615,8 @@ namespace MatrixCreator
                            const DoFHandler<dim,spacedim>    &dof,
                            const Quadrature<dim>    &q,
                            SparseMatrix<number>     &matrix,
-                           const Function<spacedim> *const coefficient)
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
   {
     Assert (matrix.m() == dof.n_dofs(),
             ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
@@ -631,6 +638,7 @@ namespace MatrixCreator
                                   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()),
@@ -648,10 +656,11 @@ namespace MatrixCreator
   void create_mass_matrix (const DoFHandler<dim,spacedim>    &dof,
                            const Quadrature<dim>    &q,
                            SparseMatrix<number>     &matrix,
-                           const Function<spacedim> *const coefficient)
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
   {
     create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof,
-                       q, matrix, coefficient);
+                       q, matrix, coefficient, constraints);
   }
 
 
@@ -663,7 +672,8 @@ namespace MatrixCreator
                            SparseMatrix<number>     &matrix,
                            const Function<spacedim>      &rhs,
                            Vector<double>           &rhs_vector,
-                           const Function<spacedim> *const coefficient)
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
   {
     Assert (matrix.m() == dof.n_dofs(),
             ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
@@ -684,6 +694,7 @@ namespace MatrixCreator
                                   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()),
@@ -703,10 +714,12 @@ namespace MatrixCreator
                            SparseMatrix<number>     &matrix,
                            const Function<spacedim>      &rhs,
                            Vector<double>           &rhs_vector,
-                           const Function<spacedim> *const coefficient)
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
   {
     create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping,
-                       dof, q, matrix, rhs, rhs_vector, coefficient);
+                       dof, q, matrix, rhs, rhs_vector, coefficient,
+                       constraints);
   }
 
 
@@ -716,7 +729,8 @@ namespace MatrixCreator
                            const hp::DoFHandler<dim,spacedim>    &dof,
                            const hp::QCollection<dim>    &q,
                            SparseMatrix<number>     &matrix,
-                           const Function<spacedim> *const coefficient)
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
   {
     Assert (matrix.m() == dof.n_dofs(),
             ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
@@ -734,6 +748,7 @@ namespace MatrixCreator
                                   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()),
@@ -748,24 +763,27 @@ namespace MatrixCreator
 
 
   template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const hp::DoFHandler<dim,spacedim>    &dof,
-                           const hp::QCollection<dim>    &q,
+  void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
+                           const hp::QCollection<dim> &q,
                            SparseMatrix<number>     &matrix,
-                           const Function<spacedim> *const coefficient)
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
   {
-    create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, matrix, coefficient);
+    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,
+  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<double>           &rhs_vector,
-                           const Function<spacedim> *const coefficient)
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
   {
     Assert (matrix.m() == dof.n_dofs(),
             ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
@@ -783,6 +801,7 @@ namespace MatrixCreator
                                   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()),
@@ -797,15 +816,16 @@ namespace MatrixCreator
 
 
   template <int dim, typename number, int spacedim>
-  void create_mass_matrix (const hp::DoFHandler<dim,spacedim>    &dof,
-                           const hp::QCollection<dim>    &q,
+  void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
+                           const hp::QCollection<dim> &q,
                            SparseMatrix<number>     &matrix,
-                           const Function<spacedim>      &rhs,
+                           const Function<spacedim> &rhs,
                            Vector<double>           &rhs_vector,
-                           const Function<spacedim> *const coefficient)
+                           const Function<spacedim> *const coefficient,
+                           const ConstraintMatrix   &constraints)
   {
     create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
-                       matrix, rhs, rhs_vector, coefficient);
+                       matrix, rhs, rhs_vector, coefficient, constraints);
   }
 
 
@@ -1706,7 +1726,8 @@ namespace MatrixCreator
                               const DoFHandler<dim,spacedim>    &dof,
                               const Quadrature<dim>    &q,
                               SparseMatrix<double>     &matrix,
-                              const Function<spacedim> *const coefficient)
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
   {
     Assert (matrix.m() == dof.n_dofs(),
             ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
@@ -1727,6 +1748,7 @@ namespace MatrixCreator
                                   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;
 
     void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &,
                                   SparseMatrix<double> *,
@@ -1749,7 +1771,8 @@ namespace MatrixCreator
   void create_laplace_matrix (const DoFHandler<dim,spacedim>    &dof,
                               const Quadrature<dim>    &q,
                               SparseMatrix<double>     &matrix,
-                              const Function<spacedim> *const coefficient)
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
   {
     create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q, matrix, coefficient);
   }
@@ -1763,7 +1786,8 @@ namespace MatrixCreator
                               SparseMatrix<double>     &matrix,
                               const Function<spacedim>      &rhs,
                               Vector<double>           &rhs_vector,
-                              const Function<spacedim> *const coefficient)
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
   {
     Assert (matrix.m() == dof.n_dofs(),
             ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
@@ -1784,6 +1808,7 @@ namespace MatrixCreator
                                   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;
 
     void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &,
                                   SparseMatrix<double> *,
@@ -1808,7 +1833,8 @@ namespace MatrixCreator
                               SparseMatrix<double>     &matrix,
                               const Function<spacedim>      &rhs,
                               Vector<double>           &rhs_vector,
-                              const Function<spacedim> *const coefficient)
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
   {
     create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q,
                           matrix, rhs, rhs_vector, coefficient);
@@ -1821,7 +1847,8 @@ namespace MatrixCreator
                               const hp::DoFHandler<dim,spacedim>    &dof,
                               const hp::QCollection<dim>    &q,
                               SparseMatrix<double>     &matrix,
-                              const Function<spacedim> *const coefficient)
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
   {
     Assert (matrix.m() == dof.n_dofs(),
             ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
@@ -1839,6 +1866,7 @@ namespace MatrixCreator
                                   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;
 
     void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &,
                                   SparseMatrix<double> *,
@@ -1861,7 +1889,8 @@ namespace MatrixCreator
   void create_laplace_matrix (const hp::DoFHandler<dim,spacedim>    &dof,
                               const hp::QCollection<dim>    &q,
                               SparseMatrix<double>     &matrix,
-                              const Function<spacedim> *const coefficient)
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
   {
     create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, matrix, coefficient);
   }
@@ -1875,7 +1904,8 @@ namespace MatrixCreator
                               SparseMatrix<double>     &matrix,
                               const Function<spacedim>      &rhs,
                               Vector<double>           &rhs_vector,
-                              const Function<spacedim> *const coefficient)
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
   {
     Assert (matrix.m() == dof.n_dofs(),
             ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
@@ -1893,6 +1923,7 @@ namespace MatrixCreator
                                   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;
 
     void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &,
                                   SparseMatrix<double> *,
@@ -1917,7 +1948,8 @@ namespace MatrixCreator
                               SparseMatrix<double>     &matrix,
                               const Function<spacedim>      &rhs,
                               Vector<double>           &rhs_vector,
-                              const Function<spacedim> *const coefficient)
+                              const Function<spacedim> *const coefficient,
+                              const ConstraintMatrix   &constraints)
   {
     create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
                           matrix, rhs, rhs_vector, coefficient);
index ad686e7118259533ad401b615e9cb5d8f0b5acde..07c6301155b9618b9b2d63d744f6acf797400f3b 100644 (file)
@@ -21,13 +21,15 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        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 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 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,
@@ -36,7 +38,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<double>     &matrix,
        const Function<deal_II_space_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient);
+       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,
@@ -44,7 +47,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<double>     &matrix,
        const Function<deal_II_space_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient);
+       const Function<deal_II_space_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
 
 
     template
@@ -53,13 +57,15 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        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 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 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,
@@ -68,7 +74,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<float>     &matrix,
        const Function<deal_II_space_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient);
+       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,
@@ -76,7 +83,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<float>     &matrix,
        const Function<deal_II_space_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_space_dimension> * const coefficient);
+       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>
@@ -121,13 +129,15 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        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 Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
     template
       void MatrixCreator::create_mass_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 Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
     template
       void MatrixCreator::create_mass_matrix<deal_II_dimension>
       (const hp::MappingCollection<deal_II_dimension>       &mapping,
@@ -136,7 +146,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<double>     &matrix,
        const Function<deal_II_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient);
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
     template
       void MatrixCreator::create_mass_matrix<deal_II_dimension>
       (const hp::DoFHandler<deal_II_dimension>    &dof,
@@ -144,7 +155,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<double>     &matrix,
        const Function<deal_II_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient);
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
 
 
     template
@@ -153,13 +165,15 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        const hp::DoFHandler<deal_II_dimension>    &dof,
        const hp::QCollection<deal_II_dimension>    &q,
        SparseMatrix<float>     &matrix,
-       const Function<deal_II_dimension> * const coefficient);
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
     template
       void MatrixCreator::create_mass_matrix<deal_II_dimension>
       (const hp::DoFHandler<deal_II_dimension>    &dof,
        const hp::QCollection<deal_II_dimension>    &q,
        SparseMatrix<float>     &matrix,
-       const Function<deal_II_dimension> * const coefficient);
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
     template
       void MatrixCreator::create_mass_matrix<deal_II_dimension>
       (const hp::MappingCollection<deal_II_dimension>       &mapping,
@@ -168,7 +182,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<float>     &matrix,
        const Function<deal_II_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient);
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
     template
       void MatrixCreator::create_mass_matrix<deal_II_dimension>
       (const hp::DoFHandler<deal_II_dimension>    &dof,
@@ -176,7 +191,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<float>     &matrix,
        const Function<deal_II_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient);
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
 
 
 #if deal_II_dimension != 1
@@ -216,14 +232,16 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
       (const DoFHandler<deal_II_dimension>    &dof,
        const Quadrature<deal_II_dimension>    &q,
        SparseMatrix<double>     &matrix,
-       const Function<deal_II_dimension> * const coefficient);
+       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 Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
     template
       void MatrixCreator::create_laplace_matrix<deal_II_dimension>
       (const Mapping<deal_II_dimension>       &mapping,
@@ -232,7 +250,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<double>     &matrix,
        const Function<deal_II_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient);
+       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,
@@ -240,7 +259,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<double>     &matrix,
        const Function<deal_II_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient);
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
 
 // hp versions of create_laplace_matrix
     template
@@ -248,14 +268,16 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
       (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 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 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,
@@ -264,7 +286,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<double>     &matrix,
        const Function<deal_II_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient);
+       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,
@@ -272,7 +295,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
        SparseMatrix<double>     &matrix,
        const Function<deal_II_dimension>      &rhs,
        Vector<double>           &rhs_vector,
-       const Function<deal_II_dimension> * const coefficient);
+       const Function<deal_II_dimension> * const coefficient,
+       const ConstraintMatrix   &constraints);
 
 
 
diff --git a/tests/deal.II/create_laplace_matrix_constraints_01.cc b/tests/deal.II/create_laplace_matrix_constraints_01.cc
new file mode 100644 (file)
index 0000000..3d1fca2
--- /dev/null
@@ -0,0 +1,139 @@
+//----------------------------  create_laplace_matrix_constraints_01.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_laplace_matrix_constraints_01.cc  ---------------------------
+
+
+// Same as create_laplace_matrix_0[1234] but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+                                  // create a system element composed
+                                  // of one Q1 and one Q2 element
+  FE_Q<dim> element(2);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1,1.),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  Functions::ExpFunction<dim> rhs_function;
+
+  Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+  MatrixTools::
+    create_laplace_matrix (mapping, dof,
+                           quadrature, matrix_ref,
+                           rhs_function, rhs_ref);
+  constraints.condense(matrix_ref, rhs_ref);
+
+  const Function<dim> *const dummy = 0;
+  MatrixTools::create_laplace_matrix (mapping, dof, quadrature, matrix,
+                                      rhs_function, rhs, dummy, constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+  rhs -= rhs_ref;
+  deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_laplace_matrix_constraints_01/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/deal.II/create_laplace_matrix_constraints_01/cmp/generic b/tests/deal.II/create_laplace_matrix_constraints_01/cmp/generic
new file mode 100644 (file)
index 0000000..f29bc07
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:1d::RHS vector error l2: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:2d::RHS vector error l2: 0
+DEAL:3d::Matrix error Frobenius: 0
+DEAL:3d::RHS vector error l2: 0
diff --git a/tests/deal.II/create_laplace_matrix_constraints_02.cc b/tests/deal.II/create_laplace_matrix_constraints_02.cc
new file mode 100644 (file)
index 0000000..a0e2e82
--- /dev/null
@@ -0,0 +1,143 @@
+//----------------------------  create_laplace_matrix_constraints_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_laplace_matrix_constraints_02.cc  ---------------------------
+
+
+// Same as create_laplace_matrix_0[1234] but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+                                  // create a system element composed
+                                  // of one Q1 and one Q2 element
+  FESystem<dim> element(FE_Q<dim>(1), 1,
+                       FE_Q<dim>(2), 1);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1.,2),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    Table<2,DoFTools::Coupling> mask (2, 2);
+    mask(0,0) = mask(1,1) = DoFTools::always;
+    mask(0,1) = mask(1,0) = DoFTools::none;
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, mask, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  Functions::ExpFunction<dim> rhs_function;
+
+  Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+  MatrixTools::
+    create_laplace_matrix (mapping, dof,
+                           quadrature, matrix_ref,
+                           rhs_function, rhs_ref);
+  constraints.condense(matrix_ref, rhs_ref);
+
+  const Function<dim> *const dummy = 0;
+  MatrixTools::create_laplace_matrix (mapping, dof, quadrature, matrix,
+                                      rhs_function, rhs, dummy, constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+  rhs -= rhs_ref;
+  deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_laplace_matrix_constraints_02/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/deal.II/create_laplace_matrix_constraints_02/cmp/generic b/tests/deal.II/create_laplace_matrix_constraints_02/cmp/generic
new file mode 100644 (file)
index 0000000..f29bc07
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:1d::RHS vector error l2: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:2d::RHS vector error l2: 0
+DEAL:3d::Matrix error Frobenius: 0
+DEAL:3d::RHS vector error l2: 0
diff --git a/tests/deal.II/create_laplace_matrix_constraints_02b.cc b/tests/deal.II/create_laplace_matrix_constraints_02b.cc
new file mode 100644 (file)
index 0000000..8e868db
--- /dev/null
@@ -0,0 +1,136 @@
+//----------------------------  create_laplace_matrix_constraints_02b.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_laplace_matrix_constraints_02b.cc  ---------------------------
+
+
+// Same as create_laplace_matrix_0[1234] but passing an additional constraint
+// matrix and comparing results without constraints (and similar to
+// create_laplace_matrix_constraints_02 but without rhs).
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+                                  // create a system element composed
+                                  // of one Q1 and one Q2 element
+  FESystem<dim> element(FE_Q<dim>(1), 1,
+                       FE_Q<dim>(2), 1);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ZeroFunction<dim>(2),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    Table<2,DoFTools::Coupling> mask (2, 2);
+    mask(0,0) = mask(1,1) = DoFTools::always;
+    mask(0,1) = mask(1,0) = DoFTools::none;
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, mask, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  MatrixTools::
+    create_laplace_matrix (mapping, dof,
+                           quadrature, matrix_ref);
+  constraints.condense(matrix_ref);
+
+  const Function<dim> *const dummy = 0;
+  MatrixTools::create_laplace_matrix (mapping, dof, quadrature, matrix,
+                                      dummy, constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_laplace_matrix_constraints_02b/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/deal.II/create_laplace_matrix_constraints_02b/cmp/generic b/tests/deal.II/create_laplace_matrix_constraints_02b/cmp/generic
new file mode 100644 (file)
index 0000000..6a3d892
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:3d::Matrix error Frobenius: 0
diff --git a/tests/deal.II/create_laplace_matrix_constraints_03.cc b/tests/deal.II/create_laplace_matrix_constraints_03.cc
new file mode 100644 (file)
index 0000000..662185c
--- /dev/null
@@ -0,0 +1,140 @@
+//----------------------------  create_laplace_matrix_constraints_03.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_laplace_matrix_constraints_03.cc  ---------------------------
+
+
+// Same as create_laplace_matrix_0[1234] but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+                                  // create a system element composed
+                                  // of one Q1 and one Q2 element
+  FE_Q<dim> element(2);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1,1.),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  Functions::ExpFunction<dim> rhs_function;
+
+  Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+  MatrixTools::
+    create_laplace_matrix (mapping, dof,
+                           quadrature, matrix_ref,
+                           rhs_function, rhs_ref,
+                           &rhs_function);
+  constraints.condense(matrix_ref, rhs_ref);
+
+  MatrixTools::create_laplace_matrix (mapping, dof, quadrature, matrix,
+                                      rhs_function, rhs, &rhs_function,
+                                      constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+  rhs -= rhs_ref;
+  deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_laplace_matrix_constraints_03/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/deal.II/create_laplace_matrix_constraints_03/cmp/generic b/tests/deal.II/create_laplace_matrix_constraints_03/cmp/generic
new file mode 100644 (file)
index 0000000..f29bc07
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:1d::RHS vector error l2: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:2d::RHS vector error l2: 0
+DEAL:3d::Matrix error Frobenius: 0
+DEAL:3d::RHS vector error l2: 0
diff --git a/tests/deal.II/create_laplace_matrix_constraints_04.cc b/tests/deal.II/create_laplace_matrix_constraints_04.cc
new file mode 100644 (file)
index 0000000..5e0237a
--- /dev/null
@@ -0,0 +1,144 @@
+//----------------------------  create_laplace_matrix_constraints_04.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_laplace_matrix_constraints_04.cc  ---------------------------
+
+
+// Same as create_laplace_matrix_0[1234] but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+                                  // create a system element composed
+                                  // of one Q1 and one Q2 element
+  FESystem<dim> element(FE_Q<dim>(1), 1,
+                       FE_Q<dim>(2), 1);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1.,2),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    Table<2,DoFTools::Coupling> mask (2, 2);
+    mask(0,0) = mask(1,1) = DoFTools::always;
+    mask(0,1) = mask(1,0) = DoFTools::none;
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, mask, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  Functions::ExpFunction<dim> rhs_function;
+
+  Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+  MatrixTools::
+    create_laplace_matrix (mapping, dof,
+                           quadrature, matrix_ref,
+                           rhs_function, rhs_ref,
+                           &rhs_function);
+  constraints.condense(matrix_ref, rhs_ref);
+
+  MatrixTools::create_laplace_matrix (mapping, dof, quadrature, matrix,
+                                      rhs_function, rhs, &rhs_function,
+                                      constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+  rhs -= rhs_ref;
+  deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_laplace_matrix_constraints_04/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/deal.II/create_laplace_matrix_constraints_04/cmp/generic b/tests/deal.II/create_laplace_matrix_constraints_04/cmp/generic
new file mode 100644 (file)
index 0000000..f29bc07
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:1d::RHS vector error l2: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:2d::RHS vector error l2: 0
+DEAL:3d::Matrix error Frobenius: 0
+DEAL:3d::RHS vector error l2: 0
diff --git a/tests/deal.II/create_mass_matrix_constraints_01.cc b/tests/deal.II/create_mass_matrix_constraints_01.cc
new file mode 100644 (file)
index 0000000..8e5682e
--- /dev/null
@@ -0,0 +1,139 @@
+//----------------------------  create_mass_matrix_constraints_01.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_mass_matrix_constraints_01.cc  ---------------------------
+
+
+// Same as create_mass_matrix_0[1234] but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+                                  // create a system element composed
+                                  // of one Q1 and one Q2 element
+  FE_Q<dim> element(2);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1,1.),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  Functions::ExpFunction<dim> rhs_function;
+
+  Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+  MatrixTools::
+    create_mass_matrix (mapping, dof,
+                       quadrature, matrix_ref,
+                       rhs_function, rhs_ref);
+  constraints.condense(matrix_ref, rhs_ref);
+
+  const Function<dim> *const dummy = 0;
+  MatrixTools::create_mass_matrix (mapping, dof, quadrature, matrix,
+                                   rhs_function, rhs, dummy, constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+  rhs -= rhs_ref;
+  deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_mass_matrix_constraints_01/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/deal.II/create_mass_matrix_constraints_01/cmp/generic b/tests/deal.II/create_mass_matrix_constraints_01/cmp/generic
new file mode 100644 (file)
index 0000000..f29bc07
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:1d::RHS vector error l2: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:2d::RHS vector error l2: 0
+DEAL:3d::Matrix error Frobenius: 0
+DEAL:3d::RHS vector error l2: 0
diff --git a/tests/deal.II/create_mass_matrix_constraints_02.cc b/tests/deal.II/create_mass_matrix_constraints_02.cc
new file mode 100644 (file)
index 0000000..14971a1
--- /dev/null
@@ -0,0 +1,143 @@
+//----------------------------  create_mass_matrix_constraints_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_mass_matrix_constraints_02.cc  ---------------------------
+
+
+// Same as create_mass_matrix_0[1234] but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+                                  // create a system element composed
+                                  // of one Q1 and one Q2 element
+  FESystem<dim> element(FE_Q<dim>(1), 1,
+                       FE_Q<dim>(2), 1);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1.,2),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    Table<2,DoFTools::Coupling> mask (2, 2);
+    mask(0,0) = mask(1,1) = DoFTools::always;
+    mask(0,1) = mask(1,0) = DoFTools::none;
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, mask, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  Functions::ExpFunction<dim> rhs_function;
+
+  Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+  MatrixTools::
+    create_mass_matrix (mapping, dof,
+                       quadrature, matrix_ref,
+                       rhs_function, rhs_ref);
+  constraints.condense(matrix_ref, rhs_ref);
+
+  const Function<dim> *const dummy = 0;
+  MatrixTools::create_mass_matrix (mapping, dof, quadrature, matrix,
+                                   rhs_function, rhs, dummy, constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+  rhs -= rhs_ref;
+  deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_mass_matrix_constraints_02/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/deal.II/create_mass_matrix_constraints_02/cmp/generic b/tests/deal.II/create_mass_matrix_constraints_02/cmp/generic
new file mode 100644 (file)
index 0000000..f29bc07
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:1d::RHS vector error l2: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:2d::RHS vector error l2: 0
+DEAL:3d::Matrix error Frobenius: 0
+DEAL:3d::RHS vector error l2: 0
diff --git a/tests/deal.II/create_mass_matrix_constraints_02b.cc b/tests/deal.II/create_mass_matrix_constraints_02b.cc
new file mode 100644 (file)
index 0000000..ad74067
--- /dev/null
@@ -0,0 +1,136 @@
+//----------------------------  create_mass_matrix_constraints_02b.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_mass_matrix_constraints_02b.cc  ---------------------------
+
+
+// Same as create_mass_matrix_0[1234]b but passing an additional constraint
+// matrix and comparing results without constraints (and similar to
+// create_mass_matrix_constraints_02 but without rhs).
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+                                  // create a system element composed
+                                  // of one Q1 and one Q2 element
+  FESystem<dim> element(FE_Q<dim>(1), 1,
+                       FE_Q<dim>(2), 1);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ZeroFunction<dim>(2),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    Table<2,DoFTools::Coupling> mask (2, 2);
+    mask(0,0) = mask(1,1) = DoFTools::always;
+    mask(0,1) = mask(1,0) = DoFTools::none;
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, mask, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  MatrixTools::
+    create_mass_matrix (mapping, dof,
+                       quadrature, matrix_ref);
+  constraints.condense(matrix_ref);
+
+  const Function<dim> *const dummy = 0;
+  MatrixTools::create_mass_matrix (mapping, dof, quadrature, matrix,
+                                   dummy, constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_mass_matrix_constraints_02b/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/deal.II/create_mass_matrix_constraints_02b/cmp/generic b/tests/deal.II/create_mass_matrix_constraints_02b/cmp/generic
new file mode 100644 (file)
index 0000000..6a3d892
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:3d::Matrix error Frobenius: 0
diff --git a/tests/deal.II/create_mass_matrix_constraints_03.cc b/tests/deal.II/create_mass_matrix_constraints_03.cc
new file mode 100644 (file)
index 0000000..0c9d3a0
--- /dev/null
@@ -0,0 +1,140 @@
+//----------------------------  create_mass_matrix_constraints_03.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_mass_matrix_constraints_03.cc  ---------------------------
+
+
+// Same as create_mass_matrix_0[1234] but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+                                  // create a system element composed
+                                  // of one Q1 and one Q2 element
+  FE_Q<dim> element(2);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1,1.),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  Functions::ExpFunction<dim> rhs_function;
+
+  Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+  MatrixTools::
+    create_mass_matrix (mapping, dof,
+                       quadrature, matrix_ref,
+                       rhs_function, rhs_ref,
+                        &rhs_function);
+  constraints.condense(matrix_ref, rhs_ref);
+
+  MatrixTools::create_mass_matrix (mapping, dof, quadrature, matrix,
+                                   rhs_function, rhs, &rhs_function,
+                                   constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+  rhs -= rhs_ref;
+  deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_mass_matrix_constraints_03/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/deal.II/create_mass_matrix_constraints_03/cmp/generic b/tests/deal.II/create_mass_matrix_constraints_03/cmp/generic
new file mode 100644 (file)
index 0000000..f29bc07
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:1d::RHS vector error l2: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:2d::RHS vector error l2: 0
+DEAL:3d::Matrix error Frobenius: 0
+DEAL:3d::RHS vector error l2: 0
diff --git a/tests/deal.II/create_mass_matrix_constraints_04.cc b/tests/deal.II/create_mass_matrix_constraints_04.cc
new file mode 100644 (file)
index 0000000..385bda3
--- /dev/null
@@ -0,0 +1,144 @@
+//----------------------------  create_mass_matrix_constraints_04.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_mass_matrix_constraints_04.cc  ---------------------------
+
+
+// Same as create_mass_matrix_0[1234] but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+                                  // create a system element composed
+                                  // of one Q1 and one Q2 element
+  FESystem<dim> element(FE_Q<dim>(1), 1,
+                       FE_Q<dim>(2), 1);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1.,2),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    Table<2,DoFTools::Coupling> mask (2, 2);
+    mask(0,0) = mask(1,1) = DoFTools::always;
+    mask(0,1) = mask(1,0) = DoFTools::none;
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, mask, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  Functions::ExpFunction<dim> rhs_function;
+
+  Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+  MatrixTools::
+    create_mass_matrix (mapping, dof,
+                       quadrature, matrix_ref,
+                       rhs_function, rhs_ref,
+                        &rhs_function);
+  constraints.condense(matrix_ref, rhs_ref);
+
+  MatrixTools::create_mass_matrix (mapping, dof, quadrature, matrix,
+                                   rhs_function, rhs, &rhs_function,
+                                   constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+  rhs -= rhs_ref;
+  deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_mass_matrix_constraints_04/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/deal.II/create_mass_matrix_constraints_04/cmp/generic b/tests/deal.II/create_mass_matrix_constraints_04/cmp/generic
new file mode 100644 (file)
index 0000000..f29bc07
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:1d::RHS vector error l2: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:2d::RHS vector error l2: 0
+DEAL:3d::Matrix error Frobenius: 0
+DEAL:3d::RHS vector error l2: 0
diff --git a/tests/hp/create_laplace_matrix_constraints_02.cc b/tests/hp/create_laplace_matrix_constraints_02.cc
new file mode 100644 (file)
index 0000000..5326b7d
--- /dev/null
@@ -0,0 +1,157 @@
+//----------------------------  create_laplace_matrix_constraints_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_laplace_matrix_constraints_02.cc  ---------------------------
+
+
+// Same as create_laplace_matrix_0[1234] but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+   hp::FECollection<dim> element;
+   element.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
+                                     FE_Q<dim>(2), 1));
+   element.push_back (FESystem<dim> (FE_Q<dim>(2), 1,
+                                     FE_Q<dim>(3), 1));
+   element.push_back (FESystem<dim> (FE_Q<dim>(3), 1,
+                                     FE_Q<dim>(4), 1));
+
+   hp::DoFHandler<dim> dof(tr);
+
+   for (typename hp::DoFHandler<dim>::active_cell_iterator
+        cell = dof.begin_active();
+        cell != dof.end(); ++cell)
+     cell->set_active_fe_index (rand() % element.size());
+
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1.,2),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    Table<2,DoFTools::Coupling> mask (2, 2);
+    mask(0,0) = mask(1,1) = DoFTools::always;
+    mask(0,1) = mask(1,0) = DoFTools::none;
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, mask, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  Functions::ExpFunction<dim> rhs_function;
+
+  Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+  MatrixTools::
+    create_laplace_matrix (hp::MappingCollection<dim>(mapping), dof,
+                           hp::QCollection<dim>(quadrature), matrix_ref,
+                           rhs_function, rhs_ref);
+  constraints.condense(matrix_ref, rhs_ref);
+
+  const Function<dim> *const dummy = 0;
+  MatrixTools::create_laplace_matrix (hp::MappingCollection<dim>(mapping), dof,
+                                      hp::QCollection<dim>(quadrature), matrix,
+                                      rhs_function, rhs, dummy, constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+  rhs -= rhs_ref;
+  deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_laplace_matrix_constraints_02/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/hp/create_laplace_matrix_constraints_02/cmp/generic b/tests/hp/create_laplace_matrix_constraints_02/cmp/generic
new file mode 100644 (file)
index 0000000..f29bc07
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:1d::RHS vector error l2: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:2d::RHS vector error l2: 0
+DEAL:3d::Matrix error Frobenius: 0
+DEAL:3d::RHS vector error l2: 0
diff --git a/tests/hp/create_laplace_matrix_constraints_02b.cc b/tests/hp/create_laplace_matrix_constraints_02b.cc
new file mode 100644 (file)
index 0000000..c348650
--- /dev/null
@@ -0,0 +1,149 @@
+//----------------------------  create_laplace_matrix_constraints_02b.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_laplace_matrix_constraints_02b.cc  ---------------------------
+
+
+// Same as create_laplace_matrix_0[1234]b but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+   hp::FECollection<dim> element;
+   element.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
+                                     FE_Q<dim>(2), 1));
+   element.push_back (FESystem<dim> (FE_Q<dim>(2), 1,
+                                     FE_Q<dim>(3), 1));
+   element.push_back (FESystem<dim> (FE_Q<dim>(3), 1,
+                                     FE_Q<dim>(4), 1));
+
+   hp::DoFHandler<dim> dof(tr);
+
+   for (typename hp::DoFHandler<dim>::active_cell_iterator
+        cell = dof.begin_active();
+        cell != dof.end(); ++cell)
+     cell->set_active_fe_index (rand() % element.size());
+
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1.,2),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    Table<2,DoFTools::Coupling> mask (2, 2);
+    mask(0,0) = mask(1,1) = DoFTools::always;
+    mask(0,1) = mask(1,0) = DoFTools::none;
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, mask, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  MatrixTools::
+    create_laplace_matrix (hp::MappingCollection<dim>(mapping), dof,
+                           hp::QCollection<dim>(quadrature), matrix_ref);
+  constraints.condense(matrix_ref);
+
+  const Function<dim> *const dummy = 0;
+  MatrixTools::create_laplace_matrix (hp::MappingCollection<dim>(mapping), dof,
+                                      hp::QCollection<dim>(quadrature), matrix,
+                                      dummy, constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_laplace_matrix_constraints_02b/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/hp/create_laplace_matrix_constraints_02b/cmp/generic b/tests/hp/create_laplace_matrix_constraints_02b/cmp/generic
new file mode 100644 (file)
index 0000000..6a3d892
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:3d::Matrix error Frobenius: 0
diff --git a/tests/hp/create_mass_matrix_constraints_02.cc b/tests/hp/create_mass_matrix_constraints_02.cc
new file mode 100644 (file)
index 0000000..8b85801
--- /dev/null
@@ -0,0 +1,157 @@
+//----------------------------  create_mass_matrix_constraints_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_mass_matrix_constraints_02.cc  ---------------------------
+
+
+// Same as create_mass_matrix_0[1234] but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+   hp::FECollection<dim> element;
+   element.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
+                                     FE_Q<dim>(2), 1));
+   element.push_back (FESystem<dim> (FE_Q<dim>(2), 1,
+                                     FE_Q<dim>(3), 1));
+   element.push_back (FESystem<dim> (FE_Q<dim>(3), 1,
+                                     FE_Q<dim>(4), 1));
+
+   hp::DoFHandler<dim> dof(tr);
+
+   for (typename hp::DoFHandler<dim>::active_cell_iterator
+        cell = dof.begin_active();
+        cell != dof.end(); ++cell)
+     cell->set_active_fe_index (rand() % element.size());
+
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1.,2),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    Table<2,DoFTools::Coupling> mask (2, 2);
+    mask(0,0) = mask(1,1) = DoFTools::always;
+    mask(0,1) = mask(1,0) = DoFTools::none;
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, mask, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  Functions::ExpFunction<dim> rhs_function;
+
+  Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+  MatrixTools::
+    create_mass_matrix (hp::MappingCollection<dim>(mapping), dof,
+                       hp::QCollection<dim>(quadrature), matrix_ref,
+                       rhs_function, rhs_ref);
+  constraints.condense(matrix_ref, rhs_ref);
+
+  const Function<dim> *const dummy = 0;
+  MatrixTools::create_mass_matrix (hp::MappingCollection<dim>(mapping), dof,
+                                   hp::QCollection<dim>(quadrature), matrix,
+                                   rhs_function, rhs, dummy, constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+  rhs -= rhs_ref;
+  deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_mass_matrix_constraints_02/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/hp/create_mass_matrix_constraints_02/cmp/generic b/tests/hp/create_mass_matrix_constraints_02/cmp/generic
new file mode 100644 (file)
index 0000000..f29bc07
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:1d::RHS vector error l2: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:2d::RHS vector error l2: 0
+DEAL:3d::Matrix error Frobenius: 0
+DEAL:3d::RHS vector error l2: 0
diff --git a/tests/hp/create_mass_matrix_constraints_02b.cc b/tests/hp/create_mass_matrix_constraints_02b.cc
new file mode 100644 (file)
index 0000000..27c8597
--- /dev/null
@@ -0,0 +1,149 @@
+//----------------------------  create_mass_matrix_constraints_02b.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  create_mass_matrix_constraints_02b.cc  ---------------------------
+
+
+// Same as create_mass_matrix_0[1234]b but passing an additional constraint
+// matrix and comparing results without constraints
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+   hp::FECollection<dim> element;
+   element.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
+                                     FE_Q<dim>(2), 1));
+   element.push_back (FESystem<dim> (FE_Q<dim>(2), 1,
+                                     FE_Q<dim>(3), 1));
+   element.push_back (FESystem<dim> (FE_Q<dim>(3), 1,
+                                     FE_Q<dim>(4), 1));
+
+   hp::DoFHandler<dim> dof(tr);
+
+   for (typename hp::DoFHandler<dim>::active_cell_iterator
+        cell = dof.begin_active();
+        cell != dof.end(); ++cell)
+     cell->set_active_fe_index (rand() % element.size());
+
+  dof.distribute_dofs(element);
+
+                                  // use a more complicated mapping
+                                  // of the domain and a quadrature
+                                  // formula suited to the elements
+                                  // we have here
+  MappingQ<dim> mapping (3);
+  QGauss<dim> quadrature(6);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  VectorTools::interpolate_boundary_values(dof, 0,
+                                           ConstantFunction<dim>(1.,2),
+                                           constraints);
+  constraints.close ();
+
+                                  // create sparsity pattern. note
+                                  // that different components should
+                                  // not couple, so use pattern
+  SparsityPattern sparsity;
+  {
+    Table<2,DoFTools::Coupling> mask (2, 2);
+    mask(0,0) = mask(1,1) = DoFTools::always;
+    mask(0,1) = mask(1,0) = DoFTools::none;
+    CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, mask, csp, constraints);
+    sparsity.copy_from(csp);
+  }
+
+  SparseMatrix<double> matrix, matrix_ref;
+  matrix.reinit (sparsity);
+  matrix_ref.reinit (sparsity);
+
+  MatrixTools::
+    create_mass_matrix (hp::MappingCollection<dim>(mapping), dof,
+                       hp::QCollection<dim>(quadrature), matrix_ref);
+  constraints.condense(matrix_ref);
+
+  const Function<dim> *const dummy = 0;
+  MatrixTools::create_mass_matrix (hp::MappingCollection<dim>(mapping), dof,
+                                   hp::QCollection<dim>(quadrature), matrix,
+                                   dummy, constraints);
+
+  // compute reference: need to cancel constrained entries as these will in
+  // general get different values
+  matrix.add(-1., matrix_ref);
+  for (unsigned int i=0; i<matrix.m(); ++i)
+    if (constraints.is_constrained(i)==true)
+      matrix.diag_element(i) = 0;
+  deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile ("create_mass_matrix_constraints_02b/output");
+  deallog << std::setprecision (3);
+  deallog.attach(logfile);
+  deallog.threshold_double(1e-10);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
diff --git a/tests/hp/create_mass_matrix_constraints_02b/cmp/generic b/tests/hp/create_mass_matrix_constraints_02b/cmp/generic
new file mode 100644 (file)
index 0000000..6a3d892
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:3d::Matrix error Frobenius: 0

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.