From ff123c1db68656885b40d3ef4d33d5d1d661ed7a Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 15 Jan 2021 20:47:58 +0100 Subject: [PATCH] Work on Work on MappingFEField for simplex meshes --- include/deal.II/fe/mapping_fe_field.h | 4 +- include/deal.II/grid/reference_cell.h | 10 ++ source/fe/mapping_fe_field.cc | 81 +++++------- source/grid/reference_cell.cc | 41 ++++++ source/grid/reference_cell.inst.in | 2 + tests/simplex/mapping_fe_fields_01.cc | 121 ++++++++++++++++++ ...e_fields_01.with_simplex_support=on.output | 30 +++++ 7 files changed, 239 insertions(+), 50 deletions(-) create mode 100644 tests/simplex/mapping_fe_fields_01.cc create mode 100644 tests/simplex/mapping_fe_fields_01.with_simplex_support=on.output diff --git a/include/deal.II/fe/mapping_fe_field.h b/include/deal.II/fe/mapping_fe_field.h index 494f5bc8ac..7f36f19cf9 100644 --- a/include/deal.II/fe/mapping_fe_field.h +++ b/include/deal.II/fe/mapping_fe_field.h @@ -13,8 +13,8 @@ // // --------------------------------------------------------------------- -#ifndef dealii_mapping_fe_h -#define dealii_mapping_fe_h +#ifndef dealii_mapping_fe_field_h +#define dealii_mapping_fe_field_h #include diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 019d41062a..2d8161feef 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -506,6 +506,16 @@ namespace ReferenceCell get_gauss_type_quadrature(const Type & reference_cell, const unsigned n_points_1D); + /** + * Return a quadrature rule with the support points of the given reference + * cell. + * + * @note The weights are not filled. + */ + template + Quadrature & + get_nodal_type_quadrature(const Type &reference_cell); + namespace internal { /** diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index 1d75c345e1..d440769a9a 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -192,28 +192,6 @@ MappingFEField::InternalData:: -namespace -{ - // Construct a quadrature formula containing the vertices of the reference - // cell in dimension dim (with invalid weights) - template - Quadrature & - get_vertex_quadrature() - { - static Quadrature quad; - if (quad.size() == 0) - { - std::vector> points(GeometryInfo::vertices_per_cell); - for (const unsigned int i : GeometryInfo::vertex_indices()) - points[i] = GeometryInfo::unit_cell_vertex(i); - quad = Quadrature(points); - } - return quad; - } -} // namespace - - - template MappingFEField::MappingFEField( const DoFHandler &euler_dof_handler, @@ -229,7 +207,8 @@ MappingFEField::MappingFEField( true)) , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int) , fe_values(this->euler_dof_handler->get_fe(), - get_vertex_quadrature(), + ReferenceCell::get_nodal_type_quadrature( + this->euler_dof_handler->get_fe().reference_cell_type()), update_values) { unsigned int size = 0; @@ -257,7 +236,8 @@ MappingFEField::MappingFEField( true)) , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int) , fe_values(this->euler_dof_handler->get_fe(), - get_vertex_quadrature(), + ReferenceCell::get_nodal_type_quadrature( + this->euler_dof_handler->get_fe().reference_cell_type()), update_values) { unsigned int size = 0; @@ -296,7 +276,8 @@ MappingFEField::MappingFEField( true)) , fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int) , fe_values(this->euler_dof_handler->get_fe(), - get_vertex_quadrature(), + ReferenceCell::get_nodal_type_quadrature( + this->euler_dof_handler->get_fe().reference_cell_type()), update_values) { unsigned int size = 0; @@ -332,7 +313,8 @@ MappingFEField::MappingFEField( , fe_mask(mapping.fe_mask) , fe_to_real(mapping.fe_to_real) , fe_values(mapping.euler_dof_handler->get_fe(), - get_vertex_quadrature(), + ReferenceCell::get_nodal_type_quadrature( + this->euler_dof_handler->get_fe().reference_cell_type()), update_values) {} @@ -372,8 +354,7 @@ MappingFEField::get_vertices( *cell, euler_dof_handler); Assert(uses_level_dofs || dof_cell->is_active() == true, ExcInactiveCell()); - AssertDimension(GeometryInfo::vertices_per_cell, - fe_values.n_quadrature_points); + AssertDimension(cell->n_vertices(), fe_values.n_quadrature_points); AssertDimension(fe_to_real.size(), euler_dof_handler->get_fe().n_components()); if (uses_level_dofs) @@ -402,7 +383,7 @@ MappingFEField::get_vertices( boost::container::small_vector, GeometryInfo::vertices_per_cell> - vertices(GeometryInfo::vertices_per_cell); + vertices(cell->n_vertices()); for (unsigned int i = 0; i < dofs_per_cell; ++i) { const unsigned int comp = fe_to_real @@ -412,7 +393,7 @@ MappingFEField::get_vertices( typename VectorType::value_type value = internal::ElementAccess::get(vector, dof_indices[i]); if (euler_dof_handler->get_fe().is_primitive(i)) - for (const unsigned int v : GeometryInfo::vertex_indices()) + for (const unsigned int v : cell->vertex_indices()) vertices[v][comp] += fe_values.shape_value(i, v) * value; else Assert(false, ExcNotImplemented()); @@ -585,23 +566,30 @@ MappingFEField::compute_face_data( data.aux.resize( dim - 1, std::vector>(n_original_q_points)); + + // TODO: only a single reference cell type possible... + const auto reference_cell_type = + this->euler_dof_handler->get_fe().reference_cell_type(); + const auto n_faces = + ReferenceCell::internal::Info::get_cell(reference_cell_type) + .n_faces(); + // Compute tangentials to the unit cell. - for (const unsigned int i : GeometryInfo::face_indices()) + for (unsigned int i = 0; i < n_faces; ++i) { data.unit_tangentials[i].resize(n_original_q_points); std::fill(data.unit_tangentials[i].begin(), data.unit_tangentials[i].end(), - GeometryInfo::unit_tangential_vectors[i][0]); + ReferenceCell::unit_tangential_vectors( + reference_cell_type, i, 0)); if (dim > 2) { - data.unit_tangentials[GeometryInfo::faces_per_cell + i] - .resize(n_original_q_points); - std::fill( - data.unit_tangentials[GeometryInfo::faces_per_cell + i] - .begin(), - data.unit_tangentials[GeometryInfo::faces_per_cell + i] - .end(), - GeometryInfo::unit_tangential_vectors[i][1]); + data.unit_tangentials[n_faces + i].resize( + n_original_q_points); + std::fill(data.unit_tangentials[n_faces + i].begin(), + data.unit_tangentials[n_faces + i].end(), + ReferenceCell::unit_tangential_vectors( + reference_cell_type, i, 1)); } } } @@ -1273,22 +1261,17 @@ namespace internal // 0. for (unsigned int d = 0; d != dim - 1; ++d) { - Assert(face_no + GeometryInfo::faces_per_cell * d < + Assert(face_no + cell->n_faces() * d < data.unit_tangentials.size(), ExcInternalError()); Assert( data.aux[d].size() <= - data - .unit_tangentials[face_no + - GeometryInfo::faces_per_cell * d] - .size(), + data.unit_tangentials[face_no + cell->n_faces() * d].size(), ExcInternalError()); mapping.transform( make_array_view( - data - .unit_tangentials[face_no + - GeometryInfo::faces_per_cell * d]), + data.unit_tangentials[face_no + cell->n_faces() * d]), mapping_contravariant, data, make_array_view(data.aux[d])); @@ -1371,6 +1354,7 @@ namespace internal if (subface_no != numbers::invalid_unsigned_int) { + // TODO const double area_ratio = GeometryInfo::subface_ratio( cell->subface_case(face_no), subface_no); @@ -2096,6 +2080,7 @@ MappingFEField::transform_real_to_unit_cell( initial_p_unit[d] = 0.5; } + // TODO initial_p_unit = GeometryInfo::project_to_unit_cell(initial_p_unit); // for (unsigned int d=0; d(); // never reached } + template + Quadrature & + get_nodal_type_quadrature(const Type &reference_cell) + { + AssertDimension(dim, get_dimension(reference_cell)); + + const auto create_quadrature = [](const Type &reference_cell) { + Triangulation tria; + make_triangulation(reference_cell, tria); + + return Quadrature(tria.get_vertices()); + }; + + if (reference_cell == get_hypercube(dim)) + { + static Quadrature quadrature = create_quadrature(reference_cell); + return quadrature; + } + else if (reference_cell == Type::Tri || reference_cell == Type::Tet) + { + static Quadrature quadrature = create_quadrature(reference_cell); + return quadrature; + } + else if (reference_cell == Type::Pyramid) + { + static Quadrature quadrature = create_quadrature(reference_cell); + return quadrature; + } + else if (reference_cell == Type::Wedge) + { + static Quadrature quadrature = create_quadrature(reference_cell); + return quadrature; + } + else + Assert(false, ExcNotImplemented()); + + static Quadrature dummy; + + return dummy; // never reached + } + #include "reference_cell.inst" } // namespace ReferenceCell diff --git a/source/grid/reference_cell.inst.in b/source/grid/reference_cell.inst.in index 198079f766..37f9309819 100644 --- a/source/grid/reference_cell.inst.in +++ b/source/grid/reference_cell.inst.in @@ -34,4 +34,6 @@ for (deal_II_dimension : DIMENSIONS) { template Quadrature get_gauss_type_quadrature( const Type &reference_cell, const unsigned n_points_1D); + template Quadrature &get_nodal_type_quadrature( + const Type &reference_cell); } diff --git a/tests/simplex/mapping_fe_fields_01.cc b/tests/simplex/mapping_fe_fields_01.cc new file mode 100644 index 0000000000..4f96c4a0fc --- /dev/null +++ b/tests/simplex/mapping_fe_fields_01.cc @@ -0,0 +1,121 @@ +// --------------------------------------------------------------------- +// +// 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 MappingFEField for simplex meshes. + +#include + +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include + +#include "../tests.h" + +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() +{ + const int dim = 2; + + Triangulation tria; + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 1); + + Simplex::FE_P fe(2); + FESystem euler_fe(fe, dim); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + DoFHandler euler_dof_handler(tria); + euler_dof_handler.distribute_dofs(euler_fe); + + 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); + + MappingFEField mapping(euler_dof_handler, euler_vector); + + Simplex::QGauss quadrature_formula(1); + + FEValues fe_values(mapping, + fe, + quadrature_formula, + update_values | update_gradients); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + } + + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handler); + + Vector solution(dof_handler.n_dofs()); + data_out.add_data_vector(solution, "solution"); + + data_out.build_patches(mapping, 2); + +#if false + std::ofstream output("test.vtk"); + data_out.write_vtk(output); +#else + data_out.write_vtk(deallog.get_file_stream()); +#endif + } +} + +int +main() +{ + initlog(); + + test(); +} \ No newline at end of file diff --git a/tests/simplex/mapping_fe_fields_01.with_simplex_support=on.output b/tests/simplex/mapping_fe_fields_01.with_simplex_support=on.output new file mode 100644 index 0000000000..146f92e589 --- /dev/null +++ b/tests/simplex/mapping_fe_fields_01.with_simplex_support=on.output @@ -0,0 +1,30 @@ + +# vtk DataFile Version 3.0 +#This file was generated +ASCII +DATASET UNSTRUCTURED_GRID + +POINTS 12 double +0.00000 0.00000 0 +1.00000 0.00000 0 +0.00000 1.00000 0 +0.500000 0.00000 0 +0.500000 0.500000 0 +0.00000 0.500000 0 +1.00000 1.00000 0 +0.00000 1.00000 0 +1.00000 0.00000 0 +0.500000 1.00000 0 +0.500000 0.500000 0 +1.00000 0.500000 0 + +CELLS 2 14 + 6 0 1 2 3 4 5 + 6 6 7 8 9 10 11 + +CELL_TYPES 2 + 22 22 +POINT_DATA 12 +SCALARS solution double 1 +LOOKUP_TABLE default +0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -- 2.39.5