From 5ae0d4c9e45cf1a114f73f4a67babcee9b1527c7 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 5 Nov 2014 14:02:16 +0100 Subject: [PATCH] Further actions to make MatrixFree/FEEvaluation work in 1D - Added test for matrix-vector product - Made compile with clang. --- include/deal.II/matrix_free/fe_evaluation.h | 14 +- .../matrix_free/mapping_info.templates.h | 2 +- tests/matrix_free/matrix_vector_18.cc | 261 ++++++++++++++++++ tests/matrix_free/matrix_vector_18.output | 13 + 4 files changed, 280 insertions(+), 10 deletions(-) create mode 100644 tests/matrix_free/matrix_vector_18.cc create mode 100644 tests/matrix_free/matrix_vector_18.output diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index c143e2f5ad..9dfa0eabcf 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -1236,12 +1236,11 @@ template class FEEvaluationAccess<1,1,Number> : public FEEvaluationBase<1,1,Number> { public: - static const int dim = 1; typedef Number number_type; typedef VectorizedArray value_type; - typedef Tensor<1,dim,VectorizedArray > gradient_type; - static const unsigned int dimension = dim; - typedef FEEvaluationBase BaseClass; + typedef Tensor<1,1,VectorizedArray > gradient_type; + static const unsigned int dimension = 1; + typedef FEEvaluationBase<1,1,Number> BaseClass; /** * Returns the value stored for the local degree of freedom with index @p @@ -4265,7 +4264,7 @@ FEEvaluationAccess<1,1,Number> template inline FEEvaluationAccess<1,1,Number> -::FEEvaluationAccess (const FEEvaluationAccess &other) +::FEEvaluationAccess (const FEEvaluationAccess<1,1,Number> &other) : FEEvaluationBase <1,1,Number>(other) {} @@ -4437,10 +4436,7 @@ FEEvaluationAccess<1,1,Number> this->cell_type == internal::MatrixFreeFunctions::general ? this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point]; - VectorizedArray new_val = jac[0][0] * grad_in[0]; - for (unsigned int e=1; egradients_quad[0][0][q_point] = new_val * JxW; + this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in[0] * JxW; } } diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index b0c3826589..d2f74ad086 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -191,7 +191,7 @@ namespace internal current_data.n_q_points.push_back (n_q_points); current_data.n_q_points_face.push_back - (Utilities::fixed_power(n_q_points_1d[q])); + (dim>1 ? Utilities::fixed_power(n_q_points_1d[q]) : 1); current_data.quadrature.push_back (Quadrature(quad[my_q][q])); current_data.face_quadrature.push_back diff --git a/tests/matrix_free/matrix_vector_18.cc b/tests/matrix_free/matrix_vector_18.cc new file mode 100644 index 0000000000..f3c5e596c0 --- /dev/null +++ b/tests/matrix_free/matrix_vector_18.cc @@ -0,0 +1,261 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +// this tests matrix-free matrix-vector products in 1D (otherwise similar as +// the other matrix-vector tests) + + +#include "../tests.h" + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + +template > +class MatrixFreeTest +{ +public: + MatrixFreeTest(const MatrixFree &data_in): + data (data_in) + {}; + + void vmult (VECTOR &dst, + const VECTOR &src) const + { + dst = 0; + data.cell_loop (&MatrixFreeTest::local_operation, + this, dst, src); + }; + +private: + const MatrixFree &data; + + void local_operation (const MatrixFree &data, + VECTOR &out, + const VECTOR &in, + const std::pair &cell_range) const + { + FEEvaluation fe_eval (data); + const unsigned int n_q_points = fe_eval.n_q_points; + + Tensor<1,dim,VectorizedArray > ones; + for (unsigned int d=0; d(3.2221)* + (fe_eval.get_hessian(q)*ones),q); + } + fe_eval.integrate (true,true); + fe_eval.distribute_local_to_global (out); + } + } +}; + + + +template +void do_test (const DoFHandler &dof, + const ConstraintMatrix &constraints, + const unsigned int parallel_option = 0) +{ + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + if (parallel_option > 0) + deallog << "Parallel option: " << parallel_option << std::endl; + + MatrixFree mf_data; + { + const QGauss<1> quad (fe_degree+1); + typename MatrixFree::AdditionalData data; + if (parallel_option == 1) + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + else if (parallel_option == 2) + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::color; + else + { + Assert (parallel_option == 0, ExcInternalError()); + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_partition; + } + data.tasks_block_size = 7; + data.mapping_update_flags |= update_second_derivatives; + + 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+1); + + FEValues fe_values (dof.get_fe(), quadrature_formula, + update_values | update_gradients | + update_JxW_values | update_second_derivatives); + + 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); + + Tensor<1,dim> ones; + for (unsigned int d=0; d::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_cube (tria); + 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_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_test (dof, constraints); +} + + + +int main () +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + + deallog << std::setprecision (3); + + { + deallog.threshold_double(1.e-9); + deallog.push("1d"); + test<1,1>(); + test<1,2>(); + test<1,3>(); + deallog.pop(); + deallog.push("2d"); + test<2,1>(); + deallog.pop(); + } +} + diff --git a/tests/matrix_free/matrix_vector_18.output b/tests/matrix_free/matrix_vector_18.output new file mode 100644 index 0000000000..35142b13b8 --- /dev/null +++ b/tests/matrix_free/matrix_vector_18.output @@ -0,0 +1,13 @@ + +DEAL:1d::Testing FE_Q<1>(1) +DEAL:1d::Norm of difference: 0 +DEAL:1d:: +DEAL:1d::Testing FE_Q<1>(2) +DEAL:1d::Norm of difference: 0 +DEAL:1d:: +DEAL:1d::Testing FE_Q<1>(3) +DEAL:1d::Norm of difference: 0 +DEAL:1d:: +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: -- 2.39.5