0

I am processing symmetric second order tensors (of stress) using numpy. In order to transform the tensors I have to generate a fully populated tensor, do the transformation and then recover the symmetric tensor in the rotated frame.

My input is a 2D numpy array of symmetric tensors (nx6). The code below works, but I'm pretty sure there must be a more efficient and/or elegant way to manipulate the arrays but I can't seem to figure it out.

I anyone can anyone suggest an improvement I'd be very grateful? The sample input is just 2 symmetric tensors but in use this could be millions of tensors, hence the concernr with efficiency

Thanks,

Doug

# Sample symmetric input (S11, S22, S33, S12, S23, S13)
sym_tens_in=np.array([[0,9], [1,10], [2,11], [3,12], [4,13], [5,14]])
   
# Expand to full tensor
tens_full=np.array([[sym_tens_in[0], sym_tens_in[3], sym_tens_in[4]],
                    [sym_tens_in[3], sym_tens_in[1], sym_tens_in[5]],
                    [sym_tens_in[4], sym_tens_in[5], sym_tens_in[2]]])

# Transpose and reshape to n x 3 x 3 
tens_full=np.transpose(tens_full, axes=(2, 0, 1))

# This where the work on the full tensor will go....

# Reshape for extraction of the symmetric tensor
tens_full=np.reshape(tens_full, (2,9))

# Create an array for the test ouput symmetric tensor
sym_tens_out=np.empty((2,6), dtype=np.int32)

# Extract the symmetric components
sym_tens_out[:,0]=tens_full[:,0]
sym_tens_out[:,1]=tens_full[:,4]
sym_tens_out[:,2]=tens_full[:,8]
sym_tens_out[:,3]=tens_full[:,2]
sym_tens_out[:,4]=tens_full[:,3]
sym_tens_out[:,5]=tens_full[:,5]

# Transpose....
sym_tens_out=np.transpose(sym_tens_out)
scotsman60
  • 251
  • 1
  • 10
  • A `scipy` module has functions to convert `vector` form to `squareform`. I haven't paid much attention to it, so can't say whether the mapping is the same as yours. https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.squareform.html#scipy.spatial.distance.squareform – hpaulj Dec 18 '22 at 22:30
  • The only reason I need to do any more than expand the symmetric tensor is because I need to transform each tensor using numpy.dot and this function needs the 3x3s to be in the second two axes and not the first two.... – scotsman60 Dec 19 '22 at 17:01

2 Answers2

0

This won't be any faster, but it's more compact:

In [166]: idx=np.array([0,3,4,3,1,5,4,5,2]).reshape(3,3)
In [167]: sym_tens_in[idx].transpose(2,0,1)
Out[167]: 
array([[[ 0,  3,  4],
        [ 3,  1,  5],
        [ 4,  5,  2]],

       [[ 9, 12, 13],
        [12, 10, 14],
        [13, 14, 11]]])

The transpose could be done first:

sym_tens_in.T[:,idx]

Similarly the reverse mapping can be done with:

In [168]: idx1 = [0,4,8,1,2,5]
In [171]: tens_full.reshape(2,-1)[:,idx1]
Out[171]: 
array([[ 0,  1,  2,  3,  4,  5],
       [ 9, 10, 11, 12, 13, 14]])

with the optional transpose.

hpaulj
  • 221,503
  • 14
  • 230
  • 353
0

OK - Based on the answers provided here I found a really cool solution. Now, I have to say that in my original question I omitted the actual reason that I was trying to get the full tensor into nx3x3 form. Basically, I'm implementing a function to rotate 2nd order stress tensors which requires solution of σ′=R⋅σ⋅RT. I was planning to use numpy.matmul for the matrix multiplication but to transform multiple stress tensors, matmul requires the 3x3 tensors to be in the last two indices of the nx3x3 matrix - hence the effort to get the data into nx3x3 from from the original 3x3xn form....

However, after I let go of numpy.matmul as my target solution and embraced numpy.einsum instead....... everything became much easier....

# Sample symmetric input (S11, S22, S33, S12, S23, S13)
sym_tens_in=np.array([[0,9], [1,10], [2,11], [3,12], [4,13], [5,14]])
idx=np.array([0,3,5,3,1,4,5,4,2]).reshape(3,3)
full=sym_tens_in[idx]
full_transformed=np.einsum('ij, jkn, lk->nil', rot_mat, full, rot_mat)

Thanks for the inspiration!!!!

scotsman60
  • 251
  • 1
  • 10