利用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)