From: Martin Kronbichler Date: Tue, 28 Feb 2017 20:59:50 +0000 (+0100) Subject: Additional tests without the template X-Git-Tag: v8.5.0-rc1~87^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F4002%2Fhead;p=dealii.git Additional tests without the template --- diff --git a/tests/matrix_free/matrix_vector_14_notempl.cc b/tests/matrix_free/matrix_vector_14_notempl.cc new file mode 100644 index 0000000000..a1952df11f --- /dev/null +++ b/tests/matrix_free/matrix_vector_14_notempl.cc @@ -0,0 +1,62 @@ +// --------------------------------------------------------------------- +// +// 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_14 but without using the template on the element +// degree + +#include "../tests.h" +#include + + +std::ofstream logfile("output"); + +#include "matrix_vector_common.h" + + + +template +void test () +{ + 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); + if (dim < 3 || fe_degree < 2) + tria.refine_global(1); + tria.begin(tria.n_levels()-1)->set_refine_flag(); + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + cell = tria.begin_active (); + for (; cell!=endc; ++cell) + if (cell->center().norm()<1e-8) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + FE_DGP fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + ConstraintMatrix constraints; + + do_test (dof, constraints); +} diff --git a/tests/matrix_free/matrix_vector_14_notempl.output b/tests/matrix_free/matrix_vector_14_notempl.output new file mode 100644 index 0000000000..2a6f032686 --- /dev/null +++ b/tests/matrix_free/matrix_vector_14_notempl.output @@ -0,0 +1,13 @@ + +DEAL:2d::Testing FE_DGP<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FE_DGP<2>(2) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FE_DGP<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_DGP<3>(2) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_stokes_qdg0b.cc b/tests/matrix_free/matrix_vector_stokes_qdg0b.cc new file mode 100644 index 0000000000..ec5facf373 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_qdg0b.cc @@ -0,0 +1,331 @@ +// --------------------------------------------------------------------- +// +// 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 +// matrix-vector products by comparing with the result of deal.II sparse +// matrix. No hanging nodes and no other constraints for a vector-valued +// problem (stokes equations). Same as matrix_vector_stokes_qdg0 but now +// without the template on the degree of the element. + +#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 "create_mesh.h" + +#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 = -velocity.get_divergence(q); + pressure.submit_value (div, q); + + // subtract p * I + for (unsigned int d=0; d &data; +}; + + + +template +void test (const unsigned int fe_degree) +{ + Triangulation triangulation; + create_mesh (triangulation); + if (fe_degree == 1) + triangulation.refine_global (4-dim); + else + triangulation.refine_global (3-dim); + + FE_Q fe_u (fe_degree+1); + FE_Q_DG0 fe_p (fe_degree); + FESystem fe (fe_u, dim, fe_p, 1); + DoFHandler dof_handler_u (triangulation); + DoFHandler dof_handler_p (triangulation); + DoFHandler dof_handler (triangulation); + + MatrixFree mf_data; + + ConstraintMatrix constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + std::vector > vec1, vec2; + + dof_handler.distribute_dofs (fe); + dof_handler_u.distribute_dofs (fe_u); + dof_handler_p.distribute_dofs (fe_p); + DoFRenumbering::component_wise (dof_handler); + + constraints.close (); + + std::vector dofs_per_block (dim+1); + DoFTools::count_dofs_per_component (dof_handler, dofs_per_block); + + //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 (dim+1,dim+1); + + for (unsigned int d=0; d 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); + } + } + + // first system_rhs with random numbers + for (unsigned int i=0; i*> dofs; + dofs.push_back(&dof_handler_u); + dofs.push_back(&dof_handler_p); + ConstraintMatrix dummy_constraints; + dummy_constraints.close(); + std::vector constraints; + constraints.push_back (&dummy_constraints); + constraints.push_back (&dummy_constraints); + QGauss<1> quad(fe_degree+2); + mf_data.reinit (dofs, constraints, quad, + typename MatrixFree::AdditionalData + (MPI_COMM_WORLD, + MatrixFree::AdditionalData::none)); + } + + system_matrix.vmult (solution, system_rhs); + + typedef std::vector > VectorType; + MatrixFreeTest mf (mf_data); + mf.vmult (vec2, vec1); + + // Verification + double error = 0.; + for (unsigned int i=0; i(1); + test<2>(2); + test<2>(3); + deallog.pop(); + deallog.push("3d"); + test<3>(1); + test<3>(2); + deallog.pop(); + } +} diff --git a/tests/matrix_free/matrix_vector_stokes_qdg0b.output b/tests/matrix_free/matrix_vector_stokes_qdg0b.output new file mode 100644 index 0000000000..db6c7de2a9 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_qdg0b.output @@ -0,0 +1,14 @@ + +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:: +DEAL:3d:: Verification fe degree 2: 0 +DEAL:3d::