# Mesh: Basic Forward 2D DC ResistivityΒΆ

2D DC forward modeling example with Tensor and Curvilinear Meshes

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 from SimPEG import Mesh, Utils, SolverLU import numpy as np import matplotlib.pyplot as plt def run(plotIt=True): """ Mesh: Basic Forward 2D DC Resistivity ===================================== 2D DC forward modeling example with Tensor and Curvilinear Meshes """ # Step1: Generate Tensor and Curvilinear Mesh sz = [40, 40] tM = Mesh.TensorMesh(sz) rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz, 'rotate')) # Step2: Direct Current (DC) operator def DCfun(mesh, pts): D = mesh.faceDiv sigma = 1e-2*np.ones(mesh.nC) MsigI = mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True) A = -D*MsigI*D.T A[-1, -1] /= mesh.vol[-1] # Remove null space rhs = np.zeros(mesh.nC) txind = Utils.meshutils.closestPoints(mesh, pts) rhs[txind] = np.r_[1, -1] return A, rhs pts = np.vstack((np.r_[0.25, 0.5], np.r_[0.75, 0.5])) # Step3: Solve DC problem (LU solver) AtM, rhstM = DCfun(tM, pts) AinvtM = SolverLU(AtM) phitM = AinvtM*rhstM ArM, rhsrM = DCfun(rM, pts) AinvrM = SolverLU(ArM) phirM = AinvrM*rhsrM if not plotIt: return # Step4: Making Figure fig, axes = plt.subplots(1, 2, figsize=(12*1.2, 4*1.2)) vmin, vmax = phitM.min(), phitM.max() dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True) dat = rM.plotImage(phirM, ax=axes[1], clim=(vmin, vmax), grid=True) cb0 = plt.colorbar(dat[0], ax=axes[0]) cb1 = plt.colorbar(dat[0], ax=axes[1]) cb0.set_label("Voltage (V)") cb1.set_label("Voltage (V)") axes[0].set_title('TensorMesh') axes[1].set_title('CurvilinearMesh') if __name__ == '__main__': run() plt.show()