## Fixed End Forces

This module computes the fixed end forces (moments and shears) due to transverse loads acting on a 2-D planar structural member.

``````

In [1]:

import numpy as np
import sys
from salib import extend

``````

### Class EF

Instances of class EF represent the 6 end-forces for a 2-D planar beam element. The forces (and local degrees of freedom) are numbered 0 through 5, and are shown here in their positive directions on a beam-element of length L. The 6 forces are labelled by prefixing the number with a letter to suggest the normal interpretation of that force: c for axial force, v for shear force, and m for moment.

For use in this module, the end forces will be fixed-end-forces.

``````

In [2]:

class EF(object):

"""Class EF represents the 6 end forces acting on a 2-D, planar, beam element."""

def __init__(self,c0=0.,v1=0.,m2=0.,c3=0.,v4=0.,m5=0.):
"""Initialize an instance with the 6 end forces.  If the first
argument is a 6-element array, initialize from a copy of that
array and ignore any other arguments."""
if np.isscalar(c0):
self.fefs = np.matrix([c0,v1,m2,c3,v4,m5],dtype=np.float64).T
else:
self.fefs = c0.copy()

def __getitem__(self,ix):
"""Retreive one of the forces by numer.  This allows allows unpacking
of all 6 end forces into 6 variables using something like:
c0,v1,m2,c3,v4,m5 = self
"""
return self.fefs[ix,0]

"""Add this set of end forces to another, returning the sum."""
assert type(self) is type(other)
new = self.__class__(self.fefs+other.fefs)
return new

def __sub__(self,other):
"""Subtract the other from this set of forces, returning the difference."""
assert type(self) is type(other)
new = self.__class__(self.fefs-other.fefs)
return new

def __mul__(self,scale):
"""Multiply this set of forces by the scalar value, returning the product."""
if scale == 1.0:
return self
return self.__class__(self.fefs*scale)

__rmul__ = __mul__

def __repr__(self):
return '{}({},{},{},{},{},{})'.format(self.__class__.__name__,*list(np.array(self.fefs.T)[0]))

``````
``````

In [3]:

##test:
f = EF(1,2,0,4,1,6)
f

``````
``````

Out[3]:

EF(1.0,2.0,0.0,4.0,1.0,6.0)

``````
``````

In [4]:

##test:
g = f+f+f
g

``````
``````

Out[4]:

EF(3.0,6.0,0.0,12.0,3.0,18.0)

``````
``````

In [5]:

##test:
f[1]

``````
``````

Out[5]:

2.0

``````
``````

In [6]:

##test:
f[np.ix_([3,0,1])]

``````
``````

Out[6]:

matrix([[ 4.,  1.,  2.]])

``````
``````

In [7]:

##test:
g[(3,0,1)]

``````
``````

Out[7]:

matrix([[ 12.],
[  3.],
[  6.]])

``````
``````

In [8]:

##test:
f0,f1,f2,f3,f4,f5 = g
f3

``````
``````

Out[8]:

12.0

``````
``````

In [9]:

##test:
g, g*5, 5*g

``````
``````

Out[9]:

(EF(3.0,6.0,0.0,12.0,3.0,18.0),
EF(15.0,30.0,0.0,60.0,15.0,90.0),
EF(15.0,30.0,0.0,60.0,15.0,90.0))

``````

Now define properties so that the individual components can be accessed like name atrributes, eg: '`ef.m3`' or '`ef.m5 = 100.`'.

``````

In [10]:

@extend
class EF:

@property
def c0(self):
return self.fefs[0,0]

@c0.setter
def c0(self,v):
self.fefs[0,0] = v

@property
def v1(self):
return self.fefs[1,0]

@v1.setter
def v1(self,v):
self.fefs[1,0] = v

@property
def m2(self):
return self.fefs[2,0]

@m2.setter
def m2(self,v):
self.fefs[2,0] = v

@property
def c3(self):
return self.fefs[3,0]

@c3.setter
def c3(self,v):
self.fefs[3,0] = v

@property
def v4(self):
return self.fefs[4,0]

@v4.setter
def v4(self,v):
self.fefs[4,0] = v

@property
def m5(self):
return self.fefs[5,0]

@m5.setter
def m5(self,v):
self.fefs[5,0] = v

``````
``````

In [11]:

##test:
f = EF(10.,11,12,13,15,15)
f, f.c0, f.v1, f.m2, f.c3, f.v4, f.m5

``````
``````

Out[11]:

(EF(10.0,11.0,12.0,13.0,15.0,15.0), 10.0, 11.0, 12.0, 13.0, 15.0, 15.0)

``````
``````

In [12]:

##test:
f.c0 *= 2
f.v1 *= 3
f.m2 *= 4
f.c3 *= 5
f.v4 *= 6
f.m5 *= 7
f

``````
``````

Out[12]:

EF(20.0,33.0,48.0,65.0,90.0,105.0)

``````

This is the base class for all the different types of member loads (point loads, UDLs, etc.) of 2D planar beam elements.

The main purpose is to calculate the fixed-end member forces, but we will also supply logic to enable calculation of internal shears and moments at any point along the span.

All types of member loads will be input using a table containing five data columns: W1, W2, A, B, and C. Each load type contains a 'TABLE_MAP' that specifies the mapping between attribute name and column name in the table.

``````

In [13]:

TABLE_MAP = {} # map from load parameter names to column names in table

def fefs(self):
"""Return the complete set of 6 fixed end forces produced by the load."""
raise NotImplementedError()

def shear(self,x):
"""Return the shear force that is in equilibrium with that
produced by the portion of the load to the left of the point at
distance 'x'.  'x' may be a scalar or a 1-dimensional array
of values."""
raise NotImplementedError()

def moment(self,x):
"""Return the bending moment that is in equilibrium with that
produced by the portion of the load to the left of the point at
distance 'x'.  'x' may be a scalar or a 1-dimensional array
of values."""
raise NotImplementedError()

``````
``````

In [14]:

@extend

@property
def vpts(self):
"""Return a descriptor of the points at which the shear force must
be evaluated in order to draw a proper shear force diagram for this
load.  The descriptor is a 3-tuple of the form: (l,r,d) where 'l'
is the leftmost point, 'r' is the rightmost point and 'd' is the
degree of the curve between.  One of 'r', 'l' may be None."""
raise NotImplementedError()

@property
def mpts(self):
"""Return a descriptor of the points at which the moment must be
evaluated in order to draw a proper bending moment diagram for this
load.  The descriptor is a 3-tuple of the form: (l,r,d) where 'l'
is the leftmost point, 'r' is the rightmost point and 'd' is the
degree of the curve between.  One of 'r', 'l' may be None."""
raise NotImplementedError()

``````

Load type PL represents a single concentrated force, of magnitude P, at a distance a from the j-end:

``````

In [15]:

TABLE_MAP = {'P':'W1','a':'A'}

def __init__(self,L,P,a):
self.L = L
self.P = P
self.a = a

def fefs(self):
P = self.P
L = self.L
a = self.a
b = L-a
m2 = -P*a*b*b/(L*L)
m5 = P*a*a*b/(L*L)
v1 = (m2 + m5 - P*b)/L
v4 = -(m2 + m5 + P*a)/L
return EF(0.,v1,m2,0.,v4,m5)

def shear(self,x):
return -self.P*(x>self.a)

def moment(self,x):
return self.P*(x-self.a)*(x>self.a)

def __repr__(self):
return '{}(L={},P={},a={})'.format(self.__class__.__name__,self.L,self.P,self.a)

``````
``````

In [16]:

##test:
p = PL(1000.,300.,400.)
p, p.fefs()

``````
``````

Out[16]:

(PL(L=1000.0,P=300.0,a=400.0), EF(0.0,-194.4,-43200.0,0.0,-105.6,28800.0))

``````
``````

In [17]:

@extend

EPSILON = 1.0E-6

@extend
class PL:

@property
def vpts(self):
return (self.a-self.EPSILON,self.a+self.EPSILON,0)

@property
def mpts(self):
return (self.a,None,1)

``````
``````

In [18]:

##test:
p = PL(1000.,300.,400.)
p.vpts

``````
``````

Out[18]:

(399.999999, 400.000001, 0)

``````
``````

In [19]:

##test:
p.mpts

``````
``````

Out[19]:

(400.0, None, 1)

``````

Load type PLA represents a single concentrated force applied parallel to the length of the segment (producing only axial forces).

``````

In [20]:

TABLE_MAP = {'P':'W1','a':'A'}

def __init__(self,L,P,a):
self.L = L
self.P = P
self.a = a

def fefs(self):
P = self.P
L = self.L
a = self.a
c0 = -P*(L-a)/L
c3 = -P*a/L
return EF(c0=c0,c3=c3)

def shear(self,x):
return 0.

def moment(self,x):
return 0.

def __repr__(self):
return '{}(L={},P={},a={})'.format(self.__class__.__name__,self.L,self.P,self.a)

``````
``````

In [21]:

##test:
p = PLA(10.,P=100.,a=4.)
p.fefs()

``````
``````

Out[21]:

EF(-60.0,0.0,0.0,-40.0,0.0,0.0)

``````
``````

In [22]:

@extend
class PLA:

@property
def vpts(self):
return (0.,self.L,0)

@property
def mpts(self):
return (0.,self.L,0)

``````

Load type UDL represents a uniformly distributed load, of magnitude w, over the complete length of the element.

``````

In [23]:

TABLE_MAP = {'w':'W1'}

def __init__(self,L,w):
self.L = L
self.w = w

def __repr__(self):
return '{}(L={},w={})'.format(self.__class__.__name__,self.L,self.w)

def fefs(self):
L = self.L
w = self.w
return EF(0.,-w*L/2., -w*L*L/12., 0., -w*L/2., w*L*L/12.)

def shear(self,x):
l = x*(x>0.)*(x<=self.L) + self.L*(x>self.L)    # length of loaded portion
return -(l*self.w)

def moment(self,x):
l = x*(x>0.)*(x<=self.L) + self.L*(x>self.L)   # length of loaded portion
d = (x-self.L)*(x>self.L)   # distance from loaded portion to x: 0 if x <= L else x-L
return self.w*l*(l/2.+d)

@property
def vpts(self):
return (0.,self.L,1)

@property
def mpts(self):
return (0.,self.L,2)

``````
``````

In [24]:

##test:
w = UDL(12,10)
w,w.fefs()

``````
``````

Out[24]:

(UDL(L=12,w=10), EF(0.0,-60.0,-120.0,0.0,-60.0,120.0))

``````

Load type LVL represents a linearly varying distributed load actiong over a portion of the span:

``````

In [25]:

TABLE_MAP = {'w1':'W1','w2':'W2','a':'A','b':'B','c':'C'}

def __init__(self,L,w1,w2=None,a=None,b=None,c=None):
if a is not None and b is not None and c is not None and L != (a+b+c):
raise Exception('Cannot specify all of a, b & c')
if a is None:
if b is not None and c is not None:
a = L - (b+c)
else:
a = 0.
if c is None:
if b is not None:
c = L - (a+b)
else:
c = 0.
if b is None:
b = L - (a+c)
if w2 is None:
w2 = w1
self.L = L
self.w1 = w1
self.w2 = w2
self.a = a
self.b = b
self.c = c

def fefs(self):
"""This mess was generated via sympy.  See:
../../examples/cive3203-notebooks/FEM-2-Partial-lvl.ipynb """
L = float(self.L)
a = self.a
b = self.b
c = self.c
w1 = self.w1
w2 = self.w2
m2 = -b*(15*a*b**2*w1 + 5*a*b**2*w2 + 40*a*b*c*w1 + 20*a*b*c*w2 + 30*a*c**2*w1 + 30*a*c**2*w2 + 3*b**3*w1 + 2*b**3*w2 + 10*b**2*c*w1 + 10*b**2*c*w2 + 10*b*c**2*w1 + 20*b*c**2*w2)/(60.*(a + b + c)**2)
m5 = b*(20*a**2*b*w1 + 10*a**2*b*w2 + 30*a**2*c*w1 + 30*a**2*c*w2 + 10*a*b**2*w1 + 10*a*b**2*w2 + 20*a*b*c*w1 + 40*a*b*c*w2 + 2*b**3*w1 + 3*b**3*w2 + 5*b**2*c*w1 + 15*b**2*c*w2)/(60.*(a + b + c)**2)
v4 = -(b*w1*(a + b/2.) + b*(a + 2*b/3.)*(-w1 + w2)/2. + m2 + m5)/L
v1 = -b*(w1 + w2)/2. - v4
return EF(0.,v1,m2,0.,v4,m5)

def __repr__(self):
return '{}(L={},w1={},w2={},a={},b={},c={})'\
.format(self.__class__.__name__,self.L,self.w1,self.w2,self.a,self.b,self.c)

def shear(self,x):
c = (x>self.a+self.b)    # 1 if x > A+B else 0
l = (x-self.a)*(x>self.a)*(1.-c) + self.b*c  # length of load portion to the left of x
return -(self.w1 + (self.w2-self.w1)*(l/self.b)/2.)*l

def moment(self,x):
c = (x>self.a+self.b)    # 1 if x > A+B else 0
#           note: ~c doesn't work if x is scalar, thus we use 1-c
l = (x-self.a)*(x>self.a)*(1.-c) + self.b*c  # length of load portion to the left of x
d = (x-(self.a+self.b))*c                  # distance from right end of load portion to x
return ((self.w1*(d+l/2.)) + (self.w2-self.w1)*(l/self.b)*(d+l/3.)/2.)*l

@property
def vpts(self):
return (self.a,self.a+self.b,1 if self.w1==self.w2 else 2)

@property
def mpts(self):
return (self.a,self.a+self.b,2 if self.w1==self.w2 else 3)

``````

Load type CM represents a single concentrated moment of magnitude M a distance a from the j-end:

``````

In [26]:

TABLE_MAP = {'M':'W1','a':'A'}

def __init__(self,L,M,a):
self.L = L
self.M = M
self.a = a

def fefs(self):
L = float(self.L)
A = self.a
B = L - A
M = self.M
m2 = B*(2.*A - B)*M/L**2
m5 = A*(2.*B - A)*M/L**2
v1 = (M + m2 + m5)/L
v4 = -v1
return EF(0,v1,m2,0,v4,m5)

def shear(self,x):
return x*0.

def moment(self,x):
return -self.M*(x>self.A)

@property
def vpts(self):
return (None,None,0)

@property
def mpts(self):
return (self.A-self.EPSILON,self.A+self.EPSILON,1)

def __repr__(self):
return '{}(L={},M={},a={})'.format(self.__class__.__name__,self.L,self.M,self.a)

``````

Finally, the function `makeMemberLoad()` will create a load object of the correct type from the data in dictionary `data`. That dictionary would normally containing the data from one row ov the input data file table.

``````

In [27]:

def all_subclasses(cls):
_all_subclasses = []
for subclass in cls.__subclasses__():
_all_subclasses.append(subclass)
_all_subclasses.extend(all_subclasses(subclass))
return _all_subclasses

if ltype is None:
ltype = data.get('TYPE',None)
if c.__name__ == ltype and hasattr(c,'TABLE_MAP'):
MAP = c.TABLE_MAP
argv = {k:data[MAP[k]] for k in MAP.keys()}
return c(L,**argv)

``````
``````

In [28]:

##test:
ml, ml.fefs()

``````
``````

Out[28]:

(UDL(L=12,w=10), EF(0.0,-60.0,-120.0,0.0,-60.0,120.0))

``````
``````

In [29]:

ans = {'TYPE':type}
return ans

``````
``````

In [30]:

##test:

``````
``````

Out[30]:

{'TYPE': 'UDL', 'W1': 10}

``````
``````

In [ ]:

``````