C     File lsqcol.f, examples from "FEA for UG", 1985


      SUBROUTINE  LSQCOL (N,NSPACE,NPROP,LEMFRE,COORD,PROP,C,S)
C     **********************************************************
C                GENERATE ELEMENT SQUARE MATRIX
C     **********************************************************
      DIMENSION  COORD(N, NSPACE), C(LEMFRE), PROP(NPROP),
     1           S(LEMFRE, LEMFRE)
C     .....................................................
C     *** PROBLEM DEPENDENT STATEMENTS FOLLOW ***
C     .....................................................
C==>>   SIMPLE SPRING, AXIAL LOADS
      EK = PROP(1)
      S(1,1) = EK
      S(2,1) = -EK
      S(1,2) = S(2,1)
      S(2,2) = EK
      RETURN
      END
      SUBROUTINE  LSQCOL (N,NSPACE,NPROP,LEMFRE,COORD,PROP,C,S)
C     **********************************************************
C                GENERATE ELEMENT SQUARE MATRIX
C     **********************************************************
      DIMENSION  COORD(N, NSPACE), C(LEMFRE), PROP(NPROP),
     1           S(LEMFRE, LEMFRE)
C     .....................................................
C     *** PROBLEM DEPENDENT STATEMENTS FOLLOW ***
C     .....................................................
C==>>  BAR WITH HEAT CONDUCTION ONLY
      DIST = COORD(2,1) - COORD(1,1)
      A = PROP(2)
      IF ( A .LE. 0.0 )  A = 1.0
      EK = PROP(1)*A/DIST
      S(1,1) = EK
      S(2,1) = -EK
      S(1,2) = S(2,1)
      S(2,2) = EK
      RETURN
      END
      SUBROUTINE  LSQCOL (N,NSPACE,NPROP,LEMFRE,COORD,PROP,C,S)
C     **********************************************************
C                GENERATE ELEMENT SQUARE MATRIX
C     **********************************************************
      DIMENSION  COORD(N, NSPACE), C(LEMFRE), PROP(NPROP),
     1           S(LEMFRE, LEMFRE)
C     .....................................................
C     *** PROBLEM DEPENDENT STATEMENTS FOLLOW ***
C     .....................................................
C==>>  HEAT CONDUCTION BAR WITH SOURCE
      DIST = COORD(2,1) - COORD(1,1)
      A = PROP(2)
      IF ( A .LE. 0.0 )  A = 1.0
      EK = PROP(1)*A/DIST
      S(1,1) = EK
      S(2,1) = -EK
      S(1,2) = S(2,1)
      S(2,2) = EK
C-->  SOURCE TERM
      Q = PROP(3)
      C(1) = 0.5*Q*DIST
      C(2) = C(1)
      RETURN
      END
      SUBROUTINE  LSQCOL (N,NSPACE,NPROP,LEMFRE,COORD,PROP,C,S)
C     **********************************************************
C                GENERATE ELEMENT SQUARE MATRIX
C     **********************************************************
      DIMENSION  COORD(N, NSPACE), C(LEMFRE), PROP(NPROP),
     1           S(LEMFRE, LEMFRE)
C     .....................................................
C     *** PROBLEM DEPENDENT STATEMENTS FOLLOW ***
C     .....................................................
C==>>   1-D QUADRATIC BAR, FEAUG PROB 6.5-3, AXIAL LOADS
C     DEFINE PROPERTIES: E-MODULUS, A-AREA, X-BODY FORCE
      E  = PROP(1)
      A  = PROP(2)
      X  = PROP(3)
C     MEMBER LENGTH
      DX = COORD(2,1) - COORD(1,1)
C     STIFFNESS
      EABYL = E*A/DX/3.
      S(1,1) = 7.*EABYL
      S(2,1) = EABYL
      S(3,1) = -8.*EABYL
      S(1,2) = S(2,1)
      S(2,2) = 7.*EABYL
      S(3,2) = -8.*EABYL
      S(1,3) = S(3,1)
      S(2,3) = S(3,2)
      S(3,3) = 16.*EABYL
C     BODY FORCE
      AXLBY6 = A*X*DX/6.
      C(1) = AXLBY6
      C(2) = AXLBY6
      C(3) = AXLBY6*4.
      RETURN
      END
      SUBROUTINE  LSQCOL (N,NSPACE,NPROP,LEMFRE,COORD,PROP,C,S)
C     **********************************************************
C                GENERATE ELEMENT SQUARE MATRIX
C     **********************************************************
      DIMENSION  COORD(N, NSPACE), C(LEMFRE), PROP(NPROP),
     1           S(LEMFRE, LEMFRE)
C     .....................................................
C     *** PROBLEM DEPENDENT STATEMENTS FOLLOW ***
C     .....................................................
C==>>  BEAM BENDING, FEAUG FIG 8.3
      DIST = COORD(2,1) - COORD(1,1)
C     PROPERTIES: 1-MODULUS 2-MOMENT OF INERTIA
      EIBYL3 = PROP(1)*PROP(2)/DIST**3
      S(1,1) = 12. * EIBYL3
      S(2,1) = 6. * DIST * EIBYL3
      S(3,1) = -12. * EIBYL3
      S(4,1) = 6. * DIST * EIBYL3
      S(1,2) = S(2,1)
      S(2,2) = 4. * DIST * DIST * EIBYL3
      S(3,2) = -6. * DIST * EIBYL3
      S(4,2) = 2. * DIST * DIST * EIBYL3
      S(1,3) = S(3,1)
      S(2,3) = S(3,2)
      S(3,3) = S(1,1)
      S(4,3) = S(3,2)
      S(1,4) = S(4,1)
      S(2,4) = S(4,2)
      S(3,4) = S(4,3)
      S(4,4) = S(2,2)
C-->   LINEAR LOADING: 4,5 LINE LOAD NODE 1,2
      P1 = PROP(3)
      P2 = PROP(4)
      DBY20 = DIST/20.
      C(1) = DBY20 * ( 7.*P1 + 3.*P2 )
      C(2) = DBY20 * ( DIST*P1 + 2.*DIST*P2/3. )
      C(3) = DBY20 * ( 3.*P1 + 7.*P2 )
      C(4) = DBY20 * ( -2.*DIST*P1/3. - DIST*P2 )
      RETURN
      END
      SUBROUTINE  LSQCOL (N,NSPACE,NPROP,LEMFRE,COORD,PROP,C,S)
C     **********************************************************
C                GENERATE ELEMENT SQUARE MATRIX
C     **********************************************************
      DIMENSION  COORD(N, NSPACE), C(LEMFRE), PROP(NPROP),
     1           S(LEMFRE, LEMFRE)
C     .....................................................
C     *** PROBLEM DEPENDENT STATEMENTS FOLLOW ***
C     .....................................................
C==>>  CYLINDER WITH PRESSURE, FEAUG  SEC. 14.5
      DATA TWOPI, T / 6.28318531, 1.0 /
C     PROPERTIES: 1-MODULUS, 2-POISSON RATIO, 3&4 PRESSURES
      E = PROP(1)
      V = PROP(2)
      P1 = PROP(3)
      P2 = PROP(4)
      D11 = E * ( 1.-V )/(1. + V)/(1. - 2.*V)
      D12 = V * E/(1. + V)/(1. - 2.*V)
      D22 = D11
      R1 = COORD(1,1)
      R2 = COORD(2,1)
      RMID = 0.5 * (R1 + R2)
      DIST = R2 - R1
      C1 = 0.5 * TWOPI * D11 * T * ( R2 + R1 )/DIST
      C2 = TWOPI * D12
      C3 = TWOPI * D22 * DIST / ( 4. * RMID )
      S(1,1) = C1 - C2 + C3
      S(2,1) = C3 - C1
      S(1,2) = C3 - C1
      S(2,2) = C1 + C2 + C3
C      PRESSURE
      C(1) = TWOPI * T * R1 * P1
      C(2) = TWOPI * T * R2 * P2
      RETURN
      END