From 718c17e3dbe895dd771acd7b9b0f3937767b923e Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 26 Dec 2020 22:55:06 +0100 Subject: [PATCH] Work on compute_embedding_matrices for simplex --- include/deal.II/fe/fe_tools.templates.h | 34 +++-- include/deal.II/simplex/fe_lib.h | 13 ++ source/simplex/fe_lib.cc | 43 ++++++ .../simplex/compute_projection_matrices_01.cc | 139 ++++++++++++++++++ ...matrices_01.with_simplex_support=on.output | 2 + 5 files changed, 221 insertions(+), 10 deletions(-) create mode 100644 tests/simplex/compute_projection_matrices_01.cc create mode 100644 tests/simplex/compute_projection_matrices_01.with_simplex_support=on.output diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h index 0d00c9062a..1443e66baa 100644 --- a/include/deal.II/fe/fe_tools.templates.h +++ b/include/deal.II/fe/fe_tools.templates.h @@ -1799,6 +1799,9 @@ namespace FETools const unsigned int n = fe.n_dofs_per_cell(); const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); + + AssertDimension(matrices.size(), nc); + for (unsigned int i = 0; i < nc; ++i) { Assert(matrices[i].n() == n, @@ -1807,18 +1810,27 @@ namespace FETools ExcDimensionMismatch(matrices[i].m(), n)); } + const auto reference_cell_type = fe.reference_cell_type(); + // Set up meshes, one with a single // reference cell and refine it once Triangulation tria; - GridGenerator::hyper_cube(tria, 0, 1); + ReferenceCell::make_triangulation(reference_cell_type, tria); tria.begin_active()->set_refine_flag(RefinementCase(ref_case)); tria.execute_coarsening_and_refinement(); const unsigned int degree = fe.degree; - QGauss q_fine(degree + 1); + + const auto &mapping = + ReferenceCell::get_default_linear_mapping( + reference_cell_type); + const auto &q_fine = + ReferenceCell::get_gauss_type_quadrature(reference_cell_type, + degree + 1); const unsigned int nq = q_fine.size(); - FEValues fine(fe, + FEValues fine(mapping, + fe, q_fine, update_quadrature_points | update_JxW_values | update_values); @@ -1865,7 +1877,10 @@ namespace FETools q_points_coarse[i](j) = q_points_fine[i](j); const Quadrature q_coarse(q_points_coarse, fine.get_JxW_values()); - FEValues coarse(fe, q_coarse, update_values); + FEValues coarse(mapping, + fe, + q_coarse, + update_values); coarse.reinit(tria.begin(0)); @@ -1921,12 +1936,11 @@ namespace FETools template void - compute_embedding_matrices(const FiniteElement &fe, - std::vector> - - > & matrices, - const bool isotropic_only, - const double threshold) + compute_embedding_matrices( + const FiniteElement & fe, + std::vector>> &matrices, + const bool isotropic_only, + const double threshold) { Threads::TaskGroup task_group; diff --git a/include/deal.II/simplex/fe_lib.h b/include/deal.II/simplex/fe_lib.h index 0d5fe5cc8e..e26dd20b20 100644 --- a/include/deal.II/simplex/fe_lib.h +++ b/include/deal.II/simplex/fe_lib.h @@ -51,6 +51,17 @@ namespace Simplex std::pair, std::vector> get_constant_modes() const override; + /** + * @copydoc dealii::FiniteElement::get_prolongation_matrix() + * + * @note Only implemented for RefinementCase::isotropic_refinement. + */ + const FullMatrix & + get_prolongation_matrix( + const unsigned int child, + const RefinementCase &refinement_case = + RefinementCase::isotropic_refinement) const override; + private: /** * @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values() @@ -59,6 +70,8 @@ namespace Simplex convert_generalized_support_point_values_to_dof_values( const std::vector> &support_point_values, std::vector & nodal_values) const override; + + mutable Threads::Mutex mutex; }; diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc index 587b91e83c..ca3141b1cd 100644 --- a/source/simplex/fe_lib.cc +++ b/source/simplex/fe_lib.cc @@ -17,6 +17,7 @@ #include #include +#include #include @@ -304,6 +305,48 @@ namespace Simplex + template + const FullMatrix & + FE_Poly::get_prolongation_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->prolongation[refinement_case - 1][child].n() == 0) + { + std::lock_guard lock(this->mutex); + + // if matrix got updated while waiting for the lock + if (this->prolongation[refinement_case - 1][child].n() == + this->n_dofs_per_cell()) + return this->prolongation[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_embedding_matrices(*this, isotropic_matrices, true); + + this_nonconst.prolongation[refinement_case - 1].swap( + isotropic_matrices.back()); + } + + // finally return the matrix + return this->prolongation[refinement_case - 1][child]; + } + + + template void FE_Poly:: diff --git a/tests/simplex/compute_projection_matrices_01.cc b/tests/simplex/compute_projection_matrices_01.cc new file mode 100644 index 0000000000..02d85f17b5 --- /dev/null +++ b/tests/simplex/compute_projection_matrices_01.cc @@ -0,0 +1,139 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test Simplex::FE_Poly::get_prolongation_matrix() +// (and indirectly FETools::compute_embedding_matrices() for simplices). + + +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +template +class RightHandSideFunction : public Function +{ +public: + RightHandSideFunction(const unsigned int n_components) + : Function(n_components) + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const + { + if (component == 0) + return p[0]; + else + return 0.0; + } +}; + +int +main() +{ + initlog(); + + const int dim = 2; + const int spacedim = 2; + + Simplex::FE_P fe(2); + MappingFE mapping(Simplex::FE_P(1)); + + const unsigned int n_refinements = 2; + + Triangulation tria_coarse, tria_fine; + GridGenerator::subdivided_hyper_cube_with_simplices(tria_coarse, 1); + tria_coarse.refine_global(n_refinements); + GridGenerator::subdivided_hyper_cube_with_simplices(tria_fine, 1); + tria_fine.refine_global(n_refinements + 1); + + DoFHandler dof_handler_coarse(tria_coarse); + dof_handler_coarse.distribute_dofs(fe); + + DoFHandler dof_handler_fine(tria_fine); + dof_handler_fine.distribute_dofs(fe); + + Vector vec_coarse(dof_handler_coarse.n_dofs()); + Vector vec_fine(dof_handler_fine.n_dofs()); + + Vector temp_coarse(fe.n_dofs_per_cell()); + Vector temp_fine(fe.n_dofs_per_cell()); + + // interpolate function onto coarse grid. + VectorTools::interpolate(mapping, + dof_handler_coarse, + RightHandSideFunction(1), + vec_coarse); + + // project the result onto fine grid (cell by cell) + for (const auto &cell_coarse : dof_handler_coarse.active_cell_iterators()) + { + DoFCellAccessor cell_coarse_on_fine_tria( + &tria_fine, + cell_coarse->level(), + cell_coarse->index(), + &dof_handler_fine); + + cell_coarse->get_dof_values(vec_coarse, temp_coarse); + + for (unsigned int c = 0; c < cell_coarse_on_fine_tria.n_children(); ++c) + { + fe.get_prolongation_matrix(c).vmult(temp_fine, temp_coarse); + cell_coarse_on_fine_tria.child(c)->set_dof_values(temp_fine, + vec_fine); + } + } + + vec_fine.print(deallog.get_file_stream()); + +#if false + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handler_coarse); + data_out.add_data_vector(vec_coarse, "solution"); + + data_out.build_patches(mapping, 2); + + std::ofstream output("test_coarse.vtk"); + data_out.write_vtk(output); + } + + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handler_fine); + data_out.add_data_vector(vec_fine, "solution"); + + data_out.build_patches(mapping, 2); + + std::ofstream output("test_fine.vtk"); + data_out.write_vtk(output); + } +#endif +} diff --git a/tests/simplex/compute_projection_matrices_01.with_simplex_support=on.output b/tests/simplex/compute_projection_matrices_01.with_simplex_support=on.output new file mode 100644 index 0000000000..d29b5b081d --- /dev/null +++ b/tests/simplex/compute_projection_matrices_01.with_simplex_support=on.output @@ -0,0 +1,2 @@ + +0.000e+00 1.250e-01 0.000e+00 6.250e-02 6.250e-02 0.000e+00 2.500e-01 1.250e-01 1.875e-01 1.875e-01 1.250e-01 0.000e+00 6.250e-02 6.250e-02 0.000e+00 3.750e-01 2.500e-01 3.125e-01 3.125e-01 2.500e-01 5.000e-01 3.750e-01 4.375e-01 4.375e-01 3.750e-01 2.500e-01 3.125e-01 3.125e-01 2.500e-01 1.250e-01 0.000e+00 6.250e-02 6.250e-02 0.000e+00 1.250e-01 1.875e-01 1.875e-01 1.250e-01 0.000e+00 6.250e-02 6.250e-02 0.000e+00 1.875e-01 1.875e-01 1.250e-01 6.250e-01 5.000e-01 5.625e-01 5.625e-01 5.000e-01 7.500e-01 6.250e-01 6.875e-01 6.875e-01 6.250e-01 5.000e-01 5.625e-01 5.625e-01 5.000e-01 8.750e-01 7.500e-01 8.125e-01 8.125e-01 7.500e-01 1.000e+00 8.750e-01 9.375e-01 9.375e-01 8.750e-01 7.500e-01 8.125e-01 8.125e-01 7.500e-01 6.250e-01 5.000e-01 5.625e-01 5.625e-01 5.000e-01 6.250e-01 6.875e-01 6.875e-01 6.250e-01 5.000e-01 5.625e-01 5.625e-01 5.000e-01 6.875e-01 6.875e-01 6.250e-01 1.250e-01 0.000e+00 6.250e-02 6.250e-02 0.000e+00 2.500e-01 1.250e-01 1.875e-01 1.875e-01 1.250e-01 0.000e+00 6.250e-02 6.250e-02 0.000e+00 3.750e-01 2.500e-01 3.125e-01 3.125e-01 2.500e-01 3.750e-01 4.375e-01 4.375e-01 3.750e-01 2.500e-01 3.125e-01 3.125e-01 2.500e-01 1.250e-01 0.000e+00 6.250e-02 6.250e-02 0.000e+00 1.250e-01 1.875e-01 1.875e-01 1.250e-01 0.000e+00 6.250e-02 6.250e-02 0.000e+00 1.875e-01 1.875e-01 1.250e-01 4.375e-01 3.750e-01 4.375e-01 4.375e-01 3.750e-01 3.125e-01 3.750e-01 4.375e-01 4.375e-01 4.375e-01 3.750e-01 3.125e-01 2.500e-01 2.500e-01 1.875e-01 2.500e-01 1.875e-01 1.250e-01 3.750e-01 3.125e-01 3.125e-01 1.000e+00 8.750e-01 1.000e+00 9.375e-01 9.375e-01 1.000e+00 7.500e-01 8.750e-01 8.125e-01 8.125e-01 8.750e-01 1.000e+00 9.375e-01 9.375e-01 1.000e+00 6.250e-01 7.500e-01 6.875e-01 6.875e-01 7.500e-01 5.000e-01 6.250e-01 5.625e-01 5.625e-01 6.250e-01 7.500e-01 6.875e-01 6.875e-01 7.500e-01 8.750e-01 1.000e+00 9.375e-01 9.375e-01 1.000e+00 8.750e-01 8.125e-01 8.125e-01 8.750e-01 1.000e+00 9.375e-01 9.375e-01 1.000e+00 8.125e-01 8.125e-01 8.750e-01 3.750e-01 5.000e-01 4.375e-01 4.375e-01 5.000e-01 2.500e-01 3.750e-01 3.125e-01 3.125e-01 3.750e-01 5.000e-01 4.375e-01 4.375e-01 5.000e-01 1.250e-01 2.500e-01 1.875e-01 1.875e-01 2.500e-01 6.250e-02 1.250e-01 1.875e-01 2.500e-01 3.750e-01 5.000e-01 4.375e-01 4.375e-01 5.000e-01 3.125e-01 3.750e-01 4.375e-01 5.000e-01 3.125e-01 3.125e-01 3.750e-01 8.750e-01 1.000e+00 9.375e-01 9.375e-01 1.000e+00 7.500e-01 8.750e-01 8.125e-01 8.125e-01 8.750e-01 1.000e+00 9.375e-01 9.375e-01 1.000e+00 6.250e-01 7.500e-01 6.875e-01 6.875e-01 7.500e-01 5.625e-01 6.250e-01 6.875e-01 7.500e-01 8.750e-01 1.000e+00 9.375e-01 9.375e-01 1.000e+00 8.125e-01 8.750e-01 9.375e-01 1.000e+00 8.125e-01 8.125e-01 8.750e-01 5.625e-01 6.250e-01 5.625e-01 5.625e-01 6.250e-01 6.875e-01 6.250e-01 5.625e-01 5.625e-01 5.625e-01 6.250e-01 6.875e-01 7.500e-01 7.500e-01 8.125e-01 7.500e-01 8.125e-01 8.750e-01 6.250e-01 6.875e-01 6.875e-01 -- 2.39.5