Surface Evolver Manual

Version 1.92 July 31, 1993

Kenneth A. Brakke The Geometry Center

1300 South Second Street Minneapolis, MN 55454

and

Susquehanna University Selinsgrove, PA 17870 brakke@geom.umn.edu

Table of Contents

1. Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2. Installation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2.1. Acquisition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2. Environment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.3. Files. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.4. Unix Compilation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.5. Customizing graphics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.5.1. Random-order interface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.5.2. Painter interface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.5.3. X-Windows graphics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.5.4. Cheyenne graphics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3. Tutorial. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.1. Basic Concepts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.2. Example 1. Cube evolving into a sphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.3. Example 2. Mound with gravity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.4. Example 3. Catenoid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.5. Example 4. Rocket fuel tank. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.5.1. Surface energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.5.2. Volume. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.5.3. Gravity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.5.4. Running. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.6. Example 5. Spherical tank. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .18 3.6.1. Surface energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.6.2. Volume. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.6.3. Gravity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.6.4. Running. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.7. Example 6. Torus partitioned into two cells. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.8. Example 7. Crystalline integrand. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

4. The Model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4.1. Geometric elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1.1. Vertices. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1.2. Edges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1.3. Facets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1.4. Bodies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

4.2. Quadratic represention. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.3. Dimension of surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.4. Simplex model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.5. Dimension of ambient space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.6. Riemannian metric. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.7. Torus domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.8. Quotient spaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.9. Symmetric surfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.10. Constraints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.11. Boundaries. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .31 4.12. Energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.13. Named quantities and methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.14. Numbered surface integrals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.15. Pressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.16. Volume or content. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .35 4.17. Diffusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.18. Motion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .35 4.19. Mobility. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .36

4.19.1. Vertex mobility. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.19.2. Area normalization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.19.3. Area normalization with effective area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.19.4. Approximate polyhedral curvature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.19.5. Approximate polyhedral curvature with effective area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.19.6. User-defined mobility. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .37

4.20. Stability. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .37 4.20.1. Zigzag string. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .37 4.20.2. Perturbed sheet with equilateral triangulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

4.21. Topology changes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.22. Refinement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .38 4.23. Adjustable constants and parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.24. The String Model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5. The Datafile. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

5.1. Lexical format. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .40 5.1.1. Lines and line splicing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .40 5.1.2. Comments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.1.3. Including files. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.1.4. Macros. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.1.5. Case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.1.6. Whitespace. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.1.7. Identifiers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.1.8. Real numbers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.1.9. Keywords. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.1.10. Colors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.1.11. Expressions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

5.2. Datafile organization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.2.1. Definitions and options. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .43 5.2.2. Vertex list. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.2.3. Edge list. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.2.4. Face list. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.2.5. Bodies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5.2.6. Commands. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .54

6. Operation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.1. System command line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.2. Initialization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.3. Error handling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.4. Commands. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

6.4.1. General language syntax. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 6.4.2. Single-letter command summary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

6.5. Alphabetical single-letter command reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.6. Command and query language. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .63 6.7. Graphics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.8. Script examples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 6.9. Interrupts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.10. Graphics output file formats. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .77

6.10.1. Pixar and Minneview. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.10.2. Geomview. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.10.3. PostScript. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.10.4. Triangle file. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 6.10.5. SoftImage file. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

7. Technical Reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .80 7.1. Notation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 7.2. Surface representation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 7.3. Energies and forces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

7.3.1. Surface tension. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .81 7.3.2. Crystalline integrand. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 7.3.3. Gravity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 7.3.4. Constraint integrals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 7.3.5. Surface integrands. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 7.3.6. Gap areas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .82 7.3.7. Ideal gas compressibility. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .82 7.3.8. Prescribed pressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 7.3.9. Squared mean curvature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .83 7.3.10. Squared Gaussian curvature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

7.4. Named quantities and methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 7.4.1. Edge length. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 7.4.2. Facet area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7.4.3. Path integrals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7.4.4. Scalar surface integral. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7.4.5. Hooke energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7.4.6. Local Hooke energy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7.4.7. Integral of mean curvature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7.4.8. Integral of gaussian curvature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 7.4.9. Gravity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 7.4.10. Average crossing number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 7.4.11. Knot energies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

7.5. Volumes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .89 7.5.1. Default facet integral. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .89 7.5.2. Symmetric content facet integral. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 7.5.3. Edge content integrals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .89 7.5.4. Volume in torus domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

7.6. Quantity integrands. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 7.6.1. Quantity of a facet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 7.6.2. Quantity of an edge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

7.7. Constraint projection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .92 7.7.1. Projection of vertex to constraints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 7.7.2. Projection of vector onto constraint tangent space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

7.8. Volume and quantity constraints. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .92 7.8.1. Volume restoring motion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 7.8.2. Motion projection in gradient mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 7.8.3. Force projection in mean curvature mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

7.9. Iteration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 7.9.1. Fixed scale motion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 7.9.2. Optimizing scale motion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .94 7.9.3. Conjugate gradient mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .94

7.10. Equiangulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 7.11. Area normalization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 7.12. Hidden surfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 7.13. Extrapolation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 7.14. Curvature test. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 7.15. Annealing (jiggling). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 7.16. Long wavelength perturbations (long jiggling). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .96 7.17. Homothety. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 7.18. Popping non-minimal edges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 7.19. Popping non-minimal vertex cones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .96 7.20. Refining. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .97 7.21. Removing tiny edges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 7.22. Weeding small triangles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 7.23. Vertex averaging. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 7.24. Zooming in on vertex. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .98 7.25. Mobility and approximate curvature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 7.26. Torus domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

8. Helpful hints and notes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 9. Bugs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 10. Version history. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 11. Bibliography. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 12. Index. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 Figures.

Figure 1. Cube for example 1 Figure 2. Mound for example 2 Figure 3. Tank for example 4 Figure 4. Sphere for example 5

Chapter 1. Introduction.

The Surface Evolver is an interactive program for the study of surfaces shaped by surface tension and other energies. A surface is implemented as a simplicial complex, that is, a union of triangles. The user defines an initial surface in a datafile. The Evolver evolves the surface toward minimal energy by a gradient descent method. The evolution is meant to be a computer model of the process of evolution by mean curvature, which was studied in [B1] for surface tension energy in the context of varifolds and geometric measure theory. The energy in the Evolver can be a combination of surface tension, gravitational energy, squared mean curvature, user-defined surface integrals, or knot energies. The Evolver can handle arbitrary topology (as seen in real soap bubble clusters), volume constraints, boundary constraints, boundary contact angles, prescribed mean curvature, crystalline integrands, gravity, and constraints expressed as surface integrals. The surface can be in an ambient space of arbitrary dimension, which can have a Riemannian metric, and the ambient space can be a quotient space under a group action. The user can interactively modify the surface to change its properties or to keep the evolution well-behaved. The Evolver was written for one and two dimensional surfaces, but it can do higher dimensional surfaces with some restrictions on the features available. Graphical output is available as screen graphics and in several file formats, including PostScript.

The Surface Evolver program is freely available (see chapter 2) and is in use by a number of researchers. Some of the applications of the Evolver so far include modelling the shape of fuel in rocket tanks in low gravity [Te], calculating areas for the Opaque Cube Problem [B4], computing capillary surfaces in cubes [MH] and in exotic containers [C], simulating grain growth, studying grain boundaries pinned by inclusions, searching for partitions of space more efficient than Kelvin’s tetrakaidecahedra, modelling the shape of molten solder on microcircuits [RSB], studying polymer chain packing, and classifying minimal surface singularities.

The strength of the Surface Evolver program is in the breadth of problems it handles, rather than optimal treatment of some specific problem. It is under continuing development, and users are invited to suggest new features.

This manual contains full operational details and mathematical descriptions. A journal article description of the Evolver is appearing in [B2].

The Evolver is written in portable C and has been run on several systems: Sun, Iris, NeXT, Xenix, Cray, MS-DOS, Macintosh, and HP 98731. It is meant to be easily portable to any system that has C.

The Evolver was written as part of the Geometry Supercomputing Project (now The Geometry Center), sponsored by the National Science Foundation, the Department of Energy, Minnesota Technology, Inc., and the University of Minnesota. The program is in the public domain.

Bug reports should be submitted by email to brakke@geom.umn.edu. Please include the Evolver version number, a description of the problem, the initial data file, and the sequence of commands necessary to reproduce the problem.

Previous users of the Evolver should consult the History chapter for the new features.

Newsletter. The group of Surface Evolver users has grown large enough that I have started a newsletter. Contents include:

• Latest features added to Evolver.

• Future plans for new features.

• Suggestions for new features, with the likelihood of their being added proportional to demand.

• Bibliography of papers using Evolver, hopefully with brief abstracts.

• Workshops, talks, etc. to be given.

• Descriptions of projects underway using Evolver.

• List of systems users have ported Evolver to.

• Anything else anybody would wish to contribute. If you would like to be on the mailing list, send your

    email address to brakke@geom.umn.edu.

Chapter 2. Installation.

This chapter explains how to get and install the Evolver. Evolver is written to be portable between systems. The only requirement is a C compiler. Graphics on different systems may need customizing; see the section on graphics porting below.

§2.1. Acquisition.

The Evolver program is public domain, and there is no charge for it.

If you can connect to Internet and can use ftp, do a ftp login to geom.umn.edu (128.101.25.31) as “anonymous” with any password. From the pub directory, get evolver.tar.Z (970Kbytes) or evolver.tar (2600Kbytes), if you don’t have uncompress. This archive contains source code, datafiles, and this manual. See below for special versions for NeXTs, Macs, and PCs.

The program is distributed in a compressed tar file evolver.tar.Z. Get this file into a working directory. Uncompress it with

    uncompress evolver.tar

(yes, leave off the.Z). Extract the files with

    tar xvf evolver.tar

If you cannot use ftp, you can get Evolver on a floppy disk. Send email to brakke@geom.umn.edu or write to

Kenneth A. Brakke

Mathematics Department Susquehanna University Selinsgrove, PA 17870 USA

Indicate what format you can handle: 5 1/4” or 3 1/2” floppy, tar, compressed, MS-DOS, etc.

This complete manual is in TeX DVI format in manual.dvi. Print it out however you print DVI files. If you can’t handle it, you can ask for a hardcopy from brakke@geom.umn.edu. The figures to accompany the tutorial examples in the manual are in PostScript files *.ps, which can be sent direct to any PostScript printer.

For those with NeXT computers, there is also a separate ftp archive pub/evolver.next.tar.Z which contains a NeXT executable Evolver and Interface Builder files. This version will only run from a NeXT console. Even though there is a command line option (-u) to not display graphics, it still insists on access to the screen. If you want a version you can run remotely, get the source code and compile vanilla version with no graphics or X windows graphics.

For those with Macintoshes, there is a Mac version available for anonymous ftp from geom.umn.edu as pub/evolver.sit.hqx. This is a binhexed Stuffit 1.5.1 archive. It includes a README file with Mac specific information and datafiles in Mac format.

For those with PC compatibles, there is a 32-bit DOS version in pub/evolver.zip. This PKZIP archive contains evolver.exe and the DOS extender program needed to run it, and a few datafiles. It was made with the Free Software Foundation’s DJGPP compiler.

Screen graphics generally uses an external graphics program. There is an X-Windows graphics interface, and the NeXT version uses NeXT graphics. The best graphics, however, are done on Silicon Graphics systems with the geomview graphics program, which is The Geometry Center’s successor to its MinneView program.

 geomview is available by anonymous ftp from geom.umn.edu. (The library pub/evolverOOGLlibs.tar.Z

is needed only if you are compiling a Minneview version of Evolver, which is obsolete).

Obsolete: The Cheyenne graphics package for Sun workstations is available via anonymous ftp from princeton.princeton.edu (128.112.128.1 or 128.112.32.1), as a compressed tar file pub/graphics/cy5.3.tar.Z. Installation instructions should be included.

§2.2. Environment.

Evolver needs to find the initial data file and sometimes other data files (Wulff vectors, color maps). If the named file is not in the current directory, then an environment variable called EVOLVERPATH will be consulted for a directory search list. The format is the same as the usual PATH variable. Set it up as usual in your system:

Unix C shell:

setenv EVOLVERPATH /usr/you/data:/usr/you

Bourne shell:

EVOLVERPATH=/usr/you/data:/usr/you
export EVOLVERPATH

MS-DOS:

set EVOLVERPATH=\evolver\data;d:\misc

For the MS-DOS Turbo C version, the graphics driver files CGA.BGI or EGAVGA.BGI should be in the current directory, or in the path given by an environment variable BGIPATH.

Output files are written to the filename given by the user in response to a prompt, so include the path you want.

§2.3. Files.

You will probably want to put the different groups of files in different directories.

Documentation:

manual.dvi
*.ps

For datafile and command parsing (user should not remake these):

datafile.lex lex source for lexyy.c
command.yac yacc source for ytab.c

For compilation:

 Makefile
 *.c
 *.h

Screen graphics interfaces:

xgraph.c  X-Windows graphics interface
geomgraph.c              geomview interface for Silicon graphics
mvgraph.c                MinneView interface for Silicon Graphics (obsolete)
iriszgraph.c  Silicon Graphics GL graphics
irisstereo.c  does stereo displays on Irises so equipped
hpzgraph.c  HP 98731 graphics

tcgraph.c  Turbo C graphics
gnugraph.c  graphics for DOS DJGPP version
mscgraph.c  MicroSoft C graphics
cheygraph.c  Cheyenne graphics for Suns
nulgraph.c  no screen graphics

Evolver data files:

*.fe (initial data files)
*.wlf (Wulff vector files)
jcmap (sample color map file)

You may put the data and colormap files in various directories. Evolver checks the environment variable EVOLVERPATH for data file paths.

§2.4. Unix Compilation.

First, you need to modify Makefile for your system. Makefile begins with sets of macro definitions for various systems. If your system is listed, remove the comment symbols ’#’ from start of your definitions. If your system is not there, use the GENERIC defines, or set up your own, and leave the graphics file as nulgraph. If you do define your own, be sure to put a corresponding section in include.h. After you get a working program you can write a screen graphics interface if you wish. Edit CFLAGS to have the proper options (optimization, floating point option, etc.). GRAPH should be the name of a screen graphics interface file (see list under Files above, or write your own as discussed below). GRAPHLIB should be the appropriate graphics library plus any other libraries needed. Check all paths and change if necessary for your system.

Silicon Graphics systems with multiple processors may compile a version that will use all processors for some calculations by including -DSGI_MULTI in CFLAGS. This version will run fine with one processor, also. Currently, the only calculations done in parallel are the “named quantities”.

The file include.h lists the include files for various systems. If your system isn’t listed, you will have to put in a list of your own. Try copying the generic list (or one of the others) and compiling. Your compiler will no doubt be very friendly and helpful in pointing out unfound header files.

include.h also has various other system-specific defines. See the GENERIC section of it for comments on these, and include the appropriate ones for your system.

Now try compiling. Hopefully, you will only have to change the system-specific parts of Makefile and include.h to get things to work. If significant changes to other files are needed, let me know. For example, the return value of sprintf() varies among systems, and that caused problems once upon a time.

Test by running on the cube sample file as described in the Tutorial chapter.

§2.5. Customizing graphics.

This section is provided for those people who need to write their own graphics interface module. All device-specific graphics output has been collected into a few basic routines. There are two styles of interface. One provides facets in random order for routines that can do their own hidden surface removal. The other provides facets in back to front order (painter algorithm).

§2.5.1. Random-order interface.

The random-order interface has several global pointers to functions which should be set to the user’s functions. This permits several sets of graphics routines in one program. The facets are generated by graphgen(), which call the user functions. The function display() is called as the overall display function; it should set the function pointers and call graphgen(). Example: iriszgraph.c

The function pointers are:

     void (*graph_start)(void)

This is called at the start of each display. It should do whatever device initialization and screen clearing is needed.

void (*graph_facet)(struct graphdata *,facet_id)

This is called to graph one triangular facet. See extern.h for the definition of graphdata. The second argument is a facet identifier for routines that are not happy with just the information in graphdata; it is for Evolver gurus only. Facets are presented in random order. The coordinates are not transformed; the current transformation matrix is in view[][], which is in homogeneous coordinates.

void (*graph_edge)(struct graphdata *)

This function is called to graph one edge in the string model.

void (*graph_end)(void)

This function is called after all data has been given.

§2.5.2. Painter interface.

The painter model is similar, except there is an extra layer of functions. Examples: cheygraph.c, tcgraph.c . In display(), the user should set

graph_start = painter_start;
graph_facet = painter_facet;
graph_end   = painter_end;

The user should also set these function pointers:

void (*init_graphics)(void)

Called by painter_start() to do device initialization.

void (*display_facet)(struct tsort *)

Called by painter_end() to graph sorted facets one at a time. See extern.h for struct

tsort.
void (*finish_graphics)(void)

Called at the end of all facets.

The user should set graph_edge as in the random interface.

§2.5.3. X-Windows graphics.

To have Evolver use X-Windows graphics for its native graphics, edit Makefile to have xgraph as the graphics module and have all the proper libraries linked in. The X window will appear when you do the “s” command. Evolver does not continually poll for window redraw events, so if you move or resize the window you will have to do “s” again to redraw the surface.

§2.5.4. Cheyenne graphics.

The graphics included for Sun workstations uses a public domain package known as Cheyenne. See section §2.1 for information. You will want to modify the file cheygraph.c to ensure that cy_open_device() gets the proper device and machine names. See the PORTING NOTES in that file. Consult Cheyenne documentation for system specific details. Note: Cheyenne graphics is obsolete. Use the X-Windows interface xgraph.c instead if possible.

Chapter 3. Tutorial.

This chapter introduces some basic concepts and then guides the user through a series of examples. This chapter can be read before the following chapters. The following chapters assume you have read the Basic Concepts section of this chapter.

§3.1. Basic Concepts.

The basic geometric elements used to represent a surface are vertices, edges, facets, and bodies. Vertices are points in Euclidean 3-space. Edges are straight line segments joining pairs of vertices. Facets are flat triangles bounded by three edges. A surface is a union of facets. (Actually, there is no separate surface entity in the program; all facets essentially belong to one logical surface.) A body is defined by giving its bounding facets.

The term “surface”, when used to refer to the entity upon which the Evolver operates, refers to all the geometric elements plus auxiliary data such as constraints, boundaries, and forces.

There are no limitations on how many edges may share a vertex nor on how many facets may share an edge. Thus arbitrary topologies are possible, including the triple junctions of surfaces characteristic of soap films.

Edges and facets are oriented for bookkeeping purposes, but there are no restrictions on the orientation of neighboring facets. Unoriented surfaces are possible.

A surface is deemed to have a total energy, arising from surface tension, gravitational energy, and possibly other sources. It is this energy which the Evolver minimizes.

No particular units of measurement are used. The program only deals with numerical values. If you wish to relate the program values to the real world, then all values should be within one consistent system, such as cgs or mks.

The initial surface is specified in a text file (hereafter referred to as the datafile) that may be created with any standard text editor. (The .fe extension I always use for datafiles stands for facet-edge, which refers to the internal data structure used to represent the surface. You may use any name you wish for a datafile.)

The basic operation of the Evolver is to read in a datafile and take commands from the user. The main command prompt is

    Enter command:

Commands are one letter (case is significant), sometimes with a numerical parameter. The most frequently used commands are:

g n do n iterations
s  show surface on screen
r  refine triangulation of surface
P  graphics output
q  quit

Commands must be followed with the ENTER key; Evolver only reads complete lines.

An iteration is one evolution step. The motion for the step is calculated as follows: First, the force on each vertex is calculated from the gradient of the total energy of the surface as a function of the position of that vertex. The force gives the direction of motion. Second, the force is made to conform to whatever constraints are applicable. Third, the actual motion is found by multiplying the force by a global scale factor.

§3.2. Example 1. Cube evolving into a sphere.

A sample datafile cube.fe comes with Evolver. The initial surface is a unit cube. The surface bounds one body, and the body is constrained to have volume 1. There is no gravity or any other force besides surface tension. Hence the minimal energy surface will turn out to be a sphere. This example illustrates the basic datafile format and some basic commands.

Let’s look at the datafile that specifies the initial unit cube (see Figure 1 for a sketch):

// cube.fe
// Evolver data for cube of prescribed volume.
vertices  /* given by coordinates */
1  0.0 0.0 0.0
2  1.0 0.0 0.0
3  1.0 1.0 0.0
4  0.0 1.0 0.0
5  0.0 0.0 1.0
6  1.0 0.0 1.0
7  1.0 1.0 1.0
8  0.0 1.0 1.0
edges  /* given by endpoints */
1   1 2
2   2 3
3   3 4
4   4 1
5   5 6
6   6 7
7   7 8
8   8 5
9   1 5
10  2 6
11  3 7
12  4 8
faces  /* given by oriented edge loop */
1   1 10 -5  -9
2   2 11 -6 -10
3   3 12 -7 -11
4   4  9 -8 -12
5   5  6  7   8
6  -4 -3 -2  -1
bodies  /* one body, defined by its oriented faces */
1   1 2 3 4 5 6  volume 1

The datafile is organized in lines, with one geometric element defined per line. Vertices must be defined first, then edges, then faces, then bodies. Each element is numbered for later reference in the datafile.

Comments are delimited by /* */ as in C, or follow // until the end of the line as in C++. Case is not significant, and all input is lower-cased immediately. Hence error messages about your datafiles will refer to items in lower case, even when you typed them in upper case.

The datafile syntax is based on keywords. The keywords VERTICES, EDGES, FACES, and BODIES signal the start of the respective sections. Note that the faces are not necessarily triangles (which is why they are called FACES and not FACETS). Any non-triangular face will be automatically triangulated by putting a vertex at its center and putting in edges to each of the original vertices. Faces don’t have to be planar. Note that a minus sign on an edge means that the edge is traversed in the opposite direction from that defined for it in the EDGES section. The face oriented normal is defined by the usual right hand rule. The cube faces all have outward normals, so they all are positive in the body list. In defining a body, the boundary faces must have outward normals. If a face as defined has an inward normal, it must be listed with a minus sign.

That the body is constrained to have a volume of 1 is indicated by the keyword VOLUME after the body definition, with the value of the volume following. In general, any attributes or properties an element has are given on the same line after its definition.

Start Evolver and load the datafile with the command line

evolver cube.fe

You should get a prompt

Enter command:

If you have geomview, do command P option 8 to get a display (option 6 if you have MinneView). It takes a couple of seconds to initialize. You can manipulate the geomview display as usual independently of the Evolver. Evolver will automatically update the image whenever the surface changes.

If you do not have geomview or MinneView and if you have screen graphics linked in, give the command s to show the surface. You should see a square divided into four triangles by diagonals. This is the front side of the cube; hidden surfaces are removed by default. The graphics command prompt is

Graphics command:

It takes strings of letters, each letter making a viewing transformation on the surface: The most used ones are

r rotate right by 6 degrees
l rotate left by 6 degrees
u rotate up by 6 degrees
d rotate down by 6 degrees
R reset to original position
q quit back to main command prompt

Try rrdd to get an oblique view of the cube. Any transformations you make will remain in effect the next time you show the surface. Now do q to get back to the main prompt.

Now do some iterations. Give the command g 10 to do 10 iterations. You should get this:

10.  area:  4.786234596177497 energy:  4.786234596177497  scale:  0.295137
9.  area:  4.987767806940851 energy:  4.987767806940851  scale:  0.380324
8.  area:  5.101912421796101 energy:  5.101912421796101  scale:  0.293334
7.  area:  5.108866386444287 energy:  5.108866386444287  scale:  0.283689
6.  area:  5.112087629027262 energy:  5.112087629027262  scale:  0.280605
5.  area:  5.112431292155377 energy:  5.112431292155377  scale:  0.280773
4.  area:  5.112484710716689 energy:  5.112484710716689  scale:  0.280554
3.  area:  5.112491934865932 energy:  5.112491934865932  scale:  0.280616
2.  area:  5.112492961234981 energy:  5.112492961234981  scale:  0.28059
1.  area:  5.112493104356962 energy:  5.112493104356962  scale:  0.280599

Note that after each iteration a line is printed with the iterations left, area, energy, and current scale factor. By default, the Evolver seeks the optimal scale factor to minimize energy. At first, there are large motions, and the volume constraint may not be exactly satisfied. The energy increase is due to the volume constraint taking hold. Volume constraints are not exactly enforced, but each iteration tries to bring the volume closer to the target. Here that results in increases in area. You can find the current volumes with the v command:

Body       target volume         actual volume        pressure
  1     1.000000000000000     0.999999779366360   3.408026016427987

Now let’s refine the triangulation with the r command. This subdivides each facet into four smaller similar facets. The printout here gives the counts of the geometric elements and the memory they take:

Vertices:  50 Edges:  144 Facets:  96  Facetedges:  288 Memory:  27554

Iterate another 10 times:

10.  area:  4.908899804670224 energy:  4.908899804670224  scale:  0.268161
9.  area:  4.909526310166165 energy:  4.909526310166165  scale:  0.204016
8.  area:  4.909119925577212 energy:  4.909119925577212  scale:  0.286541
7.  area:  4.908360229118204 energy:  4.908360229118204  scale:  0.304668
6.  area:  4.907421919968726 energy:  4.907421919968726  scale:  0.373881
5.  area:  4.906763705259419 energy:  4.906763705259419  scale:  0.261395
4.  area:  4.906032256943935 energy:  4.906032256943935  scale:  0.46086
3.  area:  4.905484754688263 energy:  4.905484754688263  scale:  0.238871
2.  area:  4.904915540917190 energy:  4.904915540917190  scale:  0.545873
1.  area:  4.904475138593070 energy:  4.904475138593070  scale:  0.227156

You can continue iterating and refining as long as you have time and memory.

Eventually, you will want to quit. So give the q command. You get

Enter new datafile name (none to continue, q to quit):

You can start a new surface by entering a datafile name (it can be the same one you just did, to start over), or continue with the present surface by hitting ENTER with no name (in case you pressed q by accident, or suddenly you remember something you didn’t do), or you can really quit with another q.

§3.3. Example 2. Mound with gravity.

This example is a mound of liquid sitting on a tabletop with gravity acting on it. The contact angle between the drop surface and the tabletop is adjustable, to simulate the different degrees to which the liquid wets the table. This example illustrates macros, adjustable constants, constraints with energy, and omitting faces from body surfaces.

The drop starts as a cube with one face (face 6 of example 1) on the tabletop (the z = 0 plane). The most straightforward way to specify a contact angle is to declare face 6 to be constrained to stay on the tabletop and give it a surface tension different than the default of 1. But this leads to problems described below. The way the contact angle is handled instead is to omit face 6 and give the edges around face 6 an energy integrand that results in the same energy we would get if we did include face 6. If we let the interface energy density for face 6 be T, then we want a vectorfield w⃗ such that

∫ ∫_{face 6} T k⃗ · d⃗S = ∫_{∂(face 6)} w⃗ · d⃗l.

So all we need is curl w⃗ = T k⃗, and I have used w⃗ = −T yi⃗.

I have chosen to parameterize the contact angle as the angle in degrees between the table and the surface on the interior of the drop. This angle can be adjusted with the A command at run time. I could have made T the parameter directly, but then I wouldn’t have had an excuse to show a macro.

Here is the datafile mound.fe (see figure 2 for a sketch):

// mound.fe

// Evolver data for drop of prescribed volume sitting on plane with gravity.
// Contact angle with plane can be varied interactively with A command.
PARAMETER angle = 90    // interior angle between plane and surface, degrees
#define T  (-cos(angle*pi/180))  // virtual tension of facet on plane
constraint 1   /* the table top */
formula:  z = 0
energy:    // for contact angle
e1:  -T*y
e2:  0
e3:  0
vertices
1   0.0  0.0 0.0  constraint 1  /* 4 vertices on plane */
2   1.0  0.0 0.0  constraint 1
3   1.0  1.0 0.0  constraint 1
4   0.0  1.0 0.0  constraint 1
5   0.0  0.0 1.0
6   1.0  0.0 1.0
7   1.0  1.0 1.0
8   0.0  1.0 1.0
9   2.0  2.0 0.0  fixed   /* for table top, display purposes only */
10  2.0 -1.0 0.0  fixed
11 -1.0 -1.0 0.0  fixed
12 -1.0  2.0 0.0  fixed
edges  /* given by endpoints and attribute */
1   1 2    constraint 1 /* 4 edges on plane */
2   2 3    constraint 1
3   3 4    constraint 1
4   4 1    constraint 1
5   5 6
6   6 7
7   7 8
8   8 5
9   1 5
10  2 6
11  3 7
12  4 8
13  9 10   fixed  /* for table top */
14 10 11   fixed
15 11 12   fixed
16 12  9   fixed
faces  /* given by oriented edge loop */
1   1 10 -5  -9
2   2 11 -6 -10
3   3 12 -7 -11
4   4  9 -8 -12
5   5  6  7   8
7  13 14 15  16  density 0 fixed /* table top for display */

    bodies  /* one body, defined by its oriented faces */
    1   1 2 3 4 5   volume 1  density 1

The mound itself was basically copied from cube.fe, but with face 6 deleted. The reason for this is that face 6 is not needed, and would actually get in the way. It is not needed for the volume calculation since it would always be at z = 0 and thus not contribute to the surface integral for volume. The bottom edges of the side faces are constrained to lie in the plane z = 0, so face 6 is not needed to keep them from catastrophically shrivelling up. We could have handled the contact angle by including face 6 with a surface tension equal to the interface energy density between the liquid and surface, but that can cause problems if the edges around face 6 try to migrate inward. After refinement a couple of times, interior vertices of face 6 have no forces acting on them, so they don’t move. Hence it would be tough for face 6 to shrink when its outer vertices ran up against its inner vertices. The tabletop face, face 7, is entirely extraneous to the calculations. Its only purpose is to make a nice display. You could remove it and all its vertices and edges without affecting the shape of the mound. It’s constraint 1 that is the tabletop as far as the mound is concerned.

Now run Evolver on mound.fe. Refine and iterate a while. You should get a nice mound. It’s not a hemisphere, since gravity is on by default with G = 1. If you use the G command to set G = 0 and iterate a while, you get a hemisphere. Try using the A command to change the contact angle, to 45 degrees or 135 degrees for example. You can also play with gravity. Set G = 10 to get a flattened drop, or G = −5 to get a drop hanging from the ceiling. G = −10 will cause the drop to break loose.

§3.4. Example 3. Catenoid.

The catenoid is the surface formed between two rings not too far apart. In cylindrical coordinates, its equation is r = (1/a) cosh(az). In cat.fe, both the upper and lower rings are given as one-parameter boundary wires. The separation and radius are adjustable constants, so you can change them during a run with the A command. The initial radius given is the minimum for which a catenoid can exist for the given separation of the rings. To get a stable catenoid, you will have to increase this value. However, if you do run with the original value, you can watch the neck pinch out. The initial surface consists of six rectangles forming a cylinder between the two circles. The vertices on the boundaries are fixed, elsewise they would slide along the boundary to short-cut the curvature; two diameters is shorter than one circumference. The boundary edges are fixed so that vertices arising from subdividing the edges are likewise fixed.

Here is the catenoid datafile:

    // cat.fe
    // Evolver data for catenoid.
    PARAMETER  RMAX = 1.5088795   // ring radius; here critical value
    PARAMETER  ZMAX = 1.0         // ring distance from midplane
    boundary 1 parameters 1     //  upper ring
    x1:    RMAX * cos(P1)
    x2:    RMAX * sin(P1)
    x3:    ZMAX
    boundary 2 parameters 1    //   lower ring
    x1:    RMAX * cos(P1)
    x2:    RMAX * sin(P1)
    x3:    -ZMAX

vertices   // given in terms of boundary parameter
1    0.00  boundary 1   fixed
2    pi/3  boundary 1   fixed
3  2*pi/3  boundary 1   fixed
4    pi    boundary 1   fixed
5  4*pi/3  boundary 1   fixed
6  5*pi/3  boundary 1   fixed
7    0.00  boundary 2   fixed
8    pi/3  boundary 2   fixed
9  2*pi/3  boundary 2   fixed
10   pi    boundary 2   fixed
11 4*pi/3  boundary 2   fixed
12 5*pi/3  boundary 2   fixed
edges  /* given by endpoints */
1    1  2  boundary 1   fixed
2    2  3  boundary 1   fixed
3    3  4  boundary 1   fixed
4    4  5  boundary 1   fixed
5    5  6  boundary 1   fixed
6    6  1  boundary 1   fixed
7    7  8  boundary 2   fixed
8    8  9  boundary 2   fixed
9    9  10 boundary 2   fixed
10   10 11 boundary 2   fixed
11   11 12 boundary 2   fixed
12   12 7  boundary 2   fixed
13   1  7
14   2  8
15   3  9
16   4  10
17   5  11
18   6  12
faces  /* given by oriented edge loops */
1   1 14 -7 -13
2   2 15 -8 -14
3   3 16 -9 -15
4   4 17 -10 -16
5   5 18 -11 -17
6   6 13 -12 -18

The parameter in a boundary definition is always P1 (and P2 in a two-parameter boundary). The Evolver can handle periodic parameterizations, as is done in this example.

Try this sequence of commands (displaying at your convenience):

r  (refine to get a crude, but workable, triangulation)
u  (equiangulation makes much better triangulation)
g 120 (takes this many iterations for neck to collapse)
t  (collapse neck to single vertex by eliminating all edges shorter than 0.05)

0.05

o  (split neck vertex to separate top and bottom surfaces)

g  (spikes collapse)

The catenoid shows some of the subtleties of evolution. Suppose the initial radius is set to RMAX = 1.0 and the initial height to ZMAX = 0.55. Fifty iterations with optimizing scale factor result in an area of 6.458483. At this point, each iteration is reducing the area by only .0000001, the triangles are all nearly equilateral, everything looks nice, and the innocent user might conclude the surface is very near its minimum. But this is really a saddle point of energy. Further iteration shows that the area change per iteration bottoms out about iteration 70, and by iteration 300 the area is down to 6.4336. The triangulation really wants to twist around so that there are edges following the lines of curvature, which are vertical meridians and horizontal circles. Hence the optimum triangulation appears to be rectangles with diagonals.

This saddle point also seems to confuse the conjugate gradient method (§7.8.3). With conjugate gradient in effect (use the U command to toggle on and off), the saddle point is reached at iteration 17 and area starts decreasing again until iteration 30, when it reaches 6.4486. But then iteration stalls out, and the conjugate gradient mode has to be turned off and on to erase the history vector. Once restarted, another 20 iterations will get the area down to 6.4334. This shows that conjugate gradient mode can work much better than ordinary mode, but it can also have problems.

Exercise for the reader: Get the Surface Evolver to display an unstable catenoid by declaring the catenoid facets to be the boundary of a body, and adjusting the body volume with the b command to get zero pressure.

§3.5. Example 4. Rocket fuel tank.

This example discusses a cylindrical flat-ended tank in low gravity with some liquid inside. Complete understanding of this example qualifies you as an Evolver Power User.

The initial configuration is shown in figure 3. The axis is horizontal with the base at the left. We will assume the fuel is in a ring around the base of the tank. The fuel comprises one body, and the boundary of the body has four components:

1. The covered part of the flat circular base at y = 0, implemented as constraint 1..

2. The cylindrical wall at x^{2} + z^{2} = R^{2}, implemented as constraint 2..

3. The triangular facets (8 of them in figure 3).

4. The gaps between the edges on the wall (the “gap edges”, edges 5,6,7,8 in figure 3) and the wall. We will assume that the gap surfaces are generated by lines radially outward from the central axis (the y axis) through the points of the gap edges. To specify contact angles, let θ_{B} be the interior contact angle with the base, and let θ_{W} be the interior

contact angle with the cylindrical wall.

§3.5.1. Surface energy.

Base.

The surface energy here comes from the contact angle. The angle is represented by pretending the uncovered base has surface tension T_{B} = cos(θ_{B}). That we handle with an energy integral on constraint 1. We want

E = ∫ ∫_{base} T_{B}j⃗ · d⃗A.  (1)

By Stokes’ theorem, this is the integral around the edges in the base (edges 1,2,3,4)

E = ∫_{edges} T_{B}zi⃗ · d⃗s,  (2)

where the edges are oriented oppositely from figure 3. Hence the energy integrand for constraint 1 should be −T_{B}zi⃗.

Wall.

The contact angle is represented by pretending the uncovered wall has surface tension T_{W} = cos(θ_{W}). The uncovered wall in infinite, so we will assume equivalently that the covered wall has tension −T_{W}. The energy is

E = ∫ ∫_{wall} −T_{W} dA = ∫ ∫_{wall} −T_{W} \frac{1}{R}(xi⃗ + zk⃗) · d⃗A.  (3)

We want to turn this into an integral over the boundary of the wall, which is made of two paths, the base circumference and the “gap trace” where the gap surfaces meet the wall. This calls for a vector potential, but our integrand above does not have zero divergence. But we could choose a divergenceless integrand that is equal on the wall:

E = ∫ ∫_{wall} F⃗ · d⃗A.  (4)

where

F⃗ = −T_{W} \frac{R}{(x^{2} + z^{2})}(xi⃗ + zk⃗).  (5)

Now we can get a vector potential w⃗

w⃗ = −T_{W} \frac{Ry}{(x^{2} + z^{2})}(−zi⃗ + xk⃗), ∇ × w⃗ = F⃗ .  (6)

Hence

E = ∫_{basecircum} w⃗ · ds + ∫_{gap trace} w⃗ · ds.  (7)

The first of these is zero since y = 0 there. The second is equal to the integral over the gap edges, since ∇ × w⃗ is parallel to the gap surfaces! So

E = ∫_{gap edges} w⃗ · ds,  (8)

where the gap edges are oppositely oriented from figure 3. Hence the constraint 2 energy integrand includes −w⃗.

The facets.

The surface tension energy here is automatically taken care of by the Evolver. We will leave the surface tension as the default value of 1.

The gaps.

We will declare constraint 2 to be CONVEX so that Evolver will include a gap surface energy for all the gap edges. We put the spring constant to be 1, as that value give the best approximation to the true area for small gaps.

§3.5.2. Volume.

The body volume V = ∫ ∫ ∫ dV can be calculated as a surface integral V = ∫ ∫ F⃗ · d⃗A where ∇ · F = 1. The integral over the facets is handled automatically, but there are two choices for F⃗. We will take F⃗ = (xi⃗ + yj⃗ + zk⃗)/3 by putting SYMMETRIC_CONTENT in the datafile. (This choice works better for the tank on its side; for an upright tank, the default F⃗ = zk⃗ is easier.) The surface integral has four components:

Base.

This part of the integral is zero, since y = 0 and the normal is j⃗.

The wall.

Here the j⃗ component of F⃗ can be dropped, and the rest put in the equivalent divergenceless form

F⃗_{W} = \frac{R^{2}}{3(x^{2} + z^{2})}(xi⃗ + zk⃗).  (9)

This has vector potential

w⃗ = \frac{R^{2}y}{3(x^{2} + z^{2})}(−zi⃗ + xk⃗),  ∇ × w⃗ = F⃗_{W} .  (10)

By the same reasoning as for wall area, this reduces to the integral of −w⃗ over the gap edges. Thus we include −w⃗ in the content integrand of constraint 2.

Facets.

Evolver takes care of these automatically.

Gaps.

We subtract some radial component from F⃗ to get a divergenceless vector F⃗_{G} with the same surface integral over the gap surfaces:

F⃗_{G} = −\frac{1}{6}xi⃗ + \frac{1}{3}yj⃗ − \frac{1}{6}zk⃗.  (11)

This has vector potential

w⃗ = \frac{1}{6}zyi⃗ − \frac{1}{6}xyk⃗.  (12)

Hence

∫ ∫_{gaps} F⃗ · dA = ∫_{gap trace} w⃗ · ds + ∫_{gap edges} w⃗ · ds.  (13)

The second part we can put directly into the constraint 2 content integrand (with opposite orientation of the edges). The gap trace part we convert to an equivalent integrand v⃗:

v⃗ = \frac{R^{2}y}{6(x^{2} + z^{2})}(zi⃗ − xk⃗).  (14)

which has radial curl. Hence the gap trace integral is the same as ∫_{gap edges} v⃗ · ds, this time with the edges in the same orientation as figure 3.

§3.5.3. Gravity.

Just for fun, we’re going to put in gravity in such a way that the user can make it act in any direction. Let the acceleration of gravity be G⃗ = g_{y}j⃗ + g_{z}k⃗. g_{y} and g_{z} will be adjustable constants. Let the fuel density be ρ. Then the gravitational potential energy is

E = − ∫ ∫ ∫_{body} G⃗ · x⃗ρdV = − ∫ ∫ ∫_{body}(g_{y}y + g_{z}z)ρdV.  (15)

First we change this to a surface integral

E = −ρ ∫ ∫_{∂body} F⃗ · d⃗A,  ∇ · F⃗ = g_{y}y + g_{z}z.  (16)

It turns out convenient to take

F⃗ = (g_{y} \frac{y^{2}}{2} + g_{z}yz)j⃗.  (17)

The surface integral has four components:

Base.

This part is zero since y = 0.

Wall.

This part is also zero since F⃗ is parallel to the wall.

Facets.

This part is taken care of by installing −ρF⃗ as the integrand for surface energy 1 for the facets.

Gaps.

We get an equivalent divergence-free integrand F⃗_{G} by adding a radial component to F⃗:

F⃗_{G} = g_{y} (\frac{y^{2}}{2} j⃗ − \frac{y}{2} (xi⃗ + zk⃗)) + g_{z} (yzj⃗ − \frac{z}{3}(xi⃗ + zk⃗)) .  (18)

This has vector potential w⃗

w⃗ = (g_{y} \frac{y^{2}}{4} + g_{z} \frac{yz}{3} ) (zi⃗ − xk⃗).  (19)

This is integrated over the gap edges with negative orientation, so ρw⃗ goes into the energy integrand of constraint 2, and over the gap trace. For the gap trace part, we construct an equivalent w⃗_{G} that has radial curl:

w⃗_{G} = (g_{y} \frac{y^{2}}{4} \frac{R^{2}}{(x^{2} + y^{2})} + g_{z} \frac{yz}{3} \frac{R^{3}}{(x^{2} + z^{2})^{(3/2)}} ) (zi⃗ − xk⃗),  (20)

so

∫_{gap trace} w⃗ · d⃗s = ∫_{gap edges} w⃗_{G} · d⃗s.  (21)

Here the integral is over the same orientation as in figure 3, so −ρw_{G} also goes into the energy integrand of constaint 2.

§3.5.4. Running.

Refine once before iterating. Iterate and refine to your heart’s content. You might try command u once in a while to regularize the triangulation. Play around with command A to change gravity and wall contact angles. Keep an eye on the scale factor. If it gets towards zero, check out the t and w histograms.

// tankex.fe
// Equilibrium shape of liquid in flat-ended cylindrical tank.
// Tank has axis along y-axis and flat bottom in x-z plane.    This
// is so gravity acting vertically draws liquid toward wall.
// Straight edges cannot conform exactly to curved walls.
// We need to give them an area so that area cannot shrink by straght edges
// pulling away from the walls.  The gaps are also accounted for
// in volume and gravitational energy.
SYMMETRIC_CONTENT     // for volume calculations
// Contact angles, initially for 45 degrees.
PARAMETER ENDT  = 0.707   /* surface tension of uncovered base   */
PARAMETER WALLT = 0.707   /* surface tension of uncovered side wall */

// Gravity components
PARAMETER GY = 0
PARAMETER GZ = -1
SPRING_CONSTANT 1  // for most accurate gap areas for constraint 2
#define  TR  1.00       /* tank radius */
#define  RHO 1.00       /* fuel density */
constraint 1  // flat base
function:    y = 0
energy:                              // for contact energy line integral
e1:      -ENDT*z
e2:      0
e3:      0
#define wstuff (WALLT*TR*y/(x^2 + z^2))   // common wall area term
#define vstuff (TR^2/3*y/(x^2 + z^2))     // common wall volume term
#define gstuff (GY*TR^2*y^2/4/(x^2 + z^2) + GZ*TR^3*y*z/3/(x^2+z^2)**1.5)
// common gap gravity term
constraint 2  CONVEX         // cylindrical wall
function:    x^2 + z^2 = TR^2
energy:                                // for contact energy and gravity
e1:    -wstuff*z  + RHO*GY*y^2*z/4 + RHO*GZ*y*z^2/3 - RHO*gstuff*z
e2:        0
e3:      wstuff*x  - RHO*GY*y^2*x/4 - RHO*GZ*y*z*x/3 + RHO*gstuff*x
content:                              // so volumes calculated correctly
c1:      vstuff*z - z*y/6 + vstuff*z/2
c2:        0
c3:    -vstuff*x + x*y/6 - vstuff*x/2
surface_energy 1   // for arbitrary direction gravity on facets
e1:    0
e2:    -RHO*GY*y^2/2 - RHO*GZ*y*z
e3:    0

// Now the specification of the initial shape
vertices
1    0.5  0.0  0.5  constraint 1
2    0.5  0.0 -0.5  constraint 1
3   -0.5  0.0 -0.5  constraint 1
4   -0.5  0.0  0.5  constraint 1
5    1.0  0.5  0.0  constraint 2
6    0.0  0.5 -1.0  constraint 2
7   -1.0  0.5  0.0  constraint 2
8    0.0  0.5  1.0  constraint 2
edges
1    2  1  constraint 1
2    1  4  constraint 1

3    4  3  constraint 1
4    3  2  constraint 1
5    5  6  constraint 2
6    6  7  constraint 2
7    7  8  constraint 2
8    8  5  constraint 2
9    1  8
10   1  5
11   2  5
12   2  6
13   3  6
14   3  7
15   4  7
16   4  8
faces
1   13   6 -14 energy 1
2    3  14 -15 energy 1
3   15   7 -16 energy 1
4    2  16  -9 energy 1
5    9   8 -10 energy 1
6    1  10 -11 energy 1
7   11   5 -12 energy 1
8    4  12 -13 energy 1
bodies
1    1 2 3 4 5 6 7 8   volume 0.6  density 0  /* no default gravity */

§3.6. Example 5. Spherical tank.

This example discusses a sperical tank in low gravity with some liquid inside. Complete understanding of this example qualifies you as an Evolver Power User.

The initial configuration is shown in figure 4. The liquid comprises one body, and the boundary of the body has three components:

1. The sphere at x^{2} + y^{2} + z^{2} = R^{2}, implemented as constraint 1.

2. The faceted surface (initially, the single face).

3. The gaps between the edges on the sphere (the “gap edges”, edges 1,2,3,4 in figure 4) and the sphere. We will assume that the gap surfaces are generated by lines radially outward from the center of the sphere through the points of the gap edges. We will often be looking for vectorfields that have radial curl, since their integral over the gap edges

will be the same as over the gap trace. As a useful fact for all the curls below, note that if

w⃗ = \frac{z^{n}}{(x^{2} + y^{2})(x^{2} + y^{2} + z^{2})^{n/2}} (yi⃗ − xj⃗)  (1)

then

∇ × w⃗ = \frac{nz^{n−1}}{(x^{2} + y^{2} + z^{2})^{1+n/2}} (xi⃗ + yj⃗ + zk⃗).  (2)

To specify contact angles, let θ_{W} be the interior contact angle with the spherical wall.

§3.6.1. Surface energy.

Wall.

The contact angle is represented by pretending the uncovered wall has surface tension T_{W} = cos(θ_{W}). We will assume equivalently that the covered wall has tension −T_{W}. The energy is

E = ∫ ∫_{wall} −T_{W} dA = ∫ ∫_{wall} −T_{W} \frac{1}{R}(xi⃗ + yj⃗ + zk⃗) · d⃗A.  (3)

We want to turn this into an integral over the boundary of the covered wall, the “gap trace” where the gap surfaces meet the wall. This calls for a vector potential, but our integrand above does not have zero divergence. But we could choose a divergenceless integrand that is equal on the wall:

E = ∫ ∫_{wall} F⃗ · d⃗A,  (4)

where

F⃗ = −T_{W} \frac{R^{2}}{(x^{2} + y^{2} + z^{2})^{3/2}} (xi⃗ + yj⃗ + zk⃗).  (5)

Now we can get a vector potential w⃗ (see useful fact above with n = 1)

w⃗ = −T_{W} \frac{R^{2}z}{(x^{2} + y^{2})(x^{2} + y^{2} + z^{2})^{1/2}} (yi⃗ − xj⃗), ∇ × w⃗ = F⃗ .  (6)

Hence

E = ∫_{gap trace} w⃗ · ds.  (7)

This is equal to the integral over the gap edges, since ∇ × w⃗ is parallel to the gap surfaces. So

E = ∫_{gap edges} w⃗ · ds,  (8)

where the gap edges are oppositely oriented from figure 4. Hence the constraint energy integrand includes −w⃗.

The facets.

The surface tension energy here is automatically taken care of by the Evolver. We will leave the surface tension as the default value of 1.

The gaps.

We will declare constraint 1 to be CONVEX so that Evolver will include a gap surface energy for all the gap edges. We put the spring constant to be 1, as that value give the best approximation to the true area for small gaps.

§3.6.2. Volume.

The body volume V = ∫ ∫ ∫ dV can be calculated as a surface integral V = ∫ ∫ F⃗ · d⃗A where ∇ · F = 1. The integral over the facets is handled automatically, but there are two choices for F⃗. We will take F⃗ = (xi⃗ + yj⃗ + zk⃗)/3 by putting SYMMETRIC_CONTENT in the datafile. Thus

V = \frac{1}{3} ∫ ∫_{body surface}(xi⃗ + yj⃗ + zk⃗) · dA.  (9)

The surface integral has three components:

The wall.

Like wall energy, this is first put in the equivalent divergenceless form

F⃗ = \frac{1}{3} \frac{R^{3}}{(x^{2} + y^{2} + z^{2})^{3/2}} (xi⃗ + yj⃗ + zk⃗).  (10)

This has vector potential (again useful fact above with n = 1)

w⃗ = \frac{1}{3} \frac{R^{3}z}{(x^{2} + y^{2})(x^{2} + y^{2} + z^{2})^{1/2}} (yi⃗ − xj⃗), ∇ × w⃗ = F⃗ .  (11)

By the same reasoning as for wall area, this reduces to the integral of −w⃗ over the gap edges. Thus we include −w⃗ in the content integrand of constraint 1. However, we also note a potential problem at z = −R where x^{2} + y^{2} = 0. We must be sure that this point always remains within the covered surface. As long as it does, it will have a constant effect. In fact, we see that at z = 0 the integrand w⃗ is zero, and the other volume contributions are zero at z = 0, so the singularity at z = −R must have a contribution to the volume of half the volume of the sphere. We can account for this by using VOLCONST 2*pi*R^3/3 with the body.

Facets.

Evolver takes care of these automatically.

Gaps.

The integrand is parallel to the gap surfaces, hence this contribution is zero.

§3.6.3. Gravity.

The sphere being very symmetric, all directions of gravity are the same. Therefore we will use Evolver’s default gravity. Let the magnitude of gravity be G (downwards). Let the liquid density be ρ. Then the gravitational potential energy is

E = G ∫ ∫ ∫_{body} zρdV.  (12)

First we change this to a surface integral

E = Gρ ∫ ∫_{∂body} F⃗ · d⃗A,  ∇ · F⃗ = z.  (13)

The default gravity mechanism uses

F⃗ = Gρ\frac{z^{2}}{2} k⃗.  (14)

The surface integral has three components:

Wall.

First, we convert to the equivalent divergence-free integrand

F⃗ = Gρ\frac{R^{4}z^{3}}{2(x^{2} + y^{2} + z^{2})^{3}} (xi⃗ + yj⃗ + zk⃗).  (15)

Note this has the same dot product with xi⃗ + yj⃗ + zk⃗ on the surface of the sphere. The vector potential is

w⃗ = Gρ\frac{R^{4}z^{4}}{8(x^{2} + y^{2})(x^{2} + y^{2} + z^{2})^{2}} (yi⃗ − xj⃗), ∇ × w⃗ = F⃗ .  (16)

So

E = ∫_{gap trace} w⃗ · d⃗s − \frac{GρR^{4}}{4} .  (17)

The constant term here is the gravitaional energy of the lower hemisphere, since w⃗ = 0 when z = 0. It will not show up in the Evolver energy; you must remember to adjust for it yourself. Due to the radial curl of w⃗, the integral is also the integral over the gap edges. The orientation is opposite that given for the gap edges, so −w⃗ goes in the energy for constraint 1.

Facets.

The integrals over facets are automatically taken care of by the Evolver.

Gaps.

We get an equivalent divergence-free integrand F⃗_{G} by adding a radial component −(z/4)(xi⃗ + yj⃗ + zk⃗) to F⃗:

F⃗_{G} = \frac{Gρ}{4}(−xzi⃗ + −yzj⃗ + z^{2}k⃗).  (18)

This has vector potential w⃗

w⃗ = −\frac{Gρz^{2}}{8} (yi⃗ − xj⃗).  (19)

This is integrated over the gap edges with opposite orientation, so −w⃗ goes into the energy integrand of constraint 1, and over the gap trace. For the gap trace part, we construct an equivalent w⃗_{G} that has radial curl (using two versions of useful fact with n = 2 and n = 4 to show radial curl):

\begin{aligned}w⃗_{G} & = \frac{Gρ}{8} \frac{z^{2}(R^{2} − z^{2})}{(x^{2} + y^{2})}(yi⃗ − xj⃗) \\ & = \frac{Gρ}{8} (\frac{z^{2}R^{4}}{(x^{2} + y^{2})(x^{2} + y^{2} + z^{2})} − \frac{z^{4}R^{4}}{(x^{2} + y^{2})(x^{2} + y^{2} + z^{2})^{2}} ) (yi⃗ − xj⃗) \\ & = \frac{Gρ}{8} \frac{z^{2}R^{4}}{(x^{2} + y^{2} + z^{2})^{2}} (yi⃗ − xj⃗) \end{aligned}  (20)

so

∫_{gap trace} w⃗ · d⃗s = ∫_{gap edges} w⃗_{G} · d⃗s.  (21)

Here the integral is over the same orientation as in figure 4, so w_{G} also goes into the energy integrand of constaint 1.

§3.6.4. Running.

Refine once before iterating. Iterate and refine to your heart’s content. You might try command u once in a while to regularize the triangulation. Play around with command A to change gravity and wall contact angles. Keep an eye on the scale factor. If it gets towards zero, check out the t and w histograms.

    // sphere.fe
    // Spherical tank.
    SYMMETRIC_CONTENT   // natural for a sphere
    /* sphere radius */
    #define R  1.00
    /* cosine of contact angle */
    PARAMETER WALLT = 0.707
    /* density of body */
    #define RHO   1.0

#define wstuff (R^2*z/(x^2 + y^2)/sqrt(x^2+y^2+z^2))
#define gstuff (G*RHO/8*R^4*z^4/(x^2 + y^2)/(x^2+y^2+z^2)^2)
#define gapstuff (G*RHO*R^4/8*z^2/(x^2+y^2+z^2)^2)
constraint 1    /* the sphere */
formula:        x^2 + y^2 + z^2 = R^2
energy:
e1:        WALLT*wstuff*y - gstuff*y + G*RHO/8*y*z^2 + gapstuff*y
e2:      -WALLT*wstuff*x + gstuff*x - G*RHO/8*x*z^2 - gapstuff*x
e3:0
content:
c1:    -wstuff*y/3
c2:      wstuff*x/3
c3:      0
vertices:
1    1  0  0  constraint 1
2    0  1  0  constraint 1
3   -1  0  0  constraint 1
4    0 -1  0  constraint 1
5    0  0  1  FIXED   // to show back hemisphere
6    0  0 -1  FIXED
7    0  1  0  FIXED
8   -1  0  0  FIXED
9    0 -1  0  FIXED
edges:
1    1  2     constraint 1
2    2  3     constraint 1
3    3  4     constraint 1
4    4  1     constraint 1
5    5  9     constraint 1  FIXED // to show back hemisphere
6    9  6     constraint 1  FIXED
7    6  7     constraint 1  FIXED
8    7  5     constraint 1  FIXED
9    5  8     constraint 1  FIXED
10   9  8     constraint 1  FIXED
11   6  8     constraint 1  FIXED
12   7  8     constraint 1  FIXED
faces:
1    1 2 3 4
2    5 10 -9  constraint 1 density 0 FIXED   // to show back hemisphere
3   -10 6 11  constraint 1 density 0 FIXED
4   -11 7 12  constraint 1 density 0 FIXED
5    8 9 -12  constraint 1 density 0 FIXED
bodies:    // start with sphere half full
1    1    volume 2*pi*R^3/3  volconst 2*pi*R^3/3

§3.7. Example 6. Torus partitioned into two cells.

This example has a flat 3-torus divided into two bodies. The unit cell is a unit cube, and the surface has the topology of Kelvin’s partitioning of space into tetrakaidecahedra, which is the least area partitioning known of space into equal volumes. The datafile handles the wrapping of edges around the torus by specifying for each direction whether an edge wraps positively (+), negatively (-), or not at all (*).

The display of a surface in a torus can be done several ways. When displaying this one, try both the connected bodies and the clipped cell options. For a screen display, you can switch between the two with the t command. Interactive geomview or tt MinneView users can use the clipped, connected, and raw cells commands to set the option. users will have to end and restart MinneView (command P).

// twointor.fe
// Two Kelvin tetrakaidecahedra in a torus.
TORUS_FILLED   // signals that domain is a torus and bodies fill it.
periods
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
vertices // values from another program
1  0.499733 0.015302 0.792314
2  0.270081 0.015548 0.500199
3  0.026251 0.264043 0.500458
4  0.755123 0.015258 0.499302
5  0.026509 0.499036 0.794636
6  0.500631 0.015486 0.293622
7  0.025918 0.750639 0.499952
8  0.499627 0.251759 0.087858
9  0.256701 0.499113 0.087842
10 0.026281 0.500286 0.292918
11 0.500693 0.765009 0.086526
12 0.770240 0.499837 0.087382
edges // with wraps in axis directions
1    1 2  * * *
2    2 3  * * *
3    1 4  * * *
4    3 5  * * *
5    2 6  * * *
6    2 7  * - *
7    1 8  * * +
8    4 6  * * *
9    5 9  * * +
10   3 10 * * *
11   3 4  - * *
12   6 8  * * *
13   6 11 * - *
14   7 4  - + *

15   8 12 * * *
16   9 8  * * *
17   9 11 * * *
18   10 7 * * *
19   11 1 * + -
20   12 5 + * -
21   5 7  * * *
22  11 12 * * *
23  10 12 - * *
24   9 10 * * *
faces
1    1 2 4 9 16 -7
2    -2 5 12 -16 24 -10
3    -4 10 18 -21
4    7 15 20 -4 11 -3
5    -1 3 8 -5
6    6 14 -11 -2
7    5 13 -17 24 18 -6
8    -12 13 19 7
9    -16 17 22 -15
10   -10 11 8 12 15 -23
11   -21 9 17 19 1 6
12   -14 -18 23 -22 -13 -8
13   -24 -9 -20 -23
14   -19 22 20 21 14 -3
bodies
1    -1 -2 -3 -4 -5 9 7 11 -9 10 12 5 14 3 volume 0.500
2     2 -6 -7 8 -10 -12 -11 -13 1 13 -14 6 4 -8 volume 0.500

§3.8. Example 7. Crystalline integrand.

This example shows how to use Wulff vectors for a crystalline integrand. In general, the Wulff vector at a point on a surface is a vector that is dotted with the surface normal to give the energy density of the surface. For ordinary area, it is the unit normal itself. But if a finite set of Wulff vectors is given, and the one used at a point is the one with the maximum dot product, then the minimal surface has only flat faces, i.e. it is crystalline. The crystal shape in this example is an octahedron. The Wulff vectors are the vertices of the octahedron, and are in octa.wlf. The datafile crystal.fe describes a cube with extra edges around the equator. These edges provide places for octahedron edges to form. Things work a lot better if the facet edges are where the crystal edges want to be.

To see the cube turn into an octahedron, refine once, use the m command to set the scale to 0.05, and iterate about 50 times. Fixed scale works much better than optimizing scale for crystalline integrands, since crystalline energy is not a smooth function of the scale factor. Sometimes you have to put in jiggling also to get over local obstacles to minimizing energy.

Here is the Wulff vector file octa.wlf:

0      0      1
0      0     -1
0.707  0.707  0

-0.707  0.707  0
0.707 -0.707  0
-0.707 -0.707  0

Here is the datafile:

// crystal.fe
// Evolver data for cube of prescribed volume with
//   crystalline surface energy.    Evolves into an octahedron.
Wulff "octa.wlf"    // Wulff vectors for octahedral crystal
vertices
1  0.0 0.0 0.0
2  1.0 0.0 0.0
3  1.0 1.0 0.0
4  0.0 1.0 0.0
5  0.0 0.0 1.0
6  1.0 0.0 1.0
7  1.0 1.0 1.0
8  0.0 1.0 1.0
9  0.0 0.0 0.5
10 1.0 0.0 0.5
11 1.0 1.0 0.5
12 0.0 1.0 0.5
edges  /* given by endpoints and attribute */
1   1 2
2   2 3
3   3 4
4   4 1
5   5 6
6   6 7
7   7 8
8   8 5
9   1 9
10  2 10
11  3 11
12  4 12
13  9 5
14  10 6
15  11 7
16  12 8
17   9 10
18  10 11
19  11 12
20  12  9
faces  /* given by oriented edge loop */
1   1 10 -17  -9
2   2 11 -18 -10

3   3 12 -19 -11
4   4  9 -20 -12
5   5  6  7   8
6  -4 -3 -2  -1
7  17 14 -5 -13
8  18 15 -6 -14
9  19 16 -7 -15
10 20 13 -8 -16
bodies  /* one body, defined by its oriented faces */
1   1 2 3 4 5 6 7 8 9 10  volume 1

Chapter 4. The Model.

This chapter describes the mathematical model of a surface that the Surface Evolver incorporates. See the Datafile section below for how to set up the initial surface and the Commands section for how to manipulate the surface. Keywords in the datafile are given in upper case, even though case is not significant.

§4.1. Geometric elements.

The surface is defined in terms of its geometric elements of each dimension. Each element has its own set of attributes. Some may be set by the user; others are set internally but may be queried by the user.

§4.1.1. Vertices.

A vertex is a point in space. The coordinates of the vertices are the parameters that determine the location of the surface. It is the coordinates that are changed when the surface evolves.

Attributes:

coordinates - either Euclidean coordinates or boundary parameters. fixed - vertex cannot be moved or eliminated. constraints - vertex may be subject to a set of constraints. boundary - vertex position is defined parametrically on a boundary. A vertex can be on at

most one boundary.

id - an identification number. original - number of the vertex in the original datafile

§4.1.2. Edges.

An edge is an oriented line segment between a tail vertex and a head vertex. Attributes:

endpoints - numbers of the tail and head vertices of the edge. length - actual length of the edge. fixed - endpoints fixed, and all vertices generated from subdividing edge will be fixed. constraints - set of constraints edge satisfies. Constraints may be used to specify volume and

energy integrands.

boundary - edge may be in one boundary. torus wrapping - whether edge crosses torus unit cell faces, and in which direction. density - energy per unit length in string model. color - for display. id - an identification number. original - number of the original datafile edge this edge is descended from.

§4.1.3. Facets.

A facet is an oriented triangle defined by its edges. Attributes:

edges - loop of three edges connected head to tail. area - actual area of the facet. fixed - fixedness inherited by all elements descended from this facet. constraints - set of constraints facet satisfies. boundary - facet may be in one boundary. density - energy per unit area; surface tension. energy - whether facet has a surface energy integrand. quantities - facet contribution to user-defined quantity integrals. displayable - whether facet should be displayed in graphics. tag - a user-defined number that can identify groups of facets. id - an identification number. original - number of the original datafile edge this edge is descended from. color - for display. frontcolor, backcolor - for coloring the two sides differently. phase - for determining boundary tension in string model.

§4.1.4. Bodies.

A body is a three-dimensional region of space. Attributes:

facets - a set of oriented facets defines the boundary of body. These facets are used for calculating body volume and gravitational energy. Only those facets needed for correct calculation need be given.

id - an identification number. density - for gravitational potential energy. volume - actual volume. target volume - constrained volume value set by user. volume constant - an adjustment used in the torus model to do correct volume calculations. pressure - for prescribed pressure, or for pressure resulting from volume constraint. phase - for determining boundary tension in soapfilm model.

§4.2. Quadratic represention.

By default, edges and facets are linear. It is possible to represent edges as quadratic curves and facets as quadratic patches by putting the QUADRATIC keyword in the datafile or by using the M command. If this is done, then each edge is endowed with a midpoint vertex. Most features are implemented for the quadratic representation. Those not are: torus domain, facet energy integrands, simplex model, squared curvature, approximate curvature.

§4.3. Dimension of surface.

In the standard Evolver surface, the surface tension resides in the two-dimensional facets (this is referred to as the soapfilm model.). However, it is possible to have the tension reside in the edges. This is referred to as the string model and is declared by the keyword STRING in the datafile. See slidestr.fe for an example. It is possible in the soapfilm model to also have edges with tension in order to represent a situation where singular curves have energy. This chapter is written for the soapfilm model. For differences in the string model, see the String Model section below. For a surface of more than two dimensions, see Simplex model.

§4.4. Simplex model.

This model represents a surface solely with vertices and facets. It was put in to enable 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. There are no edges or facet-edges. Bodies still exist for hypersurfaces.

The datafile must have the keyword SIMPLEX_REPRESENTATION in the first section, and the phrase SURFACE_DIMENSION k if k ≠ 2. If the domain is not 3-dimensional, then SPACE_DIMENSION n must also be included. There should be no EDGES section in the datafile. 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. Vertices may be FIXED. Constraints are allowed, but no integrands of any kind. No TORUS domain. No QUADRATIC representation. No surface changing except iteration and refining. Refining subdivides each simplex 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}.

§4.5. Dimension of ambient space.

By default, surfaces live in 3 dimensional space. However, the phrase SPACE_DIMENSION n in the datafile 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 extern.h. Change MAXCOORD and recompile if you want more than three dimensions.

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).

§4.6. Riemannian metric.

The ambient space can be endowed with a Riemannian metric by putting the keyword METRIC or CONFORMAL_METRIC in the datafile followed by the elements of the metric tensor. A conformal metric tensor is a multiple of the identity matrix, so only one element is needed. 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 “quantity” constraint. See quadm.fe for an example of a metric.

The conformal metric can be used to simulate a scalar surface integrand, since the only use of the metric is in calculating lengths and areas.

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.

§4.7. Torus domain.

By default, the domain of a surface is Euclidean space R^{3}. However, there are many interesting problems dealing with periodic surfaces. 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.

Certain features are not available with the torus domain, namely quadratic representation, constraints, boundaries, quantities, and gravity. (Actually, you could put them in, but they will not take into account the torus wrapping.)

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.

Vertex coordinates are given as Euclidean coordinates within the unit cell, not as linear combinations of the basis 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, 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.

There are several options for displaying a torus surface:

Raw cell - Graph the facets as they are. Good for checking your datafile if it’s all triangles. Faces that get subdivided may be unrecognizable. Connected bodies - Each body’s facets are unwrapped in the torus, so the body appears in one connected piece. Nicest option. Clipped cell - Shows the unit cell specified in the datafile. Facets are clipped on the parallelpiped faces.

§4.8. Quotient spaces.

As a generalization of the torus domain, you may declare the domain to be the quotient space of R^{n} with respect to some symmetry group. This cannot be completely specified in the datafile, but requires you to write some C routines to define group operations. Group elements are represented by integers attached to edges (like the wrap specifications in the torus model). You define the integer representation of the group elements. See the files quotient.c and symtest.fe 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 your 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.

§4.9. Symmetric surfaces.

Symmetric surfaces can often be done more easily by evolving just one fundamental domain. This is easiest if the fundamental region is bounded by planes of mirror symmetry. These planes become constraint surfaces for the fundamental domain (see Constraints below). Multiple images of the fundamental domain may be displayed by using transformation matrices (see view transforms in the Datafile chapter).

§4.10. Constraints.

A constraint is a barrier that restricts the motion of the surface. These may be represented as level sets of functions or as parametric manifolds. The term “constraint” as used in the Evolver refers to the level set formulation, and “boundary” refers to the parametric formulation. The term “volume constraint” is something entirely else, refering to bodies.

A constraint may have several roles:

1. Vertices may be required to lie on a constraint (equality constraint) or on one side (inequality constraint). A constraint may be declared GLOBAL, in which case it applies to all vertices. See mound.fe for an example.

2. A constraint may have a vectorfield associated with it that is integrated over edges lying in the constraint to give an energy. This is useful for specifying wall contact angles and for calculating gravitational energy. Integrals are not evaluated over edges that are FIXED. See mound.fe for an example.

3. A constraint may be declared CONVEX, in which case edges in the constraint have an energy associated with them that is proportional to the area between the straight edge and the curved wall. This energy (referred to as “gap energy”) is meant to compensate for the tendency for flat facets meeting a curved wall to minimize their area by lengthening some edges on the wall and shortening others, with the net effect of increasing the net gap between the edges and the wall. See tankex.fe for an example.

4. A constraint may have a vectorfield associated with it that is integrated over edges lying in the constraint to give a volume contribution to a body whose boundary facets contains the edges. This is useful for getting correct volumes for bodies without completely surrounding them with otherwise useless facets. It is important to understand how the content is added to the body in order to get the signs right. The integral is evaluated along the positive direction of the edge. If the edge is positively oriented on a facet, and the facet is positively oriented on a body, then the integral is added to the body. This may wind up giving the opposite sign to the integrand from what you think may be natural. Integrals are not evaluated over edges that are FIXED. See tankex.fe for an example.

5. A constraint may have a vectorfield associated with it that is integrated as a surface integral over facets lying in the constraint to give a contribution to the total energy. This is useful for having gravity in an arbitrary direction. Integrals are not evaluated over facets that are FIXED. See tankex.fe for an example.

Edge integrands are evaluated by gaussian quadrature. The number of points used is controlled by the INTEGRAL_ORDER option in the datafile.

§4.11. Boundaries.

Boundaries are one or two dimensional parameterized manifolds. Vertices, edges, and facets may be deemed to lie in a boundary. For a vertex, this means that the fundamental parameters of the vertex are the parameters of the boundary, and its coordinates are calculated from these. Vertices on boundaries may move during iteration, unless declared fixed. Edges on boundaries do not have energy and content integrals (yet), since one-parameter boundaries are most frequent. See cat.fe for an example.

A delicate point is how to handle wrap-arounds on a boundary such as a circle or cylinder. Subdividing a boundary edge requires a midpoint, but taking the average parameters of the endpoints can give nonsense. Therefore the average coordinates are calculated, and that point projected on the boundary as continued from one endpoint.

A general guideline is to use constraints for two-dimensional walls and boundaries for one-dimensional wires. If you are using a boundary wire, you can probably declare the vertices and edges on the boundary to be FIXED. Then the boundary becomes just a guide for refining the boundary edges.

NOTE: A vertex on a boundary cannot also have constraints.

§4.12. Energy.

The Evolver works by minimizing the total energy of the surface. This energy can have several components:

1. Surface tension. Soap films and interfaces between different fluids have an energy content proportional to their area. Hence they shrink to minimize energy. The energy per unit area can also be regarded as a surface tension, or force per unit length. Each facet has a surface tension, which is 1 unless the datafile specifies otherwise (see TENSION option for faces). Different facets may have different surface tensions. Facet tensions may be changed interactively with the ‘set facet tension ...’ command. The contribution to the total energy is the sum of all the facet areas times their respective surface tensions. The surface tension of a facet may also be specified as depending on the phases of the bodies it separates. See

PHASE.

2. Crystalline integrands. The Evolver can model energies of crystalline surfaces. These energies are proportional to the area of a facet, but they also depend on the direction of the normal. The energy is given by the largest dot product of the surface normal with a set of vectors known as the Wulff vectors. Surface area can be regarded as a crystalline integrand whose Wulff vectors are the unit sphere. See the datafile section on Wulff vectors for more. A surface has either crystalline energy or surface tension, not both.

3. Gravitational potential energy. If a body has a density (see DENSITY option for bodies), then that body contributes its gravitational energy to the total. The acceleration of gravity G is under user control. Letting ρ be the body density, the energy is defined as

E = ∫ ∫ ∫_{body} GρzdV  (1)

but is calculated by the Divergence Theorem as

E = ∫ ∫_{body surface} Gρ\frac{z^{2}}{2} k⃗ · d⃗S.  (2)

This integral is done over each facet that bounds a body. If a facet bounds two bodies of different density, then the appropriate difference in density is used. Facets lying in the z = 0 plane make no contribution, and may be omitted if they are otherwise unneeded. Facets lying in constraints may be omitted if their contributions to the gravitational energy are contained in constraint energy integrals.

4. Constraint edge integrals. An edge on a constraint may have an energy given by integrating a vectorfield F⃗ over the oriented edge:

E = ∫_{edge} F⃗ · d⃗l.  (3)

This is useful for prescribed contact angles on walls (in place of wall facets with equivalent tension) and for gravitational potential energy that would otherwise require facets in the constraint.

5. Surface energy integrals. A facet may contribute an energy resulting from the integration of a vectorfield over the facet as a surface integral. This is useful in modelling gravity from various directions, when the implicit gravity energy doesn’t apply. Because this may apply to every facet, it can slow things down. See the Surface Energies section of the datafile.

6. Convex constraints. Consider a soap film spanning a circular cylinder. The Evolver must approximate this surface with a collection of facets. The straight edges of these facets cannot conform to the curved wall, and hence the computed area of the surface leaves out the gaps between the outer edges and the wall. The Evolver will naturally try to minimize area by moving the outer vertices around so the gaps increase, ultimately resulting in a surface collapsed to a line. This is not good. Therefore there is provision for a “gap energy” to discourage this. A constraint may be declared CONVEX in the datafile. For an edge on such a constraint, an energy is calculated as

E = k ∥S⃗ × Q⃗∥ /6  (4)

where S⃗ is the edge vector and Q⃗ is the projection of the edge on the tangent plane of the constraint at the tail vertex of the edge. The constant k is a global constant called the “gap constant”. A gap constant of 1 gives the best approximation to the actual area of the gap. A larger value minimizes gaps and gets vertices nicely spread out along a constraint.

The gap energy falls off quadratically as the surface is refined. That is, refining once reduces the gap energy by a factor of four. You can see if this energy has a significant effect on the surface by changing the value of k with the k command.

7. Compressible bodies. If the ideal gas mode is in effect (see PRESSURE keyword for datafile), then each body contributes an energy

E = P_{amb}V_{0} ln(V/V_{0})  (5)

where P_{amb} is the ambient pressure, V_{0} is the volume prescribed in the datafile, and V is the actual volume of the body. To account for work done against the ambient pressure, each body also makes a negative contribution of

E = −P_{amb}V.  (6)

8. Prescribed pressure. Each body with a prescribed pressure P contributes energy

E = P V.  (7)

where V is the actual volume of the body. This can be used to generate surfaces of prescribed mean curvature, since mean curvature is proportional to pressure.

9. Squared mean curvature. An average mean curvature around each vertex can be calculated, and the integral of the square of this average counted as energy. This component is included by the SQUARE_CURVATURE phrase in the datafile. A curvature offset h_{0} (so the energy is (h−h_{0})^{2}) can be specified with H_ZERO value in the datafile. The weighting factor can be changed with the A command by changing the value of the adjustable constant square curvature modulus. For details of exactly how the integral of squared curvature is calculated, see the Technical Reference section. (This feature is not implemented for the quadratic model.) Several variants of squared mean curvature are available, toggled by effective area and normal curvature commands. See the technical reference chapter for details of the formulas. [HKS] reports some experimental results.

10. Mean curvature integral. The energy can include the integral of the signed mean curvature by including the MEAN_CURVATURE_INTEGRAL keyword with weighting factor in the datafile.

11. Squared gaussian curvature. The Gauss map of a surface maps the surface to a unit sphere by mapping a point to its normal. The Gaussian curvature at a point is the Jacobian of this map. This definition can be extended to piecewise linear surfaces. The Evolver can include the integral of the squared Gaussian curvature in the energy, but it is not yet part of the force calculation.

§4.13. Named quantities and methods.

These are an effort to provide a more systematic way of adding new energies and constraints.

A “method” is a way of calculating a scalar value from some particular type of element (vertex, edge, facet, body). Each method is implemented internally as a set of functions for calculating the value and its gradient as a function of vertex positions. Methods are referred to by their names. See the Technical Reference chapter for detailed descriptions. Currently defined methods are:

 1. edge tension or edge length: Length of edge.
 2. facet tension or facet length: Area of facet.

3. edge scalar integral: Integral of a scalar function over arclength.

4. facet scalar integral: Integral of a scalar function over facet area.

5. mean curvature integral: Integral of signed scalar mean curvature.

6. gauss curvature integral: Integral of gaussian curvature.

7. hooke energy: Energy of edges as springs with reference to a global equilibrium length.

8. local hooke energy: Energy of edges as springs with equilibrium length being average of lengths of neighbor edges. (by John Sullivan; no gradient yet, so restricted to info_only).

9. knot energy: An “electrostatic” energy in which vertices are endowed with equal charges. Power law of potential is adjustable.

10. uniform knot energy or edge knot energy: A knot energy where vertex charge is proportional to neighboring edge length.

11. uniform knot energy normalizer: Supposed to approximate the part of uniform knot energy that

is singular in the continuous limit.

12. edge edge knot energy: Between pairs of edges, energy is inverse square power of distance between midpoints of edges. (by John Sullivan)

13. facet knot energy: Charge on vertex is proportional to area of neighboring facets.

14. buck knot energy: Energy between pair of edges given by formula suggested by Greg Buck.

15. proj knot energy: Only the normal component of force between vertices is used to desingularize energy.

16. circle knot energy: Another weird way to calculate a nonsingular energy between pairs of edges.

17. sphere knot energy: Another weird way to calculate a nonsingular energy between pairs of facets.

18. curvature function: Calculates forces as function of mean and gaussian curvatures at vertices. Function may be changed by user by altering teix.c. No energy, just forces.

19. average crossings: To calculate the average crossing number in all projections of a knot. (by John Sullivan)

20. gravity method: For future use in calculating gravitational energy.

Adding a new method involves writing C routines to calculate the value and the gradient as a function of vertex coordinates, and adding a structure to the method name array in quantity.c. All the other syntax for invoking it from the datafile is already in place.

A “method instance” is a particular use of a method, with whatever particular parameters may be needed. A method instance may be created explicitly and given a name, or they may be created implicitly by defining quantities to use methods. See the Datafile chapter for details. Every method instance has a “modulus”, which is multiplied times the basic method value to give the instance value. A modulus of 0 causes the entire instance calculation to be omitted.

A “quantity” is the sum total of various method instances although usually just one instance is involved. Any quantity may be declared to be one of three types: 1) “energy” quantities are added to the overall energy of the surface; 2) “fixed” quantities that are constrained to a fixed target value (by a single Newton step at each iteration); and 3) “info only” quantities whose values are merely reported to the user. Each quantity has a “modulus”, which is just a scalar multiplier of the sum of all instance values. A modulus of 0 will turn off calculation of all the instances.

The sample datafile knotty.fe contains some examples. Quantity values can be seen with the Q or A command, and the A command can be used to change the target value of a fixed quantity, the modulus of a quantity, or some parameters associated to various methods.

It is planned that eventually all energies and global constraints will be converted to this system. However, existing syntax will remain valid wherever possible.

§4.14. Numbered surface integrals.

An integral quantity in the Evolver is a number that can be expressed as the sum of integrals of vectorfields over surfaces and edges. Quantities can be evaluated for informational purposes, or they can be used as constraints (in the mathematical sense). Body volumes and gravitational energy are examples of built-in quantities. Surface tension energy is not a quantity in this sense, since it cannot be written as a vector integral.

The Evolver has a mechanism to add user-defined quantities. These might include center of mass, moment of inertia, magnetic flux, etc. The key to many of these is the Divergence Theorem, which permits volume integrals to be written as surface integrals. You can also use a quantity in place of a body volume that is awkward to do with the built-in volume mechanisms.

Each quantity when first defined in the datafile has given the surface integrand, which may be integrated over all or some facets. For edge integrals, each constraint specification may include integrands for quantities, which are applied to edges subject to that constraint. A quantity acts as a mathematical constraint when it has the FIXED attribute and a specified value. These may be displayed and changed interactively with the Q command. Values may also be displayed with the v command (since volumes are themselves quantities, anyway). See qmound.fe for an example of body volume implemented as a quantity.

At present, these are implemented numbered “surface energy” integrals. This mechanism will be replaced by a named quantity mechanism soon.

§4.15. Pressure.

Pressure is force per unit area perpendicular to a surface. In the Evolver, pressure can arise in three ways:

1. If a body has a volume constraint, then the boundary surface is unlikely to be a minimal surface. Hence pressure is needed to counteract the desire of the surface to shrink. When there are volume constraints, the Evolver automatically calculates the pressure needed. The value of the pressure can be seen with the v command.

2. A body may have a prescribed pressure. Then the appropriate force is added to the forces on the vertices when calculating the motion of the surface. This is the way to get a surface of prescribed mean curvature, since pressure = surface tension × mean curvature. NOTE: Prescribed volume and prescribed pressure on the same body are illegal. CAUTION: Prescribed pressure can make a surface expand to infinity if the pressure is too large.

3. The Evolver can treat bodies as being made up of an isothermal ideal gas, that is, bodies can be compressible. This happens if the PRESSURE keyword comes in the first part of the datafile. The pressure given there is the ambient pressure outside all bodies. Each body must have a volume specified, which is the volume of the body at the ambient pressure.

§4.16. Volume or content.

The terms “volume” and “content” are used for the highest dimensional measure of a region in n-space: area in R^{2}, volume in R^{3}, etc. Bodies may have a volume specified in the datafile, which then becomes a volume constraint. The volume of a body can be written as

V = ∫ ∫ ∫_{body} 1dV,  (1)

which by the Divergence Theorem can be written a surface integral:

V = ∫ ∫_{body surface} zk⃗ · d⃗S.  (2)

This integral is evaluated over all the boundary facets of a body.

The part of the boundary of a body lying on a constraint need not be given as facets. In that case, Stokes’ Theorem can be used to convert the part of the surface integral on the constraint to a line integral over the edges where the body surface meets the constraint. The line integral integrands are given as constraint content integrands in the datafile.

An alternate surface integral is

V = \frac{1}{3} ∫ ∫_{body surface}(xi⃗ + yj⃗ + zk⃗) · d⃗S.  (3)

This is used if SYMMETRIC_CONTENT is specified in the datafile. It is useful if the constraint content integral (which is evaluated by an approximation) leads to asymmetric results on what should be a symmetric surface.

As with surface area, the gap between flat facets and curved constraints causes problems. You can use constraint content integrals to overcome this. See the tank and sphere examples in the tutorial section.

§4.17. 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

dm = (scale factor)(diffusion constant)(facet area)(pressure difference).  (1)

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.

§4.18. Motion.

The heart of the Evolver is the iteration step that reduces energy while obeying any constraints. The surface is changed by moving the vertices. No changes in topology or triangulation are made (unless certain options are toggled on). The idea is to calculate the velocity at each vertex and move the vertex in that direction. If the conjugate gradient option is in effect, the direction of motion is adjusted accordingly.

The force is the negative gradient of the energy. At each vertex, the gradient of the total energy (see the Energy section above) as a function of the position of that vertex is calculated. By default, the velocity is taken to be equal to the force. If the area normalization option is in effect, then the force is divided by 1/3 the area of the neighboring vertices to get an approximation to the mean curvature vector for velocity. If relevant, the gradients of quantity constraints (such as body volumes), and level set constraints are calculated at each vertex. For fixed vertices, all of these are taken as zero. Forces and quantity constraint gradients are projected onto the tangent space of the relevant level set constraints at each vertex. To enforce quantity constraints, a restoring motion along the volume gradients is calculated (there is found a constant for each quantity such that when each vertex is moved by the sum of (constant)(quantity gradient) for all quantities including that vertex, then the quantity come out right, at least in a linear approximation). Then multiples of the quantity gradients are added to the forces to keep quantities constant (for volumes, the multiplying factors are the pressures of the bodies, which are actually found as the Lagrange multipliers for the volume constraints).

The actual motion of a vertex is the quantity restoring motion plus a universal scale factor times the velocity. The scale factor can be fixed or it can be optimized (see command m). If optimized, then the Evolver calculates the energy for several values of the scale (doubling or halving each time) until the minimum energy is passed. The optimum scale is then calculated by quadratic interpolation. A Runge-Kutta step is also available.

The vertices are then moved. If jiggling is in effect, each non-fixed vertex is moved in a random direction from a gaussian distribution with standard deviation = (temperature)(characteristic length), where the characteristic length starts out as 1 and is halved at each refinement.

Finally, 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.

There is another motion available in restricted circumstances that tries to get to the minimum in one step by calculating the second derivative matrix (Hessian) for energy, and then solving for the minimum energy. At present, it can only handle surfaces whose only energy is area and whose only constraints are fixed vertices and edges. It only works sometimes, but when it does, it reaches the minimum very quickly, perhaps two or three iterations. See the hessian command.

§4.19. 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, 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.

§4.19.1. Vertex 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.

§4.19.2. Area normalization.

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 makes the resistance 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.

§4.19.3. Area normalization with effective area.

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. 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.

§4.19.4. Approximate polyhedral curvature.

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 chapter 7 for the details. S is a sparse matrix with entries corresponding to each pair of vertices joined by an edge, and M is its dense inverse.

§4.19.5. Approximate polyhedral curvature with effective area.

The previous section did 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.

§4.19.6. 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.

§4.20. 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.

§4.20.1. 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 Δt. Let V be the vertex velocity. Then the critical timestep for amplification of the perturbation is given by V Δt = −2Y, or Δt = −2Y/V.

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

Area normalization. Here the vertex star has length 2L, so V = F/L and Δ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 Δ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.

§4.20.2. 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)Δt = 2Y,

so Δ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 Δ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.

§4.21. Topology changes.

The term “topology of the surface” refers to the topology of the union of the facets. The term “triangulation” refers to the specfic way facets subdivide the surface. Some operations, such as iteration, change neither. Some, such as refinement, change the triangulation but not the topology. Some, such as short edge elimination, can change both. Some operations, such as vertex popping, are meant to change the topology.

§4.22. Refinement.

“Refinement” of a triangulation refers to creating a new triangulation by subdividing each triangle of the original triangulation. The scheme used in the Evolver is to create new vertices at the midpoints of all edges and use these to subdivide each facet into four new facets each similar to the original.

Certain attributes of new elements are inherited from the old elements in which they were created. Fixedness, constraints, and boundaries are always inherited. Torus wrapping on edges is inherited by some but not all new edges. Surface tension and displayability are inherited by the new facets.

Refinement can change surface area and energy if there are curved constraints or boundaries. Likewise for body volumes.

§4.23. Adjustable constants and parameters.

For flexibility in the interactive manipulation of the surface, there are several ways in which various constants may be changed:

1. The gravitational constant may be changed with the G command. A value of zero turns off evaluation of gravitational energy.

2. Prescribed body volumes or pressures may be changed with the b command.

3. You can declare identifiers to be adjustable constants in the datafile and use them in constraints and boundaries in their formulas and integrands. These can be changed interactively with the A command. Changing these results in automatic recalculation of the surface, unless AUTORECALC has been toggled off.

4. The diffusion constant can be changed with the f command.

5. The “gap constant” for gap energies can be changed with the k command.

6. The ambient pressure can be changed with the p command.

7. The jiggling temperature can be changed with the j or J commands.

§4.24. The String Model.

This section lists the differences you need to be aware of when using the string model. In general, a face takes on the role of a body, and an edge takes on the role of a facet. Coordinates are

still in 3D. If you want cells of prescribed volume, always put z = 0 for the vertices and for each face make a body with just that one face as boundary.

The face section of the datafile is optional. Faces are not subdivided on input. Face area is calculated solely using x,y coordinates, even though vertices are still 3D. Content refers to area. Constraint energy and content integrands are scalars (have only one component, E1 or C1 as the case

may be, and are evaluated at vertices on constraints. Tails of edges count negative, and heads positive.

Chapter 5. The Datafile.

The initial configuration of the surface is read from an ASCII datafile.

§5.1. Lexical format.

For those who know about such things, the datafile is read with a lexical analyzer generated by the lex program. The specification is in datafile.lex.

§5.1.1. Lines and line splicing.

The file is made up of lines. Line breaks are significant. The next physical line can be spliced onto the current line by having ’\’ be the last character of the current line. Line splicing is not effective in // comments. Blank lines and comment lines may be placed freely anywhere in the datafile.

§5.1.2. Comments.

Comments may be enclosed in /* */ pairs (as in C) and may span lines. // indicates the rest of the line is a comment, as in C++.

§5.1.3. Including files.

The standard C language method of including other files is available. The file name must be in quotes. If the file is not in the current directory, EVOLVERPATH will be searched. Includes may be nested to something like 10 deep. Example:

#include “common.stuff”

§5.1.4. Macros.

Simple macros (no parameters) may be defined as in C:

#DEFINE identifier string

identifier must be an identifier without other special meaning to the parser (see the keyword list below). string is the rest of the logical line, not including comments. It will be substituted for identifier whenever identifier occurs as a token subsequently. Substitutions are re-scanned. No checks for recursiveness are made. There is a maximum length (currently 500 characters) on a macro definition. Note: macro identifiers are separate tokens, so if “-M” translates into “-2”, this will be read as two tokens, not a signed number.

§5.1.5. Case.

Case is not significant. All letters are converted to lower case on input.

§5.1.6. Whitespace.

Whitespace consists of spaces, tabs, commas, colons, and semicolons. So it’s fine if you want to use commas to separate coordinate values.

§5.1.7. Identifiers.

Identifiers follow standard C rules (composed of alphanumeric characters and ’ ’ with the leading character not a digit) and must not be keywords. Identifiers are used for macros and adjustable constants.

§5.1.8. Real numbers.

Values may be in any form acceptable to atof(). This includes integers, fixed point, and scientific notation such as

    2    -3    .5    23.        5e-10    +0.7D2

§5.1.9. Keywords.

These words have special meaning to Evolver and should not be used for user-defined names. Letternumber combinations such as “Xn” include X1, X2, ...

ABS ACOS AMBIENT PRESSURE APPROX[IMATE] CURV[ATURE] AREA AREA FIXED AREA NORMALIZATION AREAWEED ASIN ATAN AUTOCHOP AUTOPOP AUTORECALC AVG BARE BACKCOLOR BACKGROUND BLUE BLACK BODY,BODIES BODY COUNT BOUNDARY,BOUNDARIES BOUNDARY CURVATURE BROWN BYE CEIL CHECK CHECK INCREASE CLEAR CLIPPED CLIPPED CELLS CLOSE SHOW COLOR COLORMAP Cn CONF EDGE CONFORMAL METRIC CONJ GRAD CONNECTED CONNECTED CELLS CONSTRAINT CONSTRAINTS CONSTRAINT TOLERANCE CONTEXT CONVEX COS COSH CYAN DARKGRAY DEBUG DELETE DENSITY DIFFUSION DIHEDRAL DIRICHLET DISSOLVE DIV NORMAL CURVATURE DO DUMP E En EDGE,EDGES EDGE COUNT EDGE DIVIDE EDGEWEED EFFECTIVE AREA ELSE ENERGY,ENERGIES EPRINT ESTIMATE EXIT EXP EXTRAPOLATE FACE,FACES FACET,FACETS FACET COLORS FACET COUNT FACET EDGES FIX FIXED FIXED AREA FLOOR FOREACH FORMULA FRONTCOLOR FUNCTION G GAP CONSTANT GLOBAL GO GRAVITY CONSTANT GREEN HELP HESSIAN HISTOGRAM HIT CONSTRAINT HOMOTHETY ID IF INTEGRAL ORDER JIGGLE LENGTH LIGHTBLUE LIGHTCYAN LIGHTGRAY LIGHTGREEN LIGHTMAGENTA LIGHTRED LINEAR LIST LOAD LOG LOGHISTOGRAM LONGJ MAGENTA MAX MEAN CURVATURE INTEGRAL MEMDEBUG MERIT FACTOR METHOD METHOD INSTANCE METRIC MIN MOBILITY MOBILITY TENSOR MODULUS NODISPLAY NONNEGATIVE NONPOSITIVE NORMAL CURVATURE NORMAL MOTION NOT NOTCH OFF OID ON ON BOUNDARY ON CONSTRAINT OPACITY OPTIMIZE OR ORIGINAL Pn PARAMETER PARAMETERS PERIODS PI PHASE PHASEFILE PINNING POST PROJECT PRESSURE PRINT PRINTF PROCEDURE,PROCEDURES QUADRATIC Qn QUANTITY QUIET QUIT RAW CELLS RAWESTV RAWV READ REBODY RECALC RED REFINE RUNGE KUTTA SCALAR˙ INTEGRAND SCALE SCALE LIMIT SET SHADING SHELL SHOW SHOW EXPR SHOW INNER SHOW OFF SHOW OUTER SHOW TRANS SHOW VOL SHOWQ SIMPLEX REPRESENTATION SIN SINH SOAPFILM SPACE DIMENSION SPRING CONSTANT SPRINTF SQR SQCURV SQGAUSS SQRT SQUARE CURVATURE STRING SUM SURFACE DIMENSION SURFACE ENERGY SYMMETRIC CONTENT SYMMETRY GROUP SYSTEM TAG TAN TEMPERATURE TENSION THEN THICKEN TOPINFO TORUS TORUS FILLED TOTAL AREA TOTAL ENERGY TOTAL LENGTH TOTAL TIME TRANSFORM DEPTH TRANSFORM EXPR TRANSFORMS UNFIX UNSET VALENCE VECTOR˙ INTEGRAND VERTEX,VERTICES VERTEX COUNT VIEW 4D VIEW MATRIX VIEW TRANSFORMS VIEW TRANSFORM GENERATORS VOLCONST VOLUME WHERE WHILE WHITE WRAP WULFF Xn X Y YELLOW Z ZOOM ZOOM RADIUS ZOOM VERTEX

§5.1.10. Colors.

The colors of edges and facets are recorded as integers. How these integers translate to colors on the screen is determined by how the graphics drivers are written. The following synonyms are supplied for the integers 0 through 15, and it is hoped that the graphics drivers will be written to display these correctly:

BLACK, BLUE, GREEN, CYAN, RED, MAGENTA, BROWN, LIGHTGRAY, DARKGRAY, LIGHTBLUE, LIGHTGREEN, LIGHTCYAN,

LIGHTRED, LIGHTMAGENTA, YELLOW, and WHITE. The special color value CLEAR (-1) makes a facet transparent. These tokens are simply translated to integer values wherever they occur, so these are reserved words.

§5.1.11. Expressions.

Variable expressions are used in constraint and boundary formulas and integrands. Constant expressions can be used wherever a real value is needed. Expressions are given in algebraic notation with the following tokens:

x1,x2,...  coordinates
x,y,z  coordinates
p1,p2  parameters for boundaries
constant  any number acceptable to atof()
G  current gravitational constant
identifier  adjustable constant
E, PI  special constants
+,-,*,/,%  arithmetic
=  treated as low-precedence -
(,)  grouping and functional notation
^n  raise to integer n power
**  raise to real power
? :  conditional expression, as in C language

functions:

sqr  square
sin,cos,tan  trig functions
acos,asin,atan inverse trig functions
sqrt  square root
log,exp  natural log, exponentiation base e
pow  pow(x,y): raise x to real power y
usrn  user-defined function

User-defined functions can be defined in C in userfunc.c. They are meant for situations where expression interpretation is too slow, or functions such as elliptic integrals are wanted. Currently, they are automatically functions of the coordinates. Do not give any arguments in the expression; for example “(usr1

+ usr3)/usr10”.

Expressions are parsed into evaluation trees, which are interpreted when needed. Constant subtrees are folded into constant nodes at parse time.

Constant expressions must evaluate to values when they are read, i.e. they have no parameters or coordinates.

In parsing an expression, the longest legal expression is used. This permits coordinates to be specified by several consecutive expressions with no special separators.

NOTE: A ’+’ or ’-’ preceded by whitespace and followed by a number is taken to be a signed number. Thus “3 - 5” and “3-5” are single expressions, but “3 -5” is not. This is for convenience in separating multiple expressions listed on the same line for vertex coordinates, metric components, etc.

The mod operator ‘%’ does a real modulus operation:

x%y = x − floor(x/y) ∗ y.

§5.2. Datafile organization.

Surface Evolver Manual

The datafile is organized into six parts:

1. Definitions and options

2. Vertices

3. Edges

4. Faces

5. Bodies

6. Commands In the syntax descriptions below, keywords will be in upper case. const expr means a constant expression,

and expr means any expression. n or k means an integer, which may be signed if it is being used as an oriented element label. Square brackets denote optional items. ’|’ means ’or’. §5.2.1. Definitions and options.

Each line starts with a keyword. The order is immaterial, except that identifiers must be defined before use and quantities must be defined before referring to them in constraints. None of these are required, but those marked as default will be assumed unless an overriding option is present.

1. Macros.

#DEFINE identifier string See Macros section above.

2. Adjustable constants.

PARAMETER identifier [= const expr] This declares identifier to be an adjustable constant with the given initial value. The default initial value

is 0. The value may be changed with the A command. Adjustable constants may be used in any subsequent expression. Changing these results in automatic recalculation of the surface, unless AUTORECALC has been toggled off.

3. Dimensionality.

    STRING

The Evolver is to use the string model (see the String Model section)

SOAPFILM (default) The Evolver is to use the standard two-dimensional surface model.

    SURFACE_DIMENSION const expr

Surface is of given dimension. Dimension over 2 is valid only in the simplex model.

4. Domain.

    EUCLIDEAN (default)

The surface lives in Euclidean space.

    SPACE_DIMENSION const expr

The surface lives in Euclidean space of the given dimension. Default is 3. The dimension must be at most the value of MAXCOORD in model.h, which is 4 in the distributed version.

    TORUS

The surface lives in a 3 dimensional flat torus. Surface is equivalent to a periodic surface.

        TORUS_FILLED

Indicates the entire torus is filled by bodies to enable the program to avoid degenerate matrices when adjusting constrained volumes.

        PERIODS

expr expr expr expr expr expr expr expr expr

Used with TORUS domain. Specifies that the three side vectors (basis vectors) of the unit cell parallelpiped follow on the next three lines. Each vector is given by its three components. Default is unit cube. Adjustable parameters may be used in the expressions, so the fundamental domain may be changed interactively with either the ’A’ command or by assigning new values to the parameters. Be sure to do a ‘recalc’ to get the period matrix re-evaluated.

        SYMMETRY_GROUP "name"

The domain is the quotient of R^{n} by a user-defined symmetry group. The user must link in C functions that define the group operations. See quotient.c for an example. “name” is a double-quoted name that is checked against a name in quotient.c to be sure the right symmetry group is linked in.

5. Representation.

LINEAR (default)

Facets are flat triangles.

        QUADRATIC

Facets are quadratic patches. See the Quadratic Representation section.

        SIMPLEX_REPRESENTATION

Facets are defined by oriented vertex list rather than edge list. See the section above on Simplex Model.

6. Surface tension energy. The surface tension energy is always included in the total energy by default. It can be turned off only by giving the facets (or edges in the string model) the “density 0” attribute.

AREA (default)

Surface energy of a facet is its area times its surface tension.

        WULFF "filename"

Specifies that a crystalline integrand is to be used. The next token should be a double-quoted filename (with path) of a file giving the Wulff vectors of the integrand. The format of the file is one Wulff vector per line with its three components in ASCII decimal format separated by spaces. The first blank line ends the specification. Some special integrands can be used by giving a special name in place of the file name. Currently, these are "hemisphere" for a Wulff shape that is an upper unit hemisphere, and "lens" for two unit spherical caps of thickness 1/2 glued together on a horizontal plane. These two don’t need separate files.

        PHASEFILE "filename"

This permits the surface tension of a grain boundary to depend on the phases or types of the adjacent grains. The information is read from an ASCII file. The first line of the file has the number of different phases. Each line after consists of two phase numbers and the surface tension between them. If a pair of phases is not mentioned, the surface tension between them is taken to be 1.0. Facets in the string model or bodies in the soapfilm model can be labelled with phases with the PHASE phrase.

7. Surface energy integrals.

    SURFACE_ENERGY n [GLOBAL]

E1: expr E2: expr E3: expr

This defines surface energy integral n. n must be between 0 and 31 (15 for DOS). See facet attributes and the Energy section. GLOBAL means this energy applies to all facets.

8. Squared mean curvature.

    SQUARE_CURVATURE const expr

This phrase indicates that the integral of squared mean curvature will be included in the energy with a weight given by const expr. The weight can be changed with the A command by changing the value of the

adjustable constant square curvature modulus.

9. Integrated mean curvature.

    MEAN_CURVATURE_INTEGRAL const expr

This phrase indicates that the integral of mean curvature will be included in the energy with a weight given by const expr. The weight can be changed with the A command by changing the value of the adjustable

constant mean curvature modulus.

10. Squared Gaussian curvature.

    SQUARE_GAUSSIAN_CURVATURE const expr

This phrase indicates that the integral of squared Gaussian curvature will be included in the energy with a weight given by const expr. The weight can be changed with the A command by changing the value of the

adjustable constant square gaussian modulus.

11. Ideal gas model.

PRESSURE const expr

This specifies that bodies are compressible and the ambient pressure is the given value. The default is that bodies with given volume are not compressible.

12. Gravity.

    GRAVITY_CONST const expr

Specifies gravitational constant G. Default 1.0.

13. Gap energy.

    GAP_CONSTANT const expr

Multiplier for convex constraint gap energy. Default 1.0.

14. Knot energy.

    INSULATING_KNOT_ENERGY const expr

The total energy will include the insulating knot energy with multiplier const expr. Abbreviation for a named quantity.

    CONDUCTING_KNOT_ENERGY const expr

The total energy will include the conducting knot energy with multiplier const expr. Abbreviation for a named quantity.

15. Mobility and motion by mean curvature.

    GRADIENT_MODE (default)

Velocity is equal to the force.

        AREA_NORMALIZATION

The velocity of a vertex is the force divided by 1/3 area of neighboring facets (or 1/2 length of neighboring edges in string model) to approximate motion by mean curvature.

        APPROXIMATE_CURVATURE

Calculates vertex velocity from force by means of polyhedral linear interpolation inner product. Do not use

together with AREA_NORMALIZATION.
        EFFECTIVE_AREA

For both area normalization and approximate curvature modes, the resistance to motion is only on the component of velocity normal to the surface.

16. Annealing.

        JIGGLE

Continuous jiggling. Default none.

        TEMPERATURE const expr

Gives temperature for jiggling. Default 0.05.

17. Diffusion.

        DIFFUSION const expr

Specifies that diffusion between bodies is in effect with diffusion constant as given. Default 0.

18. Quantities.

        QUANTITY n [GLOBAL]  [FIXED = value]

Q1 expr Q2 expr Q3 expr

This defines quantity number n, where n is between 0 and 31, inclusive. GLOBAL means that it applies to all facets; otherwise, the facet definitions have to specify which quantities they contribute to. If FIXED is present, then the quantity acts as a mathematical constraint on the surface, with the quantity being forced to the given value. The three expressions following are the components of the vector integrand over facets. The integral is done with the original orientations of the edges. See the Constraints section for more on specifying quantities.

19. Named method instances.

METHOD_INSTANCE name METHOD methodname [ MODULUS constexpr] [ GLOBAL] parameters

This is an explicit definition of a named method instance for a named quantity. Such an explicit definition is useful if your quantity has several method instances, and you want to see their individual values or apply the instances to different elements. The modulus multiplies the method value to give the instance value. The default modulus is 1. GLOBAL makes the method apply to all elements of the appropriate type. Non-global instances may be applied to elements individually.

Each method may have various parameters to specialize it to an instance. Currently the only parameters specified here are scalar integrands, with syntax

    SCALAR_INTEGRAND: expr

Some methods use global parameters (a holdover that will be done away with eventually.) The instance definition does not have to be on one line.

The current methods, their element types, and parameters are listed here:

1. edge tension or edge length: Length of edge. Type: edge. Parameters: none.

2. facet tension or facet length: Area of facet. Type: facet. Parameters: none.

3. edge scalar integral: Integral of a scalar function over arclength. Type: edge. Parameters: scalar integrand.

4. facet scalar integral: Integral of a scalar function over facet area. Type: facet. Parameters: scalar integrand.

5. mean curvature integral: Integral of signed scalar mean curvature. Type: edge. Parameters: none.

6. gauss curvature integral: Integral of gaussian curvature. Type: facet. Parameters: none.

7. hooke energy: Energy of edges as springs with reference to a global equilibrium length. Type: edge. Global parameters: hooke power, hooke length.

8. local hooke energy: Energy of edges as springs with equilibrium length being average of lengths of neighbor edges. Type: vertex. Parameters: none.

9. knot energy: An “electrostatic” energy in which vertices are endowed with equal charges. Power law of potential is adjustable. Type: vertex. Global parameters: knot power.

10. uniform knot energy or edge knot energy: A knot energy where vertex charge is proportional to neighboring edge length. Type: vertex. Global parameters: knot power.

11. uniform knot energy normalizer: Supposed to approximate the part of uniform knot energy that

is singular in the continuous limit. Type: vertex. Global parameters: knot power.

12. edge edge knot energy: Between pairs of edges, energy is inverse square power of distance between midpoints of edges. (by John Sullivan) Type: edge. Parameters: none.

13. facet knot energy: Charge on vertex is proportional to area of neighboring facets. Type: vertex. Global parameters: knot power.

14. buck knot energy: Energy between pair of edges given by formula suggested by Greg Buck. Type: edge. Global parameters: knot power.

15. proj knot energy: Only the normal component of force between vertices is used to desingularize energy. Type: edge. Parameters: knot power.

16. circle knot energy: Another weird way to calculate a nonsingular energy between pairs of edges. Type: edge. Parameters: none.

17. sphere knot energy: Another weird way to calculate a nonsingular energy between pairs of facets. Type: facet. Parameters: none.

18. curvature function: Calculates forces as function of mean and gaussian curvatures at vertices. Function may be changed by user by altering teix.c. No energy, just forces. Type: vertex. Parameters: none.

19. average crossings: To calculate the average crossing number in all projections of a knot. (by John Sullivan) Type: edge. Parameters: none.

20. gravity method: For future use in calculating gravitational energy. Type: edge. Global parameters: gravity constant. Other parameters: body densities.

20. Named quantities.

    QUANTITY name ENERGY| FIXED = value| INFO_ONLY  MODULUS constexpr methodlist

These are an effort to provide a more systematic way of adding new energies and constraints. A “method” is a way of calculating a scalar value from some particular type of element (vertex, edge, facet, body). A “quantity” is the sum total of various methods applied to various elements, although usually just one method is involved. The name is an unquoted identifier. Any quantities may be declared to be one of three types: 1) energy quantities are added to the overall energy of the surface; 2) fixed quantities that are constrained to a fixed target value (by a single Newton step at each iteration); and 3) information quantities whose values are merely reported to the user.

Each quantity has a “modulus”, which is just a scalar multiplier of the whole quantity. A modulus of 0 will turn off an energy quantity. The default modulus is 1. Adding a new method involves writing C routines to calculate the value and the gradient as a function of vertex coordinates, and adding a structure to the method name array in quantity.c.

The quantity definition may contain one or more method instances. To incorporate a previously explicitly defined instance, include

METHOD instancename To instantiate a method in the quantity definition, you essentially incorporate the instance definition, but without an instance name:

METHOD methodname [ MODULUS constexpr] [ GLOBAL] parameters See the previous subsection for method details. Usually the second, implicit definition will be more convenient, as there is usually only one method per quantity. If GLOBAL_METHOD is used instead of GLOBAL, then the method is applied to all elements of the proper type. It is equivalent to using the GLOBAL keyword in the method specification. Nonglobal instances must be applied individually to elements. That is done by simply adding the quantity or instance name to the line defining the element. If a quantity name is used, then all instances of that quantity of the appropriate type are applied to the element.

Quantity values can be seen with the Q or A command, and the A command can be used to change the target value of a fixed quantity, the modulus of a quantity, or some parameters associated to various methods.

Examples: Hooke energy:

quantity hooke ENERGY modulus 10 global method hooke energy

A little sample datafile:

// test of quantity integrands on constraints
method_instance len method edge_length
quantity lenny info_only method len
quantity sam info_only method facet_scalar_integral
scalar_integrand z

vertices
1   0 0 1 fixed
2   1 1 2 fixed
3   0 1 3 fixed
edges
1   1 2 lenny
2   2 3
3   3 1 len
faces
1   1 2 3 sam

The sample datafile knotty.fe contains more examples.

21. Constraints.

        CONSTRAINT n [GLOBAL]  [CONVEX]  [NONNEGATIVE |  NONPOSITIVE]
        EQUATION | FORMULA | FUNCTION expr
        [ ENERGY

E1 expr E2 expr E3 expr]

        [ CONTENT

C1 expr C2 expr C3 expr]

        [ QUANTITY n

Q1 expr Q2 expr Q3 expr]

        [ QUANTITY ...]

This defines constraint number n, where n is an integer between 0 and MAXCONSTR (currently 31, except 15 in MS-DOS). GLOBAL means the constraint automatically applies to all vertices (but not automatically to edges or faces). GLOBAL constraints count in the number limit. If CONVEX is given, then an additional gap energy is attributed to edges on the constraint to prevent them from trying to short-circuit a convex boundary; see the Constraints and Energy sections above and the k command. If NONNEGATIVE or NONPOSITIVE is given, then all vertices will be forced to conform appropriately to the constraint at each iteration. The EQUATION expression defines the zero level set which is the actual constraint. It may be written as an equation, since ’=’ is parsed as a low-precedence minus sign. Do not use ’>’ or ’<’ to indicate inequalities; use NONNEGATIVE or NONPOSITIVE and an expresson. Conditional expressions, as in C language, are useful for defining constraints composed of several surfaces joined smoothly, such as a cylinder with hemispherical caps.

NOTE: One-sided constraints can cause the optimal scale factor algorithm to misbehave. It may be necessary to use a fixed scale factor. See the m command below.

IMPORTANT NOTE: Do not let two constraints with the same formula apply to a vertex; that leads to a singular matrix inversion when trying to project the vertex onto the constraints. For example, do not have a vertex subject to X1 = 0 and also have a global X1 NONNEGATIVE.

ENERGY signifies that vertices or edges on the constraint are deemed to have an energy. In the SOAPFILM model, the next three lines give components of a vectorfield that will be integrated along each edge on the constraint. In the STRING model, just one component is needed, which is evaluated at each vertex on the constraint. The main purpose of this is to permit facets entirely on the constraint to be omitted. Any energy they would have had should be included here. One use is to get prescribed contact angles at a constraint. This energy should also include gravitational potential energy due to omitted facets. Integrals are not evaluated on FIXED edges.

CONTENT signifies that vertices ( STRING model ) or edges ( SOAPFILM model ) on the constraint contribute to the area or volume of bodies. If the boundary of a body that is on a constraint is not given as facets, then the body volume must get a contribution from a content integral. It is important to understand how the content is added to the body in order to get the signs right. The integral is evaluated along the positive direction of the edge. If the edge is positively oriented on a facet, and the facet is positively oriented on a body, then the integral is added to the body. This may wind up giving the opposite sign to the integrand from what you think may be natural. Always check a new datafile when you load it to be sure the integrals come out right.

QUANTITY n means that edges subject to this constraint contribute to quantity n. The three expressions that follow are the components of the integrand for the edges. The integral is done with the original orientations of the edges. Any number of quantities may be included with a constraint, as long as they have been previously defined.

        CONSTRAINT_TOLERANCE const expr

This is the tolerance within which a vertex is deemed to satisfy a constraint. Default 1e-12.

        AREA_FIXED const expr

This makes the total area of the surface a constraint with the given value. May be changed with the “ area_fixed := value” command.

22. Boundaries.

        BOUNDARY n PARAMETERS k [CONVEX]

X1 expr X2 expr X3 expr

This defines boundary number n, where n is an integer between 1 and MAXBDRY (currently 20) and k is the number of parameters (1 or 2). If CONVEX is given, then an additional energy is attributed to edges on the boundary to prevent them from trying to short-circuit a convex boundary; see the k menu option below. The following three lines have the functions for the three coordinates in terms of the parameters P1 and maybe P2. Energy and content integrals for boundaries are not implemented.

23. Integration order.

        INTEGRAL_ORDER n

Edge integrals are done by n-point gaussian quadrature. This give exact values for polynomials of degree 2n − 1. . Default is n = 2, which will do cubic polynomials exactly. No limit on order, as weights and abscissas are calculated by Evolver. Facet integrals are done with 1, 3, or 7 point integration, depending on whether n ≤ 1, n ≤ 3, or n > 3.

24. Scale factor.

        SCALE const expr [FIXED]

Sets the initial scale factor. If FIXED is present, sets fixed scale factor mode. Default is scale = 0.1 and optimizing mode.

        SCALE_LIMIT const expr

Sets upper bound on scale factor to prevent runaway motions. Default value is 1. If you use surface tensions and densities not near unity, you may have to set this value.

25. Merit factor.

        MERIT_FACTOR const expr

When the user requests a report of information (command i), Evolver also reports the area/volume merit factor, which is (total area)^{3}/(total volume)^{2} multiplied by this factor. It is useful for comparing different configurations in seeking to minimize area for given volume, as in beating Kelvin’s tetrakaidecahedron.

26. Mobility.

        MOBILITY_TENSOR

expr expr expr expr expr expr expr expr expr

or

        MOBILITY expr

The force vector is multiplied by the mobility scalar or tensor to get the velocity. Good for, say, having grain boundary mobility depend on temperature.

27. Metric.

METRIC

expr expr expr expr expr expr expr expr expr

or

CONFORMAL_METRIC

expr or

KLEIN_METRIC

The user may define a background metric for the string model only. The keyword METRIC is follwed by the N^{2} components of the metric tensor, where N is the dimension of space. The components do not have to obey any particular line layout; they may be all on one line, or each on its own line, or any combination. It is up to the user to maintain symmetry. A conformal metric is a scalar multiple of the identity matrix, and only the multiple need be given. A conformal metric will run about twice as fast. The Klein metric is a built-in metric for hyperbolic n-space modelled on the unit disk or ball.

28. Autochopping.

AUTOCHOP const expr This turns on autochopping of long edges. The constant is the maximum edge length.

29. Autopopping.

AUTOPOP

This turns on autopopping of short edges and improper vertices.

30. Total time.

TOTAL_TIME const expr

This permits setting the initial time of a surface evolving with a fixed scale. Used primarily when resuming from a dump file.

31. Runge-Kutta.

RUNGE_KUTTA const expr

This turns on doing iteration with the Runge-Kutta method.

32. Homothety scaling.

HOMOTHETY const expr

This turns on doing homothety scaling each iteration. The scaling is a uniform scaling from the origin to keep the total volume of all bodies at the given value.

33. Viewing matrix.

VIEW_MATRIX

const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr

For the specification of the initial viewing transformation matrix of the surface. The matrix is in homogeneous coordinates with translations in the fourth column. This matrix will be part of all dump files, so the view can be saved between sessions. This matrix only applies to internal graphics (Postscript, Xwindows, etc.) and not external graphics (geomview).

34. View transforms.

        VIEW_TRANSFORMS integer

[ color color]

const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr ...

For the display of several transformations of the surface simultaneously, a number of viewing transformation matrices may be given. The transforms apply to all graphics, internal and external, and are prior to the viewing transformation matrix for internal graphics. The identity transform is always done, so it does not need to be specified. The number of matrices follows the keyword VIEW_TRANSFORMS. Each matrix is in homogeneous coordinates. Individual matrices need no special separation; Evolver just goes on an expression reading frenzy until it has all the numbers it wants. Each matrix may be preceded by a color specification that applies to facets transformed by that matrix. The color applies to one transform only; it does not continue until the next color specification. Transforms may be activated or deactivated interactively with the transforms on or transforms off. See the next paragraph for a more sophisticated way to control view transforms.

35. View transform generators.

        VIEW_TRANSFORM_GENERATORS integer

const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr const expr ...

Listing all the view transforms as in the previous paragraph is tedious and inflexible. An alternative is to list just a few matrices that can generate transforms. See the transform expr command for instructions on entering the expression that generates the actual transforms.

§5.2.2. Vertex list.

The vertex list is started by the keyword VERTICES at the start of a line. It is followed by lines with one vertex specification per line in one of these formats:

k x y z [FIXED] [ CONSTRAINT c1 [ c2]] [BARE] [ quantityname ...] [ methodname ...] k p1 [ p2] [FIXED] BOUNDARY b [BARE] [ quantityname ...] [ methodname ...]

Here k is the vertex number, a positive integer. Vertices do not need to be listed in order, and there may be gaps in the numbering. However, if they are not in consecutive order, then the numbering in dump files will be different. x, y, and z are constant expressions for coordinates in the domain; p1 and p2 are constant expressions for parameter values. If FIXED is given, then the vertex never moves, except possibly for an initial projection to constraints. If CONSTRAINT is given, then one or two constraint numbers must follow. (Actually, you can list as many constraints as you want, as long as they are consistent.) The given coordinates need not lie exactly on the constraints; they will be projected onto them.

If BOUNDARY is given, then the boundary parameter values are given instead of the coordinates. The vertex coordinates will be defined by the coordinate formulas of boundary number b. A vertex may be on only one boundary.

The BARE attribute is just an instruction to the checking routines that this vertex is not supposed to have an adjacent facet in the soapfilm model, so spurious warnings will not be generated. This is useful when you want to show bare wires or outline fundamental domains.

An arbitrary number of named quantities or method instances may be listed. These add method values for this vertex to the named quantity. The named quantity or instance must have been declared in the top section of the datafile. See the Named quantity section above.

§5.2.3. Edge list.

The edge list is started by the keyword EDGES at the start of a line. It is followed by lines with one edge specification per line in this format (linespliced here):

        k v1 v2 [ s1 s2 s3] [FIXED]  [BOUNDARY b]  [CONSTRAINTS c1 [ c2]]\
        [TENSION |  DENSITY const expr]  [COLOR n]  [BARE] \

[ quantityname ...] [ methodname ...]

Here k is the edge number, with numbering following the same rules as for vertices. v1 and v2 are the numbers of the tail and head vertices of the edge. In a TORUS model, there follow three signs s1 s2 s3 indicating how the edge wraps around each unit cell direction: + for once positive, * for none, and - for once negative. FIXED means that all vertices and edges resulting from subdividing this edge will have the FIXED attribute. Likewise the BOUNDARY and CONSTRAINT attributes will be inherited by all edges and vertices derived from this edge. If a constraint has energy or content integrands, these will be done for this edge. IMPORTANT: If a constraint number is given as negative, the edge energy and content integrals will be done in the opposite orientation. In the string model, the default tension is 1, and in the soapfilm model, the default tension is 0. However, edges may be given nonzero tension in the soapfilm model, and they will contribute to the energy.

If the simplex model is in effect, there should be no EDGES section.

The BARE attribute is just an instruction to the checking routines that this ede is not supposed to have an adjacent facet in the soapfilm model, so spurious warnings will not be generated. This is useful when you want to show bare wires or outline fundamental domains.

An arbitrary number of named quantities or method instances may be listed. These add method values for this edge to the named quantity. The named quantity or instance must have been declared in the top section of the datafile. See the Named quantity section above.

§5.2.4. Face list.

The face list is started by the keyword FACES at the start of a line. It is followed by lines with one facets specification per line in this format:

        k e1 e2 ... [FIXED]  [BOUNDARY b]  [CONSTRAINTS c1 [ c2]]  [QUANTITY n]\
        [ENERGY n]  [TENSION |  DENSITY const expr]  [NODISPLAY]  [TAG const expr]\
        [COLOR n]  [FRONTCOLOR n]  [BACKCOLOR n]  [PHASE n] \

[ quantityname ...] [ methodname ...]

Here k is the face number, with numbering following the same rules as for vertices. There follows a list of oriented edge numbers in counterclockwise order around the face. A negative edge number means the opposite orientation of the edge from that defined in the edge list. The head of the last edge must be the tail of the first edge (except if you’re being tricky in the STRING model). There is no limit on the number of edges. The face will be automatically subdivided into triangles if it has more than three edges. QUANTITY n means that quantity integrand n applies to this facet. Actually, any number of quantities can apply; their numbers just need to be listed after QUANTITY or QUANTITIES. ENERGY n means that surface energy integrand n applies to this facet. Actually, any number of energies can apply; their numbers just need to be listed after ENERGY or ENERGIES. The TENSION or DENSITY value is the energy per unit area (the surface tension) of the facet; the default is 1. Density 0 facets exert no force, and can be useful to define volumes or in displays. Fractional density is useful for prescribed contact angles. NODISPLAY prevents the facet from being displayed. The TAG number is inherited by refined facets; it is useful for refering to groups of facets with the show query. The COLOR attribute applies to both sides of a facet; FRONTCOLOR applies to the positive side (edges going counterclockwise) and BACKCOLOR to the negative side. The PHASE number is used in the string model to determine the surface tension of edges between facets of different phases, if phases are used. The FIXED, BOUNDARY, CONSTRAINT, DENSITY and NODISPLAY attributes will be inherited by all facets, edges, and vertices resulting from the subdivision of the interior of the face.

If the simplex model is in effect, the edge list should be replaced by an oriented list of vertex numbers. An arbitrary number of named quantities or method instances may be listed. These add method values

for this facet to the named quantity. The named quantity or instance must have been declared in the top section of the datafile. See the Named quantity section above.

§5.2.5. Bodies.

The body list is started by the keyword BODIES at the start of a line. It is followed by lines with one body specification per line in this format:

        k f1 f2 f3 .... [VOLUME v]  [VOLCONST v]  [PRESSURE p]  [DENSITY d] \

[PHASE n] [ quantityname ...] [ methodname ...]

Here k is the body number, and f1 f2 f3 ... is an unordered list of signed facet numbers. Positive sign indicates that the facet normal (as given by the right-hand rule from the edge order in the facet list) is outward from the body and negative means the normal is inward. Giving a VOLUME value v means the body has a volume constraint, unless the ideal gas model is in effect, in which case v is the volume at the ambient pressure. VOLCONST is a value added to the volume; it is useful when the volume calculation from facet and edge integrals differs from the true volume by a constant amount. Giving a PRESSURE value p means that the body is deemed to have a constant internal pressure p; this is useful for prescribed mean curvature problems. It is incompatible with prescribed volume. Giving a DENSITY value d means that gravitational potential energy (gravitational constant G) will be included. v, p, and d are constant expressions.

To endow a facet with VOLUME, PRESSURE, or DENSITY attributes in the STRING model, define a body with just the one facet.

The PHASE number is used in the soapfilm model to determine the surface tension of facets between bodies of different phases, if phases are used.

An arbitrary number of named quantities or method instances may be listed. These add method values for this body to the named quantity. The named quantity or instance must have been declared in the top section of the datafile. See the Named quantity section above.

The BODIES section is optional.

§5.2.6. Commands.

Encountering the keyword READ in the datafile causes the Evolver to switch from datafile mode to command mode and read the rest of the datafile as command input. This feature is useful for automatic initialization of the surface with refining, iteration, defining your own commands, etc.

WARNING: Do not put the READ command in an #include file and then put more commands after the #include in the datafile. The latter set of commands will not be read.

Chapter 6. Operation.

§6.1. System command line.

Syntax:

evolver [-ffilename] [-d] [-x] datafile

The name of a datafile is a required argument. The current directory and EVOLVERPATH will be searched. If the datafile is not found, then a new search with extension “.fe” is done.

Options: -d Prints YACC debugging trace as datafile is parsed. May be helpful if you can’t figure out why your

datafile doesn’t get read in properly AND you know YACC. -f Specifies the name of a file to be used as command input. At the end of this file, input reverts to stdin.

The effect is the same as redirecting input from the file, except that -f will echo commands to the screen and revert to stdin at the end.

-t On NeXTs, sets up a terminal interface instead of menus and dialog boxes. Hence the interface is much

more like that on other machines. -u On NeXTs, prevents menus and displays. Pure tty interface. Meant for running on a NeXT from a

remote site. -x Causes Evolver to exit whenever an error occurs. Meant to be used when Evolver is run in a shell script,

so that Evolver doesn’t get hung up waiting for input in response to an error.

§6.2. Initialization.

The following steps occur when a new datafile is read in:

1. Any previous surface has all memory deallocated.

2. Defaults are initialized.

3. Datafile is read in.

4. Any non-triangular faces are divided into triangles.

5. The order of facets around each edge is determined geometrically.

6. Vertices are projected to their constraints.

7. Checks as described under the C command are run.

8. Initial areas, energies, and volumes are calculated.

9. The viewing transformation matrix is reset to center the surface on the screen.

10. The main command interface is started.

§6.3. Error handling.

There are several categories of errors:

1. WARNING. Something has happened that you should know about, but it is possible to proceed.

2. EXPRESSION ERROR. There is an error in parsing an expression.

3. PARSING ERROR. Error in datafile syntax, but parsing will continue as best it can.

4. DATAFILE ERROR. Error in datafile semantics, but parsing will continue as best it can.

5. ERROR. The current operation cannot continue. It is abandoned and you return to the command prompt.

6. FATAL ERROR. The program cannot continue. Exits immediately.

All error messages go to stderr. Errors in the datafile report line numbers and print the current line so far. Note that the real error may have been at the end of the previous line. If command input is being taken from a file at the time an error occurs, the line number of the offending command will be printed and command input will revert to stdin. If the -x option was given when Evolver was started, then Evolver will exit immediately with a nonzero error code.

§6.4. Commands.

Commands are of two main types. The first type is one letter (case is significant). These perform simple and common actions. Single-letter commands are listed below in functional groups. The second type is composed of spelled-out words, and is a sort of combination SQL type query language and a programming language.

Commands may be read from the end of the datafile, from a file given on the system command line, from stdin (the terminal), or from a file given in a READ command. The interactive command prompt is

“Enter command: ”.

§6.4.1. General language syntax.

Commands are entered one line at a time, parsed, and executed. Multi-line commands may be entered by enclosing them in braces. If a line ends while nested in braces or parenthesis, Evolver will ask for more input. It will also ask for more if the line ends with certain tokens (such as ‘+’) that cannot legally end a command. Explicit continuation to the next line may be indicated by ending a line with a backslash (linesplicing). You may want to use the read command to read long commands from a file.

Successful commands are saved in a history list. They may be accessed with !! for the last command or !string for the latest command with matching initial string. The command will be echoed. Only the first line of a command is saved. It is advised that longer commands be defined as new command words instead of being executed from the history list (see the Assignment section below).

Some single-letter commands require interactive input. For those, there are equivalent commands listed below that take input information as part of the command. This is so commands may be read from a script without having to put input data on additional lines after the command, although that can still be done for the single-letter versions.

General note: Some commands will prompt you for a value. A null response (just RETURN) will leave the old value unchanged and return you to the command prompt. On options where a zero value toggles OFF, the zero must be explicitly entered. Commands that need a real value will accept an arbitrary expression.

Many commands that change the surface or change the model will cause energies and volumes to be recalculated. If you suspect a command has not done this, the recalc command will recalculate everything. It will also update any automatic display.

General control structures.

command ; command ...

Several commands on the same line may be separated by a semicolon. A semicolon is not needed after the last command. Example:

    g 10; r; g 10; u

{command ; ... }

Curly braces group a list of commands into one command. This is useful in the various control structures. Note that ‘} ’ is not a command separator. You still need a semicolon to separate a command list from another command.

command integer

For execution of a command a fixed number of times. Be sure to leave a space between a single-letter command and the integer lest your command be interpreted as one identifier. To avoid several types of confusion, only certain types of commands are repeatable:

1. single letter commands that don’t have optional arguments (l,t,j,m,n,w have optional arguments)

2. command list in braces

3. user-defined procedure names

command | stringexpr

The output of the command is piped to the system command. The stringexpr needs to be a quoted string or a string variable. It is interpreted as a system command. Redirection with ‘>’ instead of ’—’ is not available due to the use of ‘>’ as an comparison operator. If you want output to go to a file, use "tee file" or use ordinary redirection in the system command, "cat >file". On systems (such as DOS and Macs) without true piping, the string is interpreted as a filename, so output goes directly to a file instead of a command.

Examples:

list facets | "more"
list vertices | "tee vlist" ; g 10

IF expr THEN command [ ELSE command ] For conditional execution of commands. expr is true if nonzero. Example:

if max(edges,length) > 0.02 then {r; g 100 } else g 4

WHILE expr DO command DO command WHILE expr

For iterated execution of command controlled by a logical expression. Expression is true if nonzero. Example:

while max(edges,length) > 0.02 do {r; {g 40; u} 5 }

§6.4.2. Single-letter command summary.

There are five groups of commands:

Reporting: C Run consistency checks. c Report count of elements. e Extrapolate. i Information on status. v Report volumes. z Do curvature test.

Model characteristics: A Set adjustable constants. a Toggle area normalization b Set body pressures. f Set diffusion constant. G Set gravity. J Toggle jiggling on every move. k Set boundary gap constant. M Toggle linear/quadratic model. m Toggle fixed motion scale. p Set ambient pressure. Q Report or set quantities.

U Toggle conjugate gradient method. W Homothety toggle.

Surface modification: g Go one iteration step. Often followed by a repetition count. j Jiggle once. K Skinny triangle long edge divide. l Subdivide long edges. N Set target volumes to actual. n Notch ridges and valleys. O Pop non-minimal edges. o Pop non-minimal vertices. r Refine triangulation. t Remove tiny edges. u Equiangulate. V Vertex averaging. w Weed out small triangles. y Torus duplication. Z Zoom in on vertex. 1 Minimization via Hessian.

Output: D Toggle display every iteration. d Dump surface to datafile. P Create graphics display file. s Show triangulation graphically.

Miscellaneous: F Toggle command logging.

H,h,? Help screen. q,x Exit.

The j, l,m,n,t, and w commands require real values, which may be entered on the same line, or given in response to a prompt if not.

§6.5. Alphabetical single-letter command reference.

A Adjustable constants. Displays current values and allows you to enter new values. New value is entered as the number of the constant (from the list) and the new value. Values of named quantities, their moduli, and the target values of fixed named quantities also appear here. The modulus and target value may be changed.

a Toggles area normalization of vertex forces and other gradient. Be sure you have a small enough scale factor, or else things tend to blow up. Reduce scale factor temporarily after refinement, since triangle areas are cut by a factor of 4 but the old creases remain. When this option is ON, there is an optional check that can be made for facets that move too much. This is done by computing the ratio of the length of the normal change to the length of the old normal. If this exceeds the user-specified value, then all vertices are restored to their previous position. The user should reduce the motion scale factor and iterate again.

b Permits user to change body prescribed volumes and pressures. Prints old value for each body and prompts for new.

C Runs various internal consistency checks. If no problems, just prints “Checks completed.” The checks are: Element list integrity - checks that data structures are intact. Facet-edge check - that if a facet adjoins an edge, then the edge adjoins the facet, and that the three edges around a facet link up. Facet-body check - whether adjacent facets have the same body on the same side. Collapsed elements - check if endpoints of an edge are the same, and whether neighboring facets share more than one edge and two vertices.

c Prints count of elements and memory used and prints various model parameters.

D Toggles displaying every iteration. Initially does not display.

d Dumps data to ASCII file in same format as initial data file. You will be prompted for a filename. An empty reponse will use the default dump name, which is the datafile name with a “.dmp” extension. Useful for checking your input is being read correctly, for saving current configuration, and for debugging.

e Extrapolates total energy to infinite refinement if at least two refinements have been done. Uses last energy values at three successive levels of refinement.

F Toggle logging of commands in file. If starting logging, you will be prompted for the name of a log file. Any existing file of that name will be overwritten. Logging stops automatically when the surface is exited. Only commands that change the surface are logged.

f Set diffusion constant. Prints old and prompts for new.

G Toggles gravity on or off. Gravity starts ON if any body has DENSITY; otherwise OFF. Prints old value of gravitational constant and prompts for new.

g Do iteration step. Each iteration calculates the force on each vertex due to its contribution to the total energy and moves the vertex by a multiple of the force. There is a global scale factor that multiplies the force to give the displacement. If area normaliztion is turned on, the force at each vertex is also divided by the total area of facets adjacent to the vertex to better approximate motion by mean curvature (but this seems to be often numerically badly behaved due to long skinny facets). If any bodies have prescribed volumes, the vertices are also displaced to bring the volumes back to near the prescribed values. If scale optimizing (see command m) is ON, then different global scale values will be tried until a quadratic interpolation can be done to find the optimal value. (This can blow up under certain circumstances.) The output consists of the number of iterations left (for people who wonder how close their 1000 iterations are to ending), the area and energy, and the scale factor. The user can abort repeated iterations by sending an interrupt to the process (SIGINT, to be precise; CTRL-C or whatever on your keyboard).

h,H,? Prints help screen listing these commands.

i Prints miscellaneous information:

Total energy Gap energy for convex constraints Total area of facets Count of elements and memory usage Area normalization, if on LINEAR or QUADRATIC facets Scale factor value and option (fixed or optimizing) Diffusion option and diffusion constant value Gravity option and gravitational constant value Jiggling status and temperature Gap constant (for gap energy) Ambient pressure (if ideal gas model in effect) Area-volume figure of merit, the ratio of energy cubed (usually energy is area) to volume squared. This ratio is dimensionless and useful in comparing the efficiencies of different shaped surfaces in partitioning space (trying to beat Kelvin’s tetrakaidecahedron).

J Toggles jiggling on every iteration. If jiggling gets turned on, prompts for temperature value.

j Jiggles all vertices once. Useful for shaking up surfaces that get in a rut, especially crystalline integrands. You will be prompted for a “temperature” which is used as a scaling factor, if you don’t give a temperature with the command. The actual jiggle is a random displacement of each vertex independently with a Gaussian distribution. See the longj command below for a user-definable perturbation.

K Finds skinny triangles whose smallest angle is less than a specified cutoff. You will be prompted for a value. Such triangles will have their longest edge subdivided. Should be followed with tiny edge removal (t) and equiangulation (u).

k Sets “gap constant” for gap energy for CONVEX boundaries. Adds energy roughly proportional to area between edge and boundary. You will be prompted for a value. Normal values are on the order of magnitude of unity. Value k = 1 is closest to true area. Use 0 to eliminate the energy.

l Subdivides long edges, creating new facets as necessary. You will be prompted for a cutoff edge length, if you don’t give a value with the command. Existing edges longer than the cutoff will be divided once only. Newly created edges will not be divided. Hence there may be some long edges left afterward. If you enter h, you will get a histogram of edge lengths. If you hit RETURN with no value, nothing will be done. It is much better to refine than to subdivide all edges.

M Sets model type to LINEAR or QUADRATIC. Changing from LINEAR to QUADRATIC adds vertices at the midpoints of each edge. Changing from QUADRATIC to LINEAR deletes the midpoints.

m Toggles quadratic search for optimal global motion scale factor. If search is toggled OFF, you will be prompted for a scale factor. If you give a value with the command, then you are setting a fixed scale factor. Values around 0.2 work well when the triangulation is well-behaved and area normalization (command a) is off. In optimizing mode, a scale factor getting small, say below 0.01, indicates triangulation problems. Too large a scale factor will show up as total energy increasing. If you have motion by area normalization ON (command a), use a small scale factor, like 0.001, until you get a feel for what works.

N Normalize body prescribed volumes to current actual volumes.

n Notching ridges and valleys. Finds edges that have two adjacent facets, and those facets’ normals make an angle greater than some cutoff angle. You will be prompted for the cutoff angle (radians) if you don’t give a value with the command. Qualifying edges will have the adjacent facets subdivided by putting a new vertex in the center. Should follow with equiangulation. In the string model, it will refine edges next to vertices with angle between edges (parallel orientation) exceeding the given value.

O Pop non-minimal edges. Scans for edges with more than three facets attached. Splits such edges into triple-facet edges. Splits propagate along a multiple edge until they run into some obstacle. Try octa.fe for an example.

o Pop non-minimal vertices. This command scans the surface for vertices that don’t have the topologies of one of the three minimal tangent cones. These are “popped” to proper local topologies. Not all types are properly treated as yet, and this command tends to create a lot of small edges and skinny triangles. Try popstr.fe and octa.fe for examples.

P Produce graphics output files. The view is the same as seen with the s command. Several formats are currently available, and a menu will appear. These are now Pixar, MinneView, geomview, PostScript, SoftImage, and a personal format called Triangle. If you are doing a torus surface, you will be asked for a display option, for which see the Torus section. You will be prompted for a filename. For Pixar format, you will be asked whether vertices should have normal vectors for normal interpolation (calculated as the average normal of all facets around a vertex); whether inner facets (adjacent to two bodies) outer facets (adjacent to zero or one body), or all facets are to be plotted; and whether different colors should be used for different bodies. If so, you are asked for the name of a file with a colormap in the format of RGB values, one set per line, values between 0 and 1. (This map may not be effective on all devices.)

You may also be asked if you want thickening. If you do, each facet will be recorded twice, with opposite orientations, with vertices moved from their original positions by the thickening distance (which the option lets you enter) in the normal direction. The normal used at each vertex is the same as used for normal interpolation, so all the facets around a planar vertex will have that vertex moved the same amount. Triple junctions will be separated. Thickening is good for rendering programs that insist on consistently oriented surfaces, or that can’t see the backside of a surface. The default thickening distance is one one-thousandth of the diameter of the object.

For file formats, see the section on Graphics file formats.

For those of you that have MinneView, the option for shared-memory MinneView is here. MinneView will start automatically, and will end automatically when you quit the Evolver. geomview is also here. geomview uses a named pipe interface at the moment, and may be slower than MinneView. There are options for starting a named pipe without invoking geomview. You will be told the name of the pipe, and it is up to you to start a pipe reader. Evolver blocks until a pipe reader is started. This is useful for having a second instance of Evolver feed a second surface to geomview by having geomview load the pipe.

p Sets ambient pressure in ideal gas model. Large p gives more incompressible bodies.

Q Report current values and optionally change target values of user-defined quantities.

q Exit program. You will be given a chance to have second thoughts. You may also load a new datafile. Automatically closes graphics if you really quit. Does not save anything. r Refines the triangulation. Edges are divided in two, and SOAPFILM facets are divided

into four facets with inherited attributes. Reports the number of structures and amount of memory used.

s Shows the surface on the screen. The proper display routines must have been linked in for the machine the display is on. Goes into the graphics command mode (see below). Torus surfaces have display options you will be asked for the first time, for which see the Torus section. The graphics window may be closed with the CLOSE_SHOW command.

t Eliminates tiny edges and their adjacent facets. You will be prompted for a cutoff edge length if you don’t give a value with the command. If you enter h, you will get an edge length histogram. If you hit RETURN without a value, nothing will happen. Some edges may not be eliminable due to being FIXED or endpoints having different attrtibutes from the edge.

U This toggles conjugate gradient mode, which will usually give faster convergence to the minimum energy than the default gradient descent mode. The only difference is that motion is along the conjugate gradient direction, and the maximum scale factor is not in effect. The scale factor should be in optimizing mode. The history vector is reset after every surface modification, such as refinement or equiangulation.

u This command, called “equiangulation”, tries to polish up the triangulation. Each edge that has two neighboring facets (and hence is the diagonal of a quadrilateral) is tested to see if switching the quadrilateral diagonal would make the triangles more equiangular. For a plane triangulation, a fully equiangulated triangulation is a Delauney triangulation, but the test makes sense for skew quadrilaterals in 3-space also. It may be necessary to repeat the command several times to get complete equiangulation.

V Vertex averaging. For each vertex, computes new position as area-weighted average of centroids of adjacent facets. Only adjacent facets with the same constraints and boundaries are used. Preserves volumes, at least to first order. See the rawv command below for vertex averaging without volume preservation, and rawestv for ignoring likeness of constraints. Does not move vertices on triple edges or fixed vertices. Vertices on constraints are projected back onto their constraints.

v Shows prescribed volume, actual volume, and pressure of each body. Also shows constrained quantities and information-only quantities. Pressures are really the Lagrange multipliers.

W Toggles homothety. If homothety ON, then after every iteration, the surface is scaled up so that the total volume of all bodies is 1. Meant to be used on surfaces without any constraints of any kind, to see the limiting shape as surface collapses to a point.

w Tries to weed out small triangles. You will be prompted for the cutoff area value if you don’t give a value with the command. If you enter h, you will get a histogram of areas to guide you. If you hit RETURN with no value, nothing will be done. Some small triangles may not be eliminable due to constraints or other such obstacles. The action is to eliminate an edge on the triangle, eliminating several facets in the process. WARNING: Although checks are made to see if it is reasonable to eliminate an edge, it is predicated on facets being relatively small. If you tell it to eliminate all above area 5, Evolver may eliminate your entire surface without compunction.

x Exit program. You will be given a chance to have second thoughts. You may also load a new datafile. Automatically closes graphics if you really quit. Does not save anything.

y Torus duplication. In torus model, prompts for a period number (1,2, or 3) and then doubles the torus unit cell in that direction. Useful for extending simple configurations into more extensive ones.

Z Zooms in on a vertex. Asks for vertex number and radius. Number is as given in vertex list in datafile. Beware that vertex numbers change in a dump (but correct current zoom vertex number will be recorded in dump). Eliminates all elements outside that distance from vertex 1. New edges at the radius are made FIXED. Meant to investigate tangent cones and intricate behavior, for example, where wire goes through surface in the overhand knot surface. Zooming is only implemented for surfaces without bodies.

z Do curvature test on QUADRATIC model. Supposed to be useful if you’re seeking a surface with monotone mean curvature. Currently checks angle of creases along edges and samples curvature on facet interiors. Orientation is with respect the way facets were originally defined.

§6.6. Command and query language.

The Evolver command language continues to grow by accretion, and it looks like it’s headed towards a full programming language. It has variables, expressions, subroutines, conditionals, and iteration constructs. But not subroutine arguments, local variables, structures, or arrays. Variables are either floating point, string, or subroutine names. Some longer examples of command scripts follow the language description.

Command syntax: Keywords are shown in upper case, although case is irrelevant in actual commands, except for single-letter commands. Square brackets enclose optional parts of commands.

Expressions:

Expressions are of two types: numeric and string. String expressions are either quoted strings or are created by the SPRINTF command. Numeric expressions are always floating-point, not integer. Boolean values for conditions are also floating point, 0 for false, nonzero for true. Numeric expressions are given in algebraic notation with the following terms and operators (in precedence order):

constants:

constant  an explicit number: integer, fixed point, or scientific notation
identifier  adjustable constant (listed by A command)
G  current gravitational constant
E, PI  special constants

element attributes:

X1,X2,...  coordinates
X,Y,Z  coordinates
P1,P2  parameters for boundaries
ID  element identifying number
ORIGINAL  number of parent datafile element
TAG  facet tag number from datafile
COLOR  integer representing color of facet, color of front if back different
FRONTCOLOR  color of front of facet
BACKCOLOR  color of back of facet
VALENCE  number of edges on a vertex, facets on edge, edges on a facet, or facets on a body
AREA  area of facet
LENGTH  length of edge
VOLUME  volume of body
DENSITY  density of edge, facet, or body
DIHEDRAL  dihedral angle of facets on an edge
ON CONSTRAINT n  test if element on given constraint
HIT CONSTRAINT n  test if element on a one-sided constraint
ON BOUNDARY n  test if element on given boundary

operators:

()  grouping and functional notation
^n  raise to integer n power
**  raise to real power
*,/,%  arithmetic
+,-  arithmetic

==,>,<,<=,>=,!=  comparison

NOT, !  logical NOT
AND, &&  logical AND
OR, ||  logical OR
?  :  conditional expression, as in C language

functions:  SQR, SQRT, SIN, COS, TAN, ACOS,

                                ASIN, ATAN, LOG, EXP, ABS, FLOOR, CEIL

and some internal read-only values:

vertex_count  number of vertices
edge_count  number of edges
facet_count  number of faces
body_count  number of bodies
total_energy  total energy of the surface
total_area  total area of the surface (film model)
total_length  total length of the surface (string model)
scale  current scale factor
space_dimension  dimension of ambient space
surface_dimension  dimension of surface
integration_order  points used in gaussian integration
torus  whether torus domain (Boolean)
scale  whether torus_filled specified in effect (Boolean)
symmetry_group  whether any symmetry active (Boolean)
simplex_representation  whether facets are represented as simplices (Boolean)

Also any toggle command may be used as a Boolean variable in an expression (full word toggles, not single letters).

NOTE: A ‘+’ or ‘-’ preceded by whitespace and followed by a number is taken to be a signed number. Thus “3 - 5” and “3-5” are single expressions, but “3 -5” is not. This is for convenience in separating multiple expressions listed on the same line in the datafile.

The mod operator ‘%’ does a real modulus operation:

x%y = x − floor(x/y) ∗ y,

but it works fine on integer values.

Elements inherit the original number of their parent element from the datafile; this can be referred to as original. New vertices and edges that subdivide facets have and original value of 0.

In queries and aggregate expressions, the iterating element may be named, and the name may be use to qualify an element attribute to clarify which element it refers to. The syntax is name.attribute. The following command uses a named element, although it is not necessary:

list facets ff where ff.area < 2

Element iterators. One feature different from ordinary C is the presence of iterators over geometric elements. These occur wherever an element type (vertices, edges, facets, bodies, singular or plural) appears in a command. Attributes of the iteration element may be used later in the command. The general form of an iterator is

elementtype [name] [ WHERE condition ] The iteration element may be named (useful in nested iterations), in which case the name must be used when the attribute is referenced. Element types used as attributes will just generate the elements associated with the parent element, so if ff is a facet, then ff.vertex will generate its three vertices. Currently implemented subelements are vv.EDGE, vv.FACET, ee.VERTEX, ee.FACET, ff.VERTEX, ff.EDGE, ff.BODY, and bb.FACET, where vv, ee, ff, and bb are the are names of vertices, edges, facets, and bodies respectively. But be sure to remember that in a nested iteration, an unqualified element type generates all elements of that type, not just those associated with the parent element. Also, each iterator can take a ”where” clause to limit the elements included.

Aggregate expressions. The maximum, minimum, sum, and average of an expression over a set of elements may be done with aggregate expressions. The general form is

aggregate(iterator ,expr)

where aggregate may be max, min, sum, or avg. The iterator generates elements as discussed above. Max and min over logical expressions may be used for logical any and all.

Example: This displays facet 4 and all its neighboring facets:

    show facets ff where max(ff.edge ee,max(ee.facet fff, fff.id == 4))

SQL-type queries on sets of elements:

These commands generally operate on a set of elements defined by an element iterator.

FOREACH iterator DO command

Repeat a command for each element generated by the iterator. Examples:

    foreach vertex do print x^2 + y^2 + z^2
    foreach edge ee where ee.dihedral > .4 do \
    { printf "id %g\n",id; foreach ee.vertex do printf " %g %g %g\n",x,y,z; }

LIST iterator

List elements on the screen in the same format as in the datafile. Piping to more can be used for long displays. Examples:

    list vertices | "more"
    list vertices where x < 1 and y > 2 and z >= 3 | "tee vfile"
    list edges where id == 12

REFINE iterator

Subdivides edges by putting a vertex in the middle of each edge and splitting neighboring facets in two. It is the same action as the long edge subdivide command (command l). Facets will be subdivided by putting a vertex in the center and creating edges out to the old vertices. It is strongly suggested that you follow this with equiangulation to nicen up the triangulation. Edge refinement is better than facet refinement as facet refinement can leave long edges even after equiangulation. Example:

    refine edges where not fixed and length > .1

DELETE iterator

Deletes edges by shrinking the edge to zero length (as in the tiny edge weed command t) and facets by eliminating one edge of the facet. Example:

    delete facets where area < 0.0001
DISSOLVE iterator

Removes elements from the surface without closing the gap left. The effect is the same as if the line for the element were erased from a datafile. Hence no element can be dissolved that is used by a higher dimensional element. Good for poking holes in a surface.

    dissolve facets where original == 2

SET elements [ name ] attrib expr [ WHERE expr ]

Sets an attribute to the indicated value for qualifying elements. The new value expression may refer to the element in question. SET can change the following attributes: constraint, coordinates, density, tag, target volume, pressure, color, frontcolor, backcolor, and opacity (for the appropriate type elements). For color, the color value is an integer. The integers from 0 to 15 can be referred to by the predefined names BLACK, BLUE,

GREEN, CYAN, RED, MAGENTA, BROWN, LIGHTGRAY, DARKGRAY, LIGHTBLUE, LIGHTGREEN, LIGHTCYAN, LIGHTRED,

LIGHTMAGENTA, YELLOW, and WHITE. The special color value CLEAR (-1) makes a facet transparent. On certain systems (Iris and Pixar at present) the opacity value sets the transparency of ALL facets. A value of 0.0 corresponds to transparent, and a value of 1.0 to opaque. Examples:

    set facets density 0.3 where original == 2
    set vertices x 3*x where id < 5
    set bodies volume 5 where id == 1
    set vertices constraint 1 where id == 4
    set facet color red where tag == 2
    set facet color clear where original < 5
    set facet opacity 0.5

SET elements [ name ] FIXED [ WHERE expr ] Sets the FIXED attribute for qualifying elements. Example:

    set vertices fixed where id == 4

SET name attrib expr Inside a FOREACH loop with a named element, one can set an attribute for the named element.

SET BACKGROUND color

Sets the backgound color for native graphics (not geomview). The color is an integer 0-15 or a color name listed in the general SET command below. Example:

    set background black

UNSET elements [ name ] attrib [ WHERE expr ]

Removes the attribute for qualifying elements. Attributes here are FIXED, VOLUME, PRESSURE, DENSITY, CONSTRAINT n, BOUNDARY n. A use for the last is to use a boundary or constraint to define an initial curve or surface, refine to get a decent triangulation, then use “unset vertices boundary 1” and ”unset edges boundary 1” to free the curve or surface to evolve.

UNSET elements [ name ] CONSTRAINT expr [ WHERE expr ] Removes constraint of given number from qualifying elements.

Single elements:

The verbs set, unset, list, delete, refine, fix, unfix, dissolve can be used on single named elements inside an iteration construct. Examples:

    foreach edge ee where ee.length < .01 do {
    printf "Refining edge %g.¨,ee.id; refine ee; }

Assignment:

Values can be assigned to variables. The variable names must be two or more letters, in order they not be confused with single-letter commands. Note that ’:=’ is used for assignment, not ’=’. Numeric, string, and command values may be assigned.

identifier := expr identifier ::= expr

Evaluates expression and assigns it to the named variable. If the variable does not exist, it will be created. These are the same class of variables as the adjustable parameters in the datafile, hence are all of global scope and may also be inspected and changed with the ’A’ command. If “::=” is used instead of “:=”, then the variable becomes permanent and will not be forgotten when a new surface is begun.

identifier := command identifier ::= command

Makes identifier an abbreviation for command. Subsequently, identifier may be used as a command itself. If “::=” is used instead of “:=”, then the variable becomes permanent and will not be forgotten when a new surface is begun. Example:

    gg :=  g 10; u
    gg 12
AMBIENT PRESSURE := expr

Sets the ambient pressure outside bodies.

AREA FIXED := expr

Sets the target value for the fixed area constraint. Does not activate constraint.

AUTOCHOP := expr

Sets autochop length, so edges longer than this will automatically be subdivided each iteration when autochopping is on. Does not toggle autochopping on.

COLORMAP integer := "filename" Sets the colormap file name.

DIFFUSION := expr

Sets diffusion constant. Does not toggle diffusion on.

GAP CONSTANT := expr

Sets the gap constant for assigning energy to gaps between facets and curved walls. Use 1 for best approximation to area.

GRAVITY := expr

Sets gravitational constant. Does not toggle gravity ON if it is OFF.

QUANTITY integer := expr

Sets the target value for a quantity. Does not automatically make it a constraint; use QUANTITY integer ON for that.

THICKEN := expr

Sets thickness for double layer display. Does not toggle it on.

Information:

LIST TOPINFO

Prints the first section of the datafile on the screen. This is everything before the vertices section.

PRINT expr

Prints the value of the expression. PRINT expr can also be used inside an expression, where it prints the expression and evaluates to the value of its expression. Examples:

    print max(edge,length)   print max(vertex, print (x^2+y^2+z^2) )

PRINT procedure name

Prints the commands assigned to the named procedure. These are reconstructed from the parse tree, so may not be identical in appearance to the original commands.

PRINTF string,expr,...

Prints to standard output using the standard C printf function. Only numeric expressions are valid, so the format string should use only %f and %g formats. The format string can be a string variable or a quoted string.

identifier := SPRINTF string,expr,...

Prints to a string variable using the standard C printf function. Only numeric expressions are valid, so the format string should use only %f and %g formats. The format string can be a string variable or a quoted string.

HISTOGRAM(element [ WHERE expr1 ], expr2 )

LOGHISTOGRAM(element [ WHERE expr1 ], expr2 )

Prints a histogram of the values of expr2 for the qualifying elements. The syntax is the same as the aggregate expressions. It uses 20 bins, and finds its own maximum and minimum values, so the user does not have to specify binsize. The log version will lump all zero and negative values into one bin. Example:

    histogram(edge,dihedral*180/pi)

General commands:

single-letter

Any of the single letter commands discussed in the last section may be used as a command. If commands are being read from a file, and if the command requires interactive input, each response should be on a separate line following the command.

AREAWEED expr

Weeds out facets smaller than given area. Same as ’w’ command, except does not need interactive response.

CLOSE SHOW

Closes the native graphics window started by the ‘s’ or SHOW commands. Does not affect geomview.

DIRICHLET

Does one iteration of minimizing the Dirichlet integral of the surface. The current surface is the domain, and the Dirichlet integral of the map from the current surface to the next. This is according to a scheme of Konrad Polthier and Ulrich Pinkall. At minimum Dirichlet integral, the area is minimized also. Works only on area with fixed boundary; no volume constraints or anything else. Seems to converge very slowly near minimum, so not a substitute for other iteration methods.

DUMP filename

Dumps current surface to named file in datafile format. With no filename, dumps to the default dump file, which is the datafile name with “.dmp” extension.

EDGEWEED expr

Weeds out edges shorter than given value. Same as ’t’ command, except does not need interactive response.

HESSIAN

Attempt to minimize energy in one step using Newton’s method with the Hessian matrix of the energy. Severe limitations at the present. Area can be the only energy, and fixed points and edges can be the only constraints. Even then, the surface is often not close enough to a minimum for the method to work; see the catenoid example. Does contain a safeguard to prevent the surface from blowing up.

LOAD filename

Terminates the current surface and loads a new datafile. The filename is the datafile name, and can be either a quoted string or a string variable.

LONGJ

Long jiggle once. This provides long wavelength perturbations that can test a surface for stability. The parameters are a wavevector, a phase, and a vector amplitude. The user will be prompted for values. Numbers for vectors should be entered separated by blanks, not commas. An empty reply will accept the defaults. A reply of r will generate random values. Any other will exit the command without doing a jiggle. In the random cases, a random amplitude A⃗ and a random wavelength L⃗ are chosen from a sphere whose radius is the size of the object. The wavelength is inverted to a wavevector w⃗. A random phase ψ is picked. Then each vertex v⃗ is moved by A⃗ sin(v⃗ · w⃗ + ψ).

LINEAR

Changes to linear model.

NOTCH expr

Notches ridges and valleys with dihedral angle greater than given value. Same as ’n’ command.

QUADRATIC

Changes to quadratic model.

QUIT

Exits Evolver or starts new datafile. Same as ‘q’ single-letter command.

RAWV

Does vertex averaging without conserving volume on each side of surface. Will only average vertices with those of like type of constraints. Doesn’t move vertices on boundaries.

RAWESTV

Does vertex averaging without regard for conserving volume or whether averaged vertices have like constraints. But doesn’t move vertices on boundaries.

READ filename Reads commands from the named file. The filename can be either a quoted string or a string variable.

REBODY

Recalculates connected bodies. Useful after a body has been disconnected by a neck pinching off. Facets of an old body are divided into edge-connected sets, and each set defines a new body (one of which gets the old body number). The new bodies inherit the attributes of the old body. If volume is fixed, then the target volume becomes the new actual volume. New bodies descended from an original body may be queried with the ‘original’ atttribute.

RECALC

Recalculates everything. Useful after changing some variable or something and recalculation is not automatically done.

SHELL

Invokes a system subshell for the user.

STABILITY TEST

Calculates the largest eigenvalue of the mobility matrix. Generates a random vector and applies mobility matrix 20 times. Prints out ratios of successive vector lengths. Useful in investigating stability when using approximate curvature.

SYSTEM systemcommand

Invokes a subshell to execute the given command. Command must be a quoted string or a string variable. Will wait for command to finish before resuming.

ZOOM [ integer expr ]

Zooms in on vertex whose id is the given integer, with radius the given expression. Same as ’Z’ command, but not interactive. If id and radius omitted, then uses current values.

Toggles:

Various features can be toggled on or off by giving the toggle name with ON or OFF. The default is ON. This is different from the single-letter command toggles, which always change the state. The toggle names below have brief descriptions of their actions in the ON state.

APPROX CURV

Use polyhedral curvature (linear interpolation over facets) for mean curvature vector. Actually establishes the inner product for forms or vectors to be integration over facets of euclidean inner products of linear interpolation of values at vertices.

AREA FIXED

Activation of fixed total area constraint.

AREA NORMALIZATION

Convert force on vertex to mean curvature vector by dividing force by area of star around vertex.

AUTOCHOP

Do autochopping of long edges each iteration. Use AUTOCHOP := expr to set autochop length. Each iteration, any edge that is projected to become longer than the cutoff is bisected. If any bisections are done, the motion calculation is redone.

AUTOPOP

Toggles automatic deletion of short edges and popping of improper vertices each iteration. Before each iteration, any edge projected to shorten to under the critical length is deleted by identifying its endpoints. The critical length is calculated as L_{c} = \sqrt{2Δt}, where Δt is the time step or scale factor. Hence this should be used only with a fixed scale, not optimizing scale factor. The critical length is chosen so that instabilities do not arise in motion by mean curvature. If any edges are deleted, then vertices are examined for improper vertices as in the ’o’ command. Useful in string model. Not fully implemented in soapfilm model yet, since it doesn’t check facet areas.

AUTORECALC

Toggles automatic recalculation of the surface whenever adjustable parameters or energy quantity moduli are changed. Default is ON.

BOUNDARY CURVATURE

When doing integrals of mean curvature or squared curvature, the curvature of a boundary vertex cannot be defined by its neighbor vertices, so the area of the boundary vertex star instead is counted with an adjacent interior vertex.

CHECK INCREASE

Toggles checking for increase of energy in an iteration loop. If energy increases, then the iteration loop is halted. Meant for early detection of instabilities and other problems causing the surface to misbehave. Useful in doing a multiple iteration with a fixed scale. Caution: there are circumstances where an energy increase is appropriate, for example when there are volume or quantity constraints and conforming to the constraints means an energy increase initially.

CLIPPED

Sets torus quotient space display to clip to fundamental region. Not an on-off toggle. 3-way toggle with

RAW CELLS and CONNECTED.
COLORMAP

Use colormap from file. Use COLORMAP := "filename" to set file.

CONF EDGE SQCURV

Calculation of squared curvature by fitting sphere to edge and adjacent vertices (conformal curvature).

CONJ GRAD

Uses conjugate gradient method in fixed scale iteration.

CONNECTED

Sets quotient space display to do each body as a connected, wrapped surface. Not an on-off toggle. 3-way toggle with CLIPPED and RAW CELLS.

DEBUG

Print YACC debug trace of parsing of commands.

DETURCK

Motion by unit velocity along normal, instead of by curvature vector.

DIFFUSION

Activates diffusion step each iteration.

EFFECTIVE AREA

In area normalization, the resistance factor to motion is taken to be only the projection of the vertex star area perpendicular to the motion. If squared mean curvature is being calculated, this projected area is used in calculating the curvature.

ESTIMATE

Activates estimation of energy decrease in each iteration step. Useful only for a fixed scale factor much less than optimizing, so linear approximation is good.

    noindent FACET COLORS

Enables coloring of facets in certain graphics interfaces (xgraph, cheygraph). If off, facet color is white. Default on.

GRAVITY

Includes gravitational energy in total energy.

GV BINARY

Toggles sending data to geomview in binary format, which is faster than ascii. Default is binary. Ascii is useful for debugging.

HOMOTHETY

Adjust total volume of all bodies to fixed value after each iteration by uniformly scaling entire surface.

INTERP NORMALS

Display using interpolated vertex normals for shading for those graphics interfaces that support it.

KUSNER

Calculation of squared curvature by edge formula rather than vertex formula.

NORMAL CURVATURE

Calculation of squared curvature by taking area of vertex to be the component of the volume gradient parallel to the mean curvature vector.

NORMAL MOTION

Projects motion to surface normal, defined as the volume gradient. May be useful with squared curvature if vertices tend to slither sideways into ugly patterns.

OLD AREA

In the string model with area normalization, at a triple vertex Evolver normally tries to calculate the motion so that Von Neumann’s law will be obeyed, that is, the rate of area change is proportional to the number of sides of a cell. If old area is ON, then motion is calculated simply by dividing force by star area.

QUIET

Suppresses all normal output messages. Good while running scripts.

RAW CELLS

Sets quotient space display for plain, unwrapped facets. Not an on-off toggle. 3-way toggle with CLIPPED

and CONNECTED.
RUNGE KUTTA

Use Runge-Kutta method in iteration step (fixed scale factor only).

SHADING

Toggles facet shading in certain graphics interfaces (xgraph, cheygraph). Darkness of facet depends on angle of normal from vertical, simulating a light source above surface. Default is ON.

SHOW INNER

Display interior facets, those on 2 bodies.

SHOW OUTER

Display outer facets, those on 0 or 1 body.

SQGAUSS

Include squared gaussian curvature in energy. Soapfilm model only. Included in energy only, not force. This command included just as a convenient way to calculate the integral of the squared Gaussian curvature of a surface.

THICKEN

Display thickened surface. Use THICKEN := expr to set thickness.

§6.7. Graphics.

Graphics commands control the display of the surface, both the viewing transformation and which elements are seen. Keep in mind that there are two ways to display graphics on many systems: a native or internal mode entirely under the control of the Evolver (such as the X-windows display), and an external mode where information is provided to an external program such as geomview. Remember that hardcopy output (such as Postscript) will show what the native mode shows, not what geomview shows. Some general commands that affect both modes can be made from the main prompt Enter command:. Others, that control the native graphics, are made in a special graphics command mode entered by the ‘s’ or ‘show’ commands.

Note: native mode graphics are assumed to be slow, and don’t automatically redraw by default. In particular, there usually is no polling for window redraw events from the operating system, so after moving or resizing a graphics window you may have to give an explicit ’s’ command to redraw.

The display consists entirely of facets and edges. Special edges (fixed edges, boundary edges, constraint edges, triple edges) are always shown, unless you make their color CLEAR. The individual facet edges can be toggled with the graphics command ‘e’ or the geomview command ‘ae’.

This section lists the graphics mode commands first and then the general commands.

The graphics mode prompt is “Graphics command: ”. A graphics command is a string of letters followed by RETURN. Each letter causes an action. Some commands may be preceded by an integer count of how many repetitions to do. Example command: rr15u2z Rotation commands may be preceded by a real number giving the degrees of rotation; an integer will give a repetition count with the default angle of 6 degrees. If you mean to specify the angle, be sure to include a decimal point.

Commands in graphics command mode:

Repeatable commands:

u Tip up. Rotates image six degrees about horizontal axis.

d Tip down. Rotates image six degrees other way.

r Rotate right. Rotates by six degrees about vertical axis.

l Rotate left. Rotates by six degrees the other way.

c Rotate clockwise by six degrees about center of screen.

C Rotate counterclockwise by six degrees about center of screen.

z Zoom. Expands image by factor of 1.2.

s Shrink. Contracts image by factor of 1.2.

arrow keys Move image half screen width in appropriate direction. May not work on all terminals.

Non-repeatable commands:

R Reset viewing angles to original defaults and rescale the image to fit the viewing window.

e Toggle showing all the facet edges.

h Toggle hiding hidden surfaces. When ON, takes longer to display images, but looks better.

b Toggles display of bounding box. Useful for visualizing orientation.

t Reset mode of displaying torus surfaces. Choice of raw cell, connected bodies, or clipped unit cell. See the Torus section for more explanation. w Toggles display of facets entirely on constraints. Meant to show traces where surfaces meet constraints. “w” stands for “wall”. Applies to facets run up against inequality constraints.

B Toggles display of facets on boundaries or equality constraints attribute.

v Toggles showing of convex and concave edges in different colors. “v” stands for “valleys”.

+ Increments color number used for facet edges.

- Decrements color number used for facet edges.

? Prints help screen for graphics commands.

q,x Exit from graphics mode, and return to main command mode.

Graphics commands in main command mode.

The following commands are entered at the main command prompt:

SHOW elements [ name ] [ WHERE expr ]

Shows surface with screen graphics and goes into graphics command mode. The condition will prescibe which facets (or edges in the string model) will be displayed with graphics commands. The prescription will remain in effect until the next show query. Example:

    show facets where tag <= 3
SHOWQ

Will show surface without going into graphics command mode, with whatever condition is in effect.

SHOW EXPR elements [ name ] [ WHERE expr ]

This sets the display prescription without actually doing any graphics. Useful when graphics are not available, as when working remotely, and you want to generate graphical file output with a prescription.

SHOW TRANS string

Applies string of graphics commands to the image transformation matrix without doing any graphic display. The string must be in double quotes or be a string variable, and is the same format as is accepted by the regular graphics command prompt. Example:

    show trans "rrdd5z"
TRANSFORMS [ON|OFF]

Toggles displaying of multiple viewing transforms. These are either those listed explicitly in the datafile or those generated with a previous transform expr command.

TRANSFORM EXPR string

If view transform generators were included in the datafile, then a set of view transforms may be generated by a expression with a syntax much like a regular expression. An expression generates a set of transforms, and are compounded by the following rules. Here a lower-case letter stands for one of the generators, and an upper-case letter for an expression.

a  Generates set I, a.
AB  Generates all ordered products of pairs.
nA  Generates all n-fold ordered products.
A|B  Generates union of sets A and B.
(A)  Grouping; generates same set as A.

The precedence order is that nA is higher than AB which is higher that A|B. Note that the expression string must be enclosed in double quotes or be a string variable. Examples:

transform expr "3(a|b|c)" (All products of 3 or fewer generators.) transform expr "abcd" (Generates sixteen transforms) All duplicate transforms are removed, so the growth of the sets do not get out of hand. Note that the

identity transform is always included. The letter denoting a single generator may be upper or lower case. The order is the same order as in the datafile. If 26 generators are not enough for somebody, let me know.

VIEW 4D [ON|OFF]

For geomview graphics, sends a full 4D file in 4OFF and 4VECT format. Default is OFF, so geomview gets just the 3D projection.

§6.8. Script examples.

Here are some longer examples using the command language. Since the language expects one command per line, long commands have line breaks only within braces so Evolver can tell the command continues. I tend to break long commands up into several little named commands. The best way to use these longer commands is to put their definitions at the end of your datafile in a ”read” section, or have them in a separate file you can read with the ”read” command.

1. A script to evolve and regularly save dumps in different files:

nn := 1
while nn < 100 do { g 100; ff := sprintf "stage%g.dmp",nn; dump ff; nn:=nn+1}

2. Your own iterate command to print out what you want, like energy change each step:

olde := total_energy;
gg := { g; printf "Change:  %10.5g\n",total_energy-olde; olde:=total_energy}

3. A command to print a file listing vertices and edges (happens to be the OFF format for geomview):

aa:=foreach vertex vv do { printf "%f %f %f \n",x,y,z }
bb:=foreach facet ff do { printf "3 ";
foreach ff.vertex do printf "%g ",id-1; printf "\n"}
do_off:= { printf "OFF\n%g %g %g\n",vertex_count,facet_count,edge_count;aa;bb}

This prints to stdout, so you would use it with redirection:

Enter command:  do_off | "cat >test.off"

4. A command to duplicate the fundamental region of a symmetric surface across a plane of mirror symmetry.

// Producing datafile surface duplicated across a mirror.
// This is set up to do constraint 1 only.
// Mirrored constraints given x numbers.
// Resulting datafile needs mirrored constraints added by hand.
// Usage:  mm | "cat >filename.fe"
// Mirror plane aa*x + bb*y + cc*z = dd
// You have to set these to the proper coefficients for constraint 1

// before doing mm.
aa := 1; bb := 0; cc := 0; dd := 0
ma := foreach vertex do { printf "%g  %g %g %g ",id,x,y,z;
if fixed then printf "fixed ";
if on_constraint 1 then printf "constraint 1 " ;
printf "\n";

}

mb := { voff := vertex_count }
mc := foreach vertex do if not on_constraint 1 then {
lam := 2*(dd-aa*x-bb*y-cc*z)/(aa*aa+bb*bb+cc*cc);
printf "%g  %g %g %g ",id+voff,x+lam*aa,y+lam*bb,z+lam*cc;
if fixed then printf "fixed ";
if on_constraint 1 then printf "constraint 1x " ;
printf "\n";

}

md := foreach edge ee do {
printf "%g   ",id;
foreach ee.vertex do printf "%g ",oid;
if ee.fixed then printf "fixed ";
if ee.on_constraint 1 then printf "constraint 1 " ;
printf "\n";

}

me := { eoff := edge_count }
mf := foreach edge ee do if not on_constraint 1 then {
printf "%g   ",id+eoff;
foreach ee.vertex do {
if on_constraint 1 then { printf "%g ",oid;}
else { printf "%g ", id+voff;} } ;
if ee.fixed then printf "fixed ";
if ee.on_constraint 1 then printf "constraint 1x " ;
printf "\n";

}

mg := foreach facet ff do {
printf "%g   ",id;
foreach ff.edge do printf "%g ",oid;
if ff.fixed then printf "fixed ";
if ff.on_constraint 1 then printf "constraint 1 " ;
printf "\n";

}

mh := { foff := facet_count}
mi := foreach facet ff do if not on_constraint 1 then {
printf "%g   ",id+foff;
foreach ff.edge do
{ if on_constraint 1 then { printf "%g ",id;}
else printf "%g ",(oid<0?-(id+eoff):id+eoff); } ;
if ff.fixed then printf "fixed ";
if ff.on_constraint 1 then printf "constraint 1x " ;
printf "\n";

}

mm := { list topinfo;
printf "\nVertices\n"; ma;mb;mc;
printf "\nEdges\n"; md;me;mf;

printf "\nFaces\n"; mg;mh;mi;

}

§6.9. Interrupts.

Evolver operation may be interrupted with the standard keyboard interrupt, CTRL-C usually (SIGINT for you unix gurus). During repeated iteration steps, this will set a flag which causes the repetition to cease after the current step. Otherwise, the current command is aborted and control returns to the main prompt. If Evolver receives SIGTERM (say from the unix kill command), it will dump to the default dump file and exit. This is useful for stopping background processes running on a script and saving the current state. The same thing will happen with SIGHUP, so losing a modem connection will save the current surface.

§6.10. Graphics output file formats.

This section explains the file formats produced by the P command.

§6.10.1. Pixar and Minneview.

A file extension “.quad” is automatically added to the filename. This format lists the facets as quadrilaterals, with the third vertex being repeated. Two options are

available:

1. Interpolated normals. This requires recording the normal vector at each vertex.

2. Colors. This requires recording color values (red,green,blue,alpha) at each vertex. The file starts with a keyword:

CNPOLY colors and normals
NPOLY normals
CPOLY colors
POLY  vertices only.

There follows one facet per line. Each vertex is printed as three coordinates, three normal vector components (if used), and four color values (if used).

The coordinates are the true coordinates, not transformed for viewing. Facets are not ordered back to front, since Pixar and Minneview handle viewing and hidden surfaces themselves.

§6.10.2. Geomview.

The P command menu option 8 starts the geomview program and feeds it data through a named pipe. The file fed into the named pipe follows the geomview command file format. Parts of it may be in binary format for speed. A copy of the file may be created by using the P menu option A to create a pipe, and then using cat on the pipe to read the pipe into your own file. End the file with P option B.

Some handy geomview commands (to be entered with mouse in the geomview image window): r Set mouse to rotate mode. t Set mouse to translate mode. z Set mouse to zoom mode. ab Toggle drawing bounding box. ae Toggle facet edge drawing. af Toggle drawing facets.

Remember that what you see in geomview is not what you will get with a Postscript output, since Evolver knows nothing about the commands you give geomview.

§6.10.3. PostScript.

The output file is suitable for direct input to a PostScript device. The view is from the current viewing angle as established with the s command or the show trans "string" command. Edges and facets are drawn back to front. You are given a choice whether to draw all the grid lines (all the facet edges). Even if you don’t choose to, all special edges (fixed, boundary, triple edges, bare edges) will still be drawn. You can prevent an edge from being drawn by giving it the color CLEAR. Color data is optionally included; you will be prompted for your choice.

§6.10.4. Triangle file.

This is a format I put in to get the data I wanted for a certain external application. This lists one facet per line. Coordinates are viewing transformed to the screen plane. Facets are listed back to front. The format for a line is

x1 y1 x2 y2 x3 y3 w1 w2 w3 d

Here (x1,y1), (x2,y2), and (x3,y3) are the vertex coordinates, and w1, w2, and w3 are types of edges:

type 0: ordinary edge, no constraints or boundaries, two adjacent facets. type 1: edge on constraint or boundary, but not fixed. type 3: fixed edge. type 4: edge adjacent to at least 3 facets.

These values can be used to weight different types of edges to illustrate the structure of the surface. d is the cosine of the angle of the facet normal to the screen normal. Useful for shading facets.

§6.10.5. SoftImage file.

This consists of two files, one with extension .mdl and one with .def. The extensions are added automatically.

The .mdl file has fixed content, except for the name of the picture:

SOFTIMAGE 4D Creative Environment   v 1.6  "ASCII"
MODL   "picname"
{
type    PMSH
nbdef   1
scal      1.000000  1.000000  1.000000
rot       0.000000  0.000000  0.000000
trans     0.000000  0.000000  0.000000
}

The .def file contains a vertex list followed by a quadrilateral list in terms of the vertices.

SOFTIMAGE 4D Creative Environment   v 1.6  "ASCII"
PMSH   "picfile"
{
vertex 14
1   0.000000   0.000000   0.000000
2   1.000000   0.000000   0.000000
3   1.000000   1.000000   0.000000
4   0.000000   1.000000   0.000000
5   0.000000   0.000000   1.000000

6   1.000000   0.000000   1.000000
7   1.000000   1.000000   1.000000
8   0.000000   1.000000   1.000000
9   0.500000   0.000000   0.500000
10   1.000000   0.500000   0.500000
11   0.500000   1.000000   0.500000
12   0.000000   0.500000   0.500000
13   0.500000   0.500000   1.000000
14   0.500000   0.500000   0.000000
polygon 24
5          1          9          9
6          2         10         10
7          3         11         11
8          4         12         12
8          5         13         13
2          1         14         14
9          1          2          2
9          2          6          6
9          6          5          5
10          2          3          3
10          3          7          7
10          7          6          6
11          3          4          4
11          4          8          8
11          8          7          7
12          4          1          1
12          1          5          5
12          5          8          8
13          5          6          6
13          6          7          7
13          7          8          8
14          1          4          4
14          4          3          3
14          3          2          2
}

Chapter 7. Technical Reference.

This chapter explains the mathematics and algorithms used by the Evolver. You don’t need to know this stuff to use the Evolver, but it’s here for anybody who want to know what exactly is going on and for a programmer’s reference.

This chapter currently describes only the linear soapfilm model.

§7.1. Notation.

The term “original orientation” refers to the orientation of an element as recorded in its structure. Original orientations begin with the orientations in the datafile, and are inherited thereafter during refinement and other operations.

Edges.

In referring to a edge e:

t⃗ and h⃗ are the tail and head of the side.

s⃗ is the edge vector, s⃗ = h⃗ − t⃗.

a_{i} and w_{i} are the abscissas and weights for gaussian quadrature on the interval [0,1].

Facets.

In referring to a facet f: T is the facet surface tension. v_{0}, v_{1}, and v_{2} are the vertices around a facet in counterclockwise order. s⃗_{0}, s⃗_{1}, and s⃗_{2} are the edges in counterclockwise order around a facet, and v_{0} is the tail of s⃗_{0}. a_{ij} and w_{i} are the barycentric coordinates (j is the vertex index) and weights for evaluation point i for

gaussian quadrature on a triangle.

§7.2. Surface representation.

The basic geometric elements are vertices, edges, three-sided facets, and bodies. There is a structure type for each (see skeleton.h for definitions). In addition, there is a facet-edge structure for each incident facet and edge pair. Edges, facets, and facet-edges are oriented. Vertices and bodies could in principle be oriented, but they are not. Elements can be referred to by identifiers which implicitly contain an orientation for the element. Currently these identifiers are 32-bit bitfields with one bit for an orientation bit and the rest for a list pointer. Thus one structure is used for both orientations, and the programmer doesn’t have to worry about the orientation of the stored structure; all the orientation conversions are taken care of by the structure access functions and macros.

The topology of the surface is defined by the cross-references among the elements. Each facet-edge is in two doubly-linked lists: one of facet-edges around the edge, and another of facet-edges around the facet. Each edge points to its two endpoints (called head and tail) and to one facet-edge in the facet-edge loop around it. Each facet points to one facet-edge in the facet-edge loop around it, and to the two bodies on its two sides. Each body points to one facet-edge on its boundary. The upshot of all this is that it is easy to find all the facets around an edge or all the edges around a facet, but there is no quick way to find all the elements that contain a given vertex or are on the boundary of a given body. These cases can be handled only by exhaustive search. Luckily, this can be avoided most of the time. This is the major drawback of the facet-edge representation, and something will be done about it if the Evolver gets a major re-write.

§7.3. Energies and forces.

As described in the Energies section of the Model chapter, all energies are carried by facets or edges. (Volumes for ideal gas energy are calculated from facets and edges also.) Hence it is straightforward to run through all the edges and facets, adding up their energies and accumulating their energy gradients onto the total forces on their vertices. This section gives the formulas used for all the types of energies and forces.

§7.3.1. Surface tension.

Energy of facet.

E = \frac{T}{2} ∥s⃗_{0} × s⃗_{1}∥.  (1)

Force on vertex.

F⃗ (v_{0}) = \frac{T}{2} \frac{s⃗_{1} × (s⃗_{0} × s⃗_{1})}{∥s⃗_{0} × s⃗_{1}∥}.  (2)

§7.3.2. Crystalline integrand.

W⃗ is the Wulff vector corresponding to the normal of the facet. The facet has its original orientation.

Energy of facet.

E = \frac{T}{2} W⃗ · s⃗_{0} × s⃗_{1}.  (3)

Force on vertex.

F⃗ (v_{0}) = \frac{T}{2} s⃗_{1} × W⃗ .  (4)

§7.3.3. Gravity.

Let G be the acceleration of gravity. Let ρ_{+} be the density of the body from which the facet normal is outward, and let ρ_{−} be the density of the body with inward normal. Let z_{i} be the z coordinate of v_{i}. The following formulas are exact for flat facets.

Energy of facet.

\begin{aligned}E & = G(ρ_{+} − ρ_{−}) ∫ ∫_{facet} \frac{1}{2}z^{2}k⃗ · d⃗A \\ & = G(ρ_{+} − ρ_{−})\frac{1}{12}(z_{0}^{2} + z_{1}^{2} + z_{2}^{2} + z_{0}z_{1} + z_{0}z_{2} + z_{1}z_{2})\frac{1}{2}(k⃗ · s⃗_{0} × s⃗_{1}). \end{aligned}  (5)

Force on vertex.

F⃗ (v_{0}) = −\frac{G(ρ_{+} − ρ_{−})}{24}((2z_{0} + z_{1} + z_{2})(k⃗ · s⃗_{0} × s⃗_{1})k⃗ + (z_{0}^{2} + z_{1}^{2} + z_{2}^{2} + z_{0}z_{1} + z_{0}z_{2} + z_{1}z_{2})(k⃗ × s⃗_{1})). (6)

§7.3.4. Constraint integrals.

This occurs for each constraint with energy integrand that an edge is deemed to satisfy. The edge orientation is the original one. The points evaluated for gaussian quadrature are

x⃗_{i} = a_{i}h⃗ + (1 − a_{i})t⃗.  (7)

Energy of edge energy.

E = ∫_{edge} E⃗ · d⃗s ≈ ∑_{i} w_{i}E⃗(x⃗_{i}) · s⃗.  (8)

Force on vertex.

Component notation and the Einstein summation convention are used here, as vector notation gets confusing. Comma subscript denotes partial derivative.

\begin{aligned}F_{j}(head) & = −w_{i}(a_{i}E_{k,j}(x⃗_{i})s_{k} + E_{j}(x⃗_{i})), \\ F_{j}(tail) & = −w_{i}((1 − a_{i})E_{k,j}(x⃗_{i})s_{k} − E_{j}(x⃗_{i})). \end{aligned}  (9)

§7.3.5. Surface integrands.

There is a contribution from each surface energy that a facet is deemed to have. The facet orientation is the original one. The points evaluated for gaussian quadrature are

x⃗_{i} = ∑_{j} a_{ij}v⃗_{j}.  (10)

Energy of facet.

E = ∫ ∫_{facet} E⃗ · d⃗A ≈ \frac{1}{2} ∑_{i} w_{i}E⃗(x⃗_{i}) · s⃗_{0} × s⃗_{1}.  (11)

Force on vertex.

Component notation and the Einstein summation convention are used here, as vector notation gets confusing. Comma subscript denotes partial derivative.

F_{j}(v_{m}) = −\frac{1}{2}w_{i}(a_{im}E_{k,j}(x⃗_{i})(s⃗_{0} × s⃗_{1})_{k}) + (E(x⃗_{i}) × s⃗_{1})_{j}).  (12)

§7.3.6. Gap areas.

Every edge on a CONVEX constraint contributes here. Let q⃗_{t} be the projection of s⃗ on the constraint tangent space at the tail, and q⃗_{h} likewise at the head. k is the gap constant.

Energy of edge.

E = \frac{k}{12}(\sqrt{(s^{2} − q_{t}^{2})q_{t}^{2}} + \sqrt{(s^{2} − q_{h}^{2})q_{h}^{2}}).  (13)

Force on vertex.

This energy is not the direct derivative of the previous formula. Rather it comes directly from the geometric decrease in gap area due to vertex motion. The formula is derived for the case of a flat gap surface, which should be close enough to the truth to be useful.

\begin{aligned}F⃗ (tail) & = \frac{k\sqrt{(s^{2} − q_{t}^{2})q_{t}^{2}}}{2q_{t}^{2}} q⃗_{t}, \\ F⃗ (head) & = −\frac{k\sqrt{(s^{2} − q_{h}^{2})q_{h}^{2}}}{2q_{h}^{2}} q⃗_{h}. \end{aligned}  (14)

§7.3.7. Ideal gas compressibility.

The ideal gas mode is in effect when the keyword is listed in the first part of the datafile. Then prescribed volumes become in essence moles of gas. P_{amb} is the ambient pressure. A body has prescribed volume V_{0}, actual volume V, and pressure P. The temperature is assumed constant, so P V = P_{amb}V_{0}.

Energy of body.

\begin{aligned}E & = ∫ P − P_{amb}dV = ∫ P_{amb}V_{0}/V − P_{amb}dV \\ & = P_{amb}V_{0} ln(V/V_{0}) − P_{amb}(V − V_{0}) \end{aligned}  (15)

Note the energy starts at 0 for V = V_{0}.

Force on vertex.

Let g⃗_{m} be the gradient of the volume body m as a function of the coordinates of the vertex. Let P_{m} be the pressure of body m. Then

F⃗ (v) = − ∑_{m}(P_{m} − P_{amb})g⃗_{m}.  (16)

§7.3.8. Prescribed pressure.

Ambient pressure is taken to be zero, and body pressures are constant.

Energy of body.

E = ∫ P dV = P V.  (17)

Force on vertex.

Let g⃗_{m} be the gradient of the volume body m as a function of the coordinates of the vertex. Let P_{m} be the pressure of body m. Then

F⃗ (v) = − ∑_{m} P_{m}g⃗_{m}.  (18)

§7.3.9. Squared mean curvature.

The integral of squared mean curvature in the soapfilm model is calculated as follows: Each vertex v has a star of facets around it of area A_{v}. The force due to surface tension on the vertex is

F_{v} = −\frac{∂A_{v}}{∂v} .  (19)

Since each facet has 3 vertices, the area associated with v is A_{v}/3. Hence the average mean curvature at v is

h_{v} = \frac{1}{2} \frac{F_{v}}{A_{v}/3},  (20)

and this vertex’s contribution to the total integral is

E_{v} = h_{v}^{2}A_{v}/3 = \frac{1}{4} \frac{F_{v}^{2}}{A_{v}/3}.  (21)

The corresponding calculation is done for the string model. E_{v} can be written as an exact function of the vertex cooordinates, so the gradient of E_{v} can be fed into the total force calculation.

Philosophical note: The squared mean curvature on a triangulated surface is technically infinite, so some kind of approximation scheme is needed. The alternative to locating curvature at vertices is to locate it on the edges, where it really is, and average it over the neighboring facets. But this has the problem that a least area triangulated surface would have nonzero squared curvature, whereas in the vertex formulation it would have zero squared curvature.

Practical note: The above definition of squared mean curvature seems in practice to be subject to instablities. One is that sharp corners grow sharper rather than smoothing out. Another is that some facets want to get very large at the expense of their neighbors. Hence a couple of alternate definitions have been added.

Effective area curvature: The area around a vertex is taken to be the magnitude of the gradient of the volume. This is less than the true area, so makes a larger curvature. This also eliminates the spike instability, since a spike has more area gradient but the same volume gradient. Letting N_{v} be the volume gradient at v,

h_{v} = \frac{1}{2} \frac{F_{v}}{||N_{v}||/3},  (22)

and

E_{v} = h_{v}^{2}A_{v}/3 = \frac{3}{4} \frac{F_{v}^{2}}{||N_{v}||^{2}} A_{v},  (23)

The facets of the surface must be consistently oriented for this to work, since the evolver needs an “inside” and “outside” of the surface to calculate the volume gradient. This mode is toggled by the “effective area ON | OFF” command. There are still possible instabilities where some facets grow at the expense of others.

Normal curvature: To alleviate the instability of effective area curvature, the normal curvature mode considers the area around the vertex to be the component of the volume gradient parallel to the mean curvature vector, rather than the magnitude of the volume gradient. Thus

h_{v} = \frac{1}{2} \frac{F_{v}||F_{v}||}{N_{v} · F_{v}/3},  (24)

E_{v} = \frac{3}{4} [\frac{F_{v} · F_{v}}{N_{v} · F_{v}} ]^{2} A_{v}.  (25)

This is still not perfect, but is a lot better. This mode is toggled by the “normal curvature ON | OFF” command. If you have effective area on, then normal curvature will supersede it, but effective area will be in effect when normal curvature is toggled off.

Curvature at boundary: If the edge of the surface is a free boundary on a constraint, then the above calculation gives the proper curvature under the assumption the surface is continued by reflection across the constraint. This permits symmetric surfaces to be represented by one fundamental region. If the edge of the surface is a fixed edge or on a 1-dimensional boundary, then there is no way to calculate the curvature on a boundary vertex from knowledge of neighboring facets. For example, the rings of facets around the bases of a catenoid and a spherical cap may be identical. Therefore curvature is calculated only at interior vertices, and when the surface integral is done, area along the boundary is assigned to the nearest interior vertex.

§7.3.10. Squared Gaussian curvature.

The integral of squared Gaussian curvature over a soapfilm model surface only applies where the star of facets around each vertex is planar. It is calculated as follows. The total gaussian curvature at a vertex is the angle deficit around the vertex, 2π − ∑ θ_{i}, where θ_{i} are the vertex angles of the facets adjacent to the vertex. The average Gaussian curvature is the total Gaussian curvature divided by one third of the area of the star. The average is then squared and integrated over one third of the star.

At boundary: Treated the same way squared mean curvature is.

§7.4. Named quantities and methods.

This section gives details on the calculations of the methods used in named quantities.

§7.4.1. Edge length.

Method name edge tension or edge length. This method is entirely separate from the default energy calculation for strings. For each edge, the value is the length in Euclidean coordinates. The edge density attribute, used in the default energy calculation, is not used. No metric is used.

§7.4.2. Facet area.

Method name facet tension or facet area. This method is entirely separate from the default energy calculation for surfaces. For each facet, the value is the area in Euclidean coordinates. The area density attribute, used in the default energy calculation, is not used. No metric is used.

§7.4.3. Path integrals.

Method name edge scalar integrand. A scalar integrand f(x⃗) is evaluated on a edge by gaussian quadrature:

L ∑_{i} w_{i}f(x⃗_{i})  (1)

where L is the Euclidean length of the edge, w_{i} are the gaussian weights, and x_{i} are the evaluation points.

§7.4.4. Scalar surface integral.

Method name facet scalar integral. A scalar integrand f(x⃗) is evaluated on a facet by gaussian quadrature:

A ∑_{i} w_{i}f(x⃗_{i})  (2)

where A is the Euclidean area of the facet, w_{i} are the gaussian weights, and x_{i} are the evaluation points.

§7.4.5. Hooke energy.

Method name Hooke energy. One would often like to require edges to have fixed length. The total length of some set of edges may be constrained by defining a fixed quantity. This is used to fix the total length of an evolving knot, for example. But to have one constraint for each edge would be impractical, since projecting to n constraints requires inverting an n × n matrix. Instead there is a Hooke’s Law energy available to encourage edges to have equal length. Its form per edge is

E = (L − L_{0})^{p}  (3)

where L is the edge length, L_{0} is the equilibrium length, embodied as an adjustable parameter ‘hooke length’, and the power p is an adjustable parameter ‘hooke power’. The default power is p = 2, and the default equilibrium length is the average edge length in the initial datafile. You will want to adjust this, especially if you have a total length constaint. A high modulus will decrease the hooke component of the total energy, since the restoring force is linear in displacement and the energy is quadratic (when p = 2).

§7.4.6. Local Hooke energy.

Method name local hooke energy. Energy of edges as springs with equilibrium length being average of lengths of neighbor edges. Actually, the energy is calculated per vertex,

E = (\frac{L_{1} − L_{2}}{L_{1} + L_{2}} )^{2}  (4)

where L_{1} and L_{2} are the lengths of the edges adjacent to the vertex.

§7.4.7. Integral of mean curvature.

Method name mean curvature integral. This computes the integral of the signed scalar mean curvature over a surface. The computation is exact, in the sense that for a polyhedral surface the mean curvature is concentrated on edges and singular there, but the total mean curvature for an edge is the edge length times its dihedral angle. The contribution of one edge is

L arccos (\frac{a⃗ · vecba⃗ · vecb − a⃗ · a⃗b⃗ · c⃗}{(a⃗ · a⃗b⃗ · b⃗ − (a⃗ · b⃗)^{2})^{1/2}(a⃗ · a⃗c⃗ · c⃗ − (a⃗ · c⃗)^{2})^{1/2}} ) ,  (5)

where L is the edge length, a⃗ is the edge, and b⃗, c⃗ are the two adjacent sides from the tail of a⃗.

§7.4.8. Integral of gaussian curvature.

Method name gauss curvature integral. This computes the total gaussian curvature of a surface with boundary. The gaussian curvature of a polyhedral surface may be defined at an interior vertex as the angle deficit of the adjacent angles. But as is well-known, total gaussian curvature can be computed simply in terms of the boundary vertices, which is what is done here. The total gaussian curvature is implemented as the total geodesic curvature around the boundary of the surface. The contribution of a boundary vertex is

(∑_{i} θ_{i}) − π.  (6)

The total over all boundary vertices is exactly equal to the total angle deficit of all interior vertices plus 2πχ, where χ is the Euler characteristic of the surface.

§7.4.9. Gravity.

    Method name gravity method

This does the same calculation as the usual gravity. Let G be the acceleration of gravity. Let ρ_{+} be the density of the body from which the facet normal is outward, and let ρ_{−} be the density of the body with inward normal. Let z_{i} be the z coordinate of v_{i}. The following formulas are exact for flat facets.

\begin{aligned}E & = G(ρ_{+} − ρ_{−}) ∫ ∫_{facet} \frac{1}{2}z^{2}k⃗ · d⃗A \\ & = G(ρ_{+} − ρ_{−})\frac{1}{12}(z_{0}^{2} + z_{1}^{2} + z_{2}^{2} + z_{0}z_{1} + z_{0}z_{2} + z_{1}z_{2})\frac{1}{2}(k⃗ · s⃗_{0} × s⃗_{1}). \end{aligned}  (7)

§7.4.10. Average crossing number.

Method name average crossings. Between pairs of edges, energy is inverse cube power of distance between midpoints of edges, times triple product of edge vectors and distance vector:

E = 1/d^{3} ∗ (e1, e2, d).  (8)

§7.4.11. Knot energies.

One way of smoothing a knotted curve is to put “electric charge” on it and let it seek its minimum energy position. To prevent the curve from crossing itself and unknotting, the potential energy should have a high barrier (preferably infinite) to curve crossing. This is often done by using a inverse power law potential with a higher power than the standard inverse first power law of electrostatics. A length constraint on the curve is generally necessary to prevent the charges from repelling each other off to infinity.

The Evolver implements several discrete approximations of these potentials. All use the quantity-method feature described in the next section. That is, each type of energy is a method, and these methods can be attached to user-named quantities. The power in the power law is an adjustable parameter ‘knot power’. It may be changed interactively with the ‘A’ command. The modulus is the multiple of the method added to the quantity. A modulus of 0 inactivates the method. The value of the quantity can be seen with the ’A’ command under the name of the quantity.

One of the general problems with discretization is that edges don’t know when they are crossing each other. Edges can cross without either their endpoints or midpoints getting close, especially if said edges get long. To keep edge lengths close to equal, there is also an energy called ‘hooke energy’ described in the Model chapter. Another energy of interest is total squared curvature, which is elastic bending energy. However, it does not provide any barrier to crossing.

All the knot energies proper correspond to double integrals in the continuous limit, so are double sums. Hence the computation time is quadratic in the number of vertices. Hooke energy and square curvature are linear time.

Generally, the continuous integrals are divergent, and in the literature have various regularization terms subtracted off. However, the Evolver discretizations do not have such regularization. The logic is that the discrete versions are finite, and the discete regularization terms are practically constant, so there is no sense in wasting time computing them.

A sample datafile incorporating some of these is knotty.fe.

These discrete knot energies are classified below by whether they are computed over all pairs of vertices, edges, or facets.

I. Vertex pair energies. The total energy is the sum over all pairs of vertices v_{1}, v_{2}:

E = \frac{1}{2} ∑_{v_{1}≠v_{2}} E_{v_{1}v_{2}}.

Ia. Conducting wire. There is a unit charge on each vertex, and edge lengths are not fixed.

E_{v_{1}v_{2}} = \frac{1}{|v_{1} − v_{2}|^{p}}

This is approximating the continuous integral

E = ∫ ∫ (\frac{1}{|x(s) − x(t)|^{p}} − \frac{1}{d(s, t)^{p}} ) ρ(s)ρ(t) ds dt

where s, t are arclength parameters, d(s, t) is the shortest arc distance between points, and ρ is charge density. The second term in this integral is a normalization term to prevent infinite energy. It is not present in the discrete version, since the discrete energy is finite and the normalization term is approximately constant. Note that refining doubles the total charge, so refining should be accompanied by reducing the modulus by a factor of 4.

Note that the discrete energy can apply to any dimension surface in any dimension ambient space since it does not use the regularization. It can be used as a dust energy for sphere packing (high knot power for a hard core potential) or an energy for knotted 2-spheres in 4-space.

The default power is 2.

Datafile line:

    quantity knot energy ENERGY modulus 1 global method knot energy

Ib. Insulating wire.

E_{v_{1}v_{2}} = \frac{L_{1}L_{2}}{|v_{1} − v_{2}|^{p}}

where L_{1} is the average length of the two edges adjacent to v_{1} and L_{2} is the average length of the two edges adjacent to v_{2}. This is approximating the continuous integral

E = ∫ ∫ (\frac{1}{|x(s) − x(t)|^{p}} − \frac{1}{d(s, t)^{p}} ) ds dt

where s, t are arclength parameters, d(s, t) is the shortest arc distance between points. This corresponds to a uniform charge density, which is why it is an ‘insulating wire’. The second term in the integral is a normalization term, which is not included in the discrete version since it is approximately constant. However, a discrete normalization term can be calculated as a separate quantity.

This energy assumes that each vertex has exactly 2 edges attached to it, so it is not suitable for surfaces. Datafile line:

    quantity knot energy ENERGY modulus 1 global method edge knot energy

The ‘edge’ in the method name refers to the edge-weighting of the vertex charges. The default power is 2.

Datafile line for normalization term:

    quantity norm term INFO ONLY modulus 1 global method uniform knot energy normalizer

Ic. Insulating surface.

E_{v_{1}v_{2}} = \frac{A_{1}A_{2}}{|v_{1} − v_{2}|^{p}}

where A_{1} is the area of the facets adjacent to v_{1} and A_{2} is the area of the facets adjacent to v_{2}. This corresponds to a uniform charge density, which is why it is an ‘insulating surface’.

Datafile line:

    quantity knot energy ENERGY modulus 1 global method facet knot energy

The ‘facet’ in the method name refers to the edge-weighting of the vertex charges. The default power is 4.

II. Edge pair energies. The total energy is the sum over all pairs of edges e_{1}, e_{2} with endpoints v_{11}, v_{12}, v_{21}, v_{22}:

E = \frac{1}{2} ∑_{e_{1}≠e_{2}} E_{e_{1}e_{2}}

IIa. Box energy. This energy is due to Gregory Buck. It has the form

E_{e_{1}e_{2}} = \frac{L_{1}L_{2}}{(d_{1} + d_{2} + d_{3} + d_{4} − 2(L_{1} + L_{2}))^{p}} .

Here L_{1}, L_{2} are the edge lengths and d_{1}, d_{2}, d_{3}, d_{4} are the distances between endpoints on different edges. This provides an infinite barrier to crossing in the discrete case in a simple form involving only the edge endpoints. The denominator becomes zero for parallel coincident lines or for perpendicular lines on opposite edges of a tetrahedron. This energy should not be turned on until the curve is refined enough that the denominator is always positive. The default power is 2.

Datafile line:

    quantity buck energy ENERGY modulus 1 global method buck knot energy

IIb. Normal projection energy. This energy is also due to Gregory Buck. It tries to eliminate the need for a normalization term by projecting the energy to the normal to the curve. Its form is

E_{e_{1}e_{2}} = \frac{L_{1}L_{2} cos^{p} θ}{|x_{1} − x_{2}|^{p}}

where x_{1}, x_{2} are the midpoints of the edges and θ is the angle between the normal plane of edge e_{1} and the vector x_{1} − x_{2}. The default power is 2.

Datafile line:

    quantity proj energy ENERGY modulus 1 global method proj knot energy

IIc. Conformal circle energy. This energy is due to Peter Doyle, who says it is equivalent in the continuous case to the insulating wire with power 2. Its form is

E_{e_{1}e_{2}} = \frac{L_{1}L_{2}(1 − cos α)}{|x_{1} − x_{2}|^{2}},

where x_{1}, x_{2} are the midpoints of the edges and α is the angle between edge 1 and the circle through x_{1} tangent to edge 2 at x_{2}. Only power 2 is implemented.

Datafile line:

quantity circle energy ENERGY modulus 1 global method circle knot energY

III. Facet pair energies. These energies are sums over all pairs of facets.

IIIa. Conformal sphere energy. This is the 2D surface version of the conformal circle energy. Its form is

E_{f_{1}f_{2}} = \frac{A_{1}A_{2}(1 − cos α)}{|x_{1} − x_{2}|^{4}},

where A_{1}, A_{2} are the facet areas, x_{1}, x_{2} are the barycenters of the facets, and α is the angle between f_{1} and the sphere through x_{1} tangent to f2 at x_{2}.

There is a nice expression for the general dimensional version of this energy:

4E_{f_{1}f_{2}} = \frac{A_{1}A_{2}}{r^{4}} − \frac{1}{r^{6}} det [ \begin{matrix}r · r & r · s_{1} & r · s_{2} \\ t_{1} · r & t_{1} · s_{1} & t_{1} · s_{2} \\ t_{2} · r & t_{2} · s_{1} & t_{2} · s_{2} \end{matrix} ]

where r = |x_{1} − x_{2}| and s_{1}, s_{2}, t_{1}, t_{2} are the sides of the two facets.

Datafile line:

quantity sphere energy ENERGY modulus 1 global method sphere knot energy

§7.5. Volumes.

A body can get volume from a surface integral over its boundary facets and from content integrals of boundary facets edges on constraints. Also important are the gradients of body volumes at vertices. There are two choices built-in surface integrals: default and SYMMETRIC_CONTENT.

§7.5.1. Default facet integral.

A facet oriented with normal outward from a body makes a contribution to its volume of

V = ∫ ∫_{facet} zk⃗ · d⃗A = \frac{1}{6}(z_{0} + z_{1} + z_{2})k⃗ · s⃗_{0} × s⃗_{1}.  (1)

The gradient of the volume of the body as a function of the vertex v_{0} coordinates is

g⃗ = \frac{1}{6} ((k⃗ · s⃗_{0} × s⃗_{1})k⃗ + (z_{0} + z_{1} + z_{2})k⃗ × s⃗_{1}) .  (2)

§7.5.2. Symmetric content facet integral.

A facet oriented with normal outward from a body makes a contribution to its volume of

V = \frac{1}{3} ∫ ∫_{facet}(xi⃗ + yj⃗ + zk⃗) · d⃗A = \frac{1}{6}v⃗_{0} · v⃗_{1} × v⃗_{2}.  (3)

This can be seen most easily by looking at the facet as the base of a tetrahedron with its fourth vertex at the origin. The integral over the three new faces is zero, so the integral over the original facet must be the volume of the tetrahedron.

The gradient of the volume of the body as a function of the vertex v_{0} coordinates is

g⃗ = \frac{1}{6}v⃗_{1} × v⃗_{2}.  (4)

§7.5.3. Edge content integrals.

Let body b have facet f with outward normal. Suppose f has edge s deemed to have content integrand U⃗. Then b gets volume

V = ∫_{e} U⃗ · d⃗s ≈ ∑_{i} w_{i}U⃗(x⃗_{i})s⃗.  (5)

The volume gradients for the body are (summation convention again)

\begin{aligned}g_{j}(head) & = w_{i}(a_{i}U_{k,j}(x⃗_{i})s_{k} + U_{j}(x⃗_{i})), \\ g_{j}(tail) & = w_{i}((1 − a_{i})U_{k,j}(x⃗_{i})s_{k} − U_{j}(x⃗_{i})). \end{aligned}  (6)

§7.5.4. Volume in torus domain.

The wrapping of the edges across the faces of the unit cell makes the calculation of volumes a bit tricky. Suppose body b has vertices v_{k}. Ideally, we would like to adjust the vertices by multiples of the unit cell basis vectors to get a body whose volume we could find with regular Euclidean methods. Unfortunately, all we know are the edge wraps, i.e. the differences in the adjustments to endpoints of edges. But this turns out to be enough, if we are a little careful with the initial volumes in the datafile.

Let A⃗_{k} be the vertex adjustment for vertex k, and T⃗_{j} be the wrap vector (difference in endpoint adjustments) for facet-edge j. Let m index facets. Then

\begin{aligned}V & = \frac{1}{6} ∑_{facets m}(v⃗_{m0} + A⃗_{m0}) · (v⃗_{m1} + A⃗_{m1}) × (v⃗_{m2} + A⃗_{m2}) \\ & = \frac{1}{6}(S_{1} + S_{2} + S_{3} + S_{4}) \end{aligned}  (7)

where

\begin{aligned}S_{1} & = ∑_{facets m} v⃗_{m0} · v⃗_{m1} × v⃗_{m2}, \\ S_{2} & = ∑_{facet−edges j} v⃗_{j0} · v⃗_{j1} × A⃗_{j2}, \\ S_{3} & = ∑_{facet−verts k} v⃗_{k} · A⃗_{k1} × A⃗_{k2}, \\ S_{4} & = ∑_{facets m} A⃗_{m0} · A⃗_{m1} × A⃗_{m2}. \end{aligned}  (8)

The first of these sums is straightforward. The second sum can be regrouped, pairing the two oppositely oriented facet-edges j and −j for each edge together:

\begin{aligned}S_{2} & = ∑_{edges j} v⃗_{j0} · v⃗_{j1} × A⃗_{j2} + v⃗_{−j0} · v⃗_{−j1} × A⃗_{−j2} \\ & = ∑_{edges j} v⃗_{j0} · v⃗_{j1} × (A⃗_{j2} − A⃗_{−j2}) \\ & = \frac{1}{2} ∑_{edges j} v⃗_{j0} · v⃗_{j1} × (T⃗_{−j2} + T⃗_{j1} − T⃗_{j2} − T⃗_{−j1}) \end{aligned}  (9)

which can be regrouped into a sum over facet-edges, which can be done facet by facet:

S_{2} = ∑_{facet−edges j} v⃗_{j0} · v⃗_{j1} × (T⃗_{j1} − T⃗_{j2}).

The third sum we group terms with a common vertex together, with the inner facet sum over facets around vertex k:

\begin{aligned}S_{3} & = ∑_{vertices k} (v⃗_{k} · ∑_{facetsi} A⃗_{i1} × A⃗_{i2}) \\ & = ∑_{vertices k} (v⃗_{k} · ∑_{facetsi}(A⃗_{i1} − A⃗_{k0}) × (A⃗_{i2} − A⃗_{k})) \\ & = ∑_{vertices k} (v⃗_{k} · ∑_{facets i} T⃗_{i0} × −T⃗_{i2}) \\ & = ∑_{facet−vertices k} v⃗_{k} · T⃗_{k2} × T⃗_{k0}, \end{aligned}  (10)

which again can be done facet by facet.

The fourth sum is a constant, and so only needs to be figured once. Also, it is a multiple of the unit cell volume. Therefore, if we assume the volume in the datafile is accurate to within \frac{1}{12}V_{c}, we can calculate the other sums and figure out what the fourth should be.

The body volume gradient at v_{0} will be

g⃗ = \frac{1}{6} ( ∑_{facets m on v} v⃗_{m1} × v⃗_{m2} + \frac{1}{2}v⃗_{m1} × (T⃗_{m1} − T⃗_{m2}) + T⃗_{m0} × T⃗_{m2}) .  (11)

§7.6. Quantity integrands.

Both facets and edges may contribute to quantities. If the quantities are FIXED, then the quantity gradients as functions of the vvertices must be computed. The facet orientation is the original one. The points evaluated for gaussian quadrature are

x⃗_{i} = ∑_{j} a_{ij}v⃗_{j}.  (1)

§7.6.1. Quantity of a facet.

Let Q⃗ be the facet integrand for the quantity in question. Then the contribution of a facet is

Q = ∫ ∫_{facet} Q⃗ · d⃗A ≈ \frac{1}{2} ∑_{i} w_{i}E⃗(x⃗_{i}) · s⃗_{0} × s⃗_{1}.  (2)

The gradient at vertex v_{m} follows from this. Component notation and the Einstein summation convention are used here, as vector notation gets confusing. Comma subscript denotes partial derivative.

∇Q_{j}(v_{m}) = \frac{1}{2}w_{i}(a_{im}Q_{k,j}(x⃗_{i})(s⃗_{0} × s⃗_{1})_{k}) + (Q(x⃗_{i}) × s⃗_{1})_{j}).  (3)

§7.6.2. Quantity of an edge.

Let q⃗ be a quantity integrand applicable to an edge. Then the edge’s contribution is

Q = ∫_{edge} q⃗ · d⃗s ≈ ∑_{i} w_{i}E⃗(x⃗_{i}) · s⃗.  (4)

The gradients at the endpoints follow from this. Component notation and the Einstein summation convention are used here, as vector notation gets confusing. Comma subscript denotes partial derivative.

\begin{aligned}∇Q_{j}(head) & = w_{i}(a_{i}q_{k,j}(x⃗_{i})s_{k} + q_{j}(x⃗_{i})), \\ ∇Q_{j}(tail) & = w_{i}((1 − a_{i})q_{k,j}(x⃗_{i})s_{k} − q_{j}(x⃗_{i})). \end{aligned}  (5)

§7.7. Constraint projection.

Suppose vertex v has constraints that are the zero sets of functions f_{1}, …, f_{n}.

§7.7.1. Projection of vertex to constraints.

The solution will be by Newton’s method. Several steps may be necessary to get within the tolerance desired of the constraint. Inequality constraints are included if and only if the inequality is violated.

For one step, we seek a displacement δv⃗ that is a linear combination of the constraint gradients and will cancel the difference of the vertex from the constraints:

δv⃗ = ∑_{i} c_{i}∇f_{i}  such that  δv⃗ · ∇f_{j} = −f_{j}(v)  for each j.  (1)

Thus

∑_{i} c_{i}∇f_{i} · ∇f_{j} = −f_{j}(v)  for each j.  (2)

This is a linear system that can be solved for the c_{i} and thus δv⃗.

§7.7.2. Projection of vector onto constraint tangent space.

We want to project a vector F⃗ (say, a force vector) onto the intersection of all the tangent planes of the constraints at a vertex. We project F⃗ by subtracting a linear combination of constraint gradients:

F⃗_{proj} = F⃗ − ∑_{i} c_{i}∇f_{i}  such that  F⃗_{proj} · ∇f_{j} = 0  for each j.  (3)

This leads to the system of linear equations for the c_{i}:

∑_{i} c_{i}∇f_{i} · ∇f_{j} = F⃗ · ∇f_{j}.  (4)

§7.8. Volume and quantity constraints.

This section explains how volume and quantity constraints are enforced. Henceforth in this section, everything referring to the volume of a body can also refer to the value of a quantity constraint. There are two parts to this enforcement on each iteration: a volume restoring motion that corrects for any volume deviations, and a projection of the vertex motions to the subspace of volume-preserving motions. Let g_{bv} be the gradient of the volume of body b as a function of the position of vertex v. We will also use B as a body index. If the bodies fill a torus, then one body is omitted to prevent singular systems of equations.

§7.8.1. Volume restoring motion.

Let the current excess volume of body b be δb (which may be negative, of course). We want a volume restoring motion R⃗_{v}

R⃗_{v} = ∑_{B} c_{B}g⃗_{Bv}  and  ∑_{v} R⃗_{v} · g⃗_{bv} = −δb  for each b  (1)

which leads to the linear system for the c_{B}:

∑_{B} c_{B} ∑_{v} g⃗_{Bv} · g⃗_{bv} = −δb  for each b.  (2)

§7.8.2. Motion projection in gradient mode.

The motion in this mode is the scale factor times the force on a vertex. Let F⃗_{v} be the total force at v. We want a projected force F⃗_{vproj} such that

F⃗_{vproj} = F⃗_{v} − ∑_{B} a_{B}g⃗_{Bv}  and  ∑_{v} F⃗_{vproj} · g⃗_{bv} = 0  for each b  (3)

which leads to the linear system for a_{B}

∑_{B} a_{B} ∑_{v} g⃗_{Bv} · g⃗_{bv} = ∑_{v} F⃗_{v} · g⃗_{bv}  for each b.  (4)

The coefficients a_{B} are the pressures in the bodies when the surface is in equilibrium.

§7.8.3. Force projection in mean curvature mode.

The motion in this mode is the scale factor times the force on a vertex divided by the area of the vertex. Let F⃗_{v} be the total force at v and A_{v} its area. We want a projected force F⃗_{vproj} such that

F⃗_{vproj} = F⃗_{v} − ∑_{B} a_{B}g⃗_{Bv}  and  ∑_{v} F⃗_{vproj} · g⃗_{bv}/A_{v} = 0  for each b  (5)

which leads to the linear system for a_{b}

∑_{b} a_{b} ∑_{v} g⃗_{Bv} · g⃗_{bv}/A_{v} = ∑_{v} F⃗_{v} · g⃗_{bv}/A_{v}  for each b.  (6)

The coefficients a_{b} are the pressures in the bodies even when the surface is not in equilibrium. At equilibrium, the pressures calculated with or without the mean curvature option are equal.

§7.9. Iteration.

This section explains what happens during each iteration (g command). Δt is the current scale factor.

§7.9.1. Fixed scale motion.

1. Do diffusion if called for. For each facet f, let A_{f} be the area of the facet. Let m_{1} and P_{1} be the prescribed volume (mass) of a body on one side of f, and m_{2} and P_{2} those of the body on the other side. Let κ be the diffusion constant. If there is no body on one side, the pressure there is zero or ambient pressure for the ideal gas model. The mass transfers are:

dm_{1} = −κ(P_{1} − P_{2})A_{f}Δt  and  dm_{2} = −κ(P_{2} − P_{1})A_{f}Δt.  (1)

2. All the forces on vertices are calculated, as detailed in the Energies and Forces section above. If the mean curvature option is active, each force is divided by the vertex’s associated area.

3. All forces at FIXED vertices are set to zero. Forces on constrained vertices are projected to the constraint tangent spaces. Forces on vertices on boundaries are converted to forces in boundary parameter spaces.

4. The volume restoring motion is calculated, and the forces are projected according to body volume constraints.

5. Jiggle if jiggling is toggled on.

6. Move vertices by volume restoring motion and current scale times force.

7. Enforce constraints on vertices. Do up to 10 projection steps until vertex is within tolerance of constraints.

8. Recalculate volumes, pressures, areas, and energy.

§7.9.2. Optimizing scale motion.

The idea is to find three scale factors that bracket the energy minimum and then use quadratic interpolation to estimate the optimum scale factor. The first four steps are the same as for fixed scale motion.

5. Save current coordinates.

6. Move vertices by volume restoring motion and current scale times force and enforce constraints. Recalculate energy.

7. Restore coordinates, double scale factor, move surface, and recalculate energy. If this energy is higher, keep cutting scale factor in half until energy increases again. If it is lower, keep doubling scale factor until energy increases.

8. Eventually get three scale factors s_{1}, s_{2}, s_{3} with energies e_{1}, e_{2}, e_{3} bracketing a minimum. The scale factors are successive doublings: s_{2} = 2s_{1}, s_{3} = 2s_{2}. Do quadratic interpolation to find approximate optimum:

s_{opt} = 0.75s_{2}(4e_{1} − 5e_{2} + e_{3})/(2e_{1} − 3e_{2} + e_{3}).  (2)

9. Jiggle if jiggling is on.

10. Move vertices by volume restoring motion and optimum scale times force.

11. Project to constraints.

12. Recalculate volumes, pressures, areas, and energy.

§7.9.3. Conjugate gradient mode.

The U command toggles conjugate gradient mode. The conjugate gradient method does not follow the gradient downhill, but makes an adjustment using the past history of the minimization. It usually results in much faster convergence. At iteration step i, let S_{i} be the surface, E_{i} its energy, F⃗_{i}(v) the force at vertex v, and h⃗_{i}(v) the “history vector” of v. Then

h⃗_{i}(v) = F⃗_{i}(v) + γh⃗_{i−1}  (3)

where

γ = \frac{∑_{v} F⃗_{i}(v) · F⃗_{i}(v)}{∑_{v} F⃗_{i−1}(v) · F⃗_{i−1}(v)}.  (4)

Then a one-dimensional minimization is performed in the direction of h⃗_{i}, as described in steps 6,7,8 above. It is important that all volumes and constraints be enforced during the one-dimensional minimization, or else the method can go crazy. The history vector is reset after every surface modification, such as refinement or equiangulation.

§7.10. Equiangulation.

Triangulations work best when the facets are close to equilateral (that is, equiangular). Given a set of vertices, how do you make them into a triangulation that has triangles as nearly as possible equilateral? In the plane, the answer is the Delauney triangulation, in which the circumcircle of each triangle contains no other vertex. See [S]. It is almost always unique. It can be constructed by local operations beginning with any triangulation. Consider any edge as the diagonal of the quadrilateral formed by its adjacent triangles. If the angles of the two vertices off of the edge add to more than π, then the circumcircle criterion is violated, and the diagonal should be switched to form a replacement pair of facets.

Suppose now we have a triangulation of a curved surface in space. For any edge with two adjacent facets, we switch the edge to the other diagonal of the skew quadrilateral if the sum of the angles at the off vertices is more than π. Let a be the length of the common edge, b and c the lengths of the other sides of one triangle, and d and e the lengths of the other sides of the other triangle. Let θ_{1} and θ_{2} be the off angles. Then by the law of cosines

a^{2} = b^{2} + c^{2} − 2bc cos θ_{1},  a^{2} = d^{2} + e^{2} − 2de cos θ_{2}.  (1)

The condition θ_{1} + θ_{2} > π is equivalent to cos θ_{1} + cos θ_{2} < 0. So we switch if

\frac{b^{2} + c^{2} − a^{2}}{bc} + \frac{d^{2} + e^{2} − a^{2}}{de} < 0.  (2)

This is guaranteed to terminate since a switch reduces the radii of the circumcircles and there are a finite number of triangulations on a finite number of vertices. When using the u command to equiangulate, you should use it several times until no changes happen. If you do get a genuine infinite loop, report it as a bug.

§7.11. Area normalization.

By default, the gradients of quantities like energy and volume are calculated simply as the gradient of the quantity as a function of vertex position. This gives the force on a vertex, for example. But sometimes it is desirable to have force per area instead, for example in directly estimating mean curvature. The mean curvature vectorfield h⃗ of a surface is defined to be the gradient of the area of the surface in the sense that if the surface is deformed with an instantaneous velocity u⃗ then the rate of change of surface area is

\frac{dA}{dt} = ∫ ∫_{surface} u⃗ · h⃗dA.  (1)

According to the triangulation formulation, with F⃗_{v} being the force on vertex v, the rate of change is

\frac{dA}{dt} = ∑_{v} u(v) · −F⃗_{v}.  (2)

Hence as approximation to h⃗, we take

h⃗ = −F⃗_{v}/dA_{v}  (3)

and take the area dA_{v} associated with a vertex to be one-third of the total areas of the facets surrounding the vertex. Since each facet has three vertices, this allocates all area. Command “a” can be used to toggle area normalization.

§7.12. Hidden surfaces.

This section describes the algorithm used by the painter graphics interface to hide hidden surfaces.

1. Each facet has the viewing transformation applied to it.

2. Each facet has its maximum and minimum in each dimension recorded (bounding box).

3. Facets are sorted on their maximum depth (rear of bounding box).

4. List is traversed in depth order, back to front. If a facet’s bounding box does not overlap any others’, it is drawn. If it does overlap, a linear program tries to find some point where one facet is directly in front of the other. If there is such a point and the order there is backward, the facets are switched on the list. The traversal then restarts there. If it makes over 10 switches without drawing a facet, it gives up and draws the current facet anyway.

§7.13. Extrapolation.

The final energy value at refinement step n is saved as E_{n}. Assuming the difference from minimum energy follows a power law of unknown power, three successive values are used to extrapolate to the minimum energy E_{0}:

\begin{aligned}E_{n} & = E_{0} + ac^{b}, \\ E_{n−1} & = E_{0} + a(2c)^{b}, \\ E_{n−2} & = E_{0} + a(4c)^{b}. \end{aligned}  (1)

Hence

E_{0} = E_{n} − \frac{(E_{n} − E_{n−1})^{2}}{E_{n} − 2E_{n−1} + E_{n−2}} .  (2)

§7.14. Curvature test.

This only applies to the quadratic model.

§7.15. Annealing (jiggling).

This happens at each command j, or at each iteration if the jiggling option (command J)is activated. Each coordinate of each non-FIXED vertex is moved by

δx = gT L  (1)

where g is a random value from the standard gaussian distribution (calculated from the sum of five random values from the uniform distribution on [0,1]), T is the current temperature, and L is a characteristic length that starts as the diameter of the surface and is cut in half at each refinement.

§7.16. Long wavelength perturbations (long jiggling).

This is command jj. An amplitude A⃗ and a wavelength L⃗ are chosen at random from a unit sphere by picking random vectors in a cube of side 2 centered at the origin and rejecting attempts that are outside the unit sphere. A⃗ and L⃗ are then multiplied by the size of the surface to give them the proper scale. The wavelength is inverted to a wavevector w⃗. The amplitude is multiplied by the temperature. A random phase ψ is picked. Then each vertex v⃗ is moved by A⃗ sin(v⃗ · w⃗ + ψ).

§7.17. Homothety.

This scales the total volume of all bodies to 1 by multiplying each coordinate by

( ∑_{bodies} V_{body})^{−1/3} .  (1)

In the string model, the factor is

( ∑_{bodies} V_{body})^{−1/2} .  (2)

§7.18. Popping non-minimal edges.

All facets are assumed to have equal tension. This looks for edges with more that three facets that are not fixed or on boundaries or constraints. When found, such an edge is split with a new facet in between. The two closest old facets are attached to the new edge. This is repeated until only three facets are on the original edge. Each split is propagated along the multiple junction line as far as possible. If it is impossible to propagate beyond either endpoint, the edge is subdivided to provide a vertex which can be split.

§7.19. Popping non-minimal vertex cones.

String model.

Any vertex with four or more edges is tested. The edges into the vertex are separated into two contiguous groups. The endpoints of one group are transfered to a new vertex, and a short edge is introduced between the new and old vertex. All groupings are tested to see which one pulls apart the most. This method works when the edges have non-equal tensions.

Soapfilm model.

All facets are assumed to have equal tension. This assumes that all edges have at most three facets adjoining, so edge popping should be done first. Here the facet and edge structure around each vertex is analyzed to find which have the wrong topology. The analysis is done by looking at the network formed by the intersection of the facets and edges containing the vertex with the surface of a small sphere around the vertex. The numbers of sides of the cells of the network are counted. A simple plane vertex has two cells of one side each. A triple edge vertex has three cells of two sides each. A tetrahedral point has four cells with three sides each. Any other configuration is popped.

The popping itself is done by replacing the vertex with a hollow formed by truncating each cell-cone except the cell with the largest solid angle. In case the network is disconnected, the solid angles of all the cells will add up to over 4π. Then the vertex is duplicated and the different components are assigned to different vertices. This lets necks shrink to zero thickness and pull apart.

§7.20. Refining.

First, all edges are subdivided by inserting a midpoint. Hence all facet temporarily have six sides. For edges on constraints, the midpoints get the same set of constraints, and are projected to them. For edges on boundaries, the parameters of the midpoint are calculated by projecting the vector from the edge tail to the midpoint back into the parameter space and adding that to the tail parameters. This avoids averaging parameters of endpoints, which gives bad results when done with boundaries that wrap around themselves, such as circles. The new edges and the midpoint inherits their attributes from the edge.

Second, each facet is subdivided into four facets by connecting the new midpoints into a central facet. The new edges inherit their attributes from the facet.

§7.21. Removing tiny edges.

The edge list is scanned, looking for edges shorter than the cutoff length. It is NOT removed if any of the following are true:

1. Both endpoints are FIXED.

2. The edge is FIXED.

3. The boundary and constraint attributes of the endpoints and the edge are not identical. For an edge that passes these tests, one endpoint is eliminated. If one endpoint is FIXED, that is the one kept. All the facets around the edge are eliminated (collapsing their other two edges to one), and the edge itself is eliminated.

§7.22. Weeding small triangles.

This is meant to help eliminate narrow triangles that tiny edge removal can’t because they don’t have short edges. This procedure finds the shortest edge of a triangle and eliminates it by the same process as the regular tiny edge removal.

§7.23. Vertex averaging.

For each vertex, this computes a new position as the area-weighted average of the centroids of the facets adjoining the vertex. FIXED vertices or vertices on boundaries are not moved, nor are vertices on triple or higher edges. If vertex v is connected by edges to vertices v_{1}, …, v_{n}, then the new position is calculated as

v⃗_{avg} = \frac{∑_{facets} area · x⃗_{centroid}}{∑_{facets} area} .  (1)

The volume one one side of all the facets around the vertex calculated as a cone from the vertex is

V = ∑_{facets f} v⃗ · N⃗_{f}  (2)

where N⃗_{f} is the facet normal representing its area. The total normal N⃗ is

N⃗ = ∑_{facets f} N⃗_{f}.  (3)

To preserve volume, we subtract a multiple λ of the total normal to the average position:

(v⃗_{avg} − λN⃗) · N⃗ = v⃗ · N⃗,  (4)

so

λ = \frac{v⃗_{avg} · N⃗ − v⃗ · N⃗}{N⃗ · N⃗} .  (5)

Then the new vertex position is

v⃗_{new} = (v⃗_{avg} − λN⃗).  (6)

Constrained vertices are then projected to their constraints.

The “rawv” command will do vertex averaging as just described, but without conserving volumes. The “rawestv” command will also ignore all restrictions about like constraints and triple edges, but it will not move fixed points or points on boundaries. Points on constraints will be projected back to them.

§7.24. Zooming in on vertex.

First, all vertices beyond the cutoff distance from the given vertex are deleted. Then all edges and facets containing any deleted vertices are deleted. Any remaining edge that had a facet deleted from it is made

FIXED.

§7.25. Mobility and approximate curvature.

This section derives the resistance and mobility matrices used in approximate curvature. The derivation is done in the full generality of vertex constraints, boundaries, volume and quantity constraints, and effective area.

First, we must define the inner product or metric on global vectors, which is the resistance matrix. Let W⃗ be a global vectorfield, that is, W⃗ has a component vector W⃗_{v} at each vertex v. At each vertex v, let U⃗_{v,i} be a basis for the degrees of freedom of motion of vertex v due to level set constraints or boundaries.

Consider one edge e in the string model. Let P be the projection operator onto the normal of e if effective area is in effect, and let P be the identity map if not. Then the inner product of two vector fields is

\begin{aligned}< W⃗ , V⃗ > & = L ∫_{0}^{1} P ((1 − t)W⃗_{a} + tW⃗_{b}) · P ((1 − t)V⃗_{a} + tV⃗_{b})dt \\ & = L(\frac{1}{3}P W⃗_{a} · P V⃗_{a} + \frac{1}{6}P W⃗_{a} · P V⃗_{b} + \frac{1}{6}P W⃗_{b} · P V⃗_{a} + \frac{1}{3}P W⃗_{b} · P V⃗_{b}). \end{aligned}  (1)

If the vectors are expressed in terms of the bases U_{v} with coefficient column matrices W_{v} and V_{v}, then

\begin{matrix}< W⃗ , V⃗ >= L(\frac{1}{3} ∑_{i,j} W_{ai}U⃗_{a,i} · P U⃗_{a,j}V_{aj} & + \frac{1}{6} ∑_{i,j} W_{ai}U⃗_{a,i} · P U⃗_{b,j}V_{bj} \\ +\frac{1}{6} ∑_{i,j} W_{bi}U⃗_{b,i} · P U⃗_{a,j}V_{aj} & + \frac{1}{3} ∑_{i,j} W_{bi}U⃗_{b,i} · P U⃗_{b,j}V_{bj}). \end{matrix}  (2)

So the components of the resistance matrix S are

\begin{aligned}S_{ai,aj} & = \frac{L}{3}U⃗_{a,i} · P U⃗_{a,j}, \\ S_{ai,bj} & = \frac{L}{6}U⃗_{a,i} · P U⃗_{b,j}. \end{aligned}  (3)

In the soapfilm model, for a facet of area A,

\begin{aligned}S_{ai,aj} & = \frac{A}{6}U⃗_{a,i} · P U⃗_{a,j}, \\ S_{ai,bj} & = \frac{A}{12}U⃗_{a,i} · P U⃗_{b,j}. \end{aligned}  (4)

The total matrix S is formed by adding the contributions from all edges. Note that S has a block on its diagonal for each vertex and a block for each pair of edges joined by an edge. Hence it is reasonably sparse. S is a symmetric positive semidefinite matrix, positive definite if effective area is not in effect.

Now suppose we have a form H to convert to a vector. H is represented as a covector H_{v} at each vertex v. The mobility matrix M = (M_{bj,ai}) is the inverse of S. The first step is to find the coefficients in the dual basis of U by contracting:

h_{a,i} =< U⃗_{a,i}, H_{v} > .  (5)

The U basis components of the vector dual to H are then

c_{b,j} = M_{bj,ai}h_{a,i}.  (6)

Actually, the components are found by solving the sparse system

S_{ai,bj}c_{a,i} = h_{a,i}.  (7)

The foregoing ignores volume and quantity constraints. Suppose B_{k} are gradient forms for the volume and quantity constraints. We need to project the motion in such a way as to guarantee that when the velocity is zero, the energy gradient is a linear combination of constraint gradients.

proj(H) = H − λ_{k}B_{k}.  (8)

The projected velocity is

V = M(H − λ_{k}B_{k}).  (9)

The condition for motion tangent to the constraints is

B_{m}M(H − λ_{k}B_{k}) = 0  for each m.  (10)

Therefore

B_{m}MB_{k}λ_{k} = B_{m}MH,  (11)

and so

λ_{k} = (B_{m}MB_{k})^{−1}B_{m}MH.  (12)

§7.26. Torus domain.

Special considerations needed for torus surfaces: Edge vector calculation using wrap flags. Volume and volume gradient calculation. Displaying as connected bodies or clipped unit cell. No quadratic representation. Torus duplication command.

Chapter 8. Helpful hints and notes.

When writing a new datafile, use the d command to check that it is being read as you think it should be.

Run at low resolution before refining. Having many triangles not only takes a long time to calculate, but motion can propagate only one triangle per iteration. Also run the linear model first, and use the quadratic model for final polishing.

Watch out for edges shrinking to zero length and triangles to zero area; they can stall progress. If the optimal scale factor heads toward 0, try equianglation (u), vertex averaging (V), small triangle weeder (w), or the tiny edge remover (t). Always precede and follow w with t.

Facets entirely on constraints tend not to move, and moving free boundaries tend to overrun such facets, leading to problems. Use constraint integrals rather than facets on constraints where possible.

If your edges on curved constraints try to shortcut the curve, try adding an energy integrand term to the constraint that is zero on the constraint and positive off it. Integrands are evaluted at several points in the interior of an edge, not the endpoints, so this gives an energy penalty for going away from the edge. Hopefully, it will not change the nature of your solution.

Use constraints for two-dimensional free boundaries, and use parametric boundaries for wire frames. Use the d command to save your surface at critical points. Run checks (c) after modifying surface just to make sure everything worked right. MinneView users should periodically run smem -i from the shell to clear out garbage shared memory.

Do this only when nobody is using MinneView.

Chapter 9. Bugs.

There are no known outright bugs at present. but they undoubtably exist. When you run across one, I would like to hear about it. Bug reports should be submitted by email to brakke@geom.umn.edu. Please include the Evolver version number, a description of the problem, the initial data file, and the sequence of commands necessary to reproduce the problem.

There are a few shortcomings, however: Not all features are implemented in all models. Vertex and edge popping is not elegant or complete. Zero length edges and zero area triangles stall things. Surfaces can intersect each other without knowing it. There is no save in binary format. There once was, but the ASCII dump file is shorter, and keeping two save/load formats updated got to be a pain. There are problems in detecting errors in Cheyenne initialization.

Chapter 10. Version history.

Version 1.0 August 4, 1989 First public version.

Version 1.01 August 22, 1989 Various bug fixes. Constraint integral specification changed to have just components, not density. Facet-edge specification made obsolete (but still legal). Vertex motion adjustment for quadratic model. Zoom factor changed from 2 to 1.2. Initial datafile reading revamped to use lex and yacc. Has simple macros. Comments must be delimited.

Expressions now in algebraic form, not RPN.

Version 1.10 October 18, 1989 More bug fixes. Graphics driver for HP98731 contributed by Eric Haines. Datafile reading: Equations permitted as constraints. Backslash line splicing. Constant folding in expressions.

Version 1.11 May 18, 1990 More bug fixes. x,y,z accepted as synonyms for x1,x2,x3

Version 1.12 May 25, 1990 Datafile coordinates may be given as expressions. Constraints revamped. Number of constraints per element raised to sizeof(int). One-sided constraints

may be applied selectively. GLOBAL constraints automatically apply to all vertices (they count in number limit).

Version 1.2 June 17, 1990 Comments removed from input stream before lex analyzer. So comments in macros safe now. Scale factor upper bound adjustable when toggling with m command. Adjustable constants implemented. Syntax in data file: PARAMETER name = constant-expression These

are treated as constants in expressions, but can be changed with the A command. Useful for dynamically changing contact angles and other stuff in constraint expressions.

Symmetric content evaluation added for linear model. Volumes evaluated as surface integral of (xi⃗+yj⃗+ zk⃗)/3 rather than of zk⃗. Use keyword SYMMETRIC CONTENT in datafile. Permits more accurate evaluation of lunes.

Return value of sprintf not used, as this varies among systems. Arbitrary surface energy integrands added. Vector field integrated over designated facets. Syntax

surface energy n e1: expression e2: expression e3: expression Designate facet by following facet definition with “energy n”, in same manner as giving constraints. Linear model only. This useful for changing direction of gravity, by putting in Divergence Theorem equivlent surface integrand for gravitational potential energy and putting adjustable constant in for direction.

The ’q’ or ’x’ exit command gives you a chance to load another datafile, continue with the current configuration, or exit.

You can change the initial upper limit on the ’scale’ factor by putting a line in the data file: SCALE LIMIT value value must be a number (no expression). The ’V’ command (for Vertex averaging) will replace each unconstrained vertex by the average of its

neighboring vertices (connected by edges). Good for uncramping skinny triangles sometimes. TENSION allowed as synonym for DENSITY in setting facet surface tension. Binary save/reload disabled, since it is out of date and a binary file turns out twice as large as the ascii

dump file.

Version 1.21 June 30, 1990

Shared memory interface added for MinneView (which is a public domain 3D graphics viewer for Irises written by the Geometry Project at the Minnesota Supercomputer Institute.)

Histogram added for ridge notcher (option ’n’). Ported to NeXT (no screen graphics, but PostScript output files can be displayed). Automatic energy recalculation added after all options that change energy.

Version 1.3 July 30, 1990 Reptition counts before letters in graphics commands. Clockwise (c) and counterclockwise (C) rotations added to graphics. Printout during quadratic search for optimum scale suppressed. ’i’ command added for information. ’v’ reports just volumes. Extrapolation command ’e’ fixed. Constant expressions permitted wherever real value needed. ’F’ command to log commands to file. -ffilename command line option to read commands from file. Take commands from stdin afterwards. Dump file real values printed to 15 decimal places so accuracy not lost. Adjustable accuracy on edge integrals via INTEGRAL_ORDER keyword in datafile for setting order of

gaussian quadrature.

Datafile now allows specification of constraint tolerance with CONSTRAINT_TOLERANCE. All two-word keywords made into single words by connecting the words with ’ _’. Parsing of datafile continues after errors. ’b’ command permits editing body volumes. ’v’ command prints ’none’ for prescribed volume of those bodies with no prescribed volume. Manual in TEXformat with PostScript figures. CONVEX gap energies and forces fixed up. k = 1 best approximation to area. Vertex popping fixed to handle disjoint components of tangent cones. Also not to screw up its data

structures after a modification.

Version 1.31 Added long-wavelength random jiggle (command jj). Fixed simultaneous Minneview to handle torus display options.

Simultaneous MinneView can be stopped and restarted. Fixed bug in gaussian quadrature of content integrals. Fixed equiangulation to work on arbitrarily small triangles. Histograms adjusted to current length and area order of magnitude. Zero area triangle test area cutoff adjusted to scale of surface. Vertex averaging improved to preserve volumes. Small edge removal bug fixed. In collapsing a facet, will preferentially keep an edge on a constraint. VOLCONST adjustment to body volumes added to datafile. “Quantities” added. A quantity is a sum of vector integrals over facets and edges. Can be simply tallied

for informational purposes, or can be used as mathematical constraints with fixed values. Good for center of mass, moment of inertia, magnetic flux, etc.

Conjugate gradient energy minimization added. Use the U command to toggle between gradient descent and this.

Version 1.4 posted August 20, 1990

Version 1.41 September 22, 1990 Bug fixes on vertex popping and element list management. Long jiggle command jj lets user put in own numbers, use random numbers, or use previous numbers.

Version 1.42 Bug fixes: segment violation on error message at end of data file. STRING model not graphing third dimension On SOAPFILM model, edges without facets are displayed. Useful for showing wire boundaries beyond

surface. Such edges will generate warnings when datafile is read, however.

Version 1.5 May 15, 1991 dimensional surfaces can live in N-dimensional space. See SPACE_DIMENSION. Simplex model added to represent k-dimensional surfaces in N-dimensional space. Embryonic query language added. Does list and show. Also set and unset attributes. Background metric on space added for string model only. Commands can be read from a file with read "filename"

Version 1.6 June 20, 1991 Riemannian metric extended to all dimensional surfaces. Refine, delete added to query language. List query has same format output as datafile dump. Quotient spaces added. Piping output added to commands and queries.

Version 1.63 July 21, 1991

SET FACET COLOR query added. Colors include CLEAR. Also SET FACET TRANSPARENCY for Iris and like.

Squared mean curvature added as possible energy. Queries can refer to elements by original number of parent in datafile. Adjustable constants can take values from a file according to ID number of element applied to.

Version 1.64 August 2, 1991

Revised notch command (n) to subdivide adjacent facets instead of edge itself. Supposed to be followed with equiangulation.

Put in K command to subdivide longest edges of skinny triangles (as judged by their smallest angle). Added VALENCE attribute for edges for queries. Is the number of facets on an edge.

Version 1.65 August 20, 1991

Added minimization by using Newton’s method on Hessian matrix of energy. Only for no-constraint area minimization with no other energies. Command “hessian”.

NeXT version given graphical interface. User-defined functions of coordinates added. See userfunc.c.

Version 1.76 March 20, 1992 Autopopping and autochopping in string model for automatic evolution. Phase-dependent grain boundary energies. Approximate polyhedral curvature. Stability test for approximate curvature. Squared Gaussian curvature as part of energy only, not force. “system” command to execute shell commands. “check increase” to guard against blowup during interation. “effective area” to count only resistance to motion normal to surface. Runge-Kutta iteration.

Version 1.80 July 25, 1992 Command and query language much extended. geomview interface added. Fixed area added as a constraint. Multiple viewing transforms can be specified in the datafile so one fundamental region of a symmetric

surface can be displayed as the whole surface.

Commands can be included at the end of the datafile, introduced by the keyword READ.

Version 1.83 September 9, 1992

Some alternate definitions of squared curvature added. Invoked by “effective area ON | OFF” or

“normal curvature ON | OFF”.

Version 1.84 September 10, 1992 Shaded colors added to xgraph and cheygraph.

Version 1.85 September 29, 1992 Restriction of motion to surface normal added. Toggle “tt normal motion”. Squared mean curvature, gaussian curvature, and squared gaussian curvature extended to surfaces with

boundary.

Datafile element attribute “bare” added for vertices and edges in soapfilm model so they won’t generate erroneous warnings.

Force calculation added for squared gaussian curvature, so it can be used in the energy. All prompts that require real value responses now accept arbitrary expressions.

Version 1.86 October 19, 1992 User-defined mobility added, both scalar and tensor forms. Default squared curvature works for 2-surfaces in R^{n}.

Version 1.87 October 27, 1992 “close show” command added to close show window (the native graphics window, not geomview). Graphics command checks string for illegal characters before doing any transformations. Dihedral angles now work for 2-surfaces in any dimension. Permanent variable assignments may be made with “::=” instead of “:=”. Such assignments will not

be forgotten when a new surface is begun.

Conditional expressions of C form added: expr ? expr : expr. Useful for patching constraints together.

Version 1.88 December 16, 1992 “SET BACKGROUND color” command added for native graphics. View transformation generators and expressions added. Exact bounding box calculated for PostScript files.

Version 1.88a January 6, 1993 Default constraint tolerance lowered from 1e-5 to 1e-12. Fixed bug in volume constraint calculation introduced in 1.88.

Version 1.89 February 18, 1993

Postscript draws fixed and boundary edges in interior of surface. All internal graphics should be consistent in the special edges they draw (bare edges, triple edges, etc.).

Viewing matrix can be read from datafile and will be dumped. Keyword VIEW MATRIX Mod operator ‘%’ added, and functions floor(), ceil(). ‘rebody’ command added to recalculate connected bodies after neck pinching and any other body dis-

ruption.

If squared mean curvature part of energy, then squared mean curvature at a vertex is available as a query attribute as ”sqcurve”.

These quantities can now be used in command expressions: vertex count, edge count, facet count, body count, total energy, total area, total length, scale.

Dump file records defined procedures in ‘read’ section at end. Knot energies added, both conducting and insulating wire. ‘dissolve’ command added to erase elements and leave gaps in surface, unlike delete command, which

closes gaps. Can only dissolve elements not neede by higher dimensional elements.

Command repeat numbers have been restricted to just three types of commands: 1. single letter commands that don’t have optional arguments (l,t,j,m,n,w have optional arguments) 2. command list in braces 3. user-defined procedure names This is to prevent disasters like list vertex 1293 which before would produce 1293 full lists of all vertices.

“dump” without argument will dump to default file name, which is datafile name with .dmp extension. SIGTERM is caught and causes dump to default dump file and exit. Useful for interrupting scripts

running in the background with kill -TERM. Likewise for SIGHUP.

Version 1.90 April 2, 1993

Conjugate gradient ON/OFF state saved in dump file. Note that conjugate gradient history vector is not saved.

Notching and “dihedral” attribute apply to vertices in the string model. FOREACH iterator added. Syntax: FOREACH element [ name ] [WHERE expr] DO command LOAD command added. Syntax: LOAD “filename”. Useful for starting new surfaces, especially in

scripts.

PRINTF command added for formatted printing. Syntax: PRINTF ”format string”,expr,expr,... expr is floating point, so use

String variables have been added. Can be used where quoted strings needed. Can be assigned to. SPRINTF is version of PRINTF giving string output.

A view transformation matrix in the datafile may be preceded by ”color n” to give that transform a color (overrides any facet color).

In queries, element attribute “oid” added, which returns a signed version of id. Many knot energies added. Also a “hooke energy” that keeps edges near a uniform length. PostScript output optionally includes color.

Version 1.91 May 31, 1993

Two sides of facets can have different colors. ’COLOR’ applies to both sides, ’FRONTCOLOR’ and ’BACKCOLOR’ to different sides.

Attributes of individual named elements can be set inside loops, i.e. foreach facet ff do set ff color red Every time a command changed a global variable, the surface was being recalculated. This slowed down

scripts immensely. So now the only variables that cause recalculation are 1) adjustable parameters defined in the datafile 2) quantity moduli and parameters

History commands now echoed. Surface area can be minimized by minimizing Dirichlet integral, according to scheme of Polthier and

Pinkall. Command ’dirichlet’.

To reduce need for explicit line-splicing on long commands, the parser now keeps track of depth of brace and parenthesis nesting, and will call for more input if a line ends inside nest. So if you want to type a multiline command, start with ’’ and end with ’’ many lines later. Also does auto line-splicing if certain tokens are last token in line (such as ’+’).

‘facet knot energy fix’ method added.

Command assignment fixed to assign only one command. so “ggg := g; g” will be the same as “{ggg

:= g} ; g” and not “ggg := {g;g} ”.

Queries can run through edges and facets adjacent to a vertex, and facets of a body, as in “list

vertices vv where max(vv.facet,color==red) > 0”

Improved NeXT terminal interface. -u option for no graphics, -t for terminal and graphics. view 4d command to toggle sending full 4D coordinates to geomview. Default is OFF.

Version 1.92 July 1, 1993 SGI parallelism enabled for named quantity calculations. method-instance scheme introduced. Torus periods can be specified with expressions using adjustable parameters. Verbs (list,refine,delete,dissolve) may apply to single elements: foreach edge ee do refine ee FIX and UNFIX may be used as verbs:

fix vertices where on constraint 1

Toggle command names may now be used as boolean values in expressions in commands. Also new boolean read-only variables: torus, torus filled, symmetry group, simplex representation. New numeric readonly variables: space dimension, surface dimension, integration order.

Mac version repeated commands interruptable with control-’.’.

Macintosh and Dos versions pipe to a file instead of a command:

Enter command:  list vertices | "filename"

For a torus domain, the torus periods may be specified using expressions with parameters, so the fundamental cell may be changed interactively. Do a recalc after changing such a parameter to update the torus periods.

gv binary toggle for binary/ascii data to geomview. Default ON for binary, which is faster. Ascii mode useful for debugging.

Surface Evolver Manual

Chapter 11. Bibliography.

[A] F. Almgren, “Minimal surface forms,” The Mathematical Intelligencer, vol.4, no. 4 (1982), 164–172.

[AT] F. Almgren and J. Taylor, “The geometry of soap films and soap bubbles,” Scientific American, July 1976, 82–93.

[B1] K. Brakke, The motion of a surface by its mean curvature, Princeton University Press, Princeton, NJ (1977).

[B2] K. Brakke, “The surface evolver,” Experimental Mathematics, to appear.

[B3] K. A. Brakke, “Minimal surfaces, corners, and wires,” Journal of Geometric Analysis 2 (1992) 11-36.

[B4] K. A. Brakke, “The opaque cube problem video,” Computing Optimal Geometries, J. E. Taylor, ed., American Mathematical Society, Providence RI, 1991.

[B5] K. A. Brakke, “Grain growth with the Surface Evolver,” Video Proceedings of the Workshop on Computational Crystal Growing, J. E. Taylor, ed., American Mathematical Society, Providence RI, 1992.

[C] M. Callahan, P. Concus and R. Finn, “Energy minimizing capillary surfaces for exotic containers,” Computing Optimal Geometries, American Mathematical Society, Providence RI, 1991.

[DL] D. Dobkin and M. Laszlo, “Primitives for the manipulation of three-dimensional subdivisions,” technical report CS-TR-089-87, Department of Computer Science, Princeton University, Princeton, New Jersey. (April 1987).

[HKS] L. Hsu, R. Kusner, and J. Sullivan, “Minimizing the squared mean curvature integral for surfaces in space forms,” Experimental Mathematics 1 (1991) 191–208.

[MH] H. Mittelmann, “Symmetric capillary surfaces in a cube”, Math. Comp. Simulation, to appear.

[MF1] F. Morgan, “Harnack-type mass bounds and Bernstein theorems for area-minimizing flat chains modulo ν,” Comm. Part. Diff. Eq. 11, (1986) 1257–1283.

[MF2] F. Morgan, “Surfaces minimizing area plus length of singular curves,” preprint.

[P] W. Press, B. Flannery, S. Teukolsky, and W. Vetterling, Numerical Recipes in C, Cambridge University

Press, New York, 1988.

[RSB] L. M. Racz, J. Szekely, and K. A. Brakke, “On some meniscus problems in materials processing,”

submitted to Trans. Iron and Steel Inst.

[S] R. Sibson, “Locally equiangular triangulations”, Comput. J. 21 (1978) 243-245.

[T1] J. E. Taylor, “The structure of singularities in soap-bubble-like and soap-film-like minimal surfaces,” Ann. Math. 103 (1976), 489-539.

[T2] J. E. Taylor, “Constructing crystalline minimal surfaces,” Seminar on Minimal Submanifolds, E. Bombieri, ed., Annals of Math. Studies 105 (1983), 271–288.

[T3] J. E. Taylor, “Constructions and conjectures in crystalline nondifferentiable geometry,” Proceedings of the Conference in Differential Geometry, Rio de Janiero, 1988, Pittman Publishing Ltd.

[Te] J. Tegart, “Three-dimensional fluid interfaces in cylindrical containers,” AIAA paper AIAA-91-2174, 27th Joint Propulsion Conference, Sacramento, CA, June 1991.

[TW] W. Thompson (Lord Kelvin), “On the division of space with minimum possible error,” Acta Math. 11(1887), 121–134.

Chapter 12. Index.

+ graphics command 71

- graphics command 71 0 command 59 ? graphics command 71 acquisition 2 adjustable constants 9,15,20,38,43 aggregate expressions 62 ambient pressure 32,34,39,45,59 annealing 46,92 approximate curvature 37,67

APPROXIMATE_CURVATURE 46
APPROX_CURV 68

area normalization 35,36,57,58,91

AREAWEED 66
AREA_FIXED 68
AREA_FIXED 48
AREA_NORMALIZATION 68
AREA_NORMALIZATION 45

area 61 assignment 64 attributes 8

AUTOCHOP 68

autochop 49 autochop command 68

AUTOPOP 68

autopopping 49 autopop command 68

AUTORECALC 68
AVG 62

A command 9,11,16,21,38,56 a command 57 BACKCOLOR 61 background 64 blowup 60,93 BODIES 7,52 bodies 6,28 body surface 28

body_count 61

boundaries 11,27,29,30,31,38,48,71

BOUNDARY_CURVATURE 68
BOUNDARY 50,51

boundary 27 bounding box 71 bugs 1,98 B graphics command 71 b command 38,57 b graphics command 71 case 7,40 catenoid 11

CHECK_INCREASE 68
check_increase 68

Cheyenne 2,3,5 clear 63

CLIPPED 68

clipped cell 23,30,71

CLOSE_SHOW 66

color map 4

COLORMAP 68

colors 41,59 COLOR 51,61 color 27,28,63 command line options 53 command line 53 command repitition 54 commands, graphics 70 commands 6,54,61 comments 7,40 compiling 4 compressibility 32 compressible 45 conformal metric 49

CONF_EDGE_SQCURV 69

conjugate gradient 35,60,91

CONJ_GRAD 69
CONNECTED 69

connected bodies 23,30,71 conserved quantities 28 consistency checks 57 constant expressions 42 constraint energy integral 79 constraint energy integrand 32 constraint energy 9,30 constraint volume integral 86 constraint volume integrand 14,19,31,35

CONSTRAINTS 51

constraints 27,29,30,31,36,38,47,71

CONSTRAINT_TOLERANCE 36,48
CONSTRAINT 47,50,51

constraint 88 contact angle 9,13,18,32

CONTENT 47

content 35 CONVEX 14,19,30,32,47,48,80 coordinates 27 crystalline integrands 31 crystalline integrand 24,44,79 cube 6 curvature test 61,92 C command 57 C graphics command 71 c command 57 c graphics command 71 datafile 6,40,53

DEBUG 69
DELETE 63

DENSITY 31,52

density 28,62

DETURCK 69
DIFFUSION 69
DIFFUSION 35,46

diffusion 35,38,57 dihedral angle 62 dimension 28,29

DIRICHLET 66

displayable 28,38,52 display 57,60

DISSOLVE 63

DOS 2

DO 55
DUMP 66

dump 57 D command 57 d command 57 d graphics command 8,71 edge popping 59,93 edge wraparound 29

EDGES 7,51

edges 6,27,51

EDGEWEED 66
edge_count 61

effective area 37

EFFECTIVE_AREA 69
EFFECTIVE_AREA 46
ENERGY 47,51

energy 6,31,35,78 environment 3 equiangulation 60,91 errors 53

ESTIMATE 69
EVOLVERPATH 3,53

exit 9,60,71 expressions 42,61 extrapolation 57,92 e command 57

FACES 7,51

faces 51 facet-edges 78 facet-edge 6 facets 6,27

FACET_COLORS 69
facet_count 61

file formats 74 files 3 FIXED 27,31,38,46,48,50,51,64 force 78 FOREACH 63 FRONTCOLOR 61 ftp 2 F command 57 f command 57 gap constant 58 gap energy 30,39,80 gaps 30,32,35,39,45,58,80

                                             GAP_CONSTANT 45

Gaussian curvature 33,70 gaussian quadrature 31,48

                                             GENERIC 4

geomview 2,8,59,75 GLOBAL 30,44,46,47 graphics commands 8,70 graphics interfaces 3,4

                                             GRAVITY 69
                                             GRAVITY_CONST 45

gravity 9,13,18,29,31,32,38,45,57,79 G command 11,38,57 g command 8,57,90 head 27 help 57,71 hidden surfaces 4,71,92

                                             HISTOGRAM 66

history list 54

homogeneous coordinates 4

                                             HOMOTHETY 69

homothety 49,60,93 Hooke energy 33 HP 98731 3 H command 57 h command 57 h graphics command 71 ideal gas model 39,45,52,59,81 ideal gas 32,34 identifiers 40 identifier 40 id 28,61

                                             IF 55

inequality constraint 47

information 57 initiallization 53 installation 2

                                             INTEGRAL_ORDER 31,48
                                             INTERP_NORMALS 69

interrupts 74 interrupt 57 Iris 3 iteration 6,8,35,57,90 i command 57

                                             JIGGLE 46

jiggling 24,36,39,46,58,92 jj graphics command 93 J command 39,58,92 j command 39,58,92 Kelvin tetrakaidecahedron 23 Kelvin 48 keywords 7,41 Klein hyperbolic metric 49

KLEIN_METRIC 29

knot energy 33,45,83

KUSNER 69

K command 58 k command 39,58 length 62 lex 3 line splicing 40

LINEAR 67
lines 7

LIST 63,65

LOAD 66

logging 57

LOGHISTOGRAM 66

long edges 58

LONGJ 66

l command 58 l graphics command 8,71 Macintosh 2 macros 40 macro 9 mail 2 makefile 4 MAXCOORD 29

MAX 62

mean curvature integral 33,45 mean curvature 32,34,36,89 memory 57

MERIT_FACTOR 48
METRIC 49

metric 29,49 Microsoft C 3 MinneView 8,59,74

MIN 62
MOBILITY 48

mobility 36,67 motion 35 M command 28,58 m command 24,36,58 named elements 62 NeXT 2

NODISPLAY 51
NONNEGATIVE 47
NONPOSITVE 47

normal interpolation 59

NORMAL_CURVATURE 69
NORMAL_MOTION 69

notching 59

NOTCH 67

                                             nulgraph 4

numbers 40 N command 59 n command 59 octahedron 24

                                             OLD_AREA 69

omitted faces 35 omitting faces 9 optimizing scale factor 58 orientation 6,7,46,48,51,52,78 original 28,61 o command 12,59 painter algorithm 4,92 PC 2 periodic surfaces 29

                                             PERIODS 44

perturbations 66,92

                                             PHASEFILE 44

phases 44

                                             PHASE 51

pipe 55 Pixar 59,74 PostScript 59,75 Power User 13,18 prescribed mean curvature 32,34 prescribed pressure 81 PRESSURE 34,45,52,81 pressure 28,32,34,36,38,57,60

                                             PRINT 65

P command 59,74 p command 39,59 quadratic representation 28,29,58

                                             QUADRATIC 44,67

quantities 28,29,34,46,48,59,87

                                             QUANTITY 46,47

quantity 51,52 queries 61

                                             QUIET 70
                                             QUIT 67

quit 9,59,60,71 quotient spaces 30 quotient space 44 Q command 34 q command 9,59 q graphics command 8,71

                                             RAWESTV 67
                                             RAWV 67
                                             RAW_CELLS 70
                                             READ 67
                                             REBODY 67
                                             RECALC 67

redirecting input 53

                                             REFINE 63

refining 38,52,60,93 reset 71 Riemannian metric 29 rotate 71 Runge-Kutta 49

RUNGE_KUTTA 70

R graphics command 8,71 r command 9,60 r graphics command 8,71 save 57 scale factor 6,36,48,58 scale optimization 48 scale optimizing 57,90

SCALE_LIMIT 48
SCALE 48
scale 61
SET 64
SHADING 70
SHELL 67
SHOWQ 72
SHOW_EXPR 72
SHOW_INNER 70
SHOW_OUTER 70
SHOW_TRANS 72
SHOW 72

show 60 shrink 71 signals 74 Silicon Graphics 3 simplex model 28

SIMPLEX_REPRESENTATION 44

skinny triangles 58 small facets 60 small triangles 94 soapfilm model 28,43 SoftImage 59,76 sphere example 18 sphere 6

SQGAUSS 70

squared Gaussian curvature 45 squared mean curvature 32,45,81

STABILITY_TEST 67
stability_test 67

stability 37 string model 28,39,43,47,51,52

SUM 62

Sun 2 surface energy integral 44,80 surface energy 13,19,28,31,32 surface tension 28,31,38,44,79

SURFACE_DIMENSION 43
SURFACE_ENERGY 44

surface 6

                                             SYMMETRIC_CONTENT 35,85
                                             SYMMETRY_GROUP 44

symmetry 30 s command 8,60 s graphics command 71

                                             TAG 51

tag 28,61 tail 27 tank example 13

                                             TEMPERATURE 46

temperature 36,39,58

                                             TENSION 31,51

tetrakaidecahedron 23

                                             THICKEN 70

thickening 59

time 49

tiny edges 60,94

toggles 67

topology 6,35,38,59,78 torus duplication 60 torus wrapping 27,29,38,51

                                             TORUS_FILLED 29,43
                                             TORUS_PERIODS 29
                                             TORUS 29,43

torus 23,28,29,59,86,96

                                             total_area 61
                                             total_energy 61
                                             total_length 61
                                             TRANSFORMS 72
                                             TRANSFORM_EXPR 72

translation 71 transparent 63 triangulation 38 triple junction 6 Turbo C 3 t command 12,16,21,60 t graphics command 23,71 unit cell 44 units 6

                                             UNSET 64

user-defined functions 42 U command 60,91 U graphics command 13 u command 12,16,21,60,91 u graphics command 8,71 valence 61 valleys 71 vertex averaging 94 vertex popping 11,59,93

                                             vertex_count 61
                                             VERTEX 50
                                             VERTICES 7

vertices 6,27,50 view transform generators 50 view transforms 50 viewing matrix 49

VIEW_4D 72
VOLCONST 20,52

volume constraints 36 volume constraint 9,13,18,34,35,52,89 volumes 85

VOLUME 8,52

volume 28,35,38,57,59,60,62 v command 8,34,60 v graphics command 71

WHILE 55

whitespace 40 Wulff shape 24 Wulff vectors 4,24,31

WULFF 44

W command 60 w command 16,21,60 w graphics command 71 X-Windows 3 x command 60 x graphics command 71 yacc 3 y 60

ZOOM 67

zoom 61,71,95 Z command 61 z command 61 z graphics command 71