Estimating tensor function in Fenics

I have a problem that people who know me better at Fenics can probably solve quickly, and I would be very grateful. I am trying to define a spatially dependent elastic tensor (C_ijkl). After assembling a tensor, when I draw a specific component of it (say C_1100) using the fenics plot command, it works, but if I try to evaluate it at some point in the domain, then I get an error. Code:

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

with open('Material.txt', 'r') as f:
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)

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 that, the following actions are performed:

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'm missing something rather trivial. Any light on this would be much appreciated.

+3

source to share

All Articles