Surface models.

The Surface Evolver can handle several different models of surfaces. In fact, there are several different categories of models regarding different features, so one variety of each category is active. However, some combinations are incompatible.

Incompatible models

The following combinations of surface models are incompatible:

Surface representation and combinatorics

All surfaces are simplicial complexes made of the basic elements: vertices, edges, and facets. The Evolver has three different ways of representing the combinatorics of the surface, depending on the dimension of the surface. Any of these may be used in any ambient space dimension at least as great at the surface dimension.

String model

The term "string model" means that the surface is one-dimensional. The datafile must have the keyword "STRING" or "SURFACE_DIMENSION 1" in its top section. Edges are defined in terms of their vertices, and facets by a list of boundary edges. Facets are not divided into triangles, and may have any number of edges. The edges of a facet need not form a closed loop, for example if the facet is partly bounded by a constraint. A body is defined by associating one facet to it, and the volume of the body is the area of the facet. The default energy is edge length.

Soapfilm model

The term "soapfilm model" means that the dimension of the surface is 2. This is the default model. The surface is subdivided into triangular facets, and the default energy is surface area. Edges are defined by their vertices. Facets are defined by an oriented list of three edges, which must form a closed loop. However, faces in the datafile may have more than three edges, since they are automatically refined into facets when loaded. In official Evolver-speak, a "face" is what appears in the datafile, and a "facet" is the triangle in the internal Evolver representation of the surface. Bodies are defined by a set of oriented facets, which need not form the complete boundary of the body, for example if part of the boundary is on a constraint.

Internally, the surface is held together by a set of structures called "facet-edges". There is one such structure for each incidence of an edge on a facet. There is a doubly linked list of facet-edges around each facet, so edges can be traversed in order, and there is a doubly-linked list around each edge, so the facets around an edge can be traversed in geometric order. Evolver figures out the geometric order from the geometric data in the datafile. If geometric order does not make sense, as when the space dimension is 4 or more, then the order is random.

Simplex model

The simplex model enables the representation of arbitrary dimension surfaces, but many Evolver features are not available with it. Here each facet is represented as an oriented list of k+1 vertices, where k is the dimension of the surface. Edges may be specified as k-1 dimensional simplices, but they are used only to compute constraint and named quantity integrals; a complete list of edges is not needed. Bodies are specified as lists of oriented facets.

The datafile must have the keyword SIMPLEX_REPRESENTATION in the top section, and the phrase 'SURFACE_DIMENSION k' if k isn't 2. k = 1 and k = 2 are allowed, but not very useful in comparison to the string or soapfilm models. If the domain is not 3-dimensional, then 'SPACE_DIMENSION n' must also be included. The datafile edges section is optional. Each facet should list k+1 vertex numbers. Non-simplicial facets are not allowed. See the sample datafile simplex3.fe.

Most features are not implemented. The quadratic model is not allowed, but the Lagrange model is. Vertices may be FIXED. Constraints are allowed, with energy integrands. Several basic named quantity methods work. No torus model or symmetry groups. No changing of surface topology or combinatorics is allowed except global refining with the r command. Refining subdivides each simplex 1-edge, with the edge midpoint inheriting the common attributes of the edge endpoints. Refining will increase the number of facets by a factor of 2^k.

Representation of surface elements

A surface is represented by a finite element method, where each element is a simplex. The simplest types of elements are just straight line segments and flat triangles. These combine to represent a piecewise linear surface, for which calculated quantities generally have an error of order h^2 compared to the ideal smooth surface, where h is the linear size of an element. Various ways of representing more accurate curved elements exist, and Evolver implements a couple of versions of what are called "Lagrange elements".

Linear model

In the linear model, all edges and triangular facets are flat line segments and triangles, respectively. For all calculations, an edge is defined by its two endpoints, and a facet (in the soapfilm model) is defined by its three vertices. This is the default. Quadratic or Lagrange models may be changed to linear with the M 1 or linear commands. An exception is if the spherical_arc_length method is used for length_method_name in the string model, in which case edges are computed and drawn on a sphere centered at the origin.

By default, edges and facets are linear. But it is possible to represent edges as quadratic curves and facets as quadratic patches by using the quadratic model. Each edge is endowed with a midpoint vertex. Most features are implemented for the quadratic representation, but some named quantity methods are not available, such as those involving curvature.

A special case is circular arc mode, which is effective in quadratic mode in the string model with the circular_arc_length method used for length_method_name. Then edges are calculated and drawn as exact circular arcs through their three vertices.

The smoothness of graphing of curved quadratic edges can be controlled with the internal variable `string_curve_tolerance`, which is the desired angle in degrees between successive graphed segments making up the edge.

The quadratic model can be invoked by putting the QUADRATIC keyword in the top section of the datafile or by using the quadratic or M 2 commands.

Lagrange model

The Evolver has a very limited implementation of higher-order elements. In the Lagrange model of order n, each edge is defined by interpolation on n+1 vertices evenly spaced in the parameter domain, and each facet is defined by interpolation on (n+1)(n+2)/2 vertices evenly spaced in a triangular pattern in the parameter domain. That is, the elements are Lagrange elements in the terminology of finite element analysis.

The Lagrange model is defined only for named quantities and methods, so Evolver will automatically do `convert_to_quantities` when you invoke the Lagrange model. The methods that currently accept the Lagrange model are

A surface may be converted to an order n Lagrange model with the command "lagrange n". This will convert linear or quadratic models to Lagrange, and will convert between different order Lagrange models. The commands linear and quadratic will convert Lagrange model back to the linear or quadratic models.

No triangulation manipulations are available in the Lagrange model. No refining, equiangulation, or anything. There is some vertex averaging, but just internal to edges and facets. Use the linear or quadratic model to establish your final triangulation, and just use the Lagrange model to get extra precision.

The current order can be accessed through the read-only internal variable lagrange_order. The Lagrange model can be dumped and reloaded.

As the Lagrange order is raised, calculations slow down rapidly. This is not only due to the large number of points involved, but is also due to the fact that the order of Gaussian integration is also raised.

Lagrange elements are normally plotted subdivided on their vertices, but if the smooth_graph flag is on, they are plotted with 8-fold subdivision.

The toggle command bezier_basis toggle replaces the Lagrange interpolation polynomials (which pass through the control points) with Bezier basis polynomials (which do not pass through interior control points, but have positive values, which guarantees the edge or facet is within the convex hull of the control points). This is experimental at the moment, and not all features such as graphing or refinement have been suitably adjusted.

Space dimension

By default, surfaces live in 3 dimensional space. However, the phrase "SPACE_DIMENSION n" in the datafile header sets the dimension to n. This means that all coordinates and vectors have n components. The only restriction is that Evolver has to be compiled with the MAXCOORD macro defined to be at least n in Makefile or in model.h. The default MAXCOORD is 4. Change MAXCOORD and recompile if you want more than four dimensions. The actual space dimension can be accessed in commands through the read-only variable space_dimension. Graphics will display only the first three dimensions of spaces with more than three dimensions, except for geomview, which has a four-dimensional viewer built in (although its use is awkward now).

Metric

For length and area to make sense, the ambient space must be endowed with a metric. The Evolver offers several choices, but keep in mind that they are only used to calculate default length and area. Other quantities that depend on the metric, such as volume, are up to the user to put in by hand with named quantities. All displaying is done as if the metric is Euclidean.

Euclidean metric

The default metric on the ambient space is the ordinary Euclidean metric. There are no built-in units of measurement like meters or grams, so the user should express all physical quantities in some consistent system of units, such as MKS or cgs.

Riemannian metric

The ambient space can be endowed with a general Riemannian metric by putting the keyword METRIC in the datafile followed by the elements of the metric tensor. Only one coordinate patch is allowed, but the quotient space feature makes this quite flexible. Edges and facets are linear in coordinates, they are not geodesic. The metric is used solely to calculate lengths and areas. It is not used for volume. To get a volume constraint on a body, you will have to define your own named quantity constraint. See quadm.fe for an example of a metric.

Conformal metric

The ambient space can be endowed with a conformal Riemannian metric by putting the keyword CONFORMAL_METRIC in the datafile followed by a formula for the conformal factor, i.e. the multiple of the identity matrix that gives the metric. Only one coordinate patch is allowed, but the quotient space feature makes this quite flexible. Edges and facets are linear in coordinates, they are not geodesic. The metric is used solely to calculate lengths and areas. It is not used for volume. To get a volume constraint on a body, you will have to define your own named quantity constraint. See quadm.fe for an example of a metric.

Klein hyperbolic metric

One special metric is available built-in. It is the Klein model of hyperbolic space in n dimensions. The domain is the unit disk or sphere in Euclidean coordinates. Including the keyword KLEIN_METRIC in the top section of the datafile will invoke this metric. Lengths and areas are calculated exactly, but as with other metrics you are on your own for volumes and quantities.

Symmetries

There are many interesting problems dealing with symmetric surfaces. A natural way to deal with a symmetric surface is to compute with only one fundamental domain of the symmetry, and use special boundary conditions. Some symmetries, such as mirror symmetries, can be handled with normal Evolver features. For example, a mirror can be implemented as a planar level set constraint. But symmetries such as translational or rotational symmetry require some built-in features. In any case, multiple copies of the fundamental domain may be displayed with the view_transforms command.

No symmetry

By default, the domain of a surface is Euclidean space. A symmetric surface can be done this way if its fundamental domain is bounded by mirror planes. Each mirror plane should be implemented as a linear level set constraint.

Torus model

The Evolver can take as its domain a flat torus with an arbitrary parallelpiped as its unit cell, i.e. the domain is a parallelpiped with its opposite faces identified. This is indicated by the TORUS keyword in the datafile. The defining basis vectors of the parallelpiped are given in the TORUS_PERIODS entry of the datafile. See twointor.fe for an example.

Vertex coordinates are given as Euclidean coordinates within the unit cell, not as linear combinations of the period vectors. The coordinates need not lie within the parallelpiped, as the exact shape of the unit cell is somewhat arbitrary.

The way the surface wraps around in the torus is given by saying how the edges cross the faces of the unit cell. In the datafile edges section, each edge has one symbol per dimension indicating how the edge vector crosses each identified pair of faces, and how the vector between the endpoints needs to be adjusted to get the true edge vector:

 * does not cross face + crosses in same direction as basis vector, so basis vector added to edge vector - crosses in opposite direction, so basis vector subtracted from edge vector.
Wraps are automatically maintained by the various triangulation manipulation operations.

There are several commands for ways of displaying a torus surface:
raw_cells - Graph the facets as they are, without clipping. The first vertex of a facet is used as the basepoint for any unwrapping of other vertices needed.
connected - Each body's facets are unwrapped in the torus, so the body appears in one connected piece. Nicest option, but won't show facets not on bodies.
clipped - Shows the unit cell specified in the datafile. Facets are clipped on the parallelpiped faces.

A few features are not available with the torus model, such as gravity and the simplex model. (Actually, you could put them in, but they will not take into account the torus model.)

Volumes and volume constraints are available. However, if the torus is completely partitioned into bodies of prescribed volume, then the volumes must add up to the volume of the unit cell and the TORUS_FILLED keyword must be given in the datafile. Or just don't prescribe the volume of one body.

Volumes are somewhat ambiguous. The volume calculation method is accurate only to one torus volume, so it is possible that a body whose volume is positive gets its volume calculated as negative. Evolver adjusts volumes after changes to be as continuous as possible with the previous volumes as possible, or with target volumes when available. You can also set a body's volconst attribute if you don't like the Evolver's actions.

Level set constraints can be used in the torus model, but be cautious when using them as mirror symmetry planes with volumes. The torus volume algorithm does not cope well with such partial surfaces. If you must, then use y=const symmetry planes rather than x=const, and use the -q option or do convert_to_quantities. Double-check that your volumes are turning out correctly; use volconst if necessary.

Display_periods: The displayed parallelogram unit cell can be different from the actual unit cell if you put an array called `display_periods` in the top of the datafile, in addition to the regular periods. For a string model example,

```   parameter shear = 1
torus_filled
periods
4  0
shear  4
display_periods
4 0
0 4 ```
This will always display a square, no matter how much the actual unit cell is sheared. This feature works well for shears; it may not work nicely for other kinds of deformation. Display_periods works better for the string model than the soapfilm model. For the soapfilm model, it seems to do horizontal shears best, but it can't cope with large shears, so if your shear gets too large, I advise resetting your fundamental region to less shear, say with the `unshear` command in `unshear.cmd`.

Display_origin: For a torus mode surface, if clipped mode is in effect, the center of the clip box is set with the `display_origin[]` array whose dimension is the dimension of the ambient space. This array does not exist by default, it has to be created by the user in the top of the datafile with the syntax

```  display_origin x y z
```
where x y z are the coordinates for the desired center of the clip box. At runtime, the array elements may be changed as normal:
```  display_origin[2] := 0.5
```
Changing display_origin will automatically cause the graphics to re-display.

Symmetry groups and quotient spaces

As a generalization of the torus model, you may declare the domain to be the quotient space of R^n with respect to some symmetry group. Several built-in groups are available, and ambitious users can compile C code into Evolver to define group operations. Group elements are represented by integers attached to edges (like the wrap specifications in the torus model at runtime). You define the integer representation of the group elements. See the file quotient.c for an example. See khyp.c and khyp.fe for a more intricate example modelling an octagon in Klein hyperbolic space identified into a genus 2 surface.

The datafile requires the keyword SYMMETRY_GROUP followed by the name for the group in quotes. Edges that wrap have their group element specified in the datafile by the phrase "wrap n", where n is the number of the group element. The wrap values are accessible at run time with the `wrap` attribute of edges. The group operations are accessible by the functions `wrap_inverse(w)` and `wrap_compose( w1,w2)`.

Using any Hessian commands with any symmetry group (other than the built-in torus model) will cause automatic converting to named quantities (with the " convert_to_quantities" command, since only named quantity Hessian evaluation routines have the proper symmetry transformation of the Hessian programmed in.

Volumes of bodies might not be calculated correctly with a symmetry group. The volume calculation only knows about the built-in torus model. For other symmetry groups, if you declare a body, it will use the Euclidean volume calculation. It is up to you to design an alternate volume calculation using named quantities and methods.

The currently implemented symmetry groups are:

`TORUS` symmetry group

This is the underlying symmetry for the torus model. Although the torus model has a number of special features built in to the Evolver, it can also be accessed through the general symmetry group interface. The torus group is the group on n-dimensional Euclidean space generated by n independent vectors, called the period vectors. The torus group uses the torus periods listed in the datafile.

Datafile declaration: `symmetry_group "torus"`

Group element encoding: The 32-bit code word is divided into 6-bit fields, one field for the wrap in each dimension, with low bits for the first dimension. Hence the maximum space dimension is 5. Within each bitfield, 1 codes for positive wrap and 011111 codes for negative wrap. The coding is actually a 2's complement 5-bit integer, so higher wraps could be represented.

`ROTATE` symmetry group

This is the cyclic symmetry group of rotations in the x-y plane, where the order of the group is given by the internal variable rotation_order (default value 4). There is also an internal variable generator_power (default 1) such that the angle of rotation is 2*pi*generator_power/rotation_order. Note: Since this group has fixed points, some special precautions are necessary. Vertices on the rotation axis must be labelled with the attribute axial_point in the datafile. Edges out of an axial point must have the axial point at their tail, and must have zero wrap. Facets including an axial point must have the axial point at the tail of the first edge in the facet.

Datafile declaration:

```  symmetry_group "rotate"
parameter rotation_order = 6 ```

Group element encoding: An element is encoded as the power of the group generator.

`FLIP_ROTATE` symmetry group

This is the cyclic symmetry group of rotations in the x-y plane with a flip z mapping to -z on every odd rotation, where the order of the group is given by the internal variable rotation_order, which had better be even. Note: Since this group has points that are fixed under an even number of rotations, some special precautions are necessary. Vertices on the rotation axis must be labelled with the attribute "`double_axial_point`" in the datafile. Edges out of an axial point must have the axial point at their tail, and must have zero wrap. Facets including an axial point must have the axial point at the tail of the first edge in the facet.

Datafile declaration:

```  symmetry_group "flip_rotate"
parameter rotation_order = 6 ```

Group element encoding: An element is encoded as the power of the group generator.

`CUBELAT` symmetry group

This is the full symmetry group of the unit cubic lattice.

Datafile declaration: `symmetry_group "cubelat"`

Group element encoding: wrap & {1,2,4} gives the sign changes for x,y,z respectively; (wrap&24)/8 is the power of the (xyz) permutation cycle; (wrap&32)/32 tells whether to then swap x,y. Thus wrap&63 is the same as for the cubocta symmetry group. Then wrap/64 is a torus wrap as in the torus symmetry group: three six-bit signed fields for translation in each coordinate. Translation is to be applied after other operations.

`CUBOCTA` symmetry group

This is the full symmetry group of the cube. It can be viewed as all permutations and sign changes of (x,y,z).

Datafile declaration: `symmetry_group "cubocta"`

Group element encoding: wrap & {1,2,4} gives the sign changes for x,y,z respectively; (wrap&24)/8 is the power of the (xyz) permutation cycle; (wrap&32)/32 tells whether to then swap x,y. (By John Sullivan; source in quotient.c under name pgcube)

`XYZ` symmetry group

The orientation-preserving subgroup of cubocta. Same representation.

`GENUS2` symmetry group

This is a symmetry group on the Klein model of hyperbolic space whose quotient group is a genus 2 hyperbolic manifold. The fundamental region is an octagon.

Datafile declaration:

```   Klein_metric
symmetry_group "genus2" ```
Group element encoding: There are 8 translation elements that translate the fundamental region to one of its neighbors. Translating around a vertex gives a circular string of the 8 elements. The group elements encoded are substrings of the 8, with null string being the identity. Encoding is 4 bits for each element. See khyp.fe for an example.

`DODECAHEDRON` symmetry group

This is the symmetry group of translations of hyperbolic 3 space tiled with right-angled dodecahedra. The elements of the group are represented as integers. There are 32 generators of the group so each generator is represented by five bits. Under this scheme any element that is the composition of up to five generators can be represented. If you want to use this group, you'll have to check out the source code in dodecgroup.c, since somebody else wrote this group and I don't feel like figuring it all out right now.

Datafile declaration:

```   Klein_metric
symmetry_group "dodecahedron" ```

`CENTRAL_SYMMETRY` symmetry group

This is the order 2 symmetry group of inversion through the origin, X maps to -X.

Datafile declaration:

`   symmetry_group "central_symmetry" `
Group element encoding: 0 for identity, 1 for inversion.

`SCREW_SYMMETRY` symmetry group

This is the symmetry group of screw motions along the z axis. The global parameter `screw_height` is the translation distance (default 1), and the global parameter `screw_angle` is the rotation angle in degrees (default 0).

Datafile declaration:

```  parameter screw_height = 4.0
parameter screw_angle  = 180.0
symmetry_group "screw_symmetry"
```
Group element encoding: the integer value is the power of the group generator.

quarter_turn

Symmetry group 3D torus with quarter turn in identification of top and bottom. x and y periods taken to be 1. z period is the user-defined variable quarter_turn_period. Generators x,y,z. x and y as in regular torus mode. z is vertical translation with quarter turn: (x,y,z)->(-y,x,z).
Relations: x z = z y^-1, y z = z x
Numerical representation: as in the torus symmetry group, for powers of x,y,z with generators applied in that order.

Mobility

There is a choice to be made in the conversion of the forces on vertices into velocities of vertices. Technically, force is the gradient of energy, hence a covector on the manifold of all possible configurations. In the Evolver representations of surfaces, that global covector can be represented as a covector at each vertex. The velocity is a global vector, which is represented as a vector at each vertex. Conversion from the global covector to the global vector requires multiplication by a metric tensor, i.e. singling out a particular inner product on global vectors and covectors. The tensor converting from force to velocity is the mobility tensor, represented as the mobility matrix M in some coordinate system. Its inverse, converting from velocity to force, is the resistance tensor S = M^{-1}. The same inner product has to be used in projecting the velocity tangent to the constraints, whether they be level set constraints on vertices or constraints on body volumes or quantity integrals. There are several choices implemented in the Evolver, corresponding to several different physical pictures of how the medium resists the motion of the surface through it:

Unit mobility

This is the default mobility, in which the velocity is equal to the force. Hence M and S are the identity matrices in standard coordinates. The physical interpretation of this is that there is a resistance to motion of each vertex through the medium proportional to its velocity, but not for the edges. This does not approximate motion by mean curvature, but it is very easy to calculate.

Area normalization

A type of mobility. In motion by mean curvature, the resistance to motion is really due to the surfaces, not the vertices. One way to approximate this is to say the resistance to motion of a vertex is proportional to the area associated with the vertex. So this scheme counts the area of a vertex equal to 1/3 of the area of the star of facets around it (or 1/2 the area of the star of edges in the string model). This is easy to calculate, since it is a local calculation for each vertex. S and M are diagonal matrices (see mobility). The result is that motion = force/area, which is a much better approximation to motion by mean curvature than the default of motion = force. Area normalization can be toggled with the a command or the area_normalizaton toggle.

Area normalization with effective area

A type of mobility. Simple area normalization as described in the previous paragraph isn't what's really wanted in certain circumstances. It has equal resistance for motion in all directions, both parallel and normal to the surface. If a vertex is a triple junction and migrating along the direction of one of the edges, it shouldn't matter how long that edge is. Therefore, if the effective area mode is in effect, the area associated with a vertex is the area of its star projected normal to the force at the vertex. This is a little more complicated calculation, but it is still local. S and M are block diagonal matrices, with one block for each vertex (see mobility). At a free edge not on any constraint, the force is tangent to the surface, the resistance is zero, and the mobility is infinite. But this accurately describes a popping soapfilm. Effective area can be toggled with the effective_area toggle. Note that area normalization itself must still be toggled with a or area_normalizaton.

Approximate polyhedral curvature

A type of mobility. Following a suggestion of Gerhard Dzuik and Alfred Schmidt, the inner product of global vectors is taken to be the integral of the scalar product of their linear interpolations over the facets (or edges in the string model). This has the advantage that the rate of area decrease of the surface is equal to the rate volume is swept out by the surface, which is a characteristic of motion by mean curvature. A big disadvantage is that the matrices M and S are no longer local (see mobility). S is a sparse matrix with entries corresponding to each pair of vertices joined by an edge, and M is its dense inverse. Approximate polyhedral curvature can be toggled with the approx_curv toggle command.

Approximate polyhedral curvature with effective area.

Polyhedral curvature does not make any distinction between motion parallel and perpendicular to the surface. A better approximation is to count only motion perpendicular to the surface. This can be done by projecting the interpolated vectorfields normal to the facets before integrating their scalar product. Now the rate of area decrease is equal to the rate geometric volume is swept out, as opposed to the slightly flaky way one had to calculate volume sweeping in the previous paragraph. Again S is a sparse matrix with entries corresponding to each pair of vertices joined by an edge, and M is its dense inverse. The effective area option may be toggled with effective_area.

User-defined mobility

The user may define a mobility tensor in the datafile. There is a scalar form with the keyword MOBILITY and a tensor form with MOBILITY_TENSOR. When in effect, this mobility is multiplied times the velocity to give a new velocity. This happens after any of the previous mobilities of this section have been applied and before projection to constraints. The formulas defining the mobility may include adjustable parameters, permitting the mobility to be adjusted during runtime. The formulas are evaluated at each vertex at each iteration, and so formulas may depend on vertex position and any vertex attributes.