ffc team mailing list archive

[Jake Ostien <tostien@xxxxxxxxx>]

Hi,

I have a form that takes in a symmetric stress with the intention of 
computing the deviatoric part of that stress.

It looks something like this:

el = VectorElement("DG", "tetrahedron", 0, 6)

V = TestFunction(el)
U = TrialFunction(el)

def Mat(A):
   return [ [A[0], A[3], A[5]],    \
            [A[3], A[1], A[4]],    \
            [A[5], A[4], A[2]] ]

def Vec(A):
   return [ A[0][0], A[1][1], A[2][2], A[0][1], A[1][2], A[2][0] ]

T = Mat(Function(el))

dev = T - (1.0/3.0)*mult(trace(T), Identity(3))

a = dot(V,U)*dx
L = dot(V,Vec(dev))*dx


Now if I turn debugging on in ffc with the -d1 flag, I can see that 
something is fishy with the simplified form (where as it seems OK in the 
'raw' input), and when I look at the element tensor generated for the 
linear form it looks like

   A[0] = 0.111111111111111*G0_0 - 0.0555555555555555*G1_1 - 
0.0555555555555555*G2_2;
   A[1] = 0.111111111111111*G0_1 - 0.0555555555555555*G3_0 - 
0.0555555555555555*G4_2;
   A[2] = 0.111111111111111*G0_2 - 0.0555555555555555*G5_0 - 
0.0555555555555555*G6_1;
   A[3] = 0.111111111111111*G0_3;
   A[4] = 0.111111111111111*G0_4;
   A[5] = 0.111111111111111*G0_5;

But what it _should_ look like is

   A[0] = 0.111111111111111*G0_0 - 0.0555555555555555*G1_1 - 
0.0555555555555555*G2_2;
   A[1] = 0.111111111111111*G0_1 - 0.0555555555555555*G3_0 - 
0.0555555555555555*G4_2;
   A[2] = 0.111111111111111*G0_2 - 0.0555555555555555*G5_0 - 
0.0555555555555555*G6_1;
   A[3] = G0_3;
   A[4] = G0_4;
   A[5] = G0_5;

There is essentially an extra factor of 2/3[*] multiplied on what turn 
out to be shear stresses.  I've been wading around the simplify part of 
the code for a little and I'll continue to look, but I would welcome any 
suggestions on my search direction.

Thanks,
Jake

[*] start with 0.1111111 or 1/9, multiply by through by 6 from the mass 
matrix = 2/3