--- /dev/null
+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).
+<br>
+(Peter Munch, 2020/12/12)
#include <deal.II/grid/reference_cell.h>
+#include <deal.II/hp/q_collection.h>
+
DEAL_II_NAMESPACE_OPEN
* version of this function that takes the reference cell type instead.
*/
DEAL_II_DEPRECATED static Quadrature<dim>
- project_to_all_faces(const SubQuadrature &quadrature);
+ project_to_all_faces(const Quadrature<dim - 1> &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
* each face, in order to cope possibly different orientations of the mesh.
*/
static Quadrature<dim>
- 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<dim - 1> &quadrature);
+
+ /**
+ * Like the above function, applying the same face quadrature
+ * formula on all faces.
+ */
+ static Quadrature<dim>
+ project_to_all_faces(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<dim - 1> &quadrature);
/**
* Take a face quadrature formula and generate a cell quadrature formula
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<dim - 1> &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
}
+
+template <int dim>
+Quadrature<dim> inline QProjector<dim>::project_to_all_faces(
+ const Quadrature<dim - 1> &quadrature)
+{
+ return project_to_all_faces(ReferenceCell::get_hypercube(dim), quadrature);
+}
+
+
+template <int dim>
+Quadrature<dim> inline QProjector<dim>::project_to_all_faces(
+ const ReferenceCell::Type reference_cell_type,
+ const Quadrature<dim - 1> &quadrature)
+{
+ return project_to_all_faces(reference_cell_type,
+ hp::QCollection<dim - 1>(quadrature));
+}
+
+
/* -------------- declaration of explicit specializations ------------- */
#ifndef DOXYGEN
const unsigned int face_no,
std::vector<Point<3>> & 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 <>
* it is no longer necessary.
*/
virtual std::unique_ptr<InternalDataBase>
- get_face_data(const UpdateFlags update_flags,
- const Mapping<dim, spacedim> &mapping,
- const Quadrature<dim - 1> & quadrature,
+ get_face_data(const UpdateFlags update_flags,
+ const Mapping<dim, spacedim> & mapping,
+ const hp::QCollection<dim - 1> &quadrature,
dealii::internal::FEValuesImplementation::
FiniteElementRelatedData<dim, spacedim> &output_data) const;
+ /**
+ * @deprecated Use the version taking a hp::QCollection argument.
+ */
+ virtual std::unique_ptr<InternalDataBase>
+ get_face_data(
+ const UpdateFlags update_flags,
+ const Mapping<dim, spacedim> &mapping,
+ const Quadrature<dim - 1> & quadrature,
+ internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
+ &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
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
const InternalDataBase &fe_internal,
dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
spacedim>
- &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<dim, spacedim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature<dim - 1> & quadrature,
+ const Mapping<dim, spacedim> & mapping,
+ const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
+ const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ & mapping_data,
+ const InternalDataBase &fe_internal,
+ internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
+ &output_data) const;
/**
* This function is the equivalent to FiniteElement::fill_fe_values(), but
spacedim>
&output_data) const override;
+ using FiniteElement<dim, spacedim>::fill_fe_face_values;
+
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
spacedim>
&output_data) const override;
+ using FiniteElement<dim, spacedim>::get_face_data;
+
virtual std::unique_ptr<
typename FiniteElement<dim, spacedim>::InternalDataBase>
get_face_data(
- const UpdateFlags update_flags,
- const Mapping<dim, spacedim> &mapping,
- const Quadrature<dim - 1> & quadrature,
+ const UpdateFlags update_flags,
+ const Mapping<dim, spacedim> & mapping,
+ const hp::QCollection<dim - 1> &quadrature,
dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
spacedim>
&output_data) const override;
spacedim>
&output_data) const override;
+ using FiniteElement<dim, spacedim>::fill_fe_face_values;
+
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
typename FiniteElement<1, spacedim>::InternalDataBase>();
}
+ using FiniteElement<1, spacedim>::get_face_data;
+
std::unique_ptr<typename FiniteElement<1, spacedim>::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<typename FiniteElement<1, spacedim>::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;
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
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,
dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim>
&data) const override;
+ using FiniteElement<dim, spacedim>::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
fill_fe_face_values(
const typename Triangulation<dim, dim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, dim> & mapping,
const typename Mapping<dim, dim>::InternalDataBase & mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim, dim>
spacedim>
&output_data) const override;
+ using FiniteElement<dim, spacedim>::fill_fe_face_values;
+
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
dealii::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2>
&output_data) const override;
+ using FiniteElement<2, 2>::get_face_data;
+
virtual std::unique_ptr<FiniteElement<2, 2>::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;
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.
*/
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>
spacedim>
&output_data) const override;
+ using FiniteElement<dim, spacedim>::fill_fe_face_values;
+
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
FE_Poly<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
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
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
// 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));
{
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);
}
}
return std::make_unique<InternalData>();
}
+ using FiniteElement<dim, spacedim>::get_face_data;
+
std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
get_face_data(
- const UpdateFlags update_flags,
- const Mapping<dim, spacedim> &mapping,
- const Quadrature<dim - 1> & quadrature,
+ const UpdateFlags update_flags,
+ const Mapping<dim, spacedim> & mapping,
+ const hp::QCollection<dim - 1> &quadrature,
dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
spacedim>
&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<typename FiniteElement<dim, spacedim>::InternalDataBase>
auto &data = dynamic_cast<InternalData &>(*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<double> values(0);
std::vector<double>(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,
spacedim>
&output_data) const override
{
- return get_face_data(update_flags,
- mapping,
- QProjector<dim - 1>::project_to_all_children(
- ReferenceCell::get_hypercube(dim - 1), quadrature),
- output_data);
+ return get_face_data(
+ update_flags,
+ mapping,
+ hp::QCollection<dim - 1>(QProjector<dim - 1>::project_to_all_children(
+ ReferenceCell::get_hypercube(dim - 1), quadrature)),
+ output_data);
}
virtual void
spacedim>
&output_data) const override;
+ using FiniteElement<dim, spacedim>::fill_fe_face_values;
+
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
void
FE_PolyFace<PolynomialType, dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &,
- const unsigned int face_no,
- const Quadrature<dim - 1> &quadrature,
+ const unsigned int face_no,
+ const hp::QCollection<dim - 1> &quadrature,
const Mapping<dim, spacedim> &,
const typename Mapping<dim, spacedim>::InternalDataBase &,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
ExcInternalError());
const InternalData &fe_data = static_cast<const InternalData &>(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.;
spacedim>
&output_data) const override;
+ using FiniteElement<dim, spacedim>::fill_fe_face_values;
+
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
spacedim>
&output_data) const override;
+ using FiniteElement<dim, spacedim>::get_face_data;
+
virtual std::unique_ptr<
typename FiniteElement<dim, spacedim>::InternalDataBase>
get_face_data(
- const UpdateFlags update_flags,
- const Mapping<dim, spacedim> &mapping,
- const Quadrature<dim - 1> & quadrature,
+ const UpdateFlags update_flags,
+ const Mapping<dim, spacedim> & mapping,
+ const hp::QCollection<dim - 1> &quadrature,
dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
spacedim>
&output_data) const override;
spacedim>
&output_data) const override;
+ using FiniteElement<dim, spacedim>::fill_fe_face_values;
+
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int sub_no,
- const Quadrature<dim_1> & quadrature,
+ const hp::QCollection<dim_1> & quadrature,
const CellSimilarity::Similarity cell_similarity,
const typename Mapping<dim, spacedim>::InternalDataBase &mapping_internal,
const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_data,
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/hp/q_collection.h>
+
#include <algorithm>
#include <memory>
#include <type_traits>
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
const Quadrature<dim> & 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<dim, spacedim> & mapping,
+ const FiniteElement<dim, spacedim> &fe,
+ const hp::QCollection<dim> & 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
const Quadrature<dim> & 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<dim, spacedim> &fe,
+ const hp::QCollection<dim> & 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
* is <tt>2*dim*(1<<(dim-1))</tt>, 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<dim, spacedim> & mapping,
const FiniteElement<dim, spacedim> &fe,
const Quadrature<dim - 1> & 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<dim, spacedim> & mapping,
+ const FiniteElement<dim, spacedim> &fe,
+ const hp::QCollection<dim - 1> & quadrature);
+
/**
* Boundary form of the transformation of the cell at the <tt>i</tt>th
* quadrature point. See
/**
* Store a copy of the quadrature formula here.
*/
- const Quadrature<dim - 1> quadrature;
+ const hp::QCollection<dim - 1> quadrature;
};
const Quadrature<dim - 1> & 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<dim, spacedim> & mapping,
+ const FiniteElement<dim, spacedim> &fe,
+ const hp::QCollection<dim - 1> & 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
const Quadrature<dim - 1> & 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<dim, spacedim> &fe,
+ const hp::QCollection<dim - 1> & 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.
const Quadrature<dim - 1> & 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<dim, spacedim> & mapping,
+ const FiniteElement<dim, spacedim> &fe,
+ const hp::QCollection<dim - 1> & 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
const Quadrature<dim - 1> & 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<dim, spacedim> &fe,
+ const hp::QCollection<dim - 1> & 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
inline const Quadrature<dim - 1> &
FEFaceValuesBase<dim, spacedim>::get_quadrature() const
{
- return quadrature;
+ return quadrature[0];
}
#include <deal.II/grid/tria.h>
+#include <deal.II/hp/q_collection.h>
+
#include <array>
#include <cmath>
#include <memory>
* the returned object, knowing its real (derived) type.
*/
virtual std::unique_ptr<InternalDataBase>
+ get_face_data(const UpdateFlags update_flags,
+ const hp::QCollection<dim - 1> &quadrature) const;
+
+ /**
+ * @deprecated Use the version taking a hp::QCollection argument.
+ */
+ virtual std::unique_ptr<InternalDataBase>
get_face_data(const UpdateFlags update_flags,
- const Quadrature<dim - 1> &quadrature) const = 0;
+ const Quadrature<dim - 1> &quadrature) const;
/**
* Like get_data() and get_face_data(), but in preparation for later calls
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
- &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<dim, spacedim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature<dim - 1> & quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &output_data) const;
/**
* This function is the equivalent to Mapping::fill_fe_values(), but for
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
get_data(const UpdateFlags, const Quadrature<dim> &quadrature) const override;
+ using Mapping<dim, spacedim>::get_face_data;
+
// documentation can be found in Mapping::get_face_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
- get_face_data(const UpdateFlags flags,
- const Quadrature<dim - 1> &quadrature) const override;
+ get_face_data(const UpdateFlags flags,
+ const hp::QCollection<dim - 1> &quadrature) const override;
// documentation can be found in Mapping::get_subface_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
+ using Mapping<dim, spacedim>::fill_fe_face_values;
+
// documentation can be found in Mapping::fill_fe_face_values()
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
get_data(const UpdateFlags, const Quadrature<dim> &quadrature) const override;
+ using Mapping<dim, spacedim>::get_face_data;
+
// documentation can be found in Mapping::get_face_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
- get_face_data(const UpdateFlags flags,
- const Quadrature<dim - 1> &quadrature) const override;
+ get_face_data(const UpdateFlags flags,
+ const hp::QCollection<dim - 1> &quadrature) const override;
// documentation can be found in Mapping::get_subface_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
+ using Mapping<dim, spacedim>::fill_fe_face_values;
+
// documentation can be found in Mapping::fill_fe_face_values()
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
get_data(const UpdateFlags, const Quadrature<dim> &quadrature) const override;
+ using Mapping<dim, spacedim>::get_face_data;
+
// documentation can be found in Mapping::get_face_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
- get_face_data(const UpdateFlags flags,
- const Quadrature<dim - 1> &quadrature) const override;
+ get_face_data(const UpdateFlags flags,
+ const hp::QCollection<dim - 1> &quadrature) const override;
// documentation can be found in Mapping::get_subface_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
+ using Mapping<dim, spacedim>::fill_fe_face_values;
+
// documentation can be found in Mapping::fill_fe_face_values()
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
get_data(const UpdateFlags, const Quadrature<dim> &quadrature) const override;
+ using Mapping<dim, spacedim>::get_face_data;
+
// documentation can be found in Mapping::get_face_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
- get_face_data(const UpdateFlags flags,
- const Quadrature<dim - 1> &quadrature) const override;
+ get_face_data(const UpdateFlags flags,
+ const hp::QCollection<dim - 1> &quadrature) const override;
// documentation can be found in Mapping::get_subface_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
+ using Mapping<dim, spacedim>::fill_fe_face_values;
+
// documentation can be found in Mapping::fill_fe_face_values()
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
get_data(const UpdateFlags, const Quadrature<dim> &quadrature) const override;
+ using Mapping<dim, spacedim>::get_face_data;
+
// documentation can be found in Mapping::get_face_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
- get_face_data(const UpdateFlags flags,
- const Quadrature<dim - 1> &quadrature) const override;
+ get_face_data(const UpdateFlags flags,
+ const hp::QCollection<dim - 1> &quadrature) const override;
// documentation can be found in Mapping::get_subface_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
+ using Mapping<dim, spacedim>::fill_fe_face_values;
+
// documentation can be found in Mapping::fill_fe_face_values()
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
get_data(const UpdateFlags, const Quadrature<dim> &quadrature) const override;
+ using Mapping<dim, spacedim>::get_face_data;
+
// documentation can be found in Mapping::get_face_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
- get_face_data(const UpdateFlags flags,
- const Quadrature<dim - 1> &quadrature) const override;
+ get_face_data(const UpdateFlags flags,
+ const hp::QCollection<dim - 1> &quadrature) const override;
// documentation can be found in Mapping::get_subface_data()
virtual std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
+ using Mapping<dim, spacedim>::fill_fe_face_values;
+
// documentation can be found in Mapping::fill_fe_face_values()
virtual void
fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const override;
const QCollection<q_dim> & 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<dim, FEValuesType::space_dimension>
+ &mapping_collection,
+ const FECollection<dim, FEValuesType::space_dimension> &fe_collection,
+ const std::vector<QCollection<q_dim>> & 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
const QCollection<q_dim> & 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<dim, FEValuesType::space_dimension> &fe_collection,
+ const std::vector<QCollection<q_dim>> & q_collection,
+ const UpdateFlags update_flags);
+
/**
* Copy constructor.
*/
*/
const QCollection<q_dim> 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<QCollection<q_dim>> q_collections;
+
private:
/**
* A table in which we store pointers to fe_values objects for different
const hp::QCollection<dim - 1> & 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<dim, spacedim> &mapping_collection,
+ const hp::FECollection<dim, spacedim> & fe_collection,
+ const std::vector<hp::QCollection<dim - 1>> &q_collections,
+ const UpdateFlags update_flags);
+
/**
* Constructor. This constructor is equivalent to the other one except
const hp::QCollection<dim - 1> & 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<dim, spacedim> & fe_collection,
+ const std::vector<hp::QCollection<dim - 1>> &q_collections,
+ const UpdateFlags update_flags);
+
/**
* Reinitialize the object for the given cell and face.
*
#include <deal.II/base/quadrature.h>
#include <deal.II/base/subscriptor.h>
-#include <deal.II/fe/fe.h>
-
#include <memory>
#include <vector>
}
-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;
// 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));
}
std::vector<double> 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());
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<std::pair<std::array<Point<2>, 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)},
std::vector<double> 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<Point<2>, 2> support_points;
// determine support point of the current line with the correct
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)
{
const unsigned int dim = 2;
- const unsigned int n_points = quadrature.size(),
- n_faces = GeometryInfo<dim>::faces_per_cell;
+ const unsigned int n_faces = GeometryInfo<dim>::faces_per_cell;
+
+ unsigned int n_points_total = 0;
+
+ if (quadrature.size() == 1)
+ n_points_total = quadrature[0].size() * GeometryInfo<dim>::faces_per_cell;
+ else
+ {
+ AssertDimension(quadrature.size(), GeometryInfo<dim>::faces_per_cell);
+ for (unsigned int i = 0; i < quadrature.size(); ++i)
+ n_points_total += quadrature[i].size();
+ }
// first fix quadrature points
std::vector<Point<dim>> q_points;
- q_points.reserve(n_points * n_faces);
- std::vector<Point<dim>> help(n_points);
+ q_points.reserve(n_points_total);
+ std::vector<Point<dim>> 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<double> 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<dim>(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<std::pair<std::array<Point<3>, 3>, double>, 4> faces = {
{{{{Point<3>(0.0, 0.0, 0.0),
std::vector<double> 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<Point<3>, 3> support_points;
// determine support point of the current line with the correct
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)
{
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<dim>::faces_per_cell;
+ else
+ {
+ AssertDimension(quadrature.size(), GeometryInfo<dim>::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<dim>::faces_per_cell;
+ n_points_total *= 8;
// first fix quadrature points
std::vector<Point<dim>> q_points;
- q_points.reserve(n_points * n_faces * 8);
- std::vector<Point<dim>> help(n_points);
+ q_points.reserve(n_points_total);
+ std::vector<Point<dim>> help;
+ help.reserve(quadrature.max_n_quadrature_points());
std::vector<double> 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<dim>::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<dim>(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)
+template <int dim>
+typename QProjector<dim>::DataSetDescriptor
+QProjector<dim>::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<dim - 1> &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<dim>::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<dim>::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<dim>::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(
/*------------------------------- FiniteElement ----------------------*/
+#ifndef DOXYGEN
+template <int dim, int spacedim>
+std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
+FiniteElement<dim, spacedim>::get_face_data(
+ const UpdateFlags flags,
+ const Mapping<dim, spacedim> & mapping,
+ const hp::QCollection<dim - 1> &quadrature,
+ dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
+ spacedim>
+ &output_data) const
+{
+ return get_data(flags,
+ mapping,
+ QProjector<dim>::project_to_all_faces(
+ this->reference_cell_type(), quadrature),
+ output_data);
+}
+
+
+
template <int dim, int spacedim>
std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
FiniteElement<dim, spacedim>::get_face_data(
+template <int dim, int spacedim>
+inline void
+FiniteElement<dim, spacedim>::fill_fe_face_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const hp::QCollection<dim - 1> & quadrature,
+ const Mapping<dim, spacedim> & mapping,
+ const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
+ const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
+ spacedim>
+ & mapping_data,
+ const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
+ dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
+ spacedim>
+ &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 <int dim, int spacedim>
+inline void
+FiniteElement<dim, spacedim>::fill_fe_face_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature<dim - 1> & quadrature,
+ const Mapping<dim, spacedim> & mapping,
+ const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
+ const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
+ spacedim>
+ & mapping_data,
+ const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
+ dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
+ spacedim>
+ &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 <int dim, int spacedim>
std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
FiniteElement<dim, spacedim>::get_subface_data(
FE_DGPNonparametric<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &,
const unsigned int,
- const Quadrature<dim - 1> &,
+ const hp::QCollection<dim - 1> &,
const Mapping<dim, spacedim> &,
const typename Mapping<dim, spacedim>::InternalDataBase &,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
template <int dim, int spacedim>
std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
FE_Enriched<dim, spacedim>::get_face_data(
- const UpdateFlags update_flags,
- const Mapping<dim, spacedim> &mapping,
- const Quadrature<dim - 1> & quadrature,
+ const UpdateFlags update_flags,
+ const Mapping<dim, spacedim> & mapping,
+ const hp::QCollection<dim - 1> &quadrature,
internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
&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<dim, spacedim>::InternalData>(
std::move(data)),
update_flags,
- quadrature);
+ quadrature[0]);
}
FE_Enriched<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
ExcInternalError());
const InternalData &fe_data = static_cast<const InternalData &>(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,
if (is_enriched)
multiply_by_enrichment(
- quadrature, fe_data, mapping_data, cell, output_data);
+ quadrature[0], fe_data, mapping_data, cell, output_data);
}
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,
FE_NedelecSZ<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, dim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, dim> & mapping,
const typename Mapping<dim, dim>::InternalDataBase & mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim, dim>
dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim>
&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.
// (fe_internal/fe_data) which was not filled in by get_data.
fill_edge_values(cell,
QProjector<dim>::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<dim>::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<dim>::DataSetDescriptor::face(this->reference_cell_type(),
face_no,
FE_Nothing<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &,
const unsigned int,
- const Quadrature<dim - 1> &,
+ const hp::QCollection<dim - 1> &,
const Mapping<dim, spacedim> &,
const typename Mapping<dim, spacedim>::InternalDataBase &,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
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<FiniteElement<2, 2>::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
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> &,
dealii::internal::FEValuesImplementation::FiniteElementRelatedData<2, 2>
&output_data) const
{
+ AssertDimension(quadrature.size(), 1);
+
const UpdateFlags flags(fe_internal.update_each);
// linear shape functions
// 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)
FE_PolyTensor<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
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
ExcInternalError());
const InternalData &fe_data = static_cast<const InternalData &>(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)
template <int dim, int spacedim>
std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>
FESystem<dim, spacedim>::get_face_data(
- const UpdateFlags flags,
- const Mapping<dim, spacedim> &mapping,
- const Quadrature<dim - 1> & quadrature,
+ const UpdateFlags flags,
+ const Mapping<dim, spacedim> & mapping,
+ const hp::QCollection<dim - 1> &quadrature,
dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
spacedim>
& /*output_data*/) const
internal::FEValuesImplementation::FiniteElementRelatedData<dim, spacedim>
&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));
cell,
invalid_face_number,
invalid_face_number,
- quadrature,
+ hp::QCollection<dim>(quadrature),
cell_similarity,
mapping_internal,
fe_internal,
FESystem<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const Mapping<dim, spacedim> & mapping,
const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
const dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
cell,
face_no,
sub_no,
- quadrature,
+ hp::QCollection<dim - 1>(quadrature),
CellSimilarity::none,
mapping_internal,
fe_internal,
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int sub_no,
- const Quadrature<dim_1> & quadrature,
+ const hp::QCollection<dim_1> & quadrature,
const CellSimilarity::Similarity cell_similarity,
const typename Mapping<dim, spacedim>::InternalDataBase & mapping_internal,
const typename FiniteElement<dim, spacedim>::InternalDataBase &fe_internal,
// fill_fe_face_values needs argument Quadrature<dim-1> for both cases
// dim_1==dim-1 and dim_1=dim. Hence the following workaround
- const Quadrature<dim> * cell_quadrature = nullptr;
- const Quadrature<dim - 1> *face_quadrature = nullptr;
- const unsigned int n_q_points = quadrature.size();
+ const Quadrature<dim> * cell_quadrature = nullptr;
+ const hp::QCollection<dim - 1> *face_quadrature = nullptr;
+ const Quadrature<dim - 1> * 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<dim> or Quadrature<dim-1>:
- 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<const Quadrature<dim> *>(
quadrature_base_pointer) != nullptr,
cell_quadrature =
static_cast<const Quadrature<dim> *>(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<const hp::QCollection<dim - 1> *>(
+ quadrature_base_pointer) != nullptr,
+ ExcInternalError());
+
+ face_quadrature = static_cast<const hp::QCollection<dim - 1> *>(
+ quadrature_base_pointer);
+ }
else
{
+ const Subscriptor *quadrature_base_pointer = &quadrature[0];
Assert(dim_1 == dim - 1, ExcDimensionMismatch(dim_1, dim - 1));
Assert(dynamic_cast<const Quadrature<dim - 1> *>(
quadrature_base_pointer) != nullptr,
ExcInternalError());
- face_quadrature =
+ sub_face_quadrature =
static_cast<const Quadrature<dim - 1> *>(quadrature_base_pointer);
}
base_fe.fill_fe_subface_values(cell,
face_no,
sub_no,
- *face_quadrature,
+ *sub_face_quadrature,
mapping,
mapping_internal,
mapping_data,
&shape_function_data,
std::vector<typename ProductType<Number, double>::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(),
typename ProductType<Number, dealii::Tensor<order, spacedim>>::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(),
std::vector<typename Scalar<dim, spacedim>::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(),
typename ProductType<Number, dealii::Tensor<1, 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(),
typename ProductType<Number, dealii::Tensor<order + 1, spacedim>>::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(),
&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(),
std::vector<typename Vector<dim, spacedim>::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(),
Number,
typename dealii::internal::CurlType<spacedim>::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(),
std::vector<typename Vector<dim, spacedim>::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(),
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(),
std::vector<typename SymmetricTensor<2, dim, spacedim>::
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(),
typename ProductType<Number, dealii::Tensor<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(),
std::vector<typename Tensor<2, dim, spacedim>::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(),
std::vector<typename Tensor<2, dim, spacedim>::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(),
const Mapping<dim, spacedim> & mapping,
const FiniteElement<dim, spacedim> &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())
std::vector<Number> & 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(),
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;
const dealii::Table<2, Tensor<order, spacedim>> &shape_derivatives,
std::vector<Tensor<order, spacedim, Number>> & 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(),
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;
const dealii::Table<2, Tensor<2, spacedim>> &shape_hessians,
std::vector<Number> & 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(),
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
{
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) +
+template <int dim, int spacedim>
+FEValues<dim, spacedim>::FEValues(const Mapping<dim, spacedim> & mapping,
+ const FiniteElement<dim, spacedim> &fe,
+ const hp::QCollection<dim> & q,
+ const UpdateFlags update_flags)
+ : FEValues(mapping, fe, q[0], update_flags)
+{
+ AssertDimension(q.size(), 1);
+}
+
+
+
template <int dim, int spacedim>
FEValues<dim, spacedim>::FEValues(const FiniteElement<dim, spacedim> &fe,
const Quadrature<dim> & q,
+template <int dim, int spacedim>
+FEValues<dim, spacedim>::FEValues(const FiniteElement<dim, spacedim> &fe,
+ const hp::QCollection<dim> & q,
+ const UpdateFlags update_flags)
+ : FEValues(fe, q[0], update_flags)
+{
+ AssertDimension(q.size(), 1);
+}
+
+
+
template <int dim, int spacedim>
void
FEValues<dim, spacedim>::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);
template <int dim, int spacedim>
FEFaceValuesBase<dim, spacedim>::FEFaceValuesBase(
- const unsigned int n_q_points,
+ const unsigned int dofs_per_cell,
+ const UpdateFlags flags,
+ const Mapping<dim, spacedim> & mapping,
+ const FiniteElement<dim, spacedim> &fe,
+ const Quadrature<dim - 1> & quadrature)
+ : FEFaceValuesBase<dim, spacedim>(dofs_per_cell,
+ flags,
+ mapping,
+ fe,
+ hp::QCollection<dim - 1>(quadrature))
+{}
+
+
+
+template <int dim, int spacedim>
+FEFaceValuesBase<dim, spacedim>::FEFaceValuesBase(
const unsigned int dofs_per_cell,
const UpdateFlags,
const Mapping<dim, spacedim> & mapping,
const FiniteElement<dim, spacedim> &fe,
- const Quadrature<dim - 1> & quadrature)
- : FEValuesBase<dim, spacedim>(n_q_points,
+ const hp::QCollection<dim - 1> & quadrature)
+ : FEValuesBase<dim, spacedim>(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());
+}
const FiniteElement<dim, spacedim> &fe,
const Quadrature<dim - 1> & quadrature,
const UpdateFlags update_flags)
- : FEFaceValuesBase<dim, spacedim>(quadrature.size(),
- fe.n_dofs_per_cell(),
+ : FEFaceValues<dim, spacedim>(mapping,
+ fe,
+ hp::QCollection<dim - 1>(quadrature),
+ update_flags)
+{}
+
+
+
+template <int dim, int spacedim>
+FEFaceValues<dim, spacedim>::FEFaceValues(
+ const Mapping<dim, spacedim> & mapping,
+ const FiniteElement<dim, spacedim> &fe,
+ const hp::QCollection<dim - 1> & quadrature,
+ const UpdateFlags update_flags)
+ : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
update_flags,
mapping,
fe,
const FiniteElement<dim, spacedim> &fe,
const Quadrature<dim - 1> & quadrature,
const UpdateFlags update_flags)
- : FEFaceValuesBase<dim, spacedim>(quadrature.size(),
- fe.n_dofs_per_cell(),
+ : FEFaceValues<dim, spacedim>(fe,
+ hp::QCollection<dim - 1>(quadrature),
+ update_flags)
+{}
+
+
+
+template <int dim, int spacedim>
+FEFaceValues<dim, spacedim>::FEFaceValues(
+ const FiniteElement<dim, spacedim> &fe,
+ const hp::QCollection<dim - 1> & quadrature,
+ const UpdateFlags update_flags)
+ : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
update_flags,
StaticMappingQ1<dim, spacedim>::mapping,
fe,
// 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<typename FiniteElement<dim, spacedim>::InternalDataBase> (
+ FiniteElement<dim, spacedim>::*finite_element_get_face_data)(
+ const UpdateFlags,
+ const Mapping<dim, spacedim> &,
+ const hp::QCollection<dim - 1> &,
+ dealii::internal::FEValuesImplementation::FiniteElementRelatedData<dim,
+ spacedim>
+ &) const = &FiniteElement<dim, spacedim>::get_face_data;
+
+ std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> (
+ Mapping<dim, spacedim>::*mapping_get_face_data)(
+ const UpdateFlags, const hp::QCollection<dim - 1> &) const =
+ &Mapping<dim, spacedim>::get_face_data;
+
+
Threads::Task<
std::unique_ptr<typename FiniteElement<dim, spacedim>::InternalDataBase>>
- fe_get_data =
- Threads::new_task(&FiniteElement<dim, spacedim>::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<typename Mapping<dim, spacedim>::InternalDataBase>>
mapping_get_data;
if (flags & update_mapping)
- mapping_get_data = Threads::new_task(&Mapping<dim, spacedim>::get_face_data,
+ mapping_get_data = Threads::new_task(mapping_get_face_data,
*this->mapping,
flags,
this->quadrature);
this->mapping_output,
*this->fe_data,
this->finite_element_output);
+
+ const_cast<unsigned int &>(this->n_quadrature_points) =
+ this->quadrature[this->quadrature.size() == 1 ? 0 : face_no].size();
}
const FiniteElement<dim, spacedim> &fe,
const Quadrature<dim - 1> & quadrature,
const UpdateFlags update_flags)
- : FEFaceValuesBase<dim, spacedim>(quadrature.size(),
- fe.n_dofs_per_cell(),
+ : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
update_flags,
mapping,
fe,
+template <int dim, int spacedim>
+FESubfaceValues<dim, spacedim>::FESubfaceValues(
+ const Mapping<dim, spacedim> & mapping,
+ const FiniteElement<dim, spacedim> &fe,
+ const hp::QCollection<dim - 1> & quadrature,
+ const UpdateFlags update_flags)
+ : FESubfaceValues(mapping, fe, quadrature[0], update_flags)
+{
+ AssertDimension(quadrature.size(), 1);
+}
+
+
+
template <int dim, int spacedim>
FESubfaceValues<dim, spacedim>::FESubfaceValues(
const FiniteElement<dim, spacedim> &fe,
const Quadrature<dim - 1> & quadrature,
const UpdateFlags update_flags)
- : FEFaceValuesBase<dim, spacedim>(quadrature.size(),
- fe.n_dofs_per_cell(),
+ : FEFaceValuesBase<dim, spacedim>(fe.n_dofs_per_cell(),
update_flags,
StaticMappingQ1<dim, spacedim>::mapping,
fe,
+template <int dim, int spacedim>
+FESubfaceValues<dim, spacedim>::FESubfaceValues(
+ const FiniteElement<dim, spacedim> &fe,
+ const hp::QCollection<dim - 1> & quadrature,
+ const UpdateFlags update_flags)
+ : FESubfaceValues(fe, quadrature[0], update_flags)
+{
+ AssertDimension(quadrature.size(), 1);
+}
+
+
+
template <int dim, int spacedim>
void
FESubfaceValues<dim, spacedim>::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);
*this->fe,
flags,
*this->mapping,
- this->quadrature,
+ this->quadrature[0],
this->finite_element_output);
Threads::Task<
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
Threads::new_task(&Mapping<dim, spacedim>::get_subface_data,
*this->mapping,
flags,
- this->quadrature);
+ this->quadrature[0]);
this->update_flags = flags;
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);
}
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,
return {};
}
+
+
+#ifndef DOXYGEN
+template <int dim, int spacedim>
+void
+Mapping<dim, spacedim>::fill_fe_face_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const hp::QCollection<dim - 1> & quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &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 <int dim, int spacedim>
+void
+Mapping<dim, spacedim>::fill_fe_face_values(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature<dim - 1> & quadrature,
+ const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
+ dealii::internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+ &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 <int dim, int spacedim>
+std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+Mapping<dim, spacedim>::get_face_data(
+ const UpdateFlags update_flags,
+ const hp::QCollection<dim - 1> &quadrature) const
+{
+ // base class version, implement overriden function in derived classes
+ return get_face_data(update_flags, quadrature[0]);
+}
+
+
+
+template <int dim, int spacedim>
+std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+Mapping<dim, spacedim>::get_face_data(
+ const UpdateFlags update_flags,
+ const Quadrature<dim - 1> &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<typename Mapping<dim, spacedim>::InternalDataBase>();
+}
+#endif
+
/* ---------------------------- InternalDataBase --------------------------- */
template <int dim, int spacedim>
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
MappingCartesian<dim, spacedim>::get_face_data(
- const UpdateFlags update_flags,
- const Quadrature<dim - 1> &quadrature) const
+ const UpdateFlags update_flags,
+ const hp::QCollection<dim - 1> &quadrature) const
{
+ AssertDimension(quadrature.size(), 1);
+
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(
QProjector<dim>::project_to_all_faces(ReferenceCell::get_hypercube(dim),
- quadrature));
+ quadrature[0]));
auto &data = dynamic_cast<InternalData &>(*data_ptr);
// verify that we have computed the transitive hull of the required
MappingCartesian<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
internal::FEValuesImplementation::MappingRelatedData<dim, 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
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)
maybe_compute_q_points(
const typename QProjector<dim>::DataSetDescriptor data_set,
const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
- std::vector<Point<spacedim>> &quadrature_points)
+ std::vector<Point<spacedim>> &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<spacedim> result =
maybe_update_Jacobians(
const CellSimilarity::Similarity cell_similarity,
const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
- const typename dealii::MappingFE<dim, spacedim>::InternalData &data)
+ const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
+ const unsigned int n_q_points)
{
const UpdateFlags update_flags = data.update_each;
// 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>());
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] =
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();
const CellSimilarity::Similarity cell_similarity,
const typename QProjector<dim>::DataSetDescriptor data_set,
const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
- std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
+ std::vector<DerivativeForm<2, dim, spacedim>> &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)
const CellSimilarity::Similarity cell_similarity,
const typename QProjector<dim>::DataSetDescriptor data_set,
const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
- std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
+ std::vector<Tensor<3, spacedim>> &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)
{
const CellSimilarity::Similarity cell_similarity,
const typename QProjector<dim>::DataSetDescriptor data_set,
const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
- std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
+ std::vector<DerivativeForm<3, dim, spacedim>> &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)
{
const typename QProjector<dim>::DataSetDescriptor data_set,
const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
std::vector<Tensor<4, spacedim>>
- &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)
{
const CellSimilarity::Similarity cell_similarity,
const typename QProjector<dim>::DataSetDescriptor data_set,
const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
- std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
+ std::vector<DerivativeForm<4, dim, spacedim>> &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)
{
const typename QProjector<dim>::DataSetDescriptor data_set,
const typename dealii::MappingFE<dim, spacedim>::InternalData &data,
std::vector<Tensor<5, spacedim>>
- &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)
{
template <int dim, int spacedim>
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
MappingFE<dim, spacedim>::get_face_data(
- const UpdateFlags update_flags,
- const Quadrature<dim - 1> &quadrature) const
+ const UpdateFlags update_flags,
+ const hp::QCollection<dim - 1> &quadrature) const
{
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(*this->fe);
data.initialize_face(this->requires_update_flags(update_flags),
QProjector<dim>::project_to_all_faces(
this->fe->reference_cell_type(), quadrature),
- quadrature.size());
+ quadrature.max_n_quadrature_points());
return data_ptr;
}
internal::MappingFEImplementation::maybe_compute_q_points<dim, spacedim>(
QProjector<dim>::DataSetDescriptor::cell(),
data,
- output_data.quadrature_points);
+ output_data.quadrature_points,
+ n_q_points);
internal::MappingFEImplementation::maybe_update_Jacobians<dim, spacedim>(
- computed_cell_similarity, QProjector<dim>::DataSetDescriptor::cell(), data);
+ computed_cell_similarity,
+ QProjector<dim>::DataSetDescriptor::cell(),
+ data,
+ n_q_points);
internal::MappingFEImplementation::maybe_update_jacobian_grads<dim, spacedim>(
computed_cell_similarity,
QProjector<dim>::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<dim>::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<dim>::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<dim, spacedim>(
computed_cell_similarity,
QProjector<dim>::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<dim>::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<dim, spacedim>(
computed_cell_similarity,
QProjector<dim>::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<double> &weights = quadrature.get_weights();
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());
//
// 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)
{
}
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() *
}
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<spacedim>(output_data.boundary_forms[i] /
output_data.boundary_forms[i].norm());
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data)
{
+ const unsigned int n_q_points = quadrature.size();
+
maybe_compute_q_points<dim, spacedim>(data_set,
data,
- output_data.quadrature_points);
+ output_data.quadrature_points,
+ n_q_points);
maybe_update_Jacobians<dim, spacedim>(CellSimilarity::none,
data_set,
- data);
+ data,
+ n_q_points);
maybe_update_jacobian_grads<dim, spacedim>(CellSimilarity::none,
data_set,
data,
- output_data.jacobian_grads);
+ output_data.jacobian_grads,
+ n_q_points);
maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
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<dim, spacedim>(
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<dim, spacedim>(
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<dim, spacedim>(
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<dim, spacedim>(
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);
MappingFE<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const
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);
}
template <int dim, int spacedim, typename VectorType>
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
MappingFEField<dim, spacedim, VectorType, void>::get_face_data(
- const UpdateFlags update_flags,
- const Quadrature<dim - 1> &quadrature) const
+ const UpdateFlags update_flags,
+ const hp::QCollection<dim - 1> &quadrature) const
{
+ AssertDimension(quadrature.size(), 1);
+
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
auto & data = dynamic_cast<InternalData &>(*data_ptr);
const Quadrature<dim> q(
QProjector<dim>::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;
}
MappingFEField<dim, spacedim, VectorType, void>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
internal::FEValuesImplementation::MappingRelatedData<dim, 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 possible
Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
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,
template <int dim, int spacedim>
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
MappingManifold<dim, spacedim>::get_face_data(
- const UpdateFlags update_flags,
- const Quadrature<dim - 1> &quadrature) const
+ const UpdateFlags update_flags,
+ const hp::QCollection<dim - 1> &quadrature) const
{
+ AssertDimension(quadrature.size(), 1);
+
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>();
auto &data = dynamic_cast<InternalData &>(*data_ptr);
data.initialize_face(this->requires_update_flags(update_flags),
QProjector<dim>::project_to_all_faces(
- ReferenceCell::get_hypercube(dim), quadrature),
- quadrature.size());
+ ReferenceCell::get_hypercube(dim), quadrature[0]),
+ quadrature[0].size());
return data_ptr;
}
MappingManifold<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const
{
+ AssertDimension(quadrature.size(), 1);
+
// ensure that the following cast is really correct:
Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
ExcInternalError());
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);
}
template <int dim, int spacedim>
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
MappingQ<dim, spacedim>::get_face_data(
- const UpdateFlags update_flags,
- const Quadrature<dim - 1> &quadrature) const
+ const UpdateFlags update_flags,
+ const hp::QCollection<dim - 1> &quadrature) const
{
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>();
auto &data = dynamic_cast<InternalData &>(*data_ptr);
+ std::unique_ptr<typename MappingQGeneric<dim, spacedim>::InternalDataBase> (
+ MappingQGeneric<dim, spacedim>::*mapping_get_face_data)(
+ const UpdateFlags, const hp::QCollection<dim - 1> &) const =
+ &MappingQGeneric<dim, spacedim>::get_face_data;
+
// build the Q1 and Qp internal data objects in parallel
Threads::Task<
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>>
- do_get_data =
- Threads::new_task(&MappingQGeneric<dim, spacedim>::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<
MappingQ<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const
template <int dim, int spacedim>
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
MappingQGeneric<dim, spacedim>::get_face_data(
- const UpdateFlags update_flags,
- const Quadrature<dim - 1> &quadrature) const
+ const UpdateFlags update_flags,
+ const hp::QCollection<dim - 1> &quadrature) const
{
+ AssertDimension(quadrature.size(), 1);
+
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(polynomial_degree);
auto &data = dynamic_cast<InternalData &>(*data_ptr);
data.initialize_face(this->requires_update_flags(update_flags),
QProjector<dim>::project_to_all_faces(
- ReferenceCell::get_hypercube(dim), quadrature),
- quadrature.size());
+ ReferenceCell::get_hypercube(dim), quadrature[0]),
+ quadrature[0].size());
return data_ptr;
}
MappingQGeneric<dim, spacedim>::fill_fe_face_values(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const unsigned int face_no,
- const Quadrature<dim - 1> & quadrature,
+ const hp::QCollection<dim - 1> & quadrature,
const typename Mapping<dim, spacedim>::InternalDataBase & internal_data,
internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
&output_data) const
{
+ AssertDimension(quadrature.size(), 1);
+
// ensure that the following cast is really correct:
Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
ExcInternalError());
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);
}
namespace hp
{
+ namespace
+ {
+ /**
+ * Given a quadrature collection, create a vector of quadrature collections
+ * with each containing one quadrature rule from the input collection.
+ */
+ template <int q_dim>
+ std::vector<QCollection<q_dim>>
+ translate(const QCollection<q_dim> &q_collection)
+ {
+ std::vector<QCollection<q_dim>> 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 <int q_dim>
+ QCollection<q_dim>
+ translate(const std::vector<QCollection<q_dim>> &q_collections)
+ {
+ QCollection<q_dim> result;
+
+ for (unsigned int q = 0; q < q_collections.size(); ++q)
+ result.push_back(q_collections[q][0]);
+
+ return result;
+ }
+ } // namespace
+
// -------------------------- FEValuesBase -------------------------
template <int dim, int q_dim, class FEValuesType>
: 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())
, update_flags(update_flags)
{}
-
template <int dim, int q_dim, class FEValuesType>
FEValuesBase<dim, q_dim, FEValuesType>::FEValuesBase(
+ const MappingCollection<dim, FEValuesType::space_dimension>
+ & mapping_collection,
const FECollection<dim, FEValuesType::space_dimension> &fe_collection,
- const QCollection<q_dim> & q_collection,
+ const std::vector<QCollection<q_dim>> & q_collections,
const UpdateFlags update_flags)
: fe_collection(&fe_collection)
- , mapping_collection(
- &dealii::hp::StaticMappingQ1<dim, FEValuesType::space_dimension>::
- 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)
{}
+ template <int dim, int q_dim, class FEValuesType>
+ FEValuesBase<dim, q_dim, FEValuesType>::FEValuesBase(
+ const FECollection<dim, FEValuesType::space_dimension> &fe_collection,
+ const QCollection<q_dim> & q_collection,
+ const UpdateFlags update_flags)
+ : FEValuesBase(
+ dealii::hp::StaticMappingQ1<dim, FEValuesType::space_dimension>::
+ mapping_collection,
+ fe_collection,
+ q_collection,
+ update_flags)
+ {}
+
+
+ template <int dim, int q_dim, class FEValuesType>
+ FEValuesBase<dim, q_dim, FEValuesType>::FEValuesBase(
+ const FECollection<dim, FEValuesType::space_dimension> &fe_collection,
+ const std::vector<QCollection<q_dim>> & q_collection,
+ const UpdateFlags update_flags)
+ : FEValuesBase(
+ dealii::hp::StaticMappingQ1<dim, FEValuesType::space_dimension>::
+ mapping_collection,
+ fe_collection,
+ q_collection,
+ update_flags)
+ {}
+
+
template <int dim, int q_dim, class FEValuesType>
FEValuesBase<dim, q_dim, FEValuesType>::FEValuesBase(
: 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)
{
// 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<FEValuesType>((*mapping_collection)[m_index],
(*fe_collection)[fe_index],
- q_collection[q_index],
+ q_collections[q_index],
update_flags);
});
{
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
fe_values_table(present_fe_values_index) =
std::make_unique<FEValuesType>((*mapping_collection)[mapping_index],
(*fe_collection)[fe_index],
- q_collection[q_index],
+ q_collections[q_index],
update_flags);
// now there definitely is one!
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]() {
std::make_unique<FEValuesType>(
(*mapping_collection)[mapping_index],
(*fe_collection)[fe_index],
- q_collection[q_index],
+ q_collections[q_index],
update_flags);
});
}
indices :
std::vector<unsigned int>(size, 0),
/*q_indices=*/
- (q_collection.size() > 1) ?
+ (q_collections.size() > 1) ?
indices :
std::vector<unsigned int>(size, 0));
}
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;
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());
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());
update_flags)
{}
+ template <int dim, int spacedim>
+ FEFaceValues<dim, spacedim>::FEFaceValues(
+ const hp::MappingCollection<dim, spacedim> & mapping,
+ const hp::FECollection<dim, spacedim> & fe_collection,
+ const std::vector<hp::QCollection<dim - 1>> &q_collection,
+ const UpdateFlags update_flags)
+ : hp::FEValuesBase<dim, dim - 1, dealii::FEFaceValues<dim, spacedim>>(
+ mapping,
+ fe_collection,
+ q_collection,
+ update_flags)
+ {}
+
template <int dim, int spacedim>
FEFaceValues<dim, spacedim>::FEFaceValues(
{}
+ template <int dim, int spacedim>
+ FEFaceValues<dim, spacedim>::FEFaceValues(
+ const hp::FECollection<dim, spacedim> & fe_collection,
+ const std::vector<hp::QCollection<dim - 1>> &q_collection,
+ const UpdateFlags update_flags)
+ : hp::FEValuesBase<dim, dim - 1, dealii::FEFaceValues<dim, spacedim>>(
+ fe_collection,
+ q_collection,
+ update_flags)
+ {}
+
+
template <int dim, int spacedim>
template <bool lda>
void
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;
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());
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());
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;
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());
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());
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/qprojector.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+class Fu : public Function<dim>
+{
+public:
+ Fu(const unsigned int component = numbers::invalid_unsigned_int)
+ : Function<dim>(component == numbers::invalid_unsigned_int ? dim : 1)
+ , component(component)
+ {}
+
+ virtual double
+ value(const Point<dim> &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 <int dim>
+void
+test()
+{}
+
+template <>
+void
+test<2>()
+{
+ const unsigned int dim = 2;
+
+ // test: QProjector::project_to_all_faces
+ {
+ const hp::QCollection<dim - 1> quad_ref(QGauss<dim - 1>(1),
+ QGauss<dim - 1>(2),
+ QGauss<dim - 1>(3),
+ QGauss<dim - 1>(4));
+
+ const auto quad =
+ QProjector<dim>::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<dim>::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<dim - 1> quad_ref(QGauss<dim - 1>(1),
+ QGauss<dim - 1>(2),
+ QGauss<dim - 1>(3),
+ QGauss<dim - 1>(4));
+
+ MappingFE<dim> mapping(FE_Q<dim>(1));
+ FE_Q<dim> 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<dim> data;
+ data.initialize(quad_ref.max_n_quadrature_points(), flags);
+
+ Triangulation<dim> 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<dim> 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<dim - 1> quad_ref(QGauss<dim - 1>(1),
+ QGauss<dim - 1>(2),
+ QGauss<dim - 1>(3),
+ QGauss<dim - 1>(4));
+
+ MappingFE<dim> mapping(FE_Q<dim>(1));
+ FE_Q<dim> fe(3);
+
+ const UpdateFlags flags = mapping.requires_update_flags(
+ update_values | update_quadrature_points | update_JxW_values);
+
+
+ FEFaceValues<dim> fe_face_values(mapping, fe, quad_ref, flags);
+
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ Vector<double> vector_0(dof_handler.n_dofs());
+ Vector<double> vector_1(dof_handler.n_dofs());
+
+ VectorTools::interpolate(mapping, dof_handler, Fu<dim>(0), vector_0);
+ VectorTools::interpolate(mapping, dof_handler, Fu<dim>(1), vector_1);
+
+ std::vector<double> 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<dim - 1> quad_ref(QGauss<dim - 1>(1),
+ QGauss<dim - 1>(2),
+ QGauss<dim - 1>(3),
+ QGauss<dim - 1>(4));
+
+ MappingFE<dim> mapping(FE_Q<dim>(1));
+ FESystem<dim> fe(FE_Q<dim>{3}, dim);
+
+ const UpdateFlags flags = mapping.requires_update_flags(
+ update_values | update_quadrature_points | update_JxW_values);
+
+
+ FEFaceValues<dim> fe_face_values(mapping, fe, quad_ref, flags);
+
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ Vector<double> vector_0(dof_handler.n_dofs());
+
+ VectorTools::interpolate(mapping, dof_handler, Fu<dim>(), 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<Vector<double>> values_0(quad_ref[face_no].size(),
+ Vector<double>(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<dim - 1> quad_ref(QGauss<dim - 1>(1),
+ QGauss<dim - 1>(2),
+ QGauss<dim - 1>(3),
+ QGauss<dim - 1>(4),
+ QGauss<dim - 1>(1),
+ QGauss<dim - 1>(2));
+
+ MappingFE<dim> mapping(FE_Q<dim>(1));
+ FESystem<dim> fe(FE_Q<dim>{3}, dim);
+
+ const UpdateFlags flags =
+ update_values | update_quadrature_points | update_JxW_values;
+
+
+ FEFaceValues<dim> fe_face_values(mapping, fe, quad_ref, flags);
+
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ Vector<double> vector_0(dof_handler.n_dofs());
+
+ VectorTools::interpolate(mapping, dof_handler, Fu<dim>(), 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<Vector<double>> values_0(quad_ref[face_no].size(),
+ Vector<double>(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();
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/qprojector.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/hp/fe_values.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+class Fu : public Function<dim>
+{
+public:
+ Fu(const unsigned int component = numbers::invalid_unsigned_int)
+ : Function<dim>(component == numbers::invalid_unsigned_int ? dim : 1)
+ , component(component)
+ {}
+
+ virtual double
+ value(const Point<dim> &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 <int dim>
+void
+test()
+{}
+
+template <>
+void
+test<2>()
+{
+ const unsigned int dim = 2;
+
+ // test FEFaceValues for FE_System(FE_Q)
+ {
+ const std::vector<hp::QCollection<dim - 1>> quad_ref = {
+ hp::QCollection<dim - 1>(QGauss<dim - 1>(1),
+ QGauss<dim - 1>(2),
+ QGauss<dim - 1>(3),
+ QGauss<dim - 1>(4)),
+ hp::QCollection<dim - 1>(QGauss<dim - 1>(5),
+ QGauss<dim - 1>(6),
+ QGauss<dim - 1>(7),
+ QGauss<dim - 1>(8))};
+
+ hp::MappingCollection<dim> mapping(MappingFE<dim>(FE_Q<dim>(1)));
+ hp::FECollection<dim> fe(FESystem<dim>(FE_Q<dim>{3}, dim));
+
+ const UpdateFlags flags =
+ update_values | update_quadrature_points | update_JxW_values;
+
+
+ hp::FEFaceValues<dim> fe_face_values(mapping, fe, quad_ref, flags);
+
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ Vector<double> vector_0(dof_handler.n_dofs());
+
+ VectorTools::interpolate(mapping[0 /*TODO*/],
+ dof_handler,
+ Fu<dim>(),
+ 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<Vector<double>> values_0(
+ fe_face_values.get_present_fe_values().n_quadrature_points,
+ Vector<double>(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>();
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/qprojector.h>
+
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/hp/q_collection.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+class Fu : public Function<dim>
+{
+public:
+ Fu(const unsigned int component = numbers::invalid_unsigned_int)
+ : Function<dim>(component == numbers::invalid_unsigned_int ? dim : 1)
+ , component(component)
+ {}
+
+ virtual double
+ value(const Point<dim> &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 <int dim>
+void
+test()
+{}
+
+template <>
+void
+test<2>()
+{
+ const unsigned int dim = 2;
+
+ // test FEFaceValues for FE_System(FE_P)
+ {
+ const hp::QCollection<dim - 1> quad_ref(Simplex::QGauss<dim - 1>(1),
+ Simplex::QGauss<dim - 1>(2),
+ Simplex::QGauss<dim - 1>(3));
+
+ MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+ FESystem<dim> fe(Simplex::FE_P<dim>{2}, dim);
+
+ const UpdateFlags flags = mapping.requires_update_flags(
+ update_values | update_quadrature_points | update_JxW_values);
+
+ FEFaceValues<dim> fe_face_values(mapping, fe, quad_ref, flags);
+
+ Triangulation<dim> tria;
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ Vector<double> vector_0(dof_handler.n_dofs());
+
+ VectorTools::interpolate(mapping, dof_handler, Fu<dim>(), 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<Vector<double>> values_0(quad_ref[face_no].size(),
+ Vector<double>(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<dim - 1> quad_ref(Simplex::QGauss<dim - 1>(1),
+ Simplex::QGauss<dim - 1>(2),
+ Simplex::QGauss<dim - 1>(3),
+ Simplex::QGauss<dim - 1>(1));
+
+ MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+ FESystem<dim> fe(Simplex::FE_P<dim>{2}, dim);
+
+ const UpdateFlags flags = mapping.requires_update_flags(
+ update_values | update_quadrature_points | update_JxW_values);
+
+ FEFaceValues<dim> fe_face_values(mapping, fe, quad_ref, flags);
+
+ Triangulation<dim> tria;
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1);
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ Vector<double> vector_0(dof_handler.n_dofs());
+
+ VectorTools::interpolate(mapping, dof_handler, Fu<dim>(), 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<Vector<double>> values_0(quad_ref[face_no].size(),
+ Vector<double>(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();
+}
--- /dev/null
+
+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::
projects / dealii.git / commitdiff