API — Superspheres

The superspherical and superspheroidal shape family: closed-form geometry, with no mesh and no solver involved. What turns one of these shapes into a phase of an RVE is a finite-element solve of its Eshelby problem, or a surrogate trained on those solves.

The exponent is , not — see the warning on Supersphere.

MeanFieldHomogenization.Superspheres Module

MeanFieldHomogenization.Superspheres

Superspherical and superspheroidal morphologies: shapes of the family

which interpolate between the sphere ( ), the octahedron or double cone ( ), the concave range ( ) and the cube or cylinder ( ).

These are not ellipsoids and have no closed-form Eshelby solution, which is the whole reason they are here: they are the family the package's non-ellipsoidal routes were built for. A shape becomes a phase of an RVE either through a finite-element solve of its Eshelby problem or through a surrogate trained on those solves.

This file holds the geometry alone — closed-form level set, radial map, normals, characteristic radii, and the exact volume and shadow area that any mesh is measured against. Nothing here meshes or solves.

Convention: the exponent is , after the Sevostianov–Giraud–Chen–Grgic literature. See the warning on Supersphere.

MeanFieldHomogenization.Superspheres.AbstractSuperShape Type

AbstractSuperShape{T}

A body described by a radial map: radial_distance gives, for a unit direction , the distance from the center to the surface.

Every supersphere and superspheroid is star-shaped about its center for every  , the concave ones included. That is what makes the radial map single-valued, and it is the property the whole meshing strategy rests on — a surface mesh is obtained by pushing a subdivided sphere radially onto the shape, with no risk of the map folding over.

MeanFieldHomogenization.Superspheres.Supersphere Type

Supersphere(a, p)

The body    , of cubic symmetry ().

The exponent is , not . This is the convention of the Sevostianov–Giraud–Chen–Grgic literature these shapes come from, so

	shape
	the sphere
	the regular octahedron
	concave, with conical points on the axes
	the cube

Two conventions for p

scripts/common/docviz.jl calls p the exponent itself, so its sphere sits at p = 2. The two differ by a factor of two; shape_exponent is what converts.

A Supersphere is a shape, not yet an inclusion: it has no closed-form Eshelby solution. It becomes a phase of an RVE through the finite-element or surrogate-backed inclusion built on it.

MeanFieldHomogenization.Superspheres.Superspheroid Type

Superspheroid(a, c, p)

The axisymmetric body    with   , transversely isotropic about .

Same convention for as Supersphere:   is the spheroid of semi-axes ,   the double cone,   the concave range,   the cylinder.

MeanFieldHomogenization.Superspheres.shape_exponent Function

shape_exponent(s) -> 2p

The exponent that appears in the level set, as opposed to the concavity parameter that names the shape. See the warning on Supersphere.

MeanFieldHomogenization.Superspheres.is_concave Function

is_concave(s) -> Bool
is_convex(s)  -> Bool
is_sphere(s)  -> Bool

Shape classification, from alone: concave below  , convex at or above it, and a sphere exactly at  .

These are comparisons, so they need an element type on which a comparison returns an honest Bool — every floating-point type and ForwardDiff.Dual, but not a symbolic one, where they throw rather than answer wrongly.

MeanFieldHomogenization.Superspheres.has_coordinate_mirrors Function

has_coordinate_mirrors(s) -> Bool

Whether the three coordinate planes are mirror planes of s.

This is the condition that licenses meshing one octant of a cell built around s, and it is weaker than cubic symmetry: it asks for the group {diag(±1,±1,±1)} of order 8 and nothing more. A Superspheroid has it — its level set depends only on √(x²+y²) and |z| — without having octahedral symmetry at all. What O_h adds is the symmetry class of the measured tensor, not the validity of the subdivision.

The default is false, deliberately. A new shape has to opt in by adding a method, because a default of true would turn an oversight into an answer that is wrong and plausible. check_coordinate_mirrors audits a candidate numerically before its method is written; it is not a substitute for the trait.

Two prerequisites go with it, and they are the caller's to keep: the shape is centered on the origin with the outer boundary concentric, and the reference medium shares the same mirrors. The second is free today because the corrected cell only accepts an isotropic reference; an anisotropic one would have to be orthotropic in the shape's own axes.

MeanFieldHomogenization.Superspheres.check_coordinate_mirrors Function

check_coordinate_mirrors(s; ndirs = 2000, atol = 1.0e-13) -> Bool

Numerical audit of has_coordinate_mirrors: radial_distance must be invariant under each of the three reflections, over a deterministic sweep of directions.

An audit, not a guard. A finite sample can miss a localized defect, and letting this function decide would do exactly what the trait exists to prevent: correct silently. Use it when writing a has_coordinate_mirrors method for a new shape, and as a test on the shapes already shipped.

MeanFieldHomogenization.Superspheres.level_set Function

level_set(s, x) -> f

  strictly inside,   on the surface,   outside:   .

MeanFieldHomogenization.Superspheres.radial_distance Function

radial_distance(s, n) -> r

Distance from the center to the surface along the unit direction .

No root finding is involved. The level set is positively homogeneous of degree , so    solves in closed form:

That closed form is why the radial map is cheap enough to apply at every node of a refined surface mesh, and why it differentiates cleanly with respect to the shape parameters.

MeanFieldHomogenization.Superspheres.surface_point Function

surface_point(s, n) -> x

  with from radial_distance: the radial projection of a unit direction onto the surface. This is the map that turns a subdivided sphere into a supersphere.

MeanFieldHomogenization.Superspheres.outward_normal Function

outward_normal(s, x) -> n

Unit outward normal at a surface point, from .

which is finite and non-zero as long as no coordinate vanishes. Where one does, the surface is generally not differentiable, and the two ways it fails are opposite:

• for   the vanishing coordinate's derivative vanishes too, and it contributes exactly nothing — the normal is well defined and comes from the surviving coordinates;

• for   it diverges, because the surface creases on the coordinate planes and spikes at the conical points on the axes. There the normal is not unique and no formula can invent one, so the cases are named rather than computed:

  vanishing coordinates	feature	returned
  one, say	a crease
  two or three	a conical point on an axis	the radial direction

  The crease's one-sided limits are as   and as  ; with exactly zero there is nothing to choose between them and is the convention. The conical point's normals span a cone, and the radial direction is its axis.

Writing the raw formula and normalizing afterwards does not work: at such a point it forms 0 * Inf, Julia returns NaN, and a NaN normal that reaches a mesh generator produces a negative Jacobian somewhere far away from the cause.

On a symbolic element type the case analysis is skipped — it is a set of comparisons — and the normalized gradient is returned as it stands.

MeanFieldHomogenization.Superspheres.diagonal_radius Function

diagonal_radius(s::Supersphere)

Radius along ,  .

The cubic group has three critical directions — axis, edge midpoint, body diagonal — and with   their radii are

which is monotone in . So the edge midpoint always lies between the other two, and the extremal radii of a supersphere are reached on the axis and on the diagonal, never in between — which is what makes bounding_radius and inner_radius a two-term comparison rather than an optimization.

It is also the shape's cheapest signature:   gives ,   gives , and   gives , the half-diagonal of the cube of half-side .

MeanFieldHomogenization.Superspheres.edge_radius Function

edge_radius(s::Supersphere)

Radius along ,  : the case   of diagonal_radius, always between and the diagonal radius.

MeanFieldHomogenization.Superspheres.bounding_radius Function

bounding_radius(s)

Smallest such that the body fits inside the ball of radius .

MeanFieldHomogenization.Superspheres.inner_radius Function

inner_radius(s)

Largest such that the ball of radius fits inside the body.

MeanFieldHomogenization.Superspheres.shape_volume Function

shape_volume(s) -> V

Exact volume, in closed form:

Both give at  , and at   they give (the octahedron) and (the double cone). This is the oracle a meshed volume is compared against — the mesh is what is approximate here, not this.

Uses , so it needs a floating-point or ForwardDiff.Dual element type; a symbolic one has no loggamma.

MeanFieldHomogenization.Superspheres.projected_area Function

projected_area(s::Supersphere) -> S

Exact area of the shadow on a coordinate plane, that is, the area enclosed by the superellipse   :

  gives ;   gives , the square of diagonal .

MeanFieldHomogenization.Superspheres.equivalent_sphere_radius Function

equivalent_sphere_radius(s)

Radius of the ball of the same volume, .

Two uses: normalizing contribution tensors the way the literature does, and setting a mesh size that means the same thing across — an element size tied to would be far coarser, relative to the body, at   than at  .

The triangulated surface

No CAD kernel can represent |x|^{2p} + |y|^{2p} + |z|^{2p} = a^{2p}, so the surface is discretized here and handed to a mesher as a discrete entity. What it is discretized from is a subdivided octahedron, for three reasons at once: its face edges lie exactly in the coordinate planes, so an octant is native and cubic symmetry is exploitable with no tolerance offset; for a concave   those same planes are where the surface creases, so mesh edges land on the creases rather than straddling them; and there is no polar degeneracy. And the octahedron is the   supersphere, which makes that value a bit-exact oracle rather than a converging one.

Unlike the geometry above, this part is deliberately not type-generic: a mesh exists to be handed to a mesher, which wants floating point.

MeanFieldHomogenization.Superspheres.TriSurface Type

TriSurface{T}

A triangulated surface. nodes are points, tris are 1-based node triples oriented so that     points outward.

constraint[i] records what a node may do under relaxation, and it is the field that makes octant symmetry usable:

value	meaning
0	free
1, 2, 3	confined to the plane  ,  ,
4	pinned — an axis vertex of the shape, shared by two planes

A node on a plane has the corresponding coordinate exactly zero, and stays that way through every operation here. Nothing rounds it back.

MeanFieldHomogenization.Superspheres.node_count Function

node_count(surf) -> Int
triangle_count(surf) -> Int

Sizes of a TriSurface. Spelled out rather than called nnodes and ntris, which would collide with Ferrite's own names for anyone who loads both.

MeanFieldHomogenization.Superspheres.triangle_count Function

See node_count.

MeanFieldHomogenization.Superspheres.octant_patch Function

octant_patch(level; T = Float64) -> TriSurface

The     face of the unit octahedron, subdivided level times into triangles and projected onto the unit sphere.

Nodes come from the barycentric lattice with     , normalized. A node with   therefore has exactly zero — not zero to a tolerance — which is what makes the coordinate planes exact symmetry planes of the mesh and an octant computation legitimate.

MeanFieldHomogenization.Superspheres.unit_octahedron Function

unit_octahedron(level; T = Float64) -> TriSurface

The whole subdivided octahedron projected onto the unit sphere, obtained by reflecting octant_patch through the three coordinate planes and merging the shared nodes.

Nodes on a plane merge exactly, because the reflection maps their vanishing coordinate to itself. One subtlety makes that true in floating point rather than almost true: -1.0 * 0.0 is -0.0, which compares equal to 0.0 but hashes differently, so the merge key adds 0.0 to normalize the sign of zero. Without that the three planes carry doubled nodes and the surface is not closed.

MeanFieldHomogenization.Superspheres.project_to_shape! Function

project_to_shape!(surf, shape) -> surf

Move every node radially onto the surface of shape.

Valid for every  , concave shapes included, because these bodies are star-shaped about their center and radial_distance is therefore single-valued in a direction. A projection along the normal, or any iterative snapping, would have to worry about folding; this one cannot fold.

MeanFieldHomogenization.Superspheres.shape_surface Function

shape_surface(shape, level; octant = false, relax = 0, kwargs...) -> TriSurface

The discretized surface of shape: a subdivided octahedron — or a single octant patch, when octant = true — radially mapped, and relaxed relax times. kwargs are forwarded to relax_surface!.

MeanFieldHomogenization.Superspheres.relax_surface! Function

relax_surface!(surf, shape; iterations = 20, step = 0.5, density = nothing) -> surf

Equalize the mesh by a tangential umbrella sweep followed by a radial re-projection, iterations times.

The raw radial map is far from uniform: it compresses the mesh near the axis vertices of a concave shape and stretches it near the body diagonal.

What the sweep buys, and what it costs. It equalizes edge lengths, and it does so at the expense of the minimum angle — the two are not the same quality and cannot both be maximized. Measured at level 3, 30 iterations:

shape		smallest angle
Supersphere(1, 0.30)	4.38 → 3.41	19.8° → 16.9°
Supersphere(1, 0.25)	8.07 → 6.09	12.7° → 9.4°
Superspheroid(1, 2.5, 0.35)	7.96 → 4.92	15.0° → 11.7°

So use it when element size uniformity is what matters — which is the usual case, a mesher sizing its tetrahedra from the boundary edges — and leave it off (relax = 0, the default) when the smallest angle is the binding constraint.

Two details separate a sweep that helps from one that wrecks the mesh, and both were established by getting them wrong first.

• The displacement is projected onto the tangent plane before the node moves. Its normal component is destroyed by the radial re-projection anyway, so keeping it only adds a drift that fights the curvature.

• A node confined to a coordinate plane relaxes only against its neighbors on that same plane, as a one-dimensional chain along the boundary curve. Averaging it against interior neighbors drags it along the curve unevenly and destroys an already-perfect mesh. At  , where the radial map is the identity and the raw mesh is the exact barycentric lattice, that mistake turns 60° angles into 30° ones.

With both in place the sweep is a fixed point of the   lattice, which is the property to keep and is asserted by the tests.

density(x) -> ρ biases the average toward high , so nodes cluster where is large; pass   for a target size field . The default, nothing, equalizes edge lengths.

MeanFieldHomogenization.Superspheres.mesh_area Function

mesh_area(surf)

Sum of the triangle areas.

Deliberately not compared against a closed form: a triangulation of a curved surface always under-resolves it, and unlike the volume there is no exact area for a supersphere to compare with anyway. Used for self-convergence only.

MeanFieldHomogenization.Superspheres.mesh_volume Function

mesh_volume(surf)

Volume enclosed by the surface,    over outward-oriented triangles.

For an octant patch this is the volume of the cone from the origin over the patch, which is exactly one eighth of the body — the three flat coordinate faces pass through the origin and contribute nothing to the divergence integral. So 8 * mesh_volume(patch) compares directly with shape_volume, and that identity is what says the octant and the full cell describe the same object.

MeanFieldHomogenization.Superspheres.edge_lengths Function

edge_lengths(surf) -> Vector{Float64}

Every distinct edge length, for a quality report.

MeanFieldHomogenization.Superspheres.mesh_quality Function

mesh_quality(surf) -> NamedTuple

(; nodes, triangles, area, volume, hmin, hmax, hmean, ratio, min_angle).

ratio = hmax/hmin and min_angle, in degrees, are the two numbers that decide whether a mesher will accept the surface as a volume boundary — a sliver on the boundary propagates into the tetrahedra and produces a negative Jacobian somewhere inside.

MeanFieldHomogenization.Superspheres.boundary_chains Function

boundary_chains(surf) -> NTuple{3, Vector{Int}}

The three ordered boundary chains of an octant_patch; chains[k] is the one lying in the plane  , running between the two pinned corners that bound it and including both.

Found from the mesh rather than from the subdivision lattice — a boundary edge is one belonging to a single triangle — so it survives relax_surface! and project_to_shape!, neither of which preserves the barycentric indexing.

This is what lets an octant cell be sewn together exactly: the chain is shared, node for node, between the curved patch and the flat face that meets it, so the two carry the same nodes rather than two discretizations of one curve.

MeanFieldHomogenization.Superspheres.patch_corners Function

patch_corners(surf) -> NTuple{3, Int}

The three axis vertices of an octant_patch, ordered by the axis each one lies on.

They are the nodes carrying constraint == 4: pinned, because each is shared by two coordinate planes and may move along neither.