Abaqus复杂模型六面体网格划分的Python脚本方法

来了，来了。先冒个泡，本来说写骨料的，但是那个有点难，我一时半会也搞不定，我寻思着，拖一下吧，你们肯定觉得我又没了，不太好，还是按时冒个泡，骨料的我已经在写了，写好就发。

今天写个稍微简单一点的，关于如何实现不规则几何的六面体网格画法，当然，这个不是我们常规理解的六面体网格画法，那种画复杂模型不现实，我也不会，而是像素化处理，就和玩我的世界里面的小方格一样，但是不是在3090下玩的，而是在你我的650玩的，开玩笑，在我的650，你们肯定不是。估计都上了3系的，有没有好东西的，告诉我，让我羡慕羡慕。

我记得我之前写过，类似的，那个好像是二维的Voronoi图的四边形网格，可以回顾一下前面的文章。可以对比一下，看我以前写的，与我今天写的，有没有什么不同，看我写码的能力有没有进步。那话不多说，就开始，今天的旅途吧。（我每次都是重新写，基本不看我以前写的，所以，你们帮我比比看）

+ Abaqus二维多晶体的四边形网格处理脚本

今天的目标，就是把一个钻头，应该是钻头吧。如下图，做成下下图的像素单元的形式。就是我们的目标。懂了要干什么了对吧。就是做成这种像素网格的感觉。这个模型是我找一位粉丝要的，在这里先谢过，刚好他也有需求，我也可以发点东西，算是互相帮助了，我喜欢这种简单且纯粹的合作，别让我做太难的，难的在下不会。

当时我说骨料下版让你们写，我这把文献都贴了，算法都有了，这都一年了，还在等我呢，是不是你们偷偷写好了，藏私不想拿出来，让我一个人在这里孤独落发。

开始之前，我想先说说，这种网格划分的特点，(因为这个比较简单，脚本我一说你们就会)，先浪费点时间说说别的，两种网格划分方式的对比图如下。

这种像素画的网格划分方式最大的优点，就是它可以使用六面体网格，最大的弊端，就是它的实际几何边界已经缺失了。你们会发现，这左边的这种不也是六面体网格吗，边界还顺滑，这多次一举，在这里费力不讨好，意义何在？确实，这个模型好在可以用扫掠画六面体，这里确实意义不大，但是，总有模型是复杂到画不了六面体网格的， 比如那个三维Voronoi图的，如果你想画六面体，这意义不就来了吗。是吧

这会就来说说，怎么做。先说下算法，然后再说怎么实现。

１.首先要有两个part，一个part是几何体，一个part是能覆盖到整个几何的基体，形状倒是无所谓，球的方的，都行。

２.把基体的网格画好。这个基体的网格，就是最后成品的网格大小，所以好好画。

3.遍历基体part所有的网格，计算网格中心坐标，然后判断这个中心坐标是不是在参考part(钻头)内部，如果是，把这个单元的编号存起来。最后，建立set，或者别的处理方法，或者把不在内部的单元删掉，最后怎么处理，就看个人的需求了。

4.以上就是整个实现算法的核心内容。

代码实现

下面可能就是你们感兴趣的部分了，怎么把上面的核心算法，写成python脚本。我把解释写里面，留一份未注释版本供各位复制，因为我加了中文注释abaqus就不认了，不知道你们是不是，我也不知道发生了什么。

#前处理无脑导入一下包就完事了
from abaqus import *
from abaqusConstants import *
#我需要用到numpy,拿过来用一下
import numpy as np

#定义一下model
model = mdb.models["Model-1"]
#partGear就是参考的钻头几何模型
#partBase就是几何
partGear = model.parts["Part-1"]
partBase = model.parts["Part-2"]

#存一下变量
elements = partBase.elements
nodes = partBase.nodes

#labels数组存一下在钻头内部的单元编号
labels = []
cells = partGear.cells

#这个循环，就是核心部分了
#遍历每个element
for element in elements:
    #connectivity存的是组成该单元的节点索引值，
    #注意是索引值，不是单元编号
    nodeIndex = element.connectivity
    center = np.array([0.0, 0.0, 0.0])
    #这里的循环，是把单元的8个节点坐标加起来在除8，
    #就可以计算出单元中心点坐标了
    for index in nodeIndex:
        center += nodes[index].coordinates
    center /= 8
    #这一步是未了判断单元中心点是否在参考几何的内部，
    #利用的是findAt()函数，这个函数的所有使用方法，
    #可以在abaqus帮助文档查到
    #如果findAt()找到了，会返回一个cell序列值，如果没找到，就返回一个空值
    findCell = cells.findAt((center,),printWarning=False)
    #所以，就可以下定义了，如果findCell不是空的，就表明这个单元点
    #在几何内部，把单元编号存起来
    if len(findCell):
        labels.append(element.label)

#这里主要是把内部单元存在一个set里，当然，也可以选择删除其他的单元，
#只需要最后这里改改就好了
#大家可以学学用用sequenceFromLabels()，根据单元编号找单元的序列
#这个函数非常好用，会经常用到，记住它
partBase.Set(elements=elements.sequenceFromLabels(labels=labels), name="Set-Test")

结果如下：

大体上还是对的，但是有两个瑕疵，细心的小伙伴，可能再一开始就发现了，有些单元没检测的有问题，这个我没检查原因，感觉我写的没啥问题。如果遇到这种，问题不明显的，手动去除一下就行。第二个弊端是，运行的有点慢，findAt()函数比较费时，遇到这种大基数遍历的，时间上难以招架，需要耐心等待。

没注释的我也放下面。模型，我就不放了吧，大家自己可以随便用模型测试，都是通用的。公众号可以放代码，所以原稿我是在公众号直接写的，所以排版要好一点，大家可以来公众号看，走路刷刷，每天给我涨个几分钱的广告钱，支持一下。

溜了溜了，88，下周见，88