From 615749dee5753d4cb67e49145ba9ca8b75263c6a Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 17 Sep 2013 11:26:58 +0000 Subject: [PATCH] First step towards implementing FE_FaceP in analogy to FE_FaceQ. Cleanup some code. git-svn-id: https://svn.dealii.org/trunk@30760 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/fe/fe_face.h | 107 ++++++++- deal.II/source/fe/fe_dgp.cc | 49 ++--- deal.II/source/fe/fe_face.cc | 315 ++++++++++++++++++++------- deal.II/source/fe/fe_q_base.cc | 83 +------ 4 files changed, 369 insertions(+), 185 deletions(-) diff --git a/deal.II/include/deal.II/fe/fe_face.h b/deal.II/include/deal.II/fe/fe_face.h index 3047f7a06a..e6464118ba 100644 --- a/deal.II/include/deal.II/fe/fe_face.h +++ b/deal.II/include/deal.II/fe/fe_face.h @@ -19,6 +19,7 @@ #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -105,10 +106,109 @@ public: /** * 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 &fe_other) const; + +private: + /** + * Return vector with dofs per vertex, line, quad, hex. + */ + static std::vector 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 +class FE_FaceP : public FE_PolyFace, dim, spacedim> +{ +public: + /** + * Constructor for complete basis of polynomials of degree p. The + * shape functions created using this constructor correspond to Legendre + * polynomials in each coordinate direction. + */ + FE_FaceP(unsigned int p); + + /** + * Clone method. + */ + virtual FiniteElement *clone() const; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns FE_FaceP(degree) , with dim and + * degree 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 + * source.dofs_per_face times this->dofs_per_face. This + * element only provides interpolation matrices for elements of the same + * type and FE_Nothing. For all other elements, an exception of type + * FiniteElement::ExcInterpolationNotImplemented is thrown. + */ + virtual void + get_face_interpolation_matrix (const FiniteElement &source, + FullMatrix &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 + * source.dofs_per_face times this->dofs_per_face. This + * element only provides interpolation matrices for elements of the same + * type and FE_Nothing. For all other elements, an exception of type + * FiniteElement::ExcInterpolationNotImplemented is thrown. + */ + virtual void + get_subface_interpolation_matrix (const FiniteElement &source, + const unsigned int subface, + FullMatrix &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; @@ -131,6 +231,7 @@ private: static std::vector get_dpo_vector (const unsigned int deg); }; + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/source/fe/fe_dgp.cc b/deal.II/source/fe/fe_dgp.cc index bb3f9531af..1c4d07bc3e 100644 --- a/deal.II/source/fe/fe_dgp.cc +++ b/deal.II/source/fe/fe_dgp.cc @@ -32,9 +32,8 @@ FE_DGP::FE_DGP (const unsigned int degree) std::vector(FiniteElementData( get_dpo_vector(degree), 1, degree).dofs_per_cell, std::vector(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) @@ -50,12 +49,9 @@ template std::string FE_DGP::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 << ")"; @@ -101,12 +97,10 @@ FE_DGP:: get_face_interpolation_matrix (const FiniteElement &x_source_fe, FullMatrix &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 FE; typedef FE_DGP FEDGP; AssertThrow ((x_source_fe.get_name().find ("FE_DGP<") == 0) @@ -132,12 +126,10 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe const unsigned int , FullMatrix &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 FE; typedef FE_DGP FEDGP; AssertThrow ((x_source_fe.get_name().find ("FE_DGP<") == 0) @@ -170,8 +162,7 @@ std::vector > FE_DGP:: hp_vertex_dof_identities (const FiniteElement &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*>(&fe_other) != 0) return std::vector > (); @@ -189,8 +180,7 @@ std::vector > FE_DGP:: hp_line_dof_identities (const FiniteElement &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*>(&fe_other) != 0) return std::vector > (); @@ -208,8 +198,7 @@ std::vector > FE_DGP:: hp_quad_dof_identities (const FiniteElement &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*>(&fe_other) != 0) return std::vector > (); @@ -226,8 +215,7 @@ template FiniteElementDomination::Domination FE_DGP::compare_for_face_domination (const FiniteElement &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*>(&fe_other) != 0) return FiniteElementDomination::no_requirements; @@ -243,8 +231,7 @@ bool FE_DGP::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; } diff --git a/deal.II/source/fe/fe_face.cc b/deal.II/source/fe/fe_face.cc index bbf8124ada..0478dfbad8 100644 --- a/deal.II/source/fe/fe_face.cc +++ b/deal.II/source/fe/fe_face.cc @@ -18,7 +18,9 @@ #include #include #include +#include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -73,7 +75,8 @@ FE_FaceQ::FE_FaceQ (const unsigned int degree) 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::faces_per_cell* this->unit_face_support_points.size()); const unsigned int n_face_dofs = this->unit_face_support_points.size(); @@ -94,6 +97,7 @@ FE_FaceQ::FE_FaceQ (const unsigned int degree) } + template FiniteElement * FE_FaceQ::clone() const @@ -102,17 +106,14 @@ FE_FaceQ::clone() const } + template std::string FE_FaceQ::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 << ")"; @@ -124,76 +125,11 @@ FE_FaceQ::get_name () const template void FE_FaceQ:: -get_face_interpolation_matrix (const FiniteElement &x_source_fe, +get_face_interpolation_matrix (const FiniteElement &source_fe, FullMatrix &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 *>(&x_source_fe) != 0), - (typename FiniteElement:: - 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 &source_fe - = dynamic_cast&>(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:: - ExcInterpolationNotImplemented ())); - - // generate a quadrature with the unit face support points. - const Quadrature 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 &p = face_quadrature.point (i); - - for (unsigned int j=0; jdofs_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; jdofs_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); } @@ -241,6 +177,10 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe for (unsigned int i=0; idofs_per_face; ++i) { const Point p = + subface == numbers::invalid_unsigned_int + ? + face_quadrature.point(i) + : GeometryInfo::child_to_cell_coordinates (face_quadrature.point(i), subface); @@ -342,6 +282,229 @@ compare_for_face_domination (const FiniteElement &fe_other) const return FiniteElementDomination::neither_element_dominates; } + + +// --------------------------------------- FE_FaceP -------------------------- + +template +FE_FaceP::FE_FaceP (const unsigned int degree) + : + FE_PolyFace, dim, spacedim> + (PolynomialSpace(Polynomials::Legendre::generate_complete_basis(degree)), + FiniteElementData(get_dpo_vector(degree), 1, degree, FiniteElementData::L2), + std::vector(1,true)) +{} + + +template +FiniteElement * +FE_FaceP::clone() const +{ + return new FE_FaceP(this->degree); +} + + +template +std::string +FE_FaceP::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 +bool +FE_FaceP::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 +std::vector +FE_FaceP::get_dpo_vector (const unsigned int deg) +{ + std::vector dpo(dim+1, 0U); + dpo[dim-1] = deg+1; + for (unsigned int i=1; i +bool +FE_FaceP::hp_constraints_are_implemented () const +{ + return true; +} + + + +template +FiniteElementDomination::Domination +FE_FaceP:: +compare_for_face_domination (const FiniteElement &fe_other) const +{ + if (const FE_FaceP *fe_q_other + = dynamic_cast*>(&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*>(&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 +void +FE_FaceP:: +get_face_interpolation_matrix (const FiniteElement &source_fe, + FullMatrix &interpolation_matrix) const +{ + get_subface_interpolation_matrix (source_fe, numbers::invalid_unsigned_int, + interpolation_matrix); +} + + + +template +void +FE_FaceP:: +get_subface_interpolation_matrix (const FiniteElement &x_source_fe, + const unsigned int subface, + FullMatrix &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 *source_fe + = dynamic_cast *>(&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:: + 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 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 mass (face_quadrature.size(), this->dofs_per_face); + + for (unsigned int k = 0; k < face_quadrature.size(); ++k) + { + const Point p = + GeometryInfo::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 H(mass); + Vector v_in(face_quadrature.size()); + Vector 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; idofs_per_face; ++i) + { + for (unsigned int k = 0; k < face_quadrature.size(); ++k) + { + const Point p = + GeometryInfo::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; jdofs_per_face; ++j) + { + double sum = 0.; + + for (unsigned int i=0; idofs_per_face; ++i) + sum += interpolation_matrix(j,i); + + Assert (std::fabs(sum-1) < eps, ExcInternalError()); + } + } + else if (dynamic_cast *>(&x_source_fe) != 0) + { + // nothing to do here, the FE_Nothing has no degrees of freedom anyway + } + else + AssertThrow (false,(typename FiniteElement:: + ExcInterpolationNotImplemented())); +} + + // explicit instantiations #include "fe_face.inst" diff --git a/deal.II/source/fe/fe_q_base.cc b/deal.II/source/fe/fe_q_base.cc index 6a73cb86ff..40096c09cf 100644 --- a/deal.II/source/fe/fe_q_base.cc +++ b/deal.II/source/fe/fe_q_base.cc @@ -550,83 +550,12 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, template void FE_Q_Base:: -get_face_interpolation_matrix (const FiniteElement &x_source_fe, +get_face_interpolation_matrix (const FiniteElement &source_fe, FullMatrix &interpolation_matrix) const { Assert (dim > 1, ExcImpossibleInDim(1)); - - // this is only implemented, if the source FE is also a Q element - AssertThrow ((dynamic_cast *>(&x_source_fe) != 0), - (typename FiniteElement:: - 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 &source_fe - = dynamic_cast&>(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:: - 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 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 face_quadrature - = QProjector::project_to_face (quad_face_support, 0); - for (unsigned int i=0; i &p = face_quadrature.point (i); - - for (unsigned int j=0; jdofs_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; jdofs_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); } @@ -676,7 +605,11 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe // these support points. Furthermore, check if something has to // be done for the face orientation flag in 3D. const Quadrature subface_quadrature - = QProjector::project_to_subface (quad_face_support, 0, subface); + = subface == numbers::invalid_unsigned_int + ? + QProjector::project_to_face (quad_face_support, 0) + : + QProjector::project_to_subface (quad_face_support, 0, subface); for (unsigned int i=0; idofs_per_face; ++i) { const Point &p = subface_quadrature.point (i); -- 2.39.5