1

我正在尝试使用 Rayleigh Ritz 方法生成代码以基于以下方程对有限元分析进行建模。-d/dx(a du/dx)+ cu - q 其中 U(0)=1 且du/dx=0

我已经在代码中做到了这一点，但是当它尝试集成 TypeError: unsupported operand type(s) for /: 'int' and 'function' 时出现。

L = 1.0
n = 4
dx = L/n
dofpe = 2 # Grados de libertad por elemento, por ahora solo programado 2
Uo=1 #condicion esccencial
Qo=0 #condicion natural

coords = np.linspace(0,L,n+1)
elementos = np.zeros((n,dofpe),dtype=int)

for i in range(n):
    elementos[i,0] = int(i)
    elementos[i,1] = int(i+1)

def basis(xbar,i,he):
    """Returns the ith basis funciton at x
    xbar - local element coordinate
    i - basis function index
    he - width of each element
    """
    if i == 0: return 1 - xbar/he
    if i == 1: return xbar/he
    return -1 # Este seria un codigo de error, toca mejorar

def dbasis(xbar,i,he):
    """Returns the ith basis funciton at x
    xbar - local element coordinate
    i - basis funciton index
    he - width of each element
    """
    if i == 0: return -1/he
    if i == 1: return 1/he
    return -1 # Este seria un codigo de error, toca mejorar

def a():
    return 1

def c():
    return 1

def q():
    return 1

xvec = np.linspace(0,1,10)
basisv = np.vectorize(basis)
plt.plot(xvec,basisv(xvec,0,1))
dbasisv = np.vectorize(dbasis)
plt.plot(xvec,dbasisv(xvec,0,1))




K = np.zeros((n+1,n+1))



 for e in range(n): # loop over elements
  dx = coords[elementos[e,1]] - coords[elementos[e,0]]
  for i in range(dofpe): # loop over basis i
        for j in range(dofpe): # loop over basis j
            def fun(x,i,j,dx,a,c):
             return a*dbasis(x,i,dx)*dbasis(x,j,dx)+ c*basis(x,i,dx)*basis(x,j,dx)  
            kij=scipy.integrate.quad(fun,coords[elementos[e,0]],coords[elementos[e,1]],args=(i,j,a,c,dx))
            K[elementos[e,i],elementos[e,j]] +=kij



F = np.zeros((n+1))
for e in range(n):
    dx = coords[elementos[e,1]] - coords[elementos[e,0]]
    for i in range(dofpe):
        def l(x,i,dx,q):
         return q*basis(x,i,dx)
        fij=scipy.integrate.quad(l,coords[elementos[e,0]],coords[elementos[e,1]],args=(i,dx,q))
        F[elementos[e,i]] += fij

不知道为什么，有人知道怎么解决吗？

4

0 回答 0