From 912c39151f05dacc52d6d48fecb9437e0c980333 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Sat, 2 Nov 2013 07:52:35 -0700 Subject: [PATCH 01/20] Added visulization folder and written some vtkTools. --- SimPEG/__init__.py | 1 + SimPEG/visulization/__init__.py | 1 + SimPEG/visulization/vtkTools.py | 171 ++++++++++++++++++++++++++++++++ 3 files changed, 173 insertions(+) create mode 100644 SimPEG/visulization/__init__.py create mode 100644 SimPEG/visulization/vtkTools.py diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index e8946f88..863c8ca2 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -2,3 +2,4 @@ import mesh import utils import inverse from Solver import Solver +import visulization diff --git a/SimPEG/visulization/__init__.py b/SimPEG/visulization/__init__.py new file mode 100644 index 00000000..1bd4f40a --- /dev/null +++ b/SimPEG/visulization/__init__.py @@ -0,0 +1 @@ +from vtkTools import vtkTools \ No newline at end of file diff --git a/SimPEG/visulization/vtkTools.py b/SimPEG/visulization/vtkTools.py new file mode 100644 index 00000000..2b440d79 --- /dev/null +++ b/SimPEG/visulization/vtkTools.py @@ -0,0 +1,171 @@ +import numpy as np, vtk, vtk.util.numpy_support as npsup, pdb +from SimPEG.utils import mkvc + + +class vtkTools(object): + """ + Class that interacts with VTK visulization toolkit. + + """ + + def __init__(self): + """ Initializes the VTK vtkTools. + + """ + + pass + + @staticmethod + def makeCellVTKObject(mesh,model): + """ + Make and return a cell based VTK object for a simpeg 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 + + Output: + :rtype: vtkRecilinearGrid object + :return: vtkObj + """ + + # Use rectilinear VTK grid. + # Asign the spatial information. + vtkObj = vtk.vtkRectilinearGrid() + vtkObj.SetDimensions(mesh.nNx,mesh.nNy,mesh.nNz) + vtkObj.SetXCoordinates(npsup.numpy_to_vtk(mesh.vectorNx,deep=1)) + vtkObj.SetYCoordinates(npsup.numpy_to_vtk(mesh.vectorNy,deep=1)) + vtkObj.SetZCoordinates(npsup.numpy_to_vtk(mesh.vectorNz,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) + + return vtkObj + + @staticmethod + def makeFaceVTKobject(mesh,model): + """ + Make and return a cell based VTK object for a simpeg 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). + Property array must be order hstack(Fx,Fy,Fz) + + Output: + :rtype: vtkUnstructuredGrid object + :return: vtkObj + """ + + ## Convert simpeg mesh to VTK properties + # Convert mesh nodes to vtkPoints + vtkPts = vtk.vtkPoints() + vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) + + # Define the face "cells" + # Using VTK_QUAD cell for faces (see VTK file format) + nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') + def faceR(mat,length): + return mat.T.reshape((length,1)) + # First direction + nTFx = np.prod(mesh.nFx) + FxCellBlock = np.hstack([ 4*np.ones((nTFx,1),dtype='int64'),faceR(nodeMat[:,:-1,:-1],nTFx),faceR(nodeMat[:,1: ,:-1],nTFx),faceR(nodeMat[:,1: ,1: ],nTFx),faceR(nodeMat[:,:-1,1: ],nTFx)] ) + # Second direction + nTFy = np.prod(mesh.nFy) + FyCellBlock = np.hstack([ 4*np.ones((nTFy,1),dtype='int64'),faceR(nodeMat[:-1,:,:-1],nTFy),faceR(nodeMat[1: ,:,:-1],nTFy),faceR(nodeMat[1: ,:,1: ],nTFy),faceR(nodeMat[:-1,:,1: ],nTFy)] ) + # Third direction + nTFz = np.prod(mesh.nFz) + FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] ) + # Cells -cell array + FCellArr = vtk.vtkCellArray() + FCellArr.SetCells(np.prod(mesh.nF),npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) + # Cell type + FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(np.sum(mesh.nF),dtype='uint8'),deep=1) + # Cell location + FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nF),5,dtype='int64'),deep=1) + + ## Make the object + vtkObj = vtk.vtkUnstructuredGrid() + # Set the objects properties + vtkObj.SetPoints(vtkPts) + vtkObj.SetCells(FCellType,FCellLoc,FCellArr) + + # 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) + + return vtkObj + + + # Simple write model functions. + @staticmethod + def writeVTPFile(fileName,vtkPolyObject): + '''Function to write vtk polydata file (vtp).''' + polyWriter = vtk.vtkXMLPolyDataWriter() + polyWriter.SetInput(vtkPolyObject) + polyWriter.SetFileName(fileName) + polyWriter.Update() + + @staticmethod + def writeVTUFile(fileName,vtkUnstructuredGrid): + '''Function to write vtk unstructured grid (vtu).''' + Writer = vtk.vtkXMLUnstructuredGridWriter() + Writer.SetInput(vtkUnstructuredGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def writeVTRFile(fileName,vtkRectilinearGrid): + '''Function to write vtk rectilinear grid (vtr).''' + Writer = vtk.vtkXMLRectilinearGridWriter() + Writer.SetInput(vtkRectilinearGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def writeVTSFile(fileName,vtkStructuredGrid): + '''Function to write vtk structured grid (vts).''' + Writer = vtk.vtkXMLStructuredGridWriter() + Writer.SetInput(vtkStructuredGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def readVTSFile(fileName): + '''Function to read vtk structured grid (vts) and return a grid object.''' + Reader = vtk.vtkXMLStructuredGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTUFile(fileName): + '''Function to read vtk structured grid (vtu) and return a grid object.''' + Reader = vtk.vtkXMLUnstructuredGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTRFile(fileName): + '''Function to read vtk structured grid (vtr) and return a grid object.''' + Reader = vtk.vtkXMLRectilinearGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTPFile(fileName): + '''Function to read vtk structured grid (vtp) and return a grid object.''' + Reader = vtk.vtkXMLPolyDataReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + From 304ee3ae3278ce3772b29cfdd8541b765f14ae13 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Mon, 4 Nov 2013 08:56:32 -0800 Subject: [PATCH 02/20] Renamed folder. Added vtkView Class that makes vtk figures of SimPEG mesh, models. --- SimPEG/__init__.py | 2 +- SimPEG/visulization/__init__.py | 1 - SimPEG/visulization/vtkTools.py | 171 -------------------------------- 3 files changed, 1 insertion(+), 173 deletions(-) delete mode 100644 SimPEG/visulization/__init__.py delete mode 100644 SimPEG/visulization/vtkTools.py diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 863c8ca2..8107e847 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -2,4 +2,4 @@ import mesh import utils import inverse from Solver import Solver -import visulization +import visulize diff --git a/SimPEG/visulization/__init__.py b/SimPEG/visulization/__init__.py deleted file mode 100644 index 1bd4f40a..00000000 --- a/SimPEG/visulization/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from vtkTools import vtkTools \ No newline at end of file diff --git a/SimPEG/visulization/vtkTools.py b/SimPEG/visulization/vtkTools.py deleted file mode 100644 index 2b440d79..00000000 --- a/SimPEG/visulization/vtkTools.py +++ /dev/null @@ -1,171 +0,0 @@ -import numpy as np, vtk, vtk.util.numpy_support as npsup, pdb -from SimPEG.utils import mkvc - - -class vtkTools(object): - """ - Class that interacts with VTK visulization toolkit. - - """ - - def __init__(self): - """ Initializes the VTK vtkTools. - - """ - - pass - - @staticmethod - def makeCellVTKObject(mesh,model): - """ - Make and return a cell based VTK object for a simpeg 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 - - Output: - :rtype: vtkRecilinearGrid object - :return: vtkObj - """ - - # Use rectilinear VTK grid. - # Asign the spatial information. - vtkObj = vtk.vtkRectilinearGrid() - vtkObj.SetDimensions(mesh.nNx,mesh.nNy,mesh.nNz) - vtkObj.SetXCoordinates(npsup.numpy_to_vtk(mesh.vectorNx,deep=1)) - vtkObj.SetYCoordinates(npsup.numpy_to_vtk(mesh.vectorNy,deep=1)) - vtkObj.SetZCoordinates(npsup.numpy_to_vtk(mesh.vectorNz,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) - - return vtkObj - - @staticmethod - def makeFaceVTKobject(mesh,model): - """ - Make and return a cell based VTK object for a simpeg 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). - Property array must be order hstack(Fx,Fy,Fz) - - Output: - :rtype: vtkUnstructuredGrid object - :return: vtkObj - """ - - ## Convert simpeg mesh to VTK properties - # Convert mesh nodes to vtkPoints - vtkPts = vtk.vtkPoints() - vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) - - # Define the face "cells" - # Using VTK_QUAD cell for faces (see VTK file format) - nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') - def faceR(mat,length): - return mat.T.reshape((length,1)) - # First direction - nTFx = np.prod(mesh.nFx) - FxCellBlock = np.hstack([ 4*np.ones((nTFx,1),dtype='int64'),faceR(nodeMat[:,:-1,:-1],nTFx),faceR(nodeMat[:,1: ,:-1],nTFx),faceR(nodeMat[:,1: ,1: ],nTFx),faceR(nodeMat[:,:-1,1: ],nTFx)] ) - # Second direction - nTFy = np.prod(mesh.nFy) - FyCellBlock = np.hstack([ 4*np.ones((nTFy,1),dtype='int64'),faceR(nodeMat[:-1,:,:-1],nTFy),faceR(nodeMat[1: ,:,:-1],nTFy),faceR(nodeMat[1: ,:,1: ],nTFy),faceR(nodeMat[:-1,:,1: ],nTFy)] ) - # Third direction - nTFz = np.prod(mesh.nFz) - FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] ) - # Cells -cell array - FCellArr = vtk.vtkCellArray() - FCellArr.SetCells(np.prod(mesh.nF),npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) - # Cell type - FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(np.sum(mesh.nF),dtype='uint8'),deep=1) - # Cell location - FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nF),5,dtype='int64'),deep=1) - - ## Make the object - vtkObj = vtk.vtkUnstructuredGrid() - # Set the objects properties - vtkObj.SetPoints(vtkPts) - vtkObj.SetCells(FCellType,FCellLoc,FCellArr) - - # 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) - - return vtkObj - - - # Simple write model functions. - @staticmethod - def writeVTPFile(fileName,vtkPolyObject): - '''Function to write vtk polydata file (vtp).''' - polyWriter = vtk.vtkXMLPolyDataWriter() - polyWriter.SetInput(vtkPolyObject) - polyWriter.SetFileName(fileName) - polyWriter.Update() - - @staticmethod - def writeVTUFile(fileName,vtkUnstructuredGrid): - '''Function to write vtk unstructured grid (vtu).''' - Writer = vtk.vtkXMLUnstructuredGridWriter() - Writer.SetInput(vtkUnstructuredGrid) - Writer.SetFileName(fileName) - Writer.Update() - - @staticmethod - def writeVTRFile(fileName,vtkRectilinearGrid): - '''Function to write vtk rectilinear grid (vtr).''' - Writer = vtk.vtkXMLRectilinearGridWriter() - Writer.SetInput(vtkRectilinearGrid) - Writer.SetFileName(fileName) - Writer.Update() - - @staticmethod - def writeVTSFile(fileName,vtkStructuredGrid): - '''Function to write vtk structured grid (vts).''' - Writer = vtk.vtkXMLStructuredGridWriter() - Writer.SetInput(vtkStructuredGrid) - Writer.SetFileName(fileName) - Writer.Update() - - @staticmethod - def readVTSFile(fileName): - '''Function to read vtk structured grid (vts) and return a grid object.''' - Reader = vtk.vtkXMLStructuredGridReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - - @staticmethod - def readVTUFile(fileName): - '''Function to read vtk structured grid (vtu) and return a grid object.''' - Reader = vtk.vtkXMLUnstructuredGridReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - - @staticmethod - def readVTRFile(fileName): - '''Function to read vtk structured grid (vtr) and return a grid object.''' - Reader = vtk.vtkXMLRectilinearGridReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - - @staticmethod - def readVTPFile(fileName): - '''Function to read vtk structured grid (vtp) and return a grid object.''' - Reader = vtk.vtkXMLPolyDataReader() - Reader.SetFileName(fileName) - Reader.Update() - return Reader.GetOutput() - From f04f70c455aad435c54e1fdfe49cabb7c95a9a62 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Mon, 4 Nov 2013 10:32:32 -0800 Subject: [PATCH 03/20] Committing visulize/vtk folder. Contains: vtkTools class that has abstract methods to deal with the interface between SimPEG and VTK. vtkView class that is s vtk rendering container for SimPEG --- SimPEG/visulize/__init__.py | 2 + SimPEG/visulize/vtk/__init__.py | 2 + SimPEG/visulize/vtk/vtkTools.py | 285 ++++++++++++++++++++++++++++++++ SimPEG/visulize/vtk/vtkView.py | 95 +++++++++++ 4 files changed, 384 insertions(+) create mode 100644 SimPEG/visulize/__init__.py create mode 100644 SimPEG/visulize/vtk/__init__.py create mode 100644 SimPEG/visulize/vtk/vtkTools.py create mode 100644 SimPEG/visulize/vtk/vtkView.py diff --git a/SimPEG/visulize/__init__.py b/SimPEG/visulize/__init__.py new file mode 100644 index 00000000..485f9782 --- /dev/null +++ b/SimPEG/visulize/__init__.py @@ -0,0 +1,2 @@ +import vtk +#import mpl diff --git a/SimPEG/visulize/vtk/__init__.py b/SimPEG/visulize/vtk/__init__.py new file mode 100644 index 00000000..6d60ed5c --- /dev/null +++ b/SimPEG/visulize/vtk/__init__.py @@ -0,0 +1,2 @@ +from vtkTools import vtkTools +from vtkView import vtkView \ No newline at end of file diff --git a/SimPEG/visulize/vtk/vtkTools.py b/SimPEG/visulize/vtk/vtkTools.py new file mode 100644 index 00000000..3840a1bc --- /dev/null +++ b/SimPEG/visulize/vtk/vtkTools.py @@ -0,0 +1,285 @@ +import numpy as np, vtk, vtk.util.numpy_support as npsup, pdb +from SimPEG.utils import mkvc + + +class vtkTools(object): + """ + Class that interacts with VTK visulization toolkit. + + """ + + def __init__(self): + """ Initializes the VTK vtkTools. + + """ + + pass + + @staticmethod + def makeCellVTKObject(mesh,model): + """ + Make and return a cell based VTK object for a simpeg 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 + + Output: + :rtype: vtkRecilinearGrid object + :return: vtkObj + """ + + # 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. + # Asaign 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) + + return vtkObj + + @staticmethod + def makeFaceVTKObject(mesh,model): + """ + Make and return a face based VTK object for a simpeg 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). + Property array must be order hstack(Fx,Fy,Fz) + + Output: + :rtype: vtkUnstructuredGrid object + :return: vtkObj + """ + + ## Convert simpeg mesh to VTK properties + # Convert mesh nodes to vtkPoints + vtkPts = vtk.vtkPoints() + vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) + + # Define the face "cells" + # Using VTK_QUAD cell for faces (see VTK file format) + nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') + def faceR(mat,length): + return mat.T.reshape((length,1)) + # First direction + nTFx = np.prod(mesh.nFx) + FxCellBlock = np.hstack([ 4*np.ones((nTFx,1),dtype='int64'),faceR(nodeMat[:,:-1,:-1],nTFx),faceR(nodeMat[:,1: ,:-1],nTFx),faceR(nodeMat[:,1: ,1: ],nTFx),faceR(nodeMat[:,:-1,1: ],nTFx)] ) + FyCellBlock = np.array([],dtype='int64') + FzCellBlock = np.array([],dtype='int64') + # Second direction + if mesh.dim >= 2: + nTFy = np.prod(mesh.nFy) + FyCellBlock = np.hstack([ 4*np.ones((nTFy,1),dtype='int64'),faceR(nodeMat[:-1,:,:-1],nTFy),faceR(nodeMat[1: ,:,:-1],nTFy),faceR(nodeMat[1: ,:,1: ],nTFy),faceR(nodeMat[:-1,:,1: ],nTFy)] ) + # Third direction + if mesh.dim == 3: + nTFz = np.prod(mesh.nFz) + FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] ) + # Cells -cell array + FCellArr = vtk.vtkCellArray() + FCellArr.SetNumberOfCells(np.sum(mesh.nF)) + FCellArr.SetCells(np.sum(mesh.nF)*5,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) + # Cell type + FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(np.sum(mesh.nF),dtype='uint8'),deep=1) + # Cell location + FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nF)*5,5,dtype='int64'),deep=1) + + ## Make the object + vtkObj = vtk.vtkUnstructuredGrid() + # Set the objects properties + vtkObj.SetPoints(vtkPts) + vtkObj.SetCells(FCellType,FCellLoc,FCellArr) + + # 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.Update() + return vtkObj + + @staticmethod + def makeEdgeVTKObject(mesh,model): + """ + Make and return a edge based VTK object for a simpeg 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). + Property array must be order hstack(Ex,Ey,Ez) + + Output: + :rtype: vtkUnstructuredGrid object + :return: vtkObj + """ + + ## Convert simpeg mesh to VTK properties + # Convert mesh nodes to vtkPoints + vtkPts = vtk.vtkPoints() + vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) + + # Define the face "cells" + # Using VTK_QUAD cell for faces (see VTK file format) + nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') + def edgeR(mat,length): + return mat.T.reshape((length,1)) + # First direction + nTEx = np.prod(mesh.nEx) + ExCellBlock = np.hstack([ 2*np.ones((nTEx,1),dtype='int64'),edgeR(nodeMat[:-1,:,:],nTEx),edgeR(nodeMat[1:,:,:],nTEx)]) + # Second direction + if mesh.dim >= 2: + nTEy = np.prod(mesh.nEy) + EyCellBlock = np.hstack([ 2*np.ones((nTEy,1),dtype='int64'),edgeR(nodeMat[:,:-1,:],nTEy),edgeR(nodeMat[:,1:,:],nTEy)]) + # Third direction + if mesh.dim == 3: + nTEz = np.prod(mesh.nEz) + EzCellBlock = np.hstack([ 2*np.ones((nTEz,1),dtype='int64'),edgeR(nodeMat[:,:,:-1],nTEz),edgeR(nodeMat[:,:,1:],nTEz)]) + # Cells -cell array + ECellArr = vtk.vtkCellArray() + ECellArr.SetNumberOfCells(np.sum(mesh.nE)) + ECellArr.SetCells(np.sum(mesh.nE)*3,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1)) + # Cell type + ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(np.sum(mesh.nE),dtype='uint8'),deep=1) + # Cell location + ECellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nE)*3,3,dtype='int64'),deep=1) + + ## Make the object + vtkObj = vtk.vtkUnstructuredGrid() + # Set the objects properties + vtkObj.SetPoints(vtkPts) + vtkObj.SetCells(ECellType,ECellLoc,ECellArr) + + # 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) + + return vtkObj + + @staticmethod + def makeRenderWindow(ren): + renWin = vtk.vtkRenderWindow() + renWin.AddRenderer(ren) + iren = vtk.vtkRenderWindowInteractor() + iren.SetRenderWindow(renWin) + + return iren + + + @staticmethod + def closeRenderWindow(iren): + renwin = iren.GetRenderWindow() + renwin.Finalize() + iren.TerminateApp() + del iren, renwin + + @staticmethod + def makeVTKActor(vtkObj): + """ Makes a vtk mapper and Actor""" + mapper = vtk.vtkDataSetMapper() + mapper.SetInput(vtkObj) + actor = vtk.vtkActor() + actor.SetMapper(mapper) + actor.GetProperty().SetColor(0,0,0) + actor.GetProperty().SetRepresentationToWireframe() + + return actor + + @staticmethod + def startRenderWindow(iren): + """ Start a vtk rendering window""" + iren.Initialize() + renwin = iren.GetRenderWindow() + renwin.Render() + iren.Start() + + + # Simple write/read VTK xml model functions. + @staticmethod + def writeVTPFile(fileName,vtkPolyObject): + '''Function to write vtk polydata file (vtp).''' + polyWriter = vtk.vtkXMLPolyDataWriter() + polyWriter.SetInput(vtkPolyObject) + polyWriter.SetFileName(fileName) + polyWriter.Update() + + @staticmethod + def writeVTUFile(fileName,vtkUnstructuredGrid): + '''Function to write vtk unstructured grid (vtu).''' + Writer = vtk.vtkXMLUnstructuredGridWriter() + Writer.SetInput(vtkUnstructuredGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def writeVTRFile(fileName,vtkRectilinearGrid): + '''Function to write vtk rectilinear grid (vtr).''' + Writer = vtk.vtkXMLRectilinearGridWriter() + Writer.SetInput(vtkRectilinearGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def writeVTSFile(fileName,vtkStructuredGrid): + '''Function to write vtk structured grid (vts).''' + Writer = vtk.vtkXMLStructuredGridWriter() + Writer.SetInput(vtkStructuredGrid) + Writer.SetFileName(fileName) + Writer.Update() + + @staticmethod + def readVTSFile(fileName): + '''Function to read vtk structured grid (vts) and return a grid object.''' + Reader = vtk.vtkXMLStructuredGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTUFile(fileName): + '''Function to read vtk structured grid (vtu) and return a grid object.''' + Reader = vtk.vtkXMLUnstructuredGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTRFile(fileName): + '''Function to read vtk structured grid (vtr) and return a grid object.''' + Reader = vtk.vtkXMLRectilinearGridReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + + @staticmethod + def readVTPFile(fileName): + '''Function to read vtk structured grid (vtp) and return a grid object.''' + Reader = vtk.vtkXMLPolyDataReader() + Reader.SetFileName(fileName) + Reader.Update() + return Reader.GetOutput() + diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py new file mode 100644 index 00000000..f11fcfd8 --- /dev/null +++ b/SimPEG/visulize/vtk/vtkView.py @@ -0,0 +1,95 @@ +import numpy as np, vtk +import SimPEG as simpeg +#import SimPEG.visulize.vtk.vtkTools as vtkSP # Always get an error for this import + +class vtkView(object): + """ + Class for storing and view of SimPEG models in VTK (visulization toolkit). + + Inputs: + :param mesh, SimPEG mesh. + :param propdict, dictionary of property models. + Can have these dictionary names: + 'cell' - cell model; 'face' - face model; 'edge' - edge model + The dictionary properties are given as dictionaries with: + {'NameOfThePropertyModel': np.array of the properties}. + The property array has to be ordered in compliance with SimPEG standards. + + :: + Example of usages. + + ToDo + + """ + + def __init__(self,mesh,propdict): + """ + """ + + self.name = 'VTK figure of SimPEG model' + self._mesh = mesh + self._cell = None + self._faces = None + self._edges = None + + + + self._readPropertyDictionary(propdict) + + def _readPropertyDictionary(self,propdict): + """ + Reads the property and assigns to the object + """ + import SimPEG.visulize.vtk.vtkTools as vtkSP + + # Test the property dictionary + if len(propdict) > 3: + raise(Exception,'Too many input items in the property dictionary') + for propitem in propdict.iteritems(): + if propitem[0] in ['cell','face','edge']: + if propitem[0] == 'cell': + self._cell = vtkSP.makeCellVTKObject(self._mesh,propitem[1]) + if propitem[0] == 'face': + self._face = vtkSP.makeFaceVTKObject(self._mesh,propitem[1]) + if propitem[0] == 'edge': + self._edge = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1]) + else: + raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'cell\',\'face\',\'edge\'.'.format(propitem[0])) + + def Show(self,imageType='cell'): + """ + Open the VTK figure window + """ + #vtkSP = simpeg.visulize.vtk.vtkTools + import SimPEG.visulize.vtk.vtkTools as vtkSP + + # Make a renderer + ren = vtk.vtkRenderer() + # Make renderwindow. Returns the interactor. + iren = vtkSP.makeRenderWindow(ren) + + # Sort out the actor + if imageType == 'cell': + actor = vtkSP.makeVTKActor(self._cell) + elif imageType == 'face': + actor = vtkSP.makeVTKActor(self._face) + elif imageType == 'edge': + actor = vtkSP.makeVTKActor(self._edge) + actor.GetProperty().SetRepresentationToSurface() + else: + raise(Exception,"{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) + + ren.AddActor(actor) + ren.SetBackground(.5,.5,.5) + + vtkSP.startRenderWindow(iren) + + vtkSP.closeRenderWindow(iren) + + + + + + + + From e5a3e8de7d1575123fe97803586fc751418c8b44 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Mon, 4 Nov 2013 10:35:42 -0800 Subject: [PATCH 04/20] Added more help for vtkView. --- SimPEG/visulize/vtk/vtkView.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py index f11fcfd8..20bcdaa8 100644 --- a/SimPEG/visulize/vtk/vtkView.py +++ b/SimPEG/visulize/vtk/vtkView.py @@ -58,7 +58,11 @@ class vtkView(object): def Show(self,imageType='cell'): """ - Open the VTK figure window + Open the VTK figure window and show the mesh. + + Inputs: + param: str imageType: type of image {'cell','face','edge'} + """ #vtkSP = simpeg.visulize.vtk.vtkTools import SimPEG.visulize.vtk.vtkTools as vtkSP From 47a4716d956c83efd086353fd34586ff2cc818b4 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Sun, 10 Nov 2013 23:52:39 -0800 Subject: [PATCH 05/20] vtkView works for 'cell' data, but currently not for 'face' and 'edge'. --- SimPEG/visulize/vtk/vtkTools.py | 90 ++++++++++++++++++++++++++++++--- SimPEG/visulize/vtk/vtkView.py | 60 +++++++++++++++++----- 2 files changed, 130 insertions(+), 20 deletions(-) diff --git a/SimPEG/visulize/vtk/vtkTools.py b/SimPEG/visulize/vtk/vtkTools.py index 3840a1bc..50387a14 100644 --- a/SimPEG/visulize/vtk/vtkTools.py +++ b/SimPEG/visulize/vtk/vtkTools.py @@ -55,7 +55,8 @@ class vtkTools(object): vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) vtkDoubleArr.SetName(item[0]) vtkObj.GetCellData().AddArray(vtkDoubleArr) - + + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) return vtkObj @staticmethod @@ -117,6 +118,8 @@ class vtkTools(object): 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() return vtkObj @@ -177,17 +180,18 @@ class vtkTools(object): vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) vtkDoubleArr.SetName(item[0]) vtkObj.GetCellData().AddArray(vtkDoubleArr) - + + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) return vtkObj @staticmethod def makeRenderWindow(ren): - renWin = vtk.vtkRenderWindow() - renWin.AddRenderer(ren) + renwin = vtk.vtkRenderWindow() + renwin.AddRenderer(ren) iren = vtk.vtkRenderWindowInteractor() - iren.SetRenderWindow(renWin) + iren.SetRenderWindow(renwin) - return iren + return iren, renwin @staticmethod @@ -195,6 +199,7 @@ class vtkTools(object): renwin = iren.GetRenderWindow() renwin.Finalize() iren.TerminateApp() + del iren, renwin @staticmethod @@ -206,9 +211,80 @@ class vtkTools(object): actor.SetMapper(mapper) actor.GetProperty().SetColor(0,0,0) actor.GetProperty().SetRepresentationToWireframe() - return actor + @staticmethod + def makeVTKLODActor(vtkObj,clipper): + """Make LOD vtk Actor""" + selectMapper = vtk.vtkDataSetMapper() + selectMapper.SetInputConnection(clipper.GetOutputPort()) + selectMapper.SetScalarVisibility(1) + selectMapper.SetColorModeToMapScalars() + selectMapper.SetScalarModeToUseCellData() + selectMapper.SetScalarRange(clipper.GetInputDataObject(0,0).GetCellData().GetArray(0).GetRange()) + + selectActor = vtk.vtkLODActor() + selectActor.SetMapper(selectMapper) + selectActor.GetProperty().SetEdgeColor(1,0.5,0) + selectActor.GetProperty().SetEdgeVisibility(0) + selectActor.VisibilityOn() + selectActor.SetScale(1.01, 1.01, 1.01) + return selectActor + + @staticmethod + def setScalar2View(vtkObj,scalarName): + """ Sets the sclar to view """ + useArr = vtkObj.GetCellData().GetArray(scalarName) + if useArr == None: + raise IOError('Nerty array {:s} in the vtkObject'.format(scalarName)) + vtkObj.GetCellData().SetActiveScalars(scalarName) + + @staticmethod + def makeRectiVTKVOIThres(vtkObj): + """Make volume of interest and threshold for rectilinear grid.""" + cellCore = vtk.vtkExtractRectilinearGrid() + cellCore.SetInput(vtkObj) + cellCore.SetVOI(vtkObj.GetExtent()) + cellThres = vtk.vtkThreshold() + cellThres.AllScalarsOn() + cellThres.SetInputConnection(cellCore.GetOutputPort()) + cellThres.ThresholdByUpper(-1) + cellThres.Update() + return cellThres.GetOutput(), cellCore.GetOutput() + + @staticmethod + def makePlaneClipper(vtkObj): + """Makes a plane and clipper """ + plane = vtk.vtkPlane() + clipper = vtk.vtkClipDataSet() + clipper.SetInputConnection(vtkObj.GetProducerPort()) + clipper.SetClipFunction(plane) + clipper.InsideOutOff() + return clipper, plane + + @staticmethod + def makePlaneWidget(vtkObj,iren,plane,actor): + """Make an interactive planeWidget""" + + # Callback function + def movePlane(obj, events): + obj.GetPlane(intPlane) + intActor.VisibilityOn() + + # Associate the line widget with the interactor + planeWidget = vtk.vtkImplicitPlaneWidget() + planeWidget.SetInteractor(iren) + planeWidget.SetPlaceFactor(1.25) + planeWidget.SetInput(vtkObj) + planeWidget.PlaceWidget() + #planeWidget.AddObserver("InteractionEvent", movePlane) + planeWidget.SetScaleEnabled(0) + planeWidget.SetEnabled(1) + planeWidget.SetOutlineTranslation(0) + planeWidget.GetPlaneProperty().SetOpacity(0.1) + return planeWidget + + @staticmethod def startRenderWindow(iren): """ Start a vtk rendering window""" diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py index 20bcdaa8..bb1c4f60 100644 --- a/SimPEG/visulize/vtk/vtkView.py +++ b/SimPEG/visulize/vtk/vtkView.py @@ -28,14 +28,25 @@ class vtkView(object): self.name = 'VTK figure of SimPEG model' self._mesh = mesh + # Set vtk object containers self._cell = None self._faces = None self._edges = None - - self._readPropertyDictionary(propdict) + # Setup hidden properties + self._ren = None + self._iren = None + self._renwin = None + self._core = None + self._viewobj = None + self._plane = None + self._clipper = None + self._widget = None + self._actor = None + self._lut = None + def _readPropertyDictionary(self,propdict): """ Reads the property and assigns to the object @@ -68,27 +79,50 @@ class vtkView(object): import SimPEG.visulize.vtk.vtkTools as vtkSP # Make a renderer - ren = vtk.vtkRenderer() + self._ren = vtk.vtkRenderer() # Make renderwindow. Returns the interactor. - iren = vtkSP.makeRenderWindow(ren) + self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren) + # Sort out the actor if imageType == 'cell': - actor = vtkSP.makeVTKActor(self._cell) + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell) elif imageType == 'face': - actor = vtkSP.makeVTKActor(self._face) + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._face) elif imageType == 'edge': - actor = vtkSP.makeVTKActor(self._edge) - actor.GetProperty().SetRepresentationToSurface() + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._edge) else: - raise(Exception,"{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) + raise Exception("{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) - ren.AddActor(actor) - ren.SetBackground(.5,.5,.5) - vtkSP.startRenderWindow(iren) + global intPlane, intActor + self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj) + intActor = vtkSP.makeVTKLODActor(self._vtkobj,self._clipper) + self._widget = vtkSP.makePlaneWidget(self._vtkobj,self._iren,self._clipper.GetClipFunction(),self._actor) + # Callback function + self._plane = intPlane + self._actor = intActor + def movePlane(obj, events): + global intPlane, intActor + obj.GetPlane(intPlane) + intActor.VisibilityOn() - vtkSP.closeRenderWindow(iren) + self._widget.AddObserver("InteractionEvent",movePlane) + lut = vtk.vtkLookupTable() + lut.SetNumberOfColors(256) + lut.SetHueRange(0,0.66667) + lut.Build() + self._lut = lut + self._actor.GetMapper().SetLookupTable(lut) + + # Set renderer options + self._ren.SetBackground(.5,.5,.5) + self._ren.AddActor(self._actor) + + # Start the render Window + vtkSP.startRenderWindow(self._iren) + # Close the window when exited + vtkSP.closeRenderWindow(self._iren) From 708e7a6caf7206e36eea76880f83ca27c3074e3f Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Mon, 11 Nov 2013 00:01:16 -0800 Subject: [PATCH 06/20] Working example of 3D view with VTK --- notebooks/3DRenderingWithvtkTools.ipynb | 59 +++++++++++++++++++++++++ 1 file changed, 59 insertions(+) create mode 100644 notebooks/3DRenderingWithvtkTools.ipynb diff --git a/notebooks/3DRenderingWithvtkTools.ipynb b/notebooks/3DRenderingWithvtkTools.ipynb new file mode 100644 index 00000000..7209b10e --- /dev/null +++ b/notebooks/3DRenderingWithvtkTools.ipynb @@ -0,0 +1,59 @@ +{ + "metadata": { + "name": "3D rendering with vtkTools" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "code", + "collapsed": false, + "input": "import numpy as np, vtk\nimport SimPEG as simpeg", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 1 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*5\nh2 = np.ones(10)*10\nh3 = np.ones(5)*20\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "# Make a models that correspond to the cells, faces and edges.\nmodels = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n# Make the vtk viewer object.\nvtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n ", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 3 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "# Show the image \nvtkViewer.Show(imageType='cell')\n", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": "", + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 4 + } + ], + "metadata": {} + } + ] +} From 0e03e51bd281052961c6c525cb89584903aa6989 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Mon, 11 Nov 2013 08:59:28 -0800 Subject: [PATCH 07/20] Updated the vtkView class. --- SimPEG/visulize/vtk/vtkView.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py index bb1c4f60..334561ba 100644 --- a/SimPEG/visulize/vtk/vtkView.py +++ b/SimPEG/visulize/vtk/vtkView.py @@ -27,7 +27,10 @@ class vtkView(object): """ self.name = 'VTK figure of SimPEG model' + self.extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] + self.limits = [0, 10000] self._mesh = mesh + # Set vtk object containers self._cell = None self._faces = None @@ -123,6 +126,7 @@ class vtkView(object): vtkSP.startRenderWindow(self._iren) # Close the window when exited vtkSP.closeRenderWindow(self._iren) + del self._iren, self._renwin From 2cbbb3639c5adada3b6cf707ea706c597f0e0a94 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Mon, 11 Nov 2013 09:47:35 -0800 Subject: [PATCH 08/20] Fixed SimPEG.__init__.py to include visulization (some merging issuse). --- SimPEG/__init__.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index f50e1cef..f4628b6e 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -2,10 +2,6 @@ import utils from utils import Solver import mesh import inverse -<<<<<<< HEAD -from Solver import Solver import visulize -======= import forward import regularization ->>>>>>> refs/remotes/origin/master From d7c8e723792a3a961992fca585cff926ebe4fe6c Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Thu, 14 Nov 2013 10:50:33 -0800 Subject: [PATCH 09/20] Added examples directory. --- SimPEG/__init__.py | 1 + SimPEG/examples/__init__.py | 0 2 files changed, 1 insertion(+) create mode 100644 SimPEG/examples/__init__.py diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 7f059a74..a5103cd2 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -4,3 +4,4 @@ import mesh import inverse import forward import regularization +import examples diff --git a/SimPEG/examples/__init__.py b/SimPEG/examples/__init__.py new file mode 100644 index 00000000..e69de29b From 72e844c87d2f48bf6f9789acfbdffa2eac9d2d09 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Thu, 14 Nov 2013 11:53:14 -0800 Subject: [PATCH 10/20] Added a sub-module to utils to put analytics and source functions. --- SimPEG/utils/Geophysics/__init__.py | 1 + .../MagneticDipoleVectorPotential.py | 39 +++++++++++++++++++ SimPEG/utils/Geophysics/emSources/__init__.py | 1 + SimPEG/utils/__init__.py | 1 + 4 files changed, 42 insertions(+) create mode 100644 SimPEG/utils/Geophysics/__init__.py create mode 100644 SimPEG/utils/Geophysics/emSources/MagneticDipoleVectorPotential.py create mode 100644 SimPEG/utils/Geophysics/emSources/__init__.py diff --git a/SimPEG/utils/Geophysics/__init__.py b/SimPEG/utils/Geophysics/__init__.py new file mode 100644 index 00000000..bcf01943 --- /dev/null +++ b/SimPEG/utils/Geophysics/__init__.py @@ -0,0 +1 @@ +import emSources \ No newline at end of file diff --git a/SimPEG/utils/Geophysics/emSources/MagneticDipoleVectorPotential.py b/SimPEG/utils/Geophysics/emSources/MagneticDipoleVectorPotential.py new file mode 100644 index 00000000..16c5a1fb --- /dev/null +++ b/SimPEG/utils/Geophysics/emSources/MagneticDipoleVectorPotential.py @@ -0,0 +1,39 @@ +import numpy as np +from scipy.constants import mu_0, pi + +def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)): + """ + Calculate the vector potential of a set of magnetic dipoles + at given locations 'ref. ' + + :param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z) + :param numpy.ndarray obsLoc: Where the potentials will be calculated (x, y, z) + :param str component: The component to calculate - 'x', 'y', or 'z' + :param numpy.ndarray dipoleMoment: The vector dipole moment + :rtype: numpy.ndarray + :return: The vector potential each dipole at each observation location + """ + + if component=='x': + dimInd = 0 + elif component=='y': + dimInd = 1 + elif component=='z': + dimInd = 2 + else: + raise ValueError('Invalid component') + + txLoc = np.atleast_2d(txLoc) + obsLoc = np.atleast_2d(obsLoc) + + nEdges = obsLoc.shape[0] + nTx = txLoc.shape[0] + + m = np.array(dipoleMoment).repeat(nEdges, axis=0) + A = np.empty((nEdges, nTx)) + for i in range(nTx): + dR = obsLoc - txLoc[i, np.newaxis].repeat(nEdges, axis=0) + mCr = np.cross(m, dR) + r = np.sqrt((dR**2).sum(axis=1)) + A[:, i] = -(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3) + return A \ No newline at end of file diff --git a/SimPEG/utils/Geophysics/emSources/__init__.py b/SimPEG/utils/Geophysics/emSources/__init__.py new file mode 100644 index 00000000..b4e675ef --- /dev/null +++ b/SimPEG/utils/Geophysics/emSources/__init__.py @@ -0,0 +1 @@ +from MagneticDipoleVectorPotential import MagneticDipoleVectorPotential \ No newline at end of file diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index d2170721..843e1ff0 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -10,6 +10,7 @@ from interputils import interpmat from ipythonUtils import easyAnimate as animate import Solver from Solver import Solver +import Geophysics def setKwargs(obj, **kwargs): """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist.""" From cdf1dcf79fda4e81d20191cac72be4a4685de975 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Thu, 14 Nov 2013 15:48:32 -0800 Subject: [PATCH 11/20] Reorganized. --- SimPEG/utils/Geophysics/emSources/__init__.py | 2 +- .../{MagneticDipoleVectorPotential.py => emSources.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename SimPEG/utils/Geophysics/emSources/{MagneticDipoleVectorPotential.py => emSources.py} (100%) diff --git a/SimPEG/utils/Geophysics/emSources/__init__.py b/SimPEG/utils/Geophysics/emSources/__init__.py index b4e675ef..c05fda69 100644 --- a/SimPEG/utils/Geophysics/emSources/__init__.py +++ b/SimPEG/utils/Geophysics/emSources/__init__.py @@ -1 +1 @@ -from MagneticDipoleVectorPotential import MagneticDipoleVectorPotential \ No newline at end of file +from emSources import MagneticDipoleVectorPotential \ No newline at end of file diff --git a/SimPEG/utils/Geophysics/emSources/MagneticDipoleVectorPotential.py b/SimPEG/utils/Geophysics/emSources/emSources.py similarity index 100% rename from SimPEG/utils/Geophysics/emSources/MagneticDipoleVectorPotential.py rename to SimPEG/utils/Geophysics/emSources/emSources.py From 0cafa561048556a1d315af4e62ac840a7b4a506d Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Fri, 15 Nov 2013 08:50:29 -0800 Subject: [PATCH 12/20] Fixed code so that vtkView works for cell, edge and face. extent and limits work for the rendering as well. --- SimPEG/visulize/vtk/vtkTools.py | 25 ++++++++++++++++++++++--- SimPEG/visulize/vtk/vtkView.py | 8 +++++--- 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/SimPEG/visulize/vtk/vtkTools.py b/SimPEG/visulize/vtk/vtkTools.py index 50387a14..bd6642d4 100644 --- a/SimPEG/visulize/vtk/vtkTools.py +++ b/SimPEG/visulize/vtk/vtkTools.py @@ -240,15 +240,34 @@ class vtkTools(object): vtkObj.GetCellData().SetActiveScalars(scalarName) @staticmethod - def makeRectiVTKVOIThres(vtkObj): + def makeRectiVTKVOIThres(vtkObj,VOI,limits): """Make volume of interest and threshold for rectilinear grid.""" + # Check for the input cellCore = vtk.vtkExtractRectilinearGrid() + cellCore.SetVOI(VOI) cellCore.SetInput(vtkObj) - cellCore.SetVOI(vtkObj.GetExtent()) + cellThres = vtk.vtkThreshold() cellThres.AllScalarsOn() cellThres.SetInputConnection(cellCore.GetOutputPort()) - cellThres.ThresholdByUpper(-1) + cellThres.ThresholdByUpper(limits[0]) + cellThres.ThresholdByLower(limits[1]) + cellThres.Update() + return cellThres.GetOutput(), cellCore.GetOutput() + + @staticmethod + def makeUnstructVTKVOIThres(vtkObj,extent,limits): + """Make volume of interest and threshold for rectilinear grid.""" + # Check for the input + cellCore = vtk.vtkExtractUnstructuredGrid() + cellCore.SetExtent(extent) + cellCore.SetInput(vtkObj) + + cellThres = vtk.vtkThreshold() + cellThres.AllScalarsOn() + cellThres.SetInputConnection(cellCore.GetOutputPort()) + cellThres.ThresholdByUpper(limits[0]) + cellThres.ThresholdByLower(limits[1]) cellThres.Update() return cellThres.GetOutput(), cellCore.GetOutput() diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py index 334561ba..f331e1d7 100644 --- a/SimPEG/visulize/vtk/vtkView.py +++ b/SimPEG/visulize/vtk/vtkView.py @@ -89,11 +89,13 @@ class vtkView(object): # Sort out the actor if imageType == 'cell': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell) + self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell,self.extent,self.limits) elif imageType == 'face': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._face) + extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ] + self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._face,extent,self.limits) elif imageType == 'edge': - self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._edge) + extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ] + self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._edge,extent,self.limits) else: raise Exception("{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) From c6be54712123b5933a25a8e31283eebe0de6a117 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Fri, 15 Nov 2013 09:42:06 -0800 Subject: [PATCH 13/20] Add vtkView.viewProp to indicate which model to view at any given time. No check wether the name or number exists as an array in the vtkView._vtkobj. --- SimPEG/visulize/vtk/vtkView.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py index f331e1d7..74389e05 100644 --- a/SimPEG/visulize/vtk/vtkView.py +++ b/SimPEG/visulize/vtk/vtkView.py @@ -26,11 +26,14 @@ class vtkView(object): """ """ + # ToDo: Set the properties up so that there are set/get methods self.name = 'VTK figure of SimPEG model' self.extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] - self.limits = [0, 10000] + self.limits = [0, 1e12] + self.viewProp = 0 # Int or name of the vector. self._mesh = mesh + # Set vtk object containers self._cell = None self._faces = None @@ -99,7 +102,15 @@ class vtkView(object): else: raise Exception("{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) - + # Set the active scalar. + if type(self.viewProp) == int: + actScalar = self._vtkobj.GetCellData().GetArrayName(self.viewProp) + elif type(self.viewProp) == str: + actScalar = self.viewProp + else : + raise Exception('The vtkView.viewProp has the wrong format. Has to be interger or a string.') + self._vtkobj.GetCellData().SetActiveScalars(actScalar) + # Set up the plane, clipper and the user interaction. global intPlane, intActor self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj) intActor = vtkSP.makeVTKLODActor(self._vtkobj,self._clipper) From 2a770f6fd418fbc1d1339e745d688f985952cf91 Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Fri, 15 Nov 2013 16:18:50 -0800 Subject: [PATCH 14/20] Updated nF and nE --- SimPEG/visulize/vtk/vtkTools.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/SimPEG/visulize/vtk/vtkTools.py b/SimPEG/visulize/vtk/vtkTools.py index bd6642d4..314d15e4 100644 --- a/SimPEG/visulize/vtk/vtkTools.py +++ b/SimPEG/visulize/vtk/vtkTools.py @@ -99,12 +99,12 @@ class vtkTools(object): FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] ) # Cells -cell array FCellArr = vtk.vtkCellArray() - FCellArr.SetNumberOfCells(np.sum(mesh.nF)) - FCellArr.SetCells(np.sum(mesh.nF)*5,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) + FCellArr.SetNumberOfCells(mesh.nF) + FCellArr.SetCells(mesh.nF*5,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1)) # Cell type - FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(np.sum(mesh.nF),dtype='uint8'),deep=1) + FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(mesh.nF,dtype='uint8'),deep=1) # Cell location - FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nF)*5,5,dtype='int64'),deep=1) + FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,mesh.nF*5,5,dtype='int64'),deep=1) ## Make the object vtkObj = vtk.vtkUnstructuredGrid() @@ -161,12 +161,12 @@ class vtkTools(object): EzCellBlock = np.hstack([ 2*np.ones((nTEz,1),dtype='int64'),edgeR(nodeMat[:,:,:-1],nTEz),edgeR(nodeMat[:,:,1:],nTEz)]) # Cells -cell array ECellArr = vtk.vtkCellArray() - ECellArr.SetNumberOfCells(np.sum(mesh.nE)) - ECellArr.SetCells(np.sum(mesh.nE)*3,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1)) + ECellArr.SetNumberOfCells(mesh.nE) + ECellArr.SetCells(mesh.nE*3,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1)) # Cell type - ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(np.sum(mesh.nE),dtype='uint8'),deep=1) + ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(mesh.nE,dtype='uint8'),deep=1) # Cell location - ECellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,np.sum(mesh.nE)*3,3,dtype='int64'),deep=1) + ECellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,mesh.nE*3,3,dtype='int64'),deep=1) ## Make the object vtkObj = vtk.vtkUnstructuredGrid() From 37bcf891be571bf7013be548204a77bd8701608a Mon Sep 17 00:00:00 2001 From: Gudni Karl Rosenkjaer Date: Sun, 17 Nov 2013 14:31:44 -0800 Subject: [PATCH 15/20] Fixed interface so the renderer uses Trackball interaction, extent and limits work. Added a notebook, showing a simple example of using the vtkView --- SimPEG/visulize/vtk/vtkTools.py | 1 + SimPEG/visulize/vtk/vtkView.py | 16 ++++++++-------- notebooks/3DRenderingWithvtkTools.ipynb | 20 +++++++++++++++++--- 3 files changed, 26 insertions(+), 11 deletions(-) diff --git a/SimPEG/visulize/vtk/vtkTools.py b/SimPEG/visulize/vtk/vtkTools.py index 314d15e4..f08cd345 100644 --- a/SimPEG/visulize/vtk/vtkTools.py +++ b/SimPEG/visulize/vtk/vtkTools.py @@ -189,6 +189,7 @@ class vtkTools(object): renwin = vtk.vtkRenderWindow() renwin.AddRenderer(ren) iren = vtk.vtkRenderWindowInteractor() + iren.GetInteractorStyle().SetCurrentStyleToTrackballCamera() iren.SetRenderWindow(renwin) return iren, renwin diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py index 74389e05..11e83ed4 100644 --- a/SimPEG/visulize/vtk/vtkView.py +++ b/SimPEG/visulize/vtk/vtkView.py @@ -30,7 +30,7 @@ class vtkView(object): self.name = 'VTK figure of SimPEG model' self.extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1] self.limits = [0, 1e12] - self.viewProp = 0 # Int or name of the vector. + self.viewprop = {'cell':0} # Name of the tyep and Int order of the array or name of the vector. self._mesh = mesh @@ -73,7 +73,7 @@ class vtkView(object): else: raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'cell\',\'face\',\'edge\'.'.format(propitem[0])) - def Show(self,imageType='cell'): + def Show(self): """ Open the VTK figure window and show the mesh. @@ -89,7 +89,7 @@ class vtkView(object): # Make renderwindow. Returns the interactor. self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren) - + imageType = self.viewprop.keys()[0] # Sort out the actor if imageType == 'cell': self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell,self.extent,self.limits) @@ -103,12 +103,12 @@ class vtkView(object): raise Exception("{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType)) # Set the active scalar. - if type(self.viewProp) == int: - actScalar = self._vtkobj.GetCellData().GetArrayName(self.viewProp) - elif type(self.viewProp) == str: - actScalar = self.viewProp + if type(self.viewprop.values()[0]) == int: + actScalar = self._vtkobj.GetCellData().GetArrayName(self.viewprop.values()[0]) + elif type(self.viewprop.values()[0]) == str: + actScalar = self.viewprop.values()[0] else : - raise Exception('The vtkView.viewProp has the wrong format. Has to be interger or a string.') + raise Exception('The vtkView.viewprop.values()[0] has the wrong format. Has to be interger or a string.') self._vtkobj.GetCellData().SetActiveScalars(actScalar) # Set up the plane, clipper and the user interaction. global intPlane, intActor diff --git a/notebooks/3DRenderingWithvtkTools.ipynb b/notebooks/3DRenderingWithvtkTools.ipynb index 7209b10e..c6f1663d 100644 --- a/notebooks/3DRenderingWithvtkTools.ipynb +++ b/notebooks/3DRenderingWithvtkTools.ipynb @@ -19,7 +19,7 @@ { "cell_type": "code", "collapsed": false, - "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*5\nh2 = np.ones(10)*10\nh3 = np.ones(5)*20\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", + "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*50\nh2 = np.ones(10)*100\nh3 = np.ones(5)*200\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", "language": "python", "metadata": {}, "outputs": [], @@ -43,14 +43,28 @@ "outputs": [], "prompt_number": 4 }, + { + "cell_type": "code", + "collapsed": false, + "input": "mesh.nCxx", + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "pyout", + "prompt_number": 7, + "text": "20" + } + ], + "prompt_number": 7 + }, { "cell_type": "code", "collapsed": false, "input": "", "language": "python", "metadata": {}, - "outputs": [], - "prompt_number": 4 + "outputs": [] } ], "metadata": {} From a950ba8dd1068f7107be342c11c210f8e8ef88da Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Mon, 18 Nov 2013 10:54:54 -0800 Subject: [PATCH 16/20] Bug fix in emSources. Beginning to change Cyl1DMesh to be more comparable with TensorMesh/BaseMesh. --- SimPEG/mesh/Cyl1DMesh.py | 22 +++++++++---------- .../utils/Geophysics/emSources/emSources.py | 1 + 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/SimPEG/mesh/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py index fa2bb729..28e2e1a9 100644 --- a/SimPEG/mesh/Cyl1DMesh.py +++ b/SimPEG/mesh/Cyl1DMesh.py @@ -62,11 +62,11 @@ class Cyl1DMesh(object): # Counting #################################################### - def nCr(): + def nCx(): doc = "Number of cells in the radial direction" fget = lambda self: self.hr.size return locals() - nCr = property(**nCr()) + nCx = property(**nCx()) def nCz(): doc = "Number of cells in the z direction" @@ -76,7 +76,7 @@ class Cyl1DMesh(object): def nC(): doc = "Total number of cells" - fget = lambda self: self.nCr * self.nCz + fget = lambda self: self.nCx * self.nCz return locals() nC = property(**nC()) @@ -106,7 +106,7 @@ class Cyl1DMesh(object): def nFz(): doc = "Number of z faces" - fget = lambda self: self.nNz * self.nCr + fget = lambda self: self.nNz * self.nCx return locals() nFz = property(**nFz()) @@ -242,12 +242,12 @@ class Cyl1DMesh(object): def fget(self): if self._edgeCurl is None: #1D Difference matricies - dr = sp.spdiags((np.ones((self.nCr+1, 1))*[-1, 1]).T, [-1,0], self.nCr, self.nCr, format="csr") + dr = sp.spdiags((np.ones((self.nCx+1, 1))*[-1, 1]).T, [-1,0], self.nCx, self.nCx, format="csr") dz = sp.spdiags((np.ones((self.nCz+1, 1))*[-1, 1]).T, [0,1], self.nCz, self.nCz+1, format="csr") #2D Difference matricies Dr = sp.kron(sp.eye(self.nNz), dr) - Dz = -sp.kron(dz, sp.eye(self.nCr)) #Not sure about this negative + Dz = -sp.kron(dz, sp.eye(self.nCx)) #Not sure about this negative #Edge curl operator self._edgeCurl = sp.diags(1/self.area,0)*sp.vstack((Dz, Dr))*sp.diags(self.edge,0) @@ -261,7 +261,7 @@ class Cyl1DMesh(object): def fget(self): if self._aveE2CC is None: az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr') - ar = sp.spdiags(0.5*np.ones((2, self.nCr)), [0, 1], self.nCr, self.nCr, format='csr') + ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr') ar[0,0] = 1 self._aveE2CC = sp.kron(az, ar).T return self._aveE2CC @@ -274,10 +274,10 @@ class Cyl1DMesh(object): def fget(self): if self._aveF2CC is None: az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr') - ar = sp.spdiags(0.5*np.ones((2, self.nCr)), [0, 1], self.nCr, self.nCr, format='csr') + ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr') ar[0,0] = 1 Afr = sp.kron(sp.eye(self.nCz),ar) - Afz = sp.kron(az,sp.eye(self.nCr)) + Afz = sp.kron(az,sp.eye(self.nCx)) self._aveF2CC = sp.vstack((Afr,Afz)).T return self._aveF2CC return locals() @@ -311,7 +311,7 @@ class Cyl1DMesh(object): elif type(materialProp) is float: materialProp = np.ones(self.nC)*materialProp elif materialProp.shape == (self.nCz,): - materialProp = materialProp.repeat(self.nCr) + materialProp = materialProp.repeat(self.nCx) materialProp = mkvc(materialProp) assert materialProp.shape == (self.nC,), "materialProp incorrect shape" @@ -383,7 +383,7 @@ class Cyl1DMesh(object): dFz = np.sum(dFz**2, axis=1) indBL = np.argmin(dFz) # Face below and to the left - indAL = indBL + self.nCr # Face above and to the left + indAL = indBL + self.nCx # Face above and to the left zF_BL = self.gridFz[indBL,:] zF_AL = self.gridFz[indAL,:] diff --git a/SimPEG/utils/Geophysics/emSources/emSources.py b/SimPEG/utils/Geophysics/emSources/emSources.py index 16c5a1fb..db492a94 100644 --- a/SimPEG/utils/Geophysics/emSources/emSources.py +++ b/SimPEG/utils/Geophysics/emSources/emSources.py @@ -25,6 +25,7 @@ def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0. txLoc = np.atleast_2d(txLoc) obsLoc = np.atleast_2d(obsLoc) + dipoleMoment = np.atleast_2d(dipoleMoment) nEdges = obsLoc.shape[0] nTx = txLoc.shape[0] From 8992f195a4008d96bf8d5a5a57559c4e7f6a4224 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Mon, 18 Nov 2013 11:06:35 -0800 Subject: [PATCH 17/20] Added nCv property to BaseMesh and Cyl1DMesh classes. --- SimPEG/mesh/BaseMesh.py | 11 +++++++++++ SimPEG/mesh/Cyl1DMesh.py | 6 ++++++ 2 files changed, 17 insertions(+) diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py index 2647cf44..ff804a83 100644 --- a/SimPEG/mesh/BaseMesh.py +++ b/SimPEG/mesh/BaseMesh.py @@ -228,6 +228,17 @@ class BaseMesh(object): return locals() nC = property(**nC()) + def nCv(): + doc = """ + Total number of cells in each direction + + :rtype: numpy.array (dim, ) + :return: [nCx, nCy, nCz] + """ + fget = lambda self: np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None]) + return locals() + nCv = property(**nCv()) + def nNx(): doc = """ Number of nodes in the x-direction diff --git a/SimPEG/mesh/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py index 28e2e1a9..e22e12b9 100644 --- a/SimPEG/mesh/Cyl1DMesh.py +++ b/SimPEG/mesh/Cyl1DMesh.py @@ -80,6 +80,12 @@ class Cyl1DMesh(object): return locals() nC = property(**nC()) + def nCv(): + doc = "Total number of cells in each direction" + fget = lambda self: np.array([self.nCx, self.nCz]) + return locals() + nCv = property(**nCv()) + def nNr(): doc = "Number of nodes in the radial direction" fget = lambda self: self.hr.size From 7d295dc57e9ad5a47ffd53f66c915e3677ee8c4c Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Mon, 18 Nov 2013 11:24:27 -0800 Subject: [PATCH 18/20] Fix in example notebook. --- SimPEG/visulize/vtk/vtkView.py | 6 +- notebooks/3DRenderingWithvtkTools.ipynb | 95 ++++++++++++------------- 2 files changed, 47 insertions(+), 54 deletions(-) diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visulize/vtk/vtkView.py index 11e83ed4..2ce984df 100644 --- a/SimPEG/visulize/vtk/vtkView.py +++ b/SimPEG/visulize/vtk/vtkView.py @@ -75,11 +75,7 @@ class vtkView(object): def Show(self): """ - Open the VTK figure window and show the mesh. - - Inputs: - param: str imageType: type of image {'cell','face','edge'} - + Open the VTK figure window and show the mesh. """ #vtkSP = simpeg.visulize.vtk.vtkTools import SimPEG.visulize.vtk.vtkTools as vtkSP diff --git a/notebooks/3DRenderingWithvtkTools.ipynb b/notebooks/3DRenderingWithvtkTools.ipynb index c6f1663d..b0daa408 100644 --- a/notebooks/3DRenderingWithvtkTools.ipynb +++ b/notebooks/3DRenderingWithvtkTools.ipynb @@ -1,6 +1,6 @@ { "metadata": { - "name": "3D rendering with vtkTools" + "name": "" }, "nbformat": 3, "nbformat_minor": 0, @@ -10,64 +10,61 @@ { "cell_type": "code", "collapsed": false, - "input": "import numpy as np, vtk\nimport SimPEG as simpeg", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 1 - }, - { - "cell_type": "code", - "collapsed": false, - "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*50\nh2 = np.ones(10)*100\nh3 = np.ones(5)*200\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 2 - }, - { - "cell_type": "code", - "collapsed": false, - "input": "# Make a models that correspond to the cells, faces and edges.\nmodels = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n# Make the vtk viewer object.\nvtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n ", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 3 - }, - { - "cell_type": "code", - "collapsed": false, - "input": "# Show the image \nvtkViewer.Show(imageType='cell')\n", - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 - }, - { - "cell_type": "code", - "collapsed": false, - "input": "mesh.nCxx", - "language": "python", - "metadata": {}, - "outputs": [ - { - "output_type": "pyout", - "prompt_number": 7, - "text": "20" - } + "input": [ + "import numpy as np, vtk\n", + "import SimPEG as simpeg" ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 5 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "#Make a mesh and model\n", + "x0 = np.zeros(3)\n", + "h1 = np.ones(20)*50\n", + "h2 = np.ones(10)*100\n", + "h3 = np.ones(5)*200\n", + "\n", + "mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 6 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Make a models that correspond to the cells, faces and edges.\n", + "models = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n", + "# Make the vtk viewer object.\n", + "vtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n", + " " + ], + "language": "python", + "metadata": {}, + "outputs": [], "prompt_number": 7 }, { "cell_type": "code", "collapsed": false, - "input": "", + "input": [ + "# Show the image \n", + "vtkViewer.Show()\n" + ], "language": "python", "metadata": {}, - "outputs": [] + "outputs": [], + "prompt_number": "*" } ], "metadata": {} } ] -} +} \ No newline at end of file From a2c305e6bc926952e8bbc940fbbc883f5e3a14a0 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 18 Nov 2013 11:56:35 -0800 Subject: [PATCH 19/20] Renamed module. --- SimPEG/__init__.py | 2 +- SimPEG/{visulize => visualize}/__init__.py | 0 .../{visulize => visualize}/vtk/__init__.py | 0 .../{visulize => visualize}/vtk/vtkTools.py | 0 SimPEG/{visulize => visualize}/vtk/vtkView.py | 20 +++++++++---------- notebooks/3DRenderingWithvtkTools.ipynb | 4 ++-- 6 files changed, 13 insertions(+), 13 deletions(-) rename SimPEG/{visulize => visualize}/__init__.py (100%) rename SimPEG/{visulize => visualize}/vtk/__init__.py (100%) rename SimPEG/{visulize => visualize}/vtk/vtkTools.py (100%) rename SimPEG/{visulize => visualize}/vtk/vtkView.py (90%) diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 41f7c431..df7b8b33 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -2,7 +2,7 @@ import utils from utils import Solver import mesh import inverse -import visulize +import visualize import forward import regularization import examples diff --git a/SimPEG/visulize/__init__.py b/SimPEG/visualize/__init__.py similarity index 100% rename from SimPEG/visulize/__init__.py rename to SimPEG/visualize/__init__.py diff --git a/SimPEG/visulize/vtk/__init__.py b/SimPEG/visualize/vtk/__init__.py similarity index 100% rename from SimPEG/visulize/vtk/__init__.py rename to SimPEG/visualize/vtk/__init__.py diff --git a/SimPEG/visulize/vtk/vtkTools.py b/SimPEG/visualize/vtk/vtkTools.py similarity index 100% rename from SimPEG/visulize/vtk/vtkTools.py rename to SimPEG/visualize/vtk/vtkTools.py diff --git a/SimPEG/visulize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py similarity index 90% rename from SimPEG/visulize/vtk/vtkView.py rename to SimPEG/visualize/vtk/vtkView.py index 2ce984df..6bf6dab5 100644 --- a/SimPEG/visulize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -1,6 +1,6 @@ import numpy as np, vtk import SimPEG as simpeg -#import SimPEG.visulize.vtk.vtkTools as vtkSP # Always get an error for this import +#import SimPEG.visualize.vtk.vtkTools as vtkSP # Always get an error for this import class vtkView(object): """ @@ -8,11 +8,11 @@ class vtkView(object): Inputs: :param mesh, SimPEG mesh. - :param propdict, dictionary of property models. + :param propdict, dictionary of property models. Can have these dictionary names: 'cell' - cell model; 'face' - face model; 'edge' - edge model The dictionary properties are given as dictionaries with: - {'NameOfThePropertyModel': np.array of the properties}. + {'NameOfThePropertyModel': np.array of the properties}. The property array has to be ordered in compliance with SimPEG standards. :: @@ -54,10 +54,10 @@ class vtkView(object): self._lut = None def _readPropertyDictionary(self,propdict): - """ + """ Reads the property and assigns to the object """ - import SimPEG.visulize.vtk.vtkTools as vtkSP + import SimPEG.visualize.vtk.vtkTools as vtkSP # Test the property dictionary if len(propdict) > 3: @@ -75,10 +75,10 @@ class vtkView(object): def Show(self): """ - Open the VTK figure window and show the mesh. + Open the VTK figure window and show the mesh. """ - #vtkSP = simpeg.visulize.vtk.vtkTools - import SimPEG.visulize.vtk.vtkTools as vtkSP + #vtkSP = simpeg.visualize.vtk.vtkTools + import SimPEG.visualize.vtk.vtkTools as vtkSP # Make a renderer self._ren = vtk.vtkRenderer() @@ -119,7 +119,7 @@ class vtkView(object): obj.GetPlane(intPlane) intActor.VisibilityOn() - self._widget.AddObserver("InteractionEvent",movePlane) + self._widget.AddObserver("InteractionEvent",movePlane) lut = vtk.vtkLookupTable() lut.SetNumberOfColors(256) lut.SetHueRange(0,0.66667) @@ -131,7 +131,7 @@ class vtkView(object): self._ren.SetBackground(.5,.5,.5) self._ren.AddActor(self._actor) - # Start the render Window + # Start the render Window vtkSP.startRenderWindow(self._iren) # Close the window when exited vtkSP.closeRenderWindow(self._iren) diff --git a/notebooks/3DRenderingWithvtkTools.ipynb b/notebooks/3DRenderingWithvtkTools.ipynb index b0daa408..1ad88317 100644 --- a/notebooks/3DRenderingWithvtkTools.ipynb +++ b/notebooks/3DRenderingWithvtkTools.ipynb @@ -43,7 +43,7 @@ "# Make a models that correspond to the cells, faces and edges.\n", "models = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n", "# Make the vtk viewer object.\n", - "vtkViewer = simpeg.visulize.vtk.vtkView(mesh,models) \n", + "vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models) \n", " " ], "language": "python", @@ -67,4 +67,4 @@ "metadata": {} } ] -} \ No newline at end of file +} From 1f3799caa0d1ed76a6737e2e173a00d962db77da Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 18 Nov 2013 12:19:56 -0800 Subject: [PATCH 20/20] Added try/catch on the vtk import for the documentation. Added small example at bottom of vtkView. --- SimPEG/visualize/vtk/vtkTools.py | 58 +++++++++++++++++--------------- SimPEG/visualize/vtk/vtkView.py | 24 ++++++++++--- 2 files changed, 50 insertions(+), 32 deletions(-) diff --git a/SimPEG/visualize/vtk/vtkTools.py b/SimPEG/visualize/vtk/vtkTools.py index f08cd345..9853d089 100644 --- a/SimPEG/visualize/vtk/vtkTools.py +++ b/SimPEG/visualize/vtk/vtkTools.py @@ -1,16 +1,20 @@ -import numpy as np, vtk, vtk.util.numpy_support as npsup, pdb +import numpy as np +try: + import vtk, vtk.util.numpy_support as npsup, pdb +except Exception, e: + print 'VTK import error. Please ensure you have VTK installed to use this visualization package.' from SimPEG.utils import mkvc class vtkTools(object): - """ + """ Class that interacts with VTK visulization toolkit. """ def __init__(self): """ Initializes the VTK vtkTools. - + """ pass @@ -20,7 +24,7 @@ class vtkTools(object): """ Make and return a cell based VTK object for a simpeg mesh and model. - Input: + 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 @@ -55,7 +59,7 @@ class vtkTools(object): vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) vtkDoubleArr.SetName(item[0]) vtkObj.GetCellData().AddArray(vtkDoubleArr) - + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) return vtkObj @@ -64,9 +68,9 @@ class vtkTools(object): """ Make and return a face based VTK object for a simpeg mesh and model. - Input: + Input: :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK - :param model, dictionary of numpy.array - Name('s) and array('s). + :param model, dictionary of numpy.array - Name('s) and array('s). Property array must be order hstack(Fx,Fy,Fz) Output: @@ -78,7 +82,7 @@ class vtkTools(object): # Convert mesh nodes to vtkPoints vtkPts = vtk.vtkPoints() vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) - + # Define the face "cells" # Using VTK_QUAD cell for faces (see VTK file format) nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') @@ -96,7 +100,7 @@ class vtkTools(object): # Third direction if mesh.dim == 3: nTFz = np.prod(mesh.nFz) - FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] ) + FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] ) # Cells -cell array FCellArr = vtk.vtkCellArray() FCellArr.SetNumberOfCells(mesh.nF) @@ -105,12 +109,12 @@ class vtkTools(object): FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(mesh.nF,dtype='uint8'),deep=1) # Cell location FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,mesh.nF*5,5,dtype='int64'),deep=1) - + ## Make the object vtkObj = vtk.vtkUnstructuredGrid() # Set the objects properties vtkObj.SetPoints(vtkPts) - vtkObj.SetCells(FCellType,FCellLoc,FCellArr) + vtkObj.SetCells(FCellType,FCellLoc,FCellArr) # Assign the model('s) to the object for item in model.iteritems(): @@ -118,7 +122,7 @@ class vtkTools(object): 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() return vtkObj @@ -128,9 +132,9 @@ class vtkTools(object): """ Make and return a edge based VTK object for a simpeg mesh and model. - Input: + Input: :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK - :param model, dictionary of numpy.array - Name('s) and array('s). + :param model, dictionary of numpy.array - Name('s) and array('s). Property array must be order hstack(Ex,Ey,Ez) Output: @@ -142,7 +146,7 @@ class vtkTools(object): # Convert mesh nodes to vtkPoints vtkPts = vtk.vtkPoints() vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1)) - + # Define the face "cells" # Using VTK_QUAD cell for faces (see VTK file format) nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M') @@ -158,7 +162,7 @@ class vtkTools(object): # Third direction if mesh.dim == 3: nTEz = np.prod(mesh.nEz) - EzCellBlock = np.hstack([ 2*np.ones((nTEz,1),dtype='int64'),edgeR(nodeMat[:,:,:-1],nTEz),edgeR(nodeMat[:,:,1:],nTEz)]) + EzCellBlock = np.hstack([ 2*np.ones((nTEz,1),dtype='int64'),edgeR(nodeMat[:,:,:-1],nTEz),edgeR(nodeMat[:,:,1:],nTEz)]) # Cells -cell array ECellArr = vtk.vtkCellArray() ECellArr.SetNumberOfCells(mesh.nE) @@ -172,7 +176,7 @@ class vtkTools(object): vtkObj = vtk.vtkUnstructuredGrid() # Set the objects properties vtkObj.SetPoints(vtkPts) - vtkObj.SetCells(ECellType,ECellLoc,ECellArr) + vtkObj.SetCells(ECellType,ECellLoc,ECellArr) # Assign the model('s) to the object for item in model.iteritems(): @@ -180,7 +184,7 @@ class vtkTools(object): vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) vtkDoubleArr.SetName(item[0]) vtkObj.GetCellData().AddArray(vtkDoubleArr) - + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) return vtkObj @@ -200,7 +204,7 @@ class vtkTools(object): renwin = iren.GetRenderWindow() renwin.Finalize() iren.TerminateApp() - + del iren, renwin @staticmethod @@ -239,7 +243,7 @@ class vtkTools(object): if useArr == None: raise IOError('Nerty array {:s} in the vtkObject'.format(scalarName)) vtkObj.GetCellData().SetActiveScalars(scalarName) - + @staticmethod def makeRectiVTKVOIThres(vtkObj,VOI,limits): """Make volume of interest and threshold for rectilinear grid.""" @@ -247,15 +251,15 @@ class vtkTools(object): cellCore = vtk.vtkExtractRectilinearGrid() cellCore.SetVOI(VOI) cellCore.SetInput(vtkObj) - + cellThres = vtk.vtkThreshold() - cellThres.AllScalarsOn() - cellThres.SetInputConnection(cellCore.GetOutputPort()) + cellThres.AllScalarsOn() + cellThres.SetInputConnection(cellCore.GetOutputPort()) cellThres.ThresholdByUpper(limits[0]) cellThres.ThresholdByLower(limits[1]) cellThres.Update() return cellThres.GetOutput(), cellCore.GetOutput() - + @staticmethod def makeUnstructVTKVOIThres(vtkObj,extent,limits): """Make volume of interest and threshold for rectilinear grid.""" @@ -263,10 +267,10 @@ class vtkTools(object): cellCore = vtk.vtkExtractUnstructuredGrid() cellCore.SetExtent(extent) cellCore.SetInput(vtkObj) - + cellThres = vtk.vtkThreshold() - cellThres.AllScalarsOn() - cellThres.SetInputConnection(cellCore.GetOutputPort()) + cellThres.AllScalarsOn() + cellThres.SetInputConnection(cellCore.GetOutputPort()) cellThres.ThresholdByUpper(limits[0]) cellThres.ThresholdByLower(limits[1]) cellThres.Update() diff --git a/SimPEG/visualize/vtk/vtkView.py b/SimPEG/visualize/vtk/vtkView.py index 6bf6dab5..abedd0bf 100644 --- a/SimPEG/visualize/vtk/vtkView.py +++ b/SimPEG/visualize/vtk/vtkView.py @@ -1,6 +1,10 @@ -import numpy as np, vtk +import numpy as np +try: + import vtk + #import SimPEG.visualize.vtk.vtkTools as vtkSP # Always get an error for this import +except Exception, e: + print 'VTK import error. Please ensure you have VTK installed to use this visualization package.' import SimPEG as simpeg -#import SimPEG.visualize.vtk.vtkTools as vtkSP # Always get an error for this import class vtkView(object): """ @@ -139,8 +143,18 @@ class vtkView(object): +if __name__ == '__main__': + #Make a mesh and model + x0 = np.zeros(3) + h1 = np.ones(20)*50 + h2 = np.ones(10)*100 + h3 = np.ones(5)*200 + mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0) - - - + # Make a models that correspond to the cells, faces and edges. + models = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}} + # Make the vtk viewer object. + vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models) + # Show the image + vtkViewer.Show()