From 08182a4b05e9107dba58fef96ec892eeffb6c889 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Fri, 19 Apr 2019 10:17:13 +0200 Subject: [PATCH] Mapping::get_center() --- doc/news/changes/minor/20190411LucaHeltai | 4 + include/deal.II/fe/mapping.h | 26 ++++++ source/fe/mapping.cc | 25 ++++++ tests/mappings/mapping_q_eulerian_12.cc | 96 +++++++++++++++++++++ tests/mappings/mapping_q_eulerian_12.output | 10 +++ 5 files changed, 161 insertions(+) create mode 100644 doc/news/changes/minor/20190411LucaHeltai create mode 100644 tests/mappings/mapping_q_eulerian_12.cc create mode 100644 tests/mappings/mapping_q_eulerian_12.output diff --git a/doc/news/changes/minor/20190411LucaHeltai b/doc/news/changes/minor/20190411LucaHeltai new file mode 100644 index 0000000000..ad5ad6c437 --- /dev/null +++ b/doc/news/changes/minor/20190411LucaHeltai @@ -0,0 +1,4 @@ +New: The method Mapping::get_center() allows one to retrieve the cell center of a cell when the +mapping object does not preserve vertex locations. +
+(Luca Heltai, 2019/04/11) diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index 21bec9ebc2..c1c60aac8b 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -336,6 +336,32 @@ public: get_vertices( const typename Triangulation::cell_iterator &cell) const; + /** + * Return the mapped center of a cell. + * + * If you are using a (bi-,tri-)linear mapping that preserves vertex + * locations, this function simply returns the value also produced by + * `cell->center()`. However, there are also mappings that add displacements + * or choose completely different locations, e.g., MappingQEulerian, + * MappingQ1Eulerian, or MappingFEField. + * + * By default, this function returns the average of the vertex locations + * returned by the the get_vertices() method. If the parameter + * @p map_center_of_reference_cell is set to true, than a more expensive + * algorithm is used, that returns the mapped center of the reference cell, + * obtained by the transform_unit_to_real_cell() method. + * + * @param[in] cell The cell for which you want to compute the center + * @param[in] map_center_of_reference_cell A flag that switches the algorithm + * for the computation of the cell center from vertex averages to + * transform_unit_to_real_cell() applied to the center of the reference cell + * + * @author Luca Heltai, 2019. + */ + virtual Point + get_center(const typename Triangulation::cell_iterator &cell, + const bool map_center_of_reference_cell = false) const; + /** * Return whether the mapping preserves vertex locations. In other words, * this function returns whether the mapped location of the reference cell diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc index 60d5795800..97a7a38f23 100644 --- a/source/fe/mapping.cc +++ b/source/fe/mapping.cc @@ -36,6 +36,31 @@ Mapping::get_vertices( +template +Point +Mapping::get_center( + const typename Triangulation::cell_iterator &cell, + const bool map_center_of_reference_cell) const +{ + if (map_center_of_reference_cell) + { + Point reference_center; + for (unsigned int d = 0; d < dim; ++d) + reference_center[d] = .5; + return transform_unit_to_real_cell(cell, reference_center); + } + else + { + const auto vertices = get_vertices(cell); + Point center; + for (const auto &v : vertices) + center += v; + return center / GeometryInfo::vertices_per_cell; + } +} + + + template Point Mapping::project_real_point_to_unit_point_on_face( diff --git a/tests/mappings/mapping_q_eulerian_12.cc b/tests/mappings/mapping_q_eulerian_12.cc new file mode 100644 index 0000000000..e8626585c5 --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_12.cc @@ -0,0 +1,96 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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 that MappingQEulerian knows how to compute centers of cells + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; + + +template +class Displacement : public Function +{ +public: + Displacement() + : Function(dim) + {} + + double + value(const Point &p, const unsigned int component) const + { + return p[component] * p[0] / 2; + } +}; + + +template +void +test() +{ + deallog << "dim=" << dim << std::endl; + + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation, -1, 1); + + FESystem fe(FE_Q(2), dim); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + Vector displacements(dof_handler.n_dofs()); + + VectorTools::interpolate(dof_handler, Displacement(), displacements); + + MappingQEulerian euler(2, dof_handler, displacements); + + // now the actual test + for (const auto &cell : dof_handler.active_cell_iterators()) + { + deallog << "Center: [" << cell->center() << "], with mapping [" + << euler.get_center(cell) << "]" << std::endl; + deallog << "Center: [" << cell->center() << "], with mapping + flag [" + << euler.get_center(cell, true) << "]" << std::endl; + } +} + + + +int +main() +{ + initlog(1); + + test<1>(); + test<2>(); + test<3>(); +} diff --git a/tests/mappings/mapping_q_eulerian_12.output b/tests/mappings/mapping_q_eulerian_12.output new file mode 100644 index 0000000000..eb04a22389 --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_12.output @@ -0,0 +1,10 @@ + +DEAL::dim=1 +DEAL::Center: [0.00000], with mapping [0.500000] +DEAL::Center: [0.00000], with mapping + flag [0.00000] +DEAL::dim=2 +DEAL::Center: [0.00000 0.00000], with mapping [0.500000 0.00000] +DEAL::Center: [0.00000 0.00000], with mapping + flag [0.00000 0.00000] +DEAL::dim=3 +DEAL::Center: [0.00000 0.00000 0.00000], with mapping [0.500000 0.00000 0.00000] +DEAL::Center: [0.00000 0.00000 0.00000], with mapping + flag [0.00000 0.00000 0.00000] -- 2.39.5