From b216b97109f8467151044d8f68c482a0d5dea2fe Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 28 Dec 2020 11:59:23 +0100 Subject: [PATCH] FETools::get_projection_matrix for simplices --- include/deal.II/fe/fe_tools.templates.h | 19 ++- tests/simplex/get_projection_matrix_01.cc | 139 ++++++++++++++++++ ...n_matrix_01.with_simplex_support=on.output | 3 + 3 files changed, 157 insertions(+), 4 deletions(-) create mode 100644 tests/simplex/get_projection_matrix_01.cc create mode 100644 tests/simplex/get_projection_matrix_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 7833e04cba..0d00c9062a 100644 --- a/include/deal.II/fe/fe_tools.templates.h +++ b/include/deal.II/fe/fe_tools.templates.h @@ -1598,23 +1598,34 @@ namespace FETools const unsigned int n1 = fe1.n_dofs_per_cell(); const unsigned int n2 = fe2.n_dofs_per_cell(); + const auto reference_cell_type = fe1.reference_cell_type(); + + Assert(fe1.reference_cell_type() == fe2.reference_cell_type(), + ExcNotImplemented()); + // First, create a local mass matrix for the unit cell Triangulation tr; - GridGenerator::hyper_cube(tr); + ReferenceCell::make_triangulation(reference_cell_type, tr); + + const auto &mapping = + ReferenceCell::get_default_linear_mapping( + reference_cell_type); // Choose a Gauss quadrature rule that is exact up to degree 2n-1 const unsigned int degree = std::max(fe1.tensor_degree(), fe2.tensor_degree()); Assert(degree != numbers::invalid_unsigned_int, ExcNotImplemented()); - const QGauss quadrature(degree + 1); + const auto quadrature = + ReferenceCell::get_gauss_type_quadrature(reference_cell_type, + degree + 1); // Set up FEValues. const UpdateFlags flags = update_values | update_quadrature_points | update_JxW_values; - FEValues val1(fe1, quadrature, update_values); + FEValues val1(mapping, fe1, quadrature, update_values); val1.reinit(tr.begin_active()); - FEValues val2(fe2, quadrature, flags); + FEValues val2(mapping, fe2, quadrature, flags); val2.reinit(tr.begin_active()); // Integrate and invert mass matrix. This happens in the target space diff --git a/tests/simplex/get_projection_matrix_01.cc b/tests/simplex/get_projection_matrix_01.cc new file mode 100644 index 0000000000..095fee6f5e --- /dev/null +++ b/tests/simplex/get_projection_matrix_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 FETools::get_projection_matrix 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; + } +}; + +template +void +test() +{ + Simplex::FE_P fe_coarse(1); + Simplex::FE_P fe_fine(2); + MappingFE mapping(Simplex::FE_P(1)); + + FullMatrix matrix(fe_fine.n_dofs_per_cell(), + fe_coarse.n_dofs_per_cell()); + FETools::get_projection_matrix(fe_coarse, fe_fine, matrix); + + const unsigned int n_refinements = 2; + + Triangulation tria; + GridGenerator::subdivided_hyper_cube_with_simplices( + tria, Utilities::pow(2, n_refinements)); + + DoFHandler dof_handler_coarse(tria); + dof_handler_coarse.distribute_dofs(fe_coarse); + + DoFHandler dof_handler_fine(tria); + dof_handler_fine.distribute_dofs(fe_fine); + + Vector vec_coarse(dof_handler_coarse.n_dofs()); + Vector vec_fine(dof_handler_fine.n_dofs()); + + Vector temp_coarse(fe_coarse.n_dofs_per_cell()); + Vector temp_fine(fe_fine.n_dofs_per_cell()); + + VectorTools::interpolate(mapping, + dof_handler_coarse, + RightHandSideFunction(1), + vec_coarse); + + for (const auto &cell_coarse : dof_handler_coarse.active_cell_iterators()) + { + cell_coarse->get_dof_values(vec_coarse, temp_coarse); + + matrix.vmult(temp_fine, temp_coarse); + + DoFCellAccessor(&tria, + cell_coarse->level(), + cell_coarse->index(), + &dof_handler_fine) + .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 +} + +int +main() +{ + initlog(); + + test<2>(); + test<3>(); +} diff --git a/tests/simplex/get_projection_matrix_01.with_simplex_support=on.output b/tests/simplex/get_projection_matrix_01.with_simplex_support=on.output new file mode 100644 index 0000000000..4ccd48e0a2 --- /dev/null +++ b/tests/simplex/get_projection_matrix_01.with_simplex_support=on.output @@ -0,0 +1,3 @@ + +-2.082e-17 2.500e-01 -2.082e-17 1.250e-01 1.250e-01 0.000e+00 2.500e-01 1.250e-01 2.500e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 -2.082e-17 1.250e-01 0.000e+00 2.500e-01 1.250e-01 2.500e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 -2.082e-17 1.250e-01 0.000e+00 2.500e-01 1.250e-01 2.500e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 -1.110e-16 1.250e-01 0.000e+00 2.500e-01 1.250e-01 2.500e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 +-3.331e-16 2.500e-01 -7.772e-16 -7.772e-16 1.250e-01 1.250e-01 -1.890e-16 -5.095e-17 1.250e-01 -3.454e-18 2.500e-01 2.500e-01 2.500e-01 1.250e-01 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -1.159e-16 -1.890e-16 1.250e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 3.750e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 8.750e-01 8.750e-01 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -7.772e-16 1.250e-01 -1.890e-16 -5.439e-16 1.250e-01 -5.095e-17 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -7.772e-16 1.250e-01 -1.890e-16 -3.454e-18 2.500e-01 2.500e-01 2.500e-01 1.250e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -1.159e-16 -1.890e-16 1.250e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -5.274e-16 -2.109e-15 1.250e-01 -1.890e-16 -5.439e-16 1.250e-01 -5.551e-17 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 2.500e-01 1.250e-01 2.500e-01 2.500e-01 -7.772e-16 -5.439e-16 1.250e-01 -5.095e-17 -3.331e-16 -1.159e-16 1.250e-01 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -1.159e-16 -1.890e-16 1.250e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 -7.772e-16 -5.439e-16 1.250e-01 -5.095e-17 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -5.274e-16 -3.308e-16 -1.890e-16 1.250e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 -7.772e-16 -5.095e-17 1.250e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -1.159e-16 -1.890e-16 1.250e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 8.750e-01 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 -7.772e-16 -5.439e-16 1.250e-01 -5.095e-17 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -3.331e-16 -1.159e-16 -1.890e-16 1.250e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 -2.109e-15 -5.439e-16 1.250e-01 -5.551e-17 -1.890e-16 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 2.500e-01 1.250e-01 2.500e-01 2.500e-01 -2.096e-15 -5.439e-16 1.250e-01 -5.095e-17 -8.882e-16 -1.159e-16 1.250e-01 -2.220e-17 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 8.750e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -8.882e-16 -1.159e-16 -2.220e-17 1.250e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 2.500e-01 -2.096e-15 -5.439e-16 1.250e-01 -5.095e-17 -2.220e-17 1.250e-01 2.500e-01 2.500e-01 2.500e-01 1.250e-01 5.000e-01 3.750e-01 5.000e-01 5.000e-01 3.750e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 6.250e-01 7.500e-01 7.500e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 1.000e+00 1.000e+00 8.750e-01 1.000e+00 8.750e-01 -3.454e-18 2.500e-01 1.250e-01 2.500e-01 2.500e-01 1.250e-01 2.500e-01 -1.499e-15 -3.308e-16 -2.220e-17 1.250e-01 5.000e-01 3.750e-01 3.750e-01 5.000e-01 5.000e-01 5.000e-01 3.750e-01 7.500e-01 6.250e-01 7.500e-01 7.500e-01 6.250e-01 7.500e-01 6.250e-01 1.000e+00 8.750e-01 8.750e-01 1.000e+00 1.000e+00 1.000e+00 8.750e-01 -- 2.39.5