Skip to content

Instantly share code, notes, and snippets.

@jkfurtney
Created May 25, 2021 15:29
Show Gist options
  • Save jkfurtney/7da203f7dd89dd55a5c0e3486e84ee4f to your computer and use it in GitHub Desktop.
Save jkfurtney/7da203f7dd89dd55a5c0e3486e84ee4f to your computer and use it in GitHub Desktop.
Equivalent strain tensors from displacement field using the Strain at Stations algorithm described in the Cardozo (2009), for use with DEM PFC
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment