From d0eb39dcead855fd62da184fb7edd45760ad1ea8 Mon Sep 17 00:00:00 2001
From: kronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 25 Feb 2009 09:52:20 +0000
Subject: [PATCH] Add a final test for inhomogeneous constraints.

git-svn-id: https://svn.dealii.org/trunk@18425 0785d39b-7218-0410-832d-ea1e28bc413d
---
 tests/deal.II/Makefile                        |   4 +-
 .../inhomogeneous_constraints_block.cc        | 534 ++++++++++++++++++
 .../cmp/generic                               |  25 +
 .../inhomogeneous_constraints_nonsymmetric.cc |   6 +-
 4 files changed, 565 insertions(+), 4 deletions(-)
 create mode 100644 tests/deal.II/inhomogeneous_constraints_block.cc
 create mode 100644 tests/deal.II/inhomogeneous_constraints_block/cmp/generic

diff --git a/tests/deal.II/Makefile b/tests/deal.II/Makefile
index 4e2d3e84a9..6e3edb1eea 100644
--- a/tests/deal.II/Makefile
+++ b/tests/deal.II/Makefile
@@ -24,12 +24,14 @@ default: run-tests
 tests_x = block_matrices \
 	user_data_* \
 	constraints \
-	inhomogeneous_constraints \
 	constraints_zero \
 	constraints_zero_merge \
 	constraints_zero_condense \
 	constraint_graph \
 	constraint_graph_zero \
+	inhomogeneous_constraints \
+	inhomogeneous_constraints_nonsymmetric \
+	inhomogeneous_constraints_block \
 	data_out \
 	derivative_* \
 	derivatives \
diff --git a/tests/deal.II/inhomogeneous_constraints_block.cc b/tests/deal.II/inhomogeneous_constraints_block.cc
new file mode 100644
index 0000000000..9a0db89e65
--- /dev/null
+++ b/tests/deal.II/inhomogeneous_constraints_block.cc
@@ -0,0 +1,534 @@
+//----------------  inhomogeneous_constraints_block.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_block.cc  -------------------
+
+
+// this function tests the correctness of the implementation of
+// inhomogeneous constraints on a nonsymmetric matrix that comes from a
+// discretization of the first derivative, based on block matrices instead
+// of standard matrices, by working on a vector-valued problem
+
+#include "../tests.h"
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <base/utilities.h>
+#include <lac/block_vector.h>
+#include <lac/full_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/error_estimator.h>
+#include <lac/block_sparsity_pattern.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+std::ofstream logfile("inhomogeneous_constraints_block/output");
+
+using namespace dealii;
+
+template <int dim>
+class AdvectionProblem 
+{
+  public:
+    AdvectionProblem ();
+    ~AdvectionProblem ();
+
+    void run ();
+    
+  private:
+    void setup_system ();
+    void test_equality ();
+    void assemble_reference ();
+    void assemble_test_1 ();
+    void assemble_test_2 ();
+
+    Triangulation<dim>        triangulation;
+
+    DoFHandler<dim>           dof_handler;
+    FESystem<dim>             fe;
+
+    ConstraintMatrix          hanging_nodes_only;
+    ConstraintMatrix          test_all_constraints;
+
+    BlockSparsityPattern      sparsity_pattern;
+    BlockSparseMatrix<double> reference_matrix;
+    BlockSparseMatrix<double> test_matrix;
+
+    BlockVector<double>       reference_rhs;
+    BlockVector<double>       test_rhs;
+};
+
+
+
+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>
+AdvectionProblem<dim>::AdvectionProblem ()
+		:
+                dof_handler (triangulation),
+		fe (FE_Q<dim>(2),2)
+{}
+
+
+template <int dim>
+AdvectionProblem<dim>::~AdvectionProblem () 
+{
+  dof_handler.clear ();
+}
+
+
+template <int dim>
+void AdvectionProblem<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+  hanging_nodes_only.clear ();
+  test_all_constraints.clear ();
+
+				   // add boundary conditions as
+				   // inhomogeneous constraints here.
+  {
+    std::map<unsigned int,double> boundary_values;
+    VectorTools::interpolate_boundary_values (dof_handler,
+					      0,
+					      ConstantFunction<dim>(1.,2),
+					      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 ();
+
+  BlockCompressedSimpleSparsityPattern csp (2,2);
+  {
+    const unsigned int dofs_per_block = dof_handler.n_dofs() / 2;
+    csp.block(0,0).reinit (dofs_per_block, dofs_per_block);
+    csp.block(0,1).reinit (dofs_per_block, dofs_per_block);
+    csp.block(1,0).reinit (dofs_per_block, dofs_per_block);
+    csp.block(1,1).reinit (dofs_per_block, dofs_per_block);
+    csp.collect_sizes();
+  }
+
+  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);
+
+  reference_rhs.reinit (2);
+  reference_rhs.block(0).reinit (dof_handler.n_dofs() / 2);
+  reference_rhs.block(1).reinit (dof_handler.n_dofs() / 2);
+  reference_rhs.collect_sizes();
+  test_rhs.reinit (reference_rhs);
+}
+
+
+
+				   // test whether we are equal with the
+				   // standard matrix and right hand side
+template <int dim>
+void AdvectionProblem<dim>::test_equality ()
+{
+				   // need to manually go through the
+				   // matrix, since we can have different
+				   // entries in constrained lines because
+				   // the diagonal is set differently
+
+  const BlockIndices &
+    index_mapping = sparsity_pattern.get_column_indices();
+
+  for (unsigned int i=0; i<reference_matrix.m(); ++i)
+    {
+      const unsigned int block_row = index_mapping.global_to_local(i).first;
+      const unsigned int index_in_block = index_mapping.global_to_local(i).second;
+      for (unsigned int block_col=0; block_col<sparsity_pattern.n_block_cols(); 
+	   ++block_col)
+	{
+	  SparseMatrix<double>::const_iterator reference = 
+	    reference_matrix.block(block_row,block_col).begin(index_in_block);
+	  SparseMatrix<double>::iterator test = 
+	    test_matrix.block(block_row,block_col).begin(index_in_block);
+	  if (test_all_constraints.is_constrained(i) == false)
+	    {
+	      for ( ; test != test_matrix.block(block_row,block_col).end(index_in_block); 
+		    ++test, ++reference)
+		test->value() -= reference->value();
+	    }
+	  else
+	    for ( ; test != test_matrix.block(block_row,block_col).end(index_in_block); 
+		  ++test)
+	      test->value() = 0;
+	}
+    }
+
+  double frobenius_norm = 0.;
+  for (unsigned int row=0;row<sparsity_pattern.n_block_rows(); ++row)
+    for (unsigned int col=0;col<sparsity_pattern.n_block_cols(); ++col)
+      frobenius_norm += test_matrix.block(row,col).frobenius_norm() *
+	test_matrix.block(row,col).frobenius_norm();
+  frobenius_norm = std::sqrt(frobenius_norm);
+
+  deallog << "  Matrix difference norm: " << frobenius_norm << std::endl;
+  Assert (frobenius_norm < 1e-13, ExcInternalError());
+
+				   // same here -- Dirichlet lines will have
+				   // nonzero rhs, whereas we will have zero
+				   // rhs when using inhomogeneous
+				   // constraints.
+  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 AdvectionProblem<dim>::assemble_reference () 
+{
+  reference_matrix = 0;
+  reference_rhs = 0;
+
+  QGauss<dim> quadrature_formula (3);
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+			   update_values    |  update_gradients | 
+			   update_quadrature_points  |  update_JxW_values);
+
+  const RightHandSide<dim> rhs_function;
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points = quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  std::vector<double>  rhs_values (n_q_points);
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      cell_rhs = 0;
+      fe_values.reinit (cell);
+
+      rhs_function.value_list (fe_values.get_quadrature_points(),
+			       rhs_values);
+
+      Tensor<1,dim> advection_direction;
+      advection_direction[0] = 1;
+      advection_direction[1] = 1;
+      advection_direction[dim-1] = -1;
+      
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+	for (unsigned int i=0; i<dofs_per_cell; ++i)
+	  {
+	    const unsigned int comp_i = fe.system_to_component_index(i).first;
+	    for (unsigned int j=0; j<dofs_per_cell; ++j)
+	      {
+		const unsigned int comp_j = fe.system_to_component_index(j).first;
+		if (comp_i == comp_j)
+		  cell_matrix(i,j) += (fe_values.shape_value(i,q_point) *
+				       advection_direction *
+				       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);
+
+				   // use some other rhs vector as dummy for
+				   // application of Dirichlet conditions
+  std::map<unsigned int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+					    0,
+					    ConstantFunction<dim>(1.,2),
+					    boundary_values);
+  MatrixTools::apply_boundary_values (boundary_values,
+				      reference_matrix,
+				      test_rhs,
+				      reference_rhs);
+
+  deallog << "  Reference matrix nonzeros: " << reference_matrix.n_nonzero_elements() 
+	  << ", actually: " << reference_matrix.n_actually_nonzero_elements () 
+	  << std::endl;
+}
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::assemble_test_1 () 
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+
+  QGauss<dim> quadrature_formula (3);
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+			   update_values    |  update_gradients | 
+			   update_quadrature_points  |  update_JxW_values);
+
+  const RightHandSide<dim> rhs_function;
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points = quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  std::vector<double>  rhs_values (n_q_points);
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      cell_rhs = 0;
+      fe_values.reinit (cell);
+
+      rhs_function.value_list (fe_values.get_quadrature_points(),
+			       rhs_values);
+
+      Tensor<1,dim> advection_direction;
+      advection_direction[0] = 1;
+      advection_direction[1] = 1;
+      advection_direction[dim-1] = -1;
+      
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+	for (unsigned int i=0; i<dofs_per_cell; ++i)
+	  {
+	    const unsigned int comp_i = fe.system_to_component_index(i).first;
+	    for (unsigned int j=0; j<dofs_per_cell; ++j)
+	      {
+		const unsigned int comp_j = fe.system_to_component_index(j).first;
+		if (comp_i == comp_j)
+		  cell_matrix(i,j) += (fe_values.shape_value(i,q_point) *
+				       advection_direction *
+				       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 AdvectionProblem<dim>::assemble_test_2 () 
+{
+  test_matrix = 0;
+  test_rhs = 0;
+
+  QGauss<dim> quadrature_formula (3);
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+			   update_values    |  update_gradients | 
+			   update_quadrature_points  |  update_JxW_values);
+
+  const RightHandSide<dim> rhs_function;
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points = quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  std::vector<double>  rhs_values (n_q_points);
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      cell_matrix = 0;
+      cell_rhs = 0;
+      fe_values.reinit (cell);
+
+      rhs_function.value_list (fe_values.get_quadrature_points(),
+			       rhs_values);
+
+      Tensor<1,dim> advection_direction;
+      advection_direction[0] = 1;
+      advection_direction[1] = 1;
+      advection_direction[dim-1] = -1;
+      
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+	for (unsigned int i=0; i<dofs_per_cell; ++i)
+	  {
+	    const unsigned int comp_i = fe.system_to_component_index(i).first;
+	    for (unsigned int j=0; j<dofs_per_cell; ++j)
+	      {
+		const unsigned int comp_j = fe.system_to_component_index(j).first;
+		if (comp_i == comp_j)
+		  cell_matrix(i,j) += (fe_values.shape_value(i,q_point) *
+				       advection_direction *
+				       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 AdvectionProblem<dim>::run () 
+{
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global (2);
+
+				   // manually refine the first two cells to
+				   // create some hanging nodes
+  {
+    typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active();
+    cell->set_refine_flag();
+    cell++;
+    cell->set_refine_flag();
+  }
+  triangulation.execute_coarsening_and_refinement();  
+
+  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 ();
+}
+
+
+
+int main () 
+{
+  deallog << std::setprecision (2);
+  logfile << std::setprecision (2);
+  deallog.attach(logfile);
+  //deallog.depth_console(0);
+  deallog.threshold_double(1.e-12);
+
+  {
+    AdvectionProblem<2> advection_problem;
+    advection_problem.run ();
+  }
+  {
+    AdvectionProblem<3> advection_problem;
+    advection_problem.run ();
+  }
+}
diff --git a/tests/deal.II/inhomogeneous_constraints_block/cmp/generic b/tests/deal.II/inhomogeneous_constraints_block/cmp/generic
new file mode 100644
index 0000000000..b289fdd244
--- /dev/null
+++ b/tests/deal.II/inhomogeneous_constraints_block/cmp/generic
@@ -0,0 +1,25 @@
+
+DEAL::
+DEAL::
+DEAL::  Number of active cells:       22
+DEAL::  Number of degrees of freedom: 228
+DEAL::  Number of constraints       : 18
+DEAL::  Reference matrix nonzeros: 6700, actually: 1402
+DEAL::  Test matrix 1 nonzeros: 6700, actually: 1402
+DEAL::  Matrix difference norm: 0
+DEAL::  RHS difference norm: 0
+DEAL::  Test matrix 2 nonzeros: 6700, actually: 1396
+DEAL::  Matrix difference norm: 0
+DEAL::  RHS difference norm: 0
+DEAL::
+DEAL::
+DEAL::  Number of active cells:       78
+DEAL::  Number of degrees of freedom: 1856
+DEAL::  Number of constraints       : 174
+DEAL::  Reference matrix nonzeros: 189772, actually: 28840
+DEAL::  Test matrix 1 nonzeros: 189772, actually: 28840
+DEAL::  Matrix difference norm: 0
+DEAL::  RHS difference norm: 0
+DEAL::  Test matrix 2 nonzeros: 189772, actually: 28832
+DEAL::  Matrix difference norm: 0
+DEAL::  RHS difference norm: 0
diff --git a/tests/deal.II/inhomogeneous_constraints_nonsymmetric.cc b/tests/deal.II/inhomogeneous_constraints_nonsymmetric.cc
index 1591ca73c1..58635109dd 100644
--- a/tests/deal.II/inhomogeneous_constraints_nonsymmetric.cc
+++ b/tests/deal.II/inhomogeneous_constraints_nonsymmetric.cc
@@ -130,9 +130,9 @@ void AdvectionProblem<dim>::setup_system ()
   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.
+				   // inhomogeneous constraints here. just
+				   // take the right hand side function as
+				   // boundary function
   {
     std::map<unsigned int,double> boundary_values;
     VectorTools::interpolate_boundary_values (dof_handler,
-- 
2.39.5