今天我們承接[Day13]，繼續改進box。

我們將會學習如何將box由QUAD4的shell改為HEXA8的solid，概念與[Day13]差不多，只是複雜度較高。

核心概念

[Day13]我們透過了_gen_one_layer_grids產生了一層grids。

• 同理可推，我們也可以透過呼叫兩次_gen_one_layer_grids，給予兩個不同的z_elv，產生兩層高度不同的grids，那麼這兩層grids就可以形成一層的HEXA element。
• 依此類推，如果呼叫三次_gen_one_layer_grids，可以形成兩層HEXA element
• 依此類推，如果...

透過調整不同的z_elv，並多次呼叫_gen_one_layer_grids，就可以建造任意層數的HEXA element。

Sounds good？但心裡的coding魂總是覺得哪裡不太對勁呀？

經過一番查找，發現原來是多次呼叫這個詞，觸發了警報。

觀察_gen_one_layer_grids後發現，np.meshgrid(x, y)會被多次計算，但每次的計算值都是一樣的。

針對這個情況，Python內建的functools.lru_cache function可以幫助我們。lru_cache是一個decorator，當一個function被呼叫超過一次，而所給的參數都一樣時，會直接回傳快取值，而不實際進行運算。

我們可以將np.meshgrid(x, y)前的code獨立為_get_XY function，並利用@lru_cache
進行裝飾。

# gen_seqs.py
from functools import lru_cache
from scipy.spatial.transform import Rotation as scipy_rotation


@lru_cache
def _get_XY(l, w, en1, en2):
    m1, m2 = en1+1, en2+1
    x = np.linspace(0, l, m1)
    y = np.linspace(0, w, m2)
    return np.meshgrid(x, y)


def _gen_one_layer_grids(l, w, en1, en2, *, z_elv=None, move_xy=None, rot_angle=None):
    X, Y = _get_XY(l, w, en1, en2)

    x_shape = X.shape
    z_elv = z_elv or 0
    Z = np.ones(x_shape)*z_elv

    if rot_angle is not None:
        XYZ = np.array([X.ravel(), Y.ravel(), Z.ravel()]).transpose()
        r = scipy_rotation.from_rotvec(
            rot_angle*np.array([0, 0, 1]), degrees=True)
        XYZrot = r.apply(XYZ)
        X = XYZrot[:, 0].reshape(x_shape)
        Y = XYZrot[:, 1].reshape(x_shape)
        Z = XYZrot[:, 2].reshape(x_shape)

    if move_xy is not None:
        ox, oy = move_xy
        X += ox
        Y += oy

    yield from (X.flat, Y.flat, Z.flat)

Box Grids

於gen_solid_grids我們透過一個list comprehensions收集每一層的grids到layer_grids中。接著透過zip及itertools.chain function來回傳每一個grid的座標值。

# gen_seqs.py
def gen_solid_grids(l, w, h, en1, en2, en3, *, z_elv, move_xy, rot_angle):
    layer_grids = [_gen_one_layer_grids(l, w, en1, en2,
                                        z_elv=z_elv+i*h/en3,
                                        move_xy=move_xy,
                                        rot_angle=rot_angle)
                   for i in range(en3+1)]

    x, y, z = zip(*layer_grids)
    return zip(chain.from_iterable(x),
               chain.from_iterable(y),
               chain.from_iterable(z))

Box Seqs

如果說[Day13]的Plate Seqs需要一點想像力，那麼Box Seqs看起來就需要更多一點了...

我們知道LS-DYNA的HEXA8 element是逆時針建立其四個node後，再往上一層繼續逆時針建立四個node，所以這次我們想要達成的是:

• 能夠沿著x方向，不斷逆時針寫出每個HEXA8 element的八個node id，直到x方向寫到了指定個數
• 想像自己的y座標往上走一層，繼續往x方向寫指定個數，直到y方向也寫到指定個數。
• 想像自己的z座標往上走一層，重複前面兩個步驟，直到z方向也寫到指定個數。

# gen_seqs.py
def gen_solid_seqs(en1, en2, en3, *, start=1):
    m1, m2 = en1+1, en2+1
    m12 = m1*m2
    en = en1*en2*en3

    cnt_x, cnt_y, total_cnt = 0, 0, 0
    while True:
        if total_cnt == en:
            break
        if cnt_x == en1:
            cnt_x = 0
            cnt_y += 1
            if cnt_y == en2:
                start += (en1+2)
                cnt_y = 0
            else:
                start += 1

        n1, n2, n3, n4 = start, start+1, start+en1+2, start+en1+1
        n5, n6, n7, n8 = n1+m12, n2+m12, n3+m12, n4+m12
        yield (n1, n2, n3, n4, n5, n6, n7, n8)

        total_cnt += 1
        cnt_x += 1
        start += 1

舉例來說，這邊我們想要建立一個x、y及z各有2個element的box。

en1, en2, en3 = 2, 2, 1
print(f'{list(gen_solid_seqs(en1, en2, en3))=}')

其輸出為:

list(gen_solid_seqs(en1, en2, en3))=[(1, 2, 5, 4, 10, 11, 14, 13), (2, 3, 6, 5, 11, 12, 15, 14), (4, 5, 8, 7, 13, 14, 17, 16), (5, 6, 9, 8, 14, 15, 18, 17)]

由於gen_solid_seqs是一個generator，所以需要使用list才能print出其內的元素。

Box Nodes

有了gen_solid_grids與gen_solid_seqs之後，我們可以開始來實作create_box_nodes function。

• create_box_nodes先透過gen_solid_grids取得solid_grids的generator。
• 透過一個list comprehensions配合create_node及上述generator，產生了node Entities，並收集為一list。
• 呼叫gen_solid_seqs，得到一個generator，作為下一步create_box_h8solids的準備。
• 回傳所產生的node Entities及solid_seqs。

# creators.py
def create_box_nodes(l,
                     w,
                     h,
                     en1,
                     en2,
                     en3,
                     node_start_id,
                     z_elv=None,
                     move_xy=None,
                     rot_angle=None,
                     deck=None):
    deck = deck or constants.LSDYNA
    solid_grids = gen_solid_grids(
        l, w, h, en1, en2, en3, z_elv=z_elv, move_xy=move_xy, rot_angle=rot_angle)
    keys = ('NID', 'X', 'Y', 'Z')
    nodes = [create_node(dict(zip(keys, (nid, *xyz))), deck=deck)
             for nid, xyz in enumerate(solid_grids, start=node_start_id)]
    solid_seqs = gen_solid_seqs(en1, en2, en3, start=node_start_id)
    return nodes, solid_seqs

Box Solids

create_box_h8solids的寫法與create_box_nodes差不多。我們透過create_box_nodes的第二個回傳值solid_seqs配合create_solid來產生HEAX8 Entities，並收集為一list後回傳。

# creators.py
def create_box_h8solids(solid_seqs, solid_start_id, pid, deck=None):
    deck = deck or constants.LSDYNA
    keys = ('type', 'PID', 'EID', 'N1', 'N2',
            'N3', 'N4', 'N5', 'N6', 'N7', 'N8')
    type_ = SolidType.HEXA
    return [create_solid(dict(zip(keys, (type_, pid, eid, *ns))), deck=deck)
            for eid, ns in enumerate(solid_seqs, start=solid_start_id)]

小結

至此，我們已經成功將box由shell升級為solid了，下圖是我們用LS-PrePost錄的gif動畫。

Code

本日程式碼傳送門。