* 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<Point<dim> > &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<Point<dim-1> > &points) const;
/**
* Number of base elements in a mixed
*/
FiniteElementBase (const FiniteElementData<dim> &fe_data,
const std::vector<bool> &restriction_is_additive_flags);
-
- /**
- * Compute second derivatives by
- * finite differences of
- * gradients.
- */
- void compute_2nd (const Mapping<dim> &mapping,
- const DoFHandler<dim>::cell_iterator &cell,
- const unsigned int offset,
- typename Mapping<dim>::InternalDataBase &mapping_internal,
- InternalDataBase &fe_internal,
- FEValuesData<dim> &data) const;
/**
* Projection from a fine grid
/**
* 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
* component.
*/
std::pair<unsigned int,unsigned int>
- 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
* component.
*/
std::pair<unsigned int,unsigned int>
- face_system_to_component_index (unsigned int index) const;
+ face_system_to_component_index (const unsigned int index) const;
/**
* The base element establishing a
*/
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<Point<dim> > & 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<Point<dim-1> > & 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)
*/
const std::vector<bool> 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<Point<dim> > 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<Point<dim-1> > unit_face_support_points;
+
+ /**
+ * Compute second derivatives by
+ * finite differences of
+ * gradients.
+ */
+ void compute_2nd (const Mapping<dim> &mapping,
+ const DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int offset,
+ typename Mapping<dim>::InternalDataBase &mapping_internal,
+ InternalDataBase &fe_internal,
+ FEValuesData<dim> &data) const;
+
/**
* Allow the FESystem class to
* access the restriction and
*/
~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<Point<dim> > &) 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<Point<dim-1> > &) const;
-
-
/**
* Return the polynomial degree
* of this finite element,
*/
virtual UpdateFlags update_each (UpdateFlags flags) const;
- /**
- * Compute support points, only
- * for @p{degree>0}.
- */
- static void compute_support_points (std::vector<Point<dim> > &support_points,
- const unsigned int degree);
-
/**
* Compute renumbering for rotation
* of degrees of freedom.
*/
~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<Point<dim> > &) 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<Point<dim-1> > &) const;
-
//TODO: make get_renumber private or move it some other place
/**
* Read-only access to the
std::vector<unsigned int> &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<unsigned int> &renumber,
- std::vector<Point<dim> > &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.
* 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<Point<dim> > &) 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<Point<dim-1> > &) const;
+ virtual unsigned int memory_consumption () const;
protected:
/**
typename std::vector<ElementPair> 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
+template <int dim>
+const std::vector<Point<dim> > &
+FiniteElementBase<dim>::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 <int dim>
+bool
+FiniteElementBase<dim>::has_support_points () const
+{
+ return (unit_support_points.size() != 0);
+};
+
+
+
+template <int dim>
+const std::vector<Point<dim-1> > &
+FiniteElementBase<dim>::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 <int dim>
+bool
+FiniteElementBase<dim>::has_face_support_points () const
+{
+ return (unit_face_support_points.size() != 0);
+};
+
+
+
template <int dim>
unsigned int
FiniteElementBase<dim>::memory_consumption () const
-template <int dim>
-void
-FiniteElement<dim>::get_unit_support_points (std::vector<Point<dim> > &points) const
-{
- // default implementation
- points.resize(0);
-}
-
-
-
-template <int dim>
-void
-FiniteElement<dim>::get_unit_face_support_points (std::vector<Point<dim-1> > &points) const
-{
- // default implementation
- points.resize(0);
-}
-
-
-
template <int dim>
Mapping<dim>::InternalDataBase*
FiniteElement<dim>::get_face_data (const UpdateFlags flags,
// matrix undefined, set size to zero
for (unsigned int i=0;i<GeometryInfo<dim>::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<dim; ++i)
+ unit_support_points[0](i) = 0.5;
+ }
+ else
+ {
+ // number of points: (degree+1)^dim
+ unsigned int n = degree+1;
+ for (unsigned int i=1; i<dim; ++i)
+ n *= degree+1;
+
+ unit_support_points.resize(n);
+
+ const double step = 1./degree;
+ Point<dim> 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
};
-template <int dim>
-void
-FE_DGQ<dim>::get_unit_support_points (std::vector<Point<dim> > &points) const
-{
- if (degree>0)
- compute_support_points (points, degree);
- else
- {
- // constant elements, take
- // midpoint
- points.resize(1);
- for (unsigned int i=0; i<dim; ++i)
- points[0](i)=0.5;
- }
-}
-
-
-template <int dim>
-void
-FE_DGQ<dim>::get_unit_face_support_points (std::vector<Point<dim-1> > &points) const
-{
- // no face support points
- points.resize(0);
-}
-
//----------------------------------------------------------------------
// Auxiliary functions
}
-template <int dim>
-void
-FE_DGQ<dim>::compute_support_points (std::vector<Point<dim> > &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<dim;++i)
- n *= degree+1;
-
- support_points.resize(n);
-
- const double step = 1./degree;
- Point<dim> 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 <int dim>
void
default:
Assert (false,ExcNotImplemented());
}
-}
+
+ // finally fill in support points
+ // on cell and face
+ initialize_unit_support_points ();
+ initialize_unit_face_support_points ();
+};
+
template <int dim>
-template <int dim>
-void
-FE_Q<dim>::get_unit_support_points (std::vector<Point<dim> > &points) const
-{
- compute_support_points (degree, renumber, points);
-}
+//----------------------------------------------------------------------
+// Auxiliary functions
+//----------------------------------------------------------------------
+
-
template <int dim>
-void
-FE_Q<dim>::get_unit_face_support_points (std::vector<Point<dim-1> > &points) const
+void FE_Q<dim>::initialize_unit_support_points ()
{
- FE_Q<dim-1>::compute_support_points (degree, face_renumber, points);
-}
+ // number of points: (degree+1)^dim
+ unsigned int n = degree+1;
+ for (unsigned int i=1; i<dim; ++i)
+ n *= degree+1;
+
+ unit_support_points.resize(n);
+
+ const double step = 1./degree;
+ Point<dim> 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<Point<0> > &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 <int dim>
-void
-FE_Q<dim>::compute_support_points (const unsigned int degree,
- const std::vector<unsigned int> &renumber,
- std::vector<Point<dim> > &support_points)
+void FE_Q<dim>::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<dim; ++i)
+ for (unsigned int i=1; i<codim; ++i)
n *= degree+1;
-
- support_points.resize(n);
+
+ unit_face_support_points.resize(n);
const double step = 1./degree;
- Point<dim> p;
+ Point<codim> 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 <int dim>
-template <int dim>
-void
-FESystem<dim>::
-get_unit_support_points (typename std::vector<Point<dim> > &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<Point<dim> > base_unit_support_points (base_element(0).dofs_per_cell);
- unsigned int comp = 0;
- for (unsigned int base_el=0; base_el<n_base_elements(); ++base_el)
- {
- const unsigned int
- base_element_dofs_per_cell = base_element(base_el).dofs_per_cell;
-
- base_element(base_el).get_unit_support_points (base_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
- if (base_unit_support_points.size()==0)
- {
- unit_support_points.resize(0);
- return;
- }
-
- // 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<element_multiplicity(base_el); ++n, ++comp)
- for (unsigned int i=0; i<base_element_dofs_per_cell; ++i)
- unit_support_points[component_to_system_index(comp,i)]
- = base_unit_support_points[i];
- }
-}
-
-
-
-template <int dim>
-void
-FESystem<dim>::
-get_unit_face_support_points (typename std::vector<Point<dim-1> > &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<Point<dim-1> >
- base_unit_support_points (base_element(0).dofs_per_cell);
-
- unsigned int comp = 0;
- for (unsigned int base_el=0 ; base_el<n_base_elements(); ++base_el)
- {
- const unsigned int
- base_element_dofs_per_face = base_element(base_el).dofs_per_face;
-
- base_element(base_el).get_unit_face_support_points (base_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
- if (base_unit_support_points.size()==0)
- {
- unit_support_points.resize(0);
- return;
- }
-
- // 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_face,
- ExcNotImplemented());
- for (unsigned int n=0; n<element_multiplicity(base_el); ++n, ++comp)
- for (unsigned int i=0; i<base_element_dofs_per_face; ++i)
- unit_support_points[face_component_to_system_index(comp,i)]
- = base_unit_support_points[i];
- }
-}
-
-
-
template <int dim>
UpdateFlags
FESystem<dim>::update_once (const UpdateFlags flags) const
//FESystem(FE_Q<dim> (3), 2) and for FESystem<dim>(FE_Q<dim> (1), 2,
//FE_Q<dim> (3), 1) and for FESystem<dim>(FE_Q<dim> (4), 2))
build_interface_constraints ();
+
+ // finally fill in support points
+ // on cell and face
+ initialize_unit_support_points ();
+ initialize_unit_face_support_points ();
};
+
+template <int dim>
+void
+FESystem<dim>::
+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<n_base_elements(); ++base_el)
+ if (!base_element(base_el).has_support_points())
+ {
+ unit_support_points.resize(0);
+ return;
+ };
+
+ // generate unit support points
+ // from unit support points of sub
+ // elements
+ unit_support_points.resize(dofs_per_cell);
+
+ unsigned int comp = 0;
+ for (unsigned int base_el=0; base_el<n_base_elements(); ++base_el)
+ {
+ // we know that there are
+ // support points on the cell,
+ // so collect them
+ const unsigned int
+ base_element_dofs_per_cell = base_element(base_el).dofs_per_cell;
+
+ const typename std::vector<Point<dim> >
+ & 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<element_multiplicity(base_el); ++n, ++comp)
+ for (unsigned int i=0; i<base_element_dofs_per_cell; ++i)
+ unit_support_points[component_to_system_index(comp,i)]
+ = base_unit_support_points[i];
+ }
+}
+
+
+#if deal_II_dimension == 1
+
+template <>
+void
+FESystem<1>::
+initialize_unit_face_support_points ()
+{
+ // no faces no work
+};
+
+#endif
+
+
+template <int dim>
+void
+FESystem<dim>::
+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<n_base_elements(); ++base_el)
+ if (!base_element(base_el).has_support_points())
+ {
+ unit_face_support_points.resize(0);
+ return;
+ };
+
+
+ // generate unit face support points
+ // from unit support points of sub
+ // elements
+ unit_face_support_points.resize(dofs_per_face);
+
+ unsigned int comp = 0;
+ for (unsigned int base_el=0 ; base_el<n_base_elements(); ++base_el)
+ {
+ // in some cases, finite
+ // elements have support points
+ // (we have made sure that they
+ // have above) but don't have
+ // any on the face (e.g. DG
+ // elements). in that case,
+ // don't even bother with this
+ // base element and directly go
+ // to the next one:
+ if (!base_element(base_el).has_face_support_points())
+ {
+ comp += element_multiplicity(base_el);
+ continue;
+ };
+
+ // otherwise, we know that
+ // there are support points on
+ // the face, so collect them
+ const unsigned int
+ base_element_dofs_per_face = base_element(base_el).dofs_per_face;
+
+ const typename std::vector<Point<dim-1> > &
+ 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<element_multiplicity(base_el); ++n, ++comp)
+ for (unsigned int i=0; i<base_element_dofs_per_face; ++i)
+ unit_face_support_points[face_component_to_system_index(comp,i)]
+ = base_unit_support_points[i];
+ }
+}
+
+
+
template <int dim>
FiniteElementData<dim>
FESystem<dim>::multiply_dof_numbers (const FiniteElementData<dim> &fe1,
// the unit support points of the
// fe2 element.
std::vector<double> phantom_weights(fe2.dofs_per_cell,1.);
- std::vector<Point<dim> > fe2_support_points (fe2.dofs_per_cell);
- fe2.get_unit_support_points (fe2_support_points);
+ const typename std::vector<Point<dim> > &
+ fe2_support_points = fe2.get_unit_support_points ();
Quadrature<dim> fe2_support_points_quadrature(fe2_support_points,
phantom_weights);
// avoided to evaluate
// the vectorfunction multiply at
// the same point on a cell.
- std::vector<Point<dim> > unit_support_points;
- fe.get_unit_support_points(unit_support_points);
+ const typename std::vector<Point<dim> > &
+ unit_support_points = fe.get_unit_support_points();
Assert (unit_support_points.size() != 0,
ExcNonInterpolatingFE());
// support points. this wil be used
// to obtain the quadrature points
// on the real cell's face
- typename std::vector<Point<dim-1> > unit_support_points;
- fe.get_unit_face_support_points(unit_support_points);
+ const typename std::vector<Point<dim-1> >
+ & 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