Construction of a 4th-order tensor in FEniCS


237
views
2
3 months ago by
jimmy  
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."""
    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
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 3 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 3 months ago by jimmy  

2 Answers


0
3 months ago by
pf4d  
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
3 months ago by
hsk  
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.

Similar posts:
Search »