# Zone property values not updating correctly?

I’m currently trying to assign shear-normal values based on a function for each zone, with the function based on vertical stress values. The result looks correct, but as a sanity check, I plotted the zone property shear-normal, and the values are constant for the whole model.

Is the below code correct? Or could there be an issue with how I am plotting the value?

fish define shear_stiffness
; loops through the zone list - with the current zone called adr
local syy = zone.stress.yy(adr) ; finds zone vertical stress
Ks = 10 ^ (0.093 * syy + 0.676) ; calculates shear stiffness
shear_m = 1 / ((2 * (1 + p_ratio)) / young_plane + (1 / (Ks * S))) ; calculates shear modulus
zone.prop(adr, ‘shear-normal’) = shear_m ; sets zone shear modulus value
endloop
end

Ok, I’ve narrowed it down to the following line:

local syy = zone.stress.yy(adr) ; finds zone vertical stress

When saving this as an extra (see code below), it results in 0 at every zone when plotting, whereas it should match the stress yy plot. Stay tuned while I continue troubleshooting (at least until someone can point out my mistake).

The isotropic (Mohr-coulomb) model maps correctly, matching the stress yy plot. However, for the anisotropic models, the result is zero in every zone. My current working theory is that the model properties are out. See below—stay tuned while I work through this.

[global young_plane=1.25e4, p_ratio=0.25, Kn=80, S=60]

zone cmodel assign mohr-coulomb
zone property density 2700 young [young_plane] poisson [p_ratio] …
friction 50 cohesion 10.5 tension 0

zone cmodel assign anisotropic
zone property young-plane [young_plane] young-normal [1 / (1 / young_plane + 1 / (Kn * S))]
zone property poisson-plane [p_ratio] poisson-normal [p_ratio] …
shear-normal [1 / ((2 * (1 + p_ratio)) / young_plane + (1 / (5.87 * S)))]
zone property angle 0