基於Voronoi二維多晶體的衍生腳本操作2

小夥伴們晚上好呀。2021年1月已經過完了，今天是2月的第一天，時間如風，略過我英俊的臉龐，無聲無息。過的是真快。

今天呢，還是談一談腳本。雖然我已經很久沒寫子程序了。我對腳本的興趣更大。小夥伴們期待的umeshmotion，我記著的，會寫的。

今天還是談那個二維的泰森多邊形，最近看到有粉絲的留言，大家的留言我基本都能看到，一般都會回復。有位小夥伴問我怎麼讓二維的晶粒生成的更均勻。我回覆說 「對隨機的種子點加以限制」 。因為我寫的例子，種子點都是完全隨機的，所以晶粒也是完全隨機的，所以我想的是，把種子點弄的稍微均勻一點，不就行了嗎。他回復我怎麼限制，應該怎麼做。

到此我就感覺到了大家的一些問題所在，因為有這些問題的小夥伴，我相信是寫腳本沒多久的，甚至是寫代碼沒多久的。所以各位沒有什麼大局觀，大家記住寫腳本，具有很強的目的性，我們就是通過各種手段，去實現我們的目的。不做任何無用功，不寫任何沒用的語句（注釋不算）。對於這些小夥伴，我覺得可能沒有看懂我寫的代碼，我寫的所有腳本，都不是一個整體，是一個個小塊，是可以通過各種排列組合和改寫，去實現你自己的問題的。大家下載之後，可以去嘗試，去理解，去嘗試我在幹什麼。

如果你剛學沒多久，可能對我現在寫的東西沒什麼感覺。但，如果你是一個想學，而苦於沒有好資源的時候，來看我寫的內容，那肯定是大有裨益的。所以，如果你現在不打算學，或者還不知道它的價值，沒關係，收藏一下。相信總有一天，你會用得到。

今天的內容，就是為各位有這些問題的小夥伴準備的，更偏向於新手，我想教教大家一點寫腳本的思維。

今日任務

上面的問題，想實現均勻的晶粒，核心就是讓隨機的形核點均勻。那我們現在的任務就變成了怎麼讓形核點均勻。來看一下，我之前是怎麼隨機的。

python腳本實現abaqus前處理2D多晶粒建模(附完整源碼)-Voronoi多邊形的生成

points = np.array([[random.uniform(0,length*1.2), random.uniform(0,width*1.2)] for i in range(point_number)])
這是我之前的腳本，對於形核點的處理。是個完全隨機的狀態，所以最後的結果就像下圖這種：
因為是完全隨機的，所以每個晶粒並不均勻，這是正常的。現在時代變了，咱們有新的需求了。
那咱們就對形核點進行處理，接下來就是咱們分析大佬的工作了。怎麼個均勻法，這個看個人需求。我這裡做三種，給大家展示一下做法。大家要學思路，我寫的這些話是沒什麼用的，重要是編程思維。謹記。
均勻形核點—方法一

我這第一個均勻的方法是，矩形整列固定核心點。什麼意思呢，意思就是，形核點，差不過是下面這種。
因為這些點點是我草稿紙的底圖，所以看起來不太清楚。但是，相信已經能描述我的意思了。如果，形核點是這種，那麼晶粒應該是什麼樣呢。咱們來寫一下，這些形核點改怎麼生成。
基體的長寬已知，形核點數量已知。假設橫向和縱向的形核點距離一樣。如上解法，可以把形核點的間距解出來。接下來，咱們就可以寫代碼了。
# 方法一points = []x = np.sqrt(length*width/point_number)for i in range(int(length/x)):    for j in range(int(width/x)):        points.append([i*x, j*x])
展示一下結果，夠不夠均勻，兄弟們。覺得均勻的兄弟們，素質三聯走一波，給up漲漲人氣。這個畫圖的是python3寫的。源碼我在文章末尾會分享給你們。咱們先把目標集中在，怎麼生成均勻形核點上。
又有兄弟，可能有疑問了，這個怎麼在abaqus裡生成呢。很簡單，只需要把以前的形核點替換掉即可。展示一下結果，大家自己動手試試，思路和代碼都寫了，移動一下，不是問題吧。
我做一下，證明是可以的。在abaqus裡運行結果如下。大問題沒有，但是有個小問題，看圖左下角。這個我也不知道發生了什麼，如果大家發現有少量缺胳膊少腿的情況，可以在cae裡補畫一下。
均勻形核點—方法二
方法二是什麼呢，其實是我想生成一些正六邊形的晶粒。那麼形核點的位置如何，該怎麼生成。接下來，咱們展示一下，思路。注意啊，思維是最重要的，我寫的代碼，是不重要的，重要的是算法，代碼之後的公式，思維。
把橫向和縱向的間距都求解完之後。就是寫碼的時間了。直接上，這裡注意一點就是，這個不是陣列了，每一層要錯個位，距離是x/2。
points = []x = np.sqrt(2*np.sqrt(3)*length*width/point_number)print(x,int(length/(x/2))+2)for i in range(int(length/(x/2))+2):    for j in range(int(width/(x/np.sqrt(3)))+2):        if i%2 == 0:            points.append([i*x/2, j*x/np.sqrt(3)])        else:            points.append([i*x/2, (j+0.5)*x/np.sqrt(3)])
怎麼樣，夠均勻了吧。還是老規矩，把它轉移到，我們之前的腳本裡去，測試一下成不成功。
這次沒有什麼小瑕疵了，成功運行。有沒有點意思了，只改動了很少的部分。就能有不同的結果。大家再拿到我寫的腳本之後，可以發散性的思考去嘗試，別管它有沒有用，別總是局限在你的項目，課題裡了。幹點自己感興趣的，它不香嗎？
均勻形核點—方法三
由於上面兩個，展示的都是固定形核的辦法，那麼接下來，我要展示一下，隨機形核的辦法。這個要怎麼做呢，很簡單，兄弟們。又想隨機，又想均勻。只需要如何？
只需要讓隨機的區域均勻即可，是不是這個理。啥意思，畫圖解釋一下：
兄弟們，你們懂我意思吧。什麼叫做，隨機的區域均勻。懂了不。我感覺我畫的還算清楚，你們都懂我意思的吧。有不懂的，我單獨指點一下，評論@我一下。
下面咱們就說一下，這一步，應該怎麼做。
確定算法之後，剩下的就是代碼實現了，這一步比較簡單，咱們直接上。
points = []x = np.sqrt(length*width/point_number)for i in range(int(length/x)+2):    for j in range(int(width/x)+2):        px = uniform(i*x, (i+1)*x)        py = uniform(j*x, (j+1)*x)        points.append([px, py])
這樣是不是看起來，就均勻許多了。
還是老規矩，移植到abaqus腳本裡去。看一下效果。
看起來，確實比之前的要均勻許多，晶粒大小上要均勻許多。
最後總結

以上就是我以三個例子去為大家闡述，怎麼個寫法。各位以後在寫腳本的時候，注意下。在你寫腳本之前，你就應該知道怎麼寫了。邏輯，編程的思維，才是最最要緊的。
今天的內容就到這裡，咱們下期再見。兄弟們再見。快過年了，祝願大家天天開心，水平蹭蹭的漲。
源碼

1 python3畫圖的代碼
from scipy.spatial import *from random import *import matplotlib.pyplot as pltimport numpy as np
point_number = 200length = 200width = 100


points = []x = np.sqrt(length*width/point_number)for i in range(int(length/x)+2):    for j in range(int(width/x)+2):        px = uniform(i*x, (i+1)*x)        py = uniform(j*x, (j+1)*x)        points.append([px, py])
vor = Voronoi(points)voronoi_plot_2d(vor)plt.show()
2 python腳本，或者去github取，地址你們懂得，後臺回復「腳本」
from abaqus import *from abaqusConstants import *import regionToolsetfrom scipy.spatial import *import numpy as npimport random
point_number = 250length = 150width = 100



points = []x = np.sqrt(length*width/point_number)for i in range(int(length/x)+2):    for j in range(int(width/x)+2):        px = random.uniform(i*x, (i+1)*x)        py = random.uniform(j*x, (j+1)*x)        points.append([px, py])points = np.array(points)

vor = Voronoi(points) #create instancevertices = vor.verticesedges = vor.ridge_vertices
myModel = mdb.models["Model-1"]mySketch1 = myModel.ConstrainedSketch(name='sketch1', sheetSize=200.0)mySketch1.rectangle(point1=(0.0, 0.0), point2=(length, width))myPart = myModel.Part(name='Part-voronoi', dimensionality=TWO_D_PLANAR,type=DEFORMABLE_BODY)myPart.BaseShell(sketch=mySketch1)
mySketch2 = myModel.ConstrainedSketch(name='__profile__',sheetSize=200, gridSpacing=10)mySketch2.CircleByCenterPerimeter(center=(-61.49, 2.795), point1=(-44.72, -8.385))
for edge in np.array(edges):    if np.all(edge>0):        mySketch2.Line(point1=tuple(vertices[edge[0]]), point2=tuple(vertices[edge[1]]))center = points.mean(axis=0)for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices):    simplex = np.asarray(simplex)    if np.any(simplex < 0):        i = simplex[simplex >= 0]        t = points[pointidx[1]] - points[pointidx[0]]        t = t/np.linalg.norm(t)        n = np.array([-t[1], t[0]])        midpoint = points[pointidx].mean(axis=0)        far_point = vertices[i] + np.sign(np.dot(midpoint - center, n))*n*100        mySketch2.Line(point1=tuple(vertices[i[0]]), point2=tuple(far_point[0]))
myPart.PartitionFaceBySketch(faces=myPart.faces[:], sketch=mySketch2)