Equivalent strain tensors from displacement field using the Strain at Stations algorithm described in the Cardozo (2009), for use with DEM PFC

strain_from_disp.py

import numpy as np
import pandas as pd
import matplotlib; matplotlib.rcParams["savefig.directory"] = "."
from matplotlib import pyplot as plt
from scipy.spatial import cKDTree
from numpy.linalg import lstsq

data0 = np.loadtxt("S1-1_0m.log")
data1 = np.loadtxt("S1-1_3000m.log")
assert (data0[:,0] == data1[:,0]).all()
assert (data0[:,1] == data1[:,1]).all()

p0 = data0[:,2:]
p1 = data1[:,2:]
radii = data1[:,1]
disp = p0-p1
dx, dy = disp.T

# plt.hist(dx)
# plt.hist(dy)
# plt.show()

# for pp0, pp1 in zip(p0, p1):
#     plt.plot([pp0[0], pp1[0]],[pp0[1], pp1[1]], color="black")
# plt.gca().set_aspect(1)
# plt.show()


# plt.quiver(p0[:,0], p0[:,1], dx, dy)
# plt.gca().set_aspect(1)
# plt.show()


tree = cKDTree(p1)

n = 7 # number of balls
tensors = []
dtensors = []
j2s=[]
for i, p in enumerate(p1):
    _, ind = tree.query(p, n)
    b=[]
    M=[]

    for s in ind:
        b.append(dx[s])
        b.append(dy[s])
        M.append([1,0, p1[s][0], p1[s][1], 0,         0])
        M.append([0,1, 0,        0,         p1[s][0], p1[s][1]])
    a = lstsq(M,b, rcond=None)[0]
    t1, t2, g00, g01, g10, g11 = a
    #eij = 0.5*(gij + gji - (gki gkj))
    e00 = 0.5*(g00 + g00 - (g00*g00 + g10*g10))
    e01 = 0.5*(g01 + g10 - (g00*g01 + g10*g11))
    e10 = 0.5*(g10 + g01 - (g01*g00 + g11*g10))
    e11 = 0.5*(g11 + g11 - (g01*g01 + g11*g11))
    e = [[e00, e01], [e10, e11]]
    tensors.append(e)
    isotropic = (e00+e11)/2.0
    de = [[e00-isotropic, e01], [e10, e11-isotropic]]
    d11_, d22_, d33_, d12_, d23_, d13_ = e00, e11, 0, e01, 0, 0
    a11 = (d11_ - d22_)
    a22 = (d22_ - d33_)
    a33 = (d11_ - d33_)
    j2 =  (a11*a11 + a22*a22 + a33*a33)/6.0 + d12_*d12_ + d23_*d23_ + d13_*d13_
    j2s.append(j2)
    dtensors.append(de)


np.savetxt("j2.txt", j2s)
plt.scatter(p0[:,0], p0[:,1], c=j2s, s=radii**2/100, vmin=0, vmax=1)
plt.gca().set_aspect(1)
plt.colorbar()
plt.show()