From 589854d4f5bffc952fd132d36931b919b54c8611 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 17 Jan 2021 11:22:17 +0100 Subject: [PATCH] Generalize VectorTools::get_position_vector for simplex meshes --- include/deal.II/grid/reference_cell.h | 12 ++++++++ .../vector_tools_interpolate.templates.h | 8 +++-- source/grid/reference_cell.cc | 26 +++++++++++++++++ source/grid/reference_cell.inst.in | 4 +++ tests/simplex/mapping_fe_fields_01.cc | 29 +++---------------- 5 files changed, 52 insertions(+), 27 deletions(-) diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index ca56d20f89..fe79ef35ef 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -507,6 +507,18 @@ namespace ReferenceCell make_triangulation(const Type & reference_cell, Triangulation &tria); + /** + * Return a default mapping of degree @ degree matching the given reference + * cell. If this reference cell is a hypercube, then the returned mapping is a + * MappingQGeneric; otherwise, it is an object of type MappingFE initialized + * with Simplex::FE_P (if the reference cell is a triangle and tetrahedron), + * with Simplex::FE_PyramidP (if the reference cell is a pyramid), or with + * Simplex::FE_WedgeP (if the reference cell is a wedge). + */ + template + std::unique_ptr> + get_default_mapping(const Type &reference_cell, const unsigned int degree); + /** * Return a default linear mapping matching the given reference cell. * If this reference cell is a hypercube, then the returned mapping diff --git a/include/deal.II/numerics/vector_tools_interpolate.templates.h b/include/deal.II/numerics/vector_tools_interpolate.templates.h index 7edd39de41..80138b7c19 100644 --- a/include/deal.II/numerics/vector_tools_interpolate.templates.h +++ b/include/deal.II/numerics/vector_tools_interpolate.templates.h @@ -653,8 +653,12 @@ namespace VectorTools { const Quadrature quad(fe.get_unit_support_points()); - MappingQGeneric map_q(fe.degree); - FEValues fe_v(map_q, fe, quad, update_quadrature_points); + const auto map_q = ReferenceCell::get_default_mapping( + fe.reference_cell_type(), fe.degree); + FEValues fe_v(*map_q, + fe, + quad, + update_quadrature_points); std::vector dofs(fe.n_dofs_per_cell()); AssertDimension(fe.n_dofs_per_cell(), diff --git a/source/grid/reference_cell.cc b/source/grid/reference_cell.cc index f670dde1f8..68a599be18 100644 --- a/source/grid/reference_cell.cc +++ b/source/grid/reference_cell.cc @@ -102,6 +102,32 @@ namespace ReferenceCell } + template + std::unique_ptr> + get_default_mapping(const Type &reference_cell, const unsigned int degree) + { + AssertDimension(dim, get_dimension(reference_cell)); + + if (reference_cell == get_hypercube(dim)) + return std::make_unique>(degree); + else if (reference_cell == Type::Tri || reference_cell == Type::Tet) + return std::make_unique>( + Simplex::FE_P(degree)); + else if (reference_cell == Type::Pyramid) + return std::make_unique>( + Simplex::FE_PyramidP(degree)); + else if (reference_cell == Type::Wedge) + return std::make_unique>( + Simplex::FE_WedgeP(degree)); + else + { + Assert(false, ExcNotImplemented()); + } + + return std::make_unique>(degree); + } + + template const Mapping & get_default_linear_mapping(const Type &reference_cell) diff --git a/source/grid/reference_cell.inst.in b/source/grid/reference_cell.inst.in index 37f9309819..29eec7c46d 100644 --- a/source/grid/reference_cell.inst.in +++ b/source/grid/reference_cell.inst.in @@ -20,6 +20,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const Type & reference_cell, Triangulation &tria); + template std::unique_ptr< + Mapping> + get_default_mapping(const Type &reference_cell, const unsigned int degree); + template const Mapping &get_default_linear_mapping(const Type &reference_cell); diff --git a/tests/simplex/mapping_fe_fields_01.cc b/tests/simplex/mapping_fe_fields_01.cc index 4f96c4a0fc..9168411635 100644 --- a/tests/simplex/mapping_fe_fields_01.cc +++ b/tests/simplex/mapping_fe_fields_01.cc @@ -14,7 +14,8 @@ // --------------------------------------------------------------------- -// Test MappingFEField for simplex meshes. +// Test MappingFEField and VectorTools::get_position_vector() for simplex +// meshes. #include @@ -36,21 +37,6 @@ using namespace dealii; -template -class Solution : public Function -{ -public: - Solution() - : Function(dim) - {} - - double - value(const Point &point, const unsigned int compontent) const - { - return point[compontent]; - } -}; - void test() { @@ -59,7 +45,7 @@ test() Triangulation tria; GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); - Simplex::FE_P fe(2); + Simplex::FE_P fe(1); FESystem euler_fe(fe, dim); DoFHandler dof_handler(tria); @@ -70,14 +56,7 @@ test() Vector euler_vector(euler_dof_handler.n_dofs()); - // TODO: not working (missing mapping) - // VectorTools::get_position_vector(euler_dof_handler, euler_vector); - - MappingFE mapping_interpolation(Simplex::FE_P(1)); - VectorTools::interpolate(mapping_interpolation, - euler_dof_handler, - Solution(), - euler_vector); + VectorTools::get_position_vector(euler_dof_handler, euler_vector); MappingFEField mapping(euler_dof_handler, euler_vector); -- 2.39.5