From: Marc Fehling Date: Fri, 15 Jan 2021 16:49:37 +0000 (-0700) Subject: simplex: face_interpolation for FE_Poly X-Git-Tag: v9.3.0-rc1~516^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11590%2Fhead;p=dealii.git simplex: face_interpolation for FE_Poly --- diff --git a/include/deal.II/simplex/fe_lib.h b/include/deal.II/simplex/fe_lib.h index f8b8688a89..be7ce1cbb5 100644 --- a/include/deal.II/simplex/fe_lib.h +++ b/include/deal.II/simplex/fe_lib.h @@ -62,6 +62,30 @@ namespace Simplex const RefinementCase &refinement_case = RefinementCase::isotropic_refinement) const override; + /** + * @copydoc dealii::FiniteElement::get_face_interpolation_matrix() + */ + void + get_face_interpolation_matrix(const FiniteElement &source_fe, + FullMatrix &interpolation_matrix, + const unsigned int face_no) const override; + + /** + * @copydoc dealii::FiniteElement::get_subface_interpolation_matrix() + */ + void + get_subface_interpolation_matrix( + const FiniteElement &x_source_fe, + const unsigned int subface, + FullMatrix & interpolation_matrix, + const unsigned int face_no) const override; + + /** + * @copydoc dealii::FiniteElement::hp_constraints_are_implemented() + */ + bool + hp_constraints_are_implemented() const override; + /** * @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values() */ diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index 6d05bf719b..148223d578 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -128,32 +128,57 @@ QProjector<2>::project_to_face(const ReferenceCell::Type reference_cell_type, const unsigned int face_no, std::vector> & q_points) { - Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented()); - (void)reference_cell_type; - const unsigned int dim = 2; AssertIndexRange(face_no, GeometryInfo::faces_per_cell); Assert(q_points.size() == quadrature.size(), ExcDimensionMismatch(q_points.size(), quadrature.size())); - for (unsigned int p = 0; p < quadrature.size(); ++p) - switch (face_no) - { - case 0: - q_points[p] = Point(0, quadrature.point(p)(0)); - break; - case 1: - q_points[p] = Point(1, quadrature.point(p)(0)); - break; - case 2: - q_points[p] = Point(quadrature.point(p)(0), 0); - break; - case 3: - q_points[p] = Point(quadrature.point(p)(0), 1); - break; - default: - Assert(false, ExcInternalError()); - } + if (reference_cell_type == ReferenceCell::Type::Tri) + { + // use linear polynomial to map the reference quadrature points correctly + // on faces, i.e., Simplex::ScalarPolynomial<1>(1) + for (unsigned int p = 0; p < quadrature.size(); ++p) + switch (face_no) + { + case 0: + q_points[p] = Point(quadrature.point(p)(0), 0); + break; + case 1: + q_points[p] = + Point(1 - quadrature.point(p)(0), quadrature.point(p)(0)); + break; + case 2: + q_points[p] = Point(0, 1 - quadrature.point(p)(0)); + break; + default: + Assert(false, ExcInternalError()); + } + } + else if (reference_cell_type == ReferenceCell::Type::Quad) + { + for (unsigned int p = 0; p < quadrature.size(); ++p) + switch (face_no) + { + case 0: + q_points[p] = Point(0, quadrature.point(p)(0)); + break; + case 1: + q_points[p] = Point(1, quadrature.point(p)(0)); + break; + case 2: + q_points[p] = Point(quadrature.point(p)(0), 0); + break; + case 3: + q_points[p] = Point(quadrature.point(p)(0), 1); + break; + default: + Assert(false, ExcInternalError()); + } + } + else + { + Assert(false, ExcInternalError()); + } } @@ -285,9 +310,6 @@ QProjector<2>::project_to_subface(const ReferenceCell::Type reference_cell_type, std::vector> & q_points, const RefinementCase<1> &) { - Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented()); - (void)reference_cell_type; - const unsigned int dim = 2; AssertIndexRange(face_no, GeometryInfo::faces_per_cell); AssertIndexRange(subface_no, GeometryInfo::max_children_per_face); @@ -295,65 +317,130 @@ QProjector<2>::project_to_subface(const ReferenceCell::Type reference_cell_type, Assert(q_points.size() == quadrature.size(), ExcDimensionMismatch(q_points.size(), quadrature.size())); - for (unsigned int p = 0; p < quadrature.size(); ++p) - switch (face_no) - { - case 0: - switch (subface_no) - { - case 0: - q_points[p] = Point(0, quadrature.point(p)(0) / 2); - break; - case 1: - q_points[p] = Point(0, quadrature.point(p)(0) / 2 + 0.5); - break; - default: - Assert(false, ExcInternalError()); - } - break; - case 1: - switch (subface_no) - { - case 0: - q_points[p] = Point(1, quadrature.point(p)(0) / 2); - break; - case 1: - q_points[p] = Point(1, quadrature.point(p)(0) / 2 + 0.5); - break; - default: - Assert(false, ExcInternalError()); - } - break; - case 2: - switch (subface_no) - { - case 0: - q_points[p] = Point(quadrature.point(p)(0) / 2, 0); - break; - case 1: - q_points[p] = Point(quadrature.point(p)(0) / 2 + 0.5, 0); - break; - default: - Assert(false, ExcInternalError()); - } - break; - case 3: - switch (subface_no) - { - case 0: - q_points[p] = Point(quadrature.point(p)(0) / 2, 1); - break; - case 1: - q_points[p] = Point(quadrature.point(p)(0) / 2 + 0.5, 1); - break; - default: - Assert(false, ExcInternalError()); - } - break; - - default: - Assert(false, ExcInternalError()); - } + if (reference_cell_type == ReferenceCell::Type::Tri) + { + // use linear polynomial to map the reference quadrature points correctly + // on faces, i.e., Simplex::ScalarPolynomial<1>(1) + for (unsigned int p = 0; p < quadrature.size(); ++p) + switch (face_no) + { + case 0: + switch (subface_no) + { + case 0: + q_points[p] = Point(quadrature.point(p)(0) / 2, 0); + break; + case 1: + q_points[p] = + Point(0.5 + quadrature.point(p)(0) / 2, 0); + break; + default: + Assert(false, ExcInternalError()); + } + break; + case 1: + switch (subface_no) + { + case 0: + q_points[p] = Point(1 - quadrature.point(p)(0) / 2, + quadrature.point(p)(0) / 2); + break; + case 1: + q_points[p] = Point(0.5 - quadrature.point(p)(0) / 2, + 0.5 + quadrature.point(p)(0) / 2); + break; + default: + Assert(false, ExcInternalError()); + } + break; + case 2: + switch (subface_no) + { + case 0: + q_points[p] = Point(0, 1 - quadrature.point(p)(0) / 2); + break; + case 1: + q_points[p] = + Point(0, 0.5 - quadrature.point(p)(0) / 2); + break; + default: + Assert(false, ExcInternalError()); + } + break; + default: + Assert(false, ExcInternalError()); + } + } + else if (reference_cell_type == ReferenceCell::Type::Quad) + { + for (unsigned int p = 0; p < quadrature.size(); ++p) + switch (face_no) + { + case 0: + switch (subface_no) + { + case 0: + q_points[p] = Point(0, quadrature.point(p)(0) / 2); + break; + case 1: + q_points[p] = + Point(0, quadrature.point(p)(0) / 2 + 0.5); + break; + default: + Assert(false, ExcInternalError()); + } + break; + case 1: + switch (subface_no) + { + case 0: + q_points[p] = Point(1, quadrature.point(p)(0) / 2); + break; + case 1: + q_points[p] = + Point(1, quadrature.point(p)(0) / 2 + 0.5); + break; + default: + Assert(false, ExcInternalError()); + } + break; + case 2: + switch (subface_no) + { + case 0: + q_points[p] = Point(quadrature.point(p)(0) / 2, 0); + break; + case 1: + q_points[p] = + Point(quadrature.point(p)(0) / 2 + 0.5, 0); + break; + default: + Assert(false, ExcInternalError()); + } + break; + case 3: + switch (subface_no) + { + case 0: + q_points[p] = Point(quadrature.point(p)(0) / 2, 1); + break; + case 1: + q_points[p] = + Point(quadrature.point(p)(0) / 2 + 0.5, 1); + break; + default: + Assert(false, ExcInternalError()); + } + break; + + default: + Assert(false, ExcInternalError()); + } + } + else + { + Assert(false, ExcInternalError()); + } } @@ -521,6 +608,17 @@ QProjector<2>::project_to_all_faces( { if (reference_cell_type == ReferenceCell::Type::Tri) { + const auto support_points_line = + [](const auto &face, const auto &orientation) -> std::vector> { + std::array, 2> vertices; + std::copy_n(face.first.begin(), face.first.size(), vertices.begin()); + const auto temp = + ReferenceCell::Type(ReferenceCell::Type::Line) + .permute_according_orientation(vertices, orientation); + return std::vector>(temp.begin(), + temp.begin() + face.first.size()); + }; + // reference faces (defined by its support points and arc length) const std::array, 2>, double>, 3> faces = { {{{{Point<2>(0.0, 0.0), Point<2>(1.0, 0.0)}}, 1.0}, @@ -542,21 +640,10 @@ QProjector<2>::project_to_all_faces( { const auto &face = faces[face_no]; - std::array, 2> support_points; - // determine support point of the current line with the correct // orientation - switch (orientation) - { - case 0: - support_points = {{face.first[1], face.first[0]}}; - break; - case 1: - support_points = {{face.first[0], face.first[1]}}; - break; - default: - Assert(false, ExcNotImplemented()); - } + std::vector> support_points = + support_points_line(face, orientation); // the quadrature rule to be projected ... const auto &sub_quadrature_points = diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index 77a219bca4..ddbaed6fef 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -697,6 +697,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) @@ -708,6 +709,7 @@ FE_Q_Base::get_subface_interpolation_matrix( Assert(std::fabs(sum - 1) < eps, ExcInternalError()); } +#endif } else if (dynamic_cast *>(&x_source_fe) != nullptr) { diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc index 61c14b5a43..2c31e505a8 100644 --- a/source/simplex/fe_lib.cc +++ b/source/simplex/fe_lib.cc @@ -15,6 +15,8 @@ #include +#include + #include #include #include @@ -370,6 +372,155 @@ namespace Simplex + template + void + FE_Poly::get_face_interpolation_matrix( + const FiniteElement &x_source_fe, + FullMatrix & interpolation_matrix, + const unsigned int face_no) const + { + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), + ExcDimensionMismatch(interpolation_matrix.m(), + x_source_fe.n_dofs_per_face(face_no))); + + if (const FE_Poly *source_fe = + dynamic_cast *>(&x_source_fe)) + { + const Quadrature quad_face_support( + source_fe->get_unit_face_support_points(face_no)); + + const double eps = 2e-13 * this->degree * (dim - 1); + + std::vector> face_quadrature_points( + quad_face_support.size()); + QProjector::project_to_face(this->reference_cell_type(), + quad_face_support, + face_no, + face_quadrature_points); + + 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 = + this->shape_value(this->face_to_cell_index(j, 0), + face_quadrature_points[i]); + + // 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; + } + +#ifdef DEBUG + for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j) + { + double sum = 0.; + + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) + sum += interpolation_matrix(j, i); + + Assert(std::fabs(sum - 1) < eps, ExcInternalError()); + } +#endif + } + else if (dynamic_cast *>(&x_source_fe) != nullptr) + { + // nothing to do here, the FE_Nothing has no degrees of freedom anyway + } + else + AssertThrow( + false, + (typename FiniteElement::ExcInterpolationNotImplemented())); + } + + + + template + void + FE_Poly::get_subface_interpolation_matrix( + const FiniteElement &x_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), + ExcDimensionMismatch(interpolation_matrix.m(), + x_source_fe.n_dofs_per_face(face_no))); + + if (const FE_Poly *source_fe = + dynamic_cast *>(&x_source_fe)) + { + const Quadrature quad_face_support( + source_fe->get_unit_face_support_points(face_no)); + + const double eps = 2e-13 * this->degree * (dim - 1); + + std::vector> subface_quadrature_points( + quad_face_support.size()); + QProjector::project_to_subface(this->reference_cell_type(), + quad_face_support, + face_no, + subface, + subface_quadrature_points); + + 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 = + this->shape_value(this->face_to_cell_index(j, 0), + subface_quadrature_points[i]); + + // 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; + } + +#ifdef DEBUG + for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j) + { + double sum = 0.; + + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) + sum += interpolation_matrix(j, i); + + Assert(std::fabs(sum - 1) < eps, ExcInternalError()); + } +#endif + } + else if (dynamic_cast *>(&x_source_fe) != nullptr) + { + // nothing to do here, the FE_Nothing has no degrees of freedom anyway + } + else + AssertThrow( + false, + (typename FiniteElement::ExcInterpolationNotImplemented())); + } + + + + template + bool + FE_Poly::hp_constraints_are_implemented() const + { + return true; + } + + + template void FE_Poly:: diff --git a/tests/simplex/hanging_nodes_01.cc b/tests/simplex/hanging_nodes_01.cc index 753ecdf45e..20762147c6 100644 --- a/tests/simplex/hanging_nodes_01.cc +++ b/tests/simplex/hanging_nodes_01.cc @@ -16,10 +16,23 @@ // verify hanging node constraints on locally h-refined simplex mesh +// +// dofs will be enumerated as follows +// scenario 1: +// 1-------0 +// |\ | +// | \ | +// 5---6 | +// |\ |\ | +// | \| \| +// 3---4---2 + #include #include +#include + #include #include @@ -31,6 +44,75 @@ #include "../tests.h" +// ----- diagnostics ----- + +template +void +print_dof_indices_on_faces(const DoFHandler &dofh) +{ + std::vector dof_indices; + + for (const auto &cell : dofh.active_cell_iterators()) + for (unsigned int f = 0; f < cell->n_faces(); ++f) + { + const auto &face = cell->face(f); + + if (face->has_children()) + { + for (unsigned int sf = 0; sf < face->n_children(); ++sf) + { + const auto &subface = face->child(sf); + Assert(subface->n_active_fe_indices() == 1, ExcInternalError()); + const unsigned int subface_fe_index = + subface->nth_active_fe_index(0); + const auto &subface_fe = dofh.get_fe(subface_fe_index); + + dof_indices.resize(subface_fe.n_dofs_per_face(f)); + subface->get_dof_indices(dof_indices, subface_fe_index); + + deallog << "cell:" << cell->active_cell_index() << " face:" << f + << " subface:" << sf << " dofs:"; + for (const auto &i : dof_indices) + deallog << i << " "; + deallog << std::endl; + } + } + else + { + Assert(face->n_active_fe_indices() == 1, ExcInternalError()); + const unsigned int face_fe_index = face->nth_active_fe_index(0); + const auto & face_fe = dofh.get_fe(face_fe_index); + + dof_indices.resize(face_fe.n_dofs_per_face(f)); + face->get_dof_indices(dof_indices, face_fe_index); + + deallog << "cell:" << cell->active_cell_index() << " face:" << f + << " dofs:"; + for (const auto &i : dof_indices) + deallog << i << " "; + deallog << std::endl; + } + } +} + + +template +void +print_dof_points(const DoFHandler &dofh) +{ + std::vector> points(dofh.n_dofs()); + DoFTools::map_dofs_to_support_points(MappingFE(dofh.get_fe()), + dofh, + points); + + for (unsigned int i = 0; i < dofh.n_dofs(); ++i) + deallog << "dof:" << i << " point:" << points[i] << std::endl; +} + + + +// ----- test ----- + template void test() @@ -51,9 +133,14 @@ test() dofh.distribute_dofs(Simplex::FE_P(1)); deallog << "ndofs: " << dofh.n_dofs() << std::endl; +#if false + print_dof_points(dofh); + print_dof_indices_on_faces(dofh); +#endif + // hanging node constraints AffineConstraints constraints; - // DoFTools::make_hanging_node_constraints(dofh, constraints); + DoFTools::make_hanging_node_constraints(dofh, constraints); constraints.print(deallog.get_file_stream()); deallog << "OK" << std::endl; diff --git a/tests/simplex/hanging_nodes_01.with_simplex_support=on.output b/tests/simplex/hanging_nodes_01.with_simplex_support=on.output index 1f03256372..8034a9e256 100644 --- a/tests/simplex/hanging_nodes_01.with_simplex_support=on.output +++ b/tests/simplex/hanging_nodes_01.with_simplex_support=on.output @@ -1,3 +1,5 @@ DEAL:2d::ndofs: 7 + 6 2: 0.500000 + 6 1: 0.500000 DEAL:2d::OK diff --git a/tests/simplex/hanging_nodes_02.cc b/tests/simplex/hanging_nodes_02.cc index d3816830b1..024848424e 100644 --- a/tests/simplex/hanging_nodes_02.cc +++ b/tests/simplex/hanging_nodes_02.cc @@ -16,13 +16,27 @@ // verify hanging node constraints on locally p-refined simplex mesh +// +// dofs will be enumerated as follows +// scenario 1: scenario 2: +// 6-------4 2---4---3 +// |\ | |\ | +// | \ | | \ | +// 3 2 | | 5 6 +// | \ | | \ | +// | \| | \| +// 0---1---5 0-------1 + #include #include +#include + #include #include +#include #include @@ -32,6 +46,54 @@ #include "../tests.h" +// ----- diagnostics ----- + +template +void +print_dof_indices_on_faces(const DoFHandler &dofh) +{ + std::vector dof_indices; + + for (const auto &cell : dofh.active_cell_iterators()) + for (unsigned int f = 0; f < cell->n_faces(); ++f) + { + const auto &face = cell->face(f); + + Assert(!face->has_children(), ExcInternalError()); + + const unsigned int fe_index = cell->active_fe_index(); + const auto & fe = cell->get_fe(); + + dof_indices.resize(fe.n_dofs_per_face(f)); + face->get_dof_indices(dof_indices, fe_index); + + deallog << "cell:" << cell->active_cell_index() << " face:" << f + << " dofs:"; + for (const auto &i : dof_indices) + deallog << i << " "; + deallog << std::endl; + } +} + + +template +void +print_dof_points(const DoFHandler &dofh) +{ + hp::MappingCollection mapping; + for (unsigned int i = 0; i < dofh.get_fe_collection().size(); ++i) + mapping.push_back(MappingFE(dofh.get_fe(i))); + + std::vector> points(dofh.n_dofs()); + DoFTools::map_dofs_to_support_points(mapping, dofh, points); + + for (unsigned int i = 0; i < dofh.n_dofs(); ++i) + deallog << "dof:" << i << " point:" << points[i] << std::endl; +} + + +// ----- test ----- + template void test(const hp::FECollection &fes) @@ -40,15 +102,25 @@ test(const hp::FECollection &fes) Triangulation tria; GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); +#if false + GridOut grid_out; + grid_out.write_vtk(tria, deallog.get_file_stream()); +#endif + DoFHandler dofh(tria); dofh.begin_active()->set_active_fe_index(1); dofh.distribute_dofs(fes); deallog << "ndofs: " << dofh.n_dofs() << std::endl; +#if false + print_dof_points(dofh); + print_dof_indices_on_faces(dofh); +#endif + // hanging node constraints AffineConstraints constraints; - // DoFTools::make_hanging_node_constraints(dofh, constraints); + DoFTools::make_hanging_node_constraints(dofh, constraints); constraints.print(deallog.get_file_stream()); deallog << "OK" << std::endl; diff --git a/tests/simplex/hanging_nodes_02.with_simplex_support=on.output b/tests/simplex/hanging_nodes_02.with_simplex_support=on.output index b905ef0c4e..0bc8c9bb95 100644 --- a/tests/simplex/hanging_nodes_02.with_simplex_support=on.output +++ b/tests/simplex/hanging_nodes_02.with_simplex_support=on.output @@ -1,5 +1,9 @@ DEAL:2d::ndofs: 7 + 2 6: 0.500000 + 2 5: 0.500000 DEAL:2d::OK DEAL:2d::ndofs: 7 + 5 1: 0.500000 + 5 2: 0.500000 DEAL:2d::OK