From: Peter Munch Date: Tue, 28 Jul 2020 20:04:21 +0000 (+0200) Subject: Assign faces different quadrature rules X-Git-Tag: v9.3.0-rc1~741^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d96ef0c870a3a8fbcf02fd545567fd4fc0ef4d9a;p=dealii.git Assign faces different quadrature rules --- diff --git a/doc/news/changes/incompatibilities/20201212Munch b/doc/news/changes/incompatibilities/20201212Munch new file mode 100644 index 0000000000..4ecd10bb7d --- /dev/null +++ b/doc/news/changes/incompatibilities/20201212Munch @@ -0,0 +1,6 @@ +Changed: The interface of FE::fill_fe_face_values() now accepts instead +of a Quadrature instance a hp::QCollection instance, enabling the evaluation of +shape functions for different quadrature formulas on different faces for +FiniteElement classes with supports for this feature (e.g. FE_Q, Simplex::FE_P). +
+(Peter Munch, 2020/12/12) diff --git a/include/deal.II/base/qprojector.h b/include/deal.II/base/qprojector.h index 2381799512..7ec022daf9 100644 --- a/include/deal.II/base/qprojector.h +++ b/include/deal.II/base/qprojector.h @@ -24,6 +24,8 @@ #include +#include + DEAL_II_NAMESPACE_OPEN @@ -232,12 +234,12 @@ public: * version of this function that takes the reference cell type instead. */ DEAL_II_DEPRECATED static Quadrature - project_to_all_faces(const SubQuadrature &quadrature); + project_to_all_faces(const Quadrature &quadrature); /** - * Take a face quadrature formula and generate a cell quadrature formula - * from it where the quadrature points of the given argument are projected - * on all faces. + * Take a collection of face quadrature formulas and generate a cell + * quadrature formula from it where the quadrature points of the given + * argument are projected on all faces. * * The weights of the new rule are replications of the original weights. * Thus, the sum of the weights is not one, but the number of faces, which @@ -251,8 +253,16 @@ public: * each face, in order to cope possibly different orientations of the mesh. */ static Quadrature - project_to_all_faces(const ReferenceCell::Type reference_cell_type, - const SubQuadrature & quadrature); + project_to_all_faces(const ReferenceCell::Type reference_cell_type, + const hp::QCollection &quadrature); + + /** + * Like the above function, applying the same face quadrature + * formula on all faces. + */ + static Quadrature + project_to_all_faces(const ReferenceCell::Type reference_cell_type, + const Quadrature &quadrature); /** * Take a face quadrature formula and generate a cell quadrature formula @@ -444,6 +454,18 @@ public: const bool face_rotation, const unsigned int n_quadrature_points); + /** + * Like the above function but taking a quadrature collection, enabling + * that each face might have different number of quadrature points. + */ + static DataSetDescriptor + face(const ReferenceCell::Type reference_cell_type, + const unsigned int face_no, + const bool face_orientation, + const bool face_flip, + const bool face_rotation, + const hp::QCollection &quadrature); + /** * Static function to generate an offset object for a given subface of a * cell with the given face orientation, flip and rotation. This function @@ -575,6 +597,25 @@ inline QProjector::DataSetDescriptor::operator unsigned int() const } + +template +Quadrature inline QProjector::project_to_all_faces( + const Quadrature &quadrature) +{ + return project_to_all_faces(ReferenceCell::get_hypercube(dim), quadrature); +} + + +template +Quadrature inline QProjector::project_to_all_faces( + const ReferenceCell::Type reference_cell_type, + const Quadrature &quadrature) +{ + return project_to_all_faces(reference_cell_type, + hp::QCollection(quadrature)); +} + + /* -------------- declaration of explicit specializations ------------- */ #ifndef DOXYGEN @@ -614,14 +655,11 @@ QProjector<3>::project_to_face(const ReferenceCell::Type reference_cell_type, const unsigned int face_no, std::vector> & q_points); -template <> -Quadrature<1> -QProjector<1>::project_to_all_faces(const Quadrature<0> &quadrature); template <> Quadrature<1> QProjector<1>::project_to_all_faces( const ReferenceCell::Type reference_cell_type, - const Quadrature<0> & quadrature); + const hp::QCollection<0> &quadrature); template <> diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h index d0e58d067c..318ac01493 100644 --- a/include/deal.II/fe/fe.h +++ b/include/deal.II/fe/fe.h @@ -2767,12 +2767,23 @@ protected: * it is no longer necessary. */ virtual std::unique_ptr - get_face_data(const UpdateFlags update_flags, - const Mapping &mapping, - const Quadrature & quadrature, + get_face_data(const UpdateFlags update_flags, + const Mapping & mapping, + const hp::QCollection &quadrature, dealii::internal::FEValuesImplementation:: FiniteElementRelatedData &output_data) const; + /** + * @deprecated Use the version taking a hp::QCollection argument. + */ + virtual std::unique_ptr + get_face_data( + const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature & quadrature, + internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const; + /** * Like get_data(), but return an object that will later be used for * evaluating shape function information at quadrature points on children of @@ -2964,7 +2975,7 @@ protected: fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase &mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData - &output_data) const = 0; + &output_data) const; + + /** + * @deprecated Use the version taking a hp::QCollection argument. + */ + virtual void + fill_fe_face_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature & quadrature, + const Mapping & mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const internal::FEValuesImplementation::MappingRelatedData + & mapping_data, + const InternalDataBase &fe_internal, + internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const; /** * This function is the equivalent to FiniteElement::fill_fe_values(), but diff --git a/include/deal.II/fe/fe_dgp_nonparametric.h b/include/deal.II/fe/fe_dgp_nonparametric.h index 4944c55abc..27a752a9ae 100644 --- a/include/deal.II/fe/fe_dgp_nonparametric.h +++ b/include/deal.II/fe/fe_dgp_nonparametric.h @@ -534,11 +534,13 @@ protected: spacedim> &output_data) const override; + using FiniteElement::fill_fe_face_values; + virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase &mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + using FiniteElement::get_face_data; + virtual std::unique_ptr< typename FiniteElement::InternalDataBase> get_face_data( - const UpdateFlags update_flags, - const Mapping &mapping, - const Quadrature & quadrature, + const UpdateFlags update_flags, + const Mapping & mapping, + const hp::QCollection &quadrature, dealii::internal::FEValuesImplementation::FiniteElementRelatedData &output_data) const override; @@ -610,11 +612,13 @@ protected: spacedim> &output_data) const override; + using FiniteElement::fill_fe_face_values; + virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase &mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData::InternalDataBase>(); } + using FiniteElement<1, spacedim>::get_face_data; + std::unique_ptr::InternalDataBase> get_face_data( const UpdateFlags update_flags, const Mapping<1, spacedim> & /*mapping*/, - const Quadrature<0> &quadrature, + const hp::QCollection<0> &quadrature, dealii::internal::FEValuesImplementation::FiniteElementRelatedData<1, spacedim> & /*output_data*/) const override { + AssertDimension(quadrature.size(), 1); + // generate a new data object and initialize some fields auto data_ptr = std::make_unique::InternalDataBase>(); data_ptr->update_each = requires_update_flags(update_flags); - const unsigned int n_q_points = quadrature.size(); + const unsigned int n_q_points = quadrature[0].size(); AssertDimension(n_q_points, 1); (void)n_q_points; @@ -386,7 +390,10 @@ protected: spacedim> &output_data) const override { - return get_face_data(update_flags, mapping, quadrature, output_data); + return get_face_data(update_flags, + mapping, + hp::QCollection<0>(quadrature), + output_data); } virtual void @@ -404,11 +411,13 @@ protected: spacedim> &output_data) const override; + using FiniteElement<1, spacedim>::fill_fe_face_values; + virtual void fill_fe_face_values( const typename Triangulation<1, spacedim>::cell_iterator &cell, const unsigned int face_no, - const Quadrature<0> & quadrature, + const hp::QCollection<0> & quadrature, const Mapping<1, spacedim> & mapping, const typename Mapping<1, spacedim>::InternalDataBase & mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData<1, diff --git a/include/deal.II/fe/fe_nedelec_sz.h b/include/deal.II/fe/fe_nedelec_sz.h index 7fca98c6a7..50801f6366 100644 --- a/include/deal.II/fe/fe_nedelec_sz.h +++ b/include/deal.II/fe/fe_nedelec_sz.h @@ -179,6 +179,8 @@ protected: dealii::internal::FEValuesImplementation::FiniteElementRelatedData &data) const override; + using FiniteElement::fill_fe_face_values; + /** * Compute information about the shape functions on the cell and face denoted * by the first two arguments. Note that this function must recompute the @@ -188,7 +190,7 @@ protected: fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase & mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData diff --git a/include/deal.II/fe/fe_nothing.h b/include/deal.II/fe/fe_nothing.h index 4660237f35..008bf663fe 100644 --- a/include/deal.II/fe/fe_nothing.h +++ b/include/deal.II/fe/fe_nothing.h @@ -130,11 +130,13 @@ public: spacedim> &output_data) const override; + using FiniteElement::fill_fe_face_values; + virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase &mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + using FiniteElement<2, 2>::get_face_data; + virtual std::unique_ptr::InternalDataBase> get_face_data( const UpdateFlags update_flags, const Mapping<2, 2> &, - const Quadrature<1> &quadrature, + const hp::QCollection<1> &quadrature, dealii::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2> &output_data) const override; @@ -343,6 +345,8 @@ private: internal::FEValuesImplementation::FiniteElementRelatedData<2, 2> &output_data) const override; + using FiniteElement<2, 2>::fill_fe_face_values; + /** * Compute the data on the face of the current cell. */ @@ -350,7 +354,7 @@ private: fill_fe_face_values( const Triangulation<2, 2>::cell_iterator &cell, const unsigned int face_no, - const Quadrature<1> & quadrature, + const hp::QCollection<1> & quadrature, const Mapping<2, 2> & mapping, const Mapping<2, 2>::InternalDataBase & mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData<2, 2> diff --git a/include/deal.II/fe/fe_poly.h b/include/deal.II/fe/fe_poly.h index d6ea12ba42..50d6a5652a 100644 --- a/include/deal.II/fe/fe_poly.h +++ b/include/deal.II/fe/fe_poly.h @@ -387,11 +387,13 @@ protected: spacedim> &output_data) const override; + using FiniteElement::fill_fe_face_values; + virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase &mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase & mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData::fill_fe_face_values( spacedim> &output_data) const { + const unsigned int n_q_points = + quadrature[quadrature.size() == 1 ? 0 : face_no].size(); + // convert data object to internal // data for this class. fails with // an exception if that is not @@ -365,14 +368,14 @@ FE_Poly::fill_fe_face_values( cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), - quadrature.size()); + quadrature); const UpdateFlags flags(fe_data.update_each); const bool need_to_correct_higher_derivatives = higher_derivatives_need_correcting(mapping, mapping_data, - quadrature.size(), + n_q_points, flags); // transform gradients and higher derivatives. we also have to copy @@ -380,13 +383,13 @@ FE_Poly::fill_fe_face_values( // we need to take into account the offsets if (flags & update_values) for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - for (unsigned int i = 0; i < quadrature.size(); ++i) + for (unsigned int i = 0; i < n_q_points; ++i) output_data.shape_values(k, i) = fe_data.shape_values[k][i + offset]; if (flags & update_gradients) for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) mapping.transform( - make_array_view(fe_data.shape_gradients, k, offset, quadrature.size()), + make_array_view(fe_data.shape_gradients, k, offset, n_q_points), mapping_covariant, mapping_internal, make_array_view(output_data.shape_gradients, k)); @@ -395,29 +398,26 @@ FE_Poly::fill_fe_face_values( { for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) mapping.transform( - make_array_view(fe_data.shape_hessians, k, offset, quadrature.size()), + make_array_view(fe_data.shape_hessians, k, offset, n_q_points), mapping_covariant_gradient, mapping_internal, make_array_view(output_data.shape_hessians, k)); if (need_to_correct_higher_derivatives) - correct_hessians(output_data, mapping_data, quadrature.size()); + correct_hessians(output_data, mapping_data, n_q_points); } if (flags & update_3rd_derivatives) { for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) - mapping.transform(make_array_view(fe_data.shape_3rd_derivatives, - k, - offset, - quadrature.size()), - mapping_covariant_hessian, - mapping_internal, - make_array_view(output_data.shape_3rd_derivatives, - k)); + mapping.transform( + make_array_view(fe_data.shape_3rd_derivatives, k, offset, n_q_points), + mapping_covariant_hessian, + mapping_internal, + make_array_view(output_data.shape_3rd_derivatives, k)); if (need_to_correct_higher_derivatives) - correct_third_derivatives(output_data, mapping_data, quadrature.size()); + correct_third_derivatives(output_data, mapping_data, n_q_points); } } diff --git a/include/deal.II/fe/fe_poly_face.h b/include/deal.II/fe/fe_poly_face.h index f8d9695a3a..ba4ffe6d54 100644 --- a/include/deal.II/fe/fe_poly_face.h +++ b/include/deal.II/fe/fe_poly_face.h @@ -102,17 +102,21 @@ protected: return std::make_unique(); } + using FiniteElement::get_face_data; + std::unique_ptr::InternalDataBase> get_face_data( - const UpdateFlags update_flags, - const Mapping &mapping, - const Quadrature & quadrature, + const UpdateFlags update_flags, + const Mapping & mapping, + const hp::QCollection &quadrature, dealii::internal::FEValuesImplementation::FiniteElementRelatedData &output_data) const override { (void)mapping; (void)output_data; + AssertDimension(quadrature.size(), 1); + // generate a new data object and // initialize some fields std::unique_ptr::InternalDataBase> @@ -120,7 +124,7 @@ protected: auto &data = dynamic_cast(*data_ptr); data.update_each = requires_update_flags(update_flags); - const unsigned int n_q_points = quadrature.size(); + const unsigned int n_q_points = quadrature[0].size(); // some scratch arrays std::vector values(0); @@ -139,7 +143,7 @@ protected: std::vector(n_q_points)); for (unsigned int i = 0; i < n_q_points; ++i) { - poly_space.evaluate(quadrature.point(i), + poly_space.evaluate(quadrature[0].point(i), values, grads, grad_grads, @@ -170,11 +174,12 @@ protected: spacedim> &output_data) const override { - return get_face_data(update_flags, - mapping, - QProjector::project_to_all_children( - ReferenceCell::get_hypercube(dim - 1), quadrature), - output_data); + return get_face_data( + update_flags, + mapping, + hp::QCollection(QProjector::project_to_all_children( + ReferenceCell::get_hypercube(dim - 1), quadrature)), + output_data); } virtual void @@ -192,11 +197,13 @@ protected: spacedim> &output_data) const override; + using FiniteElement::fill_fe_face_values; + virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase &mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData void FE_PolyFace::fill_fe_face_values( const typename Triangulation::cell_iterator &, - const unsigned int face_no, - const Quadrature &quadrature, + const unsigned int face_no, + const hp::QCollection &quadrature, const Mapping &, const typename Mapping::InternalDataBase &, const dealii::internal::FEValuesImplementation::MappingRelatedData::fill_fe_face_values( ExcInternalError()); const InternalData &fe_data = static_cast(fe_internal); + AssertDimension(quadrature.size(), 1); + if (fe_data.update_each & update_values) - for (unsigned int i = 0; i < quadrature.size(); ++i) + for (unsigned int i = 0; i < quadrature[0].size(); ++i) { for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k) output_data.shape_values(k, i) = 0.; diff --git a/include/deal.II/fe/fe_poly_tensor.h b/include/deal.II/fe/fe_poly_tensor.h index d12c4c9145..9ccd27b8f6 100644 --- a/include/deal.II/fe/fe_poly_tensor.h +++ b/include/deal.II/fe/fe_poly_tensor.h @@ -392,11 +392,13 @@ protected: spacedim> &output_data) const override; + using FiniteElement::fill_fe_face_values; + virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase &mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + using FiniteElement::get_face_data; + virtual std::unique_ptr< typename FiniteElement::InternalDataBase> get_face_data( - const UpdateFlags update_flags, - const Mapping &mapping, - const Quadrature & quadrature, + const UpdateFlags update_flags, + const Mapping & mapping, + const hp::QCollection &quadrature, dealii::internal::FEValuesImplementation::FiniteElementRelatedData &output_data) const override; @@ -1097,11 +1099,13 @@ protected: spacedim> &output_data) const override; + using FiniteElement::fill_fe_face_values; + virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase &mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const CellSimilarity::Similarity cell_similarity, const typename Mapping::InternalDataBase &mapping_internal, const typename FiniteElement::InternalDataBase &fe_data, diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h index 0972902873..292df40c99 100644 --- a/include/deal.II/fe/fe_values.h +++ b/include/deal.II/fe/fe_values.h @@ -38,6 +38,8 @@ #include #include +#include + #include #include #include @@ -2175,10 +2177,25 @@ public: static const unsigned int space_dimension = spacedim; /** - * Number of quadrature points. + * Number of quadrature points of the current object. Its value is + * initialized by the value of max_n_quadrature_points and is updated, + * e.g., if FEFaceValues::reinit() is called for a new cell/face. + * + * @note The default value equals to the value of max_n_quadrature_points. */ const unsigned int n_quadrature_points; + /** + * Maximum number of quadrature points. This value might be different from + * n_quadrature_points, e.g., if a QCollection with different face quadrature + * rules has been passed to initialize FEFaceValues. + * + * This is mostly useful to initialize arrays to allocate the maximum amount + * of memory that may be used when re-sizing later on to a the current + * number of quadrature points given by n_quadrature_points. + */ + const unsigned int max_n_quadrature_points; + /** * Number of shape functions per cell. If we use this base class to evaluate * a finite element on faces of cells, this is still the number of degrees @@ -3745,6 +3762,17 @@ public: const Quadrature & quadrature, const UpdateFlags update_flags); + /** + * Like the function above, but taking a collection of quadrature rules. + * + * @note We require, in contrast to FEFaceValues, that the number of quadrature + * rules in the collection is one. + */ + FEValues(const Mapping & mapping, + const FiniteElement &fe, + const hp::QCollection & quadrature, + const UpdateFlags update_flags); + /** * Constructor. This constructor is equivalent to the other one except that * it makes the object use a $Q_1$ mapping (i.e., an object of type @@ -3754,6 +3782,16 @@ public: const Quadrature & quadrature, const UpdateFlags update_flags); + /** + * Like the function above, but taking a collection of quadrature rules. + * + * @note We require, in contrast to FEFaceValues, that the number of quadrature + * rules in the collection is one. + */ + FEValues(const FiniteElement &fe, + const hp::QCollection & quadrature, + const UpdateFlags update_flags); + /** * Reinitialize the gradients, Jacobi determinants, etc for the given cell * of type "iterator into a DoFHandler object", and the finite element @@ -3865,13 +3903,24 @@ public: * is 2*dim*(1<<(dim-1)), i.e. the number of faces times the number * of subfaces per face. */ - FEFaceValuesBase(const unsigned int n_q_points, - const unsigned int dofs_per_cell, + FEFaceValuesBase(const unsigned int dofs_per_cell, const UpdateFlags update_flags, const Mapping & mapping, const FiniteElement &fe, const Quadrature & quadrature); + /** + * Like the function above, but taking a collection of quadrature rules. This + * allows to assign each face a different quadrature rule. In the case that + * the collection only contains a single face quadrature, this quadrature + * rule is use on all faces. + */ + FEFaceValuesBase(const unsigned int dofs_per_cell, + const UpdateFlags update_flags, + const Mapping & mapping, + const FiniteElement &fe, + const hp::QCollection & quadrature); + /** * Boundary form of the transformation of the cell at the ith * quadrature point. See @@ -3922,7 +3971,7 @@ protected: /** * Store a copy of the quadrature formula here. */ - const Quadrature quadrature; + const hp::QCollection quadrature; }; @@ -3967,6 +4016,17 @@ public: const Quadrature & quadrature, const UpdateFlags update_flags); + /** + * Like the function above, but taking a collection of quadrature rules. This + * allows to assign each face a different quadrature rule. In the case that + * the collection only contains a single face quadrature, this quadrature + * rule is use on all faces. + */ + FEFaceValues(const Mapping & mapping, + const FiniteElement &fe, + const hp::QCollection & quadrature, + const UpdateFlags update_flags); + /** * Constructor. This constructor is equivalent to the other one except that * it makes the object use a $Q_1$ mapping (i.e., an object of type @@ -3976,6 +4036,16 @@ public: const Quadrature & quadrature, const UpdateFlags update_flags); + /** + * Like the function above, but taking a collection of quadrature rules. This + * allows to assign each face a different quadrature rule. In the case that + * the collection only contains a single face quadrature, this quadrature + * rule is use on all faces. + */ + FEFaceValues(const FiniteElement &fe, + const hp::QCollection & quadrature, + const UpdateFlags update_flags); + /** * Reinitialize the gradients, Jacobi determinants, etc for the face with * number @p face_no of @p cell and the given finite element. @@ -4114,6 +4184,17 @@ public: const Quadrature & face_quadrature, const UpdateFlags update_flags); + /** + * Like the function above, but taking a collection of quadrature rules. + * + * @note We require, in contrast to FEFaceValues, that the number of quadrature + * rules in the collection is one. + */ + FESubfaceValues(const Mapping & mapping, + const FiniteElement &fe, + const hp::QCollection & face_quadrature, + const UpdateFlags update_flags); + /** * Constructor. This constructor is equivalent to the other one except that * it makes the object use a $Q_1$ mapping (i.e., an object of type @@ -4123,6 +4204,16 @@ public: const Quadrature & face_quadrature, const UpdateFlags update_flags); + /** + * Like the function above, but taking a collection of quadrature rules. + * + * @note We require, in contrast to FEFaceValues, that the number of quadrature + * rules in the collection is one. + */ + FESubfaceValues(const FiniteElement &fe, + const hp::QCollection & face_quadrature, + const UpdateFlags update_flags); + /** * Reinitialize the gradients, Jacobi determinants, etc for the given cell * of type "iterator into a DoFHandler object", and the finite element @@ -5810,7 +5901,7 @@ template inline const Quadrature & FEFaceValuesBase::get_quadrature() const { - return quadrature; + return quadrature[0]; } diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index 0f3c817689..f04ee0890a 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -26,6 +26,8 @@ #include +#include + #include #include #include @@ -757,8 +759,15 @@ protected: * the returned object, knowing its real (derived) type. */ virtual std::unique_ptr + get_face_data(const UpdateFlags update_flags, + const hp::QCollection &quadrature) const; + + /** + * @deprecated Use the version taking a hp::QCollection argument. + */ + virtual std::unique_ptr get_face_data(const UpdateFlags update_flags, - const Quadrature &quadrature) const = 0; + const Quadrature &quadrature) const; /** * Like get_data() and get_face_data(), but in preparation for later calls @@ -912,10 +921,22 @@ protected: fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, dealii::internal::FEValuesImplementation::MappingRelatedData - &output_data) const = 0; + &output_data) const; + + /** + * @deprecated Use the version taking a hp::QCollection argument. + */ + virtual void + fill_fe_face_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature & quadrature, + const typename Mapping::InternalDataBase & internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const; /** * This function is the equivalent to Mapping::fill_fe_values(), but for diff --git a/include/deal.II/fe/mapping_cartesian.h b/include/deal.II/fe/mapping_cartesian.h index 7826b7f02e..5fbb61f554 100644 --- a/include/deal.II/fe/mapping_cartesian.h +++ b/include/deal.II/fe/mapping_cartesian.h @@ -210,10 +210,12 @@ private: virtual std::unique_ptr::InternalDataBase> get_data(const UpdateFlags, const Quadrature &quadrature) const override; + using Mapping::get_face_data; + // documentation can be found in Mapping::get_face_data() virtual std::unique_ptr::InternalDataBase> - get_face_data(const UpdateFlags flags, - const Quadrature &quadrature) const override; + get_face_data(const UpdateFlags flags, + const hp::QCollection &quadrature) const override; // documentation can be found in Mapping::get_subface_data() virtual std::unique_ptr::InternalDataBase> @@ -230,12 +232,14 @@ private: internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + using Mapping::fill_fe_face_values; + // documentation can be found in Mapping::fill_fe_face_values() virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, internal::FEValuesImplementation::MappingRelatedData &output_data) const override; diff --git a/include/deal.II/fe/mapping_fe.h b/include/deal.II/fe/mapping_fe.h index fb26cd76da..0865491ee0 100644 --- a/include/deal.II/fe/mapping_fe.h +++ b/include/deal.II/fe/mapping_fe.h @@ -401,10 +401,12 @@ public: virtual std::unique_ptr::InternalDataBase> get_data(const UpdateFlags, const Quadrature &quadrature) const override; + using Mapping::get_face_data; + // documentation can be found in Mapping::get_face_data() virtual std::unique_ptr::InternalDataBase> - get_face_data(const UpdateFlags flags, - const Quadrature &quadrature) const override; + get_face_data(const UpdateFlags flags, + const hp::QCollection &quadrature) const override; // documentation can be found in Mapping::get_subface_data() virtual std::unique_ptr::InternalDataBase> @@ -421,12 +423,14 @@ public: dealii::internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + using Mapping::fill_fe_face_values; + // documentation can be found in Mapping::fill_fe_face_values() virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, dealii::internal::FEValuesImplementation::MappingRelatedData &output_data) const override; diff --git a/include/deal.II/fe/mapping_fe_field.h b/include/deal.II/fe/mapping_fe_field.h index 8814ae17d4..494f5bc8ac 100644 --- a/include/deal.II/fe/mapping_fe_field.h +++ b/include/deal.II/fe/mapping_fe_field.h @@ -531,10 +531,12 @@ private: virtual std::unique_ptr::InternalDataBase> get_data(const UpdateFlags, const Quadrature &quadrature) const override; + using Mapping::get_face_data; + // documentation can be found in Mapping::get_face_data() virtual std::unique_ptr::InternalDataBase> - get_face_data(const UpdateFlags flags, - const Quadrature &quadrature) const override; + get_face_data(const UpdateFlags flags, + const hp::QCollection &quadrature) const override; // documentation can be found in Mapping::get_subface_data() virtual std::unique_ptr::InternalDataBase> @@ -551,12 +553,14 @@ private: internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + using Mapping::fill_fe_face_values; + // documentation can be found in Mapping::fill_fe_face_values() virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, internal::FEValuesImplementation::MappingRelatedData &output_data) const override; diff --git a/include/deal.II/fe/mapping_manifold.h b/include/deal.II/fe/mapping_manifold.h index 7fb421cf4f..60c3f4b206 100644 --- a/include/deal.II/fe/mapping_manifold.h +++ b/include/deal.II/fe/mapping_manifold.h @@ -338,10 +338,12 @@ public: virtual std::unique_ptr::InternalDataBase> get_data(const UpdateFlags, const Quadrature &quadrature) const override; + using Mapping::get_face_data; + // documentation can be found in Mapping::get_face_data() virtual std::unique_ptr::InternalDataBase> - get_face_data(const UpdateFlags flags, - const Quadrature &quadrature) const override; + get_face_data(const UpdateFlags flags, + const hp::QCollection &quadrature) const override; // documentation can be found in Mapping::get_subface_data() virtual std::unique_ptr::InternalDataBase> @@ -358,12 +360,14 @@ public: dealii::internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + using Mapping::fill_fe_face_values; + // documentation can be found in Mapping::fill_fe_face_values() virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, dealii::internal::FEValuesImplementation::MappingRelatedData &output_data) const override; diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 6a61db592e..47ad2bf5d8 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -271,10 +271,12 @@ protected: virtual std::unique_ptr::InternalDataBase> get_data(const UpdateFlags, const Quadrature &quadrature) const override; + using Mapping::get_face_data; + // documentation can be found in Mapping::get_face_data() virtual std::unique_ptr::InternalDataBase> - get_face_data(const UpdateFlags flags, - const Quadrature &quadrature) const override; + get_face_data(const UpdateFlags flags, + const hp::QCollection &quadrature) const override; // documentation can be found in Mapping::get_subface_data() virtual std::unique_ptr::InternalDataBase> @@ -291,12 +293,14 @@ protected: internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + using Mapping::fill_fe_face_values; + // documentation can be found in Mapping::fill_fe_face_values() virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, internal::FEValuesImplementation::MappingRelatedData &output_data) const override; diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index 8a7805c860..a3000b187b 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -568,10 +568,12 @@ public: virtual std::unique_ptr::InternalDataBase> get_data(const UpdateFlags, const Quadrature &quadrature) const override; + using Mapping::get_face_data; + // documentation can be found in Mapping::get_face_data() virtual std::unique_ptr::InternalDataBase> - get_face_data(const UpdateFlags flags, - const Quadrature &quadrature) const override; + get_face_data(const UpdateFlags flags, + const hp::QCollection &quadrature) const override; // documentation can be found in Mapping::get_subface_data() virtual std::unique_ptr::InternalDataBase> @@ -588,12 +590,14 @@ public: dealii::internal::FEValuesImplementation::MappingRelatedData &output_data) const override; + using Mapping::fill_fe_face_values; + // documentation can be found in Mapping::fill_fe_face_values() virtual void fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, dealii::internal::FEValuesImplementation::MappingRelatedData &output_data) const override; diff --git a/include/deal.II/hp/fe_values.h b/include/deal.II/hp/fe_values.h index 23ce1556b3..280047e70c 100644 --- a/include/deal.II/hp/fe_values.h +++ b/include/deal.II/hp/fe_values.h @@ -80,6 +80,18 @@ namespace hp const QCollection & q_collection, const UpdateFlags update_flags); + /** + * Like the above function but taking a vector of quadrature collections. + * For hp::FEFaceValues, the ith entry of the quadrature collections are + * interpreted as the face quadrature rules to be applied the ith face. + */ + FEValuesBase( + const MappingCollection + &mapping_collection, + const FECollection &fe_collection, + const std::vector> & q_collection, + const UpdateFlags update_flags); + /** * Constructor. This constructor is equivalent to the other one except * that it makes the object use a $Q_1$ mapping (i.e., an object of type @@ -90,6 +102,16 @@ namespace hp const QCollection & q_collection, const UpdateFlags update_flags); + /** + * Like the above function but taking a vector quadrature collections. + * For hp::FEFaceValues, the ith entry of the quadrature collections are + * interpreted as the face quadrature rules to be applied the ith face. + */ + FEValuesBase( + const FECollection &fe_collection, + const std::vector> & q_collection, + const UpdateFlags update_flags); + /** * Copy constructor. */ @@ -199,6 +221,16 @@ namespace hp */ const QCollection q_collection; + /** + * Vector of quadrature collections. For hp::FEFaceValues, the ith entry of + * the quadrature collections are interpreted as the face quadrature rules + * to be applied the ith face. + * + * The variable q_collection collects the first quadrature rule of each + * quadrature collection of the vector. + */ + const std::vector> q_collections; + private: /** * A table in which we store pointers to fe_values objects for different @@ -416,6 +448,20 @@ namespace hp const hp::QCollection & q_collection, const UpdateFlags update_flags); + /** + * Like the function above, but taking a vector of collection of quadrature + * rules. This allows to assign each face a different quadrature rule: the + * ith entry of a collection is used as the face quadrature rule on the ith + * face. + * + * In the case that the collections only contains a single face quadrature, + * this quadrature rule is use on all faces. + */ + FEFaceValues(const hp::MappingCollection &mapping_collection, + const hp::FECollection & fe_collection, + const std::vector> &q_collections, + const UpdateFlags update_flags); + /** * Constructor. This constructor is equivalent to the other one except @@ -426,6 +472,19 @@ namespace hp const hp::QCollection & q_collection, const UpdateFlags update_flags); + /** + * Like the function above, but taking a vector of collection of quadrature + * rules. This allows to assign each face a different quadrature rule: the + * ith entry of a collection is used as the face quadrature rule on the ith + * face. + * + * In the case that the collections only contains a single face quadrature, + * this quadrature rule is use on all faces. + */ + FEFaceValues(const hp::FECollection & fe_collection, + const std::vector> &q_collections, + const UpdateFlags update_flags); + /** * Reinitialize the object for the given cell and face. * diff --git a/include/deal.II/hp/q_collection.h b/include/deal.II/hp/q_collection.h index 073bac7968..46339178cd 100644 --- a/include/deal.II/hp/q_collection.h +++ b/include/deal.II/hp/q_collection.h @@ -22,8 +22,6 @@ #include #include -#include - #include #include diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index e3e7b4e133..5515b54652 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -465,20 +465,13 @@ QProjector<3>::project_to_subface(const ReferenceCell::Type reference_cell_type, } -template <> -Quadrature<1> -QProjector<1>::project_to_all_faces(const Quadrature<0> &quadrature) -{ - return project_to_all_faces(ReferenceCell::Type::Line, quadrature); -} - - template <> Quadrature<1> QProjector<1>::project_to_all_faces( const ReferenceCell::Type reference_cell_type, - const Quadrature<0> & quadrature) + const hp::QCollection<0> &quadrature) { + AssertDimension(quadrature.size(), 1); Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented()); (void)reference_cell_type; @@ -496,7 +489,9 @@ QProjector<1>::project_to_all_faces( // results for (unsigned int face = 0; face < n_faces; ++face) { - project_to_face(quadrature, face, help); + project_to_face(quadrature[quadrature.size() == 1 ? 0 : face], + face, + help); std::copy(help.begin(), help.end(), std::back_inserter(q_points)); } @@ -504,9 +499,10 @@ QProjector<1>::project_to_all_faces( std::vector weights; weights.reserve(n_points * n_faces); for (unsigned int face = 0; face < n_faces; ++face) - std::copy(quadrature.get_weights().begin(), - quadrature.get_weights().end(), - std::back_inserter(weights)); + std::copy( + quadrature[quadrature.size() == 1 ? 0 : face].get_weights().begin(), + quadrature[quadrature.size() == 1 ? 0 : face].get_weights().end(), + std::back_inserter(weights)); Assert(q_points.size() == n_points * n_faces, ExcInternalError()); Assert(weights.size() == n_points * n_faces, ExcInternalError()); @@ -520,15 +516,11 @@ template <> Quadrature<2> QProjector<2>::project_to_all_faces( const ReferenceCell::Type reference_cell_type, - const SubQuadrature & quadrature) + const hp::QCollection<1> &quadrature) { if (reference_cell_type == ReferenceCell::Type::Tri) { - // the quadrature rule to be projected ... - const auto &sub_quadrature_points = quadrature.get_points(); - const auto &sub_quadrature_weights = quadrature.get_weights(); - - // ... on to the faces (defined by its support points and arc length) + // reference faces (defined by its support points and arc length) const std::array, 2>, double>, 3> faces = { {{{{Point<2>(0.0, 0.0), Point<2>(1.0, 0.0)}}, 1.0}, {{{Point<2>(1.0, 0.0), Point<2>(0.0, 1.0)}}, std::sqrt(2.0)}, @@ -543,10 +535,12 @@ QProjector<2>::project_to_all_faces( std::vector weights; // loop over all faces (lines) ... - for (const auto &face : faces) + for (unsigned int face_no = 0; face_no < faces.size(); ++face_no) // ... and over all possible orientations for (unsigned int orientation = 0; orientation < 2; ++orientation) { + const auto &face = faces[face_no]; + std::array, 2> support_points; // determine support point of the current line with the correct @@ -563,6 +557,12 @@ QProjector<2>::project_to_all_faces( Assert(false, ExcNotImplemented()); } + // the quadrature rule to be projected ... + const auto &sub_quadrature_points = + quadrature[quadrature.size() == 1 ? 0 : face_no].get_points(); + const auto &sub_quadrature_weights = + quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights(); + // loop over all quadrature points for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j) { @@ -589,60 +589,62 @@ QProjector<2>::project_to_all_faces( const unsigned int dim = 2; - const unsigned int n_points = quadrature.size(), - n_faces = GeometryInfo::faces_per_cell; + const unsigned int n_faces = GeometryInfo::faces_per_cell; + + unsigned int n_points_total = 0; + + if (quadrature.size() == 1) + n_points_total = quadrature[0].size() * GeometryInfo::faces_per_cell; + else + { + AssertDimension(quadrature.size(), GeometryInfo::faces_per_cell); + for (unsigned int i = 0; i < quadrature.size(); ++i) + n_points_total += quadrature[i].size(); + } // first fix quadrature points std::vector> q_points; - q_points.reserve(n_points * n_faces); - std::vector> help(n_points); + q_points.reserve(n_points_total); + std::vector> help; + help.reserve(quadrature.max_n_quadrature_points()); // project to each face and append // results for (unsigned int face = 0; face < n_faces; ++face) { - project_to_face(quadrature, face, help); + help.resize(quadrature[quadrature.size() == 1 ? 0 : face].size()); + project_to_face(quadrature[quadrature.size() == 1 ? 0 : face], + face, + help); std::copy(help.begin(), help.end(), std::back_inserter(q_points)); } // next copy over weights std::vector weights; - weights.reserve(n_points * n_faces); + weights.reserve(n_points_total); for (unsigned int face = 0; face < n_faces; ++face) - std::copy(quadrature.get_weights().begin(), - quadrature.get_weights().end(), - std::back_inserter(weights)); + std::copy( + quadrature[quadrature.size() == 1 ? 0 : face].get_weights().begin(), + quadrature[quadrature.size() == 1 ? 0 : face].get_weights().end(), + std::back_inserter(weights)); - Assert(q_points.size() == n_points * n_faces, ExcInternalError()); - Assert(weights.size() == n_points * n_faces, ExcInternalError()); + Assert(q_points.size() == n_points_total, ExcInternalError()); + Assert(weights.size() == n_points_total, ExcInternalError()); return Quadrature(q_points, weights); } -template <> -Quadrature<2> -QProjector<2>::project_to_all_faces(const SubQuadrature &quadrature) -{ - return project_to_all_faces(ReferenceCell::Type::Quad, quadrature); -} - - - template <> Quadrature<3> QProjector<3>::project_to_all_faces( const ReferenceCell::Type reference_cell_type, - const SubQuadrature & quadrature) + const hp::QCollection<2> &quadrature) { if (reference_cell_type == ReferenceCell::Type::Tet) { - // the quadrature rule to be projected ... - const auto &sub_quadrature_points = quadrature.get_points(); - const auto &sub_quadrature_weights = quadrature.get_weights(); - - // ... on to the faces (defined by its support points and its area) + // reference faces (defined by its support points and its area) // note: the area is later not used as a scaling factor but recomputed const std::array, 3>, double>, 4> faces = { {{{{Point<3>(0.0, 0.0, 0.0), @@ -671,10 +673,12 @@ QProjector<3>::project_to_all_faces( std::vector weights; // loop over all faces (triangles) ... - for (const auto &face : faces) + for (unsigned int face_no = 0; face_no < faces.size(); ++face_no) // ... and over all possible orientations for (unsigned int orientation = 0; orientation < 6; ++orientation) { + const auto &face = faces[face_no]; + std::array, 3> support_points; // determine support point of the current line with the correct @@ -709,6 +713,13 @@ QProjector<3>::project_to_all_faces( Assert(false, ExcNotImplemented()); } + + // the quadrature rule to be projected ... + const auto &sub_quadrature_points = + quadrature[quadrature.size() == 1 ? 0 : face_no].get_points(); + const auto &sub_quadrature_weights = + quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights(); + // loop over all quadrature points for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j) { @@ -781,68 +792,92 @@ QProjector<3>::project_to_all_faces( const unsigned int dim = 3; - SubQuadrature q_reflected = reflect(quadrature); - SubQuadrature q[8] = {quadrature, - rotate(quadrature, 1), - rotate(quadrature, 2), - rotate(quadrature, 3), - q_reflected, - rotate(q_reflected, 3), - rotate(q_reflected, 2), - rotate(q_reflected, 1)}; - + unsigned int n_points_total = 0; + if (quadrature.size() == 1) + n_points_total = quadrature[0].size() * GeometryInfo::faces_per_cell; + else + { + AssertDimension(quadrature.size(), GeometryInfo::faces_per_cell); + for (unsigned int i = 0; i < quadrature.size(); ++i) + n_points_total += quadrature[i].size(); + } - const unsigned int n_points = quadrature.size(), - n_faces = GeometryInfo::faces_per_cell; + n_points_total *= 8; // first fix quadrature points std::vector> q_points; - q_points.reserve(n_points * n_faces * 8); - std::vector> help(n_points); + q_points.reserve(n_points_total); + std::vector> help; + help.reserve(quadrature.max_n_quadrature_points()); std::vector weights; - weights.reserve(n_points * n_faces * 8); + weights.reserve(n_points_total); // do the following for all possible // mutations of a face (mutation==0 // corresponds to a face with standard // orientation, no flip and no rotation) - for (const auto &mutation : q) + for (unsigned int i = 0; i < 8; ++i) { // project to each face and append // results - for (unsigned int face = 0; face < n_faces; ++face) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) { + SubQuadrature mutation; + + const auto quadrature_f = + quadrature[quadrature.size() == 1 ? 0 : face]; + switch (i) + { + case 0: + mutation = quadrature_f; + break; + case 1: + mutation = rotate(quadrature_f, 1); + break; + case 2: + mutation = rotate(quadrature_f, 2); + break; + case 3: + mutation = rotate(quadrature_f, 3); + break; + case 4: + mutation = reflect(quadrature_f); + break; + case 5: + mutation = rotate(reflect(quadrature_f), 3); + break; + case 6: + mutation = rotate(reflect(quadrature_f), 2); + break; + case 7: + mutation = rotate(reflect(quadrature_f), 1); + break; + default: + Assert(false, ExcInternalError()) + } + + help.resize(quadrature[quadrature.size() == 1 ? 0 : face].size()); project_to_face(mutation, face, help); std::copy(help.begin(), help.end(), std::back_inserter(q_points)); - } - // next copy over weights - for (unsigned int face = 0; face < n_faces; ++face) - std::copy(mutation.get_weights().begin(), - mutation.get_weights().end(), - std::back_inserter(weights)); + std::copy(mutation.get_weights().begin(), + mutation.get_weights().end(), + std::back_inserter(weights)); + } } - Assert(q_points.size() == n_points * n_faces * 8, ExcInternalError()); - Assert(weights.size() == n_points * n_faces * 8, ExcInternalError()); + Assert(q_points.size() == n_points_total, ExcInternalError()); + Assert(weights.size() == n_points_total, ExcInternalError()); return Quadrature(q_points, weights); } -template <> -Quadrature<3> -QProjector<3>::project_to_all_faces(const SubQuadrature &quadrature) -{ - return project_to_all_faces(ReferenceCell::Type::Hex, quadrature); -} - - - template <> Quadrature<1> QProjector<1>::project_to_all_subfaces(const Quadrature<0> &quadrature) @@ -1285,6 +1320,125 @@ QProjector::DataSetDescriptor::face( +template +typename QProjector::DataSetDescriptor +QProjector::DataSetDescriptor::face( + const ReferenceCell::Type reference_cell_type, + const unsigned int face_no, + const bool face_orientation, + const bool face_flip, + const bool face_rotation, + const hp::QCollection &quadrature) +{ + if (reference_cell_type == ReferenceCell::Type::Tri || + reference_cell_type == ReferenceCell::Type::Tet) + { + unsigned int offset = 0; + + if (quadrature.size() == 1) + offset = quadrature[0].size() * face_no; + else + for (unsigned int i = 0; i < face_no; ++i) + offset += quadrature[i].size(); + + if (dim == 2) + return {2 * offset + + face_orientation * + quadrature[quadrature.size() == 1 ? 0 : face_no].size()}; + else if (dim == 3) + { + const unsigned int orientation = + (face_flip * 2 + face_rotation) * 2 + face_orientation; + + return {6 * offset + + orientation * + quadrature[quadrature.size() == 1 ? 0 : face_no].size()}; + } + } + + Assert(reference_cell_type == ReferenceCell::get_hypercube(dim), + ExcNotImplemented()); + + Assert(face_no < GeometryInfo::faces_per_cell, ExcInternalError()); + + switch (dim) + { + case 1: + case 2: + { + if (quadrature.size() == 1) + return quadrature[0].size() * face_no; + else + { + unsigned int result = 0; + for (unsigned int i = 0; i < face_no; ++i) + result += quadrature[i].size(); + return result; + } + } + case 3: + { + // in 3d, we have to account for faces that + // have non-standard face orientation, flip + // and rotation. thus, we have to store + // _eight_ data sets per face or subface + + // set up a table with the according offsets + // for non-standard orientation, first index: + // face_orientation (standard true=1), second + // index: face_flip (standard false=0), third + // index: face_rotation (standard false=0) + // + // note, that normally we should use the + // obvious offsets 0,1,2,3,4,5,6,7. However, + // prior to the changes enabling flipped and + // rotated faces, in many places of the + // library the convention was used, that the + // first dataset with offset 0 corresponds to + // a face in standard orientation. therefore + // we use the offsets 4,5,6,7,0,1,2,3 here to + // stick to that (implicit) convention + static const unsigned int offset[2][2][2] = { + {{4, 5}, // face_orientation=false; face_flip=false; + // face_rotation=false and true + {6, 7}}, // face_orientation=false; face_flip=true; + // face_rotation=false and true + {{0, 1}, // face_orientation=true; face_flip=false; + // face_rotation=false and true + {2, 3}}}; // face_orientation=true; face_flip=true; + // face_rotation=false and true + + + if (quadrature.size() == 1) + return (face_no + + offset[face_orientation][face_flip][face_rotation] * + GeometryInfo::faces_per_cell) * + quadrature[0].size(); + else + { + unsigned int n_points_i = 0; + for (unsigned int i = 0; i < face_no; ++i) + n_points_i += quadrature[i].size(); + + unsigned int n_points = 0; + for (unsigned int i = 0; i < GeometryInfo::faces_per_cell; + ++i) + n_points += quadrature[i].size(); + + return (n_points_i + + offset[face_orientation][face_flip][face_rotation] * + n_points); + } + } + + default: + Assert(false, ExcInternalError()); + } + return numbers::invalid_unsigned_int; +} + + + template <> QProjector<1>::DataSetDescriptor QProjector<1>::DataSetDescriptor::subface( diff --git a/source/fe/fe.cc b/source/fe/fe.cc index abd2bd4655..bc09e71a9d 100644 --- a/source/fe/fe.cc +++ b/source/fe/fe.cc @@ -1281,6 +1281,26 @@ FiniteElement::compute_n_nonzero_components( /*------------------------------- FiniteElement ----------------------*/ +#ifndef DOXYGEN +template +std::unique_ptr::InternalDataBase> +FiniteElement::get_face_data( + const UpdateFlags flags, + const Mapping & mapping, + const hp::QCollection &quadrature, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const +{ + return get_data(flags, + mapping, + QProjector::project_to_all_faces( + this->reference_cell_type(), quadrature), + output_data); +} + + + template std::unique_ptr::InternalDataBase> FiniteElement::get_face_data( @@ -1300,6 +1320,68 @@ FiniteElement::get_face_data( +template +inline void +FiniteElement::fill_fe_face_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const hp::QCollection & quadrature, + const Mapping & mapping, + const typename Mapping::InternalDataBase & mapping_internal, + const dealii::internal::FEValuesImplementation::MappingRelatedData + & mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const +{ + // base class version, implement overriden function in derived classes + AssertDimension(quadrature.size(), 1); + fill_fe_face_values(cell, + face_no, + quadrature[0], + mapping, + mapping_internal, + mapping_data, + fe_internal, + output_data); +} + + + +template +inline void +FiniteElement::fill_fe_face_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature & quadrature, + const Mapping & mapping, + const typename Mapping::InternalDataBase & mapping_internal, + const dealii::internal::FEValuesImplementation::MappingRelatedData + & mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const +{ + Assert(false, + ExcMessage("Use of a deprecated interface, please implement " + "fill_fe_face_values taking a hp::QCollection argument")); + (void)cell; + (void)face_no; + (void)quadrature; + (void)mapping; + (void)mapping_internal; + (void)mapping_data; + (void)fe_internal; + (void)output_data; +} +#endif + + + template std::unique_ptr::InternalDataBase> FiniteElement::get_subface_data( diff --git a/source/fe/fe_dgp_nonparametric.cc b/source/fe/fe_dgp_nonparametric.cc index 20dc3b1227..1505510170 100644 --- a/source/fe/fe_dgp_nonparametric.cc +++ b/source/fe/fe_dgp_nonparametric.cc @@ -362,7 +362,7 @@ void FE_DGPNonparametric::fill_fe_face_values( const typename Triangulation::cell_iterator &, const unsigned int, - const Quadrature &, + const hp::QCollection &, const Mapping &, const typename Mapping::InternalDataBase &, const dealii::internal::FEValuesImplementation::MappingRelatedData::setup_data( template std::unique_ptr::InternalDataBase> FE_Enriched::get_face_data( - const UpdateFlags update_flags, - const Mapping &mapping, - const Quadrature & quadrature, + const UpdateFlags update_flags, + const Mapping & mapping, + const hp::QCollection &quadrature, internal::FEValuesImplementation::FiniteElementRelatedData &output_data) const { + AssertDimension(quadrature.size(), 1); + auto data = fe_system->get_face_data(update_flags, mapping, quadrature, output_data); return setup_data(Utilities::dynamic_unique_cast< typename FESystem::InternalData>( std::move(data)), update_flags, - quadrature); + quadrature[0]); } @@ -558,7 +560,7 @@ void FE_Enriched::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase & mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData::fill_fe_face_values( ExcInternalError()); const InternalData &fe_data = static_cast(fe_internal); + AssertDimension(quadrature.size(), 1); + // call FESystem's method to fill everything without enrichment function fe_system->fill_fe_face_values(cell, face_no, @@ -584,7 +588,7 @@ FE_Enriched::fill_fe_face_values( if (is_enriched) multiply_by_enrichment( - quadrature, fe_data, mapping_data, cell, output_data); + quadrature[0], fe_data, mapping_data, cell, output_data); } diff --git a/source/fe/fe_face.cc b/source/fe/fe_face.cc index d673916566..648dbd6116 100644 --- a/source/fe/fe_face.cc +++ b/source/fe/fe_face.cc @@ -712,7 +712,7 @@ void FE_FaceQ<1, spacedim>::fill_fe_face_values( const typename Triangulation<1, spacedim>::cell_iterator &, const unsigned int face, - const Quadrature<0> &, + const hp::QCollection<0> &, const Mapping<1, spacedim> &, const typename Mapping<1, spacedim>::InternalDataBase &, const dealii::internal::FEValuesImplementation::MappingRelatedData<1, diff --git a/source/fe/fe_nedelec_sz.cc b/source/fe/fe_nedelec_sz.cc index f2a6e5d049..d319a6d184 100644 --- a/source/fe/fe_nedelec_sz.cc +++ b/source/fe/fe_nedelec_sz.cc @@ -2029,7 +2029,7 @@ void FE_NedelecSZ::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase & mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData @@ -2038,6 +2038,8 @@ FE_NedelecSZ::fill_fe_face_values( dealii::internal::FEValuesImplementation::FiniteElementRelatedData &data) const { + AssertDimension(quadrature.size(), 1); + // Note for future improvement: // We don't have the full quadrature - should use QProjector to create the 2D // quadrature. @@ -2060,18 +2062,18 @@ FE_NedelecSZ::fill_fe_face_values( // (fe_internal/fe_data) which was not filled in by get_data. fill_edge_values(cell, QProjector::project_to_all_faces( - this->reference_cell_type(), quadrature), + this->reference_cell_type(), quadrature[0]), fe_data); if (dim == 3 && this->degree > 1) { fill_face_values(cell, QProjector::project_to_all_faces( - this->reference_cell_type(), quadrature), + this->reference_cell_type(), quadrature[0]), fe_data); } const UpdateFlags flags(fe_data.update_each); - const unsigned int n_q_points = quadrature.size(); + const unsigned int n_q_points = quadrature[0].size(); const auto offset = QProjector::DataSetDescriptor::face(this->reference_cell_type(), face_no, diff --git a/source/fe/fe_nothing.cc b/source/fe/fe_nothing.cc index 721f34f8a8..e3a7b165d3 100644 --- a/source/fe/fe_nothing.cc +++ b/source/fe/fe_nothing.cc @@ -135,7 +135,7 @@ void FE_Nothing::fill_fe_face_values( const typename Triangulation::cell_iterator &, const unsigned int, - const Quadrature &, + const hp::QCollection &, const Mapping &, const typename Mapping::InternalDataBase &, const dealii::internal::FEValuesImplementation::MappingRelatedData::InternalDataBase> FE_P1NC::get_face_data( const UpdateFlags update_flags, const Mapping<2, 2> &, - const Quadrature<1> &quadrature, + const hp::QCollection<1> &quadrature, dealii::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2> &output_data) const { + AssertDimension(quadrature.size(), 1); + auto data_ptr = std::make_unique::InternalDataBase>(); data_ptr->update_each = requires_update_flags(update_flags); - const unsigned int n_q_points = quadrature.size(); + const unsigned int n_q_points = quadrature[0].size(); output_data.initialize(n_q_points, FE_P1NC(), data_ptr->update_each); // this is a linear element, so its second derivatives are zero @@ -246,7 +248,7 @@ void FE_P1NC::fill_fe_face_values( const Triangulation<2, 2>::cell_iterator &cell, const unsigned int face_no, - const Quadrature<1> & quadrature, + const hp::QCollection<1> & quadrature, const Mapping<2, 2> & mapping, const Mapping<2, 2>::InternalDataBase &, const dealii::internal::FEValuesImplementation::MappingRelatedData<2, 2> &, @@ -254,6 +256,8 @@ FE_P1NC::fill_fe_face_values( dealii::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2> &output_data) const { + AssertDimension(quadrature.size(), 1); + const UpdateFlags flags(fe_internal.update_each); // linear shape functions @@ -263,7 +267,7 @@ FE_P1NC::fill_fe_face_values( // compute on the face const Quadrature<2> quadrature_on_face = QProjector<2>::project_to_face(this->reference_cell_type(), - quadrature, + quadrature[0], face_no); if (flags & update_values) diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc index 6df444121b..c140e31765 100644 --- a/source/fe/fe_poly_tensor.cc +++ b/source/fe/fe_poly_tensor.cc @@ -945,7 +945,7 @@ void FE_PolyTensor::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase & mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData::fill_fe_face_values( spacedim> &output_data) const { + AssertDimension(quadrature.size(), 1); + // convert data object to internal // data for this class. fails with // an exception if that is not @@ -964,7 +966,7 @@ FE_PolyTensor::fill_fe_face_values( ExcInternalError()); const InternalData &fe_data = static_cast(fe_internal); - const unsigned int n_q_points = quadrature.size(); + const unsigned int n_q_points = quadrature[0].size(); // offset determines which data set // to take (all data sets for all // faces are stored contiguously) diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index de8169987e..791fcffea2 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -971,9 +971,9 @@ FESystem::get_data( template std::unique_ptr::InternalDataBase> FESystem::get_face_data( - const UpdateFlags flags, - const Mapping &mapping, - const Quadrature & quadrature, + const UpdateFlags flags, + const Mapping & mapping, + const hp::QCollection &quadrature, dealii::internal::FEValuesImplementation::FiniteElementRelatedData & /*output_data*/) const @@ -1007,7 +1007,7 @@ FESystem::get_face_data( internal::FEValuesImplementation::FiniteElementRelatedData &base_fe_output_object = data.get_fe_output_object(base_no); base_fe_output_object.initialize( - quadrature.size(), + quadrature.max_n_quadrature_points(), base_element(base_no), flags | base_element(base_no).requires_update_flags(flags)); @@ -1112,7 +1112,7 @@ FESystem::fill_fe_values( cell, invalid_face_number, invalid_face_number, - quadrature, + hp::QCollection(quadrature), cell_similarity, mapping_internal, fe_internal, @@ -1127,7 +1127,7 @@ void FESystem::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const Mapping & mapping, const typename Mapping::InternalDataBase & mapping_internal, const dealii::internal::FEValuesImplementation::MappingRelatedData::fill_fe_subface_values( cell, face_no, sub_no, - quadrature, + hp::QCollection(quadrature), CellSimilarity::none, mapping_internal, fe_internal, @@ -1191,7 +1191,7 @@ FESystem::compute_fill( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const CellSimilarity::Similarity cell_similarity, const typename Mapping::InternalDataBase & mapping_internal, const typename FiniteElement::InternalDataBase &fe_internal, @@ -1238,16 +1238,19 @@ FESystem::compute_fill( // fill_fe_face_values needs argument Quadrature for both cases // dim_1==dim-1 and dim_1=dim. Hence the following workaround - const Quadrature * cell_quadrature = nullptr; - const Quadrature *face_quadrature = nullptr; - const unsigned int n_q_points = quadrature.size(); + const Quadrature * cell_quadrature = nullptr; + const hp::QCollection *face_quadrature = nullptr; + const Quadrature * sub_face_quadrature = nullptr; + const unsigned int n_q_points = quadrature.size() == 1 ? + quadrature[0].size() : + quadrature[face_no].size(); // static cast to the common base class of quadrature being either // Quadrature or Quadrature: - const Subscriptor *quadrature_base_pointer = &quadrature; if (face_no == invalid_face_number) { + const Subscriptor *quadrature_base_pointer = &quadrature[0]; Assert(dim_1 == dim, ExcDimensionMismatch(dim_1, dim)); Assert(dynamic_cast *>( quadrature_base_pointer) != nullptr, @@ -1256,14 +1259,26 @@ FESystem::compute_fill( cell_quadrature = static_cast *>(quadrature_base_pointer); } + else if (sub_no == invalid_face_number) + { + const Subscriptor *quadrature_base_pointer = &quadrature; + Assert(dim_1 == dim - 1, ExcDimensionMismatch(dim_1, dim - 1)); + Assert(dynamic_cast *>( + quadrature_base_pointer) != nullptr, + ExcInternalError()); + + face_quadrature = static_cast *>( + quadrature_base_pointer); + } else { + const Subscriptor *quadrature_base_pointer = &quadrature[0]; Assert(dim_1 == dim - 1, ExcDimensionMismatch(dim_1, dim - 1)); Assert(dynamic_cast *>( quadrature_base_pointer) != nullptr, ExcInternalError()); - face_quadrature = + sub_face_quadrature = static_cast *>(quadrature_base_pointer); } @@ -1297,7 +1312,7 @@ FESystem::compute_fill( base_fe.fill_fe_subface_values(cell, face_no, sub_no, - *face_quadrature, + *sub_face_quadrature, mapping, mapping_internal, mapping_data, diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index ea659a9abf..4d4edceee7 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -444,10 +444,8 @@ namespace FEValuesViews &shape_function_data, std::vector::type> &values) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_values.n_cols() : values.size(); - AssertDimension(values.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = values.size(); std::fill(values.begin(), values.end(), @@ -488,10 +486,8 @@ namespace FEValuesViews typename ProductType>::type> &derivatives) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_derivatives[0].size() : derivatives.size(); - AssertDimension(derivatives.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = derivatives.size(); std::fill( derivatives.begin(), @@ -531,10 +527,8 @@ namespace FEValuesViews std::vector::template OutputType< Number>::laplacian_type> &laplacians) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_hessians[0].size() : laplacians.size(); - AssertDimension(laplacians.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = laplacians.size(); std::fill(laplacians.begin(), laplacians.end(), @@ -576,10 +570,8 @@ namespace FEValuesViews typename ProductType>::type> &values) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_values.n_cols() : values.size(); - AssertDimension(values.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = values.size(); std::fill( values.begin(), @@ -639,10 +631,8 @@ namespace FEValuesViews typename ProductType>::type> &derivatives) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_derivatives[0].size() : derivatives.size(); - AssertDimension(derivatives.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = derivatives.size(); std::fill( derivatives.begin(), @@ -708,10 +698,7 @@ namespace FEValuesViews &symmetric_gradients) { const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = dofs_per_cell > 0 ? - shape_gradients[0].size() : - symmetric_gradients.size(); - AssertDimension(symmetric_gradients.size(), n_quadrature_points); + const unsigned int n_quadrature_points = symmetric_gradients.size(); std::fill( symmetric_gradients.begin(), @@ -778,10 +765,8 @@ namespace FEValuesViews std::vector::template OutputType< Number>::divergence_type> &divergences) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_gradients[0].size() : divergences.size(); - AssertDimension(divergences.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = divergences.size(); std::fill(divergences.begin(), divergences.end(), @@ -843,10 +828,8 @@ namespace FEValuesViews Number, typename dealii::internal::CurlType::type>::type> &curls) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_gradients[0].size() : curls.size(); - AssertDimension(curls.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = curls.size(); std::fill(curls.begin(), curls.end(), @@ -1106,10 +1089,8 @@ namespace FEValuesViews std::vector::template OutputType< Number>::laplacian_type> &laplacians) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_hessians[0].size() : laplacians.size(); - AssertDimension(laplacians.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = laplacians.size(); std::fill(laplacians.begin(), laplacians.end(), @@ -1177,10 +1158,8 @@ namespace FEValuesViews dealii::SymmetricTensor<2, spacedim>>::type> &values) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_values.n_cols() : values.size(); - AssertDimension(values.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = values.size(); std::fill( values.begin(), @@ -1249,10 +1228,8 @@ namespace FEValuesViews std::vector:: template OutputType::divergence_type> &divergences) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_gradients[0].size() : divergences.size(); - AssertDimension(divergences.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = divergences.size(); std::fill(divergences.begin(), divergences.end(), @@ -1358,10 +1335,8 @@ namespace FEValuesViews typename ProductType>::type> &values) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_values.n_cols() : values.size(); - AssertDimension(values.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = values.size(); std::fill( values.begin(), @@ -1429,10 +1404,8 @@ namespace FEValuesViews std::vector::template OutputType< Number>::divergence_type> &divergences) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_gradients[0].size() : divergences.size(); - AssertDimension(divergences.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = divergences.size(); std::fill(divergences.begin(), divergences.end(), @@ -1500,10 +1473,8 @@ namespace FEValuesViews std::vector::template OutputType< Number>::gradient_type> &gradients) { - const unsigned int dofs_per_cell = dof_values.size(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_gradients[0].size() : gradients.size(); - AssertDimension(gradients.size(), n_quadrature_points); + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = gradients.size(); std::fill(gradients.begin(), gradients.end(), @@ -3076,6 +3047,7 @@ FEValuesBase::FEValuesBase( const Mapping & mapping, const FiniteElement &fe) : n_quadrature_points(n_q_points) + , max_n_quadrature_points(n_q_points) , dofs_per_cell(dofs_per_cell) , mapping(&mapping, typeid(*this).name()) , fe(&fe, typeid(*this).name()) @@ -3116,10 +3088,8 @@ namespace internal std::vector & values) { // scalar finite elements, so shape_values.size() == dofs_per_cell - const unsigned int dofs_per_cell = shape_values.n_rows(); - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_values.n_cols() : values.size(); - AssertDimension(values.size(), n_quadrature_points); + const unsigned int dofs_per_cell = shape_values.n_rows(); + const unsigned int n_quadrature_points = values.size(); // initialize with zero std::fill_n(values.begin(), @@ -3175,8 +3145,9 @@ namespace internal if (dofs_per_cell == 0) return; - const unsigned int n_quadrature_points = shape_values.n_cols(); - const unsigned int n_components = fe.n_components(); + const unsigned int n_quadrature_points = + quadrature_points_fastest ? values[0].size() : values.size(); + const unsigned int n_components = fe.n_components(); // Assert that we can write all components into the result vectors const unsigned result_components = n_components * component_multiple; @@ -3268,10 +3239,8 @@ namespace internal const dealii::Table<2, Tensor> &shape_derivatives, std::vector> & derivatives) { - const unsigned int dofs_per_cell = shape_derivatives.size()[0]; - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_derivatives[0].size() : derivatives.size(); - AssertDimension(derivatives.size(), n_quadrature_points); + const unsigned int dofs_per_cell = shape_derivatives.size()[0]; + const unsigned int n_quadrature_points = derivatives.size(); // initialize with zero std::fill_n(derivatives.begin(), @@ -3327,8 +3296,9 @@ namespace internal return; - const unsigned int n_quadrature_points = shape_derivatives[0].size(); - const unsigned int n_components = fe.n_components(); + const unsigned int n_quadrature_points = + quadrature_points_fastest ? derivatives[0].size() : derivatives.size(); + const unsigned int n_components = fe.n_components(); // Assert that we can write all components into the result vectors const unsigned result_components = n_components * component_multiple; @@ -3416,10 +3386,8 @@ namespace internal const dealii::Table<2, Tensor<2, spacedim>> &shape_hessians, std::vector & laplacians) { - const unsigned int dofs_per_cell = shape_hessians.size()[0]; - const unsigned int n_quadrature_points = - dofs_per_cell > 0 ? shape_hessians[0].size() : laplacians.size(); - AssertDimension(laplacians.size(), n_quadrature_points); + const unsigned int dofs_per_cell = shape_hessians.size()[0]; + const unsigned int n_quadrature_points = laplacians.size(); // initialize with zero std::fill_n(laplacians.begin(), @@ -3472,7 +3440,7 @@ namespace internal return; - const unsigned int n_quadrature_points = shape_hessians[0].size(); + const unsigned int n_quadrature_points = laplacians.size(); const unsigned int n_components = fe.n_components(); // Assert that we can write all components into the result vectors @@ -4209,6 +4177,7 @@ FEValuesBase::memory_consumption() const { return (sizeof(this->update_flags) + MemoryConsumption::memory_consumption(n_quadrature_points) + + MemoryConsumption::memory_consumption(max_n_quadrature_points) + sizeof(cell_similarity) + MemoryConsumption::memory_consumption(dofs_per_cell) + MemoryConsumption::memory_consumption(mapping) + @@ -4399,6 +4368,18 @@ FEValues::FEValues(const Mapping & mapping, +template +FEValues::FEValues(const Mapping & mapping, + const FiniteElement &fe, + const hp::QCollection & q, + const UpdateFlags update_flags) + : FEValues(mapping, fe, q[0], update_flags) +{ + AssertDimension(q.size(), 1); +} + + + template FEValues::FEValues(const FiniteElement &fe, const Quadrature & q, @@ -4415,6 +4396,17 @@ FEValues::FEValues(const FiniteElement &fe, +template +FEValues::FEValues(const FiniteElement &fe, + const hp::QCollection & q, + const UpdateFlags update_flags) + : FEValues(fe, q[0], update_flags) +{ + AssertDimension(q.size(), 1); +} + + + template void FEValues::initialize(const UpdateFlags update_flags) @@ -4433,8 +4425,8 @@ FEValues::initialize(const UpdateFlags update_flags) // initialize the base classes if (flags & update_mapping) - this->mapping_output.initialize(this->n_quadrature_points, flags); - this->finite_element_output.initialize(this->n_quadrature_points, + this->mapping_output.initialize(this->max_n_quadrature_points, flags); + this->finite_element_output.initialize(this->max_n_quadrature_points, *this->fe, flags); @@ -4598,20 +4590,41 @@ FEValues::memory_consumption() const template FEFaceValuesBase::FEFaceValuesBase( - const unsigned int n_q_points, + const unsigned int dofs_per_cell, + const UpdateFlags flags, + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature & quadrature) + : FEFaceValuesBase(dofs_per_cell, + flags, + mapping, + fe, + hp::QCollection(quadrature)) +{} + + + +template +FEFaceValuesBase::FEFaceValuesBase( const unsigned int dofs_per_cell, const UpdateFlags, const Mapping & mapping, const FiniteElement &fe, - const Quadrature & quadrature) - : FEValuesBase(n_q_points, + const hp::QCollection & quadrature) + : FEValuesBase(quadrature.max_n_quadrature_points(), dofs_per_cell, update_default, mapping, fe) , present_face_index(numbers::invalid_unsigned_int) , quadrature(quadrature) -{} +{ + Assert(quadrature.size() == 1 || + quadrature.size() == + ReferenceCell::internal::Info::get_cell(fe.reference_cell_type()) + .n_faces(), + ExcInternalError()); +} @@ -4654,8 +4667,21 @@ FEFaceValues::FEFaceValues( const FiniteElement &fe, const Quadrature & quadrature, const UpdateFlags update_flags) - : FEFaceValuesBase(quadrature.size(), - fe.n_dofs_per_cell(), + : FEFaceValues(mapping, + fe, + hp::QCollection(quadrature), + update_flags) +{} + + + +template +FEFaceValues::FEFaceValues( + const Mapping & mapping, + const FiniteElement &fe, + const hp::QCollection & quadrature, + const UpdateFlags update_flags) + : FEFaceValuesBase(fe.n_dofs_per_cell(), update_flags, mapping, fe, @@ -4671,8 +4697,19 @@ FEFaceValues::FEFaceValues( const FiniteElement &fe, const Quadrature & quadrature, const UpdateFlags update_flags) - : FEFaceValuesBase(quadrature.size(), - fe.n_dofs_per_cell(), + : FEFaceValues(fe, + hp::QCollection(quadrature), + update_flags) +{} + + + +template +FEFaceValues::FEFaceValues( + const FiniteElement &fe, + const hp::QCollection & quadrature, + const UpdateFlags update_flags) + : FEFaceValuesBase(fe.n_dofs_per_cell(), update_flags, StaticMappingQ1::mapping, fe, @@ -4691,27 +4728,42 @@ FEFaceValues::initialize(const UpdateFlags update_flags) // initialize the base classes if (flags & update_mapping) - this->mapping_output.initialize(this->n_quadrature_points, flags); - this->finite_element_output.initialize(this->n_quadrature_points, + this->mapping_output.initialize(this->max_n_quadrature_points, flags); + this->finite_element_output.initialize(this->max_n_quadrature_points, *this->fe, flags); // then get objects into which the FE and the Mapping can store // intermediate data used across calls to reinit. this can be done in parallel + + std::unique_ptr::InternalDataBase> ( + FiniteElement::*finite_element_get_face_data)( + const UpdateFlags, + const Mapping &, + const hp::QCollection &, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData + &) const = &FiniteElement::get_face_data; + + std::unique_ptr::InternalDataBase> ( + Mapping::*mapping_get_face_data)( + const UpdateFlags, const hp::QCollection &) const = + &Mapping::get_face_data; + + Threads::Task< std::unique_ptr::InternalDataBase>> - fe_get_data = - Threads::new_task(&FiniteElement::get_face_data, - *this->fe, - flags, - *this->mapping, - this->quadrature, - this->finite_element_output); + fe_get_data = Threads::new_task(finite_element_get_face_data, + *this->fe, + flags, + *this->mapping, + this->quadrature, + this->finite_element_output); Threads::Task< std::unique_ptr::InternalDataBase>> mapping_get_data; if (flags & update_mapping) - mapping_get_data = Threads::new_task(&Mapping::get_face_data, + mapping_get_data = Threads::new_task(mapping_get_face_data, *this->mapping, flags, this->quadrature); @@ -4831,6 +4883,9 @@ FEFaceValues::do_reinit(const unsigned int face_no) this->mapping_output, *this->fe_data, this->finite_element_output); + + const_cast(this->n_quadrature_points) = + this->quadrature[this->quadrature.size() == 1 ? 0 : face_no].size(); } @@ -4853,8 +4908,7 @@ FESubfaceValues::FESubfaceValues( const FiniteElement &fe, const Quadrature & quadrature, const UpdateFlags update_flags) - : FEFaceValuesBase(quadrature.size(), - fe.n_dofs_per_cell(), + : FEFaceValuesBase(fe.n_dofs_per_cell(), update_flags, mapping, fe, @@ -4865,13 +4919,25 @@ FESubfaceValues::FESubfaceValues( +template +FESubfaceValues::FESubfaceValues( + const Mapping & mapping, + const FiniteElement &fe, + const hp::QCollection & quadrature, + const UpdateFlags update_flags) + : FESubfaceValues(mapping, fe, quadrature[0], update_flags) +{ + AssertDimension(quadrature.size(), 1); +} + + + template FESubfaceValues::FESubfaceValues( const FiniteElement &fe, const Quadrature & quadrature, const UpdateFlags update_flags) - : FEFaceValuesBase(quadrature.size(), - fe.n_dofs_per_cell(), + : FEFaceValuesBase(fe.n_dofs_per_cell(), update_flags, StaticMappingQ1::mapping, fe, @@ -4882,6 +4948,18 @@ FESubfaceValues::FESubfaceValues( +template +FESubfaceValues::FESubfaceValues( + const FiniteElement &fe, + const hp::QCollection & quadrature, + const UpdateFlags update_flags) + : FESubfaceValues(fe, quadrature[0], update_flags) +{ + AssertDimension(quadrature.size(), 1); +} + + + template void FESubfaceValues::initialize(const UpdateFlags update_flags) @@ -4890,8 +4968,8 @@ FESubfaceValues::initialize(const UpdateFlags update_flags) // initialize the base classes if (flags & update_mapping) - this->mapping_output.initialize(this->n_quadrature_points, flags); - this->finite_element_output.initialize(this->n_quadrature_points, + this->mapping_output.initialize(this->max_n_quadrature_points, flags); + this->finite_element_output.initialize(this->max_n_quadrature_points, *this->fe, flags); @@ -4905,7 +4983,7 @@ FESubfaceValues::initialize(const UpdateFlags update_flags) *this->fe, flags, *this->mapping, - this->quadrature, + this->quadrature[0], this->finite_element_output); Threads::Task< std::unique_ptr::InternalDataBase>> @@ -4915,7 +4993,7 @@ FESubfaceValues::initialize(const UpdateFlags update_flags) Threads::new_task(&Mapping::get_subface_data, *this->mapping, flags, - this->quadrature); + this->quadrature[0]); this->update_flags = flags; @@ -5119,7 +5197,7 @@ FESubfaceValues::do_reinit(const unsigned int face_no, this->get_mapping().fill_fe_subface_values(*this->present_cell, face_no, subface_no, - this->quadrature, + this->quadrature[0], *this->mapping_data, this->mapping_output); } @@ -5127,7 +5205,7 @@ FESubfaceValues::do_reinit(const unsigned int face_no, this->get_fe().fill_fe_subface_values(*this->present_cell, face_no, subface_no, - this->quadrature, + this->quadrature[0], this->get_mapping(), *this->mapping_data, this->mapping_output, diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc index 23cbdb9aaa..288a9debcb 100644 --- a/source/fe/mapping.cc +++ b/source/fe/mapping.cc @@ -141,6 +141,76 @@ Mapping::project_real_point_to_unit_point_on_face( return {}; } + + +#ifndef DOXYGEN +template +void +Mapping::fill_fe_face_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const hp::QCollection & quadrature, + const typename Mapping::InternalDataBase & internal_data, + dealii::internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + // base class version, implement overriden function in derived classes + AssertDimension(quadrature.size(), 1); + fill_fe_face_values(cell, face_no, quadrature[0], internal_data, output_data); +} + + + +template +void +Mapping::fill_fe_face_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature & quadrature, + const typename Mapping::InternalDataBase & internal_data, + dealii::internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + Assert(false, + ExcMessage("Use of a deprecated interface, please implement " + "fill_fe_face_values taking a hp::QCollection argument")); + (void)cell; + (void)face_no; + (void)quadrature; + (void)internal_data; + (void)output_data; +} + + + +template +std::unique_ptr::InternalDataBase> +Mapping::get_face_data( + const UpdateFlags update_flags, + const hp::QCollection &quadrature) const +{ + // base class version, implement overriden function in derived classes + return get_face_data(update_flags, quadrature[0]); +} + + + +template +std::unique_ptr::InternalDataBase> +Mapping::get_face_data( + const UpdateFlags update_flags, + const Quadrature &quadrature) const +{ + Assert(false, + ExcMessage("Use of a deprecated interface, please implement " + "fill_fe_face_values taking a hp::QCollection argument")); + (void)update_flags; + (void)quadrature; + + return std::unique_ptr::InternalDataBase>(); +} +#endif + /* ---------------------------- InternalDataBase --------------------------- */ diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc index 2affce2eb5..745a44ebdb 100644 --- a/source/fe/mapping_cartesian.cc +++ b/source/fe/mapping_cartesian.cc @@ -111,13 +111,15 @@ MappingCartesian::get_data(const UpdateFlags update_flags, template std::unique_ptr::InternalDataBase> MappingCartesian::get_face_data( - const UpdateFlags update_flags, - const Quadrature &quadrature) const + const UpdateFlags update_flags, + const hp::QCollection &quadrature) const { + AssertDimension(quadrature.size(), 1); + std::unique_ptr::InternalDataBase> data_ptr = std::make_unique( QProjector::project_to_all_faces(ReferenceCell::get_hypercube(dim), - quadrature)); + quadrature[0])); auto &data = dynamic_cast(*data_ptr); // verify that we have computed the transitive hull of the required @@ -437,11 +439,13 @@ void MappingCartesian::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, internal::FEValuesImplementation::MappingRelatedData &output_data) const { + AssertDimension(quadrature.size(), 1); + // convert data object to internal // data for this class. fails with // an exception if that is not @@ -468,7 +472,7 @@ MappingCartesian::fill_fe_face_values( if (data.update_each & update_JxW_values) for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i) - output_data.JxW_values[i] = J * quadrature.weight(i); + output_data.JxW_values[i] = J * quadrature[0].weight(i); if (data.update_each & update_boundary_forms) for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i) diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index 16fa154641..a5c428ea4c 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -357,13 +357,13 @@ namespace internal maybe_compute_q_points( const typename QProjector::DataSetDescriptor data_set, const typename dealii::MappingFE::InternalData &data, - std::vector> &quadrature_points) + std::vector> &quadrature_points, + const unsigned int n_q_points) { const UpdateFlags update_flags = data.update_each; if (update_flags & update_quadrature_points) - for (unsigned int point = 0; point < quadrature_points.size(); - ++point) + for (unsigned int point = 0; point < n_q_points; ++point) { const double * shape = &data.shape(point + data_set, 0); Point result = @@ -390,7 +390,8 @@ namespace internal maybe_update_Jacobians( const CellSimilarity::Similarity cell_similarity, const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFE::InternalData &data) + const typename dealii::MappingFE::InternalData &data, + const unsigned int n_q_points) { const UpdateFlags update_flags = data.update_each; @@ -400,8 +401,6 @@ namespace internal // need to recompute jacobians... if (cell_similarity != CellSimilarity::translation) { - const unsigned int n_q_points = data.contravariant.size(); - std::fill(data.contravariant.begin(), data.contravariant.end(), DerivativeForm<1, dim, spacedim>()); @@ -442,7 +441,6 @@ namespace internal if (update_flags & update_covariant_transformation) if (cell_similarity != CellSimilarity::translation) { - const unsigned int n_q_points = data.contravariant.size(); for (unsigned int point = 0; point < n_q_points; ++point) { data.covariant[point] = @@ -453,7 +451,6 @@ namespace internal if (update_flags & update_volume_elements) if (cell_similarity != CellSimilarity::translation) { - const unsigned int n_q_points = data.contravariant.size(); for (unsigned int point = 0; point < n_q_points; ++point) data.volume_elements[point] = data.contravariant[point].determinant(); @@ -472,12 +469,13 @@ namespace internal const CellSimilarity::Similarity cell_similarity, const typename QProjector::DataSetDescriptor data_set, const typename dealii::MappingFE::InternalData &data, - std::vector> &jacobian_grads) + std::vector> &jacobian_grads, + const unsigned int n_q_points) { const UpdateFlags update_flags = data.update_each; if (update_flags & update_jacobian_grads) { - const unsigned int n_q_points = jacobian_grads.size(); + AssertIndexRange(n_q_points, jacobian_grads.size() + 1); if (cell_similarity != CellSimilarity::translation) for (unsigned int point = 0; point < n_q_points; ++point) @@ -518,13 +516,14 @@ namespace internal const CellSimilarity::Similarity cell_similarity, const typename QProjector::DataSetDescriptor data_set, const typename dealii::MappingFE::InternalData &data, - std::vector> &jacobian_pushed_forward_grads) + std::vector> &jacobian_pushed_forward_grads, + const unsigned int n_q_points) { const UpdateFlags update_flags = data.update_each; if (update_flags & update_jacobian_pushed_forward_grads) { - const unsigned int n_q_points = - jacobian_pushed_forward_grads.size(); + AssertIndexRange(n_q_points, + jacobian_pushed_forward_grads.size() + 1); if (cell_similarity != CellSimilarity::translation) { @@ -591,12 +590,13 @@ namespace internal const CellSimilarity::Similarity cell_similarity, const typename QProjector::DataSetDescriptor data_set, const typename dealii::MappingFE::InternalData &data, - std::vector> &jacobian_2nd_derivatives) + std::vector> &jacobian_2nd_derivatives, + const unsigned int n_q_points) { const UpdateFlags update_flags = data.update_each; if (update_flags & update_jacobian_2nd_derivatives) { - const unsigned int n_q_points = jacobian_2nd_derivatives.size(); + AssertIndexRange(n_q_points, jacobian_2nd_derivatives.size() + 1); if (cell_similarity != CellSimilarity::translation) { @@ -646,13 +646,15 @@ namespace internal const typename QProjector::DataSetDescriptor data_set, const typename dealii::MappingFE::InternalData &data, std::vector> - &jacobian_pushed_forward_2nd_derivatives) + & jacobian_pushed_forward_2nd_derivatives, + const unsigned int n_q_points) { const UpdateFlags update_flags = data.update_each; if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) { - const unsigned int n_q_points = - jacobian_pushed_forward_2nd_derivatives.size(); + AssertIndexRange(n_q_points, + jacobian_pushed_forward_2nd_derivatives.size() + + 1); if (cell_similarity != CellSimilarity::translation) { @@ -747,12 +749,13 @@ namespace internal const CellSimilarity::Similarity cell_similarity, const typename QProjector::DataSetDescriptor data_set, const typename dealii::MappingFE::InternalData &data, - std::vector> &jacobian_3rd_derivatives) + std::vector> &jacobian_3rd_derivatives, + const unsigned int n_q_points) { const UpdateFlags update_flags = data.update_each; if (update_flags & update_jacobian_3rd_derivatives) { - const unsigned int n_q_points = jacobian_3rd_derivatives.size(); + AssertIndexRange(n_q_points, jacobian_3rd_derivatives.size() + 1); if (cell_similarity != CellSimilarity::translation) { @@ -805,13 +808,15 @@ namespace internal const typename QProjector::DataSetDescriptor data_set, const typename dealii::MappingFE::InternalData &data, std::vector> - &jacobian_pushed_forward_3rd_derivatives) + & jacobian_pushed_forward_3rd_derivatives, + const unsigned int n_q_points) { const UpdateFlags update_flags = data.update_each; if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) { - const unsigned int n_q_points = - jacobian_pushed_forward_3rd_derivatives.size(); + AssertIndexRange(n_q_points, + jacobian_pushed_forward_3rd_derivatives.size() + + 1); if (cell_similarity != CellSimilarity::translation) { @@ -1067,8 +1072,8 @@ MappingFE::get_data(const UpdateFlags update_flags, template std::unique_ptr::InternalDataBase> MappingFE::get_face_data( - const UpdateFlags update_flags, - const Quadrature &quadrature) const + const UpdateFlags update_flags, + const hp::QCollection &quadrature) const { std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(*this->fe); @@ -1076,7 +1081,7 @@ MappingFE::get_face_data( data.initialize_face(this->requires_update_flags(update_flags), QProjector::project_to_all_faces( this->fe->reference_cell_type(), quadrature), - quadrature.size()); + quadrature.max_n_quadrature_points()); return data_ptr; } @@ -1142,51 +1147,61 @@ MappingFE::fill_fe_values( internal::MappingFEImplementation::maybe_compute_q_points( QProjector::DataSetDescriptor::cell(), data, - output_data.quadrature_points); + output_data.quadrature_points, + n_q_points); internal::MappingFEImplementation::maybe_update_Jacobians( - computed_cell_similarity, QProjector::DataSetDescriptor::cell(), data); + computed_cell_similarity, + QProjector::DataSetDescriptor::cell(), + data, + n_q_points); internal::MappingFEImplementation::maybe_update_jacobian_grads( computed_cell_similarity, QProjector::DataSetDescriptor::cell(), data, - output_data.jacobian_grads); + output_data.jacobian_grads, + n_q_points); internal::MappingFEImplementation::maybe_update_jacobian_pushed_forward_grads< dim, spacedim>(computed_cell_similarity, QProjector::DataSetDescriptor::cell(), data, - output_data.jacobian_pushed_forward_grads); + output_data.jacobian_pushed_forward_grads, + n_q_points); internal::MappingFEImplementation::maybe_update_jacobian_2nd_derivatives< dim, spacedim>(computed_cell_similarity, QProjector::DataSetDescriptor::cell(), data, - output_data.jacobian_2nd_derivatives); + output_data.jacobian_2nd_derivatives, + n_q_points); internal::MappingFEImplementation:: maybe_update_jacobian_pushed_forward_2nd_derivatives( computed_cell_similarity, QProjector::DataSetDescriptor::cell(), data, - output_data.jacobian_pushed_forward_2nd_derivatives); + output_data.jacobian_pushed_forward_2nd_derivatives, + n_q_points); internal::MappingFEImplementation::maybe_update_jacobian_3rd_derivatives< dim, spacedim>(computed_cell_similarity, QProjector::DataSetDescriptor::cell(), data, - output_data.jacobian_3rd_derivatives); + output_data.jacobian_3rd_derivatives, + n_q_points); internal::MappingFEImplementation:: maybe_update_jacobian_pushed_forward_3rd_derivatives( computed_cell_similarity, QProjector::DataSetDescriptor::cell(), data, - output_data.jacobian_pushed_forward_3rd_derivatives); + output_data.jacobian_pushed_forward_3rd_derivatives, + n_q_points); const UpdateFlags update_flags = data.update_each; const std::vector &weights = quadrature.get_weights(); @@ -1343,11 +1358,13 @@ namespace internal update_JxW_values | update_inverse_jacobians)) { if (update_flags & update_boundary_forms) - AssertDimension(output_data.boundary_forms.size(), n_q_points); + AssertIndexRange(n_q_points, + output_data.boundary_forms.size() + 1); if (update_flags & update_normal_vectors) - AssertDimension(output_data.normal_vectors.size(), n_q_points); + AssertIndexRange(n_q_points, + output_data.normal_vectors.size() + 1); if (update_flags & update_JxW_values) - AssertDimension(output_data.JxW_values.size(), n_q_points); + AssertIndexRange(n_q_points, output_data.JxW_values.size() + 1); Assert(data.aux.size() + 1 >= dim, ExcInternalError()); @@ -1413,7 +1430,7 @@ namespace internal // // to compute the cell normal, use the same method used in // fill_fe_values for cells above - AssertDimension(data.contravariant.size(), n_q_points); + AssertIndexRange(n_q_points, data.contravariant.size() + 1); for (unsigned int point = 0; point < n_q_points; ++point) { @@ -1446,8 +1463,7 @@ namespace internal } if (update_flags & update_JxW_values) - for (unsigned int i = 0; i < output_data.boundary_forms.size(); - ++i) + for (unsigned int i = 0; i < n_q_points; ++i) { output_data.JxW_values[i] = output_data.boundary_forms[i].norm() * @@ -1467,8 +1483,7 @@ namespace internal } if (update_flags & update_normal_vectors) - for (unsigned int i = 0; i < output_data.normal_vectors.size(); - ++i) + for (unsigned int i = 0; i < n_q_points; ++i) output_data.normal_vectors[i] = Point(output_data.boundary_forms[i] / output_data.boundary_forms[i].norm()); @@ -1505,47 +1520,57 @@ namespace internal internal::FEValuesImplementation::MappingRelatedData &output_data) { + const unsigned int n_q_points = quadrature.size(); + maybe_compute_q_points(data_set, data, - output_data.quadrature_points); + output_data.quadrature_points, + n_q_points); maybe_update_Jacobians(CellSimilarity::none, data_set, - data); + data, + n_q_points); maybe_update_jacobian_grads(CellSimilarity::none, data_set, data, - output_data.jacobian_grads); + output_data.jacobian_grads, + n_q_points); maybe_update_jacobian_pushed_forward_grads( CellSimilarity::none, data_set, data, - output_data.jacobian_pushed_forward_grads); + output_data.jacobian_pushed_forward_grads, + n_q_points); maybe_update_jacobian_2nd_derivatives( CellSimilarity::none, data_set, data, - output_data.jacobian_2nd_derivatives); + output_data.jacobian_2nd_derivatives, + n_q_points); maybe_update_jacobian_pushed_forward_2nd_derivatives( CellSimilarity::none, data_set, data, - output_data.jacobian_pushed_forward_2nd_derivatives); + output_data.jacobian_pushed_forward_2nd_derivatives, + n_q_points); maybe_update_jacobian_3rd_derivatives( CellSimilarity::none, data_set, data, - output_data.jacobian_3rd_derivatives); + output_data.jacobian_3rd_derivatives, + n_q_points); maybe_update_jacobian_pushed_forward_3rd_derivatives( CellSimilarity::none, data_set, data, - output_data.jacobian_pushed_forward_3rd_derivatives); + output_data.jacobian_pushed_forward_3rd_derivatives, + n_q_points); maybe_compute_face_data(mapping, cell, face_no, subface_no, - quadrature.size(), + n_q_points, data_set, data, output_data); @@ -1561,7 +1586,7 @@ void MappingFE::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, internal::FEValuesImplementation::MappingRelatedData &output_data) const @@ -1594,8 +1619,8 @@ MappingFE::fill_fe_face_values( cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), - quadrature.size()), - quadrature, + quadrature), + quadrature[quadrature.size() == 1 ? 0 : face_no], data, output_data); } diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index 2da4c2ff17..31470830ed 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -628,16 +628,18 @@ MappingFEField::get_data( template std::unique_ptr::InternalDataBase> MappingFEField::get_face_data( - const UpdateFlags update_flags, - const Quadrature &quadrature) const + const UpdateFlags update_flags, + const hp::QCollection &quadrature) const { + AssertDimension(quadrature.size(), 1); + std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(euler_dof_handler->get_fe(), fe_mask); auto & data = dynamic_cast(*data_ptr); const Quadrature q( QProjector::project_to_all_faces(ReferenceCell::get_hypercube(dim), - quadrature)); - this->compute_face_data(update_flags, q, quadrature.size(), data); + quadrature[0])); + this->compute_face_data(update_flags, q, quadrature[0].size(), data); return data_ptr; } @@ -1694,11 +1696,13 @@ void MappingFEField::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, internal::FEValuesImplementation::MappingRelatedData &output_data) const { + AssertDimension(quadrature.size(), 1); + // convert data object to internal data for this class. fails with an // exception if that is not possible Assert(dynamic_cast(&internal_data) != nullptr, @@ -1719,8 +1723,8 @@ MappingFEField::fill_fe_face_values( cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), - quadrature.size()), - quadrature, + quadrature[0].size()), + quadrature[0], data, euler_dof_handler->get_fe(), fe_mask, diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index b8dc584573..fcf99b00e9 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -282,16 +282,18 @@ MappingManifold::get_data(const UpdateFlags update_flags, template std::unique_ptr::InternalDataBase> MappingManifold::get_face_data( - const UpdateFlags update_flags, - const Quadrature &quadrature) const + const UpdateFlags update_flags, + const hp::QCollection &quadrature) const { + AssertDimension(quadrature.size(), 1); + std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(); auto &data = dynamic_cast(*data_ptr); data.initialize_face(this->requires_update_flags(update_flags), QProjector::project_to_all_faces( - ReferenceCell::get_hypercube(dim), quadrature), - quadrature.size()); + ReferenceCell::get_hypercube(dim), quadrature[0]), + quadrature[0].size()); return data_ptr; } @@ -1208,11 +1210,13 @@ void MappingManifold::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, internal::FEValuesImplementation::MappingRelatedData &output_data) const { + AssertDimension(quadrature.size(), 1); + // ensure that the following cast is really correct: Assert((dynamic_cast(&internal_data) != nullptr), ExcInternalError()); @@ -1228,8 +1232,8 @@ MappingManifold::fill_fe_face_values( cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), - quadrature.size()), - quadrature, + quadrature[0].size()), + quadrature[0], data, output_data); } diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 1f1112939a..69276944a3 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -179,21 +179,25 @@ MappingQ::get_data(const UpdateFlags update_flags, template std::unique_ptr::InternalDataBase> MappingQ::get_face_data( - const UpdateFlags update_flags, - const Quadrature &quadrature) const + const UpdateFlags update_flags, + const hp::QCollection &quadrature) const { std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(); auto &data = dynamic_cast(*data_ptr); + std::unique_ptr::InternalDataBase> ( + MappingQGeneric::*mapping_get_face_data)( + const UpdateFlags, const hp::QCollection &) const = + &MappingQGeneric::get_face_data; + // build the Q1 and Qp internal data objects in parallel Threads::Task< std::unique_ptr::InternalDataBase>> - do_get_data = - Threads::new_task(&MappingQGeneric::get_face_data, - *qp_mapping, - update_flags, - quadrature); + do_get_data = Threads::new_task(mapping_get_face_data, + *qp_mapping, + update_flags, + quadrature); if (!use_mapping_q_on_all_cells) data.mapping_q1_data = Utilities::dynamic_unique_cast< @@ -302,7 +306,7 @@ void MappingQ::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, internal::FEValuesImplementation::MappingRelatedData &output_data) const diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index f9c772999a..c481bad760 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -878,16 +878,18 @@ MappingQGeneric::get_data(const UpdateFlags update_flags, template std::unique_ptr::InternalDataBase> MappingQGeneric::get_face_data( - const UpdateFlags update_flags, - const Quadrature &quadrature) const + const UpdateFlags update_flags, + const hp::QCollection &quadrature) const { + AssertDimension(quadrature.size(), 1); + std::unique_ptr::InternalDataBase> data_ptr = std::make_unique(polynomial_degree); auto &data = dynamic_cast(*data_ptr); data.initialize_face(this->requires_update_flags(update_flags), QProjector::project_to_all_faces( - ReferenceCell::get_hypercube(dim), quadrature), - quadrature.size()); + ReferenceCell::get_hypercube(dim), quadrature[0]), + quadrature[0].size()); return data_ptr; } @@ -1140,11 +1142,13 @@ void MappingQGeneric::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, - const Quadrature & quadrature, + const hp::QCollection & quadrature, const typename Mapping::InternalDataBase & internal_data, internal::FEValuesImplementation::MappingRelatedData &output_data) const { + AssertDimension(quadrature.size(), 1); + // ensure that the following cast is really correct: Assert((dynamic_cast(&internal_data) != nullptr), ExcInternalError()); @@ -1173,8 +1177,8 @@ MappingQGeneric::fill_fe_face_values( cell->face_orientation(face_no), cell->face_flip(face_no), cell->face_rotation(face_no), - quadrature.size()), - quadrature, + quadrature[0].size()), + quadrature[0], data, output_data); } diff --git a/source/hp/fe_values.cc b/source/hp/fe_values.cc index e15ce49395..4d728b63a5 100644 --- a/source/hp/fe_values.cc +++ b/source/hp/fe_values.cc @@ -25,6 +25,41 @@ DEAL_II_NAMESPACE_OPEN namespace hp { + namespace + { + /** + * Given a quadrature collection, create a vector of quadrature collections + * with each containing one quadrature rule from the input collection. + */ + template + std::vector> + translate(const QCollection &q_collection) + { + std::vector> q_collections; + + for (unsigned int q = 0; q < q_collection.size(); ++q) + q_collections.emplace_back(q_collection[q]); + + return q_collections; + } + + /** + * Given a vector of quadrature collections, extract the first quadrature + * rule of each collection and construct with these a new collection. + */ + template + QCollection + translate(const std::vector> &q_collections) + { + QCollection result; + + for (unsigned int q = 0; q < q_collections.size(); ++q) + result.push_back(q_collections[q][0]); + + return result; + } + } // namespace + // -------------------------- FEValuesBase ------------------------- template @@ -37,6 +72,7 @@ namespace hp : fe_collection(&fe_collection) , mapping_collection(&mapping_collection) , q_collection(q_collection) + , q_collections(translate(q_collection)) , fe_values_table(fe_collection.size(), mapping_collection.size(), q_collection.size()) @@ -46,18 +82,20 @@ namespace hp , update_flags(update_flags) {} - template FEValuesBase::FEValuesBase( + const MappingCollection + & mapping_collection, const FECollection &fe_collection, - const QCollection & q_collection, + const std::vector> & q_collections, const UpdateFlags update_flags) : fe_collection(&fe_collection) - , mapping_collection( - &dealii::hp::StaticMappingQ1:: - mapping_collection) - , q_collection(q_collection) - , fe_values_table(fe_collection.size(), 1, q_collection.size()) + , mapping_collection(&mapping_collection) + , q_collection(translate(q_collections)) + , q_collections(q_collections) + , fe_values_table(fe_collection.size(), + mapping_collection.size(), + q_collection.size()) , present_fe_values_index(numbers::invalid_unsigned_int, numbers::invalid_unsigned_int, numbers::invalid_unsigned_int) @@ -65,6 +103,34 @@ namespace hp {} + template + FEValuesBase::FEValuesBase( + const FECollection &fe_collection, + const QCollection & q_collection, + const UpdateFlags update_flags) + : FEValuesBase( + dealii::hp::StaticMappingQ1:: + mapping_collection, + fe_collection, + q_collection, + update_flags) + {} + + + template + FEValuesBase::FEValuesBase( + const FECollection &fe_collection, + const std::vector> & q_collection, + const UpdateFlags update_flags) + : FEValuesBase( + dealii::hp::StaticMappingQ1:: + mapping_collection, + fe_collection, + q_collection, + update_flags) + {} + + template FEValuesBase::FEValuesBase( @@ -72,9 +138,10 @@ namespace hp : fe_collection(other.fe_collection) , mapping_collection(other.mapping_collection) , q_collection(other.q_collection) - , fe_values_table(fe_collection->size(), - mapping_collection->size(), - q_collection.size()) + , q_collections(other.q_collections) + , fe_values_table(other.fe_values_table.size(0), + other.fe_values_table.size(1), + other.fe_values_table.size(2)) , present_fe_values_index(other.present_fe_values_index) , update_flags(other.update_flags) { @@ -82,18 +149,19 @@ namespace hp // now it just contains nullptrs. Create copies of the objects that // `other.fe_values_table` stores Threads::TaskGroup<> task_group; - for (unsigned int fe_index = 0; fe_index < fe_collection->size(); + for (unsigned int fe_index = 0; fe_index < other.fe_values_table.size(0); ++fe_index) - for (unsigned int m_index = 0; m_index < mapping_collection->size(); + for (unsigned int m_index = 0; m_index < other.fe_values_table.size(1); ++m_index) - for (unsigned int q_index = 0; q_index < q_collection.size(); ++q_index) + for (unsigned int q_index = 0; q_index < other.fe_values_table.size(2); + ++q_index) if (other.fe_values_table[fe_index][m_index][q_index].get() != nullptr) task_group += Threads::new_task([&, fe_index, m_index, q_index]() { fe_values_table[fe_index][m_index][q_index] = std::make_unique((*mapping_collection)[m_index], (*fe_collection)[fe_index], - q_collection[q_index], + q_collections[q_index], update_flags); }); @@ -111,7 +179,7 @@ namespace hp { AssertIndexRange(fe_index, fe_collection->size()); AssertIndexRange(mapping_index, mapping_collection->size()); - AssertIndexRange(q_index, q_collection.size()); + AssertIndexRange(q_index, q_collections.size()); // set the triple of indices @@ -126,7 +194,7 @@ namespace hp fe_values_table(present_fe_values_index) = std::make_unique((*mapping_collection)[mapping_index], (*fe_collection)[fe_index], - q_collection[q_index], + q_collections[q_index], update_flags); // now there definitely is one! @@ -154,7 +222,7 @@ namespace hp AssertIndexRange(fe_index, fe_collection->size()); AssertIndexRange(mapping_index, mapping_collection->size()); - AssertIndexRange(q_index, q_collection.size()); + AssertIndexRange(q_index, q_collections.size()); task_group += Threads::new_task([&, fe_index, mapping_index, q_index]() { @@ -162,7 +230,7 @@ namespace hp std::make_unique( (*mapping_collection)[mapping_index], (*fe_collection)[fe_index], - q_collection[q_index], + q_collections[q_index], update_flags); }); } @@ -186,7 +254,7 @@ namespace hp indices : std::vector(size, 0), /*q_indices=*/ - (q_collection.size() > 1) ? + (q_collections.size() > 1) ? indices : std::vector(size, 0)); } @@ -238,7 +306,7 @@ namespace hp if (real_q_index == numbers::invalid_unsigned_int) { - if (this->q_collection.size() > 1) + if (this->q_collections.size() > 1) real_q_index = cell->active_fe_index(); else real_q_index = 0; @@ -256,7 +324,7 @@ namespace hp real_fe_index = cell->active_fe_index(); // some checks - AssertIndexRange(real_q_index, this->q_collection.size()); + AssertIndexRange(real_q_index, this->q_collections.size()); AssertIndexRange(real_mapping_index, this->mapping_collection->size()); AssertIndexRange(real_fe_index, this->fe_collection->size()); @@ -292,7 +360,7 @@ namespace hp real_fe_index = 0; // some checks - AssertIndexRange(real_q_index, this->q_collection.size()); + AssertIndexRange(real_q_index, this->q_collections.size()); AssertIndexRange(real_mapping_index, this->mapping_collection->size()); AssertIndexRange(real_fe_index, this->fe_collection->size()); @@ -320,6 +388,19 @@ namespace hp update_flags) {} + template + FEFaceValues::FEFaceValues( + const hp::MappingCollection & mapping, + const hp::FECollection & fe_collection, + const std::vector> &q_collection, + const UpdateFlags update_flags) + : hp::FEValuesBase>( + mapping, + fe_collection, + q_collection, + update_flags) + {} + template FEFaceValues::FEFaceValues( @@ -333,6 +414,18 @@ namespace hp {} + template + FEFaceValues::FEFaceValues( + const hp::FECollection & fe_collection, + const std::vector> &q_collection, + const UpdateFlags update_flags) + : hp::FEValuesBase>( + fe_collection, + q_collection, + update_flags) + {} + + template template void @@ -350,7 +443,7 @@ namespace hp if (real_q_index == numbers::invalid_unsigned_int) { - if (this->q_collection.size() > 1) + if (this->q_collections.size() > 1) real_q_index = cell->active_fe_index(); else real_q_index = 0; @@ -368,7 +461,7 @@ namespace hp real_fe_index = cell->active_fe_index(); // some checks - AssertIndexRange(real_q_index, this->q_collection.size()); + AssertIndexRange(real_q_index, this->q_collections.size()); AssertIndexRange(real_mapping_index, this->mapping_collection->size()); AssertIndexRange(real_fe_index, this->fe_collection->size()); @@ -421,7 +514,7 @@ namespace hp real_fe_index = 0; // some checks - AssertIndexRange(real_q_index, this->q_collection.size()); + AssertIndexRange(real_q_index, this->q_collections.size()); AssertIndexRange(real_mapping_index, this->mapping_collection->size()); AssertIndexRange(real_fe_index, this->fe_collection->size()); @@ -495,7 +588,7 @@ namespace hp if (real_q_index == numbers::invalid_unsigned_int) { - if (this->q_collection.size() > 1) + if (this->q_collections.size() > 1) real_q_index = cell->active_fe_index(); else real_q_index = 0; @@ -513,7 +606,7 @@ namespace hp real_fe_index = cell->active_fe_index(); // some checks - AssertIndexRange(real_q_index, this->q_collection.size()); + AssertIndexRange(real_q_index, this->q_collections.size()); AssertIndexRange(real_mapping_index, this->mapping_collection->size()); AssertIndexRange(real_fe_index, this->fe_collection->size()); @@ -551,7 +644,7 @@ namespace hp real_fe_index = 0; // some checks - AssertIndexRange(real_q_index, this->q_collection.size()); + AssertIndexRange(real_q_index, this->q_collections.size()); AssertIndexRange(real_mapping_index, this->mapping_collection->size()); AssertIndexRange(real_fe_index, this->fe_collection->size()); diff --git a/tests/simplex/variable_face_quadratures_01.cc b/tests/simplex/variable_face_quadratures_01.cc new file mode 100644 index 0000000000..37789de57d --- /dev/null +++ b/tests/simplex/variable_face_quadratures_01.cc @@ -0,0 +1,321 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test FEFaceValues for different face quadrature rules. + + +#include +#include + +#include +#include +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +class Fu : public Function +{ +public: + Fu(const unsigned int component = numbers::invalid_unsigned_int) + : Function(component == numbers::invalid_unsigned_int ? dim : 1) + , component(component) + {} + + virtual double + value(const Point &p, const unsigned int i = 0) const + { + if (component == numbers::invalid_unsigned_int) + return p[i]; + else + return p[component]; + } + + const unsigned int component; +}; + +template +void +test() +{} + +template <> +void +test<2>() +{ + const unsigned int dim = 2; + + // test: QProjector::project_to_all_faces + { + const hp::QCollection quad_ref(QGauss(1), + QGauss(2), + QGauss(3), + QGauss(4)); + + const auto quad = + QProjector::project_to_all_faces(ReferenceCell::Type::Quad, + quad_ref); + + const auto print = [&](const unsigned int face_no) { + deallog << "face_no=" << face_no << ":" << std::endl; + for (unsigned int q = 0, + i = QProjector::DataSetDescriptor::face( + ReferenceCell::Type::Quad, + face_no, + false, + false, + false, + quad_ref); + q < quad_ref[face_no].size(); + ++q, ++i) + { + deallog << quad.point(i) << " "; + deallog << quad.weight(i) << " "; + deallog << std::endl; + } + deallog << std::endl; + }; + + for (unsigned int i = 0; i < 4 /*TODO*/; ++i) + print(i); + + deallog << std::endl; + } + + // test: Mapping + { + const hp::QCollection quad_ref(QGauss(1), + QGauss(2), + QGauss(3), + QGauss(4)); + + MappingFE mapping(FE_Q(1)); + FE_Q fe(3); + + const UpdateFlags flags = mapping.requires_update_flags( + update_values | update_quadrature_points | update_JxW_values); + + auto data_ref = mapping.get_face_data(flags, quad_ref); + + internal::FEValuesImplementation::MappingRelatedData data; + data.initialize(quad_ref.max_n_quadrature_points(), flags); + + Triangulation tria; + GridGenerator::hyper_cube(tria); + + for (const auto &cell : tria.active_cell_iterators()) + for (const auto face_no : cell->face_indices()) + { + mapping.fill_fe_face_values(cell, face_no, quad_ref, *data_ref, data); + + deallog << "face_no=" << face_no << ":" << std::endl; + for (unsigned int q = 0; q < quad_ref[face_no].size(); ++q) + { + deallog << data.quadrature_points[q] << " "; + deallog << data.JxW_values[q] << " "; + deallog << std::endl; + } + deallog << std::endl; + } + +#if false + // not possible to test since functions are protected + internal::FEValuesImplementation::FiniteElementRelatedData data_fe; + data_fe.initialize(quad_ref.max_n_quadrature_points(), fe, flags); + + auto data_fe_ref =fe.get_face_data(flags, mapping, quad_ref, data_fe); + + for (const auto &cell : tria.active_cell_iterators()) + for (const auto face_no : cell->face_indices()) + { + fe.fill_fe_face_values(cell, face_no, quad_ref, mapping, *data_ref, data, *data_fe_ref, data_fe); + } +#endif + } + + // test FEFaceValues for FE_Q + { + const hp::QCollection quad_ref(QGauss(1), + QGauss(2), + QGauss(3), + QGauss(4)); + + MappingFE mapping(FE_Q(1)); + FE_Q fe(3); + + const UpdateFlags flags = mapping.requires_update_flags( + update_values | update_quadrature_points | update_JxW_values); + + + FEFaceValues fe_face_values(mapping, fe, quad_ref, flags); + + Triangulation tria; + GridGenerator::hyper_cube(tria); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector vector_0(dof_handler.n_dofs()); + Vector vector_1(dof_handler.n_dofs()); + + VectorTools::interpolate(mapping, dof_handler, Fu(0), vector_0); + VectorTools::interpolate(mapping, dof_handler, Fu(1), vector_1); + + std::vector values_0, values_1; + + values_0.reserve(fe_face_values.n_quadrature_points); + values_1.reserve(fe_face_values.n_quadrature_points); + + for (const auto &cell : dof_handler.active_cell_iterators()) + for (const auto face_no : cell->face_indices()) + { + fe_face_values.reinit(cell, face_no); + + values_0.resize(fe_face_values.n_quadrature_points); + values_1.resize(fe_face_values.n_quadrature_points); + + fe_face_values.get_function_values(vector_0, values_0); + fe_face_values.get_function_values(vector_1, values_1); + + for (unsigned int q = 0; q < values_0.size(); ++q) + deallog << values_0[q] << " " << values_1[q] << " " << std::endl; + + deallog << std::endl; + } + } + + // test FEFaceValues for FE_System(FE_Q) + { + const hp::QCollection quad_ref(QGauss(1), + QGauss(2), + QGauss(3), + QGauss(4)); + + MappingFE mapping(FE_Q(1)); + FESystem fe(FE_Q{3}, dim); + + const UpdateFlags flags = mapping.requires_update_flags( + update_values | update_quadrature_points | update_JxW_values); + + + FEFaceValues fe_face_values(mapping, fe, quad_ref, flags); + + Triangulation tria; + GridGenerator::hyper_cube(tria); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector vector_0(dof_handler.n_dofs()); + + VectorTools::interpolate(mapping, dof_handler, Fu(), vector_0); + + for (const auto &cell : dof_handler.active_cell_iterators()) + for (const auto face_no : cell->face_indices()) + { + fe_face_values.reinit(cell, face_no); + + std::vector> values_0(quad_ref[face_no].size(), + Vector(dim)); + fe_face_values.get_function_values(vector_0, values_0); + + deallog << "face_no=" << face_no << ":" << std::endl; + + for (unsigned int q = 0; q < values_0.size(); ++q) + deallog << values_0[q][0] << " " << values_0[q][1] << " " + << std::endl; + + deallog << std::endl; + } + } +} + +template <> +void +test<3>() +{ + const unsigned int dim = 3; + + // test FEFaceValues for FE_System(FE_Q) + { + const hp::QCollection quad_ref(QGauss(1), + QGauss(2), + QGauss(3), + QGauss(4), + QGauss(1), + QGauss(2)); + + MappingFE mapping(FE_Q(1)); + FESystem fe(FE_Q{3}, dim); + + const UpdateFlags flags = + update_values | update_quadrature_points | update_JxW_values; + + + FEFaceValues fe_face_values(mapping, fe, quad_ref, flags); + + Triangulation tria; + GridGenerator::hyper_cube(tria); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector vector_0(dof_handler.n_dofs()); + + VectorTools::interpolate(mapping, dof_handler, Fu(), vector_0); + + for (const auto &cell : dof_handler.active_cell_iterators()) + for (const auto face_no : cell->face_indices()) + { + fe_face_values.reinit(cell, face_no); + + std::vector> values_0(quad_ref[face_no].size(), + Vector(dim)); + fe_face_values.get_function_values(vector_0, values_0); + + deallog << "face_no=" << face_no << ":" << std::endl; + + for (unsigned int q = 0; q < values_0.size(); ++q) + deallog << values_0[q][0] << " " << values_0[q][1] << " " + << values_0[q][2] << " " << std::endl; + + deallog << std::endl; + } + } +} + +int +main() +{ + initlog(); + + deallog.push("dim=2"); + test<2>(); + deallog.pop(); + + deallog.push("dim=3"); + test<3>(); + deallog.pop(); +} diff --git a/tests/simplex/variable_face_quadratures_01.output b/tests/simplex/variable_face_quadratures_01.output new file mode 100644 index 0000000000..fb86f84e7c --- /dev/null +++ b/tests/simplex/variable_face_quadratures_01.output @@ -0,0 +1,117 @@ + +DEAL:dim=2::face_no=0: +DEAL:dim=2::0.00000 0.500000 1.00000 +DEAL:dim=2:: +DEAL:dim=2::face_no=1: +DEAL:dim=2::1.00000 0.211325 0.500000 +DEAL:dim=2::1.00000 0.788675 0.500000 +DEAL:dim=2:: +DEAL:dim=2::face_no=2: +DEAL:dim=2::0.112702 0.00000 0.277778 +DEAL:dim=2::0.500000 0.00000 0.444444 +DEAL:dim=2::0.887298 0.00000 0.277778 +DEAL:dim=2:: +DEAL:dim=2::face_no=3: +DEAL:dim=2::0.0694318 1.00000 0.173927 +DEAL:dim=2::0.330009 1.00000 0.326073 +DEAL:dim=2::0.669991 1.00000 0.326073 +DEAL:dim=2::0.930568 1.00000 0.173927 +DEAL:dim=2:: +DEAL:dim=2:: +DEAL:dim=2::face_no=0: +DEAL:dim=2::0.00000 0.500000 1.00000 +DEAL:dim=2:: +DEAL:dim=2::face_no=1: +DEAL:dim=2::1.00000 0.211325 0.500000 +DEAL:dim=2::1.00000 0.788675 0.500000 +DEAL:dim=2:: +DEAL:dim=2::face_no=2: +DEAL:dim=2::0.112702 0.00000 0.277778 +DEAL:dim=2::0.500000 0.00000 0.444444 +DEAL:dim=2::0.887298 0.00000 0.277778 +DEAL:dim=2:: +DEAL:dim=2::face_no=3: +DEAL:dim=2::0.0694318 1.00000 0.173927 +DEAL:dim=2::0.330009 1.00000 0.326073 +DEAL:dim=2::0.669991 1.00000 0.326073 +DEAL:dim=2::0.930568 1.00000 0.173927 +DEAL:dim=2:: +DEAL:dim=2::0.00000 0.500000 +DEAL:dim=2:: +DEAL:dim=2::1.00000 0.211325 +DEAL:dim=2::1.00000 0.788675 +DEAL:dim=2:: +DEAL:dim=2::0.112702 0.00000 +DEAL:dim=2::0.500000 0.00000 +DEAL:dim=2::0.887298 0.00000 +DEAL:dim=2:: +DEAL:dim=2::0.0694318 1.00000 +DEAL:dim=2::0.330009 1.00000 +DEAL:dim=2::0.669991 1.00000 +DEAL:dim=2::0.930568 1.00000 +DEAL:dim=2:: +DEAL:dim=2::face_no=0: +DEAL:dim=2::0.00000 0.500000 +DEAL:dim=2:: +DEAL:dim=2::face_no=1: +DEAL:dim=2::1.00000 0.211325 +DEAL:dim=2::1.00000 0.788675 +DEAL:dim=2:: +DEAL:dim=2::face_no=2: +DEAL:dim=2::0.112702 0.00000 +DEAL:dim=2::0.500000 0.00000 +DEAL:dim=2::0.887298 0.00000 +DEAL:dim=2:: +DEAL:dim=2::face_no=3: +DEAL:dim=2::0.0694318 1.00000 +DEAL:dim=2::0.330009 1.00000 +DEAL:dim=2::0.669991 1.00000 +DEAL:dim=2::0.930568 1.00000 +DEAL:dim=2:: +DEAL:dim=3::face_no=0: +DEAL:dim=3::0.00000 0.500000 0.500000 +DEAL:dim=3:: +DEAL:dim=3::face_no=1: +DEAL:dim=3::1.00000 0.211325 0.211325 +DEAL:dim=3::1.00000 0.788675 0.211325 +DEAL:dim=3::1.00000 0.211325 0.788675 +DEAL:dim=3::1.00000 0.788675 0.788675 +DEAL:dim=3:: +DEAL:dim=3::face_no=2: +DEAL:dim=3::0.112702 0.00000 0.112702 +DEAL:dim=3::0.112702 0.00000 0.500000 +DEAL:dim=3::0.112702 0.00000 0.887298 +DEAL:dim=3::0.500000 0.00000 0.112702 +DEAL:dim=3::0.500000 0.00000 0.500000 +DEAL:dim=3::0.500000 0.00000 0.887298 +DEAL:dim=3::0.887298 0.00000 0.112702 +DEAL:dim=3::0.887298 0.00000 0.500000 +DEAL:dim=3::0.887298 0.00000 0.887298 +DEAL:dim=3:: +DEAL:dim=3::face_no=3: +DEAL:dim=3::0.0694318 1.00000 0.0694318 +DEAL:dim=3::0.0694318 1.00000 0.330009 +DEAL:dim=3::0.0694318 1.00000 0.669991 +DEAL:dim=3::0.0694318 1.00000 0.930568 +DEAL:dim=3::0.330009 1.00000 0.0694318 +DEAL:dim=3::0.330009 1.00000 0.330009 +DEAL:dim=3::0.330009 1.00000 0.669991 +DEAL:dim=3::0.330009 1.00000 0.930568 +DEAL:dim=3::0.669991 1.00000 0.0694318 +DEAL:dim=3::0.669991 1.00000 0.330009 +DEAL:dim=3::0.669991 1.00000 0.669991 +DEAL:dim=3::0.669991 1.00000 0.930568 +DEAL:dim=3::0.930568 1.00000 0.0694318 +DEAL:dim=3::0.930568 1.00000 0.330009 +DEAL:dim=3::0.930568 1.00000 0.669991 +DEAL:dim=3::0.930568 1.00000 0.930568 +DEAL:dim=3:: +DEAL:dim=3::face_no=4: +DEAL:dim=3::0.500000 0.500000 0.00000 +DEAL:dim=3:: +DEAL:dim=3::face_no=5: +DEAL:dim=3::0.211325 0.211325 1.00000 +DEAL:dim=3::0.788675 0.211325 1.00000 +DEAL:dim=3::0.211325 0.788675 1.00000 +DEAL:dim=3::0.788675 0.788675 1.00000 +DEAL:dim=3:: diff --git a/tests/simplex/variable_face_quadratures_02.cc b/tests/simplex/variable_face_quadratures_02.cc new file mode 100644 index 0000000000..609372a3f8 --- /dev/null +++ b/tests/simplex/variable_face_quadratures_02.cc @@ -0,0 +1,135 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test hp::FEFaceValues for different face quadrature rules. + + +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +class Fu : public Function +{ +public: + Fu(const unsigned int component = numbers::invalid_unsigned_int) + : Function(component == numbers::invalid_unsigned_int ? dim : 1) + , component(component) + {} + + virtual double + value(const Point &p, const unsigned int i = 0) const + { + if (component == numbers::invalid_unsigned_int) + return p[i]; + else + return p[component]; + } + + const unsigned int component; +}; + +template +void +test() +{} + +template <> +void +test<2>() +{ + const unsigned int dim = 2; + + // test FEFaceValues for FE_System(FE_Q) + { + const std::vector> quad_ref = { + hp::QCollection(QGauss(1), + QGauss(2), + QGauss(3), + QGauss(4)), + hp::QCollection(QGauss(5), + QGauss(6), + QGauss(7), + QGauss(8))}; + + hp::MappingCollection mapping(MappingFE(FE_Q(1))); + hp::FECollection fe(FESystem(FE_Q{3}, dim)); + + const UpdateFlags flags = + update_values | update_quadrature_points | update_JxW_values; + + + hp::FEFaceValues fe_face_values(mapping, fe, quad_ref, flags); + + Triangulation tria; + GridGenerator::hyper_cube(tria); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector vector_0(dof_handler.n_dofs()); + + VectorTools::interpolate(mapping[0 /*TODO*/], + dof_handler, + Fu(), + vector_0); + + for (const auto &cell : dof_handler.active_cell_iterators()) + for (const auto face_no : cell->face_indices()) + for (unsigned int q_index = 0; q_index < quad_ref.size(); ++q_index) + { + fe_face_values.reinit(cell, face_no, q_index); + + std::vector> values_0( + fe_face_values.get_present_fe_values().n_quadrature_points, + Vector(dim)); + + fe_face_values.get_present_fe_values().get_function_values( + vector_0, values_0); + + deallog << "q_index=" << q_index << " face_no=" << face_no << ":" + << std::endl; + + for (unsigned int q = 0; q < values_0.size(); ++q) + deallog << values_0[q][0] << " " << values_0[q][1] << " " + << std::endl; + + deallog << std::endl; + } + } +} + +int +main() +{ + initlog(); + + test<2>(); +} diff --git a/tests/simplex/variable_face_quadratures_02.output b/tests/simplex/variable_face_quadratures_02.output new file mode 100644 index 0000000000..243a59c5cd --- /dev/null +++ b/tests/simplex/variable_face_quadratures_02.output @@ -0,0 +1,53 @@ + +DEAL::q_index=0 face_no=0: +DEAL::0.00000 0.500000 +DEAL:: +DEAL::q_index=1 face_no=0: +DEAL::0.00000 0.0469101 +DEAL::0.00000 0.230765 +DEAL::0.00000 0.500000 +DEAL::0.00000 0.769235 +DEAL::0.00000 0.953090 +DEAL:: +DEAL::q_index=0 face_no=1: +DEAL::1.00000 0.211325 +DEAL::1.00000 0.788675 +DEAL:: +DEAL::q_index=1 face_no=1: +DEAL::1.00000 0.0337652 +DEAL::1.00000 0.169395 +DEAL::1.00000 0.380690 +DEAL::1.00000 0.619310 +DEAL::1.00000 0.830605 +DEAL::1.00000 0.966235 +DEAL:: +DEAL::q_index=0 face_no=2: +DEAL::0.112702 0.00000 +DEAL::0.500000 0.00000 +DEAL::0.887298 0.00000 +DEAL:: +DEAL::q_index=1 face_no=2: +DEAL::0.0254460 0.00000 +DEAL::0.129234 0.00000 +DEAL::0.297077 0.00000 +DEAL::0.500000 0.00000 +DEAL::0.702923 0.00000 +DEAL::0.870766 0.00000 +DEAL::0.974554 0.00000 +DEAL:: +DEAL::q_index=0 face_no=3: +DEAL::0.0694318 1.00000 +DEAL::0.330009 1.00000 +DEAL::0.669991 1.00000 +DEAL::0.930568 1.00000 +DEAL:: +DEAL::q_index=1 face_no=3: +DEAL::0.0198551 1.00000 +DEAL::0.101667 1.00000 +DEAL::0.237234 1.00000 +DEAL::0.408283 1.00000 +DEAL::0.591717 1.00000 +DEAL::0.762766 1.00000 +DEAL::0.898333 1.00000 +DEAL::0.980145 1.00000 +DEAL:: diff --git a/tests/simplex/variable_face_quadratures_03.cc b/tests/simplex/variable_face_quadratures_03.cc new file mode 100644 index 0000000000..95855659f9 --- /dev/null +++ b/tests/simplex/variable_face_quadratures_03.cc @@ -0,0 +1,183 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test FEFaceValues for different face quadrature rules. + + +#include + +#include +#include +#include + +#include + +#include + +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +template +class Fu : public Function +{ +public: + Fu(const unsigned int component = numbers::invalid_unsigned_int) + : Function(component == numbers::invalid_unsigned_int ? dim : 1) + , component(component) + {} + + virtual double + value(const Point &p, const unsigned int i = 0) const + { + if (component == numbers::invalid_unsigned_int) + return p[i]; + else + return p[component]; + } + + const unsigned int component; +}; + +template +void +test() +{} + +template <> +void +test<2>() +{ + const unsigned int dim = 2; + + // test FEFaceValues for FE_System(FE_P) + { + const hp::QCollection quad_ref(Simplex::QGauss(1), + Simplex::QGauss(2), + Simplex::QGauss(3)); + + MappingFE mapping(Simplex::FE_P(1)); + FESystem fe(Simplex::FE_P{2}, dim); + + const UpdateFlags flags = mapping.requires_update_flags( + update_values | update_quadrature_points | update_JxW_values); + + FEFaceValues fe_face_values(mapping, fe, quad_ref, flags); + + Triangulation tria; + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector vector_0(dof_handler.n_dofs()); + + VectorTools::interpolate(mapping, dof_handler, Fu(), vector_0); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + for (const auto face_no : cell->face_indices()) + { + fe_face_values.reinit(cell, face_no); + + std::vector> values_0(quad_ref[face_no].size(), + Vector(dim)); + fe_face_values.get_function_values(vector_0, values_0); + + deallog << "face_no=" << face_no << ":" << std::endl; + + for (unsigned int q = 0; q < values_0.size(); ++q) + deallog << values_0[q][0] << " " << values_0[q][1] << " " + << std::endl; + + deallog << std::endl; + } + break; + } + } +} + +template <> +void +test<3>() +{ + const unsigned int dim = 3; + + // test FEFaceValues for FE_System(FE_P) + { + const hp::QCollection quad_ref(Simplex::QGauss(1), + Simplex::QGauss(2), + Simplex::QGauss(3), + Simplex::QGauss(1)); + + MappingFE mapping(Simplex::FE_P(1)); + FESystem fe(Simplex::FE_P{2}, dim); + + const UpdateFlags flags = mapping.requires_update_flags( + update_values | update_quadrature_points | update_JxW_values); + + FEFaceValues fe_face_values(mapping, fe, quad_ref, flags); + + Triangulation tria; + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector vector_0(dof_handler.n_dofs()); + + VectorTools::interpolate(mapping, dof_handler, Fu(), vector_0); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + for (const auto face_no : cell->face_indices()) + { + fe_face_values.reinit(cell, face_no); + + std::vector> values_0(quad_ref[face_no].size(), + Vector(dim)); + fe_face_values.get_function_values(vector_0, values_0); + + deallog << "face_no=" << face_no << ":" << std::endl; + + for (unsigned int q = 0; q < values_0.size(); ++q) + deallog << values_0[q][0] << " " << values_0[q][1] << " " + << " " << values_0[q][2] << " " << std::endl; + + deallog << std::endl; + } + break; + } + } +} + +int +main() +{ + initlog(); + + deallog.push("dim=2"); + test<2>(); + deallog.pop(); + + deallog.push("dim=3"); + test<3>(); + deallog.pop(); +} diff --git a/tests/simplex/variable_face_quadratures_03.with_simplex_support=on.output b/tests/simplex/variable_face_quadratures_03.with_simplex_support=on.output new file mode 100644 index 0000000000..c8b88bfefd --- /dev/null +++ b/tests/simplex/variable_face_quadratures_03.with_simplex_support=on.output @@ -0,0 +1,33 @@ + +DEAL:dim=2::face_no=0: +DEAL:dim=2::0.500000 0.00000 +DEAL:dim=2:: +DEAL:dim=2::face_no=1: +DEAL:dim=2::0.788675 0.211325 +DEAL:dim=2::0.211325 0.788675 +DEAL:dim=2:: +DEAL:dim=2::face_no=2: +DEAL:dim=2::0.00000 0.887298 +DEAL:dim=2::0.00000 0.500000 +DEAL:dim=2::0.00000 0.112702 +DEAL:dim=2:: +DEAL:dim=3::face_no=0: +DEAL:dim=3::0.333333 0.333333 0.00000 +DEAL:dim=3:: +DEAL:dim=3::face_no=1: +DEAL:dim=3::0.166667 0.00000 0.166667 +DEAL:dim=3::0.166667 0.00000 0.666667 +DEAL:dim=3::0.666667 0.00000 0.166667 +DEAL:dim=3:: +DEAL:dim=3::face_no=2: +DEAL:dim=3::0.00000 0.333333 0.333333 +DEAL:dim=3::0.00000 0.797427 0.101287 +DEAL:dim=3::0.00000 0.101287 0.797427 +DEAL:dim=3::0.00000 0.101287 0.101287 +DEAL:dim=3::0.00000 0.0597159 0.470142 +DEAL:dim=3::0.00000 0.470142 0.0597159 +DEAL:dim=3::0.00000 0.470142 0.470142 +DEAL:dim=3:: +DEAL:dim=3::face_no=3: +DEAL:dim=3::0.333333 0.333333 0.333333 +DEAL:dim=3::