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]:
from __future__ import division, print_function

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]
    
    def __add__(self,other):
        """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)

Class MemberLoad

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]:
class MemberLoad(object):
    
    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
class MemberLoad:
    
    @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

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


In [15]:
class PL(MemberLoad):
    
    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
class MemberLoad:
    
    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

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


In [20]:
class PLA(MemberLoad):
    
    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 [22]:
##test:
p = PLA(10.,P=100.,a=4.)
p.fefs()


Out[22]:
EF(-60.0,0.0,0.0,-40.0,0.0,0.0)

In [23]:
@extend
class PLA:
    
    @property
    def vpts(self):
        return (0.,self.L,0)
    
    @property
    def mpts(self):
        return (0.,self.L,0)

Load Type UDL

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


In [24]:
class UDL(MemberLoad):
    
    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 [25]:
##test:
w = UDL(12,10)
w,w.fefs()


Out[25]:
(UDL(L=12,w=10), EF(0.0,-60.0,-120.0,0.0,-60.0,120.0))

Load Type LVL

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


In [26]:
class LVL(MemberLoad):
    
    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

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


In [27]:
class CM(MemberLoad):
    
    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)

makeMemberLoad() factory function

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 [28]:
def makeMemberLoad(L,data,ltype=None):
    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)
    for c in all_subclasses(MemberLoad):
        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)
    raise Exception('Invalid load type: {}'.format(ltype))

In [29]:
##test:
ml = makeMemberLoad(12,{'TYPE':'UDL', 'W1':10})
ml, ml.fefs()


Out[29]:
(UDL(L=12,w=10), EF(0.0,-60.0,-120.0,0.0,-60.0,120.0))

In [30]:
def unmakeMemberLoad(load):
    type = load.__class__.__name__
    ans = {'TYPE':type}
    for a,col in load.TABLE_MAP.items():
        ans[col] = getattr(load,a)
    return ans

In [31]:
##test:
unmakeMemberLoad(ml)


Out[31]:
{'TYPE': 'UDL', 'W1': 10}

In [ ]: