From: Katharina Kormann Date: Wed, 25 Apr 2018 16:47:34 +0000 (+0200) Subject: Introduce the new class FEFaceEvaluation. Speed up vector access for contiguous case. X-Git-Tag: v9.0.0-rc1~93^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a1e0031ec2202cf645386ac6825e883ecd476d7e;p=dealii.git Introduce the new class FEFaceEvaluation. Speed up vector access for contiguous case. --- diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 38ab409803..9cb75a09d4 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -99,43 +99,6 @@ public: */ ~FEEvaluationBase(); - /** - * Initializes the operation pointer to the current cell. Unlike the reinit - * functions taking a cell iterator as argument below and the - * FEValues::reinit() methods, where the information related to a particular - * cell is generated in the reinit call, this function is very cheap since - * all data is pre-computed in @p matrix_free, and only a few indices have - * to be set appropriately. - */ - void reinit (const unsigned int cell); - - /** - * Initialize the data to the current cell using a TriaIterator object as - * usual in FEValues. The argument is either of type - * DoFHandler::active_cell_iterator or DoFHandler::level_cell_iterator. This - * option is only available if the FEEvaluation object was created with a - * finite element, quadrature formula and correct update flags and - * without a MatrixFree object. This initialization method loses the - * ability to use vectorization, see also the description of the - * FEEvaluation class. When this reinit method is used, FEEvaluation can - * also read from vectors (but less efficient than with data coming from - * MatrixFree). - */ - template - void reinit (const TriaIterator > &cell); - - /** - * Initialize the data to the current cell using a TriaIterator object as - * usual in FEValues. This option is only available if the FEEvaluation - * object was created with a finite element, quadrature formula and correct - * update flags and without a MatrixFree object. This initialization - * method loses the ability to use vectorization, see also the description - * of the FEEvaluation class. When this reinit method is used, FEEvaluation - * can not read from vectors because no DoFHandler information is - * available. - */ - void reinit (const typename Triangulation::cell_iterator &cell); - /** * @deprecated Use get_mapping_data_index_offset() instead. */ @@ -763,6 +726,18 @@ protected: VectorType *vectors[], const bool apply_constraints = true) const; + /** + * A unified function to read from and write into vectors based on the given + * template operation for DG-type schemes where all degrees of freedom on + * cells are contiguous. It can perform the operation for read_dof_values(), + * distribute_local_to_global(), and set_dof_values() for for several + * vectors at a time, depending on n_components. + */ + template + void + read_write_operation_contiguous (const VectorOperation &operation, + VectorType *vectors[]) const; + /** * A unified function to read from and write into vectors based on the given * template operation for the case when we do not have an underlying @@ -939,6 +914,12 @@ protected: **/ bool is_interior_face; + /** + * Stores the index an FEFaceEvaluation object is currently pointing into + * (interior face, exterior face, data associated with cell). + */ + unsigned int face_vector_access_index; + /** * Stores the current number of a face within the given cell in case * `is_face==true`, using values between `0` and `2*dim`. @@ -2046,15 +2027,74 @@ template { public: + /** + * A typedef to the base class. + */ typedef FEEvaluationAccess BaseClass; + + /** + * A underlying number type specified as template argument. + */ typedef Number number_type; + + /** + * The type of function values, e.g. `VectorizedArray` for + * `n_components=1` or `Tensor<1,dim,VectorizedArray >` for + * `n_components=dim`. + */ typedef typename BaseClass::value_type value_type; + + /** + * The type of gradients, e.g. `Tensor<1,dim,VectorizedArray>` for + * `n_components=1` or `Tensor<2,dim,VectorizedArray >` for + * `n_components=dim`. + */ typedef typename BaseClass::gradient_type gradient_type; + + /** + * The dimension given as template argument. + */ static constexpr unsigned int dimension = dim; + + /** + * The number of solution components of the evaluator given as template + * argument. + */ static constexpr unsigned int n_components = n_components_; + + /** + * The static number of quadrature points determined from the given template + * argument `n_q_points_1d`. Note that the actual number of quadrature + * points, `n_q_points`, can be different if `fe_degree=-1` is given and + * run-time loop lengths are used rather than compile time ones. + */ static constexpr unsigned int static_n_q_points = Utilities::pow(n_q_points_1d, dim); + + /** + * The static number of degrees of freedom of a scalar component determined + * from the given template argument `fe_degree`. Note that the actual number + * of degrees of freedom `dofs_per_component` can be different if + * `fe_degree=-1` is given or if the underlying is of more complicated type + * than the usual FE_Q or FE_DGQ ones, such as FE_DGP. + */ static constexpr unsigned int static_dofs_per_component = Utilities::pow(fe_degree + 1, dim); + + /** + * The static number of degrees of freedom of all components determined from + * the given template argument `fe_degree`. Note that the actual number of + * degrees of freedom `dofs_per_cell` can be different if `fe_degree=-1` is + * given or if the underlying is of more complicated type than the usual + * FE_Q or FE_DGQ ones, such as FE_DGP. + */ static constexpr unsigned int tensor_dofs_per_cell = static_dofs_per_component *n_components; + + /** + * The static number of degrees of freedom of all components determined from + * the given template argument `fe_degree`. Note that the actual number of + * degrees of freedom `dofs_per_cell` can be different if `fe_degree=-1` is + * given or if the underlying is of more complicated type than the usual + * FE_Q or FE_DGQ ones, such as FE_DGP. + */ static constexpr unsigned int static_dofs_per_cell = static_dofs_per_component *n_components; /** @@ -2160,6 +2200,43 @@ public: */ FEEvaluation &operator= (const FEEvaluation &other); + /** + * Initializes the operation pointer to the current cell batch index. Unlike + * the reinit functions taking a cell iterator as argument below and the + * FEValues::reinit() methods, where the information related to a particular + * cell is generated in the reinit call, this function is very cheap since + * all data is pre-computed in @p matrix_free, and only a few indices have + * to be set appropriately. + */ + void reinit (const unsigned int cell_batch_index); + + /** + * Initialize the data to the current cell using a TriaIterator object as + * usual in FEValues. The argument is either of type + * DoFHandler::active_cell_iterator or DoFHandler::level_cell_iterator. This + * option is only available if the FEEvaluation object was created with a + * finite element, quadrature formula and correct update flags and + * without a MatrixFree object. This initialization method loses the + * ability to use vectorization, see also the description of the + * FEEvaluation class. When this reinit method is used, FEEvaluation can + * also read from vectors (but less efficient than with data coming from + * MatrixFree). + */ + template + void reinit (const TriaIterator > &cell); + + /** + * Initialize the data to the current cell using a TriaIterator object as + * usual in FEValues. This option is only available if the FEEvaluation + * object was created with a finite element, quadrature formula and correct + * update flags and without a MatrixFree object. This initialization + * method loses the ability to use vectorization, see also the description + * of the FEEvaluation class. When this reinit method is used, FEEvaluation + * can not read from vectors because no DoFHandler information is + * available. + */ + void reinit (const typename Triangulation::cell_iterator &cell); + /** * Evaluates the function values, the gradients, and the Hessians of the * polynomial interpolation from the DoF values in the input vector to the @@ -2256,7 +2333,8 @@ public: VectorType &output_vector); /** - * Return the q-th quadrature point stored in MappingInfo. + * Return the q-th quadrature point in real coordinates stored in + * MappingInfo. */ Point > quadrature_point (const unsigned int q_point) const; @@ -2297,120 +2375,427 @@ private: -namespace internal +/** + * The class that provides all functions necessary to evaluate functions at + * quadrature points and face integrations. The design of the class is similar + * to FEEvaluation and most of the interfaces are shared with that class, in + * particular most access functions that come from the common base classes + * FEEvaluationAccess and FEEvaluatioBase. Furthermore, the relation of this + * class to FEEvaluation is similar to the relation between FEValues and + * FEFaceValues. + * + * This class has five template arguments: + * + * @tparam dim Dimension in which this class is to be used + * + * @tparam fe_degree Degree of the tensor product finite element with + * fe_degree+1 degrees of freedom per coordinate + * direction. If set to -1, the degree of the underlying + * element will be used, which acts as a run time constant + * rather than a compile time constant that slows down the + * execution. + * + * @tparam n_q_points_1d Number of points in the quadrature formula in 1D, + * usually chosen as fe_degree+1 + * + * @tparam 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 + * can be applied simultaneously with one call (and often + * more efficiently) + * + * @tparam Number Number format, usually @p double or @p float + * + * @author Katharina Kormann and Martin Kronbichler, 2018 + */ +template +class FEFaceEvaluation : public FEEvaluationAccess { - 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_component - { - // this division is always without remainder - static constexpr unsigned int value = - (DGP_dofs_per_component::value * (degree+dim)) / dim; - }; - - // base specialization: 1d elements have 'degree+1' degrees of freedom - template - struct DGP_dofs_per_component<1,degree> - { - static constexpr unsigned int value = degree+1; - }; - } -} - +public: + /** + * A typedef to the base class. + */ + typedef FEEvaluationAccess BaseClass; -/*----------------------- Inline functions ----------------------------------*/ + /** + * A underlying number type specified as template argument. + */ + typedef Number number_type; -#ifndef DOXYGEN + /** + * The type of function values, e.g. `VectorizedArray` for + * `n_components=1` or `Tensor<1,dim,VectorizedArray >` for + * `n_components=dim`. + */ + typedef typename BaseClass::value_type value_type; + /** + * The type of gradients, e.g. `Tensor<1,dim,VectorizedArray>` for + * `n_components=1` or `Tensor<2,dim,VectorizedArray >` for + * `n_components=dim`. + */ + typedef typename BaseClass::gradient_type gradient_type; + /** + * The dimension given as template argument. + */ + static constexpr unsigned int dimension = dim; -/*----------------------- FEEvaluationBase ----------------------------------*/ + /** + * The number of solution components of the evaluator given as template + * argument. + */ + static constexpr unsigned int n_components = n_components_; -template -inline -FEEvaluationBase -::FEEvaluationBase (const MatrixFree &data_in, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no_in, - const unsigned int fe_degree, - const unsigned int n_q_points, - const bool is_interior_face) - : - scratch_data_array (data_in.acquire_scratch_data()), - quad_no (quad_no_in), - n_fe_components (data_in.get_dof_info(dof_no).start_components.back()), - active_fe_index (fe_degree != numbers::invalid_unsigned_int ? - data_in.get_dof_info(dof_no).fe_index_from_degree - (first_selected_component, fe_degree) - : - 0), - active_quad_index (fe_degree != numbers::invalid_unsigned_int ? - (is_face ? - data_in.get_mapping_info().face_data[quad_no_in]. - quad_index_from_n_q_points(n_q_points) - : - data_in.get_mapping_info().cell_data[quad_no_in]. - quad_index_from_n_q_points(n_q_points)) - : - 0), - n_quadrature_points(fe_degree != numbers::invalid_unsigned_int ? n_q_points : - (is_face ? - data_in.get_shape_info - (dof_no, quad_no_in, active_fe_index, active_quad_index).n_q_points_face - : - data_in.get_shape_info - (dof_no, quad_no_in, active_fe_index, active_quad_index).n_q_points)), - matrix_info (&data_in), - dof_info (&data_in.get_dof_info(dof_no)), - mapping_data (internal::MatrixFreeFunctions::MappingInfoCellsOrFaces::get(data_in.get_mapping_info(), quad_no)), - data (&data_in.get_shape_info - (dof_no, quad_no_in, - dof_info->component_to_base_index[first_selected_component], - active_fe_index, active_quad_index)), - jacobian (nullptr), - J_value (nullptr), - normal_vectors (nullptr), - normal_x_jacobian (nullptr), - quadrature_weights (mapping_data->descriptor[active_quad_index].quadrature_weights.begin()), - cell (numbers::invalid_unsigned_int), - is_interior_face (is_interior_face), - cell_type (internal::MatrixFreeFunctions::general), - dof_values_initialized (false), - values_quad_initialized (false), - gradients_quad_initialized(false), - hessians_quad_initialized (false), - values_quad_submitted (false), - gradients_quad_submitted (false), - first_selected_component (first_selected_component) -{ - set_data_pointers(); - Assert (matrix_info->mapping_initialized() == true, - ExcNotInitialized()); - AssertDimension (matrix_info->get_size_info().vectorization_length, - VectorizedArray::n_array_elements); - AssertDimension ((is_face ? data->n_q_points_face : data->n_q_points), - n_quadrature_points); - AssertDimension (n_quadrature_points, - mapping_data->descriptor[active_quad_index].n_q_points); - Assert(dof_info->start_components.back() == 1 || - (int)n_components_ <= - (int)dof_info->start_components[dof_info->component_to_base_index[first_selected_component]+1] - first_selected_component, - ExcMessage("You tried to construct a vector-valued evaluator with " + - Utilities::to_string(n_components) + " components. However, " - "the current base element has only " + - Utilities::to_string(dof_info->start_components[dof_info->component_to_base_index[first_selected_component]+1] - first_selected_component) - + " components left when starting from local element index " + - Utilities::to_string(first_selected_component-dof_info->start_components[dof_info->component_to_base_index[first_selected_component]]) - + " (global index " + Utilities::to_string(first_selected_component) - + ")")); + /** + * The static number of quadrature points determined from the given template + * argument `n_q_points_1d` taken to the power of dim-1. Note that the actual number of quadrature + * points, `n_q_points`, can be different if `fe_degree=-1` is given and + * run-time loop lengths are used rather than compile time ones. + */ + static constexpr unsigned int static_n_q_points = Utilities::pow(n_q_points_1d, dim-1); - // do not check for correct dimensions of data fields here, should be done - // in derived classes -} + /** + * The static number of quadrature points on a cell with the same quadrature + * formula. Note that this value is only present for simpler comparison with + * the cell quadrature, as the actual number of points is given to a face by + * the `static_n_q_points` variable. + */ + static constexpr unsigned int static_n_q_points_cell = Utilities::pow(n_q_points_1d, dim); + + /** + * The static number of degrees of freedom of a scalar component determined + * from the given template argument `fe_degree`. Note that the actual number + * of degrees of freedom `dofs_per_component` can be different if + * `fe_degree=-1` is given. + */ + static constexpr unsigned int static_dofs_per_component = Utilities::pow(fe_degree + 1, dim); + + /** + * The static number of degrees of freedom of all components determined from + * the given template argument `fe_degree`. Note that the actual number of + * degrees of freedom `dofs_per_cell` can be different if `fe_degree=-1` is + * given. + */ + static constexpr unsigned int tensor_dofs_per_cell = static_dofs_per_component *n_components; + + /** + * The static number of degrees of freedom of all components determined from + * the given template argument `fe_degree`. Note that the actual number of + * degrees of freedom `dofs_per_cell` can be different if `fe_degree=-1` is + * given. + */ + static constexpr unsigned int static_dofs_per_cell = static_dofs_per_component *n_components; + + /** + * 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, the appropriate component + * can be selected by the optional arguments. + * + * @param matrix_free Data object that contains all data + * + * @param is_interior_face This selects which of the two cells of an + * internal face the current evaluator will be based upon. The interior face + * is the main face along which the normal vectors are oriented. The + * exterior face coming from the other side provides the same normal vector + * as the interior side, so if the outer normal vector to that side is + * desired, it must be multiplied by -1. + * + * @param dof_no If matrix_free was set up with multiple DoFHandler + * objects, this parameter selects to which DoFHandler/ConstraintMatrix pair + * the given evaluator should be attached to. + * + * @param quad_no If matrix_free was set up with multiple Quadrature + * objects, this parameter selects the appropriate number of the quadrature + * formula. + * + * @param first_selected_component If the dof_handler selected by dof_no + * uses an FESystem consisting of more than one base element, this parameter + * selects the number of the base element in FESystem. Note that this does + * not directly relate to the component of the respective element due to the + * possibility for a multiplicity in the element. + */ + FEFaceEvaluation (const MatrixFree &matrix_free, + const bool is_interior_face = true, + const unsigned int dof_no = 0, + const unsigned int quad_no = 0, + const unsigned int first_selected_component = 0); + + /** + * Destructor. + */ + ~FEFaceEvaluation(); + + /** + * Initializes the operation pointer to the current face. This method is the + * default choice for face integration as the data stored in MappingInfo is + * stored according to this numbering. Unlike the reinit functions taking a + * cell iterator as argument below and the FEValues::reinit() methods, where + * the information related to a particular cell is generated in the reinit + * call, this function is very cheap since all data is pre-computed in + * @p matrix_free, and only a few indices and pointers have to be set + * appropriately. + */ + void reinit (const unsigned int face_batch_number); + + /** + * As opposed to the reinit() method from the base class, this reinit() + * method initializes for a given number of cells and a face number. This + * method is less efficient than the other reinit() method taking a + * numbering of the faces because it needs to copy the data associated with + * the faces to the cells in this call. + */ + void reinit(const unsigned int cell_batch_number, + const unsigned int face_number); + + /** + * Evaluates the function values, the gradients, and the Laplacians of the + * FE function given at the DoF values stored in the internal data field + * `dof_values` (that is usually filled by the read_dof_values() method) at + * the quadrature points on the unit cell. The function arguments specify + * which parts shall actually be computed. Needs to be called before the + * functions get_value(), get_gradient() or get_normal_derivative() give + * useful information (unless these values have been set manually by + * accessing the internal data pointers). + */ + void evaluate (const bool evaluate_values, + const bool evaluate_gradients); + + /** + * Evaluates the function values, the gradients, and the Laplacians of the + * FE function given at the DoF values in the input array `values_array` at + * the quadrature points on the unit cell. If multiple components are + * involved in the current FEEvaluation object, the sorting in values_array + * is such that all degrees of freedom for the first component come first, + * then all degrees of freedom for the second, and so on. The function + * arguments specify which parts shall actually be computed. Needs to be + * called before the functions get_value(), get_gradient(), or + * get_normal_derivative() give useful information (unless these values have + * been set manually). + */ + void evaluate (const VectorizedArray *values_array, + const bool evaluate_values, + const bool evaluate_gradients); + + /** + * Reads from the input vector and evaluates the function values, the + * gradients, and the Laplacians of the FE function at the quadrature points + * on the unit cell. The function arguments specify which parts shall + * actually be computed. Needs to be called before the functions + * get_value(), get_gradient(), or get_normal_derivative() give useful + * information. + * + * This call is equivalent to calling read_dof_values() followed by + * evaluate(), but might internally use some additional optimizations. + */ + template + void gather_evaluate (const VectorType &input_vector, + const bool evaluate_values, + const bool evaluate_gradients); + + /** + * 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 + * `integrate_val` and `integrate_grad` are used to enable/disable some of + * values or gradients. The result is written into the internal data field + * `dof_values` (that is usually written into the result vector by the + * distribute_local_to_global() or set_dof_values() methods). + */ + void integrate (const bool integrate_values, + const bool integrate_gradients); + + /** + * 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 + * `integrate_val` and `integrate_grad` are used to enable/disable some of + * values or gradients. As opposed to the other integrate() method, this + * call stores the result of the testing in the given array `values_array`. + */ + void integrate (const bool integrate_values, + const bool integrate_gradients, + VectorizedArray *values_array); + + /** + * 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 + * `integrate_val` and `integrate_grad` are used to enable/disable some of + * values or gradients. + * + * This call is equivalent to calling integrate() followed by + * distribute_local_to_global(), but might internally use some additional + * optimizations. + */ + template + void integrate_scatter (const bool integrate_values, + const bool integrate_gradients, + VectorType &output_vector); + + /** + * Returns the q-th quadrature point on the face in real coordinates stored + * in MappingInfo. + */ + Point > + quadrature_point (const unsigned int q_point) const; + + /** + * The number of degrees of freedom of a single component on the cell for + * the underlying evaluation object. Usually close to + * static_dofs_per_component, but the number depends on the actual element + * selected and is thus not static. + */ + const unsigned int dofs_per_component; + + /** + * The number of degrees of freedom on the cell accumulated over all + * components in the current evaluation object. Usually close to + * static_dofs_per_cell = static_dofs_per_component*n_components, but the + * number depends on the actual element selected and is thus not static. + */ + const unsigned int dofs_per_cell; + + /** + * The number of quadrature points in use. If the number of quadrature + * points in 1d is given as a template, this number is simply the + * dim-1-th power of that value. If the element degree is set to -1 + * (dynamic selection of element degree), the static value of quadrature + * points is inaccurate and this value must be used instead. + */ + const unsigned int n_q_points; + +protected: + + /** + * For faces not oriented in the standard way, this method applies + * re-indexing on quadrature points. Called at the end of evaluate() and at + * the beginning of integrate(). + */ + void adjust_for_face_orientation(const bool integrate, + const bool values, + const bool gradients); +}; + + + +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_component + { + // this division is always without remainder + static constexpr unsigned int value = + (DGP_dofs_per_component::value * (degree+dim)) / dim; + }; + + // base specialization: 1d elements have 'degree+1' degrees of freedom + template + struct DGP_dofs_per_component<1,degree> + { + static constexpr unsigned int value = degree+1; + }; + } +} + + +/*----------------------- Inline functions ----------------------------------*/ + +#ifndef DOXYGEN + + + +/*----------------------- FEEvaluationBase ----------------------------------*/ + +template +inline +FEEvaluationBase +::FEEvaluationBase (const MatrixFree &data_in, + const unsigned int dof_no, + const unsigned int first_selected_component, + const unsigned int quad_no_in, + const unsigned int fe_degree, + const unsigned int n_q_points, + const bool is_interior_face) + : + scratch_data_array (data_in.acquire_scratch_data()), + quad_no (quad_no_in), + n_fe_components (data_in.get_dof_info(dof_no).start_components.back()), + active_fe_index (fe_degree != numbers::invalid_unsigned_int ? + data_in.get_dof_info(dof_no).fe_index_from_degree + (first_selected_component, fe_degree) + : + 0), + active_quad_index (fe_degree != numbers::invalid_unsigned_int ? + (is_face ? + data_in.get_mapping_info().face_data[quad_no_in]. + quad_index_from_n_q_points(n_q_points) + : + data_in.get_mapping_info().cell_data[quad_no_in]. + quad_index_from_n_q_points(n_q_points)) + : + 0), + n_quadrature_points(fe_degree != numbers::invalid_unsigned_int ? n_q_points : + (is_face ? + data_in.get_shape_info + (dof_no, quad_no_in, active_fe_index, active_quad_index).n_q_points_face + : + data_in.get_shape_info + (dof_no, quad_no_in, active_fe_index, active_quad_index).n_q_points)), + matrix_info (&data_in), + dof_info (&data_in.get_dof_info(dof_no)), + mapping_data (internal::MatrixFreeFunctions::MappingInfoCellsOrFaces::get(data_in.get_mapping_info(), quad_no)), + data (&data_in.get_shape_info + (dof_no, quad_no_in, + dof_info->component_to_base_index[first_selected_component], + active_fe_index, active_quad_index)), + jacobian (nullptr), + J_value (nullptr), + normal_vectors (nullptr), + normal_x_jacobian (nullptr), + quadrature_weights (mapping_data->descriptor[active_quad_index].quadrature_weights.begin()), + cell (numbers::invalid_unsigned_int), + is_interior_face (is_interior_face), + face_vector_access_index (is_face ? (is_interior_face ? 0 : 1) : 2), + cell_type (internal::MatrixFreeFunctions::general), + dof_values_initialized (false), + values_quad_initialized (false), + gradients_quad_initialized(false), + hessians_quad_initialized (false), + values_quad_submitted (false), + gradients_quad_submitted (false), + first_selected_component (first_selected_component) +{ + set_data_pointers(); + Assert (matrix_info->mapping_initialized() == true, + ExcNotInitialized()); + AssertDimension (matrix_info->get_size_info().vectorization_length, + VectorizedArray::n_array_elements); + AssertDimension ((is_face ? data->n_q_points_face : data->n_q_points), + n_quadrature_points); + AssertDimension (n_quadrature_points, + mapping_data->descriptor[active_quad_index].n_q_points); + Assert(dof_info->start_components.back() == 1 || + (int)n_components_ <= + (int)dof_info->start_components[dof_info->component_to_base_index[first_selected_component]+1] - first_selected_component, + ExcMessage("You tried to construct a vector-valued evaluator with " + + Utilities::to_string(n_components) + " components. However, " + "the current base element has only " + + Utilities::to_string(dof_info->start_components[dof_info->component_to_base_index[first_selected_component]+1] - first_selected_component) + + " components left when starting from local element index " + + Utilities::to_string(first_selected_component-dof_info->start_components[dof_info->component_to_base_index[first_selected_component]]) + + " (global index " + Utilities::to_string(first_selected_component) + + ")")); + + // do not check for correct dimensions of data fields here, should be done + // in derived classes +} @@ -2443,7 +2828,8 @@ FEEvaluationBase quadrature_weights (nullptr), cell (0), cell_type (internal::MatrixFreeFunctions::general), - is_interior_face (true), + is_interior_face (true), + face_vector_access_index (numbers::invalid_unsigned_int), dof_values_initialized (false), values_quad_initialized (false), gradients_quad_initialized(false), @@ -2508,7 +2894,8 @@ FEEvaluationBase mapping_data->descriptor[active_quad_index].quadrature_weights.begin()), cell (numbers::invalid_unsigned_int), cell_type (internal::MatrixFreeFunctions::general), - is_interior_face (other.is_interior_face), + is_interior_face (other.is_interior_face), + face_vector_access_index (other.face_vector_access_index), dof_values_initialized (false), values_quad_initialized (false), gradients_quad_initialized(false), @@ -2582,6 +2969,7 @@ FEEvaluationBase cell = numbers::invalid_unsigned_int; cell_type = internal::MatrixFreeFunctions::general; is_interior_face = other.is_interior_face; + face_vector_access_index = other.face_vector_access_index; // Create deep copy of mapped geometry for use in parallel... if (other.mapped_geometry.get() != nullptr) @@ -2668,82 +3056,12 @@ FEEvaluationBase template inline -void -FEEvaluationBase::reinit (const unsigned int cell_index) +unsigned int +FEEvaluationBase +::get_cell_data_number () const { - Assert (mapped_geometry == nullptr, - ExcMessage("FEEvaluation was initialized without a matrix-free object." - " Integer indexing is not possible")); - if (mapped_geometry != nullptr) - return; - - Assert (this->dof_info != nullptr, ExcNotInitialized()); - Assert (this->mapping_data != nullptr, ExcNotInitialized()); - this->cell = cell_index; - this->cell_type = this->matrix_info->get_mapping_info().get_cell_type(cell_index); - - const unsigned int offsets = this->mapping_data->data_index_offsets[cell_index]; - this->jacobian = &this->mapping_data->jacobians[0][offsets]; - this->J_value = &this->mapping_data->JxW_values[offsets]; - -#ifdef DEBUG - dof_values_initialized = false; - values_quad_initialized = false; - gradients_quad_initialized = false; - hessians_quad_initialized = false; -#endif -} - - - -template -template -inline -void -FEEvaluationBase -::reinit (const TriaIterator > &cell) -{ - Assert(matrix_info == nullptr, - ExcMessage("Cannot use initialization from cell iterator if " - "initialized from MatrixFree object. Use variant for " - "on the fly computation with arguments as for FEValues " - "instead")); - Assert(mapped_geometry.get() != nullptr, ExcNotInitialized()); - mapped_geometry->reinit(static_cast::cell_iterator>(cell)); - local_dof_indices.resize(cell->get_fe().dofs_per_cell); - if (level_dof_access) - cell->get_mg_dof_indices(local_dof_indices); - else - cell->get_dof_indices(local_dof_indices); -} - - - -template -inline -void -FEEvaluationBase -::reinit (const typename Triangulation::cell_iterator &cell) -{ - Assert(matrix_info == 0, - ExcMessage("Cannot use initialization from cell iterator if " - "initialized from MatrixFree object. Use variant for " - "on the fly computation with arguments as for FEValues " - "instead")); - Assert(mapped_geometry.get() != 0, ExcNotInitialized()); - mapped_geometry->reinit(cell); -} - - - -template -inline -unsigned int -FEEvaluationBase -::get_cell_data_number () const -{ - return get_mapping_data_index_offset(); -} + return get_mapping_data_index_offset(); +} @@ -2885,7 +3203,7 @@ FEEvaluationBase namespace internal { - // write access to generic vectors that have operator (). + // access to generic vectors that have operator (). template inline typename VectorType::value_type & @@ -2897,19 +3215,7 @@ namespace internal - // read access to generic vectors that have operator (). - template - inline - typename VectorType::value_type - vector_access (const VectorType &vec, - const unsigned int entry) - { - return vec(entry); - } - - - - // write access to distributed MPI vectors that have a local_element(uint) + // access to distributed MPI vectors that have a local_element(uint) // method to access data in local index space, which is what we use in // DoFInfo and hence in read_dof_values etc. template @@ -2923,20 +3229,6 @@ namespace internal - // read access to distributed MPI vectors that have a local_element(uint) - // method to access data in local index space, which is what we use in - // DoFInfo and hence in read_dof_values etc. - template - inline - Number - vector_access (const LinearAlgebra::distributed::Vector &vec, - const unsigned int entry) - { - return vec.local_element(entry); - } - - - // this is to make sure that the parallel partitioning in the // LinearAlgebra::distributed::Vector is really the same as stored in // MatrixFree @@ -2975,7 +3267,31 @@ namespace internal VectorType &vec, Number &res) const { - res = vector_access (const_cast(vec), index); + res = vector_access (vec, index); + } + + template + void process_dofs_vectorized_transpose (const unsigned int dofs_per_cell, + const unsigned int *dof_indices, + VectorType &vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + dealii::vectorized_load_and_transpose(dofs_per_cell, vec.begin(), + dof_indices, dof_values); + } + + + template + void process_dofs_vectorized_transpose (const unsigned int dofs_per_cell, + const unsigned int *dof_indices, + VectorType &vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + for (unsigned int d=0; d::n_array_elements; ++v) + dof_values[d][v] = vector_access(vec, dof_indices[v]+d); } // variant where VectorType::value_type is the same as Number -> can call @@ -2983,10 +3299,11 @@ namespace internal template void process_dof_gather (const unsigned int *indices, VectorType &vec, + const unsigned int constant_offset, VectorizedArray &res, std::integral_constant) const { - res.gather(vec.begin(), indices); + res.gather(vec.begin()+constant_offset, indices); } // variant where VectorType::value_type is not the same as Number -> must @@ -2994,11 +3311,12 @@ namespace internal template void process_dof_gather (const unsigned int *indices, VectorType &vec, + const unsigned int constant_offset, VectorizedArray &res, std::integral_constant) const { for (unsigned int v=0; v::n_array_elements; ++v) - res[v] = vector_access(const_cast(vec), indices[v]); + res[v] = vector_access(vec, indices[v]+constant_offset); } template @@ -3021,7 +3339,7 @@ namespace internal VectorType &vec, Number &res) const { - res += weight * vector_access (const_cast(vec), index); + res += weight * vector_access (vec, index); } void post_constraints (const Number &sum, @@ -3048,25 +3366,47 @@ namespace internal vector_access (vec, index) += res; } + template + void process_dofs_vectorized_transpose (const unsigned int dofs_per_cell, + const unsigned int *dof_indices, + VectorType &vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + vectorized_transpose_and_store(true, dofs_per_cell, dof_values, + dof_indices, vec.begin()); + } + + template + void process_dofs_vectorized_transpose (const unsigned int dofs_per_cell, + const unsigned int *dof_indices, + VectorType &vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + for (unsigned int d=0; d::n_array_elements; ++v) + vector_access(vec, dof_indices[v]+d) += dof_values[d][v]; + } + // variant where VectorType::value_type is the same as Number -> can call // scatter template void process_dof_gather (const unsigned int *indices, VectorType &vec, + const unsigned int constant_offset, VectorizedArray &res, std::integral_constant) const { - // TODO: enable scatter path when indices are fixed - //#if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3 -#if 1 +#if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3 for (unsigned int v=0; v::n_array_elements; ++v) - vector_access(vec, indices[v]) += res[v]; + vector_access(vec, indices[v]+constant_offset) += res[v]; #else // only use gather in case there is also scatter. VectorizedArray tmp; - tmp.gather(vec.begin(), indices); + tmp.gather(vec.begin()+constant_offset, indices); tmp += res; - tmp.scatter(indices, vec.begin()); + tmp.scatter(indices, vec.begin()+constant_offset); #endif } @@ -3075,11 +3415,12 @@ namespace internal template void process_dof_gather (const unsigned int *indices, VectorType &vec, + const unsigned int constant_offset, VectorizedArray &res, std::integral_constant) const { for (unsigned int v=0; v::n_array_elements; ++v) - vector_access(vec, indices[v]) += res[v]; + vector_access(vec, indices[v]+constant_offset) += res[v]; } template @@ -3128,23 +3469,48 @@ namespace internal vector_access (vec, index) = res; } + template + void process_dofs_vectorized_transpose (const unsigned int dofs_per_cell, + const unsigned int *dof_indices, + VectorType &vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + vectorized_transpose_and_store(false, dofs_per_cell, dof_values, + dof_indices, vec.begin()); + } + + template + void process_dofs_vectorized_transpose (const unsigned int dofs_per_cell, + const unsigned int *dof_indices, + VectorType &vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + for (unsigned int i=0; i::n_array_elements; ++v) + vector_access(vec, dof_indices[v]+i) = dof_values[i][v]; + } + template void process_dof_gather (const unsigned int *indices, VectorType &vec, + const unsigned int constant_offset, VectorizedArray &res, std::integral_constant) const { - res.scatter(indices, vec.begin()); + res.scatter(indices, vec.begin()+constant_offset); } template void process_dof_gather (const unsigned int *indices, VectorType &vec, + const unsigned int constant_offset, VectorizedArray &res, std::integral_constant) const { for (unsigned int v=0; v::n_array_elements; ++v) - vector_access(vec, indices[v]) = res[v]; + vector_access(vec, indices[v]+constant_offset) = res[v]; } template @@ -3261,7 +3627,8 @@ FEEvaluationBase const bool apply_constraints) const { // Case 1: No MatrixFree object given, simple case because we do not need to - // process constraints and need not care about vectorization + // process constraints and need not care about vectorization -> go to + // separate function if (matrix_info == nullptr) { read_write_operation_global(operation, src); @@ -3271,10 +3638,29 @@ FEEvaluationBase Assert (dof_info != nullptr, ExcNotInitialized()); Assert (matrix_info->indices_initialized() == true, ExcNotInitialized()); + if (n_fe_components == 1) + for (unsigned int comp=0; comp go to separate function + AssertIndexRange(cell, + dof_info->index_storage_variants[face_vector_access_index].size()); + if (dof_info->index_storage_variants[is_face ? face_vector_access_index : 2][cell] >= + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous) + { + read_write_operation_contiguous(operation, src); + return; + } - const unsigned int n_vectorization = VectorizedArray::n_array_elements; + // Case 3: standard operation with one index per degree of freedom -> go on + // here + + constexpr unsigned int n_vectorization = VectorizedArray::n_array_elements; const unsigned int dofs_per_component = this->data->dofs_per_component_on_cell; if (dof_info->index_storage_variants[is_face ? face_vector_access_index : 2][cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::interleaved) @@ -3284,14 +3670,14 @@ FEEvaluationBase if (n_components == 1 || n_fe_components == 1) for (unsigned int i=0; i::value>()); else for (unsigned int comp=0; comp::value>()); return; } @@ -3420,7 +3806,7 @@ FEEvaluationBase unsigned int ind_local = 0; for ( ; index_indicators != next_index_indicators; ++index_indicators) { - std::pair indicator = + const std::pair indicator = dof_info->constraint_indicator[index_indicators]; // run through values up to next constraint for (unsigned int j=0; j // check whether there is any constraint on the current cell for ( ; index_indicators != next_index_indicators; ++index_indicators) { - std::pair indicator = + const std::pair indicator = dof_info->constraint_indicator[index_indicators]; // run through values up to next constraint @@ -3550,6 +3936,71 @@ FEEvaluationBase +template +template +inline +void +FEEvaluationBase +::read_write_operation_contiguous (const VectorOperation &operation, + VectorType *src[]) const +{ + // This functions processes the functions read_dof_values, + // distribute_local_to_global, and set_dof_values with the same code for + // contiguous cell indices (DG case). The distinction between these three + // cases is made by the input VectorOperation that either reads values from + // a vector and puts the data into the local data field or write local data + // into the vector. Certain operations are no-ops for the given use case. + + std::integral_constant::value> + vector_selector; + const unsigned int ind = is_face ? face_vector_access_index : 2; + + const std::vector &dof_indices_cont + = dof_info->dof_indices_contiguous[ind]; + const unsigned int vectorization_populated = + dof_info->n_vectorization_lanes_filled[ind][this->cell]; + unsigned int dof_indices[VectorizedArray::n_array_elements]; + for (unsigned int v=0; v::n_array_elements+v] + + dof_info->component_dof_indices_offset[active_fe_index][first_selected_component]; + + // In the case with contiguous cell indices, we know that there are no + // constraints and that the indices within each element are contiguous + if (vectorization_populated == VectorizedArray::n_array_elements) + { + if (n_components == 1 || n_fe_components == 1) + for (unsigned int comp=0; compdofs_per_component_on_cell, + dof_indices, + *src[comp], values_dofs[comp], + vector_selector); + else + operation.process_dofs_vectorized_transpose(data->dofs_per_component_on_cell* + n_components, + dof_indices, *src[0], + &values_dofs[0][0], + vector_selector); + } + else + for (unsigned int comp=0; compdofs_per_component_on_cell; ++i) + operation.process_empty(values_dofs[comp][i]); + if (n_components == 1 || n_fe_components == 1) + for (unsigned int v=0; vdofs_per_component_on_cell; ++i) + operation.process_dof (dof_indices[v]+i, *src[comp], + values_dofs[comp][i][v]); + else + for (unsigned int v=0; vdofs_per_component_on_cell; ++i) + operation.process_dof (dof_indices[v]+i+comp*data->dofs_per_component_on_cell, + *src[0], values_dofs[comp][i][v]); + } +} + + + template template inline @@ -5609,6 +6060,80 @@ FEEvaluation +template +inline +void +FEEvaluation +::reinit (const unsigned int cell_index) +{ + Assert (this->mapped_geometry == nullptr, + ExcMessage("FEEvaluation was initialized without a matrix-free object." + " Integer indexing is not possible")); + if (this->mapped_geometry != nullptr) + return; + + Assert (this->dof_info != nullptr, ExcNotInitialized()); + Assert (this->mapping_data != nullptr, ExcNotInitialized()); + this->cell = cell_index; + this->cell_type = this->matrix_info->get_mapping_info().get_cell_type(cell_index); + + const unsigned int offsets = this->mapping_data->data_index_offsets[cell_index]; + this->jacobian = &this->mapping_data->jacobians[0][offsets]; + this->J_value = &this->mapping_data->JxW_values[offsets]; + +#ifdef DEBUG + this->dof_values_initialized = false; + this->values_quad_initialized = false; + this->gradients_quad_initialized = false; + this->hessians_quad_initialized = false; +#endif +} + + + +template +template +inline +void +FEEvaluation +::reinit (const TriaIterator > &cell) +{ + Assert(this->matrix_info == nullptr, + ExcMessage("Cannot use initialization from cell iterator if " + "initialized from MatrixFree object. Use variant for " + "on the fly computation with arguments as for FEValues " + "instead")); + Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized()); + this->mapped_geometry->reinit(static_cast::cell_iterator>(cell)); + this->local_dof_indices.resize(cell->get_fe().dofs_per_cell); + if (level_dof_access) + cell->get_mg_dof_indices(this->local_dof_indices); + else + cell->get_dof_indices(this->local_dof_indices); +} + + + +template +inline +void +FEEvaluation +::reinit (const typename Triangulation::cell_iterator &cell) +{ + Assert(this->matrix_info == 0, + ExcMessage("Cannot use initialization from cell iterator if " + "initialized from MatrixFree object. Use variant for " + "on the fly computation with arguments as for FEValues " + "instead")); + Assert(this->mapped_geometry.get() != 0, ExcNotInitialized()); + this->mapped_geometry->reinit(cell); +} + + + template inline @@ -5791,6 +6316,630 @@ FEEvaluation +/*-------------------------- FEFaceEvaluation ---------------------------*/ + + + +template +inline +FEFaceEvaluation +::FEFaceEvaluation (const MatrixFree &matrix_free, + const bool is_interior_face, + const unsigned int dof_no, + const unsigned int quad_no, + const unsigned int first_selected_component) + : + BaseClass(matrix_free, dof_no, first_selected_component, quad_no, fe_degree, + static_n_q_points, is_interior_face), + dofs_per_component (this->data->dofs_per_component_on_cell), + dofs_per_cell (this->data->dofs_per_component_on_cell *n_components_), + n_q_points (this->data->n_q_points_face) +{ +} + + + +template +inline +FEFaceEvaluation +::~FEFaceEvaluation () +{} + + + +template +inline +void +FEFaceEvaluation +::reinit (const unsigned int face_index) +{ + Assert (this->mapped_geometry == nullptr, + ExcMessage("FEEvaluation was initialized without a matrix-free object." + " Integer indexing is not possible")); + if (this->mapped_geometry != nullptr) + return; + + this->cell = face_index; + this->face_vector_access_index = this->is_interior_face ? 0 : 1; + if (face_index >= this->matrix_info->get_task_info().refinement_edge_face_partition_data[0]) + this->face_vector_access_index = 0; + Assert (this->mapping_data != nullptr, ExcNotInitialized()); + const unsigned int n_vectors = VectorizedArray::n_array_elements; + const internal::MatrixFreeFunctions::FaceToCellTopology &faces = + this->matrix_info->get_face_info(face_index); + if (face_index >= this->matrix_info->get_task_info().face_partition_data.back() && + face_index < this->matrix_info->get_task_info().boundary_partition_data.back()) + Assert(this->is_interior_face, ExcMessage("Boundary faces do not have a neighbor")); + + this->face_no = (this->is_interior_face ? faces.interior_face_no : faces.exterior_face_no); + this->subface_index = faces.subface_index; + if (this->is_interior_face == true) + { + this->subface_index = GeometryInfo::max_children_per_cell; + if (faces.face_orientation > 8) + this->face_orientation = faces.face_orientation - 8; + else + this->face_orientation = 0; + } + else + { + if (faces.face_orientation < 8) + this->face_orientation = faces.face_orientation; + else + this->face_orientation = 0; + } + + this->values_quad_submitted = false; + + this->cell_type = this->matrix_info->get_mapping_info().face_type[face_index]; + const unsigned int offsets = this->mapping_data->data_index_offsets[face_index]; + this->J_value = &this->mapping_data->JxW_values[offsets]; + this->normal_vectors = &this->mapping_data->normal_vectors[offsets]; + this->jacobian = &this->mapping_data->jacobians[!this->is_interior_face][offsets]; + this->normal_x_jacobian = &this->mapping_data->normals_times_jacobians[!this->is_interior_face][offsets]; + +#ifdef DEBUG + this->dof_values_initialized = false; + this->values_quad_initialized = false; + this->gradients_quad_initialized = false; + this->hessians_quad_initialized = false; +#endif +} + + + +template +inline +void +FEFaceEvaluation +::reinit (const unsigned int cell_index, + const unsigned int face_number) +{ + Assert(this->quad_no < this->matrix_info->get_mapping_info().face_data_by_cells.size(), + ExcMessage("You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method.")); + AssertIndexRange(face_number, GeometryInfo::faces_per_cell); + AssertIndexRange(cell_index, this->matrix_info->get_mapping_info().cell_type.size()); + Assert (this->mapped_geometry == nullptr, + ExcMessage("FEEvaluation was initialized without a matrix-free object." + " Integer indexing is not possible")); + Assert (this->is_interior_face==true, + ExcMessage("Cell-based FEFaceEvaluation::reinit only possible for the " + "interior face with second argument to constructor as true")); + if (this->mapped_geometry != nullptr) + return; + Assert (this->matrix_info != nullptr, ExcNotInitialized()); + + this->cell_type = this->matrix_info->get_mapping_info().cell_type[cell_index]; + this->cell = cell_index; + this->face_orientation = 0; + this->subface_index = GeometryInfo::max_children_per_cell; + this->face_no = face_number; + this->face_vector_access_index = 2; + + const unsigned int offsets = + this->matrix_info->get_mapping_info().face_data_by_cells[this->quad_no]. + data_index_offsets[cell_index*GeometryInfo::faces_per_cell+face_number]; + AssertIndexRange(offsets, this->matrix_info->get_mapping_info(). + face_data_by_cells[this->quad_no].JxW_values.size()); + this->J_value = &this->matrix_info->get_mapping_info(). + face_data_by_cells[this->quad_no].JxW_values[offsets]; + this->normal_vectors = &this->matrix_info->get_mapping_info(). + face_data_by_cells[this->quad_no].normal_vectors[offsets]; + this->jacobian = &this->matrix_info->get_mapping_info(). + face_data_by_cells[this->quad_no].jacobians[0][offsets]; + this->normal_x_jacobian = &this->matrix_info->get_mapping_info(). + face_data_by_cells[this->quad_no].normals_times_jacobians[0][offsets]; + +#ifdef DEBUG + this->dof_values_initialized = false; + this->values_quad_initialized = false; + this->gradients_quad_initialized = false; + this->hessians_quad_initialized = false; +#endif +} + + + +template +inline +void +FEFaceEvaluation +::evaluate (const bool evaluate_values, + const bool evaluate_gradients) +{ + Assert(this->dof_values_initialized, ExcNotInitialized()); + + evaluate(this->values_dofs[0], evaluate_values, evaluate_gradients); +} + + + +template +inline +void +FEFaceEvaluation +::evaluate (const VectorizedArray *values_array, + const bool evaluate_values, + const bool evaluate_gradients) +{ + if (!(evaluate_values + evaluate_gradients)) + return; + + const unsigned int static_dofs_per_face = fe_degree > -1 ? + Utilities::pow(fe_degree+1,dim-1) : numbers::invalid_unsigned_int; + const unsigned int dofs_per_face = fe_degree > -1 ? + static_dofs_per_face : + Utilities::pow(this->data->fe_degree+1, dim-1); + + // we allocate small amounts of data on the stack to signal the compiler + // that this temporary data is only needed for the calculations but the + // final results can be discarded and need not be written back to + // memory. For large sizes or when the dofs per face is not a compile-time + // constant, however, we want to go to the heap in the `scratch_data` + // variable to not risk a stack overflow. + constexpr unsigned int stack_array_size_threshold = 100; + + VectorizedArray temp_data[static_dofs_per_face < stack_array_size_threshold ? + n_components * 2 * static_dofs_per_face : 1]; + VectorizedArray *temp1; + if (static_dofs_per_face < stack_array_size_threshold) + temp1 = &temp_data[0]; + else + temp1 = this->scratch_data; + + internal::FEFaceNormalEvaluationImpl > + ::template interpolate(*this->data, values_array, temp1, + evaluate_gradients, this->face_no); + + const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0; + if (fe_degree > -1 && + this->subface_index>=GeometryInfo::max_children_per_cell && + this->data->element_type<=internal::MatrixFreeFunctions::tensor_symmetric) + internal::FEFaceEvaluationImpl > + ::evaluate_in_face(*this->data, temp1, this->begin_values(), this->begin_gradients(), + this->scratch_data+2*n_components*dofs_per_face, + evaluate_values, evaluate_gradients, this->subface_index); + else + internal::FEFaceEvaluationImpl > + ::evaluate_in_face(*this->data, temp1, this->begin_values(), this->begin_gradients(), + this->scratch_data+2*n_components*dofs_per_face, + evaluate_values, evaluate_gradients, this->subface_index); + + if (this->face_orientation) + adjust_for_face_orientation(false, evaluate_values, evaluate_gradients); + +#ifdef DEBUG + if (evaluate_values == true) + this->values_quad_initialized = true; + if (evaluate_gradients == true) + this->gradients_quad_initialized = true; +#endif +} + + + +template +inline +void +FEFaceEvaluation +::integrate (const bool integrate_values, + const bool integrate_gradients) +{ + integrate(integrate_values, integrate_gradients, this->values_dofs[0]); + +#ifdef DEBUG + this->dof_values_initialized = true; +#endif +} + + + +template +inline +void +FEFaceEvaluation +::integrate (const bool integrate_values, + const bool integrate_gradients, + VectorizedArray *values_array) +{ + if (!(integrate_values + integrate_gradients)) + return; + + if (this->face_orientation) + adjust_for_face_orientation(true, integrate_values, integrate_gradients); + + const unsigned int static_dofs_per_face = fe_degree > -1 ? + Utilities::pow(fe_degree+1,dim-1) : numbers::invalid_unsigned_int; + const unsigned int dofs_per_face = fe_degree > -1 ? + static_dofs_per_face : + Utilities::pow(this->data->fe_degree+1, dim-1); + + constexpr unsigned int stack_array_size_threshold = 100; + + VectorizedArray temp_data[static_dofs_per_face < stack_array_size_threshold ? + n_components * 2 * static_dofs_per_face : 1]; + VectorizedArray *temp1; + if (static_dofs_per_face < stack_array_size_threshold) + temp1 = &temp_data[0]; + else + temp1 = this->scratch_data; + + const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0; + if (fe_degree > -1 && + this->subface_index>=GeometryInfo::max_children_per_cell && + this->data->element_type<=internal::MatrixFreeFunctions::tensor_symmetric) + internal::FEFaceEvaluationImpl > + ::integrate_in_face(*this->data, temp1, this->begin_values(), + this->begin_gradients(), + this->scratch_data+2*n_components*dofs_per_face, + integrate_values, integrate_gradients, this->subface_index); + else + internal::FEFaceEvaluationImpl > + ::integrate_in_face(*this->data, temp1, this->begin_values(), + this->begin_gradients(), + this->scratch_data+2*n_components*dofs_per_face, + integrate_values, integrate_gradients, this->subface_index); + + internal::FEFaceNormalEvaluationImpl > + ::template interpolate(*this->data, temp1, values_array, integrate_gradients, this->face_no); +} + + + +template +template +inline +void +FEFaceEvaluation +::gather_evaluate (const VectorType &input_vector, + const bool evaluate_values, + const bool evaluate_gradients) +{ + const unsigned int side = this->face_no%2; + const unsigned int dofs_per_face = fe_degree > -1 ? + Utilities::pow(fe_degree+1,dim-1) : + Utilities::pow(this->data->fe_degree+1, dim-1); + + constexpr unsigned int stack_array_size_threshold = 100; + + VectorizedArray temp_data[dofs_per_face < stack_array_size_threshold ? + n_components_ * 2 * dofs_per_face : 1]; + VectorizedArray *__restrict temp1; + if (dofs_per_face < stack_array_size_threshold) + temp1 = &temp_data[0]; + else + temp1 = this->scratch_data; + + internal::VectorReader reader; + + if (this->dof_info->index_storage_variants[this->face_vector_access_index][this->cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous + && + this->dof_info->n_vectorization_lanes_filled[this->face_vector_access_index][this->cell] == + VectorizedArray::n_array_elements + && + ((evaluate_gradients == false && this->data->nodal_at_cell_boundaries == true) || + (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1))) + { + const unsigned int *indices = &this->dof_info->dof_indices_contiguous + [this->face_vector_access_index][this->cell*VectorizedArray::n_array_elements]; + if (evaluate_gradients == true && + this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite) + { + // we know that the gradient weights for the Hermite case on the + // right (side==1) are the negative from the value at the left + // (side==0), so we only read out one of them. + const VectorizedArray grad_weight0 = (side ? -1. : 1.) * + this->data->shape_data_on_face[0][fe_degree+1]; + const VectorizedArray grad_weight1 = (side ? -1. : 1.) * + this->data->shape_data_on_face[0][fe_degree+2]; + AssertDimension(this->data->face_to_cell_index_hermite.size(1), + 2*dofs_per_face); + + const unsigned int *index_array = &this->data->face_to_cell_index_hermite(this->face_no,0); + for (unsigned int i=0; idof_info->component_dof_indices_offset[this->active_fe_index][this->first_selected_component], + temp1[i+2*comp*dofs_per_face], + std::integral_constant::value>()); + VectorizedArray grad; + reader.process_dof_gather(indices, input_vector, + ind2+comp*static_dofs_per_component + + this->dof_info->component_dof_indices_offset[this->active_fe_index][this->first_selected_component], + grad, + std::integral_constant::value>()); + temp1[i+dofs_per_face+2*comp*dofs_per_face] = + grad_weight0 * temp1[i+2*comp*dofs_per_face] + + grad_weight1 * grad; + } + } + } + else + { + AssertDimension(this->data->face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = &this->data->face_to_cell_index_nodal(this->face_no,0); + for (unsigned int i=0; idof_info->component_dof_indices_offset[this->active_fe_index][this->first_selected_component], + temp1[i+comp*2*dofs_per_face], + std::integral_constant::value>()); + } + } + } + else + { + this->read_dof_values(input_vector); + internal::FEFaceNormalEvaluationImpl > + ::template interpolate(*this->data, this->values_dofs[0], temp1, + evaluate_gradients, this->face_no); + } + + if (fe_degree > -1 && + this->subface_index>=GeometryInfo::max_children_per_cell && + this->data->element_type <= internal::MatrixFreeFunctions::tensor_symmetric) + internal::FEFaceEvaluationImpl > + ::evaluate_in_face(*this->data, temp1, this->values_quad[0], + this->gradients_quad[0][0], this->scratch_data + + 2*n_components_*dofs_per_face, + evaluate_values, evaluate_gradients, this->subface_index); + else + internal::FEFaceEvaluationImpl > + ::evaluate_in_face(*this->data, temp1, this->values_quad[0], + this->gradients_quad[0][0], this->scratch_data + + 2*n_components_*dofs_per_face, + evaluate_values, evaluate_gradients, this->subface_index); + + if (this->face_orientation) + adjust_for_face_orientation(false, evaluate_values, evaluate_gradients); + +#ifdef DEBUG + if (evaluate_values == true) + this->values_quad_initialized = true; + if (evaluate_gradients == true) + this->gradients_quad_initialized = true; +#endif +} + + + +template +template +inline +void +FEFaceEvaluation +::integrate_scatter (const bool integrate_values, + const bool integrate_gradients, + VectorType &destination) +{ + const unsigned int side = this->face_no%2; + const unsigned int dofs_per_face = fe_degree > -1 ? + Utilities::pow(fe_degree+1,dim-1) : + Utilities::pow(this->data->fe_degree+1, dim-1); + + constexpr unsigned int stack_array_size_threshold = 100; + + VectorizedArray temp_data[dofs_per_face < stack_array_size_threshold ? + n_components_ * 2 * dofs_per_face : 1]; + VectorizedArray *__restrict temp1; + if (dofs_per_face < stack_array_size_threshold) + temp1 = &temp_data[0]; + else + temp1 = this->scratch_data; + + if (this->face_orientation) + adjust_for_face_orientation(true, integrate_values, integrate_gradients); + if (fe_degree > -1 && + this->subface_index>=GeometryInfo::max_children_per_cell && + this->data->element_type <= internal::MatrixFreeFunctions::tensor_symmetric) + internal::FEFaceEvaluationImpl > + ::integrate_in_face(*this->data, temp1, this->values_quad[0], + this->gradients_quad[0][0], this->scratch_data + + 2*n_components_*dofs_per_face, + integrate_values, integrate_gradients, this->subface_index); + else + internal::FEFaceEvaluationImpl > + ::integrate_in_face(*this->data, temp1, this->values_quad[0], + this->gradients_quad[0][0], this->scratch_data + + 2*n_components_*dofs_per_face, + integrate_values, integrate_gradients, this->subface_index); + +#ifdef DEBUG + this->dof_values_initialized = true; +#endif + + internal::VectorDistributorLocalToGlobal writer; + + if (this->dof_info->index_storage_variants[this->face_vector_access_index][this->cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous + && + this->dof_info->n_vectorization_lanes_filled[this->face_vector_access_index][this->cell] == + VectorizedArray::n_array_elements + && + ((integrate_gradients == false && this->data->nodal_at_cell_boundaries == true) || + (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1))) + { + const unsigned int *indices = &this->dof_info->dof_indices_contiguous + [this->face_vector_access_index][this->cell*VectorizedArray::n_array_elements]; + + if (integrate_gradients == true && + this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite) + { + // we know that the gradient weights for the Hermite case on the + // right (side==1) are the negative from the value at the left + // (side==0), so we only read out one of them. + const VectorizedArray grad_weight0 = (side ? -1. : 1.) * this->data->shape_data_on_face[0][fe_degree+1]; + const VectorizedArray grad_weight1 = (side ? -1. : 1.) * this->data->shape_data_on_face[0][fe_degree+2]; + AssertDimension(this->data->face_to_cell_index_hermite.size(1), + 2*dofs_per_face); + const unsigned int *index_array = &this->data->face_to_cell_index_hermite(this->face_no,0); + for (unsigned int i=0; i val = temp1[i+2*comp*dofs_per_face] + + grad_weight0 * temp1[i+dofs_per_face+2*comp*dofs_per_face]; + VectorizedArray grad = + grad_weight1 * temp1[i+dofs_per_face+2*comp*dofs_per_face]; + writer.process_dof_gather(indices, destination, + comp*static_dofs_per_component+ind1 + + this->dof_info->component_dof_indices_offset[this->active_fe_index][this->first_selected_component], + val, + std::integral_constant::value>()); + writer.process_dof_gather(indices, destination, + comp*static_dofs_per_component+ind2 + + this->dof_info->component_dof_indices_offset[this->active_fe_index][this->first_selected_component], + grad, + std::integral_constant::value>()); + } + } + } + else + { + AssertDimension(this->data->face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = &this->data->face_to_cell_index_nodal(this->face_no,0); + for (unsigned int i=0; idof_info->component_dof_indices_offset[this->active_fe_index][this->first_selected_component], + temp1[i+2*comp*dofs_per_face], + std::integral_constant::value>()); + } + } + } + else + { + internal::FEFaceNormalEvaluationImpl > + ::template interpolate(*this->data, temp1, this->values_dofs[0], + integrate_gradients, this->face_no); + this->distribute_local_to_global(destination); + } +} + + + +template +inline +void +FEFaceEvaluation +::adjust_for_face_orientation(const bool integrate, + const bool values, + const bool gradients) +{ + VectorizedArray *tmp_values = this->scratch_data; + const unsigned int *orientations = + &this->mapping_data->descriptor[this->active_fe_index].face_orientations[this->face_orientation][0]; + for (unsigned int c=0; cvalues_quad[c][q]; + else + for (unsigned int q=0; qvalues_quad[c][orientations[q]]; + for (unsigned int q=0; qvalues_quad[c][q] = tmp_values[q]; + } + if (gradients == true) + for (unsigned int d=0; dgradients_quad[c][d][q]; + else + for (unsigned int q=0; qgradients_quad[c][d][orientations[q]]; + for (unsigned int q=0; qgradients_quad[c][d][q] = tmp_values[q]; + } + } +} + + + +template +inline +Point > +FEFaceEvaluation +::quadrature_point (const unsigned int q) const +{ + AssertIndexRange (q, n_q_points); + if (this->face_vector_access_index < 2) + { + Assert(this->mapping_data->quadrature_point_offsets.empty() == false, + ExcNotImplemented()); + AssertIndexRange(this->cell, this->mapping_data->quadrature_point_offsets.size()); + return this->mapping_data->quadrature_points[this->mapping_data->quadrature_point_offsets[this->cell]+q]; + } + else + { + Assert(this->matrix_info->get_mapping_info().face_data_by_cells + [this->quad_no].quadrature_point_offsets.empty() == false, + ExcNotImplemented()); + const unsigned int index = this->cell*GeometryInfo::faces_per_cell + this->face_no; + AssertIndexRange(index, this->matrix_info->get_mapping_info().face_data_by_cells + [this->quad_no].quadrature_point_offsets.size()); + return this->matrix_info->get_mapping_info().face_data_by_cells[this->quad_no]. + quadrature_points[this->matrix_info->get_mapping_info().face_data_by_cells + [this->quad_no].quadrature_point_offsets[index]+q]; + } +} + + + +/*------------------------- end FEFaceEvaluation ------------------------- */ + + #endif // ifndef DOXYGEN