From fcab682fae83ddb18843457342e14d71c2288292 Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Mon, 1 Feb 2021 13:02:43 -0700 Subject: [PATCH] simplex: adjust get_{sub}face_interpolation_matrix() for hybrid meshes. --- source/fe/fe_q_base.cc | 25 +++++++++++++------------ source/simplex/fe_lib.cc | 38 ++++++++++++++++++++------------------ 2 files changed, 33 insertions(+), 30 deletions(-) diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index 439d3ca81f..d64cc60b12 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -621,18 +621,19 @@ FE_Q_Base::get_face_interpolation_matrix( template void FE_Q_Base::get_subface_interpolation_matrix( - const FiniteElement &x_source_fe, + const FiniteElement &source_fe, const unsigned int subface, FullMatrix & interpolation_matrix, const unsigned int face_no) const { - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), + Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face(face_no))); + source_fe.n_dofs_per_face(face_no))); - // see if source is a Q element - if (const FE_Q_Base *source_fe = - dynamic_cast *>(&x_source_fe)) + // see if source is a Q or P element + if ((dynamic_cast *>(&source_fe) != nullptr) || + (dynamic_cast *>(&source_fe) != + nullptr)) { // have this test in here since a table of size 2x0 reports its size as // 0x0 @@ -646,18 +647,18 @@ FE_Q_Base::get_subface_interpolation_matrix( // produced in that case might lead to problems in the hp-procedures, // which use this method. Assert( - this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no), (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(face_no)); + source_fe.get_unit_face_support_points(face_no)); // 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 * q_degree * (dim - 1); + double eps = 2e-13 * this->q_degree * (dim - 1); // compute the interpolation matrix by simply taking the value at the // support points. @@ -673,7 +674,7 @@ FE_Q_Base::get_subface_interpolation_matrix( quad_face_support, 0, subface); - for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i) + for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i) { const Point &p = subface_quadrature.point(i); @@ -697,7 +698,7 @@ FE_Q_Base::get_subface_interpolation_matrix( #ifdef DEBUG // 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->n_dofs_per_face(face_no); ++j) + for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j) { double sum = 0.; @@ -708,7 +709,7 @@ FE_Q_Base::get_subface_interpolation_matrix( } #endif } - else if (dynamic_cast *>(&x_source_fe) != nullptr) + else if (dynamic_cast *>(&source_fe) != nullptr) { // nothing to do here, the FE_Nothing has no degrees of freedom anyway } diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc index 9bc477ace9..05a73ed207 100644 --- a/source/simplex/fe_lib.cc +++ b/source/simplex/fe_lib.cc @@ -458,19 +458,20 @@ namespace Simplex template void FE_Poly::get_face_interpolation_matrix( - const FiniteElement &x_source_fe, + const FiniteElement &source_fe, FullMatrix & interpolation_matrix, const unsigned int face_no) const { - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), + Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face(face_no))); + source_fe.n_dofs_per_face(face_no))); - if (const FE_Poly *source_fe = - dynamic_cast *>(&x_source_fe)) + // see if source is a P or Q element + if ((dynamic_cast *>(&source_fe) != nullptr) || + (dynamic_cast *>(&source_fe) != nullptr)) { const Quadrature quad_face_support( - source_fe->get_unit_face_support_points(face_no)); + source_fe.get_unit_face_support_points(face_no)); const double eps = 2e-13 * this->degree * (dim - 1); @@ -481,7 +482,7 @@ namespace Simplex face_no, face_quadrature_points); - for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i) + for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i) for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j) { double matrix_entry = @@ -500,7 +501,7 @@ namespace Simplex } #ifdef DEBUG - for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j) + for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j) { double sum = 0.; @@ -511,7 +512,7 @@ namespace Simplex } #endif } - else if (dynamic_cast *>(&x_source_fe) != nullptr) + else if (dynamic_cast *>(&source_fe) != nullptr) { // nothing to do here, the FE_Nothing has no degrees of freedom anyway } @@ -527,20 +528,21 @@ namespace Simplex template void FE_Poly::get_subface_interpolation_matrix( - const FiniteElement &x_source_fe, + const FiniteElement &source_fe, const unsigned int subface, FullMatrix & interpolation_matrix, const unsigned int face_no) const { - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), + Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face(face_no))); + source_fe.n_dofs_per_face(face_no))); - if (const FE_Poly *source_fe = - dynamic_cast *>(&x_source_fe)) + // see if source is a P or Q element + if ((dynamic_cast *>(&source_fe) != nullptr) || + (dynamic_cast *>(&source_fe) != nullptr)) { const Quadrature quad_face_support( - source_fe->get_unit_face_support_points(face_no)); + source_fe.get_unit_face_support_points(face_no)); const double eps = 2e-13 * this->degree * (dim - 1); @@ -552,7 +554,7 @@ namespace Simplex subface, subface_quadrature_points); - for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i) + for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i) for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j) { double matrix_entry = @@ -571,7 +573,7 @@ namespace Simplex } #ifdef DEBUG - for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j) + for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j) { double sum = 0.; @@ -582,7 +584,7 @@ namespace Simplex } #endif } - else if (dynamic_cast *>(&x_source_fe) != nullptr) + else if (dynamic_cast *>(&source_fe) != nullptr) { // nothing to do here, the FE_Nothing has no degrees of freedom anyway } -- 2.39.5