From deb91748d1d3cd31f00d12436c0cca2d5eeab1b4 Mon Sep 17 00:00:00 2001 From: David Wells Date: Fri, 20 May 2022 17:02:13 -0400 Subject: [PATCH] Implement pyramid volume and fix wedge volume. --- source/grid/grid_tools_nontemplates.cc | 254 +++++++++++++------------ tests/simplex/cell_measure_01.cc | 124 ++++++++++-- tests/simplex/cell_measure_01.output | 68 +++++-- 3 files changed, 293 insertions(+), 153 deletions(-) diff --git a/source/grid/grid_tools_nontemplates.cc b/source/grid/grid_tools_nontemplates.cc index 03e867d3cc..8a32ccb5d0 100644 --- a/source/grid/grid_tools_nontemplates.cc +++ b/source/grid/grid_tools_nontemplates.cc @@ -130,126 +130,146 @@ namespace GridTools return (1.0 / 6.0) * std::abs((a - d) * cross_product_3d(b - d, c - d)); } + else if (vertex_indices.size() == 5) // pyramid + { + // This remarkably simple formula comes from Equation 4 of + // "Calculation of the volume of a general hexahedron for flow + // predictions", Davies and Salmond, AIAA Journal vol. 23 no. 6. + const auto &x0 = all_vertices[vertex_indices[0]]; + const auto &x1 = all_vertices[vertex_indices[1]]; + const auto &x2 = all_vertices[vertex_indices[2]]; + const auto &x3 = all_vertices[vertex_indices[3]]; + const auto &x4 = all_vertices[vertex_indices[4]]; + + const auto v01 = x1 - x0; + const auto v02 = x2 - x0; + const auto v03 = x3 - x0; + const auto v04 = x4 - x0; + const auto v21 = x2 - x1; + + // doing high - low consistently puts us off by -1 from the original + // paper in the first term + return -v04 * cross_product_3d(v21, v03) / 6.0 + + v03 * cross_product_3d(v01, v02) / 12.0; + } else if (vertex_indices.size() == 6) // wedge { - /* - The following python/sympy script was used: - - #!/usr/bin/env python - # coding: utf-8 - import sympy as sp - from sympy.simplify.cse_main import cse - xs = list(sp.symbols(" ".join(["x{}".format(i) for i in range(6)]))) - ys = list(sp.symbols(" ".join(["y{}".format(i) for i in range(6)]))) - zs = list(sp.symbols(" ".join(["z{}".format(i) for i in range(6)]))) - xi, eta, zeta = sp.symbols("xi eta zeta") - tphi = [(1 - xi - eta)*(1 - zeta), - (xi)*(1 - zeta), - (eta)*(1 - zeta), - (1 - xi - zeta)*(zeta), - (xi)*(zeta), - (eta)*(zeta)] - x_real = sum(xs[i]*tphi[i] for i in range(len(xs))) - y_real = sum(ys[i]*tphi[i] for i in range(len(xs))) - z_real = sum(zs[i]*tphi[i] for i in range(len(xs))) - J = sp.Matrix([[var.diff(v) for v in [xi, eta, zeta]] - for var in [x_real, y_real, z_real]]) - detJ = J.det() - detJ2 = detJ.expand().collect(zeta).collect(eta).collect(xi) - for x in xs: - detJ2 = detJ2.collect(x) - for y in ys: - detJ2 = detJ2.collect(y) - for z in zs: - detJ2 = detJ2.collect(z) - measure = sp.integrate(sp.integrate(sp.integrate(detJ2, (xi, 0, 1)), - (eta, 0, 1)), (zeta, 0, 1)) - measure2 = measure - for vs in [xs, ys, zs]: - for v in vs: - measure2 = measure2.collect(v) - pairs, expression = cse(measure2) - for pair in pairs: - print("const double " + sp.ccode(pair[0]) + " = " - + sp.ccode(pair[1]) + ";") - print("const double result = " + sp.ccode(expression[0]) + ";") - print("return result;") + /* Script used to generate volume code: + + #!/usr/bin/env python + # coding: utf-8 + import sympy as sp + from sympy.simplify.cse_main import cse + n_vertices = 6 + xs = list(sp.symbols(" ".join(["x{}".format(i) + for i in range(n_vertices)]))) + ys = list(sp.symbols(" ".join(["y{}".format(i) + for i in range(n_vertices)]))) + zs = list(sp.symbols(" ".join(["z{}".format(i) + for i in range(n_vertices)]))) + xi, eta, zeta = sp.symbols("xi eta zeta") + tphi = [(1 - xi - eta)*(1 - zeta), + (xi)*(1 - zeta), + (eta)*(1 - zeta), + (1 - xi - eta)*(zeta), + (xi)*(zeta), + (eta)*(zeta)] + x_real = sum(xs[i]*tphi[i] for i in range(n_vertices)) + y_real = sum(ys[i]*tphi[i] for i in range(n_vertices)) + z_real = sum(zs[i]*tphi[i] for i in range(n_vertices)) + J = sp.Matrix([[var.diff(v) for v in [xi, eta, zeta]] + for var in [x_real, y_real, z_real]]) + detJ = J.det() + detJ2 = detJ.expand().collect(zeta).collect(eta).collect(xi) + for x in xs: + detJ2 = detJ2.collect(x) + for y in ys: + detJ2 = detJ2.collect(y) + for z in zs: + detJ2 = detJ2.collect(z) + measure = sp.integrate(sp.integrate( + sp.integrate(detJ2, (eta, 0, 1 - xi)), + (xi, 0, 1)), (zeta, 0, 1)) + measure2 = measure + for vs in [xs, ys, zs]: + for v in vs: + measure2 = measure2.collect(v) + + pairs, expression = cse(measure2) + for vertex_no in range(n_vertices): + for (coordinate, index) in [('x', 0), ('y', 1), ('z', 2)]: + print( + "const double {}{} = all_vertices[vertex_indices[{}]][{}];" + .format(coordinate, vertex_no, vertex_no, index)) + + for pair in pairs: + print("const double " + sp.ccode(pair[0]) + " = " + + sp.ccode(pair[1]) + ";") + print("const double result = " + sp.ccode(expression[0]) + ";") + print("return result;") */ - const double x0 = all_vertices[vertex_indices[0]](0); - const double y0 = all_vertices[vertex_indices[0]](1); - const double z0 = all_vertices[vertex_indices[0]](2); - const double x1 = all_vertices[vertex_indices[1]](0); - const double y1 = all_vertices[vertex_indices[1]](1); - const double z1 = all_vertices[vertex_indices[1]](2); - const double x2 = all_vertices[vertex_indices[2]](0); - const double y2 = all_vertices[vertex_indices[2]](1); - const double z2 = all_vertices[vertex_indices[2]](2); - const double x3 = all_vertices[vertex_indices[3]](0); - const double y3 = all_vertices[vertex_indices[3]](1); - const double z3 = all_vertices[vertex_indices[3]](2); - const double x4 = all_vertices[vertex_indices[4]](0); - const double y4 = all_vertices[vertex_indices[4]](1); - const double z4 = all_vertices[vertex_indices[4]](2); - const double x5 = all_vertices[vertex_indices[5]](0); - const double y5 = all_vertices[vertex_indices[5]](1); - const double z5 = all_vertices[vertex_indices[5]](2); - - const double x6 = (1.0 / 6.0) * z5; - const double x7 = (1.0 / 12.0) * z1; - const double x8 = -x7; - const double x9 = (1.0 / 12.0) * z2; - const double x10 = -x9; - const double x11 = (1.0 / 4.0) * z5; - const double x12 = -x11; - const double x13 = (1.0 / 12.0) * z0; - const double x14 = x12 + x13; - const double x15 = (1.0 / 4.0) * z2; - const double x16 = (1.0 / 6.0) * z4; - const double x17 = (1.0 / 4.0) * z1; - const double x18 = (1.0 / 6.0) * z0; - const double x19 = x17 - x18; - const double x20 = -x6; - const double x21 = (1.0 / 4.0) * z0; - const double x22 = -x21; - const double x23 = -x17; - const double x24 = -x15; - const double x25 = (1.0 / 6.0) * z3; - const double x26 = x24 - x25; - const double x27 = x18 + x23; - const double x28 = (1.0 / 3.0) * z2; - const double x29 = (1.0 / 12.0) * z5; - const double x30 = (1.0 / 12.0) * z3; - const double x31 = -x30; - const double x32 = (1.0 / 4.0) * z4; - const double x33 = x31 + x32; - const double x34 = (1.0 / 3.0) * z1; - const double x35 = (1.0 / 12.0) * z4; - const double x36 = -x16; - const double x37 = x15 + x25; - const double x38 = -x13; - const double x39 = x11 + x38; - const double x40 = -x32; - const double x41 = x30 + x40; - const double x42 = (1.0 / 3.0) * z0; - const double x43 = (1.0 / 4.0) * z3; - const double x44 = x32 - x43; - const double x45 = x40 + x43; - return x0 * (y1 * (-x28 + x29 + x33) + y2 * (x12 + x31 + x34 - x35) + - y3 * (x20 + x7 + x9) + y4 * (x23 + x6 + x9) + - y5 * (x36 + x37 + x8)) + - x1 * (y0 * (x28 - x29 + x41) + y2 * (x11 + x33 - x42) + - y3 * (x39 + x9) + y4 * (x12 + x21 + x24) + - y5 * (x13 + x24 + x44)) + - x2 * (y0 * (x11 + x30 - x34 + x35) + y1 * (x12 + x41 + x42) + - y3 * (x39 + x8) + y4 * (x12 + x17 + x38) + - y5 * (x17 + x22 + x44)) + - x3 * (-x6 * y4 + y0 * (x10 + x6 + x8) + y1 * (x10 + x14) + - y2 * (x14 + x7) + y5 * (x15 + x16 + x19)) + - x4 * (x6 * y3 + y0 * (x10 + x17 + x20) + y1 * (x11 + x15 + x22) + - y2 * (x11 + x13 + x23) + y5 * (x26 + x27)) + - x5 * (y0 * (x16 + x26 + x7) + y1 * (x15 + x38 + x45) + - y2 * (x21 + x23 + x45) + y3 * (x24 + x27 + x36) + - y4 * (x19 + x37)); + const double x0 = all_vertices[vertex_indices[0]][0]; + const double y0 = all_vertices[vertex_indices[0]][1]; + const double z0 = all_vertices[vertex_indices[0]][2]; + const double x1 = all_vertices[vertex_indices[1]][0]; + const double y1 = all_vertices[vertex_indices[1]][1]; + const double z1 = all_vertices[vertex_indices[1]][2]; + const double x2 = all_vertices[vertex_indices[2]][0]; + const double y2 = all_vertices[vertex_indices[2]][1]; + const double z2 = all_vertices[vertex_indices[2]][2]; + const double x3 = all_vertices[vertex_indices[3]][0]; + const double y3 = all_vertices[vertex_indices[3]][1]; + const double z3 = all_vertices[vertex_indices[3]][2]; + const double x4 = all_vertices[vertex_indices[4]][0]; + const double y4 = all_vertices[vertex_indices[4]][1]; + const double z4 = all_vertices[vertex_indices[4]][2]; + const double x5 = all_vertices[vertex_indices[5]][0]; + const double y5 = all_vertices[vertex_indices[5]][1]; + const double z5 = all_vertices[vertex_indices[5]][2]; + const double x6 = (1.0 / 12.0) * z1; + const double x7 = -x6; + const double x8 = (1.0 / 12.0) * z3; + const double x9 = x7 + x8; + const double x10 = (1.0 / 12.0) * z2; + const double x11 = -x8; + const double x12 = x10 + x11; + const double x13 = (1.0 / 6.0) * z2; + const double x14 = (1.0 / 12.0) * z4; + const double x15 = (1.0 / 6.0) * z1; + const double x16 = (1.0 / 12.0) * z5; + const double x17 = -x16; + const double x18 = x16 + x7; + const double x19 = -x14; + const double x20 = x10 + x19; + const double x21 = (1.0 / 12.0) * z0; + const double x22 = x19 + x21; + const double x23 = -x10; + const double x24 = x14 + x23; + const double x25 = (1.0 / 6.0) * z0; + const double x26 = x17 + x21; + const double x27 = x23 + x8; + const double x28 = -x21; + const double x29 = x16 + x28; + const double x30 = x17 + x6; + const double x31 = x14 + x28; + const double x32 = x11 + x6; + const double x33 = (1.0 / 6.0) * z5; + const double x34 = (1.0 / 6.0) * z4; + const double x35 = (1.0 / 6.0) * z3; + const double result = + x0 * (x12 * y5 + x9 * y4 + y1 * (-x13 + x14 + x8) + + y2 * (x11 + x15 + x17) + y3 * (x18 + x20)) + + x1 * (x22 * y3 + x24 * y5 + y0 * (x11 + x13 + x19) + + y2 * (x14 + x16 - x25) + y4 * (x26 + x27)) + + x2 * (x29 * y3 + x30 * y4 + y0 * (-x15 + x16 + x8) + + y1 * (x17 + x19 + x25) + y5 * (x31 + x32)) + + x3 * (x26 * y2 + x31 * y1 + y0 * (x24 + x30) + y4 * (x28 + x33 + x7) + + y5 * (x10 + x21 - x34)) + + x4 * (x18 * y2 + x32 * y0 + y1 * (x12 + x29) + y3 * (x21 - x33 + x6) + + y5 * (x23 + x35 + x7)) + + x5 * (x20 * y1 + x27 * y0 + y2 * (x22 + x9) + y3 * (x23 + x28 + x34) + + y4 * (x10 - x35 + x6)); + return result; } AssertDimension(vertex_indices.size(), GeometryInfo<3>::vertices_per_cell); diff --git a/tests/simplex/cell_measure_01.cc b/tests/simplex/cell_measure_01.cc index 8c601aa0f5..59525e0a37 100644 --- a/tests/simplex/cell_measure_01.cc +++ b/tests/simplex/cell_measure_01.cc @@ -18,9 +18,10 @@ // Test TriaAccessor::measure(), TriaAccessor::diameter(), and // TriaAccessor::barycenter(). +#include #include #include -#include +#include #include #include @@ -37,33 +38,44 @@ process(const std::vector> &vertices, 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; - // Only implemented for tensor-product and simplices at the moment - if (cell->reference_cell() == ReferenceCells::get_simplex()) - deallog << "barycenter: " << cell->barycenter() << std::endl; - } - std::shared_ptr> mapping; - - const auto reference_cells = tria.get_reference_cells(); - + const auto reference_cells = tria.get_reference_cells(); AssertDimension(reference_cells.size(), 1); if (reference_cells[0] == ReferenceCells::get_simplex()) mapping = std::make_shared>(FE_SimplexP(1)); else if (reference_cells[0] == ReferenceCells::Wedge) mapping = std::make_shared>(FE_WedgeP(1)); + else if (reference_cells[0] == ReferenceCells::Pyramid) + mapping = std::make_shared>(FE_PyramidP(1)); else AssertThrow(false, ExcNotImplemented()); - deallog << "diameter_min: " - << GridTools::minimal_cell_diameter(tria, *mapping) << std::endl; - deallog << "diameter_max: " - << GridTools::maximal_cell_diameter(tria, *mapping) << std::endl; + const Quadrature quad = + reference_cells.front().template get_gauss_type_quadrature(2); + FE_Nothing fe(reference_cells.front()); + FEValues fe_values(*mapping, fe, quad, update_JxW_values); + + deallog << "dim=" << dim << " spacedim=" << spacedim << ':' << std::endl; + for (const auto &cell : tria.active_cell_iterators()) + { + fe_values.reinit(cell); + const double measure = std::accumulate(fe_values.get_JxW_values().begin(), + fe_values.get_JxW_values().end(), + 0.0); + + deallog << "measure: " << cell->measure() << std::endl; + deallog << "measure (via quadrature): " << measure << std::endl; + deallog << "diameter " << cell->diameter() << std::endl; + // Only implemented for tensor-product and simplices at the moment + if (cell->reference_cell() == ReferenceCells::get_simplex()) + deallog << "barycenter " << cell->barycenter() << std::endl; + } + + 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; } @@ -101,6 +113,7 @@ test<3>() const int dim = 3; const int spacedim = 3; + deallog.push("basic tetrahedron"); { std::vector> vertices; vertices.emplace_back(0, 0, 0); @@ -115,7 +128,9 @@ test<3>() process(vertices, cells); } + deallog.pop(); + deallog.push("basic wedge"); { std::vector> vertices; vertices.emplace_back(0, 0, 0); @@ -132,6 +147,79 @@ test<3>() process(vertices, cells); } + deallog.pop(); + + deallog.push("basic pyramid"); + { + std::vector> vertices; + vertices.emplace_back(0, 0, 0); + vertices.emplace_back(1, 0, 0); + vertices.emplace_back(0, 1, 0); + vertices.emplace_back(1, 1, 0); + vertices.emplace_back(0.5, 0.5, 1.0); + + std::vector> cells; + CellData cell; + cell.vertices = {0, 1, 2, 3, 4}; + cells.push_back(cell); + + process(vertices, cells); + } + deallog.pop(); + + deallog.push("general tetrahedron"); + { + std::vector> vertices; + vertices.emplace_back(0, 0, 0); + vertices.emplace_back(0.5, 0.25, 0.1); + vertices.emplace_back(0, 1, 0.2); + vertices.emplace_back(0.25, 0.5, 0.5); + + std::vector> cells; + CellData cell; + cell.vertices = {0, 1, 2, 3}; + cells.push_back(cell); + + process(vertices, cells); + } + deallog.pop(); + + deallog.push("general wedge"); + { + std::vector> vertices; + vertices.emplace_back(0, 0.2, 0); + vertices.emplace_back(1, 0, 0.1); + vertices.emplace_back(0.1, 1.3, 0.2); + vertices.emplace_back(-1, 0, 1); + vertices.emplace_back(1, 1, 1.3); + vertices.emplace_back(0, 2, 1.4); + + std::vector> cells; + CellData cell; + cell.vertices = {0, 1, 2, 3, 4, 5}; + cells.push_back(cell); + + process(vertices, cells); + } + deallog.pop(); + + deallog.push("general pyramid"); + { + std::vector> vertices; + vertices.emplace_back(-2, -1, -1.5); + vertices.emplace_back(1.2, -0.5, -1); + vertices.emplace_back(-0.3, 1.2, -0.5); + vertices.emplace_back(1, 1, 0.0); + vertices.emplace_back(0.5, 0.5, 1.0); + + std::vector> cells; + CellData cell; + cell.vertices = {0, 1, 2, 3, 4}; + cells.push_back(cell); + + process(vertices, cells); + } + deallog.pop(); } int diff --git a/tests/simplex/cell_measure_01.output b/tests/simplex/cell_measure_01.output index b079874e8b..666f623ff3 100644 --- a/tests/simplex/cell_measure_01.output +++ b/tests/simplex/cell_measure_01.output @@ -1,21 +1,53 @@ DEAL::dim=2 spacedim=2: -DEAL::measure: 0.500000 -DEAL::diameter: 1.41421 -DEAL::barycenter: 0.333333 0.333333 -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::barycenter: 0.250000 0.250000 0.250000 -DEAL::diameter_min: 1.41421 -DEAL::diameter_max: 1.41421 -DEAL:: -DEAL::dim=3 spacedim=3: -DEAL::measure: 0.500000 -DEAL::diameter: 1.73205 -DEAL::diameter_min: 1.73205 -DEAL::diameter_max: 1.73205 +DEAL::measure: 0.500000 +DEAL::measure (via quadrature): 0.500000 +DEAL::diameter 1.41421 +DEAL::barycenter 0.333333 0.333333 +DEAL::diameter_min 1.41421 +DEAL::diameter_max 1.41421 DEAL:: +DEAL:basic tetrahedron::dim=3 spacedim=3: +DEAL:basic tetrahedron::measure: 0.166667 +DEAL:basic tetrahedron::measure (via quadrature): 0.166667 +DEAL:basic tetrahedron::diameter 1.41421 +DEAL:basic tetrahedron::barycenter 0.250000 0.250000 0.250000 +DEAL:basic tetrahedron::diameter_min 1.41421 +DEAL:basic tetrahedron::diameter_max 1.41421 +DEAL:basic tetrahedron:: +DEAL:basic wedge::dim=3 spacedim=3: +DEAL:basic wedge::measure: 0.500000 +DEAL:basic wedge::measure (via quadrature): 0.500000 +DEAL:basic wedge::diameter 1.73205 +DEAL:basic wedge::diameter_min 1.73205 +DEAL:basic wedge::diameter_max 1.73205 +DEAL:basic wedge:: +DEAL:basic pyramid::dim=3 spacedim=3: +DEAL:basic pyramid::measure: 0.333333 +DEAL:basic pyramid::measure (via quadrature): 0.333333 +DEAL:basic pyramid::diameter 1.41421 +DEAL:basic pyramid::diameter_min 1.41421 +DEAL:basic pyramid::diameter_max 1.41421 +DEAL:basic pyramid:: +DEAL:general tetrahedron::dim=3 spacedim=3: +DEAL:general tetrahedron::measure: 0.0312500 +DEAL:general tetrahedron::measure (via quadrature): 0.0312500 +DEAL:general tetrahedron::diameter 1.01980 +DEAL:general tetrahedron::barycenter 0.187500 0.437500 0.200000 +DEAL:general tetrahedron::diameter_min 1.01980 +DEAL:general tetrahedron::diameter_max 1.01980 +DEAL:general tetrahedron:: +DEAL:general wedge::dim=3 spacedim=3: +DEAL:general wedge::measure: 1.13283 +DEAL:general wedge::measure (via quadrature): 1.13283 +DEAL:general wedge::diameter 2.58650 +DEAL:general wedge::diameter_min 2.58650 +DEAL:general wedge::diameter_max 2.58650 +DEAL:general wedge:: +DEAL:general pyramid::dim=3 spacedim=3: +DEAL:general pyramid::measure: 2.02375 +DEAL:general pyramid::measure (via quadrature): 2.02375 +DEAL:general pyramid::diameter 3.90512 +DEAL:general pyramid::diameter_min 3.90512 +DEAL:general pyramid::diameter_max 3.90512 +DEAL:general pyramid:: -- 2.39.5