Periodic_BC.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "f5672811-c30c-4474-a443-9d99a56da427",
      "metadata": {},
      "source": [
        "# Periodic BC\n",
        "\n",
        "## Manual solution\n",
        "\n",
        "First, we calculate an example by hand to have something to compare against."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1385,
      "id": "96785715-9ccf-4e27-bc19-2fdcbae4f70f",
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy import linalg as scpl"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1386,
      "id": "c699ab80-99c3-4e9a-abbd-1b6b47d9bfe8",
      "metadata": {},
      "outputs": [],
      "source": [
        "E = 2.1e9\n",
        "nu = 0.3"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1387,
      "id": "145034e9-bf48-4928-aea1-32d328d86247",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[2.30769231e+09, 6.92307692e+08, 0.00000000e+00],\n",
              "       [6.92307692e+08, 2.30769231e+09, 0.00000000e+00],\n",
              "       [0.00000000e+00, 0.00000000e+00, 1.61538462e+09]])"
            ]
          },
          "execution_count": 1387,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "#stiffness tensor (Kelvin mapping)\n",
        "C = E/(1-nu**2) * np.array([[1,nu,0],[nu,1,0],[0,0,1-nu]])\n",
        "C"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1388,
      "id": "703bf21f-dfe8-41a5-a7d8-d7e264246dd1",
      "metadata": {},
      "outputs": [],
      "source": [
        "#undeformed coordinates\n",
        "X_1 = np.array([0,0])\n",
        "X_2 = np.array([1,0])\n",
        "X_3 = np.array([1,1])\n",
        "X_4 = np.array([0,1])\n",
        "Xs = np.array([X_1,X_2,X_3,X_4])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1389,
      "id": "97747fa7-153d-46b0-8381-a929d2877da5",
      "metadata": {},
      "outputs": [],
      "source": [
        "#prescribed displacements\n",
        "u_1 = np.array([0,0])\n",
        "u_2 = np.array([0.04,-0.03])\n",
        "u_4 = np.array([0.05,0.04])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1390,
      "id": "38b38f82-3da3-445c-b3ac-252565cc9096",
      "metadata": {},
      "outputs": [],
      "source": [
        "#periodic bc\n",
        "u_3 = u_4 + u_2 - u_1"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1391,
      "id": "04f2590a-d1ff-481b-bc3f-77cff9d4171a",
      "metadata": {},
      "outputs": [],
      "source": [
        "#deformed coordinates\n",
        "us = np.array([u_1,u_2,u_3,u_4])\n",
        "xs = np.array([X_1+u_1,X_2+u_2,X_3+u_3,X_4+u_4])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1392,
      "id": "c2be2400-e9bf-4ddd-8945-4cb1ae11fb04",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAD4CAYAAAAHMeibAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOaElEQVR4nO3db2xdd33H8fcn9ipxBQNEDIKksbupDDpEJ2I6hvanDG0k3aQKiQct1hAVUpR2RTxst0oMCVkbDyaxqYUQVVE12SJPqFiZSqtpE/9UOupMpW1AQSG0iQlSXUCbaB50Sb57cC/tzc11fJLc42s375d0ZP9+5+dzvrLO/fjc3zm+J1WFpCvblnEXIGn8DAJJBoEkg0ASBoEkYHJcO966dWvNzMyMa/fSFenQoUMvVNXUYP/YgmBmZoalpaVx7V66IiV5bli/bw0kGQSSDAJJGASSMAgk8SoIgsVFmJmBLVu6XxcXx12RtPmM7fLhKCwuwp49cOpUt/3cc902wNzc+OqSNptNfUZwzz2vhMCvnTrV7ZfU3KYOguPHL65f0nCbOgh27Li4fknDbeogmJ+HTufcvk6n2y+puU0dBHNzsH8/TE9D0v26f78ThbpyjOqq2ZpXDZIcAP4SeL6q3jVkfYB/Am4CTgEfr6r/vrRyLt7cnC98XZlGedWsyRnBA8CuC6zfDVzbW/YAX7y4EoZbXFxkZmaGLVu2MDMzw6I3COhV5nKP8ZFeNauqNRdgBnhmlXVfAm7tax8B3rrWNnfu3FmrWVhYqE6nU8DLS6fTqYWFhVV/RtpMRnGMJ1Vw/pKs/jPAUg15PY5ijmAbcKKvvdzru2T33HMPpwai7tSpU9zjDQJ6lRjFMT7Kq2ajCIIM6Rv6sIQke5IsJVlaWVlZdYPHV7kRYLV+abMZxTE+yqtmowiCZeDqvvZ24OSwgVW1v6pmq2p2auq8T0t62Y5VIm21fmmzGcUxPsqrZqMIgoeAj6XrfcD/VNXPLmeD8/PzdAairtPpMO8NAnqVGNUxPjcHzz4LZ892v17yFbRhEwd17mTgl4GfAf9H96//J4C9wN7e+gD3AT8GngZm19pmrTFZ+OvJlOnp6UpS09PTThTqVWccxzirTBamxvTsw9nZ2fLDS6X1leRQVc0O9m/qOwsljYZBIMkgkGQQSMIgkIRBIAmDQBIGgSQMAkkYBJIwCCRhEEjCIJCEQSAJg0ASBoEkDAJJGASSMAgkYRBIwiCQhEEgCYNAEgaBJAwCSRgEkjAIJGEQSKJhECTZleRIkqNJ7h6y/vVJvpbk+0kOJ7lt9KVKasuaQZBkgu5jz3cD1wG3JrluYNhfAz+oquuBG4F/THLViGuV1JImZwQ3AEer6lhVvQQcBG4eGFPA65IEeC3wC+D0SCuV1JomQbANONHXXu719bsXeCdwEnga+FRVnR1JhZJa1yQIMqSvBtofAp4E3gb8HnBvkt88b0PJniRLSZZWVlYuslRJbWkSBMvA1X3t7XT/8ve7DXiwuo4CPwHeMbihqtpfVbNVNTs1NXWpNUsasSZB8ARwbZJrehOAtwAPDYw5DnwQIMlbgN8Bjo2yUEntmVxrQFWdTnIn8CgwARyoqsNJ9vbW7wM+CzyQ5Gm6byXuqqoXWqxb0gitGQQAVfUw8PBA376+708Cfz7a0iStF+8slGQQSDIIJGEQSMIgkIRBIAmDQBIGgSQMAkkYBJIwCCRhEEjCIJCEQSAJg0ASBoEkDAJJGASSMAgkYRBIwiCQhEEgCYNAEgaBJAwCSRgEkjAIJGEQSKJhECTZleRIkqNJ7l5lzI1JnkxyOMk3R1umpDat+TTkJBPAfcCfAcvAE0keqqof9I15A/AFYFdVHU/y5pbqldSCJmcENwBHq+pYVb0EHARuHhjzUeDBqjoOUFXPj7ZMSW1qEgTbgBN97eVeX7+3A29M8o0kh5J8bFQFSmrfmm8NgAzpqyHb2Ql8EHgN8N0kj1fVj87ZULIH2AOwY8eOi69WUiuanBEsA1f3tbcDJ4eMeaSqXqyqF4BvAdcPbqiq9lfVbFXNTk1NXWrNkkasSRA8AVyb5JokVwG3AA8NjPlX4I+STCbpAL8P/HC0pUpqy5pvDarqdJI7gUeBCeBAVR1Osre3fl9V/TDJI8BTwFng/qp6ps3CJY1Oqgbf7q+P2dnZWlpaGsu+pStVkkNVNTvY752FkgwCSQaBJAwCSRgEkjAIJGEQSMIgkIRBIAmDQBIGgSQMAkkYBJIwCCRhEEjCIJCEQSAJg0ASBoEkDAJJGASSMAgkYRBIwiCQhEEgCYNAEgaBJAwCSTQMgiS7khxJcjTJ3RcY994kZ5J8ZHQlSmrbmkGQZAK4D9gNXAfcmuS6VcZ9ju7j0yVtIk3OCG4AjlbVsap6CTgI3Dxk3CeBrwDPj7A+SeugSRBsA070tZd7fS9Lsg34MLDvQhtKsifJUpKllZWVi61VUkuaBEGG9NVA+/PAXVV15kIbqqr9VTVbVbNTU1MNS5TUtskGY5aBq/va24GTA2NmgYNJALYCNyU5XVVfHUWRktrVJAieAK5Ncg3wU+AW4KP9A6rqml9/n+QB4N8MAWnzWDMIqup0kjvpXg2YAA5U1eEke3vrLzgvIGnja3JGQFU9DDw80Dc0AKrq45dflqT15J2FkgwCSQaBJAwCSRgEkjAIJGEQSMIgkIRBIAmDQBIGgSQMAkkYBJIwCCRhEEjCIJCEQSAJg0ASBoEkDAJJGASSMAgkYRBIwiCQhEEgCYNAEgaBJAwCSTQMgiS7khxJcjTJ3UPWzyV5qrc8luT60ZcqqS1rBkGSCeA+YDdwHXBrkusGhv0E+JOqejfwWWD/qAuV1J4mZwQ3AEer6lhVvQQcBG7uH1BVj1XVL3vNx4Htoy1TUpuaBME24ERfe7nXt5pPAF8ftiLJniRLSZZWVlaaVympVU2CIEP6aujA5AN0g+CuYeuran9VzVbV7NTUVPMqJbVqssGYZeDqvvZ24OTgoCTvBu4HdlfVz0dTnqT10OSM4Ang2iTXJLkKuAV4qH9Akh3Ag8BfVdWPRl+mpDateUZQVaeT3Ak8CkwAB6rqcJK9vfX7gE8DbwK+kATgdFXNtle2pFFK1dC3+62bnZ2tpaWlsexbulIlOTTsj7R3FkoyCCQZBJIwCCRhEEjCIJCEQSAJg0ASBoEkDAJJGASSMAgkYRBIwiCQhEEgCYNAEgaBJAwCSRgEkjAIJGEQSMIgkIRBIAmDQBIGgSQMAkkYBJIwCCTRMAiS7EpyJMnRJHcPWZ8k/9xb/1SS91xuYXfc8R0mJ5dJzjI5ucwdd3zncjcpbSiLi4vMzMywZcsWZmZmWFxcHF8xVXXBhe6j0H8M/BZwFfB94LqBMTcBXwcCvA/4r7W2u3PnzlrN7bd/u+BXBdW3/Kpuv/3bq/6MtJksLCxUp9Mp4OWl0+nUwsJCq/sFlmrI63HNx6In+QPgM1X1oV77b3oB8vd9Y74EfKOqvtxrHwFurKqfrbbdCz0WfXJymTNntp/XPzGxzOnT5/dLm83MzAzPPffcef3T09M8++yzre33ch6Lvg040dde7vVd7BiS7EmylGRpZWVl1R2eOfO2i+qXNpvjx49fVH/bmgRBhvQNnkY0GUNV7a+q2aqanZqaWnWHExMnL6pf2mx27NhxUf1taxIEy8DVfe3twOArssmYxvbseRZ4caD3xV6/tPnNz8/T6XTO6et0OszPz4+noGETB3XuROAkcAy4hlcmC393YMxfcO5k4ffW2u6FJguruhOGExMnCs7UxMQJJwr1qrOwsFDT09OVpKanp8+ZKFxYqJqerkq6X0c1h8ilThYCJLkJ+DzdKwgHqmo+yd5ekOxLEuBeYBdwCritqobPBPZcaLJQupItLsKePXDq1Ct9nQ7s3w9zc5e37dUmCxsFQRsMAmm4mRkYckGB6Wm43AsKl3PVQNI6Wu3CQZsXFAwCaYNZ7cJBmxcUDAJpg5mf784J9Ot0uv1tMQikDWZurjsxOD0NSffrKCYKL2SyvU1LulRzc+2+8Ad5RiDJIJBkEEjCIJCEQSCJMd5inGQFGHIj5Xm2Ai+0XM7FsqZmrKmZ9axpuqrO+wyAsQVBU0mWht0bPU7W1Iw1NbMRavKtgSSDQNLmCIL94y5gCGtqxpqaGXtNG36OQFL7NsMZgaSWGQSSNk4QjOOxaiOoaa5Xy1NJHkty/bhr6hv33iRnknxkI9SU5MYkTyY5nOSb464pyeuTfC3J93s13dZyPQeSPJ/kmVXWr/vxfY5hn2i63gstPVZtHWp6P/DG3ve7N0JNfeP+E3gY+Mi4awLeAPwA2NFrv3kD1PS3wOd6308BvwCuarGmPwbeAzyzyvp1Pb4Hl41yRnADcLSqjlXVS8BB4OaBMTcD/1JdjwNvSPLWcdZUVY9V1S97zcfpPs+hTU1+TwCfBL4CPN9yPU1r+ijwYFUdB6iqtutqUlMBr+t9Avdr6QbB6bYKqqpv9faxmvU+vs+xUYJgZI9VW+ea+n2CbqK3ac2akmwDPgzsa7mWxjUBbwfemOQbSQ4l+dgGqOle4J10H8TzNPCpqjrbcl0Xst7H9zk2yicUjeyxaiPUeH9JPkA3CP6wxXqgWU2fB+6qqjPdP3ata1LTJLAT+CDwGuC7SR6vqh+NsaYPAU8Cfwr8NvDvSb5dVf/bUk1rWe/j+xwbJQjW/bFqI6qJJO8G7gd2V9XPW6ynaU2zwMFeCGwFbkpyuqq+OsaaloEXqupF4MUk3wKuB9oKgiY13Qb8Q3XfoB9N8hPgHcD3WqppLet9fJ9rPSckLjCR0spj1dahph3AUeD9G+X3NDD+AdqfLGzye3on8B+9sR3gGeBdY67pi8Bnet+/BfgpsLXl39UMq08WruvxPbhsiDOCqjqd5E7gUV55rNrh/seq0Z0Bv4nuC+8U3UQfd02fBt4EfKH3F/h0tfhfZA1rWldNaqqqHyZ5BHgKOAvcX1VDL6OtV03AZ4EHkjxN98V3V1W19q/ASb4M3AhsTbIM/B3wG331rOvxfV59vTSSdAXbKFcNJI2RQSDJIJBkEEjCIJCEQSAJg0AS8P/Ny6O5dFdpXAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "fig, ax = plt.subplots()\n",
        "ax.set_aspect('equal')\n",
        "for i in Xs:\n",
        "    ax.plot(i[0],i[1],marker='o',color='black',ls='')\n",
        "for i in xs:\n",
        "    ax.plot(i[0],i[1],marker='o',color='blue',ls='')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1393,
      "id": "7fad2c5f-2c5d-4d1f-954b-1a18818a489f",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.04,  0.05],\n",
              "       [-0.03,  0.04]])"
            ]
          },
          "execution_count": 1393,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "#Displacement gradient\n",
        "Grad_u = np.array([[(u_2[0]-u_1[0])/(X_2[0]-X_1[0]),(u_4[0]-u_1[0])/(X_4[1]-X_1[1])],[(u_2[1]-u_1[1])/(X_2[0]-X_1[0]),(u_4[1]-u_1[1])/(X_4[1]-X_1[1])]])\n",
        "Grad_u"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1394,
      "id": "90375b25-8dc0-4a05-87c1-fb6c1cdede1d",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1.04,  0.05],\n",
              "       [-0.03,  1.04]])"
            ]
          },
          "execution_count": 1394,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "#Deformation gradient\n",
        "F = Grad_u + np.eye(2)\n",
        "F"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1395,
      "id": "edd8359e-1132-448d-85d8-41b72da98c6c",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[0.04, 0.01],\n",
              "       [0.01, 0.04]])"
            ]
          },
          "execution_count": 1395,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "#small strain tensor\n",
        "eps = (Grad_u + Grad_u.T)/2\n",
        "eps"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1396,
      "id": "816df852-37f6-4485-8bca-d84d781c7009",
      "metadata": {},
      "outputs": [],
      "source": [
        "#polar decomposition\n",
        "R, U = scpl.polar(F)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1397,
      "id": "325d2370-5808-4e0d-a69b-450f299284ea",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([0.04      , 0.04      , 0.01414214])"
            ]
          },
          "execution_count": 1397,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "#stresses and strains as Kelvin vectors for computing analytical solution\n",
        "eps_K = np.array([eps[0,0],eps[1,1],np.sqrt(2)*eps[0,1]])\n",
        "eps_K"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1398,
      "id": "37804626-5600-4935-8d74-c4ae321a23ad",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([1.20000000e+08, 1.20000000e+08, 2.28449883e+07])"
            ]
          },
          "execution_count": 1398,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "sig_K = C.dot(eps_K)\n",
        "sig_K"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1399,
      "id": "040a2e15-b248-4757-844b-e6b167be0c5a",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[1.20000000e+08, 1.61538462e+07],\n",
              "       [1.61538462e+07, 1.20000000e+08]])"
            ]
          },
          "execution_count": 1399,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "#stress tensor of analytical solution\n",
        "sig = np.array([[sig_K[0],sig_K[2]/np.sqrt(2)],[sig_K[2]/np.sqrt(2),sig_K[1]]])\n",
        "sig"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1400,
      "id": "14d8fedc-f0ce-4781-bfe3-089ae110b0ac",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "-0.03428571428571429"
            ]
          },
          "execution_count": 1400,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "#strain in thickness direction (plane stress)\n",
        "eps_33 = -nu * (sig[0,0]+sig[1,1]) / E\n",
        "eps_33"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1401,
      "id": "01616403-4d02-4c91-839a-0ad188d42b91",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([1.23992308e+08, 1.08000000e+07])"
            ]
          },
          "execution_count": 1401,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "#Force on node 2 needed to stretch material\n",
        "l_23 = X_3+u_3 - X_2-u_2\n",
        "n_23 = l_23 / np.linalg.norm(l_23)\n",
        "n_23 = np.array([n_23[1],-n_23[0]])\n",
        "n_23\n",
        "f_2=sig.dot(n_23)*np.linalg.norm(l_23)\n",
        "f_2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1402,
      "id": "ea66901b-31c9-4279-aff1-e43f37b6282f",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([2.04000000e+07, 1.25284615e+08])"
            ]
          },
          "execution_count": 1402,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "l_43 = X_4+u_4 - X_3-u_3\n",
        "n_43 = l_43 / np.linalg.norm(l_43)\n",
        "n_43 = np.array([n_43[1],-n_43[0]])\n",
        "n_43\n",
        "f_4=sig.dot(n_43)*np.linalg.norm(l_43)\n",
        "f_4"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "835bcc15-408b-477b-922e-e19946f3bff6",
      "metadata": {},
      "source": [
        "## Calculation of corner node displacements from given deformation gradient\n",
        "\n",
        "Let's see if we recover the corner node displacements we used to calculate the analytical solution above can be reproduced assuming only $\\mathbf{F}$ is known:\n",
        "\n",
        "$$\n",
        "    [\\mathbf{u}_1\\ \\mathbf{u_2}\\ \\mathbf{u_3} ] = (\\mathbf{F} - \\mathbf{I}) [\\mathbf{X}_1\\ \\mathbf{X_2}\\ \\mathbf{X_3} ]\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1403,
      "id": "0097437d-e8d5-416a-8580-f052b35ab4a7",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Original inputs:\n",
            "[0. 0.]\n",
            "[ 0.04 -0.03]\n",
            "[0.09 0.01]\n",
            "[0.05 0.04]\n"
          ]
        }
      ],
      "source": [
        "print(\"Original inputs:\")\n",
        "for u in us:\n",
        "    print(u)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1404,
      "id": "7a553b3a-ac72-44b0-b85b-0f3b841dc815",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Calculated from F:\n",
            "[0. 0.]\n",
            "[ 0.04 -0.03]\n",
            "[0.09 0.01]\n",
            "[0.05 0.04]\n"
          ]
        }
      ],
      "source": [
        "print(\"Calculated from F:\")\n",
        "for X in Xs:\n",
        "    print((F-np.eye(2)).dot(X))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "214ca83f-426e-425b-99d0-845983062956",
      "metadata": {},
      "source": [
        "This confirms the use of the above equation."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "19b82a8b-61b5-4fba-8ba7-d1bf622a1933",
      "metadata": {},
      "source": [
        "## Averaging based on corner data\n",
        "\n",
        "The following equation is given for obtaining the average deformation gradient:\n",
        "\n",
        "$$\n",
        "    \\bar{\\mathbf{F}} = \\frac{1}{A_0} \\left[\n",
        "    (\\mathbf{x}_4 - \\mathbf{x}_1) \\otimes \\mathbf{N}_{43} \\lVert \\mathbf{X}_2 - \\mathbf{X}_1 \\rVert +\n",
        "    (\\mathbf{x}_2 - \\mathbf{x}_1) \\otimes \\mathbf{N}_{23} \\lVert \\mathbf{X}_4 - \\mathbf{X}_1 \\rVert\n",
        "    \\right]\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1405,
      "id": "8227fd58-926c-40b5-a0f8-3cc647a4ecb5",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "1.0"
            ]
          },
          "execution_count": 1405,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "N_43 = np.array([0,1])\n",
        "N_23 = np.array([1,0])\n",
        "A_0 = np.linalg.norm(X_4-X_1) * np.linalg.norm(X_2-X_1)\n",
        "A_0"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1406,
      "id": "49f3caba-f7fd-472c-a84b-0fb2163ddbde",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1.04,  0.05],\n",
              "       [-0.03,  1.04]])"
            ]
          },
          "execution_count": 1406,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "F_bar = ( np.outer(X_4+u_4 - X_1-u_1,N_43) * np.linalg.norm(X_2 - X_1) + np.outer(X_2+u_2 - X_1-u_1,N_23) * np.linalg.norm(X_4 - X_1)) / A_0\n",
        "F_bar"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1407,
      "id": "0e250c74-32e7-44e2-99c7-86dcfe50d9f4",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Difference between defined and calculated: \n",
            " [[0. 0.]\n",
            " [0. 0.]]\n"
          ]
        }
      ],
      "source": [
        "print(\"Difference between defined and calculated: \\n\", F - F_bar)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "25089c17-be6e-488e-a5a0-242c174ae3e3",
      "metadata": {},
      "source": [
        "So this equation works nicely. Now let's calculate the stresses by\n",
        "\n",
        "$$\n",
        "    \\bar{\\boldsymbol{\\sigma}} = \\frac{1}{A} \\text{sym}\\, \\left[\n",
        "    (\\mathbf{x}_4 - \\mathbf{x}_1) \\otimes \\mathbf{f}_4 + (\\mathbf{x}_2 - \\mathbf{x}_1) \\otimes \\mathbf{f}_2\n",
        "    \\right]\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1408,
      "id": "90a0c3a0-3dc9-421e-a6dd-2bbed782bcea",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "1.0832997046062554"
            ]
          },
          "execution_count": 1408,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "A = np.linalg.norm(X_4+u_4 - X_1-u_1) * np.linalg.norm(X_2+u_2 - X_1-u_1)\n",
        "A"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1409,
      "id": "e7f00df4-4dc4-4120-8cbe-4d0ae50a712d",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[1.19977878e+08, 1.61508682e+07],\n",
              "       [1.61508682e+07, 1.19977878e+08]])"
            ]
          },
          "execution_count": 1409,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "temp = np.outer(X_4+u_4 - X_1-u_1,f_4) + np.outer(X_2+u_2 - X_1-u_1,f_2)\n",
        "sig_bar = (temp + temp.T)/(2*A)\n",
        "sig_bar"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1410,
      "id": "a89d604d-81e3-476c-a4ba-090c60fc8eae",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Relative difference between defined and calculated: \n",
            " [[1.29188750e-04 1.73907933e-05]\n",
            " [1.73907933e-05 1.29188750e-04]]\n"
          ]
        }
      ],
      "source": [
        "print(\"Relative difference between defined and calculated: \\n\", (sig - sig_bar)/np.linalg.norm(sig))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c604b0c2-d01d-4ff4-876e-7d14769a607c",
      "metadata": {},
      "source": [
        "There's a very small discrepancy linked to the shear parts of the deformation, but in general it works."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8d38eebc-3bd2-40eb-be88-1b16d342526c",
      "metadata": {},
      "source": [
        "## Back-calculation of material parameters"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1411,
      "id": "e14f5818-a20c-4312-a723-510e5d55684f",
      "metadata": {},
      "outputs": [],
      "source": [
        "#hydrostatic stress\n",
        "p = np.trace(sig_bar)/3"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1412,
      "id": "d545b75c-4864-4fd9-a831-07d2087cbd0e",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.04      ,  0.01      ,  0.        ],\n",
              "       [ 0.01      ,  0.04      ,  0.        ],\n",
              "       [ 0.        ,  0.        , -0.03428571]])"
            ]
          },
          "execution_count": 1412,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "#strain from average deformation gradient\n",
        "Grad_u_bar = F_bar - np.eye(2)\n",
        "eps_bar = (Grad_u_bar + Grad_u_bar.T)/2\n",
        "eps_3D = np.append(np.append(eps_bar,[[0,0]],axis=0),[[0],[0],[0]],axis=1)\n",
        "eps_3D[2,2] = eps_33\n",
        "eps_3D"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1413,
      "id": "4a375d7f-aa62-4bc1-acd8-612cf44e33de",
      "metadata": {},
      "outputs": [],
      "source": [
        "#volume strain\n",
        "e = np.trace(eps_3D)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1414,
      "id": "3402abed-3416-4eb6-a627-52e9dbe750df",
      "metadata": {},
      "outputs": [],
      "source": [
        "#bulk modulus\n",
        "K = p/e"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1415,
      "id": "040af788-2c82-4c1b-a3ff-873165f60414",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Ratio of back-calculated and original bulk modulus: 0.9998156515640063\n"
          ]
        }
      ],
      "source": [
        "print(\"Ratio of back-calculated and original bulk modulus:\",  K/(E/(3*(1-2*nu))))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1416,
      "id": "fcaeadaa-31ce-4d28-971c-243dfd8dc9ed",
      "metadata": {},
      "outputs": [],
      "source": [
        "#shear stress (if non-zero; otherwise invariants needed)\n",
        "G = sig_bar[0,1]/(2.*eps_bar[0,1])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1417,
      "id": "da3a09a0-4f29-44fb-a714-b829f4990ef4",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Ratio of back-calculated and original shear modulus: 0.9998156515640072\n"
          ]
        }
      ],
      "source": [
        "print(\"Ratio of back-calculated and original shear modulus:\",  G/(E/(2*(1+nu))))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8d3a4c95-09f6-4799-8acc-b86f1ead68c3",
      "metadata": {},
      "source": [
        "So the material parameters can be correctly recovered."
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "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.9.13"
    },
    "widgets": {
      "application/vnd.jupyter.widget-state+json": {
        "state": {},
        "version_major": 2,
        "version_minor": 0
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}