From 57b6d20fbe6b1cb25921485a72a52724fe16a40c Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 20 Jul 2020 12:30:26 +0200 Subject: [PATCH] Implement TriaAccessor::measure() and TriaAccessor::diameter() for simplex --- .../deal.II/grid/tria_accessor.templates.h | 20 +++- source/grid/grid_tools_nontemplates.cc | 24 +++++ source/grid/tria_accessor.cc | 6 +- tests/simplex/cell_measure_01.cc | 96 +++++++++++++++++++ ..._measure_01.with_simplex_support=on.output | 9 ++ 5 files changed, 149 insertions(+), 6 deletions(-) create mode 100644 tests/simplex/cell_measure_01.cc create mode 100644 tests/simplex/cell_measure_01.with_simplex_support=on.output diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h index 1ee1438377..a8e940da23 100644 --- a/include/deal.II/grid/tria_accessor.templates.h +++ b/include/deal.II/grid/tria_accessor.templates.h @@ -1885,14 +1885,26 @@ template double TriaAccessor::diameter() const { - switch (structdim) + switch (this->reference_cell_type()) { - case 1: + case ReferenceCell::Type::Line: return (this->vertex(1) - this->vertex(0)).norm(); - case 2: + case ReferenceCell::Type::Tri: + return std::max(std::max((this->vertex(1) - this->vertex(0)).norm(), + (this->vertex(2) - this->vertex(1)).norm()), + (this->vertex(2) - this->vertex(0)).norm()); + case ReferenceCell::Type::Quad: return std::max((this->vertex(3) - this->vertex(0)).norm(), (this->vertex(2) - this->vertex(1)).norm()); - case 3: + case ReferenceCell::Type::Tet: + return std::max( + std::max(std::max((this->vertex(1) - this->vertex(0)).norm(), + (this->vertex(2) - this->vertex(0)).norm()), + std::max((this->vertex(2) - this->vertex(1)).norm(), + (this->vertex(3) - this->vertex(0)).norm())), + std::max((this->vertex(3) - this->vertex(1)).norm(), + (this->vertex(3) - this->vertex(2)).norm())); + case ReferenceCell::Type::Hex: return std::max(std::max((this->vertex(7) - this->vertex(0)).norm(), (this->vertex(6) - this->vertex(1)).norm()), std::max((this->vertex(2) - this->vertex(5)).norm(), diff --git a/source/grid/grid_tools_nontemplates.cc b/source/grid/grid_tools_nontemplates.cc index 126ea9ec8a..1bdc089ebe 100644 --- a/source/grid/grid_tools_nontemplates.cc +++ b/source/grid/grid_tools_nontemplates.cc @@ -45,6 +45,20 @@ namespace GridTools cell_measure<2>(const std::vector> & all_vertices, const ArrayView &vertex_indices) { + if (vertex_indices.size() == 3) // triangle + { + const double x[3] = {all_vertices[vertex_indices[0]](0), + all_vertices[vertex_indices[1]](0), + all_vertices[vertex_indices[2]](0)}; + + const double y[3] = {all_vertices[vertex_indices[0]](1), + all_vertices[vertex_indices[1]](1), + all_vertices[vertex_indices[2]](1)}; + + return 0.5 * std::abs((x[0] - x[2]) * (y[1] - y[0]) - + (x[0] - x[1]) * (y[2] - y[0])); + } + AssertDimension(vertex_indices.size(), GeometryInfo<2>::vertices_per_cell); /* @@ -105,6 +119,16 @@ namespace GridTools cell_measure<3>(const std::vector> & all_vertices, const ArrayView &vertex_indices) { + if (vertex_indices.size() == 4) // tetrahedron + { + const auto &a = all_vertices[vertex_indices[0]]; + const auto &b = all_vertices[vertex_indices[1]]; + const auto &c = all_vertices[vertex_indices[2]]; + const auto &d = all_vertices[vertex_indices[3]]; + + return (1.0 / 6.0) * std::abs((a - d) * cross_product_3d(b - d, c - d)); + } + AssertDimension(vertex_indices.size(), GeometryInfo<3>::vertices_per_cell); // note that this is the // cell_measure based on the new diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index a693d61ccb..d14bfcdb6e 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1269,7 +1269,8 @@ namespace vertex_indices[i] = accessor.vertex_index(i); return GridTools::cell_measure<2>( - accessor.get_triangulation().get_vertices(), vertex_indices); + accessor.get_triangulation().get_vertices(), + ArrayView(vertex_indices, accessor.n_vertices())); } @@ -1281,7 +1282,8 @@ namespace vertex_indices[i] = accessor.vertex_index(i); return GridTools::cell_measure<3>( - accessor.get_triangulation().get_vertices(), vertex_indices); + accessor.get_triangulation().get_vertices(), + ArrayView(vertex_indices, accessor.n_vertices())); } diff --git a/tests/simplex/cell_measure_01.cc b/tests/simplex/cell_measure_01.cc new file mode 100644 index 0000000000..08a57fd512 --- /dev/null +++ b/tests/simplex/cell_measure_01.cc @@ -0,0 +1,96 @@ +// --------------------------------------------------------------------- +// +// 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 TriaAccessor::measure() and TriaAccessor::diameter(). + +#include + +#include "../tests.h" + +template +void +process(const std::vector> &vertices, + const std::vector> & cells) +{ + Triangulation tria; + tria.create_triangulation(vertices, cells, SubCellData()); + + deallog << "dim=" << dim << " spacedim=" << spacedim << ":" << std::endl; + for (const auto &cell : tria.active_cell_iterators()) + { + deallog << "measure: " << cell->measure() << std::endl; + deallog << "diameter: " << cell->diameter() << std::endl; + } + deallog << std::endl; +} + +template +void +test() +{ + Assert(false, ExcNotImplemented()); +} + +template <> +void +test<2>() +{ + const int dim = 2; + const int spacedim = 2; + + std::vector> vertices; + vertices.emplace_back(0, 0); + vertices.emplace_back(1, 0); + vertices.emplace_back(0, 1); + + std::vector> cells; + CellData cell; + cell.vertices = {0, 1, 2}; + cells.push_back(cell); + + process(vertices, cells); +} + +template <> +void +test<3>() +{ + const int dim = 3; + const int spacedim = 3; + + std::vector> vertices; + vertices.emplace_back(0, 0, 0); + vertices.emplace_back(1, 0, 0); + vertices.emplace_back(0, 1, 0); + vertices.emplace_back(0, 0, 1); + + std::vector> cells; + CellData cell; + cell.vertices = {0, 1, 2, 3}; + cells.push_back(cell); + + process(vertices, cells); +} + +int +main() +{ + initlog(); + + test<2>(); + test<3>(); +} diff --git a/tests/simplex/cell_measure_01.with_simplex_support=on.output b/tests/simplex/cell_measure_01.with_simplex_support=on.output new file mode 100644 index 0000000000..a8127115b1 --- /dev/null +++ b/tests/simplex/cell_measure_01.with_simplex_support=on.output @@ -0,0 +1,9 @@ + +DEAL::dim=2 spacedim=2: +DEAL::measure: 0.500000 +DEAL::diameter: 1.41421 +DEAL:: +DEAL::dim=3 spacedim=3: +DEAL::measure: 0.166667 +DEAL::diameter: 1.41421 +DEAL:: -- 2.39.5