kip-hart / microstructpy Goto Github PK

seed_basename = 'seeds.vtk'
seed_filename = os.path.join(directory, seed_basename)
#seeds2VTK(seed_filename,seeds)
seeds.write(seed_filename)

poly_basename = 'polymesh.vtk'
poly_filename = os.path.join(directory, poly_basename)
#polymesh2VTK(poly_filename,pmesh)
pmesh.write(poly_filename)

import vtk
import vtk.util.numpy_support as vtk_np

def seeds2VTK(fname, seeds):
    numPts=len(seeds)
    geom_name = type(seeds.seeds[0].geometry).__name__.lower().strip()
    centers=np.array([s.position for s in seeds.seeds])
    
    if(centers.shape[1]==2):#extend to 3D pts
        temp = np.zeros((numPts,3))
        temp[:,:-1] = centers
        centers=temp
    
    vtkPoly=vtk.vtkPolyData()
    if(geom_name in ['circle','sphere']):
        verts = vtk.vtkPoints()
        cells = vtk.vtkCellArray()
        
        verts.SetNumberOfPoints(numPts)
        verts.SetData(vtk_np.numpy_to_vtk(centers))
        vtkPoly.SetPoints(verts)

        cells_npy = np.vstack([np.ones(numPts,dtype=np.int64),
                               np.arange(numPts,dtype=np.int64)]).T.flatten()
        cells.SetCells(numPts,vtk_np.numpy_to_vtkIdTypeArray(cells_npy))
        vtkPoly.SetVerts(cells)

        diameters=np.array([s.geometry.diameter for s in seeds.seeds])
        #print(diameters)
        diameters = vtk_np.numpy_to_vtk(diameters)
        diameters.SetName("Diameters")
        vtkPoly.GetCellData().SetScalars(diameters)
        
        
        MaterialID=np.array([s.phase for s in seeds.seeds])
        MaterialID = vtk_np.numpy_to_vtk(MaterialID)
        MaterialID.SetName("MaterialID")
        vtkPoly.GetCellData().AddArray(MaterialID)

        writer = vtk.vtkPolyDataWriter()
        writer.SetFileName(fname)
        writer.SetInputData(vtkPoly)
        writer.SetFileTypeToASCII() #We can switch to binary for large dataset
        writer.Write() 
    else:
        print('Other shape is not supported to write to VTK yet!')

def polymesh2VTK(fname,pmesh):
    ugrid = vtk.vtkUnstructuredGrid()

    # Create real polyhedron
    points = vtk.vtkPoints()

    points.SetNumberOfPoints(len(pmesh.points))
    pmesh.points=np.array(pmesh.points)
    for i in range(len(pmesh.points)):
        points.SetPoint(i, pmesh.points[i,0],pmesh.points[i,1],pmesh.points[i,2])

    ugrid.SetPoints(points)

    # List of pointIds that make up the face
    for cell in pmesh.regions:
        faceId = vtk.vtkIdList()
        faceId.InsertNextId(len(cell)) # Number faces that make up the cell.
        for face in cell:
            face_pts=pmesh.facets[face]
            faceId.InsertNextId(len(face_pts)) # Number of points in face
            [faceId.InsertNextId(i) for i in face_pts] # Insert the pointIds for the face
        ugrid.InsertNextCell(vtk.VTK_POLYHEDRON, faceId)

    MaterialID = vtk_np.numpy_to_vtk(num_array=np.array(pmesh.phase_numbers))
    MaterialID.SetName("MaterialID")
    ugrid.GetCellData().SetScalars(MaterialID)
    
    SeedsID = vtk_np.numpy_to_vtk(num_array=np.array(pmesh.seed_numbers))
    SeedsID.SetName("SeedsID")
    ugrid.GetCellData().AddArray(SeedsID)

    Vol=vtk_np.numpy_to_vtk(num_array=np.array(pmesh.volumes)) 
    Vol.SetName("Volume")
    ugrid.GetCellData().AddArray(Vol)

    #writer = vtk.vtkXMLUnstructuredGridWriter()
    writer = vtk.vtkUnstructuredGridWriter()
    writer.SetInputData(ugrid)
    writer.SetFileName(fname)
    writer.SetFileTypeToASCII()
    writer.Write()

import microstructpy


position_array = [[-0.37342, -0.289708, 0.5128352],
                  [0.204211, 0.41113, 0.4889689],
                  [-0.453932, -0.178443, 0.4169717],
                  [-0.112715, 0.434528, 0.4959116],
                  [-0.12776, -0.0181222, 0.4907776]]
seed_list = microstructpy.seeding.SeedList()
for position in position_array:
    new_seed = microstructpy.seeding.Seed.factory(seed_type="sphere", position=position)
    seed_list.append(new_seed)

domain = microstructpy.geometry.Cube(side_length=1.2)
polymesh = microstructpy.meshing.PolyMesh.from_seeds(seed_list, domain, verbose=True)

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([email protected])
    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

def getPhaseColor(i,phases):
    return phases[i].get('color', 'C' + str(i % 10))

# Create Figure
k = 1
off = 1
len_x = domain.length + 2 * off
len_y = domain.width + 2 * off
plt.figure(figsize=(k * len_x, k * len_y))

# Plot Seeds
seed_colors=[getPhaseColor(s.phase, phases) for s in seeds]
seeds.plot(color=seed_colors, alpha=0.8, edgecolor='k', linewidth=0.3)
domain.plot(facecolor='none', edgecolor='k', linewidth=0.5)

# Add legend
legend_patchs=[]
for ph in phases:
     legend_patchs.append( mpatches.Patch(color=ph["color"], ec='k',label=ph["name"]) )
plt.legend(handles=legend_patchs)
plt.axis('image')
plt.show()

# Create Figure
plt.figure(figsize=(5, 5))

seeds.plot(material_name=['Matrix','Void'],color=['C0','C1'], alpha=0.8, edgecolor='k', linewidth=0.3)
domain.plot(facecolor='none', edgecolor='k', linewidth=0.5)

plt.axis('image')
plt.show()