From: kronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 24 Feb 2009 16:33:53 +0000 (+0000)
Subject: Add a test for inhomogeneous constraints. This tests the symmetric case. Still need... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=408d62ac901bb950e3c930b31e4e539c2bb9f313;p=dealii-svn.git

Add a test for inhomogeneous constraints. This tests the symmetric case. Still need to write a test for non-symmetric matrices.

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

diff --git a/tests/deal.II/inhomogeneous_constraints.cc b/tests/deal.II/inhomogeneous_constraints.cc
new file mode 100644
index 0000000000..ea796cca55
--- /dev/null
+++ b/tests/deal.II/inhomogeneous_constraints.cc
@@ -0,0 +1,801 @@
+//------------------  inhomogeneous_constraints.cc  ------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2009 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.
+//
+//------------------  inhomogeneous_constraints.cc  ------------------------
+
+
+// this function tests the correctness of the implementation of
+// inhomogeneous constraints. The program is a modification of the step-27
+// tutorial program with hp elements and the constraints arising in that
+// situation. the idea of the test is to set up a matrix with standard tools
+// (i.e., constraints and the boundary value list), and compare that with
+// the new function.
+
+#include "../tests.h"
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <base/utilities.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_refinement.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+#include <fe/fe_q.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <numerics/error_estimator.h>
+#include <lac/compressed_simple_sparsity_pattern.h>
+#include <hp/dof_handler.h>
+#include <hp/fe_values.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+std::ofstream logfile("inhomogeneous_constraints/output");
+
+using namespace dealii;
+
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+    ~LaplaceProblem ();
+
+    void run ();
+    
+  private:
+    void setup_system ();
+    void test_equality ();
+    void assemble_reference ();
+    void assemble_test_1 ();
+    void assemble_test_2 ();
+    void solve ();
+    void create_coarse_grid ();
+    void estimate_smoothness (Vector<float> &smoothness_indicators) const;
+    void postprocess (const unsigned int cycle);
+
+    Triangulation<dim>   triangulation;
+
+    hp::DoFHandler<dim>      dof_handler;
+    hp::FECollection<dim>    fe_collection;
+    hp::QCollection<dim>     quadrature_collection;
+    hp::QCollection<dim-1>   face_quadrature_collection;
+
+    ConstraintMatrix     hanging_nodes_only;
+    ConstraintMatrix     test_all_constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> reference_matrix;
+    SparseMatrix<double> test_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       reference_rhs;
+    Vector<double>       test_rhs;
+
+    const unsigned int max_degree;
+};
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+  public:
+    RightHandSide () : Function<dim> () {}
+    
+    virtual double value (const Point<dim>   &p,
+			  const unsigned int  component) const;
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>   &p,
+			   const unsigned int  /*component*/) const
+{
+  double product = 1;
+  for (unsigned int d=0; d<dim; ++d)
+    product *= (p[d]+1);
+  return product;
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem ()
+		:
+		dof_handler (triangulation),
+		max_degree (5)
+{
+  for (unsigned int degree=2; degree<=max_degree; ++degree)
+    {
+      fe_collection.push_back (FE_Q<dim>(degree));
+      quadrature_collection.push_back (QGauss<dim>(degree+1));
+      face_quadrature_collection.push_back (QGauss<dim-1>(degree+1));
+    }
+}
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem () 
+{
+  dof_handler.clear ();
+}
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe_collection);
+
+  solution.reinit (dof_handler.n_dofs());
+  reference_rhs.reinit (dof_handler.n_dofs());
+  test_rhs.reinit (dof_handler.n_dofs());
+
+  hanging_nodes_only.clear ();
+  test_all_constraints.clear ();
+
+				   // add boundary conditions as
+				   // inhomogeneous constraints here. In
+				   // contrast to step-27, we choose a
+				   // constant function with value 1 here.
+  {
+    std::map<unsigned int,double> boundary_values;
+    VectorTools::interpolate_boundary_values (dof_handler,
+					      0,
+					      ConstantFunction<dim>(1.),
+					      boundary_values);
+    std::map<unsigned int,double>::const_iterator boundary_value = boundary_values.begin();
+    for ( ; boundary_value !=boundary_values.end(); ++boundary_value)
+      {
+	test_all_constraints.add_line(boundary_value->first);
+	test_all_constraints.set_inhomogeneity (boundary_value->first, 
+						boundary_value->second);
+      }
+  }
+  DoFTools::make_hanging_node_constraints (dof_handler,
+					   hanging_nodes_only);
+  DoFTools::make_hanging_node_constraints (dof_handler,
+					   test_all_constraints);
+  hanging_nodes_only.close ();
+  test_all_constraints.close ();
+
+  CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
+				       dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, csp,
+				   hanging_nodes_only, true);
+  sparsity_pattern.copy_from (csp);
+
+  reference_matrix.reinit (sparsity_pattern);
+  test_matrix.reinit (sparsity_pattern);
+}
+
+
+
+				   // test whether we are equal with the
+				   // standard matrix and right hand side
+template <int dim>
+void LaplaceProblem<dim>::test_equality ()
+{
+				   // need to manually go through the
+				   // matrix, since we can have different
+				   // entries in constrained lines.
+  for (unsigned int i=0; i<reference_matrix.m(); ++i)
+    {
+      SparseMatrix<double>::const_iterator reference = reference_matrix.begin(i);
+      SparseMatrix<double>::iterator test = test_matrix.begin(i);
+      if (test_all_constraints.is_constrained(i) == false)
+	{
+	  for ( ; test != test_matrix.end(i); ++test, ++reference)
+	      test->value() -= reference->value();
+	}
+      else
+	for ( ; test != test_matrix.end(i); ++test)
+	  test->value() = 0;
+    }
+
+  deallog << "Matrix difference norm: " 
+	  << test_matrix.frobenius_norm() << std::endl;
+  Assert (test_matrix.frobenius_norm() < 1e-13, ExcInternalError());
+
+				   // same here -- Dirichlet lines will have
+				   // nonzero rhs, whereas we will have
+				   // nonzero one.
+  for (unsigned int i=0; i<reference_matrix.m(); ++i)
+    if (test_all_constraints.is_constrained(i) == false)
+      test_rhs(i) -= reference_rhs(i);
+    else
+      test_rhs(i) = 0;
+
+  deallog << "rhs difference norm: " 
+	  << test_rhs.l2_norm() << std::endl;
+
+  Assert (test_rhs.l2_norm() < 1e-14, ExcInternalError());
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_reference () 
+{
+  reference_matrix = 0;
+  reference_rhs = 0;
+
+  hp::FEValues<dim> hp_fe_values (fe_collection,
+				  quadrature_collection, 
+				  update_values    |  update_gradients |
+				  update_quadrature_points  |  update_JxW_values);
+
+  const RightHandSide<dim> rhs_function;
+
+  FullMatrix<double>   cell_matrix;
+  Vector<double>       cell_rhs;
+
+  std::vector<unsigned int> local_dof_indices;
+  
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+      cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+      cell_matrix = 0;
+
+      cell_rhs.reinit (dofs_per_cell);
+      cell_rhs = 0;
+
+      hp_fe_values.reinit (cell);
+      
+      const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
+
+      std::vector<double>  rhs_values (fe_values.n_quadrature_points);
+      rhs_function.value_list (fe_values.get_quadrature_points(),
+			       rhs_values);
+      
+      for (unsigned int q_point=0;
+	   q_point<fe_values.n_quadrature_points;
+	   ++q_point)
+	for (unsigned int i=0; i<dofs_per_cell; ++i)
+	  {
+	    for (unsigned int j=0; j<dofs_per_cell; ++j)
+	      cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+				   fe_values.shape_grad(j,q_point) *
+				   fe_values.JxW(q_point));
+
+	    cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+			    rhs_values[q_point] *
+			    fe_values.JxW(q_point));
+	  }
+
+      local_dof_indices.resize (dofs_per_cell);
+      cell->get_dof_indices (local_dof_indices);
+
+      reference_matrix.add(local_dof_indices, cell_matrix);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	reference_rhs(local_dof_indices[i]) += cell_rhs(i);
+    }
+
+  hanging_nodes_only.condense (reference_matrix, reference_rhs);
+  deallog << "Reference matrix nonzeros: " << reference_matrix.n_nonzero_elements() 
+	  << ", actually: " << reference_matrix.n_actually_nonzero_elements () 
+	  << std::endl;
+
+  std::map<unsigned int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+					    0,
+					    ConstantFunction<dim>(1.),
+					    boundary_values);
+  MatrixTools::apply_boundary_values (boundary_values,
+				      reference_matrix,
+				      solution,
+				      reference_rhs);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test_1 () 
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  hp::FEValues<dim> hp_fe_values (fe_collection,
+				  quadrature_collection, 
+				  update_values    |  update_gradients |
+				  update_quadrature_points  |  update_JxW_values);
+
+  const RightHandSide<dim> rhs_function;
+  
+  FullMatrix<double>   cell_matrix;
+  Vector<double>       cell_rhs;
+
+  std::vector<unsigned int> local_dof_indices;
+  
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+      cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+      cell_matrix = 0;
+
+      cell_rhs.reinit (dofs_per_cell);
+      cell_rhs = 0;
+
+      hp_fe_values.reinit (cell);
+      
+      const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
+
+      std::vector<double>  rhs_values (fe_values.n_quadrature_points);
+      rhs_function.value_list (fe_values.get_quadrature_points(),
+			       rhs_values);
+      
+      for (unsigned int q_point=0;
+	   q_point<fe_values.n_quadrature_points;
+	   ++q_point)
+	for (unsigned int i=0; i<dofs_per_cell; ++i)
+	  {
+	    for (unsigned int j=0; j<dofs_per_cell; ++j)
+	      cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+				   fe_values.shape_grad(j,q_point) *
+				   fe_values.JxW(q_point));
+
+	    cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+			    rhs_values[q_point] *
+			    fe_values.JxW(q_point));
+	  }
+
+      local_dof_indices.resize (dofs_per_cell);
+      cell->get_dof_indices (local_dof_indices);
+
+      test_matrix.add(local_dof_indices, cell_matrix);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	test_rhs(local_dof_indices[i]) += cell_rhs(i);
+	
+    }
+
+  test_all_constraints.condense (test_matrix, test_rhs);
+  deallog << "Test matrix 1 nonzeros: " << test_matrix.n_nonzero_elements() 
+	  << ", actually: " << test_matrix.n_actually_nonzero_elements () 
+	  << std::endl;
+
+  test_equality();
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_test_2 () 
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  hp::FEValues<dim> hp_fe_values (fe_collection,
+				  quadrature_collection, 
+				  update_values    |  update_gradients |
+				  update_quadrature_points  |  update_JxW_values);
+
+  const RightHandSide<dim> rhs_function;
+  
+  FullMatrix<double>   cell_matrix;
+  Vector<double>       cell_rhs;
+
+  std::vector<unsigned int> local_dof_indices;
+  
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+      cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+      cell_matrix = 0;
+
+      cell_rhs.reinit (dofs_per_cell);
+      cell_rhs = 0;
+
+      hp_fe_values.reinit (cell);
+      
+      const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
+
+      std::vector<double>  rhs_values (fe_values.n_quadrature_points);
+      rhs_function.value_list (fe_values.get_quadrature_points(),
+			       rhs_values);
+      
+      for (unsigned int q_point=0;
+	   q_point<fe_values.n_quadrature_points;
+	   ++q_point)
+	for (unsigned int i=0; i<dofs_per_cell; ++i)
+	  {
+	    for (unsigned int j=0; j<dofs_per_cell; ++j)
+	      cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
+				   fe_values.shape_grad(j,q_point) *
+				   fe_values.JxW(q_point));
+
+	    cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+			    rhs_values[q_point] *
+			    fe_values.JxW(q_point));
+	  }
+
+      local_dof_indices.resize (dofs_per_cell);
+      cell->get_dof_indices (local_dof_indices);
+
+      test_all_constraints.distribute_local_to_global (cell_matrix,
+						       cell_rhs,
+						       local_dof_indices,
+						       test_matrix,
+						       test_rhs);
+    }
+  deallog << "Test matrix 2 nonzeros: " << test_matrix.n_nonzero_elements() 
+	  << ", actually: " << test_matrix.n_actually_nonzero_elements () 
+	  << std::endl;
+  test_equality();
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve () 
+{
+  SolverControl           solver_control (reference_rhs.size(),
+					  1e-8*reference_rhs.l2_norm());
+  SolverCG<>              cg (solver_control);
+
+  PreconditionSSOR<> preconditioner;
+  preconditioner.initialize(reference_matrix, 1.2);
+
+  cg.solve (reference_matrix, solution, reference_rhs,
+	    preconditioner);
+
+  // test also distribute function
+  Vector<double> solution_test (solution);
+
+  hanging_nodes_only.distribute (solution);
+
+  // test also distribute function
+  test_all_constraints.distribute(solution_test);
+  solution_test -= solution;
+  deallog << "Distribute error: " << solution_test.l2_norm () << std::endl;
+  Assert (solution_test.l2_norm() < 1e-8, ExcInternalError());
+}
+
+
+template <int dim>
+void LaplaceProblem<dim>::postprocess (const unsigned int cycle)
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+				      face_quadrature_collection,
+				      typename FunctionMap<dim>::type(),
+				      solution,
+				      estimated_error_per_cell);
+
+  Vector<float> smoothness_indicators (triangulation.n_active_cells());
+  estimate_smoothness (smoothness_indicators);
+
+
+  {
+    GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+						     estimated_error_per_cell,
+						     0.3, 0.03);
+
+    float max_smoothness = *std::min_element (smoothness_indicators.begin(),
+					      smoothness_indicators.end()),
+	  min_smoothness = *std::max_element (smoothness_indicators.begin(),
+					      smoothness_indicators.end());
+    {
+      typename hp::DoFHandler<dim>::active_cell_iterator
+	cell = dof_handler.begin_active(),
+	endc = dof_handler.end();
+      for (unsigned int index=0; cell!=endc; ++cell, ++index)
+	if (cell->refine_flag_set())
+	  {
+	    max_smoothness = std::max (max_smoothness,
+				       smoothness_indicators(index));
+	    min_smoothness = std::min (min_smoothness,
+				       smoothness_indicators(index));
+	  }
+    }
+    const float threshold_smoothness = (max_smoothness + min_smoothness) / 2;
+
+    {
+      typename hp::DoFHandler<dim>::active_cell_iterator
+	cell = dof_handler.begin_active(),
+	endc = dof_handler.end();
+      for (unsigned int index=0; cell!=endc; ++cell, ++index)
+	if (cell->refine_flag_set()
+	    &&
+	    (smoothness_indicators(index) > threshold_smoothness)
+	    &&
+	    (cell->active_fe_index()+1 < fe_collection.size()))
+	  {
+	    cell->clear_refine_flag();
+	    cell->set_active_fe_index (cell->active_fe_index() + 1);
+	  }
+    } 
+
+    triangulation.execute_coarsening_and_refinement ();
+  }
+}
+
+
+template <>
+void LaplaceProblem<2>::create_coarse_grid ()
+{
+  const unsigned int dim = 2;
+  
+  static const Point<2> vertices_1[]
+    = {  Point<2> (-1.,   -1.),
+         Point<2> (-1./2, -1.),
+         Point<2> (0.,    -1.),
+         Point<2> (+1./2, -1.),
+         Point<2> (+1,    -1.),
+	     
+         Point<2> (-1.,   -1./2.),
+         Point<2> (-1./2, -1./2.),
+         Point<2> (0.,    -1./2.),
+         Point<2> (+1./2, -1./2.),
+         Point<2> (+1,    -1./2.),
+	     
+         Point<2> (-1.,   0.),
+         Point<2> (-1./2, 0.),
+         Point<2> (+1./2, 0.),
+         Point<2> (+1,    0.),
+	     
+         Point<2> (-1.,   1./2.),
+         Point<2> (-1./2, 1./2.),
+         Point<2> (0.,    1./2.),
+         Point<2> (+1./2, 1./2.),
+         Point<2> (+1,    1./2.),
+	     
+         Point<2> (-1.,   1.),
+         Point<2> (-1./2, 1.),
+         Point<2> (0.,    1.),			  
+         Point<2> (+1./2, 1.),
+         Point<2> (+1,    1.)    };
+  const unsigned int
+    n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]);
+  const std::vector<Point<dim> > vertices (&vertices_1[0],
+                                           &vertices_1[n_vertices]);
+  static const int cell_vertices[][GeometryInfo<dim>::vertices_per_cell]
+    = {{0, 1, 5, 6},
+       {1, 2, 6, 7},
+       {2, 3, 7, 8},
+       {3, 4, 8, 9},
+       {5, 6, 10, 11},
+       {8, 9, 12, 13},
+       {10, 11, 14, 15},
+       {12, 13, 17, 18},
+       {14, 15, 19, 20},
+       {15, 16, 20, 21},
+       {16, 17, 21, 22},
+       {17, 18, 22, 23}};
+  const unsigned int
+    n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]);
+
+  std::vector<CellData<dim> > cells (n_cells, CellData<dim>());
+  for (unsigned int i=0; i<n_cells; ++i) 
+    {
+      for (unsigned int j=0;
+           j<GeometryInfo<dim>::vertices_per_cell;
+           ++j)
+        cells[i].vertices[j] = cell_vertices[i][j];
+      cells[i].material_id = 0;
+    }
+
+  triangulation.create_triangulation (vertices,
+                                    cells,
+                                    SubCellData());
+  triangulation.refine_global (3);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  for (unsigned int cycle=0; cycle<3; ++cycle)
+    {
+      if (cycle == 0)
+	create_coarse_grid ();
+
+      setup_system ();
+
+      deallog << std::endl << std::endl
+	      << "   Number of active cells:       "
+	      << triangulation.n_active_cells()
+	      << std::endl
+	      << "   Number of degrees of freedom: "
+	      << dof_handler.n_dofs()
+	      << std::endl
+	      << "   Number of constraints       : "
+	      << hanging_nodes_only.n_constraints()
+	      << std::endl;
+
+      assemble_reference ();
+      assemble_test_1 ();
+      assemble_test_2 ();
+
+      solve ();
+      postprocess (cycle);
+    }
+}
+
+
+				   // this function is copied verbatim from step-27
+template <int dim>
+void
+LaplaceProblem<dim>::
+estimate_smoothness (Vector<float> &smoothness_indicators) const
+{
+  const unsigned int N = max_degree;
+
+  std::vector<Tensor<1,dim> > k_vectors;
+  std::vector<unsigned int>   k_vectors_magnitude;
+  switch (dim)
+    {
+      case 2:
+      {
+	for (unsigned int i=0; i<N; ++i)
+	  for (unsigned int j=0; j<N; ++j)
+	    if (!((i==0) && (j==0))
+		&&
+		(i*i + j*j < N*N))
+	      {
+		k_vectors.push_back (Point<dim>(numbers::PI * i,
+						numbers::PI * j));
+		k_vectors_magnitude.push_back (i*i+j*j);
+	      }
+	
+	break;
+      }
+
+      case 3:
+      {
+	for (unsigned int i=0; i<N; ++i)
+	  for (unsigned int j=0; j<N; ++j)
+	    for (unsigned int k=0; k<N; ++k)
+	      if (!((i==0) && (j==0) && (k==0))
+		  &&
+		  (i*i + j*j + k*k < N*N))
+		{
+		  k_vectors.push_back (Point<dim>(numbers::PI * i,
+						  numbers::PI * j,
+						  numbers::PI * k));
+		  k_vectors_magnitude.push_back (i*i+j*j+k*k);
+	      }
+	
+	break;
+      }
+      
+      default:
+	    Assert (false, ExcNotImplemented());
+    }
+
+  const unsigned n_fourier_modes = k_vectors.size();
+  std::vector<double> ln_k (n_fourier_modes);
+  for (unsigned int i=0; i<n_fourier_modes; ++i)
+    ln_k[i] = std::log (k_vectors[i].norm());
+  
+  std::vector<Table<2,std::complex<double> > >
+    fourier_transform_matrices (fe_collection.size());
+  QGauss<1>      base_quadrature (2);
+  QIterated<dim> quadrature (base_quadrature, N);
+
+  for (unsigned int fe=0; fe<fe_collection.size(); ++fe)
+    {
+      fourier_transform_matrices[fe].reinit (n_fourier_modes,
+					     fe_collection[fe].dofs_per_cell);
+
+      for (unsigned int k=0; k<n_fourier_modes; ++k)
+	for (unsigned int j=0; j<fe_collection[fe].dofs_per_cell; ++j)
+	  {
+	    std::complex<double> sum = 0;
+	    for (unsigned int q=0; q<quadrature.size(); ++q)
+	      {
+		const Point<dim> x_q = quadrature.point(q);
+		sum += std::exp(std::complex<double>(0,1) *
+				(k_vectors[k] * x_q)) *
+		       fe_collection[fe].shape_value(j,x_q) *
+		       quadrature.weight(q);
+	      }
+	    fourier_transform_matrices[fe](k,j)
+	      = sum / std::pow(2*numbers::PI, 1.*dim/2);
+	  }
+    }
+
+  std::vector<std::complex<double> > fourier_coefficients (n_fourier_modes);
+  Vector<double>                     local_dof_values;
+
+  typename hp::DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (unsigned int index=0; cell!=endc; ++cell, ++index)
+    {
+      local_dof_values.reinit (cell->get_fe().dofs_per_cell);
+      cell->get_dof_values (solution, local_dof_values);
+
+      for (unsigned int f=0; f<n_fourier_modes; ++f)
+	{
+	  fourier_coefficients[f] = 0;
+	  
+	  for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+	    fourier_coefficients[f] += 
+	      fourier_transform_matrices[cell->active_fe_index()](f,i)
+	      *
+	      local_dof_values(i);
+	}
+
+      std::map<unsigned int, double> k_to_max_U_map;
+      for (unsigned int f=0; f<n_fourier_modes; ++f)
+	if ((k_to_max_U_map.find (k_vectors_magnitude[f]) ==
+	     k_to_max_U_map.end())
+	    ||
+	    (k_to_max_U_map[k_vectors_magnitude[f]] <
+	     std::abs (fourier_coefficients[f])))
+	  k_to_max_U_map[k_vectors_magnitude[f]]
+	    = std::abs (fourier_coefficients[f]);
+      double  sum_1           = 0,
+	      sum_ln_k        = 0,
+	      sum_ln_k_square = 0,
+	      sum_ln_U        = 0,
+	      sum_ln_U_ln_k   = 0;
+      for (unsigned int f=0; f<n_fourier_modes; ++f)
+	if (k_to_max_U_map[k_vectors_magnitude[f]] ==
+	    std::abs (fourier_coefficients[f]))
+	  {
+	    sum_1 += 1;
+	    sum_ln_k += ln_k[f];
+	    sum_ln_k_square += ln_k[f]*ln_k[f];
+	    sum_ln_U += std::log (std::abs (fourier_coefficients[f]));
+	    sum_ln_U_ln_k += std::log (std::abs (fourier_coefficients[f])) *
+			     ln_k[f];
+	  }
+
+      const double mu
+	= (1./(sum_1*sum_ln_k_square - sum_ln_k*sum_ln_k)
+	   *
+	   (sum_ln_k*sum_ln_U - sum_1*sum_ln_U_ln_k));
+
+      smoothness_indicators(index) = mu - 1.*dim/2;
+    }
+}
+
+
+int main () 
+{
+  deallog << std::setprecision (2);
+  logfile << std::setprecision (2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-12);
+
+  LaplaceProblem<2> laplace_problem;
+  laplace_problem.run ();
+}
diff --git a/tests/deal.II/inhomogeneous_constraints/cmp/generic b/tests/deal.II/inhomogeneous_constraints/cmp/generic
new file mode 100644
index 0000000000..91ffd93473
--- /dev/null
+++ b/tests/deal.II/inhomogeneous_constraints/cmp/generic
@@ -0,0 +1,46 @@
+
+DEAL::
+DEAL::
+DEAL::   Number of active cells:       768
+DEAL::   Number of degrees of freedom: 3264
+DEAL::   Number of constraints       : 0
+DEAL::Reference matrix nonzeros: 49920, actually: 49766
+DEAL::Test matrix 1 nonzeros: 49920, actually: 42514
+DEAL::Matrix difference norm: 0
+DEAL::rhs difference norm: 0
+DEAL::Test matrix 2 nonzeros: 49920, actually: 42514
+DEAL::Matrix difference norm: 0
+DEAL::rhs difference norm: 0
+DEAL:cg::Starting value 28.
+DEAL:cg::Convergence step 48 value 3.5e-07
+DEAL::Distribute error: 0
+DEAL::
+DEAL::
+DEAL::   Number of active cells:       1173
+DEAL::   Number of degrees of freedom: 5732
+DEAL::   Number of constraints       : 492
+DEAL::Reference matrix nonzeros: 98166, actually: 87212
+DEAL::Test matrix 1 nonzeros: 98166, actually: 78146
+DEAL::Matrix difference norm: 0
+DEAL::rhs difference norm: 0
+DEAL::Test matrix 2 nonzeros: 98166, actually: 78146
+DEAL::Matrix difference norm: 0
+DEAL::rhs difference norm: 0
+DEAL:cg::Starting value 31.
+DEAL:cg::Convergence step 84 value 4.4e-07
+DEAL::Distribute error: 0
+DEAL::
+DEAL::
+DEAL::   Number of active cells:       1644
+DEAL::   Number of degrees of freedom: 9764
+DEAL::   Number of constraints       : 1434
+DEAL::Reference matrix nonzeros: 200326, actually: 161794
+DEAL::Test matrix 1 nonzeros: 200326, actually: 150182
+DEAL::Matrix difference norm: 0
+DEAL::rhs difference norm: 0
+DEAL::Test matrix 2 nonzeros: 200326, actually: 150182
+DEAL::Matrix difference norm: 0
+DEAL::rhs difference norm: 0
+DEAL:cg::Starting value 35.
+DEAL:cg::Convergence step 107 value 5.2e-07
+DEAL::Distribute error: 0