2DFEM.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "2DFEM ",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/adsk2050/btp-iitg/blob/master/2DFEM.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "D4-JsKGyiWbH",
        "colab_type": "code",
        "outputId": "6ae047e8-93b9-4c4b-d27e-44a790795a45",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 122
        }
      },
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive')\n",
        "\n",
        "# https://colab.research.google.com/drive/1Uu3jHPlzxGlKdwh_bfVM6q5we9J-jvtH"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Go to this URL in a browser: https://accounts.google.com/o/oauth2/auth?client_id=947318989803-6bn6qk8qdgf4n4g3pfee6491hc0brc4i.apps.googleusercontent.com&redirect_uri=urn%3aietf%3awg%3aoauth%3a2.0%3aoob&response_type=code&scope=email%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdocs.test%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive.photos.readonly%20https%3a%2f%2fwww.googleapis.com%2fauth%2fpeopleapi.readonly\n",
            "\n",
            "Enter your authorization code:\n",
            "··········\n",
            "Mounted at /content/drive\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZniGLmMYgZtP",
        "colab_type": "code",
        "outputId": "b99b2e5d-e705-460e-823c-379916589473",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        }
      },
      "source": [
        "from google.colab import files\n",
        "import csv\n",
        "import numpy as np\n",
        "import scipy as sp\n",
        "from scipy.stats import norm\n",
        "import matplotlib as mpl\n",
        "import matplotlib.pyplot as plt\n",
        "import math\n",
        "import pandas as pd\n",
        "import shutil, os\n",
        "np.set_printoptions(precision=1)\n",
        "\n",
        "\n",
        "class Material:\n",
        "    def __init__(self, c11, c12, c21, c22, c66):\n",
        "        \"\"\" creates D matrix \"\"\"\n",
        "        self.prop = np.array([[c11, c12,   0],\n",
        "                              [c21, c22,   0],\n",
        "                              [  0,   0, c66]])\n",
        "        \n",
        "class Node:\n",
        "    def __init__(self, x, y, gi):\n",
        "        \"\"\" gi is the global index of the node \"\"\"\n",
        "        self.x = x\n",
        "        self.y = y\n",
        "        self.gi =gi\n",
        "\n",
        "class Element:\n",
        "    \"\"\" We don't have body force, so \"\"\"\n",
        "    on_boundary=False\n",
        "    nodes = np.zeros((4,3))\n",
        "    u = np.zeros((1, 4))\n",
        "    v = np.zeros((1, 4))\n",
        "    def __init__(self, node_list):\n",
        "        \"\"\" Assuming node_list contains 4 nodes in anticlockwise order \"\"\"\n",
        "        for i in range(0, 4):\n",
        "            self.nodes[i][0]=node_list[i].x\n",
        "            self.nodes[i][1]=node_list[i].y\n",
        "            self.nodes[i][2]=node_list[i].gi\n",
        "    \n",
        "    def isOnBoundary(self):\n",
        "        self.on_boundary=True\n",
        "    \n",
        "    def sizni(self, z, n):\n",
        "        \"\"\" Approximation function for x and y in terms of z and n \"\"\"\n",
        "        si1zn = 0.25*(1-z)*(1-n)\n",
        "        si2zn = 0.25*(1+z)*(1-n)\n",
        "        si3zn = 0.25*(1+z)*(1+n)\n",
        "        si4zn = 0.25*(1-z)*(1+n)\n",
        "        self.si = np.array([si1zn, si2zn, si3zn, si4zn])\n",
        "\n",
        "    def dsiidz(self, z, n):\n",
        "        \"\"\" Derivative of approximation fn wrt z \"\"\"\n",
        "        dsi1dz = -0.25*(1-n)\n",
        "        dsi2dz = 0.25*(1-n)\n",
        "        dsi3dz = 0.25*(1+n)\n",
        "        dsi4dz = -0.25*(1+n)\n",
        "        self.dsidz = np.array([dsi1dz, dsi2dz, dsi3dz, dsi4dz])\n",
        "    \n",
        "    def dsiidn(self, z, n):\n",
        "        \"\"\" Derivative of approximation fn wrt n \"\"\"\n",
        "        dsi1dn = -0.25*(1-z)\n",
        "        dsi2dn = -0.25*(1+z)\n",
        "        dsi3dn = 0.25*(1+z)\n",
        "        dsi4dn = 0.25*(1-z)\n",
        "        self.dsidn = np.array([dsi1dn, dsi2dn, dsi3dn, dsi4dn])\n",
        "\n",
        "    def Je(self, z, n):\n",
        "        \"\"\" Jacobian \"\"\"\n",
        "        x = []\n",
        "        y = []\n",
        "        for i in self.nodes:\n",
        "            x.append(i[0])\n",
        "            y.append(i[1])\n",
        "        dxdz = 0.25*(n*(x[0]-x[1]+x[2]+x[3])+(-x[0]+x[1]+x[2]-x[3]))\n",
        "        dxdn = 0.25*(z*(x[0]-x[1]+x[2]-x[3])+(-x[0]-x[1]+x[2]+x[3]))\n",
        "        dydz = 0.25*(n*(y[0]-y[1]+y[2]+y[3])+(-y[0]+y[1]+y[2]-y[3]))\n",
        "        dydn = 0.25*(z*(y[0]-y[1]+y[2]-y[3])+(-y[0]-y[1]+y[2]+y[3]))\n",
        "        self.J =  np.array([[dxdz, dxdn],\n",
        "                            [dydz, dydn]])\n",
        "    \n",
        "    def diff(self):\n",
        "        \"\"\" for each iteration  \"\"\"\n",
        "        try:\n",
        "            Jstar = np.linalg.inv(self.J)\n",
        "        except np.linalg.LinAlgError:\n",
        "            Jstar = np.zeros((2, 2))\n",
        "            print(\"The Jacobian is singular, change element size or something\")\n",
        "        tmp = [[self.dsidz[0]],\n",
        "               [self.dsidn[0]]]\n",
        "        [[dsi1dx],[dsi1dy]] = np.dot(Jstar, tmp)\n",
        "\n",
        "        tmp = [[self.dsidz[1]],\n",
        "               [self.dsidn[1]]]\n",
        "        [[dsi2dx],[dsi2dy]] = np.dot(Jstar, tmp)\n",
        "\n",
        "        tmp = [[self.dsidz[2]],\n",
        "               [self.dsidn[2]]]\n",
        "        [[dsi3dx],[dsi3dy]] = np.dot(Jstar, tmp)\n",
        "\n",
        "        tmp = [[self.dsidz[3]],\n",
        "               [self.dsidn[3]]]\n",
        "        [[dsi4dx],[dsi4dy]] = np.dot(Jstar, tmp)\n",
        "        self.dsidx = np.array([dsi1dx, dsi2dx, dsi3dx, dsi4dx])\n",
        "        self.dsidy = np.array([dsi1dy, dsi2dy, dsi3dy, dsi4dy])\n",
        "    \n",
        "    def makeB(self):\n",
        "        \"\"\" creates B matrix for [K] = integral [B]t[D][B]det([J])dzdn which is stiffness matrix (8*8) \"\"\"\n",
        "        self.B = np.array([[self.dsidx[0],             0, self.dsidx[1],             0, self.dsidx[2],             0, self.dsidx[3],             0],\n",
        "                           [            0, self.dsidy[0],             0, self.dsidy[1],             0, self.dsidy[2],             0, self.dsidy[3]], \n",
        "                           [self.dsidy[0], self.dsidx[0], self.dsidy[1], self.dsidx[1], self.dsidy[2], self.dsidx[2], self.dsidy[3], self.dsidx[3]]])\n",
        "\n",
        "    def makeK(self, z, n, mat):\n",
        "        \"\"\"\n",
        "        Will do numerical integration with only one quadrature point\n",
        "        n=z=0 and W=4\n",
        "        \"\"\"\n",
        "        pts, wts = np.polynomial.legendre.leggauss(2)\n",
        "        _K = np.zeros((8, 8))\n",
        "        for i in range(0, len(wts)):\n",
        "            for j in range(0, len(wts)):\n",
        "                z = pts[i]\n",
        "                n = pts[j]\n",
        "                wi = wts[i]\n",
        "                wj = wts[j]\n",
        "                self.sizni(z, n)\n",
        "                self.dsiidz(z, n)\n",
        "                self.dsiidn(z, n)\n",
        "                self.Je(z, n)\n",
        "                self.diff()\n",
        "                self.makeB()\n",
        "                detJ = np.linalg.det(self.J)\n",
        "                _K += wi*wj*detJ*(np.dot(np.transpose(self.B), np.dot(mat.prop, self.B))) \n",
        "        # K = np.zeros((8, 8))\n",
        "        # j=0\n",
        "        # for i in range(0,4):\n",
        "        #     K[j]=_K[2*i]\n",
        "        #     K[j+1]=_K[2*i+1]\n",
        "        #     j+=2\n",
        "        self.K=_K\n",
        "\n",
        "    # def makeFQ(self, z, n, mat, Forces):\n",
        "    #     \"\"\"\n",
        "    #     Will do numerical integration with only one quadrature point\n",
        "    #     n=z=0 and W=4\n",
        "\n",
        "    #     Also, need to find tx and ty\n",
        "    #     \"\"\"\n",
        "    #     self.sizni(z, n)\n",
        "    #     self.dsiidz(z, n)\n",
        "    #     self.dsiidn(z, n)\n",
        "    #     self.Je(z, n)\n",
        "    #     self.diff()\n",
        "    #     self.FQ = np.zeros((8,0))\n",
        "    #     if self.on_boundary==False:\n",
        "    #         pass\n",
        "    #     else:\n",
        "    #         j=0\n",
        "    #         for i in range(0,4):\n",
        "    #             FQ[i]=tx*self.si[i]\n",
        "    #             FQ[2*i]=ty*self.si[i]   \n",
        "\n",
        "def load_data():\n",
        "    node_list=[]\n",
        "    with open('/content/drive/My Drive/Colab Notebooks/nodes.csv', 'r') as f:\n",
        "        csv_reader = csv.reader(f, delimiter=',')\n",
        "        for row in csv_reader:\n",
        "            tmp = Node(float(row[1]), float(row[2]), int(row[0]))\n",
        "            node_list.append(tmp)\n",
        "        f.close()\n",
        "\n",
        "    _CM = []\n",
        "    with open('/content/drive/My Drive/Colab Notebooks/connectivity.csv', 'r') as f:\n",
        "        csv_reader = csv.reader(f, delimiter=',')\n",
        "        for row in csv_reader:\n",
        "            _CM.append([int(ri) for ri in row])\n",
        "        f.close()\n",
        "\n",
        "    no_elements = len(_CM)\n",
        "    no_nodes = len(node_list)\n",
        "    CM = np.zeros((no_elements, 4))#8))\n",
        "    for element in _CM:\n",
        "        CM[element[0]-1] = element[1:5]\n",
        "    \n",
        "    L = np.zeros((no_elements, 8))\n",
        "    a=[]\n",
        "    j=1\n",
        "    for i in range(1, no_nodes+1):\n",
        "        a.append([j, j+1])\n",
        "        j+=2\n",
        "\n",
        "    for i in range(0, no_elements):\n",
        "        k=0\n",
        "        for j in range(0, 4):\n",
        "            L[i][k:k+2]=a[int(CM[i][j])-1]\n",
        "            k+=2\n",
        "    return node_list, CM, L\n",
        "\n",
        "def Assemble(CM, L, node_list, mat, length, breadth, height):\n",
        "    \"\"\"\n",
        "    For rectangular linear element (=> with 4 nodes per element), No. of nodes = 4*ne\n",
        "    Connectivity matrix CM\n",
        "\n",
        "    \"\"\"\n",
        "    K_list = []\n",
        "    # FQ_list = []\n",
        "    # element_list = []\n",
        "    for elt in CM: #as CM=matrix, elt=row of matrix which is itself a list\n",
        "        is_onboundary=False\n",
        "        local_node_list = []\n",
        "        for gi in elt: # \n",
        "            ni = node_list[int(gi)-1]\n",
        "            if ni.x==length or ni.y==breadth:\n",
        "                is_onboundary\n",
        "            local_node_list.append(ni)\n",
        "        element = Element(local_node_list)\n",
        "        element.makeK(0, 0, mat)\n",
        "        # element.makeFQ(0, 0, mat)\n",
        "        a = element.K*height\n",
        "        K_list.append(a)\n",
        "        print(a)\n",
        "        if is_onboundary:\n",
        "            element.isOnBoundary()\n",
        "        # element_list.append(element)\n",
        "        # FQ_list.append(element.FQ)\n",
        "        \n",
        "    no_nodes = len(node_list)\n",
        "    K = np.zeros((2*no_nodes, 2*no_nodes))\n",
        "    no_elements = len(K_list)\n",
        "    for i in range(0, no_elements):\n",
        "        Ki = K_list[i]\n",
        "        for j in range(0, 8):\n",
        "            for k in range(0, 8):\n",
        "                l=int(L[i][j])-1\n",
        "                m=int(L[i][k])-1\n",
        "                K[l][m]+=Ki[j][k]\n",
        "    return K#, element_list#, FQ\n",
        "\n",
        "def ApplyBC(K, node_list, F0, height):\n",
        "    \"\"\" Apply Dirichlet (displacement) and Neuman(force) boundary conditions \"\"\"\n",
        "    no_nodes = len(node_list)\n",
        "    F = np.zeros((2*no_nodes, 1))\n",
        "    for i in range(0, no_nodes):\n",
        "        node = node_list[i]\n",
        "        print(node.x, \" \",  node.y, \"\\n\")\n",
        "        if node.x==0:\n",
        "            K[i, :]=0\n",
        "            K[:, i]=0\n",
        "            K[i, i]=1\n",
        "            K[i+1, :]=0\n",
        "            K[:, i+1]=0\n",
        "            K[i+1, i+1]=1\n",
        "        # elif node.x==0:\n",
        "        #     K[i, :]=0\n",
        "        #     K[:, i]=0\n",
        "        #     K[i, i]=1\n",
        "        # elif node.y==0:\n",
        "        #     K[i+1, :]=0\n",
        "        #     K[:, i+1]=0\n",
        "        #     K[i+1, i+1]=1\n",
        "        if node.x==10:\n",
        "            if node.y>10:\n",
        "                F[i+1]=F0\n",
        "            elif node.y<10:\n",
        "                F[i+1]=-1*F0\n",
        "    return F, K\n",
        "\n",
        "def savePlots(node_list, U):\n",
        "    x = []\n",
        "    y = []\n",
        "    u=[]\n",
        "    v=[]\n",
        "    for i in range(0, len(node_list)):\n",
        "        node = node_list[i]\n",
        "        x.append(node.x)\n",
        "        y.append(node.y)\n",
        "        u.append(U[2*i])\n",
        "        v.append(U[2*i+1])\n",
        "    # print(u, '\\n', v)\n",
        "    plt.quiver(x,y, u, v)\n",
        "    plt.savefig('test{}.jpg'.format(i))\n",
        "    files.download('test{}.jpg'.format(i))\n",
        "    plt.clf()\n",
        "\n",
        "node_list, CM, L= load_data()\n",
        "print(L)\n",
        "print(CM)\n",
        "length = 10\n",
        "breadth = 20\n",
        "height = 5e-9\n",
        "F0 = 100\n",
        "E=210e9\n",
        "nu=0.3\n",
        "# for plane stress\n",
        "c11 = E/(1-nu*nu)\n",
        "c12 = (nu*E)/(1-nu*nu)\n",
        "c21=c12\n",
        "c22 = c11\n",
        "c66 = E/(2*(1+nu))\n",
        "mat = Material(c11, c12, c21, c22, c66)\n",
        "K = Assemble(CM, L, node_list, mat, length, breadth, height)\n",
        "df1 = pd.DataFrame(K)\n",
        "no_elements = len(K)\n",
        "F, K = ApplyBC(K, node_list, F0, height)\n",
        "U = np.linalg.solve(K, F).flatten()\n",
        "savePlots(node_list, U)\n",
        "print(U)\n",
        "\n",
        "df2 = pd.DataFrame(K)\n",
        "df1.subtract(df2, fill_value=0) \n",
        "with open('dataa.csv', 'w') as f:\n",
        "    f.write(df1.to_csv(index=False))\n",
        "f.close()\n",
        "\n",
        "files.download('dataa.csv') "
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[[ 1.  2.  3.  4.  9. 10.  7.  8.]\n",
            " [ 3.  4.  5.  6. 11. 12.  9. 10.]\n",
            " [ 7.  8.  9. 10. 15. 16. 13. 14.]\n",
            " [ 9. 10. 11. 12. 31. 32. 15. 16.]\n",
            " [13. 14. 15. 16. 21. 22. 19. 20.]\n",
            " [15. 16. 17. 18. 23. 24. 21. 22.]\n",
            " [19. 20. 21. 22. 27. 28. 25. 26.]\n",
            " [21. 22. 23. 24. 29. 30. 27. 28.]]\n",
            "[[ 1.  2.  5.  4.]\n",
            " [ 2.  3.  6.  5.]\n",
            " [ 4.  5.  8.  7.]\n",
            " [ 5.  6. 16.  8.]\n",
            " [ 7.  8. 11. 10.]\n",
            " [ 8.  9. 12. 11.]\n",
            " [10. 11. 14. 13.]\n",
            " [11. 12. 15. 14.]]\n",
            "[[ 631.4  312.5 -362.2 -139.4 -282.1 -187.5   12.8   14.4]\n",
            " [ 312.5  839.7 -110.6  -70.5 -187.5 -323.7  -14.4 -445.5]\n",
            " [-362.2 -110.6  496.8  -62.5  147.4  -14.4 -282.1  187.5]\n",
            " [-139.4  -70.5  -62.5  455.1   14.4  -60.9  187.5 -323.7]\n",
            " [-282.1 -187.5  147.4   14.4  496.8   62.5 -362.2  110.6]\n",
            " [-187.5 -323.7  -14.4  -60.9   62.5  455.1  139.4  -70.5]\n",
            " [  12.8  -14.4 -282.1  187.5 -362.2  139.4  631.4 -312.5]\n",
            " [  14.4 -445.5  187.5 -323.7  110.6  -70.5 -312.5  839.7]]\n",
            "[[1168.3  500.  -899.  -326.9 -389.4 -250.   120.2   76.9]\n",
            " [ 500.  1168.3 -298.1 -399.  -250.  -389.4   48.1 -379.8]\n",
            " [-899.  -298.1 1033.7  125.   254.8   48.1 -389.4  125. ]\n",
            " [-326.9 -399.   125.   783.7   76.9    4.8  125.  -389.4]\n",
            " [-389.4 -250.   254.8   76.9  389.4  125.  -254.8   48.1]\n",
            " [-250.  -389.4   48.1    4.8  125.   389.4   76.9   -4.8]\n",
            " [ 120.2   48.1 -389.4  125.  -254.8   76.9  524.  -250. ]\n",
            " [  76.9 -379.8  125.  -389.4   48.1   -4.8 -250.   774. ]]\n",
            "[[ 833.3  437.5 -496.8 -264.4 -349.4 -187.5   12.8   14.4]\n",
            " [ 437.5 1416.7 -235.6 -455.1 -187.5 -516.   -14.4 -445.5]\n",
            " [-496.8 -235.6  564.1   62.5  282.1  -14.4 -349.4  187.5]\n",
            " [-264.4 -455.1   62.5  647.4   14.4  323.7  187.5 -516. ]\n",
            " [-349.4 -187.5  282.1   14.4  564.1  -62.5 -496.8  235.6]\n",
            " [-187.5 -516.   -14.4  323.7  -62.5  647.4  264.4 -455.1]\n",
            " [  12.8  -14.4 -349.4  187.5 -496.8  264.4  833.3 -437.5]\n",
            " [  14.4 -445.5  187.5 -516.   235.6 -455.1 -437.5 1416.7]]\n",
            "[[ 1538.5   812.5 -1201.9  -639.4  -490.4  -312.5   153.8   139.4]\n",
            " [  812.5  2226.   -610.6 -1264.4  -312.5  -677.9   110.6  -283.7]\n",
            " [-1201.9  -610.6  1269.2   437.5   423.1   110.6  -490.4    62.5]\n",
            " [ -639.4 -1264.4   437.5  1456.7   139.4   485.6    62.5  -677.9]\n",
            " [ -490.4  -312.5   423.1   139.4   423.1    62.5  -355.8   110.6]\n",
            " [ -312.5  -677.9   110.6   485.6    62.5   485.6   139.4  -293.3]\n",
            " [  153.8   110.6  -490.4    62.5  -355.8   139.4   692.3  -312.5]\n",
            " [  139.4  -283.7    62.5  -677.9   110.6  -293.3  -312.5  1254.8]]\n",
            "[[ 1125.    562.5  -721.2  -389.4  -461.5  -187.5    57.7    14.4]\n",
            " [  562.5  2250.   -360.6 -1096.2  -187.5  -836.5   -14.4  -317.3]\n",
            " [ -721.2  -360.6   721.2   187.5   461.5   -14.4  -461.5   187.5]\n",
            " [ -389.4 -1096.2   187.5  1096.2    14.4   836.5   187.5  -836.5]\n",
            " [ -461.5  -187.5   461.5    14.4   721.2  -187.5  -721.2   360.6]\n",
            " [ -187.5  -836.5   -14.4   836.5  -187.5  1096.2   389.4 -1096.2]\n",
            " [   57.7   -14.4  -461.5   187.5  -721.2   389.4  1125.   -562.5]\n",
            " [   14.4  -317.3   187.5  -836.5   360.6 -1096.2  -562.5  2250. ]]\n",
            "[[ 2.1e+03  1.1e+03 -1.7e+03 -9.5e+02 -6.6e+02 -3.8e+02  2.5e+02  2.0e+02]\n",
            " [ 1.1e+03  3.9e+03 -9.2e+02 -2.7e+03 -3.8e+02 -1.2e+03  1.7e+02  4.8e+00]\n",
            " [-1.7e+03 -9.2e+02  1.7e+03  7.5e+02  6.6e+02  1.7e+02 -6.6e+02 -5.7e-14]\n",
            " [-9.5e+02 -2.7e+03  7.5e+02  2.7e+03  2.0e+02  1.2e+03 -5.7e-14 -1.2e+03]\n",
            " [-6.6e+02 -3.8e+02  6.6e+02  2.0e+02  5.2e+02  6.0e-14 -5.2e+02  1.7e+02]\n",
            " [-3.8e+02 -1.2e+03  1.7e+02  1.2e+03  6.0e-14  7.7e+02  2.0e+02 -7.7e+02]\n",
            " [ 2.5e+02  1.7e+02 -6.6e+02 -7.6e-14 -5.2e+02  2.0e+02  9.3e+02 -3.7e+02]\n",
            " [ 2.0e+02  4.8e+00 -7.6e-14 -1.2e+03  1.7e+02 -7.7e+02 -3.7e+02  1.9e+03]]\n",
            "[[ 1506.4   687.5 -1035.3  -514.4  -618.6  -187.5   147.4    14.4]\n",
            " [  687.5  3339.7  -485.6 -1993.6  -187.5 -1285.3   -14.4   -60.9]\n",
            " [-1035.3  -485.6   967.9   312.5   685.9   -14.4  -618.6   187.5]\n",
            " [ -514.4 -1993.6   312.5  1801.3    14.4  1477.6   187.5 -1285.3]\n",
            " [ -618.6  -187.5   685.9    14.4   967.9  -312.5 -1035.3   485.6]\n",
            " [ -187.5 -1285.3   -14.4  1477.6  -312.5  1801.3   514.4 -1993.6]\n",
            " [  147.4   -14.4  -618.6   187.5 -1035.3   514.4  1506.4  -687.5]\n",
            " [   14.4   -60.9   187.5 -1285.3   485.6 -1993.6  -687.5  3339.7]]\n",
            "[[ 2884.6  1437.5 -2413.5 -1264.4  -894.2  -437.5   423.1   264.4]\n",
            " [ 1437.5  6072.1 -1235.6 -4726.   -437.5 -1831.7   235.6   485.6]\n",
            " [-2413.5 -1235.6  2346.2  1062.5   961.5   235.6  -894.2   -62.5]\n",
            " [-1264.4 -4726.   1062.5  4533.7   264.4  2024.    -62.5 -1831.7]\n",
            " [ -894.2  -437.5   961.5   264.4   692.3   -62.5  -759.6   235.6]\n",
            " [ -437.5 -1831.7   235.6  2024.    -62.5  1254.8   264.4 -1447.1]\n",
            " [  423.1   235.6  -894.2   -62.5  -759.6   264.4  1230.8  -437.5]\n",
            " [  264.4   485.6   -62.5 -1831.7   235.6 -1447.1  -437.5  2793.3]]\n",
            "0.0   0.0 \n",
            "\n",
            "5.0   0.0 \n",
            "\n",
            "10.0   0.0 \n",
            "\n",
            "0.0   5.0 \n",
            "\n",
            "5.0   5.0 \n",
            "\n",
            "10.0   5.0 \n",
            "\n",
            "0.0   10.0 \n",
            "\n",
            "5.0   10.0 \n",
            "\n",
            "10.0   10.0 \n",
            "\n",
            "0.0   15.0 \n",
            "\n",
            "5.0   15.0 \n",
            "\n",
            "10.0   15.0 \n",
            "\n",
            "0.0   20.0 \n",
            "\n",
            "5.0   20.0 \n",
            "\n",
            "10.0   20.0 \n",
            "\n",
            "10.0   10.0 \n",
            "\n",
            "[ 0.0e+00  0.0e+00  1.8e-03 -1.0e+02  0.0e+00  7.7e-04 -1.0e+02  0.0e+00\n",
            " -6.3e-04  0.0e+00  0.0e+00  1.0e-02  1.0e+02  0.0e+00 -3.9e-03  3.2e-02\n",
            " -2.8e-03  3.2e-02 -9.8e-03  1.7e-02 -9.2e-03  1.5e-02 -7.6e-03  1.5e-02\n",
            " -9.5e-03  1.6e-02 -9.6e-03  1.7e-02 -1.2e-02  1.6e-02 -1.8e-02  1.2e-02]\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 0 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xwEV3CBu57ii",
        "colab_type": "code",
        "outputId": "aa70b85c-4e9b-41d6-fcee-9aa23ff3f57e",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 680
        }
      },
      "source": [
        "print(K)\n",
        "print(F)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[[1.0e+00 0.0e+00 0.0e+00 ... 0.0e+00 0.0e+00 0.0e+00]\n",
            " [0.0e+00 1.0e+00 0.0e+00 ... 0.0e+00 0.0e+00 0.0e+00]\n",
            " [0.0e+00 0.0e+00 1.7e+03 ... 0.0e+00 0.0e+00 0.0e+00]\n",
            " ...\n",
            " [0.0e+00 0.0e+00 0.0e+00 ... 1.3e+03 0.0e+00 0.0e+00]\n",
            " [0.0e+00 0.0e+00 0.0e+00 ... 0.0e+00 4.2e+02 6.2e+01]\n",
            " [0.0e+00 0.0e+00 0.0e+00 ... 0.0e+00 6.2e+01 4.9e+02]]\n",
            "[[  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [100.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]\n",
            " [  0.]]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pBQ2GQni-COb",
        "colab_type": "code",
        "outputId": "a40cea78-a060-48a4-c9f3-12e69c4b0594",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 170
        }
      },
      "source": [
        "import numpy as np\n",
        "\n",
        "# n1 = Node(0, 0, 1)\n",
        "# n2 = Node(2, 0, 2)\n",
        "# n3 = Node(2, 2, 3)\n",
        "# n4 = Node(0, 2, 4)\n",
        "# node_list = [n1, n2, n3, n4]\n",
        "# el1 = Element(node_list)\n",
        "\n",
        "print(\"Nodes\\n\", el1.nodes)\n",
        "z=0\n",
        "n=0\n",
        "# el1.sizni(z, n)\n",
        "# el1.dsiidz(z, n)\n",
        "# el1.dsiidn(z, n)\n",
        "# el1.Je(z, n)\n",
        "# el1.diff()\n",
        "el1.makeK(z, n, mat)\n",
        "print(\"si values:\\n\", el1.si)\n",
        "print(\"d(si)/dz\\n\", el1.dsidz)\n",
        "print(\"d(si)/dn\\n\", el1.dsidn)\n",
        "print(\"d(si)/dx\\n\", el1.dsidx)\n",
        "print(\"d(si)/dy\\n\", el1.dsidy)\n",
        "print(\"Jacobian\\n\", el1.J)\n",
        "print(\"B matrix\\n\", el1.B)\n",
        "print(\"K matrix\\n\", el1.K)\n",
        "print(np.shape(el1.K))\n",
        "\n",
        "E = 210e9\n",
        "print(E)\n",
        "a = np.array([[1, -1],\n",
        "              [-1, 1]])\n",
        "try:\n",
        "    inverse = np.linalg.inv(a)\n",
        "except np.linalg.LinAlgError:\n",
        "    inverse = np.zeros((2, 2))\n",
        "print(a)\n",
        "b = np.array([[5],\n",
        "              [6]])\n",
        "count=0\n",
        "for i in range(0, 2):\n",
        "    for j in range(0, 2):\n",
        "        a[i][j]=count\n",
        "        count+=1\n",
        "print(a)\n",
        "\n",
        "q = np.zeros((4, 4))\n",
        "print(q)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "210000000000.0\n",
            "[[ 1 -1]\n",
            " [-1  1]]\n",
            "[[0 1]\n",
            " [2 3]]\n",
            "[[0. 0. 0. 0.]\n",
            " [0. 0. 0. 0.]\n",
            " [0. 0. 0. 0.]\n",
            " [0. 0. 0. 0.]]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "klxkWbOuMiqN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def setk(self, mat):\n",
        "    \"\"\" Stiffness matrix with material properties stored in mat(6x6) \"\"\"\n",
        "    k11=np.zeros((4, 4))\n",
        "    for i in range(0, 4):\n",
        "        for j in range(0, 4):\n",
        "            k[i][j] = mat.prop[0][0]*self.dsidx[i]*self.dsidx[j]+mat.prop[2][2]*self.dsidy[i]*dsidy[j]\n",
        "    \n",
        "    k12=np.zeros((4, 4))\n",
        "    for i in range(0, 4):\n",
        "        for j in range(0, 4):\n",
        "            k[i][j] = mat.prop[0][1]*self.dsidx[i]*self.dsidy[j]+mat.prop[2][2]*self.dsidy[i]*dsidy[j]\n",
        "    \n",
        "    k21=np.zeros((4, 4))\n",
        "    for i in range(0, 4):\n",
        "        for j in range(0, 4):\n",
        "            k[i][j] = mat.prop[1][0]*self.dsidy[i]*self.dsidx[j]+mat.prop[2][2]*self.dsidx[i]*dsidy[j]\n",
        "    \n",
        "    k22=np.zeros((4, 4))\n",
        "    for i in range(0, 4):\n",
        "        for j in range(0, 4):\n",
        "            k[i][j] = mat.prop[1][1]*self.dsidy[i]*self.dsidy[j]+mat.prop[2][2]*self.dsidx[i]*dsidx[j]\n",
        "\n",
        "\n",
        "    print(k11)\n",
        "    print(k12)\n",
        "    print(k21)\n",
        "    print(k22)\n",
        "\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fqsIMDx6gcU-",
        "colab_type": "code",
        "outputId": "c322ff38-a278-4cae-b4fa-29de68f0b71a",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 306
        }
      },
      "source": [
        "import numpy as np\n",
        "# import scipy as sp\n",
        "\n",
        "a = np.polynomial.hermite.hermgauss(2)\n",
        "print(a)\n",
        "a = \"Anmol{0}\".format(i)\n",
        "print(a)\n",
        "\"\"\"\n",
        "*Node\n",
        "1,           0.,           0.\n",
        "2,  0.200000003,           0.\n",
        "3,  0.400000006,           0.\n",
        "4,  0.600000024,           0.\n",
        "5,  0.800000012,           0.\n",
        "6,           1.,           0.\n",
        "7,           0.,         0.25\n",
        "8,  0.200000003,         0.25\n",
        "9,  0.400000006,         0.25\n",
        "10,  0.600000024,         0.25\n",
        "11,  0.800000012,         0.25\n",
        "12,           1.,         0.25\n",
        "*Element, type=CPS4R\n",
        "1,  1,  2,  8,  7\n",
        "2,  2,  3,  9,  8\n",
        "3,  3,  4, 10,  9\n",
        "4,  4,  5, 11, 10\n",
        "5,  5,  6, 12, 11\n",
        "\"\"\"\n",
        "# a=[]\n",
        "# j=1\n",
        "# for i in range(1, 11):\n",
        "#     a.append([j, j+1])\n",
        "#     j+=2\n",
        "\n",
        "# print(a)\n",
        "\n",
        "a = np.array([[1, 2, 0],\n",
        "              [0, 0, 0]])\n",
        "\n",
        "b = a.flatten()\n",
        "print(b)\n",
        "# print([i.all()==0 for i in a])\n",
        "\n",
        "# print(\"Connectivity matrix\\n\", CM)\n",
        "\n",
        "\n",
        "# with open('foo.txt', 'w') as f:\n",
        "#   f.write('Hello Google Drive!')\n",
        "# !cat /content/drive/My\\ Drive/foo.txt\n",
        "print(int(True))\n",
        "\n",
        "no_elements = 12\n",
        "print(no_elements)\n",
        "for k in range(0, no_elements):\n",
        "    print(k)"
      ],
      "execution_count": 0,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "(array([-0.7,  0.7]), array([0.9, 0.9]))\n",
            "Anmol11\n",
            "[1 2 0 0 0 0]\n",
            "1\n",
            "12\n",
            "0\n",
            "1\n",
            "2\n",
            "3\n",
            "4\n",
            "5\n",
            "6\n",
            "7\n",
            "8\n",
            "9\n",
            "10\n",
            "11\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}