diff --git a/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.md b/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.md index 1f71d476c4d..d32816ba60c 100644 --- a/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.md +++ b/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.md @@ -62,7 +62,7 @@ import matplotlib.pyplot as plt from matplotlib import cm # Some plot settings -plt.style.use("seaborn-deep") +plt.style.use("seaborn-v0_8-deep") plt.rcParams["lines.linewidth"] = 2.0 plt.rcParams["lines.color"] = "black" plt.rcParams["legend.frameon"] = True diff --git a/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole_convergence_analysis.ipynb b/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole_convergence_analysis.ipynb index f4cb324bfde..430a78d9518 100644 --- a/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole_convergence_analysis.ipynb +++ b/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole_convergence_analysis.ipynb @@ -34,7 +34,7 @@ "id": "4d537496-f74d-43de-8dff-b917e644a188", "metadata": {}, "source": [ - "For many typical engineering problems that need to be solved using numerical techniques, convergence plays a decisive role. It is important to obtain an accurate numerical solution, i.e., to know the mesh refinement level beyond which the results of the numerical analysis do not change significantly, anymore. \n", + "For many typical engineering problems that need to be solved using numerical techniques, convergence plays a decisive role. It is important to obtain an accurate numerical solution, i.e., to know the mesh refinement level beyond which the results of the numerical analysis do not change significantly, anymore.\n", "\n", "In this study the convergence in stress and displacement values for the disc with hole problem are analysed and compared to each other. This study is based on the accomanying Jupyter notebook \"Linear_Disc_with_hole.ipynb\", which explains the main features of the analytical and numerical solution of the disc with hole problem. For a better access to the solutions used in this document, it is recommended to read through the previously mentioned Jupyter notebook first." ] @@ -84,20 +84,20 @@ "import matplotlib.pyplot as plt\n", "from matplotlib import cm\n", "\n", - "#Some plot settings\n", - "plt.style.use('seaborn-deep')\n", - "plt.rcParams['lines.linewidth']= 2.0\n", - "plt.rcParams['lines.color']= 'black'\n", - "plt.rcParams['legend.frameon']=True\n", - "plt.rcParams['font.family'] = 'serif'\n", - "plt.rcParams['legend.fontsize']=14\n", - "plt.rcParams['font.size'] = 14\n", - "plt.rcParams['axes.spines.right'] = False\n", - "plt.rcParams['axes.spines.top'] = False\n", - "plt.rcParams['axes.spines.left'] = True\n", - "plt.rcParams['axes.spines.bottom'] = True\n", - "plt.rcParams['axes.axisbelow'] = True\n", - "plt.rcParams['figure.figsize'] = (8, 6)" + "# Some plot settings\n", + "plt.style.use(\"seaborn-v0_8-deep\")\n", + "plt.rcParams[\"lines.linewidth\"] = 2.0\n", + "plt.rcParams[\"lines.color\"] = \"black\"\n", + "plt.rcParams[\"legend.frameon\"] = True\n", + "plt.rcParams[\"font.family\"] = \"serif\"\n", + "plt.rcParams[\"legend.fontsize\"] = 14\n", + "plt.rcParams[\"font.size\"] = 14\n", + "plt.rcParams[\"axes.spines.right\"] = False\n", + "plt.rcParams[\"axes.spines.top\"] = False\n", + "plt.rcParams[\"axes.spines.left\"] = True\n", + "plt.rcParams[\"axes.spines.bottom\"] = True\n", + "plt.rcParams[\"axes.axisbelow\"] = True\n", + "plt.rcParams[\"figure.figsize\"] = (8, 6)" ] }, { @@ -112,10 +112,11 @@ "outputs": [], "source": [ "import os\n", + "\n", "# All files produced by this notebook will be put there.\n", "# ATTENTION: We assume that this notebook is executed in the directory where\n", "# it is stored. Otherwise this notebook might not work!\n", - "out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', 'out')\n", + "out_dir = os.environ.get(\"OGS_TESTRUNNER_OUT_DIR\", \"out\")\n", "os.makedirs(out_dir, exist_ok=True)" ] }, @@ -133,7 +134,7 @@ "source": [ "# This list contains all parameters for the\n", "# different meshes we create\n", - "STUDY_indices = [ 8, 16, 24, 40, 60, 80, 240 ]\n", + "STUDY_indices = [8, 16, 24, 40, 60, 80, 240]\n", "\n", "# With this parameter the length of one axis of the square plate is defined\n", "STUDY_mesh_size = 20" @@ -156,109 +157,119 @@ " # Therefore it is tracked in git and might be in a different\n", " # directory.\n", " d = \"out\" if study_idx == 240 else out_dir\n", - " \n", + "\n", " reader = pv.PVDReader(f\"{d}/disc_with_hole_idx_is_{study_idx}.pvd\")\n", - " reader.set_active_time_point(-1) # go to last timestep\n", - " \n", + " reader.set_active_time_point(-1) # go to last timestep\n", + "\n", " mesh = reader.read()[0]\n", - " \n", + "\n", " return mesh\n", "\n", + "\n", "def slice_along_line(mesh, start_point, end_point):\n", - " line = pv.Line(start_point, end_point, resolution = 2)\n", + " line = pv.Line(start_point, end_point, resolution=2)\n", " return mesh.slice_along_line(line)\n", "\n", + "\n", "def get_sigma_polar_components(mesh):\n", " sig = mesh.point_data[\"sigma\"]\n", - " \n", + "\n", " xs = mesh.points[:, 0]\n", " ys = mesh.points[:, 1]\n", " sigs_polar = vec4_to_mat3x3polar(sig, xs, ys)\n", - " \n", - " sig_rr = sigs_polar[:,0,0]\n", - " sig_tt = sigs_polar[:,1,1]\n", - " sig_rt = sigs_polar[:,0,1]\n", - " \n", + "\n", + " sig_rr = sigs_polar[:, 0, 0]\n", + " sig_tt = sigs_polar[:, 1, 1]\n", + " sig_rt = sigs_polar[:, 0, 1]\n", + "\n", " return sig_rr, sig_tt, sig_rt\n", "\n", + "\n", "def get_sort_indices_and_distances_by_distance_from_origin_2D(mesh):\n", " xs = mesh.points[:, 0]\n", " ys = mesh.points[:, 1]\n", " dist_from_origin = np.hypot(xs, ys)\n", " indices_sorted = np.argsort(dist_from_origin)\n", " dist_sorted = dist_from_origin[indices_sorted]\n", - " \n", + "\n", " return indices_sorted, dist_sorted\n", "\n", + "\n", "def compute_abs_and_rel_stress_error_rr(sigmas_rr_num, rs, theta_degree):\n", " num_points = sigmas_rr_num.shape[0]\n", " f_abs_rr = np.zeros(num_points)\n", " f_rel_rr = np.zeros(num_points)\n", - " \n", + "\n", " for pt_idx in range(num_points):\n", " r = rs[pt_idx]\n", - " \n", + "\n", " sigma_rr_ana = kirsch_sig_rr(10, r, theta_degree, 2)\n", - " \n", + "\n", " sigma_rr_num = sigmas_rr_num[pt_idx] * 1000\n", - " \n", + "\n", " f_abs_rr[pt_idx] = sigma_rr_num - sigma_rr_ana\n", - " \n", + "\n", " if sigma_rr_ana == 0:\n", " f_rel_rr[pt_idx] = f_abs_rr[pt_idx] / 1e-2\n", " else:\n", " f_rel_rr[pt_idx] = f_abs_rr[pt_idx] / sigma_rr_ana\n", - " \n", + "\n", " return f_abs_rr, f_rel_rr\n", - " \n", + "\n", + "\n", "def compute_abs_and_rel_stress_error_rt(sigmas_rt_num, rs, theta_degree):\n", " num_points = sigmas_rt_num.shape[0]\n", " f_abs_rt = np.zeros(num_points)\n", " f_rel_rt = np.zeros(num_points)\n", - " \n", + "\n", " for pt_idx in range(num_points):\n", " r = rs[pt_idx]\n", - " \n", + "\n", " sigma_rt_ana = kirsch_sig_rt(10, r, theta_degree, 2)\n", " sigma_rt_num = sigmas_rt_num[pt_idx] * 1000\n", - " \n", + "\n", " f_abs_rt[pt_idx] = sigma_rt_num - sigma_rt_ana\n", - " \n", + "\n", " if sigma_rt_ana == 0:\n", " f_rel_rt[pt_idx] = f_abs_rt[pt_idx] / 1e-2\n", " else:\n", " f_rel_rt[pt_idx] = f_abs_rt[pt_idx] / sigma_rt_ana\n", - " \n", + "\n", " return f_abs_rt, f_rel_rt\n", - " \n", + "\n", + "\n", "def compute_abs_and_rel_stress_error_tt(sigmas_tt_num, rs, theta_degree):\n", " num_points = sigmas_tt_num.shape[0]\n", " f_abs_tt = np.zeros(num_points)\n", " f_rel_tt = np.zeros(num_points)\n", - " \n", + "\n", " for pt_idx in range(num_points):\n", " r = rs[pt_idx]\n", - " \n", + "\n", " sigma_tt_ana = kirsch_sig_tt(10, r, theta_degree, 2)\n", " sigma_tt_num = sigmas_tt_num[pt_idx] * 1000\n", - " \n", + "\n", " f_abs_tt[pt_idx] = sigma_tt_num - sigma_tt_ana\n", - " \n", + "\n", " if sigma_tt_ana == 0:\n", " f_rel_tt[pt_idx] = f_abs_tt[pt_idx] / 1e-2\n", " else:\n", " f_rel_tt[pt_idx] = f_abs_tt[pt_idx] / sigma_tt_ana\n", - " \n", + "\n", " return f_abs_tt, f_rel_tt\n", "\n", + "\n", "def compute_cell_size(idx, mesh):\n", - " pt1 = (19.999,0,0)\n", - " pt2 = (19.999,20,0)\n", + " pt1 = (19.999, 0, 0)\n", + " pt2 = (19.999, 20, 0)\n", " line_mesh = slice_along_line(mesh, pt1, pt2)\n", - " number = line_mesh.points.shape[0]-1 # number of cells along the right edge of the plate\n", - " size = STUDY_mesh_size/number # height of plate divided by number of cells \n", + " number = (\n", + " line_mesh.points.shape[0] - 1\n", + " ) # number of cells along the right edge of the plate\n", + " size = STUDY_mesh_size / number # height of plate divided by number of cells\n", " return size\n", "\n", + "\n", "def resample_mesh_to_240_resolution(idx):\n", " mesh_fine = read_last_timestep_mesh(240)\n", " mesh_coarse = read_last_timestep_mesh(idx)\n", @@ -280,10 +291,10 @@ "id": "c8786de9-0b93-4e5b-b8ef-4c46d6d9d582", "metadata": {}, "source": [ - "To generate the meshes and to enable the individual parameters to be adjusted from within this notebook the script \"mesh_quarter_of_rectangle_with_hole.py\" is used. \n", + "To generate the meshes and to enable the individual parameters to be adjusted from within this notebook the script \"mesh_quarter_of_rectangle_with_hole.py\" is used.\n", "For the considered problem the parameter $a$ and $b$ which define the size of the rectangular quarter of the plate are set to a value of $20\\, \\text{cm}$. The radius of the central hole is controlled by the value of $r$, which is set at $r = 2\\, \\text{cm}$.\n", "\n", - "The arguments Nx, Ny, NR and Nr of the script describe the refinement of the mesh and are determined according to each refinement index. \n", + "The arguments Nx, Ny, NR and Nr of the script describe the refinement of the mesh and are determined according to each refinement index.\n", "\n", "Parameter $R$ describes the area for a further refinement in the vicinity of the hole. This additional refinement is needed for better capturing the stress and strain gradients near the hole. Because the fine resolution of the mesh is only needed in the area around the hole, the size of $R$ is half the size of the plate." ] @@ -519,7 +530,18 @@ " please check that, and maybe document in mesh_quarter_of_rectangle_with_hole.py\n", " \"\"\"\n", " output_file = f\"{out_dir}/disc_with_hole_idx_is_{idx}.msh\"\n", - " mesh_quarter_of_rectangle_with_hole.run(output_file, a=STUDY_mesh_size, b=STUDY_mesh_size, r=2, R=0.5*STUDY_mesh_size, Nx=idx, Ny=idx, NR=idx, Nr=idx, P=1)" + " mesh_quarter_of_rectangle_with_hole.run(\n", + " output_file,\n", + " a=STUDY_mesh_size,\n", + " b=STUDY_mesh_size,\n", + " r=2,\n", + " R=0.5 * STUDY_mesh_size,\n", + " Nx=idx,\n", + " Ny=idx,\n", + " NR=idx,\n", + " Nr=idx,\n", + " P=1,\n", + " )" ] }, { @@ -564,7 +586,7 @@ "id": "865e9ab2-4ad6-4ace-8104-d24370024ce3", "metadata": {}, "source": [ - "To get a better sense of cell sizes and the additional refinement around the hole, the meshes of refinement indices 8 and 80 are shown below. " + "To get a better sense of cell sizes and the additional refinement around the hole, the meshes of refinement indices 8 and 80 are shown below." ] }, { @@ -580,6 +602,7 @@ "outputs": [], "source": [ "import pyvista as pv\n", + "\n", "pv.set_plot_theme(\"document\")\n", "pv.set_jupyter_backend(\"static\")" ] @@ -611,18 +634,18 @@ "domain_80 = pv.read(f\"{out_dir}/disc_with_hole_idx_is_80_\" + \"domain.vtu\")\n", "\n", "p = pv.Plotter(shape=(1, 2), border=False)\n", - "p.subplot(0,0)\n", + "p.subplot(0, 0)\n", "p.add_mesh(domain_8, show_edges=True, show_scalar_bar=False, color=None, scalars=None)\n", "p.view_xy()\n", "p.show_bounds(ticks=\"outside\", xlabel=\"x / m\", ylabel=\"y / m\")\n", "p.camera.zoom(1.3)\n", "\n", - "p.subplot(0,1)\n", + "p.subplot(0, 1)\n", "p.add_mesh(domain_80, show_edges=True, show_scalar_bar=False, color=None, scalars=None)\n", "p.view_xy()\n", "p.show_bounds(ticks=\"outside\", xlabel=\"x / m\", ylabel=\"y / m\")\n", "p.camera.zoom(1.3)\n", - "p.window_size = [1000,500]\n", + "p.window_size = [1000, 500]\n", "\n", "p.show()" ] @@ -659,6 +682,7 @@ "jupyter": { "source_hidden": true }, + "lines_to_next_cell": 2, "tags": [] }, "outputs": [ @@ -682,15 +706,15 @@ } ], "source": [ - "# ATTENTION: We exclude the last study index, because its simulation takes \n", + "# ATTENTION: We exclude the last study index, because its simulation takes\n", "# too long to be included in the OGS CI pipelines.\n", "for idx in STUDY_indices[:-1]:\n", " prj_file = f\"disc_with_hole_idx_is_{idx}.prj\"\n", - " \n", + "\n", " shutil.copy2(prj_file, out_dir)\n", - " \n", + "\n", " prj_path = os.path.join(out_dir, prj_file)\n", - " \n", + "\n", " model = ogs.OGS(INPUT_FILE=prj_path, PROJECT_FILE=prj_path)\n", " model.run_model(logfile=f\"{out_dir}/out.txt\", args=f\"-o {out_dir}\")" ] @@ -716,11 +740,13 @@ { "cell_type": "markdown", "id": "b9be6c95-f016-4531-ac49-6b0b85d0d6ae", - "metadata": {}, + "metadata": { + "lines_to_next_cell": 2 + }, "source": [ "In this section we compare the numerical solution for the different refinement levels to Kirsch's analytical solution.\n", "We show plots of the radial, tangential and shear stress distribution along the x-axis.\n", - "Additionally, the absolute error of each refinement level to the analytical solution is plotted for better interpretation of variations. \n", + "Additionally, the absolute error of each refinement level to the analytical solution is plotted for better interpretation of variations.\n", "\n", "The plots show that–for the considered axis–the stress distribution around the hole converges to Kirsch´s solution with decreasing size of the mesh cells.\n", "This decrease of the extrapolation errors can also be seen in the plots for the absolute error.\n", @@ -745,14 +771,42 @@ }, "outputs": [], "source": [ - "def kirsch_sig_rr(sig,r,theta,a):\n", - " return 0.5*sig*((1-a**2/r**2)+(1+3*np.power(a,4)/np.power(r,4)-4*a**2/r**2)*np.cos(2*np.pi*theta/180)) * np.heaviside(r+1e-7-a,1) \n", - " \n", - "def kirsch_sig_tt(sig,r,theta,a):\n", - " return 0.5*sig*((1+a**2/r**2)-(1+3*np.power(a,4)/np.power(r,4))*np.cos(2*np.pi*theta/180)) * np.heaviside(r+1e-7-a,1)\n", - "\n", - "def kirsch_sig_rt(sig,r,theta,a):\n", - " return -0.5*sig*((1-3*np.power(a,4)/np.power(r,4)+2*a**2/r**2)*np.sin(2*np.pi*theta/180)) * np.heaviside(r+1e-7-a,1)" + "def kirsch_sig_rr(sig, r, theta, a):\n", + " return (\n", + " 0.5\n", + " * sig\n", + " * (\n", + " (1 - a**2 / r**2)\n", + " + (1 + 3 * np.power(a, 4) / np.power(r, 4) - 4 * a**2 / r**2)\n", + " * np.cos(2 * np.pi * theta / 180)\n", + " )\n", + " * np.heaviside(r + 1e-7 - a, 1)\n", + " )\n", + "\n", + "\n", + "def kirsch_sig_tt(sig, r, theta, a):\n", + " return (\n", + " 0.5\n", + " * sig\n", + " * (\n", + " (1 + a**2 / r**2)\n", + " - (1 + 3 * np.power(a, 4) / np.power(r, 4))\n", + " * np.cos(2 * np.pi * theta / 180)\n", + " )\n", + " * np.heaviside(r + 1e-7 - a, 1)\n", + " )\n", + "\n", + "\n", + "def kirsch_sig_rt(sig, r, theta, a):\n", + " return (\n", + " -0.5\n", + " * sig\n", + " * (\n", + " (1 - 3 * np.power(a, 4) / np.power(r, 4) + 2 * a**2 / r**2)\n", + " * np.sin(2 * np.pi * theta / 180)\n", + " )\n", + " * np.heaviside(r + 1e-7 - a, 1)\n", + " )" ] }, { @@ -768,45 +822,48 @@ "outputs": [], "source": [ "def vec4_to_mat3x3cart(vec4):\n", - " m = np.zeros((3,3))\n", - " m[0,0] = vec4[0]\n", - " m[1,1] = vec4[1]\n", - " m[2,2] = vec4[2]\n", - " m[0,1] = vec4[3]\n", - " m[1,0] = vec4[3]\n", - " \n", + " m = np.zeros((3, 3))\n", + " m[0, 0] = vec4[0]\n", + " m[1, 1] = vec4[1]\n", + " m[2, 2] = vec4[2]\n", + " m[0, 1] = vec4[3]\n", + " m[1, 0] = vec4[3]\n", + "\n", " return np.matrix(m)\n", "\n", + "\n", "def vec4_to_mat3x3cart_multi(vec4):\n", " assert vec4.shape[1] == 4\n", - " \n", + "\n", " n_pts = vec4.shape[0]\n", - " \n", + "\n", " m = np.zeros((n_pts, 3, 3))\n", - " m[:,0,0] = vec4[:,0]\n", - " m[:,1,1] = vec4[:,1]\n", - " m[:,2,2] = vec4[:,2]\n", - " m[:,0,1] = vec4[:,3]\n", - " m[:,1,0] = vec4[:,3]\n", - " \n", + " m[:, 0, 0] = vec4[:, 0]\n", + " m[:, 1, 1] = vec4[:, 1]\n", + " m[:, 2, 2] = vec4[:, 2]\n", + " m[:, 0, 1] = vec4[:, 3]\n", + " m[:, 1, 0] = vec4[:, 3]\n", + "\n", " return m\n", "\n", + "\n", "def vec4_to_mat3x3polar_single(vec4, xs, ys):\n", " m_cart = vec4_to_mat3x3cart(vec4)\n", - " \n", + "\n", " theta = np.arctan2(ys, xs)\n", "\n", " rot = np.matrix(np.eye(3))\n", - " rot[0,0] = np.cos(theta)\n", - " rot[0,1] = -np.sin(theta)\n", - " rot[1,0] = np.sin(theta)\n", - " rot[1,1] = np.cos(theta)\n", - " \n", + " rot[0, 0] = np.cos(theta)\n", + " rot[0, 1] = -np.sin(theta)\n", + " rot[1, 0] = np.sin(theta)\n", + " rot[1, 1] = np.cos(theta)\n", + "\n", " return rot.T * m_cart * rot\n", "\n", + "\n", "def vec4_to_mat3x3polar_multi(vecs4, xs, ys):\n", " \"\"\"Convert 4-vectors (Kelvin vector in 2D) to 3x3 matrices in polar coordinates at multiple points at once.\n", - " \n", + "\n", " Parameters\n", " ----------\n", " vecs4:\n", @@ -815,37 +872,38 @@ " NumPy array of x coordinates, length: N\n", " ys:\n", " NumPy array of y coordinates, length: N\n", - " \n", + "\n", " Returns\n", " -------\n", " A Numpy array of the symmetric matrices corresponding to the 4-vectors, dimensions: (N x 3 x 3)\n", " \"\"\"\n", - " \n", - " n_pts = vecs4.shape[0] \n", + "\n", + " n_pts = vecs4.shape[0]\n", " assert n_pts == xs.shape[0]\n", " assert n_pts == ys.shape[0]\n", " assert vecs4.shape[1] == 4\n", - " \n", - " m_carts = vec4_to_mat3x3cart_multi(vecs4) # vecs4 converted to symmetric matrices\n", - " \n", + "\n", + " m_carts = vec4_to_mat3x3cart_multi(vecs4) # vecs4 converted to symmetric matrices\n", + "\n", " thetas = np.arctan2(ys, xs)\n", "\n", - " rots = np.zeros((n_pts, 3, 3)) # rotation matrices at each point\n", - " rots[:,0,0] = np.cos(thetas)\n", - " rots[:,0,1] = -np.sin(thetas)\n", - " rots[:,1,0] = np.sin(thetas)\n", - " rots[:,1,1] = np.cos(thetas)\n", - " rots[:,2,2] = 1\n", - " \n", + " rots = np.zeros((n_pts, 3, 3)) # rotation matrices at each point\n", + " rots[:, 0, 0] = np.cos(thetas)\n", + " rots[:, 0, 1] = -np.sin(thetas)\n", + " rots[:, 1, 0] = np.sin(thetas)\n", + " rots[:, 1, 1] = np.cos(thetas)\n", + " rots[:, 2, 2] = 1\n", + "\n", " # rot.T * m_cart * rot for each point\n", " m_polars = np.einsum(\"...ji,...jk,...kl\", rots, m_carts, rots)\n", - " \n", + "\n", " assert m_polars.shape[0] == n_pts\n", " assert m_polars.shape[1] == 3\n", " assert m_polars.shape[2] == 3\n", - " \n", + "\n", " return m_polars\n", "\n", + "\n", "def vec4_to_mat3x3polar(vec4, xs, ys):\n", " if len(vec4.shape) == 1:\n", " # only a single 4-vector will be converted\n", @@ -871,11 +929,13 @@ "# We'll be able to reuse/read them many times\n", "STUDY_num_result_meshes_by_index = {}\n", "\n", + "\n", "def read_simulation_result_meshes():\n", " for idx in STUDY_indices:\n", " mesh = read_last_timestep_mesh(idx)\n", " STUDY_num_result_meshes_by_index[idx] = mesh\n", - " \n", + "\n", + "\n", "read_simulation_result_meshes()" ] }, @@ -893,15 +953,17 @@ "source": [ "STUDY_num_result_xaxis_meshes_by_index = {}\n", "\n", + "\n", "def compute_xaxis_meshes():\n", " for idx in STUDY_indices:\n", " mesh = STUDY_num_result_meshes_by_index[idx]\n", - " pt1 = (0,1e-6,0)\n", - " pt2 = (10,1e-6,0)\n", + " pt1 = (0, 1e-6, 0)\n", + " pt2 = (10, 1e-6, 0)\n", " line_mesh = slice_along_line(mesh, pt1, pt2)\n", "\n", " STUDY_num_result_xaxis_meshes_by_index[idx] = line_mesh\n", - " \n", + "\n", + "\n", "compute_xaxis_meshes()" ] }, @@ -919,15 +981,17 @@ "source": [ "STUDY_num_result_yaxis_meshes_by_index = {}\n", "\n", + "\n", "def compute_yaxis_meshes():\n", " for idx in STUDY_indices:\n", " mesh = STUDY_num_result_meshes_by_index[idx]\n", - " pt1 = (1e-6,0,0)\n", - " pt2 = (1e-6,10,0)\n", + " pt1 = (1e-6, 0, 0)\n", + " pt2 = (1e-6, 10, 0)\n", " line_mesh = slice_along_line(mesh, pt1, pt2)\n", "\n", " STUDY_num_result_yaxis_meshes_by_index[idx] = line_mesh\n", - " \n", + "\n", + "\n", "compute_yaxis_meshes()" ] }, @@ -939,28 +1003,33 @@ "jupyter": { "source_hidden": true }, + "lines_to_next_cell": 2, "tags": [] }, "outputs": [], "source": [ "STUDY_num_result_diagonal_meshes_by_index = {}\n", "\n", + "\n", "def compute_diagonal_meshes():\n", " for idx in STUDY_indices:\n", " mesh = STUDY_num_result_meshes_by_index[idx]\n", - " pt1 = (1e-6,1e-6,0)\n", - " pt2 = (28.28,28.28,0)\n", + " pt1 = (1e-6, 1e-6, 0)\n", + " pt2 = (28.28, 28.28, 0)\n", " line_mesh = slice_along_line(mesh, pt1, pt2)\n", "\n", " STUDY_num_result_diagonal_meshes_by_index[idx] = line_mesh\n", - " \n", + "\n", + "\n", "compute_diagonal_meshes()" ] }, { "cell_type": "markdown", "id": "4f89b0ef-cdba-4523-8f9a-fcf0d4c55dc9", - "metadata": {}, + "metadata": { + "lines_to_next_cell": 2 + }, "source": [ "### Stress distribution along the x-axis" ] @@ -992,30 +1061,33 @@ "source": [ "def plot_stress_distribution_along_xaxis():\n", " ### Step 1: Compute data ##########################################\n", - " \n", + "\n", " # These variables will hold the error data for all STUDY_indices\n", " f_abs_rr = {}\n", " f_abs_tt = {}\n", " f_abs_rt = {}\n", - " \n", - " #Plot setup\n", - " fig, ax = plt.subplots(nrows = 2, ncols=3, figsize=(22,10))\n", - " for i in range (2):\n", - " for j in range (3):\n", - " ax[i][j].axvline(2, color=\"0.6\", linestyle = \"--\", label = \"Hole Radius\")\n", + "\n", + " # Plot setup\n", + " fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(22, 10))\n", + " for i in range(2):\n", + " for j in range(3):\n", + " ax[i][j].axvline(2, color=\"0.6\", linestyle=\"--\", label=\"Hole Radius\")\n", " ax[i][j].grid(True)\n", - " ax[i][j].set(xlim=(0,STUDY_mesh_size))\n", - " ax[i][j].set_xlabel('$r$ / cm')\n", - " ax[i][1].set_ylabel('$\\Delta\\\\sigma$ / kPa')\n", - " ax[i][2].set_ylabel('$\\Delta\\\\sigma$ / $\\sigma_{\\mathrm{analytical}}$')\n", - " \n", + " ax[i][j].set(xlim=(0, STUDY_mesh_size))\n", + " ax[i][j].set_xlabel(\"$r$ / cm\")\n", + " ax[i][1].set_ylabel(\"$\\Delta\\\\sigma$ / kPa\")\n", + " ax[i][2].set_ylabel(\"$\\Delta\\\\sigma$ / $\\sigma_{\\mathrm{analytical}}$\")\n", + "\n", " for iteration, idx in enumerate(STUDY_indices):\n", " # we use the line mesh we extracted before\n", " line_mesh = STUDY_num_result_xaxis_meshes_by_index[idx]\n", "\n", " sig_rr, sig_tt, sig_rt = get_sigma_polar_components(line_mesh)\n", "\n", - " indices_sorted, dist_sorted = get_sort_indices_and_distances_by_distance_from_origin_2D(line_mesh)\n", + " (\n", + " indices_sorted,\n", + " dist_sorted,\n", + " ) = get_sort_indices_and_distances_by_distance_from_origin_2D(line_mesh)\n", "\n", " # sort sigma by distance from origin\n", " sig_rr_sorted = sig_rr[indices_sorted]\n", @@ -1024,68 +1096,149 @@ "\n", " # compute errors\n", " f_abs_rr, f_rel_rr = compute_abs_and_rel_stress_error_rr(\n", - " sig_rr_sorted, dist_sorted, -90)\n", + " sig_rr_sorted, dist_sorted, -90\n", + " )\n", " f_abs_tt, f_rel_tt = compute_abs_and_rel_stress_error_tt(\n", - " sig_tt_sorted, dist_sorted, -90)\n", + " sig_tt_sorted, dist_sorted, -90\n", + " )\n", " f_abs_rt, f_rel_rt = compute_abs_and_rel_stress_error_rt(\n", - " sig_rt_sorted, dist_sorted, -90)\n", - " \n", + " sig_rt_sorted, dist_sorted, -90\n", + " )\n", + "\n", " ### Step 2: Plot data ##############################################\n", - " \n", - " ax[0][0].set_ylabel('$\\\\sigma_{rr}$ / kPa')\n", - " ax[0][1].set_ylabel('$\\\\sigma_{\\\\theta\\\\theta}$ / kPa')\n", - " ax[0][2].set_ylabel('$\\\\sigma_{r\\\\theta}$ / kPa')\n", - " \n", + "\n", + " ax[0][0].set_ylabel(\"$\\\\sigma_{rr}$ / kPa\")\n", + " ax[0][1].set_ylabel(\"$\\\\sigma_{\\\\theta\\\\theta}$ / kPa\")\n", + " ax[0][2].set_ylabel(\"$\\\\sigma_{r\\\\theta}$ / kPa\")\n", "\n", " # analytical results\n", " if iteration == 0:\n", - " r = np.linspace(2,STUDY_mesh_size,1000)\n", - " ax[0][0].plot(r, kirsch_sig_rr(10,r,-90,2), color = \"deepskyblue\", linestyle = \":\", label = \"analytical\")\n", - " ax[0][1].plot(r, kirsch_sig_tt(10,r,-90,2), color = \"yellowgreen\", linestyle = \":\", label = \"analytical\")\n", - " ax[0][2].plot(r, kirsch_sig_rt(10,r,-90,2), color = \"orangered\", linestyle = \":\", label = \"analytical\")\n", + " r = np.linspace(2, STUDY_mesh_size, 1000)\n", + " ax[0][0].plot(\n", + " r,\n", + " kirsch_sig_rr(10, r, -90, 2),\n", + " color=\"deepskyblue\",\n", + " linestyle=\":\",\n", + " label=\"analytical\",\n", + " )\n", + " ax[0][1].plot(\n", + " r,\n", + " kirsch_sig_tt(10, r, -90, 2),\n", + " color=\"yellowgreen\",\n", + " linestyle=\":\",\n", + " label=\"analytical\",\n", + " )\n", + " ax[0][2].plot(\n", + " r,\n", + " kirsch_sig_rt(10, r, -90, 2),\n", + " color=\"orangered\",\n", + " linestyle=\":\",\n", + " label=\"analytical\",\n", + " )\n", "\n", " # numerical results\n", " cell_size = compute_cell_size(idx, STUDY_num_result_meshes_by_index[idx])\n", "\n", " if idx == 8:\n", - " ax[0][0].plot(dist_sorted, sig_rr_sorted*1000, color = \"lightskyblue\", label = f\"h = {cell_size:.3f} cm\")\n", - " ax[0][1].plot(dist_sorted, sig_tt_sorted*1000, color = \"limegreen\", label = f\"h = {cell_size:.3f} cm\")\n", - " ax[0][2].plot(dist_sorted, sig_rt_sorted*1000, color = \"lightcoral\", label = f\"h = {cell_size:.3f} cm\")\n", - " ax[1][0].plot(dist_sorted, f_abs_rr, color = \"lightskyblue\")\n", - " ax[1][1].plot(dist_sorted, f_abs_tt, color = \"limegreen\")\n", - " ax[1][2].plot(dist_sorted, f_abs_rt, color = \"lightcoral\")\n", - " \n", + " ax[0][0].plot(\n", + " dist_sorted,\n", + " sig_rr_sorted * 1000,\n", + " color=\"lightskyblue\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", + " ax[0][1].plot(\n", + " dist_sorted,\n", + " sig_tt_sorted * 1000,\n", + " color=\"limegreen\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", + " ax[0][2].plot(\n", + " dist_sorted,\n", + " sig_rt_sorted * 1000,\n", + " color=\"lightcoral\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", + " ax[1][0].plot(dist_sorted, f_abs_rr, color=\"lightskyblue\")\n", + " ax[1][1].plot(dist_sorted, f_abs_tt, color=\"limegreen\")\n", + " ax[1][2].plot(dist_sorted, f_abs_rt, color=\"lightcoral\")\n", + "\n", " if idx == 16:\n", - " ax[0][0].plot(dist_sorted, sig_rr_sorted*1000, color = \"cornflowerblue\", label = f\"h = {cell_size:.3f} cm\")\n", - " ax[0][1].plot(dist_sorted, sig_tt_sorted*1000, color = \"forestgreen\", label = f\"h = {cell_size:.3f} cm\")\n", - " ax[0][2].plot(dist_sorted, sig_rt_sorted*1000, color = \"firebrick\", label = f\"h = {cell_size:.3f} cm\")\n", - " ax[1][0].plot(dist_sorted, f_abs_rr, color = \"cornflowerblue\")\n", - " ax[1][1].plot(dist_sorted, f_abs_tt, color = \"forestgreen\")\n", - " ax[1][2].plot(dist_sorted, f_abs_rt, color = \"firebrick\")\n", - " \n", + " ax[0][0].plot(\n", + " dist_sorted,\n", + " sig_rr_sorted * 1000,\n", + " color=\"cornflowerblue\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", + " ax[0][1].plot(\n", + " dist_sorted,\n", + " sig_tt_sorted * 1000,\n", + " color=\"forestgreen\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", + " ax[0][2].plot(\n", + " dist_sorted,\n", + " sig_rt_sorted * 1000,\n", + " color=\"firebrick\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", + " ax[1][0].plot(dist_sorted, f_abs_rr, color=\"cornflowerblue\")\n", + " ax[1][1].plot(dist_sorted, f_abs_tt, color=\"forestgreen\")\n", + " ax[1][2].plot(dist_sorted, f_abs_rt, color=\"firebrick\")\n", + "\n", " if idx == 24:\n", - " ax[0][0].plot(dist_sorted, sig_rr_sorted*1000, color = \"royalblue\", label = f\"h = {cell_size:.3f} cm\")\n", - " ax[0][1].plot(dist_sorted, sig_tt_sorted*1000, color = \"darkgreen\", label = f\"h = {cell_size:.3f} cm\")\n", - " ax[0][2].plot(dist_sorted, sig_rt_sorted*1000, color = \"darkred\", label = f\"h = {cell_size:.3f} cm\")\n", - " ax[1][0].plot(dist_sorted, f_abs_rr, color = \"royalblue\")\n", - " ax[1][1].plot(dist_sorted, f_abs_tt, color = \"darkgreen\")\n", - " ax[1][2].plot(dist_sorted, f_abs_rt, color = \"darkred\")\n", - " \n", + " ax[0][0].plot(\n", + " dist_sorted,\n", + " sig_rr_sorted * 1000,\n", + " color=\"royalblue\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", + " ax[0][1].plot(\n", + " dist_sorted,\n", + " sig_tt_sorted * 1000,\n", + " color=\"darkgreen\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", + " ax[0][2].plot(\n", + " dist_sorted,\n", + " sig_rt_sorted * 1000,\n", + " color=\"darkred\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", + " ax[1][0].plot(dist_sorted, f_abs_rr, color=\"royalblue\")\n", + " ax[1][1].plot(dist_sorted, f_abs_tt, color=\"darkgreen\")\n", + " ax[1][2].plot(dist_sorted, f_abs_rt, color=\"darkred\")\n", + "\n", " if idx == 240:\n", - " ax[0][0].plot(dist_sorted, sig_rr_sorted*1000, color = \"black\", label = f\"h = {cell_size:.3f} cm\")\n", - " ax[0][1].plot(dist_sorted, sig_tt_sorted*1000, color = \"black\", label = f\"h = {cell_size:.3f} cm\")\n", - " ax[0][2].plot(dist_sorted, sig_rt_sorted*1000, color = \"black\", label = f\"h = {cell_size:.3f} cm\") \n", + " ax[0][0].plot(\n", + " dist_sorted,\n", + " sig_rr_sorted * 1000,\n", + " color=\"black\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", + " ax[0][1].plot(\n", + " dist_sorted,\n", + " sig_tt_sorted * 1000,\n", + " color=\"black\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", + " ax[0][2].plot(\n", + " dist_sorted,\n", + " sig_rt_sorted * 1000,\n", + " color=\"black\",\n", + " label=f\"h = {cell_size:.3f} cm\",\n", + " )\n", "\n", - " # final plot settings \n", - " for i in range (3):\n", + " # final plot settings\n", + " for i in range(3):\n", " ax[0][i].legend()\n", "\n", - " ax[0][0].set_title('Radial stress distribution')\n", - " ax[0][1].set_title('Tangential stress distribution')\n", - " ax[0][2].set_title('Shear stress distribution')\n", + " ax[0][0].set_title(\"Radial stress distribution\")\n", + " ax[0][1].set_title(\"Tangential stress distribution\")\n", + " ax[0][2].set_title(\"Shear stress distribution\")\n", "\n", " fig.tight_layout()\n", "\n", + "\n", "plot_stress_distribution_along_xaxis()" ] }, @@ -1143,11 +1296,11 @@ "id": "d8a28954-faff-4049-8d74-2882112dcdb7", "metadata": {}, "source": [ - "The $\\ell_{2}$ norm or Euclidean norm is the square root of the sum of the squared absolute errors at each point of the mesh. \n", + "The $\\ell_{2}$ norm or Euclidean norm is the square root of the sum of the squared absolute errors at each point of the mesh.\n", "\n", "The root mean square ($RMS$) is calculated similarly. Here, however, the influence of the number of points is taken into account by dividing by the square root of the number of points.\n", "\n", - "The $L_{2}$ norm as integral norm represents a generalization of the $\\ell_{2}$ norm for continuous functions. While the Euclidean norm only considers the values on individual mesh nodes, the integral norm considers the solution on the entire mesh. Therefore an advantage of the $L_{2}$ norm is, that big elements are considered with a higher impact than small ones, which produces more homogeneous results. " + "The $L_{2}$ norm as integral norm represents a generalization of the $\\ell_{2}$ norm for continuous functions. While the Euclidean norm only considers the values on individual mesh nodes, the integral norm considers the solution on the entire mesh. Therefore an advantage of the $L_{2}$ norm is, that big elements are considered with a higher impact than small ones, which produces more homogeneous results." ] }, { @@ -1155,8 +1308,8 @@ "id": "d5a07337-072f-4539-bf9d-7eaac25cbbe4", "metadata": {}, "source": [ - "The following plots represent the development of the Euclidean and Integral norm and $RMS$ for the refinement of the mesh. How fast the considered element converge is expressed by the slope of the lines in the plot. \n", - "First the detailed discussed error norms for the stresses are visualised and in addition to them the norms for the associated displacements to draw conclusions about the quality of convergence. \n", + "The following plots represent the development of the Euclidean and Integral norm and $RMS$ for the refinement of the mesh. How fast the considered element converge is expressed by the slope of the lines in the plot.\n", + "First the detailed discussed error norms for the stresses are visualised and in addition to them the norms for the associated displacements to draw conclusions about the quality of convergence.\n", "\n", "The main conclusion that can be drawn is that the solution for the displacements converge significantly faster than those for the stresses.\n", "As a practical consequence, it might be possible to get a sufficiently accurate displacement solution already on a relatively coarse mesh,\n", @@ -1194,7 +1347,9 @@ " integrator = vtkIntegrateAttributes()\n", " integrator.SetInputData(mesh)\n", " integrator.Update()\n", - " return pv.wrap(integrator.GetOutputDataObject(0)) # that is an entire mesh with one point and one cell" + " return pv.wrap(\n", + " integrator.GetOutputDataObject(0)\n", + " ) # that is an entire mesh with one point and one cell" ] }, { @@ -1212,15 +1367,15 @@ "def compute_ell_2_norm_sigma(idx, sigmas_test, sigmas_reference):\n", " sig_rr, sig_tt, sig_rt = sigmas_test\n", " sig_rr_240, sig_tt_240, sig_rt_240 = sigmas_reference\n", - " \n", - " list_rr = (sig_rr_240 - sig_rr)**2\n", - " list_tt = (sig_tt_240 - sig_tt)**2\n", - " list_rt = (sig_rt_240 - sig_rt)**2\n", - " \n", + "\n", + " list_rr = (sig_rr_240 - sig_rr) ** 2\n", + " list_tt = (sig_tt_240 - sig_tt) ** 2\n", + " list_rt = (sig_rt_240 - sig_rt) ** 2\n", + "\n", " l2_rr = np.sqrt(sum(list_rr))\n", " l2_tt = np.sqrt(sum(list_tt))\n", " l2_rt = np.sqrt(sum(list_rt))\n", - " \n", + "\n", " return l2_rr, l2_tt, l2_rt" ] }, @@ -1238,19 +1393,19 @@ "source": [ "def compute_ell_2_norm_displacement(idx, mesh_resampled_to_240_resolution, mesh_fine):\n", " dis = mesh_resampled_to_240_resolution.point_data[\"displacement\"]\n", - " dis_x = dis[:,0]\n", - " dis_y = dis[:,1]\n", - " \n", + " dis_x = dis[:, 0]\n", + " dis_y = dis[:, 1]\n", + "\n", " dis_240 = mesh_fine.point_data[\"displacement\"]\n", - " dis_x_240 = dis_240[:,0]\n", - " dis_y_240 = dis_240[:,1]\n", - " \n", - " list_x = (dis_x_240 - dis_x)**2\n", - " list_y = (dis_y_240 - dis_y)**2\n", - " \n", + " dis_x_240 = dis_240[:, 0]\n", + " dis_y_240 = dis_240[:, 1]\n", + "\n", + " list_x = (dis_x_240 - dis_x) ** 2\n", + " list_y = (dis_y_240 - dis_y) ** 2\n", + "\n", " l2_x = np.sqrt(sum(list_x))\n", " l2_y = np.sqrt(sum(list_y))\n", - " \n", + "\n", " return l2_x, l2_y" ] }, @@ -1269,13 +1424,13 @@ "def compute_root_mean_square_sigma(idx, sigmas_test, sigmas_reference):\n", " sig_rr, sig_tt, sig_rt = sigmas_test\n", " sig_rr_240, sig_tt_240, sig_rt_240 = sigmas_reference\n", - " \n", + "\n", " l2_rr = np.linalg.norm(sig_rr_240 - sig_rr)\n", " l2_tt = np.linalg.norm(sig_tt_240 - sig_tt)\n", " l2_rt = np.linalg.norm(sig_rt_240 - sig_rt)\n", - " \n", + "\n", " points = sig_rr.shape[0]\n", - " return l2_rr/np.sqrt(points), l2_tt/np.sqrt(points), l2_rt/np.sqrt(points)" + " return l2_rr / np.sqrt(points), l2_tt / np.sqrt(points), l2_rt / np.sqrt(points)" ] }, { @@ -1290,21 +1445,23 @@ }, "outputs": [], "source": [ - "def compute_root_mean_square_displacement(idx, mesh_resampled_to_240_resolution, mesh_fine):\n", - " points = mesh_resampled_to_240_resolution.point_data[\"sigma\"].shape[0] \n", - " \n", + "def compute_root_mean_square_displacement(\n", + " idx, mesh_resampled_to_240_resolution, mesh_fine\n", + "):\n", + " points = mesh_resampled_to_240_resolution.point_data[\"sigma\"].shape[0]\n", + "\n", " dis = mesh_resampled_to_240_resolution.point_data[\"displacement\"]\n", - " dis_x = dis[:,0]\n", - " dis_y = dis[:,1]\n", - " \n", + " dis_x = dis[:, 0]\n", + " dis_y = dis[:, 1]\n", + "\n", " dis_240 = mesh_fine.point_data[\"displacement\"]\n", - " dis_x_240 = dis_240[:,0]\n", - " dis_y_240 = dis_240[:,1]\n", - " \n", + " dis_x_240 = dis_240[:, 0]\n", + " dis_y_240 = dis_240[:, 1]\n", + "\n", " l2_x = np.linalg.norm(dis_x_240 - dis_x)\n", " l2_y = np.linalg.norm(dis_y_240 - dis_y)\n", - " \n", - " return l2_x/np.sqrt(points), l2_y/np.sqrt(points)" + "\n", + " return l2_x / np.sqrt(points), l2_y / np.sqrt(points)" ] }, { @@ -1319,27 +1476,31 @@ }, "outputs": [], "source": [ - "def compute_Ell_2_norm_sigma(idx, mesh_resampled_to_240_resolution, sigmas_test, sigmas_reference):\n", + "def compute_Ell_2_norm_sigma(\n", + " idx, mesh_resampled_to_240_resolution, sigmas_test, sigmas_reference\n", + "):\n", " sig_rr, sig_tt, sig_rt = sigmas_test\n", " sig_rr_240, sig_tt_240, sig_rt_240 = sigmas_reference\n", - " \n", - " list_rr = (sig_rr_240 - sig_rr)**2\n", - " list_tt = (sig_tt_240 - sig_tt)**2\n", - " list_rt = (sig_rt_240 - sig_rt)**2\n", - " \n", + "\n", + " list_rr = (sig_rr_240 - sig_rr) ** 2\n", + " list_tt = (sig_tt_240 - sig_tt) ** 2\n", + " list_rt = (sig_rt_240 - sig_rt) ** 2\n", + "\n", " # We add the squared differences as new point data to the mesh\n", " mesh_resampled_to_240_resolution.point_data[\"diff_rr_squared\"] = list_rr\n", " mesh_resampled_to_240_resolution.point_data[\"diff_tt_squared\"] = list_tt\n", " mesh_resampled_to_240_resolution.point_data[\"diff_rt_squared\"] = list_rt\n", - " \n", + "\n", " # this will integrate all fields at once, so you can add the tt and rt components above and call this only once.\n", - " integration_result_mesh = integrate_mesh_attributes(mesh_resampled_to_240_resolution)\n", - " \n", + " integration_result_mesh = integrate_mesh_attributes(\n", + " mesh_resampled_to_240_resolution\n", + " )\n", + "\n", " # new: integral norm, the index [0] accesses the data of the single point contained in the mesh\n", " L2_rr = np.sqrt(integration_result_mesh.point_data[\"diff_rr_squared\"][0])\n", " L2_tt = np.sqrt(integration_result_mesh.point_data[\"diff_tt_squared\"][0])\n", " L2_rt = np.sqrt(integration_result_mesh.point_data[\"diff_rt_squared\"][0])\n", - " \n", + "\n", " return L2_rr, L2_tt, L2_rt" ] }, @@ -1357,27 +1518,29 @@ "source": [ "def compute_Ell_2_norm_displacement(idx, mesh_resampled_to_240_resolution, mesh_fine):\n", " dis = mesh_resampled_to_240_resolution.point_data[\"displacement\"]\n", - " dis_x = dis[:,0]\n", - " dis_y = dis[:,1]\n", - " \n", + " dis_x = dis[:, 0]\n", + " dis_y = dis[:, 1]\n", + "\n", " dis_240 = mesh_fine.point_data[\"displacement\"]\n", - " dis_x_240 = dis_240[:,0]\n", - " dis_y_240 = dis_240[:,1]\n", + " dis_x_240 = dis_240[:, 0]\n", + " dis_y_240 = dis_240[:, 1]\n", + "\n", + " list_x = (dis_x_240 - dis_x) ** 2\n", + " list_y = (dis_y_240 - dis_y) ** 2\n", "\n", - " list_x = (dis_x_240 - dis_x)**2\n", - " list_y = (dis_y_240 - dis_y)**2\n", - " \n", " # We add the squared differences as new point data to the mesh\n", " mesh_resampled_to_240_resolution.point_data[\"diff_x_squared\"] = list_x\n", " mesh_resampled_to_240_resolution.point_data[\"diff_y_squared\"] = list_y\n", - " \n", + "\n", " # this will integrate all fields at once, so you can add the tt and rt components above and call this only once.\n", - " integration_result_mesh = integrate_mesh_attributes(mesh_resampled_to_240_resolution)\n", - " \n", + " integration_result_mesh = integrate_mesh_attributes(\n", + " mesh_resampled_to_240_resolution\n", + " )\n", + "\n", " # new: integral norm, the index [0] accesses the data of the single point contained in the mesh\n", " L2_x = np.sqrt(integration_result_mesh.point_data[\"diff_x_squared\"][0])\n", " L2_y = np.sqrt(integration_result_mesh.point_data[\"diff_y_squared\"][0])\n", - " \n", + "\n", " return L2_x, L2_y" ] }, @@ -1411,6 +1574,7 @@ "L2_x = {}\n", "L2_y = {}\n", "\n", + "\n", "def compute_error_norms():\n", " mesh_fine = STUDY_num_result_meshes_by_index[240]\n", " sigmas_reference = get_sigma_polar_components(mesh_fine)\n", @@ -1421,15 +1585,28 @@ " mesh_resampled_to_240_resolution = mesh_fine.sample(mesh_coarse)\n", " sigmas_test = get_sigma_polar_components(mesh_resampled_to_240_resolution)\n", "\n", - " l2_rr[idx], l2_tt[idx], l2_rt[idx] = compute_ell_2_norm_sigma(idx, sigmas_test, sigmas_reference)\n", - " rms_rr[idx], rms_tt[idx], rms_rt[idx] = compute_root_mean_square_sigma(idx, sigmas_test, sigmas_reference)\n", - " L2_rr[idx], L2_tt[idx], L2_rt[idx] = compute_Ell_2_norm_sigma(idx, mesh_resampled_to_240_resolution, sigmas_test, sigmas_reference)\n", - " \n", - " l2_x[idx], l2_y[idx] = compute_ell_2_norm_displacement(idx, mesh_resampled_to_240_resolution, mesh_fine)\n", - " rms_x[idx], rms_y[idx] = compute_root_mean_square_displacement(idx, mesh_resampled_to_240_resolution, mesh_fine)\n", - " L2_x[idx], L2_y[idx] = compute_Ell_2_norm_displacement(idx, mesh_resampled_to_240_resolution, mesh_fine)\n", + " l2_rr[idx], l2_tt[idx], l2_rt[idx] = compute_ell_2_norm_sigma(\n", + " idx, sigmas_test, sigmas_reference\n", + " )\n", + " rms_rr[idx], rms_tt[idx], rms_rt[idx] = compute_root_mean_square_sigma(\n", + " idx, sigmas_test, sigmas_reference\n", + " )\n", + " L2_rr[idx], L2_tt[idx], L2_rt[idx] = compute_Ell_2_norm_sigma(\n", + " idx, mesh_resampled_to_240_resolution, sigmas_test, sigmas_reference\n", + " )\n", + "\n", + " l2_x[idx], l2_y[idx] = compute_ell_2_norm_displacement(\n", + " idx, mesh_resampled_to_240_resolution, mesh_fine\n", + " )\n", + " rms_x[idx], rms_y[idx] = compute_root_mean_square_displacement(\n", + " idx, mesh_resampled_to_240_resolution, mesh_fine\n", + " )\n", + " L2_x[idx], L2_y[idx] = compute_Ell_2_norm_displacement(\n", + " idx, mesh_resampled_to_240_resolution, mesh_fine\n", + " )\n", " size[idx] = compute_cell_size(idx, mesh_coarse)\n", - " \n", + "\n", + "\n", "compute_error_norms()" ] }, @@ -1446,14 +1623,15 @@ "source": [ "def plot_slope_sketch(ax, x0, y0, slopes, xmax=None):\n", " \"\"\"Plot sketch for slopes. All slopes cross at (x0, y0)\"\"\"\n", - " if xmax is None: xmax = 2*x0\n", + " if xmax is None:\n", + " xmax = 2 * x0\n", " xs = np.linspace(x0, xmax, 20)\n", - " \n", + "\n", " for slope in slopes:\n", - " y_ = xs[0]**slope\n", - " ys = y0/y_ * xs**slope\n", - " ax.plot(xs, ys, color = \"black\")\n", - " ax.text(xs[-1]*1.05, ys[-1], slope)" + " y_ = xs[0] ** slope\n", + " ys = y0 / y_ * xs**slope\n", + " ax.plot(xs, ys, color=\"black\")\n", + " ax.text(xs[-1] * 1.05, ys[-1], slope)" ] }, { @@ -1481,47 +1659,87 @@ } ], "source": [ - "h = np.linspace(size[80],size[8],1000)\n", - "k = np.linspace(1,2, 20)\n", + "h = np.linspace(size[80], size[8], 1000)\n", + "k = np.linspace(1, 2, 20)\n", "\n", - "fig, ax = plt.subplots(ncols=3, figsize = (20,5))\n", + "fig, ax = plt.subplots(ncols=3, figsize=(20, 5))\n", "\n", - "ax[0].plot(size.values(), l2_rr.values(), color = \"firebrick\", linestyle = \":\", label = \"$\\ell_{2, rr}$\")\n", - "ax[0].plot(size.values(), l2_tt.values(), color = \"firebrick\", linestyle = \"--\", label = \"$\\ell_{2, \\\\theta\\\\theta}$\")\n", - "ax[0].plot(size.values(), l2_rt.values(), color = \"firebrick\", label = \"$\\ell_{2, r\\\\theta}$\")\n", - "ax[0].plot(size.values(), l2_x.values(), color = \"royalblue\", linestyle = \"--\", label = \"$\\ell_{2, x}$\")\n", - "ax[0].plot(size.values(), l2_y.values(), color = \"royalblue\", label = \"$\\ell_{2, y}$\")\n", + "ax[0].plot(\n", + " size.values(),\n", + " l2_rr.values(),\n", + " color=\"firebrick\",\n", + " linestyle=\":\",\n", + " label=\"$\\ell_{2, rr}$\",\n", + ")\n", + "ax[0].plot(\n", + " size.values(),\n", + " l2_tt.values(),\n", + " color=\"firebrick\",\n", + " linestyle=\"--\",\n", + " label=\"$\\ell_{2, \\\\theta\\\\theta}$\",\n", + ")\n", + "ax[0].plot(\n", + " size.values(), l2_rt.values(), color=\"firebrick\", label=\"$\\ell_{2, r\\\\theta}$\"\n", + ")\n", + "ax[0].plot(\n", + " size.values(),\n", + " l2_x.values(),\n", + " color=\"royalblue\",\n", + " linestyle=\"--\",\n", + " label=\"$\\ell_{2, x}$\",\n", + ")\n", + "ax[0].plot(size.values(), l2_y.values(), color=\"royalblue\", label=\"$\\ell_{2, y}$\")\n", "\n", - "plot_slope_sketch(ax[0], 1.5e-1, 5e-2, [1,2,3], xmax=2.5e-1)\n", + "plot_slope_sketch(ax[0], 1.5e-1, 5e-2, [1, 2, 3], xmax=2.5e-1)\n", "\n", - "ax[0].set_title('$\\ell_2$ norms')\n", + "ax[0].set_title(\"$\\ell_2$ norms\")\n", "ax[0].set_ylabel(\"$\\ell_2$ / kPa or cm\")\n", "\n", - "ax[1].plot(size.values(), rms_rr.values(), color = \"firebrick\", linestyle = \":\", label = \"$RMS_{rr}$\")\n", - "ax[1].plot(size.values(), rms_tt.values(), color = \"firebrick\", linestyle = \"--\", label = \"$RMS_{\\\\theta\\\\theta}$\")\n", - "ax[1].plot(size.values(), rms_rt.values(), color = \"firebrick\", label = \"$RMS_{r\\\\theta}$\")\n", - "ax[1].plot(size.values(), rms_x.values(), color = \"royalblue\", linestyle = \"--\", label = \"$RMS_{x}$\")\n", - "ax[1].plot(size.values(), rms_y.values(), color = \"royalblue\", label = \"$RMS_{y}$\")\n", + "ax[1].plot(\n", + " size.values(), rms_rr.values(), color=\"firebrick\", linestyle=\":\", label=\"$RMS_{rr}$\"\n", + ")\n", + "ax[1].plot(\n", + " size.values(),\n", + " rms_tt.values(),\n", + " color=\"firebrick\",\n", + " linestyle=\"--\",\n", + " label=\"$RMS_{\\\\theta\\\\theta}$\",\n", + ")\n", + "ax[1].plot(size.values(), rms_rt.values(), color=\"firebrick\", label=\"$RMS_{r\\\\theta}$\")\n", + "ax[1].plot(\n", + " size.values(), rms_x.values(), color=\"royalblue\", linestyle=\"--\", label=\"$RMS_{x}$\"\n", + ")\n", + "ax[1].plot(size.values(), rms_y.values(), color=\"royalblue\", label=\"$RMS_{y}$\")\n", "\n", - "plot_slope_sketch(ax[1], 1.5e-1, 1e-4, [1,2,3], xmax=2.5e-1)\n", + "plot_slope_sketch(ax[1], 1.5e-1, 1e-4, [1, 2, 3], xmax=2.5e-1)\n", "\n", "ax[1].set_title(\"Root-mean-square\")\n", "ax[1].set_ylabel(\"RMS / kPa or cm\")\n", "\n", - "ax[2].plot(size.values(), L2_rr.values(), color = \"firebrick\", linestyle = \":\", label = \"$L_{2, rr}$\")\n", - "ax[2].plot(size.values(), L2_tt.values(), color = \"firebrick\", linestyle = \"--\", label = \"$L_{2, \\\\theta\\\\theta}$\")\n", - "ax[2].plot(size.values(), L2_rt.values(), color = \"firebrick\", label = \"$L_{2, r\\\\theta}$\")\n", - "ax[2].plot(size.values(), L2_x.values(), color = \"royalblue\", linestyle = \"--\", label = \"$L_{2, x}$\")\n", - "ax[2].plot(size.values(), L2_y.values(), color = \"royalblue\", label = \"$L_{2, y}$\")\n", + "ax[2].plot(\n", + " size.values(), L2_rr.values(), color=\"firebrick\", linestyle=\":\", label=\"$L_{2, rr}$\"\n", + ")\n", + "ax[2].plot(\n", + " size.values(),\n", + " L2_tt.values(),\n", + " color=\"firebrick\",\n", + " linestyle=\"--\",\n", + " label=\"$L_{2, \\\\theta\\\\theta}$\",\n", + ")\n", + "ax[2].plot(size.values(), L2_rt.values(), color=\"firebrick\", label=\"$L_{2, r\\\\theta}$\")\n", + "ax[2].plot(\n", + " size.values(), L2_x.values(), color=\"royalblue\", linestyle=\"--\", label=\"$L_{2, x}$\"\n", + ")\n", + "ax[2].plot(size.values(), L2_y.values(), color=\"royalblue\", label=\"$L_{2, y}$\")\n", "\n", - "plot_slope_sketch(ax[2], 1.5e-1, 1e-3, [1,2,3], xmax=2.5e-1)\n", + "plot_slope_sketch(ax[2], 1.5e-1, 1e-3, [1, 2, 3], xmax=2.5e-1)\n", "\n", "ax[2].set_title(\"$L_2$ norms (integral norms)\")\n", "ax[2].set_ylabel(\"$L_2$ /kPa or cm\")\n", - "for i in range (3):\n", + "for i in range(3):\n", " ax[i].legend()\n", " ax[i].set_xlabel(\"h / cm\")\n", - " ax[i].loglog(base = 10)" + " ax[i].loglog(base=10)" ] }, { diff --git a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/run_script.txt b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/run_script.txt index ff5ef7e568a..b69c5909cb1 100644 --- a/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/run_script.txt +++ b/Tests/Data/Mechanics/MohrCoulombAbboSloanAnisotropic/run_script.txt @@ -34,7 +34,7 @@ import numpy as np import os as os import pyvista as pv import matplotlib.pyplot #In[2]: #Some plot settings - plt.style.use('seaborn-deep') plt.rcParams['lines.linewidth'] = 2.0 plt.rcParams['lines.color'] = 'black' plt.rcParams['legend.frameon'] = True plt.rcParams['font.family'] = 'serif' plt.rcParams['legend.fontsize'] = 14 plt.rcParams['font.size'] = 14 plt.rcParams['axes.spines.right'] = False plt.rcParams['axes.spines.top'] = False plt.rcParams['axes.spines.left'] = True plt.rcParams['axes.spines.bottom'] = True plt.rcParams['axes.axisbelow'] = True plt.rcParams['figure.figsize'] =(12, 6) + plt.style.use('seaborn-v0_8-deep') plt.rcParams['lines.linewidth'] = 2.0 plt.rcParams['lines.color'] = 'black' plt.rcParams['legend.frameon'] = True plt.rcParams['font.family'] = 'serif' plt.rcParams['legend.fontsize'] = 14 plt.rcParams['font.size'] = 14 plt.rcParams['axes.spines.right'] = False plt.rcParams['axes.spines.top'] = False plt.rcParams['axes.spines.left'] = True plt.rcParams['axes.spines.bottom'] = True plt.rcParams['axes.axisbelow'] = True plt.rcParams['figure.figsize'] =(12, 6) #In[3]: diff --git a/Tests/Data/Parabolic/LiquidFlow/BlockingConductingFracture/plot_settings.py b/Tests/Data/Parabolic/LiquidFlow/BlockingConductingFracture/plot_settings.py index 73744d12e9d..a84b0a9c920 100644 --- a/Tests/Data/Parabolic/LiquidFlow/BlockingConductingFracture/plot_settings.py +++ b/Tests/Data/Parabolic/LiquidFlow/BlockingConductingFracture/plot_settings.py @@ -1,7 +1,7 @@ import matplotlib.pyplot as plt # Some plot settings -plt.style.use("seaborn-deep") +plt.style.use("seaborn-v0_8-deep") plt.rcParams["lines.linewidth"] = 2.0 plt.rcParams["lines.color"] = "black" plt.rcParams["legend.frameon"] = True