Subject: [FFC-dev] possible bug in simplify

From: Jake Ostien <[email-masked]>
Date: Sun, 10 Feb 2008 20:36:27 -0500

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

---

From: Marie Rognes <[email-masked]>
Date: Mon, 11 Feb 2008 11:18:08 +0100

>     A[2] = 0.111111111111111*G0_2 - 0.0555555555555555*G5_0 - 

I have attempted to fix this now... The simplified and the 
non-simplified form now gives equivalent results as far as I can see.

However, are you sure that the factors should be as you have stated 
above, and not the following:

        // Compute element tensor
        A[0] = 0.166666666666666*G0_0 - 0.0555555555555555*G1_0 -
    0.0555555555555555*G2_1 - 0.0555555555555555*G3_2;
        A[1] = 0.166666666666666*G0_1 - 0.0555555555555555*G1_1 -
    0.0555555555555555*G4_0 - 0.0555555555555555*G5_2;
        A[2] = 0.166666666666666*G0_2 - 0.0555555555555555*G1_2 -
    0.0555555555555555*G6_0 - 0.0555555555555555*G7_1;
        A[3] = 0.166666666666666*G0_3;
        A[4] = 0.166666666666666*G0_4;
        A[5] = 0.166666666666666*G0_5;
      }

?

-- 
Marie E. Rognes
Ph.D Fellow, 
Centre of Mathematics for Applications, 
University of Oslo
http://folk.uio.no/meg

---

From: Jake Ostien <[email-masked]>
Date: Mon, 11 Feb 2008 09:39:08 -0500

Marie Rognes wrote:
>> 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
>>
Yes, of course, thanks for catching that.

Thanks,
Jake