Flac3D-断层带-目标点【水平、垂直应力、正应力、剪应力、剪应力和正应力比值】

import itasca as it
from math import *
it.command("python-reset-state false")
"""
1、如果只为了求解断层的正应力和剪应力比值，只需要更改zone_id即可
2、如果需要正应力和剪应力，文件夹还需要改一改
"""
#zone的id 列标配
zone_id = [3689881,3933005,3626246,3726742,3446164]

#循环每一个模型
for i in range(1,15):
    it.command(
    """
    model restore 'kaiwa-{index}.sav'
    """.format(index=i)
    )
    #循环每一个监测点
    for i in zone_id:
        #通过ID查找对应的zone
        zone = it.zone.find(i)
        #水平和垂直的应力
        sxx = zone.stress()[0][0]
        szz = zone.stress()[2][2]
        
        #正应力和剪应力求解
        zyl = (szz+sxx) * 0.5 - (sxx-szz) * cos(100*pi/180) * 0.5
        jyl = (sxx-szz) * sin(100*pi/180) * (-0.5)
        #比值
        bz = round(jyl / zyl, 3)
        #比值写入
        with open("targetZoneBz.txt", "a") as f:
            f.write(str(bz) + " ")
    with open("targetZoneBz.txt", "a") as f1:
        f1.write("\n")
        
        #正应力和剪应力写入
		with open("zyl.txt", "a") as f1:
            f1.write(str(zyl) + " ")
        with open("jyl.txt", "a") as f2:
            f2.write(str(jyl) + " ")
    with open("zyl.txt", "a") as f1_1:
        f1_1.write("\n")
    with open("jyl.txt", "a") as f2_1:
        f2_1.write("\n")