From 8586540af6efd4a62d7c4acba7f096d5c7c08bbf Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 10 Oct 2020 23:05:49 +0200 Subject: [PATCH] Extend GridTools::maximal_cell_diameter and ::minimal_cell_diameter for simplex --- .../changes/incompatibilities/20201016Munch | 4 ++++ include/deal.II/fe/mapping.h | 3 ++- include/deal.II/fe/mapping_fe_field.h | 3 ++- include/deal.II/fe/mapping_q1_eulerian.h | 3 ++- include/deal.II/fe/mapping_q_eulerian.h | 6 +++-- source/fe/mapping.cc | 14 +++++++----- source/fe/mapping_fe_field.cc | 7 ++++-- source/fe/mapping_q1_eulerian.cc | 10 +++++---- source/fe/mapping_q_eulerian.cc | 15 +++++++------ source/fe/mapping_q_generic.cc | 9 ++++++-- source/grid/grid_tools.cc | 22 +++++++++++++++---- source/opencascade/utilities.cc | 3 +-- tests/simplex/cell_measure_01.cc | 12 ++++++++++ ..._measure_01.with_simplex_support=on.output | 4 ++++ 14 files changed, 83 insertions(+), 32 deletions(-) create mode 100644 doc/news/changes/incompatibilities/20201016Munch diff --git a/doc/news/changes/incompatibilities/20201016Munch b/doc/news/changes/incompatibilities/20201016Munch new file mode 100644 index 0000000000..eeb889e55f --- /dev/null +++ b/doc/news/changes/incompatibilities/20201016Munch @@ -0,0 +1,4 @@ +Changed: The return type of the method Mapping::get_vertices() has been changed +from std::array to boost::container::small_vector. +
+(Peter Munch, 2020/10/16) diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index 8b302d8ba7..aec5119630 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -332,7 +332,8 @@ public: * information stored by the triangulation, i.e., * cell-@>vertex(v). */ - virtual std::array, GeometryInfo::vertices_per_cell> + virtual boost::container::small_vector, + GeometryInfo::vertices_per_cell> get_vertices( const typename Triangulation::cell_iterator &cell) const; diff --git a/include/deal.II/fe/mapping_fe_field.h b/include/deal.II/fe/mapping_fe_field.h index 691549b4d4..5c2eae1d94 100644 --- a/include/deal.II/fe/mapping_fe_field.h +++ b/include/deal.II/fe/mapping_fe_field.h @@ -169,7 +169,8 @@ public: * and constructs the position of the vertices according to the @p euler_vector * that was passed at construction time. */ - virtual std::array, GeometryInfo::vertices_per_cell> + virtual boost::container::small_vector, + GeometryInfo::vertices_per_cell> get_vertices(const typename Triangulation::cell_iterator &cell) const override; diff --git a/include/deal.II/fe/mapping_q1_eulerian.h b/include/deal.II/fe/mapping_q1_eulerian.h index 531f904a82..6a9d77df97 100644 --- a/include/deal.II/fe/mapping_q1_eulerian.h +++ b/include/deal.II/fe/mapping_q1_eulerian.h @@ -116,7 +116,8 @@ public: * cell but instead evaluates an externally given displacement field in * addition to the geometry of the cell. */ - virtual std::array, GeometryInfo::vertices_per_cell> + virtual boost::container::small_vector, + GeometryInfo::vertices_per_cell> get_vertices(const typename Triangulation::cell_iterator &cell) const override; diff --git a/include/deal.II/fe/mapping_q_eulerian.h b/include/deal.II/fe/mapping_q_eulerian.h index 7e5b5adab8..10bd13b13f 100644 --- a/include/deal.II/fe/mapping_q_eulerian.h +++ b/include/deal.II/fe/mapping_q_eulerian.h @@ -120,7 +120,8 @@ public: * cell but instead evaluates an externally given displacement field in * addition to the geometry of the cell. */ - virtual std::array, GeometryInfo::vertices_per_cell> + virtual boost::container::small_vector, + GeometryInfo::vertices_per_cell> get_vertices(const typename Triangulation::cell_iterator &cell) const override; @@ -204,7 +205,8 @@ private: * current cell but instead evaluates an externally given displacement * field in addition to the geometry of the cell. */ - virtual std::array, GeometryInfo::vertices_per_cell> + virtual boost::container::small_vector, + GeometryInfo::vertices_per_cell> get_vertices(const typename Triangulation::cell_iterator &cell) const override; diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc index 5433bf7471..3e10ef99e2 100644 --- a/source/fe/mapping.cc +++ b/source/fe/mapping.cc @@ -24,15 +24,17 @@ DEAL_II_NAMESPACE_OPEN template -std::array, GeometryInfo::vertices_per_cell> +boost::container::small_vector, + GeometryInfo::vertices_per_cell> Mapping::get_vertices( const typename Triangulation::cell_iterator &cell) const { - std::array, GeometryInfo::vertices_per_cell> vertices; - for (const unsigned int i : GeometryInfo::vertex_indices()) - { - vertices[i] = cell->vertex(i); - } + boost::container::small_vector, + GeometryInfo::vertices_per_cell> + vertices; + for (const unsigned int i : cell->vertex_indices()) + vertices.push_back(cell->vertex(i)); + return vertices; } diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index 163719f12e..50667b2f51 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -357,7 +357,8 @@ MappingFEField:: template -std::array, GeometryInfo::vertices_per_cell> +boost::container::small_vector, + GeometryInfo::vertices_per_cell> MappingFEField::get_vertices( const typename Triangulation::cell_iterator &cell) const { @@ -395,7 +396,9 @@ MappingFEField::get_vertices( const VectorType &vector = uses_level_dofs ? *euler_vector[cell->level()] : *euler_vector[0]; - std::array, GeometryInfo::vertices_per_cell> vertices; + boost::container::small_vector, + GeometryInfo::vertices_per_cell> + vertices(GeometryInfo::vertices_per_cell); for (unsigned int i = 0; i < dofs_per_cell; ++i) { const unsigned int comp = fe_to_real diff --git a/source/fe/mapping_q1_eulerian.cc b/source/fe/mapping_q1_eulerian.cc index 36ff806f39..3d495c32b4 100644 --- a/source/fe/mapping_q1_eulerian.cc +++ b/source/fe/mapping_q1_eulerian.cc @@ -50,11 +50,14 @@ MappingQ1Eulerian::MappingQ1Eulerian( template -std::array, GeometryInfo::vertices_per_cell> +boost::container::small_vector, + GeometryInfo::vertices_per_cell> MappingQ1Eulerian::get_vertices( const typename Triangulation::cell_iterator &cell) const { - std::array, GeometryInfo::vertices_per_cell> vertices; + boost::container::small_vector, + GeometryInfo::vertices_per_cell> + vertices(GeometryInfo::vertices_per_cell); // The assertions can not be in the constructor, since this would // require to call dof_handler.distribute_dofs(fe) *before* the mapping // object is constructed, which is not necessarily what we want. @@ -105,8 +108,7 @@ std::vector> MappingQ1Eulerian::compute_mapping_support_points( const typename Triangulation::cell_iterator &cell) const { - const std::array, GeometryInfo::vertices_per_cell> - vertices = this->get_vertices(cell); + const auto vertices = this->get_vertices(cell); std::vector> a(GeometryInfo::vertices_per_cell); for (const unsigned int i : GeometryInfo::vertex_indices()) diff --git a/source/fe/mapping_q_eulerian.cc b/source/fe/mapping_q_eulerian.cc index e10223f422..f44f0d796e 100644 --- a/source/fe/mapping_q_eulerian.cc +++ b/source/fe/mapping_q_eulerian.cc @@ -121,7 +121,8 @@ MappingQEulerian::MappingQEulerianGeneric:: // .... COMPUTE MAPPING SUPPORT POINTS template -std::array, GeometryInfo::vertices_per_cell> +boost::container::small_vector, + GeometryInfo::vertices_per_cell> MappingQEulerian::get_vertices( const typename Triangulation::cell_iterator &cell) const { @@ -130,11 +131,10 @@ MappingQEulerian::get_vertices( dynamic_cast(*this->qp_mapping) .compute_mapping_support_points(cell); - std::array, GeometryInfo::vertices_per_cell> - vertex_locations; - std::copy(a.begin(), - a.begin() + GeometryInfo::vertices_per_cell, - vertex_locations.begin()); + boost::container::small_vector, + GeometryInfo::vertices_per_cell> + vertex_locations(a.begin(), + a.begin() + GeometryInfo::vertices_per_cell); return vertex_locations; } @@ -142,7 +142,8 @@ MappingQEulerian::get_vertices( template -std::array, GeometryInfo::vertices_per_cell> +boost::container::small_vector, + GeometryInfo::vertices_per_cell> MappingQEulerian::MappingQEulerianGeneric:: get_vertices( const typename Triangulation::cell_iterator &cell) const diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 3d878bbc04..3497a2ed00 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -2518,8 +2518,13 @@ MappingQGeneric::transform_real_to_unit_cell( // cell and that value is returned. // * In 3D there is no (known to the authors) exact formula, so the Newton // algorithm is used. - const std::array, GeometryInfo::vertices_per_cell> - vertices = this->get_vertices(cell); + const auto vertices_ = this->get_vertices(cell); + + std::array, GeometryInfo::vertices_per_cell> + vertices; + for (unsigned int i = 0; i < vertices.size(); ++i) + vertices[i] = vertices_[i]; + try { switch (dim) diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 62212f8881..4b1aec7c5d 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2988,15 +2988,29 @@ namespace GridTools diameter(const typename Triangulation::cell_iterator &cell, const Mapping &mapping) { + // see also TriaAccessor::diameter() + const auto vertices = mapping.get_vertices(cell); - switch (dim) + switch (cell->reference_cell_type()) { - case 1: + case ReferenceCell::Type::Line: return (vertices[1] - vertices[0]).norm(); - case 2: + case ReferenceCell::Type::Tri: + return std::max(std::max((vertices[1] - vertices[0]).norm(), + (vertices[2] - vertices[1]).norm()), + (vertices[2] - vertices[0]).norm()); + case ReferenceCell::Type::Quad: return std::max((vertices[3] - vertices[0]).norm(), (vertices[2] - vertices[1]).norm()); - case 3: + case ReferenceCell::Type::Tet: + return std::max( + std::max(std::max((vertices[1] - vertices[0]).norm(), + (vertices[2] - vertices[0]).norm()), + std::max((vertices[2] - vertices[1]).norm(), + (vertices[3] - vertices[0]).norm())), + std::max((vertices[3] - vertices[1]).norm(), + (vertices[3] - vertices[2]).norm())); + case ReferenceCell::Type::Hex: return std::max(std::max((vertices[7] - vertices[0]).norm(), (vertices[6] - vertices[1]).norm()), std::max((vertices[2] - vertices[5]).norm(), diff --git a/source/opencascade/utilities.cc b/source/opencascade/utilities.cc index b9f11d79b9..1011d14ac6 100644 --- a/source/opencascade/utilities.cc +++ b/source/opencascade/utilities.cc @@ -619,8 +619,7 @@ namespace OpenCASCADE visited_faces[face_index] = false; // extract mapped vertex locations - std::array, GeometryInfo<2>::vertices_per_cell> - verts = mapping.get_vertices(cell); + const auto verts = mapping.get_vertices(cell); vert_to_point[v0] = verts[GeometryInfo<2>::face_to_cell_vertices( f, 0, true, false, false)]; vert_to_point[v1] = verts[GeometryInfo<2>::face_to_cell_vertices( diff --git a/tests/simplex/cell_measure_01.cc b/tests/simplex/cell_measure_01.cc index 08a57fd512..a5112d7330 100644 --- a/tests/simplex/cell_measure_01.cc +++ b/tests/simplex/cell_measure_01.cc @@ -17,8 +17,13 @@ // Test TriaAccessor::measure() and TriaAccessor::diameter(). +#include + +#include #include +#include + #include "../tests.h" template @@ -35,6 +40,13 @@ process(const std::vector> &vertices, deallog << "measure: " << cell->measure() << std::endl; deallog << "diameter: " << cell->diameter() << std::endl; } + + const MappingFE mapping(Simplex::FE_P(1)); + + deallog << "diameter_min: " << GridTools::minimal_cell_diameter(tria, mapping) + << std::endl; + deallog << "diameter_max: " << GridTools::maximal_cell_diameter(tria, mapping) + << std::endl; deallog << std::endl; } diff --git a/tests/simplex/cell_measure_01.with_simplex_support=on.output b/tests/simplex/cell_measure_01.with_simplex_support=on.output index a8127115b1..8fe42d75a2 100644 --- a/tests/simplex/cell_measure_01.with_simplex_support=on.output +++ b/tests/simplex/cell_measure_01.with_simplex_support=on.output @@ -2,8 +2,12 @@ DEAL::dim=2 spacedim=2: DEAL::measure: 0.500000 DEAL::diameter: 1.41421 +DEAL::diameter_min: 1.41421 +DEAL::diameter_max: 1.41421 DEAL:: DEAL::dim=3 spacedim=3: DEAL::measure: 0.166667 DEAL::diameter: 1.41421 +DEAL::diameter_min: 1.41421 +DEAL::diameter_max: 1.41421 DEAL:: -- 2.39.5