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

Posted at 2025-03-27

要旨

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

CAEソフトの対応状況

1. Salome-meca

Salome-mecaのSalomeにより作成したメッシュを用いて、メッシュの不具合により計算を失敗したことがありません。
内側から三角、四角、三角の平面メッシュを作り、Y軸周りに回転させたモデルです。
四面体、プリズム、ピラミッド、六面体の複合になりました。

2. Gmsh + GetDP (= Onelab)

SalomeのMED形式に対応していますので、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は二次要素に変換したい場合

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()