From: Martin Kronbichler Date: Thu, 26 Apr 2018 15:57:47 +0000 (+0200) Subject: Add tests for multiple components in single DoFHandler. X-Git-Tag: v9.0.0-rc1~93^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=34ad8e5fea44d5ed61ecae042e05f5da2334d971;p=dealii.git Add tests for multiple components in single DoFHandler. --- diff --git a/tests/matrix_free/matrix_vector_stokes_flux.cc b/tests/matrix_free/matrix_vector_stokes_flux.cc new file mode 100644 index 0000000000..4963b0882e --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_flux.cc @@ -0,0 +1,305 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2013 - 2016 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 function tests the correctness of the implementation of matrix free +// matrix-vector products by comparing with the result of a deal.II sparse +// matrix. The mesh uses a hypershell mesh with hanging nodes and constraints +// between the vector components in the form of normal flux constraints on +// the Stokes equations. Like matrix_vector_stokes_onedof except for +// different constraints + +#include "../tests.h" + +std::ofstream logfile("output"); + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + + +template +class MatrixFreeTest +{ +public: + typedef typename DoFHandler::active_cell_iterator CellIterator; + typedef double Number; + + MatrixFreeTest(const MatrixFree &data_in): + data (data_in) + {} + + void + local_apply (const MatrixFree &data, + VectorType &dst, + const VectorType &src, + const std::pair &cell_range) const + { + typedef VectorizedArray vector_t; + FEEvaluation velocity (data, 0, 0, 0); + FEEvaluation pressure (data, 0, 0, dim); + + for (unsigned int cell=cell_range.first; cell sym_grad_u = + velocity.get_symmetric_gradient (q); + vector_t pres = pressure.get_value(q); + vector_t div = -trace(sym_grad_u); + pressure.submit_value (div, q); + + // subtract p * I + for (unsigned int d=0; d::local_apply, + this, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void test () +{ + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation, 0, 1, true); + triangulation.refine_global(5-dim); + + FESystem fe (FE_Q(fe_degree+1), dim, FE_Q(fe_degree), 1); + DoFHandler dof_handler (triangulation); + + MatrixFree mf_data; + + ConstraintMatrix constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + Vector mf_solution, mf_rhs; + + dof_handler.distribute_dofs (fe); + std::vector stokes_sub_blocks (dim+1,0); + stokes_sub_blocks[dim] = 1; + DoFRenumbering::component_wise (dof_handler, stokes_sub_blocks); + + std::set no_normal_flux_boundaries; + no_normal_flux_boundaries.insert (0); + no_normal_flux_boundaries.insert (1); + DoFTools::make_hanging_node_constraints (dof_handler, + constraints); + VectorTools::compute_normal_flux_constraints (dof_handler, 0, + no_normal_flux_boundaries, + constraints); + constraints.close (); + + std::vector dofs_per_block (2); + DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, + stokes_sub_blocks); + { + BlockDynamicSparsityPattern csp (2,2); + + for (unsigned int d=0; d<2; ++d) + for (unsigned int e=0; e<2; ++e) + csp.block(d,e).reinit (dofs_per_block[d], dofs_per_block[e]); + + csp.collect_sizes(); + + DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false); + sparsity_pattern.copy_from (csp); + } + + system_matrix.reinit (sparsity_pattern); + + // this is from step-22 + { + QGauss quadrature_formula(fe_degree+2); + + FEValues fe_values (fe, quadrature_formula, + update_values | + update_JxW_values | + update_gradients); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + std::vector > phi_grads_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + local_matrix = 0; + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (local_matrix, + local_dof_indices, + system_matrix); + } + } + + + solution.reinit (2); + for (unsigned int d=0; d<2; ++d) + solution.block(d).reinit (dofs_per_block[d]); + solution.collect_sizes (); + + system_rhs.reinit (solution); + + // fill system_rhs with random numbers + for (unsigned int j=0; j quad(fe_degree+2); + // no parallelism + mf_data.reinit (dof_handler, constraints, quad, + typename MatrixFree::AdditionalData + (MatrixFree::AdditionalData::none)); + mf_data.initialize_dof_vector(mf_solution); + mf_data.initialize_dof_vector(mf_rhs); + for (unsigned int i=0; i > mf (mf_data); + mf.vmult (mf_solution, mf_rhs); + + // Verification + for (unsigned int i=0; i(); + test<2,2>(); + test<2,3>(); + deallog.pop(); + deallog.push("3d"); + test<3,1>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/matrix_vector_stokes_flux.output b/tests/matrix_free/matrix_vector_stokes_flux.output new file mode 100644 index 0000000000..0434059632 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_flux.output @@ -0,0 +1,12 @@ + +DEAL:: +DEAL::Test with doubles +DEAL:: +DEAL:2d::Verification fe degree 1: 0 +DEAL:2d:: +DEAL:2d::Verification fe degree 2: 0 +DEAL:2d:: +DEAL:2d::Verification fe degree 3: 0 +DEAL:2d:: +DEAL:3d::Verification fe degree 1: 0 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_stokes_noflux.cc b/tests/matrix_free/matrix_vector_stokes_noflux.cc index 84c9698d39..c33a019452 100644 --- a/tests/matrix_free/matrix_vector_stokes_noflux.cc +++ b/tests/matrix_free/matrix_vector_stokes_noflux.cc @@ -33,7 +33,6 @@ std::ofstream logfile("output"); #include #include #include -#include #include #include #include diff --git a/tests/matrix_free/matrix_vector_stokes_onedof.cc b/tests/matrix_free/matrix_vector_stokes_onedof.cc new file mode 100644 index 0000000000..d4576672e5 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_onedof.cc @@ -0,0 +1,306 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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. +// +// --------------------------------------------------------------------- + + + +// correctness matrix free matrix-vector products by comparing with the result +// of a deal.II sparse matrix. Similar to matrix_vector_stokes_noflux but +// putting all degrees of freedom into a single DoFHandler, where the +// selection is done through FEEvaluation + +#include "../tests.h" + +std::ofstream logfile("output"); + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + + +template +class MatrixFreeTest +{ +public: + typedef typename DoFHandler::active_cell_iterator CellIterator; + typedef double Number; + + MatrixFreeTest(const MatrixFree &data_in): + data (data_in) + {}; + + void + local_apply (const MatrixFree &data, + VectorType &dst, + const VectorType &src, + const std::pair &cell_range) const + { + typedef VectorizedArray vector_t; + FEEvaluation velocity (data, 0, 0, 0); + FEEvaluation pressure (data, 0, 0, dim); + + for (unsigned int cell=cell_range.first; cell sym_grad_u = + velocity.get_symmetric_gradient (q); + vector_t pres = pressure.get_value(q); + vector_t div = -trace(sym_grad_u); + pressure.submit_value (div, q); + + // subtract p * I + for (unsigned int d=0; d::local_apply, + this, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void test (const FESystem &fe) +{ + SphericalManifold manifold; + Triangulation triangulation; + GridGenerator::hyper_shell (triangulation, Point(), + 0.5, 1., 96, true); + triangulation.set_all_manifold_ids(0); + triangulation.set_manifold (0, manifold); + triangulation.begin_active()->set_refine_flag(); + triangulation.last()->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + triangulation.refine_global (3-dim); + triangulation.last()->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + + MappingQ mapping (3); + DoFHandler dof_handler (triangulation); + + MatrixFree mf_data; + + ConstraintMatrix constraints, constraints_u, constraints_p; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + Vector mf_solution; + + dof_handler.distribute_dofs (fe); + + std::set no_normal_flux_boundaries; + no_normal_flux_boundaries.insert (0); + no_normal_flux_boundaries.insert (1); + DoFTools::make_hanging_node_constraints (dof_handler, + constraints); + if (fe.dofs_per_vertex > 0) + VectorTools::compute_no_normal_flux_constraints (dof_handler, 0, + no_normal_flux_boundaries, + constraints, mapping); + constraints.close (); + + //std::cout << "Number of active cells: " + // << triangulation.n_active_cells() + // << std::endl + // << "Number of degrees of freedom: " + // << dof_handler.n_dofs() + // << " (" << n_u << '+' << n_p << ')' + // << std::endl; + + { + DynamicSparsityPattern dsp (dof_handler.n_dofs(),dof_handler.n_dofs()); + + DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false); + sparsity_pattern.copy_from (dsp); + } + + system_matrix.reinit (sparsity_pattern); + + // this is from step-22 + { + QGauss quadrature_formula(fe_degree+2); + + FEValues fe_values (mapping, fe, quadrature_formula, + update_values | + update_JxW_values | + update_gradients); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + std::vector > phi_grads_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + local_matrix = 0; + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (local_matrix, + local_dof_indices, + system_matrix); + } + } + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit (solution); + mf_solution.reinit (solution); + + // fill system_rhs with random numbers + for (unsigned int j=0; j quad(fe_degree+2); + // no parallelism + mf_data.reinit (mapping, dof_handler, constraints, quad, + typename MatrixFree::AdditionalData + (MatrixFree::AdditionalData::none)); + } + + system_matrix.vmult (solution, system_rhs); + + MatrixFreeTest > mf (mf_data); + mf.vmult (mf_solution, system_rhs); + + // Verification + mf_solution -= solution; + const double error = mf_solution.linfty_norm(); + const double relative = solution.linfty_norm(); + deallog << "Verification " << fe.get_name() << ": " + << error/relative << std::endl << std::endl; +} + + + +int main () +{ + initlog(); + + { + deallog << std::endl << "Test with doubles" << std::endl << std::endl; + deallog.push("2d"); + test<2,1>(FESystem<2>(FE_Q<2>(2),2, + FE_Q<2>(1),1)); + test<2,2>(FESystem<2>(FE_Q<2>(3),2, + FE_Q<2>(2),1)); + test<2,3>(FESystem<2>(FE_Q<2>(4),2, + FE_Q<2>(3),1)); + test<2,1>(FESystem<2>(FE_DGQ<2>(2),2, + FE_DGQ<2>(1),1)); + test<2,0>(FESystem<2>(FE_DGQ<2>(1),2, + FE_DGQ<2>(0),1)); + test<2,3>(FESystem<2>(FE_DGQ<2>(4),2, + FE_DGQ<2>(3),1)); + deallog.pop(); + deallog.push("3d"); + test<3,1>(FESystem<3>(FE_Q<3>(2),3, + FE_Q<3>(1),1)); + test<3,1>(FESystem<3>(FE_DGQ<3>(2),3, + FE_DGQ<3>(1),1)); + deallog.pop(); + } +} diff --git a/tests/matrix_free/matrix_vector_stokes_onedof.output b/tests/matrix_free/matrix_vector_stokes_onedof.output new file mode 100644 index 0000000000..ccae3ce6fd --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_onedof.output @@ -0,0 +1,20 @@ + +DEAL:: +DEAL::Test with doubles +DEAL:: +DEAL:2d::Verification FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(1)]: 0 +DEAL:2d:: +DEAL:2d::Verification FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(2)]: 0 +DEAL:2d:: +DEAL:2d::Verification FESystem<2>[FE_Q<2>(4)^2-FE_Q<2>(3)]: 0 +DEAL:2d:: +DEAL:2d::Verification FESystem<2>[FE_DGQ<2>(2)^2-FE_DGQ<2>(1)]: 0 +DEAL:2d:: +DEAL:2d::Verification FESystem<2>[FE_DGQ<2>(1)^2-FE_DGQ<2>(0)]: 0 +DEAL:2d:: +DEAL:2d::Verification FESystem<2>[FE_DGQ<2>(4)^2-FE_DGQ<2>(3)]: 0 +DEAL:2d:: +DEAL:3d::Verification FESystem<3>[FE_Q<3>(2)^3-FE_Q<3>(1)]: 0 +DEAL:3d:: +DEAL:3d::Verification FESystem<3>[FE_DGQ<3>(2)^3-FE_DGQ<3>(1)]: 0 +DEAL:3d::