SalomeからPyramid要素を含むメッシュをエクスポートして使用する方法

# 要旨
単純な形状であれば、計算精度の良い六面体メッシュ＋プリズムでメッシュを切ることができます。
しかし、規則性の無い形状が含まれている場合などに四面体との混合メッシュになることがあります。
Salomeによりメッシュを切ると、六面体と四面体との接続部分にピラミッド要素が作成されます（PrePoMaxではPyramidを除外する設定もできるようです）。
境界のメッシュを細かくするなど電気系での解析に必要な工夫はしますが、メッシュの種類変更などに時間を掛けたくありません。
キーワードに挙げたCAEソフトにおいて、Pyramid要素の取り込み・計算に成功しましたので、報告します。
# CAEソフトの対応状況
### 1. Salome-meca
Salome-mecaのSalomeにより作成したメッシュを用いて、メッシュの不具合により計算を失敗したことがありません。
内側から三角、四角、三角の平面メッシュを作り、Y軸周りに回転させたモデルです。
四面体、プリズム、ピラミッド、六面体の複合になりました。

![Salome-meca-Mesh.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4014034/91fc876d-779a-4bbd-bc30-df38d50a2d01.jpeg "Salomeにより作成したメッシュとInfomation")
### 2. Gmsh + GetDP (= Onelab)
SalomeのMED形式に対応していますので、Med形式によりエクスポートすれば問題なく使えます。

![Gmsh-Mesh-2枚図.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4014034/5b1ecc87-db19-4b53-b2c6-0f65bcfe22ec.png "MEDのライブラリが組み込まれています")
### 3. Elmer
二つの方法があります。変換回数が少なくて済むのは、salomeToElmerのPythonスクリプト（Warningの表示防止のため１行修正しました）です。
#### (1) salomeToElmerのPythonスクリプトを使用する方法
- SalomeをMeshメニューにする
- 変換するメッシュを選択する
- File / load Script... から、salomeToElmer.pyを選ぶ
- SalomeのプログラムのフォルダにElmerSolver用のメッシュが出力される
- ElmerGrid 2 2 Mesh_1 -autoclean (-increase) # 不要な重複の除去。increaseは二次要素に変換したい場合 

```python:salomeToElmer-r1.py
u"""
Script for exporting mesh from Salome to Elmer
"""
# ******************************************************************************
#
#    Copyright (C) 2017-, Juris Vencels
#
#    Authors: Juris Vencels
#    Email:   juris.vencels@gmail.com
#    Web:     http://vencels.com
#    Address: University of Latvia
#             Laboratory for mathematical modelling of 
#                 environmental and technological processes
#             Zellu Str. 23, Riga, LV-1002, Latvia
#
#    Original Date: 05.05.2017
#
#    This script is based on salomeToOpenFOAM.py written by Nicolas Edh
#    https://github.com/nicolasedh/salomeToOpenFOAM
#
# *****************************************************************************
#
#    This script is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    salomeToElmer is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program (in the file LICENSE); if not, write to the 
#    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 
#    Boston, MA 02110-1301, USA.
#

import sys
import SMESH
from salome.smesh import smeshBuilder
import os,time

def exportToElmer(mesh,dirname='salomeToElmer'):
    u"""
    Elmer's native mesh consists of 5 files.

    mesh.header - first line (# nodes, # all elements, # edge and boundary elements)
                  second line (# element types)
                  third and all following lines (type, # elements of this type)

    mesh.names - body and boundary names with their IDs taken from Salome mesh groups.
                 The last 'empty' belongs to edge and boundary elements that do not
                 belong to user specified groups. 

    mesh.nodes - every line represents one node (node ID, dummy, X, Y, Z)

    mesh.elements - every line represents one volume element (volume element ID,
                    body ID, type, node1, node2, node3, ...)

    mesh.boundary - every line represents one edge or boundary element (edge or 
                    boundary element ID, boundary ID, parent volume element 1, 
                    parent volume element 2, type, node1, node2, ...)
    """
    tstart=time.time()

    if not os.path.exists(dirname):
        os.makedirs(dirname)
    try:
        fileHeader   = open(dirname + "/mesh.header",'w')
        fileNodes    = open(dirname + "/mesh.nodes",'w')
        fileNames    = open(dirname + "/mesh.names",'w')
        fileElements = open(dirname + "/mesh.elements",'w')
        fileBoundary = open(dirname + "/mesh.boundary",'w')
    except Exception:
        print("ERROR: Cannot open files for writting")
        return

    meshIs3D = False
    if mesh.NbVolumes()>0:
        meshIs3D = True

    # mesh.header
    if meshIs3D:
        print("Exporting 3D mesh..\n")
        fileHeader.write("%d %d %d\n" \
            %(mesh.NbNodes(),mesh.NbVolumes(),mesh.NbEdges()+mesh.NbFaces()))
    else:
        print("Exporting 2D mesh..\n")
        fileHeader.write("%d %d %d\n" \
            %(mesh.NbNodes(),mesh.NbFaces(),mesh.NbEdges()))

    elems = {str(k): v for k, v in mesh.GetMeshInfo().items() if v}
    fileHeader.write("%d\n" %(len(elems.values())-1))

    elemTypeNames = {'202': 'Entity_Edge', '303': 'Entity_Triangle', \
                     '404': 'Entity_Quadrangle', '504': 'Entity_Tetra', \
                     '605': 'Entity_Pyramid', '706': 'Entity_Penta', \
                     '808': 'Entity_Hexa'}

    for nbr, ele in sorted(elemTypeNames.items()):
        if elems.get(ele):
            fileHeader.write("%s %d\n" %(nbr,elems.get(ele)))

    fileHeader.flush()
    fileHeader.close()

    # mesh.nodes
    points=mesh.GetElementsByType(SMESH.NODE)
    for ni in points:
        pos=mesh.GetNodeXYZ(ni)
        fileNodes.write("%d -1 %.12g %.12g %.12g\n" %(ni,pos[0],pos[1],pos[2]))
    fileNodes.flush()
    fileNodes.close()

    # initialize arrays
    invElemType = {v: k for k, v in elemTypeNames.items()}
    invElemIDs = [mesh.NbGroups()+1 for el in range(mesh.NbElements())]
    elemGrp = list(invElemIDs)

    edgeIDs = mesh.GetElementsByType(SMESH.EDGE)
    faceIDs = mesh.GetElementsByType(SMESH.FACE)
    elemIDs = edgeIDs + faceIDs
    NbBoundaryElems = mesh.NbEdges()

    if meshIs3D:
        volumeIDs = mesh.GetElementsByType(SMESH.VOLUME)
        elemIDs = elemIDs + volumeIDs
        NbBoundaryElems = NbBoundaryElems + mesh.NbFaces()

    if len(elemGrp) != max(elemIDs):
        raise Exception("ERROR: the number of elements does not match!")

    for el in range(mesh.NbElements()):
        invElemIDs[elemIDs[el]-1] = el+1

    # mesh.names
    fileNames.write("! ----- names for bodies -----\n")
    groupID = 1

    if meshIs3D:
        bodyType = SMESH.VOLUME
        boundaryType = SMESH.FACE
    else:
        bodyType = SMESH.FACE
        boundaryType = SMESH.EDGE

    for grp in mesh.GetGroups(bodyType):
        fileNames.write("$ %s = %d\n" %(grp.GetName(), groupID))
        for el in grp.GetIDs():
            elemGrp[invElemIDs[el-1]-1] = groupID
        groupID = groupID + 1

    fileNames.write("! ----- names for boundaries -----\n")

    for grp in mesh.GetGroups(boundaryType):
        fileNames.write("$ %s = %d\n" %(grp.GetName(), groupID))
        for el in grp.GetIDs():
            if elemGrp[invElemIDs[el-1]-1] > groupID:
                elemGrp[invElemIDs[el-1]-1] = groupID
        groupID = groupID + 1

    fileNames.write("$ empty = %d\n" %(mesh.NbGroups()+1))

    fileNames.flush()
    fileNames.close()

    # mesh.elements
    for el in mesh.GetElementsByType(bodyType):
        elemType = mesh.GetElementGeomType(el)
        elemTypeNbr = int(invElemType.get(str(elemType)))
        fileElements.write("%d %d %d" %(invElemIDs[el-1]-NbBoundaryElems, \
                                        elemGrp[invElemIDs[el-1]-1],elemTypeNbr))
        for nid in mesh.GetElemNodes(el):
            fileElements.write(" %d" %(nid))
        fileElements.write("\n")

    fileElements.flush()
    fileElements.close()

    # mesh.boundary
    for el in elemIDs[:NbBoundaryElems]:
        elemType = mesh.GetElementGeomType(el)
        elemTypeNbr = int(invElemType.get(str(elemType)))

        x,y,z = mesh.BaryCenter( el )
        parents = mesh.FindElementsByPoint( x,y,z, bodyType )
        # Modified due to Warning, by Akizuki 2025-03 
        # original: if len(parents) is 2 and elemTypeNbr is not 202:
        if len(parents) == 2 and elemTypeNbr != 202:
            fileBoundary.write("%d %d %d %d %d" \
                %(invElemIDs[el-1],elemGrp[invElemIDs[el-1]-1], \
                  invElemIDs[parents[0]-1]-NbBoundaryElems, \
                  invElemIDs[parents[1]-1]-NbBoundaryElems,elemTypeNbr))
        else:
            fileBoundary.write("%d %d %d 0 %d" \
                %(invElemIDs[el-1],elemGrp[invElemIDs[el-1]-1], \
                  invElemIDs[parents[0]-1]-NbBoundaryElems,elemTypeNbr))

        for nid in mesh.GetElemNodes(el):
            fileBoundary.write(" %d" %(nid))
        fileBoundary.write("\n")

    fileBoundary.flush()
    fileBoundary.close()

    print("Done exporting!\n")
    print("Total time: %0.f s\n" %(time.time()-tstart))


def findSelectedMeshes():
    meshes=list()
    smesh = smeshBuilder.New()
    nrSelected=salome.sg.SelectedCount() # total number of selected items
    
    foundMesh=False
    for i in range(nrSelected):
        selected=salome.sg.getSelected(i)
        selobjID=salome.myStudy.FindObjectID(selected)
        selobj=selobjID.GetObject()
        if selobj.__class__ == SMESH._objref_SMESH_Mesh or \
                selobj.__class__ == salome.smesh.smeshBuilder.meshProxy :
            mName=selobjID.GetName().replace(" ","_")
            foundMesh=True
            mesh=smesh.Mesh(selobj)
            meshes.append(mesh)

    if not foundMesh:
        print("ERROR: Mesh is not selected")
        return list()
    else:
        return meshes


def main():
    u""" 
    Export selected meshes
    """
    meshes=findSelectedMeshes()
    for mesh in meshes:
        if not mesh == None:
            mName=mesh.GetName()
            outdir=os.getcwd()+"/"+mName
            print("Exporting mesh to " + outdir + "\n")
            exportToElmer(mesh,outdir)
            
    
if __name__ == "__main__":
    main()
```

#### (2) Gmshにより変換する方法
- SalomeからMed形式でエクスポートする
- GmshからGmsh形式（Ver.2, ASCII）形式でエクスポートする（Ver.4形式では、拡張されたデータ？が入るようで、形状が化けました）
- ElmerGrid 14 2 Mesh_1.msh -autoclean によりElmerSolver形式に変換する

#### (3) Elmerが対応しているMeshの種類
ElmerSolver ManualのAppendix Dに記載されています
[ElmerSolver Manual](https://www.nic.funet.fi/pub/sci/physics/elmer/doc/ElmerSolverManual.pdf)
![Elmer-Element-Types.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4014034/6b88b01c-57fd-463a-9de7-58f6aabc1c12.png "Salome, Gmshと同等")

### 4. PrePoMax
- SalomeからMed形式でエクスポートする
- GmshからGmsh形式（Ver.2, ASCII）形式でエクスポートする
- 上記のファイルを再度Gmshにより取り込み、inp形式でエクスポートする

![PrePoMax-Pyramid.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4014034/249792eb-3ca4-4ffd-8eb4-49bc8e08e36a.png)

ピラミッド要素はExperimentalの位置づけです。
![PrePomax-Manual-Pyramids.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4014034/d753161c-28d3-4184-a9de-c9faa59438a7.png)

## 計算精度
電界強度は、|Grad(Potential)|です。微分ですので、誤差が誇張されます。

下図は左から、温度（電場解析でのPotential相当）、熱流速（電束密度相当）、Gradient（電界強度相当）、温度では分からない模様が明確に出ます。

![void-pot-flux-grad.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4014034/2b37769e-8f52-4b96-8fe6-f74b42e873fd.png)

誘電率5のエポキシ樹脂中の球状ボイドの電界強度は、平均電界の1.3636倍となります。
今回使用したCAEソフトの計算結果は、1.363～1.364倍で理論値通りです。

![誘電体中のボイドのフィールドマップ.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4014034/943def8b-96b6-42d8-b90d-7b195fb81f8e.png)


以下は、x方向とボイド周辺の電界強度をプロットしたものです。散布図の線が完全に重なりましたので、見やすいように上・右方法にプロットをシフトさせています。
![ボイド電界強度の比較グラフ.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4014034/a79379e5-6339-49b8-9c24-c87f1571a198.png)