From 61d9f0b3d9796cf09c9af0498ebda439cfc2fea3 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 13 Jun 2014 12:01:42 +0000 Subject: [PATCH] Implement FEEvaluation for FE_Q_DG0 elements git-svn-id: https://svn.dealii.org/trunk@33040 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 6 + deal.II/include/deal.II/base/aligned_vector.h | 1 + .../include/deal.II/fe/fe_poly.templates.h | 17 + .../deal.II/matrix_free/fe_evaluation.h | 457 ++++++++++++++---- .../include/deal.II/matrix_free/shape_info.h | 79 ++- .../matrix_free/shape_info.templates.h | 12 +- tests/matrix_free/matrix_vector_17.cc | 248 ++++++++++ tests/matrix_free/matrix_vector_17.output | 13 + .../matrix_free/matrix_vector_stokes_qdg0.cc | 333 +++++++++++++ .../matrix_vector_stokes_qdg0.output | 14 + 10 files changed, 1046 insertions(+), 134 deletions(-) create mode 100644 tests/matrix_free/matrix_vector_17.cc create mode 100644 tests/matrix_free/matrix_vector_17.output create mode 100644 tests/matrix_free/matrix_vector_stokes_qdg0.cc create mode 100644 tests/matrix_free/matrix_vector_stokes_qdg0.output diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 84a45349c6..de4f3e3f36 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -148,6 +148,12 @@ inconvenience this causes.

Specific improvements

    +
  1. New: There is now a class FEEvaluationQ_DG0 that does + optimized matrix-free evaluation for FE_Q_DG0 elements. +
    + (Martin Kronbichler, 2014/06/13) +
    2. +
  2. Bugfix: Filter libclang_rt* from the PETSc link line.
    (Matthias Maier, 2014/06/04) diff --git a/deal.II/include/deal.II/base/aligned_vector.h b/deal.II/include/deal.II/base/aligned_vector.h index 88576b5acb..f963e8e894 100644 --- a/deal.II/include/deal.II/base/aligned_vector.h +++ b/deal.II/include/deal.II/base/aligned_vector.h @@ -487,6 +487,7 @@ AlignedVector::resize_fast (const size_type size) } + template < class T > inline void diff --git a/deal.II/include/deal.II/fe/fe_poly.templates.h b/deal.II/include/deal.II/fe/fe_poly.templates.h index bfcc19d3ac..6bdbed9400 100644 --- a/deal.II/include/deal.II/fe/fe_poly.templates.h +++ b/deal.II/include/deal.II/fe/fe_poly.templates.h @@ -18,6 +18,7 @@ #include #include #include +#include #include #include @@ -480,6 +481,22 @@ namespace internal { return poly.get_numbering_inverse(); } + + template + inline + std::vector + get_poly_space_numbering (const TensorProductPolynomialsConst &poly) + { + return poly.get_numbering(); + } + + template + inline + std::vector + get_poly_space_numbering_inverse (const TensorProductPolynomialsConst &poly) + { + return poly.get_numbering_inverse(); + } } 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 02295fe322..c00b23b0b9 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -1394,13 +1394,25 @@ protected: /** * Internally stored variables for the different data fields. */ - VectorizedArray my_data_array[n_components*(dofs_per_cell+(dim*dim+2*dim+1)*n_q_points)]; + VectorizedArray my_data_array[n_components*(dofs_per_cell+1+(dim*dim+2*dim+1)*n_q_points)]; + + /** + * Checks the number of cell dofs with the value that we expect + */ + void check_dofs_per_cell (const unsigned int dofs_per_cell) const; private: /** * Sets the pointers from the data array to values_dof, etc. */ void set_data_pointers(); + + /** + * Checks that the element number of degrees of freedom given by the + * template arguments (via fe_degree) coincides with the number of degrees + * of freedom in the stored shape values. + */ + void check_template_arguments(const unsigned int fe_no); }; @@ -1576,7 +1588,7 @@ protected: private: /** - * Fills the fields shapve_???_evenodd, called in the constructor. + * Fills the fields shape_???_evenodd, called in the constructor. */ void compute_even_odd_factors(); }; @@ -1828,6 +1840,98 @@ public: }; + +/** + * The class that provides all functions necessary to evaluate functions at + * quadrature points and cell integrations. In functionality, this class is + * similar to FEValues, however, it includes a lot of specialized + * functions that make it much faster (between 5 and 500 times as fast, + * depending on the polynomial order). Access to the data fields is provided + * through functionality in the class FEEvaluationAccess. + * + * This class is an extension of FEEvaluation to work with continuous elements + * supplemented with a single discontinuous degree of freedom, i.e., FE_Q_DG0. + * + * @author Martin Kronbichler, 2014 + */ +template +class FEEvaluationQ_DG0 : + public FEEvaluation +{ +public: + typedef FEEvaluation BaseClass; + typedef Number number_type; + typedef typename BaseClass::value_type value_type; + typedef typename BaseClass::gradient_type gradient_type; + static const unsigned int dimension = dim; + static const unsigned int n_components = n_components_; + static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell + 1; + static const unsigned int n_q_points = BaseClass::n_q_points; + + /** + * Constructor. Takes all data stored in MatrixFree. If applied to problems + * with more than one finite element or more than one quadrature formula + * selected during construction of @p matrix_free, @p fe_no and @p quad_no + * allow to select the appropriate components. + */ + FEEvaluationQ_DG0 (const MatrixFree &matrix_free, + const unsigned int fe_no = 0, + const unsigned int quad_no = 0); + + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. The user has to provide a structure of type MappingFEEvaluation + * and a DoFHandler in order to allow for reading out the finite element + * data. It uses the data provided by dof_handler.get_fe(). If the element + * is vector-valued, the optional argument allows to specify the index of + * the base element (as long as the element is primitive, non-primitive are + * not supported currently). + * + * With initialization from a FEValues object, no call to a reinit method of + * this class is necessary. Instead, it is enough if the geometry is + * initialized to a given cell iterator. It can also read from or write to + * vectors in the standard way for DoFHandler::active_cell_iterator + * types (which is less efficient with MPI since index translation has to be + * done), but of course only for one cell at a time. Hence, a kernel using + * this method does not vectorize over several elements (which is most + * efficient for vector operations), but only possibly within the element if + * the evaluate/integrate routines are combined (e.g. for computing cell + * matrices). + */ + FEEvaluationQ_DG0 (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component = 0); + + /** + * Copy constructor + */ + FEEvaluationQ_DG0 (const FEEvaluationQ_DG0 &other); + + /** + * Evaluates the function values, the gradients, and the Hessians of the FE + * function given at the DoF values in the input vector at the quadrature + * points of the unit cell. The function arguments specify which parts shall + * actually be computed. Needs to be called before the functions @p + * get_value(), @p get_gradient() or @p get_laplacian give useful + * information (unless these values have been set manually). + */ + void evaluate (const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl = false); + + /** + * This function takes the values and/or gradients that are stored on + * quadrature points, tests them by all the basis functions/gradients on the + * cell and performs the cell integration. The two function arguments @p + * integrate_val and @p integrate_grad are used to enable/disable some of + * values or gradients. + */ + void integrate (const bool integrate_val, + const bool integrate_grad); +}; + + /*----------------------- Inline functions ----------------------------------*/ #ifndef DOXYGEN @@ -4355,15 +4459,53 @@ FEEvaluationGeneral const unsigned int quad_no) : BaseClass (data_in, fe_no, quad_no, dofs_per_cell, n_q_points) +{ + check_template_arguments(fe_no); + set_data_pointers(); +} + + + +template +inline +FEEvaluationGeneral +::FEEvaluationGeneral (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component) + : + BaseClass (geometry, dof_handler, first_selected_component) +{ + check_template_arguments(numbers::invalid_unsigned_int); + set_data_pointers(); +} + + + +template +inline +FEEvaluationGeneral +::FEEvaluationGeneral (const FEEvaluationGeneral &other) + : + BaseClass (other) { set_data_pointers(); +} + + +template +inline +void +FEEvaluationGeneral +::check_template_arguments(const unsigned int fe_no) +{ #ifdef DEBUG // print error message when the dimensions do not match. Propose a possible // fix - if ((dofs_per_cell != this->data->dofs_per_cell && - internal::MatrixFreeFunctions::DGP_dofs_per_cell::value != - this->data->dofs_per_cell) + if (fe_degree != this->data->fe_degree || n_q_points != this->data->n_q_points) { @@ -4376,33 +4518,36 @@ FEEvaluationGeneral message += "," + Utilities::int_to_string(n_components); message += ",Number>(data, "; message += Utilities::int_to_string(fe_no) + ", "; - message += Utilities::int_to_string(quad_no) + ")\n"; + message += Utilities::int_to_string(this->quad_no) + ")\n"; // check whether some other vector component has the correct number of // points unsigned int proposed_dof_comp = numbers::invalid_unsigned_int, proposed_quad_comp = numbers::invalid_unsigned_int; - 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) - if (this->matrix_info->get_dof_info(no).dofs_per_cell[this->active_fe_index] - == dofs_per_cell) - { - proposed_dof_comp = no; - break; - } - if (n_q_points == - this->mapping_info->mapping_data_gen[quad_no].n_q_points[this->active_quad_index]) - proposed_quad_comp = 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) - { - proposed_quad_comp = no; - break; - } + if (fe_no != numbers::invalid_unsigned_int) + { + if (fe_degree == this->data->fe_degree) + proposed_dof_comp = fe_no; + else + for (unsigned int no=0; nomatrix_info->n_components(); ++no) + if (this->matrix_info->get_shape_info(no,0,this->active_fe_index,0).fe_degree + == fe_degree) + { + proposed_dof_comp = no; + break; + } + if (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) + { + proposed_quad_comp = no; + break; + } + } if (proposed_dof_comp != numbers::invalid_unsigned_int && proposed_quad_comp != numbers::invalid_unsigned_int) { @@ -4422,7 +4567,7 @@ FEEvaluationGeneral correct_pos = " ^ "; else correct_pos = " "; - if (proposed_quad_comp != quad_no) + if (proposed_quad_comp != this->quad_no) correct_pos += " ^\n"; else correct_pos += " \n"; @@ -4430,18 +4575,17 @@ FEEvaluationGeneral } // ok, did not find the numbers specified by the template arguments in // the given list. Suggest correct template arguments - const unsigned int proposed_fe_degree = static_cast(std::pow(1.001*this->data->dofs_per_cell,1./dim))-1; const unsigned int proposed_n_q_points_1d = static_cast(std::pow(1.001*this->data->n_q_points,1./dim)); message += "Wrong template arguments:\n"; message += " Did you mean FEEvaluationdata->fe_degree) + ","; message += Utilities::int_to_string(proposed_n_q_points_1d); message += "," + Utilities::int_to_string(n_components); message += ",Number>(data, "; message += Utilities::int_to_string(fe_no) + ", "; - message += Utilities::int_to_string(quad_no) + ")?\n"; + message += Utilities::int_to_string(this->quad_no) + ")?\n"; std::string correct_pos; - if (proposed_fe_degree != fe_degree) + if (this->data->fe_degree != fe_degree) correct_pos = " ^"; else correct_pos = " "; @@ -4451,55 +4595,23 @@ FEEvaluationGeneral correct_pos += " \n"; message += " " + correct_pos; - Assert (dofs_per_cell == this->data->dofs_per_cell && + Assert (fe_degree == this->data->fe_degree && n_q_points == this->data->n_q_points, ExcMessage(message)); } - AssertDimension (n_q_points, - this->mapping_info->mapping_data_gen[this->quad_no]. - n_q_points[this->active_quad_index]); - AssertDimension (this->data->dofs_per_cell * this->n_fe_components, - this->dof_info->dofs_per_cell[this->active_fe_index]); + if (fe_no != numbers::invalid_unsigned_int) + { + AssertDimension (n_q_points, + this->mapping_info->mapping_data_gen[this->quad_no]. + n_q_points[this->active_quad_index]); + AssertDimension (this->data->dofs_per_cell * this->n_fe_components, + this->dof_info->dofs_per_cell[this->active_fe_index]); + } #endif } -template -inline -FEEvaluationGeneral -::FEEvaluationGeneral (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component) - : - 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")); -} - - -template -inline -FEEvaluationGeneral -::FEEvaluationGeneral (const FEEvaluationGeneral &other) - : - BaseClass (other) -{ - set_data_pointers(); -} - - - template inline @@ -4507,23 +4619,59 @@ void FEEvaluationGeneral ::set_data_pointers() { + AssertIndexRange(this->data->dofs_per_cell, dofs_per_cell+2); + const unsigned int desired_dofs_per_cell = this->data->dofs_per_cell; + // set the pointers to the correct position in the data array for (unsigned int c=0; cvalues_dofs[c] = &my_data_array[c*dofs_per_cell]; - this->values_quad[c] = &my_data_array[n_components*dofs_per_cell+c*n_q_points]; + this->values_dofs[c] = &my_data_array[c*desired_dofs_per_cell]; + this->values_quad[c] = &my_data_array[n_components*desired_dofs_per_cell+c*n_q_points]; for (unsigned int d=0; dgradients_quad[c][d] = &my_data_array[n_components*(dofs_per_cell+n_q_points) + this->gradients_quad[c][d] = &my_data_array[n_components*(desired_dofs_per_cell+ + n_q_points) + (c*dim+d)*n_q_points]; for (unsigned int d=0; d<(dim*dim+dim)/2; ++d) - this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+dofs_per_cell) + this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+ + desired_dofs_per_cell) + (c*(dim*dim+dim)+d)*n_q_points]; } } +template +inline +void +FEEvaluationGeneral +::check_dofs_per_cell (const unsigned int given_dofs_per_cell) const +{ +#ifdef DEBUG + if (given_dofs_per_cell != this->data->dofs_per_cell) + { + std::ostringstream str; + str << "Dofs per cell in FEEvaluation* class does not match the element (" + << given_dofs_per_cell << " != " << this->data->dofs_per_cell << "). " + << "Suggestion: "; + if (this->data->dofs_per_cell == dofs_per_cell) + str << "FEEvaluation/FEEvaluationGeneral"; + else if (this->data->dofs_per_cell == + internal::MatrixFreeFunctions::DGP_dofs_per_cell::value) + str << "FEEvaluationDGP"; + else if (this->data->dofs_per_cell == dofs_per_cell+1) + str << "FEEvaluationQ_DG0"; + else + str << "No matching suggestion found"; + + Assert(given_dofs_per_cell == this->data->dofs_per_cell, + ExcMessage(str.str().c_str())); + } +#endif +} + + namespace internal { @@ -4656,6 +4804,19 @@ namespace internal out += nn; break; case 1: + ++in; + ++out; + // faces 2 and 3 in 3D use local coordinate system zx, which + // is the other way around compared to the tensor + // product. Need to take that into account. + if (dim == 3) + { + if (dof_to_quad) + out += fe_degree; + else + in += fe_degree; + } + break; case 2: ++in; ++out; @@ -4664,10 +4825,15 @@ namespace internal Assert (false, ExcNotImplemented()); } } - if (face_direction == 1) + if (face_direction == 1 && dim == 3) { in += mm*(mm-1); out += nn*(nn-1); + // adjust for local coordinate system zx + if (dof_to_quad) + out -= (fe_degree+1)*(fe_degree+1)-1; + else + in -= (fe_degree+1)*(fe_degree+1)-1; } } } @@ -5954,6 +6120,7 @@ FEEvaluationGeneral const bool evaluate_grad, const bool evaluate_lapl) { + this->check_dofs_per_cell(dofs_per_cell); Assert (this->dof_values_initialized == true, internal::ExcAccessToUninitializedField()); internal::do_evaluate (*this, this->values_dofs, this->values_quad, @@ -5981,6 +6148,7 @@ FEEvaluationGeneral ::integrate (const bool integrate_val, const bool integrate_grad) { + this->check_dofs_per_cell(dofs_per_cell); if (integrate_val == true) Assert (this->values_quad_submitted == true, internal::ExcAccessToUninitializedField()); @@ -6244,6 +6412,7 @@ FEEvaluation const bool evaluate_grad, const bool evaluate_lapl) { + this->check_dofs_per_cell(dofs_per_cell); Assert (this->dof_values_initialized == true, internal::ExcAccessToUninitializedField()); internal::do_evaluate (*this, this->values_dofs, this->values_quad, @@ -6270,6 +6439,7 @@ void FEEvaluation ::integrate (bool integrate_val,bool integrate_grad) { + this->check_dofs_per_cell(dofs_per_cell); if (integrate_val == true) Assert (this->values_quad_submitted == true, internal::ExcAccessToUninitializedField()); @@ -6433,6 +6603,7 @@ FEEvaluationGL const bool evaluate_grad, const bool evaluate_lapl) { + this->check_dofs_per_cell(dofs_per_cell); Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); Assert (this->dof_values_initialized == true, @@ -6539,6 +6710,7 @@ void FEEvaluationGL ::integrate (const bool integrate_val, const bool integrate_grad) { + this->check_dofs_per_cell(dofs_per_cell); Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); if (integrate_val == true) @@ -6634,9 +6806,8 @@ FEEvaluationDGP : 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]; + AssertDimension(static_cast(this->values_quad[0]-this->values_dofs[0]), + n_components * dofs_per_cell); } @@ -6651,8 +6822,8 @@ FEEvaluationDGP : BaseClass (geometry, dof_handler, first_selected_component) { - for (unsigned int c=0; cvalues_dofs[c] = &this->my_data_array[c*dofs_per_cell]; + AssertDimension(static_cast(this->values_quad[0]-this->values_dofs[0]), + n_components * dofs_per_cell); } @@ -6665,8 +6836,6 @@ FEEvaluationDGP : BaseClass (other) { - for (unsigned int c=0; cvalues_dofs[c] = &this->my_data_array[c*dofs_per_cell]; } @@ -6680,6 +6849,7 @@ FEEvaluationDGP const bool evaluate_grad, const bool evaluate_lapl) { + this->check_dofs_per_cell(dofs_per_cell); Assert (this->dof_values_initialized == true, internal::ExcAccessToUninitializedField()); @@ -6733,6 +6903,7 @@ void FEEvaluationDGP ::integrate (bool integrate_val,bool integrate_grad) { + this->check_dofs_per_cell(dofs_per_cell); if (integrate_val == true) Assert (this->values_quad_submitted == true, internal::ExcAccessToUninitializedField()); @@ -6773,6 +6944,122 @@ FEEvaluationDGP +/*------------------------- FEEvaluationQ_DG0 -------------------------------*/ + +template +inline +FEEvaluationQ_DG0 +::FEEvaluationQ_DG0 (const MatrixFree &data_in, + const unsigned int fe_no, + const unsigned int quad_no) + : + BaseClass (data_in, fe_no, quad_no) +{ + AssertDimension(static_cast(this->values_quad[0]-this->values_dofs[0]), + n_components * dofs_per_cell); +} + + + +template +inline +FEEvaluationQ_DG0 +::FEEvaluationQ_DG0 (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component) + : + BaseClass (geometry, dof_handler, first_selected_component) +{ + AssertDimension(static_cast(this->values_quad[0]-this->values_dofs[0]), + n_components * dofs_per_cell); +} + + + +template +inline +FEEvaluationQ_DG0 +::FEEvaluationQ_DG0 (const FEEvaluationQ_DG0 &other) + : + BaseClass (other) +{ +} + + + +template +inline +void +FEEvaluationQ_DG0 +::evaluate (const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl) +{ + this->check_dofs_per_cell(dofs_per_cell); + Assert (this->dof_values_initialized == true, + internal::ExcAccessToUninitializedField()); + internal::do_evaluate (*this, this->values_dofs, this->values_quad, + this->gradients_quad, this->hessians_quad, + evaluate_val, evaluate_grad, evaluate_lapl, + internal::int2type()); + if (evaluate_val) + for (unsigned int c=0; cvalues_quad[c][q] += this->values_dofs[c][dofs_per_cell-1]; + +#ifdef DEBUG + if (evaluate_val == true) + this->values_quad_initialized = true; + if (evaluate_grad == true) + this->gradients_quad_initialized = true; + if (evaluate_lapl == true) + this->hessians_quad_initialized = true; +#endif +} + + + +template +inline +void +FEEvaluationQ_DG0 +::integrate (bool integrate_val,bool integrate_grad) +{ + this->check_dofs_per_cell(dofs_per_cell); + if (integrate_val == true) + Assert (this->values_quad_submitted == true, + internal::ExcAccessToUninitializedField()); + if (integrate_grad == true) + Assert (this->gradients_quad_submitted == true, + internal::ExcAccessToUninitializedField()); + + if (integrate_val) + for (unsigned int c=0; cvalues_dofs[c][dofs_per_cell-1] = this->values_quad[c][0]; + for (unsigned int q=1; qvalues_dofs[c][dofs_per_cell-1] += this->values_quad[c][q]; + } + else + for (unsigned int c=0; cvalues_dofs[c][dofs_per_cell-1] = VectorizedArray(); + + internal::do_integrate (*this, this->values_dofs, this->values_quad, + this->gradients_quad, integrate_val, integrate_grad, + internal::int2type()); + +#ifdef DEBUG + this->dof_values_initialized = true; +#endif +} + + + #endif // ifndef DOXYGEN diff --git a/deal.II/include/deal.II/matrix_free/shape_info.h b/deal.II/include/deal.II/matrix_free/shape_info.h index 3164b4679c..f9b51e21f6 100644 --- a/deal.II/include/deal.II/matrix_free/shape_info.h +++ b/deal.II/include/deal.II/matrix_free/shape_info.h @@ -72,55 +72,40 @@ namespace internal const unsigned int base_element = 0); /** - * Returns the memory consumption of this - * class in bytes. + * Returns the memory consumption of this class in bytes. */ std::size_t memory_consumption () const; /** - * Stores the shape values of the 1D finite - * element evaluated on all 1D quadrature - * points in vectorized format, i.e., as an - * array of - * VectorizedArray::n_array_elements - * equal elements. The length of this array is - * n_dofs_1d * n_q_points_1d and - * quadrature points are the index running - * fastest. + * Stores the shape values of the 1D finite element evaluated on all 1D + * quadrature points in vectorized format, i.e., as an array of + * VectorizedArray::n_array_elements equal elements. The length of + * this array is n_dofs_1d * n_q_points_1d and quadrature + * points are the index running fastest. */ AlignedVector > shape_values; /** - * Stores the shape gradients of the 1D finite - * element evaluated on all 1D quadrature - * points in vectorized format, i.e., as an - * array of - * VectorizedArray::n_array_elements - * equal elements. The length of this array is - * n_dofs_1d * n_q_points_1d and - * quadrature points are the index running - * fastest. + * Stores the shape gradients of the 1D finite element evaluated on all + * 1D quadrature points in vectorized format, i.e., as an array of + * VectorizedArray::n_array_elements equal elements. The length of + * this array is n_dofs_1d * n_q_points_1d and quadrature + * points are the index running fastest. */ AlignedVector > shape_gradients; /** - * Stores the shape Hessians of the 1D finite - * element evaluated on all 1D quadrature - * points in vectorized format, i.e., as an - * array of - * VectorizedArray::n_array_elements - * equal elements. The length of this array is - * n_dofs_1d * n_q_points_1d and - * quadrature points are the index running - * fastest. + * Stores the shape Hessians of the 1D finite element evaluated on all + * 1D quadrature points in vectorized format, i.e., as an array of + * VectorizedArray::n_array_elements equal elements. The length of + * this array is n_dofs_1d * n_q_points_1d and quadrature + * points are the index running fastest. */ AlignedVector > shape_hessians; /** - * Stores the indices from cell DoFs to face - * DoFs. The rows go through the - * 2*dim faces, and the columns the - * DoFs on the faces. + * Stores the indices from cell DoFs to face DoFs. The rows go through + * the 2*dim faces, and the columns the DoFs on the faces. */ dealii::Table<2,unsigned int> face_indices; @@ -143,15 +128,14 @@ namespace internal std::vector subface_value[2]; /** - * Non-vectorized version of shape - * values. Needed when evaluating face info. + * Non-vectorized version of shape values. Needed when evaluating face + * info. */ std::vector shape_values_number; /** - * Non-vectorized version of shape - * gradients. Needed when evaluating face - * info. + * Non-vectorized version of shape gradients. Needed when evaluating + * face info. */ std::vector shape_gradient_number; @@ -165,26 +149,28 @@ namespace internal std::vector lexicographic_numbering; /** - * Stores the number of quadrature points in - * @p dim dimensions for a cell. + * Stores the degree of the element. + */ + unsigned int fe_degree; + + /** + * Stores the number of quadrature points in @p dim dimensions for a + * cell. */ unsigned int n_q_points; /** - * Stores the number of DoFs per cell in @p - * dim dimensions. + * Stores the number of DoFs per cell in @p dim dimensions. */ unsigned int dofs_per_cell; /** - * Stores the number of quadrature points per - * face in @p dim dimensions. + * Stores the number of quadrature points per face in @p dim dimensions. */ unsigned int n_q_points_face; /** - * Stores the number of DoFs per face in @p - * dim dimensions. + * Stores the number of DoFs per face in @p dim dimensions. */ unsigned int dofs_per_face; }; @@ -200,6 +186,7 @@ namespace internal const FiniteElement &fe_in, const unsigned int base_element_number) : + fe_degree (0), n_q_points (0), dofs_per_cell (0) { diff --git a/deal.II/include/deal.II/matrix_free/shape_info.templates.h b/deal.II/include/deal.II/matrix_free/shape_info.templates.h index eb8e2bda75..826649644f 100644 --- a/deal.II/include/deal.II/matrix_free/shape_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/shape_info.templates.h @@ -22,6 +22,7 @@ #include #include #include +#include #include @@ -59,8 +60,9 @@ namespace internal Assert (fe->n_components() == 1, ExcMessage("FEEvaluation only works for scalar finite elements.")); + fe_degree = fe->degree; - const unsigned int n_dofs_1d = fe->degree+1, + const unsigned int n_dofs_1d = fe_degree+1, n_q_points_1d = quad.size(); // renumber (this is necessary for FE_Q, for example, since there the @@ -82,6 +84,8 @@ namespace internal const FE_DGP *fe_dgp = dynamic_cast*>(fe); + const FE_Q_DG0 *fe_q_dg0 = dynamic_cast*>(fe); + if (fe_poly != 0) scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse(); else if (fe_poly_piece != 0) @@ -89,9 +93,11 @@ namespace internal else if (fe_dgp != 0) { scalar_lexicographic.resize(fe_dgp->dofs_per_cell); - for (unsigned int i=0; idofs_per_cell; ++i) - scalar_lexicographic[i] = i; + for (unsigned int i=0; idofs_per_cell; ++i) + scalar_lexicographic[i] = i; } + else if (fe_q_dg0 != 0) + scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse(); else Assert(false, ExcNotImplemented()); diff --git a/tests/matrix_free/matrix_vector_17.cc b/tests/matrix_free/matrix_vector_17.cc new file mode 100644 index 0000000000..4ab84fd056 --- /dev/null +++ b/tests/matrix_free/matrix_vector_17.cc @@ -0,0 +1,248 @@ +// --------------------------------------------------------------------- +// $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. +// +// --------------------------------------------------------------------- + + + +// this test is similar to matrix_vector_06, but implements the operations for +// FE_QDG0 instead of FE_Q (where there is an additional discontinuous degree +// of freedom and different routines need to be used). The data is still not +// very useful because the matrix does not include face terms actually present +// in an approximation of the Laplacian. It only contains cell terms. + +#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_qdg0 (const MatrixFree &data, + VECTOR &dst, + const VECTOR &src, + const std::pair &cell_range) +{ + typedef typename VECTOR::value_type Number; + FEEvaluationQ_DG0 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_qdg0; + 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+1); + 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+1); + + 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_Q_DG0 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_17.output b/tests/matrix_free/matrix_vector_17.output new file mode 100644 index 0000000000..2bbf553f65 --- /dev/null +++ b/tests/matrix_free/matrix_vector_17.output @@ -0,0 +1,13 @@ + +DEAL:2d::Testing FE_Q_DG0<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q_DG0<2>(2) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FE_Q_DG0<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q_DG0<3>(2) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_stokes_qdg0.cc b/tests/matrix_free/matrix_vector_stokes_qdg0.cc new file mode 100644 index 0000000000..c746da6c99 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_qdg0.cc @@ -0,0 +1,333 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2013 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). + +#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); + FEEvaluationQ_DG0 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::local_apply, + this, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void test () +{ + 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; + + { + BlockCompressedSimpleSparsityPattern 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(); + 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_qdg0.output b/tests/matrix_free/matrix_vector_stokes_qdg0.output new file mode 100644 index 0000000000..db6c7de2a9 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_qdg0.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:: -- 2.39.5