notebook.community

Edit and run



In [1]:

    
# Import 
import SimPEG as simpeg, telluricpy
import simpegViz



In [2]:

    
ls Geological_model/









    



background.vtp          PK1_extended.vtp  topoback_vol.vtp
background.vtu          PK1.ts            topo_VOI.vtp
CDED_Lake_Coarse.ts     PK1.vtp           topo_VOI.vtu
CDED_Lake_Coarse.vtp    PK1.vtu           TsSurf_VTK_2_3DModel.py
FWR_Mag.dat             PK2.ts            viewGeology.pvsm
GKR_gocadUtils.py       PK2.vtp           VK.ts
GKR_gocadUtils.pyc      PK3.ts            VK.vtp
HK1.ts                  PK3.vtp           VOI_cliptopo.vtp
HK1.vtp                 testMesh.vtr      VOI.vtp
MEsh_TEst_extended.msh  Till.ts           VTKout.dat
MEsh_TEst.msh           till_VOI.vtu      XVK.ts
newSubset.msh           till_vol.vtp      XVK.vtp
newtopo.vtp             Till.vtp
PK1_Extended.ts         TKCtopo.dat



In [13]:

    
# 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 [14]:

    
# 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 [15]:

    
delFilt = vtk.vtkDelaunay2D()
delFilt.SetInputData(topoPolyData)
delFilt.Update()
delPoly = delFilt.GetOutput()



In [16]:

    
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 [17]:

    
telluricpy.vtkTools.io.writeVTPFile('Geological_model/newtopo.vtp',topoNorm)



In [18]:

    
# 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 [9]:

    
# 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 [22]:

    
# Test adding points for the topo and till together and make a surface from it.
# Get the till points
tillpts = npsup.vtk_to_numpy(till.GetPoints().GetData())



In [21]:

    
pts[:,2].max()









    Out[21]:





450.0418701171875



In [ ]:

    
# Make the surface
# Test deleuny
delTTFilt = vtk.vtkDelaunay2D()
delTTFilt.SetInputData(topoTillPts)
delTTFilt.Update()
delPoly = delTTFilt.GetOutput()



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 [9]:

    
# Clip the topo with the VOI
tempVTP = telluricpy.vtkTools.extraction.clipDataSetWithPolygon(topoNorm,VOI)
telluricpy.vtkTools.io.writeVTPFile('Geological_model/topo_VOI.vtp',tempVTP)



In [9]:

    
# telluricpy.vtkTools.io.writeVTPFile('Geological_model/topoback_vol.vtp',tempVTP)



In [20]:

    
tillVTP = telluricpy.vtkTools.polydata.join2Polydata(tempVTP,till,'upper','lower')
telluricpy.vtkTools.io.writeVTPFile('Geological_model/till_vol.vtp',tillVTP)









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-20-1f8afb0dc954> in <module>()
----> 1 tillVTP = telluricpy.vtkTools.polydata.join2Polydata(tempVTP,till,'upper','lower')
      2 telluricpy.vtkTools.io.writeVTPFile('Geological_model/till_vol.vtp',tillVTP)

NameError: name 'tempVTP' is not defined



In [8]:

    
# Make the till polygons
tempVTP = telluricpy.vtkTools.polydata.join2Polydata(tillVT,topoNorm,'lower','lower')



In [12]:

    
appFilt = vtk.vtkAppendPolyData()
for vtpFile in [PK1,HK1,VK]:
    appFilt.AddInputData(vtpFile)
appFilt.Update()
unstructGeo = telluricpy.vtkTools.extraction.vtp2vtuPolyhedron(appFilt.GetOutput())



In [13]:

    
# Make till
temp = telluricpy.vtkTools.polydata.join2Polydata(backVTP,PK1,'upper','upper')

telluricpy.vtkTools.io.writeVTPFile('backMinusPK1.vtp',temp)



In [14]:

    
# 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 [15]:

    
# Add data array
telluricpy.vtkTools.dataset.addNPDataArrays(unstructGeo,{'GeoType':np.arange(unstructGeo.GetNumberOfCells())},'Cell')
unstructGeo









    Out[15]:





(vtkUnstructuredGrid)0x7f1f7af943b0



In [ ]:



In [17]:

    
# 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 [19]:

    
telluricpy.vtkTools.io.writeVTUFile('BackPKHKVK1.vtu',temp,compress=False)



In [11]:

    
telluricpy.vtkTools.io.writeVTUFile('Geological_model/PK1.vtu',telluricpy.vtkTools.extraction.makePolyhedronCell(PK1,returnGrid=True),compress=False)



In [29]:

    
unstructGeo.GetNumberOfCells()









    Out[29]:





9L



In [30]:

    
telluricpy.vtkTools.extraction.makePolyhedronCell(backVTP,returnGrid=True).GetNumberOfCells()









    Out[30]:





1L



In [ ]:

    
# Convert polydata to polyhydron's