|Home | About | Journals | Submit | Contact Us | Français|
We present Laguerre Voronoi based subdivision algorithms for the quadrilateral and hexahedral meshing of particle systems within a bounded region in two and three dimensions, respectively. Particles are smooth functions over circular or spherical domains. The algorithm first breaks the bounded region containing the particles into Voronoi cells that are then subsequently decomposed into an initial quadrilateral or an initial hexahedral scaffold conforming to individual particles. The scaffolds are subsequently refined via applications of recursive subdivision (splitting and averaging rules). Our choice of averaging rules yield a particle conforming quadrilateral/hexahedral mesh, of good quality, along with being smooth and differentiable in the limit. Extensions of the basic scheme to dynamic re-meshing in the case of addition, deletion, and moving particles are also discussed. Motivating applications of the use of these static and dynamic meshes for particle systems include the mechanics of epoxy/glass composite materials, bio-molecular force field calculations, and gas hydrodynamics simulations in cosmology
Particles are smooth functions with compact support. Particle systems are used for modeling a number of physical world scenarios ranging from cosmological systems and plasma physics to molecular systems. The applications are wide and varied and include chemistry, material science, and bio-engineering. A typical astrophysics computation attempts to follow the evolution of a system of cosmological boundaries, modeled as a particle system. An important computation performed on the particle system is that of force calculations, i.e. given n particles compute the effect of the gravitational or hydrodynamic force on a given particle by the other particles. These are often termed as n-body problems. Other characteristics of the problem require tracking the dynamic structure of the particle system as particle systems move, appear and disappear.
Smooth particle hydrodynamics or SPH  is a Lagrangian numerical hydrodynamics method that uses a discrete description of gas or fluid particles in place of a continuum model. The value of a field at a point in space is given by the sum of any contributing particles, and their radially symmetric kernels, present at that point. These particles essentially act as moving centers for interpolation and carry mass, energy and the velocity of the local flow. For the meshing of particle systems, it suffices to consider particles as idealized balls, or radially symmetric domains of support of their kernels. So henceforth, we consider our particle systems to be a collection of circles in two dimensions, and spheres in three dimensions, with possibly different radii.
We present Laguerre Voronoi based subdivision algorithms for the quadrilateral and hexahedral meshing of particle systems within a bounded region, for two and three dimensions, respectively. Extensions of the basic scheme to dynamic re-meshing in the case of addition, deletion, and moving particles are also discussed.
A motivating application for finite element mesh constructions of a static particle system is for the accurate calculation of stress distributions for composite materials consisting of silica substrates with embedded particles, subjected to various loads . The re-meshing problem for time dependent particle systems arise in gas hydrodynamics simulations essential in the computational investigation of the formation of large scale structures, such as galaxies and galaxy clusters, in the universe . A third motivating application domain is the use of finite element meshes of particle systems for biophysics simulations of molecular interactions in ionic solvent .
Laguerre Voronoi diagrams [5, 6, 7] or its variants like STC , embedded Voronoi Graphs , , or medial surfaces  have been extensively used to generate hexahedral meshes. A variant of the Voronoi diagram called the Space Twist Continuum (STC)  has also been successfully used at the Sandia Labs in their meshing software CUBIT. STC is a way of representing the dual of the hexahedral mesh as a simple non-degenerate arrangement of surfaces. The STC is built in an incremental fashion using an Advancing Front Techniques (AFT) like the Whisker Weaving algorithm .
As reviewed in  , there are indirect and direct methods for unstructured quad/hex mesh generation. The indirect method is to generate triangular/tetrahedral meshes first, then convert them into quads/hexes. The direct method is to generate quads/hexes directly without first going through triangular/tetrahedral meshing . Alternatively, Pascal et al give an overview of the grid based approaches for Hexahedral meshing . Either higher order elements can be used to represent the exact boundary of the surface by intersecting it with the grid, or body-fitting techniques can be used to represent these boundary elements , . Dual contouring has also been used by Zhang et al , , for extracting hexahedral meshes from volumetric data.
Subdivision defines a smooth curve or surface as the limit of a sequence of successive refinements of some linear input . As a structured method, quad/hex mapped meshing  generates the most desirable meshes if opposite edges/faces of the domain to be meshed have equal numbers of divisions or the same surface mesh. However, it is always difficult to decompose an arbitrary geometric configuration into mapped meshable regions. In the CUBIT project  at Sandia National Labs, several techniques have been attempted to automatically recognize features and decompose geometry into mapped meshable areas or volumes.
The simplest method for post mesh quality improvement is based on Laplacian smoothing which relocates the vertex position at the average of the nodes connecting to it . Optimization-based smoothing tends to yield better results but it is more expensive than Laplacian smoothing , , 
Two dimensional particles are circles of possibly different radii. A typical input particle system along with its output conforming quadrilateral mesh, inside a rectangular bounding box, is shown in Figure 1.
Construct the 2D Laguerre Voronoi diagram of the center points of the circles within its rectangular bounding region (see figure 1). The weight of a circle center is chosen to be proportional to the radius of each circle. Each Voronoi cell is a convex polygon.
Contract the short edges into either of its endpoint vertices, or to the center point of the edge, unless boundary constraints as mentioned below are violated (see figure 2). An edge is regarded as short if its length is less than a user specified threshold value.
Note that this edge contraction may eliminate triangles (see figure 3). The aim of this small edge contraction step is to avoid producing tiny quadrilateral elements. Also note that if an edge contraction causes an intersection between the boundaries of a polygonal cell and the circle perimeter, then we do not carry out this edge contraction.
For each Voronoi cell, connect the circle center with the vertices, and the midpoint of each Voronoi cell edge. Partially delete the inner portion (i.e. within each circle) of the line segment that connects the circle center to the midpoint of each Voronoi edge, thereby creating quadrilaterals both inside and outside each circle, i.e. a quad tesselation. (see figures 5 and and66).
Next, further refinement of the quad tesselation is accomplished by first a quad mesh splitting along the radial and angular directions from each circle center, followed by an averaging/smoothing recalculation  to place all quad mesh vertices (see figure 7). Note the selective splitting of the inner most radial edges from each circle center, to preserve quad mesh topology (see figures 8. The aim of the splitting step is to achieve quad mesh adaptivity.
The averaging/smoothing based repositioning of the mesh vertices can be carried out, in a variety of ways. We choose to adopt a centroid smoothing scheme, where each vertex is repositioned in the centroid of the centroids of the maximal dimensional cells incident to it (taking care of boundary vertices as well) . While we do not have a guarantee for mesh quality, the current averaging scheme tends to produce quads with very good aspect ratio. See figures 9,,1010,,1111,,1212 for examples of our implementation of this algorithm.
In this section we present the solution in three dimensions.
Construct the 3D weighted Voronoi diagram of the center points the spheres (see figure 13). The weight of a sphere center is chosen to be proportional to the radius of the sphere.
Contract each of the short edges into single vertices. An edge is regarded as short if its length is less than a user specified threshold value. Note that this contraction may eliminate triangles (see Fig 3). Also note that this merging may lead to the intersection between the boundaries of polyhedral cells and the surface of the sphere, and so a check is made, and such edge contractions are not carried out. The aim of this step is to avoid producing very tiny hexahedral elements.
Decompose each face of the 3D Voronoi cells into quads (see Fig 14). For each Laguerre Voronoi cell, connect the particle center with the Voronoi vertices and thereby form a pyramidal partition of the cell. Each pyramid is further split into a hexahedron and another pyramid by a quad face generated by intersection with the surface of the sphere.
Adaptive subdivision of each cell is easily performed in the radial direction (see figure 15).
At this stage one has almost obtained a complete adaptive hexahedral mesh of the bounding region, conforming to the boundary of the particle spheres. The mesh elements incident to the particle centers are quad pyramids, while the rest of the mesh elements are hexahedra. By repeated applications of the radial subdivision, the pyramidal elements can be restricted to very small spheres surrounding each particle center. In many applications, where calculations are desired close to the surface of the particle spheres, or only on its exterior, this near hexahedral mesh suffices. See also figure 16 of an example mesh generated for a small portion of a heterogenous elastic material with spherical inclusions with different bulk properties to the rest of the substrate.
Further improve the hexahedral mesh quality, by applying the averaging/smoothing based repositioning of the mesh vertices using the centroid smoothing scheme of , which tends to produce hexahedra with good aspect ratio. This averaging/smoothing can also be conducted with varied crease rules, to produce user desired mesh anisotropy.
N.B. A current disadvantage of the centroid smoothing scheme  that we use for mesh quality improvement, is the lack of a guarantee of producing “good” (e.g. bounded aspect ratio, positive Poisson ratio) meshes. Furthermore, in three dimensions there is the occasional failure to maintain convexity of the subsequently partitioned hexahedra after several iterations of averaging/smoothing. While in our current implementation, we do not perform the step that causes “bad” hexahedra, we are currently working on an improved weighted averaging/smoothing scheme for hexahedral mesh quality improvement.
In the case where a complete hexahedral element mesh is required both inside and outside the spheres, one can use the following alternate hexahedral scaffolding scheme. First, split each face of the Voronoi cell into triangles. Next connect each vertex of the triangulated Voronoi faces to the center of their respective particle centers, decomposing the entire bounding region into tetrahedra. The tetrahedra are further decomposed by the (piecewise linear) surfaces of the particle spheres into triangular prisms, and tetrahedra incident to the center of the particles. Each triangular prism can be split into 3 hexahedra and the tetrahedra at the particle centers, are split into four hexahedra by adding three vertices on the edges and one vertex on the face (see figures 17 and and18).18). After this step, the entire bounding region, with the particle system, is partitioned into hexahedra.
The hexahedral mesh constructed is adaptive in the sense that it becomes denser in the regions that are close to the surfaces of the sphere. When averaging/smoothing (recursive subdivision) rules are used, the mesh is uniform in the sense that at the regions which are the same distance to a sphere, the elements are similar to each other. See figures 19, ,2020.
Our extensions to the re-meshing of dynamic particle systems is based on dynamic Laguerre Voronoi diagrams. The partitioning and recursive subdivision follows directly from the instantaneous position of the particle system and its Voronoi diagram. Dynamic Voronoi diagrams have been well studied [28, 29, 30, 31, 32, 33, 34, 35]. Dynamic Laguerre Voronoi diagrams have also been considered [36, 37].
We have incorporated our quadrilateral and hexahedral meshing schemes to two and three dimensional stress/strain finite element simulations involving epoxy/glass composite materials under static and dynamic loading scenarios in computational mechanics . We are currently also using this mesher to develop a fast non-uniform fast-Fourier transform based n-body gravitational force calculation, coupled to smooth particle gas hydrodynamics simulations in computational cosmology . Further, we hope to apply extensions of this scheme to Poisson-Boltzmann and Smoluchowski equations  as well as additional biophysics simulations.
Sincere thanks to the Dr. Guoliang Xu and Chandrasekhar Karapalem for several useful discussions and their invaluable help in the implementation of this meshing scheme. This work was supported in part by NSF grants ACI-0220037, EIA-0325550, a UT-MDACC Whitaker grant, and a subcontract from UCSD 1018140 as part of the NSF-NPACI project, Interaction Environments Thrust.