From ae986771db620b455b8ee5bac6905f966bd70408 Mon Sep 17 00:00:00 2001 From: Katharina Kormann Date: Fri, 20 Apr 2018 21:15:43 +0200 Subject: [PATCH] Introduce face functionality to fe_evaluation. Distinguish cell and face type by new template argument. --- include/deal.II/matrix_free/fe_evaluation.h | 1326 ++++++++++++------ tests/matrix_free/fe_evaluation_access_1d.cc | 2 +- 2 files changed, 883 insertions(+), 445 deletions(-) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index c30dc9faf7..c047df59e6 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -62,7 +62,7 @@ template +template class FEEvaluationBase { public: @@ -164,12 +164,6 @@ public: const internal::MatrixFreeFunctions::ShapeInfo> & get_shape_info() const; - /** - * Fills the JxW values currently used. - */ - void - fill_JxW_values(AlignedVector > &JxW_values) const; - //@} /** @@ -376,10 +370,25 @@ public: */ gradient_type get_gradient (const unsigned int q_point) const; + /** + * Return the derivative of a finite element function at quadrature point + * number @p q_point after a call to @p evaluate(...,true,...) in the + * direction normal to the face: + * $\boldsymbol \nabla u(\mathbf x_q) \cdot \mathbf n(\mathbf x_q)$ + * + * This call is equivalent to calling `get_gradient() * get_normal_vector()` + * but will use a more efficient internal representation of data. + * + * Note that the derived class FEEvaluationAccess overloads this operation + * with specializations for the scalar case (n_components == 1) and for the + * vector-valued case (n_components == dim). + */ + value_type get_normal_derivative (const unsigned int q_point) const; + /** * Write a contribution that is tested by the gradient to the field * containing the values on quadrature points with component @p q_point. - * Access to the same field as through @p get_gradient. If applied before + * Access to the same field as through get_gradient(). If applied before * the function @p integrate(...,true) is called, this specifies what is * tested by all basis function gradients on the current cell and integrated * over. @@ -391,6 +400,21 @@ public: void submit_gradient(const gradient_type grad_in, const unsigned int q_point); + /** + * Write a contribution that is tested by the gradient to the field + * containing the values on quadrature points with component @p + * q_point. Access to the same field as through get_gradient() or + * get_normal_derivative(). If applied before the function @p + * integrate(...,true) is called, this specifies what is tested by all basis + * function gradients on the current cell and integrated over. + * + * Note that the derived class FEEvaluationAccess overloads this operation + * with specializations for the scalar case (n_components == 1) and for the + * vector-valued case (n_components == dim). + */ + void submit_normal_derivative(const value_type grad_in, + const unsigned int q_point); + /** * Return the Hessian of a finite element function at quadrature point * number @p q_point after a call to @p evaluate(...,true). If only the @@ -485,10 +509,16 @@ public: * Return the determinant of the Jacobian from the unit to the real cell * times the quadrature weight. */ - VectorizedArray JxW(const unsigned int q_point) const; + VectorizedArray JxW(const unsigned int q_index) const; /** - * Gets the inverse and transposed version of Jacobian of the mapping + * Fills the JxW values currently used into the given array. + */ + void + fill_JxW_values(AlignedVector > &JxW_values) const; + + /** + * Return the inverse and transposed version of Jacobian of the mapping * between the unit to the real cell (representing the covariant * transformation). This is exactly the matrix used internally to transform * the unit cell gradients to gradients on the real cell. @@ -496,6 +526,30 @@ public: Tensor<2,dim,VectorizedArray > inverse_jacobian(const unsigned int q_index) const; + /** + * Return the unit normal vector on a face. Note that both sides of a face + * use the same orientation of the normal vector: For the faces enumerated + * as `interior` in FaceToCellTopology and selected with the + * `is_interior_face=true` flag of the constructor, this corresponds to the + * outer normal vector, whereas for faces enumerated as `exterior` in + * FaceToCellTopology and selected with the `is_interior_face=false` flag of + * the constructor, the normal points into the element as a consequence of + * the single normal vector. + * + * @note Only implemented in case `is_face == true`. + */ + Tensor<1,dim,VectorizedArray > + get_normal_vector(const unsigned int q_point) const; + + /** + * Provides a unified interface to access data in a vector of + * VectorizedArray fields of length MatrixFree::n_macro_cells() + + * MatrixFree::n_macro_ghost_cells() for both cells (plain read) and faces + * (indirect addressing). + */ + VectorizedArray + read_cell_data (const AlignedVector > &array) const; + //@} /** @@ -625,14 +679,17 @@ protected: * Constructor. Made protected to prevent users from directly using this * class. 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. + * during construction of @p matrix_free, @p dof_no, @p + * first_selected_component and @p quad_no allow to select the appropriate + * components. */ FEEvaluationBase (const MatrixFree &matrix_free, - const unsigned int fe_no, + const unsigned int dof_no, + const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, - const unsigned int n_q_points); + const unsigned int n_q_points, + const bool is_interior_face); /** * Constructor that comes with reduced functionality and works similar as @@ -802,6 +859,11 @@ protected: */ const unsigned int active_quad_index; + /** + * Stores the number of quadrature points in the current evaluation context. + */ + const unsigned int n_quadrature_points; + /** * Stores a pointer to the underlying data. */ @@ -820,7 +882,7 @@ protected: * Also contained in matrix_info, but it simplifies code if we store a * reference to it. */ - const internal::MatrixFreeFunctions::MappingInfoStorage *mapping_data; + const internal::MatrixFreeFunctions::MappingInfoStorage<(is_face?dim-1:dim),dim,Number> *mapping_data; /** * Stores a pointer to the unit cell shape data, i.e., values, gradients and @@ -844,6 +906,16 @@ protected: */ const VectorizedArray *J_value; + /** + * A pointer to the normal vectors at faces. + **/ + const Tensor<1,dim,VectorizedArray > *normal_vectors; + + /** + * A pointer to the normal vectors times the jacobian at faces. + **/ + const Tensor<1,dim,VectorizedArray > *normal_x_jacobian; + /** * A pointer to the quadrature weights of the underlying quadrature formula. */ @@ -855,6 +927,33 @@ protected: */ unsigned int cell; + /** + * Flag holding information whether a face is an interior or exterior face + * according to the defined direction of the normal. Not used for cells. + **/ + bool is_interior_face; + + /** + * Stores the current number of a face within the given cell in case + * `is_face==true`, using values between `0` and `2*dim`. + */ + unsigned int face_no; + + /** + * Stores the orientation of the given face with respect to the standard + * orientation, 0 if in standard orientation. + */ + unsigned int face_orientation; + + /** + * Stores the subface index of the given face. Usually, this variable takes + * the value numbers::invalid_unsigned_int to indicate integration over the + * full face, but in case the current physical face has a neighbor that is + * more refined, it is a subface and must scale the entries in ShapeInfo + * appropriately. + */ + unsigned int subface_index; + /** * Stores the type of the cell we are currently working with after a call to * reinit(). Valid values are @p cartesian, @p affine and @p general, which @@ -912,7 +1011,7 @@ protected: std::shared_ptr > mapped_geometry; /** - * For a FiniteElement with more than one finite element, select at which + * For a FiniteElement with more than one base element, select at which * component this data structure should start. */ const unsigned int first_selected_component; @@ -933,7 +1032,7 @@ private: /** * Make other FEEvaluationBase as well as FEEvaluation objects friends. */ - template friend class FEEvaluationBase; + template friend class FEEvaluationBase; template friend class FEEvaluation; }; @@ -946,8 +1045,8 @@ private: * * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ -template -class FEEvaluationAccess : public FEEvaluationBase +template +class FEEvaluationAccess : public FEEvaluationBase { public: typedef Number number_type; @@ -955,21 +1054,23 @@ public: typedef Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > gradient_type; static constexpr unsigned int dimension = dim; static constexpr unsigned int n_components = n_components_; - typedef FEEvaluationBase BaseClass; + typedef FEEvaluationBase BaseClass; protected: /** * Constructor. Made protected to prevent initialization in user code. 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. + * construction of @p matrix_free, @p first_selected_component and @p + * quad_no allow to select the appropriate components. */ FEEvaluationAccess (const MatrixFree &matrix_free, - const unsigned int fe_no, + const unsigned int dof_no, + const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, - const unsigned int n_q_points); + const unsigned int n_q_points, + const bool is_interior_face = true); /** * Constructor with reduced functionality for similar usage of FEEvaluation @@ -981,7 +1082,7 @@ protected: const Quadrature<1> &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBase *other); + const FEEvaluationBase *other); /** * Copy constructor @@ -1005,15 +1106,15 @@ protected: * * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ -template -class FEEvaluationAccess : public FEEvaluationBase +template +class FEEvaluationAccess : public FEEvaluationBase { public: typedef Number number_type; typedef VectorizedArray value_type; typedef Tensor<1,dim,VectorizedArray > gradient_type; static constexpr unsigned int dimension = dim; - typedef FEEvaluationBase BaseClass; + typedef FEEvaluationBase BaseClass; /** @copydoc FEEvaluationBase::get_dof_value() */ @@ -1033,15 +1134,29 @@ public: void submit_value (const value_type val_in, const unsigned int q_point); + /** @copydoc FEEvaluationBase::submit_value() + */ + void submit_value (const Tensor<1,1,VectorizedArray > val_in, + const unsigned int q_point); + /** @copydoc FEEvaluationBase::get_gradient() */ gradient_type get_gradient (const unsigned int q_point) const; + /** @copydoc FEEvaluationBase::get_normal_derivative() + */ + value_type get_normal_derivative (const unsigned int q_point) const; + /** @copydoc FEEvaluationBase::submit_gradient() */ void submit_gradient(const gradient_type grad_in, const unsigned int q_point); + /** @copydoc FEEvaluationBase::submit_normal_derivative() + */ + void submit_normal_derivative(const value_type grad_in, + const unsigned int q_point); + /** @copydoc FEEvaluationBase::get_hessian() */ Tensor<2,dim,VectorizedArray > @@ -1064,14 +1179,16 @@ protected: * Constructor. Made protected to avoid initialization in user code. 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. + * construction of @p matrix_free, @p first_selected_component and @p + * quad_no allow to select the appropriate components. */ FEEvaluationAccess (const MatrixFree &matrix_free, - const unsigned int fe_no, + const unsigned int dof_no, + const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, - const unsigned int n_q_points); + const unsigned int n_q_points, + const bool is_interior_face = true); /** * Constructor with reduced functionality for similar usage of FEEvaluation @@ -1083,7 +1200,7 @@ protected: const Quadrature<1> &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBase *other); + const FEEvaluationBase *other); /** * Copy constructor @@ -1107,8 +1224,8 @@ protected: * * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 */ -template -class FEEvaluationAccess : public FEEvaluationBase +template +class FEEvaluationAccess : public FEEvaluationBase { public: typedef Number number_type; @@ -1116,7 +1233,7 @@ public: typedef Tensor<2,dim,VectorizedArray > gradient_type; static constexpr unsigned int dimension = dim; static constexpr unsigned int n_components = dim; - typedef FEEvaluationBase BaseClass; + typedef FEEvaluationBase BaseClass; /** @copydoc FEEvaluationBase::get_gradient() */ @@ -1203,14 +1320,16 @@ protected: * Constructor. Made protected to avoid initialization in user code. 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. + * construction of @p matrix_free, @p first_selected_component and @p + * quad_no allow to select the appropriate components. */ FEEvaluationAccess (const MatrixFree &matrix_free, - const unsigned int fe_no, + const unsigned int dof_no, + const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int dofs_per_cell, - const unsigned int n_q_points); + const unsigned int n_q_points, + const bool is_interior_face = true); /** * Constructor with reduced functionality for similar usage of FEEvaluation @@ -1222,7 +1341,7 @@ protected: const Quadrature<1> &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBase *other); + const FEEvaluationBase *other); /** * Copy constructor @@ -1245,15 +1364,15 @@ protected: * @author Katharina Kormann and Martin Kronbichler, 2010, 2011, Shiva * Rudraraju, 2014 */ -template -class FEEvaluationAccess<1,1,Number> : public FEEvaluationBase<1,1,Number> +template +class FEEvaluationAccess<1,1,Number,is_face> : public FEEvaluationBase<1,1,Number,is_face> { public: typedef Number number_type; typedef VectorizedArray value_type; typedef Tensor<1,1,VectorizedArray > gradient_type; - static constexpr unsigned int dimension = 1; - typedef FEEvaluationBase<1,1,Number> BaseClass; + static constexpr unsigned int dimension = 1; + typedef FEEvaluationBase<1,1,Number,is_face> BaseClass; /** @copydoc FEEvaluationBase<1,1,Number>::get_dof_value() */ @@ -1273,15 +1392,39 @@ public: void submit_value (const value_type val_in, const unsigned int q_point); + /** @copydoc FEEvaluationBase<1,1,Number>::submit_value() + */ + void submit_value (const gradient_type val_in, + const unsigned int q_point); + /** @copydoc FEEvaluationBase<1,1,Number>::get_gradient() */ gradient_type get_gradient (const unsigned int q_point) const; + /** @copydoc FEEvaluationBase::get_normal_derivative() + */ + value_type get_normal_derivative (const unsigned int q_point) const; + /** @copydoc FEEvaluationBase<1,1,Number>::submit_gradient() */ void submit_gradient(const gradient_type grad_in, const unsigned int q_point); + /** @copydoc FEEvaluationBase<1,1,Number>::submit_gradient() + */ + void submit_gradient(const value_type grad_in, + const unsigned int q_point); + + /** @copydoc FEEvaluationBase<1,1,Number>::submit_normal_derivative() + */ + void submit_normal_derivative(const value_type grad_in, + const unsigned int q_point); + + /** @copydoc FEEvaluationBase<1,1,Number>::submit_normal_derivative() + */ + void submit_normal_derivative(const gradient_type grad_in, + const unsigned int q_point); + /** @copydoc FEEvaluationBase<1,1,Number>::get_hessian() */ Tensor<2,1,VectorizedArray > @@ -1304,14 +1447,16 @@ protected: * Constructor. Made protected to avoid initialization in user code. 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. + * construction of @p matrix_free, @p first_selected_component and @p + * quad_no allow to select the appropriate components. */ FEEvaluationAccess (const MatrixFree<1,Number> &matrix_free, - const unsigned int fe_no, + const unsigned int dof_no, + const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, - const unsigned int n_q_points); + const unsigned int n_q_points, + const bool is_interior_face = true); /** * Constructor with reduced functionality for similar usage of FEEvaluation @@ -1323,7 +1468,7 @@ protected: const Quadrature<1> &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBase<1,n_components_other,Number> *other); + const FEEvaluationBase<1,n_components_other,Number,is_face> *other); /** * Copy constructor @@ -1882,10 +2027,10 @@ protected: */ template -class FEEvaluation : public FEEvaluationAccess +class FEEvaluation : public FEEvaluationAccess { public: - typedef FEEvaluationAccess BaseClass; + typedef FEEvaluationAccess BaseClass; typedef Number number_type; typedef typename BaseClass::value_type value_type; typedef typename BaseClass::gradient_type gradient_type; @@ -1899,12 +2044,32 @@ public: /** * 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. + * 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 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 component, this parameter + * allows for selecting the component where the current evaluation routine + * should start. Note that one evaluator does not support combining + * different shape functions in different components. In other words, the + * same base element of a FESystem needs to be set for the components + * between @p first_selected_component and + * first_selected_component+n_components_. */ FEEvaluation (const MatrixFree &matrix_free, - const unsigned int fe_no = 0, - const unsigned int quad_no = 0); + const unsigned int dof_no = 0, + const unsigned int quad_no = 0, + const unsigned int first_selected_component = 0); /** * Constructor that comes with reduced functionality and works similar as @@ -2097,10 +2262,11 @@ public: const unsigned int dofs_per_cell; /** - * The number of quadrature points in use for the current evaluation - * object. It equals static_n_q_points in case the element degree is not set - * to -1, i.e., the worker kernels can be set according to the templates @p - * fe_degree and @p n_q_points_1d rather than at run time. + * 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-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; @@ -2147,38 +2313,53 @@ namespace internal /*----------------------- FEEvaluationBase ----------------------------------*/ -template +template inline -FEEvaluationBase +FEEvaluationBase ::FEEvaluationBase (const MatrixFree &data_in, - const unsigned int fe_no_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 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(fe_no_in).n_components), + n_fe_components (data_in.get_dof_info(dof_no).n_components), active_fe_index (fe_degree != numbers::invalid_unsigned_int ? - data_in.get_dof_info(fe_no_in).fe_index_from_degree(fe_degree) + data_in.get_dof_info(dof_no).fe_index_from_degree(fe_degree) : 0), active_quad_index (fe_degree != numbers::invalid_unsigned_int ? - (data_in.get_mapping_info().cell_data[quad_no_in]. + (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(fe_no_in)), - mapping_data (internal::MatrixFreeFunctions::MappingInfoCellsOrFaces::get - (data_in.get_mapping_info(), quad_no)), + 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 - (fe_no_in, quad_no_in, active_fe_index, + (dof_no, quad_no_in, 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), @@ -2186,16 +2367,16 @@ FEEvaluationBase hessians_quad_initialized (false), values_quad_submitted (false), gradients_quad_submitted (false), - first_selected_component (0) + 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 (data->dofs_per_component_on_cell*n_fe_components, - dof_info->dofs_per_cell[active_fe_index]); - AssertDimension (data->n_q_points, + 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 (n_fe_components == 1 || n_components == 1 || @@ -2211,10 +2392,10 @@ FEEvaluationBase -template +template template inline -FEEvaluationBase +FEEvaluationBase ::FEEvaluationBase (const Mapping &mapping, const FiniteElement &fe, const Quadrature<1> &quadrature, @@ -2227,16 +2408,20 @@ FEEvaluationBase n_fe_components (n_components_), active_fe_index (numbers::invalid_unsigned_int), active_quad_index (numbers::invalid_unsigned_int), + n_quadrature_points(Utilities::fixed_power(quadrature.size())), matrix_info (nullptr), dof_info (nullptr), mapping_data (nullptr), // select the correct base element from the given FE component - data (new internal::MatrixFreeFunctions::ShapeInfo>(quadrature, fe, fe.component_to_base_index(first_selected_component).first)), + data (new internal::MatrixFreeFunctions::ShapeInfo> (quadrature, fe, fe.component_to_base_index(first_selected_component).first)), jacobian (nullptr), J_value (nullptr), + normal_vectors (nullptr), + normal_x_jacobian (nullptr), quadrature_weights (nullptr), cell (0), cell_type (internal::MatrixFreeFunctions::general), + is_interior_face (true), dof_values_initialized (false), values_quad_initialized (false), gradients_quad_initialized(false), @@ -2245,10 +2430,8 @@ FEEvaluationBase gradients_quad_submitted (false), // keep the number of the selected component within the current base element // for reading dof values - first_selected_component (fe.component_to_base_index(first_selected_component).second) + first_selected_component (first_selected_component) { - const unsigned int base_element_number = - fe.component_to_base_index(first_selected_component).first; set_data_pointers(); Assert(other == nullptr || other->mapped_geometry.get() != nullptr, ExcInternalError()); @@ -2259,11 +2442,14 @@ FEEvaluationBase mapped_geometry = std::make_shared > (mapping, quadrature, update_flags); + cell = 0; mapping_data = &mapped_geometry->get_data_storage(); jacobian = mapped_geometry->get_data_storage().jacobians[0].begin(); J_value = mapped_geometry->get_data_storage().JxW_values.begin(); + const unsigned int base_element_number = + fe.component_to_base_index(first_selected_component).first; Assert(fe.element_multiplicity(base_element_number) == 1 || fe.element_multiplicity(base_element_number)-first_selected_component >= n_components_, ExcMessage("The underlying element must at least contain as many " @@ -2273,10 +2459,10 @@ FEEvaluationBase -template +template inline -FEEvaluationBase -::FEEvaluationBase (const FEEvaluationBase &other) +FEEvaluationBase +::FEEvaluationBase (const FEEvaluationBase &other) : scratch_data_array (other.matrix_info == nullptr ? new AlignedVector >() : @@ -2285,6 +2471,7 @@ FEEvaluationBase n_fe_components (other.n_fe_components), active_fe_index (other.active_fe_index), active_quad_index (other.active_quad_index), + n_quadrature_points(other.n_quadrature_points), matrix_info (other.matrix_info), dof_info (other.dof_info), mapping_data (other.mapping_data), @@ -2293,10 +2480,13 @@ FEEvaluationBase other.data), jacobian (nullptr), J_value (nullptr), + normal_vectors (nullptr), + normal_x_jacobian (nullptr), quadrature_weights (other.matrix_info == nullptr ? nullptr : 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), dof_values_initialized (false), values_quad_initialized (false), gradients_quad_initialized(false), @@ -2325,11 +2515,11 @@ FEEvaluationBase -template +template inline -FEEvaluationBase & -FEEvaluationBase -::operator= (const FEEvaluationBase &other) +FEEvaluationBase & +FEEvaluationBase +::operator= (const FEEvaluationBase &other) { AssertDimension(quad_no, other.quad_no); AssertDimension(n_fe_components, other.n_fe_components); @@ -2363,14 +2553,13 @@ FEEvaluationBase } set_data_pointers(); - jacobian = nullptr; - J_value = nullptr; quadrature_weights = (mapping_data != nullptr ? mapping_data->descriptor[active_quad_index].quadrature_weights.begin() : nullptr); cell = numbers::invalid_unsigned_int; cell_type = internal::MatrixFreeFunctions::general; + is_interior_face = other.is_interior_face; // Create deep copy of mapped geometry for use in parallel... if (other.mapped_geometry.get() != nullptr) @@ -2391,9 +2580,9 @@ FEEvaluationBase -template +template inline -FEEvaluationBase::~FEEvaluationBase () +FEEvaluationBase::~FEEvaluationBase () { if (matrix_info != nullptr) { @@ -2408,15 +2597,17 @@ FEEvaluationBase::~FEEvaluationBase () { delete scratch_data_array; delete data; + data = nullptr; } + scratch_data_array = nullptr; } -template +template inline void -FEEvaluationBase +FEEvaluationBase ::set_data_pointers() { Assert(scratch_data_array != nullptr, ExcInternalError()); @@ -2424,7 +2615,7 @@ FEEvaluationBase const unsigned int tensor_dofs_per_component = Utilities::fixed_power(this->data->fe_degree+1); const unsigned int dofs_per_component = this->data->dofs_per_component_on_cell; - const unsigned int n_quadrature_points = this->data->n_q_points; + const unsigned int n_quadrature_points = is_face ? this->data->n_q_points_face : this->data->n_q_points; const unsigned int shift = std::max(tensor_dofs_per_component+1, dofs_per_component)* n_components_*3 + 2 * n_quadrature_points; @@ -2453,10 +2644,10 @@ FEEvaluationBase -template +template inline void -FEEvaluationBase::reinit (const unsigned int cell_index) +FEEvaluationBase::reinit (const unsigned int cell_index) { Assert (mapped_geometry == nullptr, ExcMessage("FEEvaluation was initialized without a matrix-free object." @@ -2483,11 +2674,11 @@ FEEvaluationBase::reinit (const unsigned int cell_inde -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::reinit (const TriaIterator > &cell) { Assert(matrix_info == nullptr, @@ -2506,10 +2697,10 @@ FEEvaluationBase -template +template inline void -FEEvaluationBase +FEEvaluationBase ::reinit (const typename Triangulation::cell_iterator &cell) { Assert(matrix_info == 0, @@ -2523,10 +2714,10 @@ FEEvaluationBase -template +template inline unsigned int -FEEvaluationBase +FEEvaluationBase ::get_cell_data_number () const { return get_mapping_data_index_offset(); @@ -2534,10 +2725,10 @@ FEEvaluationBase -template +template inline unsigned int -FEEvaluationBase +FEEvaluationBase ::get_mapping_data_index_offset () const { if (matrix_info == 0) @@ -2551,10 +2742,10 @@ FEEvaluationBase -template +template inline internal::MatrixFreeFunctions::GeometryType -FEEvaluationBase::get_cell_type () const +FEEvaluationBase::get_cell_type () const { Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); return cell_type; @@ -2562,10 +2753,10 @@ FEEvaluationBase::get_cell_type () const -template +template inline const internal::MatrixFreeFunctions::ShapeInfo> & - FEEvaluationBase::get_shape_info() const + FEEvaluationBase::get_shape_info() const { Assert(data != nullptr, ExcInternalError()); return *data; @@ -2573,55 +2764,68 @@ const internal::MatrixFreeFunctions::ShapeInfo> & -template +template inline void -FEEvaluationBase +FEEvaluationBase ::fill_JxW_values(AlignedVector > &JxW_values) const { - AssertDimension(JxW_values.size(), data->n_q_points); - Assert (this->J_value != nullptr, ExcNotInitialized()); - if (this->cell_type == internal::MatrixFreeFunctions::cartesian || - this->cell_type == internal::MatrixFreeFunctions::affine) + AssertDimension(JxW_values.size(), n_quadrature_points); + Assert (J_value != nullptr, ExcNotInitialized()); + if (this->cell_type <= internal::MatrixFreeFunctions::affine) { - Assert (this->mapping_data != nullptr, ExcNotImplemented()); - VectorizedArray J = this->J_value[0]; - for (unsigned int q=0; qdata->n_q_points; ++q) + VectorizedArray J = J_value[0]; + for (unsigned int q=0; qn_quadrature_points; ++q) JxW_values[q] = J * this->quadrature_weights[q]; } else - for (unsigned int q=0; qn_q_points; ++q) - JxW_values[q] = this->J_value[q]; + for (unsigned int q=0; q -inline +template +inline DEAL_II_ALWAYS_INLINE +Tensor<1,dim,VectorizedArray > +FEEvaluationBase +::get_normal_vector(const unsigned int q_index) const +{ + AssertIndexRange(q_index, n_quadrature_points); + Assert(normal_vectors != nullptr, ExcMessage("Did not call reinit()!")); + if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces) + return normal_vectors[0]; + else + return normal_vectors[q_index]; +} + + + +template +inline DEAL_II_ALWAYS_INLINE VectorizedArray -FEEvaluationBase::JxW(const unsigned int q_point) const +FEEvaluationBase::JxW(const unsigned int q_index) const { - AssertIndexRange(q_point, data->n_q_points); - Assert (this->J_value != nullptr, ExcNotInitialized()); - if (this->cell_type == internal::MatrixFreeFunctions::cartesian || - this->cell_type == internal::MatrixFreeFunctions::affine) + AssertIndexRange(q_index, n_quadrature_points); + Assert (J_value != nullptr, ExcNotInitialized()); + if (this->cell_type <= internal::MatrixFreeFunctions::affine) { Assert (this->quadrature_weights != nullptr, ExcInternalError()); - return this->J_value[0] * this->quadrature_weights[q_point]; + return J_value[0] * this->quadrature_weights[q_index]; } else - return this->J_value[q_point]; + return J_value[q_index]; } -template +template inline Tensor<2,dim,VectorizedArray > -FEEvaluationBase +FEEvaluationBase ::inverse_jacobian(const unsigned int q_index) const { - AssertIndexRange(q_index, data->n_q_points); + AssertIndexRange(q_index, n_quadrature_points); Assert (this->jacobian != nullptr, ExcNotImplemented()); if (this->cell_type <= internal::MatrixFreeFunctions::affine) return jacobian[0]; @@ -2631,6 +2835,32 @@ FEEvaluationBase +template +inline +VectorizedArray +FEEvaluationBase +::read_cell_data(const AlignedVector > &array) const +{ + Assert(matrix_info != nullptr, ExcNotImplemented()); + AssertDimension(array.size(), matrix_info->get_task_info().cell_partition_data.back()); + if (is_face) + { + VectorizedArray out = make_vectorized_array(Number(1.)); + const unsigned int *cells = + is_interior_face ? + &this->matrix_info->get_face_info(cell).cells_interior[0] : + &this->matrix_info->get_face_info(cell).cells_exterior[0]; + for (unsigned int i=0; i::n_array_elements; ++i) + if (cells[i] != numbers::invalid_unsigned_int) + out[i] = array[cells[i]/VectorizedArray::n_array_elements][cells[i]%VectorizedArray::n_array_elements]; + return out; + } + else + return array[cell]; +} + + + namespace internal { // write access to generic vectors that have operator (). @@ -2999,11 +3229,11 @@ namespace internal -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_write_operation (const VectorOperation &operation, VectorType *src[]) const { @@ -3034,6 +3264,7 @@ FEEvaluationBase return; } + // Some standard checks Assert (dof_info != nullptr, ExcNotInitialized()); Assert (matrix_info->indices_initialized() == true, ExcNotInitialized()); @@ -3354,11 +3585,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_dof_values (const VectorType &src, const unsigned int first_index) { @@ -3379,11 +3610,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_dof_values_plain (const VectorType &src, const unsigned int first_index) { @@ -3399,11 +3630,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::distribute_local_to_global (VectorType &dst, const unsigned int first_index) const { @@ -3423,11 +3654,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::set_dof_values (VectorType &dst, const unsigned int first_index) const { @@ -3447,11 +3678,11 @@ FEEvaluationBase -template +template template inline void -FEEvaluationBase +FEEvaluationBase ::read_dof_values_plain (const VectorType *src[]) { // Case without MatrixFree initialization object @@ -3585,10 +3816,10 @@ FEEvaluationBase /*------------------------------ access to data fields ----------------------*/ -template +template inline const std::vector & -FEEvaluationBase:: +FEEvaluationBase:: get_internal_dof_numbering() const { return data->lexicographic_numbering; @@ -3596,10 +3827,10 @@ get_internal_dof_numbering() const -template +template inline ArrayView > -FEEvaluationBase:: +FEEvaluationBase:: get_scratch_data() const { return ArrayView >(const_cast *>(scratch_data), @@ -3609,10 +3840,10 @@ get_scratch_data() const -template +template inline const VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_dof_values () const { return &values_dofs[0][0]; @@ -3620,10 +3851,10 @@ begin_dof_values () const -template +template inline VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_dof_values () { #ifdef DEBUG @@ -3634,10 +3865,10 @@ begin_dof_values () -template +template inline const VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_values () const { Assert (values_quad_initialized || values_quad_submitted, @@ -3647,13 +3878,14 @@ begin_values () const -template +template inline VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_values () { #ifdef DEBUG + values_quad_initialized = true; values_quad_submitted = true; #endif return &values_quad[0][0]; @@ -3661,10 +3893,10 @@ begin_values () -template +template inline const VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_gradients () const { Assert (gradients_quad_initialized || gradients_quad_submitted, @@ -3674,24 +3906,25 @@ begin_gradients () const -template +template inline VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_gradients () { #ifdef DEBUG gradients_quad_submitted = true; + gradients_quad_initialized = true; #endif return &gradients_quad[0][0][0]; } -template +template inline const VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_hessians () const { Assert (hessians_quad_initialized, ExcNotInitialized()); @@ -3700,21 +3933,24 @@ begin_hessians () const -template +template inline VectorizedArray * -FEEvaluationBase:: +FEEvaluationBase:: begin_hessians () { +#ifdef DEBUG + hessians_quad_initialized = true; +#endif return &hessians_quad[0][0][0]; } -template -inline +template +inline DEAL_II_ALWAYS_INLINE Tensor<1,n_components_,VectorizedArray > -FEEvaluationBase +FEEvaluationBase ::get_dof_value (const unsigned int dof) const { AssertIndexRange (dof, this->data->dofs_per_component_on_cell); @@ -3726,15 +3962,15 @@ FEEvaluationBase -template -inline +template +inline DEAL_II_ALWAYS_INLINE Tensor<1,n_components_,VectorizedArray > -FEEvaluationBase +FEEvaluationBase ::get_value (const unsigned int q_point) const { Assert (this->values_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); Tensor<1,n_components_,VectorizedArray > return_value; for (unsigned int comp=0; compvalues_quad[comp][q_point]; @@ -3743,20 +3979,22 @@ FEEvaluationBase -template -inline +template +inline DEAL_II_ALWAYS_INLINE Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > -FEEvaluationBase +FEEvaluationBase ::get_gradient (const unsigned int q_point) const { Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); + + Assert(jacobian != nullptr, ExcNotInitialized()); Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > grad_out; // Cartesian cell - if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { for (unsigned int comp=0; comp else { const Tensor<2,dim,VectorizedArray > &jac = - this->cell_type > internal::MatrixFreeFunctions::affine ? - jacobian[q_point] : jacobian[0]; + jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ? q_point : 0]; + for (unsigned int comp=0; compgradients_quad[comp][0][q_point]; + for (unsigned int e=1; egradients_quad[comp][e][q_point]; + } + } + return grad_out; +} + + + +template +inline DEAL_II_ALWAYS_INLINE +Tensor<1,n_components_,VectorizedArray > +FEEvaluationBase +::get_normal_derivative (const unsigned int q_point) const +{ + AssertIndexRange (q_point, this->n_quadrature_points); + Assert (this->gradients_quad_initialized==true, + internal::ExcAccessToUninitializedField()); + + Assert(normal_x_jacobian != nullptr, ExcNotInitialized()); + + Tensor<1,n_components,VectorizedArray > grad_out; + if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + for (unsigned int comp=0; compgradients_quad[comp][dim-1][q_point] * + (this->normal_x_jacobian[0][dim-1]); + else + { + const unsigned int index = + this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point; for (unsigned int comp=0; compgradients_quad[comp][0][q_point]); - for (unsigned int e=1; egradients_quad[comp][e][q_point]); - } + grad_out[comp] = this->gradients_quad[comp][0][q_point] * + this->normal_x_jacobian[index][0]; + for (unsigned int d=1; dgradients_quad[comp][d][q_point] * + this->normal_x_jacobian[index][d]; } } return grad_out; @@ -3843,15 +4113,16 @@ namespace internal -template +template inline Tensor<1,n_components_,Tensor<2,dim,VectorizedArray > > -FEEvaluationBase +FEEvaluationBase ::get_hessian (const unsigned int q_point) const { + Assert(!is_face, ExcNotImplemented()); Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); Assert(jacobian != nullptr, ExcNotImplemented()); const Tensor<2,dim,VectorizedArray > &jac = @@ -3923,7 +4194,7 @@ FEEvaluationBase else { const Tensor<1,dim*(dim+1)/2,Tensor<1,dim,VectorizedArray > > &jac_grad = - mapping_data->jacobian_gradients[0][this->get_mapping_data_index_offset()+q_point]; + mapping_data->jacobian_gradients[1-this->is_interior_face][this->mapping_data->data_index_offsets[this->cell]+q_point]; for (unsigned int comp=0; comp -template +template inline Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > -FEEvaluationBase +FEEvaluationBase ::get_hessian_diagonal (const unsigned int q_point) const { + Assert(!is_face, ExcNotImplemented()); Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); Assert(jacobian != nullptr, ExcNotImplemented()); const Tensor<2,dim,VectorizedArray > &jac = @@ -4014,7 +4286,7 @@ FEEvaluationBase else { const Tensor<1,dim*(dim+1)/2,Tensor<1,dim,VectorizedArray > > &jac_grad = - mapping_data->jacobian_gradients[0][this->get_mapping_data_index_offset()+q_point]; + mapping_data->jacobian_gradients[0][this->mapping_data->data_index_offsets[this->cell]+q_point]; for (unsigned int comp=0; comp -template +template inline Tensor<1,n_components_,VectorizedArray > -FEEvaluationBase +FEEvaluationBase ::get_laplacian (const unsigned int q_point) const { + Assert (is_face == false, ExcNotImplemented()); Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); + Tensor<1,n_components_,VectorizedArray > laplacian_out; const Tensor<1,n_components_,Tensor<1,dim,VectorizedArray > > hess_diag = get_hessian_diagonal(q_point); @@ -4066,10 +4340,10 @@ FEEvaluationBase -template -inline +template +inline DEAL_II_ALWAYS_INLINE void -FEEvaluationBase +FEEvaluationBase ::submit_dof_value (const Tensor<1,n_components_,VectorizedArray > val_in, const unsigned int dof) { @@ -4083,27 +4357,29 @@ FEEvaluationBase -template -inline +template +inline DEAL_II_ALWAYS_INLINE void -FEEvaluationBase +FEEvaluationBase ::submit_value (const Tensor<1,n_components_,VectorizedArray > val_in, const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); + Assert (this->J_value != nullptr, ExcNotInitialized()); this->values_quad_submitted = true; #endif - if (this->cell_type == internal::MatrixFreeFunctions::general) + + if (this->cell_type <= internal::MatrixFreeFunctions::affine) { - const VectorizedArray JxW = J_value[q_point]; + const VectorizedArray JxW = J_value[0] * quadrature_weights[q_point]; for (unsigned int comp=0; compvalues_quad[comp][q_point] = val_in[comp] * JxW; } - else //if (this->cell_type < internal::MatrixFreeFunctions::general) + else { - const VectorizedArray JxW = J_value[0] * quadrature_weights[q_point]; + const VectorizedArray JxW = J_value[q_point]; for (unsigned int comp=0; compvalues_quad[comp][q_point] = val_in[comp] * JxW; } @@ -4111,23 +4387,23 @@ FEEvaluationBase -template -inline +template +inline DEAL_II_ALWAYS_INLINE void -FEEvaluationBase +FEEvaluationBase ::submit_gradient (const Tensor<1,n_components_, Tensor<1,dim,VectorizedArray > >grad_in, const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); this->gradients_quad_submitted = true; Assert (this->J_value != nullptr, ExcNotInitialized()); Assert (this->jacobian != nullptr, ExcNotInitialized()); #endif - if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { const VectorizedArray JxW = J_value[0] * quadrature_weights[q_point]; for (unsigned int comp=0; comp -template +template +inline DEAL_II_ALWAYS_INLINE +void +FEEvaluationBase +::submit_normal_derivative (const Tensor<1,n_components_,VectorizedArray > grad_in, + const unsigned int q_point) +{ +#ifdef DEBUG + AssertIndexRange (q_point, this->n_quadrature_points); + this->gradients_quad_submitted = true; + Assert (this->normal_x_jacobian != nullptr, ExcNotInitialized()); +#endif + + if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + for (unsigned int comp=0; compgradients_quad[comp][d][q_point] = VectorizedArray(); + this->gradients_quad[comp][dim-1][q_point] = + grad_in[comp] * (this->normal_x_jacobian[0][dim-1] * + this->J_value[0] * this->quadrature_weights[q_point]); + } + else + { + const unsigned int index = + this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point; + for (unsigned int comp=0; comp factor = grad_in[comp] * + this->J_value[index]; + if (this->cell_type <= internal::MatrixFreeFunctions::affine) + factor = factor * this->quadrature_weights[q_point]; + for (unsigned int d=0; dgradients_quad[comp][d][q_point] = factor * + this->normal_x_jacobian[index][d]; + } + } +} + + + + +template inline Tensor<1,n_components_,VectorizedArray > -FEEvaluationBase +FEEvaluationBase ::integrate_value () const { #ifdef DEBUG @@ -4170,7 +4488,7 @@ FEEvaluationBase Tensor<1,n_components_,VectorizedArray > return_value; 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->n_quadrature_points; for (unsigned int q=1; qvalues_quad[comp][q]; @@ -4182,55 +4500,58 @@ FEEvaluationBase /*----------------------- FEEvaluationAccess --------------------------------*/ -template +template inline -FEEvaluationAccess +FEEvaluationAccess ::FEEvaluationAccess (const MatrixFree &data_in, - const unsigned int fe_no, + 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 unsigned int n_q_points, + const bool is_interior_face) : - FEEvaluationBase - (data_in, fe_no, quad_no_in, fe_degree, n_q_points) + FEEvaluationBase + (data_in, dof_no, first_selected_component, quad_no_in, fe_degree, n_q_points, + is_interior_face) {} -template +template template inline -FEEvaluationAccess +FEEvaluationAccess ::FEEvaluationAccess (const Mapping &mapping, const FiniteElement &fe, const Quadrature<1> &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBase *other) + const FEEvaluationBase *other) : - FEEvaluationBase (mapping, fe, quadrature, update_flags, - first_selected_component, other) + FEEvaluationBase (mapping, fe, quadrature, update_flags, + first_selected_component, other) {} -template +template inline -FEEvaluationAccess -::FEEvaluationAccess (const FEEvaluationAccess &other) +FEEvaluationAccess +::FEEvaluationAccess (const FEEvaluationAccess &other) : - FEEvaluationBase (other) + FEEvaluationBase (other) {} -template +template inline -FEEvaluationAccess & -FEEvaluationAccess -::operator= (const FEEvaluationAccess &other) +FEEvaluationAccess & +FEEvaluationAccess +::operator= (const FEEvaluationAccess &other) { - this->FEEvaluationBase::operator=(other); + this->FEEvaluationBase::operator=(other); return *this; } @@ -4239,64 +4560,67 @@ FEEvaluationAccess /*-------------------- FEEvaluationAccess scalar ----------------------------*/ -template +template inline -FEEvaluationAccess +FEEvaluationAccess ::FEEvaluationAccess (const MatrixFree &data_in, - const unsigned int fe_no, + 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 unsigned int n_q_points, + const bool is_interior_face) : - FEEvaluationBase - (data_in, fe_no, quad_no_in, fe_degree, n_q_points) + FEEvaluationBase + (data_in, dof_no, first_selected_component, quad_no_in, fe_degree, n_q_points, + is_interior_face) {} -template +template template inline -FEEvaluationAccess +FEEvaluationAccess ::FEEvaluationAccess (const Mapping &mapping, const FiniteElement &fe, const Quadrature<1> &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBase *other) + const FEEvaluationBase *other) : - FEEvaluationBase (mapping, fe, quadrature, update_flags, - first_selected_component, other) + FEEvaluationBase (mapping, fe, quadrature, update_flags, + first_selected_component, other) {} -template +template inline -FEEvaluationAccess -::FEEvaluationAccess (const FEEvaluationAccess &other) +FEEvaluationAccess +::FEEvaluationAccess (const FEEvaluationAccess &other) : - FEEvaluationBase (other) + FEEvaluationBase (other) {} -template +template inline -FEEvaluationAccess & -FEEvaluationAccess -::operator= (const FEEvaluationAccess &other) +FEEvaluationAccess & +FEEvaluationAccess +::operator= (const FEEvaluationAccess &other) { - this->FEEvaluationBase::operator=(other); + this->FEEvaluationBase::operator=(other); return *this; } -template -inline +template +inline DEAL_II_ALWAYS_INLINE VectorizedArray -FEEvaluationAccess +FEEvaluationAccess ::get_dof_value (const unsigned int dof) const { AssertIndexRange (dof, this->data->dofs_per_component_on_cell); @@ -4305,38 +4629,49 @@ FEEvaluationAccess -template -inline +template +inline DEAL_II_ALWAYS_INLINE VectorizedArray -FEEvaluationAccess +FEEvaluationAccess ::get_value (const unsigned int q_point) const { Assert (this->values_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); return this->values_quad[0][q_point]; } -template -inline +template +inline DEAL_II_ALWAYS_INLINE +VectorizedArray +FEEvaluationAccess +::get_normal_derivative (const unsigned int q_point) const +{ + return BaseClass::get_normal_derivative(q_point)[0]; +} + + + +template +inline DEAL_II_ALWAYS_INLINE Tensor<1,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_gradient (const unsigned int q_point) const { - // could use the base class gradient, but that involves too many inefficient + // could use the base class gradient, but that involves too many expensive // initialization operations on tensors Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); Assert (this->jacobian != nullptr, ExcNotInitialized()); Tensor<1,dim,VectorizedArray > grad_out; - if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { for (unsigned int d=0; dgradients_quad[0][d][q_point] * @@ -4360,10 +4695,10 @@ FEEvaluationAccess -template +template inline Tensor<2,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_hessian (const unsigned int q_point) const { return BaseClass::get_hessian(q_point)[0]; @@ -4371,10 +4706,10 @@ FEEvaluationAccess -template +template inline Tensor<1,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_hessian_diagonal (const unsigned int q_point) const { return BaseClass::get_hessian_diagonal(q_point)[0]; @@ -4382,10 +4717,10 @@ FEEvaluationAccess -template +template inline VectorizedArray -FEEvaluationAccess +FEEvaluationAccess ::get_laplacian (const unsigned int q_point) const { return BaseClass::get_laplacian(q_point)[0]; @@ -4393,10 +4728,10 @@ FEEvaluationAccess -template +template inline -void -FEEvaluationAccess +void DEAL_II_ALWAYS_INLINE +FEEvaluationAccess ::submit_dof_value (const VectorizedArray val_in, const unsigned int dof) { @@ -4409,48 +4744,74 @@ FEEvaluationAccess -template +template inline -void -FEEvaluationAccess +void DEAL_II_ALWAYS_INLINE +FEEvaluationAccess ::submit_value (const VectorizedArray val_in, - const unsigned int q_point) + const unsigned int q_index) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_index, this->n_quadrature_points); + Assert (this->J_value != nullptr, ExcNotInitialized()); this->values_quad_submitted = true; #endif - if (this->cell_type == internal::MatrixFreeFunctions::general) + if (this->cell_type <= internal::MatrixFreeFunctions::affine) { - const VectorizedArray JxW = this->J_value[q_point]; - this->values_quad[0][q_point] = val_in * JxW; + const VectorizedArray JxW = this->J_value[0] * this->quadrature_weights[q_index]; + this->values_quad[0][q_index] = val_in * JxW; } else //if (this->cell_type < internal::MatrixFreeFunctions::general) { - const VectorizedArray JxW = this->J_value[0] * this->quadrature_weights[q_point]; - this->values_quad[0][q_point] = val_in * JxW; + this->values_quad[0][q_index] = val_in * this->J_value[q_index]; } } -template -inline +template +inline DEAL_II_ALWAYS_INLINE +void +FEEvaluationAccess +::submit_value (const Tensor<1,1,VectorizedArray > val_in, + const unsigned int q_point) +{ + submit_value(val_in[0], q_point); +} + + + +template +inline DEAL_II_ALWAYS_INLINE +void +FEEvaluationAccess +::submit_normal_derivative (const VectorizedArray grad_in, + const unsigned int q_point) +{ + Tensor<1,1,VectorizedArray > grad; + grad[0] = grad_in; + BaseClass::submit_normal_derivative(grad, q_point); +} + + + +template +inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess +FEEvaluationAccess ::submit_gradient (const Tensor<1,dim,VectorizedArray > grad_in, const unsigned int q_index) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_index, this->data->n_q_points); + AssertIndexRange (q_index, this->n_quadrature_points); this->gradients_quad_submitted = true; Assert (this->J_value != nullptr, ExcNotInitialized()); Assert (this->jacobian != nullptr, ExcNotInitialized()); #endif - if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { const VectorizedArray JxW = this->J_value[0] * this->quadrature_weights[q_index]; @@ -4480,10 +4841,10 @@ FEEvaluationAccess -template +template inline VectorizedArray -FEEvaluationAccess +FEEvaluationAccess ::integrate_value () const { return BaseClass::integrate_value()[0]; @@ -4495,64 +4856,67 @@ FEEvaluationAccess /*----------------- FEEvaluationAccess vector-valued ------------------------*/ -template +template inline -FEEvaluationAccess +FEEvaluationAccess ::FEEvaluationAccess (const MatrixFree &data_in, - const unsigned int fe_no, + 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 unsigned int n_q_points, + const bool is_interior_face) : - FEEvaluationBase - (data_in, fe_no, quad_no_in, fe_degree, n_q_points) + FEEvaluationBase + (data_in, dof_no, first_selected_component, quad_no_in, fe_degree, n_q_points, + is_interior_face) {} -template +template template inline -FEEvaluationAccess +FEEvaluationAccess ::FEEvaluationAccess (const Mapping &mapping, const FiniteElement &fe, const Quadrature<1> &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBase *other) + const FEEvaluationBase *other) : - FEEvaluationBase (mapping, fe, quadrature, update_flags, - first_selected_component, other) + FEEvaluationBase (mapping, fe, quadrature, update_flags, + first_selected_component, other) {} -template +template inline -FEEvaluationAccess -::FEEvaluationAccess (const FEEvaluationAccess &other) +FEEvaluationAccess +::FEEvaluationAccess (const FEEvaluationAccess &other) : - FEEvaluationBase (other) + FEEvaluationBase (other) {} -template +template inline -FEEvaluationAccess & -FEEvaluationAccess -::operator= (const FEEvaluationAccess &other) +FEEvaluationAccess & +FEEvaluationAccess +::operator= (const FEEvaluationAccess &other) { - this->FEEvaluationAccess::operator=(other); + this->FEEvaluationBase::operator=(other); return *this; } -template -inline +template +inline DEAL_II_ALWAYS_INLINE Tensor<2,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_gradient (const unsigned int q_point) const { return BaseClass::get_gradient (q_point); @@ -4560,21 +4924,21 @@ FEEvaluationAccess -template -inline +template +inline DEAL_II_ALWAYS_INLINE VectorizedArray -FEEvaluationAccess +FEEvaluationAccess ::get_divergence (const unsigned int q_point) const { Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); Assert (this->jacobian != nullptr, ExcNotInitialized()); VectorizedArray divergence; // Cartesian cell - if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { divergence = (this->gradients_quad[0][0][q_point] * this->jacobian[0][0][0]); @@ -4600,10 +4964,10 @@ FEEvaluationAccess -template -inline +template +inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_symmetric_gradient (const unsigned int q_point) const { // copy from generic function into dim-specialization function @@ -4636,10 +5000,10 @@ FEEvaluationAccess -template -inline +template +inline DEAL_II_ALWAYS_INLINE Tensor<1,(dim==2?1:dim),VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_curl (const unsigned int q_point) const { // copy from generic function into dim-specialization function @@ -4667,39 +5031,35 @@ FEEvaluationAccess -template -inline +template +inline DEAL_II_ALWAYS_INLINE Tensor<2,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_hessian_diagonal (const unsigned int q_point) const { - Assert (this->hessians_quad_initialized==true, - internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); - return BaseClass::get_hessian_diagonal (q_point); } -template -inline +template +inline DEAL_II_ALWAYS_INLINE Tensor<3,dim,VectorizedArray > -FEEvaluationAccess +FEEvaluationAccess ::get_hessian (const unsigned int q_point) const { Assert (this->hessians_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); return BaseClass::get_hessian(q_point); } -template -inline +template +inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess +FEEvaluationAccess ::submit_gradient (const Tensor<2,dim,VectorizedArray > grad_in, const unsigned int q_point) { @@ -4708,10 +5068,10 @@ FEEvaluationAccess -template -inline +template +inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess +FEEvaluationAccess ::submit_gradient (const Tensor<1,dim,Tensor<1,dim,VectorizedArray > > grad_in, const unsigned int q_point) @@ -4721,22 +5081,22 @@ FEEvaluationAccess -template -inline +template +inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess +FEEvaluationAccess ::submit_divergence (const VectorizedArray div_in, const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); this->gradients_quad_submitted = true; Assert (this->J_value != nullptr, ExcNotInitialized()); Assert (this->jacobian != nullptr, ExcNotInitialized()); #endif - if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { const VectorizedArray fac = this->J_value[0] * this->quadrature_weights[q_point] * div_in; @@ -4770,10 +5130,10 @@ FEEvaluationAccess -template -inline +template +inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess +FEEvaluationAccess ::submit_symmetric_gradient(const SymmetricTensor<2,dim,VectorizedArray > sym_grad, const unsigned int q_point) @@ -4783,13 +5143,13 @@ 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->n_quadrature_points); this->gradients_quad_submitted = true; Assert (this->J_value != nullptr, ExcNotInitialized()); Assert (this->jacobian != nullptr, ExcNotInitialized()); #endif - if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { const VectorizedArray JxW = this->J_value[0] * this->quadrature_weights[q_point]; for (unsigned int d=0; d -template -inline +template +inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess +FEEvaluationAccess ::submit_curl (const Tensor<1,dim==2?1:dim,VectorizedArray > curl, const unsigned int q_point) { @@ -4874,64 +5234,67 @@ FEEvaluationAccess /*-------------------- FEEvaluationAccess scalar for 1d ----------------------------*/ -template +template inline -FEEvaluationAccess<1,1,Number> +FEEvaluationAccess<1,1,Number,is_face> ::FEEvaluationAccess (const MatrixFree<1,Number> &data_in, - const unsigned int fe_no, + 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 unsigned int n_q_points, + const bool is_interior_face) : - FEEvaluationBase <1,1,Number> - (data_in, fe_no, quad_no_in, fe_degree, n_q_points) + FEEvaluationBase <1,1,Number,is_face> + (data_in, dof_no, first_selected_component, quad_no_in, fe_degree, n_q_points, + is_interior_face) {} -template +template template inline -FEEvaluationAccess<1,1,Number> +FEEvaluationAccess<1,1,Number,is_face> ::FEEvaluationAccess (const Mapping<1> &mapping, const FiniteElement<1> &fe, const Quadrature<1> &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, - const FEEvaluationBase<1,n_components_other,Number> *other) + const FEEvaluationBase<1,n_components_other,Number,is_face> *other) : - FEEvaluationBase <1,1,Number> (mapping, fe, quadrature, update_flags, - first_selected_component, other) + FEEvaluationBase <1,1,Number,is_face> (mapping, fe, quadrature, update_flags, + first_selected_component, other) {} -template +template inline -FEEvaluationAccess<1,1,Number> -::FEEvaluationAccess (const FEEvaluationAccess<1,1,Number> &other) +FEEvaluationAccess<1,1,Number,is_face> +::FEEvaluationAccess (const FEEvaluationAccess<1,1,Number,is_face> &other) : - FEEvaluationBase <1,1,Number>(other) + FEEvaluationBase <1,1,Number,is_face>(other) {} -template +template inline -FEEvaluationAccess<1,1,Number> & -FEEvaluationAccess<1,1,Number> -::operator= (const FEEvaluationAccess<1,1,Number> &other) +FEEvaluationAccess<1,1,Number,is_face> & +FEEvaluationAccess<1,1,Number,is_face> +::operator= (const FEEvaluationAccess<1,1,Number,is_face> &other) { - this->FEEvaluationBase<1,1,Number>::operator=(other); + this->FEEvaluationBase<1,1,Number,is_face>::operator=(other); return *this; } -template -inline +template +inline DEAL_II_ALWAYS_INLINE VectorizedArray -FEEvaluationAccess<1,1,Number> +FEEvaluationAccess<1,1,Number,is_face> ::get_dof_value (const unsigned int dof) const { AssertIndexRange (dof, this->data->dofs_per_component_on_cell); @@ -4940,24 +5303,24 @@ FEEvaluationAccess<1,1,Number> -template -inline +template +inline DEAL_II_ALWAYS_INLINE VectorizedArray -FEEvaluationAccess<1,1,Number> +FEEvaluationAccess<1,1,Number,is_face> ::get_value (const unsigned int q_point) const { Assert (this->values_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); return this->values_quad[0][q_point]; } -template -inline +template +inline DEAL_II_ALWAYS_INLINE Tensor<1,1,VectorizedArray > -FEEvaluationAccess<1,1,Number> +FEEvaluationAccess<1,1,Number,is_face> ::get_gradient (const unsigned int q_point) const { // could use the base class gradient, but that involves too many inefficient @@ -4965,7 +5328,7 @@ FEEvaluationAccess<1,1,Number> Assert (this->gradients_quad_initialized==true, internal::ExcAccessToUninitializedField()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); const Tensor<2,1,VectorizedArray > &jac = this->cell_type == internal::MatrixFreeFunctions::general ? @@ -4979,10 +5342,21 @@ FEEvaluationAccess<1,1,Number> -template -inline +template +inline DEAL_II_ALWAYS_INLINE +VectorizedArray +FEEvaluationAccess<1,1,Number,is_face> +::get_normal_derivative (const unsigned int q_point) const +{ + return BaseClass::get_normal_derivative(q_point)[0]; +} + + + +template +inline DEAL_II_ALWAYS_INLINE Tensor<2,1,VectorizedArray > -FEEvaluationAccess<1,1,Number> +FEEvaluationAccess<1,1,Number,is_face> ::get_hessian (const unsigned int q_point) const { return BaseClass::get_hessian(q_point)[0]; @@ -4990,10 +5364,10 @@ FEEvaluationAccess<1,1,Number> -template -inline +template +inline DEAL_II_ALWAYS_INLINE Tensor<1,1,VectorizedArray > -FEEvaluationAccess<1,1,Number> +FEEvaluationAccess<1,1,Number,is_face> ::get_hessian_diagonal (const unsigned int q_point) const { return BaseClass::get_hessian_diagonal(q_point)[0]; @@ -5001,10 +5375,10 @@ FEEvaluationAccess<1,1,Number> -template -inline +template +inline DEAL_II_ALWAYS_INLINE VectorizedArray -FEEvaluationAccess<1,1,Number> +FEEvaluationAccess<1,1,Number,is_face> ::get_laplacian (const unsigned int q_point) const { return BaseClass::get_laplacian(q_point)[0]; @@ -5012,10 +5386,10 @@ FEEvaluationAccess<1,1,Number> -template -inline -void -FEEvaluationAccess<1,1,Number> +template +inline DEAL_II_ALWAYS_INLINE +void DEAL_II_ALWAYS_INLINE +FEEvaluationAccess<1,1,Number,is_face> ::submit_dof_value (const VectorizedArray val_in, const unsigned int dof) { @@ -5028,16 +5402,16 @@ FEEvaluationAccess<1,1,Number> -template -inline +template +inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess<1,1,Number> +FEEvaluationAccess<1,1,Number,is_face> ::submit_value (const VectorizedArray val_in, const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); this->values_quad_submitted = true; #endif if (this->cell_type == internal::MatrixFreeFunctions::general) @@ -5045,7 +5419,7 @@ FEEvaluationAccess<1,1,Number> const VectorizedArray JxW = this->J_value[q_point]; this->values_quad[0][q_point] = val_in * JxW; } - else //if (this->cell_type < internal::MatrixFreeFunctions::general) + else //if (this->cell_type == internal::MatrixFreeFunctions::general) { const VectorizedArray JxW = this->J_value[0] * this->quadrature_weights[q_point]; this->values_quad[0][q_point] = val_in * JxW; @@ -5054,35 +5428,85 @@ FEEvaluationAccess<1,1,Number> -template -inline +template +inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess<1,1,Number> +FEEvaluationAccess<1,1,Number,is_face> +::submit_value (const Tensor<1,1,VectorizedArray > val_in, + const unsigned int q_point) +{ + submit_value(val_in[0], q_point); +} + + + +template +inline DEAL_II_ALWAYS_INLINE +void +FEEvaluationAccess<1,1,Number,is_face> ::submit_gradient (const Tensor<1,1,VectorizedArray > grad_in, const unsigned int q_point) +{ + submit_gradient(grad_in[0], q_point); +} + + + +template +inline DEAL_II_ALWAYS_INLINE +void +FEEvaluationAccess<1,1,Number,is_face> +::submit_gradient (const VectorizedArray grad_in, + const unsigned int q_point) { #ifdef DEBUG Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized()); - AssertIndexRange (q_point, this->data->n_q_points); + AssertIndexRange (q_point, this->n_quadrature_points); this->gradients_quad_submitted = true; #endif const Tensor<2,1,VectorizedArray > &jac = - this->cell_type > internal::MatrixFreeFunctions::affine ? + this->cell_type == internal::MatrixFreeFunctions::general ? this->jacobian[q_point] : this->jacobian[0]; const VectorizedArray JxW = - this->cell_type > internal::MatrixFreeFunctions::affine ? + this->cell_type == internal::MatrixFreeFunctions::general ? this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point]; - this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in[0] * JxW; + this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in * JxW; +} + + + +template +inline DEAL_II_ALWAYS_INLINE +void +FEEvaluationAccess<1,1,Number,is_face> +::submit_normal_derivative (const VectorizedArray grad_in, + const unsigned int q_point) +{ + Tensor<1,1,VectorizedArray > grad; + grad[0] = grad_in; + BaseClass::submit_normal_derivative(grad, q_point); } -template +template +inline DEAL_II_ALWAYS_INLINE +void +FEEvaluationAccess<1,1,Number,is_face> +::submit_normal_derivative (const Tensor<1,1,VectorizedArray > grad_in, + const unsigned int q_point) +{ + BaseClass::submit_normal_derivative(grad_in, q_point); +} + + + +template inline VectorizedArray -FEEvaluationAccess<1,1,Number> +FEEvaluationAccess<1,1,Number,is_face> ::integrate_value () const { return BaseClass::integrate_value()[0]; @@ -5100,9 +5524,10 @@ inline FEEvaluation ::FEEvaluation (const MatrixFree &data_in, const unsigned int fe_no, - const unsigned int quad_no) + const unsigned int quad_no, + const unsigned int first_selected_component) : - BaseClass (data_in, fe_no, quad_no, fe_degree, static_n_q_points), + BaseClass (data_in, fe_no, first_selected_component, quad_no, fe_degree, static_n_q_points), 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) @@ -5124,7 +5549,7 @@ FEEvaluation : BaseClass (mapping, fe, quadrature, update_flags, first_selected_component, - static_cast*>(nullptr)), + static_cast*>(nullptr)), 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) @@ -5145,7 +5570,7 @@ FEEvaluation : BaseClass (StaticMappingQ1::mapping, fe, quadrature, update_flags, first_selected_component, - static_cast*>(nullptr)), + static_cast*>(nullptr)), 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) @@ -5212,18 +5637,19 @@ template -::check_template_arguments(const unsigned int fe_no, +::check_template_arguments(const unsigned int dof_no, const unsigned int first_selected_component) { - (void)fe_no; + (void)dof_no; (void)first_selected_component; #ifdef DEBUG // print error message when the dimensions do not match. Propose a possible // fix - if ((fe_degree != -1 && static_cast(fe_degree) != this->data->fe_degree) + if ((static_cast(fe_degree) != numbers::invalid_unsigned_int && + static_cast(fe_degree) != this->data->fe_degree) || - (fe_degree != -1 && static_n_q_points != this->data->n_q_points)) + n_q_points != this->n_quadrature_points) { std::string message = "-------------------------------------------------------\n"; @@ -5233,27 +5659,33 @@ FEEvaluation message += Utilities::int_to_string(n_q_points_1d); message += "," + Utilities::int_to_string(n_components); message += ",Number>(data"; - if (fe_no != numbers::invalid_unsigned_int) + if (first_selected_component != numbers::invalid_unsigned_int) { - message += ", " + Utilities::int_to_string(fe_no) + ", "; - message += Utilities::int_to_string(this->quad_no); + message += ", " + Utilities::int_to_string(dof_no) + ", "; + message += Utilities::int_to_string(this->quad_no) + ", "; + message += Utilities::int_to_string(first_selected_component); } message += ")\n"; // check whether some other vector component has the correct number of // points - unsigned int proposed_dof_comp = numbers::invalid_unsigned_int, + unsigned int proposed_dof_comp = numbers::invalid_unsigned_int, + proposed_fe_comp = numbers::invalid_unsigned_int, proposed_quad_comp = numbers::invalid_unsigned_int; - if (fe_no != numbers::invalid_unsigned_int) + if (dof_no != numbers::invalid_unsigned_int) { if (static_cast(fe_degree) == this->data->fe_degree) - proposed_dof_comp = fe_no; + { + proposed_dof_comp = dof_no; + proposed_fe_comp = first_selected_component; + } else for (unsigned int no=0; nomatrix_info->n_components(); ++no) if (this->matrix_info->get_shape_info(no,0,this->active_fe_index,0).fe_degree == static_cast(fe_degree)) { proposed_dof_comp = no; + proposed_fe_comp = 0; break; } if (n_q_points == @@ -5271,7 +5703,7 @@ FEEvaluation if (proposed_dof_comp != numbers::invalid_unsigned_int && proposed_quad_comp != numbers::invalid_unsigned_int) { - if (proposed_dof_comp != fe_no) + if (proposed_dof_comp != first_selected_component) message += "Wrong vector component selection:\n"; else message += "Wrong quadrature formula selection:\n"; @@ -5280,18 +5712,23 @@ FEEvaluation message += Utilities::int_to_string(n_q_points_1d); message += "," + Utilities::int_to_string(n_components); message += ",Number>(data"; - if (fe_no != numbers::invalid_unsigned_int) + if (dof_no != numbers::invalid_unsigned_int) { message += ", " + Utilities::int_to_string(proposed_dof_comp) + ", "; - message += Utilities::int_to_string(proposed_quad_comp); + message += Utilities::int_to_string(proposed_quad_comp) + ", "; + message += Utilities::int_to_string(proposed_fe_comp); } message += ")?\n"; std::string correct_pos; - if (proposed_dof_comp != fe_no) + if (proposed_dof_comp != dof_no) correct_pos = " ^ "; else correct_pos = " "; if (proposed_quad_comp != this->quad_no) + correct_pos += " ^ "; + else + correct_pos += " "; + if (proposed_fe_comp != first_selected_component) correct_pos += " ^\n"; else correct_pos += " \n"; @@ -5299,17 +5736,18 @@ FEEvaluation } // ok, did not find the numbers specified by the template arguments in // the given list. Suggest correct template arguments - 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_n_q_points_1d = static_cast(std::pow(1.001*this->n_quadrature_points,1./dim)); message += "Wrong template arguments:\n"; message += " Did you mean FEEvaluationdata->fe_degree) + ","; message += Utilities::int_to_string(proposed_n_q_points_1d); message += "," + Utilities::int_to_string(n_components); message += ",Number>(data"; - if (fe_no != numbers::invalid_unsigned_int) + if (dof_no != numbers::invalid_unsigned_int) { - message += ", " + Utilities::int_to_string(fe_no) + ", "; + message += ", " + Utilities::int_to_string(dof_no) + ", "; message += Utilities::int_to_string(this->quad_no); + message += ", " + Utilities::int_to_string(first_selected_component); } message += ")?\n"; std::string correct_pos; @@ -5324,10 +5762,10 @@ FEEvaluation message += " " + correct_pos; Assert (static_cast(fe_degree) == this->data->fe_degree && - static_n_q_points == this->data->n_q_points, + n_q_points == this->n_quadrature_points, ExcMessage(message)); } - if (fe_no != numbers::invalid_unsigned_int) + if (dof_no != numbers::invalid_unsigned_int) { AssertDimension (n_q_points, this->mapping_data->descriptor[this->active_quad_index].n_q_points); @@ -5359,8 +5797,6 @@ FEEvaluation } AssertIndexRange (q, n_q_points); - const Point > *quadrature_points = &this->mapping_data-> - quadrature_points[this->mapping_data->quadrature_point_offsets[this->cell]]; const unsigned int n_q_points_1d_actual = fe_degree == -1 ? this->data->n_q_points_1d : n_q_points_1d; @@ -5368,6 +5804,8 @@ FEEvaluation // Cartesian mesh: not all quadrature points are stored, only the // diagonal. Hence, need to find the tensor product index and retrieve the // value from that + const Point > *quadrature_points = + &this->mapping_data->quadrature_points[this->mapping_data->quadrature_point_offsets[this->cell]]; if (this->cell_type == internal::MatrixFreeFunctions::cartesian) { Point > point; diff --git a/tests/matrix_free/fe_evaluation_access_1d.cc b/tests/matrix_free/fe_evaluation_access_1d.cc index 852c38d112..16dbc86c93 100644 --- a/tests/matrix_free/fe_evaluation_access_1d.cc +++ b/tests/matrix_free/fe_evaluation_access_1d.cc @@ -26,6 +26,6 @@ int main () { initlog(); - FEEvaluationAccess<1,1,double> *test; // didn't compile before + FEEvaluationAccess<1,1,double,false> *test; // didn't compile before deallog << "OK" << std::endl; } -- 2.39.5