This notebook reads files describing a structure, and the files output by Frame2D after an analysis, and checks that the member end forces computed here from the displacements and member loads agree with those computed by Frame2D.
It does this in the simplest way possible, using quite different logic than Frame2D, resulting in a higher degree of confidence in the results. It would have been better had someone else programmed it, but oh well ...
In [1]:
ds = 'KG82'
lcase = 'all'
#ds = 'l22x6'
#lcase = 'Case-2b'
def filename(basename,lc=None):
if lc is not None:
basename = lc + '/' + basename
return 'data/' + ds + '.d/' + basename + '.csv'
def Warn(msg):
print('!!!!! Warning: {}'.format(msg))
In [2]:
import pandas as pd
import math
In [3]:
class Node(object):
def __init__(self,id,x,y):
self.id = id
self.x = x
self.y = y
self.deltaX = 0.
self.deltaY = 0.
self.thetaZ = 0.
In [4]:
ntable = pd.read_csv(filename('nodes'))
NODES = {}
for i,n in ntable.iterrows():
if n.NODEID in NODES:
Warn("Node '{}' is multiply defined.".format(n.NODEID))
NODES[n.NODEID] = Node(n.NODEID,float(n.X),float(n.Y))
#ntable
In [5]:
dtable = pd.read_csv(filename('node_displacements',lcase))
for i,n in dtable.iterrows():
node = NODES[n.NODEID]
node.deltaX = float(n.DX)
node.deltaY = float(n.DY)
node.thetaZ = float(n.RZ)
dtable
Out[5]:
In [6]:
pd.DataFrame([vars(v) for v in NODES.values()]).set_index('id')
Out[6]:
In [7]:
class Member(object):
E = 200000.
def __init__(self,id,nodej,nodek):
self.id = id
self.nodej = nodej
self.nodek = nodek
dx = nodek.x - nodej.x
dy = nodek.y - nodej.y
self.L = L = math.sqrt(dx*dx + dy*dy)
self.cosx = dx/L
self.cosy = dy/L
self.Ix = 0.
self.A = 0.
self.loads = []
self.releases = set()
for a in 'FXJ FXK FYJ FYK MZJ MZK'.split():
setattr(self,a,0.)
In [8]:
table = pd.read_csv(filename('members'))
MEMBERS = {}
for i,m in table.iterrows():
if m.MEMBERID in MEMBERS:
Warn("Member '{}' is multiply defined.".format(m.MEMBERID))
MEMBERS[m.MEMBERID] = Member(m.MEMBERID,NODES[m.NODEJ],NODES[m.NODEK])
In [9]:
import sst
SST = sst.SST()
table = pd.read_csv(filename('properties'))
defIx = defA = None
for i,row in table.iterrows():
if not pd.isnull(row.SIZE):
defIx,defA = SST.section(row.SIZE,'Ix,A')
memb = MEMBERS[row.MEMBERID]
memb.Ix = float(defIx if pd.isnull(row.IX) else row.IX)
memb.A = float(defA if pd.isnull(row.A) else row.A)
if not pd.isnull(row.IX):
defIx = row.IX
if not pd.isnull(row.A):
defA = row.A
In [10]:
try:
lctable = pd.read_csv(filename('load_combinations'))
use_all = False
COMBO = {}
for i,row in lctable.iterrows():
if row.CASE == lcase:
COMBO[row.LOAD.lower()] = row.FACTOR
except OSError:
use_all = True
COMBO = None
COMBO
In [11]:
table = pd.read_csv(filename('member_loads'))
for i,row in table.iterrows():
memb = MEMBERS[row.MEMBERID]
typ = row.TYPE
f = 1.0 if use_all else COMBO.get(row.LOAD.lower(),0.)
if f != 0.:
w1 = None if pd.isnull(row.W1) else (float(row.W1)*f)
w2 = None if pd.isnull(row.W2) else (float(row.W2)*f)
a = None if pd.isnull(row.A) else float(row.A)
b = None if pd.isnull(row.B) else float(row.B)
c = None if pd.isnull(row.C) else float(row.C)
memb.loads.append((typ,w1,w2,a,b,c))
#MEMBERS['LC'].loads
In [12]:
table = pd.read_csv(filename('releases'))
for i,row in table.iterrows():
memb = MEMBERS[row.MEMBERID]
memb.releases.add(row.RELEASE.upper())
In [13]:
t = pd.DataFrame([vars(v) for v in MEMBERS.values()]).set_index('id')
del t['nodej']
del t['nodek']
del t['loads']
t
Out[13]:
In [14]:
MEFS = pd.read_csv(filename('member_end_forces',lcase)).set_index('MEMBERID')
MEFS
Out[14]:
In [15]:
cols = 'FXJ FXK FYJ FYK MZJ MZK'.split()
for m in MEMBERS.values():
for a in cols:
setattr(m,a,0.)
# difference in end displacements, global coords
dX = m.nodek.deltaX - m.nodej.deltaX
dY = m.nodek.deltaY - m.nodej.deltaY
# axial deformation / force:
ldX = dX*m.cosx + dY*m.cosy
T = m.E*m.A*ldX/m.L
m.FXK += T
m.FXJ += -T
#print(m.id,ldX,T)
# shear deformation / force:
vdY = dY*m.cosx - dX*m.cosy
M = -6.*m.E*m.Ix*vdY/(m.L*m.L)
V = 2.*M/m.L
m.MZJ += M
m.MZK += M
m.FYJ += V
m.FYK += -V
#print(m.id,vdY,M,V)
# end rotations / moments:
MJ = (m.E*m.Ix/m.L)*(4.*m.nodej.thetaZ + 2.*m.nodek.thetaZ)
MK = (m.E*m.Ix/m.L)*(2.*m.nodej.thetaZ + 4.*m.nodek.thetaZ)
VJ = (MJ+MK)/m.L
m.MZJ += MJ
m.MZK += MK
m.FYJ += VJ
m.FYK += -VJ
#print(m.id,m.nodej.thetaZ,m.nodek.thetaZ,MJ,MK,VJ)
# applied loads: fixed-end moments and shears:
for ltype,w1,w2,a,b,c in m.loads:
mj = mk = 0.
vj = vk = 0.
if ltype == 'PL':
b = m.L - a
P = w1
mj = -P*a*b*b/(m.L*m.L)
mk = P*b*a*a/(m.L*m.L)
vj = (-P*b + nj + mk)/m.L
vk = -P - vj
elif ltype == 'UDL':
mj = -w1*m.L**2/12.
mk = -mj
vj = -w1*m.L/2.
vk = vj
else:
Warn("Load type '{}' not implemented here ...".format(ltype))
continue
m.MZJ += mj
m.MZK += mk
m.FYJ += vj
m.FYK += vk
# member end moment releases:
relc = m.releases.copy()
if 'MZJ' in m.releases:
mj = -m.MZJ
mk = 0. if 'MZK' in m.releases else 0.5*mj
vj = (mj + mk)/m.L
##print(m.id,'MZJ',m.MZJ,m.MZK,mj,mk,rel)
m.MZJ += mj
m.MZK += mk
m.FYJ += vj
m.FYK += -vj
relc.remove('MZJ')
if 'MZK' in m.releases:
mk = -m.MZK
mj = 0. if 'MZJ' in m.releases else 0.5*mk
vj = (mj + mk)/m.L
##print(m.id,'MZK',m.MZJ,m.MZK,mj,mk,rel)
m.MZJ += mj
m.MZK += mk
m.FYJ += vj
m.FYK += -vj
relc.remove('MZK')
if relc:
Warn("Member end-releases not processed: {}".format(relc))
In [16]:
computed = pd.DataFrame([{k:getattr(m,k) for k in ['id']+cols}
for m in MEMBERS.values()]).set_index('id')
diff = (computed - MEFS[cols])
lim = 1E-12
for c in cols:
biggest = MEFS[c].abs().max()
diff[c][diff[c].abs() < biggest*lim] = 0
diff
Out[16]:
In [17]:
diff.abs().max()
Out[17]:
In [18]:
diff = (computed - MEFS[cols])
for c in cols:
biggest = MEFS[c].abs().max()
r = (diff[c]/biggest)
idr = r.abs().idxmax()
print(c,idr,r[idr])
In [ ]: