### Construction of a 4th-order tensor in FEniCS

588
views
2
11 months ago by
Hello everyone. I am solving an elasticity problem, and I need to construct a 4th-order tensor P and afterward take the dot product with a 2nd-order strain tensor E. The exact operation that I need to perform is,

Tij = Pijkl Ekl

where T is also a second-order tensor. Now here's the catch: the tensor P is not uniform throughout the domain. It varies node-by-node. I've already written the code to compute P's components at each node in NumPy. Let's call the Numpy array p with shape (2, 2, 2, 2, k) where k is the number of nodes, and p stores all the values that need to go in P. My question boils down to this: how can I define P and assign the values of p into it in a consistent manner?

Here is a sample of my code:
...

V = VectorFunctionSpace(mesh, "CG", 1)  # vector space

def strain(u):
"""Small strain operator."""

u = TrialFunction(V)  # trial function for displacement
v = TestFunction(V)  # test function for displacement

P = ???

a = inner(strain(v), dot(P, strain(u)))*dx
L = inner(v, Constant((0, -9.81)))*dx

displacement = Function(V)
solve(a == L, displacement, BC)

...​
Community: FEniCS Project
You can do this using ufl.as_tensor and many sleepless nights wishing you had a smooth field for P rather than definitions at each mesh vertex.
written 11 months ago by Nate
How can I avoid shape mismatches when using ufl.as_tensor, particularly when assigning values into it? I am not very familiar with that function. Thank you.
written 11 months ago by jimmy

0
11 months ago by
As that bloke Nate said, use ufl.as_tensor e.g., here.

If you have numpy arrays, set them locally to Functions (via the Function::vector().set_local(array)) and use these as the components to as_tensor.
0
11 months ago by
import numpy
def VoigtToTensor (A) :
A11 , A12 , A13 , A14 , A15 , A16 = A[ 0 , 0 ] , A[ 0 , 1 ] , A[ 0 , 2 ] , A[ 0 , 3 ] , A[ 0 , 4 ] , A[ 0 , 5 ]
A22 , A23 , A24 , A25 , A26 = A[ 1 , 1 ] , A[ 1 , 2 ] , A[ 1 , 3 ] , A[ 1 , 4 ] , A[ 1 , 5 ]
A33 , A34 , A35 , A36 = A[ 2 , 2 ] , A[ 2 , 3 ] , A[ 2 , 4 ] , A[ 2 , 5 ]
A44 , A45 , A46 = A[ 3 , 3 ] , A[ 3 , 4 ] , A[ 3 , 5 ]
A55 , A56 = A[ 4 , 4 ] , A[ 4 , 5 ]
A66 = A[ 5 , 5 ]
A21 , A31 , A41 , A51 , A61 = A12 , A13 , A14 , A15 , A16
A32 , A42 , A52 , A62 = A23 , A24 , A25 , A26
A43 , A53 , A63 = A34 , A35 , A36
A54 , A64 = A45 , A46
A65 = A56
return as_tensor( [ \
[ \
[ [ A11 , A16 , A15 ] , [ A16 , A12 , A14 ] , [ A15 , A14 , A13 ] ] , \
[ [ A61 , A66 , A65 ] , [ A66 , A62 , A64 ] , [ A65 , A64 , A63 ] ] , \
[ [ A51 , A56 , A55 ] , [ A56 , A52 , A54 ] , [ A55 , A54 , A53 ] ] \
] , [ \
[ [ A61 , A66 , A65 ] , [ A66 , A62 , A64 ] , [ A65 , A64 , A63 ] ] , \
[ [ A21 , A26 , A25 ] , [ A26 , A22 , A24 ] , [ A25 , A24 , A23 ] ] , \
[ [ A41 , A46 , A45 ] , [ A46 , A42 , A44 ] , [ A45 , A44 , A43 ] ] \
] , [ \
[ [ A51 , A56 , A55 ] , [ A56 , A52 , A54 ] , [ A55 , A54 , A53 ] ] , \
[ [ A41 , A46 , A45 ] , [ A46 , A42 , A44 ] , [ A45 , A44 , A43 ] ] , \
[ [ A31 , A36 , A35 ] , [ A36 , A32 , A34 ] , [ A35 , A34 , A33 ] ] ] \
])

Cvoigt = numpy.array( [ \
[ la +2.*mu, la , la , 0 , 0 , 0 ] , \
[ la , la +2.*mu, la , 0 , 0 , 0 ] , \
[ la , la , la +2.*mu, 0 , 0 , 0 ] , \
[ 0 , 0 , 0 , mu, 0 , 0 ] , \
[ 0 , 0 , 0 , 0 , mu, 0 ] , \
[ 0 , 0 , 0 , 0 , 0 , mu ]])
C = VoigtToTensor ( Cvoigt )

nu = 0.3
E = 210e3
G = E/(2.0*(1+nu))
la = 2.0*G*nu/(1.0-2.0*nu)
mu = G