ffc team mailing list archive

[Jake Ostien <tostien@xxxxxxxxx>]

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

Thanks,
Jake