In [1]:
%load_ext autoreload
%autoreload 2
In [2]:
import numpy as np
from scipy.sparse.linalg import cg, lsqr, gmres
import matplotlib.pyplot as plt
import matplotlib.tri
import okada_wrapper
import tectosaur as tct
In [28]:
TCTN = 121
CURVE = 1.0
corners = [[-1, 0, -1], [1, 0, -1], [1, 0, 1], [-1, 0, 1]]
src_mesh = tct.make_rect(TCTN, TCTN, corners)
src_mesh[0][:,1] = CURVE * np.abs(src_mesh[0][:,0])
#src_mesh[0][:,1] = CURVE * (src_mesh[0][:,0] ** 2 + src_mesh[0][:,2] ** 2)
def gauss_slip_fnc(x, z):
r2 = x ** 2 + z ** 2
R = 0.6
out = (np.cos(np.sqrt(r2) * np.pi / R) + 1) / 2.0
out[np.sqrt(r2) > R] = 0.0
return out
dof_pts = src_mesh[0][src_mesh[1]]
x = dof_pts[:,:,0]
z = dof_pts[:,:,2]
slip = np.zeros((src_mesh[1].shape[0], 3, 3)).astype(np.float32)
slip[:,:,0] = gauss_slip_fnc(x, z)
In [29]:
plt.tripcolor(src_mesh[0][:,0], src_mesh[0][:,2], src_mesh[1], src_mesh[0][:,1])
plt.show()
In [30]:
H = tct.RegularizedSparseIntegralOp(
6, 6, 6, 2, 5, 2.5,
'elasticRH3', 'elasticRH3', [1.0, 0.25], src_mesh[0], src_mesh[1], np.float32,
farfield_op_type = tct.TriToTriDirectFarfieldOp
)
traction_mass_op = tct.MassOp(4, src_mesh[0], src_mesh[1])
In [31]:
dof_pts = src_mesh[0][src_mesh[1]]
from tectosaur.util.geometry import unscaled_normals
ns = unscaled_normals(dof_pts)
ns /= np.linalg.norm(ns, axis = 1)[:, np.newaxis]
ns = np.repeat(ns, 3, axis = 0)
edge_dofs = tct.free_edge_dofs(src_mesh[1], tct.find_free_edges(src_mesh[1]))
obs_pts = dof_pts.reshape((-1,3))[edge_dofs]
obs_ns = ns[edge_dofs]
interior_op = tct.InteriorOp(obs_pts.copy(), obs_ns.copy(), src_mesh, 'elasticH3', 4.0, 11, 7, 7, [1.0, 0.25], np.float32)
C_ik = 0.5 * np.ones
trac_edges = -interior_op.dot(slip.flatten())# / 4.0
vec_edge_dofs = np.tile(3 * np.array(edge_dofs)[:,np.newaxis], (1,3))
vec_edge_dofs[:,1] += 1
vec_edge_dofs[:,2] += 2
vec_edge_dofs = vec_edge_dofs.flatten()
edge_cs = tct.simple_constraints(vec_edge_dofs, trac_edges)
In [32]:
cs = tct.continuity_constraints(src_mesh[0], src_mesh[1], src_mesh[1].shape[0])
cs.extend(edge_cs)
cmU, _, _ = tct.build_constraint_matrix(cs, src_mesh[1].shape[0] * 9)
In [33]:
csT1, csT2 = tct.traction_admissibility_constraints(src_mesh[0], src_mesh[1])
csT1.extend(edge_cs)
#csT1.extend(csT2)
cmT, c_rhs, _ = tct.build_constraint_matrix(csT1, src_mesh[1].shape[0] * 9)
import scipy.sparse
rows = []
cols = []
data = []
cmT2_rhs = np.zeros((len(csT2)))
for i in range(len(csT2)):
c = csT2[i]
for j in range(len(c.terms)):
rows.append(i)
cols.append(c.terms[j].dof)
data.append(c.terms[j].val)
cmT2_rhs[i] = c.rhs
cmT2 = scipy.sparse.csr_matrix((data, (rows, cols)), shape = (len(csT2), cmT.shape[0]))
constrained_traction_mass_op = cmU.T.dot(traction_mass_op.mat.dot(cmT))
full_lhs = scipy.sparse.vstack((constrained_traction_mass_op, cmT2.dot(cmT)))
rhs = -cmU.T.dot(H.dot(slip.flatten())) - cmU.T.dot(traction_mass_op.mat.dot(c_rhs))
full_rhs = np.concatenate((rhs, cmT2_rhs - cmT2.dot(c_rhs)))
# full_lhs.shape, full_rhs.shape, cmT.shape, cmU.shape, cmT2.shape, constrained_traction_mass_op.shape
#soln = cg(constrained_traction_mass_op, rhs)
soln = scipy.sparse.linalg.lsmr(full_lhs, full_rhs, show = True, atol = 1e-8)
out = cmT.dot(soln[0]) + c_rhs
pts = src_mesh[0][src_mesh[1]].reshape((-1,3))
idxs = np.where(np.abs(pts[:,2] - 0.115) < 0.005)[0]
#idxs = np.where(np.abs(pts[:,0]) < 0.01)[0]
plt.figure(figsize = (18, 6))
for d in range(3):
plt.subplot(1,3,d+1)
plt.plot(pts[idxs][:,0], out.reshape((-1,3))[idxs][:,d], '-*')
plt.show()
pts = src_mesh[0][src_mesh[1]].reshape((-1,3))
idxs = np.where(np.abs(pts[:,0]) < 0.01)[0]
#idxs = np.where(np.abs(pts[:,0]) < 0.01)[0]
plt.figure(figsize = (18, 6))
for d in range(3):
plt.subplot(1,3,d+1)
plt.plot(pts[idxs][:,2], out.reshape((-1,3))[idxs][:,d], '.')
plt.show()
In [34]:
def plotter(fields):
levels = []
for f in fields:
minf = np.min(f)
maxf = np.max(f)
levels.append(np.linspace(minf - 1e-12, maxf + 1e-12, 21))
#levels[2] = np.linspace(-0.15, 0.15, 11)
# levels = [
# np.linspace(-0.1, 0.6, 21),
# np.linspace(-1.9, 1.9, 21),
# np.linspace(-0.2, 0.2, 21)
# ]
tri_centers = np.mean(src_mesh[0][src_mesh[1]], axis = 1)
side = np.zeros(src_mesh[1].shape[0]).astype(np.int)
side[tri_centers[:,0] > 0] = 1
n_sides = 2
cmap = 'RdBu_r'
plt.figure(figsize = (16, 6))
for d in range(len(fields)):
plt.subplot(1,len(fields),d+1)
field = fields[d]
for i in range(n_sides):
which_tris = np.where(side == i)[0]
reduced_m = tct.mesh.modify.remove_unused_pts((src_mesh[0], src_mesh[1][which_tris]))
soln_vals = np.empty(reduced_m[0].shape[0])
soln_vals[reduced_m[1]] = field[which_tris]
cntf = plt.tricontourf(
reduced_m[0][:,0], reduced_m[0][:,2], reduced_m[1], soln_vals,
cmap = cmap, linewidth = 0.5, levels = levels[d], extend = 'both'
)
if levels[d].shape[0] < 25:
plt.tricontour(
reduced_m[0][:,0], reduced_m[0][:,2], reduced_m[1], soln_vals,
linestyles = 'solid', colors='k', linewidths = 0.7, levels = levels[d]
)
plt.colorbar(cntf)
plt.tight_layout()
plt.show()
zvec = np.cross(ns, [0,1,0])
xvec = np.cross(zvec, ns)
normal_trac = np.sum(out.reshape((-1,3)) * ns, axis = 1).reshape((-1,3))
horizontal_shear = np.sum(out.reshape((-1,3)) * xvec, axis = 1).reshape((-1,3))
vertical_shear = np.sum(out.reshape((-1,3)) * zvec, axis = 1).reshape((-1,3))
# plotter([
# horizontal_shear,
# normal_trac,
# vertical_shear
# ])
plotter([
out.reshape((-1,3,3))[:,:,d]
for d in range(3)
])
In [35]:
cmT.shape, cmT2.shape, cmU.shape, cmT.shape[1] - cmT2.shape[0]
Out[35]: