From c83b460672df464598600f7b793927d26e4a4605 Mon Sep 17 00:00:00 2001
From: Rowan Cockett <rowan@3ptscience.com>
Date: Fri, 29 Apr 2016 11:43:31 -0700
Subject: [PATCH 1/3] Surface to Indices (GoCAD and VTK)

---
 SimPEG/Utils/io_utils.py | 132 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 132 insertions(+)
 create mode 100644 SimPEG/Utils/io_utils.py

diff --git a/SimPEG/Utils/io_utils.py b/SimPEG/Utils/io_utils.py
new file mode 100644
index 00000000..f68ee68c
--- /dev/null
+++ b/SimPEG/Utils/io_utils.py
@@ -0,0 +1,132 @@
+from SimPEG import np, Mesh
+import time as tm
+import vtk, vtk.util.numpy_support as npsup
+import re
+
+def read_GOCAD_ts(tsfile):
+    """
+
+    Read GOCAD triangulated surface (*.ts) file
+    INPUT:
+    tsfile: Triangulated surface
+
+    OUTPUT:
+    vrts : Array of vertices in XYZ coordinates [n x 3]
+    trgl : Array of index for triangles [m x 3]. The order of the vertices
+            is important and describes the normal
+            n = cross( (P2 - P1 ) , (P3 - P1) )
+
+    Author: @fourndo
+
+    """
+
+
+    fid = open(tsfile,'r')
+    line = fid.readline()
+
+    # Skip all the lines until the vertices
+    while re.match('TFACE',line)==None:
+        line = fid.readline()
+
+    line = fid.readline()
+    vrtx = []
+
+    # Run down all the vertices and save in array
+    while re.match('VRTX',line):
+        l_input  = re.split('[\s*]',line)
+        temp = np.array(l_input[2:5])
+        vrtx.append(temp.astype(np.float))
+
+        # Read next line
+        line = fid.readline()
+
+    vrtx = np.asarray(vrtx)
+
+    # Skip lines to the triangles
+    while re.match('TRGL',line)==None:
+        line = fid.readline()
+
+    # Run down the list of triangles
+    trgl = []
+
+    # Run down all the vertices and save in array
+    while re.match('TRGL',line):
+        l_input  = re.split('[\s*]',line)
+        temp = np.array(l_input[1:4])
+        trgl.append(temp.astype(np.int))
+
+        # Read next line
+        line = fid.readline()
+
+    trgl = np.asarray(trgl)
+
+    return vrtx, trgl
+
+def surface2inds(vrtx, trgl, mesh, boundaries=True, internal=True):
+    """"
+    Function to read gocad polystructure file and output indexes of mesh with in the structure.
+
+    """
+    # Adjust the index
+    trgl = trgl - 1
+
+    # Make vtk pts
+    ptsvtk = vtk.vtkPoints()
+    ptsvtk.SetData(npsup.numpy_to_vtk(vrtx,deep=1))
+
+    # Make the polygon connection
+    polys = vtk.vtkCellArray()
+    for face in trgl:
+        poly = vtk.vtkPolygon()
+        poly.GetPointIds().SetNumberOfIds(len(face))
+        for nrv, vert in enumerate(face):
+            poly.GetPointIds().SetId(nrv,vert)
+        polys.InsertNextCell(poly)
+
+    # Make the polydata, structure of connections and vrtx
+    polyData = vtk.vtkPolyData()
+    polyData.SetPoints(ptsvtk)
+    polyData.SetPolys(polys)
+
+    # Make implicit func
+    ImpDistFunc = vtk.vtkImplicitPolyDataDistance()
+    ImpDistFunc.SetInput(polyData)
+
+    # Convert the mesh
+    vtkMesh = vtk.vtkRectilinearGrid()
+    vtkMesh.SetDimensions(mesh.nNx,mesh.nNy,mesh.nNz)
+    vtkMesh.SetXCoordinates(npsup.numpy_to_vtk(mesh.vectorNx, deep=1))
+    vtkMesh.SetYCoordinates(npsup.numpy_to_vtk(mesh.vectorNy, deep=1))
+    vtkMesh.SetZCoordinates(npsup.numpy_to_vtk(mesh.vectorNz, deep=1))
+    # Add indexes
+    vtkInd = npsup.numpy_to_vtk(np.arange(mesh.nC), deep=1)
+    vtkInd.SetName('Index')
+    vtkMesh.GetCellData().AddArray(vtkInd)
+
+    extractImpDistRectGridFilt = vtk.vtkExtractGeometry() # Object constructor
+    extractImpDistRectGridFilt.SetImplicitFunction(ImpDistFunc) #
+    extractImpDistRectGridFilt.SetInputData(vtkMesh)
+
+    if boundaries is True:
+        extractImpDistRectGridFilt.ExtractBoundaryCellsOn()
+
+    else:
+        extractImpDistRectGridFilt.ExtractBoundaryCellsOff()
+
+    if internal=True is True:
+        extractImpDistRectGridFilt.ExtractInsideOn()
+
+    else:
+        extractImpDistRectGridFilt.ExtractInsideOff()
+
+    print "Extracting indices from grid..."
+    # Executing the pipe
+    extractImpDistRectGridFilt.Update()
+
+    # Get index inside
+    insideGrid = extractImpDistRectGridFilt.GetOutput()
+    insideGrid = npsup.vtk_to_numpy(insideGrid.GetCellData().GetArray('Index'))
+
+
+    # Return the indexes inside
+    return insideGrid

From 028a16a45a3f0306e4911f675e29f5acee120de1 Mon Sep 17 00:00:00 2001
From: Rowan Cockett <rowan@3ptscience.com>
Date: Fri, 29 Apr 2016 11:44:42 -0700
Subject: [PATCH 2/3] Syntax bug.

---
 SimPEG/Utils/io_utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/SimPEG/Utils/io_utils.py b/SimPEG/Utils/io_utils.py
index f68ee68c..e46c13db 100644
--- a/SimPEG/Utils/io_utils.py
+++ b/SimPEG/Utils/io_utils.py
@@ -113,7 +113,7 @@ def surface2inds(vrtx, trgl, mesh, boundaries=True, internal=True):
     else:
         extractImpDistRectGridFilt.ExtractBoundaryCellsOff()
 
-    if internal=True is True:
+    if internal is True:
         extractImpDistRectGridFilt.ExtractInsideOn()
 
     else:

From 00db6746d4272dd953773300dc815c0a0ee0008e Mon Sep 17 00:00:00 2001
From: Rowan Cockett <rowan@3ptscience.com>
Date: Fri, 29 Apr 2016 11:50:56 -0700
Subject: [PATCH 3/3] Add a warnign about mesh attributes

---
 SimPEG/Utils/io_utils.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/SimPEG/Utils/io_utils.py b/SimPEG/Utils/io_utils.py
index e46c13db..323feb0c 100644
--- a/SimPEG/Utils/io_utils.py
+++ b/SimPEG/Utils/io_utils.py
@@ -18,6 +18,11 @@ def read_GOCAD_ts(tsfile):
 
     Author: @fourndo
 
+
+    .. note::
+
+        Remove all attributes from the GoCAD surface before exporting it!
+
     """