As we have seen, seismic wave propagation problems place special demands on mesh generators, including the need for tight control over mesh resolution and aspect ratio, and the need to support extremely large problem sizes. We have developed a fast, stable, and efficient meshing algorithm for generating very large scale meshes, suitable for the large basins we target. Since repeated computations will be performed with a single mesh (one or two dozen earthquake scenarios, each involving thousands of time steps), we have decided to generate and partition each mesh sequentially. However, care must be taken in designing and implementing efficient algorithms for these steps, lest they become bottlenecks.
Mesh generation begins with a database of the material properties of the soil and rock within and around a particular basin. The material properties--the shear-wave velocity, the dilatational-wave velocity, and the density--are estimated throughout the basin from soil borings, from geological records, and from seismic prospecting studies.
The meshing algorithm comprises two steps. First, we generate an octree that resolves the local wavelength of shear waves. The wavelength is known from the shear-wave velocity and the frequency of excitation. We have seen in the previous section that 8-10 nodes per wavelength is sufficient for ``engineering'', or 95%, accuracy when using linear finite elements. When constructing the octree, we enforce the rule that adjacent cells may not differ in edge length by more than a factor of two, producing a balanced octree. This is crucial for producing elements with bounded aspect ratios, since aspect ratios far from one lead to poorly conditioned stiffness matrices, which, in turn, can lead to instability in time integration.
Once a balanced octree is created such that no cell is wider than one-tenth the length of the wave that passes through it, a finite element node is placed at each cell vertex. This set of nodes is then tetrahedralized according to the Delaunay criterion. Delaunay tetrahedralization is performed by a straightforward implementation of the Bowyer/Watson incremental algorithm [5, 25], which constructs the tetrahedralization by adding one node at a time and locally adjusting the mesh to maintain the Delaunay criterion.
We have found that the Bowyer/Watson algorithm is occasionally sensitive to floating-point roundoff error; tetrahedral mesh generation can fail dramatically because of roundoff when processing near-degenerate geometric features. Such failures became increasingly common for us as the size of our meshes grew. To overcome this problem, we have developed a method for fast exact arithmetic that is particularly well-suited for certain tests that arise in computational geometry codes [20]. Our method is used to construct predicates that determine whether a point falls to the left or right side of a line, or whether a point falls inside or outside a sphere. These predicates are adaptive in the sense that they only use exact arithmetic to the extent it is needed to ensure a correct answer. Hence, if a point falls very close to a line, high precision arithmetic may be needed to resolve which side of the line it falls on; if a point is far from a line, approximate arithmetic will suffice, so the test can be performed quickly. Because the latter case is far more common, our exact arithmetic predicates are on average only slightly slower than ordinary, non-robust floating-point predicates, and our Delaunay tetrahedralization code runs quickly while ensuring the integrity of its results.
Our use of the Delaunay tetrahedralization of the vertices of a balanced octree guarantees that element aspect ratios are bounded, and that element sizes are chosen appropriately so that wavelengths are sufficiently resolved without unnecessary resolution (provided the material properties do not vary too rapidly).