From: kormann Date: Wed, 13 Jun 2012 19:52:16 +0000 (+0000) Subject: Fig bug and augment FEEValuation. Rename MatrixFree::get_quadrature_formula into... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fe962b376bb331e812c26401dc63dad9654b3e80;p=dealii-svn.git Fig bug and augment FEEValuation. Rename MatrixFree::get_quadrature_formula into MatrixFree::get_quadrature. Fix typos in tutorial programs. Cleanup in matrix free class for hp case. git-svn-id: https://svn.dealii.org/trunk@25626 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-37/doc/intro.dox b/deal.II/examples/step-37/doc/intro.dox index fd54c93359..19e8a07e9f 100644 --- a/deal.II/examples/step-37/doc/intro.dox +++ b/deal.II/examples/step-37/doc/intro.dox @@ -8,7 +8,7 @@ The algorithm for the matrix-vector product is based on the article A generic interface for parallel cell-based finite element operator application by Martin Kronbichler and Katharina Kormann, Computers -and Fluids, 2012, and the paper "Parallel finite element operator +and Fluids 63:135–147, 2012, and the paper "Parallel finite element operator application: Graph partitioning and coloring" by Katharina Kormann and Martin Kronbichler in: Proceedings of the 7th IEEE International Conference on e-Science, 2011. diff --git a/deal.II/examples/step-37/doc/results.dox b/deal.II/examples/step-37/doc/results.dox index 8c269d575b..ab9a9b0648 100644 --- a/deal.II/examples/step-37/doc/results.dox +++ b/deal.II/examples/step-37/doc/results.dox @@ -204,8 +204,9 @@ factor 12. This is the threshold when the cache in the processor can no longer hold all data necessary for the matrix-vector products and all matrix elements must be fetched from main memory. The second deviation is the times for the matrix-free solve which increase by less than a factor 8. This is because of -more parallelism from more cells, exploited by the (involved) dynamic tasks -scheduling approach CellFEOperator. Note that about 30% of the time in the +more parallelism from more cells, exploited by the (involved) dynamic task +scheduling approach taken in the cell loop of the MatrixFree class. Note +that about 30% of the time in the matrix-free solver is spent on restriction and prolongation, which use sparse matrices. So the speedup could be even better if all parts where done efficiently. diff --git a/deal.II/examples/step-37/step-37.cc b/deal.II/examples/step-37/step-37.cc index bc24a553e6..9f6e4b8987 100644 --- a/deal.II/examples/step-37/step-37.cc +++ b/deal.II/examples/step-37/step-37.cc @@ -1,8 +1,9 @@ -/* Author: Katharina Kormann, Martin Kronbichler, Uppsala University, 2009-2011 */ +/* $Id$ */ +/* Author: Katharina Kormann, Martin Kronbichler, Uppsala University, 2009-2012 */ -/* $Id$ */ +/* $Id$ */ /* */ -/* Copyright (C) 2009, 2010, 2011, 2012 by the deal.II authors */ +/* Copyright (C) 2009, 2010, 2011, 2012 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -189,7 +190,7 @@ namespace Step37 // writing to vectors), and hence SIMD can // be used to process several cells at // once. In all what follows, you can think - // of an AlignedVector to hold data from + // of a VectorizedArray to hold data from // several cells. For example, we evaluate // the coefficient shown here not on a // simple point as usually done, but we @@ -206,7 +207,7 @@ namespace Step37 // contains the x-coordinate for the first // point and the second point. You may // access the coordinates individually - // using e.g. p[0][j], j=1,2, + // using e.g. p[0][j], j=0,1, // but it is recommended to define // operations on a VectorizedArray as much // as possible in order to make use of diff --git a/deal.II/examples/step-48/doc/intro.dox b/deal.II/examples/step-48/doc/intro.dox index b1ceccc553..04d9ee7e4e 100644 --- a/deal.II/examples/step-48/doc/intro.dox +++ b/deal.II/examples/step-48/doc/intro.dox @@ -7,7 +7,7 @@ The algorithm for the matrix-vector product is based on the article A generic interface for parallel cell-based finite element operator application by Martin Kronbichler and Katharina Kormann, Computers -and Fluids, 2012, and the paper "Parallel finite element operator +and Fluids 63:135–147, 2012, and the paper "Parallel finite element operator application: Graph partitioning and coloring" by Katharina Kormann and Martin Kronbichler in: Proceedings of the 7th IEEE International Conference on e-Science, 2011. @@ -23,10 +23,12 @@ MatrixFree class handles constraints and how it can be parallelized over distributed nodes. Finally, we will use an explicit time-stepping method to solve the problem and introduce Gauss-Lobatto finite elements that are very convenient in this case since they have -a diagonally, and thus trivially invertable, mass matrix. Moreover, +a diagonal, and thus trivially invertible, mass matrix. Moreover, this type of elements clusters the nodes towards the element -boundaries which is why it has good properties for high-order -discretization methods. +boundaries which is why they have good properties for high-order +discretization methods. Indeed, the condition number of standard FE_Q +elements with equidistant nodes grows exponentially with the degree, +which destroys any benefit for orders of about five and higher.

Problem statement and discretization

diff --git a/deal.II/examples/step-48/doc/results.dox b/deal.II/examples/step-48/doc/results.dox index 34e7cf6dbd..52308a1754 100644 --- a/deal.II/examples/step-48/doc/results.dox +++ b/deal.II/examples/step-48/doc/results.dox @@ -2,7 +2,7 @@

Comparison with a sparse matrix

-In order to demonstrate the gain in using the CellFEOperator class instead of +In order to demonstrate the gain in using the MatrixFree class instead of the standard deal.II assembly routines for evaluating the information from old time steps, we study a simple serial run of the code on a nonadaptive mesh. Since much time is spent on evaluating the sine function, we @@ -59,7 +59,7 @@ following table. -It is apparent that the CellFEOperator outperforms the standard assembly +It is apparent that the matrix-free code outperforms the standard assembly routines in deal.II by far. In 3D and for fourth order elements, one operator application is also almost ten times as fast as a sparse matrix-vector product. diff --git a/deal.II/examples/step-48/step-48.cc b/deal.II/examples/step-48/step-48.cc index ddd3174156..451281dc28 100644 --- a/deal.II/examples/step-48/step-48.cc +++ b/deal.II/examples/step-48/step-48.cc @@ -1,9 +1,9 @@ /* $Id$ */ -/* Author: Katharina Kormann, Martin Kronbichler, Uppsala University, 2011 */ +/* Author: Katharina Kormann, Martin Kronbichler, Uppsala University, 2011-2012 */ /* $Id$ */ /* */ -/* Copyright (C) 2011, 2012 by the deal.II authors */ +/* Copyright (C) 2011, 2012 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -96,8 +96,6 @@ namespace Step48 class SineGordonOperation { public: - typedef VectorizedArray vector_t; - SineGordonOperation(const MatrixFree &data_in, const double time_step); @@ -105,12 +103,12 @@ namespace Step48 const std::vector*> &src) const; private: - const MatrixFree &data; - const vector_t delta_t_sqr; + const MatrixFree &data; + const VectorizedArray delta_t_sqr; parallel::distributed::Vector inv_mass_matrix; - void local_apply (const MatrixFree &data, - parallel::distributed::Vector &dst, + void local_apply (const MatrixFree &data, + parallel::distributed::Vector &dst, const std::vector*>&src, const std::pair &cell_range) const; }; @@ -141,12 +139,12 @@ namespace Step48 template SineGordonOperation:: SineGordonOperation(const MatrixFree &data_in, - const double time_step) + const double time_step) : data(data_in), delta_t_sqr(make_vectorized_array(time_step*time_step)) { - vector_t one = make_vectorized_array (1.); + VectorizedArray one = make_vectorized_array (1.); data.initialize_dof_vector (inv_mass_matrix); @@ -242,8 +240,8 @@ namespace Step48 for (unsigned int q=0; q current_value = current.get_value(q); + const VectorizedArray old_value = old.get_value(q); current.submit_value (2.*current_value - old_value - delta_t_sqr * std::sin(current_value),q); @@ -269,9 +267,9 @@ namespace Step48 // the cell loop is implemented in the cell // finite element operator class. On each cell // it applies the routine defined as the - // operator () method of the + // local_apply() method of the // class SineGordonOperation, - // i.e., *this. One could also + // i.e., this. One could also // provide a function with the same signature // that is not part of a class. template @@ -468,7 +466,7 @@ namespace Step48 // to use shared-memory parallelization (hence // one would use multithreading for intra-node // parallelism and not MPI; note that we here - // choose the standard option &mdash if we + // choose the standard option — if we // wanted to disable shared memory // parallelization, we would choose @p // none). Finally, three solution vectors are diff --git a/deal.II/include/deal.II/matrix_free/dof_info.templates.h b/deal.II/include/deal.II/matrix_free/dof_info.templates.h index bba6378016..8c680c9646 100644 --- a/deal.II/include/deal.II/matrix_free/dof_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/dof_info.templates.h @@ -566,6 +566,8 @@ namespace MatrixFreeFunctions std::vector &renumbering, std::vector &irregular_cells) { + if (max_fe_index < 2) + return; const unsigned int n_active_cells = size_info.n_active_cells; const unsigned int vectorization_length = size_info.vectorization_length; irregular_cells.resize (0); 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 7e8a42ed09..580042ef36 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -38,6 +38,11 @@ namespace parallel namespace internal { DeclException0 (ExcAccessToUninitializedField); + + template + void do_evaluate (FEEval &, const bool, const bool, const bool); + template + void do_integrate (FEEval &, const bool, const bool); } @@ -74,13 +79,15 @@ namespace internal * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ template + int n_components_, typename Number> class FEEvaluationBase { public: - typedef Tensor<1,n_components,VectorizedArray > value_type; - typedef Tensor<1,n_components,Tensor<1,dim,VectorizedArray > > gradient_type; + typedef Number number_type; + typedef Tensor<1,n_components_,VectorizedArray > value_type; + typedef Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > gradient_type; static const unsigned int dimension = dim; + static const unsigned int n_components = n_components_; static const unsigned int dofs_per_cell = dofs_per_cell_; static const unsigned int n_q_points = n_q_points_; @@ -175,7 +182,7 @@ public: */ template void read_dof_values (const std::vector &src, - const unsigned int first_index=0); + const unsigned int first_index=0); /** * Reads data from several vectors. Same as @@ -184,21 +191,7 @@ public: */ template void read_dof_values (const std::vector &src, - const unsigned int first_index=0); - - /** - * For a collection of several vector @p src, - * read out the values on the degrees of - * freedom of the current cell for @p - * n_components (template argument), and store - * them internally. Similar functionality as - * the function - * ConstraintMatrix::read_dof_values. Note - * that if vectorization is enabled, the DoF - * values for several cells are set. - */ - template - void read_dof_values (const VectorType * src_data[]); + const unsigned int first_index=0); /** * For the vector @p src, read out the values @@ -248,20 +241,6 @@ public: void read_dof_values_plain (const std::vector &src, const unsigned int first_index=0); - /** - * For a collection of several vector @p src, - * read out the values on the degrees of - * freedom of the current cell for @p - * n_components (template argument), and store - * them internally. Similar functionality as - * the function - * DoFAccessor::read_dof_values. Note - * that if vectorization is enabled, the DoF - * values for several cells are set. - */ - template - void read_dof_values_plain (const VectorType * src_data[]); - /** * Takes the values stored internally on dof * values of the current cell and sums them @@ -304,23 +283,6 @@ public: void distribute_local_to_global (std::vector &dst, const unsigned int first_index=0) const; - /** - * Takes the values stored internally on dof - * values of the current cell for a - * vector-valued problem consisting of @p - * n_components (template argument) and sums - * them into the collection of vectors 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. - */ - template - void distribute_local_to_global (VectorType * dst_data[]) const; - /** * Takes the values stored internally on dof * values of the current cell and sums them @@ -363,23 +325,6 @@ public: void set_dof_values (std::vector &dst, const unsigned int first_index=0) const; - /** - * Takes the values stored internally on dof - * values of the current cell for a - * vector-valued problem consisting of @p - * n_components (template argument) and sums - * them into the collection of vectors 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. - */ - template - void set_dof_values (VectorType * dst_data[]) const; - //@} /** @@ -518,7 +463,7 @@ public: * (n_components == 1) and for the * vector-valued case (n_components == dim). */ - Tensor<1,n_components,Tensor<2,dim,VectorizedArray > > + Tensor<1,n_components_,Tensor<2,dim,VectorizedArray > > get_hessian (const unsigned int q_point) const; /** @@ -569,6 +514,132 @@ public: //@} + /** + * @name 4: Access to internal data + */ + //@{ + /** + * Returns a read-only pointer to the first + * field of function values on quadrature + * points. First come the function values on + * all quadrature points for the first + * component, then all values for the second + * component, and so on. This is related to + * the internal data structures used in this + * class. The raw data after a call to @p + * evaluate only contains unit cell + * operations, so possible transformations, + * quadrature weights etc. must be applied + * manually. In general, it is safer to use + * the get_value() function instead, which + * does all the transformation internally. + */ + const VectorizedArray * begin_values () const; + + /** + * Returns a read and write pointer to the + * first field of function values on + * quadrature points. First come the function + * values on all quadrature points for the + * first component, then all values for the + * second component, and so on. This is + * related to the internal data structures + * used in this class. The raw data after a + * call to @p evaluate only contains unit + * cell operations, so possible + * transformations, quadrature weights + * etc. must be applied manually. In general, + * it is safer to use the get_value() function + * instead, which does all the transformation + * internally. + */ + VectorizedArray * begin_values (); + + /** + * Returns a read-only pointer to the first + * field of function gradients on quadrature + * points. First comes the x-component of the + * gradient for the first component on all + * quadrature points, then the y-component, + * and so on. Next comes the x-component of + * the second component, and so on. This is + * related to the internal data structures + * used in this class. The raw data after a + * call to @p evaluate only contains unit + * cell operations, so possible + * transformations, quadrature weights + * etc. must be applied manually. In general, + * it is safer to use the get_gradient() function + * instead, which does all the transformation + * internally. + */ + const VectorizedArray * begin_gradients () const; + + /** + * Returns a read and write pointer to the + * first field of function gradients on + * quadrature points. First comes the + * x-component of the gradient for the first + * component on all quadrature points, then + * the y-component, and so on. Next comes the + * x-component of the second component, and so + * on. This is related to the internal data + * structures used in this class. The raw data + * after a call to @p evaluate only + * contains unit cell operations, so possible + * transformations, quadrature weights + * etc. must be applied manually. In general, + * it is safer to use the get_gradient() + * function instead, which does all the + * transformation internally. + */ + VectorizedArray * begin_gradients (); + + /** + * Returns a read-only pointer to the first + * field of function hessians on quadrature + * points. First comes the xx-component of the + * hessian for the first component on all + * quadrature points, then the yy-component, + * zz-component in (3D), then the + * xy-component, and so on. Next comes the + * xx-component of the second component, and + * so on. This is related to the internal data + * structures used in this class. The raw data + * after a call to @p evaluate only + * contains unit cell operations, so possible + * transformations, quadrature weights + * etc. must be applied manually. In general, + * it is safer to use the get_laplacian() or + * get_hessian() functions instead, which does + * all the transformation internally. + */ + const VectorizedArray * begin_hessians () const; + + /** + * Returns a read and write pointer to the + * first field of function hessians on + * quadrature points. First comes the + * xx-component of the hessian for the first + * component on all quadrature points, then + * the yy-component, zz-component in (3D), + * then the xy-component, and so on. Next + * comes the xx-component of the second + * component, and so on. This is related to + * the internal data structures used in this + * class. The raw data after a call to @p + * evaluate only contains unit cell + * operations, so possible transformations, + * quadrature weights etc. must be applied + * manually. In general, it is safer to use + * the get_laplacian() or get_hessian() + * functions instead, which does all the + * transformation internally. + */ + VectorizedArray * begin_hessians (); + + //@} + protected: /** @@ -585,6 +656,33 @@ protected: const unsigned int fe_no = 0, const unsigned int quad_no = 0); + /** + * A unified function to read from and write + * into vectors based on the given template + * operation. It can perform the operation for + * @p read_dof_values, @p + * distribute_local_to_global, and @p + * set_dof_values. It performs the operation + * for several vectors at a time. + */ + template + void read_write_operation (const VectorOperation &operation, + VectorType *vectors[]) const; + + /** + * For a collection of several vector @p src, + * read out the values on the degrees of + * freedom of the current cell for @p + * n_components (template argument), and store + * them internally. Similar functionality as + * the function + * DoFAccessor::read_dof_values. Note + * that if vectorization is enabled, the DoF + * values for several cells are set. + */ + template + void read_dof_values_plain (const VectorType * src_data[]); + /** * Internal data fields that store the * values. Since all array lengths are known @@ -698,52 +796,6 @@ protected: */ const internal::MatrixFreeFunctions::ShapeInfo &data; - /** - * After a call to reinit(), stores the number - * of the cell we are currently working with. - */ - unsigned int cell; - - /** - * Stores the type of the cell we are - * currently working with after a call to - * reinit(). Valid values are @p cartesian, @p - * affine and @p general, which have different - * implications on how the Jacobian - * transformations are stored internally in - * MappingInfo. - */ - internal::MatrixFreeFunctions::CellType cell_type; - - /** - * The stride to access the correct data in - * MappingInfo. - */ - unsigned int cell_data_number; - - /** - * Stores whether the present cell chunk used - * in vectorization is not completely filled - * up with physical cells. E.g. if - * vectorization dictates that four cells - * should be worked with but only three - * physical cells are left, this flag will be - * set to true, otherwise to false. Mainly - * used for internal checking when reading - * from vectors or writing to vectors. - */ - bool at_irregular_cell; - - /** - * If the present cell chunk for vectorization - * is not completely filled up with data, this - * field stores how many physical cells are - * underlying. Is between 1 and - * VectorizedArray::n_array_elements-1 - * (inclusive). - */ - unsigned int n_irreg_components_filled; - /** * A pointer to the Cartesian Jacobian * information of the present cell. Only set @@ -797,13 +849,59 @@ protected: */ const Tensor<1,(dim>1?dim*(dim-1)/2:1),Tensor<1,dim,VectorizedArray > > * jacobian_grad_upper; + /** + * After a call to reinit(), stores the number + * of the cell we are currently working with. + */ + unsigned int cell; + + /** + * Stores the type of the cell we are + * currently working with after a call to + * reinit(). Valid values are @p cartesian, @p + * affine and @p general, which have different + * implications on how the Jacobian + * transformations are stored internally in + * MappingInfo. + */ + internal::MatrixFreeFunctions::CellType cell_type; + + /** + * The stride to access the correct data in + * MappingInfo. + */ + unsigned int cell_data_number; + + /** + * If the present cell chunk for vectorization + * is not completely filled up with data, this + * field stores how many physical cells are + * underlying. Is between 1 and + * VectorizedArray::n_array_elements-1 + * (inclusive). + */ + unsigned int n_irreg_components_filled; + + /** + * Stores whether the present cell chunk used + * in vectorization is not completely filled + * up with physical cells. E.g. if + * vectorization dictates that four cells + * should be worked with but only three + * physical cells are left, this flag will be + * set to true, otherwise to false. Mainly + * used for internal checking when reading + * from vectors or writing to vectors. + */ + bool at_irregular_cell; + /** * Debug information to track whether dof * values have been initialized before * accessed. Used to control exceptions when * uninitialized data is used. */ - bool dof_values_initialized; + bool dof_values_initialized; /** * Debug information to track whether values @@ -811,7 +909,7 @@ protected: * before accessed. Used to control exceptions * when uninitialized data is used. */ - bool values_quad_initialized; + bool values_quad_initialized; /** * Debug information to track whether @@ -820,7 +918,7 @@ protected: * control exceptions when uninitialized data * is used. */ - bool gradients_quad_initialized; + bool gradients_quad_initialized; /** * Debug information to track whether @@ -829,7 +927,7 @@ protected: * control exceptions when uninitialized data * is used. */ - bool hessians_quad_initialized; + bool hessians_quad_initialized; /** * Debug information to track whether values @@ -838,7 +936,7 @@ protected: * actually stared. Used to control exceptions * when uninitialized data is used. */ - bool values_quad_submitted; + bool values_quad_submitted; /** * Debug information to track whether @@ -848,7 +946,7 @@ protected: * control exceptions when uninitialized data * is used. */ - bool gradients_quad_submitted; + bool gradients_quad_submitted; }; @@ -862,17 +960,19 @@ protected: * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ template + int n_components_, typename Number> class FEEvaluationAccess : - public FEEvaluationBase + public FEEvaluationBase { public: - typedef Tensor<1,n_components,VectorizedArray > value_type; - typedef Tensor<1,n_components,Tensor<1,dim,VectorizedArray > > gradient_type; + typedef Number number_type; + typedef Tensor<1,n_components_,VectorizedArray > value_type; + typedef Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > gradient_type; static const unsigned int dimension = dim; + static const unsigned int n_components = n_components_; static const unsigned int dofs_per_cell = dofs_per_cell_; static const unsigned int n_q_points = n_q_points_; - typedef FEEvaluationBase BaseClass; protected: @@ -907,6 +1007,7 @@ class FEEvaluationAccess : public FEEvaluationBase { public: + typedef Number number_type; typedef VectorizedArray value_type; typedef Tensor<1,dim,VectorizedArray > gradient_type; static const unsigned int dimension = dim; @@ -1072,9 +1173,11 @@ class FEEvaluationAccess : public FEEvaluationBase { public: - typedef Tensor<1,dim,VectorizedArray > value_type; - typedef Tensor<2,dim,VectorizedArray > gradient_type; + typedef Number number_type; + typedef Tensor<1,dim,VectorizedArray > value_type; + typedef Tensor<2,dim,VectorizedArray > gradient_type; static const unsigned int dimension = dim; + static const unsigned int n_components = dim; static const unsigned int dofs_per_cell = dofs_per_cell_; static const unsigned int n_q_points = n_q_points_; typedef FEEvaluationBase BaseClass; @@ -1246,21 +1349,24 @@ protected: * * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ -template +template class FEEvaluationGeneral : public FEEvaluationAccess::value, Utilities::fixed_int_power::value, - n_components,Number> + n_components_,Number> { public: typedef FEEvaluationAccess::value, Utilities::fixed_int_power::value, - n_components, Number> BaseClass; + n_components_, Number> 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; static const unsigned int n_q_points = BaseClass::n_q_points; @@ -1288,8 +1394,9 @@ class FEEvaluationGeneral : * @p get_gradient() or @p get_laplacian * return useful information. */ - void evaluate (bool evaluate_val, bool evaluate_grad, - bool evaluate_hess=false); + void evaluate (const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_hess = false); /** * This function takes the values and/or @@ -1301,7 +1408,8 @@ class FEEvaluationGeneral : * integrate_grad are used to enable/disable * some of values or gradients. */ - void integrate (bool integrate_val, bool integrate_grad); + void integrate (const bool integrate_val, + const bool integrate_grad); /** * Returns the q-th quadrature point stored in @@ -1313,21 +1421,63 @@ class FEEvaluationGeneral : protected: /** - * Internal function that applies the shape - * function data of the tensor product in a - * given coordinate direction (first template + * Internal function that applies the function + * values of the tensor product in a given + * coordinate direction (first template + * argument), from polynomials to values on + * quadrature points (second flag set to true) + * or in an integration loop from values on + * quadrature points to values tested by + * different test function (second flag set to + * false), and if the result is to be added to + * previous content in the data fields or + * not. + */ + template + void apply_values (const VectorizedArray in [], + VectorizedArray out []); + + /** + * Internal function that applies the gradient + * operation of the tensor product in a given + * coordinate direction (first template * argument), from polynomials to values on * quadrature points (second flag set to true) * or in an integration loop from values on * quadrature points to values tested by * different test function (second flag set to * false), and if the result is to be added to - * some previous results or not. + * previous content in the data fields or + * not. + */ + template + void apply_gradients (const VectorizedArray in [], + VectorizedArray out []); + + /** + * Internal function that applies the second + * derivative operation (Hessian) of the + * tensor product in a given coordinate + * direction (first template argument), from + * polynomials to values on quadrature points + * (second flag set to true) or in an + * integration loop from values on quadrature + * points to values tested by different test + * function (second flag set to false), and if + * the result is to be added to previous + * content in the data fields or not. */ template - void apply_tensor_prod (const VectorizedArray * shape_data, - const VectorizedArray in [], - VectorizedArray out []); + void apply_hessians (const VectorizedArray in [], + VectorizedArray out []); + + /** + * Friend declaration. + */ + template friend void + internal::do_evaluate (FEEval &, const bool, const bool, const bool); + template friend void + internal::do_integrate (FEEval &, const bool, const bool); }; @@ -1373,15 +1523,18 @@ protected: * * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ -template +template class FEEvaluation : - public FEEvaluationGeneral + public FEEvaluationGeneral { public: - typedef FEEvaluationGeneral BaseClass; + typedef FEEvaluationGeneral BaseClass; + typedef Number number_type; typedef typename BaseClass::value_type value_type; typedef typename BaseClass::gradient_type gradient_type; + 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; static const unsigned int n_q_points = BaseClass::n_q_points; @@ -1411,8 +1564,9 @@ class FEEvaluation : * (unless these values have been set * manually). */ - void evaluate (bool evaluate_val, bool evaluate_grad, - bool evaluate_hess=false); + void evaluate (const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_hess = false); /** * This function takes the values and/or @@ -1424,9 +1578,11 @@ class FEEvaluation : * integrate_grad are used to enable/disable * some of values or gradients. */ - void integrate (bool integrate_val, bool integrate_grad); + void integrate (const bool integrate_val, + const bool integrate_grad); protected: + /** * Internal function that applies the function * values of the tensor product in a given @@ -1477,6 +1633,14 @@ protected: template void apply_hessians (const VectorizedArray in [], VectorizedArray out []); + + /** + * Friend declarations. + */ + template friend void + internal::do_evaluate (FEEval &, const bool, const bool, const bool); + template friend void + internal::do_integrate (FEEval &, const bool, const bool); }; @@ -1518,14 +1682,17 @@ protected: * * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ -template +template class FEEvaluationGL : - public FEEvaluation + public FEEvaluation { public: - typedef FEEvaluation BaseClass; + 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; static const unsigned int n_q_points = BaseClass::n_q_points; @@ -1555,8 +1722,9 @@ class FEEvaluationGL : * (unless these values have been set * manually). */ - void evaluate (bool evaluate_val, bool evaluate_grad, - bool evaluate_lapl=false); + void evaluate (const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl = false); /** * This function takes the values and/or @@ -1568,7 +1736,8 @@ class FEEvaluationGL : * integrate_grad are used to enable/disable * some of values or gradients. */ - void integrate (bool integrate_val, bool integrate_grad); + void integrate (const bool integrate_val, + const bool integrate_grad); protected: /** @@ -1596,15 +1765,15 @@ protected: #ifndef DOXYGEN -/*----------------------- FEEvaluationBase -------------------------------*/ +/*----------------------- FEEvaluationBase ----------------------------------*/ template + int n_components_, typename Number> inline -FEEvaluationBase:: -FEEvaluationBase (const MatrixFree &data_in, - const unsigned int fe_no_in, - const unsigned int quad_no_in) +FEEvaluationBase +::FEEvaluationBase (const MatrixFree &data_in, + const unsigned int fe_no_in, + const unsigned int quad_no_in) : quad_no (quad_no_in), n_fe_components (data_in.get_dof_info(fe_no_in).n_components), @@ -1619,11 +1788,6 @@ FEEvaluationBase (const MatrixFree &data_in, data (data_in.get_shape_info (fe_no_in, quad_no_in, active_fe_index, active_quad_index)), - cell (numbers::invalid_unsigned_int), - cell_type (internal::MatrixFreeFunctions::undefined), - cell_data_number (0), - at_irregular_cell (false), - n_irreg_components_filled (0), cartesian_data (0), jacobian (0), J_value (0), @@ -1631,7 +1795,12 @@ FEEvaluationBase (const MatrixFree &data_in, quadrature_weights[active_quad_index].begin()), quadrature_points (0), jacobian_grad (0), - jacobian_grad_upper(0) + jacobian_grad_upper(0), + cell (numbers::invalid_unsigned_int), + cell_type (internal::MatrixFreeFunctions::undefined), + cell_data_number (0), + n_irreg_components_filled (0), + at_irregular_cell (false) { Assert (matrix_info.indices_initialized() == true, ExcNotInitialized()); @@ -1654,11 +1823,11 @@ FEEvaluationBase (const MatrixFree &data_in, template + int n_components_, typename Number> inline void -FEEvaluationBase:: -reinit (const unsigned int cell_in) +FEEvaluationBase +::reinit (const unsigned int cell_in) { AssertIndexRange (cell_in, dof_info.row_starts.size()-1); AssertDimension (((dof_info.cell_active_fe_index.size() > 0) ? @@ -1732,11 +1901,11 @@ reinit (const unsigned int cell_in) template + int n_components_, typename Number> inline unsigned int -FEEvaluationBase:: -get_cell_data_number () const +FEEvaluationBase +::get_cell_data_number () const { Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); return cell_data_number; @@ -1745,11 +1914,11 @@ get_cell_data_number () const template + int n_components_, typename Number> inline internal::MatrixFreeFunctions::CellType -FEEvaluationBase:: -get_cell_type () const +FEEvaluationBase +::get_cell_type () const { Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); return cell_type; @@ -1839,82 +2008,148 @@ namespace internal "Use MatrixFree::initialize_dof_vector to get a " "compatible vector.")); } -} + // A class to use the same code to read from + // and write to vector + template + struct VectorReader + { + template + void process_dof (const unsigned int index, + VectorType &vec, + Number &res) const + { + res = vector_access (vec, index); + } + void pre_constraints (const Number &, + Number &res) const + { + res = Number(); + } -template -template -inline -void -FEEvaluationBase:: -read_dof_values (const VectorType &src) -{ - AssertDimension (n_components, n_fe_components); - // only need one component, but to avoid - // compiler warnings, use n_components copies - // here (but these will not be used) - const VectorType * src_data[n_components]; - for (unsigned int d=0; d + void process_constraint (const unsigned int index, + const Number weight, + VectorType &vec, + Number &res) const + { + res += weight * vector_access (vec, index); + } + void post_constraints (const Number &sum, + Number &write_pos) const + { + write_pos = sum; + } + void process_empty (Number &res) const + { + res = Number(); + } + }; -template -template -inline -void -FEEvaluationBase:: -read_dof_values (const std::vector &src, - const unsigned int first_index) -{ - AssertIndexRange (first_index, src.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= src.size()) : true), - ExcIndexRange (first_index + n_components, 0, src.size())); - const VectorType * src_data [n_components]; - for (unsigned int comp=0; comp + struct VectorDistributorLocalToGlobal + { + template + void process_dof (const unsigned int index, + VectorType &vec, + Number &res) const + { + vector_access (vec, index) += res; + } + void pre_constraints (const Number &input, + Number &res) const + { + res = input; + } + template + void process_constraint (const unsigned int index, + const Number weight, + VectorType &vec, + Number &res) const + { + vector_access (vec, index) += weight * res; + } + + void post_constraints (const Number &, + Number &) const + { + } + + void process_empty (Number &) const + { + } + }; + + + // A class to use the same code to read from + // and write to vector + template + struct VectorSetter + { + template + void process_dof (const unsigned int index, + VectorType &vec, + Number &res) const + { + vector_access (vec, index) = res; + } + + void pre_constraints (const Number &, + Number &) const + { + } + + template + void process_constraint (const unsigned int, + const Number, + VectorType&, + Number &) const + { + } + + void post_constraints (const Number &, + Number &) const + { + } + + void process_empty (Number &) const + { + } + }; -template -template -inline -void -FEEvaluationBase:: -read_dof_values (const std::vector &src, - const unsigned int first_index) -{ - AssertIndexRange (first_index, src.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= src.size()) : true), - ExcIndexRange (first_index + n_components, 0, src.size())); - const VectorType * src_data [n_components]; - for (unsigned int comp=0; comp -template + int n_components_, typename Number> +template inline void -FEEvaluationBase:: -read_dof_values (const VectorType * src[]) +FEEvaluationBase +::read_write_operation (const VectorOperation &operation, + VectorType *src[]) const { + // This functions processes all the functions + // read_dof_values, + // distribute_local_to_global, and + // set_dof_values with the same code. The + // distinction between these three cases is + // made by the input VectorOperation that + // either reads values from a vector and puts + // the data into the local data field or write + // local data into the vector. Certain + // operations are no-ops for the given use + // case. + Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); // loop over all local dofs. ind_local holds @@ -1938,9 +2173,10 @@ read_dof_values (const VectorType * src[]) VectorizedArray::n_array_elements * dofs_per_cell; for (unsigned int comp=0; comp(values_dofs[comp]); + local_data[comp] = + const_cast(reinterpret_cast(values_dofs[comp])); // standard case where there are sufficiently // many cells to fill all vectors @@ -1955,8 +2191,9 @@ read_dof_values (const VectorType * src[]) // run through values up to next constraint for (unsigned int j=0; jfirst; ++j) for (unsigned int comp=0; compfirst; dof_indices += indicators->first; @@ -1965,19 +2202,22 @@ read_dof_values (const VectorType * src[]) // according to constraints Number value [n_components]; for (unsigned int comp=0; compsecond); const Number * end_pool = matrix_info.constraint_pool_end(indicators->second); for ( ; data_val != end_pool; ++data_val, ++dof_indices) for (unsigned int comp=0; comp(n_local_dofs)); for (unsigned int j=0; j= n_irreg_components_filled) { for (unsigned int comp=0; compsecond); const Number * end_pool = matrix_info.constraint_pool_end(indicators->second); + for ( ; data_val != end_pool; ++data_val, ++dof_indices) for (unsigned int comp=0; comp::n_array_elements >= n_irreg_components_filled) { for (unsigned int comp=0; comp::n_array_elements >= n_irreg_components_filled) { for (unsigned int comp=0; comp::n_array_elements * n_components; - Number * local_src_number = reinterpret_cast(values_dofs[0]); + Number * local_data = + const_cast(reinterpret_cast(values_dofs[0])); if (at_irregular_cell == false) { // check whether there is any constraint on @@ -2105,33 +2351,35 @@ read_dof_values (const VectorType * src[]) { // run through values up to next constraint for (unsigned int j=0; jfirst; ++j) - local_src_number[ind_local+j] = - internal::vector_access (*src[0], dof_indices[j]); + operation.process_dof (dof_indices[j], *src[0], + local_data[ind_local+j]); ind_local += indicators->first; dof_indices += indicators->first; // constrained case: build the local value as // a linear combination of the global value // according to constraints - Number value = 0; + Number value; + operation.pre_constraints (local_data[ind_local], value); + const Number * data_val = matrix_info.constraint_pool_begin(indicators->second); const Number * end_pool = matrix_info.constraint_pool_end(indicators->second); + for ( ; data_val != end_pool; ++data_val, ++dof_indices) - value += - (internal::vector_access (*src[0], *dof_indices) * - (*data_val)); + operation.process_constraint (*dof_indices, *data_val, + *src[0], value); - local_src_number[ind_local] = value; + operation.post_constraints (value, local_data[ind_local]); ind_local++; } // get the dof values past the last // constraint for(; ind_local(n_local_dofs)); for (unsigned int j=0; j::n_array_elements >= n_irreg_components_filled) { - local_src_number[ind_local] = 0.; + operation.process_empty (local_data[ind_local]); ++ind_local; } } @@ -2179,20 +2427,24 @@ read_dof_values (const VectorType * src[]) // constrained case: build the local value as // a linear combination of the global value // according to constraint - Number value = 0; + Number value; + operation.pre_constraints (local_data[ind_local], value); + const Number * data_val = matrix_info.constraint_pool_begin(indicators->second); const Number * end_pool = matrix_info.constraint_pool_end(indicators->second); + for ( ; data_val != end_pool; ++data_val, ++dof_indices) - value += - internal::vector_access (*src[0], *dof_indices) * (*data_val); - local_src_number[ind_local] = value; + operation.process_constraint (*dof_indices, *data_val, + *src[0], value); + + operation.post_constraints (value, local_data[ind_local]); ind_local++; while (ind_local % VectorizedArray::n_array_elements >= n_irreg_components_filled) { - local_src_number[ind_local] = 0.; + operation.process_empty (local_data[ind_local]); ++ind_local; } } @@ -2204,18 +2456,41 @@ read_dof_values (const VectorType * src[]) // non-constrained case: copy the data from // the global vector, src, to the local one, // local_dst. - local_src_number[ind_local] = - internal::vector_access (*src[0], *dof_indices); + operation.process_dof (*dof_indices, *src[0], + local_data[ind_local]); ++ind_local; while (ind_local % VectorizedArray::n_array_elements >= n_irreg_components_filled) { - local_src_number[ind_local] = 0.; + operation.process_empty (local_data[ind_local]); ++ind_local; } } } } +} + + + +template +template +inline +void +FEEvaluationBase +::read_dof_values (const VectorType &src) +{ + AssertDimension (n_components_, n_fe_components); + + // only need one component, but to silent + // compiler warnings, use n_components copies + // here (but these will not be used) + VectorType * src_data[n_components]; + for (unsigned int d=0; d(&src); + + internal::VectorReader reader; + read_write_operation (reader, src_data); #ifdef DEBUG dof_values_initialized = true; @@ -2224,18 +2499,73 @@ read_dof_values (const VectorType * src[]) +template +template +inline +void +FEEvaluationBase +::read_dof_values (const std::vector &src, + const unsigned int first_index) +{ + AssertIndexRange (first_index, src.size()); + Assert (n_fe_components == 1, ExcNotImplemented()); + Assert ((n_fe_components == 1 ? + (first_index+n_components <= src.size()) : true), + ExcIndexRange (first_index + n_components_, 0, src.size())); + + VectorType * src_data [n_components]; + for (unsigned int comp=0; comp(&src[comp+first_index]); + + internal::VectorReader reader; + read_write_operation (reader, src_data); + +#ifdef DEBUG + dof_values_initialized = true; +#endif +} template + int n_components_, typename Number> template inline void -FEEvaluationBase:: -read_dof_values_plain (const VectorType &src) +FEEvaluationBase +::read_dof_values (const std::vector &src, + const unsigned int first_index) +{ + AssertIndexRange (first_index, src.size()); + Assert (n_fe_components == 1, ExcNotImplemented()); + Assert ((n_fe_components == 1 ? + (first_index+n_components <= src.size()) : true), + ExcIndexRange (first_index + n_components_, 0, src.size())); + + const VectorType * src_data [n_components]; + for (unsigned int comp=0; comp(src[comp+first_index]); + + internal::VectorReader reader; + read_write_operation (reader, src_data); + +#ifdef DEBUG + dof_values_initialized = true; +#endif +} + + + +template +template +inline +void +FEEvaluationBase +::read_dof_values_plain (const VectorType &src) { - AssertDimension (n_components, n_fe_components); + AssertDimension (n_components_, n_fe_components); // only need one component, but to avoid // compiler warnings, use n_components copies // here (but these will not be used) @@ -2248,19 +2578,19 @@ read_dof_values_plain (const VectorType &src) template + int n_components_, typename Number> template inline void -FEEvaluationBase:: -read_dof_values_plain (const std::vector &src, - const unsigned int first_index) +FEEvaluationBase +::read_dof_values_plain (const std::vector &src, + const unsigned int first_index) { AssertIndexRange (first_index, src.size()); Assert (n_fe_components == 1, ExcNotImplemented()); Assert ((n_fe_components == 1 ? (first_index+n_components <= src.size()) : true), - ExcIndexRange (first_index + n_components, 0, src.size())); + ExcIndexRange (first_index + n_components_, 0, src.size())); const VectorType * src_data [n_components]; for (unsigned int comp=0; comp &src, template + int n_components_, typename Number> template inline void -FEEvaluationBase:: -read_dof_values_plain (const std::vector &src, - const unsigned int first_index) +FEEvaluationBase +::read_dof_values_plain (const std::vector &src, + const unsigned int first_index) { AssertIndexRange (first_index, src.size()); Assert (n_fe_components == 1, ExcNotImplemented()); Assert ((n_fe_components == 1 ? (first_index+n_components <= src.size()) : true), - ExcIndexRange (first_index + n_components, 0, src.size())); + ExcIndexRange (first_index + n_components_, 0, src.size())); const VectorType * src_data [n_components]; for (unsigned int comp=0; comp &src, template + int n_components_, typename Number> template inline void -FEEvaluationBase:: -read_dof_values_plain (const VectorType * src[]) +FEEvaluationBase +::distribute_local_to_global (VectorType &dst) const { - Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - Assert (dof_info.store_plain_indices == true, ExcNotInitialized()); + AssertDimension (n_components_, n_fe_components); + Assert (dof_values_initialized==true, + internal::ExcAccessToUninitializedField()); - // loop over all local dofs. ind_local holds - // local number on cell, index iterates over - // the elements of index_local_to_global and - // dof_indices points to the global indices - // stored in index_local_to_global - const unsigned int * dof_indices = dof_info.begin_indices_plain(cell); + // only need one component, but to avoid + // compiler warnings, use n_components copies + // here (but these will not be used) + VectorType * dst_data [n_components]; + for (unsigned int d=0; d::n_array_elements * dofs_per_cell; - for (unsigned int comp=0; comp(values_dofs[comp]); + internal::VectorDistributorLocalToGlobal distributor; + read_write_operation (distributor, dst_data); +} - // standard case where there are sufficiently - // many cells to fill all vectors - if (at_irregular_cell == false) - { - for (unsigned int j=0; j 0, ExcInternalError()); - for(unsigned int ind_local=0; ind_local::n_array_elements >= n_irreg_components_filled) - { - for (unsigned int comp=0; comp::n_array_elements * n_components; - Number * local_src_number = reinterpret_cast(values_dofs[0]); - if (at_irregular_cell == false) - { - for (unsigned int j=0; j 0, ExcInternalError()); - for(unsigned int ind_local=0; ind_local::n_array_elements >= n_irreg_components_filled) - { - local_src_number[ind_local] = 0.; - ++ind_local; - } - } - } - } +template +template +inline +void +FEEvaluationBase +::distribute_local_to_global (std::vector &dst, + const unsigned int first_index) const +{ + AssertIndexRange (first_index, dst.size()); + Assert (n_fe_components == 1, ExcNotImplemented()); + Assert ((n_fe_components == 1 ? + (first_index+n_components <= dst.size()) : true), + ExcIndexRange (first_index + n_components_, 0, dst.size())); + Assert (dof_values_initialized==true, + internal::ExcAccessToUninitializedField()); -#ifdef DEBUG - dof_values_initialized = true; -#endif + VectorType * dst_data [n_components]; + for (unsigned int comp=0; comp distributor; + read_write_operation (distributor, dst_data); } template + int n_components_, typename Number> template inline void -FEEvaluationBase:: -distribute_local_to_global (VectorType &dst) const +FEEvaluationBase +::distribute_local_to_global (std::vector &dst, + const unsigned int first_index) const { - AssertDimension (n_components, n_fe_components); + AssertIndexRange (first_index, dst.size()); + Assert (n_fe_components == 1, ExcNotImplemented()); + Assert ((n_fe_components == 1 ? + (first_index+n_components <= dst.size()) : true), + ExcIndexRange (first_index + n_components_, 0, dst.size())); + Assert (dof_values_initialized==true, + internal::ExcAccessToUninitializedField()); + + VectorType * dst_data [n_components]; + for (unsigned int comp=0; comp distributor; + read_write_operation (distributor, dst_data); +} + + + +template +template +inline +void +FEEvaluationBase +::set_dof_values (VectorType &dst) const +{ + AssertDimension (n_components_, n_fe_components); + Assert (dof_values_initialized==true, + internal::ExcAccessToUninitializedField()); + // only need one component, but to avoid // compiler warnings, use n_components copies // here (but these will not be used) VectorType * dst_data [n_components]; for (unsigned int d=0; d setter; + read_write_operation (setter, dst_data); } template + int n_components_, typename Number> template inline void -FEEvaluationBase:: -distribute_local_to_global (std::vector &dst, - const unsigned int first_index) const +FEEvaluationBase +::set_dof_values (std::vector &dst, + const unsigned int first_index) const { AssertIndexRange (first_index, dst.size()); Assert (n_fe_components == 1, ExcNotImplemented()); Assert ((n_fe_components == 1 ? (first_index+n_components <= dst.size()) : true), - ExcIndexRange (first_index + n_components, 0, dst.size())); + ExcIndexRange (first_index + n_components_, 0, dst.size())); + + Assert (dof_values_initialized==true, + internal::ExcAccessToUninitializedField()); VectorType * dst_data [n_components]; for (unsigned int comp=0; comp setter; + read_write_operation (setter, dst_data); } template + int n_components_, typename Number> template inline void -FEEvaluationBase:: -distribute_local_to_global (std::vector &dst, - const unsigned int first_index) const +FEEvaluationBase +::set_dof_values (std::vector &dst, + const unsigned int first_index) const { AssertIndexRange (first_index, dst.size()); Assert (n_fe_components == 1, ExcNotImplemented()); Assert ((n_fe_components == 1 ? (first_index+n_components <= dst.size()) : true), - ExcIndexRange (first_index + n_components, 0, dst.size())); + ExcIndexRange (first_index + n_components_, 0, dst.size())); + + Assert (dof_values_initialized==true, + internal::ExcAccessToUninitializedField()); VectorType * dst_data [n_components]; for (unsigned int comp=0; comp setter; + read_write_operation (setter, dst_data); } template + int n_components_, typename Number> template inline void -FEEvaluationBase:: -distribute_local_to_global (VectorType * dst[]) const +FEEvaluationBase +::read_dof_values_plain (const VectorType * src[]) { + // this is different from the other three + // operations because we do not use + // constraints here, so this is a separate + // function. Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - Assert (dof_values_initialized==true, - internal::ExcAccessToUninitializedField()); + Assert (dof_info.store_plain_indices == true, ExcNotInitialized()); // loop over all local dofs. ind_local holds // local number on cell, index iterates over // the elements of index_local_to_global and // dof_indices points to the global indices // stored in index_local_to_global - const unsigned int * dof_indices = dof_info.begin_indices(cell); - const std::pair * indicators = - dof_info.begin_indicators(cell); - const std::pair * indicators_end = - dof_info.end_indicators(cell); - unsigned int ind_local = 0; + const unsigned int * dof_indices = dof_info.begin_indices_plain(cell); // scalar case (or case when all components // have the same degrees of freedom and sit on @@ -2508,111 +2811,45 @@ distribute_local_to_global (VectorType * dst[]) const const unsigned int n_local_dofs = VectorizedArray::n_array_elements * dofs_per_cell; for (unsigned int comp=0; comp(values_dofs[comp]); + local_src_number[comp] = reinterpret_cast(values_dofs[comp]); + + // standard case where there are sufficiently + // many cells to fill all vectors if (at_irregular_cell == false) { - // check whether there is no constraint at all - if (indicators != indicators_end) - { - // run from one constraint to the next - for ( ; indicators != indicators_end; ++indicators) - { - // distribute values up to the constraint - // (values not constrained) - for (unsigned int j=0; jfirst; ++j) - for (unsigned int comp=0; compfirst; - ind_local += indicators->first; - - // constrained case: build the local value as - // a linear combination of the global value - // according to constraint - const Number * data_val = - matrix_info.constraint_pool_begin(indicators->second); - const Number * end_pool = - matrix_info.constraint_pool_end(indicators->second); - for ( ; data_val != end_pool; ++data_val, ++dof_indices) - for (unsigned int comp=0; comp(n_local_dofs)); - for (unsigned int j=0; j 0, ExcInternalError()); - for ( ; indicators != indicators_end; ++indicators) + // non-standard case: need to fill in zeros + // for those components that are not present + // (a bit more expensive), but there is not + // more than one such cell + else { - for(unsigned int j=0; jfirst; ++j) + Assert (n_irreg_components_filled > 0, ExcInternalError()); + for(unsigned int ind_local=0; ind_local::n_array_elements == n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements-n_irreg_components_filled; - } - dof_indices += indicators->first; - - // constrained case: distribute according to - // the constraint - const Number * data_val = - matrix_info.constraint_pool_begin(indicators->second); - const Number * end_pool = - matrix_info.constraint_pool_end(indicators->second); - for ( ; data_val != end_pool; ++data_val, ++dof_indices) - { - for (unsigned int comp=0; comp::n_array_elements >= n_irreg_components_filled) + { + for (unsigned int comp=0; comp::n_array_elements == - n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements- - n_irreg_components_filled; - } - for(; ind_local::n_array_elements == - n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements-n_irreg_components_filled; } } else @@ -2624,408 +2861,144 @@ distribute_local_to_global (VectorType * dst[]) const // is ensured by the way MatrixFree reads // out the indices. { - internal::check_vector_compatibility (*dst[0], dof_info); - Assert (n_fe_components == n_components, ExcNotImplemented()); + internal::check_vector_compatibility (*src[0], dof_info); + Assert (n_fe_components == n_components_, ExcNotImplemented()); const unsigned int n_local_dofs = dofs_per_cell * VectorizedArray::n_array_elements * n_components; - const Number * local_dst_number = - reinterpret_cast(values_dofs[0]); + Number * local_src_number = reinterpret_cast(values_dofs[0]); if (at_irregular_cell == false) { - // check whether there is no constraint at all - if (indicators != indicators_end) - { - // run from one constraint to the next - for ( ; indicators != indicators_end; ++indicators) - { - // distribute values up to the constraint - // (values not constrained) - for (unsigned int j=0; jfirst; ++j) - internal::vector_access (*dst[0], dof_indices[j]) - += local_dst_number[ind_local+j]; - dof_indices += indicators->first; - ind_local += indicators->first; - - // constrained case: build the local value as - // a linear combination of the global value - // according to constraint - const Number * data_val = - matrix_info.constraint_pool_begin(indicators->second); - const Number * end_pool = - matrix_info.constraint_pool_end(indicators->second); - for ( ; data_val != end_pool; ++data_val, ++dof_indices) - internal::vector_access (*dst[0], *dof_indices) - += local_dst_number[ind_local] * (*data_val); - ++ind_local; - } - // distribute values after the last constraint - // (values not constrained) - for(; ind_local(n_local_dofs)); - for (unsigned int j=0; j 0, ExcInternalError()); - for ( ; indicators != indicators_end; ++indicators) + // non-standard case: need to fill in zeros + // for those components that are not present + // (a bit more expensive), but there is not + // more than one such cell + else { - for(unsigned int j=0; jfirst; ++j) + Assert (n_irreg_components_filled > 0, ExcInternalError()); + for(unsigned int ind_local=0; ind_local::n_array_elements == n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements-n_irreg_components_filled; - } - dof_indices += indicators->first; - - // constrained case: distribute according to - // the constraint - const Number * data_val = - matrix_info.constraint_pool_begin(indicators->second); - const Number * end_pool = - matrix_info.constraint_pool_end(indicators->second); - for ( ; data_val != end_pool; ++data_val, ++dof_indices) - { - internal::vector_access (*dst[0], *dof_indices) - += local_dst_number[ind_local] * (*data_val); + while (ind_local % VectorizedArray::n_array_elements >= n_irreg_components_filled) + { + local_src_number[ind_local] = 0.; + ++ind_local; + } } - ++ind_local; - if (ind_local % VectorizedArray::n_array_elements == n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements-n_irreg_components_filled; } - for(; ind_local::n_array_elements == n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements-n_irreg_components_filled; - } - Assert (dof_indices == dof_info.end_indices(cell), - ExcInternalError()); } + +#ifdef DEBUG + dof_values_initialized = true; +#endif } + +/*------------------------------ access to data fields ----------------------*/ + template -template inline -void +const VectorizedArray * FEEvaluationBase:: -set_dof_values (VectorType &dst) const +begin_values () const { - AssertDimension (n_components, n_fe_components); - // only need one component, but to avoid - // compiler warnings, use n_components copies - // here (but these will not be used) - VectorType * dst_data [n_components]; - for (unsigned int d=0; d -template inline -void +VectorizedArray * FEEvaluationBase:: -set_dof_values (std::vector &dst, - const unsigned int first_index) const +begin_values () { - AssertIndexRange (first_index, dst.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= dst.size()) : true), - ExcIndexRange (first_index + n_components, 0, dst.size())); - - VectorType * dst_data [n_components]; - for (unsigned int comp=0; comp -template inline -void +const VectorizedArray * FEEvaluationBase:: -set_dof_values (std::vector &dst, - const unsigned int first_index) const +begin_gradients () const { - AssertIndexRange (first_index, dst.size()); - Assert (n_fe_components == 1, ExcNotImplemented()); - Assert ((n_fe_components == 1 ? - (first_index+n_components <= dst.size()) : true), - ExcIndexRange (first_index + n_components, 0, dst.size())); - - VectorType * dst_data [n_components]; - for (unsigned int comp=0; comp -template inline -void +VectorizedArray * FEEvaluationBase:: -set_dof_values (VectorType * dst[]) const +begin_gradients () { - Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - Assert (dof_values_initialized==true, - internal::ExcAccessToUninitializedField()); +#ifdef DEBUG + gradients_quad_submitted = true; +#endif + return &gradients_quad[0][0][0]; +} - // loop over all local dofs. ind_local holds - // local number on cell, index iterates over - // the elements of index_local_to_global and - // glob_indices points to the global indices - // stored in index_local_to_global - const unsigned int * dof_indices = dof_info.begin_indices(cell); - const std::pair * indicators = - dof_info.begin_indicators(cell); - const std::pair * indicators_end = - dof_info.end_indicators(cell); - unsigned int ind_local = 0; - if (n_fe_components == 1) - { - const unsigned int n_local_dofs = - VectorizedArray::n_array_elements * dofs_per_cell; - for (unsigned int comp=0; compsize(), - dof_info.vector_partitioner->size()); - const Number * local_dst_number [n_components]; - for (unsigned int comp=0; comp(values_dofs[comp]); - if (at_irregular_cell == false) - { - // check whether there is no constraint at all - if (indicators != indicators_end) - { - // run from one constraint to the next - for ( ; indicators != indicators_end; ++indicators) - { - // distribute values up to the constraint - // (values not constrained) - for (unsigned int j=0; jfirst; ++j) - for (unsigned int comp=0; compfirst; - ind_local += indicators->first; - - // jump over constraints - const unsigned int row_length = - matrix_info.constraint_pool_end(indicators->second)- - matrix_info.constraint_pool_begin(indicators->second); - dof_indices += row_length; - ++ind_local; - } - // distribute values after the last constraint - // (values not constrained) - for(; ind_local +inline +const VectorizedArray * +FEEvaluationBase:: +begin_hessians () const +{ + Assert (hessians_quad_initialized, ExcNotInitialized()); + return &hessians_quad[0][0][0]; +} - // irregular case - Assert (n_irreg_components_filled > 0, ExcInternalError()); - for ( ; indicators != indicators_end; ++indicators) - { - for(unsigned int j=0; jfirst; ++j) - { - // non-constrained case - for (unsigned int comp=0; comp::n_array_elements == n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements-n_irreg_components_filled; - } - dof_indices += indicators->first; - - // jump over constraint - const unsigned int row_length = - matrix_info.constraint_pool_end(indicators->second)- - matrix_info.constraint_pool_begin(indicators->second); - dof_indices += row_length; - ++ind_local; - if (ind_local % VectorizedArray::n_array_elements == - n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements - - n_irreg_components_filled; - } - for(; ind_local::n_array_elements == n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements-n_irreg_components_filled; - } - } - else - // case with vector-valued finite elements - // where all components are included in one - // single vector. Assumption: first come all - // entries to the first component, then all - // entries to the second one, and so on. This - // is ensured by the way MatrixFree reads - // out the indices. - { - AssertDimension (dst[0]->size(), - dof_info.vector_partitioner->size()); - Assert (n_fe_components == n_components, ExcNotImplemented()); - const unsigned int n_local_dofs = - dofs_per_cell * VectorizedArray::n_array_elements * n_components; - const Number * local_dst_number = - reinterpret_cast(values_dofs[0]); - if (at_irregular_cell == false) - { - // check whether there is no constraint at all - if (indicators != indicators_end) - { - // run from one constraint to the next - for ( ; indicators != indicators_end; ++indicators) - { - // distribute values up to the constraint - // (values not constrained) - for (unsigned int j=0; jfirst; ++j) - internal::vector_access (*dst[0], dof_indices[j]) - = local_dst_number[ind_local+j]; - dof_indices += indicators->first; - ind_local += indicators->first; - - // jump over constraints - const unsigned int row_length = - matrix_info.constraint_pool_end(indicators->second) - - matrix_info.constraint_pool_begin(indicators->second); - dof_indices += row_length; - ++ind_local; - } - // distribute values after the last constraint - // (values not constrained) - for(; ind_local 0, ExcInternalError()); - for ( ; indicators != indicators_end; ++indicators) - { - for(unsigned int j=0; jfirst; ++j) - { - // non-constrained case - internal::vector_access (*dst[0], dof_indices[j]) - = local_dst_number[ind_local]; - ++ind_local; - if (ind_local % VectorizedArray::n_array_elements == n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements-n_irreg_components_filled; - } - dof_indices += indicators->first; - - // jump over constraint - const unsigned int row_length = - matrix_info.constraint_pool_end(indicators->second)- - matrix_info.constraint_pool_begin(indicators->second); - dof_indices += row_length; - ++ind_local; - if (ind_local % VectorizedArray::n_array_elements == n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements-n_irreg_components_filled; - } - for(; ind_local::n_array_elements == n_irreg_components_filled) - ind_local += VectorizedArray::n_array_elements-n_irreg_components_filled; - } - Assert (dof_indices == dof_info.end_indices (cell), - ExcInternalError()); - } +template +inline +VectorizedArray * +FEEvaluationBase:: +begin_hessians () +{ + return &hessians_quad[0][0][0]; } -// ------------------------------ access to data fields --------------------- - template + int n_components_, typename Number> inline -Tensor<1,n_components,VectorizedArray > -FEEvaluationBase:: -get_dof_value (const unsigned int dof) const +Tensor<1,n_components_,VectorizedArray > +FEEvaluationBase +::get_dof_value (const unsigned int dof) const { AssertIndexRange (dof, dofs_per_cell); - Tensor<1,n_components,VectorizedArray > return_value (false); + Tensor<1,n_components_,VectorizedArray > return_value (false); for(unsigned int comp=0;compvalues_dofs[comp][dof]; return return_value; @@ -3034,16 +3007,16 @@ get_dof_value (const unsigned int dof) const template + int n_components_, typename Number> inline -Tensor<1,n_components,VectorizedArray > -FEEvaluationBase:: -get_value (const unsigned int q_point) const +Tensor<1,n_components_,VectorizedArray > +FEEvaluationBase +::get_value (const unsigned int q_point) const { Assert (this->values_quad_initialized==true, internal::ExcAccessToUninitializedField()); AssertIndexRange (q_point, n_q_points); - Tensor<1,n_components,VectorizedArray > return_value (false); + Tensor<1,n_components_,VectorizedArray > return_value (false); for(unsigned int comp=0;compvalues_quad[comp][q_point]; return return_value; @@ -3052,17 +3025,17 @@ get_value (const unsigned int q_point) const template + int n_components_, typename Number> inline -Tensor<1,n_components,Tensor<1,dim,VectorizedArray > > -FEEvaluationBase:: -get_gradient (const unsigned int q_point) const +Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > +FEEvaluationBase +::get_gradient (const unsigned int q_point) const { Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); AssertIndexRange (q_point, n_q_points); - Tensor<1,n_components,Tensor<1,dim,VectorizedArray > > grad_out (false); + Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > grad_out (false); // Cartesian cell if (this->cell_type == internal::MatrixFreeFunctions::cartesian) @@ -3072,33 +3045,20 @@ get_gradient (const unsigned int q_point) const grad_out[comp][d] = (this->gradients_quad[comp][d][q_point] * cartesian_data[0][d]); } - // cell with general Jacobian - else if (this->cell_type == internal::MatrixFreeFunctions::general) - { - for(unsigned int comp=0;compgradients_quad[comp][0][q_point]); - for (unsigned e=1; egradients_quad[comp][e][q_point]); - } - } - } - // cell with general Jacobian, but constant - // within the cell - else // if (this->cell_type == internal::MatrixFreeFunctions::affine) + // cell with general/affine Jacobian + else { + const Tensor<2,dim,VectorizedArray > & jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + jacobian[q_point] : jacobian[0]; for(unsigned int comp=0;compgradients_quad[comp][0][q_point]); for (unsigned e=1; egradients_quad[comp][e][q_point]); } } @@ -3108,12 +3068,57 @@ get_gradient (const unsigned int q_point) const +namespace internal +{ + // compute tmp = hess_unit(u) * J^T. do this + // manually because we do not store the lower + // diagonal because of symmetry + template + inline + void + hessian_unit_times_jac (const Tensor<2,dim,VectorizedArray > &jac, + const VectorizedArray hessians_quad[][n_q_points], + const unsigned int q_point, + VectorizedArray tmp[dim][dim]) + { + for (unsigned int d=0; d + int n_components_, typename Number> inline -Tensor<1,n_components,Tensor<2,dim,VectorizedArray > > -FEEvaluationBase:: -get_hessian (const unsigned int q_point) const +Tensor<1,n_components_,Tensor<2,dim,VectorizedArray > > +FEEvaluationBase +::get_hessian (const unsigned int q_point) const { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); @@ -3167,37 +3172,9 @@ get_hessian (const unsigned int q_point) const // because it needs to access unscaled // gradient data VectorizedArray tmp[dim][dim]; + internal::hessian_unit_times_jac (jac, this->hessians_quad[comp], + q_point, tmp); - // compute tmp = hess_unit(u) * J^T. do this - // manually because we do not store the lower - // diagonal because of symmetry - for (unsigned int d=0; dhessians_quad[comp][0][q_point]; - break; - case 2: - tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] + - jac[d][1] * this->hessians_quad[comp][2][q_point]); - tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] + - jac[d][1] * this->hessians_quad[comp][1][q_point]); - break; - case 3: - tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] + - jac[d][1] * this->hessians_quad[comp][3][q_point] + - jac[d][2] * this->hessians_quad[comp][4][q_point]); - tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] + - jac[d][1] * this->hessians_quad[comp][1][q_point] + - jac[d][2] * this->hessians_quad[comp][5][q_point]); - tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] + - jac[d][1] * this->hessians_quad[comp][5][q_point] + - jac[d][2] * this->hessians_quad[comp][2][q_point]); - break; - default: Assert (false, ExcNotImplemented()); - } - } // compute first part of hessian, // J * tmp = J * hess_unit(u) * J^T for (unsigned int d=0; d tmp[dim][dim]; + internal::hessian_unit_times_jac (jac, this->hessians_quad[comp], + q_point, tmp); - // compute tmp = hess_unit(u) * J^T. do this - // manually because we do not store the lower - // diagonal because of symmetry - for (unsigned int d=0; dhessians_quad[comp][0][q_point]; - break; - case 2: - tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] + - jac[d][1] * this->hessians_quad[comp][2][q_point]); - tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] + - jac[d][1] * this->hessians_quad[comp][1][q_point]); - break; - case 3: - tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] + - jac[d][1] * this->hessians_quad[comp][3][q_point] + - jac[d][2] * this->hessians_quad[comp][4][q_point]); - tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] + - jac[d][1] * this->hessians_quad[comp][1][q_point] + - jac[d][2] * this->hessians_quad[comp][5][q_point]); - tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] + - jac[d][1] * this->hessians_quad[comp][5][q_point] + - jac[d][2] * this->hessians_quad[comp][2][q_point]); - break; - default: Assert (false, ExcNotImplemented()); - } - } // compute first part of hessian, // J * tmp = J * hess_unit(u) * J^T for (unsigned int d=0; d > >(hessian_out); + return Tensor<1,n_components_,Tensor<2,dim,VectorizedArray > >(hessian_out); } template + int n_components_, typename Number> inline -Tensor<1,n_components,Tensor<1,dim,VectorizedArray > > -FEEvaluationBase:: -get_hessian_diagonal (const unsigned int q_point) const +Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > +FEEvaluationBase +::get_hessian_diagonal (const unsigned int q_point) const { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); AssertIndexRange (q_point, n_q_points); - Tensor<1,n_components,Tensor<1,dim,VectorizedArray > > hessian_out (false); + Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > hessian_out (false); // Cartesian cell if (this->cell_type == internal::MatrixFreeFunctions::cartesian) @@ -3329,37 +3278,9 @@ get_hessian_diagonal (const unsigned int q_point) const // because it needs to access unscaled // gradient data VectorizedArray tmp[dim][dim]; + internal::hessian_unit_times_jac (jac, this->hessians_quad[comp], + q_point, tmp); - // compute tmp = hess_unit(u) * J^T. do this - // manually because we do not store the lower - // diagonal because of symmetry - for (unsigned int d=0; dhessians_quad[comp][0][q_point]; - break; - case 2: - tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] + - jac[d][1] * this->hessians_quad[comp][2][q_point]); - tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] + - jac[d][1] * this->hessians_quad[comp][1][q_point]); - break; - case 3: - tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] + - jac[d][1] * this->hessians_quad[comp][3][q_point] + - jac[d][2] * this->hessians_quad[comp][4][q_point]); - tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] + - jac[d][1] * this->hessians_quad[comp][1][q_point] + - jac[d][2] * this->hessians_quad[comp][5][q_point]); - tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] + - jac[d][1] * this->hessians_quad[comp][5][q_point] + - jac[d][2] * this->hessians_quad[comp][2][q_point]); - break; - default: Assert (false, ExcNotImplemented()); - } - } // compute only the trace part of hessian, // J * tmp = J * hess_unit(u) * J^T for (unsigned int d=0; d tmp[dim][dim]; + internal::hessian_unit_times_jac (jac, this->hessians_quad[comp], + q_point, tmp); - // compute tmp = hess_unit(u) * J^T. do this - // manually because we do not store the lower - // diagonal because of symmetry - for (unsigned int d=0; dhessians_quad[comp][0][q_point]; - break; - case 2: - tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] + - jac[d][1] * this->hessians_quad[comp][2][q_point]); - tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][2][q_point] + - jac[d][1] * this->hessians_quad[comp][1][q_point]); - break; - case 3: - tmp[0][d] = (jac[d][0] * this->hessians_quad[comp][0][q_point] + - jac[d][1] * this->hessians_quad[comp][3][q_point] + - jac[d][2] * this->hessians_quad[comp][4][q_point]); - tmp[1][d] = (jac[d][0] * this->hessians_quad[comp][3][q_point] + - jac[d][1] * this->hessians_quad[comp][1][q_point] + - jac[d][2] * this->hessians_quad[comp][5][q_point]); - tmp[2][d] = (jac[d][0] * this->hessians_quad[comp][4][q_point] + - jac[d][1] * this->hessians_quad[comp][5][q_point] + - jac[d][2] * this->hessians_quad[comp][2][q_point]); - break; - default: Assert (false, ExcNotImplemented()); - } - } // compute only the trace part of hessian, // J * tmp = J * hess_unit(u) * J^T for (unsigned int d=0; d + int n_components_, typename Number> inline -Tensor<1,n_components,VectorizedArray > -FEEvaluationBase:: -get_laplacian (const unsigned int q_point) const +Tensor<1,n_components_,VectorizedArray > +FEEvaluationBase +::get_laplacian (const unsigned int q_point) const { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); AssertIndexRange (q_point, n_q_points); - Tensor<1,n_components,VectorizedArray > laplacian_out (false); - const Tensor<1,n_components,Tensor<1,dim,VectorizedArray > > hess_diag + Tensor<1,n_components_,VectorizedArray > laplacian_out (false); + const Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > hess_diag = get_hessian_diagonal(q_point); for (unsigned int comp=0; comp + int n_components_, typename Number> inline void -FEEvaluationBase:: -submit_dof_value (const Tensor<1,n_components,VectorizedArray > val_in, - const unsigned int dof) +FEEvaluationBase +::submit_dof_value (const Tensor<1,n_components_,VectorizedArray > val_in, + const unsigned int dof) { #ifdef DEBUG this->dof_values_initialized = true; @@ -3475,12 +3368,12 @@ submit_dof_value (const Tensor<1,n_components,VectorizedArray > val_in, template + int n_components_, typename Number> inline void -FEEvaluationBase:: -submit_value (const Tensor<1,n_components,VectorizedArray > val_in, - const unsigned int q_point) +FEEvaluationBase +::submit_value (const Tensor<1,n_components_,VectorizedArray > val_in, + const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); @@ -3504,13 +3397,13 @@ submit_value (const Tensor<1,n_components,VectorizedArray > val_in, template + int n_components_, typename Number> inline void -FEEvaluationBase:: -submit_gradient (const Tensor<1,n_components, - Tensor<1,dim,VectorizedArray > > grad_in, - const unsigned int q_point) +FEEvaluationBase +::submit_gradient (const Tensor<1,n_components_, + Tensor<1,dim,VectorizedArray > >grad_in, + const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); @@ -3525,26 +3418,20 @@ submit_gradient (const Tensor<1,n_components, this->gradients_quad[comp][d][q_point] = (grad_in[comp][d] * cartesian_data[0][d] * JxW); } - else if (this->cell_type == internal::MatrixFreeFunctions::general) - { - for (unsigned int comp=0; comp new_val = jacobian[q_point][0][d] * grad_in[comp][0]; - for (unsigned e=1; egradients_quad[comp][d][q_point] = new_val * J_value[q_point]; - } - } - else //if (this->cell_type == internal::MatrixFreeFunctions::affine) + else { - const VectorizedArray JxW = J_value[0] * quadrature_weights[q_point]; + const Tensor<2,dim,VectorizedArray > &jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + jacobian[q_point] : jacobian[0]; + const VectorizedArray JxW = + this->cell_type == internal::MatrixFreeFunctions::general ? + J_value[q_point] : J_value[0] * quadrature_weights[q_point]; for (unsigned int comp=0; comp new_val = jacobian[0][0][d] * grad_in[comp][0]; + VectorizedArray new_val = jac[0][d] * grad_in[comp][0]; for (unsigned e=1; egradients_quad[comp][d][q_point] = new_val * JxW; } } @@ -3553,21 +3440,21 @@ submit_gradient (const Tensor<1,n_components, template + int n_components_, typename Number> inline -Tensor<1,n_components,VectorizedArray > -FEEvaluationBase:: -integrate_value () const +Tensor<1,n_components_,VectorizedArray > +FEEvaluationBase +::integrate_value () const { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); Assert (this->values_quad_submitted == true, internal::ExcAccessToUninitializedField()); #endif - Tensor<1,n_components,VectorizedArray > return_value (false); + Tensor<1,n_components_,VectorizedArray > return_value (false); for (unsigned int comp=0; compvalues_quad[comp][0]; - for (unsigned int q=0; qvalues_quad[comp][q]; return (return_value); @@ -3575,33 +3462,33 @@ integrate_value () const -/*----------------------- FEEvaluationAccess -------------------------------*/ +/*----------------------- FEEvaluationAccess --------------------------------*/ template + int n_components_, typename Number> inline -FEEvaluationAccess:: -FEEvaluationAccess (const MatrixFree &data_in, - const unsigned int fe_no, - const unsigned int quad_no_in) +FEEvaluationAccess +::FEEvaluationAccess (const MatrixFree &data_in, + const unsigned int fe_no, + const unsigned int quad_no_in) : - FEEvaluationBase + FEEvaluationBase (data_in, fe_no, quad_no_in) {} -/*-------------------- FEEvaluationAccess scalar --------------------------*/ +/*-------------------- FEEvaluationAccess scalar ----------------------------*/ template inline -FEEvaluationAccess:: -FEEvaluationAccess (const MatrixFree &data_in, - const unsigned int fe_no, - const unsigned int quad_no_in) +FEEvaluationAccess +::FEEvaluationAccess (const MatrixFree &data_in, + const unsigned int fe_no, + const unsigned int quad_no_in) : FEEvaluationBase (data_in, fe_no, quad_no_in) @@ -3612,8 +3499,8 @@ FEEvaluationAccess (const MatrixFree &data_in, template inline VectorizedArray -FEEvaluationAccess:: -get_dof_value (const unsigned int dof) const +FEEvaluationAccess +::get_dof_value (const unsigned int dof) const { AssertIndexRange (dof, dofs_per_cell); return this->values_dofs[0][dof]; @@ -3624,8 +3511,8 @@ get_dof_value (const unsigned int dof) const template inline VectorizedArray -FEEvaluationAccess:: -get_value (const unsigned int q_point) const +FEEvaluationAccess +::get_value (const unsigned int q_point) const { Assert (this->values_quad_initialized==true, internal::ExcAccessToUninitializedField()); @@ -3638,8 +3525,8 @@ get_value (const unsigned int q_point) const template inline Tensor<1,dim,VectorizedArray > -FEEvaluationAccess:: -get_gradient (const unsigned int q_point) const +FEEvaluationAccess +::get_gradient (const unsigned int q_point) const { // could use the base class gradient, but that // involves too many inefficient @@ -3658,29 +3545,17 @@ get_gradient (const unsigned int q_point) const grad_out[d] = (this->gradients_quad[0][d][q_point] * this->cartesian_data[0][d]); } - // cell with general Jacobian - else if (this->cell_type == internal::MatrixFreeFunctions::general) - { - for (unsigned int d=0; djacobian[q_point][d][0] * - this->gradients_quad[0][0][q_point]); - for (unsigned e=1; ejacobian[q_point][d][e] * - this->gradients_quad[0][e][q_point]); - } - } - // cell with general Jacobian, but constant - // within the cell - else // if (this->cell_type == internal::MatrixFreeFunctions::affine) + // cell with general/constant Jacobian + else { + const Tensor<2,dim,VectorizedArray > &jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->jacobian[q_point] : this->jacobian[0]; for (unsigned int d=0; djacobian[0][d][0] * - this->gradients_quad[0][0][q_point]); + grad_out[d] = (jac[d][0] * this->gradients_quad[0][0][q_point]); for (unsigned e=1; ejacobian[0][d][e] * - this->gradients_quad[0][e][q_point]); + grad_out[d] += (jac[d][e] * this->gradients_quad[0][e][q_point]); } } return grad_out; @@ -3691,8 +3566,8 @@ get_gradient (const unsigned int q_point) const template inline Tensor<2,dim,VectorizedArray > -FEEvaluationAccess:: -get_hessian (const unsigned int q_point) const +FEEvaluationAccess +::get_hessian (const unsigned int q_point) const { return BaseClass::get_hessian(q_point)[0]; } @@ -3702,8 +3577,8 @@ get_hessian (const unsigned int q_point) const template inline Tensor<1,dim,VectorizedArray > -FEEvaluationAccess:: -get_hessian_diagonal (const unsigned int q_point) const +FEEvaluationAccess +::get_hessian_diagonal (const unsigned int q_point) const { return BaseClass::get_hessian_diagonal(q_point)[0]; } @@ -3713,8 +3588,8 @@ get_hessian_diagonal (const unsigned int q_point) const template inline VectorizedArray -FEEvaluationAccess:: -get_laplacian (const unsigned int q_point) const +FEEvaluationAccess +::get_laplacian (const unsigned int q_point) const { return BaseClass::get_laplacian(q_point)[0]; } @@ -3724,8 +3599,8 @@ get_laplacian (const unsigned int q_point) const template inline void -FEEvaluationAccess:: -submit_dof_value (const VectorizedArray val_in, +FEEvaluationAccess +::submit_dof_value (const VectorizedArray val_in, const unsigned int dof) { #ifdef DEBUG @@ -3740,9 +3615,9 @@ submit_dof_value (const VectorizedArray val_in, template inline void -FEEvaluationAccess:: -submit_value (const VectorizedArray val_in, - const unsigned int q_point) +FEEvaluationAccess +::submit_value (const VectorizedArray val_in, + const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); @@ -3766,9 +3641,9 @@ submit_value (const VectorizedArray val_in, template inline void -FEEvaluationAccess:: -submit_gradient (const Tensor<1,dim,VectorizedArray > grad_in, - const unsigned int q_point) +FEEvaluationAccess +::submit_gradient (const Tensor<1,dim,VectorizedArray > grad_in, + const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); @@ -3783,24 +3658,20 @@ submit_gradient (const Tensor<1,dim,VectorizedArray > grad_in, this->cartesian_data[0][d] * JxW); } - else if (this->cell_type == internal::MatrixFreeFunctions::general) - { - for (unsigned int d=0; d new_val = this->jacobian[q_point][0][d] * grad_in[0]; - for (unsigned e=1; ejacobian[q_point][e][d] * grad_in[e]; - this->gradients_quad[0][d][q_point] = new_val * this->J_value[q_point]; - } - } - else //if (this->cell_type == internal::MatrixFreeFunctions::affine) + // general/affine cell type + else { - const VectorizedArray JxW = this->J_value[0] * this->quadrature_weights[q_point]; + const Tensor<2,dim,VectorizedArray > & jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->jacobian[q_point] : this->jacobian[0]; + const VectorizedArray JxW = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point]; for (unsigned int d=0; d new_val = this->jacobian[0][0][d] * grad_in[0]; + VectorizedArray new_val = jac[0][d] * grad_in[0]; for (unsigned e=1; ejacobian[0][e][d] * grad_in[e]; + new_val += jac[e][d] * grad_in[e]; this->gradients_quad[0][d][q_point] = new_val * JxW; } } @@ -3811,8 +3682,8 @@ submit_gradient (const Tensor<1,dim,VectorizedArray > grad_in, template inline VectorizedArray -FEEvaluationAccess:: -integrate_value () const +FEEvaluationAccess +::integrate_value () const { return BaseClass::integrate_value()[0]; } @@ -3820,15 +3691,15 @@ integrate_value () const -/*----------------- FEEvaluationAccess vector-valued ----------------------*/ +/*----------------- FEEvaluationAccess vector-valued ------------------------*/ template inline -FEEvaluationAccess:: -FEEvaluationAccess (const MatrixFree &data_in, - const unsigned int fe_no, - const unsigned int quad_no_in) +FEEvaluationAccess +::FEEvaluationAccess (const MatrixFree &data_in, + const unsigned int fe_no, + const unsigned int quad_no_in) : FEEvaluationBase (data_in, fe_no, quad_no_in) @@ -3839,8 +3710,8 @@ FEEvaluationAccess (const MatrixFree &data_in, template inline Tensor<2,dim,VectorizedArray > -FEEvaluationAccess:: -get_gradient (const unsigned int q_point) const +FEEvaluationAccess +::get_gradient (const unsigned int q_point) const { return BaseClass::get_gradient (q_point); } @@ -3850,8 +3721,8 @@ get_gradient (const unsigned int q_point) const template inline VectorizedArray -FEEvaluationAccess:: -get_divergence (const unsigned int q_point) const +FEEvaluationAccess +::get_divergence (const unsigned int q_point) const { Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); @@ -3868,32 +3739,18 @@ get_divergence (const unsigned int q_point) const divergence += (this->gradients_quad[d][d][q_point] * this->cartesian_data[0][d]); } - // cell with general Jacobian - else if (this->cell_type == internal::MatrixFreeFunctions::general) - { - divergence = (this->jacobian[q_point][0][0] * - this->gradients_quad[0][0][q_point]); - for (unsigned e=1; ejacobian[q_point][0][e] * - this->gradients_quad[0][e][q_point]); - for (unsigned int d=1; djacobian[q_point][d][e] * - this->gradients_quad[d][e][q_point]); - } - // cell with general Jacobian, but constant - // within the cell - else // if (this->cell_type == internal::MatrixFreeFunctions::affine) + // cell with general/constant Jacobian + else { - divergence = (this->jacobian[0][0][0] * - this->gradients_quad[0][0][q_point]); + const Tensor<2,dim,VectorizedArray > &jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->jacobian[q_point] : this->jacobian[0]; + divergence = (jac[0][0] * this->gradients_quad[0][0][q_point]); for (unsigned e=1; ejacobian[0][0][e] * - this->gradients_quad[0][e][q_point]); + divergence += (jac[0][e] * this->gradients_quad[0][e][q_point]); for (unsigned int d=1; djacobian[0][d][e] * - this->gradients_quad[d][e][q_point]); + divergence += (jac[d][e] * this->gradients_quad[d][e][q_point]); } return divergence; } @@ -3903,8 +3760,8 @@ get_divergence (const unsigned int q_point) const template inline SymmetricTensor<2,dim,VectorizedArray > -FEEvaluationAccess:: -get_symmetric_gradient (const unsigned int q_point) const +FEEvaluationAccess +::get_symmetric_gradient (const unsigned int q_point) const { // copy from generic function into // dim-specialization function @@ -3940,8 +3797,8 @@ get_symmetric_gradient (const unsigned int q_point) const template inline Tensor<1,dim==2?1:dim,VectorizedArray > -FEEvaluationAccess:: -get_curl (const unsigned int q_point) const +FEEvaluationAccess +::get_curl (const unsigned int q_point) const { // copy from generic function into // dim-specialization function @@ -3972,8 +3829,8 @@ get_curl (const unsigned int q_point) const template inline Tensor<2,dim,VectorizedArray > -FEEvaluationAccess:: -get_hessian_diagonal (const unsigned int q_point) const +FEEvaluationAccess +::get_hessian_diagonal (const unsigned int q_point) const { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); @@ -3987,8 +3844,8 @@ get_hessian_diagonal (const unsigned int q_point) const template inline Tensor<3,dim,VectorizedArray > -FEEvaluationAccess:: -get_hessian (const unsigned int q_point) const +FEEvaluationAccess +::get_hessian (const unsigned int q_point) const { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); @@ -4001,9 +3858,9 @@ get_hessian (const unsigned int q_point) const template inline void -FEEvaluationAccess:: -submit_gradient (const Tensor<2,dim,VectorizedArray > grad_in, - const unsigned int q_point) +FEEvaluationAccess +::submit_gradient (const Tensor<2,dim,VectorizedArray > grad_in, + const unsigned int q_point) { BaseClass::submit_gradient (grad_in, q_point); } @@ -4013,9 +3870,10 @@ submit_gradient (const Tensor<2,dim,VectorizedArray > grad_in, template inline void -FEEvaluationAccess:: -submit_gradient (const Tensor<1,dim,Tensor<1,dim,VectorizedArray > > grad_in, - const unsigned int q_point) +FEEvaluationAccess +::submit_gradient (const Tensor<1,dim,Tensor<1,dim,VectorizedArray > > + grad_in, + const unsigned int q_point) { BaseClass::submit_gradient(grad_in, q_point); } @@ -4025,10 +3883,10 @@ submit_gradient (const Tensor<1,dim,Tensor<1,dim,VectorizedArray > > gra template inline void -FEEvaluationAccess:: -submit_symmetric_gradient (const SymmetricTensor<2,dim,VectorizedArray > - sym_grad, - const unsigned int q_point) +FEEvaluationAccess +::submit_symmetric_gradient(const SymmetricTensor<2,dim,VectorizedArray > + sym_grad, + const unsigned int q_point) { // could have used base class operator, but // that involves some overhead which is @@ -4057,52 +3915,31 @@ submit_symmetric_gradient (const SymmetricTensor<2,dim,VectorizedArray > this->cartesian_data[0][e]); } } - else if (this->cell_type == internal::MatrixFreeFunctions::general) + // general/affine cell type + else { + const VectorizedArray JxW = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point]; + const Tensor<2,dim,VectorizedArray > &jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->jacobian[q_point] : this->jacobian[0]; VectorizedArray weighted [dim][dim]; - { - const VectorizedArray JxW = this->J_value[q_point]; - for (unsigned int i=0; i value = sym_grad.access_raw_entry(counter) * JxW; - weighted[i][j] = value; - weighted[j][i] = value; - } - } - for (unsigned int comp=0; comp new_val = this->jacobian[q_point][0][d] * weighted[comp][0]; - for (unsigned e=1; ejacobian[q_point][e][d] * weighted[comp][e]; - this->gradients_quad[comp][d][q_point] = new_val; + const VectorizedArray value = sym_grad.access_raw_entry(counter) * JxW; + weighted[i][j] = value; + weighted[j][i] = value; } - } - else //if (this->cell_type == internal::MatrixFreeFunctions::affine) - { - VectorizedArray weighted [dim][dim]; - { - const VectorizedArray JxW = (this->J_value[0] * - this->quadrature_weights[q_point]); - for (unsigned int i=0; i value = sym_grad.access_raw_entry(counter) * JxW; - weighted[i][j] = value; - weighted[j][i] = value; - } - } for (unsigned int comp=0; comp new_val = this->jacobian[q_point][0][d] * weighted[comp][0]; + VectorizedArray new_val = jac[0][d] * weighted[comp][0]; for (unsigned e=1; ejacobian[q_point][e][d] * weighted[comp][e]; + new_val += jac[e][d] * weighted[comp][e]; this->gradients_quad[comp][d][q_point] = new_val; } } @@ -4113,9 +3950,9 @@ submit_symmetric_gradient (const SymmetricTensor<2,dim,VectorizedArray > template inline void -FEEvaluationAccess:: -submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray > curl, - const unsigned int q_point) +FEEvaluationAccess +::submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray > curl, + const unsigned int q_point) { Tensor<2,dim,VectorizedArray > grad; switch (dim) @@ -4146,13 +3983,13 @@ submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray > curl, /*----------------------- FEEvaluationGeneral -------------------------------*/ -template inline -FEEvaluationGeneral:: -FEEvaluationGeneral (const MatrixFree &data_in, - const unsigned int fe_no, - const unsigned int quad_no_in) +FEEvaluationGeneral +::FEEvaluationGeneral (const MatrixFree &data_in, + const unsigned int fe_no, + const unsigned int quad_no_in) : BaseClass (data_in, fe_no, quad_no_in) { @@ -4245,259 +4082,440 @@ FEEvaluationGeneral (const MatrixFree &data_in, -template -inline -void -FEEvaluationGeneral:: -evaluate (bool evaluate_val, bool evaluate_grad, bool evaluate_lapl) +namespace internal { - Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - Assert (this->dof_values_initialized == true, - internal::ExcAccessToUninitializedField()); + // evaluates the given shape data in 1d-3d + // using the tensor product form. does not use + // the tensor product form and corresponds to + // a usual matrix-matrix product + template + inline + void + apply_tensor_product (const VectorizedArray*shape_data, + const VectorizedArray in [], + VectorizedArray out []) + { + AssertIndexRange (direction, dim); + const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, + nn = dof_to_quad ? n_q_points_1d : (fe_degree+1); - const VectorizedArray * val = this->data.shape_values.begin(); - const VectorizedArray * grad = this->data.shape_gradients.begin(); - const VectorizedArray * hess = this->data.shape_hessians.begin(); + const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); + const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); + const int stride = Utilities::fixed_int_power::value; - for(unsigned int comp=0;comp val0; + if (dof_to_quad == true) + val0 = shape_data[col]; + else + val0 = shape_data[col*n_q_points_1d]; + VectorizedArray res0 = val0 * in[0]; + for (int ind=1; ind + inline + void + do_evaluate (FEEval &fe_eval, + const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl) + { + Assert (fe_eval.cell != numbers::invalid_unsigned_int, + ExcNotInitialized()); + Assert (fe_eval.dof_values_initialized == true, + internal::ExcAccessToUninitializedField()); + + const unsigned int temp_size = fe_eval.dofs_per_cell > fe_eval.n_q_points ? + fe_eval.dofs_per_cell : fe_eval.n_q_points; + const unsigned int n_components = fe_eval.n_components; + const unsigned int dim = fe_eval.dimension; + + for(unsigned int c=0; c temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points]; - VectorizedArray temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points]; + VectorizedArray temp1[temp_size]; + VectorizedArray temp2[temp_size]; - if (dim == 3) - { - if (evaluate_grad == true) - { - // grad x - apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1); - apply_tensor_prod<1,true,false> (val, temp1, temp2); - apply_tensor_prod<2,true,false> (val, temp2, this->gradients_quad[comp][0]); - } + switch (dim) + { + case 3: - if (evaluate_lapl == true) - { - // grad xz - if (evaluate_grad == false) - { - apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1); - apply_tensor_prod<1,true,false> (val, temp1, temp2); - } - apply_tensor_prod<2,true,false>(grad, temp2, this->hessians_quad[comp][4]); + if (evaluate_grad == true) + { + // grad x + fe_eval.template apply_gradients<0,true,false> + (fe_eval.values_dofs[c], temp1); + fe_eval.template apply_values<1,true,false> + (temp1, temp2); + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.gradients_quad[c][0]); + } - // grad xy - apply_tensor_prod<1,true,false>(grad, temp1, temp2); - apply_tensor_prod<2,true,false> (val, temp2, this->hessians_quad[comp][3]); + if (evaluate_lapl == true) + { + // grad xz + if (evaluate_grad == false) + { + fe_eval.template apply_gradients<0,true,false> + (fe_eval.values_dofs[c], temp1); + fe_eval.template apply_values<1,true,false> + (temp1, temp2); + } + fe_eval.template apply_gradients<2,true,false> + (temp2, fe_eval.hessians_quad[c][4]); + + // grad xy + fe_eval.template apply_gradients<1,true,false> + (temp1, temp2); + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.hessians_quad[c][3]); + + // grad xx + fe_eval.template apply_hessians<0,true,false> + (fe_eval.values_dofs[c], temp1); + fe_eval.template apply_values<1,true,false> + (temp1, temp2); + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.hessians_quad[c][0]); + } - // grad xx - apply_tensor_prod<0,true,false> (hess, this->values_dofs[comp], temp1); - apply_tensor_prod<1,true,false> (val, temp1, temp2); - apply_tensor_prod<2,true,false> (val, temp2, this->hessians_quad[comp][0]); - } + // grad y + fe_eval.template apply_values<0,true,false> + (fe_eval.values_dofs[c], temp1); + if (evaluate_grad == true) + { + fe_eval.template apply_gradients<1,true,false> + (temp1, temp2); + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.gradients_quad[c][1]); + } - // grad y - apply_tensor_prod<0,true,false> (val, this->values_dofs[comp], temp1); - if (evaluate_grad == true) - { - apply_tensor_prod<1,true,false> (grad, temp1, temp2); - apply_tensor_prod<2,true,false> (val, temp2, this->gradients_quad[comp][1]); - } + if (evaluate_lapl == true) + { + // grad yz + if (evaluate_grad == false) + fe_eval.template apply_gradients<1,true,false> + (temp1, temp2); + fe_eval.template apply_gradients<2,true,false> + (temp2, fe_eval.hessians_quad[c][5]); + + // grad yy + fe_eval.template apply_hessians<1,true,false> + (temp1, temp2); + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.hessians_quad[c][1]); + } - if (evaluate_lapl == true) - { - // grad yz - if (evaluate_grad == false) - apply_tensor_prod<1,true,false> (grad, temp1, temp2); - apply_tensor_prod<2,true,false> (grad, temp2, this->hessians_quad[comp][5]); + // grad z: can use the values applied in x direction stored in temp1 + fe_eval.template apply_values<1,true,false> + (temp1, temp2); + if (evaluate_grad == true) + fe_eval.template apply_gradients<2,true,false> + (temp2, fe_eval.gradients_quad[c][2]); - // grad yy - apply_tensor_prod<1,true,false> (hess, temp1, temp2); - apply_tensor_prod<2,true,false> (val, temp2, this->hessians_quad[comp][1]); - } + // grad zz: can use the values applied in x and y direction stored + // in temp2 + if (evaluate_lapl == true) + fe_eval.template apply_hessians<2,true,false> + (temp2, fe_eval.hessians_quad[c][2]); - // grad z: can use the values applied in x direction stored in temp1 - apply_tensor_prod<1,true,false> (val, temp1, temp2); - if (evaluate_grad == true) - apply_tensor_prod<2,true,false> (grad, temp2, this->gradients_quad[comp][2]); + // val: can use the values applied in x & y direction stored in temp2 + if (evaluate_val == true) + fe_eval.template apply_values<2,true,false> + (temp2, fe_eval.values_quad[c]); - // grad zz: can use the values applied in x and y direction stored in temp2 - if (evaluate_lapl == true) - apply_tensor_prod<2,true,false> (hess, temp2, this->hessians_quad[comp][2]); + break; - // val: can use the values applied in x & y direction stored in temp2 - if (evaluate_val == true) - apply_tensor_prod<2,true,false> (val, temp2, this->values_quad[comp]); - } - else if (dim == 2) - { - // grad x - if (evaluate_grad == true) - { - apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1); - apply_tensor_prod<1,true,false> (val, temp1, this->gradients_quad[comp][0]); - } - if (evaluate_lapl == true) - { - // grad xy - if (evaluate_grad == false) - apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], temp1); - apply_tensor_prod<1,true,false> (grad, temp1, this->hessians_quad[comp][2]); + case 2: + + // grad x + if (evaluate_grad == true) + { + fe_eval.template apply_gradients<0,true,false> + (fe_eval.values_dofs[c], temp1); + fe_eval.template apply_values<1,true,false> + (temp1, fe_eval.gradients_quad[c][0]); + } + if (evaluate_lapl == true) + { + // grad xy + if (evaluate_grad == false) + fe_eval.template apply_gradients<0,true,false> + (fe_eval.values_dofs[c], temp1); + fe_eval.template apply_gradients<1,true,false> + (temp1, fe_eval.hessians_quad[c][2]); + + // grad xx + fe_eval.template apply_hessians<0,true,false> + (fe_eval.values_dofs[c], temp1); + fe_eval.template apply_values<1,true,false> + (temp1, fe_eval.hessians_quad[c][0]); + } + + // grad y + fe_eval.template apply_values<0,true,false> + (fe_eval.values_dofs[c], temp1); + if (evaluate_grad == true) + fe_eval.template apply_gradients<1,true,false> + (temp1, fe_eval.gradients_quad[c][1]); + + // grad yy + if (evaluate_lapl == true) + fe_eval.template apply_hessians<1,true,false> + (temp1, fe_eval.hessians_quad[c][1]); + + // val: can use values applied in x + if (evaluate_val == true) + fe_eval.template apply_values<1,true,false> + (temp1, fe_eval.values_quad[c]); + + break; + + case 1: + if (evaluate_val == true) + fe_eval.template apply_values<0,true,false> + (fe_eval.values_dofs[c], fe_eval.values_quad[c]); + if (evaluate_grad == true) + fe_eval.template apply_gradients<0,true,false> + (fe_eval.values_dofs[c], fe_eval.gradients_quad[c][0]); + if (evaluate_lapl == true) + fe_eval.template apply_hessians<0,true,false> + (fe_eval.values_dofs[c], fe_eval.hessians_quad[c][0]); + break; + + default: + Assert (false, ExcNotImplemented()); + } + } + +#ifdef DEBUG + if (evaluate_val == true) + fe_eval.values_quad_initialized = true; + if (evaluate_grad == true) + fe_eval.gradients_quad_initialized = true; + if (evaluate_lapl == true) + fe_eval.hessians_quad_initialized = true; +#endif + } + + + + template + inline + void + do_integrate (FEEval &fe_eval, + const bool integrate_val, + const bool integrate_grad) + { + Assert (fe_eval.cell != numbers::invalid_unsigned_int, ExcNotInitialized()); + if (integrate_val == true) + Assert (fe_eval.values_quad_submitted == true, + ExcAccessToUninitializedField()); + if (integrate_grad == true) + Assert (fe_eval.gradients_quad_submitted == true, + ExcAccessToUninitializedField()); + + const unsigned int temp_size = fe_eval.dofs_per_cell > fe_eval.n_q_points ? + fe_eval.dofs_per_cell : fe_eval.n_q_points; + const unsigned int n_components = fe_eval.n_components; + const unsigned int dim = fe_eval.dimension; + + + for(unsigned int c=0; c temp1[temp_size]; + VectorizedArray temp2[temp_size]; + + switch (dim) + { + case 3: + + if (integrate_val == true) + { + // val + fe_eval.template apply_values<0,false,false> + (fe_eval.values_quad[c], temp1); + } + if (integrate_grad == true) + { + // grad x: can sum to temporary value in temp1 + if (integrate_val == true) + fe_eval.template apply_gradients<0,false,true> + (fe_eval.gradients_quad[c][0], temp1); + else + fe_eval.template apply_gradients<0,false,false> + (fe_eval.gradients_quad[c][0], temp1); + } + fe_eval.template apply_values<1,false,false> + (temp1, temp2); + if (integrate_grad == true) + { + // grad y: can sum to temporary x value in temp2 + fe_eval.template apply_values<0,false,false> + (fe_eval.gradients_quad[c][1], temp1); + fe_eval.template apply_gradients<1,false,true> + (temp1, temp2); + } + fe_eval.template apply_values<2,false,false> + (temp2, fe_eval.values_dofs[c]); + if (integrate_grad == true) + { + // grad z: can sum to temporary x and y value in output + fe_eval.template apply_values<0,false,false> + (fe_eval.gradients_quad[c][2], temp1); + fe_eval.template apply_values<1,false,false> + (temp1, temp2); + fe_eval.template apply_gradients<2,false,true> + (temp2, fe_eval.values_dofs[c]); + } + + break; + + case 2: + + // val + if (integrate_val == true) + fe_eval.template apply_values<0,false,false> + (fe_eval.values_quad[c], temp1); + if (integrate_grad == true) + { + //grad x + if (integrate_val == true) + fe_eval.template apply_gradients<0,false,true> + (fe_eval.gradients_quad[c][0], temp1); + else + fe_eval.template apply_gradients<0,false,false> + (fe_eval.gradients_quad[c][0], temp1); + } + fe_eval.template apply_values<1,false,false> + (temp1, fe_eval.values_dofs[c]); + if (integrate_grad == true) + { + // grad y + fe_eval.template apply_values<0,false,false> + (fe_eval.gradients_quad[c][1], temp1); + fe_eval.template apply_gradients<1,false,true> + (temp1, fe_eval.values_dofs[c]); + } - // grad xx - apply_tensor_prod<0,true,false> (hess, this->values_dofs[comp], temp1); - apply_tensor_prod<1,true,false> (val, temp1, this->hessians_quad[comp][0]); - } + break; - // grad y - apply_tensor_prod<0,true,false> (val, this->values_dofs[comp], temp1); - if (evaluate_grad == true) - apply_tensor_prod<1,true,false> (grad, temp1, this->gradients_quad[comp][1]); + case 1: - // grad yy - if (evaluate_lapl == true) - apply_tensor_prod<1,true,false> (hess, temp1, this->hessians_quad[comp][1]); + if (integrate_grad == true) + fe_eval.template apply_gradients<0,false,false> + (fe_eval.gradients_quad[c][0], fe_eval.values_dofs[c]); + if (integrate_val == true) + { + if (integrate_grad == true) + fe_eval.template apply_values<0,false,true> + (fe_eval.values_quad[c], fe_eval.values_dofs[c]); + else + fe_eval.template apply_values<0,false,false> + (fe_eval.values_quad[c], fe_eval.values_dofs[c]); + } + break; - // val: can use values applied in x - if (evaluate_val == true) - apply_tensor_prod<1,true,false> (val, temp1, this->values_quad[comp]); - } - else if (dim == 1) - { - if (evaluate_val == true) - apply_tensor_prod<0,true,false> (val, this->values_dofs[comp], - this->values_quad[comp]); - if (evaluate_grad == true) - apply_tensor_prod<0,true,false> (grad, this->values_dofs[comp], - this->gradients_quad[comp][0]); - if (evaluate_lapl == true) - apply_tensor_prod<0,true,false> (hess, this->values_dofs[comp], - this->hessians_quad[comp][0]); - } + default: + Assert (false, ExcNotImplemented()); + } } #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; + fe_eval.dof_values_initialized = true; #endif + } } -template inline void -FEEvaluationGeneral:: -integrate (bool integrate_val,bool integrate_grad) +FEEvaluationGeneral +::evaluate (const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl) { -#ifdef DEBUG - Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - if (integrate_val == true) - Assert (this->values_quad_submitted == true, - internal::ExcAccessToUninitializedField()); - if (integrate_grad == true) - Assert (this->gradients_quad_submitted == true, - internal::ExcAccessToUninitializedField()); -#endif - - const VectorizedArray * val = this->data.shape_values.begin(); - const VectorizedArray * grad = this->data.shape_gradients.begin(); + internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl); +} - for(unsigned int comp=0;comp temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points]; - VectorizedArray temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points]; - if (dim == 3) - { - if (integrate_val == true) - { - // val - apply_tensor_prod<0,false,false> (val, this->values_quad[comp], temp1); - } - if (integrate_grad == true) - { - // grad x: can sum to temporary value in temp1 - if (integrate_val == true) - apply_tensor_prod<0,false,true> - (grad, this->gradients_quad[comp][0],temp1); - else - apply_tensor_prod<0,false,false> - (grad, this->gradients_quad[comp][0],temp1); - } - apply_tensor_prod<1,false,false> (val, temp1, temp2); - if (integrate_grad == true) - { - // grad y: can sum to temporary x value in temp2 - apply_tensor_prod<0,false,false> (val, this->gradients_quad[comp][1], temp1); - apply_tensor_prod<1,false,true> (grad, temp1, temp2); - } - apply_tensor_prod<2,false,false> (val, temp2, this->values_dofs[comp]); - if (integrate_grad == true) - { - // grad z: can sum to temporary x and y value in output - apply_tensor_prod<0,false,false> (val, this->gradients_quad[comp][2], temp1); - apply_tensor_prod<1,false,false> (val, temp1, temp2); - apply_tensor_prod<2,false,true> (grad, temp2, this->values_dofs[comp]); - } - } - else if (dim == 2) - { - // val - if (integrate_val == true) - apply_tensor_prod<0,false,false> (val, this->values_quad[comp], temp1); - if (integrate_grad == true) - { - //grad x - if (integrate_val == true) - apply_tensor_prod<0,false,true> - (grad, this->gradients_quad[comp][0],temp1); - else - apply_tensor_prod<0,false,false> - (grad, this->gradients_quad[comp][0],temp1); - } - apply_tensor_prod<1,false,false> (val, temp1, this->values_dofs[comp]); - if (integrate_grad == true) - { - // grad y - apply_tensor_prod<0,false,false> (grad, this->gradients_quad[comp][1], temp1); - apply_tensor_prod<1,false,true> (val, temp1, this->values_dofs[comp]); - } - } - else if (dim == 1) - { - if (integrate_grad == true) - apply_tensor_prod<0,false,false> (grad, this->gradients_quad[comp][0], - this->values_dofs[comp]); - if (integrate_val == true) - { - if (integrate_grad == true) - apply_tensor_prod<0,false,true> (val, this->values_quad[comp], - this->values_dofs[comp]); - else - apply_tensor_prod<0,false,false> (val, this->values_quad[comp], - this->values_dofs[comp]); - } - } - } -#ifdef DEBUG - this->dof_values_initialized = true; -#endif +template +inline +void +FEEvaluationGeneral +::integrate (const bool integrate_val, + const bool integrate_grad) +{ + internal::do_integrate (*this, integrate_val, integrate_grad); } -template inline Point > -FEEvaluationGeneral:: -quadrature_point (const unsigned int q) const +FEEvaluationGeneral +::quadrature_point (const unsigned int q) const { Assert (this->mapping_info.quadrature_points_initialized == true, ExcNotInitialized()); @@ -4535,95 +4553,64 @@ quadrature_point (const unsigned int q) const } - // General tensor product application for up - // to three spatial dimensions. Does not - // assume any symmetry in the shape values - // field -template template inline void -FEEvaluationGeneral:: -apply_tensor_prod (const VectorizedArray*shape_data, - const VectorizedArray input [], - VectorizedArray output []) +FEEvaluationGeneral +::apply_values(const VectorizedArray in [], + VectorizedArray out []) { - AssertIndexRange (direction, dim); - const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, - nn = dof_to_quad ? n_q_points_1d : (fe_degree+1); + internal::apply_tensor_product + (this->data.shape_values.begin(), in, out); +} - const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); - const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); - const int stride = ((direction > 0 ? nn : 1 ) * - (direction > 1 ? nn : 1)); - const VectorizedArray * in = &input[0]; - VectorizedArray * out = &output[0]; - for (int i2=0; i2 val0; - if (dof_to_quad == true) - val0 = shape_data[col]; - else - val0 = shape_data[col*n_q_points_1d]; - VectorizedArray res0 = val0 * in[0]; - for (int ind=1; ind +template +inline +void +FEEvaluationGeneral +::apply_gradients(const VectorizedArray in [], + VectorizedArray out []) +{ + internal::apply_tensor_product + (this->data.shape_gradients.begin(), in, out); +} + + + +template +template +inline +void +FEEvaluationGeneral +::apply_hessians(const VectorizedArray in [], + VectorizedArray out []) +{ + internal::apply_tensor_product + (this->data.shape_hessians.begin(), in, out); } -/*----------------------- FEEvaluation -------------------------------*/ +/*-------------------------- FEEvaluation -----------------------------------*/ -template inline -FEEvaluation:: -FEEvaluation (const MatrixFree &data_in, - const unsigned int fe_no, - const unsigned int quad_no) +FEEvaluation +::FEEvaluation (const MatrixFree &data_in, + const unsigned int fe_no, + const unsigned int quad_no) : BaseClass (data_in, fe_no, quad_no) { @@ -4685,257 +4672,42 @@ FEEvaluation (const MatrixFree &data_in, -template inline void -FEEvaluation:: -evaluate (bool evaluate_val, bool evaluate_grad, bool evaluate_lapl) +FEEvaluation +::evaluate (const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl) { - Assert (this->cell != numbers::invalid_unsigned_int, - ExcNotInitialized()); - Assert (this->dof_values_initialized == true, - internal::ExcAccessToUninitializedField()); - - for(unsigned int comp=0;comp temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points]; - VectorizedArray temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points]; - - if (dim == 3) - { - if (evaluate_grad == true) - { - // grad x - apply_gradients<0,true,false> (this->values_dofs[comp], temp1); - apply_values<1,true,false> (temp1, temp2); - apply_values<2,true,false> (temp2, this->gradients_quad[comp][0]); - } - - if (evaluate_lapl == true) - { - // grad xz - if (evaluate_grad == false) - { - apply_gradients<0,true,false> (this->values_dofs[comp], temp1); - apply_values<1,true,false> (temp1, temp2); - } - apply_gradients<2,true,false>(temp2, this->hessians_quad[comp][4]); - - // grad xy - apply_gradients<1,true,false>(temp1, temp2); - apply_values<2,true,false> (temp2, this->hessians_quad[comp][3]); - - // grad xx - apply_hessians<0,true,false> (this->values_dofs[comp], temp1); - apply_values<1,true,false> (temp1, temp2); - apply_values<2,true,false> (temp2, this->hessians_quad[comp][0]); - } - - // grad y - apply_values<0,true,false> (this->values_dofs[comp], temp1); - if (evaluate_grad == true) - { - apply_gradients<1,true,false> (temp1, temp2); - apply_values<2,true,false> (temp2, this->gradients_quad[comp][1]); - } - - if (evaluate_lapl == true) - { - // grad yz - if (evaluate_grad == false) - apply_gradients<1,true,false> (temp1, temp2); - apply_gradients<2,true,false> (temp2, this->hessians_quad[comp][5]); - - // grad yy - apply_hessians<1,true,false> (temp1, temp2); - apply_values<2,true,false> (temp2, this->hessians_quad[comp][1]); - } - - // grad z: can use the values applied in x direction stored in temp1 - apply_values<1,true,false> (temp1, temp2); - if (evaluate_grad == true) - apply_gradients<2,true,false> (temp2, this->gradients_quad[comp][2]); - - // grad zz: can use the values applied in x and y direction stored in temp2 - if (evaluate_lapl == true) - apply_hessians<2,true,false> (temp2, this->hessians_quad[comp][2]); - - // val: can use the values applied in x & y direction stored in temp2 - if (evaluate_val == true) - apply_values<2,true,false> (temp2, this->values_quad[comp]); - } - else if (dim == 2) - { - // grad x - if (evaluate_grad == true) - { - apply_gradients<0,true,false> (this->values_dofs[comp], temp1); - apply_values<1,true,false> (temp1, this->gradients_quad[comp][0]); - } - if (evaluate_lapl == true) - { - // grad xy - if (evaluate_grad == false) - apply_gradients<0,true,false> (this->values_dofs[comp], temp1); - apply_gradients<1,true,false> (temp1, this->hessians_quad[comp][2]); - - // grad xx - apply_hessians<0,true,false> (this->values_dofs[comp], temp1); - apply_values<1,true,false> (temp1, this->hessians_quad[comp][0]); - } - - // grad y - apply_values<0,true,false> (this->values_dofs[comp], temp1); - if (evaluate_grad == true) - apply_gradients<1,true,false> (temp1, this->gradients_quad[comp][1]); - - // grad yy - if (evaluate_lapl == true) - apply_hessians<1,true,false> (temp1, this->hessians_quad[comp][1]); - - // val: can use values applied in x - if (evaluate_val == true) - apply_values<1,true,false> (temp1, this->values_quad[comp]); - } - else if (dim == 1) - { - if (evaluate_val == true) - apply_values<0,true,false> (this->values_dofs[comp], - this->values_quad[comp]); - if (evaluate_grad == true) - apply_gradients<0,true,false> (this->values_dofs[comp], - this->gradients_quad[comp][0]); - if (evaluate_lapl == true) - apply_hessians<0,true,false> (this->values_dofs[comp], - this->hessians_quad[comp][0]); - } - } - -#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 + internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl); } -template inline void -FEEvaluation:: -integrate (bool integrate_val,bool integrate_grad) +FEEvaluation +::integrate (bool integrate_val,bool integrate_grad) { -#ifdef DEBUG - Assert (this->cell != numbers::invalid_unsigned_int, - ExcNotInitialized()); - if (integrate_val == true) - Assert (this->values_quad_submitted == true, - internal::ExcAccessToUninitializedField()); - if (integrate_grad == true) - Assert (this->gradients_quad_submitted == true, - internal::ExcAccessToUninitializedField()); -#endif - - for(unsigned int comp=0;comp temp1[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points]; - VectorizedArray temp2[fe_degree >= n_q_points_1d ? dofs_per_cell : n_q_points]; - - if (dim == 3) - { - if (integrate_val == true) - { - // val - apply_values<0,false,false> (this->values_quad[comp], temp1); - } - if (integrate_grad == true) - { - // grad x: can sum to temporary value in temp1 - if (integrate_val == true) - apply_gradients<0,false,true> (this->gradients_quad[comp][0], - temp1); - else - apply_gradients<0,false,false> (this->gradients_quad[comp][0], - temp1); - } - apply_values<1,false,false> (temp1, temp2); - if (integrate_grad == true) - { - // grad y: can sum to temporary x value in temp2 - apply_values<0,false,false> (this->gradients_quad[comp][1], temp1); - apply_gradients<1,false,true> (temp1, temp2); - } - apply_values<2,false,false> (temp2, this->values_dofs[comp]); - if (integrate_grad == true) - { - // grad z: can sum to temporary x and y value in output - apply_values<0,false,false> (this->gradients_quad[comp][2], temp1); - apply_values<1,false,false> (temp1, temp2); - apply_gradients<2,false,true> (temp2, this->values_dofs[comp]); - } - } - else if (dim == 2) - { - // val - if (integrate_val == true) - apply_values<0,false,false> (this->values_quad[comp], temp1); - if (integrate_grad == true) - { - //grad x - if (integrate_val == true) - apply_gradients<0,false,true> (this->gradients_quad[comp][0], - temp1); - else - apply_gradients<0,false,false> (this->gradients_quad[comp][0], - temp1); - } - apply_values<1,false,false> (temp1, this->values_dofs[comp]); - if (integrate_grad == true) - { - // grad y - apply_values<0,false,false> (this->gradients_quad[comp][1], temp1); - apply_gradients<1,false,true> (temp1, this->values_dofs[comp]); - } - } - else if (dim == 1) - { - if (integrate_grad == true) - apply_gradients<0,false,false> (this->gradients_quad[comp][0], - this->values_dofs[comp]); - if (integrate_val == true) - { - if (integrate_grad == true) - apply_values<0,false,true> (this->values_quad[comp], - this->values_dofs[comp]); - else - apply_values<0,false,false> (this->values_quad[comp], - this->values_dofs[comp]); - } - } - } -#ifdef DEBUG - this->dof_values_initialized = true; -#endif + internal::do_integrate (*this, integrate_val, integrate_grad); } -// ----------------- optimized implementation tensor product symmetric case +/*----------------- optimized implementation tensor product symmetric case --*/ -template template inline void -FEEvaluation:: -apply_values (const VectorizedArray input [], - VectorizedArray output []) +FEEvaluation +::apply_values (const VectorizedArray in [], + VectorizedArray out []) { AssertIndexRange (direction, dim); const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, @@ -4945,50 +4717,78 @@ apply_values (const VectorizedArray input [], const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); - const int stride = ((direction > 0 ? nn : 1 ) * - (direction > 1 ? nn : 1)); - - const VectorizedArray * in = &input[0]; - VectorizedArray * out = &output[0]; + const int stride = Utilities::fixed_int_power::value; + + // This loop specializes the general + // application of tensor-product based + // elements for "symmetric" finite elements, + // i.e., when the shape functions are + // symmetric about 0.5 and the quadrature + // points are, too. In that case, the 1D shape + // values read (sorted lexicographically, rows + // run over 1D dofs, columns over quadrature + // points): + // Q2 --> [ 0.687 0 -0.087 ] + // [ 0.4 1 0.4 ] + // [-0.087 0 0.687 ] + // Q3 --> [ 0.66 0.003 0.002 0.049 ] + // [ 0.521 1.005 -0.01 -0.230 ] + // [-0.230 -0.01 1.005 0.521 ] + // [ 0.049 0.002 0.003 0.66 ] + // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ] + // [ 0.608 1.059 0 0.039 0.176 ] + // [-0.409 -0.113 1 -0.113 -0.409 ] + // [ 0.176 0.039 0 1.059 0.608 ] + // [-0.032 -0.007 0 0.022 0.658 ] + // + // In these matrices, we want to use avoid + // computations involving zeros and ones and + // in addition use the symmetry in entries to + // reduce the number of read operations. + const VectorizedArray * shape_values = this->data.shape_values.begin(); for (int i2=0; i2 val0, val1, res0, res1; + VectorizedArray val0, val1, in0, in1, res0, res1; if (dof_to_quad == true) { - val0 = this->data.shape_values[col]; - val1 = this->data.shape_values[nn-1-col]; + val0 = shape_values[col]; + val1 = shape_values[nn-1-col]; } else { - val0 = this->data.shape_values[col*n_q_points_1d]; - val1 = this->data.shape_values[(col+1)*n_q_points_1d-1]; + val0 = shape_values[col*n_q_points_1d]; + val1 = shape_values[(col+1)*n_q_points_1d-1]; } if (mid > 0) { - res0 = val0 * in[0]; - res1 = val1 * in[0]; - res0 += val1 * in[stride*(mm-1)]; - res1 += val0 * in[stride*(mm-1)]; + in0 = in[0]; + in1 = in[stride*(mm-1)]; + res0 = val0 * in0; + res1 = val1 * in0; + res0 += val1 * in1; + res1 += val0 * in1; for (int ind=1; inddata.shape_values[ind*n_q_points_1d+col]; - val1 = this->data.shape_values[ind*n_q_points_1d+nn-1-col]; + val0 = shape_values[ind*n_q_points_1d+col]; + val1 = shape_values[ind*n_q_points_1d+nn-1-col]; } else { - val0 = this->data.shape_values[col*n_q_points_1d+ind]; - val1 = this->data.shape_values[(col+1)*n_q_points_1d-1-ind]; + val0 = shape_values[col*n_q_points_1d+ind]; + val1 = shape_values[(col+1)*n_q_points_1d-1-ind]; } - res0 += val0 * in[stride*ind]; - res1 += val1 * in[stride*ind]; - res0 += val1 * in[stride*(mm-1-ind)]; - res1 += val0 * in[stride*(mm-1-ind)]; + in0 = in[stride*ind]; + in1 = in[stride*(mm-1-ind)]; + res0 += val0 * in0; + res1 += val1 * in0; + res0 += val1 * in1; + res1 += val0 * in1; } } else @@ -4997,7 +4797,7 @@ apply_values (const VectorizedArray input [], { if (mm % 2 == 1) { - val0 = this->data.shape_values[mid*n_q_points_1d+col]; + val0 = shape_values[mid*n_q_points_1d+col]; val1 = val0 * in[stride*mid]; res0 += val1; res1 += val1; @@ -5007,7 +4807,7 @@ apply_values (const VectorizedArray input [], { if (mm % 2 == 1 && nn % 2 == 0) { - val0 = this->data.shape_values[col*n_q_points_1d+mid]; + val0 = shape_values[col*n_q_points_1d+mid]; val1 = val0 * in[stride*mid]; res0 += val1; res1 += val1; @@ -5034,14 +4834,14 @@ apply_values (const VectorizedArray input [], else if (dof_to_quad == true && nn%2==1) { VectorizedArray res0; - VectorizedArray val0 = this->data.shape_values[n_cols]; + VectorizedArray val0 = shape_values[n_cols]; if (mid > 0) { res0 = in[0] + in[stride*(mm-1)]; res0 *= val0; for (int ind=1; inddata.shape_values[ind*n_q_points_1d+n_cols]; + val0 = shape_values[ind*n_q_points_1d+n_cols]; VectorizedArray val1 = in[stride*ind] + in[stride*(mm-1-ind)]; val1 *= val0; res0 += val1; @@ -5051,7 +4851,7 @@ apply_values (const VectorizedArray input [], res0 = VectorizedArray(); if (mm % 2 == 1) { - val0 = this->data.shape_values[mid*n_q_points_1d+n_cols]; + val0 = shape_values[mid*n_q_points_1d+n_cols]; res0 += val0 * in[stride*mid]; } if (add == false) @@ -5064,12 +4864,12 @@ apply_values (const VectorizedArray input [], VectorizedArray res0; if (mid > 0) { - VectorizedArray val0 = this->data.shape_values[n_cols*n_q_points_1d]; + VectorizedArray val0 = shape_values[n_cols*n_q_points_1d]; res0 = in[0] + in[stride*(mm-1)]; res0 *= val0; for (int ind=1; inddata.shape_values[n_cols*n_q_points_1d+ind]; + val0 = shape_values[n_cols*n_q_points_1d+ind]; VectorizedArray val1 = in[stride*ind] + in[stride*(mm-1-ind)]; val1 *= val0; res0 += val1; @@ -5114,14 +4914,14 @@ apply_values (const VectorizedArray input [], -template template inline void -FEEvaluation:: -apply_gradients (const VectorizedArray input [], - VectorizedArray output []) +FEEvaluation +::apply_gradients (const VectorizedArray in [], + VectorizedArray out []) { AssertIndexRange (direction, dim); const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, @@ -5131,50 +4931,74 @@ apply_gradients (const VectorizedArray input [], const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); - const int stride = ((direction > 0 ? nn : 1 ) * - (direction > 1 ? nn : 1)); + const int stride = Utilities::fixed_int_power::value; - const VectorizedArray * in = &input[0]; - VectorizedArray * out = &output[0]; + const VectorizedArray * shape_gradients = this->data.shape_gradients.begin(); for (int i2=0; i2 [-2.549 -1 0.549 ] + // [ 3.098 0 -3.098 ] + // [-0.549 1 2.549 ] + // Q3 --> [-4.315 -1.03 0.5 -0.44 ] + // [ 6.07 -1.44 -2.97 2.196 ] + // [-2.196 2.97 1.44 -6.07 ] + // [ 0.44 -0.5 1.03 4.315 ] + // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ] + // [10.111 -2.76 -2.667 2.066 -2.306 ] + // [-5.688 5.773 0 -5.773 5.688 ] + // [ 2.306 -2.066 2.667 2.76 -10.111 ] + // [-0.413 0.353 -0.333 -0.353 0.413 ] + // + // In these matrices, we want to use avoid + // computations involving zeros and ones and + // in addition use the symmetry in entries to + // reduce the number of read operations. for (int col=0; col val0, val1, res0, res1; + VectorizedArray val0, val1, in0, in1, res0, res1; if (dof_to_quad == true) { - val0 = this->data.shape_gradients[col]; - val1 = this->data.shape_gradients[nn-1-col]; + val0 = shape_gradients[col]; + val1 = shape_gradients[nn-1-col]; } else { - val0 = this->data.shape_gradients[col*n_q_points_1d]; - val1 = this->data.shape_gradients[(nn-col-1)*n_q_points_1d]; + val0 = shape_gradients[col*n_q_points_1d]; + val1 = shape_gradients[(nn-col-1)*n_q_points_1d]; } if (mid > 0) { - res0 = val0 * in[0]; - res1 = val1 * in[0]; - res0 -= val1 * in[stride*(mm-1)]; - res1 -= val0 * in[stride*(mm-1)]; + in0 = in[0]; + in1 = in[stride*(mm-1)]; + res0 = val0 * in0; + res1 = val1 * in0; + res0 -= val1 * in1; + res1 -= val0 * in1; for (int ind=1; inddata.shape_gradients[ind*n_q_points_1d+col]; - val1 = this->data.shape_gradients[ind*n_q_points_1d+nn-1-col]; + val0 = shape_gradients[ind*n_q_points_1d+col]; + val1 = shape_gradients[ind*n_q_points_1d+nn-1-col]; } else { - val0 = this->data.shape_gradients[col*n_q_points_1d+ind]; - val1 = this->data.shape_gradients[(nn-col-1)*n_q_points_1d+ind]; + val0 = shape_gradients[col*n_q_points_1d+ind]; + val1 = shape_gradients[(nn-col-1)*n_q_points_1d+ind]; } - res0 += val0 * in[stride*ind]; - res1 += val1 * in[stride*ind]; - res0 -= val1 * in[stride*(mm-1-ind)]; - res1 -= val0 * in[stride*(mm-1-ind)]; + in0 = in[stride*ind]; + in1 = in[stride*(mm-1-ind)]; + res0 += val0 * in0; + res1 += val1 * in0; + res0 -= val1 * in1; + res1 -= val0 * in1; } } else @@ -5182,9 +5006,9 @@ apply_gradients (const VectorizedArray input [], if (mm % 2 == 1) { if (dof_to_quad == true) - val0 = this->data.shape_gradients[mid*n_q_points_1d+col]; + val0 = shape_gradients[mid*n_q_points_1d+col]; else - val0 = this->data.shape_gradients[col*n_q_points_1d+mid]; + val0 = shape_gradients[col*n_q_points_1d+mid]; val1 = val0 * in[stride*mid]; res0 += val1; res1 -= val1; @@ -5204,17 +5028,17 @@ apply_gradients (const VectorizedArray input [], { VectorizedArray val0, res0; if (dof_to_quad == true) - val0 = this->data.shape_gradients[n_cols]; + val0 = shape_gradients[n_cols]; else - val0 = this->data.shape_gradients[n_cols*n_q_points_1d]; + val0 = shape_gradients[n_cols*n_q_points_1d]; res0 = in[0] - in[stride*(mm-1)]; res0 *= val0; for (int ind=1; inddata.shape_gradients[ind*n_q_points_1d+n_cols]; + val0 = shape_gradients[ind*n_q_points_1d+n_cols]; else - val0 = this->data.shape_gradients[n_cols*n_q_points_1d+ind]; + val0 = shape_gradients[n_cols*n_q_points_1d+ind]; VectorizedArray val1 = in[stride*ind] - in[stride*(mm-1-ind)]; val1 *= val0; res0 += val1; @@ -5261,14 +5085,14 @@ apply_gradients (const VectorizedArray input [], // same symmetry relations hold. However, it // is not possible to omit some values that // are zero for the values -template template inline void -FEEvaluation:: -apply_hessians (const VectorizedArray input [], - VectorizedArray output []) +FEEvaluation +::apply_hessians (const VectorizedArray in [], + VectorizedArray out []) { AssertIndexRange (direction, dim); const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, @@ -5278,18 +5102,15 @@ apply_hessians (const VectorizedArray input [], const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); - const int stride = ((direction > 0 ? nn : 1 ) * - (direction > 1 ? nn : 1)); + const int stride = Utilities::fixed_int_power::value; - const VectorizedArray * in = &input[0]; - VectorizedArray * out = &output[0]; for (int i2=0; i2 val0, val1, res0, res1; + VectorizedArray val0, val1, in0, in1, res0, res1; if (dof_to_quad == true) { val0 = this->data.shape_hessians[col]; @@ -5302,10 +5123,12 @@ apply_hessians (const VectorizedArray input [], } if (mid > 0) { - res0 = val0 * in[0]; - res1 = val1 * in[0]; - res0 += val1 * in[stride*(mm-1)]; - res1 += val0 * in[stride*(mm-1)]; + in0 = in[0]; + in1 = in[stride*(mm-1)]; + res0 = val0 * in0; + res1 = val1 * in0; + res0 += val1 * in1; + res1 += val0 * in1; for (int ind=1; ind input [], val0 = this->data.shape_hessians[col*n_q_points_1d+ind]; val1 = this->data.shape_hessians[(col+1)*n_q_points_1d-1-ind]; } - res0 += val0 * in[stride*ind]; - res1 += val1 * in[stride*ind]; - res0 += val1 * in[stride*(mm-1-ind)]; - res1 += val0 * in[stride*(mm-1-ind)]; + in0 = in[stride*ind]; + in1 = in[stride*(mm-1-ind)]; + res0 += val0 * in0; + res1 += val1 * in0; + res0 += val1 * in1; + res1 += val0 * in1; } } else @@ -5413,13 +5238,13 @@ apply_hessians (const VectorizedArray input [], } -/*----------------------- FEEvaluationGL -------------------------------*/ +/*------------------------- FEEvaluationGL ----------------------------------*/ -template +template inline -FEEvaluationGL:: -FEEvaluationGL (const MatrixFree &data_in, +FEEvaluationGL +::FEEvaluationGL (const MatrixFree &data_in, const unsigned int fe_no, const unsigned int quad_no) : @@ -5429,6 +5254,7 @@ FEEvaluationGL (const MatrixFree &data_in, std::string error_mess = "FEEvaluationGL not appropriate. It assumes:\n"; error_mess += " - identity operation for shape values\n"; error_mess += " - zero diagonal at interior points for gradients\n"; + error_mess += " - gradient equal to unity at element boundary\n"; error_mess += "Try FEEvaluation<...> instead!"; const double zero_tol = @@ -5451,16 +5277,21 @@ FEEvaluationGL (const MatrixFree &data_in, for (unsigned int i=1; idata.shape_gradients[i*n_points_1d+i][0])data.shape_gradients[n_points_1d-1][0]- + (n_points_1d%2==0 ? -1. : 1.)) < zero_tol, + ExcMessage (error_mess.c_str())); #endif } -template +template inline void -FEEvaluationGL:: -evaluate (bool evaluate_val,bool evaluate_grad,bool evaluate_lapl) +FEEvaluationGL +::evaluate (const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl) { Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); @@ -5476,6 +5307,9 @@ evaluate (bool evaluate_val,bool evaluate_grad,bool evaluate_lapl) this->values_quad_initialized = true; #endif } + // separate implementation here compared to + // the general case because the values are an + // identity operation if (evaluate_grad == true) { for(unsigned int comp=0;comp +template inline void -FEEvaluationGL:: -integrate (bool integrate_val, bool integrate_grad) +FEEvaluationGL +::integrate (const bool integrate_val, const bool integrate_grad) { Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); @@ -5634,13 +5468,13 @@ integrate (bool integrate_val, bool integrate_grad) -template +template template inline void -FEEvaluationGL:: -apply_gradients (const VectorizedArray input [], - VectorizedArray output []) +FEEvaluationGL +::apply_gradients (const VectorizedArray in [], + VectorizedArray out []) { AssertIndexRange (direction, dim); const int mm = fe_degree+1; @@ -5650,34 +5484,75 @@ apply_gradients (const VectorizedArray input [], const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); - const int stride = ((direction > 0 ? nn : 1 ) * - (direction > 1 ? nn : 1)); - - const VectorizedArray * in = &input[0]; - VectorizedArray * out = &output[0]; + const int stride = Utilities::fixed_int_power::value; + + // This loop specializes the application of + // the tensor product loop for Gauss-Lobatto + // elements which are symmetric about 0.5 just + // as the general class of elements treated by + // FEEvaluation, have diagonal shape matrices + // for the values and have the following + // gradient matrices (notice the zeros on the + // diagonal in the interior points, which is + // due to the construction of Legendre + // polynomials): + // Q2 --> [-3 -1 1 ] + // [ 4 0 -4 ] + // [-1 1 3 ] + // Q3 --> [-6 -1.618 0.618 -1 ] + // [ 8.09 0 -2.236 3.09 ] + // [-3.09 2.236 0 -8.09 ] + // [ 1 -0.618 1.618 6 ] + // Q4 --> [-10 -2.482 0.75 -0.518 1 ] + // [ 13.51 0 -2.673 1.528 -2.82 ] + // [-5.333 3.491 0 -3.491 5.333 ] + // [ 2.82 -1.528 2.673 0 -13.51 ] + // [-1 0.518 -0.75 2.482 10 ] for (int i2=0; i2 val0, val1, res0, res1; - if (dof_to_quad == true) - { - val0 = this->data.shape_gradients[col]; - val1 = this->data.shape_gradients[nn-1-col]; - } - else - { - val0 = this->data.shape_gradients[col*mm]; - val1 = this->data.shape_gradients[(nn-col-1)*mm]; - } + VectorizedArray val0, val1, in0, in1, res0, res1; if (mid > 0) { - res0 = val0 * in[0]; - res1 = val1 * in[0]; - res0 -= val1 * in[stride*(mm-1)]; - res1 -= val0 * in[stride*(mm-1)]; + if (dof_to_quad == true) + { + val0 = this->data.shape_gradients[col]; + val1 = this->data.shape_gradients[nn-1-col]; + } + else + { + val0 = this->data.shape_gradients[col*mm]; + val1 = this->data.shape_gradients[(nn-col-1)*mm]; + } + in0 = in[0]; + in1 = in[stride*(mm-1)]; + if (col == 0) + { + if ((mm+dof_to_quad)%2 == 1) + { + res0 = val0 * in0; + res1 = -in0; + res0 += in1; + res1 -= val0 * in1; + } + else + { + res0 = val0 * in0; + res0 -= in1; + res1 = in0; + res1 -= val0 * in1; + } + } + else + { + res0 = val0 * in0; + res1 = val1 * in0; + res0 -= val1 * in1; + res1 -= val0 * in1; + } for (int ind=1; ind input [], // at inner points, the gradient is zero for // ind==col + in0 = in[stride*ind]; + in1 = in[stride*(mm-1-ind)]; if (ind == col) { - res1 += val1 * in[stride*ind]; - res0 -= val1 * in[stride*(mm-1-ind)]; + res1 += val1 * in0; + res0 -= val1 * in1; } else { - res0 += val0 * in[stride*ind]; - res1 += val1 * in[stride*ind]; - res0 -= val1 * in[stride*(mm-1-ind)]; - res1 -= val0 * in[stride*(mm-1-ind)]; + res0 += val0 * in0; + res1 += val1 * in0; + res0 -= val1 * in1; + res1 -= val0 * in1; } } } diff --git a/deal.II/include/deal.II/matrix_free/mapping_info.h b/deal.II/include/deal.II/matrix_free/mapping_info.h index 1848193b73..d8913efc2b 100644 --- a/deal.II/include/deal.II/matrix_free/mapping_info.h +++ b/deal.II/include/deal.II/matrix_free/mapping_info.h @@ -261,14 +261,15 @@ namespace MatrixFreeFunctions * underlying the problem (constructed from a * 1D tensor product quadrature formula). */ - dealii::hp::QCollection quadrature_formula; + dealii::hp::QCollection quadrature; /** * The (dim-1)-dimensional quadrature formula - * underlying the problem (constructed from a - * 1D tensor product quadrature formula). + * corresponding to face evaluation + * (constructed from a 1D tensor product + * quadrature formula). */ - dealii::hp::QCollection quadrature_formula_faces; + dealii::hp::QCollection face_quadrature; /** * The number of quadrature points for the diff --git a/deal.II/include/deal.II/matrix_free/mapping_info.templates.h b/deal.II/include/deal.II/matrix_free/mapping_info.templates.h index c607df65e6..9d1626f45a 100644 --- a/deal.II/include/deal.II/matrix_free/mapping_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/mapping_info.templates.h @@ -130,11 +130,11 @@ namespace MatrixFreeFunctions template void MappingInfo::initialize - (const dealii::Triangulation &tria, + (const dealii::Triangulation &tria, const std::vector > &cells, const std::vector &active_fe_index, const Mapping &mapping, - const std::vector > &quad, + const std::vector > &quad, const UpdateFlags update_flags_input) { clear(); @@ -204,16 +204,16 @@ namespace MatrixFreeFunctions current_data.n_q_points_face.push_back (Utilities::fixed_power(n_q_points_1d[q])); - current_data.quadrature_formula.push_back + current_data.quadrature.push_back (Quadrature(quad[my_q][q])); - current_data.quadrature_formula_faces.push_back + current_data.face_quadrature.push_back (Quadrature(quad[my_q][q])); // set quadrature weights in vectorized form current_data.quadrature_weights[q].resize(n_q_points); for (unsigned int i=0; i 1) current_data.quad_index_conversion[q] = n_q_points; @@ -260,7 +260,7 @@ namespace MatrixFreeFunctions // finite element, so just hold a vector of // FEValues std::vector > > - fe_values (current_data.quadrature_formula.size()); + fe_values (current_data.quadrature.size()); UpdateFlags update_flags_feval = (update_flags & update_inverse_jacobians ? update_jacobians : update_default) | (update_flags & update_jacobian_grads ? update_jacobian_grads : update_default) | @@ -307,7 +307,7 @@ namespace MatrixFreeFunctions if (fe_values[fe_index].get() == 0) fe_values[fe_index].reset (new FEValues (mapping, dummy_fe, - current_data.quadrature_formula[fe_index], + current_data.quadrature[fe_index], update_flags_feval)); FEValues &fe_val = *fe_values[fe_index]; data.resize (n_q_points); @@ -912,8 +912,8 @@ namespace MatrixFreeFunctions memory += MemoryConsumption::memory_consumption (jacobians_grad_upper); memory += MemoryConsumption::memory_consumption (rowstart_q_points); memory += MemoryConsumption::memory_consumption (quadrature_points); - memory += MemoryConsumption::memory_consumption (quadrature_formula); - memory += MemoryConsumption::memory_consumption (quadrature_formula_faces); + memory += MemoryConsumption::memory_consumption (quadrature); + memory += MemoryConsumption::memory_consumption (face_quadrature); memory += MemoryConsumption::memory_consumption (quadrature_weights); memory += MemoryConsumption::memory_consumption (n_q_points); memory += MemoryConsumption::memory_consumption (n_q_points_face); diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.h b/deal.II/include/deal.II/matrix_free/matrix_free.h index 8c7a87b2e9..bacc960834 100644 --- a/deal.II/include/deal.II/matrix_free/matrix_free.h +++ b/deal.II/include/deal.II/matrix_free/matrix_free.h @@ -932,8 +932,16 @@ public: * given hp index. */ const Quadrature & - get_quad (const unsigned int quad_index = 0, - const unsigned int hp_active_fe_index = 0) const; + get_quadrature (const unsigned int quad_index = 0, + const unsigned int hp_active_fe_index = 0) const; + + /** + * Returns the quadrature rule for + * given hp index. + */ + const Quadrature & + get_face_quadrature (const unsigned int quad_index = 0, + const unsigned int hp_active_fe_index = 0) const; /** * Queries whether or not the @@ -1352,8 +1360,15 @@ MatrixFree::create_cell_subrange_hp const unsigned int degree, const unsigned int vector_component) const { - if (dof_info[vector_component].cell_active_fe_index.size() == 0) - return range; + AssertIndexRange (vector_component, dof_info.size()); + if (dof_info[vector_component].cell_active_fe_index.empty()) + { + AssertDimension (dof_info[vector_component].fe_index_conversion.size(),1); + if (dof_info[vector_component].fe_index_conversion[0].first == degree) + return range; + else + return std::pair (range.second,range.second); + } const unsigned int fe_index = dof_info[vector_component].fe_index_from_degree(degree); @@ -1691,12 +1706,25 @@ MatrixFree::get_shape_info (const unsigned int index_fe, template inline const Quadrature & -MatrixFree::get_quad (const unsigned int quad_index, - const unsigned int active_fe_index) const +MatrixFree::get_quadrature (const unsigned int quad_index, + const unsigned int active_fe_index) const +{ + AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size()); + return mapping_info.mapping_data_gen[quad_index]. + quadrature[active_fe_index]; +} + + + +template +inline +const Quadrature & +MatrixFree::get_face_quadrature (const unsigned int quad_index, + const unsigned int active_fe_index) const { AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size()); return mapping_info.mapping_data_gen[quad_index]. - quadrature_formula[active_fe_index]; + face_quadrature[active_fe_index]; } 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 675f5679ae..1fc7d26320 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 @@ -184,7 +184,12 @@ internal_reinit(const Mapping &mapping, // cell, Jacobian determinants, quadrature // points in real space, based on the ordering // of the cells determined in @p - // extract_local_to_global_indices. + // extract_local_to_global_indices. The + // algorithm assumes that the active FE index + // for the transformations is given the active + // FE index in the zeroth DoFHandler. TODO: + // how do things look like in the more general + // case? if(additional_data.initialize_mapping == true) { mapping_info.initialize (dof_handler[0]->get_tria(), cell_level_index, @@ -486,14 +491,15 @@ void MatrixFree::initialize_indices for (unsigned int f=0; fget_fe().size(); ++ind) dof_info[no].fe_index_conversion[ind] = std::pair(fe[ind].degree, fe[ind].dofs_per_cell); + if (fe.size() > 1) + dof_info[no].cell_active_fe_index.resize(n_active_cells, + numbers::invalid_unsigned_int); } else { @@ -502,6 +508,10 @@ void MatrixFree::initialize_indices &*dof_handlers.dof_handler[no] : &*dof_handlers.mg_dof_handler[no]; fes.push_back (&dofh->get_fe()); dof_info[no].max_fe_index = 1; + dof_info[no].fe_index_conversion.resize (1); + dof_info[no].fe_index_conversion[0] = + std::pair(fes.back()->degree, + fes.back()->dofs_per_cell); } lexicographic_inv[no].resize (fes.size()); @@ -654,8 +664,9 @@ void MatrixFree::initialize_indices cell_level_index[counter].first, cell_level_index[counter].second, dofh); - dof_info[no].cell_active_fe_index[counter] = - cell_it->active_fe_index(); + if (dofh->get_fe().size() > 1) + dof_info[no].cell_active_fe_index[counter] = + cell_it->active_fe_index(); 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, diff --git a/deal.II/source/numerics/matrix_free.inst.in b/deal.II/source/numerics/matrix_free.inst.in index 913f3220b7..5c9d3882e1 100644 --- a/deal.II/source/numerics/matrix_free.inst.in +++ b/deal.II/source/numerics/matrix_free.inst.in @@ -17,36 +17,36 @@ for (deal_II_dimension : DIMENSIONS) template class MatrixFree; // reinit for DoFHandler - template void MatrixFree::reinit + template void MatrixFree::internal_reinit (const Mapping &, const std::vector*> &, const std::vector &, const std::vector &, - const std::vector > &, + const std::vector > &, const MatrixFree::AdditionalData); - template void MatrixFree::reinit + template void MatrixFree::internal_reinit (const Mapping &, const std::vector*> &, const std::vector &, const std::vector &, - const std::vector > &, + const std::vector > &, const MatrixFree::AdditionalData); // reinit for MGDoFHandler - template void MatrixFree::reinit + template void MatrixFree::internal_reinit (const Mapping &, const std::vector*> &, const std::vector &, const std::vector &, - const std::vector > &, + const std::vector > &, const MatrixFree::AdditionalData); - template void MatrixFree::reinit + template void MatrixFree::internal_reinit (const Mapping &, const std::vector*> &, const std::vector &, const std::vector &, - const std::vector > &, + const std::vector > &, const MatrixFree::AdditionalData); template void MatrixFree:: diff --git a/tests/matrix_free/get_functions_common.h b/tests/matrix_free/get_functions_common.h index 122db07661..f942ffd026 100644 --- a/tests/matrix_free/get_functions_common.h +++ b/tests/matrix_free/get_functions_common.h @@ -42,12 +42,10 @@ template ::n_vectors; - MatrixFreeTest(const MatrixFree &data_in): data (data_in), fe_val (data.get_dof_handler().get_fe(), - Quadrature(data.get_quad(0)), + Quadrature(data.get_quadrature(0)), update_values | update_gradients | update_hessians) {}; @@ -55,7 +53,7 @@ class MatrixFreeTest const Mapping &mapping): data (data_in), fe_val (mapping, data.get_dof_handler().get_fe(), - Quadrature(data.get_quad(0)), + Quadrature(data.get_quadrature(0)), update_values | update_gradients | update_hessians) {}; @@ -70,7 +68,6 @@ class MatrixFreeTest const Vector &src, const std::pair &cell_range) const { - typedef VectorizedArray vector_t; FEEvaluation fe_eval (data); std::vector reference_values (fe_eval.n_q_points); @@ -198,9 +195,6 @@ template class MatrixFreeTest { public: - typedef VectorizedArray vector_t; - static const std::size_t n_vectors = VectorizedArray::n_vectors; - MatrixFreeTest(const MatrixFree &) {}; diff --git a/tests/matrix_free/get_functions_gl.cc b/tests/matrix_free/get_functions_gl.cc index 3d80aa9182..62a7d65f2c 100644 --- a/tests/matrix_free/get_functions_gl.cc +++ b/tests/matrix_free/get_functions_gl.cc @@ -23,8 +23,6 @@ template class MatrixFreeTestGL : public MatrixFreeTest { public: - typedef VectorizedArray vector_t; - static const std::size_t n_vectors = VectorizedArray::n_array_elements; MatrixFreeTestGL(const MatrixFree &data, const Mapping &mapping): diff --git a/tests/matrix_free/get_functions_multife.cc b/tests/matrix_free/get_functions_multife.cc index 7ecb02f458..dce1bb24e1 100644 --- a/tests/matrix_free/get_functions_multife.cc +++ b/tests/matrix_free/get_functions_multife.cc @@ -38,17 +38,15 @@ template vector_t; typedef std::vector > VectorType; - static const std::size_t n_vectors = VectorizedArray::n_array_elements; MatrixFreeTest(const MatrixFree &data_in): data (data_in), fe_val0 (data.get_dof_handler(0).get_fe(), - Quadrature(data.get_quad(0)), + Quadrature(data.get_quadrature(0)), update_values | update_gradients | update_hessians), fe_val1 (data.get_dof_handler(1).get_fe(), - Quadrature(data.get_quad(1)), + Quadrature(data.get_quadrature(1)), update_values | update_gradients | update_hessians) {}; diff --git a/tests/matrix_free/get_functions_multife2.cc b/tests/matrix_free/get_functions_multife2.cc index 61499d7c9b..8c26693f5f 100644 --- a/tests/matrix_free/get_functions_multife2.cc +++ b/tests/matrix_free/get_functions_multife2.cc @@ -40,20 +40,18 @@ template vector_t; typedef std::vector > VectorType; - static const std::size_t n_vectors = VectorizedArray::n_array_elements; MatrixFreeTest(const MatrixFree &data_in): data (data_in), fe_val0 (data.get_dof_handler(0).get_fe(), - Quadrature(data.get_quad(0)), + Quadrature(data.get_quadrature(0)), update_values | update_gradients | update_hessians), fe_val1 (data.get_dof_handler(1).get_fe(), - Quadrature(data.get_quad(1)), + Quadrature(data.get_quadrature(1)), update_values | update_gradients | update_hessians), fe_val2 (data.get_dof_handler(2).get_fe(), - Quadrature(data.get_quad(1)), + Quadrature(data.get_quadrature(1)), update_values | update_gradients | update_hessians) {}; diff --git a/tests/matrix_free/get_functions_q_hierarchical.cc b/tests/matrix_free/get_functions_q_hierarchical.cc index b1ca17475b..9e2042debf 100644 --- a/tests/matrix_free/get_functions_q_hierarchical.cc +++ b/tests/matrix_free/get_functions_q_hierarchical.cc @@ -25,9 +25,6 @@ template class MatrixFreeTestGen : public MatrixFreeTest { public: - typedef VectorizedArray vector_t; - static const std::size_t n_vectors = VectorizedArray::n_array_elements; - MatrixFreeTestGen(const MatrixFree &data, const Mapping &mapping): MatrixFreeTest(data, mapping) diff --git a/tests/matrix_free/get_functions_rect.cc b/tests/matrix_free/get_functions_rect.cc index 5cb26618e4..86b2ec2dc0 100644 --- a/tests/matrix_free/get_functions_rect.cc +++ b/tests/matrix_free/get_functions_rect.cc @@ -51,7 +51,7 @@ void test () tria.begin(tria.n_levels()-1)->set_refine_flag(); tria.last()->set_refine_flag(); tria.execute_coarsening_and_refinement(); - tria.refine_global (4-dim); + tria.refine_global (1); FE_Q fe (fe_degree); DoFHandler dof (tria); diff --git a/tests/matrix_free/get_functions_variants.cc b/tests/matrix_free/get_functions_variants.cc index 5dd4c9d91c..68035a30b8 100644 --- a/tests/matrix_free/get_functions_variants.cc +++ b/tests/matrix_free/get_functions_variants.cc @@ -37,9 +37,7 @@ template class MatrixFreeTest { public: - typedef VectorizedArray vector_t; typedef Vector VectorType; - static const std::size_t n_vectors = VectorizedArray::n_array_elements; MatrixFreeTest(const MatrixFree &data_in): data (data_in) @@ -121,7 +119,7 @@ operator () (const MatrixFree &data, // FEEvaluations. Those are tested in other // functions and seen as reference here for (unsigned int q=0; q::n_array_elements; ++j) { errors[0] += std::fabs(fe_eval.get_value(q)[j]- fe_eval2.get_value(q)[j]); diff --git a/tests/matrix_free/get_values_plain.cc b/tests/matrix_free/get_values_plain.cc index 709f985288..5c97f82bb2 100644 --- a/tests/matrix_free/get_values_plain.cc +++ b/tests/matrix_free/get_values_plain.cc @@ -38,8 +38,6 @@ template ::n_array_elements; - MatrixFreeTest(const MatrixFree &data_in): data (data_in) {}; @@ -52,8 +50,6 @@ class MatrixFreeTest const Vector &src, const std::pair &cell_range) const { - typedef VectorizedArray vector_t; - const unsigned int n_vectors = sizeof(vector_t)/sizeof(Number); FEEvaluation fe_eval (data); FEEvaluation fe_eval_plain (data); for(unsigned int cell=cell_range.first;cell::n_array_elements; ++j) { error += std::fabs(fe_eval.get_dof_value(i)[j]- fe_eval_plain.get_dof_value(i)[j]); diff --git a/tests/matrix_free/integrate_functions.cc b/tests/matrix_free/integrate_functions.cc index 12dce33f76..73cf3a85da 100644 --- a/tests/matrix_free/integrate_functions.cc +++ b/tests/matrix_free/integrate_functions.cc @@ -36,14 +36,12 @@ template class MatrixFreeTest { public: - typedef VectorizedArray vector_t; typedef std::vector*> VectorType; - static const std::size_t n_vectors = VectorizedArray::n_array_elements; MatrixFreeTest(const MatrixFree &data_in): data (data_in), fe_val (data.get_dof_handler().get_fe(), - Quadrature(data.get_quad(0)), + Quadrature(data.get_quadrature(0)), update_values | update_gradients | update_JxW_values) {}; @@ -83,8 +81,8 @@ operator () (const MatrixFree &data, FEEvaluation fe_eval (data); const unsigned int n_q_points = fe_eval.n_q_points; const unsigned int dofs_per_cell = fe_eval.dofs_per_cell; - AlignedVector values (n_q_points); - AlignedVector gradients (dim*n_q_points); + AlignedVector > values (n_q_points); + AlignedVector > gradients (dim*n_q_points); std::vector dof_indices (dofs_per_cell); for(unsigned int cell=cell_range.first;cell &data, for (unsigned int q=0; q submit (false); + Tensor<1,dim,VectorizedArray > submit (false); for (unsigned int d=0; d class MatrixFreeTest { public: - typedef VectorizedArray vector_t; typedef std::vector > VectorType; - static const std::size_t n_vectors = VectorizedArray::n_array_elements; MatrixFreeTest(const MatrixFree &data_in): data (data_in), fe_val0 (data.get_dof_handler(0).get_fe(), - Quadrature(data.get_quad(0)), + Quadrature(data.get_quadrature(0)), update_values | update_gradients | update_JxW_values), fe_val01 (data.get_dof_handler(0).get_fe(), - Quadrature(data.get_quad(1)), + Quadrature(data.get_quadrature(1)), update_values | update_gradients | update_JxW_values), fe_val1 (data.get_dof_handler(1).get_fe(), - Quadrature(data.get_quad(1)), + Quadrature(data.get_quadrature(1)), update_values | update_gradients | update_JxW_values) {}; @@ -92,10 +90,10 @@ operator () (const MatrixFree &data, const unsigned int n_q_points1 = fe_eval1.n_q_points; const unsigned int dofs_per_cell0 = fe_eval0.dofs_per_cell; const unsigned int dofs_per_cell1 = fe_eval1.dofs_per_cell; - AlignedVector values0 (n_q_points0); - AlignedVector gradients0 (dim*n_q_points0); - AlignedVector values1 (n_q_points1); - AlignedVector gradients1 (dim*n_q_points1); + AlignedVector > values0 (n_q_points0); + AlignedVector > gradients0 (dim*n_q_points0); + AlignedVector > values1 (n_q_points1); + AlignedVector > gradients1 (dim*n_q_points1); std::vector dof_indices0 (dofs_per_cell0); std::vector dof_indices1 (dofs_per_cell1); for(unsigned int cell=cell_range.first;cell &data, for (unsigned int q=0; q submit (false); + Tensor<1,dim,VectorizedArray > submit (false); for (unsigned int d=0; d &data, for (unsigned int q=0; q submit (false); + Tensor<1,dim,VectorizedArray > submit (false); for (unsigned int d=0; d &data, for (unsigned int q=0; q submit (false); + Tensor<1,dim,VectorizedArray > submit (false); for (unsigned int d=0; d #include +#include "create_mesh.h" + #include std::ofstream logfile("integrate_functions_multife2/output"); @@ -38,20 +40,18 @@ template class MatrixFreeTest { public: - typedef VectorizedArray vector_t; typedef std::vector > VectorType; - static const std::size_t n_vectors = VectorizedArray::n_array_elements; MatrixFreeTest(const MatrixFree &data_in): data (data_in), fe_val0 (data.get_dof_handler(0).get_fe(), - Quadrature(data.get_quad(0)), + Quadrature(data.get_quadrature(0)), update_values | update_gradients | update_JxW_values), fe_val01 (data.get_dof_handler(0).get_fe(), - Quadrature(data.get_quad(1)), + Quadrature(data.get_quadrature(1)), update_values | update_gradients | update_JxW_values), fe_val1 (data.get_dof_handler(1).get_fe(), - Quadrature(data.get_quad(1)), + Quadrature(data.get_quadrature(1)), update_values | update_gradients | update_JxW_values) {}; @@ -93,10 +93,10 @@ operator () (const MatrixFree &data, const unsigned int n_q_points1 = fe_eval1.n_q_points; const unsigned int dofs_per_cell0 = fe_eval0.dofs_per_cell; const unsigned int dofs_per_cell1 = fe_eval1.dofs_per_cell; - AlignedVector values0 (n_q_points0); - AlignedVector gradients0 (dim*n_q_points0); - AlignedVector values1 (n_q_points1); - AlignedVector gradients1 (dim*n_q_points1); + AlignedVector > values0 (n_q_points0); + AlignedVector > gradients0 (dim*n_q_points0); + AlignedVector > values1 (n_q_points1); + AlignedVector > gradients1 (dim*n_q_points1); std::vector dof_indices0 (dofs_per_cell0); std::vector dof_indices1 (dofs_per_cell1); for(unsigned int cell=cell_range.first;cell &data, for (unsigned int q=0; q submit (false); + Tensor<1,dim,VectorizedArray > submit (false); for (unsigned int d=0; d &data, for (unsigned int q=0; q submit (false); + Tensor<1,dim,VectorizedArray > submit (false); for (unsigned int d=0; d &data, for (unsigned int q=0; q submit (false); + Tensor<1,dim,VectorizedArray > submit (false); for (unsigned int d=0; d tria; - GridGenerator::hyper_ball (tria); - static const HyperBallBoundary boundary; - tria.set_boundary (0, boundary); - 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(); + create_mesh (tria); + tria.begin_active ()->set_refine_flag(); tria.execute_coarsening_and_refinement(); + typename Triangulation::active_cell_iterator cell, endc; cell = tria.begin_active (); + endc = tria.end(); for (; cell!=endc; ++cell) - if (cell->center().norm()<0.2) + if (cell->center().norm()<0.5) cell->set_refine_flag(); tria.execute_coarsening_and_refinement(); - 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(); + tria.refine_global(1); cell = tria.begin_active (); - for (unsigned int i=0; i<7-2*dim; ++i) + for (unsigned int i=0; i<10-3*dim; ++i) { cell = tria.begin_active (); + endc = tria.end(); unsigned int counter = 0; for (; cell!=endc; ++cell, ++counter) if (counter % (7-i) == 0) @@ -324,7 +319,7 @@ int main () deallog << std::setprecision (3); { - deallog.threshold_double(1.e-12); + deallog.threshold_double(1.e-11); deallog.push("2d"); test<2,1,double>(); test<2,2,double>(); diff --git a/tests/matrix_free/matrix_vector_06.cc b/tests/matrix_free/matrix_vector_06.cc index fc0d926027..1708f87255 100644 --- a/tests/matrix_free/matrix_vector_06.cc +++ b/tests/matrix_free/matrix_vector_06.cc @@ -22,6 +22,8 @@ std::ofstream logfile("matrix_vector_06/output"); template void test () { + if (fe_degree > 1) + return; Triangulation tria; create_mesh (tria); tria.begin_active ()->set_refine_flag(); @@ -33,8 +35,6 @@ void test () if (cell->center().norm()<0.5) cell->set_refine_flag(); tria.execute_coarsening_and_refinement(); - 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(); diff --git a/tests/matrix_free/matrix_vector_06/cmp/generic b/tests/matrix_free/matrix_vector_06/cmp/generic index c8d1cd1e9f..e9e176a4c1 100644 --- a/tests/matrix_free/matrix_vector_06/cmp/generic +++ b/tests/matrix_free/matrix_vector_06/cmp/generic @@ -2,12 +2,6 @@ DEAL:2d::Testing FE_Q<2>(1) DEAL:2d::Norm of difference: 0 DEAL:2d:: -DEAL:2d::Testing FE_Q<2>(2) -DEAL:2d::Norm of difference: 0 -DEAL:2d:: 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_08.cc b/tests/matrix_free/matrix_vector_08.cc index 3511799f87..1f468c9af8 100644 --- a/tests/matrix_free/matrix_vector_08.cc +++ b/tests/matrix_free/matrix_vector_08.cc @@ -20,6 +20,9 @@ std::ofstream logfile("matrix_vector_08/output"); template void test () { + if (fe_degree > 1) + return; + Triangulation tria; create_mesh (tria, 1e20); tria.begin_active ()->set_refine_flag(); @@ -31,8 +34,6 @@ void test () if (cell->center().norm()<0.5*1e20) cell->set_refine_flag(); tria.execute_coarsening_and_refinement(); - 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(); diff --git a/tests/matrix_free/matrix_vector_08/cmp/generic b/tests/matrix_free/matrix_vector_08/cmp/generic index c8d1cd1e9f..e9e176a4c1 100644 --- a/tests/matrix_free/matrix_vector_08/cmp/generic +++ b/tests/matrix_free/matrix_vector_08/cmp/generic @@ -2,12 +2,6 @@ DEAL:2d::Testing FE_Q<2>(1) DEAL:2d::Norm of difference: 0 DEAL:2d:: -DEAL:2d::Testing FE_Q<2>(2) -DEAL:2d::Norm of difference: 0 -DEAL:2d:: 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_09.cc b/tests/matrix_free/matrix_vector_09.cc index 20ca29bab8..556f6940fe 100644 --- a/tests/matrix_free/matrix_vector_09.cc +++ b/tests/matrix_free/matrix_vector_09.cc @@ -20,6 +20,9 @@ std::ofstream logfile("matrix_vector_09/output"); template void test () { + if (fe_degree > 1) + return; + Triangulation tria; create_mesh (tria, 1e-20); tria.begin_active ()->set_refine_flag(); @@ -31,12 +34,9 @@ void test () if (cell->center().norm()<0.5*1e-20) cell->set_refine_flag(); tria.execute_coarsening_and_refinement(); - 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(); - tria.refine_global(1); cell = tria.begin_active (); for (unsigned int i=0; i<10-3*dim; ++i) { diff --git a/tests/matrix_free/matrix_vector_09/cmp/generic b/tests/matrix_free/matrix_vector_09/cmp/generic index c8d1cd1e9f..e9e176a4c1 100644 --- a/tests/matrix_free/matrix_vector_09/cmp/generic +++ b/tests/matrix_free/matrix_vector_09/cmp/generic @@ -2,12 +2,6 @@ DEAL:2d::Testing FE_Q<2>(1) DEAL:2d::Norm of difference: 0 DEAL:2d:: -DEAL:2d::Testing FE_Q<2>(2) -DEAL:2d::Norm of difference: 0 -DEAL:2d:: 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_hp.cc b/tests/matrix_free/matrix_vector_hp.cc index 0577fe6a37..fa7e807efa 100644 --- a/tests/matrix_free/matrix_vector_hp.cc +++ b/tests/matrix_free/matrix_vector_hp.cc @@ -25,9 +25,6 @@ template class MatrixFreeTestHP { public: - typedef VectorizedArray vector_t; - static const std::size_t n_vectors = VectorizedArray::n_array_elements; - MatrixFreeTestHP(const MatrixFree &data_in): data (data_in) {}; diff --git a/tests/matrix_free/matrix_vector_stokes.cc b/tests/matrix_free/matrix_vector_stokes.cc index f1caacaf72..5c81f3effa 100644 --- a/tests/matrix_free/matrix_vector_stokes.cc +++ b/tests/matrix_free/matrix_vector_stokes.cc @@ -7,7 +7,7 @@ // this function tests the correctness of the implementation of matrix free // matrix-vector products by comparing with the result of deal.II sparse -// matrix. The mesh uses a hypercube mesh with no hanging nodes and no other +// matrix. No hanging nodes and no other // constraints for a vector-valued problem (stokes equations). #include "../tests.h" @@ -34,6 +34,8 @@ std::ofstream logfile("matrix_vector_stokes/output"); #include #include +#include "create_mesh.h" + #include #include #include @@ -113,23 +115,11 @@ template void test () { Triangulation triangulation; - { - std::vector subdivisions (dim, 1); - subdivisions[0] = 4; - - const Point bottom_left = (dim == 2 ? - Point(-2,-1) : - Point(-2,0,-1)); - const Point top_right = (dim == 2 ? - Point(2,0) : - Point(2,1,0)); - - GridGenerator::subdivided_hyper_rectangle (triangulation, - subdivisions, - bottom_left, - top_right); - } - triangulation.refine_global (4-dim); + 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 fe_p (fe_degree); diff --git a/tests/matrix_free/matrix_vector_stokes_noflux.cc b/tests/matrix_free/matrix_vector_stokes_noflux.cc index e9c8909c60..6c58266926 100644 --- a/tests/matrix_free/matrix_vector_stokes_noflux.cc +++ b/tests/matrix_free/matrix_vector_stokes_noflux.cc @@ -341,11 +341,9 @@ int main () test<2,1>(); test<2,2>(); test<2,3>(); - test<2,4>(); deallog.pop(); deallog.push("3d"); test<3,1>(); - test<3,2>(); deallog.pop(); } } diff --git a/tests/matrix_free/matrix_vector_stokes_noflux/cmp/generic b/tests/matrix_free/matrix_vector_stokes_noflux/cmp/generic index afab4527f1..0434059632 100644 --- a/tests/matrix_free/matrix_vector_stokes_noflux/cmp/generic +++ b/tests/matrix_free/matrix_vector_stokes_noflux/cmp/generic @@ -8,9 +8,5 @@ DEAL:2d::Verification fe degree 2: 0 DEAL:2d:: DEAL:2d::Verification fe degree 3: 0 DEAL:2d:: -DEAL:2d::Verification fe degree 4: 0 -DEAL:2d:: DEAL:3d::Verification fe degree 1: 0 DEAL:3d:: -DEAL:3d::Verification fe degree 2: 0 -DEAL:3d:: diff --git a/tests/matrix_free/quadrature_points.cc b/tests/matrix_free/quadrature_points.cc index ba030b88e9..99c3b98af8 100644 --- a/tests/matrix_free/quadrature_points.cc +++ b/tests/matrix_free/quadrature_points.cc @@ -63,9 +63,9 @@ void test () } double error_points = 0, abs_points = 0; - const unsigned int n_cells = mf_data.get_size_info().n_macro_cells; + const unsigned int n_cells = mf_data.n_macro_cells(); FEEvaluation fe_eval (mf_data); - FEValues fe_values (mapping, fe, mf_data.get_quad(), + FEValues fe_values (mapping, fe, mf_data.get_quadrature(), update_quadrature_points); typedef VectorizedArray vector_t; diff --git a/tests/matrix_free/thread_correctness_hp.cc b/tests/matrix_free/thread_correctness_hp.cc index 84b58e609a..9b4940e2ab 100644 --- a/tests/matrix_free/thread_correctness_hp.cc +++ b/tests/matrix_free/thread_correctness_hp.cc @@ -26,9 +26,6 @@ template class MatrixFreeTestHP { public: - typedef VectorizedArray vector_t; - static const std::size_t n_vectors = VectorizedArray::n_array_elements; - MatrixFreeTestHP(const MatrixFree &data_in): data (data_in) {}; @@ -132,11 +129,11 @@ void do_test (const unsigned int parallel_option) dof.distribute_dofs(fe_collection); ConstraintMatrix constraints; DoFTools::make_hanging_node_constraints (dof, - constraints); + constraints); VectorTools::interpolate_boundary_values (dof, - 0, - ZeroFunction(), - constraints); + 0, + ZeroFunction(), + constraints); constraints.close (); //std::cout << "Number of cells: " << dof.get_tria().n_active_cells() << std::endl; @@ -151,48 +148,55 @@ void do_test (const unsigned int parallel_option) mf_data.reinit (dof, constraints, quadrature_collection_mf, data); MatrixFreeTestHP mf (mf_data); - MatrixFree mf_data_par; - if (parallel_option == 0) - { - data.tasks_parallel_scheme = - MatrixFree::AdditionalData::partition_partition; - deallog << "Parallel option partition/partition" << std::endl; - } - else + // test different block sizes, starting from + // auto setting (= 0) + for (unsigned int block_size = 0; block_size < 5; ++block_size) { - data.tasks_parallel_scheme = - MatrixFree::AdditionalData::partition_color; - deallog << "Parallel option partition/color" << std::endl; - } - data.tasks_block_size = 1; - mf_data_par.reinit (dof, constraints, quadrature_collection_mf, data); - MatrixFreeTestHP mf_par(mf_data_par); + deallog.push ("blk_" + Utilities::int_to_string(block_size,1)); + MatrixFree mf_data_par; + if (parallel_option == 0) + { + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_partition; + deallog << "Parallel option partition/partition" << std::endl; + } + else + { + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + deallog << "Parallel option partition/color" << std::endl; + } + data.tasks_block_size = 1; + mf_data_par.reinit (dof, constraints, quadrature_collection_mf, data); + MatrixFreeTestHP mf_par(mf_data_par); // fill a right hand side vector with random // numbers in unconstrained degrees of freedom - Vector src (dof.n_dofs()); - Vector result_ref(src), result_mf (src); + Vector src (dof.n_dofs()); + Vector result_ref(src), result_mf (src); - for (unsigned int i=0; i