From b3569a82d914372ed04083382233300ce201c213 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 21 Jul 2015 09:06:28 -0500 Subject: [PATCH] Split FEValuesData into two base classes. There is now one base class that stores Mapping-related data, and one that stores FiniteElement-related data. These are declared in fe_update_flags() so that they can be referenced from the finite element and mapping classes without having to know about fe_values.h. --- include/deal.II/fe/fe_update_flags.h | 191 ++++++++++++++++++ include/deal.II/fe/fe_values.h | 170 +--------------- .../matrix_free/mapping_data_on_the_fly.h | 6 +- include/deal.II/matrix_free/mapping_info.h | 2 +- .../matrix_free/mapping_info.templates.h | 12 +- source/fe/fe_values.cc | 134 +++++++----- 6 files changed, 283 insertions(+), 232 deletions(-) diff --git a/include/deal.II/fe/fe_update_flags.h b/include/deal.II/fe/fe_update_flags.h index 0a934f287f..4622f309a3 100644 --- a/include/deal.II/fe/fe_update_flags.h +++ b/include/deal.II/fe/fe_update_flags.h @@ -18,9 +18,19 @@ #include +#include +#include +#include +#include + +#include + DEAL_II_NAMESPACE_OPEN +template class FiniteElement; + + /*!@addtogroup feaccess */ /*@{*/ @@ -336,6 +346,187 @@ namespace CellSimilarity } +namespace internal +{ + namespace FEValues + { + /** + * A class that stores all of the mapping related data used in + * dealii::FEValues, dealii::FEFaceValues, and dealii::FESubfaceValues + * objects. Objects of this kind will be given + * as output argument when dealii::FEValues::reinit() + * calls Mapping::fill_fe_values() for a given cell, face, or subface. + * + * The data herein will then be provided as input argument in + * the following call to FiniteElement::fill_fe_values(). + * + * @ingroup feaccess + */ + template + class MappingRelatedData + { + public: + /** + * Initialize all vectors to correct size. + */ + void initialize (const unsigned int n_quadrature_points, + const UpdateFlags flags); + + /** + * Store an array of weights times the Jacobi determinant at the quadrature + * points. This function is reset each time reinit() is called. The Jacobi + * determinant is actually the reciprocal value of the Jacobi matrices + * stored in this class, see the general documentation of this class for + * more information. + * + * However, if this object refers to an FEFaceValues or FESubfaceValues + * object, then the JxW_values correspond to the Jacobian of the + * transformation of the face, not the cell, i.e. the dimensionality is that + * of a surface measure, not of a volume measure. In this case, it is + * computed from the boundary forms, rather than the Jacobian matrix. + */ + std::vector JxW_values; + + /** + * Array of the Jacobian matrices at the quadrature points. + */ + std::vector< DerivativeForm<1,dim,spacedim> > jacobians; + + /** + * Array of the derivatives of the Jacobian matrices at the quadrature + * points. + */ + std::vector > jacobian_grads; + + /** + * Array of the inverse Jacobian matrices at the quadrature points. + */ + std::vector > inverse_jacobians; + + /** + * Array of quadrature points. This array is set up upon calling reinit() + * and contains the quadrature points on the real element, rather than on + * the reference element. + */ + std::vector > quadrature_points; + + /** + * List of outward normal vectors at the quadrature points. This field is + * filled in by the finite element class. + */ + std::vector > normal_vectors; + + /** + * List of boundary forms at the quadrature points. This field is filled in + * by the finite element class. + */ + std::vector > boundary_forms; + }; + + + /** + * A class that stores all of the shape function related data used in + * dealii::FEValues, dealii::FEFaceValues, and dealii::FESubfaceValues + * objects. Objects of this kind will be given + * as output argument when dealii::FEValues::reinit() + * calls FiniteElement::fill_fe_values(). + * + * @ingroup feaccess + */ + template + class FiniteElementRelatedData + { + public: + /** + * Initialize all vectors to correct size. + */ + void initialize (const unsigned int n_quadrature_points, + const FiniteElement &fe, + const UpdateFlags flags); + + /** + * Storage type for shape values. Each row in the matrix denotes the values + * of a single shape function at the different points, columns are for a + * single point with the different shape functions. + * + * If a shape function has more than one non-zero component (in deal.II + * diction: it is non-primitive), then we allocate one row per non-zero + * component, and shift subsequent rows backward. Lookup of the correct row + * for a shape function is thus simple in case the entire finite element is + * primitive (i.e. all shape functions are primitive), since then the shape + * function number equals the row number. Otherwise, use the + * #shape_function_to_row_table array to get at the first row that belongs + * to this particular shape function, and navigate among all the rows for + * this shape function using the FiniteElement::get_nonzero_components() + * function which tells us which components are non-zero and thus have a row + * in the array presently under discussion. + */ + typedef dealii::Table<2,double> ShapeVector; + + /** + * Storage type for gradients. The layout of data is the same as for the + * #ShapeVector data type. + */ + typedef std::vector > > GradientVector; + + /** + * Likewise for second order derivatives. + */ + typedef std::vector > > HessianVector; + + /** + * Store the values of the shape functions at the quadrature points. See the + * description of the data type for the layout of the data in this field. + */ + ShapeVector shape_values; + + /** + * Store the gradients of the shape functions at the quadrature points. See + * the description of the data type for the layout of the data in this + * field. + */ + GradientVector shape_gradients; + + /** + * Store the 2nd derivatives of the shape functions at the quadrature + * points. See the description of the data type for the layout of the data + * in this field. + */ + HessianVector shape_hessians; + + /** + * When asked for the value (or gradient, or Hessian) of shape function i's + * c-th vector component, we need to look it up in the #shape_values, + * #shape_gradients and #shape_hessians arrays. The question is where in + * this array does the data for shape function i, component c reside. This + * is what this table answers. + * + * The format of the table is as follows: - It has dofs_per_cell times + * n_components entries. - The entry that corresponds to shape function i, + * component c is i * n_components + c. - The value stored at + * this position indicates the row in #shape_values and the other tables + * where the corresponding datum is stored for all the quadrature points. + * + * In the general, vector-valued context, the number of components is larger + * than one, but for a given shape function, not all vector components may + * be nonzero (e.g., if a shape function is primitive, then exactly one + * vector component is non-zero, while the others are all zero). For such + * zero components, #shape_values and friends do not have a row. + * Consequently, for vector components for which shape function i is zero, + * the entry in the current table is numbers::invalid_unsigned_int. + * + * On the other hand, the table is guaranteed to have at least one valid + * index for each shape function. In particular, for a primitive finite + * element, each shape function has exactly one nonzero component and so for + * each i, there is exactly one valid index within the range + * [i*n_components, (i+1)*n_components). + */ + std::vector shape_function_to_row_table; + }; + } +} + + /*@}*/ diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h index dd8afbea03..f65a767a3a 100644 --- a/include/deal.II/fe/fe_values.h +++ b/include/deal.II/fe/fe_values.h @@ -1243,38 +1243,19 @@ namespace internal //TODO: Add access to mapping values to FEValuesBase -//TODO: Several FEValuesBase of a system should share Mapping /** * A class that contains all data vectors for FEValues, FEFaceValues, and * FESubfaceValues. * - * This class has been extracted from FEValuesBase to encapsulate in one - * place all of the data, independent of the functions thater later - * access this data in the public interfaces of the FEValues and related - * classes. Consequently, this base class is protected in FEValuesBase. - * - * The second reason is because in FEValuesBase::reinit, we first need to - * call Mapping::fill_fe_values() to compute mapping related data, and later - * call FiniteElement::fill_fe_values() to compute shape function related - * data. In the first step, Mapping::fill_fe_values() gets a pointer to - * its own internal data structure and a pointer to the FEValuesData base - * object of FEValuesBase, and the mapping then places the computed data - * into the data fields that pertain to the mapping below. In the second - * step, the finite element receives a pointer to its own internal object, - * and to the current object, and from both of these computes the shape - * function related information and, again, places it into the current - * FEValuesData object. - * * More information can be found on the page on * @ref UpdateFlagsEssay. * * @ingroup feaccess - * @author Guido Kanschat - * @date 2000 */ template -class FEValuesData +class FEValuesData : public internal::FEValues::MappingRelatedData, + public internal::FEValues::FiniteElementRelatedData { public: /** @@ -1284,153 +1265,6 @@ public: const FiniteElement &fe, const UpdateFlags flags); - /** - * @name Fields filled by the finite element - * @{ - */ - - /** - * Storage type for shape values. Each row in the matrix denotes the values - * of a single shape function at the different points, columns are for a - * single point with the different shape functions. - * - * If a shape function has more than one non-zero component (in deal.II - * diction: it is non-primitive), then we allocate one row per non-zero - * component, and shift subsequent rows backward. Lookup of the correct row - * for a shape function is thus simple in case the entire finite element is - * primitive (i.e. all shape functions are primitive), since then the shape - * function number equals the row number. Otherwise, use the - * #shape_function_to_row_table array to get at the first row that belongs - * to this particular shape function, and navigate among all the rows for - * this shape function using the FiniteElement::get_nonzero_components() - * function which tells us which components are non-zero and thus have a row - * in the array presently under discussion. - */ - typedef Table<2,double> ShapeVector; - - /** - * Storage type for gradients. The layout of data is the same as for the - * #ShapeVector data type. - */ - typedef std::vector > > GradientVector; - - /** - * Likewise for second order derivatives. - */ - typedef std::vector > > HessianVector; - - /** - * Store the values of the shape functions at the quadrature points. See the - * description of the data type for the layout of the data in this field. - */ - ShapeVector shape_values; - - /** - * Store the gradients of the shape functions at the quadrature points. See - * the description of the data type for the layout of the data in this - * field. - */ - GradientVector shape_gradients; - - /** - * Store the 2nd derivatives of the shape functions at the quadrature - * points. See the description of the data type for the layout of the data - * in this field. - */ - HessianVector shape_hessians; - - /** - * @} - */ - - /** - * @name Fields filled by the mapping - * @{ - */ - - /** - * Store an array of weights times the Jacobi determinant at the quadrature - * points. This function is reset each time reinit() is called. The Jacobi - * determinant is actually the reciprocal value of the Jacobi matrices - * stored in this class, see the general documentation of this class for - * more information. - * - * However, if this object refers to an FEFaceValues or FESubfaceValues - * object, then the JxW_values correspond to the Jacobian of the - * transformation of the face, not the cell, i.e. the dimensionality is that - * of a surface measure, not of a volume measure. In this case, it is - * computed from the boundary forms, rather than the Jacobian matrix. - */ - std::vector JxW_values; - - /** - * Array of the Jacobian matrices at the quadrature points. - */ - std::vector< DerivativeForm<1,dim,spacedim> > jacobians; - - /** - * Array of the derivatives of the Jacobian matrices at the quadrature - * points. - */ - std::vector > jacobian_grads; - - /** - * Array of the inverse Jacobian matrices at the quadrature points. - */ - std::vector > inverse_jacobians; - - /** - * Array of quadrature points. This array is set up upon calling reinit() - * and contains the quadrature points on the real element, rather than on - * the reference element. - */ - std::vector > quadrature_points; - - /** - * List of outward normal vectors at the quadrature points. This field is - * filled in by the finite element class. - */ - std::vector > normal_vectors; - - /** - * List of boundary forms at the quadrature points. This field is filled in - * by the finite element class. - */ - std::vector > boundary_forms; - - /** - * @} - */ - - /** - * When asked for the value (or gradient, or Hessian) of shape function i's - * c-th vector component, we need to look it up in the #shape_values, - * #shape_gradients and #shape_hessians arrays. The question is where in - * this array does the data for shape function i, component c reside. This - * is what this table answers. - * - * The format of the table is as follows: - It has dofs_per_cell times - * n_components entries. - The entry that corresponds to shape function i, - * component c is i * n_components + c. - The value stored at - * this position indicates the row in #shape_values and the other tables - * where the corresponding datum is stored for all the quadrature points. - * - * In the general, vector-valued context, the number of components is larger - * than one, but for a given shape function, not all vector components may - * be nonzero (e.g., if a shape function is primitive, then exactly one - * vector component is non-zero, while the others are all zero). For such - * zero components, #shape_values and friends do not have a row. - * Consequently, for vector components for which shape function i is zero, - * the entry in the current table is numbers::invalid_unsigned_int. - * - * On the other hand, the table is guaranteed to have at least one valid - * index for each shape function. In particular, for a primitive finite - * element, each shape function has exactly one nonzero component and so for - * each i, there is exactly one valid index within the range - * [i*n_components, (i+1)*n_components). - */ - std::vector shape_function_to_row_table; - /** * Original update flags handed to the constructor of FEValues. */ diff --git a/include/deal.II/matrix_free/mapping_data_on_the_fly.h b/include/deal.II/matrix_free/mapping_data_on_the_fly.h index 1f812b4170..12b187c14d 100644 --- a/include/deal.II/matrix_free/mapping_data_on_the_fly.h +++ b/include/deal.II/matrix_free/mapping_data_on_the_fly.h @@ -96,7 +96,7 @@ namespace internal * mapped quadrature points are accessible, as no finite element data is * actually used). */ - const FEValues &get_fe_values () const; + const dealii::FEValues &get_fe_values () const; /** * Return a vector of inverse transpose Jacobians. For compatibility @@ -152,7 +152,7 @@ namespace internal /** * An underlying FEValues object that performs the (scalar) evaluation. */ - FEValues fe_values; + dealii::FEValues fe_values; /** * Get 1D quadrature formula to be used for reinitializing shape info. @@ -272,7 +272,7 @@ namespace internal template inline - const FEValues & + const dealii::FEValues & MappingDataOnTheFly::get_fe_values() const { return fe_values; diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h index a5293eab9d..906e2ae6ad 100644 --- a/include/deal.II/matrix_free/mapping_info.h +++ b/include/deal.II/matrix_free/mapping_info.h @@ -328,7 +328,7 @@ namespace internal const unsigned int my_q, CellType (&cell_t_prev)[n_vector_elements], CellType (&cell_t)[n_vector_elements], - FEValues &fe_values, + dealii::FEValues &fe_values, CellData &cell_data) const; }; diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index d2f74ad086..1c68ab104e 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -240,7 +240,7 @@ namespace internal // hp::DoFHandler::active_cell_iterator, we need to manually // select the correct finite element, so just hold a vector of // FEValues - std::vector > > + std::vector > > fe_values (current_data.quadrature.size()); UpdateFlags update_flags_feval = (update_flags & update_inverse_jacobians ? update_jacobians : update_default) | @@ -290,10 +290,10 @@ namespace internal const unsigned int n_q_points = current_data.n_q_points[fe_index]; if (fe_values[fe_index].get() == 0) fe_values[fe_index].reset - (new FEValues (mapping, dummy_fe, - current_data.quadrature[fe_index], - update_flags_feval)); - FEValues &fe_val = *fe_values[fe_index]; + (new dealii::FEValues (mapping, dummy_fe, + current_data.quadrature[fe_index], + update_flags_feval)); + dealii::FEValues &fe_val = *fe_values[fe_index]; data.resize (n_q_points); // if the fe index has changed from the previous cell, set the @@ -608,7 +608,7 @@ namespace internal const unsigned int my_q, CellType (&cell_t_prev)[n_vector_elements], CellType (&cell_t)[n_vector_elements], - FEValues &fe_val, + dealii::FEValues &fe_val, CellData &data) const { const unsigned int n_q_points = fe_val.n_quadrature_points; diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index 909f6d11d3..49792db89b 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -2073,67 +2073,93 @@ get_interpolated_dof_values (const IndexSet &, /* --------------------- FEValuesData ----------------- */ +namespace internal +{ + namespace FEValues + { + template + void + MappingRelatedData::initialize (const unsigned int n_quadrature_points, + const UpdateFlags flags) + { + if (flags & update_quadrature_points) + this->quadrature_points.resize(n_quadrature_points); + + if (flags & update_JxW_values) + this->JxW_values.resize(n_quadrature_points); + + if (flags & update_jacobians) + this->jacobians.resize(n_quadrature_points); + + if (flags & update_jacobian_grads) + this->jacobian_grads.resize(n_quadrature_points); + + if (flags & update_inverse_jacobians) + this->inverse_jacobians.resize(n_quadrature_points); + + if (flags & update_boundary_forms) + this->boundary_forms.resize(n_quadrature_points); + + if (flags & update_normal_vectors) + this->normal_vectors.resize(n_quadrature_points); + } + + + template + void + FiniteElementRelatedData::initialize (const unsigned int n_quadrature_points, + const FiniteElement &fe, + const UpdateFlags flags) + { + + // initialize the table mapping + // from shape function number to + // the rows in the tables storing + // the data by shape function and + // nonzero component + this->shape_function_to_row_table + = make_shape_function_to_row_table (fe); + + // count the total number of non-zero + // components accumulated over all shape + // functions + unsigned int n_nonzero_shape_components = 0; + for (unsigned int i=0; i= fe.dofs_per_cell, + ExcInternalError()); + + // with the number of rows now + // known, initialize those fields + // that we will need to their + // correct size + if (flags & update_values) + this->shape_values.reinit(n_nonzero_shape_components, + n_quadrature_points); + + if (flags & update_gradients) + this->shape_gradients.resize (n_nonzero_shape_components, + std::vector > (n_quadrature_points)); + + if (flags & update_hessians) + this->shape_hessians.resize (n_nonzero_shape_components, + std::vector > (n_quadrature_points)); + } + } +} + + template void FEValuesData::initialize (const unsigned int n_quadrature_points, const FiniteElement &fe, const UpdateFlags flags) { - this->update_flags = flags; + // initialize the base classes + internal::FEValues::MappingRelatedData::initialize(n_quadrature_points, flags); + internal::FEValues::FiniteElementRelatedData::initialize(n_quadrature_points, fe, flags); - // initialize the table mapping - // from shape function number to - // the rows in the tables storing - // the data by shape function and - // nonzero component - this->shape_function_to_row_table - = make_shape_function_to_row_table (fe); - - // count the total number of non-zero - // components accumulated over all shape - // functions - unsigned int n_nonzero_shape_components = 0; - for (unsigned int i=0; i= fe.dofs_per_cell, - ExcInternalError()); - - // with the number of rows now - // known, initialize those fields - // that we will need to their - // correct size - if (flags & update_values) - this->shape_values.reinit(n_nonzero_shape_components, - n_quadrature_points); - - if (flags & update_gradients) - this->shape_gradients.resize (n_nonzero_shape_components, - std::vector > (n_quadrature_points)); - - if (flags & update_hessians) - this->shape_hessians.resize (n_nonzero_shape_components, - std::vector > (n_quadrature_points)); - - if (flags & update_quadrature_points) - this->quadrature_points.resize(n_quadrature_points); - - if (flags & update_JxW_values) - this->JxW_values.resize(n_quadrature_points); - - if (flags & update_jacobians) - this->jacobians.resize(n_quadrature_points); - - if (flags & update_jacobian_grads) - this->jacobian_grads.resize(n_quadrature_points); - - if (flags & update_inverse_jacobians) - this->inverse_jacobians.resize(n_quadrature_points); - - if (flags & update_boundary_forms) - this->boundary_forms.resize(n_quadrature_points); - - if (flags & update_normal_vectors) - this->normal_vectors.resize(n_quadrature_points); + this->update_flags = flags; } -- 2.39.5