From f5f1dce43774e5174fc2b25a885bc1dda9482122 Mon Sep 17 00:00:00 2001 From: kayser-herold Date: Thu, 27 Jul 2006 16:32:30 +0000 Subject: [PATCH] The FiniteElement class was extended by two methods, which create interpolation matrices between element faces. These functions will be needed for the hp FE method. git-svn-id: https://svn.dealii.org/trunk@13450 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 58 +++++++- deal.II/deal.II/include/fe/fe_base.h | 4 +- deal.II/deal.II/include/fe/fe_q.h | 58 +++++++- deal.II/deal.II/source/dofs/dof_tools.cc | 4 + deal.II/deal.II/source/fe/fe.cc | 33 +++++ deal.II/deal.II/source/fe/fe_q.cc | 161 +++++++++++++++++++++++ 6 files changed, 315 insertions(+), 3 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 7dd6811cdf..37aaef223e 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -896,7 +896,7 @@ class FiniteElement : public Subscriptor, */ virtual void get_interpolation_matrix (const FiniteElement &source, - FullMatrix &matrix) const; + FullMatrix &matrix) const; //@} /** @@ -904,6 +904,62 @@ class FiniteElement : public Subscriptor, * @{ */ + + /** + * 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 @p dofs_per_face times + * source.dofs_per_face. + * + * Derived elements will have to + * implement this function. They + * may only provide interpolation + * matrices for certain source + * finite elements, for example + * those from the same family. If + * they don't implement + * interpolation from a given + * element, then they must throw + * an exception of type + * FiniteElement::ExcInterpolationNotImplemented. + */ + 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 subface of + * the neighboring element. + * The size of the matrix is + * then @p dofs_per_face times + * source.dofs_per_face. + * + * Derived elements will have to + * implement this function. They + * may only provide interpolation + * matrices for certain source + * finite elements, for example + * those from the same family. If + * they don't implement + * interpolation from a given + * element, then they must throw + * an exception of type + * FiniteElement::ExcInterpolationNotImplemented. + */ + virtual void + get_subface_interpolation_matrix (const FiniteElement &source, + const unsigned int subface, + FullMatrix &matrix) const; + //@} + + /** * If, on a vertex, several * finite elements are active, diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index e48f7c0d2c..929528a24d 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -560,7 +560,9 @@ face_to_equivalent_cell_index (const unsigned int index) const } - +//TODO: This method produces results, that differ from those which +// are returned bz "face_to_cell_index (index, 0)". This has to be +// verified. template <> inline unsigned int diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index 9ebe9dab1d..a3e6915032 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -271,8 +271,64 @@ class FE_Q : public FE_Poly,dim> */ virtual void get_interpolation_matrix (const FiniteElement &source, - FullMatrix &matrix) const; + 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 @p dofs_per_face times + * source.dofs_per_face. + * + * Derived elements will have to + * implement this function. They + * may only provide interpolation + * matrices for certain source + * finite elements, for example + * those from the same family. If + * they don't implement + * interpolation from a given + * element, then they must throw + * an exception of type + * FiniteElement::ExcInterpolationNotImplemented. + */ + 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 @p dofs_per_face times + * source.dofs_per_face. + * + * Derived elements will have to + * implement this function. They + * may only provide interpolation + * matrices for certain source + * finite elements, for example + * those from the same family. If + * they don't implement + * interpolation from a given + * element, then they must throw + * an exception of type + * FiniteElement::ExcInterpolationNotImplemented. + */ + virtual void + get_subface_interpolation_matrix (const FiniteElement &source, + const unsigned int subface, + FullMatrix &matrix) const; + //@} + + /** * Check for non-zero values on a face. * diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 7f2f251e39..59534efa5c 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -1704,6 +1704,8 @@ namespace internal = this_face->child(child)->dof_index(dof, fe_index); Assert (next_index == dofs_on_children.size(), ExcInternalError()); + +// TODO[OKH]: FE::get_subface_interpolation_matrix () // for each row in the constraint // matrix for this line: @@ -1729,6 +1731,8 @@ namespace internal Assert (cell->face(face) ->fe_index_is_active(cell->active_fe_index()) == true, ExcInternalError()); + +// TODO[OKH]: FE::get_face_interpolation_matrix () } } #endif diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 13b1c8ea87..dcc9cc824b 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -442,6 +442,39 @@ get_interpolation_matrix (const FiniteElement &, +template +void +FiniteElement:: +get_face_interpolation_matrix (const FiniteElement &, + FullMatrix &) const +{ + // by default, no interpolation + // implemented. so throw exception, + // as documentation says + AssertThrow (false, + typename FiniteElement:: + ExcInterpolationNotImplemented()); +} + + + +template +void +FiniteElement:: +get_subface_interpolation_matrix (const FiniteElement &, + const unsigned int, + FullMatrix &) const +{ + // by default, no interpolation + // implemented. so throw exception, + // as documentation says + AssertThrow (false, + typename FiniteElement:: + ExcInterpolationNotImplemented()); +} + + + template std::vector > FiniteElement:: diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 3b2fbf8f0a..d25befe7a4 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -324,6 +324,167 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, } +template +void +FE_Q:: +get_face_interpolation_matrix (const FiniteElement &x_source_fe, + FullMatrix &interpolation_matrix) const +{ + // this is only implemented, if the + // source FE is also a + // Q element + AssertThrow ((x_source_fe.get_name().find ("FE_Q<") == 0) + || + (dynamic_cast*>(&x_source_fe) != 0), + typename FiniteElement:: + ExcInterpolationNotImplemented()); + + // ok, source is a Q element, so + // we will be able to do the work + const FE_Q &source_fe + = dynamic_cast&>(x_source_fe); + + Assert (interpolation_matrix.m() == this->dofs_per_face, + ExcDimensionMismatch (interpolation_matrix.m(), + this->dofs_per_face)); + Assert (interpolation_matrix.n() == source_fe.dofs_per_face, + ExcDimensionMismatch (interpolation_matrix.m(), + source_fe.dofs_per_face)); + + const std::vector &index_map= + this->poly_space.get_numbering(); + + + // generate a point on this + // cell and evaluate the + // shape functions there + Quadrature quad_face_support (source_fe.get_unit_face_support_points ()); + + + // compute the interpolation + // matrix by simply taking the + // value at the support points. + for (unsigned int i=0; i p = QProjector::project_to_face (quad_face_support, 0).point (i); + + for (unsigned int j=0; jdofs_per_face; ++j) + interpolation_matrix(j,i) = this->shape_value (this->face_to_cell_index(j, 0), p); + } + + // cut off very small values + for (unsigned int i=0; idofs_per_face; ++i) + for (unsigned int j=0; jdofs_per_face; ++i) + sum += interpolation_matrix(i,j); + + Assert (std::fabs(sum-1) < 2e-14*this->degree*(dim-1), + ExcInternalError()); + } +} + + +template +void +FE_Q:: +get_subface_interpolation_matrix (const FiniteElement &x_source_fe, + const unsigned int subface, + FullMatrix &interpolation_matrix) const +{ + // this is only implemented, if the + // source FE is also a + // Q element + AssertThrow ((x_source_fe.get_name().find ("FE_Q<") == 0) + || + (dynamic_cast*>(&x_source_fe) != 0), + typename FiniteElement:: + ExcInterpolationNotImplemented()); + + // ok, source is a Q element, so + // we will be able to do the work + const FE_Q &source_fe + = dynamic_cast&>(x_source_fe); + + Assert (interpolation_matrix.m() == this->dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + this->dofs_per_cell)); + Assert (interpolation_matrix.n() == source_fe.dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + source_fe.dofs_per_cell)); + + const std::vector &index_map= + this->poly_space.get_numbering(); + + // compute the interpolation + // matrices in much the same way as + // we do for the embedding matrices + // from mother to child. + FullMatrix cell_interpolation (this->dofs_per_cell, + this->dofs_per_cell); + FullMatrix source_interpolation (this->dofs_per_cell, + source_fe.dofs_per_cell); + FullMatrix tmp (this->dofs_per_cell, + source_fe.dofs_per_cell); + for (unsigned int j=0; jdofs_per_cell; ++j) + { + // generate a point on this + // cell and evaluate the + // shape functions there + const Point + p = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell, + FE_Q_Helper::int2type()); + for (unsigned int i=0; idofs_per_cell; ++i) + cell_interpolation(j,i) = this->poly_space.compute_value (i, p); + + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++i) + { + double sum = 0.; + for (unsigned int j=0; jdegree*dim, + ExcInternalError()); + } +} + template std::vector > -- 2.39.5