In [1]:
import numpy as np
import Triangulation as tr
import HierarchicalSmooth as hs
import Base as base
from scipy.sparse import csc_matrix
import time
In [2]:
import plotly as py
from plotly.tools import FigureFactory as FF
from plotly.graph_objs import graph_objs
py.offline.init_notebook_mode()
from scipy.spatial import Delaunay
import Base as base
In [3]:
filePrefix = '/home/siddharth/CMU/intro_to_research/intro/DREAM3D/DREAM3D-6.1.3-Linux-x86_64/Data/Mesh'
P = np.loadtxt( filePrefix + '/SharedVertexList.txt' ).T
tri = np.loadtxt( filePrefix + '/SharedTriList.txt', dtype=int )
nFaceLabels = np.loadtxt( filePrefix + '/FaceLabels.txt', dtype=int )
nType = np.loadtxt( filePrefix + '/NodeType.txt', dtype=int )
In [4]:
BoundaryDict = {} # hashing surface elements takes a few minutes; better now than later!
for i in range( len ( nFaceLabels ) ):
if nFaceLabels[i,0] > nFaceLabels[i,1]:
nFaceLabels[i,0], nFaceLabels[i,1] = nFaceLabels[i,1], nFaceLabels[i,0]
tri[i,0], tri[i,1] = tri[i,1], tri[i,0]
try:
BoundaryDict[ str( nFaceLabels[i] ) ].append( i )
except KeyError:
BoundaryDict[ str( nFaceLabels[i] ) ] = [ i ] # start a new list
nUniqFaces = np.vstack( { tuple( row ) for row in nFaceLabels} )
In [5]:
print 'Number of grain boundaries = %d' % len( BoundaryDict )
print 'Number of unique grains = %d' % np.unique( nUniqFaces ).size
In [ ]:
nGB = 1
# triGB = tri[ np.where( base.ismember( nFaceLabels, nUniqFaces[nGB].reshape( 1, -1 ), 'rows' )[0] )[0], : ]
triGB = tri[ BoundaryDict[ str( nUniqFaces[nGB] ) ], : ]
pointGB = np.unique( triGB )
triGB = np.array( base.ismember( triGB, pointGB )[1] ).reshape( -1, 3 )
xGB = P[ :, pointGB ]
typeGB = nType[ pointGB ]
t0 = time.time()
FB, FBsec = hs.DifferentiateFaces( tr.Triangulation( triGB, xGB ) )
print 'Differentiating faces %s seconds' % str( time.time() - t0 )
In [14]:
%autoreload
nGrain = np.array( [ 182, 146, 23 ] )
nGrainPatches = np.any( base.ismember( nFaceLabels, nGrain )[0], axis=1 )
nGrainFaces = nFaceLabels[ nGrainPatches ]
triGrain = tri[ nGrainPatches, : ]
pointGrain = np.unique( triGrain )
triGrain = np.array( base.ismember( triGrain, pointGrain )[1] ).reshape( -1, 3 )
PGrain = P[ :, pointGrain ]
nGrainType = nType[ pointGrain ]
print 'Patches = %d' % len( triGrain )
print 'nodes = %d' % PGrain.shape[1]
# nFaceLabels[ np.any( nFaceLabels==nGrain, axis=1 ) ]
In [ ]:
fig1 = FF.create_trisurf(x=PGrain[0], y=PGrain[1], z=PGrain[2],
simplices=triGrain,
title="System of %d grains" % nGrain.size,
aspectratio=dict(
x = np.ptp( PGrain[0] ),
y = np.ptp( PGrain[1] ),
z = np.ptp( PGrain[2] )
)
)
py.offline.iplot(fig1 )
In [12]:
%autoreload
xS, bPS = hs.HierarchicalSmooth( PGrain, triGrain, nGrainFaces, nGrainType )
In [13]:
np.savetxt( filePrefix + '/SmoothVertexList.txt', xS.T, fmt='%f' )
np.savetxt( filePrefix + '/SmoothTriList.txt', triGrain, fmt='%d')
np.savetxt( filePrefix + '/BoolSmoothed.txt', bPS, fmt='%d' )
np.savetxt( filePrefix + '/TypeSmoothed.txt', nGrainType, fmt='%d' )
In [6]:
import HierarchicalSmooth2 as hs2
In [7]:
%autoreload
triBad = tri[ BoundaryDict[ '[126 716]' ], : ]
pBad = np.unique( triBad )
triBad = np.array( base.ismember( triBad, pBad )[1] ).reshape( -1, 3 )
pointBad = P[ :, pBad ]
typeBad = nType[ pBad ]
T = tr.Triangulation( triBad, pointBad )
FB, FBsec = hs.DifferentiateFaces( T )
thisFB = FBsec[0]
thisRange = FB[ thisFB[0]:thisFB[1]+1, 0 ]
# pointPerimeter = pointBad[ :, thisRange ]
# fixPerimeter = np.where( ( typeBad[ thisRange ] )%10==4 )[0]
# xS = hs.Smooth( pointPerimeter.T, 1.e-4, 2000,
# nFixed=fixPerimeter
# )[0]
# np.savetxt( filePrefix + '/Perimeter.txt', pointPerimeter, fmt='%f' )
# np.savetxt( filePrefix + '/Fixed.txt', fixPerimeter, fmt='%d' )
# np.savetxt( filePrefix + '/PSmooth.txt', xS.T, fmt='%f' )
xS, PS = hs.HierarchicalSmooth(
xPoints=pointBad,
tri=triBad,
nFaceLabels=np.array( [ 182, 667 ] ).reshape( 1, -1 ).repeat( len( triBad ), axis=0 ),
nNodeType=typeBad,
fThreshold=1.e-5
)
np.savetxt( filePrefix + '/BadVertex.txt', pointBad.T, fmt='%f' )
np.savetxt( filePrefix + '/BadTri.txt', triBad, fmt='%d' )
np.savetxt( filePrefix + '/BadType.txt', typeBad, fmt='%d' )
np.savetxt( filePrefix + '/BadSmooth.txt', xS.T, fmt='%f' )
np.savetxt( filePrefix + '/BadBool.txt', PS, fmt='%d' )
# print triBad.shape
# print pBad.shape
# # print nType[ pBad ]
# T = tr.Triangulation( triBad, pointBad )
# print hs.DifferentiateFaces( T )
In [11]:
fig1 = FF.create_trisurf(x=pointBad[0], y=pointBad[1], z=pointBad[2],
simplices=triBad,
title="Misbehaving boundary",
aspectratio=dict(
x = np.ptp( pointBad[0] ),
y = np.ptp( pointBad[1] ),
z = np.ptp( pointBad[2] )
)
)
py.offline.iplot(fig1 )
In [14]:
T = tr.Triangulation( triBad, pointBad )
print T.freeBoundary()
In [15]:
str( [ 2, 3 ] )
Out[15]: