From 75a865f099f3ea7303739084841a417a39ab0654 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 8 Mar 2001 15:18:12 +0000 Subject: [PATCH] Move unit support points on cells and faces to base class. Initialize them in derived classes if necessary. git-svn-id: https://svn.dealii.org/trunk@4173 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 36 --- deal.II/deal.II/include/fe/fe_base.h | 170 +++++++++++++-- deal.II/deal.II/include/fe/fe_dgq.h | 41 ---- deal.II/deal.II/include/fe/fe_q.h | 58 ++--- deal.II/deal.II/include/fe/fe_system.h | 60 ++--- deal.II/deal.II/source/fe/fe.cc | 70 ++++-- deal.II/deal.II/source/fe/fe_dgq.cc | 96 ++++---- deal.II/deal.II/source/fe/fe_q.cc | 95 ++++---- deal.II/deal.II/source/fe/fe_system.cc | 241 ++++++++++++--------- deal.II/deal.II/source/fe/fe_tools.cc | 4 +- deal.II/deal.II/source/numerics/vectors.cc | 9 +- 11 files changed, 482 insertions(+), 398 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 2a86519a75..6e4edfc9f7 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -67,42 +67,6 @@ class FiniteElement : public FiniteElementBase * are deleted properly. */ virtual ~FiniteElement (); - - /** - * Return the support points of the - * trial functions on the unit cell. - * - * The order of points in the - * array matches that returned by - * the #cell->get_dof_indices# - * function, but: - * - * If the shape functions are not - * Lagrangian interpolants at - * some points, the size of the - * array will be zero after - * calling this function. This is - * the standard behavior, if the - * function is not overloaded. - */ - virtual void get_unit_support_points (std::vector > &points) const; - - /** - * Return the support points of - * the trial functions on the - * first face of the unit cell. - * - * The order of points in - * the array matches that returned by - * the #cell->get_dof_indices# function. - * - * If the shape functions are not - * Lagrangian interpolants at some - * points, the size of the array - * will be zero. This is the standard behavior, - * if the function is not overloaded. - */ - virtual void get_unit_face_support_points (std::vector > &points) const; /** * Number of base elements in a mixed diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index c8cca72bf2..518c86c781 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -375,18 +375,6 @@ class FiniteElementBase : public Subscriptor, */ FiniteElementBase (const FiniteElementData &fe_data, const std::vector &restriction_is_additive_flags); - - /** - * Compute second derivatives by - * finite differences of - * gradients. - */ - void compute_2nd (const Mapping &mapping, - const DoFHandler::cell_iterator &cell, - const unsigned int offset, - typename Mapping::InternalDataBase &mapping_internal, - InternalDataBase &fe_internal, - FEValuesData &data) const; /** * Projection from a fine grid @@ -503,8 +491,8 @@ class FiniteElementBase : public Subscriptor, /** * Compute system index from components. */ - unsigned int component_to_system_index (unsigned int component, - unsigned int component_index) const; + unsigned int component_to_system_index (const unsigned int component, + const unsigned int component_index) const; /** * Compute component and index from @@ -515,13 +503,13 @@ class FiniteElementBase : public Subscriptor, * component. */ std::pair - system_to_component_index (unsigned int index) const; + system_to_component_index (const unsigned int index) const; /** * Compute system index from components on a face. */ - unsigned int face_component_to_system_index (unsigned int component, - unsigned int component_index) const; + unsigned int face_component_to_system_index (const unsigned int component, + const unsigned int component_index) const; /** * Compute component and index from system @@ -532,7 +520,7 @@ class FiniteElementBase : public Subscriptor, * component. */ std::pair - face_system_to_component_index (unsigned int index) const; + face_system_to_component_index (const unsigned int index) const; /** * The base element establishing a @@ -556,6 +544,109 @@ class FiniteElementBase : public Subscriptor, */ bool restriction_is_additive (const unsigned int component) const; + /** + * Return the support points of + * the trial functions on the + * unit cell, if the derived + * finite element defines some. + * Finite elements that allow + * some kind of interpolation + * operation usually have support + * points. On the other hand, + * elements that define their + * degrees of freedom by, for + * example, moments on faces, or + * as derivatives, don't have + * support points. In that case, + * the returned field is empty. + * + * If the finite element defines + * support points, then their + * number equals the number of + * degrees of freedom of the + * element. The order of points + * in the array matches that + * returned by the + * @p{cell->get_dof_indices} + * function. + */ + const std::vector > & get_unit_support_points () const; + + /** + * Return whether a finite + * element has defined support + * points. If the result is true, + * then a call to the + * @p{get_unit_support_points} + * yields a non-empty array. + */ + bool has_support_points () const; + + /** + * Return the support points of + * the trial functions on the + * unit face, if the derived + * finite element defines some. + * Finite elements that allow + * some kind of interpolation + * operation usually have support + * points. On the other hand, + * elements that define their + * degrees of freedom by, for + * example, moments on faces, or + * as derivatives, don't have + * support points. In that case, + * the returned field is empty + * + * Note that elements that have + * support points need not + * necessarily have some on the + * faces, even if the + * interpolation points are + * located physically on a + * face. For example, the + * discontinuous elements have + * interpolation points on the + * vertices, and for higher + * degree elements also on the + * faces, but they are not + * defined to be on faces since + * in that case degrees of + * freedom from both sides of a + * face (or from all adjacent + * elements to a vertex) would be + * identified with each other, + * which is not what we would + * like to have). Logically, + * these degrees of freedom are + * therefore defined to belong to + * the cell, rather than the face + * or vertex. In that case, the + * returned element would + * therefore have length zero. + * + * If the finite element defines + * support points, then their + * number equals the number of + * degrees of freedom on the face + * (@p{dofs_per_face}). The order + * of points in the array matches + * that returned by the + * @p{cell->get_dof_indices} + * function. + */ + const std::vector > & get_unit_face_support_points () const; + + /** + * Return whether a finite + * element has defined support + * points on faces. If the result + * is true, then a call to the + * @p{get_unit_support_points} + * yields a non-empty array. + */ + bool has_face_support_points () const; + /** * Determine an estimate for the * memory consumption (in bytes) @@ -721,6 +812,49 @@ class FiniteElementBase : public Subscriptor, */ const std::vector restriction_is_additive_flags; + /** + * List of support points on the + * unit cell, in case the finite + * element has any. The + * constructor leaves this field + * empty, derived classes may + * write in some contents. + * + * Finite elements that allow + * some kind of interpolation + * operation usually have support + * points. On the other hand, + * elements that define their + * degrees of freedom by, for + * example, moments on faces, or + * as derivatives, don't have + * support points. In that case, + * this field remains empty. + */ + std::vector > unit_support_points; + + /** + * Same for the faces. See the + * description of the + * @p{get_unit_face_support_points} + * function for a discussion of + * what contributes a face + * support point. + */ + std::vector > unit_face_support_points; + + /** + * Compute second derivatives by + * finite differences of + * gradients. + */ + void compute_2nd (const Mapping &mapping, + const DoFHandler::cell_iterator &cell, + const unsigned int offset, + typename Mapping::InternalDataBase &mapping_internal, + InternalDataBase &fe_internal, + FEValuesData &data) const; + /** * Allow the FESystem class to * access the restriction and diff --git a/deal.II/deal.II/include/fe/fe_dgq.h b/deal.II/deal.II/include/fe/fe_dgq.h index aa3ba2f922..dedbb0b76d 100644 --- a/deal.II/deal.II/include/fe/fe_dgq.h +++ b/deal.II/deal.II/include/fe/fe_dgq.h @@ -38,40 +38,6 @@ class FE_DGQ : public FiniteElement */ ~FE_DGQ (); - /** - * Return the support points of the - * trial functions on the unit cell. - * - * The order of points in - * the array matches that returned by - * the #cell->get_dof_indices# function, but: - * - * If the shape functions are not - * Lagrangian interpolants at some - * points, the size of the array - * will be zero. This is the standard behavior, - * if the function is not overloaded. - */ - virtual void get_unit_support_points (std::vector > &) const; - - /** - * Return the support points of - * the trial functions on the - * first face of the unit cell. - * - * The order of points in - * the array matches that returned by - * the #cell->get_dof_indices# function. - * - * If the shape functions are not - * Lagrangian interpolants at some - * points, the size of the array - * will be zero. This is the standard behavior, - * if the function is not overloaded. - */ - virtual void get_unit_face_support_points (std::vector > &) const; - - /** * Return the polynomial degree * of this finite element, @@ -227,13 +193,6 @@ class FE_DGQ : public FiniteElement */ virtual UpdateFlags update_each (UpdateFlags flags) const; - /** - * Compute support points, only - * for @p{degree>0}. - */ - static void compute_support_points (std::vector > &support_points, - const unsigned int degree); - /** * Compute renumbering for rotation * of degrees of freedom. diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index 237ad8a641..04206bf33a 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -40,44 +40,6 @@ class FE_Q : public FiniteElement */ ~FE_Q (); - /** - * Return the support points of - * the trial functions on the - * unit cell. - * - * The order of points in the - * array matches that returned by - * the @p{cell->get_dof_indices} - * function, but: - * - * If the shape functions are not - * Lagrangian interpolants at - * some points, the size of the - * array will be zero. This is - * the standard behavior, if the - * function is not overloaded. - */ - virtual void get_unit_support_points (std::vector > &) const; - - /** - * Return the support points of - * the trial functions on the - * first face of the unit cell. - * - * The order of points in the - * array matches that returned by - * the @p{cell->get_dof_indices} - * function. - * - * If the shape functions are not - * Lagrangian interpolants at - * some points, the size of the - * array will be zero. This is - * the standard behavior, if the - * function is not overloaded. - */ - virtual void get_unit_face_support_points (std::vector > &) const; - //TODO: make get_renumber private or move it some other place /** * Read-only access to the @@ -304,13 +266,23 @@ class FE_Q : public FiniteElement std::vector &numbering); /** - * Compute support points. + * Initialize the + * @p{unit_support_points} field + * of the @ref{FiniteElementBase} + * class. Called from the + * constructor. */ - static void compute_support_points (const unsigned int degree, - const std::vector &renumber, - std::vector > &support_points); - + void initialize_unit_support_points (); + /** + * Initialize the + * @p{unit_face_support_points} field + * of the @ref{FiniteElementBase} + * class. Called from the + * constructor. + */ + void initialize_unit_face_support_points (); + /** * Compute flags for initial * update only. diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index a287fab700..a6246becc2 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -170,47 +170,7 @@ class FESystem : public FiniteElement * pointers to their base class, * rather than the class itself. */ - virtual unsigned int memory_consumption () const; - - /** - * Return the support points of - * the trial functions on the - * unit cell. - * - * The order of points in the - * array matches that returned by - * the @p{cell->get_dof_indices} - * function, but: - * - * If one of the base elements - * has no support points, then it - * makes no sense to define - * support points for the - * composed element, so return an - * empty array to demonstrate - * that fact. - */ - virtual void get_unit_support_points (typename std::vector > &) const; - - /** - * Return the support points of - * the trial functions on the - * first face of the unit cell. - * - * The order of points in the - * array matches that returned by - * the @p{cell->get_dof_indices} - * function, but: - * - * If one of the base elements - * has no support points, then it - * makes no sense to define - * support points for the - * composed element, so return an - * empty array to demonstrate - * that fact. - */ - virtual void get_unit_face_support_points (typename std::vector > &) const; + virtual unsigned int memory_consumption () const; protected: /** @@ -340,6 +300,24 @@ class FESystem : public FiniteElement typename std::vector base_elements; + /** + * Initialize the + * @p{unit_support_points} field + * of the @ref{FiniteElementBase} + * class. Called from the + * constructor. + */ + void initialize_unit_support_points (); + + /** + * Initialize the + * @p{unit_face_support_points} field + * of the @ref{FiniteElementBase} + * class. Called from the + * constructor. + */ + void initialize_unit_face_support_points (); + /** * Helper function used in the constructor: * take a @p{FiniteElementData} object diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 01eda10290..67abddfe97 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -220,6 +220,56 @@ bool FiniteElementBase::operator == (const FiniteElementBase &f) const +template +const std::vector > & +FiniteElementBase::get_unit_support_points () const +{ + // a finite element may define + // support points, but only if + // there are as many as there are + // degrees of freedom + Assert ((unit_support_points.size() == 0) || + (unit_support_points.size() == dofs_per_cell), + ExcInternalError()); + return unit_support_points; +}; + + + +template +bool +FiniteElementBase::has_support_points () const +{ + return (unit_support_points.size() != 0); +}; + + + +template +const std::vector > & +FiniteElementBase::get_unit_face_support_points () const +{ + // a finite element may define + // support points, but only if + // there are as many as there are + // degrees of freedom on a face + Assert ((unit_face_support_points.size() == 0) || + (unit_face_support_points.size() == dofs_per_face), + ExcInternalError()); + return unit_face_support_points; +}; + + + +template +bool +FiniteElementBase::has_face_support_points () const +{ + return (unit_face_support_points.size() != 0); +}; + + + template unsigned int FiniteElementBase::memory_consumption () const @@ -333,26 +383,6 @@ FiniteElement::~FiniteElement () -template -void -FiniteElement::get_unit_support_points (std::vector > &points) const -{ - // default implementation - points.resize(0); -} - - - -template -void -FiniteElement::get_unit_face_support_points (std::vector > &points) const -{ - // default implementation - points.resize(0); -} - - - template Mapping::InternalDataBase* FiniteElement::get_face_data (const UpdateFlags flags, diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index e9d98c7a3f..b9da2f48ef 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -161,6 +161,46 @@ FE_DGQ::FE_DGQ (unsigned int degree) // matrix undefined, set size to zero for (unsigned int i=0;i::children_per_cell;++i) restriction[i].reinit(0); + + + // finally fill in support points + if (degree == 0) + { + // constant elements, take + // midpoint + unit_support_points.resize(1); + for (unsigned int i=0; i p; + + unsigned int k=0; + for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz) + for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy) + for (unsigned int ix=0; ix<=degree; ++ix) + { + p(0) = ix * step; + if (dim>1) + p(1) = iy * step; + if (dim>2) + p(2) = iz * step; + + unit_support_points[k++] = p; + }; + }; + + // note: no face support points for + // DG elements }; @@ -185,31 +225,6 @@ FE_DGQ::clone() const -template -void -FE_DGQ::get_unit_support_points (std::vector > &points) const -{ - if (degree>0) - compute_support_points (points, degree); - else - { - // constant elements, take - // midpoint - points.resize(1); - for (unsigned int i=0; i -void -FE_DGQ::get_unit_face_support_points (std::vector > &points) const -{ - // no face support points - points.resize(0); -} - //---------------------------------------------------------------------- // Auxiliary functions @@ -257,37 +272,6 @@ FE_DGQ::update_each (const UpdateFlags flags) const } -template -void -FE_DGQ::compute_support_points (std::vector > &support_points, - const unsigned int degree) -{ - Assert(degree>0, ExcInternalError()); - // number of points: (degree+1)^dim - unsigned int n = degree+1; - for (unsigned int i=1;i p; - - unsigned int k=0; - for (unsigned int iz=0;iz <= ((dim>2) ? degree : 0) ; ++iz) - for (unsigned int iy=0;iy <= ((dim>1) ? degree : 0) ; ++iy) - for (unsigned int ix=0;ix<=degree;++ix) - { - p(0) = ix * step; - if (dim>1) - p(1) = iy * step; - if (dim>2) - p(2) = iz * step; - - support_points[k++] = p; - } -} - template void diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 90b17638d2..24a8ab12f8 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -570,7 +570,13 @@ FE_Q::FE_Q (const unsigned int degree) default: Assert (false,ExcNotImplemented()); } -} + + // finally fill in support points + // on cell and face + initialize_unit_support_points (); + initialize_unit_face_support_points (); +}; + template @@ -593,71 +599,82 @@ FE_Q::clone() const -template -void -FE_Q::get_unit_support_points (std::vector > &points) const -{ - compute_support_points (degree, renumber, points); -} +//---------------------------------------------------------------------- +// Auxiliary functions +//---------------------------------------------------------------------- + - template -void -FE_Q::get_unit_face_support_points (std::vector > &points) const +void FE_Q::initialize_unit_support_points () { - FE_Q::compute_support_points (degree, face_renumber, points); -} + // number of points: (degree+1)^dim + unsigned int n = degree+1; + for (unsigned int i=1; i p; + + unsigned int k=0; + for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz) + for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy) + for (unsigned int ix=0; ix<=degree; ++ix) + { + p(0) = ix * step; + if (dim>1) + p(1) = iy * step; + if (dim>2) + p(2) = iz * step; + + unit_support_points[renumber[k++]] = p; + }; +}; #if deal_II_dimension == 1 template <> -void -FE_Q<1>::get_unit_face_support_points (std::vector > &points) const +void FE_Q<1>::initialize_unit_face_support_points () { - points.resize(0); -} + // no faces in 1d, so nothing to do +}; #endif - -//---------------------------------------------------------------------- -// Auxiliary functions -//---------------------------------------------------------------------- - - template -void -FE_Q::compute_support_points (const unsigned int degree, - const std::vector &renumber, - std::vector > &support_points) +void FE_Q::initialize_unit_face_support_points () { - // number of points: (degree+1)^dim + const unsigned int codim = dim-1; + + // number of points: (degree+1)^codim unsigned int n = degree+1; - for (unsigned int i=1; i p; + Point p; unsigned int k=0; - for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz) - for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy) + for (unsigned int iz=0; iz <= ((codim>2) ? degree : 0) ; ++iz) + for (unsigned int iy=0; iy <= ((codim>1) ? degree : 0) ; ++iy) for (unsigned int ix=0; ix<=degree; ++ix) - { + { p(0) = ix * step; - if (dim>1) + if (codim>1) p(1) = iy * step; - if (dim>2) + if (codim>2) p(2) = iz * step; - support_points[renumber[k++]] = p; - } -} + unit_face_support_points[face_renumber[k++]] = p; + }; +}; + template diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 140c00d0ab..e7f923b145 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -230,104 +230,6 @@ FESystem::clone() const -template -void -FESystem:: -get_unit_support_points (typename std::vector > &unit_support_points) const -{ - // generate unit support points - // from unit support points of sub - // elements - unit_support_points.resize(dofs_per_cell); - - typename std::vector > base_unit_support_points (base_element(0).dofs_per_cell); - unsigned int comp = 0; - for (unsigned int base_el=0; base_el -void -FESystem:: -get_unit_face_support_points (typename std::vector > &unit_support_points) const -{ - // generate unit face support points - // from unit support points of sub - // elements - unit_support_points.resize(dofs_per_face); - - typename std::vector > - base_unit_support_points (base_element(0).dofs_per_cell); - - unsigned int comp = 0; - for (unsigned int base_el=0 ; base_el UpdateFlags FESystem::update_once (const UpdateFlags flags) const @@ -1181,6 +1083,11 @@ void FESystem::initialize () //FESystem(FE_Q (3), 2) and for FESystem(FE_Q (1), 2, //FE_Q (3), 1) and for FESystem(FE_Q (4), 2)) build_interface_constraints (); + + // finally fill in support points + // on cell and face + initialize_unit_support_points (); + initialize_unit_face_support_points (); }; @@ -1201,6 +1108,144 @@ FESystem::multiply_dof_numbers (const FiniteElementData &fe_data, + +template +void +FESystem:: +initialize_unit_support_points () +{ + // if one of the base elements + // has no support points, then + // it makes no sense to define + // support points for the + // composed element, so return + // an empty array to + // demonstrate that + // fact + for (unsigned int base_el=0 ; base_el > + & base_unit_support_points = base_element(base_el).get_unit_support_points (); + + // otherwise distribute the + // support points of this base + // element to all degrees of + // freedom contributed by this + // base element + Assert(base_unit_support_points.size()==base_element_dofs_per_cell, + ExcInternalError()); + for (unsigned int n=0; n +void +FESystem<1>:: +initialize_unit_face_support_points () +{ + // no faces no work +}; + +#endif + + +template +void +FESystem:: +initialize_unit_face_support_points () +{ + // if one of the base elements + // has no support points, then + // it makes no sense to define + // support points for the + // composed element, so return + // an empty array to + // demonstrate that fact (note + // that we ask whether the base + // element has no support + // points at all, not only none + // on the face!) + for (unsigned int base_el=0 ; base_el > & + base_unit_support_points = base_element(base_el).get_unit_face_support_points (); + + // distribute the support + // points of this base element + // to all degrees of freedom + // contributed by this base + // element + Assert(base_unit_support_points.size()==base_element_dofs_per_face, + ExcNotImplemented()); + for (unsigned int n=0; n FiniteElementData FESystem::multiply_dof_numbers (const FiniteElementData &fe1, diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 4d7ce88639..09b897d8a7 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -50,8 +50,8 @@ void FETools::get_interpolation_matrix(const FiniteElement &fe1, // the unit support points of the // fe2 element. std::vector phantom_weights(fe2.dofs_per_cell,1.); - std::vector > fe2_support_points (fe2.dofs_per_cell); - fe2.get_unit_support_points (fe2_support_points); + const typename std::vector > & + fe2_support_points = fe2.get_unit_support_points (); Quadrature fe2_support_points_quadrature(fe2_support_points, phantom_weights); diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 69f15b5e67..98864bfb46 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -91,8 +91,8 @@ void VectorTools::interpolate (const DoFHandler &dof, // avoided to evaluate // the vectorfunction multiply at // the same point on a cell. - std::vector > unit_support_points; - fe.get_unit_support_points(unit_support_points); + const typename std::vector > & + unit_support_points = fe.get_unit_support_points(); Assert (unit_support_points.size() != 0, ExcNonInterpolatingFE()); @@ -630,8 +630,9 @@ VectorTools::interpolate_boundary_values (const DoFHandler &dof, // support points. this wil be used // to obtain the quadrature points // on the real cell's face - typename std::vector > unit_support_points; - fe.get_unit_face_support_points(unit_support_points); + const typename std::vector > + & unit_support_points = fe.get_unit_face_support_points(); + // check whether there are support // points on the face, if not, then // this FE does not allow to be -- 2.39.5