]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement pyramid volume and fix wedge volume. 13766/head
authorDavid Wells <drwells@email.unc.edu>
Fri, 20 May 2022 21:02:13 +0000 (17:02 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sat, 21 May 2022 21:26:10 +0000 (17:26 -0400)
source/grid/grid_tools_nontemplates.cc
tests/simplex/cell_measure_01.cc
tests/simplex/cell_measure_01.output

index 03e867d3cc5bfc401110d709d2c47d71ade4aca8..8a32ccb5d02f98cd7cd9851257079a79890472cf 100644 (file)
@@ -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);
index 8c601aa0f5f75093d1f39e4eda7cedc8703803d3..59525e0a371d4f453d9a120de5e4c2668ac53f33 100644 (file)
 // Test TriaAccessor::measure(), TriaAccessor::diameter(), and
 // TriaAccessor::barycenter().
 
+#include <deal.II/fe/fe_nothing.h>
 #include <deal.II/fe/fe_pyramid_p.h>
 #include <deal.II/fe/fe_simplex_p.h>
-#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/fe_wedge_p.h>
 #include <deal.II/fe/mapping_fe.h>
 
@@ -37,33 +38,44 @@ process(const std::vector<Point<spacedim>> &vertices,
   Triangulation<dim, spacedim> 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<dim>())
-        deallog << "barycenter: " << cell->barycenter() << std::endl;
-    }
-
   std::shared_ptr<MappingFE<dim>> 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<dim>())
     mapping = std::make_shared<MappingFE<dim>>(FE_SimplexP<dim>(1));
   else if (reference_cells[0] == ReferenceCells::Wedge)
     mapping = std::make_shared<MappingFE<dim>>(FE_WedgeP<dim>(1));
+  else if (reference_cells[0] == ReferenceCells::Pyramid)
+    mapping = std::make_shared<MappingFE<dim>>(FE_PyramidP<dim>(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<dim> quad =
+    reference_cells.front().template get_gauss_type_quadrature<dim>(2);
+  FE_Nothing<dim, spacedim> fe(reference_cells.front());
+  FEValues<dim, spacedim>   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<dim>())
+        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<Point<spacedim>> vertices;
     vertices.emplace_back(0, 0, 0);
@@ -115,7 +128,9 @@ test<3>()
 
     process(vertices, cells);
   }
+  deallog.pop();
 
+  deallog.push("basic wedge");
   {
     std::vector<Point<spacedim>> vertices;
     vertices.emplace_back(0, 0, 0);
@@ -132,6 +147,79 @@ test<3>()
 
     process(vertices, cells);
   }
+  deallog.pop();
+
+  deallog.push("basic pyramid");
+  {
+    std::vector<Point<spacedim>> 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<CellData<dim>> cells;
+    CellData<dim>              cell;
+    cell.vertices = {0, 1, 2, 3, 4};
+    cells.push_back(cell);
+
+    process(vertices, cells);
+  }
+  deallog.pop();
+
+  deallog.push("general tetrahedron");
+  {
+    std::vector<Point<spacedim>> 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<CellData<dim>> cells;
+    CellData<dim>              cell;
+    cell.vertices = {0, 1, 2, 3};
+    cells.push_back(cell);
+
+    process(vertices, cells);
+  }
+  deallog.pop();
+
+  deallog.push("general wedge");
+  {
+    std::vector<Point<spacedim>> 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<CellData<dim>> cells;
+    CellData<dim>              cell;
+    cell.vertices = {0, 1, 2, 3, 4, 5};
+    cells.push_back(cell);
+
+    process(vertices, cells);
+  }
+  deallog.pop();
+
+  deallog.push("general pyramid");
+  {
+    std::vector<Point<spacedim>> 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<CellData<dim>> cells;
+    CellData<dim>              cell;
+    cell.vertices = {0, 1, 2, 3, 4};
+    cells.push_back(cell);
+
+    process(vertices, cells);
+  }
+  deallog.pop();
 }
 
 int
index b079874e8b6339198377efeb69a92ee3f1a5df03..666f623ff337465e08ecb755511a1b5089461d11 100644 (file)
@@ -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::

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.