2022-03-15 07:04:08 +00:00
{
"cells": [
{
"cell_type": "code",
2022-03-16 14:44:20 +00:00
"execution_count": 147,
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:20 +00:00
"execution_count": 148,
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:20 +00:00
"execution_count": 149,
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:20 +00:00
"execution_count": 150,
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:20 +00:00
"execution_count": 151,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [
{
"data": {
2022-03-16 14:44:20 +00:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPYAAADyCAYAAAB+gzjeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAB8OUlEQVR4nO29d3hcZ5n+/zln+qj3almyZcu922kLIbRkAymUkCyQQqgBwpJAQkL4sYFdyhdIWMouWSAQajokgYQkEEiANHdLtmxLVi+WRm2k6eWc9/fHaEZTzhTJki3bc19XrlgzZ06ZOfd5nvcp9yMJIcgiiyzOLMin+gSyyCKL+UeW2FlkcQYiS+wssjgDkSV2FlmcgcgSO4sszkBkiZ1FFmcg9Gnez+bCsshi4SHN9w6zFjuLLM5AZImdRRZnILLEziKLMxBZYmeRxRmILLGzyOIMRJbYWWRxBiJL7CyyOAORJXYWWZyByBI7iyzOQGSJnUUWZyCyxM4iizMQWWJnkcUZiCyxs8jiDESW2FlkcQYiS+wssjgDkSX2KYAQAr/fTzAYJCv/nMVCIJ3QQhbzDFVV8fv9eL3eyGs6nQ6DwYBer0en0yFJ8953n8VZBimNxciak3mCEIJgMEgwGESSJAKBQOR1IQSqqkYI7fP5yMvLw2g0Zol+dmDef+CsxT4JCLve0eQNQ5IkJElCluXIth0dHdTX12O1WoGsRc9i9sgSe4ERDAbp7+9HURRqamqQJClipbUIGia6TqdDp9NFrLnH44lsr9frI/9liZ6FFrLEXiBEu96qqkZc8NlCy6IrikIwGIxso9frIxZdluUs0bPIEnshoKoqgUAg4nqHrXSmSLV9eH9hxBNdkqQYi54l+tmJLLHnEWGShQNjYSubjKjJ3PHZQIvowWAwcg5ha6/X6zEajVminyXIEnueIIQgEAigKEoC2eKJnc6Kz9bCx382nuj9/f0AVFVVZS36WYIssecB4dx02AJrRb5PVSFK9PmEg3GBQCDGoofX6DqdLkv0MwRZYp8A4nPTYdc7HvO5xj5RhCPuYWgRPRyI0+v1mg+qLBY/ssSeI+Jz06lu/rkQdb6Jnez8tIju9/vx+XxAKE5gMBgiFj1L9NMDWWLPAeEAWTLXOx5zsdinCqmIHh2Ii3bds1h8yP4qs0DYbd2/fz8+ny/j9eipdsVPZF9hoocDbQB+v59du3Zht9uZmprC7XZH0ntZLA5kLXaGiM5NhwNlmeJUBs/mE9FBOK/XG7Hsfr8fv98PkLXoiwRZYqdBfG467I4uJLFPlwdBdPkrzDS0xBM9us49S/STgyyxUyBZblqSpFm5nacLUWeLZA0tYYSJ7vP5EoJxWaIvLLLEToJUuemFtsBnyoMgHdGFEDFuezi9lsWJI0vsOES73sly07IsZy32HKBFdFVVI6ITg4OD1NbWYjQas51rJ4gssaOQaW56PizwQpWUnk6I/46Hhoaorq7OqsvMA7LEnka6stBozEfw7Gwg7lwQn0OP70XPEj0znPXEzrQsNBonGjwbHR2ltbUVnU5HQUEBRUVFFBUVodfrNbc/W6HVi54lemY4q4mtqio2mw0hBEVFRRnfFHN1xYUQHDt2DLvdzpYtW5BlmampKSYmJujp6UGSJAoLC/F6vfNe7HEm3PCZED2rLhPCWUns6ACZw+FACEFxcXHGn59L8CwYDLJ7926KiorYtm1bRFmluLg4cuxAIIDdbmd8fJxjx47R399PUVERxcXF5OXlzfkmPVOtfzJ1mXDgMwyj0YjJZDqrOtfOOmLHu946nS5SfJIpZmux7XY7drudzZs3U1pamnQ7g8FAWVkZDoeD/Px88vLyGB8fp7+/H4fDgcViobi4mKKiIqxW61lzk2YKLaL39vZiMBgoLy+PaVE903vRzypia0kWzdb6QubBs7Di6OjoKPn5+SlJHY3wg8NkMlFVVUVVVRVCCDweD+Pj43R2duJ2u8nLy4usz81m86yu4WxA+HsMu+Znk4zUWUHsVLnp2Ua4IbPgmd/vp7m5mfz8fDZt2kRLS8us9q+VHrNarVitVmpraxFC4HA4mJiY4MiRI/j9/phAnMFgmNU1nalQVTVGoiqdjNSZQvQzntipJItg9hHu8GdSPQzGx8c5fPgwK1eupKysLNLiOZ+QJIn8/Hzy8/NZunQpqqoyOTnJxMQEfX19CCEoLCykqKjorO66iiZ2PLSIfqaoy5zRxM4kNz1XV1zrM0IIurq6GBkZYcuWLVgsFuDklJTKshyx1hDSMw8H4mw2W6QjKxyIW2w12gsV4EtF7HjMRl1msRP9jCT2bHLTcyG2FvH8fj8tLS3k5OSwffv2mGOeiry0Xq+ntLSU0tLSyPrbaDQyODiIw+HAZDJFAnE5OTmn/AadD8VWLcyG2PFIJToxODhIRUUFVqt1UcpInXHEno1kEcyd2IqiRP6emJigtbWVFStWUF5errl9PLFTHXchHgR6vZ6KigoqKioAIoG47u5uXC4XOTk5EaKHPY2TicVI7HhEE318fJyKiooYdZmwRV8MvehnFLHDblP4S87kRjmRqLgQgu7uboaHh9m8eXNk1lY8FmNJqcVioaamhpqaGoQQuFwuxsfHaWtrw+fzkZ+fH3HtjUbjgp/PfBLwZOxXUZTIGCaY+T3Dveif/vSn+dKXvsSqVavm/diZ4Iwgdtj1PnLkCIWFhZpWMxnmarGDwSD79u3DYrGwY8eOeb15FsJip3rISZJEbm4uubm51NXVoapqpCJuYGAARVEigbjCwsJI6et84nSw2PH7jXbTo0tcIWTRT4XnE8ZpT+zo3PRC5qSj4fF46OvrY82aNVRWVqbdfrGsuzKFLMsUFhZSWFhIQ0MDiqJgt9uZmJigu7sbSZLw+XxMTExQUFAwL8QRQpxWFhtS/64ul4vc3NwFOW4mOG2JrTVOR6fTzcn6ZvqZcCXTwMAAlZWVGZE63f60bo75ttgnui+dTkdJSQklJSVAaH733r17sdlstLe3YzQaI6Wvubm5c3qQaY0Yng8sFLHTfadutztL7NkiWW56PlNX8QgEAhw8eBCTycTKlSux2+1zOfUzAuEAUVNTEwBer5eJiQl6e3txOp1YrdaYQFwmhD3dXPF0CAQCJyU2kQynHbFT5aZ1Ol1MtDoTZELsqakpDh48SENDA1VVVYyNjS1o8Guxt23Gk9BsNseUvrrd7kgji9frJTc3N0J0k8mkuc+FtKwLsd/Fvrw6bYidSW5aluVZN3SkIrYQgr6+PgYGBti4cSM5OTnA4ifeQiOVdZUkiZycHHJycliyZAmqquJ0OhkfH6e1tZVgMKjZg75QFnshkO63XwzXcloQO9Pc9Hy64sFgkEOHDqHT6dixY0dMBHQux5kN5lLmejIxm4eaLMuR0tf6+noURYmUvvb29gJQWFiIyWQ65WTIFOFUVyqcanIvemLPZpzOXINn8XA4HLS0tFBfX091dbXmZxbaFV/smOs56nQ6zR70oaEh7HY7Ho9nXnrQFxLpiB0MBtMSf6GxaIk9F8kiWZZnvcaOP+bAwAC9vb1s2LAhaVRzPix2uht2rg8Oj8eLxbKwLZzzaY3CPeh6vR6z2UxdXR0TExMxPehhoi+WHnRFUVLej+FKvlOJRUlsRVFwOp0YjcZZ1d+eCOGCwSCtra0A7NixI2URxsmw2LPd/6TLy03/+wwHu4fIMelYUZ7HxVsaefeFm+b9/BbCzQwHz0wmUySVGO5Bn5iYiPSgRwfiTlUPenxxSjxOdQ4bFhmxw7lpn8/Hvn37OPfcc2d1A82V2IqisHPnTurq6qitrU27/WILnv3mpUN86/GX8QRUQMbr8DPmmeK1nr189Xe7ydcFWL+kiA++fQcXbW6aF1LON7G1HhbRPejh0len05lxD/pC/UbpXPGsxY5CdG46XA0225tnLumuwcFBPB4P5513Hnl5eRl95kRdca/Xy4EDByI3ZtgChW/MTB8ctkkXH7/vLxzsHECKcg2FJMH09yckmUnVxD973PzzJy8iKX+mOt/AG9bU8sG3bmdlXcWsz38hCJPJ7y1JEnl5eeTl5UV
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:20 +00:00
"execution_count": 152,
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:20 +00:00
"execution_count": 153,
2022-03-15 07:04:08 +00:00
"metadata": {},
"outputs": [
{
2022-03-15 07:05:49 +00:00
"data": {
2022-03-16 14:44:20 +00:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPcAAADyCAYAAACRQVPgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAADDBklEQVR4nOz9eZRdWV7fiX723me8U8yDpFBISs1K5ZwSVbYfuO2Fn/08NWv5YeN+xnTZbUwDnvADm3bTLrDBmMk2NGAoJmO7sB9gJkM1ZqgyDZWVykwplZrHkBQRCsV4405n3Hu/P07cIKTUEKEMKYeK71q5MjPinnNunLO/5/fbv+H7E9ZatrCFLXz0IN/vL7CFLWzhyWCL3FvYwkcUW+TewhY+otgi9xa28BHFFrm3sIWPKLbIvYUtfEThPOL3W3myLWzhyUM8iZNuWe4tbOEjii1yb2ELH1FskXsLW/iIYovcW9jCRxRb5N7CFj6i2CL3FrbwEcUWubewhY8otsi9hS18RLFF7i1s4SOKLXJvYQsfUWyRewtb+Ihii9xb2MJHFFvk3sIWPqLYIvcWtvARxRa5t7CFjyi2yP0+wFpLmqbkec6WtPQWnhQeJdawhU2GMYY0TYnjePVnSilc18VxHJRSCPFEeve38EUG8QjLsWVWNgnWWvI8J89zhBBkWbb6c2stxphVUidJQrVaxfO8LbJ/ceCJPOAty/0U0HXD1xK4CyEEQgiklKufvXr1Krt376ZUKgFbln0Lj4ctcj9h5HnO5OQkWmt27NiBEGLVWt+PpF2yK6VQSq1a9SiKVj/vOM7qP1tk38KDsEXuJ4S1brgxZtUd3yjuZ9m11uR5vvoZx3FWLbuUcovsWwC2yP1EYIwhy7JVN7xrrdeLh32+e74u7iW7EOIuy75F9i9ebJF7E9ElWjdY1rW2DyLrg1zzjeB+ZM/zfPU7dK2+4zh4nrdF9i8ibJF7k2CtJcsytNbvIty95H6UNd+opb/32HvJPjk5CcC2bdu2LPsXEbbIvQno5q67lvh+EfH3q1hl7ffpBuiyLLvLsnf37EqpLbJ/hLBF7veAe3PXXTf8Xmzmnvu9ohuJ7+J+ZO8G5xzHue/LagsfDmyR+zFxb+76YQR4HLJuNrkf9P3uR/Y0TUmSBCjiBq7rrlr2LbJ/eLBF7sdAN2j2IDf8XjyO5X6/8DCyrw3OrXXjt/DBxNaT2QC6LuypU6dIkmTd+9P32y1/L+fqkr0bfANI05QTJ05Qr9dpNBp0Op3V1N8WPjjYstzrxNrcdTd4tl68nwG1zcTawFwcx6sWPk1T0jQF2LLsHyBskfsRuDd33XVNnyS5Pywvg7WlsvCHTTD3kn1tXfwW2Z8etsj9EDwody2E2JAL+mEh60bxoCaYLrpkT5LkXQG6LbI/eWyR+wF4WO76SVvij8rL4FFkt9be5cJ3U29b2BxskfserHXDH5S7llJuWe7HwP3IboxZFa6Ynp5mbGwMz/O2Ot42AVvkXoP15q7fqyU2xjA5OYnruvT19eF53ns6/4cV997jmZkZtm/fvqVSs0nYIvcKHlVCuhbvJaAWRRGnT5+mt7eXNE2ZmprCGENPTw99fX309va+1z/lQ417c+z39rJvkX39+KIn93pLSNficQNqc3NzXLp0iSNHjlCpVFY9BK019XqdpaUlJiYmVvekQgh6enq+aINO9+tl3yL7+vFFTW5jDLOzs1hr6evrW/fCeBy3udFokOc5x44dw/O81dQaFAt0YGCAgYEBAC5duoRSitnZWa5cubLqvvf19VGtVh9b9OHDjvWQfUul5g/xRUnutUGzZrOJtZb+/v51H7+RgFocx5w+fRohBK+88sq6FpvjONRqNQYHB1fPsbS0xOTkJK1WizAMV8leKpUeec6P6v79QSo13WBoF57n4fv+F13H2xcdue91w5VSd1nR9WC9lrvrhu/du5epqanHXlhBELBt2za2bduGtZZOp8PS0hLXrl0jiiIqlQp9fX309/fj+/5jXeOjgPuR/ebNm7iuy/Dw8F3trV8MvexfVOS+n/zRRtNa8OiAmjGGK1eu0Gg0ePXVVwFWBRPWg0cJOZTLZcrlMmNjY1hraTabLC0tcf78ebIso1arrVp213U39Ld9lNC9j103/YtNkuqLgtwPy11vNPINDw+odd3w/v7+VTf8SdaiCyGo1WrUajV27dqFMYbl5WWWlpa4desW1lqUUlQqFbTWd0WjvxhgjLlL7upRklQfJbJ/5Mn9MPkj2Hjku3vM/cg3Pz/PxYsXOXTo0Gpw7GGffxKQUq5abSikla9cuUK73eatt95CKbXqwler1Y98JH4tue/F/cj+UVKp+UiTez2568d1y9ceY63lypUr1Ot1Xn311Xfte9/P8lPHcSiXy/T09LBt2zbSNGVpaYnbt29z8eJFfN9ffRlUKpX3bfE+qZffw8h9LzaiUvNhIPtHktwbyV0/DrnXki9JktWilFdfffW+D/uDVHHmeR4jIyOMjIwARVHN0tISN2/epNVqUS6XV8kehuFTW7yboQR7P2yE3PfiYcIV09PTjIyMUCqVPrCSVB85cm9E/ggen9xaaxYWFrhw4QIHDx5cTVs96PP3kvth132aL4MwDAnDkO3bt2Otpd1us7S0xJUrV4jjmGq1ukr2JxmJ/yCS+16sJfvi4iIjIyN3qdR0LfsHpZf9I0XurgvVvdHrWSyPS+65uTlmZ2d55ZVXCILgkZ+/l6wfFEu+FkIIKpUKlUqFnTt3YoxZjcSfO3eOPM9Xy2Q3+/tvJgmfxnm7wcm1vezwh8IV3/AN38A/+Sf/hEOHDm36tdeLjwS5u274hQsX6O3tZXh4eN3HbpTcSZIwMTGB67q8+uqrT2ThPAnL/ThWUUpJT08PPT097N69G631aiS+0+nw5ptv0tvbS19fHz09Pe8pEv9hsNz3nnft37u2HBYKyx6G4aZfdyP40JN7be76SeSs12JxcZHz58+vFkRsJFDzUYBSiv7+fvr6+lhcXOT555+nXq8zPz/P1atXcRznrkj8Rv5ua+2HynLDw59ru92mUqk8keuuFx9act9vdI9S6rFc7EcdY63l2rVrLCws8Morr9BoNFheXn7s7772vE8jALfZXkD3e7uuy9DQEENDQ0Dh1SwtLTE1NUWz2SQIgtX9erlcfigZ7jfeeDPwpMj9qHva6XS2yP04eFDuejPSWvciTVNOnz5NtVpddcObzeaW0ud94Ps+o6OjjI6OYq1djcRPTEysWrK1kfi1+LC55Y9ClmXv6tN/2vjQkfthuWulFFrrDZ3vYeTuBpIOHDiwap26x7wXa2itZWpqiizLGBgYeFfzxwcpdXY/rIeIQghKpRKlUokdO3ZgraXVarG0tMSlS5dIkuSuMtknaWGfVFzkg44PDbnXk7uWUm64CeR+5LbWcv36debm5nj55ZffZWXeC/nyPOfMmTM4jkOpVLqr+aO/v39D3WnvFx7HygohqFarVKtVxsfHMcbQaDRW3fjuc5ufn6e3txfH+eAuzUc9+yflhWwUH9w7uAbrzV1vhluepinvvPMO5XKZY8eObYqGWhetVovTp0+ze/duhoeHyfP8ruaPxcVFzpw5QxzHhGFIEAT09vZ+4OrBN8OrkFLS29tLb28ve/bsWS2kqdfrTExMIIRYteofNMGK9dTofxAI/oEn90ZG9zxuQK2Lrhu+f//+h6bTHsdyZ1nG6dOnee6556hWq6udSd3zdZs/du/ezdTUFK1Wi8XFRa5du/aeotBPCpv9HYQQhGHIvn37gOJ+LS0tvUuwor+//30tk4VHkzvP8w/EC/kDS+7HkT+SUm54z9291vXr15mdnb2vG36/66z3JWKM4eLFi2RZxh/9o3901d182OKUUlIqldi5cydQRKEXFxeZnJyk2Wyuloj29/e/L7nUJ2GV7j1ntwe7+5LtClbcunXrLsGK7j14mmTXWj90Pbbbbcrl8lP7Pg/CB5LcWmtarRae522oXvdx3OU0TYmiiDiOH+iG34v1Wu44jnn77bcZHh4mDMN17yPvPb/v+3eJNXRLRLuBqW7V2NP
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:20 +00:00
"execution_count": 154,
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:20 +00:00
"execution_count": 155,
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]])"
]
},
2022-03-16 14:44:20 +00:00
"execution_count": 155,
2022-03-16 14:44:07 +00:00
"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:20 +00:00
"execution_count": 156,
2022-03-15 07:04:08 +00:00
"metadata": {},
2022-03-16 14:44:07 +00:00
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
2022-03-16 14:44:20 +00:00
"[-0.19625549 1.80841983 -1.29709571 3.27946985 -1.31278846 1.81960843\n",
" -0.19625549 -0.58876648 -0.25543315 -0.19625549 -0.25543315 -0.58876648\n",
" 2.4657089 2.90036722 2.46946706 -0.14351276]\n"
2022-03-16 14:44:07 +00:00
]
},
{
"data": {
2022-03-16 14:44:20 +00:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAADyCAYAAABgSghtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABirklEQVR4nO2dd3hb9dn+P0eSJe894tiJneVMx3ESJyEUwoYGCNCyWyDQlsILdNBBaftr6QLeTtrSwlugpRQIIwHCTBmFQIEkTuLEK7EdO96x5W1Z1jzn/P5wjpBtjSNZsp1E93X5gkhnSTr3eZ7vM+5HkGWZCCKI4OSEZqovIIIIIggfIgSPIIKTGBGCRxDBSYwIwSOI4CRGhOARRHASI0LwCCI4iaHz834khxZBBOGHEK4DRyx4BBGcxIgQPIIITmJECB5BBCcxIgSPIIKTGBGCRxDBSYwIwSOI4CRGhOARRHASI0LwCCI4iREheAQRnMSIEDyCCE5iRAgeQQQnMSIEjyCCkxgRgkcQwUmMCMEjiOAkRoTgEURwEiNC8CmALMvY7XacTicR2eoIwgl/gg8RhBiSJGG327Fara7XtFotUVFR6HQ6tFotghC2/v8ITjEIfixIxLyECLIs43Q6cTqdCIKAw+FwvS7LMpIkuYhts9lISEhAr9dHCH9qIGw/cMSCTwIUl9ydxAoEQUAQBDQajWvb+vp68vPziY2NBSIWPoLgESF4mOF0OmltbUUURXJychAEwWW1PRFVIbxWq0Wr1bqsu8VicW2v0+lcfxHCR+ALEYKHCe4uuSRJLtc8UHiy8KIo4nQ6XdvodDqXhddoNBHCR+BChOBhgCRJOBwOl0uuWG0FZrOZuro6YmNjSUlJITEx0UVgYNz27lCOp2As4QVBGGXhI4Q/tREheAihkE0JoCmkdSdse3s7jY2NzJ07F7vdTkdHB7W1tej1elJSUkhJSQkodeaJ8E6n03UNivXX6XTo9foI4U8xRKLoIYIsyzgcDkRRHEe6jo4OTCYTFosFgCVLliBJ0qh1uNVqpa+vj76+Prq7u4mPjycjI4OUlBTi4uKCJqUsy7S1tQGQnZ0dsfDTE5Eo+nSGkttWCDuWMBaLhdbWVgoKCsjJyQHAbreP2iY6Oprs7Gyys7OpqqoiKysLu91OU1MTZrOZmJgYl4WPjY1VTUr361GCdg6HY5SFV9bwWq02QviTDBGCTwBjc9vu62jl/ZaWFlpaWsjMzHSR2x80Gg0Gg4H09HRmzpyJLMsMDw/T19dHQ0MDw8PDxMfHuwgfHR0dEOG1Wu2oaxxLeCVgp9PpPD6wIjhxECF4kBib2x5LAofDQWVlJQaDgYULFzIwMBDw8RUIgkBcXBxxcXHk5uYiyzJms5m+vj7q6uqwWq3jCD8W3kjqifB2ux2bzQaMPGyioqJcFj5C+BMLEYIHASWQ5s0l7+vro7q6mnnz5jFjxgy6u7sDDpz5ez8+Pp74+HhmzZqFLMuYTCb6+vo4fPgwdrudxMREF+EDgS/Cuwfs3F36CKYvIgQPAIpLXl5ezvz588dZSlmWOXr0KF1dXRQXF7sq0XylvTwhmO0TExNJTEwkLy8PSZIwmUz09vbS3t7O8PAwMTEx6HQ6UlJSiIqKCujYCuGVa7Lb7Rw4cIBFixa5XPkI4acnIgRXCffcthJQc4fNZqOiooKEhARKSkpU57XDAY1GQ1JSEklJSQA0Nzdjs9kYGhqitbUVSZJISkoiJSWF5ORkdDp1t4F7sM5qtbqIb7fbXUHDiIWfXogQ3A/G5rYVN9WdsD09PRw+fJiCggIyMjLGHSPcFtwfNBoNcXFxzJw5ExhZYvT399PX10dTUxMAycnJrj93F93fdY618IpL70549zr6COEnFxGC+4C33LYgCEiShCRJHDlyhIGBAVatWuUxuKVsP536vrVaLWlpaaSlpQEj9fL9/f309vZy9OhRNBoNycnJrio7b4T31jijQCG8zWYbF7SLEH5yECG4F/jKbQuCgNVqpaqqivT0dFavXu0zMDbVFtwfdDod6enppKenAyMud39/P0ajkSNHjrjW7ikpKSQkJKgmpT/Cy7I8yp1X0nIRhA4Rgo+Bu0vuKbcNI+vt6upqCgsLVUWpp5sF9we9Xk9mZiaZmZnAyOft6+vj2LFj1NTUYDAYsNvtDA4OkpCQEFTRDeDqlFPEL9rb28nNzUWv10c65UKECMHd4C+3LUkShw8fZnh4mKVLl6pOQXkiuL+Gkun0QDAYDMyYMYMZM2YAI2W1+/fvp7W1laGhIaKjo10WPpCyWk8lvTNnzoyo3YQQEYIfh79yU7PZTHl5uStQpTYQBdOPsBNFdHQ0UVFRLFmyBFmWsVgsroDd0NAQcXFxrjV8IGW1wLgc/Nhe+AjhA8MpT3B/5abwWQfY0qVLSUpK4vDhw0iSpPoc030NPhEIgkBsbCyxsbHk5OT4LauNiYkJ6Nhje+EjhA8MpzTBJUnCaDQiyzIpKSnjbg6n00l1dTUAa9asceWLT0TCTtaN76msdmhoiL6+Pmpra116c4qF95Z58HZsf4SPqN2MxilJcPdAmslkQpZlUlNTR20zODhIZWUleXl5zJw5c9SNotFoTigLPpUPF0EQSEhIICEhgdmzZyNJkovwY8tqA/lOlWN7UrtRAqQK9Ho9BoPhlOyUO+UIPtYl12q1riIW5f2Wlhba2tpYvnw58fHx444x1YQ9kaHRaMaV1Q4ODtLX14fVaqW0tHRUlV2gZbVjCd/c3ExUVBSZmZmjWmNPlV74U4rgnqSU3K2x0gGm1+tZs2aN10Da2Eo2f4g8ELxDKapJTk6mu7ublStXMjAwQH9/P83NzciyHFRZLXz2PSou+6kob3VKENxXblsh69gOMF9QKtnUwhNhh4aG0Gg0HoNOpxLBx0Kr1ZKamupaMjmdTgYGBujr66OxsRFBEFzr96SkJL/ZDEmSRkln+ZO3OtkIf9IT3JeUkoKenh66u7tHdYD5wkQssizLNDQ0YDQa0Wq12O12l4UKtNPrVIBOpxtVVutwOOjv76enp4f6+nq0Wu0owo/NgrgTfCw8Ef5kU7s5qQnuL7dts9mor69HEATWrl2rugQz2CCb3W6noqKC+Ph4Vq5cCYzcVIqFamlpQZZlNBoNoiiSlpYWUL79RIbaB2ZUVBQZGRmuph73stq6ujqioqJGldX6IvhYBKJ2c6IQ/qQkuJrcttIBlp2djc1mC6jpIRgL7nA4KC0tZcGCBWRmZrpiARqNZpQwg9PppLa2FpPJxP79+9HpdKSmprpu2Ol+QwULb4Mg/MFbWW17ezsmk2mUwm2g358v8Yv29naysrKIjY2d1vJWJx3B1ZSbuneAWa1Wl+qoWgiCgCiKqq+npaUFi8XC6aef7loCePMCdDod8fHx6PV6ZsyYgc1mo7e3l9bWVkwmE7GxsS7Cx8TETLsbKlgES/CxGFtWW1FRQVRUlOv7cxevDFSt1p3wvb29LmFMRe1GsfDTqRf+pCK44k4pX/bYH89isVBeXj6qA0x5GAQCtVF0p9NJZWUlUVFRrmovBWo9AIPB4FJbda8SO3LkCFarlYSEBFJSUkhNTUWv1wf0OaYTAnGlA4EgCC5L615W29jYiNlsJi4ublSVnVrCi6LoGi8Fo9Vu7HY7d955Jz/+8Y9ZtGhRyD9TIDgpCK645IcPHyY5Odnlrrmjs7OTI0eOsGTJklFNIoGup0Gdi24ymaioqCA/P5/s7Gw+/fTTCR9/bJWYIs3U19dHZWUloii6Ak5jU0rT3dKHyoKPxdgo+tiyWrPZTH9/P/X19QGV1UqSNMp9dy+dhRELH0hZbrhwwhPcPbftiaxKB5jVaqWkpGSclQs0p63s4+uh0NbWRlNTk9dCmVDBXZopPz9/lFJLY2Oja33vcDgCyh/7QzhSeEpwMdTwF0VXxCt9ldUqhDcYDOP29waz2RzW314tTliCexoTpNVqRxHPvQNs8eLFHn+QQHPayj6ebnJRFDl06BCSJI2qXff1GbxdUzAkGqvUokSYm5ub6e7uxmg0utbvE52WEmpr62m0cqiOG0gUfWx
2022-03-16 14:44:07 +00:00
"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
}
