python_flac3d.py


import itasca as it
it.command("python-reset-state false")
it.command("""
program load module 'zone'
program load guimodule 'zone'
program load module 'body'
program load guimodule 'body'
program load module 'extruder'
program load guimodule 'extruder'
program load module 'sel'
program load module 'zonesel'
program load module 'selpython'
program load module 'wallsel'
program load guimodule 'sel'
program load module 'pfcsel'
program load module 'rblocksel'
program load module 'dfnzone'
program load module 'ballzone'
program load module 'zonepython'
program load module 'wallzone'
""")



it.command("""
model new
model large-strain off
zone create brick size 10 10 10
zone cmodel assign elastic
zone property density 2950 young 12e9 poisson 0.25
cycle 1
""")



it.zone.count()



z = it.zone.find(1)
print(z)



z.pos()



print("zone", z.id(), "with model", z.model(), "at", z.pos())



volume_sum = 0.0
for z in it.zone.list():
    volume_sum += z.vol()



print(volume_sum)
print(z.vol() * it.zone.count())
assert volume_sum == z.vol() * it.zone.count()



z = it.zone.near((5,5,5))



z.pos()



z.props()



z.props()['bulk']




z.prop('shear')



z.set_prop('bulk', 8.5e9)
print(z.prop('bulk'))




gp = it.gridpoint.near((2,2,2))
print(gp.pos())



total_mass = 0
for gp in it.gridpoint.list():
    total_mass += gp.mass_gravity()

print(total_mass, z.vol()*1000*z.density())



it.command("""
struct cable create by-line (0.0, 5, 5) (10, 5, 5) id=1 segments=4
zone mechanical damping combined
struct cable property cross-sectional-area=2e-3 young=200e9 ...
                      yield-tension=1.1e20 grout-stiffness=1e10 ...
                      grout-cohesion=1e20 grout-friction=0.123
struct beam create by-line (0.0, 2, 2) (10, 2, 2) id=2 segments=4
struct beam property young=1.234e9 poisson=0.223 plastic-moment=23.4e4 ...
                     cross-sectional-area=1.234 moi-y=55.66 moi-z=3458.33
""")



for s in it.structure.list():
    if type(s) is it.structure.Cable:
        print("Cable element with id: {} and component id: {}"
              .format(s.id(), s.component_id()))
    elif type(s) is it.structure.Beam:
        print("Beam element with Young's modulus: {}".format(s.young()))



s1 = it.structure.find(1)
s2 = it.structure.near((0,2,2))

print(s1)
print(s2)



s_node = it.structure.node.find(1)
s_link = s_node.links()[0]
print("SEL node {} is attached to zone {}"
      .format(s_node.id(), s_link.target().id()))

it.command("struct cable delete")
it.command("struct beam delete")



z.set_extra(1, 1.23)
z.set_extra(2, "a text string")

print(z.extra(1))
print(z.extra(2))



gp.set_extra(1, gp.pos())
gp.set_extra(2, 86856527)

print(gp.extra(1))
print(gp.extra(2))



it.command("""
zone group "lower" range position-z 0 5
zone group "upper" range position-z 5 10
""")




lower_zones, upper_zones = 0, 0

for z in it.zone.list():
    if z.group("default") == "lower":
        lower_zones += 1
        z.set_prop("young", 1.6e10)
    elif z.group("default") == "upper":
        upper_zones += 1

print(lower_zones, upper_zones)



lower_gridpoints, upper_gridpoints = 0, 0

for gp in it.gridpoint.list():
    if gp.pos_z() == 0:
        lower_gridpoints += 1
        gp.set_fix(0, True)
        gp.set_fix(1, True)
        gp.set_fix(2, True)
    if gp.pos_z() == 10:
        upper_gridpoints += 1
        gp.set_force_load((1e6,2e6,1e6))

print(lower_gridpoints, upper_gridpoints)



it.command("""
model save "before_cycling"
model solve
""")


gp_top_middle = it.gridpoint.near((5,5,10))
gp_middle = it.gridpoint.near((5,5,5))
print(gp_top_middle.disp())
print(gp_middle.disp())




central_zone = it.zone.near((5,5,5))
stress = central_zone.stress()
strain = central_zone.strain()

print(stress)
print()
print(strain)



for modulus in [6e9, 8e9, 10e9, 12e9]:
    it.command("""
    model restore 'before_cycling'
    zone prop young {}
    model solve
    """.format(modulus))
    vertical_disp = it.gridpoint.near((5,5,10)).disp_z()
    print("Modulus: {:.01e} Pa, vertical displacement: {:.03e}"
          .format(modulus, vertical_disp))




try:
    it.command("""
    zone prop young 0
    model solve
    """)
except Exception as data:
    print("an error occurred", data)

print("Python program execution continues...")

it.command("zone prop young 8e9")  # fix the problem



z = it.zone.near((5,5,5))
len(z.gridpoints())



for gp in z.gridpoints():
    print("gridpoint with id: {} at {}".format(gp.id(), gp.pos()))



z = it.zone.near((5,5,5))
print("central zone id: {}, position: {}".format(z.id(), z.pos()))
print()
for i, neighbor in enumerate(z.adjacent_zones()):
    print("neighbor zone {} id: {}, position: {}"
          .format(i, neighbor.id(), neighbor.pos()))



z = it.zone.near((0,0,0))
for neighbor in  z.adjacent_zones():
    print(neighbor)



for i,face in enumerate(z.faces()):
    print("face {} contains the gridpoints with IDs: {}"
          .format(i, [gp.id() for gp in face]))



for i,normal in enumerate(z.face_normals()):
    print("face {} has normal vector {}".format(i, normal))




gp = it.gridpoint.near((2,2,2))
print(gp.pos())
for z in gp.zones():
    print("gp id: {} is connected to zone id: {}, position: {}"
          .format(gp.id(), z.id(), z.pos()))



print(type(1.0) is int)

print(type(1.0) is float)

print(type("Itasca") is str)

print(type(gp) is it.gridpoint.Gridpoint)

print(type(gp) is it.zone.Zone)

print(type(z) is it.gridpoint.Gridpoint)

print(type(z) is it.zone.Zone)


i=0
def my_callback(*args):
    global i
    i += 1
    print("in Python callback function.")

it.set_callback("my_callback", -1)
it.command("cycle 5")
print("The Python callback function was called {} times".format(i))

it.remove_callback("my_callback",-1)
i=0
it.command("cycle 5")
print("The Python callback function was called {} times".format(i))