diff --git a/simpegPF/BaseMag.py b/simpegPF/BaseMag.py index 2aacc9c4..c7519038 100644 --- a/simpegPF/BaseMag.py +++ b/simpegPF/BaseMag.py @@ -105,6 +105,54 @@ class BaseMagSurvey(Survey.BaseSurvey): return np.r_[bfx, bfy, bfz] +class LinearSurvey(Survey.BaseSurvey): + """Base Magnetics Survey""" + + rxLoc = None #: receiver locations + rxType = None #: receiver type + + def __init__(self, srcField, **kwargs): + self.srcField = srcField + Survey.BaseSurvey.__init__(self, **kwargs) + + def eval(self, u): + return u + + @property + def nD(self): + return self.prob.G.shape[0] + + @property + def nRx(self): + return self.srcField.rxList[0].locs.shape[0] + # def setBackgroundField(self, SrcField): + + # if getattr(self, 'B0', None) is None: + # self._B0 = SrcField.param[0] * dipazm_2_xyz(SrcField.param[1],SrcField.param[2]) + + # return self._B0 + +class SrcField(Survey.BaseSrc): + """ Define the inducing field """ + + param = None #: Inducing field param (Amp, Incl, Decl) + + def __init__(self, rxList, A, I, D, **kwargs): + self.param = (A,I,D) + super(SrcField, self).__init__(rxList, **kwargs) + +class RxObs(Survey.BaseRx): + """A station location must have be located in 3-D""" + def __init__(self, locsXYZ, **kwargs): + locs = locsXYZ + assert locsXYZ.shape[1] == 3, 'locs must in 3-D (x,y,z).' + super(RxObs, self).__init__(locs, 'tmi', storeProjections=False, **kwargs) + + @property + def nD(self): + """Number of data in the receiver.""" + return self.locs[0].shape[0] + class MagSurveyBx(object): """docstring for MagSurveyBx""" def __init__(self, **kwargs): @@ -144,352 +192,3 @@ class WeightMap(Maps.IdentityMap): return Utils.sdiag(self.weight) - - -def readUBCmagObs(obs_file): - - """ - Read and write UBC mag file format - - INPUT: - :param fileName, path to the UBC obs mag file - - OUTPUT: - :param dobs, observation in (x y z [data] [wd]) - :param B, primary field information (BI, BD, B0) - :param M, magnetization orentiaton (MI, MD) - - """ - - fid = open(obs_file,'r') - - # First line has the inclination,declination and amplitude of B0 - line = fid.readline() - B = np.array(line.split(),dtype=float) - - # Second line has the magnetization orientation and a flag - line = fid.readline() - M = np.array(line.split(),dtype=float) - - # Third line has the number of rows - line = fid.readline() - ndat = np.array(line.split(),dtype=int) - - # Pre-allocate space for obsx, obsy, obsz, data, uncert - line = fid.readline() - temp = np.array(line.split(),dtype=float) - - dobs = np.zeros((ndat,len(temp)), dtype=float) - - - for ii in range(ndat): - - dobs[ii,:] = np.array(line.split(),dtype=float) - line = fid.readline() - - return B, M, dobs - -def read_MAGfwr_inp(input_file): - - """Read input files for forward modeling MAG data with integral form - INPUT: - input_file: File name containing the forward parameter - - OUTPUT: - mshfile - obsfile - modfile - magfile - topofile - # All files should be in the working directory, otherwise the path must - # be specified. - - Created on Jul 17, 2013 - - @author: dominiquef - """ - - - fid = open(input_file,'r') - - line = fid.readline() - l_input = line.split('!') - mshfile = l_input[0].rstrip() - - line = fid.readline() - l_input = line.split('!') - obsfile = l_input[0].rstrip() - - line = fid.readline() - l_input = line.split('!') - modfile = l_input[0].rstrip() - - line = fid.readline() - l_input = line.split('!') - if l_input=='null': - magfile = [] - - else: - magfile = l_input[0].rstrip() - - - line = fid.readline() - l_input = line.split('!') - if l_input=='null': - topofile = [] - - else: - topofile = l_input[0].rstrip() - - return mshfile, obsfile, modfile, magfile, topofile - -def read_MAGinv_inp(input_file): - """Read input files for forward modeling MAG data with integral form - INPUT: - input_file: File name containing the forward parameter - - OUTPUT: - mshfile - obsfile - topofile - start model - ref model - mag model - weightfile - chi_target - as, ax ,ay, az - upper, lower bounds - lp, lqx, lqy, lqz - - # All files should be in the working directory, otherwise the path must - # be specified. - - Created on Dec 21th, 2015 - - @author: dominiquef - """ - - - fid = open(input_file,'r') - - # Line 1 - line = fid.readline() - l_input = line.split('!') - mshfile = l_input[0].rstrip() - - # Line 2 - line = fid.readline() - l_input = line.split('!') - obsfile = l_input[0].rstrip() - - # Line 3 - line = fid.readline() - l_input = re.split('[!\s]',line) - if l_input=='null': - topofile = [] - - else: - topofile = l_input[0].rstrip() - - - # Line 4 - line = fid.readline() - l_input = re.split('[!\s]',line) - if l_input[0]=='VALUE': - mstart = float(l_input[1]) - - else: - mstart = l_input[0].rstrip() - - # Line 5 - line = fid.readline() - l_input = re.split('[!\s]',line) - if l_input[0]=='VALUE': - mref = float(l_input[1]) - - else: - mref = l_input[0].rstrip() - - - # Line 6 - line = fid.readline() - l_input = re.split('[!\s]',line) - if l_input=='DEFAULT': - magfile = [] - - else: - magfile = l_input[0].rstrip() - - # Line 7 - line = fid.readline() - l_input = re.split('[!\s]',line) - if l_input=='DEFAULT': - wgtfile = [] - - else: - wgtfile = l_input[0].rstrip() - - # Line 8 - line = fid.readline() - l_input = re.split('[!\s]',line) - chi = float(l_input[0]) - - # Line 9 - line = fid.readline() - l_input = re.split('[!\s]',line) - val = np.array(l_input[0:4]) - alphas = val.astype(np.float) - - # Line 10 - line = fid.readline() - l_input = re.split('[!\s]',line) - if l_input[0]=='VALUE': - val = np.array(l_input[1:3]) - bounds = val.astype(np.float) - - else: - bounds = l_input[0].rstrip() - - # Line 11 - line = fid.readline() - l_input = re.split('[!\s]',line) - if l_input[0]=='VALUE': - val = np.array(l_input[1:6]) - lpnorms = val.astype(np.float) - - else: - lpnorms = l_input[0].rstrip() - - return mshfile, obsfile, topofile, mstart, mref, magfile, wgtfile, chi, alphas, bounds, lpnorms - -def read_GOCAD_ts(tsfile): - """Read GOCAD triangulated surface (*.ts) file - INPUT: - tsfile: Triangulated surface - - OUTPUT: - vrts : Array of vertices in XYZ coordinates [n x 3] - trgl : Array of index for triangles [m x 3]. The order of the vertices - is important and describes the normal - n = cross( (P2 - P1 ) , (P3 - P1) ) - - - Created on Jan 13th, 2016 - - Author: @fourndo - """ - - - fid = open(tsfile,'r') - line = fid.readline() - - # Skip all the lines until the vertices - while re.match('TFACE',line)==None: - line = fid.readline() - - line = fid.readline() - vrtx = [] - - # Run down all the vertices and save in array - while re.match('VRTX',line): - l_input = re.split('[\s*]',line) - temp = np.array(l_input[2:5]) - vrtx.append(temp.astype(np.float)) - - # Read next line - line = fid.readline() - - vrtx = np.asarray(vrtx) - - # Skip lines to the triangles - while re.match('TRGL',line)==None: - line = fid.readline() - - # Run down the list of triangles - trgl = [] - - # Run down all the vertices and save in array - while re.match('TRGL',line): - l_input = re.split('[\s*]',line) - temp = np.array(l_input[1:4]) - trgl.append(temp.astype(np.int)) - - # Read next line - line = fid.readline() - - trgl = np.asarray(trgl) - - return vrtx, trgl - -def gocad2vtk(gcFile,mesh,bcflag,inflag): - """" - Function to read gocad polystructure file and output indexes of mesh with in the structure. - - """ - import vtk, vtk.util.numpy_support as npsup - - print "Reading GOCAD ts file..." - vrtx, trgl = read_GOCAD_ts(gcFile) - # Adjust the index - trgl = trgl - 1 - - # Make vtk pts - ptsvtk = vtk.vtkPoints() - ptsvtk.SetData(npsup.numpy_to_vtk(vrtx,deep=1)) - - # Make the polygon connection - polys = vtk.vtkCellArray() - for face in trgl: - poly = vtk.vtkPolygon() - poly.GetPointIds().SetNumberOfIds(len(face)) - for nrv, vert in enumerate(face): - poly.GetPointIds().SetId(nrv,vert) - polys.InsertNextCell(poly) - - # Make the polydata, structure of connections and vrtx - polyData = vtk.vtkPolyData() - polyData.SetPoints(ptsvtk) - polyData.SetPolys(polys) - - # Make implicit func - ImpDistFunc = vtk.vtkImplicitPolyDataDistance() - ImpDistFunc.SetInput(polyData) - - # Convert the mesh - vtkMesh = vtk.vtkRectilinearGrid() - vtkMesh.SetDimensions(mesh.nNx,mesh.nNy,mesh.nNz) - vtkMesh.SetXCoordinates(npsup.numpy_to_vtk(mesh.vectorNx,deep=1)) - vtkMesh.SetYCoordinates(npsup.numpy_to_vtk(mesh.vectorNy,deep=1)) - vtkMesh.SetZCoordinates(npsup.numpy_to_vtk(mesh.vectorNz,deep=1)) - # Add indexes - vtkInd = npsup.numpy_to_vtk(np.arange(mesh.nC),deep=1) - vtkInd.SetName('Index') - vtkMesh.GetCellData().AddArray(vtkInd) - - extractImpDistRectGridFilt = vtk.vtkExtractGeometry() # Object constructor - extractImpDistRectGridFilt.SetImplicitFunction(ImpDistFunc) # - extractImpDistRectGridFilt.SetInputData(vtkMesh) - - if bcflag is True: - extractImpDistRectGridFilt.ExtractBoundaryCellsOn() - - else: - extractImpDistRectGridFilt.ExtractBoundaryCellsOff() - - if inflag is True: - extractImpDistRectGridFilt.ExtractInsideOn() - - else: - extractImpDistRectGridFilt.ExtractInsideOff() - - print "Extracting indices from grid..." - # Executing the pipe - extractImpDistRectGridFilt.Update() - - # Get index inside - insideGrid = extractImpDistRectGridFilt.GetOutput() - insideGrid = npsup.vtk_to_numpy(insideGrid.GetCellData().GetArray('Index')) - - - # Return the indexes inside - return insideGrid \ No newline at end of file diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py index efe91ce1..4caa384b 100644 --- a/simpegPF/Magnetics.py +++ b/simpegPF/Magnetics.py @@ -1,18 +1,43 @@ from SimPEG import * -import BaseMag +import BaseMag as MAG from scipy.constants import mu_0 from MagAnalytics import spheremodel, CongruousMagBC +import re class MagneticIntegral(Problem.BaseProblem): - surveyPair = Survey.LinearSurvey + #surveyPair = Survey.LinearSurvey - def __init__(self, mesh, G, mapping=None, **kwargs): + storeG = True #: Store the forward matrix by default, otherwise just compute d + actInd = None #: Active cell indices provided + M = None #: Magnetization matrix provided, otherwise all induced + def __init__(self, mesh, mapping=None, **kwargs): Problem.BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs) - self.G = G + + + def fwr_ind(self): + # Add forward function + # kappa = self.curModel.kappa TODO + kappa = self.mapping*self.curModel + return self.G.dot(kappa) + + def fwr_rem(self): + #TODO check if we are inverting for M + return self.G.dot(self.mapping(m)) + def fields(self, m): - - return self.G.dot(self.mapping*(m)) + self.curModel = m + total = np.zeros(self.survey.nRx) + induced = self.fwr_ind() + # rem = self.rem + + if induced is not None: + total += induced + + return total + + + # return self.G.dot(self.mapping*(m)) def Jvec(self, m, v, u=None): dmudm = self.mapping.deriv(m) @@ -22,15 +47,171 @@ class MagneticIntegral(Problem.BaseProblem): dmudm = self.mapping.deriv(m) return dmudm.T * (self.G.T.dot(v)) + @property + def G(self): + if not self.ispaired: + raise Exception('Need to pair!') + if getattr(self,'_G', None) is None: + self._G = self.Intrgl_Fwr_Op( 'ind' ) + + + return self._G + + # @property + # def Grem(self): + # if not self.ispaired: + # raise Exception('Need to pair!') + + # if getattr(self,'_Grem', None) is None: + # self._Grem = Intrgl_Fwr_Op('full') + + + # return self._Grem + + def Intrgl_Fwr_Op(self, flag): + + """ + + Magnetic forward operator in integral form + + flag = 'ind' | 'full' + + 1- ind : Magnetization fixed by user + + 3- full: Full tensor matrix stored with shape([3*ndata, 3*nc]) + + Return + _G = Linear forward modeling operation + + Created on March, 13th 2016 + + @author: dominiquef + + """ + # Find non-zero cells + #inds = np.nonzero(actv)[0] + if getattr(self, 'actInd', None) is not None: + if self.actInd.dtype=='bool': + inds = np.asarray([inds for inds, elem in enumerate(self.actInd, 1) if elem], dtype = int) - 1 + else: + inds = self.actInd + + + + else: + + inds = np.asarray(range(self.mesh.nC)) + + nC = len(inds) + + # Create active cell projector + P = sp.csr_matrix((np.ones(nC),(inds, range(nC))), + shape=(self.mesh.nC, nC)) + + + + + # Create vectors of nodal location (lower and upper coners for each cell) + xn = self.mesh.vectorNx; + yn = self.mesh.vectorNy; + zn = self.mesh.vectorNz; + + yn2,xn2,zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:]) + yn1,xn1,zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1]) + + Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)] + Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)] + Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)] + + rxLoc = self.survey.srcField.rxList[0].locs + ndata = rxLoc.shape[0] + + survey = self.survey + + # Pre-allocate space and create magnetization matrix if required + if (flag=='ind'): + + # # If assumes uniform magnetization direction + # if M.shape != (nC,3): + + # print 'Magnetization vector must be Nc x 3' + # return + if getattr(self, 'M', None) is None: + M = dipazm_2_xyz(np.ones(nC) * survey.srcField.param[1],np.ones(nC) * survey.srcField.param[2]) + + + Mx = Utils.sdiag(M[:,0]*survey.srcField.param[0]) + My = Utils.sdiag(M[:,1]*survey.srcField.param[0]) + Mz = Utils.sdiag(M[:,2]*survey.srcField.param[0]) + + Mxyz = sp.vstack((Mx,My,Mz)) + + + + if survey.srcField.rxList[0].rxType == 'tmi': + G = np.zeros((ndata, nC)) + + # Convert Bdecination from north to cartesian + D = (450.-float(survey.srcField.param[2]))%360. + I = survey.srcField.param[1] + # Projection matrix + Ptmi = mkvc(np.r_[np.cos(np.deg2rad(I))*np.cos(np.deg2rad(D)), + np.cos(np.deg2rad(I))*np.sin(np.deg2rad(D)), + np.sin(np.deg2rad(I))],2).T; + + elif survey.srcField.rxList[0].rxType == 'xyz': + + G = np.zeros((int(3*ndata), nC)) + + elif flag == 'full': + G = np.zeros((int(3*ndata), int(3*nC))) + + + else: + print """Flag must be either 'tmi' | 'xyz' | 'full', please revised""" + return + + + # Loop through all observations and create forward operator (ndata-by-nC) + print "Begin calculation of forward operator: " + flag + + # Add counter to dsiplay progress. Good for large problems + count = -1; + for ii in range(ndata): + + + tx, ty, tz = get_T_mat(Xn,Yn,Zn,rxLoc[ii,:]) + if flag == 'ind': + + if survey.srcField.rxList[0].rxType =='tmi': + G[ii,:] = Ptmi.dot(np.vstack((tx,ty,tz)))*Mxyz + + elif survey.srcField.rxList[0].rxType =='xyz': + G[ii,:] = tx*Mxyz + G[ii+ndata,:] = ty*Mxyz + G[ii+2*ndata,:] = tz*Mxyz + + elif flag == 'full': + G[ii,:] = tx + G[ii+ndata,:] = ty + G[ii+2*ndata,:] = tz + + + # Display progress + count = progress(ii,count,ndata) + + print "Done 100% ...forward operator completed!!\n" + + return G class MagneticsDiffSecondary(Problem.BaseProblem): """ Secondary field approach using differential equations! """ - surveyPair = BaseMag.BaseMagSurvey - modelPair = BaseMag.BaseMagMap + surveyPair = MAG.BaseMagSurvey + modelPair = MAG.BaseMagMap def __init__(self, model, mapping=None, **kwargs): Problem.BaseProblem.__init__(self, model, mapping=mapping, **kwargs) @@ -823,9 +1004,9 @@ def dipazm_2_xyz(dip,azm_N): M = np.zeros((nC,3)) # Modify azimuth from North to Cartesian-X - azm_X = (450.- azm_N) % 360. + azm_X = (450.- np.asarray(azm_N)) % 360. - D = np.deg2rad(dip) + D = np.deg2rad(np.asarray(dip)) I = np.deg2rad(azm_X) M[:,0] = np.cos(D) * np.cos(I) ; @@ -926,7 +1107,7 @@ def get_dist_wgt(mesh,rxLoc,actv,R,R0): return wr -def writeUBCobs(filename,B,M,rxLoc,d,wd): +def writeUBCobs(filename,survey,d): """ writeUBCobs(filename,B,M,rxLoc,d,wd) @@ -934,11 +1115,8 @@ def writeUBCobs(filename,B,M,rxLoc,d,wd): INPUT filename : Name of out file including directory - B : Inducing field parameters [Inc, Decl, Intensity] - M : Magnetization orientation [Inc, Decl, dtype] - rxLoc : Observation locations [obsx, obsy, obsz] - d : Data vector - wd : Uncertainty vector + survey + flag : dobs | dpred OUTPUT Obsfile @@ -948,11 +1126,17 @@ def writeUBCobs(filename,B,M,rxLoc,d,wd): @author: dominiquef """ + B = survey.srcField.param + + rxLoc = survey.srcField.rxList[0].locs + + wd = survey.std + data = np.c_[rxLoc , d , wd] with file(filename,'w') as fid: - fid.write('%6.2f %6.2f %6.2f\n' %(B[0], B[1], B[2]) ) - fid.write('%6.2f %6.2f %6.2f\n' %(M[0], M[1], 1) ) + fid.write('%6.2f %6.2f %6.2f\n' %(B[2], B[1], B[0]) ) + fid.write('%6.2f %6.2f %6.2f\n' %(B[2], B[1], 1) ) fid.write('%i\n' %len(d) ) np.savetxt(fid, data, fmt='%e',delimiter=' ',newline='\n') @@ -1048,4 +1232,359 @@ def plot_obs_2D(rxLoc,d,wd,varstr): plt.scatter(rxLoc[:,0],rxLoc[:,1], c=d, s=20) plt.title(varstr) plt.gca().set_aspect('equal', adjustable='box') - \ No newline at end of file + +def readUBCmagObs(obs_file): + + """ + Read and write UBC mag file format + + INPUT: + :param fileName, path to the UBC obs mag file + + OUTPUT: + :param survey + :param M, magnetization orentiaton (MI, MD) + + """ + + fid = open(obs_file,'r') + + # First line has the inclination,declination and amplitude of B0 + line = fid.readline() + B = np.array(line.split(),dtype=float) + + # Second line has the magnetization orientation and a flag + line = fid.readline() + M = np.array(line.split(),dtype=float) + + # Third line has the number of rows + line = fid.readline() + ndat = np.array(line.split(),dtype=int) + + # Pre-allocate space for obsx, obsy, obsz, data, uncert + line = fid.readline() + temp = np.array(line.split(),dtype=float) + + d = np.zeros(ndat, dtype=float) + wd = np.zeros(ndat, dtype=float) + locXYZ = np.zeros( (ndat,3), dtype=float) + + for ii in range(ndat): + + temp = np.array(line.split(),dtype=float) + locXYZ[ii,:] = temp[:3] + d[ii] = temp[3] + wd[ii] = temp[4] + line = fid.readline() + + rxLoc = MAG.RxObs(locXYZ) + srcField = MAG.SrcField([rxLoc],B[2],B[0],B[1]) + survey = MAG.LinearSurvey(srcField) + survey.dobs = d + survey.std = wd + return survey + +def read_MAGfwr_inp(input_file): + + """Read input files for forward modeling MAG data with integral form + INPUT: + input_file: File name containing the forward parameter + + OUTPUT: + mshfile + obsfile + modfile + magfile + topofile + # All files should be in the working directory, otherwise the path must + # be specified. + + Created on Jul 17, 2013 + + @author: dominiquef + """ + + + fid = open(input_file,'r') + + line = fid.readline() + l_input = line.split('!') + mshfile = l_input[0].rstrip() + + line = fid.readline() + l_input = line.split('!') + obsfile = l_input[0].rstrip() + + line = fid.readline() + l_input = line.split('!') + modfile = l_input[0].rstrip() + + line = fid.readline() + l_input = line.split('!') + if l_input=='null': + magfile = [] + + else: + magfile = l_input[0].rstrip() + + + line = fid.readline() + l_input = line.split('!') + if l_input=='null': + topofile = [] + + else: + topofile = l_input[0].rstrip() + + return mshfile, obsfile, modfile, magfile, topofile + +def read_MAGinv_inp(input_file): + """Read input files for forward modeling MAG data with integral form + INPUT: + input_file: File name containing the forward parameter + + OUTPUT: + mshfile + obsfile + topofile + start model + ref model + mag model + weightfile + chi_target + as, ax ,ay, az + upper, lower bounds + lp, lqx, lqy, lqz + + # All files should be in the working directory, otherwise the path must + # be specified. + + Created on Dec 21th, 2015 + + @author: dominiquef + """ + + + fid = open(input_file,'r') + + # Line 1 + line = fid.readline() + l_input = line.split('!') + mshfile = l_input[0].rstrip() + + # Line 2 + line = fid.readline() + l_input = line.split('!') + obsfile = l_input[0].rstrip() + + # Line 3 + line = fid.readline() + l_input = re.split('[!\s]',line) + if l_input=='null': + topofile = [] + + else: + topofile = l_input[0].rstrip() + + + # Line 4 + line = fid.readline() + l_input = re.split('[!\s]',line) + if l_input[0]=='VALUE': + mstart = float(l_input[1]) + + else: + mstart = l_input[0].rstrip() + + # Line 5 + line = fid.readline() + l_input = re.split('[!\s]',line) + if l_input[0]=='VALUE': + mref = float(l_input[1]) + + else: + mref = l_input[0].rstrip() + + + # Line 6 + line = fid.readline() + l_input = re.split('[!\s]',line) + if l_input=='DEFAULT': + magfile = [] + + else: + magfile = l_input[0].rstrip() + + # Line 7 + line = fid.readline() + l_input = re.split('[!\s]',line) + if l_input=='DEFAULT': + wgtfile = [] + + else: + wgtfile = l_input[0].rstrip() + + # Line 8 + line = fid.readline() + l_input = re.split('[!\s]',line) + chi = float(l_input[0]) + + # Line 9 + line = fid.readline() + l_input = re.split('[!\s]',line) + val = np.array(l_input[0:4]) + alphas = val.astype(np.float) + + # Line 10 + line = fid.readline() + l_input = re.split('[!\s]',line) + if l_input[0]=='VALUE': + val = np.array(l_input[1:3]) + bounds = val.astype(np.float) + + else: + bounds = l_input[0].rstrip() + + # Line 11 + line = fid.readline() + l_input = re.split('[!\s]',line) + if l_input[0]=='VALUE': + val = np.array(l_input[1:6]) + lpnorms = val.astype(np.float) + + else: + lpnorms = l_input[0].rstrip() + + return mshfile, obsfile, topofile, mstart, mref, magfile, wgtfile, chi, alphas, bounds, lpnorms + +def read_GOCAD_ts(tsfile): + """Read GOCAD triangulated surface (*.ts) file + INPUT: + tsfile: Triangulated surface + + OUTPUT: + vrts : Array of vertices in XYZ coordinates [n x 3] + trgl : Array of index for triangles [m x 3]. The order of the vertices + is important and describes the normal + n = cross( (P2 - P1 ) , (P3 - P1) ) + + + Created on Jan 13th, 2016 + + Author: @fourndo + """ + + + fid = open(tsfile,'r') + line = fid.readline() + + # Skip all the lines until the vertices + while re.match('TFACE',line)==None: + line = fid.readline() + + line = fid.readline() + vrtx = [] + + # Run down all the vertices and save in array + while re.match('VRTX',line): + l_input = re.split('[\s*]',line) + temp = np.array(l_input[2:5]) + vrtx.append(temp.astype(np.float)) + + # Read next line + line = fid.readline() + + vrtx = np.asarray(vrtx) + + # Skip lines to the triangles + while re.match('TRGL',line)==None: + line = fid.readline() + + # Run down the list of triangles + trgl = [] + + # Run down all the vertices and save in array + while re.match('TRGL',line): + l_input = re.split('[\s*]',line) + temp = np.array(l_input[1:4]) + trgl.append(temp.astype(np.int)) + + # Read next line + line = fid.readline() + + trgl = np.asarray(trgl) + + return vrtx, trgl + +def gocad2vtk(gcFile,mesh,bcflag,inflag): + """" + Function to read gocad polystructure file and output indexes of mesh with in the structure. + + """ + import vtk, vtk.util.numpy_support as npsup + + print "Reading GOCAD ts file..." + vrtx, trgl = read_GOCAD_ts(gcFile) + # Adjust the index + trgl = trgl - 1 + + # Make vtk pts + ptsvtk = vtk.vtkPoints() + ptsvtk.SetData(npsup.numpy_to_vtk(vrtx,deep=1)) + + # Make the polygon connection + polys = vtk.vtkCellArray() + for face in trgl: + poly = vtk.vtkPolygon() + poly.GetPointIds().SetNumberOfIds(len(face)) + for nrv, vert in enumerate(face): + poly.GetPointIds().SetId(nrv,vert) + polys.InsertNextCell(poly) + + # Make the polydata, structure of connections and vrtx + polyData = vtk.vtkPolyData() + polyData.SetPoints(ptsvtk) + polyData.SetPolys(polys) + + # Make implicit func + ImpDistFunc = vtk.vtkImplicitPolyDataDistance() + ImpDistFunc.SetInput(polyData) + + # Convert the mesh + vtkMesh = vtk.vtkRectilinearGrid() + vtkMesh.SetDimensions(mesh.nNx,mesh.nNy,mesh.nNz) + vtkMesh.SetXCoordinates(npsup.numpy_to_vtk(mesh.vectorNx,deep=1)) + vtkMesh.SetYCoordinates(npsup.numpy_to_vtk(mesh.vectorNy,deep=1)) + vtkMesh.SetZCoordinates(npsup.numpy_to_vtk(mesh.vectorNz,deep=1)) + # Add indexes + vtkInd = npsup.numpy_to_vtk(np.arange(mesh.nC),deep=1) + vtkInd.SetName('Index') + vtkMesh.GetCellData().AddArray(vtkInd) + + extractImpDistRectGridFilt = vtk.vtkExtractGeometry() # Object constructor + extractImpDistRectGridFilt.SetImplicitFunction(ImpDistFunc) # + extractImpDistRectGridFilt.SetInputData(vtkMesh) + + if bcflag is True: + extractImpDistRectGridFilt.ExtractBoundaryCellsOn() + + else: + extractImpDistRectGridFilt.ExtractBoundaryCellsOff() + + if inflag is True: + extractImpDistRectGridFilt.ExtractInsideOn() + + else: + extractImpDistRectGridFilt.ExtractInsideOff() + + print "Extracting indices from grid..." + # Executing the pipe + extractImpDistRectGridFilt.Update() + + # Get index inside + insideGrid = extractImpDistRectGridFilt.GetOutput() + insideGrid = npsup.vtk_to_numpy(insideGrid.GetCellData().GetArray('Index')) + + + # Return the indexes inside + return insideGrid \ No newline at end of file diff --git a/simpegPF/notebooks/SimPEG Tutorial - MAG Linear Problem.ipynb b/simpegPF/notebooks/SimPEG Tutorial - MAG Linear Problem.ipynb index 0abe8ef6..169b3cfc 100644 --- a/simpegPF/notebooks/SimPEG Tutorial - MAG Linear Problem.ipynb +++ b/simpegPF/notebooks/SimPEG Tutorial - MAG Linear Problem.ipynb @@ -43,7 +43,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 4, "metadata": { "collapsed": false }, @@ -53,19 +53,15 @@ "output_type": "stream", "text": [ "Using matplotlib backend: nbAgg\n", - "Populating the interactive namespace from numpy and matplotlib\n", - "Efficiency Warning: Interpolation will be slow, use setup.py!\n", - "\n", - " python setup.py build_ext --inplace\n", - " \n" + "Populating the interactive namespace from numpy and matplotlib\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ - "C:\\Users\\dominiquef.MIRAGEOSCIENCE\\AppData\\Local\\Continuum\\Anaconda\\lib\\site-packages\\IPython\\kernel\\__init__.py:13: ShimWarning: The `IPython.kernel` package has been deprecated. You should import from ipykernel or jupyter_client instead.\n", - " \"You should import from ipykernel or jupyter_client instead.\", ShimWarning)\n" + "WARNING: pylab import has clobbered these variables: ['linalg']\n", + "`%matplotlib` prevents importing * from pylab and numpy\n" ] } ], @@ -78,7 +74,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 5, "metadata": { "collapsed": false, "scrolled": true @@ -168,7 +164,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 6, "metadata": { "collapsed": false }, @@ -205,7 +201,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 7, "metadata": { "collapsed": false }, @@ -992,903 +988,6 @@ "plt.gca().set_aspect('equal', adjustable='box')" ] }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "application/javascript": [ - "/* Put everything inside the global mpl namespace */\n", - "window.mpl = {};\n", - "\n", - "mpl.get_websocket_type = function() {\n", - " if (typeof(WebSocket) !== 'undefined') {\n", - " return WebSocket;\n", - " } else if (typeof(MozWebSocket) !== 'undefined') {\n", - " return MozWebSocket;\n", - " } else {\n", - " alert('Your browser does not have WebSocket support.' +\n", - " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", - " 'Firefox 4 and 5 are also supported but you ' +\n", - " 'have to enable WebSockets in about:config.');\n", - " };\n", - "}\n", - "\n", - "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", - " this.id = figure_id;\n", - "\n", - " this.ws = websocket;\n", - "\n", - " this.supports_binary = (this.ws.binaryType != undefined);\n", - "\n", - " if (!this.supports_binary) {\n", - " var warnings = document.getElementById(\"mpl-warnings\");\n", - " if (warnings) {\n", - " warnings.style.display = 'block';\n", - " warnings.textContent = (\n", - " \"This browser does not support binary websocket messages. \" +\n", - " \"Performance may be slow.\");\n", - " }\n", - " }\n", - "\n", - " this.imageObj = new Image();\n", - "\n", - " this.context = undefined;\n", - " this.message = undefined;\n", - " this.canvas = undefined;\n", - " this.rubberband_canvas = undefined;\n", - " this.rubberband_context = undefined;\n", - " this.format_dropdown = undefined;\n", - "\n", - " this.image_mode = 'full';\n", - "\n", - " this.root = $('
');\n", - " this._root_extra_style(this.root)\n", - " this.root.attr('style', 'display: inline-block');\n", - "\n", - " $(parent_element).append(this.root);\n", - "\n", - " this._init_header(this);\n", - " this._init_canvas(this);\n", - " this._init_toolbar(this);\n", - "\n", - " var fig = this;\n", - "\n", - " this.waiting = false;\n", - "\n", - " this.ws.onopen = function () {\n", - " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", - " fig.send_message(\"send_image_mode\", {});\n", - " fig.send_message(\"refresh\", {});\n", - " }\n", - "\n", - " this.imageObj.onload = function() {\n", - " if (fig.image_mode == 'full') {\n", - " // Full images could contain transparency (where diff images\n", - " // almost always do), so we need to clear the canvas so that\n", - " // there is no ghosting.\n", - " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", - " }\n", - " fig.context.drawImage(fig.imageObj, 0, 0);\n", - " fig.waiting = false;\n", - " };\n", - "\n", - " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", - " }\n", - "\n", - " this.ws.onmessage = this._make_on_message_function(this);\n", - "\n", - " this.ondownload = ondownload;\n", - "}\n", - "\n", - "mpl.figure.prototype._init_header = function() {\n", - " var titlebar = $(\n", - " '
');\n", - " var titletext = $(\n", - " '
');\n", - " titlebar.append(titletext)\n", - " this.root.append(titlebar);\n", - " this.header = titletext[0];\n", - "}\n", - "\n", - "\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._init_canvas = function() {\n", - " var fig = this;\n", - "\n", - " var canvas_div = $('
');\n", - "\n", - " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", - "\n", - " function canvas_keyboard_event(event) {\n", - " return fig.key_event(event, event['data']);\n", - " }\n", - "\n", - " canvas_div.keydown('key_press', canvas_keyboard_event);\n", - " canvas_div.keyup('key_release', canvas_keyboard_event);\n", - " this.canvas_div = canvas_div\n", - " this._canvas_extra_style(canvas_div)\n", - " this.root.append(canvas_div);\n", - "\n", - " var canvas = $('');\n", - " canvas.addClass('mpl-canvas');\n", - " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", - "\n", - " this.canvas = canvas[0];\n", - " this.context = canvas[0].getContext(\"2d\");\n", - "\n", - " var rubberband = $('');\n", - " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", - "\n", - " var pass_mouse_events = true;\n", - "\n", - " canvas_div.resizable({\n", - " start: function(event, ui) {\n", - " pass_mouse_events = false;\n", - " },\n", - " resize: function(event, ui) {\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " stop: function(event, ui) {\n", - " pass_mouse_events = true;\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " });\n", - "\n", - " function mouse_event_fn(event) {\n", - " if (pass_mouse_events)\n", - " return fig.mouse_event(event, event['data']);\n", - " }\n", - "\n", - " rubberband.mousedown('button_press', mouse_event_fn);\n", - " rubberband.mouseup('button_release', mouse_event_fn);\n", - " // Throttle sequential mouse events to 1 every 20ms.\n", - " rubberband.mousemove('motion_notify', mouse_event_fn);\n", - "\n", - " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", - " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", - "\n", - " canvas_div.on(\"wheel\", function (event) {\n", - " event = event.originalEvent;\n", - " event['data'] = 'scroll'\n", - " if (event.deltaY < 0) {\n", - " event.step = 1;\n", - " } else {\n", - " event.step = -1;\n", - " }\n", - " mouse_event_fn(event);\n", - " });\n", - "\n", - " canvas_div.append(canvas);\n", - " canvas_div.append(rubberband);\n", - "\n", - " this.rubberband = rubberband;\n", - " this.rubberband_canvas = rubberband[0];\n", - " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", - " this.rubberband_context.strokeStyle = \"#000000\";\n", - "\n", - " this._resize_canvas = function(width, height) {\n", - " // Keep the size of the canvas, canvas container, and rubber band\n", - " // canvas in synch.\n", - " canvas_div.css('width', width)\n", - " canvas_div.css('height', height)\n", - "\n", - " canvas.attr('width', width);\n", - " canvas.attr('height', height);\n", - "\n", - " rubberband.attr('width', width);\n", - " rubberband.attr('height', height);\n", - " }\n", - "\n", - " // Set the figure to an initial 600x600px, this will subsequently be updated\n", - " // upon first draw.\n", - " this._resize_canvas(600, 600);\n", - "\n", - " // Disable right mouse context menu.\n", - " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", - " return false;\n", - " });\n", - "\n", - " function set_focus () {\n", - " canvas.focus();\n", - " canvas_div.focus();\n", - " }\n", - "\n", - " window.setTimeout(set_focus, 100);\n", - "}\n", - "\n", - "mpl.figure.prototype._init_toolbar = function() {\n", - " var fig = this;\n", - "\n", - " var nav_element = $('
')\n", - " nav_element.attr('style', 'width: 100%');\n", - " this.root.append(nav_element);\n", - "\n", - " // Define a callback function for later on.\n", - " function toolbar_event(event) {\n", - " return fig.toolbar_button_onclick(event['data']);\n", - " }\n", - " function toolbar_mouse_event(event) {\n", - " return fig.toolbar_button_onmouseover(event['data']);\n", - " }\n", - "\n", - " for(var toolbar_ind in mpl.toolbar_items) {\n", - " var name = mpl.toolbar_items[toolbar_ind][0];\n", - " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", - " var image = mpl.toolbar_items[toolbar_ind][2];\n", - " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", - "\n", - " if (!name) {\n", - " // put a spacer in here.\n", - " continue;\n", - " }\n", - " var button = $('');\n", + " button.click(method_name, toolbar_event);\n", + " button.mouseover(tooltip, toolbar_mouse_event);\n", + " nav_element.append(button);\n", + " }\n", + "\n", + " // Add the status bar.\n", + " var status_bar = $('');\n", + " nav_element.append(status_bar);\n", + " this.message = status_bar[0];\n", + "\n", + " // Add the close button to the window.\n", + " var buttongrp = $('
');\n", + " var button = $('');\n", + " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", + " button.mouseover('Close figure', toolbar_mouse_event);\n", + " buttongrp.append(button);\n", + " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", + " titlebar.prepend(buttongrp);\n", + "}\n", + "\n", + "\n", + "mpl.figure.prototype._canvas_extra_style = function(el){\n", + " // this is important to make the div 'focusable\n", + " el.attr('tabindex', 0)\n", + " // reach out to IPython and tell the keyboard manager to turn it's self\n", + " // off when our div gets focus\n", + "\n", + " // location in version 3\n", + " if (IPython.notebook.keyboard_manager) {\n", + " IPython.notebook.keyboard_manager.register_events(el);\n", + " }\n", + " else {\n", + " // location in version 2\n", + " IPython.keyboard_manager.register_events(el);\n", + " }\n", + "\n", + "}\n", + "\n", + "mpl.figure.prototype._key_event_extra = function(event, name) {\n", + " var manager = IPython.notebook.keyboard_manager;\n", + " if (!manager)\n", + " manager = IPython.keyboard_manager;\n", + "\n", + " // Check for shift+enter\n", + " if (event.shiftKey && event.which == 13) {\n", + " this.canvas_div.blur();\n", + " event.shiftKey = false;\n", + " // Send a \"J\" for go to next cell\n", + " event.which = 74;\n", + " event.keyCode = 74;\n", + " manager.command_mode();\n", + " manager.handle_keydown(event);\n", + " }\n", + "}\n", + "\n", + "mpl.figure.prototype.handle_save = function(fig, msg) {\n", + " fig.ondownload(fig, null);\n", + "}\n", + "\n", + "\n", + "mpl.find_output_cell = function(html_output) {\n", + " // Return the cell and output element which can be found *uniquely* in the notebook.\n", + " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", + " // IPython event is triggered only after the cells have been serialised, which for\n", + " // our purposes (turning an active figure into a static one), is too late.\n", + " var cells = IPython.notebook.get_cells();\n", + " var ncells = cells.length;\n", + " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", + " data = data.data;\n", + " }\n", + " if (data['text/html'] == html_output) {\n", + " return [cell, data, j];\n", + " }\n", + " }\n", + " }\n", + " }\n", + "}\n", + "\n", + "// Register the function which deals with the matplotlib target/channel.\n", + "// The kernel may be null if the page has been refreshed.\n", + "if (IPython.notebook.kernel != null) {\n", + " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", + "}\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "" ], "text/plain": [ "" @@ -3431,10 +3459,10 @@ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 9, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" }