torsional rigidity of (-2,2)^2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9a568a9e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "656dc504078d4185839415298290215f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(value={'ngsolve_version': '6.2.2203', 'mesh_dim': 2, 'order2d': 2, 'order3d': 2, 'draw_vol': Fals…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqzklEQVR4nO3deXhU5d3/8fc3+wJJCCQhGwlh38IWIEBZpNYdUUGFCm5Yxa3aaqutfXxqtX1s+2hbt1pxF0UUUXGroqCCrGGHsAVCQkjIQkLIQvb790fG55diAgmZmTPL93VduZjMOTPz4WT4cHLmPvcRYwxKKaXcn4/VAZRSStmHFrpSSnkILXSllPIQWuhKKeUhtNCVUspD+Fn1wj169DDJyclWvbxSSrmlzZs3lxhjolpbZlmhJycnk5GRYdXLK6WUWxKRnLaW6SEXpZTyEFroSinlIbTQlVLKQ2ihK6WUh9BCV0opD6GFrpRSHkILXSmlPIRl49CVsooxhuNVdeSfOEVpVR0nqus5UV1HbUMTdQ1N1DcZ/HyEAD8fAnx9CAv2p1uIP91CA4gNDyK6axC+PmL1X0OpH9BCVx6t8GQNW3PL2Husgn3HKjhQVEleWTU19U3n/Jz+vkJcRDB9orrQP6YrA3t2ZWh8OH2iQhHRolfW0UJXHuVwSRXrDx1n4+FSNh0u5UjpKQBEILl7KH2juzC1fxQJ3YKJiwime5cAIkICiAj2J8jflwA/H/x8hMYmQ32jobahkfJT9ZRV11NWVUd++Snyyk5xpLSarKJKVh8opr6x+SIx4cH+jOwVwYQ+3Zk6IJp+0V204JVTiVVXLEpLSzN66r/qrKYmw46j5Xyx+xgrMgs5UFQJQPfQANKSuzEmOZK05EgGxHQlOMDX7q9f39jEweJKth85wZacE2TklHKwuAqA+IhgpgyI4rwB0Uzu34NAP/u/vvI+IrLZGJPW6jItdOWODhRWsGTTET7akU/hyVp8fYSxyZFcMCSGyf2jSOlh3eGP/BOn+HpfMV/vK+K7rBKq6hoJD/bn0tRYrhwZT1pSN91zV+dMC115hMraBj7ens+SjCNszT2Bn48wbWA0Fw3tybSB0USEBFgd8QfqGpr47mAJH249yue7CzlV30hCt2CuHp3Idem96NEl0OqIys1ooSu3VlB+ipfXZPPWhlyq6hrpG92Fa9MSuXJUvFsVYlVtA5/vPsayLUdZk1VCgK8Pl4+I46aJyQyJC7c6nnITWujKLe07VsEL3x7iw21HMcBlqbHcMCGZkYkRbn/I4mBxJa9+d5ilm/M4Vd9Iekokd0/rx4Q+3d3+76YcSwtduZWsokr++vlePt9dSLC/L9eOSWT+j3qTGBlidTS7K6+uZ0lGLi+vOcyxkzWM6x3JfRcMYGzvSKujKRelha7cwrHyGv7+5X7eyThCsL8vt0xK4cYJyXQLdb1j4/ZWU9/I2xtzefbrgxRX1DKpXw/uv2AAwxMjrI6mXIwWunJplbUNPLMyi1e+y6bJGOamJ3HXeX3p7kbHx+3lVF0jb6w/zPPfHKK0qo5ZoxP49UUDiO4aZHU05SK00JVLMsbw6c5j/OHj3RSerOXKkfH88if9PfLQSkd9/5/cS2sOEejnyz0/7scNE5IJ8NPpl7ydFrpyOYdLqnh4+W6+3V/M4NgwHrtyKKN6dbM6lsvJLqni0Y8zWbm3iJSoUP505TDSU7pbHUtZSAtduYyGxib+9e0h/vHVAQJ8fbjvgv7MS0/Cz1f3PM9k5d5CHvkok5zj1Vw/PokHLhpIaKDO3OGNzlTo+o5QTnOouJJfvrOdbUdOcOmwWB6ePpiYMD023B7TBsYwPqUHf/18H6+szWbl3iL+MjOVCX17WB1NuRDdLVIO19RkeH3dYS55ajXZJVU8PWckz143Ssu8g4IDfHl4+mDevW08/r4+/PTFDfzug52cqmu0OppyEbqHrhyq6GQN9727ndUHSpg6IIo/z0zVIu+ktORIPv35JJ74Yh8vrslmY3Ypz/50FP1iulodTVlM99CVw6w9WMIlT60h43AZf7xyKK/cOEbL3E6CA3z53WWDef3msZRW1TH9mTUs2ZSLVZ+JKdegha7srqnJ8OyqLOa+uIHwYD8+vGsi141L0lPaHWBy/yg+/fkkRvXqxgPv7eSet7dRUVNvdSxlET3kouzqRHUdv1iyjVX7ipk+PI7HrxqmozEcLDosiDfmj+O5VVn87cv9ZBacZOH1afTuEWp1NOVkZ91DF5FEEVklIpkisltE7mllnRkiskNEtolIhoj8yDFxlSvbX1jB9GfWsCarhEdnDOGp2SO0zJ3E10e4+8f9WHTLOI5X1jLjmTV8u7/Y6ljKydpzyKUBuM8YMxhIB+4UkcGnrfMVMNwYMwK4GXjRrimVy1u1r4irnltLTX0T79w2nnnjk/UQiwUm9OnB8rt+RFxEMDe+spEXVx/S4+pe5KyFbowpMMZssd2uAPYA8aetU2n+/7smFNB3kJcwxvDKd9nMf3UTvSJD+PDOiYzUMz4tlRgZwnu3T+Ang2N47JM9/GrpDuoazv2i2Mp9dOhDURFJBkYCG1pZdqWI7AU+oXkvvbXH32o7JJNRXKy/Drq7+sYmfvfBLh75KJMfD4rh3QXjiYsItjqWAkID/fjndaO59/x+LN2cx02vbtQPS71Au0/9F5EuwDfAH40xy86w3mTgYWPM+Wd6Pj31372dqmvkzre2sHJvEbdNSeGBCwfi46OHWFzRe5vzeOC9HfSL6cqrN+nQUXd3plP/27WHLiL+wHvAm2cqcwBjzLdAiojoOckeqry6nnkvbWDVviIeu2Iov7l4kJa5C5s5OoGXbhxDzvEqrnpuLVlFlVZHUg7SnlEuArwE7DHGPNnGOn1t6yEio4BA4Lg9gyrXUHiyhmtfWMf2vBM8M2cUc9OTrI6k2mFK/yiW3Dqe2oZGZj2/ls05pVZHUg7Qnj30icA8YJptWOI2EblERBaIyALbOjOBXSKyDXgWuNboR+se53BJFbOeX0tuaTWv3DiWS1NjrY6kOmBYQjjLbp9It5AA5r20kXUHdZ/L0+j0uapdDhRWMGfhBpqM4ZUbx+il0dxYUUUN1y3cQG5pNQuvT2Ny/yirI6kO6PQxdOXd9hdWMGfhekTgndvStczdXHTXIN6+NZ2UqC7c8loGX+0ptDqSshMtdHVG+wsrmPPCenxEePvWdPpG64x+nqB7l0AW/2wcA2O7smDRZv69q8DqSMoOtNBVm/Yday5zP9/mMu8T1cXqSMqOIkICWHTLOFITIrjzra1a6h5AC121at+x5sMszWU+nhQtc48UFuTP6zePZURiBHcv3sqqfUVWR1KdoIWufuBwSRVzX9qAv63MddY+zxYa6MfLN45hQM+uLHhjs45+cWNa6Oo/HCuvYe5LG2hobGLR/HFa5l4iPNif128eR1L3EOa/toktuWVWR1LnQAtd/Z+yqjrmvbSBsqo6Xr1prF7SzMtEhgawaP44orsGcsPLG9l1tNzqSKqDtNAVAJW1Ddz46iZySqtZeEOaDk30UtFhQbz5s3TCgvy58ZWN5B6vtjqS6gAtdEVtQyO3vZHBrqPlPDNnJBP66DQ83iw+IpjX54+lsclwwysbOV5Za3Uk1U5a6F7OGMMDS3fwXdZx/jIzlQuG9LQ6knIBfaK68OINY8g/cYr5r2Vwqq7R6kiqHbTQvdyTK/bzwbZ87r+gPzNHJ1gdR7mQ0UndeGrOSHbkneDuxVtoaNSLZLg6LXQv9s6mIzy9Motr0xK587y+VsdRLujCIT15ZMZQvtxTxH99uFsvZ+fi9Aq+Xmr1gWJ++/5OJvXrwWNXDtXrf6o2zUtPouDEKZ77+iC9IkO4fWofqyOpNmihe6G9x05yx6It9I3uwnPXjcLfV39RU2f2qwsHcKTsFH/5fC99o7vwk8ExVkdSrdB/yV7meGUt81/NICTQl5dvHEPXIH+rIyk3ICL8dVYqw+LDufftrew9dtLqSKoVWuhepL6xiTve3EJJZS0Lr0/TCzqrDgny92Xh9Wl0CfJj/qsZlOhwRpejhe5FHv04kw3ZpTw+cxipCRFWx1FuKCYsiIXXp1FSWcvtizZT26DDGV2JFrqXWLwxl9fX5XDr5BSuHKnDE9W5S02I4H+vHs6mw2X81we7dOSLC9FC9wIZh0t5+MNdTOrXgwcuGmh1HOUBpg+P4+fT+vJORh5vbcy1Oo6y0UL3cMfKa1iwaAvxEcE8M2cUvj46PFHZx73n92fqgCgeWZ7JtiMnrI6j0EL3aPWNTdz11haq6xpYeH0a4SE6okXZj4+P8PdrRxAdFsgdizbrnC8uQAvdg/3l33vJyCnj8ZmpOhWucoiIkACenzuakqo67nl7G41NejzdSlroHurz3cdYuDqbeelJXD48zuo4yoMNjQ/nsRlDWZNVwpMr9lkdx6tpoXugnONV3P/udlITwvndZYOsjqO8wDVjEpk9JpFnVx3ky8xCq+N4LS10D1NT38gdb27BR4RnfzqKQD9fqyMpL/H7y4cwJC6M+5dup6D8lNVxvJIWuof5w8eZ7M4/yZPXDCcxMsTqOMqLBPn78vSckdQ1NHHP4m063a4FtNA9yGc7C3hrQy63TUnhx4N08iTlfClRXXjsiqFsPFzK0yuzrI7jdbTQPUT+iVM8uGwnwxPCuf+CAVbHUV7sqlEJXDUqnqdXHmD9oeNWx/EqWugeoLHJ8Isl26hvbOIfs0fqdLjKco/OGEpS91DufXsbZVV1VsfxGmf9ly8iiSKySkQyRWS3iNzTyjrXicgOEdkpImtFZLhj4qrWPP/NQTZkl/KHGUNJ7hFqdRylCA304+k5IymtquNXS7frfC9O0p5duQbgPmPMYCAduFNEBp+2TjYwxRgzDHgUeMG+MVVbtuaW8eSK/UwfHsfMUfFWx1Hq/wyND+fBiwfy5Z4ine/FSc5a6MaYAmPMFtvtCmAPEH/aOmuNMWW2b9cDOp2fE1TU1HPP29voGRbEY1foZeSU67lxQjI/6tuDxz7ew+GSKqvjeLwOHWwVkWRgJLDhDKvNBz5r4/G3ikiGiGQUFxd35KVVKx79OJO8smr+MXsE4cE6T4tyPT4+wl+vTsXfV/jlOzqU0dHaXegi0gV4D7jXGNPq9adE5DyaC/2B1pYbY14wxqQZY9KioqLOJa+y+WpPIe9k5LFgSh/SkiOtjqNUm2LDg3n0iqFsyT3B898ctDqOR2tXoYuIP81l/qYxZlkb66QCLwIzjDE6VsmByqrqeHDZTgb27Mo95/ezOo5SZzVjRDyXpcby9y8PsOtoudVxPFZ7RrkI8BKwxxjzZBvr9AKWAfOMMfvtG1Gd7r+X76asqo4nrhmup/Yrt/HYFUPp3iWAXyzZRk29XrrOEdqzhz4RmAdME5Fttq9LRGSBiCywrfMw0B14zrY8w1GBvd0nOwpYvj2fe37cjyFx4VbHUardIkIC+Mus4RwoquR/P9dZGR3B72wrGGPWAGccPmGMuQW4xV6hVOuKK2r53QfNZ4PePrWP1XGU6rAp/aO4blwvXvoum0tSYxnVq5vVkTyKnlLoJowx/Pb9nVTVNfLENcPx07NBlZt68OKBxIYF8eulO6ht0EMv9qSt4CY+3lHAisxC7r+gP32j9epDyn11DfLnT1cNI6uokqe/0gm87EkL3Q2UVdXx++W7SU0I5+aJva2Oo1SnTR0QzcxRCfzzm4M66sWOtNDdwGOf7KH8VD1/npmqh1qUx/ivywYRGRrAr5buoF5POLILbQcX9+3+Yt7b0nwC0aDYMKvjKGU3ESEBPHbFUPYUnOT5r/WEI3vQQndhVbUN/Pb9naREhXLXtL5Wx1HK7i4c0pNLU2N5auUBsooqrI7j9rTQXdgTX+wnr+wUj1+VSpC/nkCkPNMjlw8hJMCP376/S6fZ7SQtdBe17cgJXlmbzdz0XoztrXO1KM/Vo0sgv7l4IBuzS3l3c57VcdyaFroLamhs4rfLdhLTNYgHLhpodRylHO6atETSkrrxP5/uoVSvcHTOtNBd0OvrcsgsOMnD0wfTNUinxVWez8dH+NNVw6ioaeBPn+6xOo7b0kJ3MYUna3hyxX6m9I/i4qE9rY6jlNP0j+nKzyansHRznl5c+hxpobuYRz/OpK6xiUcuH6JXIFJe5+fT+pEYGcxD7+/UaQHOgRa6C1l9oJiPdxRw59S+erFn5ZWCA3x5dMZQDhZX8cI3h6yO43a00F1ETX0jD3+4m+TuIdw2JcXqOEpZZuqAaC4dFsszq7I4UlptdRy3ooXuIl749hDZJVU8esVQHXOuvN5Dlw7CR4Q/fqIfkHaEFroLOFJazbOrsrgsNZZJ/fRaq0rFRQRz17S+/Hv3MVYf0AvKt5cWugv44yd78BHhoUsHWR1FKZdxy6TeJHcP4ffLd1PXoJN3tYcWusXWZpXw793HuPO8PsSGB1sdRymXEejny8PTB3OwuIrX1h62Oo5b0EK3UENjE498lEliZDC3TNIPQpU63bSBMUwbGM3fv9xP0ckaq+O4PC10C725IZd9hRU8dMlg/SBUqTY8fNlg6hsNj3+21+ooLk8L3SJlVXU8uWI/E/t258IhMVbHUcplJfcI5WeTe7Ns61E255RaHcelaaFb5IkV+6isbeC/p+sZoUqdzZ3n9SUmLJA/fLyHpiadYrctWugWyMw/yVsbcpmXnkT/GL3gs1JnExLgx68uHMj2Iyf4aEe+1XFclha6kxljeOyTTMKD/fnF+f2tjqOU27hqZDxD48P482d7qanXeV5ao4XuZKv2FbH24HHuPb8/4SE6Na5S7eXjIzx0yWDyy2t4aU221XFckha6EzU0NvGnT/fSu0coPx3Xy+o4Srmd8X26c8HgGJ5blUVRhQ5jPJ0WuhMtyThCVlElD148EH9f3fRKnYvfXDKIusYmnvxiv9VRXI62ipNU1jbwtxX7GZscyQWDdZiiUueqd49Qrh+fzJKMI2Tmn7Q6jks5a6GLSKKIrBKRTBHZLSL3tLLOQBFZJyK1InK/Y6K6t399c5CSyjp+e+kgHaaoVCf9fFo/woP99XJ1p2nPHnoDcJ8xZjCQDtwpIoNPW6cU+Dnwv3bO5xEKyk+xcPUhLh8ex4jECKvjKOX2wkP8uXtaP9ZklehsjC2ctdCNMQXGmC222xXAHiD+tHWKjDGbgHqHpHRzT3yxn6Ym+NWFA6yOopTHmJvei4RuwTz+2V492cimQ8fQRSQZGAlscEgaD7Sn4CTvbcnjponJJEaGWB1HKY8R6OfLfRf0Z3f+ST3ZyKbdhS4iXYD3gHuNMef0SYSI3CoiGSKSUVzsHb8m/fXzfXQN9OOOqX2tjqKUx5kxPJ5BsWE88cV+nTOddha6iPjTXOZvGmOWneuLGWNeMMakGWPSoqI8/8o8mw6XsnJvEQum9tGTiJRyAB8f4dcXDSC3tJrFG3OtjmO59oxyEeAlYI8x5knHR/IMxhj+/NleorsGctOE3lbHUcpjTe0fRXpKJE99dYDK2gar41iqPXvoE4F5wDQR2Wb7ukREFojIAgAR6SkiecAvgd+JSJ6IhDkwt8tbta+IjJwyfv7jfgQH6FznSjmKiPDgxYM4XlXHi6sPWR3HUn5nW8EYswY448BpY8wxIMFeodxdU5PhL//eR1L3EK4dk2h1HKU83ojECC4e2pOF3x7iunFJRHUNtDqSJfRMUQf4aEc+e49V8Muf9NdT/JVykvsvHEBNQxPPrsqyOopltG3srK6hiSe+2M+g2DCmp8ZZHUcpr9EnqgszR8Xz1oZc8k+csjqOJbTQ7WxJxhFyS6v59YUD8PHRU/yVcqa7p/XDYHjGS/fStdDtqKa+kWdWHmBMcjemDvD8YZlKuZrEyBBmj+nFO5uOcKS02uo4TqeFbkdvbcil8GQt910wQCfgUsoid03ri6+P8I+vDlgdxem00O3kVF0jz319kPEp3UlP6W51HKW8VkxYEHPTk1i2JY+DxZVWx3EqLXQ7WbQ+h5LKWn7xE71OqFJWu31qHwL9fPnHl961l66FbgfVdQ08/81BJvXrwdjekVbHUcrr9egSyE0Tk21DiL3nIhha6Hbw+rocjlfVce/5uneulKu4dXIKXQL8vOpSdVronVRZ28C/vjnIlP5RjE7qZnUcpZRNREgA8yf15ovMQnbnl1sdxym00DvptbWHKauu12PnSrmgmyb2pmugH8+s9I5x6VronVBRU88L3x5i2sBovbScUi4oPNifmyYm89muY+w7VmF1HIfTQu+E19flUH6qnnvP72d1FKVUG27+UW9CA3x5eqXnj3jRQj9HVbUNvLj6EOcNiCI1IcLqOEqpNkSEBHDDhGQ+2VnAgULP3kvXQj9Hb23Ipay6nrum6d65Uq7ulkkpBPv7evwcL1ro56CmvpEXVh9iYt/uOrJFKTcQGRrAvPQkPtqe79Fnj2qhn4Mlm45QXFHL3bp3rpTbuGVSCgF+Ph49X7oWegfVNTTx/DcHGZPcjXF6VqhSbiOqayDXjUviw235HC6psjqOQ2ihd9CyLXkUlNdw17R+OqOiUm7mtskp+PkI//z6oNVRHEILvQMaGpt47uuDDE8IZ3K/HlbHUUp1UHRYENekJbJsax4F5Z53VSMt9A5Yvj2f3NJq3TtXyo3dOjmFJgMvrs62OordaaG3U1OT4dlVWQzs2ZXzB0VbHUcpdY4SI0O4fHgcizfmUlZVZ3Ucu9JCb6cVewo5WFzFHef11b1zpdzcgil9qK5r5LV1h62OYlda6O1gjOH5bw6SGBnMJUN7Wh1HKdVJA2y/ab+69jBVtQ1Wx7EbLfR2yMgpY2vuCX42KQU/X91kSnmC26f25UR1PYs35lodxW60ndrh+a8PEhkawNWjE62OopSyk9FJzeeSvLg6m9qGRqvj2IUW+lnsL6zgq71F3DA+meAAX6vjKKXs6I7z+nLsZA0fbD1qdRS70EI/i399c4hgf1+uH59kdRSllJ1N7teDIXFhPP/NIRqbjNVxOk0L/QwKyk/x4bajXDsmkW6hAVbHUUrZmYhw+9Q+ZJdUsSKz0Oo4nXbWQheRRBFZJSKZIrJbRO5pZR0RkadEJEtEdojIKMfEda6X12RjgPk/6m11FKWUg1w0pCeJkcEsXH3I6iid1p499AbgPmPMYCAduFNEBp+2zsVAP9vXrcA/7ZrSAuXV9by1IZfpqbEkRoZYHUcp5SB+vj7Mn9ibzTllbM4ptTpOp5y10I0xBcaYLbbbFcAeIP601WYAr5tm64EIEYm1e1onWrQhh6q6Rm6d3MfqKEopB7s6LZHwYH8Wfuve0wF06Bi6iCQDI4ENpy2KB460+D6PH5Y+InKriGSISEZxcXEHozpPTX0jr3x3mCn9oxgcF2Z1HKWUg4UG+jE3vRefZx5z66l1213oItIFeA+41xhz8lxezBjzgjEmzRiTFhUVdS5P4RTLthylpLKW26akWB1FKeUkN4xPxt/Hh5fWuO9eersKXUT8aS7zN40xy1pZ5SjQ8qybBNt9bqepybBw9SFSE8IZn9Ld6jhKKSeJDgviipFxvLv5CKVuOmlXe0a5CPASsMcY82Qbqy0HrreNdkkHyo0xBXbM6TQr9xaRXVLFrZNTdBIupbzMzyalUFPfxKL1OVZHOSft2UOfCMwDponINtvXJSKyQEQW2Nb5FDgEZAELgTscE9fxXlqTTVx4EBcN0Um4lPI2/WK6Mm1gNK+tPUxNvftNB+B3thWMMWuAM+6qGmMMcKe9QlklM/8k6w4d5zcXD9RJuJTyUj+blMKchet5f+tR5oztZXWcDtHWauGV77IJ9vdl9hj3+iEqpewnPSWSofFhvLj6EM37qu5DC92mpLKWD7flM2t0AuEh/lbHUUpZRES4eWJvDhZXsfpAidVxOkQL3ebN9bnUNTZx48Rkq6MopSx2aWosPboE8sp37jWEUQsdqG1o5I31OZw3IIo+UV2sjqOUsligny9z03uxal8xh4orrY7TblrowEfbCyiprGX+j/REIqVUs5+O64W/r/Da2sNWR2k3ry90Ywwvr8mmf0wXJvbVE4mUUs2iuwYxPTWOpZvzOFlTb3WcdvH6Qt+QXUpmwUlunthbTyRSSv2Hmyb2pqqukXcz8qyO0i5eX+gvr8mmW4g/V4z8wVxiSikvNywhnLSkbry29rBbXNHIqws953gVK/YUct24JIL89XqhSqkfumlib3JLq1m5t8jqKGfl1YX+6trD+PkI8/R6oUqpNlwwJIbY8CC3GMLotYVeUVPPuxl5XJYaR0xYkNVxlFIuyt/Xh3njk1h78Dh7j53TzOFO47WF/v7Wo1TWNnDDhGSroyilXNycMb0I8vdx+SGMXlnoxhheX5dDakI4IxIjrI6jlHJx3UIDuHJkPMu2HKXMhedK98pCX3+olKyiSual67FzpVT73DihN7UNTSzelGt1lDZ5ZaG/sf4wESH+TB8eZ3UUpZSbGNCzKxP6dGfRuhwaGpusjtMqryv0wpM1fL67kGvSEnWoolKqQ64fn0R+eQ2r9rnmRe69rtDf2pBLkzFcN07nPFdKdcz5g2KICQt02UvUeVWh1zc2sXhjLlP6R5HUPdTqOEopN+Pn68PsMb349kAxOcerrI7zA15V6F/sLqSoopbr9UQipdQ5mjO2Fz4ivLXB9T4c9apCf33dYRK6BTOlf7TVUZRSbqpneBA/GRTDOxlHXO5C0l5T6PsLK9iQXcrc9CR8fXRWRaXUuZubnkRZdT2f7SqwOsp/8JpCf2NdDgF+PlyTlmh1FKWUm5vQpzu9e4SyaL1rHXbxikKvqKln2ZY8LkuNJTI0wOo4Sik35+MjXDeuF5tzysjMd535Xbyi0D/YepSqukauH59sdRSllIeYNTqBQD8fFm1wnSGMHl/o38/bMiw+nOEJ4VbHUUp5iIiQAKYPj+ODrUepcJFL1Hl8oW/MLuWAbd4WvcScUsqe5qYnUV3XyAdbj1odBfCCQn970xG6BvnpvC1KKbsbnhDOsPhw3lifgzHWX6LOowv9RHUdn+ws4IoR8QQH6LwtSin7EhHmpvdif2Elmw6XWR3Hswv9/a1HqWtoYs5YnbdFKeUY04fH0TXIzyXmdzlroYvIyyJSJCK72ljeTUTeF5EdIrJRRIbaP2bHGWNYvDGX4QnhDI4LszqOUspDhQT4MXNUAv/edYxSiy9+0Z499FeBi86w/LfANmNMKnA98A875Oq0Lbkn2F9YyWzdO1dKOdjssYnUNTbxvsUfjp610I0x3wKlZ1hlMLDStu5eIFlEYuwT79wt3phLaICvfhiqlHK4gT3DGJEYwdsbcy39cNQex9C3A1cBiMhYIAlIaG1FEblVRDJEJKO42HETxJ+sqefjHflcPiKOLoF+DnsdpZT63uwxiRwoqmRL7gnLMtij0B8HIkRkG3A3sBVodQoyY8wLxpg0Y0xaVFSUHV66dR9uPUpNvX4YqpRynsuGxxES4MsSC6852ulCN8acNMbcZIwZQfMx9CjgUGeftxN5eGvjEQbHhjEsXs8MVUo5R5dAPy4fHsdH2wssO3O004UuIhEi8v2MV7cA3xpjLJutZkdeOXsKTjJnXC89M1Qp5VTXjknkVH0jH223Zlrd9gxbXAysAwaISJ6IzBeRBSKywLbKIGCXiOwDLgbucVzcs3t7Uy7B/r7MGKEfhiqlnGtEYgQDYrpadtjlrJ8YGmPmnGX5OqC/3RJ1QmVtAx9uy+ey1FjCgvytjqOU8jIiwrVjEvnDx5lk5p90+jkwHnWm6Efb86mua9Sx50opy1w1Kp4APx/eyTji9Nf2qEJ/e2MuA2K6MqpXhNVRlFJeKiIkgIuG9GTZljynX3PUYwp9d3452/PKmT02UT8MVUpZavaYRE7WNPDvXcec+roeU+jvZuQR4OvDlSPjrY6ilPJy6Snd6RUZwttO/nDUIwq9rqGJD7cd5SeDY4gI0WuGKqWs5ePT/OHo+kOlZJdUOe91nfZKDrRybyFl1fXMSmt1xgGllHK6WaMT8PURlmxy3oejHlHo72bkERMWyOR+jptOQCmlOiImLIjzBkSzdHMe9Y1NTnlNty/0oooavt5fzJUjm/83VEopV3FNWgIllbV8u99xkxG25PaF/sHWozQ2Ga7Wwy1KKRdz3sBouocGsHRznlNez60L3RjDuxl5jOoVQZ+oLlbHUUqp/+Dv68OMEfF8uaeQMidczcitC31HXjkHiiqZNTrR6ihKKdWqq9MSqG80LN+e7/DXcutCf3fzEYL8fbhseKzVUZRSqlWDYsMYEhfmlMMublvoNfWNLN+Wz0VDeupEXEoplzZrdAI7j5az95hjZxZ320JfkVnIyZoGPdyilHJ5M0bE4+8rvOfgvXS3LfR3N+cRHxHMhD7drY6ilFJnFBkawLSB0by/Nd+hY9LdstALyk+x+kAxM0fF46Njz5VSbmDW6ESHj0l3y0JftuUoxsDM0Tr2XCnlHqYOiHL4mHS3K3RjDEs35zG2dyRJ3UOtjqOUUu3i7+vDFSMdOybd7Qp9c04Z2SVVXK1750opNzNrtGPHpLtdoYvA5P5RXDJMx54rpdzLoNgwZoyIIyLEMUOtxRjjkCc+m7S0NJORkWHJayullLsSkc3GmLTWlrndHrpSSqnWaaErpZSH0EJXSikPoYWulFIeQgtdKaU8hBa6Ukp5CC10pZTyEFroSinlISw7sUhEioGcc3x4D6DEjnHsxVVzgetm01wdo7k6xhNzJRljolpbYFmhd4aIZLR1ppSVXDUXuG42zdUxmqtjvC2XHnJRSikPoYWulFIewl0L/QWrA7TBVXOB62bTXB2juTrGq3K55TF0pZRSP+Sue+hKKaVOo4WulFIewmULXUSuFpHdItIkIm0O7xGRi0Rkn4hkiciDLe7vLSIbbPcvEZEAO+WKFJEVInLA9me3VtY5T0S2tfiqEZErbMteFZHsFstGOCuXbb3GFq+9vMX9Vm6vESKyzvbz3iEi17ZYZtft1db7pcXyQNvfP8u2PZJbLPuN7f59InJhZ3KcQ65fikimbft8JSJJLZa1+jN1Uq4bRaS4xevf0mLZDbaf+wERucHJuf7WItN+ETnRYpkjt9fLIlIkIrvaWC4i8pQt9w4RGdViWee3lzHGJb+AQcAA4GsgrY11fIGDQAoQAGwHBtuWvQPMtt1+HrjdTrn+Ajxou/0g8OezrB8JlAIhtu9fBWY5YHu1KxdQ2cb9lm0voD/Qz3Y7DigAIuy9vc70fmmxzh3A87bbs4ElttuDbesHAr1tz+PrxFzntXgP3f59rjP9TJ2U60bgmVYeGwkcsv3ZzXa7m7Nynbb+3cDLjt5etueeDIwCdrWx/BLgM0CAdGCDPbeXy+6hG2P2GGP2nWW1sUCWMeaQMaYOeBuYISICTAOW2tZ7DbjCTtFm2J6vvc87C/jMGFNtp9dvS0dz/R+rt5cxZr8x5oDtdj5QBLR6Jlwntfp+OUPepcCPbdtnBvC2MabWGJMNZNmezym5jDGrWryH1gPOuEp6e7ZXWy4EVhhjSo0xZcAK4CKLcs0BFtvptc/IGPMtzTtwbZkBvG6arQciRCQWO20vly30dooHjrT4Ps92X3fghDGm4bT77SHGGFNgu30MiDnL+rP54Zvpj7Zft/4mIoFOzhUkIhkisv77w0C40PYSkbE073UdbHG3vbZXW++XVtexbY9ymrdPex7ryFwtzad5L+97rf1MnZlrpu3ns1REEjv4WEfmwnZoqjewssXdjtpe7dFWdrtsL79OReskEfkS6NnKooeMMR86O8/3zpSr5TfGGCMibY77tP3POwz4vMXdv6G52AJoHov6APAHJ+ZKMsYcFZEUYKWI7KS5tM6ZnbfXG8ANxpgm293nvL08kYjMBdKAKS3u/sHP1BhzsPVnsLuPgMXGmFoRuY3m326mOem122M2sNQY09jiPiu3l0NZWujGmPM7+RRHgcQW3yfY7jtO868yfra9rO/v73QuESkUkVhjTIGtgIrO8FTXAO8bY+pbPPf3e6u1IvIKcL8zcxljjtr+PCQiXwMjgfeweHuJSBjwCc3/ma9v8dznvL1a0db7pbV18kTEDwin+f3Unsc6Mhcicj7N/0lOMcbUfn9/Gz9TexTUWXMZY463+PZFmj8z+f6xU0977Nd2yNSuXC3MBu5seYcDt1d7tJXdLtvL3Q+5bAL6SfMIjQCaf3jLTfOnDKtoPn4NcANgrz3+5bbna8/z/uDYna3Uvj9ufQXQ6qfhjsglIt2+P2QhIj2AiUCm1dvL9rN7n+Zji0tPW2bP7dXq++UMeWcBK23bZzkwW5pHwfQG+gEbO5GlQ7lEZCTwL+ByY0xRi/tb/Zk6MVdsi28vB/bYbn8OXGDL1w24gP/8TdWhuWzZBtL8AeO6Fvc5cnu1x3Lgettol3Sg3LbTYp/t5ahPezv7BVxJ83GkWqAQ+Nx2fxzwaYv1LgH20/w/7EMt7k+h+R9cFvAuEGinXN2Br4ADwJdApO3+NODFFusl0/y/rs9pj18J7KS5mBYBXZyVC5hge+3ttj/nu8L2AuYC9cC2Fl8jHLG9Wnu/0HwI53Lb7SDb3z/Ltj1SWjz2Idvj9gEX2/n9frZcX9r+HXy/fZaf7WfqpFz/A+y2vf4qYGCLx95s245ZwE3OzGX7/vfA46c9ztHbazHNo7Tqae6v+cACYIFtuQDP2nLvpMUIPntsLz31XymlPIS7H3JRSillo4WulFIeQgtdKaU8hBa6Ukp5CC10pZTyEFroSinlIbTQlVLKQ/w/uNwC5YACoKcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torsional rigidity = 35.98771580341495\n"
     ]
    }
   ],
   "source": [
    "# torsional rigidity\n",
    "# solve the Poisson equation -Delta u = f\n",
    "# with Dirichlet boundary condition u = 0\n",
    "\n",
    "from ngsolve import *\n",
    "import netgen.geom2d as geom2d\n",
    "from ngsolve.webgui import Draw\n",
    "\n",
    "ngsglobals.msg_level = 1\n",
    "\n",
    "geo = geom2d.SplineGeometry()\n",
    "p1,p2,p3,p4 = [ geo.AppendPoint(x,y) for x,y in [(-2,-2), (2,-2), (2,2), (-2,2)] ]\n",
    "geo.Append ([\"line\", p1, p2], leftdomain=1, rightdomain=0)\n",
    "geo.Append ([\"line\", p2, p3], leftdomain=1, rightdomain=0)\n",
    "geo.Append ([\"line\", p3, p4], leftdomain=1, rightdomain=0)\n",
    "geo.Append ([\"line\", p4, p1], leftdomain=1, rightdomain=0)\n",
    "        \n",
    "# generate a triangular mesh of mesh-size maxh\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.05))\n",
    "\n",
    "# H1-conforming finite element space\n",
    "fes = H1(mesh, order=3, dirichlet=[1,2,3,4])\n",
    "\n",
    "# define trial- and test-functions\n",
    "u = fes.TrialFunction()\n",
    "v = fes.TestFunction()\n",
    "\n",
    "# define the coefficient function\n",
    "b = 1\n",
    "c = 2\n",
    "\n",
    "# the right hand side\n",
    "f = LinearForm(fes)\n",
    "f +=  c *v * dx\n",
    "\n",
    "# the bilinear-form \n",
    "a = BilinearForm(fes, symmetric=True)\n",
    "a += b*grad(u)*grad(v)*dx\n",
    "\n",
    "a.Assemble()\n",
    "f.Assemble()\n",
    "\n",
    "# the solution field \n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = a.mat.Inverse(fes.FreeDofs(), inverse=\"sparsecholesky\") * f.vec\n",
    "# print (u.vec)\n",
    "\n",
    "\n",
    "# plot the solution (netgen-gui only)\n",
    "Draw (gfu)\n",
    "#Draw (-grad(gfu), mesh, \"Flux\")\n",
    "\n",
    "# 2D plot\n",
    "import matplotlib.pyplot as plt\n",
    "N = 100\n",
    "ypnts = [i/N for i in range(-N, N+1, 1)]\n",
    "vals = [gfu(mesh(0,yi)) for yi in ypnts]\n",
    "\n",
    "plt.plot(ypnts,vals)\n",
    "plt.show(block=False)\n",
    "\n",
    "print (\"torsional rigidity =\", Integrate ( grad(gfu)*grad(gfu), mesh))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "3d644e1e",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[35.9936, 35.5328, 0.9900990099009901]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "M1 = 0.1406 *  4 ** 4\n",
    "M2 = 0.1388* 4**4\n",
    "t = 1/1.01\n",
    "print([M1, M2, t])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3896ab6b",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c4be6612565444de8186b1d2dcef373e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(value={'ngsolve_version': '6.2.2203', 'mesh_dim': 2, 'order2d': 2, 'order3d': 2, 'draw_vol': Fals…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0abdc3a0931d44d1996fceb4d812aa8c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(value={'ngsolve_version': '6.2.2203', 'mesh_dim': 2, 'order2d': 2, 'order3d': 2, 'draw_vol': Fals…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1c4c72c6c3a34f1d9889e67aa953a9b6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(value={'ngsolve_version': '6.2.2203', 'mesh_dim': 2, 'order2d': 2, 'order3d': 2, 'draw_vol': Fals…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torsional rigidity = 35.972111800921475\n"
     ]
    }
   ],
   "source": [
    "# cell problem\n",
    "from ngsolve import *\n",
    "import netgen.geom2d as geom2d\n",
    "from ngsolve.webgui import Draw\n",
    "\n",
    "geo = geom2d.SplineGeometry()\n",
    "geo.AddCircle(c=(-0,0), r=1, bc=\"c\", leftdomain=2, rightdomain=1)\n",
    "geo.AddRectangle((-2,-2), (2,2), bcs=[\"b\",\"r\",\"t\",\"l\"], leftdomain=1)\n",
    "\n",
    "geo.SetMaterial(1, \"outer\")\n",
    "geo.SetMaterial(2, \"inner\")\n",
    "\n",
    "# generate a triangular mesh of mesh-size maxh\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.05))\n",
    "\n",
    "fes1 = H1(mesh, definedon=\"inner\")\n",
    "u1 = GridFunction(fes1, \"u1\")\n",
    "u1.Set (1.01)\n",
    "\n",
    "fes2 = H1(mesh, definedon=\"outer\")\n",
    "u2 = GridFunction(fes2, \"u2\")\n",
    "u2.Set (1)\n",
    "\n",
    "\n",
    "fes = H1(mesh, order=3,dirichlet=\"b|l|r|t\")\n",
    "u = fes.TrialFunction()\n",
    "v = fes.TestFunction()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += 2*v*dx\n",
    "f.Assemble()\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += SymbolicBFI (u1*grad(u)*grad(v), definedon=mesh.Materials(\"inner\"))\n",
    "a += SymbolicBFI (u2*grad(u)*grad(v), definedon=mesh.Materials(\"outer\"))\n",
    "a.Assemble()\n",
    "\n",
    "gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec\n",
    "\n",
    "sts = [-grad(gfu)[1], -grad(gfu)[0]]\n",
    "\n",
    "Draw (u1); Draw (u2)\n",
    "Draw (gfu)\n",
    "#Draw (sts, mesh, \"Flux\")\n",
    "\n",
    "print (\"torsional rigidity =\", Integrate ( u1*grad(gfu)*grad(gfu), mesh)+Integrate ( u2*grad(gfu)*grad(gfu), mesh))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d0e6ff6f",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'u1' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp\\ipykernel_15196\\1772886189.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     31\u001b[0m \u001b[0msts\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mgrad\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgfu\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m-\u001b[0m\u001b[0mgrad\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgfu\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\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[0;32m     32\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 33\u001b[1;33m \u001b[0mDraw\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mu1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m;\u001b[0m \u001b[0mDraw\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mu2\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     34\u001b[0m \u001b[0mDraw\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mgfu\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     35\u001b[0m \u001b[1;31m#Draw (sts, mesh, \"Flux\")\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'u1' is not defined"
     ]
    }
   ],
   "source": [
    "from ngsolve import *\n",
    "import netgen.geom2d as geom2d\n",
    "from ngsolve.webgui import Draw\n",
    "\n",
    "geo = geom2d.SplineGeometry()\n",
    "geo.AddCircle(c=(-0,0), r=1, bc=\"c\", leftdomain=2, rightdomain=1)\n",
    "geo.AddRectangle((-2,-2), (2,2), bcs=[\"b\",\"r\",\"t\",\"l\"], leftdomain=1)\n",
    "geo.SetMaterial(1, \"outer\")\n",
    "\n",
    "# generate a triangular mesh of mesh-size maxh\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.05))\n",
    "\n",
    "fes = H1(mesh, order=3,dirichlet=\"b|l|r|t|c\", definedon=\"outer\")\n",
    "u = fes.TrialFunction()\n",
    "v = fes.TestFunction()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += 2*v*dx\n",
    "f.Assemble()\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += u * v * dx\n",
    "a.Assemble()\n",
    "\n",
    "gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec\n",
    "\n",
    "sts = [-grad(gfu)[1], -grad(gfu)[0]]\n",
    "\n",
    "Draw (u1); Draw (u2)\n",
    "Draw (gfu)\n",
    "#Draw (sts, mesh, \"Flux\")\n",
    "\n",
    "print (\"torsional rigidity =\", Integrate ( grad(gfu)*grad(gfu), mesh))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "647d4cb9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}