From 51ff8f24790765de48405a546fb97eb3efb21aa5 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 25 Feb 2009 08:52:46 +0000 Subject: [PATCH] Add another test for inhomogeneous constraints. git-svn-id: https://svn.dealii.org/trunk@18424 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/deal.II/inhomogeneous_constraints.cc | 24 +- .../inhomogeneous_constraints/cmp/generic | 70 +-- .../inhomogeneous_constraints_nonsymmetric.cc | 488 ++++++++++++++++++ .../cmp/generic | 25 + 4 files changed, 560 insertions(+), 47 deletions(-) create mode 100644 tests/deal.II/inhomogeneous_constraints_nonsymmetric.cc create mode 100644 tests/deal.II/inhomogeneous_constraints_nonsymmetric/cmp/generic diff --git a/tests/deal.II/inhomogeneous_constraints.cc b/tests/deal.II/inhomogeneous_constraints.cc index ea796cca55..d0d58c7474 100644 --- a/tests/deal.II/inhomogeneous_constraints.cc +++ b/tests/deal.II/inhomogeneous_constraints.cc @@ -41,7 +41,6 @@ #include #include #include -#include #include #include #include @@ -214,20 +213,21 @@ void LaplaceProblem::test_equality () test->value() = 0; } - deallog << "Matrix difference norm: " + 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. + // nonzero rhs, whereas we will have zero + // rhs when using inhomogeneous + // constraints. for (unsigned int i=0; i::assemble_reference () } 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 boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 0, @@ -312,6 +308,10 @@ void LaplaceProblem::assemble_reference () reference_matrix, solution, reference_rhs); + + deallog << " Reference matrix nonzeros: " << reference_matrix.n_nonzero_elements() + << ", actually: " << reference_matrix.n_actually_nonzero_elements () + << std::endl; } @@ -380,7 +380,7 @@ void LaplaceProblem::assemble_test_1 () } test_all_constraints.condense (test_matrix, test_rhs); - deallog << "Test matrix 1 nonzeros: " << test_matrix.n_nonzero_elements() + deallog << " Test matrix 1 nonzeros: " << test_matrix.n_nonzero_elements() << ", actually: " << test_matrix.n_actually_nonzero_elements () << std::endl; @@ -452,7 +452,7 @@ void LaplaceProblem::assemble_test_2 () test_matrix, test_rhs); } - deallog << "Test matrix 2 nonzeros: " << test_matrix.n_nonzero_elements() + deallog << " Test matrix 2 nonzeros: " << test_matrix.n_nonzero_elements() << ", actually: " << test_matrix.n_actually_nonzero_elements () << std::endl; test_equality(); @@ -612,7 +612,7 @@ void LaplaceProblem<2>::create_coarse_grid () triangulation.create_triangulation (vertices, cells, SubCellData()); - triangulation.refine_global (3); + triangulation.refine_global (2); } diff --git a/tests/deal.II/inhomogeneous_constraints/cmp/generic b/tests/deal.II/inhomogeneous_constraints/cmp/generic index 91ffd93473..d85cc69215 100644 --- a/tests/deal.II/inhomogeneous_constraints/cmp/generic +++ b/tests/deal.II/inhomogeneous_constraints/cmp/generic @@ -1,46 +1,46 @@ DEAL:: DEAL:: -DEAL:: Number of active cells: 768 -DEAL:: Number of degrees of freedom: 3264 +DEAL:: Number of active cells: 192 +DEAL:: Number of degrees of freedom: 864 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:: Reference matrix nonzeros: 12672, actually: 8974 +DEAL:: Test matrix 1 nonzeros: 12672, actually: 8974 +DEAL:: Matrix difference norm: 0 +DEAL:: RHS difference norm: 0 +DEAL:: Test matrix 2 nonzeros: 12672, actually: 8974 +DEAL:: Matrix difference norm: 0 +DEAL:: RHS difference norm: 0 +DEAL:cg::Starting value 20. +DEAL:cg::Convergence step 25 value 2.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:: Number of active cells: 309 +DEAL:: Number of degrees of freedom: 1575 +DEAL:: Number of constraints : 187 +DEAL:: Reference matrix nonzeros: 25901, actually: 17333 +DEAL:: Test matrix 1 nonzeros: 25901, actually: 17333 +DEAL:: Matrix difference norm: 0 +DEAL:: RHS difference norm: 0 +DEAL:: Test matrix 2 nonzeros: 25901, actually: 17333 +DEAL:: Matrix difference norm: 0 +DEAL:: RHS difference norm: 0 +DEAL:cg::Starting value 22. +DEAL:cg::Convergence step 42 value 3.2e-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:: Number of active cells: 450 +DEAL:: Number of degrees of freedom: 2692 +DEAL:: Number of constraints : 475 +DEAL:: Reference matrix nonzeros: 51910, actually: 33908 +DEAL:: Test matrix 1 nonzeros: 51910, actually: 33908 +DEAL:: Matrix difference norm: 0 +DEAL:: RHS difference norm: 0 +DEAL:: Test matrix 2 nonzeros: 51910, actually: 33908 +DEAL:: Matrix difference norm: 0 +DEAL:: RHS difference norm: 0 +DEAL:cg::Starting value 25. +DEAL:cg::Convergence step 61 value 3.8e-07 DEAL::Distribute error: 0 diff --git a/tests/deal.II/inhomogeneous_constraints_nonsymmetric.cc b/tests/deal.II/inhomogeneous_constraints_nonsymmetric.cc new file mode 100644 index 0000000000..1591ca73c1 --- /dev/null +++ b/tests/deal.II/inhomogeneous_constraints_nonsymmetric.cc @@ -0,0 +1,488 @@ +//---------------- inhomogeneous_constraints_nonsymmetric.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_nonsymmetric.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 + +#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_nonsymmetric/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; + FE_Q fe; + + ConstraintMatrix hanging_nodes_only; + ConstraintMatrix test_all_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix reference_matrix; + SparseMatrix test_matrix; + + Vector reference_rhs; + Vector 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 (2) +{} + + +template +AdvectionProblem::~AdvectionProblem () +{ + dof_handler.clear (); +} + + +template +void AdvectionProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + + 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 boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + RightHandSide(), + 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 (); + + 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 +void AdvectionProblem::test_equality () +{ + // need to manually go through the + // matrix, since we can have different + // entries in constrained lines. + for (unsigned int i=0; i::const_iterator reference = reference_matrix.begin(i); + SparseMatrix::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 zero + // rhs when using inhomogeneous + // constraints. + for (unsigned int i=0; i +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, + RightHandSide(), + 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_nonsymmetric/cmp/generic b/tests/deal.II/inhomogeneous_constraints_nonsymmetric/cmp/generic new file mode 100644 index 0000000000..d6a10009e4 --- /dev/null +++ b/tests/deal.II/inhomogeneous_constraints_nonsymmetric/cmp/generic @@ -0,0 +1,25 @@ + +DEAL:: +DEAL:: +DEAL:: Number of active cells: 22 +DEAL:: Number of degrees of freedom: 114 +DEAL:: Number of constraints : 9 +DEAL:: Reference matrix nonzeros: 1618, actually: 701 +DEAL:: Test matrix 1 nonzeros: 1618, actually: 701 +DEAL:: Matrix difference norm: 0 +DEAL:: RHS difference norm: 0 +DEAL:: Test matrix 2 nonzeros: 1618, actually: 698 +DEAL:: Matrix difference norm: 0 +DEAL:: RHS difference norm: 0 +DEAL:: +DEAL:: +DEAL:: Number of active cells: 78 +DEAL:: Number of degrees of freedom: 928 +DEAL:: Number of constraints : 87 +DEAL:: Reference matrix nonzeros: 46980, actually: 14420 +DEAL:: Test matrix 1 nonzeros: 46980, actually: 14420 +DEAL:: Matrix difference norm: 0 +DEAL:: RHS difference norm: 0 +DEAL:: Test matrix 2 nonzeros: 46980, actually: 14416 +DEAL:: Matrix difference norm: 0 +DEAL:: RHS difference norm: 0 -- 2.39.5