API reference

Every public class and function, from the docstrings. For how the pieces fit together, start with Getting started (the conventions) and Relaxing a surface.

The engine

There is one Surface Evolver engine per process; every Evolver is a handle to it.

class pysurfaceevolver.Evolver(*args, process=False, **kwargs)

A handle to the Surface Evolver engine of this process.

Parameters:

  • datafile – Optional path of a .fe datafile to load right away.

  • echo – If true, print Evolver’s output live while commands run. It’s always returned from command() either way.

  • input – Optional callable input(prompt) -> str | None that answers Evolver’s interactive prompts, for example in hessian_menu. None (the default) answers every prompt with end-of-file.

  • process (bool) – If true, the surface lives in a worker process of its own, with an engine of its own: several surfaces at once (see Notes).

  • threads – With process=True: the worker’s threads (default: as set_threads() would choose; several busy workers share the cores).

  • args (Any)

  • kwargs (Any)

Return type:

Any

Notes

Surface Evolver keeps its state in C globals, so there is one engine per process, and every Evolver object is a handle to it: they all see the same surface, and loading through one changes it for all. echo and input belong to each handle. Use save() and restore() to keep a surface around, and process=True to work on several at once.

Evolver(process=True) starts a worker process (a new Python interpreter) with an engine of its own and returns a handle on it that works like any other: attribute reads, assignments and method calls go to the worker, views (ev.body(1), ev.parameters) are handles too, and results (meshes, IterationResults, arrays, snapshots) come back as copies. Each call costs a round trip (a fraction of a millisecond) plus copying its arguments and result. Callables given to methods (relax(callback=...), undo_if=) are copied to the worker and run there (sending functions defined in a script or notebook needs cloudpickle); input= isn’t available. plot() and live_view() draw in this process (a live view updates when you call update(); a callback running in the worker can’t reach it). Prints and warnings from the worker show here, and Ctrl-C interrupts the worker’s call. ev.close() (or with, or the last handle going away) ends the worker; ev.pid is its process id. If it dies, the call raises WorkerCrashed. POSIX systems only.

Calls are serialized: a call made while another thread is running one raises RuntimeError instead of waiting.

Ctrl-C during a call (from the main thread) stops Evolver at the next iteration or statement; pressing it again aborts the operation. Both raise KeyboardInterrupt.

adapt(max_turn=15.0, *, coarsen=True)

Fit the mesh to the curvature: split the edges along which the surface turns by more than max_turn degrees (the angle between the vertex normals at its ends), such as on the strongly curved rim of a drop at a high contact angle, then equiangulate. One level per call.

With coarsen (default), first merge the ends of short interior edges (under 3/4 of the bulk edge length, the 75th percentile) along which it turns by less than a third of that, where the surface has flattened since it was refined. A merge is refused if it would tilt a facet by more than max_turn/2, flip one, or make an edge longer than 4/3 of the bulk length (checked in Evolver’s edge deletion; collapse_max_tilt and collapse_max_edge). Boundary edges (contact lines, wires) are never merged.

Edges flagged no_refine are kept whole, and edges touching a junction of three or more facets (a triple line) are not judged (no normal there). Relax afterwards; relax(adapt=True) does both.

Linear model only. Returns the number of edges split plus merged.

Parameters:
Return type:

int

bodies()

Return volumes, target volumes and pressures of all bodies (a snapshot).

Return type:

Bodies

body(number)

A live handle on one body (1-based, as in the datafile).

ev.body(1).target = 0.05 sets the prescribed volume, and ev.body(1).pressure reads the current pressure. For all bodies at once, as arrays, use bodies().

Parameters:

number (int)

Return type:

BodyView

body_surfaces(*, curved='tessellate', n=None, cap=False)

The closed surface around each body, with outward normals.

See Mesh.body_surfaces(). With cap=True, an opening whose boundary vertices all lie on a constraint (a drop on a curved solid, say) is closed by a cap on that constraint: Evolver projects its points, so any constraint formula works (3D models; elsewhere caps are flat). cap_constraints says which constraint each cap is on (None: not on one; a flat cap).

Projecting creates and dissolves temporary vertices: the surface is left as it was, but Evolver counts it as changed (the next Newton step rebuilds its matrix pattern, and later new elements may get other numbers than without the export).

Parameters:
Return type:

Dict[int, BodySurface]

check()

Evolver’s consistency check of the surface’s topology (check): the problems found, one per line; empty when it is sound.

Return type:

List[str]

command(text)

Run Evolver commands, such as "g 10; r; hessian".

Returns everything Evolver printed to standard output. Raises EvolverError if Evolver reports an error.

Parameters:

text (str)

Return type:

str

contact_angles(constraint, body=None)

The contact angle (degrees, through the liquid) at each vertex where the surface meets a constraint (number or name): the angle between the surface’s normal out of body (area-weighted over its flat facets) and the constraint’s normal, the gradient of its formula. Vertices also on another constraint (corners) are left out.

body: the liquid; needed when the surface has several bodies (without bodies, facet normals are taken as they are). The gradient is taken to point out of the solid, as with the walls of pysurfaceevolver.constraints (a plane’s normal into the liquid; spheres and cylinders solid); otherwise the angles come out as 180 degrees minus the contact angle.

Parameters:
Return type:

ndarray

contact_lines(constraint)

Where the surface meets a constraint (number or name): the edges on it that bound the surface (on one facet only; for strings, the vertices’ chain on it), joined into polylines of vertex positions, longest first; a closed line repeats its first point. Curved models: through the corner vertices. For a drop’s contact radius, say, np.linalg.norm(line[:-1, :2] - line[:-1, :2].mean(axis=0), axis=1).mean().

Parameters:

constraint (int | str)

Return type:

List[ndarray]

property counts: dict

Number of vertices, edges, facets and bodies.

property datafile: str

Name of the loaded datafile, or "" if none.

define(name, value)

Make a new variable (or set one), for use in commands and expressions: ev.define("h", 0.1), then ev.command("t h").

Parameters:
Return type:

None

define_attribute(element, name, dtype='real')

Define an extra per-element attribute (real or integer).

Read it with values() and write it with set_values().

Parameters:
Return type:

None

derivative(name, step=None, *, move=None, method='newton')

The derivative of the equilibrium energy (energy) with respect to a parameter, by central differences: the parameter set to its value plus and minus step (default 1e-5 times its size, at least 1e-5; 1e-3 with method="relax"), the energy found each time, the surface put back as it was. With a parameter that moves a wall or a pad, minus this is the force on it (Brakke’s bga-4 tutorial compares seven ways).

method:

  • "newton" (default): safeguarded Newton steps reconverge the surface each time, as Brakke’s zforce5. Newton moves vertices along the normal only, so the mesh doesn’t drift between the two and the difference is clean even for tiny steps.

  • "fixed": no reconverging, the energy corrected for the volume and quantity changes by their multipliers (minus pressure times change), as Brakke’s preferred zforce4: as accurate, and far cheaper on a large surface, but it needs move.

  • "relax": relax() each time; it tidies the mesh too, which changes the energy by more than a small step does, hence its larger default step.

move(xyz): how the vertices follow a unit change of the parameter, an (n, sdim) array (a pad at height h: lambda x: np.c_[0*x[:, :2], x[:, 2]/h], everything stretched with it); the surface is moved that way before the energy is taken, which also gives Newton a close start. An optimizing parameter is held fixed meanwhile: at its optimum the derivative is zero (with a load in the energy, the joint’s force balances it). (A parameter that only enters a quantity’s modulus, which Evolver evaluates once, gives 0.)

Parameters:
Return type:

float

dump(path)

Save the current surface as a datafile (Evolver’s dump).

Numbers are written with 17 significant digits, so loading the file gives back the same coordinates exactly.

Parameters:

path (str | PathLike)

Return type:

None

eigen_counts(shift=0.0)

How many Hessian eigenvalues lie below, at and above shift (Evolver’s eigenprobe; costs one factorization).

At an equilibrium, negative > 0 means it is unstable (a saddle): some motion lowers the energy. With constraints, only motions that keep them count.

Parameters:

shift (float)

Return type:

EigenCounts

energies()

What the energy is made of, term by term: {name: energy} for every energy method instance once Evolver has converted everything to named quantities (default_area_inst for the facets’ tension, constraint_1_energy_inst for a wall’s line integral, gravity_inst, …). Brakke’s way of checking a datafile (convert_to_quantities; show_all_quantities; Q): a term with the wrong sign or size points at the mistake. The terms add up to total_energy. Works on a copy: the surface is left as it was.

Return type:

Dict[str, float]

property energy: float

The total energy including what the walls’ line integrals leave out (pysurfaceevolver.constraints.Constraint.energy_constant, added by make_datafile(); Evolver can’t add a constant): total_energy plus that constant.

equiangulate()

Flip edges to improve triangle shapes (Evolver’s u).

Return type:

None

eval(expr)

Evaluate a numeric Evolver expression, like "body[1].volume".

Nothing gets defined in Evolver’s symbol table. Results are exact doubles, including infinities and NaN.

Parameters:

expr (str)

Return type:

float

fix(element, where=None)

Fix vertices, edges or facets (all, or those selected by where).

Parameters:

element (str)

Return type:

None

property gravity: float

no gravity). Bodies weigh their density times it.

Type:

The gravitational constant (Evolver’s gravity_constant; 0

health(tol=1e-08, stability=False)

A check-up of the surface: is it in equilibrium (residual() below tol), stable (with stability=True, stability()), are its facets in shape, does it come close to (or through) a wall or mirror it isn’t attached to, or to itself?

The gaps are in units of the median edge length. wall_gap leaves out the vertices within two edges of a constraint’s own vertices (next to a contact line the surface is close to its wall by design); self_gap measures between vertices more than two edges apart. junctions counts free junctions that break Plateau’s laws (edges where four or more films meet, vertices whose tangent cone Evolver’s o pops, such as more than four triple edges, or, in a string network, four or more strings meeting): the surface is in equilibrium for its topology, but popping them (Evolver’s o) lowers the energy. target_gap is the largest relative gap between a fixed body volume or quantity and its target (Evolver’s projection onto them can fail on a degenerate mesh; an issue above 1e-4). crossed counts vertices on the far side of a constraint: the forbidden side of a one-sided one, or a few vertices (at most 5%) against all the others for an equality constraint. issues puts what looks wrong in words: unconverged, unstable, facets below 1 degree, crossings, gaps under 0.1, targets missed. (For a crystalline energy the residual is left out of issues: it doesn’t vanish at the Wulff shape.) relax() attaches one to its result. Cost: about a gradient step (the stability check: two factorizations).

Parameters:
Return type:

Health

iterate(n=1, *, callback=None, every=1)

Run n gradient-descent iterations (Evolver’s g).

Returns the energy, area and scale factor after each iteration. callback(ev, i) is called after every every-th iteration and after the last one, for example to update a LiveView. Ctrl-C stops between iterations with KeyboardInterrupt.

Parameters:

  • n (int)

  • callback (Optional[Callable[['Evolver', int], Any]])

  • every (int)

Return type:

IterationResult

property lagrange_order: int
live_view(scalars=None, **kwargs)

A PyVista view that redraws on LiveView.update(); see LiveView.

Typical use: ev.iterate(200, callback=ev.live_view().update, every=10).

Parameters:

  • scalars (Union[None, str, np.ndarray])

  • kwargs (Any)

Return type:

LiveView

load(datafile)

Load a .fe datafile, replacing the current surface.

Evolver looks for the file in the current directory, then in the directories listed in the EVOLVERPATH environment variable. Returns the output printed while loading.

If the datafile has errors, this raises EvolverError, and there is no valid surface until a datafile loads successfully.

Parameters:

datafile (str | PathLike)

Return type:

str

load_arrays(vertices, faces=None, **kwargs)

Load a surface given as arrays; see make_datafile() for options. The walls among constraints (from pysurfaceevolver.constraints) are remembered for plot(walls=True).

Parameters:

kwargs (Any)

Return type:

str

load_mesh_file(path, *, volume=None, fix_boundary=True, merge_tolerance=1e-09, **kwargs)

Load a triangle surface from a mesh file (STL, OBJ, PLY, VTK, …).

Duplicate points (closer than merge_tolerance times the bounding box size) are merged, and quads are split into triangles.

volume:

None: no body. A number: one body enclosed by the whole surface, with that target volume. "current": the same, with the volume the surface encloses now. The surface must be closed; its orientation is fixed up so the volume is positive.

fix_boundary:

Fix the vertices and edges on open boundaries, so a film spanning a wire keeps its wire.

Other keyword arguments go to make_datafile().

Parameters:
Return type:

str

load_string(text, name='surface.fe')

Load a datafile given as a string.

The text goes to a temporary file first, so #include and read paths are resolved from the current directory.

Parameters:
Return type:

str

mesh()

Return the current geometry and connectivity as NumPy arrays.

Faces exist only in the soapfilm representation. For quadratic and Lagrange models the high-order node layout is included too; see Mesh.tessellate().

The result is cached until the surface changes (any command, load or coordinate write), so its arrays are read-only; copy them to modify.

Return type:

Mesh

mesh_quality(skinny_angle=15.0)

Edge lengths and facet shapes; see Mesh.quality(). In a torus or quotient space, from Evolver’s own edge lengths (vertex coordinates are wrapped, so facets across the cell’s sides would look torn).

Parameters:

skinny_angle (float)

Return type:

MeshQuality

property model: str

"linear", "quadratic" or "lagrange".

modes(n=1)

The n lowest Hessian eigenvalues and their modes: the motions of the vertices, with the constraints (volumes, quantities, walls) held.

values are ascending, with Evolver’s linear_metric as in stability() (the continuous surface’s eigenvalues: a negative one is a motion that lowers the energy). vectors[k] is mode k as a displacement per vertex (0 for fixed vertices), scaled to a root mean square of 1 over the vertices and signed so that its largest component is positive (among components within 1% of the largest, the first in vertex order), so a simple eigenvalue’s mode comes out the same each time; a repeated eigenvalue’s (a symmetry’s rotations or translations) is any combination of its eigenspace. ev.vertices + 0.01*m.vectors[0] steps along the first. Unlike saddle(), the surface is left as it is. Fewer than n come back when the surface has fewer degrees of freedom. A mode that moves optimizing parameters moves them too; vectors has only the vertices. Uses Evolver’s ritz, slow on large surfaces.

Parameters:

n (int)

Return type:

Modes

newton(steps=1, *, seek=None, tol=None, undo_if=None, threads=None)

Up to steps safeguarded Newton steps.

Each step is a plain Newton step (hessian), or, if that one is rejected, one with a line search (hessian_seek); seek=True or False uses only the one kind. A step is rejected (undone) if it makes the energy or the residual() non-finite, raises both the energy and the residual, raises the energy by more than 10% plus what meeting the volume and quantity targets costs (pressure times the change), or makes undo_if(ev) true; a step that brings an unmet target (say a new body volume) at least halfway may raise both. With no step accepted, Newton stops. It also stops when neither the residual nor an unmet target is improving, or, with tol, once the residual is below tol with the targets met.

Newton converges to the equilibrium nearest in its own sense, saddles included: after a large jump of a target (a volume made 16 times larger in one go) it can end at an unstable one. Check with stability(), or let relax() (gradient steps first) or a continuation in smaller steps take such jumps. Returns an IterationResult: newton_steps (the steps kept), residual, and converged: the residual below tol (without one, relax()’s default 1e-8) and the targets met. Stability isn’t checked (stability()). (hessian() takes one step with no checks.)

Parameters:

  • steps (int)

  • seek (Optional[bool])

  • tol (Optional[float])

  • undo_if (Optional[Callable[['Evolver'], bool]])

  • threads (Optional[int])

Return type:

IterationResult

on_constraint(constraint, element='vertex')

Which elements are on a constraint (number or name): one bool per mesh() row of that element type.

Parameters:
Return type:

ndarray

property parameters: Parameters

The datafile’s parameters, as a live dict-like view.

ev.parameters["angle"] = 60 assigns, and Evolver recalculates whatever depends on it.

plot(scalars=None, *, element=None, n=None, off_screen=False, screenshot=None, mirror=None, walls=None, **kwargs)

Show the surface with PyVista.

scalars colors it: an Evolver expression evaluated per facet (or per vertex with element="vertex"), or an array with one value per facet or vertex row. Curved elements are tessellated with n subdivisions. mirror adds mirror images, for a symmetric piece of a surface: a list of planes, each "x" (x = 0), ("x", c) (x = c) or (normal, point), applied in order, each doubling what is there (an eighth of a cell: ["z", "y", "x"]). walls draws solids and walls with it: constraints from pysurfaceevolver.constraints (C.sphere(...), C.plane(...); planes as translucent patches around the surface) or PyVista meshes; walls=True draws those given to load_arrays(). Other keyword arguments go to Plotter.add_mesh.

Parameters:
Return type:

Any

pop()

Pop the junctions a minimal surface can’t have (Evolver’s o: four films on an edge, a cube’s cone at a vertex, four strings at a point, …), so they can split into Plateau’s junctions as the surface relaxes. Returns how many were popped. relax() does this by default in the linear model.

Return type:

int

quantities()

Named quantities from the datafile, recalculated now.

Return type:

Dict[str, Quantity]

quantity(name)

A live handle on a named quantity from the datafile.

ev.quantity("stretch").modulus = 100 changes a modulus, ev.quantity("vol").target = 2 sets (and fixes) a target and = None frees it. For all quantities at once (a snapshot), use quantities().

Parameters:

name (str)

Return type:

QuantityView

recalc()

Recalculate energies and quantities.

Return type:

None

refine(times=1)

Refine the surface: split every edge and facet (Evolver’s r).

Parameters:

times (int)

Return type:

None

relax(tol=1e-08, max_iter=1000, *, levels=0, newton=None, stability=None, adapt=False, tidy=10, seek=None, undo_if=None, callback=None, every=1, threads=None)

Relax the surface to an equilibrium: gradient steps that keep the mesh in shape, then safeguarded Newton steps, until residual() is below tol. The defaults are meant to work as they are; the result says whether it converged, whether the equilibrium is stable, and how the surface looks (health).

Common options: levels=n refines and relaxes again n times; newton=0 uses gradient steps only; stability=False skips the stability check (about two factorizations); adapt=True (linear model) then refines where the surface curves too much for its edges and coarsens where it has flattened (adapt(); or give the turning angle in degrees, default 15), up to three passes or four times the facets the call began with. callback(ev, i) runs every every-th gradient step and at the end of each level; threads sets the threads for this call only (see pysurfaceevolver.threads_limit()).

How it works:

  • With levels=0, a surface already in equilibrium (residual and the gaps to its volume and quantity targets below tol, no degenerate facet) is left as it is: only the stability check below (and adapt, if given) runs. Gradient steps on it would be mesh motion only, which next to a shallow contact line slides vertices through the wall.

  • In the quadratic and Lagrange models (a final stage from a mesh settled in the linear model) with levels=0, five gradient steps and then safeguarded Newton steps come first (as the manual’s g 5; hessian), Newton as long as each step cuts the residual by 10%; if that converges, that is all. If not, the surface goes back and the steps below follow; right after such a miss, the next call skips the attempt.

  • A body volume or fixed quantity whose target is more than 25% from its value: if the first gradient step toward it leaves degenerate facets (below 1 degree) the surface didn’t have, it is undone and the target is approached in stages of at most 10%, each followed by 20 gradient steps (tidied as below): in one go, Evolver’s projection onto a far target can fold facets.

  • Gradient steps (Evolver’s g) come in rounds of tidy steps, each followed by equiangulation and vertex averaging (u; V; only V for strings and the Lagrange model, only u in the simplex representation; tidy=0: none), until the energy changes by less than 1e-5 (relative) over three rounds when Newton follows (it finishes faster and more reliably, and long gradient runs let slow mesh drift grow), 1e-9 with newton=0, or max_iter steps. When Newton follows, a round that more than doubles the residual from its lowest so far is undone and Newton takes over (on a coarse mesh a contact line can collapse within a few steps). Rounds that have used half of max_iter and keep slowing without settling (the residual falling by less than 5% a round, three rounds running: a contact line creeping over a curved wall) get a trial of Newton steps, kept if they converge or cut the residual a hundredfold; otherwise the surface goes back and the rounds go on (a neck shrinking toward a pinch-off creeps the same way).

  • Then up to newton safeguarded Newton steps (see newton(); seek and undo_if as there), until the residual is below tol or stops falling.

  • If that leaves the surface unconverged, one more try with conjugate gradients (far faster on flat or stiff shapes, but they wear the mesh down where contact lines travel), then Newton, with up to 200 of the steps the level left of max_iter; kept only if that converges. Otherwise plain gradient steps (your own conj_grad setting is restored after).

  • With levels=n, if the last level still ends unconverged: the levels again from the call’s start refined once first, kept only if that converges. A coarse start can hold equilibria the shape doesn’t have, and which one its gradient rounds find can turn on the last digit of a volume (a bridge at 20 degrees from a tube 12 facets around).

  • A surface that ends unconverged with vertices sunk into a wall from pysurfaceevolver.constraints whose contact angle is below 15 degrees: the levels again from the call’s start with the small angles raised to 20 degrees, then lowered in stages (by 0.7 each), relaxed at each, kept only if that converges at their own angles. At a small angle the liquid meets the solid almost tangentially, gradient steps from a crude start slide the vertices next to the contact line through it, and the wall’s integrals reward that (a trimer between three beads at 5 degrees). Walls with such angles carry a datafile parameter pyse_angle_floor (0) for this.

  • Once converged (linear model), edges leaving a contact line that are longer than 1.5 times the contact-line edges there are split and the surface relaxed again, kept only if that converges: a shrinking contact circle (a drop at a high contact angle) would otherwise end up smaller than the facets next to it.

  • Then (linear model, films and strings) junctions that break Plateau’s laws are popped (Evolver’s o: four films on an edge, four strings at a point, …) and the surface relaxed again, up to five times. If it is still unconverged, a last resort: up to 300 gradient steps without tidying (which can hold a shrinking neck open), then edges that have collapsed below a 20th of the median edge are deleted (Evolver’s t: a film or neck shrinking to nothing, as before a foam’s T1 change or a pinch-off), junctions popped and the surface relaxed, kept only if that converges, or lowers the energy on a sound mesh without raising the residual. A converged surface never loses edges. After k failed attempts in a row, the next 2**k - 1 calls skip the last resort. popped and collapsed count them.

  • Once converged (or within 100 times tol), with stability, stability() sets stable and negative_modes and warns (UnstableEquilibriumWarning) at a saddle, such as a liquid column past its Rayleigh-Plateau limit: the solver finds saddles too, so the warning is the only sign.

  • Without a Hessian (gap energies, simplex bodies outside the Lagrange model, …) the Newton steps give way to gradient rounds settled as with newton=0, while no facet is degenerate (below 1 degree); a crystalline (Wulff) energy takes gradient rounds only and counts as converged once its energy settles, since its residual doesn’t vanish at the Wulff shape. settled says the energy settled with the residual above tol.

If the energy becomes non-finite, the surface is put back to the start of the gradient round where it happened and EvolverError is raised. A gradient round that pushes a met body volume or fixed quantity off its target (Evolver’s projection can fail on a degenerate mesh) is undone and ends the rounds.

Returns an IterationResult: one entry per gradient step in energy, area, scale and level (adaptive passes count as further levels), plus newton_steps, residual, converged, settled, stable, negative_modes, health, popped and collapsed.

Parameters:

  • tol (float)

  • max_iter (int)

  • levels (int)

  • newton (Optional[int])

  • stability (Optional[bool])

  • adapt (Union[bool, float])

  • tidy (int)

  • seek (Optional[bool])

  • undo_if (Optional[Callable[['Evolver'], bool]])

  • callback (Optional[Callable[['Evolver', int], Any]])

  • every (int)

  • threads (Optional[int])

Return type:

IterationResult

remesh(target=None, *, max_edge=None, min_edge=None, equiangulate=True, average=False, protect=None, keep_boundary=True)

Even out edge lengths: delete edges shorter than min_edge (Evolver’s t), split edges longer than max_edge (refine edge where length > max_edge), then equiangulate (u) and, with average=True, average the vertices (V). Deleting first: merging vertices lengthens the edges around them. One call splits each long edge once and doesn’t delete every short edge (Evolver refuses unsafe merges); call it every few steps rather than expecting one call to finish the job.

target=h sets max_edge=1.6*h and min_edge=0.5*h (explicit values win). Call it every few steps when the surface stretches or shrinks a lot, for example while a contact line moves.

protect (an edge mask, aligned with mesh() rows) keeps edges from being split, by setting no_refine on them for the split; their own flags are restored after. Edges flagged no_refine already are never split (unlike Evolver’s l, which ignores the flag). Edges leading from a contact line on a curved constraint into the surface are good candidates: their midpoint would land off the constraint, e.g. inside a solid sphere (mesh.edges_touching(ev.on_constraint(k), "one")).

keep_boundary (default): edges on the surface’s boundary (contact lines, wires) are never deleted: they are short where the geometry is small, and deleting them would coarsen it.

Linear model only (Evolver deletes edges only there). Returns the numbers of edges deleted (about: each deletion removes three), split and switched.

Parameters:
Return type:

Dict[str, int]

property representation: str

"string", "soapfilm" or "simplex".

residual()

How far the surface is from equilibrium, as a dimensionless number (0 at an equilibrium; nothing moves).

From the vertex velocities a gradient step would use (forces projected on the constraints, with the volume and quantity multipliers applied), the part that changes the shape: along each vertex’s normal, projected into its constraints (at a contact line: the imbalance of the contact angle). Tangential parts, and freedoms that only slide a vertex along a wire, are mesh motion: Newton steps leave them to the gradient steps and so does the residual. Each velocity is about the vertex’s share of the area times the local pressure imbalance, so the root mean square times N/sqrt(area) measures that imbalance relative to the surface’s size: comparable across meshes and scales. With an elastic energy (Evolver’s elastic methods: in-plane forces are physical there) all directions count, not only the normal. Optimizing parameters count too: the force on each (minus the derivative of the energy, with the volume and quantity multipliers applied), as if spread over the vertices. (For a string: the root mean square times N, the imbalance times the length; in general the imbalance times the surface measure to the power 1/dimension.)

Return type:

float

restore(snapshot)

Replace the surface with a snapshot from save().

Parameters:

snapshot (Snapshot)

Return type:

None

saddle()

Step off a saddle point (Evolver’s saddle): along the eigenvector of the lowest Hessian eigenvalue, if it is negative, as far as the energy falls. Safe far from equilibrium too (Brakke’s eigenvalue tutorial). Returns how much the energy fell, 0 at a local minimum; follow with relax() to settle into the minimum it leads to. (Evolver needs pressures from a gradient step first: one is taken if none was.)

Return type:

float

save()

Take a snapshot of the surface, to restore() later.

The snapshot is an exact dump: coordinates and energies come back bit for bit. Settings that dump doesn’t record (most display options) are not part of it.

Return type:

Snapshot

property sdim: int

Dimension of the ambient space.

set_constraint(element, constraint, where=None, on=True)

Put elements on a constraint (or take them off with on=False).

Parameters:
Return type:

None

set_flag(element, flag, where=None, on=True)

Set (or with on=False clear) a yes/no attribute such as fixed or no_refine on all elements of a type, or those selected by where (a boolean mask or rows, aligned with mesh() rows).

Prefer this over Evolver’s foreach ... where ... in commands: there, a bare element type inside an aggregate (max(vertex, ...)) means all elements of the surface, and vertex[2] means global vertex 2, not an edge’s second vertex (that is ee.vertex[2] with a named element).

Parameters:
Return type:

None

set_model(model, order=None)

Switch to the "linear", "quadratic" or "lagrange" model.

Lagrange needs an order.

Parameters:

  • model (str)

  • order (int | None)

Return type:

None

set_values(element, attribute, values, *, where=None)

Set a numeric attribute on every element of a type.

values is one number per element, aligned with mesh() / values() rows, or a single number for all. where (a boolean mask or rows) limits which elements change. Works for any attribute Evolver’s set command accepts: coordinates (x, y, …), density, tension, target (bodies), extra attributes from define_attribute(), …

Coordinates and scalar extra attributes are written directly in C; other attributes through one Evolver set command (the values and selection go in scratch attributes, pyse_<type>_value and pyse_<type>_mark).

Parameters:

  • element (str)

  • attribute (str)

Return type:

None

stability(nearest=6, *, eigenvalues=False)

Whether the surface sits at a stable equilibrium (a minimum, as far as second order tells; call it once converged, see residual()).

Counts the Hessian eigenvalues below zero (one factorization); if there are any, counts again below a small negative threshold: -1% of the scale of the lowest ones (the median magnitude of the positive ones among the nearest to zero, from Evolver’s ritz, which is slow on large surfaces). Exact zero modes from symmetries (a drop sliding on a plane, a barrel along its fibre) come out slightly above or below zero numerically; the threshold keeps them from reading as instabilities, at the price of noticing a real one slightly past its onset.

lowest and threshold are eigenvalues with Evolver’s linear_metric (as its eigenvalue tutorial recommends): those of the continuous surface, which don’t shrink as the mesh is refined (the counts are the same either way). On a stable surface lowest is empty unless eigenvalues=True (then ritz runs anyway: the lowest eigenvalue falling toward zero over a continuation shows a fold or instability coming).

Parameters:
Return type:

Stability

property torus: bool

Whether the domain is periodic (torus model).

property total_area: float

Total area as of the last iteration or recalculation.

property total_energy: float

Total energy as of the last iteration or recalculation.

unfix(element, where=None)

Unfix vertices, edges or facets (all, or those selected by where).

Parameters:

element (str)

Return type:

None

property valid: bool

Whether there is a usable surface (see InvalidSurfaceError).

values(element, expr)

Evaluate an expression for every element of one type.

element is "vertex", "edge", "facet" or "body". Inside expr, attributes refer to the element at hand, for example ev.values("vertex", "x^2 + y^2") or ev.values("facet", "area").

The result is aligned with the rows of mesh() (vertices, edges, facets) and bodies().

Parameters:
Return type:

ndarray

vertex_average()

Move vertices to the average of their neighbors (Evolver’s V).

Return type:

None

property vertices: ndarray

Vertex coordinates as an (n, sdim) array (a copy).

Assigning an array of the same shape moves the vertices and recalculates energies. Constraints aren’t re-projected until the next iteration.

write(path, *, curved='tessellate', n=None, file_format=None)

Write the surface to a mesh file with meshio (STL, OBJ, PLY, VTU, Gmsh .msh, XDMF, …); the format follows the file extension.

curved="tessellate" writes flat triangles sampled on curved elements (n per facet edge); curved="native" writes the high-order elements themselves: Gmsh up to order 10, VTU/VTK any order, XDMF order 2; STL, OBJ, PLY etc. only take flat triangles. XDMF needs h5py. .msh means Gmsh format 2.2 (meshio alone would pick ANSYS, and its Gmsh 4.1 output can’t be opened by Gmsh). See Mesh.to_meshio() for cell data.

Scenes for viewers and 3D printing software go through PyVista: .gltf/.glb and .wrl (VRML; flat triangles, no cell data).

Parameters:
Return type:

None

write_bodies(pattern='body_{id}.stl', *, curved='tessellate', n=None, cap=False, require_watertight=True, file_format=None)

Write one closed surface per body, for volume meshing.

pattern is formatted with the body id. Raises ValueError for a body whose surface isn’t watertight (for example a drop that Evolver closes with a constraint plane) unless cap=True closes it or require_watertight=False. Returns {body id: path}.

Parameters:

  • pattern (str)

  • curved (str)

  • n (int | None)

  • cap (bool)

  • require_watertight (bool)

  • file_format (str | None)

Return type:

Dict[int, str]

class pysurfaceevolver.IterationResult(energy, area, scale, output, converged=None, newton_steps=0, level=None, residual=nan, stable=None, negative_modes=None, health=None, popped=0, collapsed=0, settled=None)

What Evolver.iterate(), Evolver.relax() or Evolver.newton() did.

energy, area, scale and level (the refinement level, from 0) have one entry per gradient iteration (none for newton()). converged, newton_steps and residual are set by Evolver.relax() and Evolver.newton().

Parameters:
energy: ndarray
area: ndarray
scale: ndarray
output: str
converged: bool | None = None
newton_steps: int = 0
level: ndarray | None = None
residual: float = nan
stable: bool | None = None
negative_modes: int | None = None
health: Health | None = None
popped: int = 0
collapsed: int = 0
settled: bool | None = None
class pysurfaceevolver.Health(residual, stable, negative_modes, angle_min, skinny, wall_gap, wall, crossed, self_gap, junctions, target_gap, issues)

See Evolver.health(). Gaps are in units of the median edge length.

Parameters:
residual: float
stable: bool | None
negative_modes: int | None
angle_min: float
skinny: int
wall_gap: float
wall: int | None
crossed: Dict[int, int]
self_gap: float
junctions: int
target_gap: float
issues: List[str]
property ok: bool
summary()

"ok", or the issues in a line.

Return type:

str

class pysurfaceevolver.BodyView(ev, number)

A live handle on one body; see Evolver.body().

Reading an attribute asks Evolver now; setting one changes the surface.

Parameters:
property volume: float

The current volume.

property target: float | None

The prescribed volume, or None if the volume is free. Setting a number fixes the volume; setting None frees it.

property volconst: float

A constant added to the computed volume (for parts of the body the surface doesn’t enclose, such as a solid inside it).

property pressure: float

the Lagrange multiplier of the volume constraint.

Type:

The pressure

property density: float

The density, for gravity (see Evolver.gravity).

property fixed: bool

Whether the volume is prescribed.

class pysurfaceevolver.QuantityView(ev, name)

A live handle on one named quantity; see Evolver.quantity().

Reading an attribute asks Evolver now; setting one changes the surface.

Parameters:
property value: float

The current value.

property kind: str

"energy", "fixed", "info" or "conserved".

property target: float | None

The prescribed value, or None unless the quantity is fixed. Setting a number fixes it at that value; setting None frees it (it becomes an info quantity, as Evolver’s unfix).

property modulus: float

The factor its value enters the energy with (energy quantities).

property pressure: float

The Lagrange multiplier of a fixed quantity.

class pysurfaceevolver.EigenCounts(negative, zero, positive)

Hessian eigenvalues below, at and above a shift; see Evolver.eigen_counts().

Parameters:
class pysurfaceevolver.Stability(stable, negative, lowest, threshold)

See Evolver.stability().

Parameters:
stable: bool
negative: int
lowest: ndarray
threshold: float
class pysurfaceevolver.Modes(values, vectors)

See Evolver.modes().

Parameters:
values: ndarray
vectors: ndarray
class pysurfaceevolver.Snapshot(text, datafile, walls=<factory>, generation=-1)

A saved surface from Evolver.save(), as exact datafile text.

Parameters:
text: str
datafile: str
walls: Dict[object, object]
generation: int = -1
write(path)

Write the snapshot as a .fe datafile.

Parameters:

path (str | PathLike)

Return type:

None

class pysurfaceevolver.Parameters(evolver)

Live view of the parameters declared in the loaded datafile.

Only existing parameters can be assigned, and none can be deleted.

Parameters:

evolver (Evolver)

property optimizing: frozenset

ones the energy minimization moves too, like vertex coordinates (Evolver’s optimizing_parameter; see optimize()).

Type:

Names of the optimizing parameters

optimize(name, on=True)

Make a parameter optimizing (Evolver’s unfix): Evolver.relax() then finds its value too, say the height a loaded joint settles at. on=False fixes it again at its current value.

Parameters:
Return type:

None

Surface data

class pysurfaceevolver.Mesh(vertices, edges, facets, vertex_ids, edge_ids, facet_ids, facet_bodies, fixed, order=1, bezier=False, edge_nodes=None, edge_node_index=None, facet_nodes=None, facet_node_index=None)

A snapshot of the surface geometry.

edges and facets hold row indices into vertices. The *_ids arrays hold Evolver’s own 1-based element numbers, as used in commands such as vertex[5].x.

In quadratic and Lagrange models, vertices also holds the extra nodes, facets and edges use the corner vertices only, and facet_nodes and edge_nodes list every node of each element. Use tessellate() for flat triangles that follow the curved facets.

Parameters:
vertices: ndarray
edges: ndarray
facets: ndarray | None
vertex_ids: ndarray
edge_ids: ndarray
facet_ids: ndarray | None
facet_bodies: ndarray | None
fixed: ndarray
order: int = 1
bezier: bool = False
edge_nodes: ndarray | None = None
edge_node_index: ndarray | None = None
facet_nodes: ndarray | None = None
facet_node_index: ndarray | None = None
edge_lengths()

The length of each edge (straight, between its end vertices), aligned with edges. Read-only (kept with the mesh).

Return type:

ndarray

edge_valence()

The number of facets on each edge (all 0 without facets), aligned with edges: 1 on a boundary or contact line, 3 or more at a junction of films. Read-only.

Return type:

ndarray

facet_areas()

The area of each facet, flat between its corner vertices (see tessellate() for curved ones), aligned with facets. Read-only.

Return type:

ndarray

mirrored(planes)

The mesh with its mirror images (see mirror_matrices() for planes), as one mesh: for showing or exporting a full cell from a symmetric piece. Facets of odd images are reversed, so normals keep pointing the same way relative to the surface. Element ids repeat (one copy per image). Linear meshes only; plots mirror curved ones too (ev.plot(mirror=...)).

Return type:

Mesh

quality(skinny_angle=15.0)

Edge lengths and facet shapes, from the corner vertices.

skinny counts facets with an angle below skinny_angle degrees, degenerate those whose area is below 1e-12 of the median.

Parameters:

skinny_angle (float)

Return type:

MeshQuality

edges_touching(vertex_mask, how='any')

Edges by their end vertices, as a mask over edges rows.

vertex_mask is one bool per vertex row (for example from Evolver.on_constraint()). how: "any" (at least one end selected), "all" (both ends) or "one" (exactly one end: the edges leading off a contact line into the surface).

Parameters:

how (str)

Return type:

ndarray

slice(normal=(0.0, 0.0, 1.0), point=(0.0, 0.0, 0.0), *, n=None)

Where the surface crosses a plane: a list of polylines, (k, 3) arrays of points in order along the curve, longest first (a closed curve repeats its first point at the end). Curves end where the surface does, and at junctions of several films.

The plane goes through point with the given normal; for example mesh.slice((0, 1, 0)) is the section y = 0. Curved (quadratic, Lagrange) facets are cut through tessellate() (with its n). Vertices lying in the plane belong to the curve, so a mirror plane’s own vertices come out as its curve.

Parameters:

n (int | None)

Return type:

List[ndarray]

tessellate(n=None, *, merge=True, values=None)

Subdivide every facet into n*n flat triangles, for plotting or export.

Points are sampled on the curved (quadratic or Lagrange) facets, so the triangles follow the actual surface. n defaults to 1 for linear elements and 2*order otherwise. Triangle t belongs to facet t // n**2, and is oriented like that facet’s facets row. More than pyse.tessellation_limit triangles (default 10 million) gives a LargeTessellationWarning.

With merge=True (the default), points shared by neighboring facets appear once, so the result is a connected surface. With merge=False, each facet gets its own (n+1)(n+2)/2 points.

values, one per vertex row (shape (n_vertices,) or (n_vertices, c)), are interpolated like the coordinates and returned as a third item.

In torus models, facets that cross the periodic boundary aren’t unwrapped.

Parameters:
tessellate_edges(n=None, *, merge=True, values=None)

Subdivide every edge into n straight segments, for plotting.

Returns (points, segments) (and interpolated values if given), sampled on the curved edges like tessellate(). Segment s belongs to edge s // n. n defaults to 1 for linear elements and 2*order otherwise.

Parameters:
native_cells(reorient=True)

Facets as high-order triangles in Gmsh/VTK node order.

Returns (cell_type, cells) with meshio cell type names ("triangle", "triangle6", "triangle10", …) and cells indexing vertices. With reorient, cells follow the facet orientation of facets.

Parameters:

reorient (bool)

Return type:

Tuple[str, ndarray]

native_edge_cells()

Edges as high-order lines in Gmsh/VTK order (ends, then interior).

Return type:

Tuple[str, ndarray]

body_surfaces(curved='tessellate', n=None, *, cap=False, project=None)

The surface around each body, with outward normals.

Facets with the body in front keep their orientation; facets with it behind are flipped. A film between two bodies appears in both.

Bodies that Evolver closes off with a constraint (for instance a drop on a plane) have open boundary loops. cap=True closes each loop with a cap: rings of triangles between the loop and its centre, about as fine as the loop (curved="tessellate" only). The cap is flat, which is exact for planar loops; project(points, loop_points), returning the points moved onto the surface the loop lies on, makes it follow a curved one (Evolver.body_surfaces() does this for loops on a constraint). cap_ids tells the caps apart. Check watertight.

Parameters:
Return type:

Dict[int, BodySurface]

to_meshio(curved='tessellate', n=None, flavor='gmsh')

The whole surface as a meshio mesh.

Soapfilm surfaces give triangles (or native high-order triangles with curved="native"); the string model gives lines. Cell data: facet_id (or edge_id), front_body and back_body, and Gmsh tags (gmsh:physical = front body, gmsh:geometrical = Evolver element id). flavor ("gmsh", "vtk", "xdmf" or "linear") picks native cell names for the target format.

Parameters:

  • curved (str)

  • n (Optional[int])

  • flavor (str)

Return type:

meshio.Mesh

to_pyvista(n=None, *, point_values=None, cell_values=None)

The surface as a PyVista PolyData, tessellated for curved elements.

Soapfilm surfaces give triangles (cell data facet_id, front_body, back_body); the string model gives lines (cell data edge_id). point_values (one value per vertex row) are interpolated onto the tessellation; cell_values (one per facet, or per edge in the string model) are repeated for its sub-cells.

Parameters:

  • n (Optional[int])

  • point_values (Optional[Dict[str, np.ndarray]])

  • cell_values (Optional[Dict[str, np.ndarray]])

Return type:

pyvista.PolyData

class pysurfaceevolver.MeshQuality(edge_min, edge_median, edge_max, angle_min, angle_median, skinny, degenerate, skinny_angle)

Shape statistics of a mesh’s facets and edges (see Mesh.quality()).

Parameters:
class pysurfaceevolver.Bodies(ids, volume, target, pressure, fixed)

Body data, one entry per body.

Parameters:
class pysurfaceevolver.Quantity(name, value, target, modulus, pressure, kind)

A named quantity declared in the datafile.

Parameters:
class pysurfaceevolver.BodySurface(body, points, cells, cell_type, facet_ids, watertight, open_loops, order=1, bezier=False, cap_ids=None, cap_constraints=<factory>)

The closed surface around one body, oriented with outward normals.

cells are triangles (cell_type == "triangle") or native high-order triangles in Gmsh/VTK node order ("triangle6", "triangle10", …).

Parameters:
body: int
points: ndarray
cells: ndarray
cell_type: str
facet_ids: ndarray
watertight: bool
open_loops: int
order: int = 1
bezier: bool = False
cap_ids: ndarray | None = None
cap_constraints: Dict[int, int | None]
property area: float

The area of the surface (caps included; see areas()).

areas()

The area of each part: index 0 the body’s own facets, k cap k.

Return type:

ndarray

property volume: float

The volume enclosed, from the outward-oriented triangles; NaN unless the surface is closed (watertight, e.g. with caps).

to_meshio(flavor='gmsh')

As a meshio mesh; flavor picks high-order cell names (see native_cell_name()).

Parameters:

flavor (str)

Return type:

meshio.Mesh

volume_mesh(size=None, path=None, *, algorithm=10, names=None)

Tetrahedra filling this closed surface, made with Gmsh (pip install gmsh); the surface triangles are kept as they are.

Physical groups: 1 “surface” for the Evolver facets, 1 + k “cap k” for cap k (see cap_ids), and the body id for the tetrahedra. The result is a meshio mesh with tetrahedra and boundary triangles, their groups in cell_data["gmsh:physical"]. size is the largest element size (default: the mean edge length of the surface); path also writes the mesh (.msh through Gmsh, with the group names; other formats through meshio). algorithm is Gmsh’s Mesh.Algorithm3D (10: HXT, 1: Delaunay). names overrides group names, keyed by cap id (0: the Evolver facets).

Parameters:

  • size (Optional[float])

  • path (Optional[str])

  • algorithm (int)

  • names (Optional[Dict[int, str]])

Return type:

meshio.Mesh

to_pyvista()
Return type:

pyvista.PolyData

pysurfaceevolver.is_watertight(triangles)

True if every edge is shared by exactly two triangles, with opposite directions (a closed, consistently oriented surface).

Parameters:

triangles (ndarray)

Return type:

bool

Building surfaces

pysurfaceevolver.make_datafile(vertices, faces=None, *, edges=None, bodies=None, fixed=None, constraints=None, vertex_constraints=None, edge_fixed=None, edge_constraints=None, parameters=None, gravity=None, boundaries=None, vertex_boundaries=None, string=False, header='', commands='')

Write Surface Evolver datafile text for a surface given as arrays.

Parameters:

  • vertices – (n, sdim) coordinates. A space_dimension line is added when sdim != 3.

  • faces (ndarray | Sequence[Sequence[int]] | None) – Faces as vertex-row loops: an (k, 3) array of triangles, or a list of polygons. Edges and the faces’ oriented edge loops are derived from them. Optional in the string model.

  • edges – Extra (m, 2) edges, for example the curves of a string model.

  • bodies (Iterable[Body | Mapping] | None) – Body objects (or dicts with the same keys).

  • fixed – Vertex rows (or a boolean mask) to fix.

  • constraints (Mapping[Any, Any] | None) – {number: formula}, e.g. {1: "z = 0"}; a key may also be a name ({"floor": ...}). A value that already contains formula: (or function:) is copied as is, so energy and content integrals can be included. A value from pysurfaceevolver.constraints (planes, mirrors, spheres and cylinders with contact angles) writes its integrals itself, and its constants go where they belong: volconst to the bodies that touch it, the left-out wetting energy to Evolver.energy.

  • vertex_constraints (None | Mapping[Any, Iterable] | Sequence[Any]) – {constraint: vertex rows or mask} putting vertices on constraints, or one entry per vertex: the constraint (or a set of them, or None) each vertex is on.

  • edge_fixed – Same for edges, as rows into the derived edge list (see below). By default, a boundary edge (used by at most one face) is fixed, or on a constraint, when both of its vertices are.

  • edge_constraints (Mapping[int, Iterable] | None) – Same for edges, as rows into the derived edge list (see below). By default, a boundary edge (used by at most one face) is fixed, or on a constraint, when both of its vertices are.

  • parameters (Mapping[str, float] | None) – {name: value}, written as parameter name = value.

  • gravity (float | None) – The gravitational constant (gravity_constant); bodies weigh their density times it, along -z.

  • boundaries (Mapping[int, Any] | None) – Curves vertices slide along (Evolver’s boundary; a curved wire for a film), {number: Boundary} from pysurfaceevolver.constraints.circle() or curve, and {number: vertex rows or mask}. A boundary edge (on one face at most) whose ends are both on a boundary is on it too. Vertices on a boundary are fixed (free, they would bunch up along it: that lowers the area of the discrete film); refining puts the new ones on the curve, so the wire stays round.

  • vertex_boundaries (Mapping[int, Iterable] | None) – Curves vertices slide along (Evolver’s boundary; a curved wire for a film), {number: Boundary} from pysurfaceevolver.constraints.circle() or curve, and {number: vertex rows or mask}. A boundary edge (on one face at most) whose ends are both on a boundary is on it too. Vertices on a boundary are fixed (free, they would bunch up along it: that lowers the area of the discrete film); refining puts the new ones on the curve, so the wire stays round.

  • string (bool) – Use the string model (curves; faces become 2D cells).

  • header (str) – Extra datafile lines placed before vertices (gravity, quantities, options, …).

  • commands (str) – Commands placed after read at the end of the file.

  • order (Edges are numbered in this)

  • use. (derived from faces in order of first)

Return type:

str

class pysurfaceevolver.Body(volume=None, faces=None, orientation=None, density=None, volconst=None)

A body for make_datafile().

volume fixes the volume; None leaves it free. faces are 0-based face rows; None (default) means all faces. orientation (+1/-1 per face) says whether each face’s normal, by its vertex order, points out of the body (+1) or into it (-1). By default it is worked out when the faces close a surface (each edge on two of them): consistent across the edges and outward; otherwise all +1, and the body’s starting volume is checked when it loads (see Evolver.load()). density sets the body’s weight density (for gravity). volconst is added to the computed volume (for example a solid inside the body); by default it is the sum of pysurfaceevolver.constraints.Constraint.volconst over the walls from pysurfaceevolver.constraints that the body’s vertices touch, and, with a single body, over the walls no vertex touches too (a solid immersed in the liquid).

Parameters:

Starting shapes

Starting shapes for pysurfaceevolver.make_datafile().

Each builder returns a Shape: vertices, triangles (normals outward, away from the axis or up), and its rims by name (rows of the vertices along each open edge, in order around it), ready for vertex_constraints:

from pysurfaceevolver import shapes, constraints as C
drop = shapes.cap(radius=1, height=0.6)
ev.load_arrays(drop.vertices, drop.faces,
               constraints={1: C.plane((0, 0, 1), 0.0, contact_angle=60)},
               vertex_constraints={1: drop.rims["rim"]},
               bodies=[pyse.Body(volume=0.5)])

The start needn’t have the target volume: relaxing meets it (an open cube of volume 1 relaxes to targets from 0.2 to 10).

class pysurfaceevolver.shapes.Shape(vertices, faces, rims=<factory>)

A triangulated starting shape: vertices (n, 3), faces (k, 3) vertex rows, rims: name -> vertex rows along an open edge.

Parameters:
vertices: ndarray
faces: ndarray
rims: Dict[str, List[int]]
pysurfaceevolver.shapes.cap(radius=1.0, height=1.0, *, nr=8, nt=32, z0=0.0)

A spherical cap of sphere radius and height standing on the plane z = z0 (a drop’s start); rim "rim" on that plane.

Parameters:
Return type:

Shape

pysurfaceevolver.shapes.disc(radius=1.0, *, nr=4, nt=24, z=0.0)

A flat disc in the plane z = z (a film’s start), normal +z; rim "rim".

Parameters:
Return type:

Shape

pysurfaceevolver.shapes.tube(radius=1.0, z0=0.0, z1=1.0, *, nz=8, nt=24)

A cylinder r = radius from z0 to z1 (a liquid bridge’s start), normals outward; rims "bottom" and "top".

Parameters:
Return type:

Shape

pysurfaceevolver.shapes.revolve(r, z, nt=24)

A surface of revolution about the z axis through the meridian points (r_i, z_i), from the bottom ring to the top ring; normals away from the axis. Rims "bottom" and "top". A point with r = 0 at either end closes it there (a pole) instead of a ring.

Parameters:

nt (int)

Return type:

Shape

pysurfaceevolver.shapes.box(size=(1.0, 1.0, 1.0), *, open_bottom=False, n=1)

An axis-aligned box from the origin to size, each face split into n x n squares of two triangles, normals outward. With open_bottom, no face at z = 0, and the rim there is "rim" (a drop or a film on a table).

Parameters:
Return type:

Shape

Walls, mirrors and wires

Constraint builders: planes, mirrors, spheres and cylinders, with contact angles.

A solid wall wetted by the liquid contributes -tension*cos(contact_angle) times its wetted area to the energy, and the wetted part of the wall closes the liquid’s volume. Evolver sees neither directly: both are line integrals along the contact line (a constraint’s energy and content integrands). These builders write those integrands, so a datafile needs only:

from pysurfaceevolver import constraints as C
make_datafile(..., constraints={1: C.plane((0, 0, 1), 0.0, contact_angle=60),
                                2: C.sphere((0, 0, 0), 0.48, contact_angle=40,
                                            span=np.pi/2, wet_poles=("north",))})

Conventions (the usual ones for bodies):

  • facets are oriented with their normals pointing out of the liquid;

  • a plane’s normal points into the liquid (away from the solid); spheres and cylinders are solid, with the liquid outside;

  • contact_angle is in degrees, measured through the liquid; None (or a mirror) means no wetting energy. It may be a datafile parameter’s name or an expression ("theta"), so ev.parameters["theta"] = 50 changes it; then Constraint.energy_constant is NaN and the left-out energy is -tension*cos(theta)*area_constant. An angle given as a parameter, or a number below 20 degrees, enters as maximum(angle, pyse_angle_floor) with a datafile parameter pyse_angle_floor (0) that relax() raises for a while when a crude start at a small angle fails (see there), and Constraint.energy_constant is NaN for it too;

  • volumes are Evolver’s default (z dx dy over facets), not symmetric_content.

The line integrals only see the contact line. Keep mirror planes vertical or in z = 0, where they need no volume closure; a tilted mirror face bounded partly by other walls can’t be closed by line integrals. Faces of the liquid’s region that no contact line touches (a floor under the whole film, say) aren’t closed by any of them: add their part of Evolver’s volume (z dx dy) to the body’s volconst yourself (it is zero for a face in the plane z = 0 or parallel to the z axis). Where the wetted region’s boundary also runs along a mirror plane, the integrands are chosen so that those parts contribute nothing (for spheres and cylinders: mirror planes through the axis, and for spheres the plane through the centre normal to the axis). A sphere’s wetted region may contain a pole (where the azimuth is undefined); say so with wet_poles. pysurfaceevolver.make_datafile() then adds Constraint.volconst to the volume of the body touching the wall (or of the only body), and Constraint.energy_constant to pysurfaceevolver.Evolver.energy (Evolver can’t add a constant energy; pressures don’t need it).

class pysurfaceevolver.constraints.Constraint(formula, energy=None, content=None, volconst=0.0, energy_constant=0.0, nonnegative=False, area_constant=0.0, energy_constant_text=None, parameters=(), shape=None, contact_angle=None)

A constraint for make_datafile(): its datafile text, and the constants its line integrals leave out (see the module notes).

Parameters:
formula: str
energy: Sequence[str] | None = None
content: Sequence[str] | None = None
volconst: float = 0.0
energy_constant: float = 0.0
nonnegative: bool = False
area_constant: float = 0.0
energy_constant_text: str | None = None
parameters: Tuple[str, ...] = ()
shape: Tuple[Any, ...] | None = None
contact_angle: Any = None
property floored: bool

Whether its contact angle takes ANGLE_FLOOR as its least (see relax()).

text()

The constraint’s body in a datafile (after constraint n).

Return type:

str

to_pyvista(bounds=None, resolution=60)

The wall as a PyVista surface, for drawing next to the liquid (see walls in Evolver.plot()): a sphere, a cylinder cut to bounds (xmin, xmax, ymin, ymax, zmin, zmax) or a plane patch covering them. None for a constraint without known geometry.

Parameters:

resolution (int)

Return type:

Any

class pysurfaceevolver.constraints.Boundary(formulas, parameter, content=None)

A curve for make_datafile() (Evolver’s boundary): refining puts new vertices on it, so a film keeps a curved wire however much it is refined (a wire of fixed vertices alone stays the polygon it started as).

formulas give x, y, z in terms of the parameter p1; parameter(xyz) gives each vertex’s parameter from its coordinates (an (n, 3) array in, n values out).

content, three expressions in x, y, z, is a vector field integrated along the body’s edges on the wire and added to its volume (Evolver’s content integral, oriented by the facets as on a constraint): the volume of a face that bounds the body but isn’t in the surface. For a flat face at height h spanning the wire, ("-h*y/2", "h*x/2", "0") adds its z dx dy (a bridge pinned on the rim of a disk at z = h).

Parameters:
formulas: Tuple[str, str, str]
parameter: Any
content: Tuple[str, str, str] | None = None
text()

The boundary’s body in a datafile (after boundary n parameters 1).

Return type:

str

pysurfaceevolver.constraints.plane(normal, offset=0.0, contact_angle=None, *, point=None, ref=None, origin=None, tension=1.0)

The plane normal . x = offset (or through point), the liquid on the side normal points to. With contact_angle, the wetted part adds -tension*cos(contact_angle) per area; the volume is closed over the wetted part.

In-plane coordinates: s along ref (default: x, or y for planes normal to x), t across it, both from origin (default: the point of the plane nearest the coordinate origin). The integrands are t ds forms: they vanish along lines s = const and along t = 0. So where the wetted face’s boundary runs along mirror planes, choose ref and origin so those lines are s = const or t = 0 (a ceiling z = 0.5 cut by mirrors x = 0, x = c and y = 0: the defaults).

Parameters:
Return type:

Constraint

pysurfaceevolver.constraints.mirror(axis, at=0.0)

A mirror plane (no wetting energy): "x" for x = at (also "y", "z"), or (normal, point). The liquid may lie on either side; the volume is closed over the face on the mirror.

Parameters:

at (float)

Return type:

Constraint

pysurfaceevolver.constraints.sphere(center, radius, contact_angle=None, *, wet_poles=(), span=6.283185307179586, tension=1.0, nonnegative=False)

A solid sphere, the liquid outside. The azimuth runs about the z axis through the centre; its poles are center +- radius*z.

wet_poles: the poles (“north”, “south”) inside the wetted region, and span: the azimuthal extent of the computed piece (2 pi for a whole surface, pi/2 for a quarter cut by mirrors through the axis). They set Constraint.volconst and Constraint.energy_constant. nonnegative makes it one-sided (vertices may not enter the sphere).

A coordinate of center may be a datafile parameter’s name or an expression ((0, 0, "1 + gap/2")), so ev.parameters["gap"] = 0.5 moves the sphere (and pysurfaceevolver.recipes.arclength() can follow the parameter). The volume over a wet pole is then in the line integral itself rather than in volconst, which can’t follow the parameter; with both poles wet the centre must be numbers.

Parameters:
Return type:

Constraint

pysurfaceevolver.constraints.cylinder(point, direction, radius, contact_angle=None, *, gauge='azimuthal', ref=None, tension=1.0, nonnegative=False)

A solid cylinder (a fibre) of radius about the line through point along direction, the liquid outside.

Where the wetted region’s boundary runs along a mirror plane, the line integrals must vanish there; pick the form that does:

  • gauge="azimuthal" (default): integrals of l dphi (l along the axis, from point; phi around it). Vanishes on mirrors through the axis and on the mirror normal to the axis through point. For example a vertical fibre through a film, cut by mirrors through its axis.

  • gauge="axial": integrals of phi dl, phi measured from ref (a direction normal to the axis). Vanishes on every mirror normal to the axis and on the mirror through the axis containing ref. For example a fibre along x, cut by mirrors x = const and by z = 0 (ref along y). The wetted region must not cross phi = +-pi (opposite ref).

Parameters:
Return type:

Constraint

pysurfaceevolver.constraints.circle(center=(0.0, 0.0, 0.0), radius=1.0, normal=(0.0, 0.0, 1.0), content=None)

A circular wire of radius about center in the plane normal to normal, parametrized by the angle p1 (radians), with an optional content integral (see Boundary).

Parameters:

radius (float)

Return type:

Boundary

pysurfaceevolver.constraints.curve(x, y, z, parameter, content=None)

A curve x(p1), y(p1), z(p1) given as datafile expressions, with parameter(xyz) mapping vertex coordinates to p1, and an optional content integral (see Boundary).

Parameters:
Return type:

Boundary

Visualization

class pysurfaceevolver.LiveView(ev, scalars=None, element=None, n=None, plotter=None, off_screen=False, mirror=None, walls=None, **mesh_kwargs)

A PyVista window (or notebook widget) that follows the surface.

Call update() to redraw, for example after every few iterations with ev.iterate(100, callback=view.update, every=5). While the facets are unchanged, an update only moves the points (and refreshes the scalars); after refinement or other topology changes the surface is rebuilt. The camera is kept.

Parameters:

  • ev (Evolver) – The Evolver to show.

  • scalars (Scalars) – As for Evolver.plot(); an expression is re-evaluated on every update.

  • element (Optional[str]) – As for Evolver.plot(); an expression is re-evaluated on every update.

  • n (Optional[int]) – As for Evolver.plot(); an expression is re-evaluated on every update.

  • plotter (Optional['pyvista.Plotter']) – An existing pyvista.Plotter; by default a new one is made.

  • off_screen (bool) – Render without a window (for tests and image files).

  • mirror – Planes to show mirror images in (see Evolver.plot()).

  • walls – Walls to draw with it (see Evolver.plot()).

  • **mesh_kwargs (Any) – Passed to Plotter.add_mesh (cmap, show_edges, …).

show()

Open the window (or display the notebook widget) without blocking.

Return type:

Any

update(*_)

Redraw with the current surface. Accepts and ignores callback args.

Parameters:

_ (Any)

Return type:

None

screenshot(filename=None, **kwargs)

Save (or return) an image of the current view.

Parameters:

  • filename (str | None)

  • kwargs (Any)

close()
Return type:

None

Recipes

Workflow recipes built on the core API.

continuation() follows a family of equilibria: it steps a value (a body’s volume, a parameter, …), relaxes after each step, can checkpoint to disk, and yields a record per step:

from pysurfaceevolver import recipes

for step in recipes.continuation(ev, ev.body(1), np.linspace(1.0, 0.2, 50),
                                 checkpoint="runs/drain"):
    print(step.value, step.pressures[0])
    if something_happened(ev):
        ev.load_string(rebuilt_surface)   # same Evolver: the loop goes on with it

A run stopped for any reason resumes from its last finished step with resume="runs/drain".

arclength() follows a branch of equilibria through a fold (where continuation() stops) and on along the unstable branch past it.

class pysurfaceevolver.recipes.Step(index, value, energy, volumes, pressures, result=None, extra=<factory>, substeps=1, predicted=False)

One step of continuation(), after relaxing.

Parameters:
index: int
value: Any
energy: float
volumes: ndarray
pressures: ndarray
result: Any = None
extra: dict
substeps: int = 1
predicted: bool = False
pysurfaceevolver.recipes.arclength(ev, what, ds, steps, *, bounds=None, tol=1e-08, tidy=True, halvings=6, along=None, remesh=None)

Follow a branch of equilibria through folds: pseudo-arclength continuation in what (a parameter or variable name, ev.body(i)’s target volume or ev.quantity(name)’s target).

Where continuation() steps the value and stops at a fold (no equilibrium past it), this steps along the branch itself: each point is a distance ds from the last, measured as the change of the value and the root mean square displacement of the vertices together. So it turns at a fold and goes on along the other branch, typically unstable (a catenoid past its widest span, a drop or bridge past its largest volume). Yields a Step per point, starting with the first two; result.negative_modes counts the Hessian’s negative eigenvalues there with the constraints (volumes, …) held, from Evolver.eigen_counts() (0 on a stable branch, 1 past a simple fold; a symmetry’s zero modes can come out slightly negative, which Evolver.stability() sorts out, at the cost of Evolver’s ritz).

It starts by relaxing the surface at the current value, then takes Newton steps to the value plus ds (the sign sets the direction). From there each point is a prediction along the secant of the last two (where that fails, along the branch’s tangent at the last one), corrected by Newton’s method on the equilibrium together with the condition of being ds along the branch (each iteration: two Newton steps of Evolver’s). No gradient steps: those would leave an unstable branch. Between points (tidy=True) the vertices are spread along the surface and edges flipped, which Newton steps alone don’t do (as Evolver.relax() tidies: no flips in the Lagrange model or for strings; the vertices are spread in the linear model only, since in the others the energy depends on where the nodes sit along the surface and Newton steps don’t move them back).

along (a displacement per vertex, such as ev.modes(1).vectors[0]) starts on another branch: the first step goes a distance ds along it at the starting value instead of in the value, onto a branch that crosses this one there, such as a less symmetric one past a symmetric bifurcation (start just past it, where the mode’s eigenvalue has turned negative: from the point a first run stopped at). The sign of ds picks the direction along it; from the second step on the run follows the new branch, either way in the value.

remesh keeps the edges’ lengths in range on a branch that stretches or shrinks the surface: h (Evolver.remesh()’s target) or a dict of its arguments, called between points when an edge has left the range they set (target=h: 0.5h to 1.6h; edges on the boundary are never deleted), or a function of the surface, called between points. The remeshed surface is brought back to the branch at the same point, and the last step carried over to its new vertices (a split edge’s midpoint takes the mean of its neighbours’); if it doesn’t get back, the surface is put back as it was and the run goes on. Without it the vertices stay the same (refine before). Next to a contact line that shrinks, let it merge contact-line edges too (remesh=dict(target=h, keep_boundary=False)).

Stops after steps points, at the first point outside bounds ((low, high) on the value; the surface is put back at the last one inside), or, with a warning, when a step (the first one too) doesn’t converge, or doesn’t settle again after tidying, even at ds/2**halvings (put back at the last point) or the facets degenerate (below 1 degree: refine and start again from there). A step is also cut when its point lands farther from the prediction than the step itself: on the branch’s plane, but on some other branch or none. tol bounds the residual of each point, as in Evolver.relax().

Parameters:
Return type:

Iterator[Step]

pysurfaceevolver.recipes.continuation(ev, what, values, *, relax=None, checkpoint=None, resume=None, stop_on=None, halvings=0, predictor=True)

Step through values: for each, set what to it, relax, and yield a Step.

what: a parameter or variable name ("theta"), a body’s target volume (ev.body(1)), a fixed quantity’s target (ev.quantity("vol")), or a callable (ev, value) for anything else. relax: options for Evolver.relax() (a dict), a callable relax(ev), or None for ev.relax().

checkpoint: a path prefix; after every step the surface goes to <prefix>.dmp and the steps so far to <prefix>.steps (both written atomically). resume: such a prefix; the surface is loaded from it and the values it finished are skipped (values must be the same sequence). A step’s extra dict is yours to fill in the loop; it is checkpointed with the next step.

Changing ev inside the loop (refining, loading a rebuilt surface) is fine: the next step continues from whatever ev holds.

stop_on="unconverged" (or "unstable": also a converged but unstable equilibrium) ends the run at the first step whose relax result says so, or whose relax pinched the surface off (collapsed edges: a film or bridge breaking past a fold, then converged on as pieces), with the surface put back to the last good step and a warning. halvings=k first retries such a step from the last good one in 2, 4, … 2**k smaller steps (numeric values only); Step.substeps says how many it took.

predictor: start each relax from the secant prediction, every vertex extrapolated from the last two steps (from the third step on, for numeric values and steps at most twice the last, and only while the vertices stay the same: after a refinement, remeshing or pop a step starts from the last equilibrium). It saves a fifth to a half of the gradient steps on smooth continuations. A predicted step that comes out unconverged is relaxed again from the last equilibrium’s positions before anything else; Step.predicted says whether the prediction was kept.

Parameters:
Return type:

Iterator[Step]

pysurfaceevolver.recipes.load_steps(prefix)

The steps a checkpointed continuation() has finished.

Parameters:

prefix (str)

Return type:

List[Step]

Parameter sweeps

pysurfaceevolver.map(fn, params, *, processes=None, threads=None, errors='raise', callback=None)

Run fn(param) for every param in worker processes, in parallel.

fn typically builds or loads a surface with Evolver, runs it, and returns numbers, arrays or a Snapshot. Results come back in the order of params.

processes:

Number of worker processes (default: number of CPUs, at most the number of jobs).

threads:

Threads per worker for Evolver’s parallel loops and Newton steps (default 1: for sweeps, more processes beat threads); sets OMP_NUM_THREADS in the workers.

errors:

"raise" (default): after all jobs finish, raise the first failure as JobError (its results attribute has all results, with failures as JobError objects). "return": put JobError objects in the result list instead.

callback:

callback(index, result) is called in this process as each job finishes.

fn must be picklable: a module-level function, or anything when cloudpickle is installed (which also covers functions defined in notebooks). A job that crashes its worker gives a JobError wrapping WorkerCrashed; the other jobs are unaffected. Workers that die while starting raise WorkerStartError right away (see there for the usual causes).

Parameters:
Return type:

List[Any]

class pysurfaceevolver.JobError(index, param, error, traceback='')

A job raised an exception (or crashed its worker).

index and param identify the job; error is the exception, traceback the worker’s formatted traceback (if any).

Parameters:
exception pysurfaceevolver.WorkerCrashed(exitcode)

A job’s worker process died (for example a segfault) while running it.

Parameters:

exitcode (Optional[int])

exception pysurfaceevolver.WorkerStartError(exitcode)

Worker processes died before they could take a job.

Workers are started with the spawn method, which re-imports the main script in every worker. That fails for a script read from stdin or python -c, and a script whose top-level code calls map() without an if __name__ == "__main__": guard starts workers recursively.

Parameters:

exitcode (Optional[int])

Settings

pysurfaceevolver.set_threads(n)

Threads for Evolver’s parallel loops and the Newton-step solver; n <= 0 restores the default: the physical cores (or OMP_NUM_THREADS / PYSE_THREADS when set).

Parameters:

n (int)

Return type:

None

pysurfaceevolver.threads()

Threads the parallel loops and Newton steps use (1 without OpenMP).

Return type:

int

pysurfaceevolver.threads_limit(n)

Use n threads inside the with block, then restore the previous setting (including “the default”):

with pyse.threads_limit(4):
    ev.relax(newton=10)
Parameters:

n (int)

Return type:

Generator[None, None, None]

pysurfaceevolver.set_solver(name)

Sparse factoring for Newton steps (hessian): "mumps" (the default when pySE is built with MUMPS) or "evolver" (Evolver’s own minimal-degree factoring). Also set by the PYSE_SOLVER environment variable at start-up.

Parameters:

name (str)

Return type:

None

pysurfaceevolver.solver()

The sparse factoring used for Newton steps (see set_solver()).

Return type:

str

pysurfaceevolver.tessellation_limit: int | None = 10000000

Tessellations (plots, live views, export of curved surfaces) with more triangles than this give a LargeTessellationWarning; None: no check.

pysurfaceevolver.busy_timeout: float | None = None

Seconds a call waits while Evolver runs a call from another thread, then EvolverBusyError; None: fail at once; float(“inf”): wait as long as needed.

Errors and warnings

exception pysurfaceevolver.EvolverError(message, errnum=0, output='')

An error reported by Surface Evolver. The surface is still usable.

command is the command text that failed (set by Evolver.command()).

Parameters:
exception pysurfaceevolver.InvalidSurfaceError(message, errnum=0, output='')

There is no valid surface: the last datafile failed to load, or an unrecoverable error happened. Only Evolver.load() works until a datafile loads successfully.

Parameters:
exception pysurfaceevolver.EvolverFatalError(message, errnum=0, output='')

An unrecoverable Evolver error. Load a new datafile before going on.

Parameters:
exception pysurfaceevolver.EvolverExit(code, output='')

Evolver tried to exit, for example because of the quit command.

The process keeps running and the engine is still usable.

Parameters:
exception pysurfaceevolver.EvolverBusyError

Evolver is running another call: from another thread (and pyse.busy_timeout is None, or ran out), or the running call itself (a re-entrant call from a callback). A subclass of RuntimeError.

exception pysurfaceevolver.EvolverWarning

A warning printed by Surface Evolver.

exception pysurfaceevolver.LargeTessellationWarning

A tessellation has more triangles than pyse.tessellation_limit.

exception pysurfaceevolver.UnstableEquilibriumWarning

relax() converged to an equilibrium that some motion makes lower in energy (a saddle): see Evolver.stability().

exception pysurfaceevolver.SharedEngineWarning

An Evolver handle was used after another handle loaded a different surface: every handle is a view of the one engine in the process.

Sample files

Surface Evolver’s sample datafiles and command scripts, installed with the package.

The directory is also appended to EVOLVERPATH when the package is imported, so Evolver("cube.fe") and commands like read "obj.cmd" find these files from any working directory (files in the working directory, and directories already on EVOLVERPATH, come first).

>>> from pysurfaceevolver import Evolver, examples
>>> examples.names()[:3]
['100grain.fe', 'addload_example.fe', 'cat.fe']
>>> ev = Evolver(examples.path("cube.fe"))
pysurfaceevolver.examples.directory()

The directory holding the sample files.

Return type:

Path

pysurfaceevolver.examples.names(pattern='*.fe')

Names of the sample files matching a glob pattern (datafiles by default; "*.cmd" lists the command scripts).

Parameters:

pattern (str)

Return type:

List[str]

pysurfaceevolver.examples.path(name)

Full path of a sample file, e.g. path("cube.fe").

Parameters:

name (str)

Return type:

str