flac3d_array_interface.py


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)