Rayleigh–Ritz method for FEM code

Question

I´m trying to generate a code to model a finite elements analysis based on the following equation, using the Rayleigh Ritz Method. -d/dx(adu/dx)+ c u - q with U(0)=1 and adu/dx=0

i´ve got this far in the code but when it tries to integrate a TypeError: unsupported operand type(s) for /: 'int' and 'function' appears.

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

I don´t know why, does any one knows how to fix it?

Rayleigh–Ritz method for FEM code

0 Answers0