2022-03-15 07:04:08 +00:00
{
"cells": [
{
"cell_type": "code",
2022-03-16 14:44:07 +00:00
"execution_count": 125,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import mpl_toolkits.mplot3d\n",
"import numpy as np\n",
"import scipy.sparse as scps\n",
"import scipy.sparse.linalg as ssl\n",
"import math"
]
},
{
"cell_type": "code",
2022-03-16 14:44:07 +00:00
"execution_count": 126,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
"def maillage_carre(n: int):\n",
" \"\"\"\n",
" Une discrétisation possible d'une EDP elliptique sur le domaine ]0,1[ x ]0,1[.\n",
" Le carre [0,1]x[0,1] est maille uniquement avec des triangles.\n",
" Les conditions limites sont de type Dirichlet uniquement -> `neumann=[]`.\n",
"\n",
" Args:\n",
" n (int): nombre de points par cote du care => Npts points de discretisation au total\n",
"\n",
" Returns:\n",
" coordinates : matrice a deux colonnes. Chaque ligne contient les coordonnes 2D d'un des points de la discretisation. Ces sommets seront identifies a l'indice de la ligne correspondante dans la matrice coordinates.\n",
" elements3 : matrice a trois colonnes. Chaque ligne contient les indices des sommets d'un element triangle, dans le sens antihoraire.\n",
" dirichlet : vecteur colonne des indices des sommets de la frontiere de Dirichlet.\n",
" neumann : matrice a deux colonnes. Chaque ligne contient les indices des deux sommets d'une arete de la frontiere de Neumann. (neumann est vide sur cet exemple)\n",
" \"\"\"\n",
"\n",
" h = 1 / (n - 1)\n",
" n_pts = n * n\n",
" n_elm = 2 * (n - 1) * (n - 1)\n",
" coordinates = np.zeros((n_pts, 2))\n",
" elements3 = np.zeros((n_elm, 3), dtype=int)\n",
" neumann = []\n",
" dirichlet = np.zeros((4 * n - 4, 1), dtype=int)\n",
"\n",
" # Coordonnees et connectivites :\n",
" e = -1\n",
" p = -1\n",
" x = np.zeros((n + 1, 1))\n",
" x[n, 0] = 1.0\n",
"\n",
" for l in range(n + 1):\n",
" x[l, 0] = l * h\n",
"\n",
" for j in range(n):\n",
" for i in range(n):\n",
" p = p + 1\n",
" coordinates[p, 0] = x[i, 0]\n",
" coordinates[p, 1] = x[j, 0]\n",
" if (i != n - 1) & (j != n - 1):\n",
" p1 = p\n",
" p2 = p1 + 1\n",
" p3 = p1 + n\n",
" p4 = p2 + n\n",
" e = e + 1\n",
" elements3[e, 0] = p1\n",
" elements3[e, 1] = p2\n",
" elements3[e, 2] = p3\n",
" e = e + 1\n",
" elements3[e, 0] = p4\n",
" elements3[e, 1] = p3\n",
" elements3[e, 2] = p2\n",
"\n",
" # Liste des sommets de la frontiere de Dirichlet:\n",
" p = -1\n",
" for j in range(n):\n",
" p = p + 1\n",
" dirichlet[p, 0] = j\n",
"\n",
" for j in range(n * 2 - 1, n * (n - 1), n):\n",
" p = p + 1\n",
" dirichlet[p, 0] = j\n",
"\n",
" for j in range(n * n - 1, n * n - n - 1, -1):\n",
" p = p + 1\n",
" dirichlet[p, 0] = j\n",
"\n",
" for j in range(n * n - 2 * n, n - 1, -n):\n",
" p = p + 1\n",
" dirichlet[p, 0] = j\n",
"\n",
" return coordinates, elements3, dirichlet, neumann\n"
]
},
{
"cell_type": "code",
2022-03-16 14:44:07 +00:00
"execution_count": 127,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
"def show(coordinates, u) -> None:\n",
" \"\"\"Fonction d'affichage de la solution u sur le maillage defini par elements3, coordinates.\n",
"\n",
" Args:\n",
" elements3 : matrice a trois colonnes contenant les elements triangles de la discretisation, identifies par les indices de leurs trois sommets.\n",
" coordinates : matrice a deux colonnes contenant les coordonnes 2D des points de la discretisation.\n",
" u : vecteur colonne de longueur egale au nombre de lignes de coordinates contenant les valeurs de la solution a afficher aux points de la discretisation.\n",
"\n",
" Returns:\n",
" None, plots a figure\n",
" \"\"\"\n",
"\n",
" ax = plt.figure().add_subplot(projection=\"3d\")\n",
" ax.plot_trisurf(\n",
" coordinates[:, 0], coordinates[:, 1], u, linewidth=0.2, antialiased=True\n",
" )\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Partie I : maillage triangulaire et conditions de Dirichlet**\n"
]
},
{
"cell_type": "code",
2022-03-16 14:44:07 +00:00
"execution_count": 128,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
"def f(x, y):\n",
2022-03-15 07:05:49 +00:00
" return 2 * np.pi ** 2 * np.sin(np.pi * x) * np.sin(np.pi * y)\n",
2022-03-16 14:44:07 +00:00
" # try:\n",
" # return np.ones(x.shape[0])\n",
" # except:\n",
" # return 1\n",
2022-03-15 07:04:08 +00:00
"\n",
"\n",
"def u_ex(x, y):\n",
" return np.sin(np.pi * x) * np.sin(np.pi * y)\n",
"\n",
"\n",
"def u_d(x, y):\n",
2022-03-16 14:44:07 +00:00
" # return np.zeros(x.shape[0])\n",
" # return 1\n",
" return np.ones(x.shape[0])\n",
2022-03-15 07:04:08 +00:00
"\n",
"def g(x):\n",
2022-03-16 14:44:07 +00:00
" # return np.cos(x)\n",
" return 1"
2022-03-15 07:04:08 +00:00
]
},
{
"cell_type": "code",
2022-03-16 14:44:07 +00:00
"execution_count": 129,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [
{
"data": {
2022-03-16 14:44:07 +00:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAADyCAYAAACLfbNuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACEsElEQVR4nO29d3hkd3X//7rTm0Yzo96llbb36oZtSgBj3Kg2JAZTY0I3JIGEJ3EK5QGSfPPFlN8XnJhQjLEN2GCzAQw2GBt7d70rrXZXq97LjKb3en9/jO7dGWmatJJ21zvv5/EDq7lz752Zez7nfM55n/cRRFGkjDLKePlDcaFvoIwyylgflI29jDIuE5SNvYwyLhOUjb2MMi4TlI29jDIuE5SNvYwyLhOoirxersuVUcbaQ1iPi5Q9exllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigbexllXCYoG3sZZVwmKBt7GWVcJigb+wWAKIrEYjESiQRlKe8y1gvFxCvKWGWkUilisRiRSET+m1KpRK1Wo1KpUCqVCMK6aBmUcZlBKOJZym5nlSCKIolEgkQigSAIxONx+e+iKJJKpWQjj0ajVFRUoNFoysZ/eWBdfuCyZ18HSGF7pkFLEAQBQRBQKBTysUNDQ7S3t2MwGICy5y9jdVA29jVGIpFgcnKSZDJJU1MTgiDI3jyX0UrGr1QqUSqVstcPh8Py8SqVSv6vbPxllIqysa8RMsP2VColh+/LRS7Pn0wmSSQS8jEqlUr2/AqFomz8ZeRE2djXAKlUing8LoftkjcvFYWOl84nYbHxC4KQ5fnLxl+GhLKxryIkw5OSb5I3zme8+UL55SCX8ScSCfkepKhApVKh0WjKxn8Zo2zsqwRRFInH4ySTySUGuNjYi3n75UYCi9+72PgnJycBaGhoKHv+yxhlY18FSLVzyVPnyrhfKPJM5v1ICb94PJ7l+aU9v1KpLBv/yxhlYz8PLK6dS2H7Yqzmnv18IWX6JeQyfinZp1Kpci5eZVyaKBv7CrG4dl7IIBYb7+zsLP39/ajVaiwWC1arlcrKyiVGuJrId3+5jD8WixGNRoF03kGtVsuev2z8ly7Kxr4CSEm4fGH7YkjGnkwm6evrIxaLsX//fkRRxOv1Mj8/z9DQEEqlEqvVKmfyLwQKGX9msi8z7C/j0kCZLrsMSGF7T08PXV1d6HS6kt7ndDqZnp4mEAjQ1NRES0tL1mIhIRqN4vF4GBkZQRRFDAYDVqsVq9WKyWRasUednJxEoVDQ2Ni4ovdLkJ6VVCrF6dOn2bJlixzul43/vFCmy15MyKydS8m4UuF0OnE4HBw4cACz2Zz3OK1WS11dHR6Ph/r6ejQaDR6Ph4mJCfx+f5bxGwyGdQ+nMxN9kUhEjgBisRixWAyg7PkvYpSNvQgW186lULYUY08kEpw5c4ZIJEJtbW1BQ8+EFPbr9Xr0ej0NDQ2Iokg4HMbtdjMyMkIwGMRkMsnGr9frz+tzrgSZ1F4419Sz2Pgzef1l479wKBt7AeSrnQuCUHRP7ff7OXnyJK2trRiNRqanp8/rXgRBwGAwYDAYaGpqQhRFgsEgbrebgYEBIpEIFRUVsvFrtdrzul6p97T434tr/KIoEo1GlyT8ysa//igbex4Uqp0XKo2JosjU1BTj4+Ps3LmTiooKvF7vqpfeBEHAZDJhMploaWkhlUrh9/txu92cPn2aeDxOZWUlVquVZDJ5QYyqmPGLopgV8kulvjLWBmVjX4TMsD1f7VyhUOT07IlEglOnTqFQKDh06BAqVfrrXQ9SjUKhoLKyksrKStrb20mlUni9XtxuN3Nzc3IkYLVasVgs8r2tJ3IZfyqVkoU8pqenaW5uRqPRlDv61gBlY89AqbXzXMbr8/no7e2lvb19SdY71/FrRZeVoFAo5JBeo9EAoNfrcbvdjI6OIghC3hr/emHxdzw7O0tjY2NZxWeNUDb2BRSjvGYiM0EniiITExNMTU2xe/dujEbjkuMvJF1WgkKhoKqqiqqqKgDi8Tgej2dJjd9qtWI2my/YXnpxjX9xL3/Z+FeOy97YS6W8ZkJK0MXjcXp7e9FqtRw6dCivd7yY6LIS1Go1NTU11NTUAOnymdvtltl9Wq12VWr854Ncvfxl4185LmtjT6VS2O12RFHEarWW/KAIgoDf76evr4/Ozk7q6+uLHn+hPXuxz6bRaKirq6Ourg6ASCSC2+3OW+O/ECjF+MsqPvlxWRp7ZhLO7/cjiiI2m63k9/p8PiKRCPv37y/pwb/Qnn0l59LpdDQ0NOSt8UejUaanp7Fareh0uovG82cmVyVoNBq0Wu1l39F32Rn74rBdqVTKhJliiMVinDx5kkQiQVdXV8ke7mLw7OeDXDX+F154gWQyecFq/Pnuc7Hxj4+Po1arqa2tzWrnvRx7+S8rY88lF5WvjLYYUv1648aNBAKBZT0kF9qzrzakRbKlpUWu8QcCAVwu15Iav8VikasBF+I+pVq+SqW67CW8LgtjL1Q7L0Z9FUWR4eFh5ufn2bdvH3q9nmAwuKyutJebsS+GQqHAbDZjNpuX1PglZV2pzLfeNf5UKpUlD1ZMwuvlbPwve2MvJBcFhamv0WiUkydPUlFRwcGDB4tqyuXDpWa854vMGj+kW4I9Hs8FqfFnGvti5DL+l7OKz8va2EupnecL451OJ319fWzatEkuTxV7Tz4sNvZgMEhPTw+CIGC1WrHZbFRUVKx4MbnYoVQqS6rxS7Lbq1njX875lqPicyka/8vS2JdTO19suKIoMjg4iMfjYf/+/Tl71s/Hs8/OzjI8PMyWLVtQq9V4PB6mp6fx+/3odDqsViuxWOyCJbnWA7lq/C6Xi0QiwdGjR9FoNHJkUFFRcV4GdT6LRyEhj+npaerq6jAYDJeMhNfLztiXIxcF2cYeiUTo6enBZrNx4MCBgnTZZDJZ8j1JWwWp3VXaEiQSCbm2nVnecrlc2O12XC6X7Plfzsav0WioqalhamqK/fv3yzX+yclJAoGAvAharVaMRuOyDGo1I4VM43e5XNTV1WWp+Eie/2Lt5X9ZGbsUcklffCkPhWTsDoeD/v5+tmzZIoebhd6zHM8eiUQIBoM0NjayZcsWBEHImugC2eWtZDKJWq3GZDLJGe5EIiFnuK1W6wVpZFlLZBplvhr/6OgowWAQo9GY1cdf6Hde7W2BhGQyKY/ognNcBqmX/yMf+Qif+9zn2LJly6pfe6V4WTwxUtje19eHxWKhtrZ2We/3eDxEo1EOHDhQkgddThgvLSJarZb29vaS3iM9vBUVFVRUVNDW1kYymZQz3GNjY/J+X0pyFXugL+bwEvIPzCjUxz84OEgkEpFFPHJFQGtl7KlUKivEz6TvQtrzXwhBkUK45I09s3a+3MRZOBzm1KlTCILA/v37SzaIUq4j7f29Xi8HDx7k6NGjJd9XLiiVSmw2m8z0i8fjuN1u7HY7AwMDaDQabDbbunDZ1yJ5KIpiyX0Ji/v4A4EAbrebM2fOEIvFsmr8a2Xs0r3kg6QkdDHhkjX2XKOWlEplycYuGUlnZyfT09OrSpKJRqP09PRgtVoLLiKFvFkxg5JYYVIUE4lEcLlcjI+PEwgEMBqNsvGvtnGuxtiqxcg1zroUZNb429raSKVS+Hw+ec/v9/sZGRmhqqpqVWv8xb7TUChUNvbVQL7aeSkeN5VKcfbsWcLhMAcPHswaj1QqCtXmJabd5s2bqa6uXtZ5zwc6nY7GxkYaGxvlUNflctHf34/f78doNMr1b7VavW73VSpWawFRKBRYLBYsFgsdHR0cOXKEqqoqvF4vY2NjAOtS44/H4xeMOZgPl5yxF6qdK5XKglnyUChET08P9fX1cqJsJRrtuRJ0oigyOjqK3W7PW7IrFedbZ88MdVtbWxkbGyORSBAIBJiYmEAURSwWCzabbUUP/Fp59rUKt6urq+WFN5FI4PF4cDqd59XHf7HnQHLhkjH2UmrnCoUib1PL7OwsQ0ND7Ni
2022-03-15 07:04:08 +00:00
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"n = 10\n",
"coords, elems3, dirichlet, neumann = maillage_carre(n)\n",
"show(coords, f(coords[:, 0], coords[:, 1]))\n",
"# show(coords, u_ex(coords[:, 0], coords[:, 1]))\n",
"# print(maillage_carre(3))"
]
},
{
"cell_type": "code",
2022-03-16 14:44:07 +00:00
"execution_count": 130,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
"def raideur(triangle):\n",
" M = np.zeros((3, 3))\n",
" x = triangle[:, 0]\n",
" y = triangle[:, 1]\n",
"\n",
" # calcul de alpha\n",
" mat_alpha = np.array(\n",
" [\n",
" [x[1] - x[0], x[2] - x[0]],\n",
" [y[1] - y[0], y[2] - y[0]]\n",
" ]\n",
" )\n",
" alpha = np.linalg.det(mat_alpha)\n",
"\n",
"\n",
" for i in range(3):\n",
" grad_eta_i = np.array(\n",
" [\n",
" y[(i+1)%3] - y[(i+2)%3],\n",
" x[(i+2)%3] - x[(i+1)%3]\n",
" ]\n",
" )\n",
" for j in range(3):\n",
" grad_eta_j = np.array(\n",
" [\n",
" y[(j+1)%3] - y[(j+2)%3],\n",
" x[(j+2)%3] - x[(j+1)%3]\n",
" ]\n",
" )\n",
"\n",
" M[i, j] = np.dot(grad_eta_i, grad_eta_j)\n",
"\n",
" return M / alpha / 2"
]
},
{
"cell_type": "code",
2022-03-16 14:44:07 +00:00
"execution_count": 131,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [
{
2022-03-15 07:05:49 +00:00
"data": {
2022-03-16 14:44:07 +00:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAADyCAYAAACLfbNuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAC+1UlEQVR4nOz9eZRdWXbWi/7W7k8fcaIPhfoulcpeKdku7AIuPC7w7oDBfdQDc9+zDeaCMfajuxgwxmBz7QsYYxv7YsC4KRqXwQbcFlUu92W70qnMlJTqUr0Uir45/Tm7X+v9sfY5pVSqiYiUlMrK+MbIkZkRp9lxzp5rzTXnN79PKKXYxja28aUP44O+gG1sYxtPBtvBvo1tfESwHezb2MZHBNvBvo1tfESwHezb2MZHBNvBvo1tfERgPeT32325bWzj8UM8iTfZ3tm3sY2PCLaDfRvb+IhgO9i3sY2PCLaDfRvb+IhgO9i3sY2PCLaDfRvb+IhgO9i3sY2PCLaDfRvb+IhgO9i3sY2PCLaDfRvb+IhgO9i3sY2PCLaDfRvb+IhgO9i3sY2PCLaDfRvb+IhgO9i3sY2PCLaD/QOAUoooikiShG0p7208KTxMvGIbjxhSSqIoIgiCwc9M08S2bSzLwjRNhHgiWgbb+IhBPGRn2d52HhGUUiRJQpIkCCGI43jwc6UUUspBkIdhSKlUwnGc7eD/aOCJfMHbO/sTQD9tvzOg+xBCIITAMIzBY69du8aePXvI5/PA9s6/jUeD7WB/zEiShLm5OdI0ZceOHQghBrv5vYK2H/ymaWKa5mDX931/8HjLsgb/bAf/NjaK7WB/TLgzbZdSDtL3zeJeO3+apiRJMniMZVmDnd8wjO3g38Y9sR3sjwFSSuI4HqTt/d18o3jQ4/uv18fdwS+EeNfOvx382+hjO9gfIfqB1y++9Xfj+wXv/VL5zeBewZ8kyeAa+lmBZVk4jrMd/B9hbAf7I4JSijiOSdP0PQF4d7A/bLffbCZw93PvDv65uTkApqamtnf+jzC2g/0RoN877+/U96q4f1DkmTuvp1/wi+P4XTt//8xvmuZ28H8JYzvY3wfu7p330/a78SjP7O8X/Up/H/cK/n6xz7Ksey5e2/hwYjvYt4i7e+cPCoitBO+jDvb7Xd+9gj+KIsIwBHTdwbbtwc6/HfwfXmwH+xbQL8LdL22/G1vZ2T8oPCj47yz23Zn2b+PDge1vahPop7ynT58mDMMNn28/6DT+/bxWP/j7xTyAKIo4efIkjUaDVqtFr9cbtBq38fRie2ffIO7snfeLcRvFB1mge5S4s9AXBMEgA4iiiCiKALZ3/qcY28H+ENzdO++nso8z2D8si8Od1F744lDP3cF/J69/O/g/OGwH+wNwv965EGJTKeuHJXg3i/sN9fTRD/4wDN9T8NsO/ieP7WC/Dx7UO3/cO/WXyuLwsOBXSr0r5e+3+rbxeLAd7HfhzrT9fr1zwzDe184ehiEXL17EcRyq1SqVSuVdFfAvVdwr+KWUAyGPhYUFZmZmcBxne6LvMWA72O/ARnvn72enXl9f55133mHfvn1IKVlfX+fatWuYpsnw8DDVanVwLV/quPszXlpaYnp6elvF5zFhO9gzPIzyeie2UqCTUnL16lXq9TrHjh3DNE2klIyNjQG6ol2v11lYWGBtbY1ms0mn06FarVIoFD4yN/ndPf67Z/m3g3/r+MgH+0Ypr3diswW6OI6p1WoUCgVeffXVd8lS9eE4DhMTE0xMTHDt2jXy+TxKKW7dukWn06FQKDA8PMzw8DC5XO4jcZPfa5Z/O/i3jo90sEspWVlZQSnF8PDwhm+UzaTx6+vrXLhwgXw+z8GDBzd8bY7jMDIywvT0NEoput0u9Xqdq1evEgQBpVJpEPyu627omj/s2Ejwb6v43B8fyWC/swjXbrdRSg3OyhvBRgp0SimuX7/O+vo6zz//PNevX9/w699rJLZYLFIsFtm5cydSStrtNvV6nQsXLpAkCZVKZRD8lvXur/VL9fx/PxWffnG1D8dxcF33Iz/R95EL9rvTdtM035NSPwwP29nDMOTs2bOUy2VeffXVTTPuHgbDMKhUKlQqFfbs2UOapjSbTer1OrOzswAMDQ0xPDxMpVJ5ZO/7tONewT87O4tt24yPj79rnPejOMv/kQr2e8lFbbaNBg8u0NVqNS5evMihQ4cGxbfH3Wc3TZNqtTrITpIkoV6vs7a2xrVr10iShGKxSD6fp1QqfWSILP3PsZ/Wf9QlvD4Swf6g3vlmK+tw7wLdnWn7sWPH8DzvXY9/kqQay7IYGxsbLDY3btwgDEMWFhZot9u4rku1WmV4ePhLvtIvpXyXPNjDJLy+lIP/Sz7YHyQXBZuvrPefc2cwRlHE22+/PUjb7945P2hGnGVZeJ7H1NQUAL7vU6/XuXnzJt1ul0KhMAj+XC73gV3n48CdwX437hX8X8oqPl/Swb6R3vlW0/j+c+6Vtt+Np40um8vlyOVy76n0X758eeBG0yf4OI7z2K7jTjyuv/dBwX43NqPi82EM/i/JYN9M73wrwd7PBq5fv87q6up70vZ7Pf5prYg/qNJ//vz5h1b6HxUehdLuvbCZYL8bDxLyWFhYYGJignw+/6GR8PqSC/bNyEXB1oI9TVNu377N+Pg4x48ff+jNdK9gf9D7fpCLw0Yr/UmSkKbpI+P0P43BfjfuDP5arcbExMS7VHz6O//TOsv/JRXs/ZSr/8Fv5ObZbLD3z7ojIyMcPnx4Q8+5V/A+rTv93bi70h/HMY1Gg4WFBd566y0syxrs+u+n0v8og/JJvG5/obtzlh++KOTxTd/0TXzbt30bzzzzzCN/763iSyLY+2n7O++8w9DQEOPj4xt+7kaDXSnFjRs3WF1dZe/evaRp+n4u+YF4HDv7o9o1bdtmbGyMmzdvcvz4ccIwHHD62+02nucNgn8zlf4Pw85+9+vemdXcSd8FvfM/bcXOD32w39k7f9Q98z6iKOLs2bMUi0WOHz/O6uoq7XZ7w+/xNJ/jtoI7Py/XdZmcnGRychKl1Hsq/cVi8V2c/ge95odpZ4cHf6/9v/1pwoc22O9ltdSfJNsMHtZ661NSDx48OMgYHtXO+yAn16dFcPJ+r3e/687n8+TzeXbs2DGo9NdqtXdV+vttvjsr/feys34UeFzB/rDPtNfrbQf7o8D9eufvt41293vcvHmTlZUVXnnllXftSlvpzX8UcWelf9euXYNKf61WY35+njRNB7Rex3E+VMH+MMRx/MTalhvFhy7YH9Q7N01z02fpewV7P20vFAr3rLZvhXW3GTzNrTrY+vn6zkp/v+7RbDap1WrUajXCMOTatWtUq1XK5fIjqfQ/ruPBh/Fo9qEJ9o30zg3D2PRQy93B3mg0OH/+/LvS9rvxtAfj48ajKqbdWelvt9vMzs5SLpdZXV3l6tWrg0p/tVqlVCo9NQH2sO/+cRUb3y8+FMG+0d75+0nj70zbX375ZfL5/EOf87jwtB8THsdC1x9YuZPT36/0z8/Pv69K/6PGRvgFT2PAP/XBvhmrpa0W6JRSnDp1ilwut2WSzKPE03aT3AuP+hrvVaB7VJX+R42HBXuSJE+lgOhTG+xbkYsyDGPTZ/ZGo0G32+XAgQNMTExs6DmPYmd/WLA8zceEx7FrPew171Xp73Q67+L0l8vlQfA/zuJYmqYPvB/7w0VPG57KYE/TlE6nM6jQbvTG2kwQ9vXdlpaWyOfzGw50eP87u5SSy5cvEwQBIyMj75GWetprAo8j2DdbNRdCUCqVKJVKg0p/q9UapP39Sn9/w3iUnP67CTV342nsscNTFuz93nkYhpw6dYov//Iv39RNtdFgj+OYs2fPksvlOHHiBK+99tqmrvP9BGMQBJw5c4axsTFGRkYGBcH+zVmtVreUNcyvtviNczdZqHdotAOCOMa1TfKeg0gint01wh+rjuA59pau+2486Z39YTAMg6GhIYaGhgaV/kajwfLyMqdPn0YIMdj1K5XK+6rQPyyN397ZH4I7e+f91tZmv/yNtN6azSbnzp3bVNp+N7aaxvc1448cOUKlUiG
2022-03-15 07:05:49 +00:00
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPcAAADyCAYAAACRQVPgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAADENElEQVR4nOz9d5heWX7Xi37WWju9sXJQlVTKsdXqoFbPjDn3mnuI9nNteLgc7oXLwVwbzsUH+8A55gAGXxwwDjgQbGODx8YmGYPHOHucxvYAnh6pW2rlHEqqKlUOb9pxrXX/2PXWlKolVWhJ091T3+fR0+rSu/d+a+/13b/f+oXvT1hr2cY2tvHRg/xif4FtbGMbzwfb5N7GNj6i2Cb3NrbxEcU2ubexjY8otsm9jW18RLFN7m1s4yMKZ51/386TbWMbzx/ieZx023JvYxsfUWyTexvb+Ihim9zb2MZHFNvk3sY2PqLYJvc2tvERxTa5t7GNjyi2yb2NbXxEsU3ubWzjI4ptcm9jGx9RbJN7G9v4iGKb3NvYxkcU2+TexjY+otgm9za28RHFNrm3sY2PKLbJvY1tfESxTe4vAqy1JElClmVsS0tv43lhPbGGbTxjGGNIkoQoilZ+ppTCdV0cx0EphRDPpXd/G19iEOtYjm2z8oxgrSXLMrIsQwhBmqYrP7fWYoxZIXUcx1QqFTzP2yb7lwaeywPettwvAG03fDWB2xBCIIRASrny2du3b7Nnzx6KxSKwbdm3sTVsk/s5I8syxsbG0FozPDyMEGLFWj+OpG2yK6VQSq1Y9TAMVz7vOM7Kn22yb+NJ2Cb3c8JqN9wYs+KObxaPs+xaa7IsW/mM4zgrll1KuU32bQDb5H4uMMaQpumKG9621hvF0z7fPl8ba8kuhHjEsm+T/UsX2+R+hmgTrR0sa1vbJ5H1Sa75ZvA4smdZtvId2lbfcRw8z9sm+5cQtsn9jGCtJU1TtNbvIdxacq9nzTdr6dceu5bsY2NjAOzYsWPbsn8JYZvczwDt3HXbEj8uIv7FKlZZ/X3aAbo0TR+x7O09u1Jqm+wfIWyT+31gbe667YavxbPcc79ftCPxbTyO7O3gnOM4j31ZbePDgW1ybxFrc9dPI8BWyPqsyf2k7/c4sidJQhzHQB43cF13xbJvk/3Dg21ybwHtoNmT3PC12Irl/mLhaWRfHZxb7cZv44OJ7SezCbRd2HfffZc4jje8P/1iu+Xv51xtsreDbwBJknDmzBkWFxep1Wq0Wq2V1N82PjjYttwbxOrcdTt4tlF8MQNqzxKrA3NRFK1Y+CRJSJIEYNuyf4CwTe51sDZ33XZNnye5Pywvg9WlsvCFJpi1ZF9dF79N9heHbXI/BU/KXQshNuWCfljIulk8qQmmjTbZ4zh+T4Bum+zPH9vkfgKelrt+3pb4o/IyWI/s1tpHXPh26m0bzwbb5F6D1W74k3LXUspty70FPI7sxpgV4YqJiQl27tyJ53nbHW/PANvkXoWN5q7fryW21vLw4UMcx6GrqwvHcZ76+Y8q1t7jyclJhoaGtlVqnhG2yb2M9UpIV+P9BNTiOObChQuUSiUARkdHEULQ1dVFV1cXHR0d7+8X+ZBjbY59bS/7Ntk3ji95cm+0hHQ1thpQm5ub49q1axw+fJiOjo4VDyFNUxYXF5mZmeHWrVukaUqSJDiOQ6VS+ZJdwI/rZd8m+8bxJU1uYwzT09NYa+nq6trwwtiK21yv17l9+zYnT54kCIKV1BqA67r09fXR19cHwPXr11FKMTY2RqPRIAgCuru76erqolgsbln04cOOjZB9W6XmC/iSJPfqoFm9XsdaS3d394aP30xALY5jLl26hLWWN954Y0Oegeu6VKtVent7sdYShiELCwvcuXOHMAwplUorZA+CYN3zfVT3709SqWkHQ9vwPA/f97/kOt6+5Mi91g1XSj1iRTeCjVru+fl5rl69yt69e3n48OGWcrpCCIrFIsVikeHhYay1NBoNFhYWuHbtGkmS0NHRsbJnd11309f4qOBxZL9//z6u69Lf3/9Ie+uXQi/7lxS5Hyd/tNm0FqwfULPWcufOHebm5jh58iRCCCYmJjZ8/vWEHCqVCpVKhZGREYwx1Go15ufnGRsbwxhDZ2cnXV1ddHZ2PhKg+lJD+z623fQvNUmqLwlyPy13vdnINzw9oJYkCRcuXKBSqay44c+zFl1KSWdnJ52dnUCutrq0tMT8/Dx3795FSomUkkqlgjHmS64ibPXvvBFJqo8S2T/y5H6a/BFsPvLdPuZx5FtYWODKlSscOnRoJTj2tM8/DziOQ09PDz09PUD+srl16xZLS0u8/fbb+L6/4sKXy+UP9eLdCJ72Qnsc2T9KKjUfaXJvJHe9Vbd89THWWu7evcvMzAyvv/46hULhkc9/MctPPc+jUqnQ1dXFjh07iKKI+fl57t+/T6PRoFgsrgTnCoXCF23xPq+X32a8lc2o1HwYyP6RJPdmctdbIfdq8iVJwsWLFymXy5w6deqx1/ogVZwFQcDQ0BBDQ0NYa2m1WiwsLHDr1i2iKKJcLq+Q3ff9F/a9noUS7OPwfrYiTxOumJiYYGBggGKx+IGVpPrIkXsz8kewdXJrrVfc8IMHD9Lf3//Uz68l99Ou+6JeBkIISqUSpVKJnTt3Yq2lXq8zPz/PlStXyLJsJRLf2dn5XCPxH0Ryr8Vqss/PzzMwMPCISk3bsn9Qetk/UuRuu1DtG72RxbJVcs/NzTE1NfVYN/xxn19L1g+KJV8NIQTVapVqtcqePXvQWq9E4u/fv79S7NPV1fXMv//zCvY9r/NqrVdGPgGPeHJJkvAN3/ANfMu3fAtHjhx55tfeKD4S5G674deuXaOzs/OpVnQtNkvuNE0ZHR1FSsmbb775XBbO87DcW7GKSqkVMkMeiV9cXGR2dpZWq8W5c+dW/r1Sqbyve/FhsNxrz7vaZV9dDgu5ZV/vpf+88aEn9+rc9fPIWa/G4uIily9fpre3d1Nu1wdpH/Z+4DgOvb299PT0sLS0xEsvvcT8/DwTExPU63WCIFghe6lU2tTvba39UFluePpzbTablMvl53LdjeJDS+7Hje5RSm3JxV7vGGsto6OjTE5O8tprr9FoNFhaWtryd1993sctkGdtuZ+1F9D+3p7nMTg4yODgIABhGDI/P8+9e/dWFneb7OtZsceNN34WeF7kXu+etlqtbXJvBU/KXT+LtNZapGnKpUuXCIJgxQ1vNpvbSp+PQaFQYHh4eKVMttlsMj8/z40bN4jjmEqlshKJ9zzvkWM/bG75ekjT9D2/44vGh47cT8tdK6XQWm/qfE8j99LSEpcvX2bfvn0r1ql9zPuxhm2xhjRN6enped958ReNjRBRCEG5XKZcLq+UybYj8ePj42itHymTfZ4W9nnFRT7o+NCQeyO5aynlpptAHkfudsPBw4cPefXVVykWi4/8+/shn9aay5cvA1AsFlesWrVaXbFqH3RsxcpKKeno6KCjo4O9e/eitWZxcZGFhQXu3bu38gwWFhbo6Oj4oqeRnob1nv3z8kI2iw8FuTeau34Wbnmaply+fBnP8zh16tRjGy+2ch3IgywXLlxg165dDA4OkmXZe5o/Hjx4QJIkFAoFisXiB3KhPwuvQin1SJlsO902PT3NrVu3cBxn5WX3QROsaKfBnoYPAsE/8OTezOierQbU2mi74Xv37mXHjh1PPWazCzxNU86fP8/x48epVqsrnUnw3uaPtkjD9PQ0N2/eXKkH7+7u3nQU+nnhWX8HKSXFYpFDhw4BeR98u9Ot0WhQKBRWgnNbFax4VliP3FmWfSC68T6w5N6K/JGUctN77va17t+/z/j4OK+88sqKvtnTrrPRl4gxhps3b5KmKV/2ZV+2UuW13kuqVCqxa9cugBWxhtVR6O7ubrq7u19oiWgbz8MqrT2n7/vs2LGDHTt2rAhWzM/PrwhWrI7Eb0Sw4llCa/3U9dhsNtddQy8CH0hya61pNBp4nrepet2tuMtZlhGGIfV6nTfffHNDb9yNWu62GGJ3dzeFQmHD5Ztrz18oFCgUCiv14I1G45ES0dW
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
2022-03-15 07:04:08 +00:00
}
],
"source": [
"def assemblage(coords, elems3):\n",
" Ns = len(coords)\n",
" A = np.zeros((Ns, Ns))\n",
" for triangle in elems3:\n",
" M = raideur(coords[triangle])\n",
" for i, a in enumerate(triangle):\n",
" for j, b in enumerate(triangle):\n",
" A[a, b] += M[i, j]\n",
" return A\n",
"\n",
"\n",
"def second_membre(coords, elem3, f):\n",
" Ns = len(coords)\n",
" b = np.zeros(Ns)\n",
" for triangle in elem3:\n",
" coords_triangle = coords[triangle]\n",
" centre = np.mean(coords_triangle, 0)\n",
"\n",
" # calcul de alpha\n",
" x = coords_triangle[:, 0]\n",
" y = coords_triangle[:, 1]\n",
" mat_alpha = np.array([[x[1] - x[0], x[2] - x[0]], [y[1] - y[0], y[2] - y[0]]])\n",
" alpha = np.linalg.det(mat_alpha)\n",
"\n",
" b[triangle] += alpha / 6 * f(centre[0], centre[1])\n",
"\n",
" return b\n",
"\n",
"\n",
"def calcul_Ud(coords, dirichlet):\n",
" Ns = len(coords)\n",
" U = np.zeros(Ns)\n",
" # for d in dirichlet:\n",
" # x, y = coords[d].flatten()\n",
" # U[d] = u_d(x, y)\n",
" U[dirichlet.T] = u_d(coords[dirichlet, 0], coords[dirichlet, 1])\n",
"\n",
" return U\n",
"\n",
"\n",
"def tildage(A, b, coords, dirichlet):\n",
" A_tild = np.delete(A, dirichlet, 0)\n",
" A_tild = np.delete(A_tild, dirichlet, 1)\n",
" b_tild = np.delete(b, dirichlet, 0)\n",
" coords_tild = np.delete(coords, dirichlet, 0)\n",
"\n",
" return A_tild, b_tild, coords_tild\n",
"\n",
"\n",
"def untildage(x, dirichlet, U_d):\n",
" x_untild = np.zeros(U_d.shape[0])\n",
2022-03-16 14:44:07 +00:00
" not_dirichlet = np.setdiff1d(range(U_d.shape[0]), dirichlet)\n",
"\n",
" # print(x_untild)\n",
" # print(not_dirichlet)\n",
" # print(dirichlet)\n",
" # print(x)\n",
2022-03-15 07:04:08 +00:00
"\n",
" x_untild[dirichlet] = U_d[dirichlet]\n",
" x_untild[not_dirichlet] = x\n",
"\n",
" return x_untild\n",
"\n",
"n = 50\n",
"coords, elems3, dirichlet, neumann = maillage_carre(n)\n",
"\n",
"A = assemblage(coords, elems3)\n",
"b = second_membre(coords, elems3, f)\n",
"U_d = calcul_Ud(coords, dirichlet)\n",
"b -= np.dot(A, U_d)\n",
"\n",
"A_tild, b_tild, coords_tild = tildage(A, b, coords, dirichlet)\n",
"\n",
"x = np.linalg.solve(A_tild, b_tild)\n",
"x_untild = untildage(x, dirichlet, U_d)\n",
"\n",
"# show(coords_tild, x)\n",
"# print(coords.shape, x_untild.shape)\n",
"# print(coords, x_untild)\n",
"show(coords, x_untild)\n",
"show(coords, u_ex(coords[:, 0], coords[:, 1]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Partie II : maillage mixte et ajoût des conditions de Neumann**\n"
]
},
{
"cell_type": "code",
2022-03-16 14:44:07 +00:00
"execution_count": 132,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [],
"source": [
"e3 = np.array(\n",
" [[1, 2, 12], [2, 3, 12], [3, 4, 14], [4, 5, 14], [2, 15, 3], [3, 15, 4]]\n",
").astype(int)\n",
"\n",
"e4 = np.array(\n",
" [\n",
" [0, 1, 12, 11],\n",
" [11, 12, 13, 10],\n",
" [12, 3, 14, 13],\n",
" [10, 13, 8, 9],\n",
" [13, 14, 7, 8],\n",
" [14, 5, 6, 7],\n",
" ]\n",
").astype(int)\n",
"\n",
"dds = np.array([2, 15, 4, 6, 7, 8, 9, 10, 11, 0]).astype(int)\n",
"\n",
"nns = np.array([[4, 5], [5, 6], [0, 1], [1, 2]]).astype(int)\n",
"\n",
"ccs = np.array(\n",
" [\n",
" [0.0, 0.0],\n",
" [1 / 3, 0],\n",
" [0.53333333333333, 0.0],\n",
" [2 / 3, 1 / 3],\n",
" [1.0, 0.47],\n",
" [1, 2 / 3],\n",
" [1.0, 1.0],\n",
" [2 / 3, 1.0],\n",
" [1 / 3, 1.0],\n",
" [0.0, 1.0],\n",
" [0.0, 2 / 3],\n",
" [0.0, 1 / 3],\n",
" [1 / 3, 1 / 3],\n",
" [1 / 3, 2 / 3],\n",
" [2 / 3, 2 / 3],\n",
" [1.0, 0.0],\n",
" ]\n",
")"
]
},
{
"cell_type": "code",
2022-03-16 14:44:07 +00:00
"execution_count": 133,
2022-03-15 07:04:08 +00:00
"metadata": {},
2022-03-16 14:44:07 +00:00
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.66666667, -0.16666667, -0.33333333, -0.16666667],\n",
" [-0.16666667, 0.66666667, -0.16666667, -0.33333333],\n",
" [-0.33333333, -0.16666667, 0.66666667, -0.16666667],\n",
" [-0.16666667, -0.33333333, -0.16666667, 0.66666667]])"
]
},
"execution_count": 133,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def raideur_quadrangle(quadrangle):\n",
" x = quadrangle[:, 0]\n",
" y = quadrangle[:, 1]\n",
"\n",
" # calcul de la jacobienne\n",
" J_kk = np.array([[x[1] - x[0], x[3] - x[0]], [y[1] - y[0], y[3] - y[0]]])\n",
"\n",
" det_J_kk = np.linalg.det(J_kk)\n",
"\n",
" # on récupère les coefficients\n",
" coeffs = np.linalg.inv(np.matmul(J_kk.T, J_kk))\n",
" a = coeffs[0, 0]\n",
" b = coeffs[0, 1]\n",
" c = coeffs[1, 1]\n",
"\n",
" # on calcul M (on a calculé toutes les intégrales au préalable)\n",
" M = np.array(\n",
" [\n",
" [2 * a + 3 * b + 2 * c, -2 * a + c, -a - 3 * b - c, a - 2 * c],\n",
" [-2 * a + c, 2 * a - 3 * b + 2 * c, a - 2 * c, -a + 3 * b - c],\n",
" [-a - 3 * b - c, a - 2 * c, 2 * a + 3 * b + 2 * c, -2 * a + c],\n",
" [a - 2 * c, -a + 3 * b - c, -2 * a + c, 2 * a - 3 * b + 2 * c],\n",
" ]\n",
" )\n",
"\n",
" return det_J_kk / 6 * M\n",
"\n",
"raideur_quadrangle(ccs[e4[0]])"
]
2022-03-15 07:04:08 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Compléments d'analyse du système**\n"
]
},
{
"cell_type": "code",
2022-03-16 14:44:07 +00:00
"execution_count": 146,
2022-03-15 07:04:08 +00:00
"metadata": {},
2022-03-16 14:44:07 +00:00
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[-0.30555556 1.63888889 -1.67777778 2.11685185 -1.69356874 1.65116133\n",
" -0.30555556 -0.94444444 -0.61111111 -0.30555556 -0.61111111 -0.94444444\n",
" 1.31296296 1.77777778 1.31777778 -0.35296296]\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAADyCAYAAABkv9hQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABerElEQVR4nO2dd3xb9b3330eSJe+97XgkjrOcYSdOApSE0RIaSpLSAiktu6VwGS20t5S2Ty/tvaU83YMWnkvTlrJHCqSMUHaBQpb3iPdesi0P2dbWef5wpMiyxpEs2U6s9+sVSOyjc45sfc7v+/tOQRRFQoQIcXYjW+gbCBEiRPAJCT1EiCVASOghQiwBQkIPEWIJEBJ6iBBLgJDQQ4RYAii8fD8UewsRIvgIwb5AaEUPEWIJEBJ6iBBLgJDQQ4RYAoSEHiLEEiAk9BAhlgAhoYcIsQQICT1EiCVASOghQiwBQkIPEWIJEBJ6iBBLgJDQQ4RYAoSEHiLEEiAk9BAhlgAhoYcIsQQICT1EiCVASOgLgCiKGI1GzGYzoXbbIeYDb40nQgQYq9WK0WhEr9fbvyaXywkLC0OhUCCXyxGEoPchCLHEELysKKHlJkCIoojZbMZsNiMIAiaTyf51URSxWq12gRsMBmJiYlAqlSHhLw2C/gsOrejzgM1UdxSzDUEQEAQBmUxmP7alpYW8vDwiIyOB0IofYu6EhB5kzGYz3d3dWCwWsrKyEATBvoq7EqxN+HK5HLlcbl/tdTqd/XiFQmH/ExJ+CCmEhB4kHE11q9VqN9l9xdWKb7FYMJvN9mMUCoV9xZfJZCHhh5hFSOhBwGq1YjKZ7Ka6bRWXiqfjbeez4Sx8QRBmrPgh4YeAkNADik10NkebbRV2J1x35rsvuBK+2Wy234PNGlAoFCiVypDwlyghoQcIURQxmUxYLJZZ4nMWurdV3lcLwPm1zsLv7u4GICMjI7TiL1FCQg8Atti4bYV25VlfqMQYx/uxOfdMJtOMFd+2x5fL5SHhn6WEhD4HnGPjNlPdmUDu0eeKzaNvw5XwbY49hULh8sEV4swjJHQ/cY6NexKDP8INtNDd3Z8r4RuNRgwGAzDtZwgLC7Ov+CHhn5mEhO4HNoebO1PdGX9W9IXCk/AdHXuOpn6IxU/ot+QDNjO3oqICg8EgeT+70Kb7XM5lE77NcQdgNBo5duwYo6OjjI+PMzU1ZQ8nhlichFZ0iTjGxm2ON6kspDMukDg69fR6vX3lNxqNGI1GgNCKv0gJCd0LzrFxm/kaTKGfKQ8Gx3RdOF2g4yx8xzz9kPAXhpDQPeAuNi4Igk9m6pkiXF9xV6BjwyZ8g8Ewy7kXEv78EhK6GzzFxoO9Qp8tDwZvwhdFcYaZbwvnhQg8IaE74Wiqu4uNy2Sy0IruB66Eb7Va7U04ent7yc7ORqlUhirzAkxI6A5IjY0HYoUOVgrsmYTzz7i/v5/MzMxQ950gEBL6KbylsToSbGfcUsY5hu9cix8Svn8seaFLTWN1JNjOuNCDYRpXtfgh4fvHkha61WpFrVYjiiIJCQmSPyRnonDPBgFIEX6o+45rlqTQHR1uWq0WURRJTEyU/Pq5OuNGRkaora1FEAQSEhJISEggPj7ebrYupsy4xYy77js2R6oNpVKJSqVa0pV5S07ozqa6XC63J8NIxd8VXRRF2tvbUavVbNy4EZlMxtjYGBqNhtbWVuRyOQkJCej1emJiYnx9a0seV8Lv7OwkLCyM1NTUGSW5S60Wf0kJ3VWLJ19XZ/DPGWexWCgvLycyMpLS0lIsFgtWq5Xk5GSSk5OB6VTSkZERBgcHGR0dpbe3l8TERBISEoiKiloyH8pAYXvA2kz5pdx2a0kI3VNs3FfRgu/OOK1Wy9jYGBs2bCAtLQ2YroBzRqlUkpaWhk6nIzIykpiYGDQaDe3t7UxOThIdHW039SMiIny656WK1Wqd0dLLW9uts1X4Z73QPbV4At9Fa3uNlIeDzXTs7e0lOjraLnKpREREkJWVRVZWFqIoMjExwcjICI2NjRgMBmJjY+0rflhYmE/nXio4Ct0ZV8I/W7vvnNVClxIb99d09/Yas9lMTU0NYWFhlJSUUFFRIfn87hJsYmJiiImJIScnB6vVyvj4OBqNhq6uLkRRJD4+fpZj70whWA5DT0J3xpfuO2ea8M9KofsSG/dH6N5W9PHxcWpqasjPzycjIwOLxRLwD7JMJiM+Pp74+Hhg+sEyOjo6y7GXmJh4RnjdA9ER1xW+CN0ZT004ent7SUtLIzIy8oxou3XWCd2XFk/gv9Bd7bFFUaSnp4euri42bNhAdHS0/XhnsXm6rj/hNYVC4dKx19vby/DwMEqlErPZvGgde4tR6M44Cl+j0ZCWljaj+45txV+MtfhnldBtZpbthy7lgxMor7vZbKaurg5BENi6deuMlcCVcIO9ytoce2lpaXR1ddk/8IvVsRdIQc7HeS0Wi31sFpz+fdpq8e+44w5+8IMfsHr16oBf2x/OCqHbTPWTJ08SHx9Pamqq5NcGwnSfmJigurqaZcuWkZ2d7dO5pJw/EKhUKtLT07069uLj41EqlQG9thTOhBXd+bzOD3Ngxoq/GB6gNs54oTvGxucjJm57je06vb29tLe3s379erdJLovNTPbk2Ovu7sZqtXp07AXDGhFF8Yxa0cHz79VmNS0Wzlihuxp/JJfL/Vqd/d2j19TUYLFY2Lp1KwqF7z9KTxNV5zMF1hfHnu1hFuiHl6uR0oE6bzCE7u1nOjU1FRL6XHEXGw9WqMwZm9d1xYoVLFu2bNGt2HPFk2NPq9WiUqkwGAxMTEwEzLF3ppnu3jCZTAuyBXLHGSd0T7FxuVzu0hvuCV+F3t/fT1tbG4mJieTk5Ph0Lakshmo3Rxwde3DaJxFIx14wV95gnPdMe7ifMUKXEhuXyWQ+F6hIFbrVaqWhoQG9Xs/q1asZHh726TpnE+Hh4ahUKoqKigLm2AvWih4MvD2EF+N7OSOELjU2HizTfWpqiqqqKtLT01m9ejVjY2NBHVbgj99gPnGeDDsXx57jORdT3NkTttCaJxab2Be90H0Zf+SvM84TarWapqYm1q1bZ3dWBdu0XkwfEHd4etj64tiziTtYzrhg4E3oZrN50aUgL1qh+9PiSSaT+bxHd4fVaqWpqYmJiQlKS0tnmKD+WA7OePtQL6Y9ujO+rFZSHHuJiYmLThiesFgsHj+Pk5OTREVFzeMdeWdRCt1isTAxMYFSqfQpfzgQAgTQ6/VUVVWRnJxMSUmJy4q3YK/oZ4vQnXF27Ol0OjQaDf39/UxOTmIymRZVxp4rnJNlnFlsMXRYZEK3xcYNBgPl5eVs377dpw9UIIQ+NDREQ0MDa9ascdtearELcT4IlJltK8VVqVSMjY2RmpoasIy9YP2OvJnuoRXdA46xcVu2mq8fJn/Ca47Xb25uZnR0lC1btqBSqdweGyjLwR3BeJAEcv8bzMw4V469kZERnxx7NoKd5+6OiYmJ0IruCl96qnvCXwFarVaOHz9OQkICW7Zs8Xr9uQrR1pBiamrKvloFc48aaGEGw6PsSpSOjr38/HzMZrPLHnvOjj1P5wwE3vbotkSixcSCCl3K+CNf8Efow8PDTE1NsXr1arvDKBjXsWE2m6murkalUhEfH2//0CoUCpKSkkhISADOHmdcIM9p+xklJSUBpx17fX19NDQ02B17tlLchVrRF1v6Kyyg0H2tG5eCL6a7KIq0trYyPDxMZGSk/cMjBX9XdK1WS3V1NXl5eaSlpWEymUhJSQGm02o1Gg2dnZ2Mjo4SFhZm/+B62kYsBMF4CPkjSleOvZGREXvGXnh4OAaDAZ1OF1DHnhRnXGhFJ3CmujNSV1qj0Uh1dTXR0dFs2bKFI0eO+LRK+SN0k8lEdXW1vcrN+YGkUqnIyMggIyMDtVqNRqPBaDRSV1eH2WwmPj5+Xsx8qSzEiu6NiIgIIiIiyMzMRBRFhoaGaGtrm+HYs3n055KHbrFYPBYxLXmvu6+xcV9/+VKEPjIyQl1dHStXrrTXrdteJ3VF8cV0t8XjTSYT27dvl7Q6C4KASqUiNze
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def assemblage_quadrangle(coords, elems4):\n",
" Ns = len(coords)\n",
" A = np.zeros((Ns, Ns))\n",
" for quadrangle in elems4:\n",
" M = raideur_quadrangle(coords[quadrangle])\n",
" for i, a in enumerate(quadrangle):\n",
" for j, b in enumerate(quadrangle):\n",
" A[a, b] += M[i, j]\n",
" return A\n",
"\n",
"def second_membre_quadrangle(coords, elem4, f):\n",
" Ns = len(coords)\n",
" b = np.zeros(Ns)\n",
" for quadrangle in elem4:\n",
" coords_triangle = coords[quadrangle]\n",
" centre = np.mean(coords_triangle, 0)\n",
"\n",
" # calcul de alpha\n",
" x = coords_triangle[:, 0]\n",
" y = coords_triangle[:, 1]\n",
" mat_alpha = np.array([[x[1] - x[0], x[2] - x[0]], [y[1] - y[0], y[2] - y[0]]])\n",
" alpha = np.linalg.det(mat_alpha)\n",
"\n",
" b[quadrangle] += alpha / 4 * f(centre[0], centre[1])\n",
"\n",
" return b\n",
"\n",
"def condition_neumann(ccs, nns):\n",
" Ns = len(ccs)\n",
" coeffs = np.zeros(Ns)\n",
" for i, j in nns:\n",
" point1 = ccs[i]\n",
" point2 = ccs[j]\n",
" \n",
" valeur = np.linalg.norm(point1 - point2) / 2 * g((point1 + point2) / 2)\n",
" coeffs[i] += valeur\n",
" coeffs[j] += valeur\n",
"\n",
" return coeffs\n",
"\n",
"#------------------------------------------------------------------------------------------------------------------------\n",
"\n",
"A3 = assemblage(ccs, e3)\n",
"A4 = assemblage_quadrangle(ccs, e4)\n",
"\n",
"b3 = second_membre(ccs, e3, f)\n",
"b4 = second_membre_quadrangle(ccs, e4, f)\n",
"\n",
"A = A3 + A4\n",
"b = b3 + b4\n",
"\n",
"# print(A)\n",
"\n",
"U_d = calcul_Ud(ccs, dds)\n",
"b -= np.dot(A, U_d)\n",
"b += condition_neumann(ccs, nns)\n",
"print(b)\n",
"\n",
"# print(b3)\n",
"# print(b4)\n",
"\n",
"# A_tild, b_tild, ccs_tild = tildage_quadrangle(A, b, ccs, dds, nns)\n",
"A_tild, b_tild, ccs_tild = tildage(A, b, ccs, dds)\n",
"\n",
"x = np.linalg.solve(A_tild, b_tild)\n",
"# print(x)\n",
"x_untild = untildage(x, dds, U_d)\n",
"\n",
"# print(A)\n",
"# print(b)\n",
"# print(U_d)\n",
"show(ccs, x_untild)\n",
"# show(ccs, u_ex(ccs[:, 0], ccs[:, 1]))"
]
2022-03-15 07:04:08 +00:00
}
],
"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.10.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
