Most mesh generation software seems to be aimed at building nicely shaped elements for FEM. I'm curious about a different situation:

I need to numerically integrate over an irregular region. I don't need to solve for any unknowns in this region -- just integrate a known function. And I need to integrate over many such regions as fast as possible. My current approach is to mesh the region (using either Triangle or Tetgen, depending on 2d vs 3d) and then perform Gaussian quadrature over each individual tri/tet.

However, I suspect that these meshing algorithms (which are a bottleneck) are doing far more work than necessary for this use case. Differences:

  • Small internal angles are fine.
  • The occasional zero volume element would be fine -- integrating over a zero volume element will give just zero.
  • Nonconforming elements are not a problem.

Essentially, anything goes as long as the mesh covers the whole volume and there is no overlap between cells.

Does anyone know of research that addresses this situation?

  • 1
    $\begingroup$ Your function to integrate has which form ? If it is simple and can be integrated in closed form over a simplex, then you can take an arbitrary point, connect it to all the triangles, and compute the sum of the integrals over the so-generated tets. Some of them go outside the volume, but signs will cancel-out. $\endgroup$ – BrunoLevy Sep 29 '16 at 7:20
  • $\begingroup$ @BrunoLevy Cool idea! Note: one needs to multiply the tetrahedron integrals by $\pm 1$ depending on whether the current face normal points towards or away from the chosen center point $\endgroup$ – Nick Alger Sep 29 '16 at 8:06

Since you don't care about the quality of the tetrahedra (or triangles), you could quickly tetrahedralize your shape with the following 2-step procedure that is commonly used in computer graphics and animation:

  1. Partition the input nonconvex polyhedron into a union of convex polyhedra organized in a binary space partition tree (BSP tree), by means of the standard recursive cut plane algorithm (successively cutting the interior volume into 2 pieces with successive facets).

  2. Further break down each of these convex polyhedron into tetrahedra by choosing any vertex and connecting it with edges to all other vertices.

The ordering of the cut facets in step 1 and the choice of vertex in step 2 are can be chosen however you like (this will affect the quality of the resulting triangulation); random is usually a good choice.

Once your shape is tetrahedralized, you can do standard quadrature on each tetrahedron in the usual way.

For a 2D example see the picture below.

bsp based triangulation

The extension to 3D is straightforward but it's more difficult to draw. In 3D the cut edges in the BSP step become cut planes, but the edges in the tetrahedralizing step are still edges.

There are lots of good tutorials about this online that were written for the Quake engine (A popular 3D computer game from the 1990's on which many other subsequent 3d engines were based).

  • 1
    $\begingroup$ I think that a BSP does more than what's necessary in the present context. One can use the (simpler) "ear clipping" algorithm, that does not generate additional vertices (en.wikipedia.org/wiki/Polygon_triangulation). $\endgroup$ – BrunoLevy Sep 29 '16 at 6:21
  • 2
    $\begingroup$ @BrunoLevy I don't think ear clipping works in 3D. In 3D there exist nonconvex polyhedra that are not tetrahedralizable, so any extension of the algorithm would need to create new points somehow. The 2D version only works with the points that are there already. For an example of a nontetrahedralizable polyhedron, imagine extruding a triangle to make a triangular cylinder, triangulating the quadrilateral faces, then twisting the top face relative to the bottom one so that the previously vertical edges nearly touch at the center. There does not exist an initial ear for this shape. $\endgroup$ – Nick Alger Sep 29 '16 at 6:53
  • 1
    $\begingroup$ yes, this is only for the 2D case (forgot to say it in the comment). For the 3D case, BSP would work, but would generate many intersection points (so I'd rather use a constrained Delaunay algorithm such as the one in tetgen, while switching-off points insertion / quality optimization). $\endgroup$ – BrunoLevy Sep 29 '16 at 7:18
  • $\begingroup$ @BenThompson Yeah, hope it helps. By the way, although this answer is how I would solve the problem and should be reasonably fast, also take seriously the suggestions by BrunoLevy. He is a legitimate expert in the field, whereas I am not, and his way might turn out to be better. $\endgroup$ – Nick Alger Oct 1 '16 at 23:23

In 3D, if you are using tetgen, you can deactivate insertion of additional points and optimization of mesh quality. When I want to to that, I am using the following flags:

-p (input data is surfacic)
-O0 (do not optimize mesh)
-YY (do not insert Steiner points on boundary)
-AA (generete region tags for each shell) (if needed)

tetgen command line switches are documented here [1].

However, I am unsure that Gaussian quadrature on very skinny elements will be stable enough. If you do that, I'd recommend to evaluate the quality by running the algorithm on both an optimized mesh and a non-optimized one and compare the results.

[1] http://wias-berlin.de/software/tetgen/switches.html

  • $\begingroup$ Thanks! I can modify the quadrature rule to have more points along the longer dimensions of a skinny tri/tet. Perhaps this would be called an anisotropic tensor product quadrature rule... $\endgroup$ – Ben Thompson Sep 29 '16 at 13:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.