【分享】DFN程序-python库

说明

本程序基于开源程序修改，检查了源代码中存在的一些错误，做了修改，编译成window系统下可用，支持python3.10.x。使用者请自行验证，no warranty！来源论文:《UnblockGen – A Python library for 3D rock mass generation and analysis》

一个开源的DFN生成程序–UnblockGen

在源代码基础上添加了：裂隙面属性输出,法向量、中心点坐标、顶点，倾向/倾角（可以使用mlsteronet绘制裂隙倾向/倾角分布结果）添加了椭圆形裂隙面的生成倾向/倾角基于高斯分布的椭圆形裂隙面生成添加获取裂隙面数量、块体数量的函数，实际上程序也有返回相应列表的功能，使用python中的len()也能获取添加获取生成的block的id、顶点、面、边、中心点、几何体Polygon等几何参数的函数添加了生成一个指定参数（法向量、中心点、宽度、高度）的平面计算裂隙面与目标平面的交点get_IntersectLineToPlane，获取目标平面与裂隙面的交线，目前仅使用圆形裂隙面

生成结果

程序生成的dfn为vtk格式，可使用paraview进行浏览，当然你也可以自己写一个可视化界面，将源代码编译成库进行使用。

扩展

如果你熟悉DDA,获取可以基于该程序进行扩展，它已提供了基本的几何信息，如何让块体的”动起来”，需要你的贡献。DDA原理：https://www.researchgate.net/publication/355696579_Discontinuous_Deformation_Analysis_in_Rock_Mechanics_and_Rock_Engineering

基于物理引擎库的仿真，例如Bullet(pyBullet)、web端的canonjs等，做一些动画模拟，隧洞、边坡滑塌，期待你的贡献。

示例程序

from unblocks import

 *

import plotTools as

 plotTools

import matplotlib.pyplot as

 plt

import numpy as

 np

import

 mplstereonet

#定义裂隙网络模型

dfn = DFN()

dfn.set_RandomSeed(100)#设置随机数-则相同随机数下生成的模型是相同的dfn.set_RegionMaxCorner([15,15,150

])

dfn.add_FractureSet()#添加裂隙集1dfn.add_FractureSet()#添加裂隙集2#dfn.add_FractureSet()#添加裂隙集3dfn.add_LineMapping([10,30,3], [0.766582,4.631392,24.30794]) #线映射#dfn.add_QuadrilateralMapping([0,0,20],[100,0,20],[100,100,20],[0,100,20])#面映射dfn.add_VolumeMapping()#体积映射#生成裂隙集while (dfn.linesMapping[0].get_P10(0) < 0.05

):

 dfn.fractureSets[0].add_BaecherFracture(110,48,50000,“det”,200,0

)

# while (dfn.surfacesMapping[0].get_P21(0) < 0.2):#  dfn.fractureSets[1].add_BaecherFracture(90, 45, 20, “exp”, 40, 0);while (dfn.volumesMapping[0].get_P30(1) < 0.0001

):

 dfn.fractureSets[1].add_BaecherFracture(70, 45, 100, “exp”, 5, 0

);

while (dfn.volumesMapping[0].get_P32(1) < 0.002

):

 dfn.fractureSets[1].add_BaecherFracture(180, 65, 100, “exp”, 4, 1

);

#向裂隙面集合1添加椭圆形裂隙面fracCenter=[7.5,7.5,10

]

fracDirection=175fracdipAngle=90fracRadius=5dfn.fractureSets[0].add_EllipseFracture(fracCenter,fracDirection,fracdipAngle,fracRadius,fracRadius*2

)

#添加一个倾向/倾角高斯分布的椭圆形裂隙面meanDipDirection=75meanDipAngle=60sigmaDipDirection=5sigmaDipAngle=3sizeDistribution=exp #det 、logmeanFractureSize=5sigmaFractureSize=2for i  in range(10

):

 dfn.fractureSets[0

].add_GaussDistFracture(meanDipDirection, meanDipAngle, sigmaDipDirection, sigmaDipAngle, sizeDistribution, meanFractureSize, sigmaFractureSize)

#导出区域dfn.export_RegionVtk(“./cube/region”

)

#获取裂隙集数量print(“FractureSetNum:”,dfn.get_FractureSetsNum(),dfn.FractureSetsNum)#两个一样tpCenter=[15,15,75

]

tpNormal=[0,0,1

]

dipDirectionsAndAngle=[]

#遍历每个生成的裂隙集及其裂隙面for i in

 range(dfn.FractureSetsNum):

 fs=dfn.fractureSets[i]

 n=fs.get_FractureNum() #集合中的裂隙数量 for k in

 range(n):

  f=fs.fractures[k]#获取裂隙面  #打印裂隙面的属性  print(“id”,f.id,“Center:”,f.get_Center(),“normal:”,f.get_UnitVector(),“radius:”,f.get_BoundingSphereRadius(),“area:”

,f.get_Area())

  #print(“direction/angle”,f.get_DipDirectionDipAngle())

  dipDirectionsAndAngle.append(f.get_DipDirectionDipAngle())

  #计算裂隙面与目标平面的交点

  intersection_points=f.get_IntersectLineToPlane(tpCenter,tpNormal)

  print(“intersection_points”

,intersection_points)

dfn.export_DFNVtk(“./cube/dfn”

)

#块体生成

generator = Generator()

#参数设置generator.set_MinInscribedSphereRadius(0.05

) 

generator.set_MaxAspectRatio(30

)

#执行生成

generator.generate_RockMass(dfn)

#生成一个平面，输入平面中心点、法向量、宽度、高度—作用可自行探索generator.generate_PlaneVTK(“./cube/plane1”,tpCenter,tpNormal,50,50

)

#访问块体信息#块体数量

blocksNum=generator.get_BlocksNum()

for i in

 range(blocksNum):

 block0=generator.blocks[i]

 print(“#”*20,“block-“,i,“块体信息”,“#”*20

)

 print(“Id:”,block0.get_Id(),“Volume:”

,block0.get_Volume())

 print(“Alpha:”,block0.get_Alpha(),“Beta:”,block0.get_Beta(),“AspectRatio:”

,block0.get_AspectRatio())

 print(“*”*20,“block-“,i,“几何信息”,“*”*20

)

 print(“Center:”

,block0.get_Center())

 print(“Vetices:”

,block0.get_Vertices())

 print(“Edges:”

,block0.get_Edges())

 print(“Planes:”

,block0.get_Planes())

 print(“Polygons:”

,block0.get_Polygons())

 block0.export_BlockVtk(“./cube/blocks/b”+str(i))#导出单个块体generator.export_BlocksVtk(“./cube/block”

)

dipDirectionsAndAngle=np.array(dipDirectionsAndAngle)

#绘制块体统计结果#统计plotTools.blockVolumeDistribution(generator.get_Volumes(True

))

plotTools.blockShapeDiagram(generator.get_AlphaValues(True), generator.get_BetaValues(True), generator.get_Volumes(True), 0.05

)

plotTools.BlockShapeDistribution(generator.get_AlphaValues(True), generator.get_BetaValues(True), generator.get_Volumes(True

))

plotTools.showPlots()

#Stero统计

fig, ax = mplstereonet.subplots()

cax = ax.density_contourf(dipDirectionsAndAngle[:,0], dipDirectionsAndAngle[:,1], measurement=poles

)

ax.pole(dipDirectionsAndAngle[:,0], dipDirectionsAndAngle[:,1

])

ax.grid(True

)

fig.colorbar(cax)

plt.show()

stereonet

私信回复UBG获取下载链接