From 8960a464202f69ea24befb48133f7d8f074e5a0b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 24 Feb 2011 02:43:43 +0000 Subject: [PATCH] Handle constraints between FE_Q and FE_Nothing. git-svn-id: https://svn.dealii.org/trunk@23444 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/fe/fe_q.cc | 167 +++++++++++++++++++------------------- 1 file changed, 85 insertions(+), 82 deletions(-) diff --git a/deal.II/source/fe/fe_q.cc b/deal.II/source/fe/fe_q.cc index bd1751d11b..f08b5fc4dd 100644 --- a/deal.II/source/fe/fe_q.cc +++ b/deal.II/source/fe/fe_q.cc @@ -973,103 +973,106 @@ 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 - typedef FE_Q FEQ; - typedef FiniteElement FEL; - AssertThrow ((x_source_fe.get_name().find ("FE_Q<") == 0) - || - (dynamic_cast(&x_source_fe) != 0), - typename FEL:: - 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 &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 FEL:: - ExcInterpolationNotImplemented ()); - - // generate a point on this - // cell and evaluate the - // shape functions there - const 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. - double eps = 2e-13*this->degree*(dim-1); - - // compute the interpolation - // matrix by simply taking the - // value at the support points. + // see if source is a Q element + if (const FE_Q *source_fe + = dynamic_cast *>(&x_source_fe)) + { + // have this test in here since + // a table of size 2x0 reports + // its size as 0x0 + Assert (interpolation_matrix.n() == this->dofs_per_face, + ExcDimensionMismatch (interpolation_matrix.n(), + this->dofs_per_face)); + + // 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 point on this + // cell and evaluate the + // shape functions there + const 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. + 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 subface_quadrature - = QProjector::project_to_subface (quad_face_support, 0, subface); - for (unsigned int i=0; i &p = subface_quadrature.point (i); + const Quadrature subface_quadrature + = 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); - for (unsigned int j=0; jdofs_per_face; ++j) + for (unsigned int j=0; jdofs_per_face; ++j) + { + double matrix_entry = this->shape_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; ++j) { - double matrix_entry = this->shape_value (this->face_to_cell_index(j, 0), p); + double sum = 0.; - // 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; + for (unsigned int i=0; idofs_per_face; ++i) + sum += interpolation_matrix(j,i); - interpolation_matrix(i,j) = matrix_entry; + Assert (std::fabs(sum-1) < 2e-13*this->degree*this->degree*dim, + ExcInternalError()); } } - - // 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 *>(&x_source_fe) != 0) { - double sum = 0.; - - for (unsigned int i=0; idofs_per_face; ++i) - sum += interpolation_matrix(j,i); - - Assert (std::fabs(sum-1) < 2e-13*this->degree*this->degree*dim, - ExcInternalError()); + // nothing to do here, the + // FE_Nothing has no degrees of + // freedom anyway } + else + AssertThrow (false, + (typename FiniteElement:: + ExcInterpolationNotImplemented())); } -- 2.39.5