• Cylinders are pretty common in FEM analysis and would be a nice additional benchmark for 3D.

# Mesh a cylinder
from mpi4py import MPI
import meshio

import SeismicMesh

comm = MPI.COMM_WORLD

hmin = 0.10

cylinder = SeismicMesh.Cylinder(h=1.0, r=0.5)

points, cells = SeismicMesh.generate_mesh(
    domain=cylinder,
    edge_length=hmin,
)

points, cells = SeismicMesh.sliver_removal(
    points=points,
    domain=cylinder,
    edge_length=hmin,
)

if comm.rank == 0:
    meshio.write_points_cells(
        "Cylinder.vtk",
        points,
        [("tetra", cells)],
        file_format="vtk",
    )

Adding L-shaped domains for SeismicMesh 3.0.5. The 3D L-shaped domain indeed has some trouble with the pockets between SDFs so that remains to be improved but it's a valid mesh with a bounded minimum dihedral angle.

Note: min() and max() of two SDFs do not always produce a valid SDF. The result is actually not a distance field, but just a lower bound to the real distance to the resulting surface. Here I use the union to approximate the domain which results in errors "deep" in the interior of the geometry but this okay for meshing in these examples...

• 2D

 import numpy as np
 import meshio
 from SeismicMesh import Rectangle, generate_mesh, geometry
 
 hmin = 0.05
 bbox = (0.0, 1.0, 0.0, 1.0)
 
 
 bbox0 = (0.0, 1.0, 0.0, 0.5)
 rect0 = Rectangle(bbox0)
 
 bbox1 = (0.0, 0.5, 0.0, 1.0)
 rect1 = Rectangle(bbox1)
 
 
 corners = geometry.corners
 
 
 def union(x):
     return np.minimum.reduce([rect0.eval(x), rect1.eval(x)])
 
 
 pfix = np.vstack((corners(bbox0), corners(bbox1)))
 
 points, cells = generate_mesh(
     bbox=bbox,
     domain=union,
     edge_length=hmin,
     pfix=pfix,
     verbose=2,
 )
 
 meshio.write_points_cells(
     "L_shape_2d.vtk",
     points,
     [("triangle", cells)],
 )

• 3D

 from SeismicMesh import Cube, generate_mesh, sliver_removal, geometry
 
 hmin = 0.05
 bbox = (0.0, 1.0, 0.0, 1.0, 0.0, 1.0)
 
 tol = hmin/10
 
 
 bbox0 = (0.0, 1.0, 0.0, 1.0, 0.0, 0.50 + tol)
 cube0 = Cube(bbox0)
 
 bbox1 = (0.0, 1.0, 0.5, 1.0, 0.50 - tol, 1.0)
 cube1 = Cube(bbox1)
 
 bbox2 = (0.5 - tol, 1.0, 0.0, 1.0, 0.50 - tol, 1.0)
 cube2 = Cube(bbox2)
 
 corners = geometry.corners
 
 
 def union(x):
     return np.minimum.reduce([cube0.eval(x), cube1.eval(x), cube2.eval(x)])
 
 
 pfix = np.vstack((corners(bbox0), corners(bbox1), corners(bbox2)))
 
 points, cells = generate_mesh(
     bbox=bbox,
     domain=union,
     edge_length=hmin,
     pfix=pfix,
     verbose=2,
 )
 
 points, cells = sliver_removal(
     points=points,
     bbox=bbox,
     domain=union,
     edge_length=hmin,
     verbose=2,
 )
 
 meshio.write_points_cells(
     "L_shape.vtk",
     points,
     [("tetra", cells)],
 )

One more FEMish benchmark that I've used a bunch:

 import meshio
 import SeismicMesh
 
 h = 0.025
 
 rect = SeismicMesh.Rectangle((0.0, 1.0, 0.0, 1.0))
 disk0 = SeismicMesh.Disk([0.0, 0.0], 0.25)
 diff0 = SeismicMesh.Difference([rect, disk0])
 
 disk1 = SeismicMesh.Disk([0.0, 0.0], 1.0)
 quarter = SeismicMesh.Intersection([diff0, disk1])
 
 
 def sizes(x):
     return h + 0.10 * abs(disk0.eval(x))
 
 
 points, cells = SeismicMesh.generate_mesh(domain=quarter, edge_length=sizes, h0=h)
 meshio.Mesh(points, {"triangle": cells}).write("out.vtk")

Feel free to reject, I don't want to add more work for you but just think it would bolster this neat page.

You mentioned being willing to add things...

The graphs are lovely, but someone who is new to FE might not be able to interpret them as to when would be a good reason to use one generator over another, which is a good multi-purpose generator (i.e. if you only learn one, this is the most useful...) You can certainly preface it with "IMHO", but a Recommendation section at the end of the ReadMe would certainly help the utility of the comparison.

• A good way I can estimate the numerically stable timestep and the "goodness" of a 2D/3D simplical mesh is by measuring the spectral radius of the scalar wave operator.
• dt < 2/sqrt(spectral_radius)
• Here I used KMV elements which admit diagonal mass matrices for spatial orders < 5.
• I've confirmed through ad hoc tests that the min. dihedral angle of the mesh reduces the spectral radius, thus increasing the maximum stable timestep.
• The spectral radius can be calculated like this (in Firedrake):

import scipy
import numpy as np

import firedrake as fd
from firedrake import dot, grad
import finat


def estimate_timestep(mesh, V, c, estimate_max_eigenvalue=True):
    """Estimate the maximum stable timestep based on the spectral radius
    using optionally the Gershgorin Circle Theorem to estimate the
    maximum generalized eigenvalue. Otherwise computes the maximum
    generalized eigenvalue exactly
    Currently only works with KMV elements.
    """

    u, v = fd.TrialFunction(V), fd.TestFunction(V)
    quad_rule = finat.quadrature.make_quadrature(
        V.finat_element.cell, V.ufl_element().degree(), "KMV"
    )
    dxlump = fd.dx(rule=quad_rule)
    A = fd.assemble(u * v * dxlump)
    ai, aj, av = A.petscmat.getValuesCSR()
    av_inv = []
    for value in av:
        if value == 0:
            av_inv.append(0.0)
        else:
            av_inv.append(1 / value)
    Asp = scipy.sparse.csr_matrix((av, aj, ai))
    Asp_inv = scipy.sparse.csr_matrix((av_inv, aj, ai))

    K = fd.assemble(c*c*dot(grad(u), grad(v)) * dxlump)
    ai, aj, av = K.petscmat.getValuesCSR()
    Ksp = scipy.sparse.csr_matrix((av, aj, ai))

    # operator
    Lsp = Asp_inv.multiply(Ksp)
    if estimate_max_eigenvalue:
        # absolute maximum of diagonals
        max_eigval = np.amax(np.abs(Lsp.diagonal()))
    else:
        print(
            "Computing exact eigenvalues is extremely computationally demanding!",
            flush=True,
        )
        max_eigval = scipy.sparse.linalg.eigs(
            Ksp, M=Asp, k=1, which="LM", return_eigenvectors=False
        )[0]

    # print(max_eigval)
    max_dt = np.float(2 / np.sqrt(max_eigval))
    print(
        f"Maximum stable timestep should be about: {np.float(2 / np.sqrt(max_eigval))} seconds",
        flush=True,
    )
    return max_dt

• The minimum dihedral angle would be a useful diagnostic to plot for 3D (in additional too the existing cell-quality measure).

• The element can still appear to be well-shaped but have a near zero dihedral angle leading to failure in FEM. I suppose this is already captured implicitly in the number of CG iterations.