#include <deal.II/base/config.h>
#include <deal.II/base/tensor_product_polynomials.h>
+#include <deal.II/base/polynomial_space.h>
#include <deal.II/fe/fe_poly_face.h>
DEAL_II_NAMESPACE_OPEN
/**
* Return whether this element implements its hanging node constraints in
* the new way, which has to be used to make elements "hp compatible".
+ */
+ virtual bool hp_constraints_are_implemented () const;
+
+ /**
+ * Return whether this element dominates the one given as argument when they
+ * meet at a common face, whether it is the other way around, whether
+ * neither dominates, or if either could dominate.
+ *
+ * For a definition of domination, see FiniteElementBase::Domination and in
+ * particular the @ref hp_paper "hp paper".
+ */
+ virtual
+ FiniteElementDomination::Domination
+ compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const;
+
+private:
+ /**
+ * Return vector with dofs per vertex, line, quad, hex.
+ */
+ static std::vector<unsigned int> get_dpo_vector (const unsigned int deg);
+};
+
+
+
+
+/**
+ * A finite element, which is a Legendre on each face (i.e., FE_DGP)
+ * and undefined in the interior of the cells. The basis functions on
+ * the faces are from Polynomials::Legendre.
+ *
+ * @note Since these are only finite elements on faces, only
+ * FEFaceValues and FESubfaceValues will be able to extract reasonable
+ * values from any face polynomial. In order to make the use of
+ * FESystem simpler, FEValues objects will not fail using this finite
+ * element space, but all shape function values extracted will equal
+ * to zero.
+ *
+ * @ingroup fe
+ * @author Martin Kronbichler
+ * @date 2013
+ */
+template <int dim, int spacedim=dim>
+class FE_FaceP : public FE_PolyFace<PolynomialSpace<dim-1>, dim, spacedim>
+{
+public:
+ /**
+ * Constructor for complete basis of polynomials of degree <tt>p</tt>. The
+ * shape functions created using this constructor correspond to Legendre
+ * polynomials in each coordinate direction.
+ */
+ FE_FaceP(unsigned int p);
+
+ /**
+ * Clone method.
+ */
+ virtual FiniteElement<dim,spacedim> *clone() const;
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>FE_FaceP<dim>(degree)</tt> , with <tt>dim</tt> and
+ * <tt>degree</tt> replaced by appropriate values.
+ */
+ virtual std::string get_name () const;
+
+ /**
+ * Return the matrix interpolating from a face of of one element to the face
+ * of the neighboring element. The size of the matrix is then
+ * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+ * element only provides interpolation matrices for elements of the same
+ * type and FE_Nothing. For all other elements, an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+ */
+ virtual void
+ get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+ FullMatrix<double> &matrix) const;
+
+ /**
+ * Return the matrix interpolating from a face of of one element to the face
+ * of the neighboring element. The size of the matrix is then
+ * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+ * element only provides interpolation matrices for elements of the same
+ * type and FE_Nothing. For all other elements, an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+ */
+ virtual void
+ get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+ const unsigned int subface,
+ FullMatrix<double> &matrix) const;
+
+ /**
+ * Check for non-zero values on a face.
*
- * For the FE_Q class the result is always true (independent of the degree
- * of the element), as it implements the complete set of functions necessary
- * for hp capability.
+ * This function returns @p true, if the shape function @p shape_index has
+ * non-zero values on the face @p face_index.
+ *
+ * Implementation of the interface in FiniteElement
+ */
+ virtual bool has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const;
+
+ /**
+ * Return whether this element implements its hanging node constraints in
+ * the new way, which has to be used to make elements "hp compatible".
*/
virtual bool hp_constraints_are_implemented () const;
static std::vector<unsigned int> get_dpo_vector (const unsigned int deg);
};
+
DEAL_II_NAMESPACE_CLOSE
#endif
std::vector<ComponentMask>(FiniteElementData<dim>(
get_dpo_vector(degree), 1, degree).dofs_per_cell, std::vector<bool>(1,true)))
{
- // Reinit the vectors of
- // restriction and prolongation
- // matrices to the right sizes
+ // Reinit the vectors of restriction and prolongation matrices to the right
+ // sizes
this->reinit_restriction_and_prolongation_matrices();
// Fill prolongation matrices with embedding operators
if (dim == spacedim)
std::string
FE_DGP<dim,spacedim>::get_name () const
{
- // note that the
- // FETools::get_fe_from_name
- // function depends on the
- // particular format of the string
- // this function returns, so they
- // have to be kept in synch
+ // note that the FETools::get_fe_from_name function depends on the
+ // particular format of the string this function returns, so they have to be
+ // kept in synch
std::ostringstream namebuf;
namebuf << "FE_DGP<" << dim << ">(" << this->degree << ")";
get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
FullMatrix<double> &interpolation_matrix) const
{
- // this is only implemented, if the source
- // FE is also a DGP element. in that case,
- // both elements have no dofs on their
- // faces and the face interpolation matrix
- // is necessarily empty -- i.e. there isn't
- // much we need to do here.
+ // this is only implemented, if the source FE is also a DGP element. in that
+ // case, both elements have no dofs on their faces and the face
+ // interpolation matrix is necessarily empty -- i.e. there isn't much we
+ // need to do here.
typedef FiniteElement<dim,spacedim> FE;
typedef FE_DGP<dim,spacedim> FEDGP;
AssertThrow ((x_source_fe.get_name().find ("FE_DGP<") == 0)
const unsigned int ,
FullMatrix<double> &interpolation_matrix) const
{
- // this is only implemented, if the source
- // FE is also a DGP element. in that case,
- // both elements have no dofs on their
- // faces and the face interpolation matrix
- // is necessarily empty -- i.e. there isn't
- // much we need to do here.
+ // this is only implemented, if the source FE is also a DGP element. in that
+ // case, both elements have no dofs on their faces and the face
+ // interpolation matrix is necessarily empty -- i.e. there isn't much we
+ // need to do here.
typedef FiniteElement<dim,spacedim> FE;
typedef FE_DGP<dim,spacedim> FEDGP;
AssertThrow ((x_source_fe.get_name().find ("FE_DGP<") == 0)
FE_DGP<dim,spacedim>::
hp_vertex_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
{
- // there are no such constraints for DGP
- // elements at all
+ // there are no such constraints for DGP elements at all
if (dynamic_cast<const FE_DGP<dim,spacedim>*>(&fe_other) != 0)
return
std::vector<std::pair<unsigned int, unsigned int> > ();
FE_DGP<dim,spacedim>::
hp_line_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
{
- // there are no such constraints for DGP
- // elements at all
+ // there are no such constraints for DGP elements at all
if (dynamic_cast<const FE_DGP<dim,spacedim>*>(&fe_other) != 0)
return
std::vector<std::pair<unsigned int, unsigned int> > ();
FE_DGP<dim,spacedim>::
hp_quad_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
{
- // there are no such constraints for DGP
- // elements at all
+ // there are no such constraints for DGP elements at all
if (dynamic_cast<const FE_DGP<dim,spacedim>*>(&fe_other) != 0)
return
std::vector<std::pair<unsigned int, unsigned int> > ();
FiniteElementDomination::Domination
FE_DGP<dim,spacedim>::compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
{
- // check whether both are discontinuous
- // elements, see the description of
+ // check whether both are discontinuous elements, see the description of
// FiniteElementDomination::Domination
if (dynamic_cast<const FE_DGP<dim,spacedim>*>(&fe_other) != 0)
return FiniteElementDomination::no_requirements;
FE_DGP<dim,spacedim>::has_support_on_face (const unsigned int,
const unsigned int) const
{
- // all shape functions have support on all
- // faces
+ // all shape functions have support on all faces
return true;
}
#include <deal.II/fe/fe_face.h>
#include <deal.II/fe/fe_poly_face.templates.h>
#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_tools.h>
#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/householder.h>
#include <sstream>
DEAL_II_NAMESPACE_OPEN
AssertDimension (k, this->unit_face_support_points.size());
}
- // initialize unit support points
+ // initialize unit support points (this makes it possible to assign initial
+ // values to FE_FaceQ)
this->unit_support_points.resize(GeometryInfo<dim>::faces_per_cell*
this->unit_face_support_points.size());
const unsigned int n_face_dofs = this->unit_face_support_points.size();
}
+
template <int dim, int spacedim>
FiniteElement<dim,spacedim> *
FE_FaceQ<dim,spacedim>::clone() const
}
+
template <int dim, int spacedim>
std::string
FE_FaceQ<dim,spacedim>::get_name () const
{
- // note that the
- // FETools::get_fe_from_name
- // function depends on the
- // particular format of the string
- // this function returns, so they
- // have to be kept in synch
-
+ // note that the FETools::get_fe_from_name function depends on the
+ // particular format of the string this function returns, so they have to be
+ // kept in synch
std::ostringstream namebuf;
namebuf << "FE_FaceQ<" << dim << ">(" << this->degree << ")";
template <int dim, int spacedim>
void
FE_FaceQ<dim,spacedim>::
-get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source_fe,
FullMatrix<double> &interpolation_matrix) const
{
- // this function is similar to the respective method in FE_Q
-
- // this is only implemented, if the source FE is also a FE_FaceQ element
- AssertThrow ((dynamic_cast<const FE_FaceQ<dim,spacedim> *>(&x_source_fe) != 0),
- (typename FiniteElement<dim,spacedim>::
- ExcInterpolationNotImplemented()));
-
- Assert (interpolation_matrix.n() == this->dofs_per_face,
- ExcDimensionMismatch (interpolation_matrix.n(),
- this->dofs_per_face));
- Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
- ExcDimensionMismatch (interpolation_matrix.m(),
- x_source_fe.dofs_per_face));
-
- // ok, source is a FaceQ element, so we will be able to do the work
- const FE_FaceQ<dim,spacedim> &source_fe
- = dynamic_cast<const FE_FaceQ<dim,spacedim>&>(x_source_fe);
-
- // Make sure that the element for which the DoFs should be constrained is
- // the one with the higher polynomial degree. Actually the procedure will
- // work also if this assertion is not satisfied. But the matrices produced
- // in that case might lead to problems in the hp procedures, which use this
- // method.
- Assert (this->dofs_per_face <= source_fe.dofs_per_face,
- (typename FiniteElement<dim,spacedim>::
- ExcInterpolationNotImplemented ()));
-
- // generate a quadrature with the unit face support points.
- const Quadrature<dim-1> face_quadrature (source_fe.get_unit_face_support_points ());
-
- // Rule of thumb for FP accuracy, that can be expected for a given
- // polynomial degree. This value is used to cut off values close to zero.
- const double eps = 2e-13*(this->degree+1)*(dim-1);
-
- // compute the interpolation matrix by simply taking the value at the
- // support points.
- for (unsigned int i=0; i<source_fe.dofs_per_face; ++i)
- {
- const Point<dim-1> &p = face_quadrature.point (i);
-
- for (unsigned int j=0; j<this->dofs_per_face; ++j)
- {
- double matrix_entry = this->poly_space.compute_value (j, p);
-
- // Correct the interpolated value. I.e. if it is close to 1 or 0,
- // make it exactly 1 or 0. Unfortunately, this is required to avoid
- // problems with higher order elements.
- if (std::fabs (matrix_entry - 1.0) < eps)
- matrix_entry = 1.0;
- if (std::fabs (matrix_entry) < eps)
- matrix_entry = 0.0;
-
- interpolation_matrix(i,j) = matrix_entry;
- }
- }
-
- // make sure that the row sum of each of the matrices is 1 at this
- // point. this must be so since the shape functions sum up to 1
- for (unsigned int j=0; j<source_fe.dofs_per_face; ++j)
- {
- double sum = 0.;
-
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- sum += interpolation_matrix(j,i);
-
- Assert (std::fabs(sum-1) < eps, ExcInternalError());
- }
+ get_subface_interpolation_matrix (source_fe, numbers::invalid_unsigned_int,
+ interpolation_matrix);
}
for (unsigned int i=0; i<source_fe->dofs_per_face; ++i)
{
const Point<dim-1> p =
+ subface == numbers::invalid_unsigned_int
+ ?
+ face_quadrature.point(i)
+ :
GeometryInfo<dim-1>::child_to_cell_coordinates (face_quadrature.point(i),
subface);
return FiniteElementDomination::neither_element_dominates;
}
+
+
+// --------------------------------------- FE_FaceP --------------------------
+
+template <int dim, int spacedim>
+FE_FaceP<dim,spacedim>::FE_FaceP (const unsigned int degree)
+ :
+ FE_PolyFace<PolynomialSpace<dim-1>, dim, spacedim>
+ (PolynomialSpace<dim-1>(Polynomials::Legendre::generate_complete_basis(degree)),
+ FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
+ std::vector<bool>(1,true))
+{}
+
+
+template <int dim, int spacedim>
+FiniteElement<dim,spacedim> *
+FE_FaceP<dim,spacedim>::clone() const
+{
+ return new FE_FaceP<dim,spacedim>(this->degree);
+}
+
+
+template <int dim, int spacedim>
+std::string
+FE_FaceP<dim,spacedim>::get_name () const
+{
+ // note that the FETools::get_fe_from_name function depends on the
+ // particular format of the string this function returns, so they have to be
+ // kept in synch
+ std::ostringstream namebuf;
+ namebuf << "FE_FaceP<" << dim << ">(" << this->degree << ")";
+
+ return namebuf.str();
+}
+
+
+
+template <int dim, int spacedim>
+bool
+FE_FaceP<dim,spacedim>::has_support_on_face (
+ const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ return (face_index == (shape_index/this->dofs_per_face));
+}
+
+
+
+template <int dim, int spacedim>
+std::vector<unsigned int>
+FE_FaceP<dim,spacedim>::get_dpo_vector (const unsigned int deg)
+{
+ std::vector<unsigned int> dpo(dim+1, 0U);
+ dpo[dim-1] = deg+1;
+ for (unsigned int i=1; i<dim-1; ++i)
+ {
+ dpo[dim-1] *= deg+1+i;
+ dpo[dim-1] /= i+1;
+ }
+ return dpo;
+}
+
+
+
+
+template <int dim, int spacedim>
+bool
+FE_FaceP<dim,spacedim>::hp_constraints_are_implemented () const
+{
+ return true;
+}
+
+
+
+template <int dim, int spacedim>
+FiniteElementDomination::Domination
+FE_FaceP<dim,spacedim>::
+compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
+{
+ if (const FE_FaceP<dim,spacedim> *fe_q_other
+ = dynamic_cast<const FE_FaceP<dim,spacedim>*>(&fe_other))
+ {
+ if (this->degree < fe_q_other->degree)
+ return FiniteElementDomination::this_element_dominates;
+ else if (this->degree == fe_q_other->degree)
+ return FiniteElementDomination::either_element_can_dominate;
+ else
+ return FiniteElementDomination::other_element_dominates;
+ }
+ else if (dynamic_cast<const FE_Nothing<dim>*>(&fe_other) != 0)
+ {
+ // the FE_Nothing has no degrees of freedom and it is typically used in
+ // a context where we don't require any continuity along the interface
+ return FiniteElementDomination::no_requirements;
+ }
+
+ Assert (false, ExcNotImplemented());
+ return FiniteElementDomination::neither_element_dominates;
+}
+
+
+
+
+template <int dim, int spacedim>
+void
+FE_FaceP<dim,spacedim>::
+get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source_fe,
+ FullMatrix<double> &interpolation_matrix) const
+{
+ get_subface_interpolation_matrix (source_fe, numbers::invalid_unsigned_int,
+ interpolation_matrix);
+}
+
+
+
+template <int dim, int spacedim>
+void
+FE_FaceP<dim,spacedim>::
+get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+ const unsigned int subface,
+ FullMatrix<double> &interpolation_matrix) const
+{
+ // this function is similar to the respective method in FE_Q
+
+ Assert (interpolation_matrix.n() == this->dofs_per_face,
+ ExcDimensionMismatch (interpolation_matrix.n(),
+ this->dofs_per_face));
+ Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ x_source_fe.dofs_per_face));
+
+ // see if source is a FaceP element
+ if (const FE_FaceP<dim,spacedim> *source_fe
+ = dynamic_cast<const FE_FaceP<dim,spacedim> *>(&x_source_fe))
+ {
+ // Make sure that the element for which the DoFs should be constrained
+ // is the one with the higher polynomial degree. Actually the procedure
+ // will work also if this assertion is not satisfied. But the matrices
+ // produced in that case might lead to problems in the hp procedures,
+ // which use this method.
+ Assert (this->dofs_per_face <= source_fe->dofs_per_face,
+ (typename FiniteElement<dim,spacedim>::
+ ExcInterpolationNotImplemented ()));
+
+ // do this as in FETools by solving a least squares problem where we
+ // force the source FE polynomial to be equal the given FE on all
+ // quadrature points
+ const QGauss<dim-1> face_quadrature (source_fe->degree+1);
+
+ // Rule of thumb for FP accuracy, that can be expected for a given
+ // polynomial degree. This value is used to cut off values close to
+ // zero.
+ const double eps = 2e-13*(this->degree+1)*(dim-1);
+
+ FullMatrix<double> mass (face_quadrature.size(), this->dofs_per_face);
+
+ for (unsigned int k = 0; k < face_quadrature.size(); ++k)
+ {
+ const Point<dim-1> p =
+ GeometryInfo<dim-1>::child_to_cell_coordinates (face_quadrature.point(k),
+ subface);
+
+ for (unsigned int j = 0; j < this->dofs_per_face; ++j)
+ mass (k , j) = this->poly_space.compute_value(j, p);
+ }
+
+ Householder<double> H(mass);
+ Vector<double> v_in(face_quadrature.size());
+ Vector<double> v_out(this->dofs_per_face);
+
+
+ // compute the interpolation matrix by evaluating on the fine side and
+ // then solving the least squares problem
+ for (unsigned int i=0; i<source_fe->dofs_per_face; ++i)
+ {
+ for (unsigned int k = 0; k < face_quadrature.size(); ++k)
+ {
+ const Point<dim-1> p =
+ GeometryInfo<dim-1>::child_to_cell_coordinates (face_quadrature.point(k),
+ subface);
+ v_in(k) = source_fe->poly_space.compute_value(i, p);
+ }
+ const double result = H.least_squares(v_out, v_in);
+ Assert(result < 1e-12, FETools::ExcLeastSquaresError (result));
+
+ for (unsigned int j = 0; j < this->dofs_per_face; ++j)
+ {
+ double matrix_entry = v_out(j);
+
+ // Correct the interpolated value. I.e. if it is close to 1 or 0,
+ // make it exactly 1 or 0. Unfortunately, this is required to avoid
+ // problems with higher order elements.
+ if (std::fabs (matrix_entry - 1.0) < eps)
+ matrix_entry = 1.0;
+ if (std::fabs (matrix_entry) < eps)
+ matrix_entry = 0.0;
+
+ interpolation_matrix(i,j) = matrix_entry;
+ }
+ }
+
+ // make sure that the row sum of each of the matrices is 1 at this
+ // point. this must be so since the shape functions sum up to 1
+ for (unsigned int j=0; j<source_fe->dofs_per_face; ++j)
+ {
+ double sum = 0.;
+
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ sum += interpolation_matrix(j,i);
+
+ Assert (std::fabs(sum-1) < eps, ExcInternalError());
+ }
+ }
+ else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != 0)
+ {
+ // nothing to do here, the FE_Nothing has no degrees of freedom anyway
+ }
+ else
+ AssertThrow (false,(typename FiniteElement<dim,spacedim>::
+ ExcInterpolationNotImplemented()));
+}
+
+
// explicit instantiations
#include "fe_face.inst"
template <class POLY, int dim, int spacedim>
void
FE_Q_Base<POLY,dim,spacedim>::
-get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source_fe,
FullMatrix<double> &interpolation_matrix) const
{
Assert (dim > 1, ExcImpossibleInDim(1));
-
- // this is only implemented, if the source FE is also a Q element
- AssertThrow ((dynamic_cast<const FE_Q_Base<POLY,dim,spacedim> *>(&x_source_fe) != 0),
- (typename FiniteElement<dim,spacedim>::
- ExcInterpolationNotImplemented()));
-
- Assert (interpolation_matrix.n() == this->dofs_per_face,
- ExcDimensionMismatch (interpolation_matrix.n(),
- this->dofs_per_face));
- Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
- ExcDimensionMismatch (interpolation_matrix.m(),
- x_source_fe.dofs_per_face));
-
- // ok, source is a Q element, so we will be able to do the work
- const FE_Q_Base<POLY,dim,spacedim> &source_fe
- = dynamic_cast<const FE_Q_Base<POLY,dim,spacedim>&>(x_source_fe);
-
- // Make sure that the element for which the DoFs should be constrained is
- // the one with the higher polynomial degree. Actually the procedure will
- // work also if this assertion is not satisfied. But the matrices produced
- // in that case might lead to problems in the hp procedures, which use this
- // method.
- Assert (this->dofs_per_face <= source_fe.dofs_per_face,
- (typename FiniteElement<dim,spacedim>::
- ExcInterpolationNotImplemented ()));
-
- // generate a quadrature with the unit support points. This is later based
- // as a basis for the QProjector, which returns the support points on the
- // face.
- Quadrature<dim-1> quad_face_support (source_fe.get_unit_face_support_points ());
-
- // Rule of thumb for FP accuracy, that can be expected for a given
- // polynomial degree. This value is used to cut off values close to zero.
- const double eps = 2e-13*this->degree*(dim-1);
-
- // compute the interpolation matrix by simply taking the value at the
- // support points.
-//TODO: Verify that all faces are the same with respect to
-// these support points. Furthermore, check if something has to
-// be done for the face orientation flag in 3D.
- const Quadrature<dim> face_quadrature
- = QProjector<dim>::project_to_face (quad_face_support, 0);
- for (unsigned int i=0; i<source_fe.dofs_per_face; ++i)
- {
- const Point<dim> &p = face_quadrature.point (i);
-
- for (unsigned int j=0; j<this->dofs_per_face; ++j)
- {
- double matrix_entry = this->poly_space.compute_value (this->face_to_cell_index(j, 0), p);
-
- // Correct the interpolated value. I.e. if it is close to 1 or 0,
- // make it exactly 1 or 0. Unfortunately, this is required to avoid
- // problems with higher order elements.
- if (std::fabs (matrix_entry - 1.0) < eps)
- matrix_entry = 1.0;
- if (std::fabs (matrix_entry) < eps)
- matrix_entry = 0.0;
-
- interpolation_matrix(i,j) = matrix_entry;
- }
- }
-
- // make sure that the row sum of each of the matrices is 1 at this
- // point. this must be so since the shape functions sum up to 1
- for (unsigned int j=0; j<source_fe.dofs_per_face; ++j)
- {
- double sum = 0.;
-
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- sum += interpolation_matrix(j,i);
-
- Assert (std::fabs(sum-1) < eps, ExcInternalError());
- }
+ get_subface_interpolation_matrix (source_fe, numbers::invalid_unsigned_int,
+ interpolation_matrix);
}
// these support points. Furthermore, check if something has to
// be done for the face orientation flag in 3D.
const Quadrature<dim> subface_quadrature
- = QProjector<dim>::project_to_subface (quad_face_support, 0, subface);
+ = subface == numbers::invalid_unsigned_int
+ ?
+ QProjector<dim>::project_to_face (quad_face_support, 0)
+ :
+ QProjector<dim>::project_to_subface (quad_face_support, 0, subface);
for (unsigned int i=0; i<source_fe->dofs_per_face; ++i)
{
const Point<dim> &p = subface_quadrature.point (i);
projects / dealii-svn.git / commitdiff