python - 类对象内嵌套 for 循环和 3D 绘图

Question

我确信这是一个很容易处理的问题，但我无法弄清楚。我创建了一个钻孔类，想要计算每个钻孔/井周围的孔隙压力。沿着单个轴，我的代码如下所示:

from scipy.special import *
import matplotlib.pyplot as plt
import numpy as np
from math import *


## Globale Variablen ##

rhof = 1000                                     # Dichte Flüssigkeit [kg/m³]
lameu = 11.2*10**9                              # Lamé-Parameter, undrained [GPa]
lame = 8.4*10**9                                # Lamé-Parameter, drained [GPa]
pi                                              # durch Pythonmodul "math" gegeben
alpha = 0.65                                    # Biot-Willis-Koeffizient
G = 8.4*10**9                                   # Schermodul [GPa]
k = 1.0e-15                                     # Permeabilität [m²] bzw. [Darcy] 
eta = 0.001                                     # Viskosität des Fluids [Pa*s]

## Berechnung der Parameter ##

kappa = k/eta                                                    
c = ((kappa*(lameu-lame)*(lame+2*G))/((alpha**2)*(lameu+2*G)))   

## Wertebereich ##

xmin = 0
xmax = 100
xsteps = 1.0
x = np.arange(xmin, xmax, xsteps)

## Class ##

class Bohrloch(object):
    loch_zaehler = 0

    def __init__(self, xlage, tstart, q):    # Funktion, um BL zu erzeugen
        self.xlage = xlage
        #self.ylage = ylage                  # Lage der Bohrung
        self.tstart = tstart                 # Start der Injektion/Produktion
        self.q = q                           # Fluidmenge


## Druck ##     

    def getPressure(self, t): # gibt nach Zeit t die zugehörigen Druckwerte aus
        if (t-self.tstart<0): # Fehlermeldung, falls Startpunkt nach t liegt
            return () 
            print "Startpunkt liegt außerhalb des Förderzeitraumes!"
        else:
            self.r = np.sqrt((x-self.xlage)**2)
            self.P = (self.q/(rhof*4*pi*kappa))*(expn(1,self.r**2/(4*c*(t-self.tstart))))
            #self.P[self.xlage] = 0         # gibt Bohrlochlage wieder
            self.z = self.P/1e6
            return self.z                   # Druckwerte in [MPa]


    def pressureTable (self, t, xschritt):       # erstellt Wertetabelle
        self.getPressure(t)
        for i in range (xmin, xmax, xschritt):
            print i, "  ", self.z[i]

t = 1000*24*3600
b1 = Bohrloch(50,0*24*3600,6.0/1000)
b1.pressureTable(t,1)

通过这种方法我得到了我想要的压力表。现在我想要一个 x 和 y 值的压力表，包括 3D 图。这是我到目前为止的代码:

from mpl_toolkits.mplot3d import Axes3D 
from matplotlib import cm
from scipy.special import *
import matplotlib.pyplot as plt
import numpy as np
from math import *


## Globale Variablen ##

rhof = 1000                                     # Dichte Flüssigkeit [kg/m³]
lameu = 11.2*10**9                              # Lamé-Parameter, undrained [GPa]
lame = 8.4*10**9                                # Lamé-Parameter, drained [GPa]
pi                                              # durch Pythonmodul "math" gegeben
alpha = 0.65                                    # Biot-Willis-Koeffizient
G = 8.4*10**9                                   # Schermodul [GPa]
k = 1.0e-15                                     # Permeabilität [m²] bzw. [Darcy] 
eta = 0.001                                     # Viskosität des Fluids [Pa*s]

## Berechnung der Parameter ##

kappa = k/eta                                                     
c = ((kappa*(lameu-lame)*(lame+2*G))/((alpha**2)*(lameu+2*G)))   

## Wertebereich ##

xmin = 0
xmax = 100
xsteps = 1.0
x = np.arange(xmin,xmax,xsteps)
ymin = 0
ymax = 100
ysteps = 1.0
y = np.arange(ymin,ymax,ysteps)

## Klassendefinition ## 

class Bohrloch(object):
    loch_zaehler = 0

    def __init__(self, xlage, ylage, tstart, q):     # Funktion, um BL zu erzeugen
        self.xlage = xlage                   # x-Lage der Bohrung
        self.ylage = ylage                   # y-Lage der Bohrung
        self.tstart = tstart                 # Start der Injektion/Produktion
        self.q = q                           # Fluidmenge


## Druck ##      

    def getPressure(self, t):                
        if (t-self.tstart<0):                
            return () 
            print "Startpunkt liegt außerhalb des Förderzeitraumes!"
        else:
            self.r = np.sqrt((x-self.xlage)**2+(y-self.ylage)**2)
            self.P = (self.q/(rhof*4*pi*kappa))*(expn(1,self.r**2/(4*c*(t-self.tstart))))
            self.P[self.xlage] = np.nan
            self.P[self.ylage] = np.nan            
            self.z = self.P/1e6
            return self.z                    # Druckwerte in [MPa]

    def pressureTable (self, t, xschritt, yschritt):
        self.getPressure(t)
        for k in range (xmin, xmax, xschritt):
            for l in range (ymin, ymax, yschritt):
                # my mistake should be here?
                print k, "  ",  l, "  ",  self.z[k][l]                   

    def pressurePlot3D (self, t):
        self.getPressure(t)
        Z = self.z
        X, Y = np.meshgrid(x,y)
        Z[Z == np.inf] = np.nan                              
        fig = plt.figure()
        ax = fig.gca(projection='3d')
        surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet, linewidth=0,
                               antialiased=False, vmin=np.nanmin(Z), vmax=np.nanmax(Z))
        fig.colorbar(surf, shrink=0.5, aspect=5)
        ax.set_xlim(xmin,xmax)      # x-Achsenskala vorgeben
        ax.set_ylim(ymin,ymax)      # y-Achsenskala vorgeben
        ax.set_title('Druckverteilung')
        ax.set_xlabel('x-Richtung [m]')
        ax.set_ylabel('y-Richtung Well [m]')
        ax.set_zlabel('Druck in [MPa]')
        plt.show()

t = 1000*24*3600
b1 = Bohrloch(50,50,0*24*3600,6.0/1000)
b1.pressureTable(t,1)
b1.pressurePlot3D(t)

不幸的是，我的 table 不起作用，并且所需的 3D 绘图看起来很奇怪。我仍然是 Python 的初学者，需要一些建议。

有人可以帮忙吗？

Answer 1

问题是 self.z 不是二维数组/列表。因此，尝试访问 self.z[k][l] 会导致 IndexError:标量变量的索引无效。

<小时/>

• 我不太明白你想如何实现第二个维度。您引入了 y 位置，但是随后，您只需使用

  中的 x 和 y 位置来计算一维半径数组

  self.r = np.sqrt((x-self.xlage)**2+(y-self.ylage)**2)
  

• 下一个问题是，您打算做什么:

  self.P[self.xlage] = np.nan
  self.P[self.ylage] = np.nan
  

  如果将xsteps和ysteps更改为10，并调用:

  b1 = Bohrloch(2,3,0*24*3600,6.0/1000)
  print b1.getPressure(t)  
  

  您的输出将是:

  [ 5.44152501 4.40905986        nan        nan  2.87481753  2.64950827
    2.46756653 2.31503845 2.18379093 2.06866598]
  

  为什么要用 nan 替换第三个和第四个元素？

<小时/>

这些问题也是您绘图例程的基础。由于您的数组中现在有 np.nan 值，因此这些值不会显示在图中。由于 self.z 不是二维的，因此您可能无法获得您期望的表面:

<小时/>

这是一种实现 2D 实现的简单方法。我不太熟悉你想要做什么，但它传达了这个想法:

    def getPressure(self, t):                
        if (t-self.tstart<0):                
            return () 
            print "Startpunkt liegt außerhalb des Förderzeitraumes!"
        else:

            # you need to initialize r, P and z as list of lists
            # make this dependent on your x coordinates
            # the second dimension will grow dynamically

            self.r = [[] for ri in range(len(x))]
            self.P = [[] for ri in range(len(x))]
            self.z = [[] for ri in range(len(x))]

            # iterate through both x and y independently

            for ii in range(len(x)):
                for jj in range(len(y)):

            # append to the list that corresponds to the current x -value
            # also, use x[ii] and y[jj] to call one x-, y-value at a time  

                    self.r[ii].append(np.sqrt((x[ii]-self.xlage)**2+(y[jj]-self.ylage)**2))

            # calling r[ii][-1] ensures you are using the value that was last added to the list:

                    self.P[ii].append((self.q/(rhof*4*pi*kappa))*(expn(1,self.r[ii][-1]**2/(4*c*(t-self.tstart)))))

                    self.z[ii].append(self.P[ii][-1]/1e6)

            # now, you can use xlage and ylage to blank one value
            # do this for both P and z, because z is now calculated inside the loop

            self.P[self.xlage][self.ylage] = np.nan
            self.z[self.xlage][self.ylage] = np.nan

            return self.z

从绘图例程中删除此行:Z[Z == np.inf] = np.nan，使用原始命令:

b1 = Bohrloch(50,50,0*24*3600,6.0/1000)  
b1.pressurePlot3D(t)  

现在你将得到这个图: