FEA Code to solve simple Beam Problems
Prof. P. Labossiere November 19, 2007

Off[General::spell1];
Off[General::spell];
deleteColumn[matrix_, column_] := 
        Map[Delete[#, column] &, matrix];
deleteRow[matrix_, row_] := 
        Transpose[Map[Delete[#, row] &, Transpose[matrix]]];

INPUT

Nodal coordinates {nodenumber, xcoord, ycoord}

MatrixForm[nodeinf = {{1, 0, 0}, {2, 1, 0}, {3, 1, 1}}]

Connectivity {elementnumber, Material/geomnumber, node1, node2}

MatrixForm[
      elcon = {{1, 1, 1, 2}, {2, 1, 2, 3}, {3, 2, 1, 3}}]

Applied loads

Fg = Transpose[{{0, 0, 0, 0, 0, 0, 2, -1, 0}}];

Fixed displacement boundary conditions (node number, dof 1=x, 2=y)

fixnode = {{1, 1}, {1, 2}, {2, 2}};
numfix = Length[fixnode];

Material Properties, (matnumber,E,Area,Momment of Inertia)

MatrixForm[
      Mat = {{1, 100,  .1,  .01}, {2, 
            100,  .1,  .01}}]

Problem Formation

Construct and print the local stiffness matrices

numnode = Length[nodeinf];
numel = Length[elcon];
Kg = 
        Table[0. , {i, 3 numnode}, {j, 
            3 numnode}];
Clear[x, y]
\(\(Do[{x1 = nodeinf[\([elcon[\([elnum, 3]\)], 2]\)]; 
	
          x2 = nodeinf[\([elcon[\([elnum, 4]\)], 2]\)]; 
	
          y1 = nodeinf[\([elcon[\([elnum, 3]\)], 3]\)]; 
	
          y2 = nodeinf[\([elcon[\([elnum, 4]\)], 
                3]\)]; 
L1 = 
            Sqrt[\((x2 - x1)\)^2 + \((y2 - 
                      y1)\)^2]; 
Angle = 
            ArcTan[x2 - x1, y2 - y1]; 
c = Cos[Angle]; 
          s = Sin[Angle]; 
matno = 
            elcon[\([elnum, 2]\)]; 
E1 = 
            Mat[\([matno, 2]\)]; 
A1 = 
            Mat[\([matno, 3]\)]; 
I1 = 
            Mat[\([matno, 4]\)]; 
T1 = 
            E1\ A1/\ L1; 
B1 = 
            E1\ I1/L1; 
Tmat = {{c, s, 0, 0, 0, 
                0}, {\(-s\), c, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0}, {0, 0, 0, c, 
                s, 0}, {0, 0, 0, \(-s\), c, 0}, {0, 0, 0, 0, 0, 
                1}}; 
KmatE = {{T1, 0, 0, \(-T1\), 0, 
                0}, {0, 12\ B1, 6\ L1\ B1, 0, \(-12\)\ B1, \ 6\ L1\ B1}, {0, 
                6  L1\ \ B1, 4\ L1^2\ B1, 0, \(-6\)\ L1\ B1, 
                2\ L1^2\ B1}, {\(-T1\), 0, 0, T1, 0, 
                0}, {0, \(-12\)\ B1, \(-6\)\ L1\ B1, 0, 
                12\ B1, \(-6\)\ L1\ B1}, {0, 6\ L1\ B1, 2\ L1^2\ B1, 
                0, \(-6\)\ L1\ B1, 4\ L1^2\ B1}}; 
Kmat = 
            Transpose[Tmat] . KmatE . Tmat; 
Do[
            Do[Do[\(Do[\(Kg[\([i, j]\)] = 
                        Kg[\([i, j]\)] + 
                          Kmat[\([i - 3\ elcon[\([elnum, k + 2]\)] + \ 3\ k, 
                              j - 3\ elcon[\([elnum, l + 2]\)] + 
                                3\ l]\)];\)
, {i, 
                      3\ elcon[\([elnum, k + 2]\)] - 2, 
                      3\ elcon[\([elnum, k + 2]\)]}];\)
					, {j, 
                  3\ elcon[\([elnum, l + 2]\)] - 2, 
                  3\ elcon[\([elnum, l + 2]\)]}], {k, 1, 2}], {l, 1, 
              2}]; 
Print[
            NumberForm[MatrixForm[Kmat], 5]]}, {elnum, 1, 
          numel}];\)\)
Print[NumberForm[MatrixForm[Kg], 5]]

Apply the fixed displacment boundary condition by deleting the appropriate rows and columns

Kred = Kg;
Do[Kred = 
        deleteColumn[
          deleteRow[Kred, 3 fixnode[[k, 1]] - 3 + fixnode[[k, 2]]], 
          3 fixnode[[k, 1]] - 3 + 
            fixnode[[k, 2]]]; 
Print[
        3 fixnode[[k, 1]] - 3 + fixnode[[k, 2]]], {k, numfix, 
        1, -1}]
NumberForm[MatrixForm[Kred], 5]

Fred = Fg;
Do[Fred = 
        deleteRow[Fred, 
          3 fixnode[[k, 1]] - 3 + fixnode[[k, 2]]], {k, numfix, 
        1, -1}]
MatrixForm[Fred]

Solution

MatrixForm[Ured = Inverse[Kred] . Fred]

Postprocessing...