From: kronbichler Date: Mon, 17 Feb 2014 11:22:08 +0000 (+0000) Subject: Make FEEvaluation collaborate with some data from FEValues. Cleanup of some internal... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=12289cdc4dc745dfc07aee55b66804b84c71c8e8;p=dealii-svn.git Make FEEvaluation collaborate with some data from FEValues. Cleanup of some internal data structures to make this possible. Implement FEEvaluationDGP for DGP elements by truncating the tensor product. Shorten a few matrix-free stress tests in debug mode. git-svn-id: https://svn.dealii.org/trunk@32496 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/matrix_free/evaluated_geometry.h b/deal.II/include/deal.II/matrix_free/evaluated_geometry.h new file mode 100644 index 0000000000..f4ec60c484 --- /dev/null +++ b/deal.II/include/deal.II/matrix_free/evaluated_geometry.h @@ -0,0 +1,259 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2011 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +#ifndef __deal2__matrix_free_evaluated_geometry_h +#define __deal2__matrix_free_evaluated_geometry_h + + +#include +#include +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +/** + * The class makes FEEvaluation with the mapping information generated by + * FEValues. + */ +template +class EvaluatedGeometry : public Subscriptor +{ +public: + /** + * Constructor, similar to FEValues. + */ + EvaluatedGeometry (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags); + + /** + * Constructor. Instead of providing a mapping, use MappingQ1. + */ + EvaluatedGeometry (const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags); + + /** + * Initialize with the given cell iterator. + */ + template + void reinit(ITERATOR &cell); + + /** + * Return a vector of inverse transpose Jacobians. For compatibility with + * FEEvaluation, it returns tensors of vectorized arrays, even though all + * components are equal. + */ + const AlignedVector > >& + get_inverse_jacobians() const; + + /** + * Return a vector of quadrature weights times the Jacobian determinant + * (JxW). For compatibility with FEEvaluation, it returns tensors of + * vectorized arrays, even though all components are equal. + */ + const AlignedVector >& + get_JxW_values() const; + + /** + * Return a vector of quadrature points in real space on the given + * cell. For compatibility with FEEvaluation, it returns tensors of + * vectorized arrays, even though all components are equal. + */ + const AlignedVector > >& + get_quadrature_points() const; + + /** + * Return a vector of quadrature points in real space on the given + * cell. For compatibility with FEEvaluation, it returns tensors of + * vectorized arrays, even though all components are equal. + */ + const AlignedVector > >& + get_normal_vectors() const; + + /** + * Return a reference to 1D quadrature underlying this object. + */ + const Quadrature<1>& + get_quadrature () const; + +private: + /** + * An underlying FEValues object that performs the (scalar) evaluation. + */ + FEValues fe_values; + + /** + * Get 1D quadrature formula to be used for reinitializing shape info. + */ + const Quadrature<1> quadrature_1d; + + /** + * Inverse Jacobians, stored in vectorized array form. + */ + AlignedVector > > inverse_jacobians; + + /** + * Stored Jacobian determinants and quadrature weights + */ + AlignedVector > jxw_values; + + /** + * Stored quadrature points + */ + AlignedVector > > quadrature_points; + + /** + * Stored normal vectors (for face integration) + */ + AlignedVector > > normal_vectors; +}; + + +/*----------------------- Inline functions ----------------------------------*/ + +template +inline +EvaluatedGeometry::EvaluatedGeometry (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags) + : + fe_values(mapping, fe, Quadrature(quadrature), + internal::MatrixFreeFunctions::MappingInfo::compute_update_flags(update_flags)), + quadrature_1d(quadrature), + inverse_jacobians(fe_values.get_quadrature().size()), + jxw_values(fe_values.get_quadrature().size()), + quadrature_points(fe_values.get_quadrature().size()), + normal_vectors(fe_values.get_quadrature().size()) +{ + Assert(!(fe_values.get_update_flags() & update_jacobian_grads), + ExcNotImplemented()); +} + + + +template +inline +EvaluatedGeometry::EvaluatedGeometry (const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags) + : + fe_values(fe, Quadrature(quadrature), + internal::MatrixFreeFunctions::MappingInfo::compute_update_flags(update_flags)), + quadrature_1d(quadrature), + inverse_jacobians(fe_values.get_quadrature().size()), + jxw_values(fe_values.get_quadrature().size()), + quadrature_points(fe_values.get_quadrature().size()), + normal_vectors(fe_values.get_quadrature().size()) +{ + Assert(!(fe_values.get_update_flags() & update_jacobian_grads), + ExcNotImplemented()); +} + + + +template +template +inline +void +EvaluatedGeometry::reinit(ITERATOR &cell) +{ + fe_values.reinit(cell); + for (unsigned int q=0; q +inline +const AlignedVector > >& +EvaluatedGeometry::get_inverse_jacobians() const +{ + return inverse_jacobians; +} + + + +template +inline +const AlignedVector > >& +EvaluatedGeometry::get_normal_vectors() const +{ + return normal_vectors; +} + + + +template +inline +const AlignedVector > >& +EvaluatedGeometry::get_quadrature_points() const +{ + return quadrature_points; +} + + + +template +inline +const AlignedVector >& +EvaluatedGeometry::get_JxW_values() const +{ + return jxw_values; +} + + + +template +inline +const Quadrature<1>& +EvaluatedGeometry::get_quadrature() const +{ + return quadrature_1d; +} + + +#ifndef DOXYGEN + + +#endif // ifndef DOXYGEN + + +DEAL_II_NAMESPACE_CLOSE + +#endif 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 38ddc1bd49..d39ba8dce3 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2011 - 2013 by the deal.II authors +// Copyright (C) 2011 - 2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -24,7 +24,10 @@ #include #include #include +#include #include +#include +#include DEAL_II_NAMESPACE_OPEN @@ -43,44 +46,24 @@ namespace parallel namespace internal { DeclException0 (ExcAccessToUninitializedField); - - template - void do_evaluate (FEEval &, const bool, const bool, const bool, int2type<1>); - template - void do_evaluate (FEEval &, const bool, const bool, const bool, int2type<2>); - template - void do_evaluate (FEEval &, const bool, const bool, const bool, int2type<3>); - template - void do_integrate (FEEval &, const bool, const bool, int2type<1>); - template - void do_integrate (FEEval &, const bool, const bool, int2type<2>); - template - void do_integrate (FEEval &, const bool, const bool, int2type<3>); } /** * This is the base class for the FEEvaluation classes. This class is a base - * class and needs usually not be called in user code. Use one of the derived - * classes FEEvaluationGeneral, FEEvaluation or FEEvaluationGL instead. It - * implements a reinit method that is used to set pointers so that operations - * on quadrature points can be performed quickly, access functions to vectors - * for the @p read_dof_values, @p set_dof_values, and @p - * distributed_local_to_global functions, as well as methods to access values - * and gradients of finite element functions. + * class and needs usually not be called in user code. It does not have any + * public constructor. Use one of the derived classes FEEvaluationGeneral, + * FEEvaluation or FEEvaluationGL instead. It implements a reinit method that + * is used to set pointers so that operations on quadrature points can be + * performed quickly, access functions to vectors for the @p read_dof_values, + * @p set_dof_values, and @p distributed_local_to_global functions, as well as + * methods to access values and gradients of finite element functions. * - * This class has five template arguments: + * This class has three template arguments: * * @param dim Dimension in which this class is to be used * - * @param dofs_per_cell Number of degrees of freedom of the FE per cell, - * usually (fe_degree+1)^dim for elements based on a tensor - * product - * - * @param n_q_points Number of points in the quadrature formula, usually - * (fe_degree+1)^dim for tensor-product quadrature formulas - * * @param n_components Number of vector components when solving a system of * PDEs. If the same operation is applied to several * components of a PDE (e.g. a vector Laplace equation), they @@ -114,6 +97,14 @@ public: */ void reinit (const unsigned int cell); + /** + * Initializes the operation pointer to the current cell. This is a reinit + * call similar to FEValues where the necessary information is computed on + * the fly. + */ + template + void reinit (const ITERATOR &cell_iterator); + /** * For the transformation information stored in MappingInfo, this function * returns the index which belongs to the current cell as specified in @p @@ -527,6 +518,26 @@ protected: const unsigned int dofs_per_cell, const unsigned int n_q_points); + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. The user has to provide a structure of type EvaluatedGeometry + * and a DoFHandler in order to allow for reading out the finite element + * data. It uses the data provided by dof_handler.get_fe(). If the element + * is vector-valued, the optional argument allows to specify the index of + * the base element (as long as the element is primitive, non-primitive are + * not supported currently). + * + * With this initialization, no call to a reinit method of this + * class. Instead, it is enough if the geometry is initialized to a given + * cell iterator. Moreover, beware that a kernel using this method does not + * vectorize over several elements (which is most efficient for vector + * operations), but only possibly within the element if the + * evaluate/integrate routines are combined (e.g. for matrix assembly). + */ + FEEvaluationBase (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element = 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, @@ -548,11 +559,12 @@ protected: void read_dof_values_plain (const VectorType *src_data[]); /** - * Internal data fields that store the values. Since all array lengths are - * known at compile time and since they are rarely more than a few - * kilobytes, allocate them on the stack. This makes it possible to cheaply - * set up a FEEvaluation object and write thread-safe programs by letting - * each thread own a private object of this type. + * Internal data fields that store the values. Derived classes will know the + * length of all arrays at compile time and allocate the memory on the + * stack. This makes it possible to cheaply set up a FEEvaluation object and + * write thread-safe programs by letting each thread own a private object of + * this type. In this base class, only pointers to the actual data are + * stored. * * This field stores the values for local degrees of freedom (e.g. after * reading out from a vector but before applying unit cell transformations @@ -609,24 +621,24 @@ protected: const unsigned int active_quad_index; /** - * Stores a reference to the underlying data. + * Stores a pointer to the underlying data. */ - const MatrixFree &matrix_info; + const MatrixFree *matrix_info; /** - * Stores a reference to the underlying DoF indices and constraint + * Stores a pointer to the underlying DoF indices and constraint * description for the component specified at construction. Also contained * in matrix_info, but it simplifies code if we store a reference to it. */ - const internal::MatrixFreeFunctions::DoFInfo &dof_info; + const internal::MatrixFreeFunctions::DoFInfo *dof_info; /** - * Stores a reference to the underlying transformation data from unit to + * Stores a pointer to the underlying transformation data from unit to * real cells for the given quadrature formula specified at construction. * Also contained in matrix_info, but it simplifies code if we store a * reference to it. */ - const internal::MatrixFreeFunctions::MappingInfo &mapping_info; + const internal::MatrixFreeFunctions::MappingInfo *mapping_info; /** * Stores a reference to the unit cell data, i.e., values, gradients and @@ -634,7 +646,7 @@ protected: * product. Also contained in matrix_info, but it simplifies code if we * store a reference to it. */ - const internal::MatrixFreeFunctions::ShapeInfo &data; + std_cxx1x::shared_ptr > data; /** * A pointer to the Cartesian Jacobian information of the present cell. Only @@ -740,6 +752,16 @@ protected: * stared. Used to control exceptions when uninitialized data is used. */ bool gradients_quad_submitted; + + /** + * Geometry data generated by FEValues on the fly. + */ + SmartPointer > evaluated_geometry; + + /** + * A pointer to the underlying DoFHandler. + */ + const DoFHandler *dof_handler; }; @@ -776,6 +798,27 @@ protected: const unsigned int quad_no, const unsigned int dofs_per_cell, const unsigned int n_q_points); + + + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. The user has to provide a structure of type EvaluatedGeometry + * and a DoFHandler in order to allow for reading out the finite element + * data. It uses the data provided by dof_handler.get_fe(). If the element + * is vector-valued, the optional argument allows to specify the index of + * the base element (as long as the element is primitive, non-primitive are + * not supported currently). + * + * With this initialization, no call to a reinit method of this + * class. Instead, it is enough if the geometry is initialized to a given + * cell iterator. Moreover, beware that a kernel using this method does not + * vectorize over several elements (which is most efficient for vector + * operations), but only possibly within the element if the + * evaluate/integrate routines are combined (e.g. for matrix assembly). + */ + FEEvaluationAccess (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element = 0); }; @@ -898,6 +941,26 @@ protected: const unsigned int quad_no, const unsigned int dofs_per_cell, const unsigned int n_q_points); + + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. The user has to provide a structure of type EvaluatedGeometry + * and a DoFHandler in order to allow for reading out the finite element + * data. It uses the data provided by dof_handler.get_fe(). If the element + * is vector-valued, the optional argument allows to specify the index of + * the base element (as long as the element is primitive, non-primitive are + * not supported currently). + * + * With this initialization, no call to a reinit method of this + * class. Instead, it is enough if the geometry is initialized to a given + * cell iterator. Moreover, beware that a kernel using this method does not + * vectorize over several elements (which is most efficient for vector + * operations), but only possibly within the element if the + * evaluate/integrate routines are combined (e.g. for matrix assembly). + */ + FEEvaluationAccess (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element = 0); }; @@ -1029,6 +1092,26 @@ protected: const unsigned int quad_no, const unsigned int dofs_per_cell, const unsigned int n_q_points); + + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. The user has to provide a structure of type EvaluatedGeometry + * and a DoFHandler in order to allow for reading out the finite element + * data. It uses the data provided by dof_handler.get_fe(). If the element + * is vector-valued, the optional argument allows to specify the index of + * the base element (as long as the element is primitive, non-primitive are + * not supported currently). + * + * With this initialization, no call to a reinit method of this + * class. Instead, it is enough if the geometry is initialized to a given + * cell iterator. Moreover, beware that a kernel using this method does not + * vectorize over several elements (which is most efficient for vector + * operations), but only possibly within the element if the + * evaluate/integrate routines are combined (e.g. for matrix assembly). + */ + FEEvaluationAccess (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element = 0); }; @@ -1096,6 +1179,26 @@ public: const unsigned int fe_no = 0, const unsigned int quad_no = 0); + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. The user has to provide a structure of type EvaluatedGeometry + * and a DoFHandler in order to allow for reading out the finite element + * data. It uses the data provided by dof_handler.get_fe(). If the element + * is vector-valued, the optional argument allows to specify the index of + * the base element (as long as the element is primitive, non-primitive are + * not supported currently). + * + * With this initialization, no call to a reinit method of this + * class. Instead, it is enough if the geometry is initialized to a given + * cell iterator. Moreover, beware that a kernel using this method does not + * vectorize over several elements (which is most efficient for vector + * operations), but only possibly within the element if the + * evaluate/integrate routines are combined (e.g. for matrix assembly). + */ + FEEvaluationGeneral (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element = 0); + /** * Evaluates the function values, the gradients, and the Laplacians of the * FE function given at the DoF values in the input vector at the quadrature @@ -1123,8 +1226,6 @@ public: Point > quadrature_point (const unsigned int q_point) const; -protected: - /** * Internal function that applies the function values of the tensor product * in a given coordinate direction (first template argument), from @@ -1162,13 +1263,7 @@ protected: void apply_hessians (const VectorizedArray in [], VectorizedArray out []); - /** - * Friend declaration. - */ - template friend void - internal::do_evaluate (FEEval &, const bool, const bool, const bool, internal::int2type); - template friend void - internal::do_integrate (FEEval &, const bool, const bool, internal::int2type); +protected: /** * Internally stored variables for the different data fields. @@ -1185,6 +1280,16 @@ protected: * functions that make it much faster (between 5 and 500, depending on the * polynomial order). * + * This class can be used in two different ways. The first way is to + * initialize it from a MatrixFree object that caches everything related to + * the degrees of freedom and the mapping information. This way, it is + * possible to use vectorization for applying a vector operation for several + * cells at once. The second form of usage is to initialize it from geometry + * information generated by FEValues, which is stored in the class + * EvaluatedGeometry. Here, the operations can only work on a single cell, but + * possibly be vectorized by combining several operations (e.g. when + * performing matrix assembly). + * * This class is a specialization of FEEvaluationGeneral designed for standard * FE_Q or FE_DGQ elements and quadrature points symmetric around 0.5 (like * Gauss quadrature), and hence the most common situation. Note that many of @@ -1244,6 +1349,26 @@ public: const unsigned int fe_no = 0, const unsigned int quad_no = 0); + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. The user has to provide a structure of type EvaluatedGeometry + * and a DoFHandler in order to allow for reading out the finite element + * data. It uses the data provided by dof_handler.get_fe(). If the element + * is vector-valued, the optional argument allows to specify the index of + * the base element (as long as the element is primitive, non-primitive are + * not supported currently). + * + * With this initialization, no call to a reinit method of this + * class. Instead, it is enough if the geometry is initialized to a given + * cell iterator. Moreover, beware that a kernel using this method does not + * vectorize over several elements (which is most efficient for vector + * operations), but only possibly within the element if the + * evaluate/integrate routines are combined (e.g. for matrix assembly). + */ + FEEvaluation (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element = 0); + /** * Evaluates the function values, the gradients, and the Laplacians of the * FE function given at the DoF values in the input vector at the quadrature @@ -1266,8 +1391,6 @@ public: 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 coordinate direction (first template argument), from @@ -1305,17 +1428,16 @@ protected: void apply_hessians (const VectorizedArray in [], VectorizedArray out []); +protected: VectorizedArray shape_val_evenodd[fe_degree+1][(n_q_points_1d+1)/2]; VectorizedArray shape_gra_evenodd[fe_degree+1][(n_q_points_1d+1)/2]; VectorizedArray shape_hes_evenodd[fe_degree+1][(n_q_points_1d+1)/2]; +private: /** - * Friend declarations. + * Fills the fields shapve_???_evenodd, called in the constructor. */ - template friend void - internal::do_evaluate (FEEval &, const bool, const bool, const bool, internal::int2type); - template friend void - internal::do_integrate (FEEval &, const bool, const bool, internal::int2type); + void compute_even_odd_factors(); }; @@ -1381,6 +1503,26 @@ public: const unsigned int fe_no = 0, const unsigned int quad_no = 0); + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. The user has to provide a structure of type EvaluatedGeometry + * and a DoFHandler in order to allow for reading out the finite element + * data. It uses the data provided by dof_handler.get_fe(). If the element + * is vector-valued, the optional argument allows to specify the index of + * the base element (as long as the element is primitive, non-primitive are + * not supported currently). + * + * With this initialization, no call to a reinit method of this + * class. Instead, it is enough if the geometry is initialized to a given + * cell iterator. Moreover, beware that a kernel using this method does not + * vectorize over several elements (which is most efficient for vector + * operations), but only possibly within the element if the + * evaluate/integrate routines are combined (e.g. for matrix assembly). + */ + FEEvaluationGL (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element = 0); + /** * Evaluates the function values, the gradients, and the Hessians of the FE * function given at the DoF values in the input vector at the quadrature @@ -1403,7 +1545,6 @@ public: void integrate (const bool integrate_val, const bool integrate_grad); -protected: /** * Internal function that applies the gradient operation of the tensor * product in a given coordinate direction (first template argument), from @@ -1419,12 +1560,148 @@ protected: +namespace internal +{ + namespace MatrixFreeFunctions + { + // a helper function to compute the number of DoFs of a DGP element at compile + // time, depending on the degree + template + struct DGP_dofs_per_cell + { + // this division is always without remainder + static const unsigned int value = + (DGP_dofs_per_cell::value * (degree+dim)) / dim; + }; + + // base specialization: 1d elements have 'degree+1' degrees of freedom + template + struct DGP_dofs_per_cell<1,degree> + { + static const unsigned int value = degree+1; + }; + } +} + + + +/** + * The class that provides all functions necessary to evaluate functions at + * quadrature points and cell integrations. In functionality, this class is + * similar to FEValues, however, it includes a lot of specialized + * functions that make it much faster (between 5 and 500 times as fast, + * depending on the polynomial order). Access to the data fields is provided + * through functionality in the class FEEvaluationAccess. + * + * This class is an extension of FEEvaluationGeneral to work with elements of + * complete polynomial degree p, FE_DGP. In this case, the polynomial basis is + * a truncated tensor product, so the evaluate and integrate routines use a + * truncation. + * + * @author Martin Kronbichler, 2014 + */ +template +class FEEvaluationDGP : + public FEEvaluationGeneral +{ +public: + 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 = internal::MatrixFreeFunctions::DGP_dofs_per_cell::value; + static const unsigned int n_q_points = BaseClass::n_q_points; + + /** + * Constructor. Takes all data stored in MatrixFree. If applied to problems + * with more than one finite element or more than one quadrature formula + * selected during construction of @p matrix_free, @p fe_no and @p quad_no + * allow to select the appropriate components. + */ + FEEvaluationDGP (const MatrixFree &matrix_free, + const unsigned int fe_no = 0, + const unsigned int quad_no = 0); + + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. The user has to provide a structure of type EvaluatedGeometry + * and a DoFHandler in order to allow for reading out the finite element + * data. It uses the data provided by dof_handler.get_fe(). If the element + * is vector-valued, the optional argument allows to specify the index of + * the base element (as long as the element is primitive, non-primitive are + * not supported currently). + * + * With this initialization, no call to a reinit method of this + * class. Instead, it is enough if the geometry is initialized to a given + * cell iterator. Moreover, beware that a kernel using this method does not + * vectorize over several elements (which is most efficient for vector + * operations), but only possibly within the element if the + * evaluate/integrate routines are combined (e.g. for matrix assembly). + */ + FEEvaluationDGP (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element = 0); + + /** + * Evaluates the function values, the gradients, and the Hessians of the FE + * function given at the DoF values in the input vector at the quadrature + * points of the unit cell. The function arguments specify which parts shall + * actually be computed. Needs to be called before the functions @p + * get_value(), @p get_gradient() or @p get_laplacian give useful + * information (unless these values have been set manually). + */ + void evaluate (const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl = false); + + /** + * This function takes the values and/or gradients that are stored on + * quadrature points, tests them by all the basis functions/gradients on the + * cell and performs the cell integration. The two function arguments @p + * integrate_val and @p integrate_grad are used to enable/disable some of + * values or gradients. + */ + void integrate (const bool integrate_val, + const bool integrate_grad); +}; + /*----------------------- Inline functions ----------------------------------*/ #ifndef DOXYGEN +namespace internal +{ + namespace MatrixFreeFunctions + { + // a small class that gives control over the delete behavior of + // std::shared_ptr: we need to disable it when we initialize a pointer + // from another structure. + template + struct DummyDeleter + { + DummyDeleter (const bool do_delete = false) + : + do_delete(do_delete) + {} + + void operator () (CLASS *pointer) + { + if (do_delete) + delete pointer; + } + + const bool do_delete; + }; + } +} + + + /*----------------------- FEEvaluationBase ----------------------------------*/ template @@ -1443,16 +1720,18 @@ FEEvaluationBase active_quad_index (data_in.get_mapping_info(). mapping_data_gen[quad_no_in]. quad_index_from_n_q_points(n_q_points)), - matrix_info (data_in), - dof_info (data_in.get_dof_info(fe_no_in)), - mapping_info (data_in.get_mapping_info()), - data (data_in.get_shape_info + matrix_info (&data_in), + dof_info (&data_in.get_dof_info(fe_no_in)), + mapping_info (&data_in.get_mapping_info()), + data (&data_in.get_shape_info (fe_no_in, quad_no_in, active_fe_index, - active_quad_index)), + active_quad_index), + internal::MatrixFreeFunctions::DummyDeleter + >(false)), cartesian_data (0), jacobian (0), J_value (0), - quadrature_weights (mapping_info.mapping_data_gen[quad_no]. + quadrature_weights (mapping_info->mapping_data_gen[quad_no]. quadrature_weights[active_quad_index].begin()), quadrature_points (0), jacobian_grad (0), @@ -1470,14 +1749,14 @@ FEEvaluationBase for (unsigned int d=0; d<(dim*dim+dim)/2; ++d) hessians_quad[c][d] = 0; } - Assert (matrix_info.mapping_initialized() == true, + Assert (matrix_info->mapping_initialized() == true, ExcNotInitialized()); - AssertDimension (matrix_info.get_size_info().vectorization_length, + AssertDimension (matrix_info->get_size_info().vectorization_length, VectorizedArray::n_array_elements); - AssertDimension (data.dofs_per_cell, - dof_info.dofs_per_cell[active_fe_index]/n_fe_components); - AssertDimension (data.n_q_points, - mapping_info.mapping_data_gen[quad_no].n_q_points[active_quad_index]); + AssertDimension (data->dofs_per_cell, + dof_info->dofs_per_cell[active_fe_index]/n_fe_components); + AssertDimension (data->n_q_points, + mapping_info->mapping_data_gen[quad_no].n_q_points[active_quad_index]); Assert (n_fe_components == 1 || n_components == 1 || n_components == n_fe_components, @@ -1492,66 +1771,116 @@ FEEvaluationBase +template +inline +FEEvaluationBase +::FEEvaluationBase (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler_in, + const unsigned int base_element) + : + quad_no (-1), + n_fe_components (n_components_), + active_fe_index (-1), + active_quad_index (-1), + matrix_info (0), + dof_info (0), + mapping_info (0), + data (new internal::MatrixFreeFunctions::ShapeInfo(geometry.get_quadrature(), dof_handler_in.get_fe(), base_element)), + cartesian_data (0), + jacobian (geometry.get_inverse_jacobians().begin()), + J_value (geometry.get_JxW_values().begin()), + quadrature_weights (0), + quadrature_points (geometry.get_quadrature_points().begin()), + jacobian_grad (0), + jacobian_grad_upper(0), + cell (0), + cell_type (internal::MatrixFreeFunctions::general), + cell_data_number (0), + evaluated_geometry (&geometry), + dof_handler (&dof_handler_in) +{ + for (unsigned int c=0; cget_fe().element_multiplicity(base_element) == 1 || + dof_handler->get_fe().element_multiplicity(base_element) >= n_components_, + ExcMessage("The underlying element must at least contain as many " + "components as requested by this class")); +} + + + template inline void 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) ? - dof_info.cell_active_fe_index[cell_in] : 0), + Assert (evaluated_geometry == 0, ExcMessage("FEEvaluation was initialized without a matrix-free object. Integer indexing is not possible")); + if (evaluated_geometry != 0) + return; + Assert (dof_info != 0, ExcNotInitialized()); + Assert (mapping_info != 0, ExcNotInitialized()); + AssertIndexRange (cell_in, dof_info->row_starts.size()-1); + AssertDimension (((dof_info->cell_active_fe_index.size() > 0) ? + dof_info->cell_active_fe_index[cell_in] : 0), active_fe_index); cell = cell_in; - cell_type = mapping_info.get_cell_type(cell); - cell_data_number = mapping_info.get_cell_data_index(cell); + cell_type = mapping_info->get_cell_type(cell); + cell_data_number = mapping_info->get_cell_data_index(cell); - if (mapping_info.quadrature_points_initialized == true) + if (mapping_info->quadrature_points_initialized == true) { - AssertIndexRange (cell_data_number, mapping_info. + AssertIndexRange (cell_data_number, mapping_info-> mapping_data_gen[quad_no].rowstart_q_points.size()); - const unsigned int index = mapping_info.mapping_data_gen[quad_no]. + const unsigned int index = mapping_info->mapping_data_gen[quad_no]. rowstart_q_points[cell]; - AssertIndexRange (index, mapping_info.mapping_data_gen[quad_no]. + AssertIndexRange (index, mapping_info->mapping_data_gen[quad_no]. quadrature_points.size()); quadrature_points = - &mapping_info.mapping_data_gen[quad_no].quadrature_points[index]; + &mapping_info->mapping_data_gen[quad_no].quadrature_points[index]; } if (cell_type == internal::MatrixFreeFunctions::cartesian) { - cartesian_data = &mapping_info.cartesian_data[cell_data_number].first; - J_value = &mapping_info.cartesian_data[cell_data_number].second; + cartesian_data = &mapping_info->cartesian_data[cell_data_number].first; + J_value = &mapping_info->cartesian_data[cell_data_number].second; } else if (cell_type == internal::MatrixFreeFunctions::affine) { - jacobian = &mapping_info.affine_data[cell_data_number].first; - J_value = &mapping_info.affine_data[cell_data_number].second; + jacobian = &mapping_info->affine_data[cell_data_number].first; + J_value = &mapping_info->affine_data[cell_data_number].second; } else { - const unsigned int rowstart = mapping_info. + const unsigned int rowstart = mapping_info-> mapping_data_gen[quad_no].rowstart_jacobians[cell_data_number]; - AssertIndexRange (rowstart, mapping_info. + AssertIndexRange (rowstart, mapping_info-> mapping_data_gen[quad_no].jacobians.size()); jacobian = - &mapping_info.mapping_data_gen[quad_no].jacobians[rowstart]; - if (mapping_info.JxW_values_initialized == true) + &mapping_info->mapping_data_gen[quad_no].jacobians[rowstart]; + if (mapping_info->JxW_values_initialized == true) { - AssertIndexRange (rowstart, mapping_info. + AssertIndexRange (rowstart, mapping_info-> mapping_data_gen[quad_no].JxW_values.size()); - J_value = &(mapping_info.mapping_data_gen[quad_no]. + J_value = &(mapping_info->mapping_data_gen[quad_no]. JxW_values[rowstart]); } - if (mapping_info.second_derivatives_initialized == true) + if (mapping_info->second_derivatives_initialized == true) { - AssertIndexRange(rowstart, mapping_info. + AssertIndexRange(rowstart, mapping_info-> mapping_data_gen[quad_no].jacobians_grad_diag.size()); - jacobian_grad = &mapping_info.mapping_data_gen[quad_no]. + jacobian_grad = &mapping_info->mapping_data_gen[quad_no]. jacobians_grad_diag[rowstart]; - AssertIndexRange(rowstart, mapping_info. + AssertIndexRange(rowstart, mapping_info-> mapping_data_gen[quad_no].jacobians_grad_upper.size()); - jacobian_grad_upper = &mapping_info.mapping_data_gen[quad_no]. + jacobian_grad_upper = &mapping_info->mapping_data_gen[quad_no]. jacobians_grad_upper[rowstart]; } } @@ -1828,22 +2157,24 @@ FEEvaluationBase // into the local data field or write local data into the vector. Certain // operations are no-ops for the given use case. - Assert (matrix_info.indices_initialized() == true, + Assert (matrix_info != 0, ExcNotInitialized()); + Assert (dof_info != 0, ExcNotInitialized()); + Assert (matrix_info->indices_initialized() == true, ExcNotInitialized()); Assert (cell != numbers::invalid_unsigned_int, 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 unsigned int *dof_indices = dof_info->begin_indices(cell); const std::pair *indicators = - dof_info.begin_indicators(cell); + dof_info->begin_indicators(cell); const std::pair *indicators_end = - dof_info.end_indicators(cell); + dof_info->end_indicators(cell); unsigned int ind_local = 0; - const unsigned int dofs_per_cell = this->data.dofs_per_cell; + const unsigned int dofs_per_cell = this->data->dofs_per_cell; - const unsigned int n_irreg_components_filled = dof_info.row_starts[cell][2]; + const unsigned int n_irreg_components_filled = dof_info->row_starts[cell][2]; const bool at_irregular_cell = n_irreg_components_filled > 0; // scalar case (or case when all components have the same degrees of freedom @@ -1853,7 +2184,7 @@ FEEvaluationBase const unsigned int n_local_dofs = VectorizedArray::n_array_elements * dofs_per_cell; for (unsigned int comp=0; comp value[comp]); const Number *data_val = - matrix_info.constraint_pool_begin(indicators->second); + matrix_info->constraint_pool_begin(indicators->second); const Number *end_pool = - matrix_info.constraint_pool_end(indicators->second); + matrix_info->constraint_pool_end(indicators->second); for ( ; data_val != end_pool; ++data_val, ++dof_indices) for (unsigned int comp=0; comp { // no constraint at all: loop bounds are known, compiler can // unroll without checks - AssertDimension (dof_info.end_indices(cell)-dof_indices, + AssertDimension (dof_info->end_indices(cell)-dof_indices, static_cast(n_local_dofs)); for (unsigned int j=0; j value[comp]); const Number *data_val = - matrix_info.constraint_pool_begin(indicators->second); + matrix_info->constraint_pool_begin(indicators->second); const Number *end_pool = - matrix_info.constraint_pool_end(indicators->second); + matrix_info->constraint_pool_end(indicators->second); for ( ; data_val != end_pool; ++data_val, ++dof_indices) for (unsigned int comp=0; comp } for (; ind_localend_indices(cell), ExcInternalError()); // non-constrained case: copy the data from the global vector, @@ -2005,7 +2336,7 @@ FEEvaluationBase // the first component, then all entries to the second one, and so // on. This is ensured by the way MatrixFree reads out the indices. { - internal::check_vector_compatibility (*src[0], dof_info); + 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; @@ -2031,9 +2362,9 @@ FEEvaluationBase operation.pre_constraints (local_data[ind_local], value); const Number *data_val = - matrix_info.constraint_pool_begin(indicators->second); + matrix_info->constraint_pool_begin(indicators->second); const Number *end_pool = - matrix_info.constraint_pool_end(indicators->second); + matrix_info->constraint_pool_end(indicators->second); for ( ; data_val != end_pool; ++data_val, ++dof_indices) operation.process_constraint (*dof_indices, *data_val, @@ -2047,14 +2378,14 @@ FEEvaluationBase for (; ind_localend_indices(cell), ExcInternalError()); } else { // no constraint at all: loop bounds are known, compiler can // unroll without checks - AssertDimension (dof_info.end_indices(cell)-dof_indices, + AssertDimension (dof_info->end_indices(cell)-dof_indices, static_cast(n_local_dofs)); for (unsigned int j=0; j operation.pre_constraints (local_data[ind_local], value); const Number *data_val = - matrix_info.constraint_pool_begin(indicators->second); + matrix_info->constraint_pool_begin(indicators->second); const Number *end_pool = - matrix_info.constraint_pool_end(indicators->second); + matrix_info->constraint_pool_end(indicators->second); for ( ; data_val != end_pool; ++data_val, ++dof_indices) operation.process_constraint (*dof_indices, *data_val, @@ -2113,7 +2444,7 @@ FEEvaluationBase } for (; ind_localend_indices(cell), ExcInternalError()); // non-constrained case: copy the data from the global vector, @@ -2436,18 +2767,20 @@ FEEvaluationBase { // this is different from the other three operations because we do not use // constraints here, so this is a separate function. - Assert (matrix_info.indices_initialized() == true, + Assert (matrix_info != 0, ExcNotInitialized()); + Assert (dof_info != 0, ExcNotInitialized()); + Assert (matrix_info->indices_initialized() == true, ExcNotInitialized()); Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - Assert (dof_info.store_plain_indices == true, ExcNotInitialized()); + 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_plain(cell); - const unsigned int dofs_per_cell = this->data.dofs_per_cell; + const unsigned int *dof_indices = dof_info->begin_indices_plain(cell); + const unsigned int dofs_per_cell = this->data->dofs_per_cell; - const unsigned int n_irreg_components_filled = dof_info.row_starts[cell][2]; + const unsigned int n_irreg_components_filled = dof_info->row_starts[cell][2]; const bool at_irregular_cell = n_irreg_components_filled > 0; // scalar case (or case when all components have the same degrees of freedom @@ -2457,7 +2790,7 @@ FEEvaluationBase const unsigned int n_local_dofs = VectorizedArray::n_array_elements * dofs_per_cell; for (unsigned int comp=0; comp // the first component, then all entries to the second one, and so // on. This is ensured by the way MatrixFree reads out the indices. { - internal::check_vector_compatibility (*src[0], dof_info); + 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; @@ -2655,7 +2988,7 @@ Tensor<1,n_components_,VectorizedArray > FEEvaluationBase ::get_dof_value (const unsigned int dof) const { - AssertIndexRange (dof, this->data.dofs_per_cell); + AssertIndexRange (dof, this->data->dofs_per_cell); Tensor<1,n_components_,VectorizedArray > return_value (false); for (unsigned int comp=0; compvalues_dofs[comp][dof]; @@ -2672,7 +3005,7 @@ FEEvaluationBase { Assert (this->values_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); Tensor<1,n_components_,VectorizedArray > return_value (false); for (unsigned int comp=0; compvalues_quad[comp][q_point]; @@ -2689,7 +3022,7 @@ FEEvaluationBase { Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > grad_out (false); @@ -2789,7 +3122,7 @@ FEEvaluationBase { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); Tensor<2,dim,VectorizedArray > hessian_out [n_components]; @@ -2828,7 +3161,7 @@ FEEvaluationBase // cell with general Jacobian else if (this->cell_type == internal::MatrixFreeFunctions::general) { - Assert (this->mapping_info.second_derivatives_initialized == true, + Assert (this->mapping_info->second_derivatives_initialized == true, ExcNotInitialized()); const Tensor<2,dim,VectorizedArray > &jac = jacobian[q_point]; const Tensor<2,dim,VectorizedArray > &jac_grad = jacobian_grad[q_point]; @@ -2914,7 +3247,7 @@ FEEvaluationBase { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > hessian_out (false); @@ -2930,7 +3263,7 @@ FEEvaluationBase // cell with general Jacobian else if (this->cell_type == internal::MatrixFreeFunctions::general) { - Assert (this->mapping_info.second_derivatives_initialized == true, + Assert (this->mapping_info->second_derivatives_initialized == true, ExcNotInitialized()); const Tensor<2,dim,VectorizedArray > &jac = jacobian[q_point]; const Tensor<2,dim,VectorizedArray > &jac_grad = jacobian_grad[q_point]; @@ -2992,7 +3325,7 @@ FEEvaluationBase { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); Tensor<1,n_components_,VectorizedArray > laplacian_out (false); const Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > hess_diag = get_hessian_diagonal(q_point); @@ -3017,7 +3350,7 @@ FEEvaluationBase #ifdef DEBUG this->dof_values_initialized = true; #endif - AssertIndexRange (dof, this->data.dofs_per_cell); + AssertIndexRange (dof, this->data->dofs_per_cell); for (unsigned int comp=0; compvalues_dofs[comp][dof] = val_in[comp]; } @@ -3033,7 +3366,7 @@ FEEvaluationBase { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); this->values_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::general) @@ -3062,7 +3395,7 @@ FEEvaluationBase { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); this->gradients_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::cartesian) @@ -3108,7 +3441,7 @@ FEEvaluationBase Tensor<1,n_components_,VectorizedArray > return_value (false); for (unsigned int comp=0; compvalues_quad[comp][0]; - const unsigned int n_q_points = this->data.n_q_points; + const unsigned int n_q_points = this->data->n_q_points; for (unsigned int q=1; qvalues_quad[comp][q]; @@ -3135,6 +3468,18 @@ FEEvaluationAccess +template +inline +FEEvaluationAccess +::FEEvaluationAccess (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element) + : + FEEvaluationBase (geometry, dof_handler, base_element) +{} + + + /*-------------------- FEEvaluationAccess scalar ----------------------------*/ @@ -3154,13 +3499,25 @@ FEEvaluationAccess +template +inline +FEEvaluationAccess +::FEEvaluationAccess (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element) + : + FEEvaluationBase (geometry, dof_handler, base_element) +{} + + + template inline VectorizedArray FEEvaluationAccess ::get_dof_value (const unsigned int dof) const { - AssertIndexRange (dof, this->data.dofs_per_cell); + AssertIndexRange (dof, this->data->dofs_per_cell); return this->values_dofs[0][dof]; } @@ -3174,7 +3531,7 @@ FEEvaluationAccess { Assert (this->values_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); return this->values_quad[0][q_point]; } @@ -3191,7 +3548,7 @@ FEEvaluationAccess Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); Tensor<1,dim,VectorizedArray > grad_out (false); @@ -3262,7 +3619,7 @@ FEEvaluationAccess { #ifdef DEBUG this->dof_values_initialized = true; - AssertIndexRange (dof, this->data.dofs_per_cell); + AssertIndexRange (dof, this->data->dofs_per_cell); #endif this->values_dofs[0][dof] = val_in; } @@ -3278,7 +3635,7 @@ FEEvaluationAccess { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); this->values_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::general) @@ -3304,7 +3661,7 @@ FEEvaluationAccess { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); this->gradients_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::cartesian) @@ -3366,6 +3723,18 @@ FEEvaluationAccess +template +inline +FEEvaluationAccess +::FEEvaluationAccess (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element) + : + FEEvaluationBase (geometry, dof_handler, base_element) +{} + + + template inline Tensor<2,dim,VectorizedArray > @@ -3385,7 +3754,7 @@ FEEvaluationAccess { Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); VectorizedArray divergence; @@ -3491,7 +3860,7 @@ FEEvaluationAccess { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); return BaseClass::get_hessian_diagonal (q_point); } @@ -3506,7 +3875,7 @@ FEEvaluationAccess { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); return BaseClass::get_hessian(q_point); } @@ -3546,7 +3915,7 @@ FEEvaluationAccess { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); this->gradients_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::cartesian) @@ -3596,7 +3965,7 @@ FEEvaluationAccess // that saves some operations #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data.n_q_points); + AssertIndexRange (q_point, this->data->n_q_points); this->gradients_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::cartesian) @@ -3711,8 +4080,11 @@ FEEvaluationGeneral #ifdef DEBUG // print error message when the dimensions do not match. Propose a possible // fix - if (dofs_per_cell != this->data.dofs_per_cell || - n_q_points != this->data.n_q_points) + if ((dofs_per_cell != this->data->dofs_per_cell && + internal::MatrixFreeFunctions::DGP_dofs_per_cell::value != + this->data->dofs_per_cell) + || + n_q_points != this->data->n_q_points) { std::string message = "-------------------------------------------------------\n"; @@ -3729,22 +4101,22 @@ FEEvaluationGeneral // points unsigned int proposed_dof_comp = numbers::invalid_unsigned_int, proposed_quad_comp = numbers::invalid_unsigned_int; - if (dofs_per_cell == this->matrix_info.get_dof_info(fe_no).dofs_per_cell[this->active_fe_index]) + if (dofs_per_cell == this->matrix_info->get_dof_info(fe_no).dofs_per_cell[this->active_fe_index]) proposed_dof_comp = fe_no; else - for (unsigned int no=0; nomatrix_info.n_components(); ++no) - if (this->matrix_info.get_dof_info(no).dofs_per_cell[this->active_fe_index] + for (unsigned int no=0; nomatrix_info->n_components(); ++no) + if (this->matrix_info->get_dof_info(no).dofs_per_cell[this->active_fe_index] == dofs_per_cell) { proposed_dof_comp = no; break; } if (n_q_points == - this->mapping_info.mapping_data_gen[quad_no].n_q_points[this->active_quad_index]) + this->mapping_info->mapping_data_gen[quad_no].n_q_points[this->active_quad_index]) proposed_quad_comp = quad_no; else - for (unsigned int no=0; nomapping_info.mapping_data_gen.size(); ++no) - if (this->mapping_info.mapping_data_gen[no].n_q_points[this->active_quad_index] + for (unsigned int no=0; nomapping_info->mapping_data_gen.size(); ++no) + if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index] == n_q_points) { proposed_quad_comp = no; @@ -3777,8 +4149,8 @@ FEEvaluationGeneral } // ok, did not find the numbers specified by the template arguments in // the given list. Suggest correct template arguments - const unsigned int proposed_fe_degree = static_cast(std::pow(1.001*this->data.dofs_per_cell,1./dim))-1; - const unsigned int proposed_n_q_points_1d = static_cast(std::pow(1.001*this->data.n_q_points,1./dim)); + const unsigned int proposed_fe_degree = static_cast(std::pow(1.001*this->data->dofs_per_cell,1./dim))-1; + const unsigned int proposed_n_q_points_1d = static_cast(std::pow(1.001*this->data->n_q_points,1./dim)); message += "Wrong template arguments:\n"; message += " Did you mean FEEvaluation correct_pos += " \n"; message += " " + correct_pos; - Assert (dofs_per_cell == this->data.dofs_per_cell && - n_q_points == this->data.n_q_points, + Assert (dofs_per_cell == this->data->dofs_per_cell && + n_q_points == this->data->n_q_points, ExcMessage(message)); } AssertDimension (n_q_points, - this->mapping_info.mapping_data_gen[this->quad_no]. + this->mapping_info->mapping_data_gen[this->quad_no]. n_q_points[this->active_quad_index]); - AssertDimension (dofs_per_cell * this->n_fe_components, - this->dof_info.dofs_per_cell[this->active_fe_index]); + AssertDimension (this->data->dofs_per_cell * this->n_fe_components, + this->dof_info->dofs_per_cell[this->active_fe_index]); #endif } +template +inline +FEEvaluationGeneral +::FEEvaluationGeneral (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element) + : + BaseClass (geometry, dof_handler, base_element) +{ + // set the pointers to the correct position in the data array + for (unsigned int c=0; cvalues_dofs[c] = &my_data_array[c*dofs_per_cell]; + this->values_quad[c] = &my_data_array[n_components*dofs_per_cell+c*n_q_points]; + for (unsigned int d=0; dgradients_quad[c][d] = &my_data_array[n_components*(dofs_per_cell+n_q_points) + + + (c*dim+d)*n_q_points]; + for (unsigned int d=0; d<(dim*dim+dim)/2; ++d) + this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+dofs_per_cell) + + + (c*(dim*dim+dim)+d)*n_q_points]; + } +} + + + namespace internal { // evaluates the given shape data in 1d-3d using the tensor product @@ -4886,40 +5286,29 @@ namespace internal inline void do_evaluate (FEEval &fe_eval, + VectorizedArray* values_dofs[], + VectorizedArray* values_quad[], + VectorizedArray* gradients_quad[][1], + VectorizedArray* hessians_quad[][1], const bool evaluate_val, const bool evaluate_grad, const bool evaluate_lapl, internal::int2type<1>) { - AssertDimension(FEEval::dimension, 1); - Assert (fe_eval.cell != numbers::invalid_unsigned_int, - ExcNotInitialized()); - Assert (fe_eval.dof_values_initialized == true, - internal::ExcAccessToUninitializedField()); - const unsigned int n_components = FEEval::n_components; for (unsigned int c=0; c - (fe_eval.values_dofs[c], fe_eval.values_quad[c]); + (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]); + (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]); + (values_dofs[c], fe_eval.hessians_quad[c][0]); } - -#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 } @@ -4927,17 +5316,15 @@ namespace internal inline void do_evaluate (FEEval &fe_eval, + VectorizedArray* values_dofs[], + VectorizedArray* values_quad[], + VectorizedArray* gradients_quad[][2], + VectorizedArray* hessians_quad[][3], const bool evaluate_val, const bool evaluate_grad, const bool evaluate_lapl, internal::int2type<2>) { - AssertDimension(FEEval::dimension, 2); - Assert (fe_eval.cell != numbers::invalid_unsigned_int, - ExcNotInitialized()); - Assert (fe_eval.dof_values_initialized == true, - internal::ExcAccessToUninitializedField()); - const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ? FEEval::dofs_per_cell : FEEval::n_q_points; const unsigned int n_components = FEEval::n_components; @@ -4951,69 +5338,58 @@ namespace internal if (evaluate_grad == true) { fe_eval.template apply_gradients<0,true,false> - (fe_eval.values_dofs[c], temp1); + (values_dofs[c], temp1); fe_eval.template apply_values<1,true,false> - (temp1, fe_eval.gradients_quad[c][0]); + (temp1, 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); + (values_dofs[c], temp1); fe_eval.template apply_gradients<1,true,false> - (temp1, fe_eval.hessians_quad[c][2]); + (temp1, hessians_quad[c][2]); // grad xx fe_eval.template apply_hessians<0,true,false> - (fe_eval.values_dofs[c], temp1); + (values_dofs[c], temp1); fe_eval.template apply_values<1,true,false> - (temp1, fe_eval.hessians_quad[c][0]); + (temp1, hessians_quad[c][0]); } // grad y fe_eval.template apply_values<0,true,false> - (fe_eval.values_dofs[c], temp1); + (values_dofs[c], temp1); if (evaluate_grad == true) fe_eval.template apply_gradients<1,true,false> - (temp1, fe_eval.gradients_quad[c][1]); + (temp1, 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]); + (temp1, 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]); + (temp1, values_quad[c]); } - -#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_evaluate (FEEval &fe_eval, + VectorizedArray* values_dofs[], + VectorizedArray* values_quad[], + VectorizedArray* gradients_quad[][3], + VectorizedArray* hessians_quad[][6], const bool evaluate_val, const bool evaluate_grad, const bool evaluate_lapl, internal::int2type<3>) { - AssertDimension(FEEval::dimension, 3); - Assert (fe_eval.cell != numbers::invalid_unsigned_int, - ExcNotInitialized()); - Assert (fe_eval.dof_values_initialized == true, - internal::ExcAccessToUninitializedField()); - const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ? FEEval::dofs_per_cell : FEEval::n_q_points; const unsigned int n_components = FEEval::n_components; @@ -5027,11 +5403,11 @@ namespace internal { // grad x fe_eval.template apply_gradients<0,true,false> - (fe_eval.values_dofs[c], temp1); + (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]); + (temp2, gradients_quad[c][0]); } if (evaluate_lapl == true) @@ -5040,37 +5416,37 @@ namespace internal if (evaluate_grad == false) { fe_eval.template apply_gradients<0,true,false> - (fe_eval.values_dofs[c], temp1); + (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]); + (temp2, 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]); + (temp2, hessians_quad[c][3]); // grad xx fe_eval.template apply_hessians<0,true,false> - (fe_eval.values_dofs[c], temp1); + (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]); + (temp2, hessians_quad[c][0]); } // grad y fe_eval.template apply_values<0,true,false> - (fe_eval.values_dofs[c], temp1); + (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]); + (temp2, gradients_quad[c][1]); } if (evaluate_lapl == true) @@ -5080,13 +5456,13 @@ namespace internal 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]); + (temp2, 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]); + (temp2, hessians_quad[c][1]); } // grad z: can use the values applied in x direction stored in temp1 @@ -5094,28 +5470,19 @@ namespace internal (temp1, temp2); if (evaluate_grad == true) fe_eval.template apply_gradients<2,true,false> - (temp2, fe_eval.gradients_quad[c][2]); + (temp2, gradients_quad[c][2]); // 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]); + (temp2, hessians_quad[c][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]); + (temp2, values_quad[c]); } - -#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 } @@ -5124,57 +5491,43 @@ namespace internal inline void do_integrate (FEEval &fe_eval, + VectorizedArray* values_dofs[], + VectorizedArray* values_quad[], + VectorizedArray* gradients_quad[][1], const bool integrate_val, const bool integrate_grad, internal::int2type<1>) { - 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 n_components = FEEval::n_components; for (unsigned int c=0; c - (fe_eval.gradients_quad[c][0], fe_eval.values_dofs[c]); + (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]); + (values_quad[c], values_dofs[c]); else fe_eval.template apply_values<0,false,false> - (fe_eval.values_quad[c], fe_eval.values_dofs[c]); + (values_quad[c], values_dofs[c]); } } - -#ifdef DEBUG - fe_eval.dof_values_initialized = true; -#endif } template inline void do_integrate (FEEval &fe_eval, + VectorizedArray* values_dofs[], + VectorizedArray* values_quad[], + VectorizedArray* gradients_quad[][2], const bool integrate_val, const bool integrate_grad, internal::int2type<2>) { - 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 = FEEval::dofs_per_cell > FEEval::n_q_points ? FEEval::dofs_per_cell : FEEval::n_q_points; const unsigned int n_components = FEEval::n_components; @@ -5187,51 +5540,42 @@ namespace internal // val if (integrate_val == true) fe_eval.template apply_values<0,false,false> - (fe_eval.values_quad[c], temp1); + (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); + (gradients_quad[c][0], temp1); else fe_eval.template apply_gradients<0,false,false> - (fe_eval.gradients_quad[c][0], temp1); + (gradients_quad[c][0], temp1); } if (integrate_val || integrate_grad) fe_eval.template apply_values<1,false,false> - (temp1, fe_eval.values_dofs[c]); + (temp1, values_dofs[c]); if (integrate_grad == true) { // grad y fe_eval.template apply_values<0,false,false> - (fe_eval.gradients_quad[c][1], temp1); + (gradients_quad[c][1], temp1); fe_eval.template apply_gradients<1,false,true> - (temp1, fe_eval.values_dofs[c]); + (temp1, values_dofs[c]); } } - -#ifdef DEBUG - fe_eval.dof_values_initialized = true; -#endif } template inline void do_integrate (FEEval &fe_eval, + VectorizedArray* values_dofs[], + VectorizedArray* values_quad[], + VectorizedArray* gradients_quad[][3], const bool integrate_val, const bool integrate_grad, internal::int2type<3>) { - 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 = FEEval::dofs_per_cell > FEEval::n_q_points ? FEEval::dofs_per_cell : FEEval::n_q_points; const unsigned int n_components = FEEval::n_components; @@ -5245,17 +5589,17 @@ namespace internal { // val fe_eval.template apply_values<0,false,false> - (fe_eval.values_quad[c], temp1); + (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); + (gradients_quad[c][0], temp1); else fe_eval.template apply_gradients<0,false,false> - (fe_eval.gradients_quad[c][0], temp1); + (gradients_quad[c][0], temp1); } if (integrate_val || integrate_grad) fe_eval.template apply_values<1,false,false> @@ -5264,28 +5608,24 @@ namespace internal { // 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); + (gradients_quad[c][1], temp1); fe_eval.template apply_gradients<1,false,true> (temp1, temp2); } if (integrate_val || integrate_grad) fe_eval.template apply_values<2,false,false> - (temp2, fe_eval.values_dofs[c]); + (temp2, 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); + (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]); + (temp2, values_dofs[c]); } } - -#ifdef DEBUG - fe_eval.dof_values_initialized = true; -#endif } } // end of namespace internal @@ -5301,8 +5641,21 @@ FEEvaluationGeneral const bool evaluate_grad, const bool evaluate_lapl) { - internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl, + Assert (this->dof_values_initialized == true, + internal::ExcAccessToUninitializedField()); + internal::do_evaluate (*this, this->values_dofs, this->values_quad, + this->gradients_quad, this->hessians_quad, + evaluate_val, evaluate_grad, evaluate_lapl, internal::int2type()); + +#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 } @@ -5315,8 +5668,20 @@ FEEvaluationGeneral ::integrate (const bool integrate_val, const bool integrate_grad) { - internal::do_integrate (*this, integrate_val, integrate_grad, + if (integrate_val == true) + Assert (this->values_quad_submitted == true, + internal::ExcAccessToUninitializedField()); + if (integrate_grad == true) + Assert (this->gradients_quad_submitted == true, + internal::ExcAccessToUninitializedField()); + + internal::do_integrate (*this, this->values_dofs, this->values_quad, + this->gradients_quad, integrate_val, integrate_grad, internal::int2type()); + +#ifdef DEBUG + this->dof_values_initialized = true; +#endif } @@ -5328,7 +5693,7 @@ Point > FEEvaluationGeneral ::quadrature_point (const unsigned int q) const { - Assert (this->mapping_info.quadrature_points_initialized == true, + Assert (this->mapping_info->quadrature_points_initialized == true, ExcNotInitialized()); AssertIndexRange (q, n_q_points); @@ -5374,7 +5739,7 @@ FEEvaluationGeneral { internal::apply_tensor_product, direction, dof_to_quad, add> - (this->data.shape_values.begin(), in, out); + (this->data->shape_values.begin(), in, out); } @@ -5390,7 +5755,7 @@ FEEvaluationGeneral { internal::apply_tensor_product, direction, dof_to_quad, add> - (this->data.shape_gradients.begin(), in, out); + (this->data->shape_gradients.begin(), in, out); } @@ -5406,7 +5771,7 @@ FEEvaluationGeneral { internal::apply_tensor_product, direction, dof_to_quad, add> - (this->data.shape_hessians.begin(), in, out); + (this->data->shape_hessians.begin(), in, out); } @@ -5422,6 +5787,33 @@ FEEvaluation const unsigned int quad_no) : BaseClass (data_in, fe_no, quad_no) +{ + compute_even_odd_factors(); +} + + + +template +inline +FEEvaluation +::FEEvaluation (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element) + : + BaseClass (geometry, dof_handler, base_element) +{ + compute_even_odd_factors(); +} + + + +template +inline +void +FEEvaluation +::compute_even_odd_factors() { // check whether element is appropriate #ifdef DEBUG @@ -5436,8 +5828,8 @@ FEEvaluation const unsigned int n_dofs_1d = fe_degree + 1; for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) for (unsigned int j=0; jdata.shape_values[i*n_q_points_1d+j][0] - - this->data.shape_values[(n_dofs_1d-i)*n_q_points_1d + Assert (std::fabs(this->data->shape_values[i*n_q_points_1d+j][0] - + this->data->shape_values[(n_dofs_1d-i)*n_q_points_1d -j-1][0]) < zero_tol, ExcMessage(error_message)); @@ -5446,10 +5838,10 @@ FEEvaluation if (n_q_points_1d%2 == 1 && n_dofs_1d%2 == 1) { for (int i=0; i(n_dofs_1d/2); ++i) - Assert (std::fabs(this->data.shape_values[i*n_q_points_1d+ + Assert (std::fabs(this->data->shape_values[i*n_q_points_1d+ n_q_points_1d/2][0]) < zero_tol, ExcMessage(error_message)); - Assert (std::fabs(this->data.shape_values[(n_dofs_1d/2)*n_q_points_1d+ + Assert (std::fabs(this->data->shape_values[(n_dofs_1d/2)*n_q_points_1d+ n_q_points_1d/2][0]-1.)< zero_tol, ExcMessage(error_message)); } @@ -5458,12 +5850,12 @@ FEEvaluation // quadrature point for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) for (unsigned int j=0; jdata.shape_gradients[i*n_q_points_1d+j][0] + - this->data.shape_gradients[(n_dofs_1d-i)*n_q_points_1d- + Assert (std::fabs(this->data->shape_gradients[i*n_q_points_1d+j][0] + + this->data->shape_gradients[(n_dofs_1d-i)*n_q_points_1d- j-1][0]) < zero_tol, ExcMessage(error_message)); if (n_dofs_1d%2 == 1 && n_q_points_1d%2 == 1) - Assert (std::fabs(this->data.shape_gradients[(n_dofs_1d/2)*n_q_points_1d+ + Assert (std::fabs(this->data->shape_gradients[(n_dofs_1d/2)*n_q_points_1d+ (n_q_points_1d/2)][0]) < zero_tol, ExcMessage(error_message)); @@ -5471,8 +5863,8 @@ FEEvaluation // symmetry for Laplacian for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) for (unsigned int j=0; jdata.shape_hessians[i*n_q_points_1d+j][0] - - this->data.shape_hessians[(n_dofs_1d-i)*n_q_points_1d- + Assert (std::fabs(this->data->shape_hessians[i*n_q_points_1d+j][0] - + this->data->shape_hessians[(n_dofs_1d-i)*n_q_points_1d- j-1][0]) < zero_tol, ExcMessage(error_message)); #endif @@ -5483,35 +5875,35 @@ FEEvaluation for (unsigned int q=0; q<(n_q_points_1d+1)/2; ++q) { shape_val_evenodd[i][q] = - 0.5 * (this->data.shape_values[i*n_q_points_1d+q] + - this->data.shape_values[i*n_q_points_1d+n_q_points_1d-1-q]); + 0.5 * (this->data->shape_values[i*n_q_points_1d+q] + + this->data->shape_values[i*n_q_points_1d+n_q_points_1d-1-q]); shape_val_evenodd[fe_degree-i][q] = - 0.5 * (this->data.shape_values[i*n_q_points_1d+q] - - this->data.shape_values[i*n_q_points_1d+n_q_points_1d-1-q]); + 0.5 * (this->data->shape_values[i*n_q_points_1d+q] - + this->data->shape_values[i*n_q_points_1d+n_q_points_1d-1-q]); shape_gra_evenodd[i][q] = - 0.5 * (this->data.shape_gradients[i*n_q_points_1d+q] + - this->data.shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]); + 0.5 * (this->data->shape_gradients[i*n_q_points_1d+q] + + this->data->shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]); shape_gra_evenodd[fe_degree-i][q] = - 0.5 * (this->data.shape_gradients[i*n_q_points_1d+q] - - this->data.shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]); + 0.5 * (this->data->shape_gradients[i*n_q_points_1d+q] - + this->data->shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]); shape_hes_evenodd[i][q] = - 0.5 * (this->data.shape_hessians[i*n_q_points_1d+q] + - this->data.shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]); + 0.5 * (this->data->shape_hessians[i*n_q_points_1d+q] + + this->data->shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]); shape_hes_evenodd[fe_degree-i][q] = - 0.5 * (this->data.shape_hessians[i*n_q_points_1d+q] - - this->data.shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]); + 0.5 * (this->data->shape_hessians[i*n_q_points_1d+q] - + this->data->shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]); } if (fe_degree % 2 == 0) for (unsigned int q=0; q<(n_q_points_1d+1)/2; ++q) { shape_val_evenodd[fe_degree/2][q] = - this->data.shape_values[(fe_degree/2)*n_q_points_1d+q]; + this->data->shape_values[(fe_degree/2)*n_q_points_1d+q]; shape_gra_evenodd[fe_degree/2][q] = - this->data.shape_gradients[(fe_degree/2)*n_q_points_1d+q]; + this->data->shape_gradients[(fe_degree/2)*n_q_points_1d+q]; shape_hes_evenodd[fe_degree/2][q] = - this->data.shape_hessians[(fe_degree/2)*n_q_points_1d+q]; + this->data->shape_hessians[(fe_degree/2)*n_q_points_1d+q]; } } @@ -5526,8 +5918,21 @@ FEEvaluation const bool evaluate_grad, const bool evaluate_lapl) { - internal::do_evaluate (*this, evaluate_val, evaluate_grad, evaluate_lapl, + Assert (this->dof_values_initialized == true, + internal::ExcAccessToUninitializedField()); + internal::do_evaluate (*this, this->values_dofs, this->values_quad, + this->gradients_quad, this->hessians_quad, + evaluate_val, evaluate_grad, evaluate_lapl, internal::int2type()); + +#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 } @@ -5539,8 +5944,20 @@ void FEEvaluation ::integrate (bool integrate_val,bool integrate_grad) { - internal::do_integrate (*this, integrate_val, integrate_grad, + if (integrate_val == true) + Assert (this->values_quad_submitted == true, + internal::ExcAccessToUninitializedField()); + if (integrate_grad == true) + Assert (this->gradients_quad_submitted == true, + internal::ExcAccessToUninitializedField()); + + internal::do_integrate (*this, this->values_dofs, this->values_quad, + this->gradients_quad, integrate_val, integrate_grad, internal::int2type()); + +#ifdef DEBUG + this->dof_values_initialized = true; +#endif } @@ -5564,7 +5981,7 @@ FEEvaluation else internal::apply_tensor_product_values, direction, dof_to_quad, add> - (this->data.shape_values.begin(), in, out); + (this->data->shape_values.begin(), in, out); } @@ -5585,7 +6002,7 @@ FEEvaluation else internal::apply_tensor_product_gradients, direction, dof_to_quad, add> - (this->data.shape_gradients.begin(), in, out); + (this->data->shape_gradients.begin(), in, out); } @@ -5609,7 +6026,7 @@ FEEvaluation else internal::apply_tensor_product_hessians, direction, dof_to_quad, add> - (this->data.shape_hessians.begin(), in, out); + (this->data->shape_hessians.begin(), in, out); } @@ -5640,19 +6057,19 @@ FEEvaluationGL for (unsigned int j=0; jdata.shape_values[i*n_points_1d+j][0])data->shape_values[i*n_points_1d+j][0])data.shape_values[i*n_points_1d+ + Assert (std::fabs(this->data->shape_values[i*n_points_1d+ j][0]-1.)data.shape_gradients[i*n_points_1d+i][0])data->shape_gradients[i*n_points_1d+i][0])data.shape_gradients[n_points_1d-1][0]- + Assert (std::fabs(this->data->shape_gradients[n_points_1d-1][0]- (n_points_1d%2==0 ? -1. : 1.)) < zero_tol, ExcMessage (error_mess.c_str())); #endif @@ -5660,6 +6077,18 @@ FEEvaluationGL +template +inline +FEEvaluationGL +::FEEvaluationGL (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element) + : + BaseClass (geometry, dof_handler, base_element) +{} + + + template inline void @@ -5852,10 +6281,142 @@ FEEvaluationGL { internal::apply_tensor_product_gradients_gl, direction, dof_to_quad, add> - (this->data.shape_gradients.begin(), in, out); + (this->data->shape_gradients.begin(), in, out); } + +/*------------------------- FEEvaluationDGP ---------------------------------*/ + +template +inline +FEEvaluationDGP +::FEEvaluationDGP (const MatrixFree &data_in, + const unsigned int fe_no, + const unsigned int quad_no) + : + BaseClass (data_in, fe_no, quad_no) +{} + + + +template +inline +FEEvaluationDGP +::FEEvaluationDGP (const EvaluatedGeometry &geometry, + const DoFHandler &dof_handler, + const unsigned int base_element) + : + BaseClass (geometry, dof_handler, base_element) +{} + + +template +inline +void +FEEvaluationDGP +::evaluate (const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl) +{ + Assert (this->dof_values_initialized == true, + internal::ExcAccessToUninitializedField()); + + // expand dof_values to tensor product + VectorizedArray data_array[n_components*Utilities::fixed_int_power::value]; + VectorizedArray *expanded_dof_values[n_components]; + for (unsigned int c=0; c::value]; + + unsigned int count_p = 0, count_q = 0; + for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i) + { + for (int j=0; j<(dim>1?fe_degree+1-i:1); ++j) + { + for (int k=0; kvalues_dofs[c][count_p]; + } + for (int k=fe_degree+1-j-i; k(); + } + for (unsigned int j=(dim>1?fe_degree+1-i:1); j<(dim>1?fe_degree+1:1); ++j) + for (unsigned int k=0; k(); + } + AssertDimension(count_q, BaseClass::dofs_per_cell); + AssertDimension(count_p, dofs_per_cell); + + internal::do_evaluate (*this, expanded_dof_values, this->values_quad, + this->gradients_quad, this->hessians_quad, + evaluate_val, evaluate_grad, evaluate_lapl, + internal::int2type()); + +#ifdef DEBUG + if (evaluate_val == true) + this->values_quad_initialized = true; + if (evaluate_grad == true) + this->gradients_quad_initialized = true; + if (evaluate_lapl == true) + this->hessians_quad_initialized = true; +#endif +} + + + +template +inline +void +FEEvaluationDGP +::integrate (bool integrate_val,bool integrate_grad) +{ + if (integrate_val == true) + Assert (this->values_quad_submitted == true, + internal::ExcAccessToUninitializedField()); + if (integrate_grad == true) + Assert (this->gradients_quad_submitted == true, + internal::ExcAccessToUninitializedField()); + + VectorizedArray data_array[n_components*Utilities::fixed_int_power::value]; + VectorizedArray *expanded_dof_values[n_components]; + for (unsigned int c=0; c::value]; + internal::do_integrate (*this, expanded_dof_values, this->values_quad, + this->gradients_quad, integrate_val, integrate_grad, + internal::int2type()); + + // truncate tensor product + unsigned int count_p = 0, count_q = 0; + for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i) + { + for (int j=0; j<(dim>1?fe_degree+1-i:1); ++j) + { + for (int k=0; kvalues_dofs[c][count_p] = expanded_dof_values[c][count_q]; + } + count_q += j+i; + } + count_q += i*(fe_degree+1); + } + AssertDimension(count_q, BaseClass::dofs_per_cell); + AssertDimension(count_p, dofs_per_cell); + +#ifdef DEBUG + this->dof_values_initialized = true; +#endif +} + + + #endif // ifndef DOXYGEN 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 ee1f8c1961..d2d616bb75 100644 --- a/deal.II/include/deal.II/matrix_free/mapping_info.h +++ b/deal.II/include/deal.II/matrix_free/mapping_info.h @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2011 - 2013 by the deal.II authors +// Copyright (C) 2011 - 2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -86,9 +86,10 @@ namespace internal * Helper function to determine which update flags must be set in the * internal functions to initialize all data as requested by the user. */ - UpdateFlags + static UpdateFlags compute_update_flags (const UpdateFlags update_flags, - const std::vector > &quad) const; + const std::vector > &quad = + std::vector >()); /** * Returns the type of a given cell as detected during initialization. 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 5d093bc971..74051a1c85 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 @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2011 - 2013 by the deal.II authors +// Copyright (C) 2011 - 2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -62,7 +62,7 @@ namespace internal UpdateFlags MappingInfo:: compute_update_flags (const UpdateFlags update_flags, - const std::vector > &quad) const + const std::vector > &quad) { // this class is build around the evaluation this class is build around // the evaluation of inverse gradients, so compute them in any case @@ -89,22 +89,23 @@ namespace internal // one quadrature point on the first component, but more points on later // components, we need to have Jacobian gradients anyway in order to // determine whether the Jacobian is constant throughout a cell - bool formula_with_one_point = false; - for (unsigned int i=0; i 1) + if (quad.empty() == false) + { + bool formula_with_one_point = false; + for (unsigned int i=0; i 1) + { + new_flags |= update_jacobian_grads; + } + } return new_flags; } 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 dfdc478210..4b84ba49d0 100644 --- a/deal.II/include/deal.II/matrix_free/matrix_free.h +++ b/deal.II/include/deal.II/matrix_free/matrix_free.h @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2011 - 2013 by the deal.II authors +// Copyright (C) 2011 - 2014 by the deal.II authors // // This file is part of the deal.II library. // 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 25dab92538..b50d23bc01 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 @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2011 - 2013 by the deal.II authors +// Copyright (C) 2011 - 2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -526,38 +526,22 @@ void MatrixFree::initialize_indices // there the vertex DoFs come first, which is incompatible with the // lexicographic ordering necessary to apply tensor products // efficiently) - const FE_Poly,dim,dim> *fe_poly = - dynamic_cast,dim,dim>*> - (&fe.base_element(0)); - const FE_Poly >,dim,dim> *fe_poly_piece = - dynamic_cast >,dim,dim>*> - (&fe.base_element(0)); - - // This class currently only works for elements derived from - // FE_Poly,dim,dim> or piecewise - // polynomials. For any other element, the dynamic casts above will - // fail and give fe_poly == 0. - Assert (fe_poly != 0 || fe_poly_piece != 0, ExcNotImplemented()); + if (n_fe_components == 1) { - if (fe_poly != 0) - lexicographic_inv[no][fe_index] = - fe_poly->get_poly_space_numbering_inverse(); - else - lexicographic_inv[no][fe_index] = - fe_poly_piece->get_poly_space_numbering_inverse(); + lexicographic_inv[no][fe_index] = + internal::MatrixFreeFunctions::get_lexicographic_numbering_inverse(fe); AssertDimension (lexicographic_inv[no][fe_index].size(), dof_info[no].dofs_per_cell[fe_index]); } else { - // ok, we have more than one component + // ok, we have more than one component, get the inverse + // permutation, invert it, sort the components one after one, + // and invert back Assert (n_fe_components > 1, ExcInternalError()); std::vector scalar_lex = - fe_poly != 0 ? fe_poly->get_poly_space_numbering() : - fe_poly_piece->get_poly_space_numbering(); + Utilities::invert_permutation(internal::MatrixFreeFunctions::get_lexicographic_numbering_inverse(fe.base_element(0))); AssertDimension (scalar_lex.size() * n_fe_components, dof_info[no].dofs_per_cell[fe_index]); std::vector lexicographic (dof_info[no].dofs_per_cell[fe_index]); diff --git a/deal.II/include/deal.II/matrix_free/shape_info.h b/deal.II/include/deal.II/matrix_free/shape_info.h index 8bbae186b8..24c5e74ae8 100644 --- a/deal.II/include/deal.II/matrix_free/shape_info.h +++ b/deal.II/include/deal.II/matrix_free/shape_info.h @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2011 - 2013 by the deal.II authors +// Copyright (C) 2011 - 2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -49,6 +49,14 @@ namespace internal */ ShapeInfo (); + /** + * Constructor that initializes the data fields using the reinit method. + */ + template + ShapeInfo (const Quadrature<1> &quad, + const FiniteElement &fe, + const unsigned int base_element = 0); + /** * Initializes the data fields. Takes a one-dimensional quadrature * formula and a finite element as arguments and evaluates the shape @@ -59,7 +67,8 @@ namespace internal */ template void reinit (const Quadrature<1> &quad, - const FiniteElement &fe_dim); + const FiniteElement &fe_dim, + const unsigned int base_element = 0); /** * Returns the memory consumption of this @@ -170,6 +179,35 @@ namespace internal unsigned int dofs_per_face; }; + + /** + * Extracts the mapping between the actual numbering inside a scalar + * element and lexicographic numbering as required by the evaluation + * routines. + */ + template + std::vector + get_lexicographic_numbering_inverse (const FiniteElement &fe); + + + + // ------------------------------------------ inline functions + + template + template + inline + ShapeInfo::ShapeInfo (const Quadrature<1> &quad, + const FiniteElement &fe_in, + const unsigned int base_element_number) + : + n_q_points (0), + dofs_per_cell (0) + { + reinit (quad, fe_in, base_element_number); + } + + + } // end of namespace MatrixFreeFunctions } // end of namespace internal diff --git a/deal.II/include/deal.II/matrix_free/shape_info.templates.h b/deal.II/include/deal.II/matrix_free/shape_info.templates.h index 852a6369a0..d11c8654e0 100644 --- a/deal.II/include/deal.II/matrix_free/shape_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/shape_info.templates.h @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2011 - 2013 by the deal.II authors +// Copyright (C) 2011 - 2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -21,6 +21,7 @@ #include #include #include +#include #include @@ -33,6 +34,43 @@ namespace internal namespace MatrixFreeFunctions { + // helper function + template + std::vector + get_lexicographic_numbering_inverse(const FiniteElement &fe) + { + Assert(fe.n_components() == 1, + ExcMessage("Expected a scalar element")); + + const FE_Poly,dim,dim> *fe_poly = + dynamic_cast,dim,dim>*>(&fe); + + const FE_Poly >,dim,dim> *fe_poly_piece = + dynamic_cast >,dim,dim>*> (&fe); + + const FE_DGP *fe_dgp = dynamic_cast*>(&fe); + + std::vector lexicographic; + if (fe_poly != 0) + lexicographic = fe_poly->get_poly_space_numbering_inverse(); + else if (fe_poly_piece != 0) + lexicographic = fe_poly_piece->get_poly_space_numbering_inverse(); + else if (fe_dgp != 0) + { + lexicographic.resize(fe_dgp->dofs_per_cell); + for (unsigned int i=0; idofs_per_cell; ++i) + lexicographic[i] = i; + } + else + Assert(false, ExcNotImplemented()); + + return lexicographic; + } + + + // ----------------- actual ShapeInfo functions -------------------- template @@ -48,43 +86,39 @@ namespace internal template void ShapeInfo::reinit (const Quadrature<1> &quad, - const FiniteElement &fe) + const FiniteElement &fe_in, + const unsigned int base_element_number) { - Assert (fe.n_components() == 1, + const FiniteElement *fe = &fe_in; + if (fe_in.n_components() > 1) + fe = &fe_in.base_element(base_element_number); + + Assert (fe->n_components() == 1, ExcMessage("FEEvaluation only works for scalar finite elements.")); - const unsigned int n_dofs_1d = fe.degree+1, + + const unsigned int n_dofs_1d = fe->degree+1, n_q_points_1d = quad.size(); - AssertDimension(fe.dofs_per_cell, Utilities::fixed_power(n_dofs_1d)); - std::vector lexicographic (fe.dofs_per_cell); + std::vector lexicographic (fe->dofs_per_cell); // renumber (this is necessary for FE_Q, for example, since there the // vertex DoFs come first, which is incompatible with the lexicographic // ordering necessary to apply tensor products efficiently) { - const FE_Poly,dim,dim> *fe_poly = - dynamic_cast,dim,dim>*>(&fe); - const FE_Poly >,dim,dim> *fe_poly_piece = - dynamic_cast >,dim,dim>*> (&fe); - Assert (fe_poly != 0 || fe_poly_piece, ExcNotImplemented()); - lexicographic = fe_poly != 0 ? - fe_poly->get_poly_space_numbering_inverse() : - fe_poly_piece->get_poly_space_numbering_inverse(); + lexicographic = get_lexicographic_numbering_inverse(*fe); // to evaluate 1D polynomials, evaluate along the line where y=z=0, // assuming that shape_value(0,Point()) == 1. otherwise, need // other entry point (e.g. generating a 1D element by reading the // name, as done before r29356) - Assert(std::fabs(fe.shape_value(lexicographic[0], Point())-1) < 1e-13, + Assert(std::fabs(fe->shape_value(lexicographic[0], Point())-1) < 1e-13, ExcInternalError()); } n_q_points = Utilities::fixed_power(n_q_points_1d); - dofs_per_cell = Utilities::fixed_power(n_dofs_1d); + dofs_per_cell = fe->dofs_per_cell; n_q_points_face = dim>1?Utilities::fixed_power(n_q_points_1d):1; - dofs_per_face = dim>1?Utilities::fixed_power(n_dofs_1d):1; + dofs_per_face = fe->dofs_per_face; const unsigned int array_size = n_dofs_1d*n_q_points_1d; this->shape_gradients.resize_fast (array_size); @@ -113,25 +147,25 @@ namespace internal // non-vectorized fields Point q_point; q_point[0] = quad.get_points()[q][0]; - shape_values_number[i*n_q_points_1d+q] = fe.shape_value(my_i,q_point); - shape_gradient_number[i*n_q_points_1d+q] = fe.shape_grad (my_i,q_point)[0]; + shape_values_number[i*n_q_points_1d+q] = fe->shape_value(my_i,q_point); + shape_gradient_number[i*n_q_points_1d+q] = fe->shape_grad (my_i,q_point)[0]; shape_values [i*n_q_points_1d+q] = shape_values_number [i*n_q_points_1d+q]; shape_gradients[i*n_q_points_1d+q] = shape_gradient_number[i*n_q_points_1d+q]; shape_hessians[i*n_q_points_1d+q] = - fe.shape_grad_grad(my_i,q_point)[0][0]; + fe->shape_grad_grad(my_i,q_point)[0][0]; q_point[0] *= 0.5; - subface_value[0][i*n_q_points_1d+q] = fe.shape_value(my_i,q_point); + subface_value[0][i*n_q_points_1d+q] = fe->shape_value(my_i,q_point); q_point[0] += 0.5; - subface_value[1][i*n_q_points_1d+q] = fe.shape_value(my_i,q_point); + subface_value[1][i*n_q_points_1d+q] = fe->shape_value(my_i,q_point); } Point q_point; - this->face_value[0][i] = fe.shape_value(my_i,q_point); - this->face_gradient[0][i] = fe.shape_grad(my_i,q_point)[0]; + this->face_value[0][i] = fe->shape_value(my_i,q_point); + this->face_gradient[0][i] = fe->shape_grad(my_i,q_point)[0]; q_point[0] = 1; - this->face_value[1][i] = fe.shape_value(my_i,q_point); - this->face_gradient[1][i] = fe.shape_grad(my_i,q_point)[0]; + this->face_value[1][i] = fe->shape_value(my_i,q_point); + this->face_gradient[1][i] = fe->shape_grad(my_i,q_point)[0]; } // face information @@ -157,7 +191,7 @@ namespace internal } case 2: { - for (unsigned int i=0; idofs_per_face; i++) { this->face_indices(0,i) = n_dofs_1d*i; this->face_indices(1,i) = n_dofs_1d*i + n_dofs_1d-1; diff --git a/tests/matrix_free/matrix_vector_07.cc b/tests/matrix_free/matrix_vector_07.cc index 663b35ff98..3735885557 100644 --- a/tests/matrix_free/matrix_vector_07.cc +++ b/tests/matrix_free/matrix_vector_07.cc @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2013 by the deal.II authors +// Copyright (C) 2013-2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -52,12 +52,14 @@ void test () if (cell->center().norm()<0.2) cell->set_refine_flag(); tria.execute_coarsening_and_refinement(); +#ifndef DEBUG 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(4-dim); +#endif cell = tria.begin_active (); for (unsigned int i=0; i<10-3*dim; ++i) { diff --git a/tests/matrix_free/matrix_vector_14.cc b/tests/matrix_free/matrix_vector_14.cc new file mode 100644 index 0000000000..7f0bd5ac0a --- /dev/null +++ b/tests/matrix_free/matrix_vector_14.cc @@ -0,0 +1,254 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// this test is similar to matrix_vector_07, but implements the operations for +// FE_DGP instead of FE_DGQ (where there is no complete tensor product and +// different routines need to be used). The data is still not very useful +// because the matrix does not include face terms actually present in an +// approximation of the Laplacian. It only contains cell terms. + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +std::ofstream logfile("output"); + + + +template +void +helmholtz_operator_dgp (const MatrixFree &data, + VECTOR &dst, + const VECTOR &src, + const std::pair &cell_range) +{ + typedef typename VECTOR::value_type Number; + FEEvaluationDGP fe_eval (data); + const unsigned int n_q_points = fe_eval.n_q_points; + + for (unsigned int cell=cell_range.first; cell > +class MatrixFreeTest +{ +public: + typedef VectorizedArray vector_t; + + MatrixFreeTest(const MatrixFree &data_in): + data (data_in) + {}; + + void vmult (VECTOR &dst, + const VECTOR &src) const + { + dst = 0; + const std_cxx1x::function &, + VECTOR &, + const VECTOR &, + const std::pair &)> + wrap = helmholtz_operator_dgp; + data.cell_loop (wrap, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void do_test (const DoFHandler &dof, + const ConstraintMatrix &constraints, + const unsigned int parallel_option = 0) +{ + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + if (parallel_option > 0) + deallog << "Parallel option: " << parallel_option << std::endl; + //std::cout << "Number of cells: " << dof.get_tria().n_active_cells() << std::endl; + //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; + //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl; + + MatrixFree mf_data; + { + const QGauss<1> quad (fe_degree+1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + data.tasks_block_size = 3; + + mf_data.reinit (dof, constraints, quad, data); + } + + MatrixFreeTest mf (mf_data); + Vector in (dof.n_dofs()), out (dof.n_dofs()); + Vector in_dist (dof.n_dofs()); + Vector out_dist (in_dist); + + for (unsigned int i=0; i sparse_matrix (sparsity); + { + QGauss quadrature_formula(fe_degree+1); + + FEValues fe_values (dof.get_fe(), quadrature_formula, + update_values | update_gradients | + update_JxW_values); + + const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + fe_values.reinit (cell); + + for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, + local_dof_indices, + sparse_matrix); + } + } + + sparse_matrix.vmult (out, in); + out -= out_dist; + const double diff_norm = out.linfty_norm() / out_dist.linfty_norm(); + + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} + + + +template +void test () +{ + Triangulation tria; + GridGenerator::hyper_ball (tria); + static const HyperBallBoundary boundary; + tria.set_boundary (0, boundary); + if (dim < 3 || fe_degree < 2) + tria.refine_global(1); + tria.begin(tria.n_levels()-1)->set_refine_flag(); + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + typename Triangulation::active_cell_iterator + cell = tria.begin_active (), + endc = tria.end(); + for (; cell!=endc; ++cell) + if (cell->center().norm()<1e-8) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + FE_DGP fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + ConstraintMatrix constraints; + + do_test (dof, constraints); +} + + + +int main () +{ + deallog.attach(logfile); + deallog.depth_console(0); + + deallog << std::setprecision (3); + + { + deallog.threshold_double(5.e-11); + deallog.push("2d"); + test<2,1>(); + test<2,2>(); + deallog.pop(); + deallog.push("3d"); + test<3,1>(); + test<3,2>(); + deallog.pop(); + } +} + + diff --git a/tests/matrix_free/matrix_vector_14.output b/tests/matrix_free/matrix_vector_14.output new file mode 100644 index 0000000000..2a6f032686 --- /dev/null +++ b/tests/matrix_free/matrix_vector_14.output @@ -0,0 +1,13 @@ + +DEAL:2d::Testing FE_DGP<2>(1) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:2d::Testing FE_DGP<2>(2) +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Testing FE_DGP<3>(1) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: +DEAL:3d::Testing FE_DGP<3>(2) +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/matrix_free/thread_correctness.cc b/tests/matrix_free/thread_correctness.cc index edbd65c2ac..35cc44fbce 100644 --- a/tests/matrix_free/thread_correctness.cc +++ b/tests/matrix_free/thread_correctness.cc @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2013 by the deal.II authors +// Copyright (C) 2013-2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -42,6 +42,7 @@ void sub_test() if (cell->center().norm()<0.5) cell->set_refine_flag(); tria.execute_coarsening_and_refinement(); +#ifndef DEBUG if (dim < 3 || fe_degree < 2) tria.refine_global(1); tria.begin(tria.n_levels()-1)->set_refine_flag(); @@ -51,6 +52,7 @@ void sub_test() tria.refine_global(2); else tria.refine_global(1); +#endif FE_Q fe (fe_degree); DoFHandler dof (tria); diff --git a/tests/matrix_free/thread_correctness_hp.cc b/tests/matrix_free/thread_correctness_hp.cc index 1cb0a8cb88..0d9015d9bb 100644 --- a/tests/matrix_free/thread_correctness_hp.cc +++ b/tests/matrix_free/thread_correctness_hp.cc @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2013 by the deal.II authors +// Copyright (C) 2013-2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -98,7 +98,10 @@ void do_test (const unsigned int parallel_option) { Triangulation tria; create_mesh (tria); + // larger mesh in release mode +#ifndef DEBUG tria.refine_global(2); +#endif // refine a few cells for (unsigned int i=0; i<11-3*dim; ++i) @@ -171,12 +174,18 @@ void do_test (const unsigned int parallel_option) MatrixFree::AdditionalData::partition_partition; deallog << "Parallel option partition/partition" << std::endl; } - else + else if (parallel_option == 1) { data.tasks_parallel_scheme = MatrixFree::AdditionalData::partition_color; deallog << "Parallel option partition/color" << std::endl; } + else + { + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::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);