I have seen other people asking something similar but I can not figure out the problem anyway. I am trying to translate a Matlab code to Python and I have a problem after the following line in the FOR loop:
dx = abs(np.diff(g_coord(num)))
Below you have the code up to that loop. Any help will be appreciated. I really tried to fix it by myself but unsuccessfully. Sorry if it is a stupid mistake. The MATLAB lines are kept as Python comments in case it helps.
import numpy as np
from scipy.sparse import lil_matrix
# physical parameters
seconds_per_yr = 60*60*24*365; # number of seconds in one year
lx = 10000 ; #length of spatial domain (m)
Cp = 1e3 ; # rock heat capacity (J/kg/K)
rho = 2700 ; # rock density (kg/mˆ3)
K = 3.3 ; # bulk thermal conductivity (W/m/K)
kappa = K/(Cp*rho); # thermal diffusivity (mˆ2/s)
Tb = 0 ; # temperatures at boundaries (o C)
A = 2.6e-6 ; # heat production (W/mˆ3)
H = A/(rho*Cp); # heat source term (o K/s) % numerical parameters
dt = 1000*seconds_per_yr ; # time step (s)
ntime = 5000 ; # number of time steps
nels = 40 ; # total number of elements
nod = 2 ; # number of nodes per element
nn = nels+1 # total number of nodes
dx = lx/nels ; # element size
g_coord = np.arange(0, lx+1, dx)#[0:dx:lx]
bcdof = np.array([1, nn]); #[ 1 nn ] ; boundary nodes
bcval = np.array([Tb, Tb]); #[ Tb Tb ] ; # boudary values
g_num = np.zeros((nod, nels), float); #zeros(nod,nels) ;
g_num[0,:]=np.arange(1, nn); #g_num(1,:) = [1:nn-1] ;
g_num[1,:]=np.arange(2, nn+1); #g_num(2,:) = [2:nn] ;
# initialise matrices and vectors
ff = np.zeros((nn,1), float); # system load vector
b = np.zeros((nn,1), float); # system rhs vector
lhs=lil_matrix((nn, nn)) #lhs = sparse(nn,nn); system lhs matrix
rhs=lil_matrix((nn, nn)) #rhs = sparse(nn,nn); system rhs matrix
displ = np.zeros((nn,1), float); # initial temperature (o C)
#-----------------------------------------------------
# matrix assembly
#-----------------------------------------------------
# Matlab version of the loop
#-----------------------------------------------------
#for iel=1:nels # loop over all elements
# num = g_num(:,iel) ; # retrieve equation number
# dx = abs(diff(g_coord(num))) ; # length of element
# MM = dx*[1/3 1/6 ; 1/6 1/3 ] ;# mass matrix
# KM = [kappa/dx -kappa/dx ; -kappa/dx kappa/dx ]; #diffn matrix
# F = dx*H*[1/2 ; 1/2] ; # load vector
# lhs(num,num) = lhs(num,num) + MM/dt + KM ; # assemble lhs
# rhs(num,num) = rhs(num,num) + MM/dt ; # assemble rhs
# ff(num) = ff(num) + F ; # assemble load
#end # end of element loop
#Python version of the loop
#-----------------------------------------------------
for iel in range(0, nels): # loop over all elements
num = g_num[:,iel] # retrieve equation number
#print(num)
dx = abs(np.diff(g_coord[num])) # length of element
MM = dx*(np.array([[1/3, 1/6],[1/6, 1/3]])) # mass matrix
KM = np.array([[kappa/dx, -kappa/dx],[-kappa/dx, kappa/dx]])
F = dx*H*(np.array([1/2, 1/2])).reshape(-1,1) # load vector
lhs[num,num] = lhs[num,num] + MM/dt + KM # assemble lhs
rhs[num,num] = rhs[num,num] + MM/dt # assemble rhs
ff[num] = ff[num] + F # assemble load
callableerror occurred when you tried to index an array with()(matlab style), right? You correct that, and you get thisintegererror? New versions ofnumpyhave gotten picky about the index type. Why did you intializeg_numasfloat?lilsparse matrix - this is a start. But collecting element values incoostyle inputs is probably faster. At least that's how I used to do with in MATLAB. Duplicates incooinputs are summed when converted tocsrfor calculation purposes.