# Surface Evolver Documentation

Basic surface evolution is done by gradient descent iteration. Relevant topics are:

Each g command iteration does the following:
• Calculates the force vector at each vertex as the gradient of the total energy.
• Calculates the gradients at each vertex of constrained body volumes and named quantities. Also calculates the multiple of gradient motion needed to restore target values of the constraints.
• Projects the forces to be tangent to level set constraints, volume constraints (pressure is calculated here), and fixed quantity constraints.
• If any options are in effect, the proper force to velocity conversion is done.
• Does the conjugate gradient projection, if in effect.
• The current vertex positions are saved.
• If the optimizing mode is on, then trial motions are made to find the optimal scale factor.
• The scale factor is multiplied by the value of the scale_scale internal variable. (Useful only if you want to muck around with modifying the optimal scale factor for some strange reason.)
• Each vertex is moved by the scale factor times the velocity, plus the volume and fixed quantity restoring motions. If runge_kutta is in effect, then a fourth-order Runge-Kutta step is done instead of a simple Euler step.
• All level set constraints are enforced. Vertices violating an inequality or equality constraint are projected to the constraint (Newton's method). Several projection steps may be needed, until the violation is less that a certain tolerance or a certain number of steps are done (which generates a warning message). The default constraint tolerance is 1e-12, but it can be set with the CONSTRAINT_TOLERANCE option in the datafile, or setting the constraint_tolerance variable.
• If jiggling is on, the surface is randomly perturbed.
• If autopop or autochop are on, then the appropriate edge deletion or division is done.
• New energies, volumes, and quantities are calculated.
• If check_increase is on and the surface energy has increased, then the vertices are restored to their old values.
• If estimating is on, then a linear estimate of the energy change is printed. This is calculated from the velocities, gradients, and scale factor.
• The number of iterations left, the new area, energy, and scale factor are printed.
• If a graphics display is active and set to autodisplay, then graphics are redrawn.

## Scale factor

Once a direction of motion is found, the direction must be multiplied by a scale factor to compute the actual motion. If one interprets the direction of motion as a velocity (as in motion by mean curvature), then the scale factor becomes a time step. The scale factor may be fixed with the m command, or it may be in optimizing mode. The default is to start in optimizing mode.

## Optimizing scale

In using gradient descent to seek a minimum energy, one finds a direction of motion and does a line search along that direction to find the minimum energy. Evolver will do that in optimizing scale mode. The line search consists of halving or doubling the current scale factor until an energy minimum is bracketed; then quadratic interpolation is used to estimate the optimum scale. Optimizing scale is the default; it also may be turned on with the m command or the optimizing command.

For safety, there is an upper bound to the scale; it defaults to 1 but may be changed with the optimizing command. There is also a lower bound; if Evolver gets a scale below 1e-12 of the scale bound when attempting to find a minimum, it gives up and just uses scale 0. Scale 0 is not a null operation since it still projects to constraints, if they are not exactly satisfied.

In general, a good scale factor depends on the type of energy being minimized and the level of refinement. However, for minimizing area, when the triangulation is well-behaved and area normalization is off, the best scale factor is usually around 0.2, independent of refinement. In optimizing mode, a scale factor getting small, say below 0.01, indicates triangulation problems. Too large a fixed scale factor will show up as total energy increasing. If you have motion by area normalization ON use a small scale factor, like 0.001, until you get a feel for what works.

If check_increase is toggled on, then the motion is not done if it would increase energy. But be aware that energy sometimes may have to increase in order to satisfy constraints.

"Conjugate gradient" is a method of accelerating gradient descent. In ordinary gradient descent, one uses the gradient of energy to find the steepest downhill direction, then moves along that line to the minimum energy in that direction. Hence successive steps are at right angles. However, this can be very inefficient, as you can spend a lot of time zigzagging across an energy "valley" without making much progress "downstream". With conjugate gradient, the search direction is chosen to be in a "conjugate" direction to the previous direction. For a mathematical explanation, see any decent book in numerical analysis, such as [P]. In practice, the conjugate gradient method remembers a cumulative "history vector", which it combines with the ordinary gradient to figure out the conjugate gradient direction. The upshot is that conjugate gradient can converge much faster than ordinary gradient descent.

Conjugate gradient can be toggled with the U command, or with the conj_grad toggle. It should always be used with optimizing scale.

Notes: The conjugate gradient method is designed for quadratic energy functions. As long as the energy function is nearly quadratic, as it should be near an energy minimum, conjugate gradient works well. Otherwise, it may misbehave, either by taking too big steps or by getting stalled. Both effects are due to the history vector being misleading. To prevent too big steps, one should iterate without conjugate gradient for a few steps whenever significant changes are made to the surface (refining, changing a constraint, etc.). On the other hand, if it looks like conjugate gradient is converging, it may have simply become confused by its own history. See the catenoid example for a case in point. A danger signal is the scale factor going to zero. If you are suspicious, toggle conjugate gradient off and on ("U 2" does nicely) to erase the history vector and start over.

## Diffusion

The Evolver can simulate the real-life phenomenon of gas diffusion between neighboring bubbles. This diffusion is driven by the pressure difference across a surface. This is invoked by the keyword DIFFUSION in the first part of the datafile, followed by the value of the diffusion constant. The amount diffused across a facet during an iteration is calculated as scale*diffusion_constant*facet_area*pressure_difference. The scale factor is included as the time step of an iteration. The amount is added to or subtracted from the prescribed volumes of the bodies on either side of the facet.

## Stability

The timestep of an iteration should not be so large as to amplify perturbations of the surface. Short wavelength perturbations are most prone to amplification. This section contains a sketch of the stability characteristics of the various mobility modes, enough to let the user relate the maximum timestep to the minimum facet or edge size. Two examples are discussed: a zigzag string and a nearly flat surface with equilateral triangulation. Effective area is not included, as it is an insignificant correction for nearly flat surfaces. The general moral of this section is that the maximum time step in iteration is limited by the length of the shortest edge or the area of the smallest facet, except in one case.

### Zigzag string

Let the amplitude of the perturbation about the midline be Y and the edge length L. Then the force on a vertex is F = -4Y/L for small Y. Let the timestep (the Evolver scale factor) be \Delta t. Let V be the vertex velocity. Then the critical timestep for amplification of the perturbation is given by V\Delta t = -2Y, or \Delta t = -2Y/V.

Vertex mobility. Here V = F, so \Delta t = L/2.

Area normalization. Here the vertex star has length 2L, so V = F/L and \Delta t = L^2/2.

Approximate curvature. It turns out that the zigzag is an eigenvector of the mobility matrix M for the largest eigenvalue 3/L, so V = 3F/L and \Delta t = L^2/6. This is a major disadvantage of approximate curvature. If perturbation instability is the limitation on the timestep, it will take three times as many iterations as with area normalization to do the same evolution.

### Perturbed sheet with equilateral triangulation.

Consider a plane surface triangulated with equilateral triangles of area A. The perturbation consists of a tiling of hexagonal dimples with their centers of height Y above their peripheries. The force at a central vertex turns out to be -sqrt(3)Y and the force at a peripheral vertex sqrt(3)Y/2.

Vertex mobility. The critical time step is given by
(\sqrt(3)Y + sqrt(3)Y/2)\Delta t = 2Y,
so \Delta t = 4/3*sqrt(3). Note that this is independent of the triangle size. This is consistent with experience in evolving with optimizing scale factor, where the optimum time step is in the range 0.2 - 0.3 indenpendent of triangle size. This is a definite advantage of this version of mobility, since different parts of the surface can have different size triangulations and one size time step can work for all.

Area normalization. The star area of each vertex is 6A, so the velocities become -sqrt(3)Y/2A and sqrt(3)Y/4A, and the critical time step \Delta t = 4/6*sqrt(3)A. Hence on a surface with varying size triangles, the timestep is limited by the area of the smallest triangles.

Approximate area. This force turns out to be an eigenvector of the mobility matrix with eigenvalue 2/A. Hence the velocity is four times that of the area normalization, and the critical time step four times shorter.