Evaluating a tensor valued function at a point

Question

Hi, I have a problem which people better versed than me in Fenics can probably solve quickly, and I'd appreciate it very much. I am trying to define a spatially dependent elasticity tensor (C_ijkl). After assembling the tensor, when I plot a particular component of it (let's say C_1100) using the fenics plot command then it works but if I try to evaluate it at some point within the domain then I get an error. The code is:

mesh = Mesh("geometry.xml")
cd = MeshFunction('size_t',mesh,"geometry_physical_region.xml")

def readMP2():
 with open('Material.txt', 'r') as f:
    N=([int(x) for x in f.readline().split()])[0];
    rhoL=[];

    for i, line in enumerate(f):
        if i==N-1:
            break
        rhoL.append(([float(x) for x in line.split()])[0])

    rhoL.append(([float(x) for x in line.split()])[0])
    rho=np.asarray(rhoL)

    lmL=[];

    for i, line in enumerate(f):
        lmL.append(([float(x) for x in line.split()])[:2])
        if i==N-1:
            break

    lm=np.asarray(lmL)

    return (rho,lm)

r,lm = readMP2()

V0=FunctionSpace(mesh, 'DG', 0)
M0=TensorFunctionSpace(mesh, 'DG', 0, shape=(2,2,2,2))
rho,lam,mu=Function(V0),Function(V0),Function(V0)
C=Function(M0)
i=Index()
j=Index()
k=Index()
l=Index()

delta=Identity(2)

rho.vector()[:] = numpy.choose(numpy.asarray(cd.array(), dtype=numpy.int32), r)
lam.vector()[:] = numpy.choose(numpy.asarray(cd.array(), dtype=numpy.int32), lm[:,0])
mu.vector()[:] = numpy.choose(numpy.asarray(cd.array(), dtype=numpy.int32), lm[:,1])

C=as_tensor((lam*(delta[i,j]*delta[k,l])+mu*(delta[i,k]*delta[j,l]+delta[i,l]*delta[j,k])),[i,j,k,l])

After this the following works:

plot(C[1,1,0,0])
interactive()

But if I try to do the following:

C1=C[0,0,0,0]
print C1(.0001,.0001)

then I get the following error:

ufl.log.UFLException: Expecting dim to match the geometric dimension, got dim=1 and gdim=2.

I feel like I am missing something rather trivial. Any light on this would be very appreciated

Answer 1

Please post a complete, but minimal code example. There's a lot of junk here.

To evaluate the first component you could try either one of

print C.sub(0)(0.01, 0.01)
print C(0.01, 0.01)[0]