Projection and symmetry detection

Finding the closest tensor of a given material symmetry, and identifying which symmetry a tensor has. The algorithms are on Projection onto a symmetry class.

proj_tens

B, d, drel = proj_tens(sym, A)              # orientation optimized — needs NLopt
B, d, drel = proj_tens(sym, A, n_or_frame)  # orientation given

with sym ∈ (:ISO, :TI, :ORTHO, :CUBIC), A an order-2 or order-4 array or tensor, and

Returned	Meaning
B	the projection, in the storage type of the class
d	absolute Frobenius distance
drel	relative distance

Always look at drel, not d. The relative distance is dimensionless, so a tolerance on it is independent of the units of the moduli.

using TensND, LinearAlgebra

n = [0.0, 0.0, 1.0]
C = get_array(tens_TI(10.0, 3.0, 2.5, 12.0, 2.0, n))

B, d, drel = proj_tens(:TI, C, n)
(typeof(B), drel)

(TensTI{4, Float64, 5}, 8.709305949236423e-17)

Projecting onto a class the tensor already belongs to is exact. Projecting about the wrong axis is not:

[round(proj_tens(:TI, C, [sind(α), 0.0, cosd(α)])[3], digits = 5) for α in (0, 15, 30, 45, 90)]

5-element Vector{Float64}:
 0.0
 0.14676
 0.20301
 0.1887
 0.18989

That function of the orientation is what the optimizer minimizes.

Orientation given, or found

Without NLopt, the axis or frame must be supplied. With it, the two-argument form searches:

using NLopt

θ, ϕ = 0.6, 1.1
n_tilt = [sin(θ)cos(ϕ), sin(θ)sin(ϕ), cos(θ)]
C_tilt = get_array(tens_TI(10.0, 3.0, 2.5, 12.0, 2.0, n_tilt))

B_opt, _, drel_opt = proj_tens(:TI, C_tilt)
(round.(collect(axis(B_opt)), digits = 8), drel_opt)

([0.25611964, 0.50321353, 0.82533561], 2.169465605562266e-15)

The axis is recovered up to a sign, which is immaterial: and define the same symmetry. Similarly an orthotropic frame is determined only up to the reorderings and sign flips that map the three orthogonal planes onto themselves.

The search is a deterministic multi-start, so repeated calls return bit-identical results:

length(unique([proj_tens(:TI, C_tilt)[3] for _ in 1:5]))

1

best_sym_tens

Tries the classes from the most restrictive to the least and returns the first whose relative error falls below ε:

B, d, drel, sym = best_sym_tens(t; proj = (:ISO, :CUBIC, :TI, :ORTHO), ε = 1e-6,
                                optimize_angles = false)

sym is one of :ISO, :CUBIC, :TI, :ORTHO, :ANISO. The default order is by number of constants — 2, 3, 5, 9 — so the tightest class that fits within ε is the one reported. :CUBIC and :TI are incomparable, neither containing the other, so a tensor satisfying both is reported cubic, having the fewer constants. The argument must be an AbstractTens, not a bare array — wrap a raw array with Tens first.

best_sym_tens(Tens(C_tilt))[4]

:TI

optimize_angles	axis / frame	needs NLopt
false (default)	from the tensor if it is a structured container, otherwise from the Kelvin–Mandel eigenstructure	no
true	found by multi-start optimization	yes

The cheap path is usually enough

The eigenstructure candidate is exact whenever the tensor genuinely has the symmetry sought, so the default path identifies tilted TI and rotated orthotropic tensors correctly and roughly an order of magnitude faster. Reach for optimize_angles = true when the tensor is only approximately of the class and the best orientation is itself the question.

Cubic symmetry

:CUBIC takes a cube frame, like :ORTHO, and returns a TensCubic — or a TensISO at order 2, the two classes coinciding there.

ℬ = CanonicalBasis{3, Float64}()
cub = tens_cubic(10.0, 4.0, 2.0, ℬ)
B, d, drel = proj_tens(:CUBIC, get_array(cub), ℬ)
drel

7.30078163725866e-17

Omitting the frame optimizes over cube orientations, exactly as for a TI axis or an orthotropic frame, and needs NLopt in the same way. The octahedral group is discrete but the orientation of the cube is not — it is an ordinary rotation — so the objective is smooth in the Euler angles; the group shows up only as a 24-fold degeneracy of the minimum, and all 24 frames describe the same tensor with the same coefficients.

The residual is the number worth reporting. When the morphology and the medium leave the octahedral group invariant, the answer belongs to the class by group theory, so drel is bounded by the discretization error and by nothing else — an error estimate with no reference solution in it. What lives inside the class, and is therefore not measured by drel, is cubic_anisotropy: the two together separate a real material anisotropy from a numerical artifact, since an artifact breaks the symmetry while a real anisotropy does not.

Predicates

is_ISO, is_TI, is_ORTHO and is_CUBIC are the same computation with a boolean answer, and they respect the hierarchy ISO ⊂ TI ⊂ ORTHO:

𝕀, 𝕁, 𝕂 = ISO(Val(3), Val(Float64))
(is_ISO(𝕀), is_TI(𝕀), is_ORTHO(𝕀))

(true, false, false)

(is_ISO(C), is_TI(C), is_ORTHO(C))

(false, true, true)

Each accepts an optional axis or frame, and the optimize_angles keyword.

Two things worth stating in a report

Which tensor was projected. The Euclidean distance is not invariant under inversion, so projecting a stiffness and projecting the corresponding compliance give different materials:

ort = TensOrtho(10.0, 8.0, 9.0, 3.0, 2.0, 4.0, 2.5, 3.0, 1.5, CanonicalBasis{3, Float64}())
Biso, _, _ = proj_tens(:ISO, get_array(ort))
Siso, _, _ = proj_tens(:ISO, get_array(inv(ort)))
round(norm(get_array(inv(Biso)) - get_array(Siso)) / norm(get_array(Siso)), digits = 4)

0.0951

Alternatives invariant under inversion are discussed in Isotropic tensors.

Whether you wanted a projection at all. Projecting onto the five-parameter major-symmetric TI subspace is a best fit: it forces   and discards . If the object is not major-symmetric — a strain-concentration tensor, typically — an exact average over the rotation group preserves that content where a best fit destroys it. The distinction is developed on The extended Walpole algebra and Projection onto a symmetry class.