在类Obj中嵌套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绘图。这是我目前为止的代码：

^{pr2}$

不幸的是，我的表格不工作，而且所需的三维绘图看起来很奇怪。我仍然是Python的初学者，需要一些建议。在

有人能帮忙吗？在

Answer 1

问题是self.z不是二维数组/列表。因此，尝试访问self.z[k][l]会导致IndexError: invalid index to scalar variable.

---

• 我不太明白你想如何实现第二个维度。您引入了y位置，但是，您只需使用中的x和y位置来计算一维半径数组

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

• 下一个问题是，你打算做什么：

  ^{pr2}$

  如果您将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)  

现在你将得到这个情节：