From 4f987bdfdfc5be8239ee0899faf7e10eb09ca58a Mon Sep 17 00:00:00 2001 From: Magdalena Schreter Date: Thu, 18 Jan 2024 22:15:10 +0100 Subject: [PATCH] allow interpolate for combination FENothing/FE_DGQ --- doc/news/changes/minor/20240118Schreter | 4 + source/fe/fe_dgq.cc | 145 +++++++++++++----------- tests/fe/dgq_1.cc | 30 +++++ tests/fe/dgq_1.output | 6 + tests/hp/solution_transfer_03.cc | 2 +- 5 files changed, 121 insertions(+), 66 deletions(-) create mode 100644 doc/news/changes/minor/20240118Schreter diff --git a/doc/news/changes/minor/20240118Schreter b/doc/news/changes/minor/20240118Schreter new file mode 100644 index 0000000000..3021b4753a --- /dev/null +++ b/doc/news/changes/minor/20240118Schreter @@ -0,0 +1,4 @@ +Changed: +FE_DGQ::get_interpolation_matrix() now also takes as a source element FENothing. +
+(Magdalena Schreter, Peter Munch, 2024/01/18) diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc index 3f24eb3875..f1d15b8667 100644 --- a/source/fe/fe_dgq.cc +++ b/source/fe/fe_dgq.cc @@ -274,78 +274,93 @@ FE_DGQ::get_interpolation_matrix( const FiniteElement &x_source_fe, FullMatrix &interpolation_matrix) const { - // this is only implemented, if the - // source FE is also a - // DGQ element - using FE = FiniteElement; - AssertThrow((dynamic_cast *>(&x_source_fe) != - nullptr), - typename FE::ExcInterpolationNotImplemented()); - - // ok, source is a Q element, so - // we will be able to do the work - const FE_DGQ &source_fe = - dynamic_cast &>(x_source_fe); - - Assert(interpolation_matrix.m() == this->n_dofs_per_cell(), - ExcDimensionMismatch(interpolation_matrix.m(), - this->n_dofs_per_cell())); - Assert(interpolation_matrix.n() == source_fe.n_dofs_per_cell(), - ExcDimensionMismatch(interpolation_matrix.n(), - source_fe.n_dofs_per_cell())); - - - // compute the interpolation - // matrices in much the same way as - // we do for the embedding matrices - // from mother to child. - FullMatrix cell_interpolation(this->n_dofs_per_cell(), - this->n_dofs_per_cell()); - FullMatrix source_interpolation(this->n_dofs_per_cell(), - source_fe.n_dofs_per_cell()); - FullMatrix tmp(this->n_dofs_per_cell(), source_fe.n_dofs_per_cell()); - for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j) + // go through the list of elements we can interpolate from + if (const FE_DGQ *source_fe = + dynamic_cast *>(&x_source_fe)) { - // generate a point on this - // cell and evaluate the - // shape functions there - const Point p = this->unit_support_points[j]; - for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i) - cell_interpolation(j, i) = this->poly_space->compute_value(i, p); - - for (unsigned int i = 0; i < source_fe.n_dofs_per_cell(); ++i) - source_interpolation(j, i) = source_fe.poly_space->compute_value(i, p); - } + Assert(interpolation_matrix.m() == this->n_dofs_per_cell(), + ExcDimensionMismatch(interpolation_matrix.m(), + this->n_dofs_per_cell())); + Assert(interpolation_matrix.n() == source_fe->n_dofs_per_cell(), + ExcDimensionMismatch(interpolation_matrix.n(), + source_fe->n_dofs_per_cell())); + + + // compute the interpolation + // matrices in much the same way as + // we do for the embedding matrices + // from mother to child. + FullMatrix cell_interpolation(this->n_dofs_per_cell(), + this->n_dofs_per_cell()); + FullMatrix source_interpolation(this->n_dofs_per_cell(), + source_fe->n_dofs_per_cell()); + FullMatrix tmp(this->n_dofs_per_cell(), + source_fe->n_dofs_per_cell()); + for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j) + { + // generate a point on this + // cell and evaluate the + // shape functions there + const Point p = this->unit_support_points[j]; + for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i) + cell_interpolation(j, i) = this->poly_space->compute_value(i, p); + + for (unsigned int i = 0; i < source_fe->n_dofs_per_cell(); ++i) + source_interpolation(j, i) = + source_fe->poly_space->compute_value(i, p); + } - // then compute the - // interpolation matrix matrix - // for this coordinate - // direction - cell_interpolation.gauss_jordan(); - cell_interpolation.mmult(interpolation_matrix, source_interpolation); + // then compute the + // interpolation matrix matrix + // for this coordinate + // direction + cell_interpolation.gauss_jordan(); + cell_interpolation.mmult(interpolation_matrix, source_interpolation); - // cut off very small values - for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i) - for (unsigned int j = 0; j < source_fe.n_dofs_per_cell(); ++j) - if (std::fabs(interpolation_matrix(i, j)) < 1e-15) - interpolation_matrix(i, j) = 0.; + // cut off very small values + for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i) + for (unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j) + if (std::fabs(interpolation_matrix(i, j)) < 1e-15) + interpolation_matrix(i, j) = 0.; #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 i = 0; i < this->n_dofs_per_cell(); ++i) - { - double sum = 0.; - for (unsigned int j = 0; j < source_fe.n_dofs_per_cell(); ++j) - sum += interpolation_matrix(i, j); + // 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 i = 0; i < this->n_dofs_per_cell(); ++i) + { + double sum = 0.; + for (unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j) + sum += interpolation_matrix(i, j); - Assert(std::fabs(sum - 1) < 5e-14 * std::max(this->degree, 1U) * dim, - ExcInternalError()); - } + Assert(std::fabs(sum - 1) < 5e-14 * std::max(this->degree, 1U) * dim, + ExcInternalError()); + } #endif + } + else if (dynamic_cast *>(&x_source_fe)) + { + // the element we want to interpolate from is an FE_Nothing. this + // element represents a function that is constant zero and has no + // degrees of freedom, so the interpolation is simply a multiplication + // with a n_dofs x 0 matrix. there is nothing to do here + + // we would like to verify that the number of rows and columns of + // the matrix equals this->n_dofs_per_cell() and zero. unfortunately, + // whenever we do FullMatrix::reinit(m,0), it sets both rows and + // columns to zero, instead of m and zero. thus, only test the + // number of columns + Assert(interpolation_matrix.n() == x_source_fe.n_dofs_per_cell(), + ExcDimensionMismatch(interpolation_matrix.m(), + x_source_fe.n_dofs_per_cell())); + } + else + AssertThrow( + false, + (typename FiniteElement::ExcInterpolationNotImplemented())); } diff --git a/tests/fe/dgq_1.cc b/tests/fe/dgq_1.cc index 879f0967c1..974ad0fea1 100644 --- a/tests/fe/dgq_1.cc +++ b/tests/fe/dgq_1.cc @@ -21,6 +21,7 @@ // result doesn't change #include +#include #include #include @@ -56,6 +57,31 @@ test(const unsigned int degree1, const unsigned int degree2) } +template +void +test_fe_nothing() +{ + deallog << "FE_DGQ<" << dim << ">(1)" + << " to FE_Nothing<" << dim << ">()" << std::endl; + + FE_DGQ fe1(1); + FE_Nothing fe2; + + FullMatrix m; + fe1.get_interpolation_matrix(fe2, m); + + for (unsigned int i = 0; i < m.m(); ++i) + { + for (unsigned int j = 0; j < m.n(); ++j) + deallog << m(i, j) << ' '; + + deallog << std::endl; + } + + deallog << std::endl; +} + + int main() { @@ -74,5 +100,9 @@ main() for (unsigned int degree2 = 0; degree2 <= 2; ++degree2) test<3>(degree1, degree2); + test_fe_nothing<1>(); + test_fe_nothing<2>(); + test_fe_nothing<3>(); + return 0; } diff --git a/tests/fe/dgq_1.output b/tests/fe/dgq_1.output index 8686f9a711..1b6ca44b96 100644 --- a/tests/fe/dgq_1.output +++ b/tests/fe/dgq_1.output @@ -402,3 +402,9 @@ DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0. DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 DEAL:: +DEAL::FE_DGQ<1>(1) to FE_Nothing<1>() +DEAL:: +DEAL::FE_DGQ<2>(1) to FE_Nothing<2>() +DEAL:: +DEAL::FE_DGQ<3>(1) to FE_Nothing<3>() +DEAL:: diff --git a/tests/hp/solution_transfer_03.cc b/tests/hp/solution_transfer_03.cc index f2f1c52ae6..d4e94555b4 100644 --- a/tests/hp/solution_transfer_03.cc +++ b/tests/hp/solution_transfer_03.cc @@ -91,7 +91,7 @@ main() data_out.write_vtu(deallog.get_file_stream()); - // Interpoalte solution + // Interpolate solution SolutionTransfer<2, Vector> solultion_trans(dof_handler); solultion_trans.prepare_for_coarsening_and_refinement(solution); -- 2.39.5