From 25ad1488f5c617a2f9d24aab1af93db76d372595 Mon Sep 17 00:00:00 2001
From: GudniRos <gkrosen@gmail.com>
Date: Wed, 2 Dec 2015 19:20:21 -0800
Subject: [PATCH] Added a function to read UBC octree mesh. Updated __init__ to
 import the new functions.

---
 SimPEG/Utils/__init__.py  |  2 +-
 SimPEG/Utils/meshutils.py | 66 +++++++++++++++++++++++++++++++++++++--
 2 files changed, 65 insertions(+), 3 deletions(-)

diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py
index 5280ae79..4d5c8e2e 100644
--- a/SimPEG/Utils/__init__.py
+++ b/SimPEG/Utils/__init__.py
@@ -1,6 +1,6 @@
 from matutils import *
 from codeutils import *
-from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile
+from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCocTreeFiles, readUBCocTreeFiles, writeUBCTensorModel, readVTRFile, writeVTRFile, writeVTUFile
 from curvutils import volTetra, faceInfo, indexCube
 from interputils import interpmat
 from ipythonutils import easyAnimate as animate
diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py
index 81b15566..3537f0b5 100644
--- a/SimPEG/Utils/meshutils.py
+++ b/SimPEG/Utils/meshutils.py
@@ -205,11 +205,13 @@ def writeUBCTensorModel(fileName, mesh, model):
 
     np.savetxt(fileName, modelMatTR.ravel())
 
-def writeUBCocTreeFiles(fileName,mesh,modelDict=None,):
+def writeUBCocTreeFiles(fileName,mesh,modelDict=None):
     '''
         Write UBC ocTree mesh and model files from a simpeg ocTree mesh and model.
 
-        :param
+        :param str fileName: File to write to
+        :param simpeg.Mesh.TreeMesh mesh: The mesh
+        :param dictionary modelDict: The models in a dictionary, where the keys is the name of the of the model file
 
     '''
 
@@ -256,6 +258,66 @@ def writeUBCocTreeFiles(fileName,mesh,modelDict=None,):
             # Save the data
             np.savetxt(item[0],item[1][ubcReorder],fmt='%3.5e')
 
+def readUBCocTreeFiles(meshFile,modelFiles=None):
+    """
+        Read UBC 3D OcTree mesh and/or modelFiles
+
+        Input:
+        :param str meshFile: path to the UBC GIF OcTree mesh file to read
+        :param list of str modelFiles: list of paths modelFiles
+
+        Output:
+        :return SimPEG.Mesh.TreeMesh mesh: The octree mesh
+        :return list of ndarray's: models as a list of numpy array's
+    """
+
+    ## Read the file lines
+    fileLines = np.genfromtxt(meshFile,dtype=str,delimiter='\n')
+    # Extract the data
+    nCunderMesh = np.array(fileLines[0].split(),dtype=float)
+    # I think this is the case?
+    if np.unique(nCunderMesh).size >1:
+        raise Exception('SimPEG TreeMeshes have the same number of cell in all directions')
+    tswCorn = np.array(fileLines[1].split(),dtype=float)
+    smallCell = np.array(fileLines[2].split(),dtype=float)
+    nrCells = np.array(fileLines[3].split(),dtype=float)
+    # Read the index array
+    indArr = np.genfromtxt(fileLines[4::],dtype=np.int)
+
+    ## Calculate simpeg parameters
+    h1,h2,h3 = [np.ones(nr)*sz for nr,sz in zip(nCunderMesh,smallCell)]
+    x0 = tswCorn - np.sum(h3)
+    # Need to convert the index array to a points list that complies with SimPEG TreeMesh.
+    # Shift to start at 0
+    simpegCellPt = indArr[:,0:-1].copy() - np.array([1.,1.,1.])
+    # Need reindex the z index to be from the bottom-left-close corner and to be from the global bottom.
+    simpegCellPt[:,2] = ( nCunderMesh[-1] + 2) - (simpegCellPt[:,2] - indArr[:,3])
+    # Figure out the reordering
+    simpegReorder = np.argsort(simpegCellPt.view(','.join(3*['float'])),axis=0,order=['f2','f1','f0'])[:,0]
+
+    # Calculate the cell level
+    simpegLevel = np.log2(np.min(nCunderMesh)) - np.log2(indArr[:,3])
+    # Make a pointer matrix
+    simpegPointers = np.concatenate((simpegCellPt[simpegReorder,:],simpegLevel[simpegReorder].reshape((-1,1))),axis=1)
+    # Make an index set
+
+    ## Make the tree mesh
+    from SimPEG.Mesh import TreeMesh
+    mesh = TreeMesh([h1,h2,h3],x0)
+    mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()])
+
+    if modelFiles is None:
+        return mesh
+    else:
+        modList = []
+        for modFile in modelFiles:
+            modArr = np.loadtxt(modFile)
+            if len(modArr.shape) == 1:
+                modList.append(modArr[simpegReorder])
+            else:
+                modList.append(modArr[simpegReorder,:])
+        return mesh, modList
+
 def readVTRFile(fileName):
     """
         Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model