In [2]:
# Import
import SimPEG as simpeg, telluricpy
import simpegViz
In [3]:
def simpeg2vtk(mesh,modDict):
from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter, VTK_VERSION
from vtk.util.numpy_support import numpy_to_vtk
# Deal with dimensionalities
if mesh.dim >= 1:
vX = mesh.vectorNx
xD = mesh.nNx
yD,zD = 1,1
vY, vZ = np.array([0,0])
if mesh.dim >= 2:
vY = mesh.vectorNy
yD = mesh.nNy
if mesh.dim == 3:
vZ = mesh.vectorNz
zD = mesh.nNz
# Use rectilinear VTK grid.
# Assign the spatial information.
vtkObj = rectGrid()
vtkObj.SetDimensions(xD,yD,zD)
vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
# Assign the model('s) to the object
if modDict is not None:
for item in modDict.iteritems():
# Convert numpy array
vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
vtkDoubleArr.SetName(item[0])
vtkObj.GetCellData().AddArray(vtkDoubleArr)
# Set the active scalar
vtkObj.GetCellData().SetActiveScalars(modDict.keys()[0])
return vtkObj
In [4]:
# Load all the data
PK1 = telluricpy.vtkTools.io.readVTPFile('Geological_model/PK1.vtp')
HK1 = telluricpy.vtkTools.io.readVTPFile('Geological_model/HK1.vtp')
VK = telluricpy.vtkTools.io.readVTPFile('Geological_model/VK.vtp')
VOI = telluricpy.vtkTools.polydata.triangulatePolyData(telluricpy.vtkTools.io.readVTPFile('Geological_model/VOI.vtp'))
topo = telluricpy.vtkTools.polydata.cleanPolyData(telluricpy.vtkTools.io.readVTPFile('Geological_model/topo_VOI.vtp'))
till = telluricpy.vtkTools.extraction.vtu2vtp(telluricpy.vtkTools.io.readVTUFile('Geological_model/till_VOI.vtu'))
In [5]:
# Make the background mesh
x0 = np.array([557100.0,7133350.0,-100])
mesh = simpeg.Mesh.TensorMesh([[(25,22)],[(25,24)],[(25,24)]],x0 = x0)
# Mesh
meshVTR = simpeg2vtk(mesh,{'id':np.arange(mesh.nC)})
In [6]:
# Clip the vtr mesh with the topo
groundClip = telluricpy.vtkTools.extraction.clipDataSetWithPolygon(meshVTR,topo,insideOut=False,extractBounds=False)
In [7]:
# telluricpy.vtkTools.dataset.addNPDataArrays(groundClip,{'GeoType':np.ones(groundClip.GetNumberOfCells())*0})
In [8]:
telluricpy.vtkTools.io.writeVTUFile('groundMesh.vtu',groundClip)
In [9]:
# Till volume
tillClip = telluricpy.vtkTools.extraction.clipDataSetWithPolygon(groundClip,till,insideOut=True,extractBounds=False)
In [10]:
telluricpy.vtkTools.dataset.addNPDataArrays(tillClip,{'GeoType':np.ones(groundClip.GetNumberOfCells())*1.})
In [11]:
telluricpy.vtkTools.io.writeVTUFile('tillMesh.vtu',tillClip)
In [12]:
# Background
backClip = telluricpy.vtkTools.extraction.clipDataSetWithPolygon(groundClip,till,insideOut=False,extractBounds=False)
In [13]:
telluricpy.vtkTools.dataset.addNPDataArrays(backClip,{'GeoType':np.ones(groundClip.GetNumberOfCells())*0.})
In [14]:
telluricpy.vtkTools.io.writeVTUFile('backMesh.vtu',backClip)
In [15]:
# Append to a single mesh
appFilt = vtk.vtkAppendFilter()
appFilt.AddInputData(tillClip)
appFilt.AddInputData(backClip)
appFilt.Update()
In [16]:
tillbackMesh = appFilt.GetOutput()
In [17]:
telluricpy.vtkTools.io.writeVTUFile('tillbackMesh.vtu',tillbackMesh,compress=False)
In [18]:
# Cut all the diamond structures
temp1 = telluricpy.vtkTools.extraction.clipDataSetWithPolygon(tillbackMesh,PK1,insideOut=False,extractBounds=False)
In [19]:
temp2 = telluricpy.vtkTools.extraction.clipDataSetWithPolygon(temp1,HK1,insideOut=False,extractBounds=False)
tillbackNoStruct = telluricpy.vtkTools.extraction.clipDataSetWithPolygon(temp2,VK,insideOut=False,extractBounds=False)
In [20]:
telluricpy.vtkTools.io.writeVTUFile('tillback_noDiaStMesh.vtu',tillbackNoStruct)
In [21]:
# Make the
appFilt = vtk.vtkAppendFilter()
appFilt.AddInputData(tillbackNoStruct)
for nr, vtpObj in enumerate([PK1,HK1,VK]):
tempObj = telluricpy.vtkTools.extraction.clipDataSetWithPolygon(groundClip,vtpObj,insideOut=True,extractBounds=True)
telluricpy.vtkTools.dataset.addNPDataArrays(tempObj,{'GeoType':np.ones(groundClip.GetNumberOfCells())* (nr+2.)})
appFilt.AddInputData(tempObj)
appFilt.Update()
# All the geo Unitis
geoMesh = appFilt.GetOutput()
# Save
telluricpy.vtkTools.io.writeVTUFile('geoMesh_allUnits.vtu',geoMesh)
In [ ]:
# Something is curupt in the topo structure, restructure the polydata
topopts = npsup.vtk_to_numpy(topo.GetPoints().GetData())
ptopopts = vtk.vtkPoints()
ptopopts.SetData(npsup.numpy_to_vtk(topopts,deep=1))
# Make the cells array of polygons
polys = vtk.vtkCellArray()
poly = vtk.vtkTriangle()
poly.GetPointIds().SetNumberOfIds(3)
for i in range(topo.GetNumberOfCells()):
for j in range(topo.GetCell(i).GetPointIds().GetNumberOfIds()):
poly.GetPointIds().SetId(j,topo.GetCell(i).GetPointId(j))
polys.InsertNextCell(poly)
topoPolyData = vtk.vtkPolyData()
topoPolyData.SetPoints(ptopopts)
# topoPolyData.SetPolys(polys)
topoNorm = telluricpy.vtkTools.polydata.normFilter(topoPolyData)
In [ ]:
delFilt = vtk.vtkDelaunay2D()
delFilt.SetInputData(topoPolyData)
delFilt.Update()
delPoly = delFilt.GetOutput()
In [ ]:
topoExt = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(delPoly,VOI)
topoNorm = telluricpy.vtkTools.polydata.normFilter(delFilt.GetOutput())#(topoExt)
# Reverse the order
reverseFilt = vtk.vtkReverseSense()
reverseFilt.SetInputData(topoNorm);
reverseFilt.ReverseCellsOn();
reverseFilt.ReverseNormalsOn()
reverseFilt.Update()
topoNorm = reverseFilt.GetOutput()
In [ ]:
telluricpy.vtkTools.io.writeVTPFile('Geological_model/newtopo.vtp',topoNorm)
In [ ]:
# Make the background polygon
backVTP = telluricpy.vtkTools.polydata.join2Polydata(VOI,till,'upper','lower')
backVTU = telluricpy.vtkTools.extraction.makePolyhedronCell(backVTP,returnGrid=True)
telluricpy.vtkTools.io.writeVTPFile('Geological_model/background.vtp',backVTP)
telluricpy.vtkTools.io.writeVTUFile('Geological_model/background.vtu',backVTU,compress=False)
In [ ]:
# Clip the topo with the VOI
temptopoVTP = telluricpy.vtkTools.extraction.clipDataSetWithPolygon(topoNorm,VOI)
telluricpy.vtkTools.io.writeVTPFile('Geological_model/topo_VOI.vtp',temptopoVTP)
# Clip the VOI with the Topo
tempVOIVTP = telluricpy.vtkTools.extraction.clipDataSetWithPolygon(VOI,topoNorm,insideOut=False)
telluricpy.vtkTools.io.writeVTPFile('Geological_model/VOI_cliptopo.vtp',tempVOIVTP)
In [ ]:
print topopts[topopts[:,1]<7134300].shape
print topopts.shape
In [ ]:
# Test adding points for the topo and till together and make a surface from it.
# Get the till points
topofpts = topopts[topopts[:,1]<7134300].copy()
tillpts = npsup.vtk_to_numpy(till.GetPoints().GetData())
topotillpts = np.vstack((topofpts,tillpts))
ttpts = vtk.vtkPoints()
ttpts.SetData(npsup.numpy_to_vtk(topotillpts,deep=1))
# Make a polydata
topotillPolyData = vtk.vtkPolyData()
topotillPolyData.SetPoints(ttpts)
In [ ]:
telluricpy.vtkTools.io.writeVTPFile('topotillpts.vtp',telluricpy.dataFiles.XYZtools.makeCylinderPtsVTP(topotillpts,height=5,radius=5))
In [ ]:
# surfFilt = vtk.vtkSurfaceReconstructionFilter()
# surfFilt.SetInputData(topotillPolyData)
# contFilt = vtk.vtkContourFilter()
# contFilt.SetInputConnection(surfFilt.GetOutputPort())
# contFilt.SetValue(0,-0.01)
# reverseFilt = vtk.vtkReverseSense()
# reverseFilt.SetInputConnection(contFilt.GetOutputPort())
# reverseFilt.ReverseCellsOn()
# reverseFilt.ReverseNormalsOn()
# reverseFilt.Update()
# # TopoTill volume
# ttvolVTP = reverseFilt.GetOutput()
# telluricpy.vtkTools.io.writeVTPFile('till_surfRec.vtp',ttvolVTP)
In [ ]:
# Make the surface
# Test deleuny
delTTFilt = vtk.vtkDelaunay2D()
delTTFilt.SetInputData(topotillPolyData)
delTTFilt.Update()
delPoly = delTTFilt.GetOutput()
telluricpy.vtkTools.io.writeVTPFile('till_surfDel.vtp',delPoly)
# Get the polydata and add normals
ttNorm = telluricpy.vtkTools.polydata.normFilter(delPoly)
In [ ]:
# Join the 2D del and the till surface to form a volume vtp
ttVTP = telluricpy.vtkTools.polydata.join2Polydata(ttNorm,till,'lower','lower')
telluricpy.vtkTools.io.writeVTPFile('till_TTvol.vtp',ttVTP)
In [ ]:
telluricpy.vtkTools.io.writeVTUFile('till_TTvol.vtu',telluricpy.vtkTools.extraction.makePolyhedronCell(ttVTP,returnGrid=True),compress=False)
In [ ]:
# Join with the VOI
tillAbove = telluricpy.vtkTools.polydata.join2Polydata(VOI,ttVTP,'lower','lower')
telluricpy.vtkTools.io.writeVTPFile('till_TTVOIvol.vtp',ttVTP)
In [ ]:
# Make the till polygons
tillAbove = telluricpy.vtkTools.polydata.join2Polydata(VOI,till,'lower','lower')
tillvolVTP = telluricpy.vtkTools.polydata.join2Polydata(tillAbove,topoNorm,'upper','lower')
telluricpy.vtkTools.io.writeVTPFile('Geological_model/till_vol.vtp',tillvolVTP)
In [ ]:
# Clip the topo with the VOI
tempVTP = telluricpy.vtkTools.extraction.clipDataSetWithPolygon(topoNorm,VOI)
telluricpy.vtkTools.io.writeVTPFile('Geological_model/topo_VOI.vtp',tempVTP)
In [ ]:
# telluricpy.vtkTools.io.writeVTPFile('Geological_model/topoback_vol.vtp',tempVTP)
In [ ]:
tillVTP = telluricpy.vtkTools.polydata.join2Polydata(tempVTP,till,'upper','lower')
telluricpy.vtkTools.io.writeVTPFile('Geological_model/till_vol.vtp',tillVTP)
In [ ]:
# Make the till polygons
tempVTP = telluricpy.vtkTools.polydata.join2Polydata(tillVT,topoNorm,'lower','lower')
In [ ]:
appFilt = vtk.vtkAppendPolyData()
for vtpFile in [PK1,HK1,VK]:
appFilt.AddInputData(vtpFile)
appFilt.Update()
unstructGeo = telluricpy.vtkTools.extraction.vtp2vtuPolyhedron(appFilt.GetOutput())
In [ ]:
# Make till
temp = telluricpy.vtkTools.polydata.join2Polydata(backVTP,PK1,'upper','upper')
telluricpy.vtkTools.io.writeVTPFile('backMinusPK1.vtp',temp)
In [ ]:
# unstructGeo = telluricpy.vtkTools.extraction.vtp2vtuPolyhedron(appFilt.GetOutput())
appVTU = vtk.vtkAppendFilter()
appVTU.AddInputData(backVTU)
appVTU.AddInputData(unstructGeo)
appFilt.Update()
unstrBackGeo = appFilt.GetOutput()
telluricpy.vtkTools.dataset.addNPDataArrays(unstrBackGeo,{'GeoType':np.arange(unstrBackGeo.GetNumberOfCells())},'Cell')
telluricpy.vtkTools.io.writeVTUFile('BackPKHKVK1.vtu',unstrBackGeo,compress=False)
In [ ]:
# Add data array
telluricpy.vtkTools.dataset.addNPDataArrays(unstructGeo,{'GeoType':np.arange(unstructGeo.GetNumberOfCells())},'Cell')
unstructGeo
In [ ]:
In [ ]:
# Make till
temp = telluricpy.vtkTools.polydata.join2Polydata(backVTP,PK1,'upper','upper')
temp = telluricpy.vtkTools.polydata.join2Polydata(temp,HK1,'upper','upper')
temp = telluricpy.vtkTools.polydata.join2Polydata(temp,VK,'upper','upper')
telluricpy.vtkTools.io.writeVTPFile('backMinusPKHK1VK.vtp',temp)
In [ ]:
telluricpy.vtkTools.io.writeVTUFile('BackPKHKVK1.vtu',temp,compress=False)
In [ ]:
telluricpy.vtkTools.io.writeVTUFile('Geological_model/PK1.vtu',telluricpy.vtkTools.extraction.makePolyhedronCell(PK1,returnGrid=True),compress=False)
In [ ]:
unstructGeo.GetNumberOfCells()
In [ ]:
telluricpy.vtkTools.extraction.makePolyhedronCell(backVTP,returnGrid=True).GetNumberOfCells()
In [ ]:
# Convert polydata to polyhydron's