From f791b78827fa0444093b31dcb2f48ba6e1a32d82 Mon Sep 17 00:00:00 2001 From: David Wells <drwells@email.unc.edu> Date: Sat, 31 Aug 2019 17:21:23 -0400 Subject: [PATCH] Permit reinitialization with iterators instead of (sub)face numbers. This is convenient when paired with the new cell->face_iterators() function. --- doc/news/changes/minor/20191027DavidWells | 4 ++ include/deal.II/fe/fe_values.h | 66 ++++++++++++++++++ source/fe/fe_values.cc | 54 +++++++++++++++ source/fe/fe_values.inst.in | 26 ++++++- tests/fe/shapes.h | 82 ++++++++++++++++++----- 5 files changed, 212 insertions(+), 20 deletions(-) create mode 100644 doc/news/changes/minor/20191027DavidWells diff --git a/doc/news/changes/minor/20191027DavidWells b/doc/news/changes/minor/20191027DavidWells new file mode 100644 index 0000000000..fb70fa6b38 --- /dev/null +++ b/doc/news/changes/minor/20191027DavidWells @@ -0,0 +1,4 @@ +New: Added two new overloads of FEValues::reinit() that take face and subface +iterators instead of face and subface numbers. +<br> +(David Wells, 2019/10/27) diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h index 8229d7b760..da9b132a57 100644 --- a/include/deal.II/fe/fe_values.h +++ b/include/deal.II/fe/fe_values.h @@ -3757,6 +3757,18 @@ public: level_dof_access>> &cell, const unsigned int face_no); + /** + * Reinitialize the gradients, Jacobi determinants, etc for face @p face + * and cell @p cell. + * + * @note @p face must be one of @p cell's face iterators. + */ + template <template <int, int> class DoFHandlerType, bool level_dof_access> + void + reinit(const TriaIterator<DoFCellAccessor<DoFHandlerType<dim, spacedim>, + level_dof_access>> & cell, + const typename Triangulation<dim, spacedim>::face_iterator &face); + /** * Reinitialize the gradients, Jacobi determinants, etc for the given face * on a given cell of type "iterator into a Triangulation object", and the @@ -3774,6 +3786,25 @@ public: reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell, const unsigned int face_no); + /* + * Reinitialize the gradients, Jacobi determinants, etc for the given face + * on a given cell of type "iterator into a Triangulation object", and the + * given finite element. Since iterators into a triangulation alone only + * convey information about the geometry of a cell, but not about degrees of + * freedom possibly associated with this cell, you will not be able to call + * some functions of this class if they need information about degrees of + * freedom. These functions are, above all, the + * <tt>get_function_value/gradients/hessians/third_derivatives</tt> + * functions. If you want to call these functions, you have to call the @p + * reinit variants that take iterators into DoFHandler or other DoF handler + * type objects. + * + * @note @p face must be one of @p cell's face iterators. + */ + void + reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell, + const typename Triangulation<dim, spacedim>::face_iterator &face); + /** * Return a reference to this very object. * @@ -3877,6 +3908,17 @@ public: const unsigned int face_no, const unsigned int subface_no); + /** + * Alternative reinitialization function that takes, as arguments, iterators + * to the face and subface instead of their numbers. + */ + template <template <int, int> class DoFHandlerType, bool level_dof_access> + void + reinit(const TriaIterator<DoFCellAccessor<DoFHandlerType<dim, spacedim>, + level_dof_access>> & cell, + const typename Triangulation<dim, spacedim>::face_iterator &face, + const typename Triangulation<dim, spacedim>::face_iterator &subface); + /** * Reinitialize the gradients, Jacobi determinants, etc for the given * subface on a given cell of type "iterator into a Triangulation object", and @@ -3895,6 +3937,30 @@ public: const unsigned int face_no, const unsigned int subface_no); + /** + * Reinitialize the gradients, Jacobi determinants, etc for the given + * subface on a given cell of type "iterator into a Triangulation object", and + * the given finite element. Since iterators into a triangulation alone only + * convey information about the geometry of a cell, but not about degrees of + * freedom possibly associated with this cell, you will not be able to call + * some functions of this class if they need information about degrees of + * freedom. These functions are, above all, the + * <tt>get_function_value/gradients/hessians/third_derivatives</tt> + * functions. If you want to call these functions, you have to call the @p + * reinit variants that take iterators into DoFHandler or other DoF handler + * type objects. + * + * This does the same thing as the previous function but takes iterators + * instead of numbers as arguments. + * + * @note @p face and @p subface must correspond to a face (and a subface of + * that face) of @p cell. + */ + void + reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell, + const typename Triangulation<dim, spacedim>::face_iterator &face, + const typename Triangulation<dim, spacedim>::face_iterator &subface); + /** * Return a reference to this very object. * diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index 4208da6368..914733969f 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -4784,6 +4784,19 @@ FEFaceValues<dim, spacedim>::reinit( +template <int dim, int spacedim> +template <template <int, int> class DoFHandlerType, bool lda> +void +FEFaceValues<dim, spacedim>::reinit( + const TriaIterator<DoFCellAccessor<DoFHandlerType<dim, spacedim>, lda>> &cell, + const typename Triangulation<dim, spacedim>::face_iterator & face) +{ + const auto face_n = cell->face_iterator_to_index(face); + reinit(cell, face_n); +} + + + template <int dim, int spacedim> void FEFaceValues<dim, spacedim>::reinit( @@ -4806,6 +4819,18 @@ FEFaceValues<dim, spacedim>::reinit( +template <int dim, int spacedim> +void +FEFaceValues<dim, spacedim>::reinit( + const typename Triangulation<dim, spacedim>::cell_iterator &cell, + const typename Triangulation<dim, spacedim>::face_iterator &face) +{ + const auto face_n = cell->face_iterator_to_index(face); + reinit(cell, face_n); +} + + + template <int dim, int spacedim> void FEFaceValues<dim, spacedim>::do_reinit(const unsigned int face_no) @@ -4980,6 +5005,21 @@ FESubfaceValues<dim, spacedim>::reinit( +template <int dim, int spacedim> +template <template <int, int> class DoFHandlerType, bool lda> +void +FESubfaceValues<dim, spacedim>::reinit( + const TriaIterator<DoFCellAccessor<DoFHandlerType<dim, spacedim>, lda>> &cell, + const typename Triangulation<dim, spacedim>::face_iterator & face, + const typename Triangulation<dim, spacedim>::face_iterator &subface) +{ + reinit(cell, + cell->face_iterator_to_index(face), + face->child_iterator_to_index(subface)); +} + + + template <int dim, int spacedim> void FESubfaceValues<dim, spacedim>::reinit( @@ -5005,6 +5045,20 @@ FESubfaceValues<dim, spacedim>::reinit( +template <int dim, int spacedim> +void +FESubfaceValues<dim, spacedim>::reinit( + const typename Triangulation<dim, spacedim>::cell_iterator &cell, + const typename Triangulation<dim, spacedim>::face_iterator &face, + const typename Triangulation<dim, spacedim>::face_iterator &subface) +{ + reinit(cell, + cell->face_iterator_to_index(face), + face->child_iterator_to_index(subface)); +} + + + template <int dim, int spacedim> void FESubfaceValues<dim, spacedim>::do_reinit(const unsigned int face_no, diff --git a/source/fe/fe_values.inst.in b/source/fe/fe_values.inst.in index 5eca3410a2..b2cfec8f4f 100644 --- a/source/fe/fe_values.inst.in +++ b/source/fe/fe_values.inst.in @@ -67,13 +67,35 @@ for (dof_handler : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS; DoFCellAccessor<dof_handler<deal_II_dimension, deal_II_space_dimension>, lda>> &, unsigned int); + + template void + FEFaceValues<deal_II_dimension, deal_II_space_dimension>::reinit( + const TriaIterator< + DoFCellAccessor<dof_handler<deal_II_dimension, deal_II_space_dimension>, + lda>> &, + const TriaIterator<TriaAccessor<deal_II_dimension - 1, + deal_II_dimension, + deal_II_space_dimension>> &); + template void FESubfaceValues<deal_II_dimension, deal_II_space_dimension>::reinit( const TriaIterator< DoFCellAccessor<dof_handler<deal_II_dimension, deal_II_space_dimension>, lda>> &, - unsigned int, - unsigned int); + const unsigned int, + const unsigned int); + + template void + FESubfaceValues<deal_II_dimension, deal_II_space_dimension>::reinit( + const TriaIterator< + DoFCellAccessor<dof_handler<deal_II_dimension, deal_II_space_dimension>, + lda>> &, + const TriaIterator<TriaAccessor<deal_II_dimension - 1, + deal_II_dimension, + deal_II_space_dimension>> &, + const TriaIterator<TriaAccessor<deal_II_dimension - 1, + deal_II_dimension, + deal_II_space_dimension>> &); #endif } diff --git a/tests/fe/shapes.h b/tests/fe/shapes.h index dbe5cda337..d61631e06d 100644 --- a/tests/fe/shapes.h +++ b/tests/fe/shapes.h @@ -157,6 +157,17 @@ plot_face_shape_functions(Mapping<dim> & mapping, q, UpdateFlags(uflags | update_quadrature_points)); + // Test the iterator-based initialization routines too: + FEFaceValues<dim> fe2(mapping, + finel, + q, + UpdateFlags(uflags | update_quadrature_points)); + FESubfaceValues<dim> sub2(mapping, + finel, + q, + UpdateFlags(uflags | update_quadrature_points)); + + sprintf(fname, "Face%dd-%s", dim, name); deallog.push(fname); @@ -165,6 +176,7 @@ plot_face_shape_functions(Mapping<dim> & mapping, if (!c->face(f)->has_children()) { fe.reinit(c, f); + fe2.reinit(c, c->face(f)); unsigned int k = 0; for (unsigned int my = 0; my <= ((dim > 2) ? div : 0); ++my) @@ -187,23 +199,47 @@ plot_face_shape_functions(Mapping<dim> & mapping, c) { if (uflags & update_values) - AssertThrow((fe.shape_value(i, k) == - fe.shape_value_component(i, k, c)), - ExcInternalError()); + { + AssertThrow( + (fe.shape_value(i, k) == + fe.shape_value_component(i, k, c)), + ExcInternalError()); + AssertThrow((fe.shape_value(i, k) == + fe2.shape_value(i, k)), + ExcInternalError()); + } if (uflags & update_gradients) - AssertThrow((fe.shape_grad(i, k) == - fe.shape_grad_component(i, k, c)), - ExcInternalError()); + { + AssertThrow( + (fe.shape_grad(i, k) == + fe.shape_grad_component(i, k, c)), + ExcInternalError()); + AssertThrow((fe.shape_grad(i, k) == + fe2.shape_grad(i, k)), + ExcInternalError()); + } if (uflags & update_hessians) - AssertThrow( - (fe.shape_hessian(i, k) == - fe.shape_hessian_component(i, k, c)), - ExcInternalError()); + { + AssertThrow( + (fe.shape_hessian(i, k) == + fe.shape_hessian_component(i, k, c)), + ExcInternalError()); + AssertThrow((fe.shape_hessian(i, k) == + fe2.shape_hessian(i, k)), + ExcInternalError()); + } if (uflags & update_3rd_derivatives) - AssertThrow( - (fe.shape_3rd_derivative(i, k) == - fe.shape_3rd_derivative_component(i, k, c)), - ExcInternalError()); + { + AssertThrow( + (fe.shape_3rd_derivative(i, k) == + fe.shape_3rd_derivative_component(i, + k, + c)), + ExcInternalError()); + AssertThrow((fe.shape_3rd_derivative(i, k) == + fe2.shape_3rd_derivative(i, k)), + ExcInternalError()); + } } else { @@ -239,6 +275,7 @@ plot_face_shape_functions(Mapping<dim> & mapping, ++s) { sub.reinit(c, f, s); + sub2.reinit(c, c->face(f), c->face(f)->child(s)); unsigned int k = 0; for (unsigned int my = 0; my <= ((dim > 2) ? div : 0); ++my) @@ -267,23 +304,29 @@ plot_face_shape_functions(Mapping<dim> & mapping, const double v1 = sub.shape_value(i, k), v2 = sub.shape_value_component( - i, k, c); + i, k, c), + v3 = sub2.shape_value(i, k); Assert(v1 == v2, ExcInternalError()); + Assert(v1 == v3, ExcInternalError()); } if (uflags & update_gradients) { const Tensor<1, dim> g1 = sub.shape_grad(i, k), - g2 = sub.shape_grad_component(i, k, c); + g2 = sub.shape_grad_component(i, k, c), + g3 = sub2.shape_grad(i, k); Assert(g1 == g2, ExcInternalError()); + Assert(g1 == g3, ExcInternalError()); } if (uflags & update_hessians) { const Tensor<2, dim> s1 = sub.shape_hessian(i, k), s2 = - sub.shape_hessian_component(i, k, c); + sub.shape_hessian_component(i, k, c), + s3 = sub2.shape_hessian(i, k); Assert(s1 == s2, ExcInternalError()); + Assert(s1 == s3, ExcInternalError()); } if (uflags & update_3rd_derivatives) { @@ -292,8 +335,11 @@ plot_face_shape_functions(Mapping<dim> & mapping, t2 = sub.shape_3rd_derivative_component(i, k, - c); + c), + t3 = sub2.shape_3rd_derivative(i, k); + Assert(t1 == t2, ExcInternalError()); + Assert(t1 == t3, ExcInternalError()); } } else -- 2.39.5