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<dim, spacedim> tr;
- GridGenerator::hyper_cube(tr);
+ ReferenceCell::make_triangulation(reference_cell_type, tr);
+
+ const auto &mapping =
+ ReferenceCell::get_default_linear_mapping<dim, spacedim>(
+ 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<dim> quadrature(degree + 1);
+ const auto quadrature =
+ ReferenceCell::get_gauss_type_quadrature<dim>(reference_cell_type,
+ degree + 1);
// Set up FEValues.
const UpdateFlags flags =
update_values | update_quadrature_points | update_JxW_values;
- FEValues<dim> val1(fe1, quadrature, update_values);
+ FEValues<dim> val1(mapping, fe1, quadrature, update_values);
val1.reinit(tr.begin_active());
- FEValues<dim> val2(fe2, quadrature, flags);
+ FEValues<dim> val2(mapping, fe2, quadrature, flags);
val2.reinit(tr.begin_active());
// Integrate and invert mass matrix. This happens in the target space
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/lac/householder.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+class RightHandSideFunction : public Function<dim>
+{
+public:
+ RightHandSideFunction(const unsigned int n_components)
+ : Function<dim>(n_components)
+ {}
+
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const
+ {
+ if (component == 0)
+ return p[0];
+ else
+ return 0.0;
+ }
+};
+
+template <int dim, int spacedim = dim>
+void
+test()
+{
+ Simplex::FE_P<dim, spacedim> fe_coarse(1);
+ Simplex::FE_P<dim, spacedim> fe_fine(2);
+ MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+
+ FullMatrix<double> 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<dim, spacedim> tria;
+ GridGenerator::subdivided_hyper_cube_with_simplices(
+ tria, Utilities::pow(2, n_refinements));
+
+ DoFHandler<dim, spacedim> dof_handler_coarse(tria);
+ dof_handler_coarse.distribute_dofs(fe_coarse);
+
+ DoFHandler<dim, spacedim> dof_handler_fine(tria);
+ dof_handler_fine.distribute_dofs(fe_fine);
+
+ Vector<double> vec_coarse(dof_handler_coarse.n_dofs());
+ Vector<double> vec_fine(dof_handler_fine.n_dofs());
+
+ Vector<double> temp_coarse(fe_coarse.n_dofs_per_cell());
+ Vector<double> temp_fine(fe_fine.n_dofs_per_cell());
+
+ VectorTools::interpolate(mapping,
+ dof_handler_coarse,
+ RightHandSideFunction<dim>(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<dim, spacedim, false>(&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<dim> 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<dim> 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>();
+}
projects / dealii.git / commitdiff