Ph or ? for anisotropic property by region

Dear freefem community,

I would like to have your advices for the best implementation of anisotropic physical property for each material.

In fact I imagine this for attributing physical property to each material, see below, is it correct ?

real [int, int][int] z(4)
for (i=1;i<=3;i++) {
z[i].resize(4,4); z[i] = 0.; // for instance 3 materials (I don’t use 0 index)
}

// Def for a physical property to whatever i material
macro Z(i) [[ z[i] (1,1), z[i] (1,2), z[i] (1,3) ],
[ z[i] (2,1), z[i] (2,2), z[i] (2,3) ],
[ z[i] (3,1), z[i] (3,2), z[i] (3,3) ]] // EOM

// and the varf for each material…1,2,3
varf VZ(u,v) = int3d(Th,1)( grad(v)’ *(Z(1)*grad(u)) )
+ int3d(Th,2)( grad(v)’ *(Z(2)*grad(u)) )
+ int3d(Th,3)( grad(v)’ *(Z(3)*grad(u)) );
Is the proposal of my previous post correct ?

Or should I use definition with

Ph Phy = Phy1*(region==1)+Phy2*(region==2) + etc. But the latter definition is well documented only for scalar values even for elasticity in isotropic conditions with E, nu or lambda, mu.

In fact, I am searching for the best way to properly model these anisotropic properties, especially when some solution fields are dependent of such region property (for instance stress field).

Best regards.