From 3f06b79376f9f342994b0c0a686b31015a621e36 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 16 Feb 2017 12:18:03 +0100 Subject: [PATCH] Tests for FEEvaluation without template on polynomial degree --- include/deal.II/matrix_free/fe_evaluation.h | 10 +- tests/matrix_free/get_functions_notempl.cc | 115 ++++++ .../matrix_free/get_functions_notempl.output | 147 ++++++++ tests/matrix_free/matrix_vector_06_notempl.cc | 71 ++++ .../matrix_vector_06_notempl.output | 13 + tests/matrix_free/matrix_vector_common.h | 2 +- .../matrix_vector_stokes_notempl.cc | 347 ++++++++++++++++++ .../matrix_vector_stokes_notempl.output | 12 + 8 files changed, 711 insertions(+), 6 deletions(-) create mode 100644 tests/matrix_free/get_functions_notempl.cc create mode 100644 tests/matrix_free/get_functions_notempl.output create mode 100644 tests/matrix_free/matrix_vector_06_notempl.cc create mode 100644 tests/matrix_free/matrix_vector_06_notempl.output create mode 100644 tests/matrix_free/matrix_vector_stokes_notempl.cc create mode 100644 tests/matrix_free/matrix_vector_stokes_notempl.output diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 7a95e3c881..3ff41dcc65 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -7404,9 +7404,9 @@ FEEvaluation #ifdef DEBUG // print error message when the dimensions do not match. Propose a possible // fix - if (fe_degree != -1 && static_cast(fe_degree) != this->data->fe_degree + if ((fe_degree != -1 && static_cast(fe_degree) != this->data->fe_degree) || - n_q_points != this->data->n_q_points) + (fe_degree != -1 && static_n_q_points != this->data->n_q_points)) { std::string message = "-------------------------------------------------------\n"; @@ -7439,13 +7439,13 @@ FEEvaluation proposed_dof_comp = no; break; } - if (n_q_points == + if (static_n_q_points == this->mapping_info->mapping_data_gen[this->quad_no].n_q_points[this->active_quad_index]) proposed_quad_comp = this->quad_no; else for (unsigned int no=0; nomapping_info->mapping_data_gen.size(); ++no) if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index] - == n_q_points) + == static_n_q_points) { proposed_quad_comp = no; break; @@ -7507,7 +7507,7 @@ FEEvaluation message += " " + correct_pos; Assert (static_cast(fe_degree) == this->data->fe_degree && - n_q_points == this->data->n_q_points, + static_n_q_points == this->data->n_q_points, ExcMessage(message)); } if (fe_no != numbers::invalid_unsigned_int) diff --git a/tests/matrix_free/get_functions_notempl.cc b/tests/matrix_free/get_functions_notempl.cc new file mode 100644 index 0000000000..418cd4a27e --- /dev/null +++ b/tests/matrix_free/get_functions_notempl.cc @@ -0,0 +1,115 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 +// operations in getting the function values, the function gradients, and the +// function Laplacians on a hyperball mesh with different sizes in the number +// of degrees of freedom per cell and quadrature points per cell. This is the +// same test as get_functions_rect but without using a template parameter on +// the degree for FEEvaluation. + +#include "../tests.h" + +std::ofstream logfile("output"); + +#include "get_functions_common.h" + + +template +void sub_test (const DoFHandler &dof, + const ConstraintMatrix &constraints, + MatrixFree &mf_data, + Vector &solution) +{ + deallog << "Test with fe_degree " << fe_degree + << ", n_q_points_1d: " << (n_q_points_1d) << std::endl; + const QGauss<1> quad (n_q_points_1d); + MappingQ mapping (2); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_gradients | update_hessians; + mf_data.reinit (mapping, dof, constraints, quad, data); + MatrixFreeTest mf (mf_data, mapping); + mf.test_functions (solution); +} + + + +template +void test () +{ + typedef double number; + const SphericalManifold manifold; + Triangulation tria; + GridGenerator::hyper_ball (tria); + typename Triangulation::active_cell_iterator + cell = tria.begin_active (), + endc = tria.end(); + for (; cell!=endc; ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->at_boundary(f)) + cell->face(f)->set_all_manifold_ids(0); + tria.set_manifold (0, manifold); + + // refine first and last cell + tria.begin(tria.n_levels()-1)->set_refine_flag(); + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + tria.refine_global (1); + + FE_Q fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + constraints.close(); + + + // in the other functions, use do_test in + // get_functions_common, but here we have to + // manually choose non-rectangular tests. + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + //std::cout << "Number of cells: " << dof.get_triangulation().n_active_cells() + // << std::endl; + //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; + //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl; + + Vector solution (dof.n_dofs()); + + // create vector with random entries + for (unsigned int i=0; i mf_data; + if (fe_degree > 1) + sub_test (dof, constraints, mf_data, + solution); + sub_test (dof, constraints, mf_data, + solution); + sub_test (dof, constraints, mf_data, + solution); + if (dim == 2) + sub_test (dof, constraints, mf_data, + solution); +} diff --git a/tests/matrix_free/get_functions_notempl.output b/tests/matrix_free/get_functions_notempl.output new file mode 100644 index 0000000000..bcea3cad30 --- /dev/null +++ b/tests/matrix_free/get_functions_notempl.output @@ -0,0 +1,147 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Test with fe_degree 1, n_q_points_1d: 1 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Test with fe_degree 1, n_q_points_1d: 3 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Test with fe_degree 1, n_q_points_1d: 4 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(2) +DEAL:2d::Test with fe_degree 2, n_q_points_1d: 1 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Test with fe_degree 2, n_q_points_1d: 2 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Test with fe_degree 2, n_q_points_1d: 4 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Test with fe_degree 2, n_q_points_1d: 5 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(3) +DEAL:2d::Test with fe_degree 3, n_q_points_1d: 2 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Test with fe_degree 3, n_q_points_1d: 3 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Test with fe_degree 3, n_q_points_1d: 5 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Test with fe_degree 3, n_q_points_1d: 6 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(4) +DEAL:2d::Test with fe_degree 4, n_q_points_1d: 3 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Test with fe_degree 4, n_q_points_1d: 4 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Test with fe_degree 4, n_q_points_1d: 6 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Test with fe_degree 4, n_q_points_1d: 7 +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Test with fe_degree 1, n_q_points_1d: 1 +DEAL:3d::Error function values: 0 +DEAL:3d::Error function gradients: 0 +DEAL:3d::Error function Laplacians: 0 +DEAL:3d::Error function diagonal of Hessian: 0 +DEAL:3d::Error function Hessians: 0 +DEAL:3d:: +DEAL:3d::Test with fe_degree 1, n_q_points_1d: 3 +DEAL:3d::Error function values: 0 +DEAL:3d::Error function gradients: 0 +DEAL:3d::Error function Laplacians: 0 +DEAL:3d::Error function diagonal of Hessian: 0 +DEAL:3d::Error function Hessians: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Test with fe_degree 2, n_q_points_1d: 1 +DEAL:3d::Error function values: 0 +DEAL:3d::Error function gradients: 0 +DEAL:3d::Error function Laplacians: 0 +DEAL:3d::Error function diagonal of Hessian: 0 +DEAL:3d::Error function Hessians: 0 +DEAL:3d:: +DEAL:3d::Test with fe_degree 2, n_q_points_1d: 2 +DEAL:3d::Error function values: 0 +DEAL:3d::Error function gradients: 0 +DEAL:3d::Error function Laplacians: 0 +DEAL:3d::Error function diagonal of Hessian: 0 +DEAL:3d::Error function Hessians: 0 +DEAL:3d:: +DEAL:3d::Test with fe_degree 2, n_q_points_1d: 4 +DEAL:3d::Error function values: 0 +DEAL:3d::Error function gradients: 0 +DEAL:3d::Error function Laplacians: 0 +DEAL:3d::Error function diagonal of Hessian: 0 +DEAL:3d::Error function Hessians: 0 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_06_notempl.cc b/tests/matrix_free/matrix_vector_06_notempl.cc new file mode 100644 index 0000000000..1772e8db0b --- /dev/null +++ b/tests/matrix_free/matrix_vector_06_notempl.cc @@ -0,0 +1,71 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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. +// +// --------------------------------------------------------------------- + + + +// same as matrix_vector_06 but without setting a template parameter on the +// polynomial degree, which requires to read it as a run time parameter in a +// separate code path in FEEvaluation + +#include "../tests.h" +#include +#include "create_mesh.h" + +std::ofstream logfile("output"); + +#include "matrix_vector_common.h" + +template +void test () +{ + Triangulation tria; + create_mesh (tria); + tria.begin_active ()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + typename Triangulation::active_cell_iterator cell, endc; + cell = tria.begin_active (); + endc = tria.end(); + for (; cell!=endc; ++cell) + if (cell->center().norm()<0.5) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + tria.begin(tria.n_levels()-1)->set_refine_flag(); + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + tria.refine_global(1); + cell = tria.begin_active (); + for (unsigned int i=0; i<10-3*dim; ++i) + { + cell = tria.begin_active (); + endc = tria.end(); + unsigned int counter = 0; + for (; cell!=endc; ++cell, ++counter) + if (counter % (7-i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + FE_Q fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction(), + constraints); + constraints.close(); + + // do not enable templates in FEEvaluation + do_test (dof, constraints); +} diff --git a/tests/matrix_free/matrix_vector_06_notempl.output b/tests/matrix_free/matrix_vector_06_notempl.output new file mode 100644 index 0000000000..c8d1cd1e9f --- /dev/null +++ b/tests/matrix_free/matrix_vector_06_notempl.output @@ -0,0 +1,13 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(2) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_common.h b/tests/matrix_free/matrix_vector_common.h index 87c97352c3..1837813c78 100644 --- a/tests/matrix_free/matrix_vector_common.h +++ b/tests/matrix_free/matrix_vector_common.h @@ -61,7 +61,7 @@ void do_test (const DoFHandler &dof, //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl; - MappingQGeneric mapping(fe_degree); + MappingQGeneric mapping(dof.get_fe().degree); MatrixFree mf_data; { const QGauss<1> quad (n_q_points_1d); diff --git a/tests/matrix_free/matrix_vector_stokes_notempl.cc b/tests/matrix_free/matrix_vector_stokes_notempl.cc new file mode 100644 index 0000000000..0e17648af1 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_notempl.cc @@ -0,0 +1,347 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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. +// +// --------------------------------------------------------------------- + + + +// same as matrix_vector_stokes_noflux but no template parameter on the +// polynomial degree + +#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 +#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); + FEEvaluation pressure (data, 1); + + 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 unsigned int fe_degree) +{ + SphericalManifold manifold; + HyperShellBoundary boundary; + Triangulation triangulation; + GridGenerator::hyper_shell (triangulation, Point(), + 0.5, 1., 96, true); + triangulation.set_all_manifold_ids(0); + triangulation.set_all_manifold_ids_on_boundary(1); + triangulation.set_manifold (0, manifold); + triangulation.set_manifold (1, boundary); + + 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); + FE_Q fe_u_scal (fe_degree+1); + FESystem fe_u (fe_u_scal,dim); + FE_Q fe_p (fe_degree); + FESystem fe (fe_u_scal, dim, fe_p, 1); + DoFHandler dof_handler_u (triangulation); + DoFHandler dof_handler_p (triangulation); + DoFHandler dof_handler (triangulation); + + MatrixFree mf_data; + + ConstraintMatrix constraints, constraints_u, constraints_p; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + BlockVector mf_solution; + + dof_handler.distribute_dofs (fe); + dof_handler_u.distribute_dofs (fe_u); + dof_handler_p.distribute_dofs (fe_p); + 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_no_normal_flux_constraints (dof_handler, 0, + no_normal_flux_boundaries, + constraints, mapping); + constraints.close (); + DoFTools::make_hanging_node_constraints (dof_handler_u, + constraints_u); + VectorTools::compute_no_normal_flux_constraints (dof_handler_u, 0, + no_normal_flux_boundaries, + constraints_u, mapping); + constraints_u.close (); + DoFTools::make_hanging_node_constraints (dof_handler_p, + constraints_p); + constraints_p.close (); + + std::vector dofs_per_block (2); + DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, + stokes_sub_blocks); + + //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; + + { + 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 (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 (2); + for (unsigned int d=0; d<2; ++d) + solution.block(d).reinit (dofs_per_block[d]); + solution.collect_sizes (); + + system_rhs.reinit (solution); + mf_solution.reinit (solution); + + // fill system_rhs with random numbers + for (unsigned int j=0; j*> dofs; + dofs.push_back(&dof_handler_u); + dofs.push_back(&dof_handler_p); + std::vector constraints; + constraints.push_back (&constraints_u); + constraints.push_back (&constraints_p); + QGauss<1> quad(fe_degree+2); + // no parallelism + mf_data.reinit (mapping, dofs, constraints, quad, + typename MatrixFree::AdditionalData + (MPI_COMM_WORLD, + 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 degree " << fe_degree << ": " + << error/relative << std::endl << std::endl; +} + + + +int main () +{ + deallog.attach(logfile); + + deallog << std::setprecision (3); + + { + deallog << std::endl << "Test with doubles" << std::endl << std::endl; + deallog.threshold_double(1.e-12); + deallog.push("2d"); + test<2>(1); + 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_notempl.output b/tests/matrix_free/matrix_vector_stokes_notempl.output new file mode 100644 index 0000000000..0434059632 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_notempl.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:: -- 2.39.5