import numpy as np
import A
nels=4
np_types=2
loaded_nodes=4
fixed_freedoms=2
nod=2
nprops=1
nn=nels+1
nr=2
nodof=2
ndof=nod*nodof
ell=np.array([[2.5],[2.5],[3.0],[2.0]])
# 初始化定义数组
nf=np.ones((nodof,nn),dtype=np.int64)
g=np.ones((ndof,1),dtype=np.int64)
num=np.ones((nod,1),dtype=np.int64)
etype=np.ones((nels,1),dtype=np.int64)
eld=np.ones((ndof,1))
km=np.ones((ndof,ndof))
mm=np.zeros((ndof,ndof))
action=np.ones((ndof,1))
g_g=np.ones((ndof,nels))
prop=np.ones((nprops,np_types))
if np_types==1:
etype[:,0]=1
else:
etype[:,0]=(1,1,2,2)
if np_types==1:
prop[:nprops,np_types-1]=5e5
else:
prop[0,:]=(4.0e4,2.0e4)
# 读取 nr
if nr!=0:
Dim_1=[1,4]
nf_value=np.array([[0,0],[1,1]])
m=0
for i in Dim_1:
for j in range(1,nodof+1):
nf[j-1,i-1]=nf_value[j-1,m]
m=m+1
#form nf
A.formnf(nf)
neq=int(max(nf.reshape(nf.shape[0]*nf.shape[1],1)))
kdiag=np.zeros((neq,1),dtype=np.int64)
loads=np.zeros((neq+1,1))
## 累加单元去发现全局尺寸
for iel in range(1,nels+1):
num[:nod,0]=(iel,iel+1)
A.num_to_g(num,nf,g)
g_g[:,iel-1]=g[:,0]
#call fkidiag
A.fkdiag(kdiag,g)
print(neq)
for i in range(1,neq):
kdiag[i]=kdiag[i]+kdiag[i-1]
kv=np.zeros((kdiag[neq-1,0],1),dtype=float)
print(' 一共有 '+str(neq)外汇跟单gendan5.com+' 等式和天际线存储个数为 '+str(kdiag[neq-1]))
# 全局刚度矩阵安装
for iel in range(1,nels+1):
#call pin_jointed
A.beam_km(km,prop[0,etype[iel-1,0]-1],ell[iel-1,0])
g[:,0]=g_g[:,iel-1]
if nprops>1:
A.beam_mm(mm,prop[1,etype[iel-1,0]-1],ell[iel-1,0])
#call fsparv
A.fsparv(kv,km+mm,g,kdiag)
### 读荷载和位移
if loaded_nodes!=0:
Dim_2 = [2,3,4,5]
load_value=np.array([[-20.0,-6.0,-8.8,-1.2],[0.0,-3.0,2.2,0.5333]])
m=0
for i in Dim_2:
for j in range(1,nodof+1):
loads[nf[j-1,i-1]-1]=load_value[j-1,m]
m=m+1
if fixed_freedoms!=0:
node=np.ones((fixed_freedoms,1),dtype=np.int64)
sense=np.ones((fixed_freedoms,1),dtype=np.int64)
value=np.ones((fixed_freedoms,1))
no=np.ones((fixed_freedoms,1),dtype=np.int64)
node[:,0]=(1,3)
sense[:,0]=(2,1)
value[:,0]=(-0.001,-0.005)
#### 位移赋值
for i in range(1,fixed_freedoms+1):
no[i-1]=nf[sense[i-1]-1,node[i-1]-1]
kv[kdiag[no-1]-1]=kv[kdiag[no-1]-1]+1e20
loads[no-1,0]=kv[kdiag[no-1,0]-1,0]*value
##### 方程求解
#call sparin
A.sparin(kv,kdiag)
##call spabac
A.spabac(kv,loads,kdiag)
loads[neq]=0
print(' 节点 位移 转角 ')
for k in range(1,nn+1):
print(k,end=' ')
for m in range(1,nodof+1):
print('{:9.4e}'.format(loads[nf[m-1,k-1]-1,0]),end=' ')
print( )
## 重新取得单元和节点力
print(' 单元 力 力矩 力 力矩 ')
axial=[0]
for iel in range(1,nels+1):
#call pin_jointed
#call pin_jointed
A.beam_km(km,prop[0,etype[iel-1,0]-1],ell[iel-1,0])
g[:,0]=g_g[:,iel-1]
if nprops>1:
A.beam_mm(mm,prop[1,etype[iel-1,0]-1],ell[iel-1,0])
##########
eld=loads[g-1,0]
action=np.dot(km+mm,eld)
print(iel,end=' ')
for i in range(1,ndof+1):
print('{:9.4e}'.format(action[i-1,0]),end=' ')
print( )