利用python在abaqus中画Voronoi多面体简单示例
利用scipy.spatial库得到Voronoi多面体顶点坐标
简单例子：
from scipy.spatial import Voronoi
import numpy as np
import csv
aa=np.array([[0,0,0],[10,0,0],[-10,0,0],[0,10,0],[0,-10,0],[0,0,10],[0,0,-10]])
vor = Voronoi(aa)
with open("vertices.csv", "w", encoding="utf-8", newline="") as f:
csv_writer = csv.writer(f)
csv_writer.writerows(vor.vertices)
with open("ridge_vertices.csv", "w", encoding="utf-8", newline="") as f:
csv_writer = csv.writer(f)
for vorRidge in vor.ridge_vertices:
if (vorRidge[0]>=0 and vorRidge[1]>=0 and vorRidge[2]>=0):
csv_writer.writerow(vorRidge)
aa为种子
vor.vertices为Voronoi多面体顶点坐标
vor.ridge_vertices为Voronoi多面体各面上的顶点序号
abaqus中绘制多面体
CAE操作得到相应rpy文件
好习惯先 输入 session.journalOptions.setValues(replayGeometry=COORDINATE,recoverGeometry=COORDINATE)
0、 将vertices.csv和ridge_vertices.csv导入abaqus
import csv
csv_reader1 = csv.reader(open("ridge_vertices.csv"))
ridgeList=[]
for row in csv_reader1:
ridgeList0=[]
for floatRow in row:
ridgeList0=ridgeList0+[float(floatRow),]
ridgeList=ridgeList+[ridgeList0]
csv_reader2 = csv.reader(open("vertices.csv"))
verticesList=[]
for row in csv_reader2:
verticesList0=[]
for floatRow in row:
verticesList0=verticesList0+[float(floatRow),]
verticesList=verticesList+[verticesList0]
1、 新建一个part
在脚本文件中输入p = mdb.models['Model-1'].Part(name='Part-1', dimensionality=THREE_D, type=DEFORMABLE_BODY)
或 CAE 新一个part 后删除几何
2、创建点
从rpy文件读取响应指令
p = mdb.models['Model-1'].parts['Part-1']
p.DatumPointByCoordinate(coords=(0.0, 0.0, 0.0))
修改后：
# create points
p = mdb.models['Model-1'].parts['Part-1']
for vertices in verticesList:
p.DatumPointByCoordinate(coords=(vertices[0], vertices[1], vertices[2]))
3、画线
从rpy文件读取响应指令
p = mdb.models['Model-1'].parts['Part-1']
d = p.datums
p.WirePolyLine(points=((d[0], d[1])), mergeType=IMPRINT, meshable=ON)
修改后：
# create line
p = mdb.models['Model-1'].parts['Part-1']
d = p.datums
for ridge in ridgeList:
for i in range(len(ridge)):
p.WirePolyLine(points=((d[int(ridge[i])+1], d[int(ridge[i-1])+1])), mergeType=IMPRINT, meshable=ON)
4、画面
从rpy文件读取响应指令
p = mdb.models['Model-1'].parts['Part-1']
e = p.edges
p.CoverEdges(edgeList=(e.findAt(coordinates=(5.0, 5.0, 2.5)), e.findAt(
coordinates=(-2.5, 5.0, 5.0)), e.findAt(coordinates=(2.5, 5.0, -5.0)),
e.findAt(coordinates=(-5.0, 5.0, -2.5))), tryAnalytical=False)
修改后：
# create plane
p = mdb.models['Model-1'].parts['Part-1']
e = p.edges
for ridge in ridgeList:
xyzLine=[]
lineList=[]
for i in range(len(ridge)):
xLine=(verticesList[int(ridge[i])][0]+verticesList[int(ridge[i-1])][0])/2.0
yLine=(verticesList[int(ridge[i])][1]+verticesList[int(ridge[i-1])][1])/2.0
zLine=(verticesList[int(ridge[i])][2]+verticesList[int(ridge[i-1])][2])/2.0
xyzLine=xyzLine+[[xLine,yLine,zLine],]
lineList=lineList+[e.findAt(coordinates=[xLine,yLine,zLine]),]
p.CoverEdges(edgeList=lineList, tryAnalytical=False)
完整代码
# -*- coding: mbcs -*-
#
# Abaqus/CAE Release 2020 replay file
# Internal Version: 2019_09_14-01.49.31 163176
# Run by lichen on Tue Dec 19 10:55:39 2023
#
# from driverUtils import executeOnCaeGraphicsStartup
# executeOnCaeGraphicsStartup()
#: Executing "onCaeGraphicsStartup()" in the site directory ...
from abaqus import *
from abaqusConstants import *
session.Viewport(name='Viewport: 1', origin=(0.0, 0.0), width=209.183197021484,
height=206.733337402344)
session.viewports['Viewport: 1'].makeCurrent()
session.viewports['Viewport: 1'].maximize()
from caeModules import *
from driverUtils import executeOnCaeStartup
executeOnCaeStartup()
session.viewports['Viewport: 1'].partDisplay.geometryOptions.setValues(
referenceRepresentation=ON)
cliCommand("""session.journalOptions.setValues(replayGeometry=COORDINATE,recoverGeometry=COORDINATE)""")
import csv
csv_reader1 = csv.reader(open("ridge_vertices.csv"))
ridgeList=[]
for row in csv_reader1:
ridgeList0=[]
for floatRow in row:
ridgeList0=ridgeList0+[float(floatRow),]
ridgeList=ridgeList+[ridgeList0]
csv_reader2 = csv.reader(open("vertices.csv"))
verticesList=[]
for row in csv_reader2:
verticesList0=[]
for floatRow in row:
verticesList0=verticesList0+[float(floatRow),]
verticesList=verticesList+[verticesList0]
p = mdb.models['Model-1'].Part(name='Part-1', dimensionality=THREE_D,
type=DEFORMABLE_BODY)
p = mdb.models['Model-1'].parts['Part-1']
# create points
p = mdb.models['Model-1'].parts['Part-1']
for vertices in verticesList:
p.DatumPointByCoordinate(coords=(vertices[0], vertices[1], vertices[2]))
# create line
p = mdb.models['Model-1'].parts['Part-1']
d = p.datums
for ridge in ridgeList:
for i in range(len(ridge)):
p.WirePolyLine(points=((d[int(ridge[i])+1], d[int(ridge[i-1])+1])), mergeType=IMPRINT, meshable=ON)
# create plane
p = mdb.models['Model-1'].parts['Part-1']
e = p.edges
for ridge in ridgeList:
xyzLine=[]
lineList=[]
for i in range(len(ridge)):
xLine=(verticesList[int(ridge[i])][0]+verticesList[int(ridge[i-1])][0])/2.0
yLine=(verticesList[int(ridge[i])][1]+verticesList[int(ridge[i-1])][1])/2.0
zLine=(verticesList[int(ridge[i])][2]+verticesList[int(ridge[i-1])][2])/2.0
xyzLine=xyzLine+[[xLine,yLine,zLine],]
lineList=lineList+[e.findAt(coordinates=[xLine,yLine,zLine]),]
p.CoverEdges(edgeList=lineList, tryAnalytical=False)