From 408d62ac901bb950e3c930b31e4e539c2bb9f313 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Tue, 24 Feb 2009 16:33:53 +0000 Subject: [PATCH] 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 --- tests/deal.II/inhomogeneous_constraints.cc | 801 ++++++++++++++++++ .../inhomogeneous_constraints/cmp/generic | 46 + 2 files changed, 847 insertions(+) create mode 100644 tests/deal.II/inhomogeneous_constraints.cc create mode 100644 tests/deal.II/inhomogeneous_constraints/cmp/generic 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 +#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 + +std::ofstream logfile("inhomogeneous_constraints/output"); + +using namespace dealii; + +template +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 &smoothness_indicators) const; + void postprocess (const unsigned int cycle); + + Triangulation triangulation; + + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; + hp::QCollection quadrature_collection; + hp::QCollection face_quadrature_collection; + + ConstraintMatrix hanging_nodes_only; + ConstraintMatrix test_all_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix reference_matrix; + SparseMatrix test_matrix; + + Vector solution; + Vector reference_rhs; + Vector test_rhs; + + const unsigned int max_degree; +}; + + + +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 +LaplaceProblem::LaplaceProblem () + : + dof_handler (triangulation), + max_degree (5) +{ + for (unsigned int degree=2; degree<=max_degree; ++degree) + { + fe_collection.push_back (FE_Q(degree)); + quadrature_collection.push_back (QGauss(degree+1)); + face_quadrature_collection.push_back (QGauss(degree+1)); + } +} + + +template +LaplaceProblem::~LaplaceProblem () +{ + dof_handler.clear (); +} + + +template +void LaplaceProblem::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 boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ConstantFunction(1.), + 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 LaplaceProblem::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 + // nonzero one. + for (unsigned int i=0; i +void LaplaceProblem::assemble_reference () +{ + reference_matrix = 0; + reference_rhs = 0; + + hp::FEValues hp_fe_values (fe_collection, + quadrature_collection, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const RightHandSide rhs_function; + + FullMatrix cell_matrix; + Vector cell_rhs; + + std::vector local_dof_indices; + + typename hp::DoFHandler::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 &fe_values = hp_fe_values.get_present_fe_values (); + + std::vector 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_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.), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + reference_matrix, + solution, + reference_rhs); +} + + + +template +void LaplaceProblem::assemble_test_1 () +{ + test_matrix = 0; + test_rhs = 0; + + hp::FEValues hp_fe_values (fe_collection, + quadrature_collection, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const RightHandSide rhs_function; + + FullMatrix cell_matrix; + Vector cell_rhs; + + std::vector local_dof_indices; + + typename hp::DoFHandler::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 &fe_values = hp_fe_values.get_present_fe_values (); + + std::vector 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_pointget_dof_indices (local_dof_indices); + + test_matrix.add(local_dof_indices, cell_matrix); + for (unsigned int i=0; i +void LaplaceProblem::assemble_test_2 () +{ + test_matrix = 0; + test_rhs = 0; + + hp::FEValues hp_fe_values (fe_collection, + quadrature_collection, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const RightHandSide rhs_function; + + FullMatrix cell_matrix; + Vector cell_rhs; + + std::vector local_dof_indices; + + typename hp::DoFHandler::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 &fe_values = hp_fe_values.get_present_fe_values (); + + std::vector 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_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 LaplaceProblem::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 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 +void LaplaceProblem::postprocess (const unsigned int cycle) +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + KellyErrorEstimator::estimate (dof_handler, + face_quadrature_collection, + typename FunctionMap::type(), + solution, + estimated_error_per_cell); + + Vector 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::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::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 > vertices (&vertices_1[0], + &vertices_1[n_vertices]); + static const int cell_vertices[][GeometryInfo::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 > cells (n_cells, CellData()); + for (unsigned int i=0; i::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 +void LaplaceProblem::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 +void +LaplaceProblem:: +estimate_smoothness (Vector &smoothness_indicators) const +{ + const unsigned int N = max_degree; + + std::vector > k_vectors; + std::vector k_vectors_magnitude; + switch (dim) + { + case 2: + { + for (unsigned int i=0; i(numbers::PI * i, + numbers::PI * j)); + k_vectors_magnitude.push_back (i*i+j*j); + } + + break; + } + + case 3: + { + for (unsigned int i=0; i(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 ln_k (n_fourier_modes); + for (unsigned int i=0; i > > + fourier_transform_matrices (fe_collection.size()); + QGauss<1> base_quadrature (2); + QIterated quadrature (base_quadrature, N); + + for (unsigned int fe=0; fe sum = 0; + for (unsigned int q=0; q x_q = quadrature.point(q); + sum += std::exp(std::complex(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 > fourier_coefficients (n_fourier_modes); + Vector local_dof_values; + + typename hp::DoFHandler::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; fget_fe().dofs_per_cell; ++i) + fourier_coefficients[f] += + fourier_transform_matrices[cell->active_fe_index()](f,i) + * + local_dof_values(i); + } + + std::map k_to_max_U_map; + for (unsigned int f=0; f 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 -- 2.39.5