From d0a50788c1deddf023ec77024f67f7adac38cae4 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +std::ofstream logfile("inhomogeneous_constraints_block/output"); + +using namespace dealii; + +template +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 triangulation; + + DoFHandler dof_handler; + FESystem fe; + + ConstraintMatrix hanging_nodes_only; + ConstraintMatrix test_all_constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix reference_matrix; + BlockSparseMatrix test_matrix; + + BlockVector reference_rhs; + BlockVector test_rhs; +}; + + + +template +class RightHandSide : public Function +{ + public: + RightHandSide () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double product = 1; + for (unsigned int d=0; d +AdvectionProblem::AdvectionProblem () + : + dof_handler (triangulation), + fe (FE_Q(2),2) +{} + + +template +AdvectionProblem::~AdvectionProblem () +{ + dof_handler.clear (); +} + + +template +void AdvectionProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + + hanging_nodes_only.clear (); + test_all_constraints.clear (); + + // add boundary conditions as + // inhomogeneous constraints here. + { + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ConstantFunction(1.,2), + boundary_values); + std::map::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 +void AdvectionProblem::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::const_iterator reference = + reference_matrix.block(block_row,block_col).begin(index_in_block); + SparseMatrix::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 +void AdvectionProblem::assemble_reference () +{ + reference_matrix = 0; + reference_rhs = 0; + + QGauss quadrature_formula (3); + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const RightHandSide rhs_function; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + std::vector rhs_values (n_q_points); + + typename DoFHandler::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_pointget_dof_indices (local_dof_indices); + + reference_matrix.add(local_dof_indices, cell_matrix); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ConstantFunction(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 +void AdvectionProblem::assemble_test_1 () +{ + test_matrix = 0; + test_rhs = 0; + + + QGauss quadrature_formula (3); + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const RightHandSide rhs_function; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + std::vector rhs_values (n_q_points); + + typename DoFHandler::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_pointget_dof_indices (local_dof_indices); + + test_matrix.add(local_dof_indices, cell_matrix); + for (unsigned int i=0; i +void AdvectionProblem::assemble_test_2 () +{ + test_matrix = 0; + test_rhs = 0; + + QGauss quadrature_formula (3); + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const RightHandSide rhs_function; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + std::vector rhs_values (n_q_points); + + typename DoFHandler::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_pointget_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 +void AdvectionProblem::run () +{ + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (2); + + // manually refine the first two cells to + // create some hanging nodes + { + typename DoFHandler::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::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 boundary_values; VectorTools::interpolate_boundary_values (dof_handler, -- 2.39.5