Численные расчеты с python

def DNx(self, Xi, Eta, Zeta, i):
        '''DNx - производная от функция формы fNi по x в лок. коорд.'''
        return self.JcbInv[0, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[0, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[0, 2] * self.DNZeta(Xi, Eta, i)

    def DNy(self, Xi, Eta, Zeta, i):
        '''DNy - производная от функция формы fNi по y в лок. коорд.'''
        return self.JcbInv[1, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[1, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[1, 2] * self.DNZeta(Xi, Eta, i)

    def DNz(self, Xi, Eta, Zeta, i):
        '''DNz - производная от функция формы fNi по z в лок. коорд.'''
        return self.JcbInv[2, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[2, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[2, 2] * self.DNZeta(Xi, Eta, i)

    def BMat(self, Xi, Eta, Zeta):
        '''Подпрограмма для вычисления матрицы связи (B) между перемещениями и
        деформациями (eps=B*u).'''
        nFree = self.nFree
        self.B.fill(0)
        for i in xrange(self.eNode):
            dnx = self.DNx(Xi, Eta, Zeta, i)
            dny = self.DNy(Xi, Eta, Zeta, i)
            dnz = self.DNz(Xi, Eta, Zeta, i)
            self.B[0, nFree * i] = dnx
            self.B[1, nFree * i + 1] = dny
            self.B[2, nFree * i + 2] = dnz
            self.B[3, nFree * i] = dny
            self.B[3, nFree * i + 1] = dnx
            self.B[4, nFree * i + 1] = dnz
            self.B[4, nFree * i + 2] = dny
            self.B[5, nFree * i] = dnz
            self.B[5, nFree * i + 2] = dnx

# Скользящее среднее с окном размером 2.

a = numpy.random.random(10000)
res = numpy.empty(len(a) - 1)

# Медленно:
for i in range(len(a) - 1):
    res[i] = 0.5*(a[i] + a[i+1])
# Быстрee:
res = 0.5*(a[1:] + a[:-1])

class Solid8():
    def __init__(self):
        self.nDef = 6
        self.nFree = 3
        self.eNode = 8
        self.eInt = 3
        self.Xi0 = [-1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0]
        self.Eta0 = [-1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0]
        self.Zeta0 = [-1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0]
        self.Jcb = np.empty((self.nFree, self.nFree))
        self.JcbInv = np.empty((self.nFree, self.nFree))
        self.B = np.zeros((self.nDef, self.nFree * self.eNode))

    def fNi(self, Xi, Eta, Zeta, i):
        return (1.0 + Xi * self.Xi0[i]) * (1.0 + Eta * self.Eta0[i])\
            * (1.0 + Zeta * self.Zeta0[i]) / 8.0

    def DNXi(self, Eta, Zeta, i):
        return self.Xi0[i] * (1.0 + Eta * self.Eta0[i])\
            * (1.0 + Zeta * self.Zeta0[i]) / 8.0

    def DNEta(self, Xi, Zeta, i):
        return self.Eta0[i] * (1.0 + Xi * self.Xi0[i])\
            * (1.0 + Zeta * self.Zeta0[i]) / 8.0

    def DNZeta(self, Xi, Eta, i):
        return self.Zeta0[i] * (1.0 + Xi * self.Xi0[i])\
            * (1.0 + Eta * self.Eta0[i]) / 8.0

    def Jacobian(self, Xi, Eta, Zeta, x0, y0, z0):
        self.Jcb = np.zeros((self.nFree, self.nFree))

        # Вычисление Jcb
        for i in xrange(self.eNode):
            self.Jcb[0, 0] += x0[i] * self.DNXi(Eta, Zeta, i)
            self.Jcb[0, 1] += y0[i] * self.DNXi(Eta, Zeta, i)
            self.Jcb[0, 2] += z0[i] * self.DNXi(Eta, Zeta, i)
            self.Jcb[1, 0] += x0[i] * self.DNEta(Xi, Zeta, i)
            self.Jcb[1, 1] += y0[i] * self.DNEta(Xi, Zeta, i)
            self.Jcb[1, 2] += z0[i] * self.DNEta(Xi, Zeta, i)
            self.Jcb[2, 0] += x0[i] * self.DNZeta(Xi, Eta, i)
            self.Jcb[2, 1] += y0[i] * self.DNZeta(Xi, Eta, i)
            self.Jcb[2, 2] += z0[i] * self.DNZeta(Xi, Eta, i)

        # Вычисление определителя
        det = np.linalg.det(self.Jcb)

        # Вычисление обратной матрицы JcbInv
        self.JcbInv = np.linalg.inv(self.Jcb)
        return det

    def DNx(self, Xi, Eta, Zeta, i):
        return self.JcbInv[0, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[0, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[0, 2] * self.DNZeta(Xi, Eta, i)

    def DNy(self, Xi, Eta, Zeta, i):
        return self.JcbInv[1, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[1, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[1, 2] * self.DNZeta(Xi, Eta, i)

    def DNz(self, Xi, Eta, Zeta, i):
        return self.JcbInv[2, 0] * self.DNXi(Eta, Zeta, i)\
            + self.JcbInv[2, 1] * self.DNEta(Xi, Zeta, i)\
            + self.JcbInv[2, 2] * self.DNZeta(Xi, Eta, i)

    def BMat(self, Xi, Eta, Zeta):
        nFree = self.nFree
        self.B.fill(0)
        for i in xrange(self.eNode):
            dnx = self.DNx(Xi, Eta, Zeta, i)
            dny = self.DNy(Xi, Eta, Zeta, i)
            dnz = self.DNz(Xi, Eta, Zeta, i)
            self.B[0, nFree * i] = dnx
            self.B[1, nFree * i + 1] = dny
            self.B[2, nFree * i + 2] = dnz
            self.B[3, nFree * i] = dny
            self.B[3, nFree * i + 1] = dnx
            self.B[4, nFree * i + 1] = dnz
            self.B[4, nFree * i + 2] = dny
            self.B[5, nFree * i] = dnz
            self.B[5, nFree * i + 2] = dnx

elem = Solid8()
x0 = [0, 1, 1, 0, 0, 1, 1, 0]
y0 = [0, 0, 1, 1, 0, 0, 1, 1]
z0 = [0, 0, 0, 0, 1, 1, 1, 1]
det = elem.Jacobian(0, 0, 0, x0, y0, z0)
elem.BMat(1, 2, 3)
print(det)
print(elem.B)