From: Peter Munch Date: Sun, 28 Feb 2021 06:13:53 +0000 (+0100) Subject: Implement FE_SimplexPoly::get_restriction_matrix() X-Git-Tag: v9.3.0-rc1~397^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11819%2Fhead;p=dealii.git Implement FE_SimplexPoly::get_restriction_matrix() --- diff --git a/include/deal.II/fe/fe_simplex_p.h b/include/deal.II/fe/fe_simplex_p.h index 43fb6b08a3..2d240489c6 100644 --- a/include/deal.II/fe/fe_simplex_p.h +++ b/include/deal.II/fe/fe_simplex_p.h @@ -60,6 +60,17 @@ public: const RefinementCase &refinement_case = RefinementCase::isotropic_refinement) const override; + /** + * @copydoc dealii::FiniteElement::get_restriction_matrix() + * + * @note Only implemented for RefinementCase::isotropic_refinement. + */ + virtual const FullMatrix & + get_restriction_matrix( + const unsigned int child, + const RefinementCase &refinement_case = + RefinementCase::isotropic_refinement) const override; + /** * @copydoc dealii::FiniteElement::get_face_interpolation_matrix() */ diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h index f96914254c..48512f06fc 100644 --- a/include/deal.II/fe/fe_tools.templates.h +++ b/include/deal.II/fe/fe_tools.templates.h @@ -2164,9 +2164,15 @@ namespace FETools const unsigned int nd = fe.n_components(); const unsigned int degree = fe.degree; + const ReferenceCell reference_cell = fe.reference_cell(); + + const auto &mapping = + reference_cell.template get_default_linear_mapping(); + const auto &q_fine = + reference_cell.get_gauss_type_quadrature(degree + 1); + // prepare FEValues, quadrature etc on // coarse cell - QGauss q_fine(degree + 1); const unsigned int nq = q_fine.size(); // create mass matrix on coarse cell. @@ -2174,9 +2180,10 @@ namespace FETools { // set up a triangulation for coarse cell Triangulation tr; - GridGenerator::hyper_cube(tr, 0, 1); + GridGenerator::reference_cell(reference_cell, tr); - FEValues coarse(fe, + FEValues coarse(mapping, + fe, q_fine, update_JxW_values | update_values); @@ -2212,9 +2219,10 @@ namespace FETools const auto compute_one_case = - [&fe, &q_fine, n, nd, nq](const unsigned int ref_case, - const FullMatrix &inverse_mass_matrix, - std::vector> &matrices) { + [&reference_cell, &mapping, &fe, &q_fine, n, nd, nq]( + const unsigned int ref_case, + const FullMatrix & inverse_mass_matrix, + std::vector> &matrices) { const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); @@ -2228,11 +2236,12 @@ namespace FETools // create a respective refinement on the triangulation Triangulation tr; - GridGenerator::hyper_cube(tr, 0, 1); + GridGenerator::reference_cell(reference_cell, tr); tr.begin_active()->set_refine_flag(RefinementCase(ref_case)); tr.execute_coarsening_and_refinement(); - FEValues fine(fe, + FEValues fine(mapping, + fe, q_fine, update_quadrature_points | update_JxW_values | update_values); diff --git a/source/fe/fe_simplex_p.cc b/source/fe/fe_simplex_p.cc index 9f32a2f293..fe9b39cab5 100644 --- a/source/fe/fe_simplex_p.cc +++ b/source/fe/fe_simplex_p.cc @@ -354,6 +354,48 @@ FE_SimplexPoly::get_prolongation_matrix( +template +const FullMatrix & +FE_SimplexPoly::get_restriction_matrix( + const unsigned int child, + const RefinementCase &refinement_case) const +{ + Assert(refinement_case == RefinementCase::isotropic_refinement, + ExcNotImplemented()); + AssertDimension(dim, spacedim); + + // initialization upon first request + if (this->restriction[refinement_case - 1][child].n() == 0) + { + std::lock_guard lock(this->mutex); + + // if matrix got updated while waiting for the lock + if (this->restriction[refinement_case - 1][child].n() == + this->n_dofs_per_cell()) + return this->restriction[refinement_case - 1][child]; + + // now do the work. need to get a non-const version of data in order to + // be able to modify them inside a const function + auto &this_nonconst = const_cast &>(*this); + + std::vector>> isotropic_matrices( + RefinementCase::isotropic_refinement); + isotropic_matrices.back().resize( + GeometryInfo::n_children(RefinementCase(refinement_case)), + FullMatrix(this->n_dofs_per_cell(), this->n_dofs_per_cell())); + + FETools::compute_projection_matrices(*this, isotropic_matrices, true); + + this_nonconst.restriction[refinement_case - 1].swap( + isotropic_matrices.back()); + } + + // finally return the matrix + return this->restriction[refinement_case - 1][child]; +} + + + template void FE_SimplexPoly::get_face_interpolation_matrix(