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

588

views

2

Hello everyone. I am solving an elasticity problem, and I need to construct a 4th-order tensor

where

Here is a sample of my code:

How can I avoid shape mismatches when using

*and afterward take the dot product with a 2nd-order strain tensor***P***. The exact operation that I need to perform is,***E***T*_{ij}=*P*_{ijkl}*E*_{kl}where

*is also a second-order tensor. Now here's the catch: the tensor***T***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**`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 *. My question boils down to this: how can I define***P***and assign the values of***P**`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."""
return sym(nabla_grad(u))
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

`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

### 2 Answers

0

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.

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

```
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
```

Please login to add an answer/comment or follow this question.

`ufl.as_tensor`

and many sleepless nights wishing you had a smooth field for P rather than definitions at each mesh vertex.