import functools
import numpy as np
np.set_printoptions(threshold=20)
import itasca as it
from itasca import zonearray as za
from itasca import gridpointarray as gpa
it.command("""
model reset
model large-strain off
zone create brick size 10 10 10
zone cmodel assign elastic
zone property density 2950 young 12e9 poisson 0.25
model cycle 1
""")
za.pos()
p = za.pos()
type(p)
p.shape
it.zone.count()
x,y,z = p.T
print(x)
print(x.shape)
z > 5
upper_zones = z > 5
print((upper_zones).sum())
modulus = np.zeros_like(x)
modulus[upper_zones] = 1.2e10
lower_zones = np.logical_not(upper_zones)
modulus[lower_zones] = 1.6e10
za.set_prop_scalar("young", modulus)
it.command("zone group \"lower\" range position-z 0 5")
za.in_group("lower")
za.in_group("lower").sum(), "zones in lower group."
from functools import reduce
corner_mask = reduce(np.logical_and, (x<3, y<3, z<3))
za.set_group(corner_mask, "corner", "geometry")
print(za.in_group("corner", "geometry").sum(), "zones in corner group.")
gpos = gpa.pos()
gx, gy, gz = gpos.T
print(gz)
f = gpa.fixity()
print(f)
f[:,][gz==0] = True, True, True
print(f)
gpa.set_fixity(f)
top_gridpoints = gz==10
radial_distance = np.sqrt((gx-5)**2+(gy-5)**2)
central_gridpoints = radial_distance < 5
mask = np.logical_and(top_gridpoints, central_gridpoints)
print("boundary load applied to {} gridpoints".format(mask.sum()))
fapp = gpa.force_app()
print(fapp)
fapp[:,2] = mask*1e6*(5.0-radial_distance)/5.0
gpa.set_force_app(fapp)
it.command("model solve")
print("gridpoint displacements:")
print(gpa.disp())
print("gridpoint displacement magnitudes: ")
mag = np.linalg.norm(gpa.disp(), axis=1)
print(mag)
max_index = np.argmax(mag)
print("Maximum displacement: {} at location {}".format(gpa.disp()[max_index],
gpa.pos()[max_index]))
print("Vertical displacement along vertical line x=5, y=5: from z=0 to z=10")
print(gpa.disp()[np.logical_and(gx==5, gy==5)][:,2])
za.stress()
za.stress_flat()
za.ids()
gpa.ids()
za.gridpoints()
print("average gridpoint locations: ")
print(gpa.pos()[za.gridpoints()].sum(axis=1)/8.0)
print("zone centroids: ")
print(za.pos())
za.gridpoints()
za.faces()
print(za.neighbors())
print()
exterior_zone_mask = [(v==-1).any() for v in za.neighbors()]
n_exterior = sum(exterior_zone_mask)
print(n_exterior, "zones on the model exterior")
N = 10
assert n_exterior == 2*N*N + 2*(N-2)*N + 2*(N-2)*(N-2)
