## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# The Voronoi Implicit Interface Method for computing multiphase physics

Edited by Leslie Greengard, New York University, New York, NY, and approved October 11, 2011 (received for review July 21, 2011)

## Abstract

We introduce a numerical framework, the Voronoi Implicit Interface Method for tracking multiple interacting and evolving regions (phases) whose motion is determined by complex physics (fluids, mechanics, elasticity, etc.), intricate jump conditions, internal constraints, and boundary conditions. The method works in two and three dimensions, handles tens of thousands of interfaces and separate phases, and easily and automatically handles multiple junctions, triple points, and quadruple points in two dimensions, as well as triple lines, etc., in higher dimensions. Topological changes occur naturally, with no surgery required. The method is first-order accurate at junction points/lines, and of arbitrarily high-order accuracy away from such degeneracies. The method uses a single function to describe all phases simultaneously, represented on a fixed Eulerian mesh. We test the method’s accuracy through convergence tests, and demonstrate its applications to geometric flows, accurate prediction of von Neumann’s law for multiphase curvature flow, and robustness under complex fluid flow with surface tension and large shearing forces.

There are a host of physical problems that involve interconnected moving interfaces, including dry foams, crystal grain growth, mixing of multiple materials, and multicellular structures. These interfaces are the boundaries of individual regions/cells, which we refer to as phases. The physics, chemistry, and mechanics that drive the motion of these interfaces are often complex, and include topological challenges when interfaces are destroyed and created.

Producing good mathematical models that capture the motion of these interfaces, especially at degeneracies, such as triple points and triple lines where multiple interfaces meet, is challenging. Building robust numerical methods to tackle these problems is equally difficult, requiring numerical resolution of sharp corners and singularities, and recharacterization of domains when topologies change. A variety of methods have been proposed to handle these problems, including (*i*) front tracking methods, which explicitly track the interface, modeled as moving segments in two dimensions and moving triangles in three dimensions, (*ii*) volume of fluid methods, which use fixed Eulerian cells and assign a volume fraction for each phase within a cell, (*iii*) level set methods (1), which use an implicit formulation to represent the interface, and treat each region/phase separately, followed by a repair procedure which reattaches the evolving regions to each other (1⇓–3), and (*iv*) diffusion generated motion which combine diffusion via convolution with reconstruction procedures to simulate multiphase mean curvature flow (4). Although there are advantages and disadvantages to each of these approaches, it has remained a challenge to robustly and accurately handle the wide range of possible motions of evolving, highly complex interconnected interfaces separating a large number of phases under time-resolved physics.

In this paper, we present a numerical method for tracking the interface in general multiphase problems. We call this method the Voronoi Implicit Interface Method (VIIM), because it relies on a robust interaction between Voronoi diagrams and implicit interface methods. The method has a variety of important features, including:

Efficiency and consistency: The method uses a fixed Eulerian mesh, and a single function plus an indicator function to track the entire multiphase system. Phases are coupled together in a consistent fashion, with no gaps, overlaps, or ambiguities.

Multiple junctions and topological change: Triple points and triple lines, where more than two phases touch, as well as breakage, merger, creation, and disappearance of phases, are all handled naturally. Transition events occur automatically, with no special attention paid to discontinuous topological change.

Coupling with time-dependent physics: The method uses a physical time step, and complex physics may be solved at each time step and correctly incorporated into the interface evolution. Feedback from the physics affects the interface, and changes in the interface affects the physics.

Accurate calculation of geometric quantities and extension to any dimension: The method allows one to accurately calculate curvature and geometric constraints as part of the interface evolution, and is fundamentally unchanged, regardless of the dimension of the problems.

After a short background about implicit interface level set methods, we introduce the basic VIIM. We next summarize some of the convergence tests, and then demonstrate the basic method on a collection of geometrical and fluid mechanics problems.

## The Voronoi Implicit Interface Method

### Background.

We begin by recalling the basics of level set methods, introduced by Osher and Sethian (5), which rely in part on the theory of curve and surface evolution given in refs. 6 and 7, and the link between front propagation and hyperbolic conservation laws discussed in ref. 8. Start with an interface Γ separating two phases *A* and *B*, and moving with a given speed *F* in its normal direction. We can interpret this interface as an *n* - 1 dimensional hypersurface in , and then embed this interface as the zero level set of a signed distance function *ϕ*(*x*), so that *ϕ*(*x*) is the initial (signed) distance to the interface Γ, chosen to be positive inside one phase and negative inside the other. The zero level set {*ϕ* = 0} corresponds to the original interface Γ separating *A* and *B*. Adding time, we can ensure that the zero level set of this function *ϕ* always corresponds to the evolving interface Γ through the initial value partial-differential equation given by Here, the speed function *F* may depend on a variety of factors, including local geometry (normal direction and curvature), integral properties (enclosed area/volume), as well as the solution of complex partial differential equations (PDEs), such as fluid flow, material elasticity, diffusion, etc., with jump conditions and source terms supplied by the interface position.

This initial value PDE is approximated using upwind finite difference schemes: In the original formulation presented in ref. 5, the interface is embedded throughout the entire computational domain and hence adds unnecessary labor. More sophisticated versions employ the Narrow Band Level Set Method (9), which adaptively focuses computation in a small thin tube around the moving interfaces, and hence reduces the complexity of the algorithm to roughly the number of elements on the front. For details and a comprehensive review on the level set method, see ref. 1.

Consider now three different regions, or phases, A, B, and C. The standard level set construction, which relies on a signed distance function whose zero level set is the interface of interest, breaks down at triple points where three interfaces meet, because there is no longer an inside and an outside. Various authors have approached this problem by using multiple level set functions, in most cases by simply dividing the region into three separate level set functions, namely, inside A and not inside A, inside B and not inside B, etc., advancing each level set function separately. Although more economical classifications can be formed, they all boil down to using multiple functions in multiple regions.

### Voronoi Sets.

Given a set of *m* nodes in , the Voronoi diagram is the decomposition of into different cells such that all the points in a given cell are closer to one node than any other. The boundaries of these cells are then a set of points, which are equidistant from at least two nodes. This construction may be extended to a set of regions *A*_{1},…,*A*_{m} instead of nodes, which may be curves, surfaces, or any other objects: The tessellation of the domain into cells still provides a Voronoi subdivision such that all points in a particular cell are closer to region *A*_{i} than to any other, and the boundaries of these cells are equidistant from at least two regions.

Using this idea, we now formulate the VIIM, which is a PDE-based method based on a single function, defined on a fixed Eulerian background mesh, and tracks evolving interface problems of multiple regions in two and three dimensions, regardless of the complexity of the multiply-connected junctions.

### The Voronoi Implicit Interface Method.

Given a domain in , we imagine a collection of phases such that each point *x* in the domain is either in a unique phase, or on the boundary between two or more phases: This boundary (which is really a collection of boundaries) will be known as the interface. These are the initial conditions for the problem; we further assume that we are given a speed *F* (or velocity field **u**) defined on the interface.

At each point *x* of the domain, we let *ϕ*(*x*) be the distance from *x* to the closest point on the interface: *ϕ*(*x*) is the embedding of the interface as the zero level set of an unsigned distance function. Furthermore, we create an extension velocity, which defines a speed function *F*_{ext}(*x*) in all of the domain. This extension velocity can be constructed so that it smoothly varies away from the interface itself.

The key idea is to rely on a main feature of level set methods: The level set corresponding to the interface is embedded in a family of nearby level sets. Thus, the motion of the zero level set corresponding to the interface is bracketed by the motion of surrounding level sets. This property is an immediate consequence of a comparison theorem under suitable restrictions on the speed function that moves the neighboring level sets, and these restrictions are often satisfied through the construction of extension velocities as outlined in ref. 10.

We utilize this property to construct the VIIM. The motion of the interface in a multiphase system is determined by nearby level sets, given by the *ϵ*-level sets of *ϕ*, where *ϵ* > 0. Although the interface where multiple phases touch may have high-order junctions such as triple points and triple lines, the *ϵ*-level sets are hypersurfaces, which exist solely in a single phase and do not contain such degeneracies. The motion of these nearby interfaces controls the evolution of the multiphase interface through the Voronoi construction, as follows.

The VIIM method consists of three steps. First, we evolve the unsigned distance function a small time step Δ*t* by solving the governing PDE. For a short period of time, the neighboring level sets will remain nice hypersurfaces. Next, we reconstruct the interface after time Δ*t* as the Voronoi interface of the nearby *ϵ*-level sets. Obviously, Δ*t* and the choice of the value of *ϵ* are linked: We need to choose *ϵ* large enough (or conversely, Δ*t* small enough) so that the *ϵ*-level sets remain hypersurfaces during the time step. Finally, we use the Voronoi interface to construct a new unsigned distance function at time Δ*t*.

### Two Illustrative Examples.

We illustrate with two simple examples. Imagine an interface (Fig. 1) in one spatial dimension, which is a single point located at *x* = *X*(*t*) at time *t*, moving with speed *F*(*x*,*t*), thus . At *t* = 0, the zero level set of the unsigned distance function *ϕ* corresponds to the initial location of the point, and is at an extremum of *ϕ*. However, at the nearby level sets with value *ϵ*, the distance function is smooth. Hence, we can update the distance function everywhere away from the zero level set in a straightforward manner. What remains is to obtain a suitable definition of the zero level set at time Δ*t* in such a manner that it corresponds to the location of the front at this time. We can do so by defining the zero level set as the point equidistant from the *ϵ* level set from each side: This single point is the Voronoi set.

In higher dimensions, and in the presence of high-order junctions, this same construction works. In Fig. 2 (*Left*) we see a triple point of an unsigned distance function, which is nonnegative everywhere. In Fig. 2 (*Middle*) we show the *ϵ*-level sets, and, Fig. 2 (*Right*), the Voronoi interface reconstructed from those *ϵ*-level sets.

## Algorithm and Implementation

Given a multiphase problem with *N* different phases, our goal is to transform the above into a numerical method. Begin by defining the unsigned distance function *ϕ*(*x*), which is zero on the boundary of each phase, and measures the positive distance to the boundary everywhere else. Let there also be an indicator function *χ*(*x*), which indicates the phase in which the point *x* is located. Let us further suppose that we are given a speed function *F*, defined in all of space (and which can depend on position, normal, and curvature, as well as associated physics.) We fix *ϵ* and Δ*t*. Our goal is to execute the following sequence:

Advance the unsigned level set distance function in time by solving the standard level set equation, using forward Euler with a time step Δ

*t*depending on*F*.Reconstruct the unsigned distance function

*ϕ*at the new time step, using the Voronoi interface reconstructed from the*ϵ*-level sets and the old indicator function*χ*.Update the speed function

*F*to reflect the appropriate physics, mechanics, etc., by solving the relevant equations of motion. This calculation may involve computing geometric quantities from the Voronoi interface, as well as extracting jump conditions, source terms, etc., from the Voronoi interface as input to PDE solvers, as well as building extension velocities to define the speed within the domain. Loop to 1.

This formulation, by construction, cannot create overlaps or vacuums: Each point in the domain is part of a unique phase at all times. Our reconstruction differs from existing algorithms in that there is neither a no overlap/vacuum condition, nor is there a penalty term that penalizes region of vacuums or overlap. The interface is always well-defined at every time step. Furthermore, the time *t* in the above algorithm has physical meaning: The interface moves under physical speed *F* during the time step Δ*t*.

The key to the above is an accurate and fast way to execute step 2, namely to construct the new unsigned distance function from the *ϵ*-level sets. This procedure involves the computation of distances from the *ϵ*-level sets on an Eulerian mesh, and so the core tool is an Eikonal equation solver. We make extensive use of the accurate and fast Eikonal reinitialization algorithm developed by Chopp (11), which utilizes bicubic [tricubic in three dimensions (3D)] interpolation to accurately initialize an initial band. In more detail, if a grid cell is identified to contain the zero level set of a function *ψ* (which need not be a distance function), then *ψ* is interpolated in that cell using a bicubic/tricubic patch. A Newton solver is then used to find the closest point on the interpolated interface, from which the exact distance to the interpolated interface is computed at the nodes of all such grid cells. Note that this construction does not require explicit construction of the zero level set of *ψ*. This procedure creates a small initial band that can be input to the efficient Fast Marching Method (12), which is a Dijkstra-like ordered upwind finite difference scheme for solving the Eikonal equation outside this initial band. A different Dijkstra-like control theoretic discretization of the Eikonal equation stemming from optimal control was developed in ref. 13, and we refer the reader to ref. 14 for a detailed discussion and extensions of Fast Marching Methods to more general front propagation problems.

VIIM uses this methodology as follows. First, without explicit construction, by locating cells containing the *ϵ*-level sets and building bi/tricubic patches, we solve the Eikonal equation with zero boundary value on these sets. With this solution, we locate cells containing the Voronoi set, build bi/tricubic patches, and again solve the Eikonal equation to provide the new unsigned distance function.

## Convergence Tests

An extensive convergence analysis of the VIIM has been performed, analyzing convergence for various values of time step Δ*t*, grid size *h*, and choice of parameter *ϵ*. One may fix *ϵ* from the outset, and study convergence as the grid size *h* vanishes. Here, we instead couple *ϵ* to the grid size and choose *ϵ* = *αh*, where *α* is a constant. Generally speaking, larger values of *ϵ* give larger numerical errors, but with judicious programming, α can be taken as any nonnegative number. Here, we summarize some of the key findings.

### Convergence Tests for Smooth Interfaces.

Consider a problem involving a single interface separating two phases, for example, a circle collapsing under its curvature. In this case, no corners occur in the moving interface. Using second-order finite difference schemes in space and second-order in time, the VIIM gives second-order convergence to the known exact solution.

### Convergence Test of T-Junction to Y-Junction.

For multiphase problems with more than two phases, and in the presence of triple points, etc., exact solutions are less well-known, and hence we rely on grid convergence. We first analyze the motion of a single T-junction that evolves into a Y-junction under curvature flow, that is, *F* = *κ*. Suitably interpreted, it is well-known that curvature flow applied to a triple point results in 120° angles: This property is known as Young’s law, and arises as a natural consequence of interpreting curvature flow as minimizing length (surface area in 3D). We consider Dirichlet boundary conditions, in which the junction is anchored at the boundary of the domain [0,1]^{2} in two dimensions, and evolve over a time interval 0 ≤ *t* ≤ *T*. In Fig. 3, we show a snapshot of the evolution, which shows that the junction develops 120° angles, and ultimately converges to a triple point with three straight lines. In Table 1, we show convergence results using grid refinement for various grid sizes, coupling the choice of *ϵ* to *h*, where *h* is the size of one grid cell of a uniform Cartesian grid. Here, is a metric using the Hausdorff distance *d*_{H} in space and measures the convergence of the interface evolution in both space and time. The results show first-order convergence. It is important to note that first-order convergence at triple points is probably all that can be expected in this framework.

### Convergence Test Using von Neumann–Mullins’ Law.

Imagine now a large number of phases with a connected network of interfaces and many triple points. von Neumann and Mullins, under the assumption of Young’s law (15), showed that if the speed of the interface is a constant multiple γ of curvature,* then the rate of change of area of a particular phase depends on its number of sides, namely We use this result to further test the convergence of the VIIM. Start with 25 randomly positioned phases and apply curvature flow (*γ* = 1, grid size 256 × 256) with periodic boundary conditions (Fig. 4). Each phase grows, shrinks, or conserves its area, with the rate of area change a constant given by the above law. When phases change neighborship with other phases or collapse, the topology of the network changes and may alter the number of sides of each phase. Throughout the entire evolution, it follows that the area of each phase is a piecewise linear function of time. In Fig. 4 we plot the area of a selected set of nine phases. Our results show correct match with von Neumann–Mullins’ law throughout the evolution.

### Three-Dimensional Convergence Tests.

Under curvature flow with Neumann boundary conditions, Fig. 5 shows the evolution of a 3D analogue of a T-junction, which has four different triple lines and a quadruple point. Qualitatively, we see that the surfaces make 120° angles at triple lines, which is one of Plateau’s laws on the shapes of soap bubbles in a foam. Using grid refinement to measure convergence, Table 2 shows that the VIIM converges at first order, in both space and time.

### Efficiency.

Formally, the operation count on the method per time step is , where *N* is the number of grid cells containing the interface, and *k* is the size of the narrow band. The Voronoi reconstruction takes twice as long as the reinitialization step in the narrow band level set method, because there are two Eikonal solves.

## Geometric Flows with Constraints

By adding a discontinuous source term to the right-hand side of the mean curvature flow equation, we can simulate mean curvature flow with volume conservation. The modified forward Euler step is where Here denotes the volume [area in two dimensions (2D)] of phase *i* at time step *n*, is the initial volume (area) of phase *i*, and is its surface area (perimeter in 2D) at time step *n*. In effect, the source term *s*^{n}|∇*ϕ*^{n}| grows or shrinks each phase equally around its boundary by an amount that corrects for any mass loss/gain.

Despite each phase potentially growing or shrinking at different rates, the VIIM robustly and smoothly handles this discontinuity and conserves volume almost exactly. Fig. 6 demonstrates the method in 2D and 3D on a set of 100 initially random phases (using zero Neumann boundary conditions). Mean curvature flow minimizes the total length (surface area in 3D) of all interfaces, which, subject to the constraint of area (volume) conservation, eventually attains an equilibrium (Fig. 6). Various topological changes occur, and at all times, triple points (lines) have 120° angles.

## Fluid Flow with Permeability

We now use the VIIM to incorporate fluid dynamics into the dynamics of dry foams. A foam is a collection of gas bubbles separated by a liquid component and is considered dry when the liquid makes up less than 10% of the total volume. In the multiphase evolution of dry foams, the following occurs:

Fluid mechanics of the gas drives the flexible membranes. The flow is incompressible, with jumps in pressure across the membrane due to surface tension effects. Membranes meet at 120° junctions: triple points in 2D and Plateau borders in 3D.

Membranes may be permeable to gas, which diffuses through the boundary across phases, leading to diffusive coarsening. In general, phases with a large number of faces grow whereas those with a small number of phases shrink. Once these phases disappear, the topology of the interconnected network changes, leading to different growth laws for each phase.

We consider a collection of membranes, moving under the combined effects of surface tension, fluid mechanics, and permeability. We assume that membranes are massless and thin, and follow the model used in ref. 18, in which the gas is modeled as an incompressible Newtonian fluid, satisfying the Navier–Stokes equations with density *ρ* and viscosity *μ*. Surface tension at the interface induces a fluid flow in addition to any external forcing flow, which in turn moves the interface, providing a feedback mechanism between membrane dynamics and fluid physics. In more detail, interface surface tension induces a pressure jump of [*p*] = *σκ*, where *σ* is the surface tension coefficient and *κ* is the mean curvature of the interface. In the absence of permeability, the interface is advected by the velocity **u** of the gas, taken as continuous across the interface. With permeability, the diffusion rate of gas per unit length is proportional to the pressure difference (15). As in ref. 18, this permeability is modeled with a slip of the interface in the normal direction relative to **u**, hence **u** - *Mσκ***n**, where *M*≥0 is a constant, denoting the amount of permeability and **n** is the normal unit of the interface (with the same orientation as that used to calculate *κ*). We adopt a continuum surface tension model, whereby surface tension at the interface becomes a body force through the use of a Dirac delta function. The governing equations of motion are thereforewhere **F** is any additional body forces (such as gravity) and **st** is a surface tension body force appropriate to a multiphase system. For two-phase fluid flow, this term takes the form **st** = -*σκδ*(*ϕ*)∇*ϕ*, where *ϕ* is a signed level set function. In our case of multiple phases which meet at triple junctions (and quad points in 3D), curvature needs to be suitably defined. It is both physically and mathematically natural to take the same definition applied to each phase (gas bubble) separately, sum all of these, and normalize by a factor of two. This formulation is physically consistent, because each bubble is separated from the others by a thin membrane, and it is mathematically natural because the resulting formula effectively enforces Young’s law at triple points, i.e., triple junctions instantaneously obtain 120° angles. Using this definition, we have , where *ϕ*_{i} is a signed level set function for phase *i* and *N* is the total number of phases. Although the normal is not well-defined at triple junctions, the curvature times the normal, *κ*(*ϕ*_{i})∇*ϕ*_{i}, is well-defined as a distribution, and is a Dirac delta function with magnitude related to the angle of the corner.

A careful analysis and detailed computation of the equations of motion were first performed in ref. 18, in which an explicit 2D Lagrangian-based front tracking method together with an immersed boundary method was used to couple gas fluid mechanics to membrane dynamics under the effects of large shear forces, showing a variety of phenomena, including the balance between evolution toward structures satisfying Young’s law and the evolving fluid mechanics. Here, the VIIM allows topological changes in an implicit fashion, and can track 3D flow. Numerical details will be presented elsewhere. Briefly, we use a fixed Cartesian mesh and a second-order projection method, using upwinding for the advection term and Crank–Nicholson for the diffusion term, with the interface surface tension term mollified to the background grid.

### Two-Dimensional Flow With and Without Permeability.

Consider a foam being agitated by a strong external forcing which first spins counterclockwise, then settles, and then reverses. We take a unit square [0,1]^{2} with periodic boundary conditions, with *ρ* = 1, *σ* = 1, and *μ* = 0.005, chosen so that effects of inertia, viscosity, and surface tension are of similar importance. The agitator is an external force **F** in the Navier–Stokes momentum equations: **F**(*x*,*y*,*t*) = 15(sin *πx* sin 2*πy*,- sin 2*πx* sin *πy*) sin *πt*, corresponding to a force which first spins in the counterclockwise direction and then reverses. The factor of 15 in **F** was chosen to give a relatively strong shearing/spinning force that dominates the stabilization effects of surface tension. Starting with 25 random phases, we evolve the Navier–Stokes equations with the agitator forcing. First, consider *M* = 0, i.e., no permeability (Fig. 7). Flow undergoes significant shearing, causing considerable rearrangement of the phases and topological change. Stream function plots computed from the velocity field **u** show flow strongly affected by agitator forcing, but localized in nature due to effects of surface tension.

Fig. 8 shows the effect of permeability (*M* = 0.05). Diffusion of gas across the membranes cause phases to collapse, others to grow, all of which occurs through large shearing forces.

### Three-Dimensional Flow With and Without Permeability.

We repeat the above in 3D, further demonstrating the capabilities of VIIM. The domain is a unit cube [0,1]^{3} and the external force is the same but with no forcing in the *z*-direction. We start with 125 random phases and evolve the Navier–Stokes equations with the agitator forcing. Once again, the flow undergoes significant shearing, causing considerable rearrangement of the phases and many topological changes. In the case of no permeability (Fig. 9, *Top*), all phases conserve their volume, whereas in the case of permeability (Fig. 9, *Bottom*), diffusion causes some phases to collapse and others to grow under the large shear. In 3D, the rate of volume change matches with a recently found generalization of von Neumann’s law (19, 20).

## Conclusion

The Voronoi Implicit Interface Method is a robust, accurate, and efficient numerical method to track a large number of evolving interfaces moving under coupled complex interactions of geometry, physics, constraints, and internal boundary conditions. The accuracy of the method has been demonstrated under both convergence studies and verification of von Neumann–Mullins’ law, and applications have been presented for fluid flow in the presence of permeability and large shearing forces.

## Acknowledgments

This research was supported in part by the Applied Mathematical Sciences subprogram of the Office of Energy Research, US Department of Energy, under contract number DE-AC02-05CH11231, by the Division of Mathematical Sciences of the National Science Foundation, and by National Cancer Institute U54CA143833. J.A.S. was also supported by the Miller Foundation at University of California, Berkeley, and as an Einstein Visiting Fellow of the Einstein Foundation, Berlin. R.I.S. was also supported by an American Australian Association Sir Keith Murdoch Fellowship.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: sethian{at}math.berkeley.edu.

Author contributions: R.I.S. and J.A.S. designed research, performed research, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

↵

^{*}von Neumann’s 1951 paper (16) discussed the rate of change of area as a function of the number of sides due to gas diffusion, and because pressure is taken as constant, the boundaries become circular arcs. Mullins’ 1956 paper (17) studied the issue for metal grains, in which the movement of boundaries is governed by local conditions, and produced the same law, without the assumption of circular arcs.

## References

- ↵
- Sethian JA

- ↵
- ↵
- Smith KA,
- Solis FJ,
- Chopp DL

- ↵
- ↵
- ↵
- Sethian JA

- ↵
- ↵
- Sethian JA

- ↵
- ↵
- ↵
- ↵
- Sethian JA

- ↵
- ↵
- ↵
- Weaire D,
- Hutzler S

- ↵
- von Neumann J

- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

### More Articles of This Classification

### Physical Sciences

### Related Content

- No related articles found.