fem_frame_complete_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.set_printoptions(suppress=True)\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import sympy\n",
    "from matplotlib.patches import Ellipse, Arc\n",
    "from matplotlib import cm, colors, patches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Draw_truss(Nodes,Element,TotalNodes,TotalElements,x_min,x_max,y_min,y_max,R,mr):\n",
    "    fig = plt.figure(figsize=(12,9))\n",
    "    ax = fig.add_subplot(1,1,1)\n",
    "    L1,L2 = (x_max-x_min)/6,(y_max-y_min)/6\n",
    "    L = min(L1,L2)\n",
    "    \n",
    "    for i in range(TotalElements):\n",
    "        n1,n2 = Element[i]['From_To'][0], Element[i]['From_To'][1]\n",
    "        X = [Nodes[n1]['coordinates'][0],Nodes[n2]['coordinates'][0]]\n",
    "        Y = [Nodes[n1]['coordinates'][1],Nodes[n2]['coordinates'][1]]\n",
    "        ax.plot(X,Y,linewidth=1,color='b')\n",
    "        \n",
    "        ### ELEMENT ###\n",
    "        x0,y0 = (X[0]+X[1])/2,(Y[0]+Y[1])/2\n",
    "        ax.text(x0-0.1*L1, y0+0.1*L2, \"E\"+str(i+1))\n",
    "    \n",
    "    if mr==1:\n",
    "        print(\"\\nExternal Forces\\n\")\n",
    "    if mr==2:\n",
    "        print(\"\\nReaction Forces\\n\")\n",
    "    \n",
    "    for i in range(TotalNodes):\n",
    "        x,y = Nodes[i]['coordinates'][0],Nodes[i]['coordinates'][1]\n",
    "        \n",
    "        if Nodes[i]['joint'] == 'h':\n",
    "            ax.plot(x,y,marker='^',color='orange',markersize=18)\n",
    "        elif Nodes[i]['joint'] == 'r':\n",
    "            ax.plot(x,y,marker='o',color='orange',markersize=15)\n",
    "        elif Nodes[i]['joint'] == 'f':\n",
    "            ax.plot(x,y,marker='s',color='orange',markersize=15)\n",
    "        else:\n",
    "            ax.plot(x,y,marker='s',color='k',markersize=6)\n",
    "        if mr==1:\n",
    "            Fx,Fy,Mz = round(Nodes[i]['external_force'][0],3),round(Nodes[i]['external_force'][1],3),round(Nodes[i]['external_force'][2],3)\n",
    "        if mr==2:\n",
    "            Fx,Fy,Mz = round(R[3*i,0],3),round(R[3*i+1,0],3),round(R[3*i+2,0],3)\n",
    "        dx,dy = 0.5*L1,0.5*L2\n",
    "        head_width = 0.2*L\n",
    "        x1,y1 = head_width+1.1*dx,head_width+1.1*dy\n",
    "        if Fx<0:\n",
    "            ax.arrow(x+x1, y, -dx, 0, head_width = head_width, width = 0.03*L2, ec ='green')\n",
    "            ax.text(x+x1, y+0.1*L2, \"Fx = \"+str(Fx))\n",
    "        if Fx>0:\n",
    "            ax.arrow(x-x1, y, dx, 0, head_width = head_width, width = 0.03*L2, ec ='green')\n",
    "            ax.text(x-2*x1, y+0.1*L2, \"Fx = \"+str(Fx))\n",
    "        if Fy<0:\n",
    "            ax.arrow(x, y+y1, 0, -dy, head_width = head_width, width = 0.03*L1, ec ='green')\n",
    "            ax.text(x, y+y1*1.1, \"Fy = \"+str(Fy))\n",
    "        if Fy>0:\n",
    "            ax.arrow(x, y-y1, 0, dy, head_width = head_width, width = 0.03*L1, ec ='green')\n",
    "            ax.text(x, y-y1*1.1, \"Fy = \"+str(Fy))\n",
    "        r1,r2,r3 = 0.9*L1*0.8,0.9*L2*0.8,0.45*L2*0.8\n",
    "        if Mz>0 and Fy<=0:\n",
    "            ellipse = Arc([x,y],r1,r2,0,0,270,color='green', linewidth=1)\n",
    "            ax.add_patch(ellipse)\n",
    "            ax.plot(x,y-r3,marker='>',color='green',markersize=8)\n",
    "            ax.text(x, y-(r3*1.3), \"Mz = \"+str(Mz))\n",
    "        if Mz>0 and Fy>0:\n",
    "            ellipse = Arc([x,y],r1,r2,0,180,450,color='green', linewidth=1)\n",
    "            ax.add_patch(ellipse)\n",
    "            ax.plot(x,y+r3,marker='<',color='green',markersize=8)\n",
    "            ax.text(x, y+(r3*1.3), \"Mz = \"+str(Mz))\n",
    "        if Mz<0 and Fy<=0:\n",
    "            ellipse = Arc([x,y],r1,r2,0,-90,180,color='green', linewidth=1)\n",
    "            ax.add_patch(ellipse)\n",
    "            ax.plot(x,y-r3,marker='<',color='green',markersize=8)\n",
    "            ax.text(x, y-(r3*1.3), \"Mz = \"+str(Mz))\n",
    "        if Mz<0 and Fy>0:\n",
    "            ellipse = Arc([x,y],r1,r2,0,90,360,color='green', linewidth=1)\n",
    "            ax.add_patch(ellipse)\n",
    "            ax.plot(x,y+r3,marker='>',color='green',markersize=8)\n",
    "            ax.text(x, y+(r3*1.3), \"Mz = \"+str(Mz))\n",
    "            \n",
    "        ### NODE ###\n",
    "        ax.text(x+0.2*L1, y+0.2*L2, \"N\"+str(i+1))\n",
    "        \n",
    "        ### FORCE ###\n",
    "            \n",
    "    ax.set_ylim(y_min-2*L2, y_max+2*L2)\n",
    "    ax.set_xlim(x_min-2*L1,x_max+2*L1)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enter total number of nodes: c\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "invalid literal for int() with base 10: 'c'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-26-4a30cd02b5ab>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m    212\u001b[0m     \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"\\nAnalysis ends here!\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    213\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 214\u001b[1;33m \u001b[0mmain\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m<ipython-input-26-4a30cd02b5ab>\u001b[0m in \u001b[0;36mmain\u001b[1;34m()\u001b[0m\n\u001b[0;32m      2\u001b[0m     \u001b[1;31m##### 2D FEM FOR TRUSS #####\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 4\u001b[1;33m     \u001b[0mTotalNodes\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0minput\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Enter total number of nodes: \"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      5\u001b[0m     \u001b[0mdof\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m3\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mTotalNodes\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      6\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: invalid literal for int() with base 10: 'c'"
     ]
    }
   ],
   "source": [
    "def main():\n",
    "    ##### 2D FEM FOR TRUSS #####\n",
    "\n",
    "    TotalNodes = int(input(\"Enter total number of nodes: \"))\n",
    "    dof = 3*TotalNodes\n",
    "\n",
    "    ### DEFINING NODE QUANTITIES ###\n",
    "\n",
    "    Nodes = {}\n",
    "    freeNode = []\n",
    "    F = []\n",
    "    x_min,x_max,y_min,y_max = 1e10,-1e10,1e10,-1e10\n",
    "\n",
    "    for i in range(TotalNodes):\n",
    "        print()\n",
    "        X = input(\"Enter (x,y) coordinates of node-\"+str(i+1)+\" (Example: 2,-3.7): \").split(',')\n",
    "        support = input(\"Enter support type at node-\"+str(i+1)+\" ---> write 'P/p' for Planar and 'I/i' for Inclined: \")\n",
    "        if support.upper() == 'P':\n",
    "            angle = 0\n",
    "        elif support.upper() == 'I':\n",
    "            angle = float(input(\"Enter inclination angle in degrees (Example: 45.3): \"))\n",
    "        joint = input(\"Enter constraint type at node-\"+str(i+1)+\" ---> write 'f' for Fixed, 'fr' for Free, 'h' for hinged and 'r' for Roller support: \")\n",
    "        if joint == 'fr':\n",
    "            freeNode.append(3*i)\n",
    "            freeNode.append(3*i+1)\n",
    "            freeNode.append(3*i+2)\n",
    "            Ext_F = input(\"Enter external force (Fx,Fy,Mz) at node-\"+str(i+1)+\" (Example: 1.2,-2,4.76 or 0,2.75,0): \").split(',')\n",
    "            F.append(float(Ext_F[0]))\n",
    "            F.append(float(Ext_F[1]))\n",
    "            F.append(float(Ext_F[2]))\n",
    "        elif joint == 'r':\n",
    "            freeNode.append(3*i)\n",
    "            freeNode.append(3*i+2)\n",
    "            Ext_Fx = input(\"Enter external force (Fx,Mz) at node-\"+str(i+1)+\" (Example: 1.2,-3.7): \").split(',')\n",
    "            Ext_F = [Ext_Fx[0],'0',Ext_Fx[1]]\n",
    "            F.append(float(Ext_F[0]))\n",
    "            F.append(float(Ext_F[2]))\n",
    "        elif joint == 'h':\n",
    "            freeNode.append(3*i+2)\n",
    "            Ext_Fx = input(\"Enter external force Mz at node-\"+str(i+1)+\" (Example: 0.75): \")\n",
    "            Ext_F = ['0','0',Ext_Fx]\n",
    "            F.append(float(Ext_F[2]))\n",
    "        else:\n",
    "            Ext_F = ['0','0','0']\n",
    "\n",
    "        Nodes[i] = {'coordinates':[float(X[0]),float(X[1])],\n",
    "                    'angle': angle,\n",
    "                    'joint': joint,\n",
    "                    'external_force': [float(Ext_F[0]),float(Ext_F[1]),float(Ext_F[2])]}\n",
    "\n",
    "        if x_min>Nodes[i]['coordinates'][0]:\n",
    "            x_min = Nodes[i]['coordinates'][0]\n",
    "        if y_min>Nodes[i]['coordinates'][1]:\n",
    "            y_min = Nodes[i]['coordinates'][1]\n",
    "        if x_max<Nodes[i]['coordinates'][0]:\n",
    "            x_max = Nodes[i]['coordinates'][0]\n",
    "        if y_max<Nodes[i]['coordinates'][1]:\n",
    "            y_max = Nodes[i]['coordinates'][1]\n",
    "\n",
    "    print()\n",
    "\n",
    "    ### NODE ELEMENT RELATIONSHIP TABLE ###\n",
    "\n",
    "    TotalElements = int(input(\"Enter total number of elements: \"))\n",
    "\n",
    "    Element = {}\n",
    "\n",
    "    for i in range(TotalElements):\n",
    "        print()\n",
    "        FromTo = input(\"Enter (From_node, To_node) for element-\"+str(i+1)+\" (Example: 2,3): \").split(',')\n",
    "        E = sympy.sympify(input(\"Enter Young's modulus of element-\"+str(i+1)+\" (Example: 2.27e10 if E = 2.27*10^10 or 2.7e10*(2**(0.5)) if E = 2.7*sqrt(2)*10^10): \"))\n",
    "        A = sympy.sympify(input(\"Enter area of element-\"+str(i+1)+\" (Example: 2.27e-3 if A = 2.27*10^-3 or 2.7*(2**(0.5)) if A = 2.7*sqrt(2)): \"))\n",
    "        I = sympy.sympify(input(\"Enter Izz of element-\"+str(i+1)+\" (Example: 2.27e-3 if I = 2.27*10^-3 or 2.7*(2**(0.5)) if I = 2.7*sqrt(2)): \"))\n",
    "        X,Y = np.array(Nodes[int(FromTo[0])-1]['coordinates']), np.array(Nodes[int(FromTo[1])-1]['coordinates'])\n",
    "        L,Z = np.linalg.norm(X-Y),Y-X\n",
    "        theta = np.arctan2(Z[1], Z[0])\n",
    "        s,c = math.sin(theta),math.cos(theta)\n",
    "\n",
    "        ### LOCAL STIFNESS MATRIX ###\n",
    "\n",
    "        c1 = A*c*c + (12*I*s*s/(L*L))\n",
    "        c2 = (A - (12*I/(L*L)))*c*s\n",
    "        c3 = 6*I*s/L\n",
    "        c4 = A*s*s + (12*I*c*c/(L*L))\n",
    "        c5 = 6*I*c/L\n",
    "\n",
    "        k = np.array([[c1, c2, -c3, -c1, -c2, -c3],\n",
    "                      [c2, c4, c5, -c2, -c4, c5],\n",
    "                      [-c3, c5, 4*I, c3, -c5, 2*I],\n",
    "                      [-c1, -c2, c3, c1, c2, c3],\n",
    "                      [-c2, -c4, -c5, c2, c4, -c5],\n",
    "                      [-c3, c5, 2*I, c3, -c5, 4*I]])\n",
    "        k = k*(E/L)\n",
    "        Element[i] = {'s':s,\n",
    "                      'c':c,\n",
    "                      'L':L,\n",
    "                      'A':A,\n",
    "                      'E':E,\n",
    "                      'I':I,\n",
    "                      'k':k,\n",
    "                      'From_To':[int(FromTo[0])-1,int(FromTo[1])-1]}\n",
    "\n",
    "\n",
    "    ### ASSEMBLE GLOBAL STIFFNESS MATRIX ###    \n",
    "    print()\n",
    "\n",
    "    K = np.zeros((dof,dof))\n",
    "    for i in range(TotalElements):\n",
    "        n1,n2 = Element[i]['From_To'][0], Element[i]['From_To'][1]\n",
    "        k = Element[i]['k']\n",
    "        print('LOCAL STIFFNESS MATRIX k'+str(i+1)+':\\n')\n",
    "        print(np.around(k.astype(np.double),4))\n",
    "        print()\n",
    "        map_local_to_global = {0:3*n1,\n",
    "                               1:3*n1+1,\n",
    "                               2:3*n1+2,\n",
    "                               3:3*n2,\n",
    "                               4:3*n2+1,\n",
    "                               5:3*n2+2}\n",
    "        for i in range(6):\n",
    "            for j in range(6):\n",
    "                K[map_local_to_global[i],map_local_to_global[j]] += k[i,j]\n",
    "\n",
    "    print('GLOBAL STIFFNESS MATRIX K:\\n')\n",
    "    print(np.around(K,4))\n",
    "    print()\n",
    "\n",
    "    ### TRANSFORM GLOBAL STIFFNESS MATRIX ###\n",
    "\n",
    "    T = np.zeros((dof,dof))\n",
    "    for i in range(TotalNodes):\n",
    "        alpha = (Nodes[i]['angle'])*np.pi/180\n",
    "        s,c = math.sin(alpha),math.cos(alpha)\n",
    "        t = np.array([[c,s,0],\n",
    "                      [-s,c,0],\n",
    "                      [0,0,1]])\n",
    "        T[3*i:3*i+3,3*i:3*i+3] = t\n",
    "\n",
    "    K = (T@K)@(T.T)\n",
    "    print('TRANSFORMED GLOBAL STIFFNESS MATRIX K:\\n')\n",
    "    print(np.around(K,4))\n",
    "    print()\n",
    "\n",
    "    ### SOLVING KQ = F ###\n",
    "\n",
    "    K_small = np.zeros((len(freeNode),len(freeNode)))\n",
    "\n",
    "    map_to_small = {}\n",
    "    for i in range(len(freeNode)):\n",
    "        map_to_small[freeNode[i]] = i\n",
    "\n",
    "    for i in range(dof):\n",
    "        if i not in freeNode:\n",
    "            continue\n",
    "        for j in range(dof):\n",
    "            if j not in freeNode:\n",
    "                continue\n",
    "            K_small[map_to_small[i],map_to_small[j]] = K[i,j]\n",
    "\n",
    "    F = np.array(F).reshape(len(freeNode),1)\n",
    "    q = (np.linalg.inv(K_small))@F\n",
    "\n",
    "    Q = np.zeros((dof,1))\n",
    "    for i in range(dof):\n",
    "        if i not in freeNode:\n",
    "            continue\n",
    "        Q[i,0] = q[map_to_small[i],0]\n",
    "\n",
    "    print(\"DISPLACEMENT MATRIX Q:\\n\")\n",
    "    print(np.around(Q,4))\n",
    "    print()\n",
    "\n",
    "    ### REACTION FORCES AND STRESSES ###\n",
    "\n",
    "    Z = K@Q\n",
    "    R = Z[:,:]\n",
    "    for i in range(TotalNodes):\n",
    "        if Nodes[i]['joint']=='fr':\n",
    "            R[3*i,0],R[3*i+1,0],R[3*i+2,0] = 0,0,0\n",
    "        elif Nodes[i]['joint']=='r':\n",
    "            R[3*i,0],R[3*i+2,0] = 0,0\n",
    "        elif Nodes[i]['joint']=='h':\n",
    "            R[3*i+2,0] = 0\n",
    "        else:\n",
    "            pass\n",
    "\n",
    "    print(\"REACTION FORCE R:\\n\")\n",
    "    print(np.around(R,4))\n",
    "    print()\n",
    "\n",
    "    '''\n",
    "    for i in range(TotalElements):\n",
    "        n1,n2 = Element[i]['From_To'][0], Element[i]['From_To'][1]\n",
    "        q = np.array([Q[2*n1,0],Q[2*n1+1,0],Q[2*n2,0],Q[2*n2+1,0]])\n",
    "        c,s = Element[i]['c'],Element[i]['s']\n",
    "        E,L = Element[i]['E'],Element[i]['L']\n",
    "        C = np.array([-c,-s,c,s])\n",
    "        S = np.dot(C,q)\n",
    "        S = S*(E/L)\n",
    "        if S<0:\n",
    "            print(\"ELEMENT-\"+str(i+1)+\" IS UNDER COMPRESSIVE STRESS OF\",S,\"UNITS\")\n",
    "            print()\n",
    "        elif S>0:\n",
    "            print(\"ELEMENT-\"+str(i+1)+\" IS UNDER TENSILE STRESS OF\",S,\"UNITS\")\n",
    "            print()\n",
    "        else:\n",
    "            print(\"ELEMENT-\"+str(i+1)+\" IS STRESS FREE\")\n",
    "            print()\n",
    "    '''\n",
    "    Draw_truss(Nodes,Element,TotalNodes,TotalElements,x_min,x_max,y_min,y_max,R,1)\n",
    "    Draw_truss(Nodes,Element,TotalNodes,TotalElements,x_min,x_max,y_min,y_max,R,2)\n",
    "    print(\"\\nAnalysis ends here!\")\n",
    "    \n",
    "main()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}