Gmsh mesh generation engine support

[BinWang0213]

Is your feature request related to a problem? Please describe.
Gmsh mesh generator is support by many open-source simulation code. It will be nice to add a Gmsh meshing interface.

• Also, meshio library support Gmsh very well. Once we have gmsh file, we can use meshio to convert the mesh into any mesh format, such as

Abaqus, ANSYS msh, AVS-UCD, CGNS, DOLFIN XML, Exodus, FLAC3D, H5M, Kratos/MDPA, Medit, MED/Salome, Nastran (bulk data), Neuroglancer precomputed format, Gmsh (format versions 2.2, 4.0, and 4.1), OBJ, OFF, PERMAS, PLY, STL, Tecplot .dat, TetGen .node/.ele, SVG (2D only, output only), UGRID, VTK, VTU (not raw binary data), WKT (TIN), XDMF.

Describe the solution you'd like
The working code is attached, please feel free to integrate it into the MicroStructPy. It support to mesh any polymesh.region we want.

Noted that Gmsh is intefaced by generating its built-in script file *.geo. We need to use Gmsh to load the geo file and generate mesh.

pmesh = msp.meshing.PolyMesh.from_seeds(seeds, domain)

tri_basename = 'polymesh_gmsh.geo'
tri_filename = os.path.join(directory, tri_basename)

Material_number=0
MaterialIDs=[i for i,p in enumerate(pmesh.phase_numbers) if p==Material_number]

#Current function usage
#geo=pgm.Mesh.mspPolymesh2geo(pmesh,regionIDs=MaterialIDs,max_edge_size=0.3)
#with open(tri_filename, 'w') as f:
#    f.write(geo)

#Expected msp usage
tmesh = msp.meshing.TriMesh.from_polymesh(pmesh, phases, max_edge_length=0.3)
tmesh.write(tri_filename, "geo", seeds, pmesh)

def mspPolymesh2geo(polymesh,regionIDs=[],max_edge_size=-1):
    """Convert MicroStructPy polymesh to gmsh geo script

    Method
    --------
    Chose the regionID you want to output

    Author: Bin Wang(binwang.0213@gmail.com)
    Date: April. 2020
    """
    NumRegions=len(polymesh.regions)
    if(len(regionIDs)==0): regionIDs=[n for n in range(NumRegions)] #All region
    #Collect Pts,facets and volumes for requested material
    facets_out=[]
    regionIDs_out=[]
    for ri,region in enumerate(polymesh.regions):
        for mID in regionIDs:
            if(ri==mID):
                facets_out+=region
                regionIDs_out+=[ri]
    pts_out=list(set([pi for fi in facets_out for pi in polymesh.facets[fi]]))

    string=''
    string += '/* -------Genetrated by PyGeoMesh------- */\n\n'
    string += 'SetFactory("Built-in");\n' #This is not working well for fixing orientation
    string +='\n'

    #* 1. -----------Collect Points------------
    string += '/* --Model_Points-- */\n'
    for pi, pts in enumerate(polymesh.points):#Loop over regions
        if(pi not in pts_out): continue
        string += "p%d = newp; Point(p%d) = {%.20f,%.20f,%.20f};\n" % (pi,pi, pts[0], pts[1], pts[2])
    string +='\n'

    #* 2. -----------Collect Polyhedron Facet Edges------------
    Line_id=0
    string += '/* --Polyhedron Edges/Facets-- */\n'
    for fi,facet in enumerate(polymesh.facets):
        if(fi not in facets_out): continue
        NumPts=len(facet)
        temp ="clp%d = newll; Curve Loop(clp%d) = {" % (fi,fi)
        for pi in range(NumPts):
            id0,id1 = facet[pi],facet[(pi+1)%NumPts] #round end connect pair
            string += "lp%d = newl; Line(lp%d) = {p%d,p%d};\n" % (Line_id, Line_id, id0, id1)
            temp += "lp%d," % (Line_id)
            Line_id = Line_id + 1
        temp=temp[:-1] + "};\n"
        #handle boundary inward normal vector
        temp += "sp%d = news;  Plane Surface(sp%d) = {clp%d};\n" % (fi,fi,fi)
        string += temp + '\n'
    string += '\n'

    #* 3. -----------Collect Polyhedron------------
    string += '/* --Polyhedron Cells-- */\n'
    for ri,region in enumerate(polymesh.regions):
        if(ri not in regionIDs_out): continue
        string+='vlp%d = newsl; Surface Loop(vlp%d)={' %(ri,ri)
        for fi in region:
            facet_neigh=polymesh.facet_neighbors[fi]
            if(facet_neigh[1]==ri):  string+='-sp%d,'%(fi) #Backside of the facet
            else: string+='sp%d,'%(fi)
        string=string[:-1]+'};\n'
        string+='vp%d = newv; Volume(vp%d)={vlp%d};\n' %(ri,ri,ri)
        string +='\n'

    string += '/* --Combine Duplicated Entites, Polyhedron Edges here-- */\n'
    string+='Coherence;\n'
    
    #* 4. ------------Reverse boundary facet orientation-----------
    inward_facet=[]
    regionIDs_wall=list( set(range(len(polymesh.regions))) - set(regionIDs_out) )
    id2tagName={-1:'X-',-2:'X+',-3:'Y-',-4:'Y+',-5:'Z-',-6:'Z+'}
    boundary_facets={'X-':[],'X+':[],'Y-':[],'Y+':[],'Z-':[],'Z+':[],'Wall':[]}
    for fi,facet in enumerate(polymesh.facets):
        if(fi not in facets_out): continue
        facet_neigh=np.array(polymesh.facet_neighbors[fi])
        check1=facet_neigh[0]<0 or facet_neigh[1]<0 #boundary facet X-,X+...
        check2=facet_neigh[0] in regionIDs_wall or facet_neigh[1] in regionIDs_wall #material interface wall
        #print('facet',fi,'neigh',facet_neigh,check1,check2)
        
        if(check1 or check2):#This is a boundary facet
            if(facet_neigh[1] in regionIDs_out):#Normal always point toward facet_neigh[1] cell
                inward_facet+=[fi]
        if(check1):
            boundary_id=facet_neigh[facet_neigh<0][0]
            boundary_name=id2tagName[boundary_id]
            boundary_facets[boundary_name].append(fi)
        if(check2):
            boundary_facets['Wall'].append(fi)
    string +='\n'

    #print(inward_facet,regionIDs_wall)
    string += '/* --Fix boundary surface orientation-- */\n'
    string += 'ReverseMesh Surface{%s};\n' % ','.join('sp'+str(x) for x in inward_facet)
    string +='\n'

    string += '/* --Physical Tags-- */\n'
    print('boundary facets Tag=',boundary_facets)
    for name,facet_list in boundary_facets.items():
        if(len(facet_list)>0):
            string+='Physical Surface("%s") = {%s};\n' %(name,','.join('sp'+str(f) for f in facet_list))
    
    for ri,region_geo in enumerate(polymesh.regions):
        if(ri not in regionIDs_out): continue
        string+='Physical Volume("Material%d") = {vp%d};\n' %(ri,ri)
    string+='\n'

    string += '/* --Mesh Size Setup-- */\n'
    if(max_edge_size==-1):#Auto mesh size
        max_edge_size=np.sum(polymesh.volumes)**(1/3.0)/8.0

    string +='//Mesh.CharacteristicLengthMin = %s;\n' %(max_edge_size)
    string +='Mesh.CharacteristicLengthMax = %s;\n' %(max_edge_size)

    return string

[kip-hart]

Hi Bin,

I'm planning on inheriting PolyMesh and TriMesh from the meshio.Mesh class in v2.0 of the code. At that point, I'll integrate gmsh support and include pygmsh so the PolyMesh can be saved to a .geo for externally calling gmsh or the TriMesh can be created with pygmsh like TriMesh.from_poly(..., method='gmsh') then written to a .msh file like TriMesh.write(..., format='msh'). That will also allow for the meshes to be written to a variety of formats. I'll update once I get closer to tackling this.

Best,
Kip

[BinWang0213]

Hi Kip,

Sounds good. Please be careful the normal orientation issue. All boundary faces should be outward normal.

You can see my effort @ #16 (comment) to adapt this for Gmsh.

Thanks,
Bin

[kip-hart]

Thanks, that's a good point. The exterior faces can be oriented positive outward, while leaving the interior faces arbitrary.

Since the facets are created by the PolyMesh.from_seeds method (https://github.com/kip-hart/MicroStructPy/blob/v1.2.2/src/microstructpy/meshing/polymesh.py#L603), I can make that modification before adding gmsh.

[kip-hart]

I realize it's been several months, but I've added gmsh to the code MicroStructPy code. It can be specified using Trimesh.from_polymesh(..., mesher='gmsh'). I'm planning on a v1.4.0 release soon and at that point it'll be part of the PyPI release. Thanks for the feature suggestion!

[MayaStr]

Hey Kip,

I have been playing around with your package to import a mesh into FENICS using meshio. I was able to create the triangular mesh like you described. Trimesh.from_polymesh(..., mesher='gmsh'). Is there a possibility now to export the mesh with facets, so it can be imported by meshio? What is the workflow? Or use Bin's script?