大家好啊，long time no see。如果把东西切片，每次写的少一点，简单一点，循序渐进，是不是就能改了时不时消失的缺点了。
今天讲一个遗留的问题，三维骨料的建模问题的核心问题—多面体的凸性判定。故事的起源源自于这个：骨料上。当时做了8个面的基体，为啥只做了这8个面，因为当时我只会做这个基体。理论上，基体做完以后，还会涉及点的延拓，再重新生产面，即下图的过程。
在原有的8面体的基础上，又随机一个点，跟基体的点连一连，好像看起来，不是很难的样子，这也是目前生成随机骨料的主流思想，应该是吧。引用一下（本文均指的是凸壳多面体的生成）：
三维凸壳的算法通常有：卷包裹法、分治算法、Z3-3算法和增量算法。
以上的算法，怎么做的，我统统都不知道
，我不知道哪个好，我也不知道怎么实现，但是，大概在写这篇文章之前，我也查阅了一些中文的文献，关于骨料的生产算法，其生成方法有不少，最难实现的也是最核心的部分，我觉得就是怎么判断多面体是不是凸的，所以我今天就写最核心的。
问题来了，什么叫凸性判断？
假设现在空间有不共面的4点，毫无疑问，这个一定可以连成一个四面体：
如果空间中有不共面的5点，怎么连？
现在就有两种连接方法了，左图是凹面体，右图是凸面体。上图是我人连接的，我可以很容易判断凸面体的连接方法，但是脚本并不知道，这这是5点6个面的多面体，如果空间有25点呢，如下图，怎么连？才能保证是凸面体？
上面25个空间点积是我在一个球面上随机生成的，多面体的随机顶点很好生成，不同的随机顶点就对应了随机的多面体骨料，但是，麻烦的在于，怎么判断点与点之间的连接，使得整个多面体保持凸性。我猜，不同的骨料生成算法，其实质就是通过不同的手段保持骨料整体的凸性。
进行凸性判断的算法如下[1]，这是我参考一篇论文的：
设新生成的多面体为凸多面体，令面 CiCjCk (i，j，k=1~7，且 i≠j≠k)为新多面体的任意一个面，Ci、Cj、Ck分别为该面的 3 个顶点，Cp (p=1~7，p≠i，j，k)为除此 3 个顶点之外的任一顶点，当该多面体为凸多面体时，点 Cp 必位于面 CiCjCk 的同一侧。因此，设Vi 为面 CiCjCk 的外法线向量，Vip 为经过点Ci 指向点 Cp 的向量，根据向量运算法则，必有Vi·Vip的值为负。
感兴趣的，可以看一下原文，参考文献我附在文末，我解释一下这段话的意思，它意思说凸性多面体的有一个特别重要的特点，多面体的所有顶点，不会出现在多面体任一个面的两侧，只会在同一侧，我们可以利用这个特点，去判断当前的平面连接是否合理，满足这个条件，说明当前的连接方法合理。
比如下图，之前随机了25个多面体顶点，假定下图的连接方法，可以保证多面体的所有其他顶点均在该面的同一侧，图中的x负方向。这种方式，即所谓的符合条件。
如果下图这种连接方式，顶点分布在面的左右两侧，这种方式连成一定是凹面，甚至有可能连不成封闭空间多面体。这种连接方式，算不符合条件。
程序判断方法，选择一种连接方式，如图中平面abc，定义该面的外法相方向向量为n，遍历所有多面体顶点，定义该点指向与平面abc任意点的连接方向向量为V，如果n·V的值小于0(向量点积，如果值<0，表明两个向量的夹角大于90°；=0，向量垂直；>0，夹角小于90°)，遍历顶点后，如果有值小于0，表明有顶点在面外，该面不符合凸性条件，舍去该面，继续下一个判断。
凸性条件的判断法则即介绍完了，接下来，咱们来实现：
step1 - 生成随机顶点
顶点的随机性决定了生成后的多面体骨料的随机性，这个说简单也很简单，说复杂也有点复杂，看个人需求，如果想最后的骨料长的稍微规矩一点，那得在随机顶点这里下手，不能太随机，我是测试的，所以随机在直径为1的球面上生成，随机两个角度就好了，可以计算出顶点坐标，储存一下，最后封装成一个函数，我们把它叫做createVertices函数吧。
def createVertices(verticeNumber=25):
&#39;&#39;&#39;
:param verticeNumber:number of vertices
:return: coordinates of vertices
&#39;&#39;&#39;
vertices = [] # 储存顶点坐标
for i in range(verticeNumber):
angle1 = np.random.uniform(0, pi*2)
angle2 = np.random.uniform(0, pi*2)
# 随机生成顶点坐标 x,y,z
z = cos(angle1)
x = sin(angle1)*cos(angle2)
y = sin(angle1)*sin(angle2)
vertices.append([x, y, z])
return verticesstep2 - 排列组合定连接方案
这一步的目的，以顶点为总体点集，3个顶点为一组，即一个面，排列出所有的3点的组合，python里有自带排列的工具箱itertools，我百度的，还挺好用的
# 排列组合
surfaces = itertools.combinations(vertices, 3)打印一下结果给各位看看，这个surfaces里是什么样的，其实就是三点的坐标，这里面每一组的三个点就对应了一个平面，接下来，只需要判断，这里面的哪一些平面符合凸性条件，挑出来即可。
[([-0.09482537560790799, 0.9389803435438531, 0.33064189477301104], [0.5710388473688103, -0.3661841964994861, -0.7347270030628562], [-0.12507170945502372, 0.25491098039746707, -0.9588417281109538]),
([-0.09482537560790799, 0.9389803435438531, 0.33064189477301104], [0.5710388473688103, -0.3661841964994861, -0.7347270030628562], [0.7771115167398512, -0.41105885822058597, 0.47657979985378757]),
([-0.09482537560790799, 0.9389803435438531, 0.33064189477301104], [0.5710388473688103, -0.3661841964994861, -0.7347270030628562], [0.11278485993333813, 0.4566570606769418, -0.8824646759523653]),
([-0.09482537560790799, 0.9389803435438531, 0.33064189477301104], [0.5710388473688103, -0.3661841964994861, -0.7347270030628562], [-0.071483670472563, 0.05190044054175497, 0.996090572752971]),
([-0.09482537560790799, 0.9389803435438531, 0.33064189477301104], [0.5710388473688103, -0.3661841964994861, -0.7347270030628562], [-0.736367506151871, -0.5372447547541759, -0.4112553578651586]),
......step3 - 凸性判断
这里就用到了上述的的原理部分了，这里首先要建立一个工具，求解平面的外法向向量，对吧，建立一个名为directN的函数，输入值为三点坐标，返回值为一个外法向向量，大家可以看一下思路，这种求解出的，是外法向的向量。
# 方向向量计算
def directN(coords):
coord1, coord2, coord3 = coords[0],coords[1],coords[2]
x1, y1, z1 = coord1[0], coord1[1], coord1[2]
x2, y2, z2 = coord2[0], coord2[1], coord2[2]
x3, y3, z3 = coord3[0], coord3[1], coord3[2]
a = (y2-y1)*(z3-z1) - (y3-y1)*(z2-z1)
b = (z2-z1)*(x3-x1) - (z3-z1)*(x2-x1)
c = (x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)
n = np.array([a, b, c])
n = n/np.linalg.norm(n)
d = arccos(np.dot(n, coord1)/np.linalg.norm(coord1))
if d&gt;pi/2.:
n = -n
return n外法向方向向量的工具求解完成后，接下来就是找出凸性面了，整体思路就是遍历所有顶点，与平面的任一个顶点连接的向量V，计算该向量与平面法向之间的角度，向量点积就行了，numpy里的dot函数，根据值的正负取舍，这里是夹角要小于90°，点积后的值应该大于0。所以一旦有小于0的值出现，表明该平面一定不是凸性面，就不用接下来的判断了，直接break就行。
我们以随机的顶点vertices为输入值，封装成如下的函数，这个函数大家可以直接拿去用，把随机后的顶点坐标丢给它，它会帮你判断出哪些是可以用于建模的凸性面，然后提供给你选择后的返回值，这就是今天的重点部分。
def chosePlane(vertices):
vertices = [np.array(vertice) for vertice in vertices]
threePoints = list(itertools.combinations(vertices, 3))
chosenPlane = []
for threePoint in threePoints:
n = directN(threePoint)
sign = True
for vertice in vertices:
vector1 = threePoint[0]-vertice
result = np.dot(vector1, n)
if result&lt;-1e-5:
sign = False
break
if sign:
chosenPlane.append(list(threePoint))
return chosenPlane
step4 - 工具使用及建模
最后一步，建完了工具，当然是测试使用情况了。
首先使用上面三个函数，把凸性面的数据选择处理，这部分很简单。把凸性面的数据储存在chosePlane变量里。
vertices = createVertices(verticeNumber=25)
chosePlane = chosePlane(vertices)获得了数据后，放在abaqus里建模就好了，这一部分也很简单，直接丢源码就行，各位自己尝试
# create model
#
if mdb.models.has_key(&quot;model-test&quot;):
del mdb.models[&quot;model-test&quot;]
model = mdb.Model(name=&quot;model-test&quot;, modelType=STANDARD_EXPLICIT)
part = model.Part(name=&quot;part-test&quot;, dimensionality=THREE_D, type=DEFORMABLE_BODY)
# create point
for vertice in vertices:
part.DatumPointByCoordinate(coords=tuple(vertice))
# create plane
for coords in chosePlane:
coords.append(coords[0])
wire = part.WirePolyLine(mergeType=SEPARATE, meshable=ON, points=(coords))
face_edge = part.getFeatureEdges(name=wire.name)
part.CoverEdges(edgeList = face_edge, tryAnalytical=True)
结果展示
必不可少的结果展示部分
25顶点的结果如下，大家可以看到，基本上保持了它整体的凸性，但是因为点的随机性比较大，所以骨料不会显得那么圆润，对于这一点，可以从随机点的生成下手，把点随机的的分布弄的更均匀，结果自然会好一点。
对于我这种随机点的生成方法，由于我是在圆上随机产生的，所以点越多，结果肯定越圆，我试一下，50，100，150和200点的结果如何(点越多，越费时间，我是算着玩的，正常顶点应该不会超过50个吧，太多了有限元网格画的费劲)：
明显的越来越圆，这也与我们预测的一致，当然其实很明显结果
源码自取，我还是会上传到github里的，我都忘了怎么提交了，又得复习一下命令
# encoding=utf-8
# code by wuhuiguo
#
from abaqus import *
from abaqusConstants import *
import numpy as np
import itertools
from numpy import pi,sin,cos,arccos
def createVertices(verticeNumber=25):
&#39;&#39;&#39;
:param verticeNumber:number of vertices
:return: coordinates of vertices
&#39;&#39;&#39;
vertices = [] # 储存顶点坐标
for i in range(verticeNumber):
angle1 = np.random.uniform(0, pi*2)
angle2 = np.random.uniform(0, pi*2)
# 随机生成顶点坐标 x,y,z
z = cos(angle1)
x = sin(angle1)*cos(angle2)
y = sin(angle1)*sin(angle2)
vertices.append([x, y, z])
return vertices
# 方向向量计算
def directN(coords):
coord1, coord2, coord3 = coords[0],coords[1],coords[2]
x1, y1, z1 = coord1[0], coord1[1], coord1[2]
x2, y2, z2 = coord2[0], coord2[1], coord2[2]
x3, y3, z3 = coord3[0], coord3[1], coord3[2]
a = (y2-y1)*(z3-z1) - (y3-y1)*(z2-z1)
b = (z2-z1)*(x3-x1) - (z3-z1)*(x2-x1)
c = (x2-x1)*(y3-y1) - (x3-x1)*(y2-y1)
n = np.array([a, b, c])
n = n/np.linalg.norm(n)
d = arccos(np.dot(n, coord1)/np.linalg.norm(coord1))
if d&gt;pi/2.:
n = -n
return n
# 平面判断
def chosenPlane(vertices):
vertices = [np.array(vertice) for vertice in vertices]
threePoints = list(itertools.combinations(vertices, 3))
chosenPlane = []
for threePoint in threePoints:
n = directN(threePoint)
sign = True
for vertice in vertices:
vector1 = threePoint[0]-vertice
result = np.dot(vector1, n)
if result&lt;-1e-5:
sign = False
break
if sign:
chosenPlane.append(list(threePoint))
return chosenPlane
#
if mdb.models.has_key(&quot;model-test&quot;):
del mdb.models[&quot;model-test&quot;]
model = mdb.Model(name=&quot;model-test&quot;, modelType=STANDARD_EXPLICIT)
part = model.Part(name=&quot;part-test&quot;, dimensionality=THREE_D, type=DEFORMABLE_BODY)
vertices = createVertices(verticeNumber=20)
chosePlane = chosenPlane(vertices)
# create point
for vertice in vertices:
part.DatumPointByCoordinate(coords=tuple(vertice))
# create plane
for coords in chosePlane:
coords.append(coords[0])
wire = part.WirePolyLine(mergeType=SEPARATE, meshable=ON, points=(coords))
face_edge = part.getFeatureEdges(name=wire.name)
part.CoverEdges(edgeList = face_edge, tryAnalytical=True)
核心思想我已经摆出来了，其实三维骨料的建模下集篇的所有内容，基本上我今天都写完了，算是还愿了。但是考虑到，大家肯定还有不知道怎么下手的，肯定还有，所以我会把完整三维骨料建模小脚本下次发出来。估计那才是你们需要的对不对。
再来个预告，这个完了以后，我想做一个关于voronoi多边形界面建模的教学，这个其实我以前做过一次，但是是二维的，就讲个二维的吧，有小伙伴提到过这个界面的建模问题。先把二维的搞出来，然后再试试三维的。88，溜了溜了
参考文献
【1】方秦,张锦华,还毅,张亚栋.全级配混凝土三维细观模型的建模方法研究[J].工程力学.2013:(30)14-21.
【2】www.baidu.com，我参考了一大堆百度的资料，这次我没记录给忘记了，其实里面有很多东西值得反复观看观看的，下次写的时候我留个心眼，都记录一下
