From 01974dca0b1f6887ad40848dd7b6d662a19a3f64 Mon Sep 17 00:00:00 2001 From: kayser-herold Date: Thu, 27 Jul 2006 19:38:45 +0000 Subject: [PATCH] Fixed a bug in the face/subface interpolation procedures. These do not have any useful meaning in 1D. Now there is an empty implementation for the 1D case. git-svn-id: https://svn.dealii.org/trunk@13458 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_q.cc | 143 +++++++++++++++++++----------- 1 file changed, 93 insertions(+), 50 deletions(-) diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 62e1dedaa5..c3a75618f4 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -279,9 +279,9 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, 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 + // 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()); @@ -323,6 +323,34 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, } +#if deal_II_dimension == 1 + +template <> +void +FE_Q<1>:: +get_face_interpolation_matrix (const FiniteElement<1> &/*x_source_fe*/, + FullMatrix &/*interpolation_matrix*/) const +{ + Assert (false, + FiniteElement<1>:: + ExcInterpolationNotImplemented ()); +} + + +template <> +void +FE_Q<1>:: +get_subface_interpolation_matrix (const FiniteElement<1> &/*x_source_fe*/, + const unsigned int /*subface*/, + FullMatrix &/*interpolation_matrix*/) const +{ +} + +#endif + + +#if deal_II_dimension > 1 + template void FE_Q:: @@ -343,6 +371,21 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, const FE_Q &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 ()); + Assert (interpolation_matrix.m() == this->dofs_per_face, ExcDimensionMismatch (interpolation_matrix.m(), this->dofs_per_face)); @@ -350,13 +393,12 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, 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 + // 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 ()); @@ -419,71 +461,72 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, const FE_Q &source_fe = dynamic_cast&>(x_source_fe); - Assert (interpolation_matrix.m() == this->dofs_per_cell, + // 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 ()); + + Assert (interpolation_matrix.m() == this->dofs_per_face, ExcDimensionMismatch (interpolation_matrix.m(), - this->dofs_per_cell)); - Assert (interpolation_matrix.n() == source_fe.dofs_per_cell, + this->dofs_per_face)); + Assert (interpolation_matrix.n() == source_fe.dofs_per_face, ExcDimensionMismatch (interpolation_matrix.m(), - source_fe.dofs_per_cell)); - - const std::vector &index_map= - this->poly_space.get_numbering(); + source_fe.dofs_per_face)); - // 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); + Quadrature quad_face_support (source_fe.get_unit_face_support_points ()); - for (unsigned int i=0; i p = QProjector::project_to_subface (quad_face_support, 0, subface).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_cell; ++i) - for (unsigned int j=0; jdofs_per_face; ++i) + for (unsigned int j=0; jdofs_per_cell; ++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, + Assert (std::fabs(sum-1) < 2e-14*this->degree*(dim-1), ExcInternalError()); } } +#endif template std::vector > -- 2.39.5