From 94b8776ac2f58a2fd82547622b2d7f1ce3746d57 Mon Sep 17 00:00:00 2001 From: Gudni UBC-Talva Date: Thu, 5 Feb 2015 12:57:10 -0800 Subject: [PATCH 1/5] Added a function to write VTK rectilinear file (vtr) for a Tensor mesh and model (optional) for view in an external viewer. --- SimPEG/Utils/meshutils.py | 56 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 55 insertions(+), 1 deletion(-) diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 42c9d479..b953f6c2 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -142,7 +142,6 @@ def readUBCTensorMesh(fileName): tensMsh = Mesh.TensorMesh([h1,h2,h3],x0) return tensMsh - def readUBCTensorModel(fileName, mesh): """ ReadUBC 3DTensor mesh model and generate 3D Tensor mesh model in simpegTD @@ -198,6 +197,61 @@ def writeUBCTensorModel(mesh, model, fileName): np.savetxt(fileName, modelMatTR.ravel()) +def writeVTRFile(mesh,model=None,fileName): + """ + Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model. + + Input: + :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK + :param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells + :param model, path to the output vtk file + + """ + + # 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 = vtk.vtkRectilinearGrid() + vtkObj.SetDimensions(xD,yD,zD) + vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1)) + vtkObj.SetYCoordinates(npsup.numpy_to_vtk(vY,deep=1)) + vtkObj.SetZCoordinates(npsup.numpy_to_vtk(vZ,deep=1)) + + # Assign the model('s) to the object + for item in model.iteritems(): + # Convert numpy array + vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtkObj.GetCellData().AddArray(vtkDoubleArr) + + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + vtkObj.Update() + + # Write the file. + # Check the extension of the fileName + os.path.splitext(fileName)[1] + if ext is '': + fileName = fileName + '.vtr' + elif ext is not '.vtr': + raise IOError('{:s} is an incorrect extension, has to be .vtr') + vtrWriteFilter = vtk.vtkXMLRectilinearGridWriter() + vtrWriteFilter.SetInput(vtkObj) + vtrWriteFilter.SetFileName(fileName) + vtrWriteFilter.Update() + + + if __name__ == '__main__': from SimPEG import Mesh From a119ba30834e4fe17b181539bb672fff2483fe7d Mon Sep 17 00:00:00 2001 From: Gudni UBC-Talva Date: Thu, 5 Feb 2015 13:05:08 -0800 Subject: [PATCH 2/5] Added readVTRFile to import a VTK rectilinear XML model file as SimPEG tensor mesh and models --- SimPEG/Utils/meshutils.py | 68 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 64 insertions(+), 4 deletions(-) diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index b953f6c2..f11ff8cb 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -197,6 +197,62 @@ def writeUBCTensorModel(mesh, model, fileName): np.savetxt(fileName, modelMatTR.ravel()) + +def readVTRFile(vtrFileName): + """ + Read VTK vtr (Rectilinear xml file) and return SimPEG Tensor mesh + + Input: + :param vtrFileName, path to the vtr model file + + Output: + :param SimPEG TensorMesh object + :param SimPEG model dictionary + :return + """ + # Import + from vtk import vtkXMLRectilinearGridReader as vtrFileReader + from vtk.util.numpy_support import vtk_to_numpy + + # Read the file + vtrReader = vtrFileReader() + vtrReader.SetFileName(vtrFileName) + vtrReader.Update() + vtrGrid = vtrReader.GetOutput() + # Sort information + hx = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetXCoordinates()))) + xR = vtk_to_numpy(vtrGrid.GetXCoordinates())[0] + hy = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetYCoordinates()))) + yR = vtk_to_numpy(vtrGrid.GetYCoordinates())[0] + zD = np.diff(vtk_to_numpy(vtrGrid.GetZCoordinates())) + # Check the direction of hz + if np.all(zD < 0): + hz = np.abs(zD[::-1]) + zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[-1] + else: + hz = np.abs(zD) + zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[0] + x0 = np.array([xR,yR,zR]) + + # Make the SimPEG object + from SimPEG import Mesh + tensMsh = Mesh.TensorMesh([hx,hy,hz],x0) + + # Grap the models + modelDict = {} + for i in np.arange(vtrGrid.GetCellData().GetNumberOfArrays()): + modelName = vtrGrid.GetCellData().GetArrayName(i) + if np.all(zD < 0): + modFlip = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i)) + tM = tensMsh.r(modFlip,'CC','CC','M') + modArr = tensMsh.r(tM[:,:,::-1],'CC','CC','V') + else: + modArr = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i)) + modelDict[modelName] = modArr + + # Return the data + return tensMsh, modelDict + def writeVTRFile(mesh,model=None,fileName): """ Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model. @@ -207,6 +263,9 @@ def writeVTRFile(mesh,model=None,fileName): :param model, path to the output vtk file """ + # Import + from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter + from vtk.util.numpy_support import vtk_to_numpy # Deal with dimensionalities if mesh.dim >= 1: @@ -222,7 +281,7 @@ def writeVTRFile(mesh,model=None,fileName): zD = mesh.nNz # Use rectilinear VTK grid. # Assign the spatial information. - vtkObj = vtk.vtkRectilinearGrid() + vtkObj = rectGrid() vtkObj.SetDimensions(xD,yD,zD) vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1)) vtkObj.SetYCoordinates(npsup.numpy_to_vtk(vY,deep=1)) @@ -234,18 +293,19 @@ def writeVTRFile(mesh,model=None,fileName): vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) vtkDoubleArr.SetName(item[0]) vtkObj.GetCellData().AddArray(vtkDoubleArr) - + # Set the active scalar vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) vtkObj.Update() - # Write the file. + # Check the extension of the fileName os.path.splitext(fileName)[1] if ext is '': fileName = fileName + '.vtr' elif ext is not '.vtr': raise IOError('{:s} is an incorrect extension, has to be .vtr') - vtrWriteFilter = vtk.vtkXMLRectilinearGridWriter() + # Write the file. + vtrWriteFilter = rectWriter() vtrWriteFilter.SetInput(vtkObj) vtrWriteFilter.SetFileName(fileName) vtrWriteFilter.Update() From 529fed6e6efe2c71db993cf43e2d1a6cccbd4bc3 Mon Sep 17 00:00:00 2001 From: SEOGI KANG Date: Thu, 5 Feb 2015 13:59:50 -0800 Subject: [PATCH 3/5] Modify some bugs of vtk import --- SimPEG/Utils/meshutils.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index f11ff8cb..21c47642 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -3,6 +3,7 @@ from scipy import sparse as sp from matutils import mkvc, ndgrid, sub2ind, sdiag from codeutils import asArray_N_x_Dim from codeutils import isScalar +import os def exampleLrmGrid(nC, exType): assert type(nC) == list, "nC must be a list containing the number of nodes" @@ -169,7 +170,7 @@ def writeUBCTensorMesh(mesh, fileName): s += '%i %i %i\n' %tuple(mesh.vnC) origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated. origin.dtype = float - + s += '%.2f %.2f %.2f\n' %tuple(origin) s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx) s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy) @@ -230,7 +231,7 @@ def readVTRFile(vtrFileName): hz = np.abs(zD[::-1]) zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[-1] else: - hz = np.abs(zD) + hz = np.abs(zD) zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[0] x0 = np.array([xR,yR,zR]) @@ -253,7 +254,7 @@ def readVTRFile(vtrFileName): # Return the data return tensMsh, modelDict -def writeVTRFile(mesh,model=None,fileName): +def writeVTRFile(mesh,fileName,model=None): """ Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model. @@ -263,9 +264,9 @@ def writeVTRFile(mesh,model=None,fileName): :param model, path to the output vtk file """ - # Import + # Import from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter - from vtk.util.numpy_support import vtk_to_numpy + from vtk.util.numpy_support import numpy_to_vtk # Deal with dimensionalities if mesh.dim >= 1: @@ -283,26 +284,26 @@ def writeVTRFile(mesh,model=None,fileName): # Assign the spatial information. vtkObj = rectGrid() vtkObj.SetDimensions(xD,yD,zD) - vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1)) - vtkObj.SetYCoordinates(npsup.numpy_to_vtk(vY,deep=1)) - vtkObj.SetZCoordinates(npsup.numpy_to_vtk(vZ,deep=1)) + 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 for item in model.iteritems(): # Convert numpy array - vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) + vtkDoubleArr = numpy_to_vtk(item[1],deep=1) vtkDoubleArr.SetName(item[0]) vtkObj.GetCellData().AddArray(vtkDoubleArr) # Set the active scalar vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) vtkObj.Update() - + # Check the extension of the fileName - os.path.splitext(fileName)[1] + ext = os.path.splitext(fileName)[1] if ext is '': fileName = fileName + '.vtr' - elif ext is not '.vtr': + elif ext not in '.vtr': raise IOError('{:s} is an incorrect extension, has to be .vtr') # Write the file. vtrWriteFilter = rectWriter() From c6c2713461343d9684c0d79bdaf18166df21cf7a Mon Sep 17 00:00:00 2001 From: Gudni UBC-Talva Date: Thu, 5 Feb 2015 14:56:47 -0800 Subject: [PATCH 4/5] Fixed the __init__.py for meshutils to include read/writeVTRFile functions. Fixed order of inputs in write mesh functions to all be the same. --- SimPEG/Utils/__init__.py | 2 +- SimPEG/Utils/meshutils.py | 19 ++++++++++--------- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 85effc55..6b04e042 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, writeUBCTensorMesh, writeUBCTensorModel +from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile from lrmutils 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 21c47642..98627f34 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -158,12 +158,13 @@ def readUBCTensorModel(fileName, mesh): return model -def writeUBCTensorMesh(mesh, fileName): +def writeUBCTensorMesh(fileName, mesh): """ Writes a SimPEG TensorMesh to a UBC-GIF format mesh file. - :param simpeg.Mesh.TensorMesh mesh: The mesh :param str fileName: File to write to + :param simpeg.Mesh.TensorMesh mesh: The mesh + """ assert mesh.dim == 3 s = '' @@ -179,14 +180,14 @@ def writeUBCTensorMesh(mesh, fileName): f.write(s) f.close() -def writeUBCTensorModel(mesh, model, fileName): +def writeUBCTensorModel(fileName, mesh, model): """ Writes a model associated with a SimPEG TensorMesh to a UBC-GIF format model file. + :param str fileName: File to write to :param simpeg.Mesh.TensorMesh mesh: The mesh :param numpy.ndarray model: The model - :param str fileName: File to write to """ # Reshape model to a matrix @@ -199,9 +200,9 @@ def writeUBCTensorModel(mesh, model, fileName): np.savetxt(fileName, modelMatTR.ravel()) -def readVTRFile(vtrFileName): +def readVTRFile(fileName): """ - Read VTK vtr (Rectilinear xml file) and return SimPEG Tensor mesh + Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model Input: :param vtrFileName, path to the vtr model file @@ -217,7 +218,7 @@ def readVTRFile(vtrFileName): # Read the file vtrReader = vtrFileReader() - vtrReader.SetFileName(vtrFileName) + vtrReader.SetFileName(fileName) vtrReader.Update() vtrGrid = vtrReader.GetOutput() # Sort information @@ -254,14 +255,14 @@ def readVTRFile(vtrFileName): # Return the data return tensMsh, modelDict -def writeVTRFile(mesh,fileName,model=None): +def writeVTRFile(fileName,mesh,model=None): """ Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model. Input: + :param str, path to the output vtk file :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK :param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells - :param model, path to the output vtk file """ # Import From 234cad848b376887d05e9c5c0351cc6e519536e2 Mon Sep 17 00:00:00 2001 From: Gudni UBC-Talva Date: Thu, 5 Feb 2015 15:09:37 -0800 Subject: [PATCH 5/5] Updated the function help text --- SimPEG/Utils/meshutils.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 98627f34..85316869 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -145,8 +145,14 @@ def readUBCTensorMesh(fileName): def readUBCTensorModel(fileName, mesh): """ - ReadUBC 3DTensor mesh model and generate 3D Tensor mesh model in simpegTD + Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg + Input: + :param fileName, path to the UBC GIF mesh file to read + :param mesh, TensorMesh object, mesh that coresponds to the model + + Output: + :return numpy array, model with TensorMesh ordered """ f = open(fileName, 'r') model = np.array(map(float, f.readlines())) @@ -205,12 +211,12 @@ def readVTRFile(fileName): Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model Input: - :param vtrFileName, path to the vtr model file + :param vtrFileName, path to the vtr model file to write to Output: - :param SimPEG TensorMesh object - :param SimPEG model dictionary - :return + :return SimPEG TensorMesh object + :return SimPEG model dictionary + """ # Import from vtk import vtkXMLRectilinearGridReader as vtrFileReader