From 651053649030b203abe1152fdf1b667993dfb496 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Fri, 28 Feb 2014 16:05:32 +0000 Subject: [PATCH] Fix bug git-svn-id: https://svn.dealii.org/trunk@32590 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/matrix_free/fe_evaluation.h | 48 ++-- tests/matrix_free/matrix_vector_14b.cc | 245 ++++++++++++++++++ tests/matrix_free/matrix_vector_14b.output | 13 + 3 files changed, 291 insertions(+), 15 deletions(-) create mode 100644 tests/matrix_free/matrix_vector_14b.cc create mode 100644 tests/matrix_free/matrix_vector_14b.output diff --git a/deal.II/include/deal.II/matrix_free/fe_evaluation.h b/deal.II/include/deal.II/matrix_free/fe_evaluation.h index 05f544e929..3d7aed87b9 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -1753,10 +1753,10 @@ namespace internal template class FEEvaluationDGP : - public FEEvaluationGeneral + public FEEvaluationGeneral { public: - typedef FEEvaluationGeneral BaseClass; + typedef FEEvaluationGeneral BaseClass; typedef Number number_type; typedef typename BaseClass::value_type value_type; typedef typename BaseClass::gradient_type gradient_type; @@ -2397,6 +2397,8 @@ FEEvaluationBase mapped_geometry->get_cell()->level(), mapped_geometry->get_cell()->index(), dof_handler); + Assert (sizeof(src) == sizeof(VectorType *), + ExcMessage("Only one single vector allowed")); local_dof_indices.resize(dof_handler->get_fe().dofs_per_cell); cell->get_dof_indices(local_dof_indices); @@ -4382,7 +4384,7 @@ FEEvaluationGeneral // points unsigned int proposed_dof_comp = numbers::invalid_unsigned_int, proposed_quad_comp = numbers::invalid_unsigned_int; - if (dofs_per_cell == this->matrix_info->get_dof_info(fe_no).dofs_per_cell[this->active_fe_index]) + if (dofs_per_cell == this->data->dofs_per_cell) proposed_dof_comp = fe_no; else for (unsigned int no=0; nomatrix_info->n_components(); ++no) @@ -4476,6 +4478,14 @@ FEEvaluationGeneral BaseClass (geometry, dof_handler, first_selected_component) { set_data_pointers(); + + Assert ((dofs_per_cell == this->data->dofs_per_cell || + internal::MatrixFreeFunctions::DGP_dofs_per_cell::value == + this->data->dofs_per_cell + ) + && + n_q_points == this->data->n_q_points, + ExcMessage("Underlying element and template arguments do not match")); } @@ -6625,7 +6635,11 @@ FEEvaluationDGP const unsigned int quad_no) : BaseClass (data_in, fe_no, quad_no) -{} +{ + // reset values_dofs pointer as it has wider gaps in the base class + for (unsigned int c=0; cvalues_dofs[c] = &this->my_data_array[c*dofs_per_cell]; +} @@ -6638,7 +6652,10 @@ FEEvaluationDGP const unsigned int first_selected_component) : BaseClass (geometry, dof_handler, first_selected_component) -{} +{ + for (unsigned int c=0; cvalues_dofs[c] = &this->my_data_array[c*dofs_per_cell]; +} @@ -6649,7 +6666,10 @@ FEEvaluationDGP ::FEEvaluationDGP (const FEEvaluationDGP &other) : BaseClass (other) -{} +{ + for (unsigned int c=0; cvalues_dofs[c] = &this->my_data_array[c*dofs_per_cell]; +} @@ -6666,10 +6686,10 @@ FEEvaluationDGP internal::ExcAccessToUninitializedField()); // expand dof_values to tensor product - VectorizedArray data_array[n_components*Utilities::fixed_int_power::value]; + VectorizedArray data_array[n_components*BaseClass::dofs_per_cell]; VectorizedArray *expanded_dof_values[n_components]; for (unsigned int c=0; c::value]; + expanded_dof_values[c] = &data_array[c*BaseClass::dofs_per_cell]; unsigned int count_p = 0, count_q = 0; for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i) @@ -6677,15 +6697,13 @@ FEEvaluationDGP for (unsigned int j=0; j<(dim>1?fe_degree+1-i:1); ++j) { for (unsigned int k=0; kvalues_dofs[c][count_p]; - } + for (unsigned int c=0; cvalues_dofs[c][count_p]; for (unsigned int k=fe_degree+1-j-i; k(); } - for (unsigned int j=(dim>1?fe_degree+1-i:1); j<(dim>1?fe_degree+1:1); ++j) + for (unsigned int j=fe_degree+1-i; j(); @@ -6724,10 +6742,10 @@ FEEvaluationDGP Assert (this->gradients_quad_submitted == true, internal::ExcAccessToUninitializedField()); - VectorizedArray data_array[n_components*Utilities::fixed_int_power::value]; + VectorizedArray data_array[n_components*BaseClass::dofs_per_cell]; VectorizedArray *expanded_dof_values[n_components]; for (unsigned int c=0; c::value]; + expanded_dof_values[c] = &data_array[c*BaseClass::dofs_per_cell]; internal::do_integrate (*this, expanded_dof_values, this->values_quad, this->gradients_quad, integrate_val, integrate_grad, internal::int2type()); diff --git a/tests/matrix_free/matrix_vector_14b.cc b/tests/matrix_free/matrix_vector_14b.cc new file mode 100644 index 0000000000..04e33dd37f --- /dev/null +++ b/tests/matrix_free/matrix_vector_14b.cc @@ -0,0 +1,245 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2014 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 using more quadrature points than the FE +// degree would suggest + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +std::ofstream logfile("output"); + + + +template +void +helmholtz_operator_dgp (const MatrixFree &data, + VECTOR &dst, + const VECTOR &src, + const std::pair &cell_range) +{ + typedef typename VECTOR::value_type Number; + FEEvaluationDGP fe_eval (data); + const unsigned int n_q_points = fe_eval.n_q_points; + + for (unsigned int cell=cell_range.first; cell > +class MatrixFreeTest +{ +public: + typedef VectorizedArray vector_t; + + MatrixFreeTest(const MatrixFree &data_in): + data (data_in) + {}; + + void vmult (VECTOR &dst, + const VECTOR &src) const + { + dst = 0; + const std_cxx1x::function &, + VECTOR &, + const VECTOR &, + const std::pair &)> + wrap = helmholtz_operator_dgp; + data.cell_loop (wrap, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void do_test (const DoFHandler &dof, + const ConstraintMatrix &constraints) +{ + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + MatrixFree mf_data; + { + const QGauss<1> quad (fe_degree+2); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + data.tasks_block_size = 3; + + mf_data.reinit (dof, constraints, quad, data); + } + + MatrixFreeTest mf (mf_data); + Vector in (dof.n_dofs()), out (dof.n_dofs()); + Vector in_dist (dof.n_dofs()); + Vector out_dist (in_dist); + + for (unsigned int i=0; i sparse_matrix (sparsity); + { + QGauss quadrature_formula(fe_degree+2); + + FEValues fe_values (dof.get_fe(), quadrature_formula, + update_values | update_gradients | + update_JxW_values); + + const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + fe_values.reinit (cell); + + for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, + local_dof_indices, + sparse_matrix); + } + } + + sparse_matrix.vmult (out, in); + out -= out_dist; + const double diff_norm = out.linfty_norm() / out_dist.linfty_norm(); + + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} + + + +template +void test () +{ + Triangulation tria; + GridGenerator::hyper_ball (tria); + static const HyperBallBoundary boundary; + tria.set_boundary (0, boundary); + 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(); + typename Triangulation::active_cell_iterator + cell = tria.begin_active (), + endc = tria.end(); + 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); +} + + + +int main () +{ + deallog.attach(logfile); + deallog.depth_console(0); + + deallog << std::setprecision (3); + + { + deallog.threshold_double(5.e-11); + deallog.push("2d"); + test<2,1>(); + test<2,2>(); + deallog.pop(); + deallog.push("3d"); + test<3,1>(); + test<3,2>(); + deallog.pop(); + } +} + + diff --git a/tests/matrix_free/matrix_vector_14b.output b/tests/matrix_free/matrix_vector_14b.output new file mode 100644 index 0000000000..2a6f032686 --- /dev/null +++ b/tests/matrix_free/matrix_vector_14b.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:: -- 2.39.5