From ebc1ee3f357caa1677c27e187cf6698421de9efd Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Sat, 24 Oct 2015 12:27:04 -0500 Subject: [PATCH] Update test and rename to constraint_linear_operator_01 --- tests/lac/constraint_linear_operator_01.cc | 261 ++++++++++ ...nt_linear_operator_01.with_cxx11=on.output | 10 + tests/lac/linear_operator_09.cc | 482 ------------------ .../linear_operator_09.with_cxx11=on.output | 16 - 4 files changed, 271 insertions(+), 498 deletions(-) create mode 100644 tests/lac/constraint_linear_operator_01.cc create mode 100644 tests/lac/constraint_linear_operator_01.with_cxx11=on.output delete mode 100644 tests/lac/linear_operator_09.cc delete mode 100644 tests/lac/linear_operator_09.with_cxx11=on.output diff --git a/tests/lac/constraint_linear_operator_01.cc b/tests/lac/constraint_linear_operator_01.cc new file mode 100644 index 0000000000..188f8baa4d --- /dev/null +++ b/tests/lac/constraint_linear_operator_01.cc @@ -0,0 +1,261 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Set up a laplace problem with inhomogeneity 1 on the boundary, and run a +// cycle with local refinement to have hanging nodes. +// +// This test compares the solution vector obtained by the classical method +// of applying the constraints directly to the solution obtained by solving +// the system: +// (C^T * A * C + Id_c) * x = C^t * (b - A * k) +// with the help of constrainted_linear_operator and +// constrained_right_hand_side + +#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 + +using namespace dealii; + +template +class Step6 +{ +public: + Step6 (); + + void run (); + +private: + void setup_system (); + + void assemble_system (); + + void solve (); + + void refine_grid (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + ConstraintMatrix constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + SparseMatrix system_matrix_lo; + + Vector solution; + Vector solution_lo; + Vector system_rhs; + Vector system_rhs_lo; +}; + + +template +Step6::Step6 () + : + fe (2), + dof_handler (triangulation) +{} + + +template +void Step6::setup_system () +{ + dof_handler.distribute_dofs (fe); + + solution.reinit (dof_handler.n_dofs()); + solution_lo.reinit (dof_handler.n_dofs()); + + system_rhs.reinit (dof_handler.n_dofs()); + system_rhs_lo.reinit (dof_handler.n_dofs()); + + constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + constraints); + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ConstantFunction(1.), + constraints); + constraints.close (); + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit (sparsity_pattern); + system_matrix_lo.reinit (sparsity_pattern); +} + + +template +void Step6::assemble_system () +{ + const QGauss quadrature_formula(3); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + 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); + + 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); + + for (unsigned int q_index=0; q_indexget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + cell->distribute_local_to_global(cell_matrix, system_matrix_lo); + cell->distribute_local_to_global(cell_rhs, system_rhs_lo); + } +} + + +template +void Step6::solve () +{ + SolverControl solver_control (1000, 1e-8); + SolverCG<> solver (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + 10, + 60); + constraints.distribute (solution); + + const auto A = linear_operator(system_matrix_lo); + const auto M = constrained_linear_operator(constraints, A); + Vector rhs = + constrained_right_hand_side(constraints, A, system_rhs_lo); + + check_solver_within_range( + solver.solve(M, solution_lo, rhs, preconditioner), + solver_control.last_step(), + 10, + 60); + constraints.distribute(solution_lo); +} + + +template +void Step6::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::estimate (dof_handler, + QGauss(3), + typename FunctionMap::type(), + solution, + estimated_error_per_cell); + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.03); + + triangulation.execute_coarsening_and_refinement (); +} + + +template +void Step6::run () +{ + for (unsigned int cycle=0; cycle<3; ++cycle) + { + if (cycle == 0) + { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(3); + } + else + refine_grid(); + + setup_system(); + assemble_system(); + solve(); + + Vector diff = solution - solution_lo; + if (diff.l2_norm() < 1e-10) + deallog << "OK" << std::endl; + else + deallog << "ERROR! Output does not match!" << std::endl; + } +} + + +int main() +{ + initlog(); + deallog << std::setprecision(10); + + Step6<2> laplace_problem_2d; + laplace_problem_2d.run (); + + return 0; +} diff --git a/tests/lac/constraint_linear_operator_01.with_cxx11=on.output b/tests/lac/constraint_linear_operator_01.with_cxx11=on.output new file mode 100644 index 0000000000..30e009ec6f --- /dev/null +++ b/tests/lac/constraint_linear_operator_01.with_cxx11=on.output @@ -0,0 +1,10 @@ + +DEAL::Solver stopped within 10 - 60 iterations +DEAL::Solver stopped within 10 - 60 iterations +DEAL::OK +DEAL::Solver stopped within 10 - 60 iterations +DEAL::Solver stopped within 10 - 60 iterations +DEAL::OK +DEAL::Solver stopped within 10 - 60 iterations +DEAL::Solver stopped within 10 - 60 iterations +DEAL::OK diff --git a/tests/lac/linear_operator_09.cc b/tests/lac/linear_operator_09.cc deleted file mode 100644 index 6369b81617..0000000000 --- a/tests/lac/linear_operator_09.cc +++ /dev/null @@ -1,482 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2015 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - -// This test is based on step-6 -// The aim is to test constrained linear operators -// and compare them to the old implementation. - -#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 -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace dealii; - -template -class Step6 -{ -public: - Step6 (); - ~Step6 (); - - void run (); - -private: - void setup_system (); - - void assemble_system (); - void assemble_system_lo (); - - void solve (); - void solve_lo (); - - void refine_grid (); - - Triangulation triangulation; - Triangulation triangulation_lo; - - DoFHandler dof_handler; - DoFHandler dof_handler_lo; - - FE_Q fe; - FE_Q fe_lo; - - ConstraintMatrix constraints; - ConstraintMatrix constraints_lo; - - SparsityPattern sparsity_pattern; - SparsityPattern sparsity_pattern_lo; - - SparseMatrix system_matrix; - SparseMatrix system_matrix_lo; - - Vector solution; - Vector solution_lo; - - Vector system_rhs; - Vector system_rhs_lo; -}; - - - - -template -class Coefficient : public Function -{ -public: - Coefficient () : Function() {} - - virtual double value (const Point &p, - const unsigned int component = 0) const; - - virtual void value_list (const std::vector > &points, - std::vector &values, - const unsigned int component = 0) const; -}; - - - -template -double Coefficient::value (const Point &p, - const unsigned int) const -{ - if (p.square() < 0.5*0.5) - return 20; - else - return 1; -} - - - -template -void Coefficient::value_list (const std::vector > &points, - std::vector &values, - const unsigned int component) const -{ - const unsigned int n_points = points.size(); - - Assert (values.size() == n_points, - ExcDimensionMismatch (values.size(), n_points)); - - Assert (component == 0, - ExcIndexRange (component, 0, 1)); - - for (unsigned int i=0; i -Step6::Step6 () - : - dof_handler (triangulation), - dof_handler_lo (triangulation_lo), - fe (2), - fe_lo (2) -{} - -template -Step6::~Step6 () -{ - dof_handler.clear (); -} - - - -template -void Step6::setup_system () -{ - dof_handler.distribute_dofs (fe); - - solution.reinit (dof_handler.n_dofs()); - solution_lo.reinit (dof_handler.n_dofs()); - - system_rhs.reinit (dof_handler.n_dofs()); - system_rhs_lo.reinit (dof_handler.n_dofs()); - - constraints.clear (); - DoFTools::make_hanging_node_constraints (dof_handler, - constraints); - - - VectorTools::interpolate_boundary_values (dof_handler, - 0, - ConstantFunction(1.), - constraints); - - - constraints.close (); - - DynamicSparsityPattern dsp(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern(dof_handler, - dsp, - constraints, - /*keep_constrained_dofs = */ false); - - sparsity_pattern.copy_from(dsp); - - system_matrix.reinit (sparsity_pattern); - - constraints_lo.clear (); - constraints_lo.close (); - - DynamicSparsityPattern dsp_lo(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern(dof_handler, - dsp_lo, - constraints_lo, - /*keep_constrained_dofs = */ false); - - sparsity_pattern_lo.copy_from(dsp_lo); - - system_matrix_lo.reinit (sparsity_pattern_lo); -} - - - -template -void Step6::assemble_system () -{ - const QGauss quadrature_formula(3); - - FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - - 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); - - const Coefficient coefficient; - std::vector coefficient_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); - - coefficient.value_list (fe_values.get_quadrature_points(), - coefficient_values); - - for (unsigned int q_index=0; q_indexget_dof_indices (local_dof_indices); - constraints.distribute_local_to_global (cell_matrix, - cell_rhs, - local_dof_indices, - system_matrix, - system_rhs); - } -} - -template -void Step6::assemble_system_lo () -{ - const QGauss quadrature_formula(3); - - FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | - update_quadrature_points | update_JxW_values); - - 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); - - const Coefficient coefficient; - std::vector coefficient_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); - - coefficient.value_list (fe_values.get_quadrature_points(), - coefficient_values); - - for (unsigned int q_index=0; q_indexget_dof_indices (local_dof_indices); - constraints_lo.distribute_local_to_global (cell_matrix, - cell_rhs, - local_dof_indices, - system_matrix_lo, - system_rhs_lo); - } - -} - - -template -void Step6::solve () -{ - // constraints.condense(system_matrix, system_rhs); - unsigned int n_iterations = 0; - auto M = linear_operator(system_matrix); - - SolverControl solver_control (10000, 1e-12); - SolverCG<> solver (solver_control); - - PreconditionSSOR<> preconditioner; - preconditioner.initialize(system_matrix, 1.2); - - auto M_inv = inverse_operator( M, - solver, - preconditioner); - - - - M_inv.vmult(solution, system_rhs); - n_iterations = solver_control.last_step(); - constraints.distribute (solution); -} - - -template -void Step6::solve_lo () -{ - unsigned int n_iterations = 0; - - auto new_system_rhs_lo = constrained_rhs< >( - constraints, system_matrix_lo, system_rhs_lo);; - auto M = constrained_linear_operator< >( constraints, system_matrix_lo ); - - SolverControl solver_control (10000, 1e-12); - SolverCG<> solver (solver_control); - - PreconditionSSOR<> preconditioner; - preconditioner.initialize(system_matrix, 1.2); - - auto M_inv = inverse_operator( M, - solver, - preconditioner); - - M_inv.vmult(solution_lo, new_system_rhs_lo); - n_iterations = solver_control.last_step(); - - constraints.distribute(solution_lo); -} - - - -template -void Step6::refine_grid () -{ - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - KellyErrorEstimator::estimate (dof_handler, - QGauss(3), - typename FunctionMap::type(), - solution, - estimated_error_per_cell); - - GridRefinement::refine_and_coarsen_fixed_number (triangulation, - estimated_error_per_cell, - 0.3, 0.03); - - triangulation.execute_coarsening_and_refinement (); -} - -template -void Step6::run () -{ - for (unsigned int cycle=0; cycle<3; ++cycle) - { - if (cycle == 0) - { - GridGenerator::hyper_ball (triangulation); - - static const SphericalManifold boundary; - triangulation.set_all_manifold_ids_on_boundary(0); - triangulation.set_manifold (0, boundary); - - triangulation.refine_global (1); - } - else - refine_grid (); - - setup_system (); - - assemble_system (); - solve (); - - assemble_system_lo (); - solve_lo (); - - // auto rhs_tmp = system_rhs; - // new_system_rhs_lo.apply(rhs_tmp); - // compare_solution(system_rhs, rhs_tmp, dof_handler.n_dofs(), 1e-6, false, false, "RHS"); - Vector diff = solution-solution_lo; - if (diff.l2_norm() < 1e-10) - deallog << "OK" << std::endl; - } -} - -int main() -{ - initlog(); - deallog << std::setprecision(10); - - { - try - { - Step6<2> laplace_problem_2d; - laplace_problem_2d.run (); - } - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - return 0; - - deallog << "OK" << std::endl; - } -} diff --git a/tests/lac/linear_operator_09.with_cxx11=on.output b/tests/lac/linear_operator_09.with_cxx11=on.output deleted file mode 100644 index 80feeb7ee1..0000000000 --- a/tests/lac/linear_operator_09.with_cxx11=on.output +++ /dev/null @@ -1,16 +0,0 @@ - -DEAL:cg::Starting value 13.25516422 -DEAL:cg::Convergence step 24 value 1.400407060e-13 -DEAL:cg::Starting value 13.25516422 -DEAL:cg::Convergence step 24 value 1.399854635e-13 -DEAL::OK -DEAL:cg::Starting value 17.40605136 -DEAL:cg::Convergence step 33 value 3.786179069e-13 -DEAL:cg::Starting value 17.40605136 -DEAL:cg::Convergence step 33 value 3.789771041e-13 -DEAL::OK -DEAL:cg::Starting value 17.40587637 -DEAL:cg::Convergence step 56 value 4.293614863e-13 -DEAL:cg::Starting value 17.40587637 -DEAL:cg::Convergence step 56 value 4.295516840e-13 -DEAL::OK -- 2.39.5