From 30ee7565f32ca9479abadf588f5653480d82ea02 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 24 Feb 2014 22:09:03 +0000 Subject: [PATCH] Change name of a class to better represent what it is doing, add inverse mass operator for DG with matrix-free. Test all those new things. git-svn-id: https://svn.dealii.org/trunk@32546 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/matrix_free/fe_evaluation.h | 355 ++++++++++++------ ...ted_geometry.h => mapping_fe_evaluation.h} | 103 +++-- .../matrix_free/matrix_free.templates.h | 123 +++--- .../include/deal.II/matrix_free/operators.h | 224 +++++++++++ .../include/deal.II/matrix_free/shape_info.h | 19 +- .../matrix_free/shape_info.templates.h | 93 +++-- tests/matrix_free/inverse_mass_01.cc | 216 +++++++++++ tests/matrix_free/inverse_mass_01.output | 41 ++ tests/matrix_free/inverse_mass_02.cc | 205 ++++++++++ tests/matrix_free/inverse_mass_02.output | 26 ++ tests/matrix_free/inverse_mass_03.cc | 220 +++++++++++ tests/matrix_free/inverse_mass_03.output | 26 ++ tests/matrix_free/matrix_vector_14.cc | 8 +- tests/matrix_free/matrix_vector_15.cc | 242 ++++++++++++ tests/matrix_free/matrix_vector_15.output | 25 ++ tests/matrix_free/matrix_vector_16.cc | 237 ++++++++++++ tests/matrix_free/matrix_vector_16.output | 16 + 17 files changed, 1900 insertions(+), 279 deletions(-) rename deal.II/include/deal.II/matrix_free/{evaluated_geometry.h => mapping_fe_evaluation.h} (65%) create mode 100644 deal.II/include/deal.II/matrix_free/operators.h create mode 100644 tests/matrix_free/inverse_mass_01.cc create mode 100644 tests/matrix_free/inverse_mass_01.output create mode 100644 tests/matrix_free/inverse_mass_02.cc create mode 100644 tests/matrix_free/inverse_mass_02.output create mode 100644 tests/matrix_free/inverse_mass_03.cc create mode 100644 tests/matrix_free/inverse_mass_03.output create mode 100644 tests/matrix_free/matrix_vector_15.cc create mode 100644 tests/matrix_free/matrix_vector_15.output create mode 100644 tests/matrix_free/matrix_vector_16.cc create mode 100644 tests/matrix_free/matrix_vector_16.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 f49488d2c7..bd755ee0ac 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -27,7 +27,7 @@ #include #include #include -#include +#include DEAL_II_NAMESPACE_OPEN @@ -115,6 +115,18 @@ public: */ internal::MatrixFreeFunctions::CellType get_cell_type() const; + /** + * Returns a reference to the ShapeInfo object currently in use. + */ + const internal::MatrixFreeFunctions::ShapeInfo & + get_shape_info() const; + + /** + * Fills the JxW values currently used. + */ + void + fill_JxW_values(AlignedVector > &JxW_values) const; + //@} /** @@ -136,6 +148,14 @@ public: * ConstraintMatrix::distribute. When accessing vector entries during the * solution of linear systems, the temporary solution should always have * homogeneous constraints and this method is the correct one. + * + * If the class was constructed through a MappingFEEvaluation object, only + * one single cell is used by this class and this function extracts the + * values of the underlying components on this cell. This call is slower + * than the ones done through a MatrixFree object and lead to a structure + * that does not effectively use vectorization in the evaluate routines + * based on these values (instead, VectorizedArray::n_array_elements + * same copies are worked on). */ template void read_dof_values (const VectorType &src); @@ -171,6 +191,15 @@ public: * function is also necessary when inhomogeneous constraints are to be used, * as MatrixFree can only handle homogeneous constraints. Note that if * vectorization is enabled, the DoF values for several cells are set. + * + * If the class was constructed through a MappingFEEvaluation object, only + * one single cell is used by this class and this function extracts the + * values of the underlying components on this cell. This call is slower + * than the ones done through a MatrixFree object and lead to a structure + * that does not effectively use vectorization in the evaluate routines + * based on these values (instead, VectorizedArray::n_array_elements + * same copies are worked on). In that case, no constraints can be + * processed as these are not available here. */ template void read_dof_values_plain (const VectorType &src); @@ -199,8 +228,17 @@ public: * Takes the values stored internally on dof values of the current cell and * sums them into the vector @p dst. The function also applies constraints * during the write operation. The functionality is hence similar to the - * function ConstraintMatrix::distribute_local_to_global. Note that if - * vectorization is enabled, the DoF values for several cells are used. + * function ConstraintMatrix::distribute_local_to_global. If vectorization + * is enabled, the DoF values for several cells are used. + * + * If the class was constructed through a MappingFEEvaluation object, only + * one single cell is used by this class and this function extracts the + * values of the underlying components on this cell. This call is slower + * than the ones done through a MatrixFree object and lead to a structure + * that does not effectively use vectorization in the evaluate routines + * based on these values (instead, VectorizedArray::n_array_elements + * same copies are worked on). In that case, no constraints can be + * processed as these are not available here. */ template void distribute_local_to_global (VectorType &dst) const; @@ -211,8 +249,8 @@ public: * and sums them into the collection of vectors vector @p dst, starting at * index @p first_index. The function also applies constraints during the * write operation. The functionality is hence similar to the function - * ConstraintMatrix::distribute_local_to_global. Note that if vectorization - * is enabled, the DoF values for several cells are used. + * ConstraintMatrix::distribute_local_to_global. If vectorization is + * enabled, the DoF values for several cells are used. */ template void distribute_local_to_global (std::vector &dst, @@ -232,6 +270,15 @@ public: * during the write operation. The functionality is hence similar to the * function ConstraintMatrix::distribute_local_to_global. Note that if * vectorization is enabled, the DoF values for several cells are used. + * + * If the class was constructed through a MappingFEEvaluation object, only + * one single cell is used by this class and this function extracts the + * values of the underlying components on this cell. This call is slower + * than the ones done through a MatrixFree object and lead to a structure + * that does not effectively use vectorization in the evaluate routines + * based on these values (instead, VectorizedArray::n_array_elements + * same copies are worked on). In that case, no constraints can be + * processed as these are not available here. */ template void set_dof_values (VectorType &dst) const; @@ -512,7 +559,7 @@ protected: /** * Constructor that comes with reduced functionality and works similar as - * FEValues. The user has to provide a structure of type EvaluatedGeometry + * 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 @@ -526,9 +573,9 @@ protected: * operations), but only possibly within the element if the * evaluate/integrate routines are combined (e.g. for matrix assembly). */ - FEEvaluationBase (const EvaluatedGeometry &geometry, - const DoFHandler &dof_handler, - const unsigned int base_element = 0); + FEEvaluationBase (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component = 0); /** * Copy constructor @@ -638,13 +685,18 @@ protected: const internal::MatrixFreeFunctions::MappingInfo *mapping_info; /** - * Stores a reference to the unit cell data, i.e., values, gradients and + * In case the class is initialized from MappingFEEvaluation instead of + * MatrixFree, this data structure holds the evaluated shape data. + */ + std_cxx1x::shared_ptr > stored_shape_info; + + /** + * Stores a pointer to the unit cell shape data, i.e., values, gradients and * Hessians in 1D at the quadrature points that constitute the tensor * product. Also contained in matrix_info, but it simplifies code if we - * store a reference to it. If the object is initialized without MatrixFree - * object, the constructor creates this data structure. + * store a reference to it. */ - std_cxx1x::shared_ptr > data; + const internal::MatrixFreeFunctions::ShapeInfo *data; /** * A pointer to the Cartesian Jacobian information of the present cell. Only @@ -754,12 +806,24 @@ protected: /** * Geometry data generated by FEValues on the fly. */ - SmartPointer > evaluated_geometry; + SmartPointer > mapped_geometry; /** * A pointer to the underlying DoFHandler. */ const DoFHandler *dof_handler; + + /** + * For a DoFHandler with more than one finite element, select at which + * component this data structure should start. + */ + const unsigned int first_selected_component; + + /** + * A temporary data structure necessary to read degrees of freedom when no + * MatrixFree object was given at initialization. + */ + mutable std::vector local_dof_indices; }; @@ -799,7 +863,7 @@ protected: /** * Constructor that comes with reduced functionality and works similar as - * FEValues. The user has to provide a structure of type EvaluatedGeometry + * 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 @@ -813,9 +877,9 @@ protected: * operations), but only possibly within the element if the * evaluate/integrate routines are combined (e.g. for matrix assembly). */ - FEEvaluationAccess (const EvaluatedGeometry &geometry, - const DoFHandler &dof_handler, - const unsigned int base_element = 0); + FEEvaluationAccess (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component = 0); /** * Copy constructor @@ -946,7 +1010,7 @@ protected: /** * Constructor that comes with reduced functionality and works similar as - * FEValues. The user has to provide a structure of type EvaluatedGeometry + * 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 @@ -960,9 +1024,9 @@ protected: * operations), but only possibly within the element if the * evaluate/integrate routines are combined (e.g. for matrix assembly). */ - FEEvaluationAccess (const EvaluatedGeometry &geometry, - const DoFHandler &dof_handler, - const unsigned int base_element = 0); + FEEvaluationAccess (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component = 0); /** * Copy constructor @@ -1102,7 +1166,7 @@ protected: /** * Constructor that comes with reduced functionality and works similar as - * FEValues. The user has to provide a structure of type EvaluatedGeometry + * 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 @@ -1116,9 +1180,9 @@ protected: * operations), but only possibly within the element if the * evaluate/integrate routines are combined (e.g. for matrix assembly). */ - FEEvaluationAccess (const EvaluatedGeometry &geometry, - const DoFHandler &dof_handler, - const unsigned int base_element = 0); + FEEvaluationAccess (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component = 0); /** * Copy constructor @@ -1193,7 +1257,7 @@ public: /** * Constructor that comes with reduced functionality and works similar as - * FEValues. The user has to provide a structure of type EvaluatedGeometry + * 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 @@ -1207,9 +1271,9 @@ public: * operations), but only possibly within the element if the * evaluate/integrate routines are combined (e.g. for matrix assembly). */ - FEEvaluationGeneral (const EvaluatedGeometry &geometry, - const DoFHandler &dof_handler, - const unsigned int base_element = 0); + FEEvaluationGeneral (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component = 0); /** * Copy constructor @@ -1309,7 +1373,7 @@ private: * possible to use vectorization for applying a vector operation for several * cells at once. The second form of usage is to initialize it from geometry * information generated by FEValues, which is stored in the class - * EvaluatedGeometry. Here, the operations can only work on a single cell, but + * MappingFEEvaluation. Here, the operations can only work on a single cell, but * possibly be vectorized by combining several operations (e.g. when * performing matrix assembly). * @@ -1379,7 +1443,7 @@ public: /** * Constructor that comes with reduced functionality and works similar as - * FEValues. The user has to provide a structure of type EvaluatedGeometry + * 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 @@ -1393,9 +1457,9 @@ public: * operations), but only possibly within the element if the * evaluate/integrate routines are combined (e.g. for matrix assembly). */ - FEEvaluation (const EvaluatedGeometry &geometry, - const DoFHandler &dof_handler, - const unsigned int base_element = 0); + FEEvaluation (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component = 0); /** * Evaluates the function values, the gradients, and the Laplacians of the @@ -1533,7 +1597,7 @@ public: /** * Constructor that comes with reduced functionality and works similar as - * FEValues. The user has to provide a structure of type EvaluatedGeometry + * 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 @@ -1547,9 +1611,9 @@ public: * operations), but only possibly within the element if the * evaluate/integrate routines are combined (e.g. for matrix assembly). */ - FEEvaluationGL (const EvaluatedGeometry &geometry, - const DoFHandler &dof_handler, - const unsigned int base_element = 0); + FEEvaluationGL (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component = 0); /** * Copy constructor @@ -1660,7 +1724,7 @@ public: /** * Constructor that comes with reduced functionality and works similar as - * FEValues. The user has to provide a structure of type EvaluatedGeometry + * 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 @@ -1674,9 +1738,9 @@ public: * operations), but only possibly within the element if the * evaluate/integrate routines are combined (e.g. for matrix assembly). */ - FEEvaluationDGP (const EvaluatedGeometry &geometry, - const DoFHandler &dof_handler, - const unsigned int base_element = 0); + FEEvaluationDGP (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component = 0); /** * Copy constructor @@ -1712,33 +1776,6 @@ public: #ifndef DOXYGEN -namespace internal -{ - namespace MatrixFreeFunctions - { - // a small class that gives control over the delete behavior of - // std::shared_ptr: we need to disable it when we initialize a pointer - // from another structure. - template - struct DummyDeleter - { - DummyDeleter (const bool do_delete = false) - : - do_delete(do_delete) - {} - - void operator () (CLASS *pointer) - { - if (do_delete) - delete pointer; - } - - const bool do_delete; - }; - } -} - - /*----------------------- FEEvaluationBase ----------------------------------*/ @@ -1763,9 +1800,7 @@ FEEvaluationBase mapping_info (&data_in.get_mapping_info()), data (&data_in.get_shape_info (fe_no_in, quad_no_in, active_fe_index, - active_quad_index), - internal::MatrixFreeFunctions::DummyDeleter - >(false)), + active_quad_index)), cartesian_data (0), jacobian (0), J_value (0), @@ -1776,7 +1811,9 @@ FEEvaluationBase jacobian_grad_upper(0), cell (numbers::invalid_unsigned_int), cell_type (internal::MatrixFreeFunctions::undefined), - cell_data_number (0) + cell_data_number (0), + dof_handler (0), + first_selected_component (0) { for (unsigned int c=0; c template inline FEEvaluationBase -::FEEvaluationBase (const EvaluatedGeometry &geometry, - const DoFHandler &dof_handler_in, - const unsigned int base_element) +::FEEvaluationBase (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler_in, + const unsigned int first_selected_component) : quad_no (-1), n_fe_components (n_components_), @@ -1823,7 +1860,8 @@ FEEvaluationBase matrix_info (0), dof_info (0), mapping_info (0), - data (new internal::MatrixFreeFunctions::ShapeInfo(geometry.get_quadrature(), dof_handler_in.get_fe(), base_element)), + stored_shape_info (new internal::MatrixFreeFunctions::ShapeInfo(geometry.get_quadrature(), dof_handler_in.get_fe(), 0)), + data (stored_shape_info.get()), cartesian_data (0), jacobian (geometry.get_inverse_jacobians().begin()), J_value (geometry.get_JxW_values().begin()), @@ -1834,8 +1872,9 @@ FEEvaluationBase cell (0), cell_type (internal::MatrixFreeFunctions::general), cell_data_number (0), - evaluated_geometry (&geometry), - dof_handler (&dof_handler_in) + mapped_geometry (&geometry), + dof_handler (&dof_handler_in), + first_selected_component (first_selected_component) { for (unsigned int c=0; c for (unsigned int d=0; d<(dim*dim+dim)/2; ++d) hessians_quad[c][d] = 0; } - Assert(dof_handler->get_fe().element_multiplicity(base_element) == 1 || - dof_handler->get_fe().element_multiplicity(base_element) >= n_components_, + Assert(dof_handler->get_fe().element_multiplicity(0) == 1 || + dof_handler->get_fe().element_multiplicity(0)-first_selected_component >= n_components_, ExcMessage("The underlying element must at least contain as many " "components as requested by this class")); } @@ -1877,8 +1916,9 @@ FEEvaluationBase cell (other.cell), cell_type (other.cell_type), cell_data_number (other.cell_data_number), - evaluated_geometry (other.evaluated_geometry), - dof_handler (other.dof_handler) + mapped_geometry (other.mapped_geometry), + dof_handler (other.dof_handler), + first_selected_component (other.first_selected_component) { for (unsigned int c=0; c - template inline void -FEEvaluationBase -::reinit (const unsigned int cell_in) +FEEvaluationBase::reinit (const unsigned int cell_in) { - Assert (evaluated_geometry == 0, ExcMessage("FEEvaluation was initialized without a matrix-free object. Integer indexing is not possible")); - if (evaluated_geometry != 0) + Assert (mapped_geometry == 0, ExcMessage("FEEvaluation was initialized without a matrix-free object. Integer indexing is not possible")); + if (mapped_geometry != 0) return; Assert (dof_info != 0, ExcNotInitialized()); Assert (mapping_info != 0, ExcNotInitialized()); @@ -1988,8 +2026,7 @@ FEEvaluationBase template inline internal::MatrixFreeFunctions::CellType -FEEvaluationBase -::get_cell_type () const +FEEvaluationBase::get_cell_type () const { Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); return cell_type; @@ -1997,6 +2034,40 @@ FEEvaluationBase +template +inline +const internal::MatrixFreeFunctions::ShapeInfo & +FEEvaluationBase::get_shape_info() const +{ + Assert(data != 0, ExcInternalError()); + return *data; +} + + + +template +inline +void +FEEvaluationBase +::fill_JxW_values(AlignedVector > &JxW_values) const +{ + AssertDimension(JxW_values.size(), data->n_q_points); + Assert (this->J_value != 0, ExcNotImplemented()); + if (this->cell_type == internal::MatrixFreeFunctions::cartesian || + this->cell_type == internal::MatrixFreeFunctions::affine) + { + Assert (this->mapping_info != 0, ExcNotImplemented()); + VectorizedArray J = this->J_value[0]; + for (unsigned int q=0; qdata->n_q_points; ++q) + JxW_values[q] = J * this->quadrature_weights[q]; + } + else + for (unsigned int q=0; qn_q_points; ++q) + JxW_values[q] = this->J_value[q]; +} + + + namespace internal { // write access to generic vectors that have operator (). @@ -2086,6 +2157,14 @@ namespace internal res = vector_access (const_cast(vec), index); } + template + void process_dof_global (const types::global_dof_index index, + VectorType &vec, + Number &res) const + { + res = const_cast(vec)(index); + } + void pre_constraints (const Number &, Number &res) const { @@ -2125,6 +2204,14 @@ namespace internal vector_access (vec, index) += res; } + template + void process_dof_global (const types::global_dof_index index, + VectorType &vec, + Number &res) const + { + vec(index) += res; + } + void pre_constraints (const Number &input, Number &res) const { @@ -2163,6 +2250,14 @@ namespace internal vector_access (vec, index) = res; } + template + void process_dof_global (const types::global_dof_index index, + VectorType &vec, + Number &res) const + { + vec(index) = res; + } + void pre_constraints (const Number &, Number &) const { @@ -2235,7 +2330,32 @@ FEEvaluationBase // into the local data field or write local data into the vector. Certain // operations are no-ops for the given use case. - Assert (matrix_info != 0, ExcNotInitialized()); + // Case 1: No MatrixFree object given, simple case because we do not need to + // process constraints and need not care about vectorization + if (matrix_info == 0) + { + Assert (dof_handler != 0, ExcNotInitialized()); + typename DoFHandler::cell_iterator cell (&dof_handler->get_tria(), + mapped_geometry->get_cell()->level(), + mapped_geometry->get_cell()->index(), + dof_handler); + local_dof_indices.resize(dof_handler->get_fe().dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + + unsigned int index = first_selected_component * this->data->dofs_per_cell; + for (unsigned int comp = 0; compdata->dofs_per_cell; ++i, ++index) + { + operation.process_dof_global(local_dof_indices[this->data->lexicographic_numbering[index]], + *src[0], values_dofs[comp][i][0]); + for (unsigned int v=1; v::n_array_elements; ++v) + operation.process_empty(values_dofs[comp][i][v]); + } + } + return; + } + Assert (dof_info != 0, ExcNotInitialized()); Assert (matrix_info->indices_initialized() == true, ExcNotInitialized()); @@ -2843,9 +2963,16 @@ void FEEvaluationBase ::read_dof_values_plain (const VectorType *src[]) { + // Case without MatrixFree initialization object + if (matrix_info == 0) + { + internal::VectorReader reader; + read_write_operation (reader, src); + return; + } + // this is different from the other three operations because we do not use // constraints here, so this is a separate function. - Assert (matrix_info != 0, ExcNotInitialized()); Assert (dof_info != 0, ExcNotInitialized()); Assert (matrix_info->indices_initialized() == true, ExcNotInitialized()); @@ -3549,11 +3676,11 @@ FEEvaluationAccess template inline FEEvaluationAccess -::FEEvaluationAccess (const EvaluatedGeometry &geometry, +::FEEvaluationAccess (const MappingFEEvaluation &geometry, const DoFHandler &dof_handler, - const unsigned int base_element) + const unsigned int first_selected_component) : - FEEvaluationBase (geometry, dof_handler, base_element) + FEEvaluationBase (geometry, dof_handler, first_selected_component) {} @@ -3590,11 +3717,11 @@ FEEvaluationAccess template inline FEEvaluationAccess -::FEEvaluationAccess (const EvaluatedGeometry &geometry, +::FEEvaluationAccess (const MappingFEEvaluation &geometry, const DoFHandler &dof_handler, - const unsigned int base_element) + const unsigned int first_selected_component) : - FEEvaluationBase (geometry, dof_handler, base_element) + FEEvaluationBase (geometry, dof_handler, first_selected_component) {} @@ -3824,11 +3951,11 @@ FEEvaluationAccess template inline FEEvaluationAccess -::FEEvaluationAccess (const EvaluatedGeometry &geometry, +::FEEvaluationAccess (const MappingFEEvaluation &geometry, const DoFHandler &dof_handler, - const unsigned int base_element) + const unsigned int first_selected_component) : - FEEvaluationBase (geometry, dof_handler, base_element) + FEEvaluationBase (geometry, dof_handler, first_selected_component) {} @@ -4284,11 +4411,11 @@ template inline FEEvaluationGeneral -::FEEvaluationGeneral (const EvaluatedGeometry &geometry, +::FEEvaluationGeneral (const MappingFEEvaluation &geometry, const DoFHandler &dof_handler, - const unsigned int base_element) + const unsigned int first_selected_component) : - BaseClass (geometry, dof_handler, base_element) + BaseClass (geometry, dof_handler, first_selected_component) { set_data_pointers(); } @@ -5917,11 +6044,11 @@ template inline FEEvaluation -::FEEvaluation (const EvaluatedGeometry &geometry, +::FEEvaluation (const MappingFEEvaluation &geometry, const DoFHandler &dof_handler, - const unsigned int base_element) + const unsigned int first_selected_component) : - BaseClass (geometry, dof_handler, base_element) + BaseClass (geometry, dof_handler, first_selected_component) { compute_even_odd_factors(); } @@ -6213,11 +6340,11 @@ FEEvaluationGL template inline FEEvaluationGL -::FEEvaluationGL (const EvaluatedGeometry &geometry, +::FEEvaluationGL (const MappingFEEvaluation &geometry, const DoFHandler &dof_handler, - const unsigned int base_element) + const unsigned int first_selected_component) : - BaseClass (geometry, dof_handler, base_element) + BaseClass (geometry, dof_handler, first_selected_component) {} @@ -6448,11 +6575,11 @@ template inline FEEvaluationDGP -::FEEvaluationDGP (const EvaluatedGeometry &geometry, +::FEEvaluationDGP (const MappingFEEvaluation &geometry, const DoFHandler &dof_handler, - const unsigned int base_element) + const unsigned int first_selected_component) : - BaseClass (geometry, dof_handler, base_element) + BaseClass (geometry, dof_handler, first_selected_component) {} diff --git a/deal.II/include/deal.II/matrix_free/evaluated_geometry.h b/deal.II/include/deal.II/matrix_free/mapping_fe_evaluation.h similarity index 65% rename from deal.II/include/deal.II/matrix_free/evaluated_geometry.h rename to deal.II/include/deal.II/matrix_free/mapping_fe_evaluation.h index 28d486ca53..2c42f1e98f 100644 --- a/deal.II/include/deal.II/matrix_free/evaluated_geometry.h +++ b/deal.II/include/deal.II/matrix_free/mapping_fe_evaluation.h @@ -15,8 +15,8 @@ // --------------------------------------------------------------------- -#ifndef __deal2__matrix_free_evaluated_geometry_h -#define __deal2__matrix_free_evaluated_geometry_h +#ifndef __deal2__matrix_free_mapping_fe_evaluation_h +#define __deal2__matrix_free_mapping_fe_evaluation_h #include @@ -27,39 +27,62 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN /** - * The class makes FEEvaluation with the mapping information generated by - * FEValues. + * This class provides evaluated mapping information using standard deal.II + * information in a form that FEEvaluation and friends can use for vectorized + * access. Since no vectorization over cells is available with the + * DoFHandler/Triangulation cell iterators, the interface to FEEvaluation's + * vectorization model is to use @p VectorizedArray::n_array_element copies of + * the same element. This interface is thus primarily useful for evaluating + * several operators on the same cell, e.g., when assembling cell matrices. + * + * @author Martin Kronbichler */ template -class EvaluatedGeometry : public Subscriptor +class MappingFEEvaluation : public Subscriptor { public: /** - * Constructor, similar to FEValues. + * Constructor, similar to FEValues. Since this class only evaluates the + * geometry, no finite element has to be specified and the simplest element, + * FE_Nothing, is used internally for the underlying FEValues object. */ - EvaluatedGeometry (const Mapping &mapping, - const FiniteElement &fe, + MappingFEEvaluation (const Mapping &mapping, const Quadrature<1> &quadrature, const UpdateFlags update_flags); /** * Constructor. Instead of providing a mapping, use MappingQ1. */ - EvaluatedGeometry (const FiniteElement &fe, - const Quadrature<1> &quadrature, + MappingFEEvaluation (const Quadrature<1> &quadrature, const UpdateFlags update_flags); /** - * Initialize with the given cell iterator. + * Initialize with the given cell iterator. This works for all type of + * iterators that can be converted into a Triangulation::cell_iterator + * (DoFHandler::cell_iterator or Triangulation::cell_iterator). */ template - void reinit(ITERATOR &cell); + void reinit(ITERATOR cell); + + /** + * Return a triangulation iterator to the current cell. + */ + typename Triangulation::cell_iterator get_cell () const; + + /** + * Return a reference to the underlying FEValues object that evaluates + * certain quantities (only mapping-related ones like Jacobians or mapped + * quadrature points are accessible, as no finite element data is actually + * used). + */ + const FEValues& get_fe_values () const; /** * Return a vector of inverse transpose Jacobians. For compatibility with @@ -100,6 +123,12 @@ public: get_quadrature () const; private: + /** + * Dummy finite element object necessary for initializing the FEValues + * object. + */ + FE_Nothing fe_dummy; + /** * An underlying FEValues object that performs the (scalar) evaluation. */ @@ -136,12 +165,11 @@ private: template inline -EvaluatedGeometry::EvaluatedGeometry (const Mapping &mapping, - const FiniteElement &fe, - const Quadrature<1> &quadrature, - const UpdateFlags update_flags) +MappingFEEvaluation::MappingFEEvaluation (const Mapping &mapping, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags) : - fe_values(mapping, fe, Quadrature(quadrature), + fe_values(mapping, fe_dummy, Quadrature(quadrature), internal::MatrixFreeFunctions::MappingInfo::compute_update_flags(update_flags)), quadrature_1d(quadrature), inverse_jacobians(fe_values.get_quadrature().size()), @@ -157,11 +185,10 @@ EvaluatedGeometry::EvaluatedGeometry (const Mapping &mapping, template inline -EvaluatedGeometry::EvaluatedGeometry (const FiniteElement &fe, - const Quadrature<1> &quadrature, - const UpdateFlags update_flags) +MappingFEEvaluation::MappingFEEvaluation (const Quadrature<1> &quadrature, + const UpdateFlags update_flags) : - fe_values(fe, Quadrature(quadrature), + fe_values(fe_dummy, Quadrature(quadrature), internal::MatrixFreeFunctions::MappingInfo::compute_update_flags(update_flags)), quadrature_1d(quadrature), inverse_jacobians(fe_values.get_quadrature().size()), @@ -179,9 +206,9 @@ template template inline void -EvaluatedGeometry::reinit(ITERATOR &cell) +MappingFEEvaluation::reinit(ITERATOR cell) { - fe_values.reinit(cell); + fe_values.reinit(static_cast::cell_iterator>(cell)); for (unsigned int q=0; q::reinit(ITERATOR &cell) +template +inline +typename Triangulation::cell_iterator +MappingFEEvaluation::get_cell() const +{ + return fe_values.get_cell(); +} + + + +template +inline +const FEValues& +MappingFEEvaluation::get_fe_values() const +{ + return fe_values; +} + + + template inline const AlignedVector > >& -EvaluatedGeometry::get_inverse_jacobians() const +MappingFEEvaluation::get_inverse_jacobians() const { return inverse_jacobians; } @@ -214,7 +261,7 @@ EvaluatedGeometry::get_inverse_jacobians() const template inline const AlignedVector > >& -EvaluatedGeometry::get_normal_vectors() const +MappingFEEvaluation::get_normal_vectors() const { return normal_vectors; } @@ -224,7 +271,7 @@ EvaluatedGeometry::get_normal_vectors() const template inline const AlignedVector > >& -EvaluatedGeometry::get_quadrature_points() const +MappingFEEvaluation::get_quadrature_points() const { return quadrature_points; } @@ -234,7 +281,7 @@ EvaluatedGeometry::get_quadrature_points() const template inline const AlignedVector >& -EvaluatedGeometry::get_JxW_values() const +MappingFEEvaluation::get_JxW_values() const { return jxw_values; } @@ -244,7 +291,7 @@ EvaluatedGeometry::get_JxW_values() const template inline const Quadrature<1>& -EvaluatedGeometry::get_quadrature() const +MappingFEEvaluation::get_quadrature() const { return quadrature_1d; } diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h index b50d23bc01..ad08fa64aa 100644 --- a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h +++ b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h @@ -111,6 +111,21 @@ internal_reinit(const Mapping &mapping, const std::vector > &quad, const typename MatrixFree::AdditionalData additional_data) { + + // Reads out the FE information and stores the shape function values, + // gradients and Hessians for quadrature points. + { + const unsigned int n_fe = dof_handler.size(); + const unsigned int n_quad = quad.size(); + shape_info.reinit (TableIndices<4>(n_fe, n_quad, 1, 1)); + for (unsigned int no=0; noget_fe()); + } + } + if (additional_data.initialize_indices == true) { clear(); @@ -191,21 +206,6 @@ internal_reinit(const Mapping &mapping, } } - // Reads out the FE information and stores the shape function values, - // gradients and Hessians for quadrature points. - const unsigned int n_fe = dof_handler.size(); - const unsigned int n_quad = quad.size(); - shape_info.reinit (TableIndices<4>(n_fe, n_quad, 1, 1)); - for (unsigned int no=0; no &fe = dof_handler[no]->get_fe(); - for (unsigned int nq =0; nq &mapping, const std::vector > &quad, const typename MatrixFree::AdditionalData additional_data) { + // Reads out the FE information and stores the shape function values, + // gradients and Hessians for quadrature points. + { + const unsigned int n_components = dof_handler.size(); + const unsigned int n_quad = quad.size(); + unsigned int n_fe_in_collection = 0; + for (unsigned int i=0; iget_fe().size()); + unsigned int n_quad_in_collection = 0; + for (unsigned int q=0; q(n_components, n_quad, + n_fe_in_collection, + n_quad_in_collection)); + for (unsigned int no=0; noget_fe().size(); ++fe_no) + for (unsigned int nq =0; nqget_fe()[fe_no]); + } + if (additional_data.initialize_indices == true) { clear(); @@ -314,30 +337,6 @@ internal_reinit(const Mapping &mapping, } } - // Reads out the FE information and stores the shape function values, - // gradients and Hessians for quadrature points. - const unsigned int n_components = dof_handler.size(); - const unsigned int n_quad = quad.size(); - unsigned int n_fe_in_collection = 0; - for (unsigned int i=0; iget_fe().size()); - unsigned int n_quad_in_collection = 0; - for (unsigned int q=0; q(n_components, n_quad, - n_fe_in_collection, - n_quad_in_collection)); - for (unsigned int no=0; noget_fe().size(); ++fe_no) - { - const FiniteElement &fe = dof_handler[no]->get_fe()[fe_no]; - for (unsigned int nq =0; nq*> &dof_hand template void MatrixFree::initialize_indices (const std::vector &constraint, - const std::vector &locally_owned_set) + const std::vector &locally_owned_set) { const unsigned int n_fe = dof_handlers.n_dof_handlers; const unsigned int n_active_cells = cell_level_index.size(); @@ -507,12 +506,11 @@ void MatrixFree::initialize_indices fes.back()->dofs_per_cell); } - lexicographic_inv[no].resize (fes.size()); for (unsigned int fe_index = 0; fe_index &fe = *fes[fe_index]; Assert (fe.n_base_elements() == 1, - ExcMessage ("MatrixFree only works for DoFHandler with one base element")); + ExcMessage ("MatrixFree currently only works for DoFHandler with one base element")); const unsigned int n_fe_components = fe.element_multiplicity (0); // cache number of finite elements and dofs_per_cell @@ -521,40 +519,7 @@ void MatrixFree::initialize_indices dof_info[no].dimension = dim; dof_info[no].n_components = n_fe_components; - // get permutation that gives lexicographic renumbering of the cell - // dofs renumber (this is necessary for FE_Q, for example, since - // there the vertex DoFs come first, which is incompatible with the - // lexicographic ordering necessary to apply tensor products - // efficiently) - - if (n_fe_components == 1) - { - lexicographic_inv[no][fe_index] = - internal::MatrixFreeFunctions::get_lexicographic_numbering_inverse(fe); - AssertDimension (lexicographic_inv[no][fe_index].size(), - dof_info[no].dofs_per_cell[fe_index]); - } - else - { - // ok, we have more than one component, get the inverse - // permutation, invert it, sort the components one after one, - // and invert back - Assert (n_fe_components > 1, ExcInternalError()); - std::vector scalar_lex = - Utilities::invert_permutation(internal::MatrixFreeFunctions::get_lexicographic_numbering_inverse(fe.base_element(0))); - AssertDimension (scalar_lex.size() * n_fe_components, - dof_info[no].dofs_per_cell[fe_index]); - std::vector lexicographic (dof_info[no].dofs_per_cell[fe_index]); - for (unsigned int comp=0; comp::initialize_indices local_dof_indices.resize (dof_info[no].dofs_per_cell[0]); cell_it->get_dof_indices(local_dof_indices); dof_info[no].read_dof_indices (local_dof_indices, - lexicographic_inv[no][0], + shape_info(no,0,0,0).lexicographic_numbering, *constraint[no], counter, constraint_values, cell_at_boundary); @@ -623,7 +588,7 @@ void MatrixFree::initialize_indices local_dof_indices.resize (dof_info[no].dofs_per_cell[0]); cell_it->get_mg_dof_indices(local_dof_indices); dof_info[no].read_dof_indices (local_dof_indices, - lexicographic_inv[no][0], + shape_info(no,0,0,0).lexicographic_numbering, *constraint[no], counter, constraint_values, cell_at_boundary); @@ -643,7 +608,7 @@ void MatrixFree::initialize_indices local_dof_indices.resize (cell_it->get_fe().dofs_per_cell); cell_it->get_dof_indices(local_dof_indices); dof_info[no].read_dof_indices (local_dof_indices, - lexicographic_inv[no][cell_it->active_fe_index()], + shape_info(no,0,cell_it->active_fe_index(),0).lexicographic_numbering, *constraint[no], counter, constraint_values, cell_at_boundary); diff --git a/deal.II/include/deal.II/matrix_free/operators.h b/deal.II/include/deal.II/matrix_free/operators.h new file mode 100644 index 0000000000..7e6821e344 --- /dev/null +++ b/deal.II/include/deal.II/matrix_free/operators.h @@ -0,0 +1,224 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2011 - 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. +// +// --------------------------------------------------------------------- + + +#ifndef __deal2__matrix_free_operators_h +#define __deal2__matrix_free_operators_h + + +#include +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace MatrixFreeOperators +{ + /** + * This class implements the operation of the action of the inverse of a mass + * matrix on an element for the special case of an evaluation object with as + * many quadrature points as there are cell degrees of freedom. It uses + * algorithms from FEEvaluation and produces the exact mass matrix for DGQ + * elements. This algorithm uses tensor products of inverse 1D shape matrices + * over quadrature points, so the inverse operation is exactly as expensive as + * applying the forward operator on each cell. Of course, for continuous + * finite elements this operation does not produce the inverse of a mass + * operation as the coupling between the elements cannot be considered by this + * operation. + * + * The equation may contain variable coefficients, so the user is required to + * provide an array for the inverse of the local coefficient (this class + * provide a helper method 'fill_inverse_JxW_values' to get the inverse of a + * constant-coefficient operator). + * + * @author Martin Kronbichler, 2014 + */ + template + class CellwiseInverseMassMatrix + { + public: + /** + * Constructor. Initializes the shape information from the ShapeInfo field + * in the class FEEval. + */ + CellwiseInverseMassMatrix (const FEEvaluationBase &fe_eval); + + /** + * Applies the inverse mass matrix operation on an input array. It is + * assumed that the passed input and output arrays are of correct size, + * namely FEEval::dofs_per_cell * n_components long. The inverse of the + * local coefficient (also containing the inverse JxW values) must be passed + * as first argument. Passing more than one component in the coefficient is + * allowed. + */ + void apply(const AlignedVector > &inverse_coefficient, + const unsigned int n_actual_components, + const VectorizedArray *in_array, + VectorizedArray *out_array) const; + + /** + * Fills the given array with the inverse of the JxW values, i.e., a mass + * matrix with coefficient 1. Non-unit coefficients must be multiplied (in + * inverse form) to this array. + */ + void fill_inverse_JxW_values(AlignedVector > &inverse_jxw) const; + + private: + /** + * A reference to the FEEvaluation object for getting the JxW_values. + */ + const FEEvaluationBase &fe_eval; + + /** + * A structure to hold inverse shape functions, put into an aligned vector + * in order to make sure it can be allocated. + */ + struct Inverse1DShape + { + VectorizedArray evenodd[fe_degree+1][fe_degree/2+1]; + }; + AlignedVector inverse_shape; + }; + + + + // ------------------------------------ inline functions --------------------- + + template + inline + CellwiseInverseMassMatrix + ::CellwiseInverseMassMatrix (const FEEvaluationBase &fe_eval) + : + fe_eval (fe_eval) + { + FullMatrix shapes_1d(fe_degree+1, fe_degree+1); + for (unsigned int i=0, c=0; i + inline + void + CellwiseInverseMassMatrix + ::fill_inverse_JxW_values(AlignedVector > &inverse_jxw) const + { + const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; + Assert(inverse_jxw.size() > 0 && + inverse_jxw.size() % dofs_per_cell == 0, + ExcMessage("Expected diagonal to be a multiple of scalar dof per cells")); + + // temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values + // from fe_eval (will not reallocate any memory) + const unsigned int previous_size = inverse_jxw.size(); + inverse_jxw.resize(dofs_per_cell); + fe_eval.fill_JxW_values(inverse_jxw); + + // invert + inverse_jxw.resize_fast(previous_size); + for (unsigned int q=0; q + inline + void + CellwiseInverseMassMatrix + ::apply(const AlignedVector > &inverse_coefficients, + const unsigned int n_actual_components, + const VectorizedArray *in_array, + VectorizedArray *out_array) const + { + const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; + Assert(inverse_coefficients.size() > 0 && + inverse_coefficients.size() % dofs_per_cell == 0, + ExcMessage("Expected diagonal to be a multiple of scalar dof per cells")); + if (inverse_coefficients.size() != dofs_per_cell) + AssertDimension(n_actual_components * dofs_per_cell, inverse_coefficients.size()); + + Assert(dim == 2 || dim == 3, ExcNotImplemented()); + + const unsigned int shift_coefficient = + inverse_coefficients.size() > dofs_per_cell ? dofs_per_cell : 0; + const VectorizedArray *inv_coefficient = &inverse_coefficients[0]; + VectorizedArray temp_data_field[dofs_per_cell]; + for (unsigned int d=0; d* in = in_array+d*dofs_per_cell; + VectorizedArray* out = out_array+d*dofs_per_cell; + internal::apply_tensor_product_evenodd,0,false,false,2> + (inverse_shape[0].evenodd, in, temp_data_field); + internal::apply_tensor_product_evenodd,1,false,false,2> + (inverse_shape[0].evenodd, temp_data_field, out); + + if (dim == 3) + { + internal::apply_tensor_product_evenodd,2,false,false,2> + (inverse_shape[0].evenodd, out, temp_data_field); + for (unsigned int q=0; q,2,true,false,2> + (inverse_shape[0].evenodd, temp_data_field, out); + } + else if (dim == 2) + for (unsigned int q=0; q,1,true,false,2> + (inverse_shape[0].evenodd, out, temp_data_field); + internal::apply_tensor_product_evenodd,0,true,false,2> + (inverse_shape[0].evenodd, temp_data_field, out); + + inv_coefficient += shift_coefficient; + } + } + +} // end of namespace MatrixFreeOperators + + +DEAL_II_NAMESPACE_CLOSE + +#endif 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 7b6057a43b..3164b4679c 100644 --- a/deal.II/include/deal.II/matrix_free/shape_info.h +++ b/deal.II/include/deal.II/matrix_free/shape_info.h @@ -155,6 +155,15 @@ namespace internal */ std::vector shape_gradient_number; + /** + * Renumbering from deal.II's numbering of cell degrees of freedom to + * lexicographic numbering used inside the FEEvaluation schemes of the + * underlying element in the DoFHandler. For vector-valued elements, the + * renumbering starts with a lexicographic numbering of the first + * component, then everything of the second component, and so on. + */ + std::vector lexicographic_numbering; + /** * Stores the number of quadrature points in * @p dim dimensions for a cell. @@ -181,16 +190,6 @@ namespace internal }; - /** - * Extracts the mapping between the actual numbering inside a scalar - * element and lexicographic numbering as required by the evaluation - * routines. - */ - template - std::vector - get_lexicographic_numbering_inverse (const FiniteElement &fe); - - // ------------------------------------------ inline functions 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 d11c8654e0..eb8e2bda75 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 @@ -34,43 +34,6 @@ namespace internal namespace MatrixFreeFunctions { - // helper function - template - std::vector - get_lexicographic_numbering_inverse(const FiniteElement &fe) - { - Assert(fe.n_components() == 1, - ExcMessage("Expected a scalar element")); - - const FE_Poly,dim,dim> *fe_poly = - dynamic_cast,dim,dim>*>(&fe); - - const FE_Poly >,dim,dim> *fe_poly_piece = - dynamic_cast >,dim,dim>*> (&fe); - - const FE_DGP *fe_dgp = dynamic_cast*>(&fe); - - std::vector lexicographic; - if (fe_poly != 0) - lexicographic = fe_poly->get_poly_space_numbering_inverse(); - else if (fe_poly_piece != 0) - lexicographic = fe_poly_piece->get_poly_space_numbering_inverse(); - else if (fe_dgp != 0) - { - lexicographic.resize(fe_dgp->dofs_per_cell); - for (unsigned int i=0; idofs_per_cell; ++i) - lexicographic[i] = i; - } - else - Assert(false, ExcNotImplemented()); - - return lexicographic; - } - - - // ----------------- actual ShapeInfo functions -------------------- template @@ -99,19 +62,67 @@ namespace internal const unsigned int n_dofs_1d = fe->degree+1, n_q_points_1d = quad.size(); - std::vector lexicographic (fe->dofs_per_cell); // renumber (this is necessary for FE_Q, for example, since there the // vertex DoFs come first, which is incompatible with the lexicographic // ordering necessary to apply tensor products efficiently) + std::vector scalar_lexicographic; { - lexicographic = get_lexicographic_numbering_inverse(*fe); + // find numbering to lexicographic + Assert(fe->n_components() == 1, + ExcMessage("Expected a scalar element")); + + const FE_Poly,dim,dim> *fe_poly = + dynamic_cast,dim,dim>*>(fe); + + const FE_Poly >,dim,dim> *fe_poly_piece = + dynamic_cast >,dim,dim>*> (fe); + + const FE_DGP *fe_dgp = dynamic_cast*>(fe); + + if (fe_poly != 0) + scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse(); + else if (fe_poly_piece != 0) + scalar_lexicographic = fe_poly_piece->get_poly_space_numbering_inverse(); + 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; + } + else + Assert(false, ExcNotImplemented()); + + // Finally store the renumbering into the member variable of this + // class + if (fe_in.n_components() == 1) + lexicographic_numbering = scalar_lexicographic; + else + { + // have more than one component, get the inverse + // permutation, invert it, sort the components one after one, + // and invert back + std::vector scalar_inv = + Utilities::invert_permutation(scalar_lexicographic); + std::vector lexicographic (fe_in.dofs_per_cell); + for (unsigned int comp=0; comp()) == 1. otherwise, need // other entry point (e.g. generating a 1D element by reading the // name, as done before r29356) - Assert(std::fabs(fe->shape_value(lexicographic[0], Point())-1) < 1e-13, + Assert(std::fabs(fe->shape_value(scalar_lexicographic[0], + Point())-1) < 1e-13, ExcInternalError()); } @@ -138,7 +149,7 @@ namespace internal { // need to reorder from hierarchical to lexicographic to get the // DoFs correct - const unsigned int my_i = lexicographic[i]; + const unsigned int my_i = scalar_lexicographic[i]; for (unsigned int q=0; q +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +std::ofstream logfile("output"); + + + +template > +class MatrixFreeTest +{ +public: + + MatrixFreeTest(const MatrixFree &data_in): + data (data_in) + {}; + + void + local_mass_operator (const MatrixFree &data, + VECTOR &dst, + const VECTOR &src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval (data); + const unsigned int n_q_points = fe_eval.n_q_points; + + for (unsigned int cell=cell_range.first; cell &data, + VECTOR &dst, + const VECTOR &src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval (data); + MatrixFreeOperators::CellwiseInverseMassMatrix mass_inv(fe_eval); + const unsigned int n_q_points = fe_eval.n_q_points; + AlignedVector > inverse_coefficients(n_q_points); + + for (unsigned int cell=cell_range.first; cell::local_mass_operator, + this, dst, src); + }; + + void apply_inverse (VECTOR &dst, + const VECTOR &src) const + { + dst = 0; + data.cell_loop (&MatrixFreeTest::local_inverse_mass_operator, + this, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void do_test (const DoFHandler &dof) +{ + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + MappingQ mapping(4); + + 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; + ConstraintMatrix constraints; + + mf_data.reinit (mapping, dof, constraints, quad, data); + } + + MatrixFreeTest mf (mf_data); + Vector in (dof.n_dofs()), inverse (dof.n_dofs()), reference(dof.n_dofs()); + + for (unsigned int i=0; i > solver(control); + solver.solve (mf, reference, in, PreconditionIdentity()); + + inverse -= reference; + const double diff_norm = inverse.linfty_norm() / reference.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_DGQ fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + + do_test (dof); + + if (dim == 2) + { + deallog.push("float"); + deallog.threshold_double(1.e-6); + do_test (dof); + deallog.threshold_double(5.e-11); + deallog.pop(); + } +} + + + +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>(); + test<2,4>(); + deallog.pop(); + deallog.push("3d"); + test<3,1>(); + test<3,2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/inverse_mass_01.output b/tests/matrix_free/inverse_mass_01.output new file mode 100644 index 0000000000..98c2aaf8e9 --- /dev/null +++ b/tests/matrix_free/inverse_mass_01.output @@ -0,0 +1,41 @@ + +DEAL:2d::Testing FE_DGQ<2>(1) +DEAL:2d:cg::Starting value 6.32 +DEAL:2d:cg::Convergence step 51 value 0 +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d:float::Testing FE_DGQ<2>(1) +DEAL:2d:float:cg::Starting value 6.15 +DEAL:2d:float:cg::Convergence step 66 value 0 +DEAL:2d:float::Norm of difference: 0 +DEAL:2d:float:: +DEAL:2d::Testing FE_DGQ<2>(2) +DEAL:2d:cg::Starting value 8.62 +DEAL:2d:cg::Convergence step 116 value 0 +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d:float::Testing FE_DGQ<2>(2) +DEAL:2d:float:cg::Starting value 9.05 +DEAL:2d:float:cg::Convergence step 158 value 0 +DEAL:2d:float::Norm of difference: 0 +DEAL:2d:float:: +DEAL:2d::Testing FE_DGQ<2>(4) +DEAL:2d:cg::Starting value 14.6 +DEAL:2d:cg::Convergence step 354 value 0 +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d:float::Testing FE_DGQ<2>(4) +DEAL:2d:float:cg::Starting value 14.7 +DEAL:2d:float:cg::Convergence step 437 value 0 +DEAL:2d:float::Norm of difference: 0 +DEAL:2d:float:: +DEAL:3d::Testing FE_DGQ<3>(1) +DEAL:3d:cg::Starting value 13.7 +DEAL:3d:cg::Convergence step 131 value 0 +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_DGQ<3>(2) +DEAL:3d:cg::Starting value 13.5 +DEAL:3d:cg::Convergence step 255 value 0 +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/matrix_free/inverse_mass_02.cc b/tests/matrix_free/inverse_mass_02.cc new file mode 100644 index 0000000000..e3efb75c44 --- /dev/null +++ b/tests/matrix_free/inverse_mass_02.cc @@ -0,0 +1,205 @@ +// --------------------------------------------------------------------- +// $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. +// +// --------------------------------------------------------------------- + + + +// Tests CellwiseInverseMassMatrix on vector DG elements, otherwise the same +// as inverse_mass_01.cc + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +std::ofstream logfile("output"); + + + +template > +class MatrixFreeTest +{ +public: + + MatrixFreeTest(const MatrixFree &data_in): + data (data_in) + {}; + + void + local_mass_operator (const MatrixFree &data, + VECTOR &dst, + const VECTOR &src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval (data); + const unsigned int n_q_points = fe_eval.n_q_points; + + for (unsigned int cell=cell_range.first; cell &data, + VECTOR &dst, + const VECTOR &src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval (data); + MatrixFreeOperators::CellwiseInverseMassMatrix mass_inv(fe_eval); + const unsigned int n_q_points = fe_eval.n_q_points; + AlignedVector > inverse_coefficients(n_q_points); + + for (unsigned int cell=cell_range.first; cell::local_mass_operator, + this, dst, src); + }; + + void apply_inverse (VECTOR &dst, + const VECTOR &src) const + { + dst = 0; + data.cell_loop (&MatrixFreeTest::local_inverse_mass_operator, + this, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void do_test (const DoFHandler &dof) +{ + + 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; + ConstraintMatrix constraints; + + mf_data.reinit (dof, constraints, quad, data); + } + + MatrixFreeTest mf (mf_data); + Vector in (dof.n_dofs()), inverse (dof.n_dofs()), reference(dof.n_dofs()); + + for (unsigned int i=0; i > solver(control); + solver.solve (mf, reference, in, PreconditionIdentity()); + + inverse -= reference; + const double diff_norm = inverse.linfty_norm() / reference.linfty_norm(); + + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} + + + +template +void test () +{ + Triangulation tria; + GridGenerator::hyper_cube (tria, -1, 1); + tria.refine_global(1); + 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(); + + FESystem fe (FE_DGQ(fe_degree), dim); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + + do_test (dof); +} + + + +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>(); + test<2,4>(); + deallog.pop(); + deallog.push("3d"); + test<3,1>(); + test<3,2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/inverse_mass_02.output b/tests/matrix_free/inverse_mass_02.output new file mode 100644 index 0000000000..67b9f69785 --- /dev/null +++ b/tests/matrix_free/inverse_mass_02.output @@ -0,0 +1,26 @@ + +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^2] +DEAL:2d:cg::Starting value 7.90 +DEAL:2d:cg::Convergence step 7 value 0 +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2] +DEAL:2d:cg::Starting value 11.6 +DEAL:2d:cg::Convergence step 15 value 0 +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^2] +DEAL:2d:cg::Starting value 19.1 +DEAL:2d:cg::Convergence step 46 value 0 +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d:cg::Starting value 25.0 +DEAL:3d:cg::Convergence step 10 value 0 +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d:cg::Starting value 24.3 +DEAL:3d:cg::Convergence step 32 value 0 +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/matrix_free/inverse_mass_03.cc b/tests/matrix_free/inverse_mass_03.cc new file mode 100644 index 0000000000..02abac7109 --- /dev/null +++ b/tests/matrix_free/inverse_mass_03.cc @@ -0,0 +1,220 @@ +// --------------------------------------------------------------------- +// $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. +// +// --------------------------------------------------------------------- + + + +// Tests CellwiseInverseMassMatrix on vector DG elements, similar test as +// inverse_mass_02 but using different coefficients on different components + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +std::ofstream logfile("output"); + + + +template > +class MatrixFreeTest +{ +public: + + MatrixFreeTest(const MatrixFree &data_in): + data (data_in) + {}; + + void + local_mass_operator (const MatrixFree &data, + VECTOR &dst, + const VECTOR &src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval (data); + const unsigned int n_q_points = fe_eval.n_q_points; + + for (unsigned int cell=cell_range.first; cell > val = fe_eval.get_value(q); + val[0] *= make_vectorized_array(0.8314159); + val[1] *= make_vectorized_array(2.3); + val[2] *= make_vectorized_array(1.98); + fe_eval.submit_value (val, q); + } + fe_eval.integrate (true, false); + fe_eval.distribute_local_to_global (dst); + } + } + + void + local_inverse_mass_operator (const MatrixFree &data, + VECTOR &dst, + const VECTOR &src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval (data); + MatrixFreeOperators::CellwiseInverseMassMatrix mass_inv(fe_eval); + const unsigned int n_q_points = fe_eval.n_q_points; + AlignedVector > inverse_coefficients(3*n_q_points); + + for (unsigned int cell=cell_range.first; cell::local_mass_operator, + this, dst, src); + }; + + void apply_inverse (VECTOR &dst, + const VECTOR &src) const + { + dst = 0; + data.cell_loop (&MatrixFreeTest::local_inverse_mass_operator, + this, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void do_test (const DoFHandler &dof) +{ + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + MappingQ mapping(4); + 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; + ConstraintMatrix constraints; + + mf_data.reinit (mapping, dof, constraints, quad, data); + } + + MatrixFreeTest mf (mf_data); + Vector in (dof.n_dofs()), inverse (dof.n_dofs()), reference(dof.n_dofs()); + + for (unsigned int i=0; i > solver(control); + solver.solve (mf, reference, in, PreconditionIdentity()); + + inverse -= reference; + const double diff_norm = inverse.linfty_norm() / reference.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(); + + FESystem fe (FE_DGQ(fe_degree),3); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + + do_test (dof); +} + + + +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>(); + test<2,4>(); + deallog.pop(); + deallog.push("3d"); + test<3,1>(); + test<3,2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/inverse_mass_03.output b/tests/matrix_free/inverse_mass_03.output new file mode 100644 index 0000000000..35799127ea --- /dev/null +++ b/tests/matrix_free/inverse_mass_03.output @@ -0,0 +1,26 @@ + +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^3] +DEAL:2d:cg::Starting value 10.6 +DEAL:2d:cg::Convergence step 108 value 0 +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^3] +DEAL:2d:cg::Starting value 15.1 +DEAL:2d:cg::Convergence step 252 value 0 +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^3] +DEAL:2d:cg::Starting value 25.4 +DEAL:2d:cg::Convergence step 699 value 0 +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3] +DEAL:3d:cg::Starting value 23.7 +DEAL:3d:cg::Convergence step 412 value 0 +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3] +DEAL:3d:cg::Starting value 23.7 +DEAL:3d:cg::Convergence step 1019 value 0 +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_14.cc b/tests/matrix_free/matrix_vector_14.cc index 7f0bd5ac0a..80126fc3ea 100644 --- a/tests/matrix_free/matrix_vector_14.cc +++ b/tests/matrix_free/matrix_vector_14.cc @@ -104,16 +104,10 @@ private: template void do_test (const DoFHandler &dof, - const ConstraintMatrix &constraints, - const unsigned int parallel_option = 0) + const ConstraintMatrix &constraints) { deallog << "Testing " << dof.get_fe().get_name() << std::endl; - if (parallel_option > 0) - deallog << "Parallel option: " << parallel_option << std::endl; - //std::cout << "Number of cells: " << dof.get_tria().n_active_cells() << std::endl; - //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; - //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl; MatrixFree mf_data; { diff --git a/tests/matrix_free/matrix_vector_15.cc b/tests/matrix_free/matrix_vector_15.cc new file mode 100644 index 0000000000..0c8f9aa4c6 --- /dev/null +++ b/tests/matrix_free/matrix_vector_15.cc @@ -0,0 +1,242 @@ +// --------------------------------------------------------------------- +// $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 on +// the fly from an FEValues implementation instead of data cached in +// MatrixFree. + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +std::ofstream logfile("output"); + + + +template > +class MatrixFreeTest +{ +public: + MatrixFreeTest(const DoFHandler &dof_handler, + const ConstraintMatrix &constraints) + : + dof_handler (dof_handler), + constraints (constraints) + {} + + void vmult (VECTOR &dst, + const VECTOR &src) const + { + MappingFEEvaluation mapping(QGauss<1>(fe_degree+1), + update_gradients | update_values | + update_JxW_values); + VECTOR src_cpy = src; + constraints.distribute(src_cpy); + FEEvaluation fe_eval(mapping, dof_handler); + dst = 0; + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for ( ; cell != endc; ++cell) + { + mapping.reinit(cell); + fe_eval.read_dof_values (src_cpy); + fe_eval.evaluate (true, true, false); + for (unsigned int q=0; q &dof_handler; + const ConstraintMatrix &constraints; +}; + + + +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; + //std::cout << "Number of cells: " << dof.get_tria().n_active_cells() << std::endl; + //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; + //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl; + + MatrixFreeTest mf (dof, constraints); + 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 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); + + if (dim == 2) + { + deallog.push("float"); + deallog.threshold_double(1.e-6); + do_test (dof, constraints); + deallog.threshold_double(5.e-11); + deallog.pop(); + } +} + + + +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>(); + test<2,4>(); + deallog.pop(); + deallog.push("3d"); + test<3,1>(); + test<3,2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/matrix_vector_15.output b/tests/matrix_free/matrix_vector_15.output new file mode 100644 index 0000000000..0cd94749bd --- /dev/null +++ b/tests/matrix_free/matrix_vector_15.output @@ -0,0 +1,25 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d:float::Testing FE_Q<2>(1) +DEAL:2d:float::Norm of difference: 0 +DEAL:2d:float:: +DEAL:2d::Testing FE_Q<2>(2) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d:float::Testing FE_Q<2>(2) +DEAL:2d:float::Norm of difference: 0 +DEAL:2d:float:: +DEAL:2d::Testing FE_Q<2>(4) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d:float::Testing FE_Q<2>(4) +DEAL:2d:float::Norm of difference: 0 +DEAL:2d:float:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_16.cc b/tests/matrix_free/matrix_vector_16.cc new file mode 100644 index 0000000000..5b2e37ff55 --- /dev/null +++ b/tests/matrix_free/matrix_vector_16.cc @@ -0,0 +1,237 @@ +// --------------------------------------------------------------------- +// $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. +// +// --------------------------------------------------------------------- + + + +// evaluating the geometry on the fly as in matrix_vector_15, but using +// vector-valued finite elements. + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +std::ofstream logfile("output"); + + + +template > +class MatrixFreeTest +{ +public: + MatrixFreeTest(const DoFHandler &dof_handler, + const ConstraintMatrix &constraints) + : + dof_handler (dof_handler), + constraints (constraints) + {} + + void vmult (VECTOR &dst, + const VECTOR &src) const + { + MappingFEEvaluation mapping(QGauss<1>(fe_degree+1), + update_gradients | update_values | + update_JxW_values); + VECTOR src_cpy = src; + constraints.distribute(src_cpy); + FEEvaluation fe_eval(mapping, dof_handler); + dst = 0; + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for ( ; cell != endc; ++cell) + { + mapping.reinit(cell); + fe_eval.read_dof_values (src_cpy); + fe_eval.evaluate (true, true, false); + for (unsigned int q=0; q(10.) * + fe_eval.get_value(q),q); + fe_eval.submit_gradient (fe_eval.get_gradient(q),q); + } + fe_eval.integrate (true,true); + fe_eval.distribute_local_to_global (dst); + } + constraints.condense(dst); + }; + +private: + const DoFHandler &dof_handler; + const ConstraintMatrix &constraints; +}; + + + +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; + //std::cout << "Number of cells: " << dof.get_tria().n_active_cells() << std::endl; + //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; + //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl; + + MatrixFreeTest mf (dof, constraints); + 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(); + + FESystem fe(FE_Q(fe_degree), dim); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction(dim), + constraints); + constraints.close(); + + 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>(); + test<2,4>(); + deallog.pop(); + deallog.push("3d"); + test<3,1>(); + test<3,2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/matrix_vector_16.output b/tests/matrix_free/matrix_vector_16.output new file mode 100644 index 0000000000..b1dbf5d0cd --- /dev/null +++ b/tests/matrix_free/matrix_vector_16.output @@ -0,0 +1,16 @@ + +DEAL:2d::Testing FESystem<2>[FE_Q<2>(1)^2] +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_Q<2>(2)^2] +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FESystem<2>[FE_Q<2>(4)^2] +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FESystem<3>[FE_Q<3>(1)^3] +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FESystem<3>[FE_Q<3>(2)^3] +DEAL:3d::Norm of difference: 0 +DEAL:3d:: -- 2.39.5