gmsh のチェックについて

  • 負の体積の要素を検出するプログラム

    gmsh チェックプログラム
     1import numpy as np
 2import os, sys, subprocess
 3
 4# ========================================================= #
 5# ===  check__gmshConnectivity.py                       === #
 6# ========================================================= #
 7
 8def check__gmshConnectivity():
 9
10    inpFile  = "msh/model.msh"
11    
12    # ------------------------------------------------- #
13    # --- [1] convert into elmer format             --- #
14    # ------------------------------------------------- #
15    cmd = "ElmerGrid 14 2 {0}".format( inpFile )
16    print( cmd )
17    subprocess.call( cmd.split() )
18
19    # ------------------------------------------------- #
20    # --- [2] load element / node file              --- #
21    # ------------------------------------------------- #
22    
23    dirname  = "/".join( ( inpFile.split( "/" ) )[:-1] )
24    FileBase = ( ( inpFile.split( "/" ) )[-1] ).split( ".msh" )[0]
25    
26    elemFile = dirname + "/" + FileBase + "/" + "mesh.elements"
27    nodeFile = dirname + "/" + FileBase + "/" + "mesh.nodes"
28
29    with open( elemFile, "r" ) as f:
30        elems = np.loadtxt( f, dtype=np.int64 )
31        elems = elems[:,3:]
32    with open( nodeFile, "r" ) as f:
33        nodes = np.loadtxt( f, dtype=np.float )
34        nodes = nodes[:,2:]
35
36    # ------------------------------------------------- #
37    # --- [3] calculate volume of the element       --- #
38    # ------------------------------------------------- #
39    nElems   = elems.shape[0]
40    onesixth = 1.0 / 6.0
41    volumes  = np.zeros( (nElems) )
42    for iv,vert in enumerate( elems ):
43        iv0,iv1,iv2,iv3 =    vert[0]-1,    vert[1]-1,    vert[2]-1,    vert[3]-1
44        nd0,nd1,nd2,nd3 = nodes[iv0,:], nodes[iv1,:], nodes[iv2,:], nodes[iv3,:]
45        vc1,vc2,vc3     =    nd1 - nd0,    nd2 - nd0,    nd3 - nd0
46        matrix          = np.concatenate( [vc1[:,None],vc2[:,None],vc3[:,None]], axis=1 )
47        volumes[iv]     = onesixth * np.linalg.det( matrix )
48
49    # ------------------------------------------------- #
50    # --- [4] output                                --- #
51    # ------------------------------------------------- #
52    #  -- [4-1] statistics                          --  #
53    avgvol = np.average( volumes )
54    stddev = np.std    ( volumes )
55    print()
56    print( "[check__gmshConnectivity.py]  --- check of the gmsh File --- " )
57    print( "[check__gmshConnectivity.py]    #.of elements     :: {0}".format( nElems ) )
58    print( "[check__gmshConnectivity.py]   average( volumes ) :: {0}".format( avgvol ) )
59    print( "[check__gmshConnectivity.py]    stddev( volumes ) :: {0}".format( stddev ) )
60    print()
61
62    #  -- [4-2] output data                         --  #
63    Data        = np.zeros( (nElems,6) )
64    Data[:,0  ] = np.arange( 1,nElems+1 )
65    Data[:,1:5] = np.copy( elems[:,:] )
66    Data[:,5  ] = np.copy( volumes    )
67
68    #  -- [4-3] negative volume elements            --  #
69    index       = np.where( Data[:,5] <= 0.0 )
70    illegal     = Data[index]
71
72    #  -- [4-4] save in a file                      --  #
73    import nkUtilities.save__pointFile as spf
74    outFile1    = "dat/gmsh_check.dat"
75    outFile2    = "dat/illegal_elements.dat"
76    spf.save__pointFile( outFile=outFile1, Data=Data )
77    spf.save__pointFile( outFile=outFile2, Data=illegal, \
78                         fmt=["%12d","%12d","%12d","%12d","%12d","%15.8e"] )
79    return()
80
81    
82
83# ========================================================= #
84# ===   実行部                                          === #
85# ========================================================= #
86
87if ( __name__=="__main__" ):
88    check__gmshConnectivity()