今天我們來學習怎麼樣能夠針對1st order solid elements，建立segment contact。

Why Segment Contact

可能有人會問為什麼要特別用segment contact呢？Part set to Part set不就可以一招打天下了嗎？

原來是有些情況，你已經事先知道contact的區域，可能僅佔該part的一小部份。如果為了方便設定，使用part set來設定contact，會包含過多無用的elements，會降低contact的計算效率。但是偏偏這類型的題目所contact的區域，大多是幾何外型不容易選取的地方，相比於選取的難度，很多人會傾向寧願選擇多花費計算資源。

不容易操作的繁複工作，這正正是咱們二次開發發揮的好機會呀！

我們的目標是做出讓使用者點選兩組Property(每組內可含1~多個Property)，就能建立在某一tolerance下的segment contact。

我們拍了一小段影片，大家可以先睹為快(此範例tolerance=50)，請點選下圖。

Sounds interesting, right? Let’s move on!

解題思路

這個問題的關鍵在於，如何快速找出兩邊符合tolerance下的solidfacet。

經過一番努力(google...)之後，我們發現scipy.spatial.KDTree有一個query_ball_point相當符合我們的需求。這真是太好了，因為scipy正是ANSA附的五個third-party package之一。

關於KDTree的詳細說明，可能需要參考scipy的文件說明或machine learning相關的資料，例如scikit-learn。如果只是想簡單了解，我們發現這篇知乎的文章是個不錯的參考。

我們總共有以下幾個function需要講解：

• create_shells
• _get_s
• _query_ball_set
• _get_candidate
• remove_skin_related
• main

至於create_set、create_segment及create_contact及schemas部份，這些都是老朋友了，不佔篇幅了。

create_shells

• 透過base.CollectEntities收集prop中的所有solid element。
• 透過base.CreateShellsFromSolidFacets搭配skin exclude internal bounds建立skin shell並回傳。

#quick_1st_solid_ct.py
def create_shells(prop):
    solids = base.CollectEntities(deck, prop, LSDYNAType.SOLID)
    return base.CreateShellsFromSolidFacets(
        "skin exclude internal bounds",
        ret_ents=True,
        solids=solids)

_get_s

• _get_s除收集各shell的node.position，還收集了各shell node.position的中點，可以使得搜尋結果比較smooth，不會東一塊西一塊。

#quick_1st_solid_ct.py
def _get_s(shells):
    fields = ('N1', 'N2', 'N3', 'N4')  # can cover both tri and quad
    s = []
    for shell in shells:
        result_dict = shell.get_entity_values(deck, fields)
        pts = [node.position for node in result_dict.values()]
        for pt in pts:
            s.append((shell, pt))
        s.append((shell, np.mean(pts, axis=0)))
    return s

_query_ball_set

• spatial.KDTree.query_ball_point可以找出ptas內所有與ptbs在seartol範圍內的點。比較特別的是其回傳的是index，所以我們需要用一個set comprehensions取出所有index，再用另一個set comprehensions由shs取出相對應index的shell後並回傳。

#quick_1st_solid_ct.py
def _query_ball_set(ptas, ptbs, shs, seartol):
    # 1->2 : pt1s, pt2s, sh1s
    # 2->1 : pt2s, pt1s, sh2s
    tree = spatial.KDTree(ptas)
    balls = tree.query_ball_point(ptbs, seartol)
    idx_set = {idx
               for ball in balls
               for idx in ball}
    return {shs[idx] for idx in idx_set}

_get_candidate

• _get_candidate算是一個utility function，幫我們分別取得兩邊的shell set。

#quick_1st_solid_ct.py
def _get_candidate(shells1, shells2, seartol):
    sh1s, pt1s = zip(*_get_s(shells1))
    sh2s, pt2s = zip(*_get_s(shells2))

    set1 = _query_ball_set(pt1s, pt2s, sh1s, seartol)
    set2 = _query_ball_set(pt2s, pt1s, sh2s, seartol)

    return set1, set2

remove_skin_related

• 於create_shells建立shell時，會產生ansapart及property Entity，其名都會是以Skin of開頭，所以我們可以透過base.NameToEnts一次收集以Skin of開頭的Entity，然後透過一個list comprehensions過濾出ansapart及property Entity。最後透過base.DeleteEntity刪除這些多餘的Entity。

#quick_1st_solid_ct.py
def remove_skin_related():
    default_entities = base.NameToEnts("^Skin\sof\s.*")
    if default_entities:
        to_dels = [ent
                   for ent in default_entities
                   if ent.ansa_type(deck) == LSDYNAType.PART
                   or ent.ansa_type(deck) == SecType.SECTION_SHELL]
        base.DeleteEntity(to_dels, force=False)

main

• 分別透過create_shells對兩組property建立shells1及shells2。
• 透過_get_candidate取得ball_set1及ball_set2兩個shell set。
• 透過create_segment分別對ball_set1及ball_set2內的shell建立segment Entity。
• 透過create_set分別建立seg1及seg1。
• 透過create_contact建立contact Entity。
• 呼叫remove_skin_related刪除多餘Entity。

#quick_1st_solid_ct.py
def main(prop1, prop2, seartol=1.0):
    print(f'contact creation is initialized')
    shells1 = create_shells(prop1)
    shells2 = create_shells(prop2)

    ball_set1, ball_set2 = _get_candidate(shells1, shells2, seartol)

    seg1 = create_segment(ball_set1)
    seg_set1 = create_set(seg1)

    seg2 = create_segment(ball_set2)
    seg_set2 = create_set(seg2)

    contact = create_contact(ssid=seg_set2._id,
                             msid=seg_set1._id,
                             sstyp=ContactType.TYPE0_SEGMENT_SET.value,
                             mstyp=ContactType.TYPE0_SEGMENT_SET.value)

    # remember to delete
    remove_skin_related()
    print(f'contact creation is done')

if name == 'main'

• 透過base.PickEntities讓使用者可以選擇兩組property。
• 檢查確實選擇了兩組property後，呼叫main function。

#quick_1st_solid_ct.py
if __name__ == '__main__':
    deck = constants.LSDYNA
    type_ = LSDYNAType.PROPERTY
    seartol = 50
    prop1 = base.PickEntities(deck, [type_], initial_type=type_)
    if prop1:
        prop2 = base.PickEntities(deck, [type_], initial_type=type_)
        if prop2:
            main(prop1, prop2, seartol)

小結

實作完成這個小script後，發現真的要寫的code並不多，這都給歸功於scipy有提供好用的function呀！

說實話，我們一開始實作的時候，是自己土炮的去打迴圈計算距離，效率實在不佳，還好google及stackoverflow又再次拉了我們一把。

可能有先進會問2nd solid elements可不可以也這樣處理呢？

其實是可以的，但是還是那句，留一點給我們賺嘛XD

Code

本日程式碼傳送門。