From a28f720d1a8432305522f2dd6b37a772b67f6a61 Mon Sep 17 00:00:00 2001 From: Katharina Kormann Date: Sun, 24 Jun 2012 14:36:54 +0000 Subject: [PATCH] Fix bug pointed out by Saul Teukolsky and add submit_divergence function. git-svn-id: https://svn.dealii.org/trunk@25639 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/matrix_free/fe_evaluation.h | 144 ++++++--- .../include/deal.II/matrix_free/matrix_free.h | 2 +- tests/matrix_free/matrix_vector_div.cc | 289 ++++++++++++++++++ .../matrix_free/matrix_vector_div/cmp/generic | 12 + 4 files changed, 410 insertions(+), 37 deletions(-) create mode 100644 tests/matrix_free/matrix_vector_div.cc create mode 100644 tests/matrix_free/matrix_vector_div/cmp/generic 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 580042ef36..f68e429ea4 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -429,14 +429,16 @@ public: gradient_type get_gradient (const unsigned int q_point) const; /** - * Write a gradient to the field containing - * the values on quadrature points with - * component @p q_point. Access to the same - * field as through @p get_gradient. If + * Write a contribution that is tested + * by the gradient to the field + * containing the values on quadrature + * points with component @p + * q_point. Access to the same field + * as through @p get_gradient. If * applied before the function @p - * integrate(...,true) is called, - * this specifies the gradient which is tested - * by all basis function gradients on the + * integrate(...,true) is called, this + * specifies what is tested by all + * basis function gradients on the * current cell and integrated over. * * Note that the derived class @@ -1081,17 +1083,17 @@ class FEEvaluationAccess : gradient_type get_gradient (const unsigned int q_point) const; /** - * Write a gradient to the field + * Write a contribution that is tested + * by the gradient to the field * containing the values on quadrature * points with component @p * q_point. Access to the same field * as through @p get_gradient. If * applied before the function @p * integrate(...,true) is called, this - * specifies the gradient which is - * tested by all basis function - * gradients on the current cell and - * integrated over. + * specifies what is tested by all + * basis function gradients on the + * current cell and integrated over. */ void submit_gradient(const gradient_type grad_in, const unsigned int q_point); @@ -1241,45 +1243,70 @@ class FEEvaluationAccess : gradient_type get_hessian_diagonal (const unsigned int q_point) const; /** - * Write a gradient to the field containing - * the values on quadrature points with - * component @p q_point. Access to the same - * field as through @p get_gradient. If + * Write a contribution that is tested + * by the gradient to the field + * containing the values on quadrature + * points with component @p + * q_point. Access to the same field + * as through @p get_gradient. If * applied before the function @p - * integrate(...,true) is called, - * this specifies the gradient which is tested - * by all basis function gradients on the + * integrate(...,true) is called, this + * specifies what is tested by all + * basis function gradients on the * current cell and integrated over. */ void submit_gradient(const gradient_type grad_in, const unsigned int q_point); /** - * Write a gradient to the field containing - * the values on quadrature points with - * component @p q_point. This function is an - * alternative to the other submit_gradient - * function when using a system of fixed - * number of equations which happens to - * coincide with the dimension for some - * dimensions, but not all. To allow for - * dimension-independent programming, this - * function can be used instead. + * Write a contribution that is tested + * by the gradient to the field + * containing the values on quadrature + * points with component @p + * q_point. This function is an + * alternative to the other + * submit_gradient function when using + * a system of fixed number of + * equations which happens to coincide + * with the dimension for some + * dimensions, but not all. To allow + * for dimension-independent + * programming, this function can be + * used instead. */ void submit_gradient(const Tensor<1,dim,Tensor<1,dim,VectorizedArray > > grad_in, const unsigned int q_point); /** - * Write a gradient to the field containing - * the values on quadrature points with - * component @p q_point. Access to the same - * field as through @p get_gradient. If + * Write a constribution that is + * tested by the divergence to the field + * containing the values on quadrature + * points with component @p + * q_point. Access to the same field + * as through @p get_gradient. If * applied before the function @p - * integrate(...,true) is called, - * this specifies the gradient which is tested - * by all basis function gradients on the + * integrate(...,true) is called, this + * specifies what is tested by all + * basis function gradients on the * current cell and integrated over. */ + void submit_divergence (const VectorizedArray div_in, + const unsigned int q_point); + + /** + * Write a contribution that is tested + * by the gradient to the field + * containing the values on quadrature + * points with component @p + * q_point. Access to the same field + * as through @p get_gradient. If + * applied before the function @p + * integrate(...,true) is called, this + * specifies the gradient which is + * tested by all basis function + * gradients on the current cell and + * integrated over. + */ void submit_symmetric_gradient(const SymmetricTensor<2,dim,VectorizedArray > grad_in, const unsigned int q_point); @@ -3878,6 +3905,51 @@ FEEvaluationAccess BaseClass::submit_gradient(grad_in, q_point); } +template +inline +void +FEEvaluationAccess +::submit_divergence (const VectorizedArray div_in, + const unsigned int q_point) +{ +#ifdef DEBUG + Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); + AssertIndexRange (q_point, n_q_points); + this->gradients_quad_submitted = true; +#endif + if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + const VectorizedArray fac = this->J_value[0] * + this->quadrature_weights[q_point] * div_in; + for (unsigned int d=0; dgradients_quad[d][d][q_point] = (fac * + this->cartesian_data[0][d]); + for(unsigned int e=d+1;egradients_quad[d][e][q_point] = VectorizedArray(); + this->gradients_quad[e][d][q_point] = VectorizedArray(); + } + } + } + else + { + const Tensor<2,dim,VectorizedArray > &jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->jacobian[q_point] : this->jacobian[0]; + const VectorizedArray fac = + (this->cell_type == internal::MatrixFreeFunctions::general ? + this->J_value[q_point] : this->J_value[0] * + this->quadrature_weights[q_point]) * div_in; + for (unsigned int d=0; dgradients_quad[d][e][q_point] = jac[d][e] * fac; + } + } +} + template diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.h b/deal.II/include/deal.II/matrix_free/matrix_free.h index bacc960834..864307f03f 100644 --- a/deal.II/include/deal.II/matrix_free/matrix_free.h +++ b/deal.II/include/deal.II/matrix_free/matrix_free.h @@ -1218,7 +1218,7 @@ MatrixFree::initialize_dof_vector(VectorType &vec, const unsigned int comp) const { AssertIndexRange (comp, n_components()); - vec.reinit(dof_info[comp].vector_partitioner->global_size); + vec.reinit(dof_info[comp].vector_partitioner->size()); } diff --git a/tests/matrix_free/matrix_vector_div.cc b/tests/matrix_free/matrix_vector_div.cc new file mode 100644 index 0000000000..4ecf10ca7d --- /dev/null +++ b/tests/matrix_free/matrix_vector_div.cc @@ -0,0 +1,289 @@ +//------------------ matrix_vector_div.cc ------------------------ +// $Id$ +// Version: $Name$ +// +//------------------ matrix_vector_div.cc ------------------------ + + +// this function tests the correctness of the implementation of matrix free +// matrix-vector products by comparing with the result of deal.II sparse +// matrix. The mesh uses a hyperball mesh with hanging nodes for a +// vector-valued problem (div-div operator which does not really make a lot +// of sense from a problem point of view, though). + +#include "../tests.h" + +std::ofstream logfile("matrix_vector_div/output"); + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "create_mesh.h" + +const double global_coefficient = 0.1; + + +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 phi (data); + vector_t coeff = make_vectorized_array(global_coefficient); + + for(unsigned int cell=cell_range.first;cell::local_apply, + this, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void test () +{ + Triangulation tria; + create_mesh (tria); + tria.refine_global(4-dim); + + // refine a few cells + for (unsigned int i=0; i<10-3*dim; ++i) + { + typename Triangulation::active_cell_iterator + 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_sca (QGaussLobatto<1>(fe_degree+1)); + FESystem fe (fe_sca, dim); + DoFHandler dof_handler_sca (tria); + DoFHandler dof_handler (tria); + + 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_sca.distribute_dofs (fe_sca); + DoFRenumbering::component_wise (dof_handler); + + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + const unsigned int dofs_per_block = dof_handler_sca.n_dofs(); + { + BlockCompressedSimpleSparsityPattern csp (dim,dim); + for (unsigned int d=0; d quadrature_formula(fe_degree+1); + + 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 sc (0); + + std::vector phi_div (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; q phi_grad = fe_values[sc].gradient(k,q); + phi_div[k] = trace(phi_grad); + } + + for (unsigned int i=0; iget_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 quad(fe_degree+1); + mf_data.reinit (dof_handler_sca, 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(); + test<2,2>(); + deallog.pop(); + deallog.push("3d"); + test<3,1>(); + test<3,2>(); + deallog.pop(); + } +} + diff --git a/tests/matrix_free/matrix_vector_div/cmp/generic b/tests/matrix_free/matrix_vector_div/cmp/generic new file mode 100644 index 0000000000..619c272047 --- /dev/null +++ b/tests/matrix_free/matrix_vector_div/cmp/generic @@ -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:3d:: Verification fe degree 1: 0 +DEAL:3d:: +DEAL:3d:: Verification fe degree 2: 0 +DEAL:3d:: -- 2.39.5