const unsigned int n = fe.n_dofs_per_cell();
const unsigned int nc =
GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
+
+ AssertDimension(matrices.size(), nc);
+
for (unsigned int i = 0; i < nc; ++i)
{
Assert(matrices[i].n() == n,
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<dim, spacedim> tria;
- GridGenerator::hyper_cube(tria, 0, 1);
+ ReferenceCell::make_triangulation(reference_cell_type, tria);
tria.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
tria.execute_coarsening_and_refinement();
const unsigned int degree = fe.degree;
- QGauss<dim> q_fine(degree + 1);
+
+ const auto &mapping =
+ ReferenceCell::get_default_linear_mapping<dim, spacedim>(
+ reference_cell_type);
+ const auto &q_fine =
+ ReferenceCell::get_gauss_type_quadrature<dim>(reference_cell_type,
+ degree + 1);
const unsigned int nq = q_fine.size();
- FEValues<dim, spacedim> fine(fe,
+ FEValues<dim, spacedim> fine(mapping,
+ fe,
q_fine,
update_quadrature_points |
update_JxW_values | update_values);
q_points_coarse[i](j) = q_points_fine[i](j);
const Quadrature<dim> q_coarse(q_points_coarse,
fine.get_JxW_values());
- FEValues<dim, spacedim> coarse(fe, q_coarse, update_values);
+ FEValues<dim, spacedim> coarse(mapping,
+ fe,
+ q_coarse,
+ update_values);
coarse.reinit(tria.begin(0));
template <int dim, typename number, int spacedim>
void
- compute_embedding_matrices(const FiniteElement<dim, spacedim> &fe,
- std::vector<std::vector<FullMatrix<number>>
-
- > & matrices,
- const bool isotropic_only,
- const double threshold)
+ compute_embedding_matrices(
+ const FiniteElement<dim, spacedim> & fe,
+ std::vector<std::vector<FullMatrix<number>>> &matrices,
+ const bool isotropic_only,
+ const double threshold)
{
Threads::TaskGroup<void> task_group;
std::pair<Table<2, bool>, std::vector<unsigned int>>
get_constant_modes() const override;
+ /**
+ * @copydoc dealii::FiniteElement::get_prolongation_matrix()
+ *
+ * @note Only implemented for RefinementCase::isotropic_refinement.
+ */
+ const FullMatrix<double> &
+ get_prolongation_matrix(
+ const unsigned int child,
+ const RefinementCase<dim> &refinement_case =
+ RefinementCase<dim>::isotropic_refinement) const override;
+
private:
/**
* @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values()
convert_generalized_support_point_values_to_dof_values(
const std::vector<Vector<double>> &support_point_values,
std::vector<double> & nodal_values) const override;
+
+ mutable Threads::Mutex mutex;
};
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_tools.h>
#include <deal.II/simplex/fe_lib.h>
+ template <int dim, int spacedim>
+ const FullMatrix<double> &
+ FE_Poly<dim, spacedim>::get_prolongation_matrix(
+ const unsigned int child,
+ const RefinementCase<dim> &refinement_case) const
+ {
+ Assert(refinement_case == RefinementCase<dim>::isotropic_refinement,
+ ExcNotImplemented());
+ AssertDimension(dim, spacedim);
+
+ // initialization upon first request
+ if (this->prolongation[refinement_case - 1][child].n() == 0)
+ {
+ std::lock_guard<std::mutex> 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<FE_Poly<dim, spacedim> &>(*this);
+
+ std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
+ RefinementCase<dim>::isotropic_refinement);
+ isotropic_matrices.back().resize(
+ GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)),
+ FullMatrix<double>(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 <int dim, int spacedim>
void
FE_Poly<dim, spacedim>::
--- /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 Simplex::FE_Poly::get_prolongation_matrix()
+// (and indirectly FETools::compute_embedding_matrices() 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;
+ }
+};
+
+int
+main()
+{
+ initlog();
+
+ const int dim = 2;
+ const int spacedim = 2;
+
+ Simplex::FE_P<dim, spacedim> fe(2);
+ MappingFE<dim> mapping(Simplex::FE_P<dim>(1));
+
+ const unsigned int n_refinements = 2;
+
+ Triangulation<dim, spacedim> 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<dim, spacedim> dof_handler_coarse(tria_coarse);
+ dof_handler_coarse.distribute_dofs(fe);
+
+ DoFHandler<dim, spacedim> dof_handler_fine(tria_fine);
+ dof_handler_fine.distribute_dofs(fe);
+
+ Vector<double> vec_coarse(dof_handler_coarse.n_dofs());
+ Vector<double> vec_fine(dof_handler_fine.n_dofs());
+
+ Vector<double> temp_coarse(fe.n_dofs_per_cell());
+ Vector<double> temp_fine(fe.n_dofs_per_cell());
+
+ // interpolate function onto coarse grid.
+ VectorTools::interpolate(mapping,
+ dof_handler_coarse,
+ RightHandSideFunction<dim>(1),
+ vec_coarse);
+
+ // project the result onto fine grid (cell by cell)
+ for (const auto &cell_coarse : dof_handler_coarse.active_cell_iterators())
+ {
+ DoFCellAccessor<dim, spacedim, false> 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<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
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff