From: David Wells Date: Thu, 9 Feb 2023 18:21:15 +0000 (-0500) Subject: Make GridTools::volume() work with simplices. X-Git-Tag: v9.5.0-rc1~555^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14780%2Fhead;p=dealii.git Make GridTools::volume() work with simplices. --- diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index b5344d499f..e706ce5566 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -146,7 +146,34 @@ namespace GridTools * Compute the volume (i.e. the dim-dimensional measure) of the * triangulation. We compute the measure using the integral $\sum_K \int_K 1 * \; dx$ where $K$ are the cells of the given triangulation. The integral - * is approximated via quadrature for which we need the mapping argument. + * is approximated via quadrature. This version of the function uses a + * linear mapping to compute the JxW values on each cell. + * + * If the triangulation is a dim-dimensional one embedded in a higher + * dimensional space of dimension spacedim, then the value returned is the + * dim-dimensional measure. For example, for a two-dimensional triangulation + * in three-dimensional space, the value returned is the area of the surface + * so described. (This obviously makes sense since the spacedim-dimensional + * measure of a dim-dimensional triangulation would always be zero if dim @< + * spacedim). + * + * This function also works for objects of type + * parallel::distributed::Triangulation, in which case the function is a + * collective operation. + * + * @param tria The triangulation. + * @return The dim-dimensional measure of the domain described by the + * triangulation, as discussed above. + */ + template + double + volume(const Triangulation &tria); + + /** + * Compute the volume (i.e. the dim-dimensional measure) of the + * triangulation. We compute the measure using the integral $\sum_K \int_K 1 + * \; dx$ where $K$ are the cells of the given triangulation. The integral + * is approximated via quadrature for which we use the mapping argument. * * If the triangulation is a dim-dimensional one embedded in a higher * dimensional space of dimension spacedim, then the value returned is the @@ -161,24 +188,18 @@ namespace GridTools * collective operation. * * @param tria The triangulation. - * @param mapping An optional argument used to denote the mapping that - * should be used when describing whether cells are bounded by straight or - * curved faces. The default is to use a $Q_1$ mapping, which corresponds to - * straight lines bounding the cells. + * @param mapping The Mapping which computes the Jacobians used to + * approximate the volume via quadrature. Explicitly using a higher-order + * Mapping (i.e., instead of using the other version of this function) will + * result in a more accurate approximation of the volume on Triangulations + * with curvature described by Manifold objects. * @return The dim-dimensional measure of the domain described by the * triangulation, as discussed above. */ template double volume(const Triangulation &tria, - const Mapping & mapping = - (ReferenceCells::get_hypercube() -#ifndef _MSC_VER - .template get_default_linear_mapping() -#else - .ReferenceCell::get_default_linear_mapping() -#endif - )); + const Mapping & mapping); /** * Return an approximation of the diameter of the smallest active cell of a diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index f5d400f7b3..a558986d80 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -32,6 +32,7 @@ #include #include #include +#include #include #include @@ -135,6 +136,20 @@ namespace GridTools + template + double + volume(const Triangulation &triangulation) + { + Assert(triangulation.get_reference_cells().size() == 1, + ExcNotImplemented()); + const ReferenceCell reference_cell = triangulation.get_reference_cells()[0]; + return volume( + triangulation, + reference_cell.template get_default_linear_mapping()); + } + + + template double volume(const Triangulation &triangulation, @@ -145,30 +160,31 @@ namespace GridTools if (const auto *p = dynamic_cast *>(&mapping)) mapping_degree = p->get_degree(); else if (const auto *p = - dynamic_cast *>(&mapping)) + dynamic_cast *>(&mapping)) mapping_degree = p->get_degree(); // then initialize an appropriate quadrature formula - const QGauss quadrature_formula(mapping_degree + 1); + Assert(triangulation.get_reference_cells().size() == 1, + ExcNotImplemented()); + const ReferenceCell reference_cell = triangulation.get_reference_cells()[0]; + const Quadrature quadrature_formula = + reference_cell.template get_gauss_type_quadrature(mapping_degree + + 1); const unsigned int n_q_points = quadrature_formula.size(); // we really want the JxW values from the FEValues object, but it // wants a finite element. create a cheap element as a dummy // element - FE_Nothing dummy_fe; + FE_Nothing dummy_fe(reference_cell); FEValues fe_values(mapping, dummy_fe, quadrature_formula, update_JxW_values); - typename Triangulation::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); - double local_volume = 0; // compute the integral quantities by quadrature - for (; cell != endc; ++cell) + for (const auto &cell : triangulation.active_cell_iterators()) if (cell->is_locally_owned()) { fe_values.reinit(cell); diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 052630b44b..4d22dc0060 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -219,6 +219,9 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) diameter( const Triangulation &); + template double + volume(const Triangulation &); + template double volume(const Triangulation &, const Mapping &); diff --git a/tests/grid/grid_tools_03.cc b/tests/grid/grid_tools_03.cc index cb1aeeb6b2..96159b55f8 100644 --- a/tests/grid/grid_tools_03.cc +++ b/tests/grid/grid_tools_03.cc @@ -20,7 +20,6 @@ #include #include #include -#include #include #include "../tests.h" @@ -42,8 +41,13 @@ test1() deallog << dim << "d, " << "hypercube volume, " << i * 2 << " refinements: " << GridTools::volume(tria) << std::endl; - }; - }; + } + + Triangulation simplex_tria; + GridGenerator::convert_hypercube_to_simplex_mesh(tria, simplex_tria); + deallog << dim << "d, simplex hypercube volume: " + << GridTools::volume(simplex_tria) << std::endl; + } // test 2: hyperball if (dim >= 2) @@ -51,9 +55,6 @@ test1() Triangulation tria; GridGenerator::hyper_ball(tria, Point(), 1); - static const SphericalManifold boundary; - tria.set_manifold(0, boundary); - for (unsigned int i = 0; i < 4; ++i) { tria.refine_global(1); @@ -61,6 +62,13 @@ test1() << "hyperball volume, " << i << " refinements: " << GridTools::volume(tria) << std::endl; } + + Triangulation simplex_tria; + GridGenerator::convert_hypercube_to_simplex_mesh(tria, simplex_tria); + simplex_tria.set_all_manifold_ids(numbers::flat_manifold_id); + deallog << dim << "d, simplex hyperball volume: " + << GridTools::volume(simplex_tria) << std::endl; + deallog << "exact value=" << (dim == 2 ? numbers::PI : 4. / 3. * numbers::PI) << std::endl; } diff --git a/tests/grid/grid_tools_03.output b/tests/grid/grid_tools_03.output index f785507852..27857a71b9 100644 --- a/tests/grid/grid_tools_03.output +++ b/tests/grid/grid_tools_03.output @@ -1,17 +1,22 @@ DEAL::1d, hypercube volume, 0 refinements: 1.000 DEAL::1d, hypercube volume, 2 refinements: 1.000 +DEAL::1d, simplex hypercube volume: 1.000 DEAL::2d, hypercube volume, 0 refinements: 1.000 DEAL::2d, hypercube volume, 2 refinements: 1.000 +DEAL::2d, simplex hypercube volume: 1.000 DEAL::2d, hyperball volume, 0 refinements: 2.828 DEAL::2d, hyperball volume, 1 refinements: 3.061 DEAL::2d, hyperball volume, 2 refinements: 3.121 DEAL::2d, hyperball volume, 3 refinements: 3.137 +DEAL::2d, simplex hyperball volume: 3.137 DEAL::exact value=3.142 DEAL::3d, hypercube volume, 0 refinements: 1.000 DEAL::3d, hypercube volume, 2 refinements: 1.000 +DEAL::3d, simplex hypercube volume: 1.000 DEAL::3d, hyperball volume, 0 refinements: 3.210 DEAL::3d, hyperball volume, 1 refinements: 3.917 DEAL::3d, hyperball volume, 2 refinements: 4.119 DEAL::3d, hyperball volume, 3 refinements: 4.171 +DEAL::3d, simplex hyperball volume: 4.171 DEAL::exact value=4.189