{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Mesh basics\n", "\n", "In order to solve a model numerically in a region, we have to discretise it. There are two main ways of discretising the space: finite-difference and finite-element discretisation. `discretisedfield` deals only with finite-difference discretisation at the moment. This means that we are dividing our cubic region into smaller \"chunks\" - small cubes. We refer to the discretised region as a **mesh**:\n", "\n", "$$\\text{MESH} = \\text{REGION} + \\text{DISCRETISATION}$$\n", "\n", "In this tutorial, we show how to define it, as well as some basic operations we can perform with meshes.\n", "\n", "As we showed in previous tutorials, region is always cuboidal and it is defined by any two diagonally opposite corner points. We are going to use the same region as before, defined by the following two diagonally opposite points\n", "\n", "$$\\mathbf{p}_{1} = (0, 0, 0)$$\n", "$$\\mathbf{p}_{2} = (l_{x}, l_{y}, l_{z})$$\n", "\n", "with $l_{x} = 100 \\,\\text{nm}$, $l_{y} = 50 \\,\\text{nm}$, and $l_{z} = 20 \\,\\text{nm}$.\n", "\n", "So, let us start by defining the region:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import discretisedfield as df\n", "\n", "p1 = (0, 0, 0)\n", "p2 = (100e-9, 50e-9, 20e-9)\n", "\n", "region = df.Region(p1=p1, p2=p2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The region is now defined. Another missing piece is the discretisation and we need to decide how we are going to discretise the region. In other words, we need to decide into what size \"chunks\" we are going to discretise our region in. We refer to the \"chunk\" as the **discretisation cell**. In `discretisedfield`, there are two ways how we can define the discretisation. We can define either:\n", "\n", "1. The number of discretisation cells in all geometric directions, or\n", "2. The size of a single discretisation cell.\n", "\n", "Let us start with the first case. The number of discretisation cells in all geometric directions can be passed using `n` argument, which is, for a three-dimensional geometry, a length-3 tuple:\n", "\n", "$$n = (n_{x}, n_{y}, n_{z})$$\n", "\n", "For instance, we want to discretise our region in 5 cells in the x-direction, 2 in the y-direction and 1 cell in the z-direction. Therefore, knowing the region as well as the discretisation `n`, we pass them both to `Mesh`:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "n = (5, 2, 1)\n", "\n", "mesh = df.Mesh(region=region, n=n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The mesh is defined. Based on the region dimensions and the number of discretisation cells, we can ask the mesh to give us the size of a single discretisation cell:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2.0e-08, 2.5e-08, 2.0e-08])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh.cell" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Knowing this value, we could have defined the mesh passing this value using `cell` argument, and we would have got exactly the same mesh." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "cell = (20e-9, 25e-9, 20e-9)\n", "\n", "mesh = df.Mesh(region=region, cell=cell)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we now ask our new mesh about the number of discretisation cells:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([5, 2, 1])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh.n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is no difference whatsoever how we are going to define the mesh. However, defining the mesh with `cell` can result in an error, if the region cannot be divided into chunks of that size. For instance:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exception raised.\n" ] } ], "source": [ "try:\n", " mesh = df.Mesh(region=region, cell=(3e-9, 3e-9, 3e-9))\n", "except ValueError:\n", " print(\"Exception raised.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us now have a look at some basic properties we can ask the mesh object for. First of all, region object is a part of the mesh object:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Region\n", "" ], "text/plain": [ "Region(pmin=[0.0, 0.0, 0.0], pmax=[1e-07, 5e-08, 2e-08], dims=['x', 'y', 'z'], units=['m', 'm', 'm'])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh.region" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Therefore, we can perform all the operations on the region we saw previously, but now through the mesh object (`mesh.region`). For instance: " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0., 0., 0.])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh.region.pmin # minimum point" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1.e-07, 5.e-08, 2.e-08])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh.region.edges # edge lenghts" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([5.0e-08, 2.5e-08, 1.0e-08])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh.region.centre # centre point" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By asking the mesh object directly, we can get the number of discretisation cells in all three directions $n = (n_{x}, n_{y}, n_{z})$:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([5, 2, 1])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh.n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and the size of a single discretisation cell:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2.0e-08, 2.5e-08, 2.0e-08])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh.cell" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The total number of discretisation cells is:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(mesh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This number is simply $n_{x}n_{y}n_{z}$. We can conclude that the entire region is now divided into 10 small cubes (discretisation cells). Each cell in the mesh has its index and its coordinate. We can get indices of all discretisation cells:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# NBVAL_IGNORE_OUTPUT\n", "mesh.indices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This gives us a generator object which we can use as an iterable in different Pyhton contexts. For instance, we can give it to the `list`." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(0, 0, 0),\n", " (1, 0, 0),\n", " (2, 0, 0),\n", " (3, 0, 0),\n", " (4, 0, 0),\n", " (0, 1, 0),\n", " (1, 1, 0),\n", " (2, 1, 0),\n", " (3, 1, 0),\n", " (4, 1, 0)]" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(mesh.indices)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "List function now \"unpacks\" the generator and gives us a list of tuples. Each tuple has three unsigned integers. For instance, we can interpret index `(2, 1, 0)` as an index which belongs to the third cell in the x-direction, second in the y, and the first in the z direction. Please note that indexing in Python starts from 0. Therefore, we say that the \"fifth element\" has index 4.\n", "\n", "Another thing we can associate to every discretisation cell is its coordinate. The coordinate of the cell is the coordinate of its centre point. So, the coordinate of index `(2, 1, 0)` cell is:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([5.00e-08, 3.75e-08, 1.00e-08])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "index = (2, 1, 0)\n", "\n", "mesh.index2point(index)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is very often the case we need to iterate through all discretisation cells and use their coordinates. For that, we can use the mesh object itself, which is also an iterable:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[array([1.00e-08, 1.25e-08, 1.00e-08]),\n", " array([3.00e-08, 1.25e-08, 1.00e-08]),\n", " array([5.00e-08, 1.25e-08, 1.00e-08]),\n", " array([7.00e-08, 1.25e-08, 1.00e-08]),\n", " array([9.00e-08, 1.25e-08, 1.00e-08]),\n", " array([1.00e-08, 3.75e-08, 1.00e-08]),\n", " array([3.00e-08, 3.75e-08, 1.00e-08]),\n", " array([5.00e-08, 3.75e-08, 1.00e-08]),\n", " array([7.00e-08, 3.75e-08, 1.00e-08]),\n", " array([9.00e-08, 3.75e-08, 1.00e-08])]" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(mesh)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since mesh object is an iterator itself, we can use it, for example, in for loops:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1.00e-08 1.25e-08 1.00e-08]\n", "[3.00e-08 1.25e-08 1.00e-08]\n", "[5.00e-08 1.25e-08 1.00e-08]\n", "[7.00e-08 1.25e-08 1.00e-08]\n", "[9.00e-08 1.25e-08 1.00e-08]\n", "[1.00e-08 3.75e-08 1.00e-08]\n", "[3.00e-08 3.75e-08 1.00e-08]\n", "[5.00e-08 3.75e-08 1.00e-08]\n", "[7.00e-08 3.75e-08 1.00e-08]\n", "[9.00e-08 3.75e-08 1.00e-08]\n" ] } ], "source": [ "for point in mesh:\n", " print(point)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A function, which is opposite to `index2point`, is `point2index`. This function takes any point in the region and returns the index of a cell it belongs to:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2, 1, 0)" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "point = (41.6e-9, 35.2e-9, 4.71e-9)\n", "\n", "mesh.point2index(point)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also ask the mesh to give us the midpoints along a certain axis:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1e-08, 3.0000000000000004e-08, 5e-08, 7e-08, 9e-08]" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(mesh.cells.x)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1.25e-08, 3.75e-08]" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(mesh.cells.y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, we can get the vertices of the cells along a certain axis:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0.0, 2e-08, 4e-08, 6.000000000000001e-08, 8e-08, 1e-07]" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(mesh.vertices.x)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0.0, 2.5e-08, 5e-08]" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(mesh.vertices.y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can compare meshes using `==` and `!=` relational operators. Let us define two meshes and compare them to the one we have already defined:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh_same = df.Mesh(region=region, n=(5, 2, 1))\n", "mesh_different = df.Mesh(region=region, n=(10, 5, 7))\n", "\n", "mesh == mesh_same" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh == mesh_different" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh != mesh_different" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, mesh has its representation string:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Mesh\n", "" ], "text/plain": [ "Mesh(Region(pmin=[0.0, 0.0, 0.0], pmax=[1e-07, 5e-08, 2e-08], dims=['x', 'y', 'z'], units=['m', 'm', 'm']), n=[5, 2, 1])" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the representation string, we see `pmin`, `pmax`, and `n` we discussed earlier, but there are also `bc` and `subregions` we did not and we will look at some more advanced mesh properties in the next tutorials." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Scale and translate\n", "Similarly to regions, a mesh has two methods `scale` and `translate`. Both optionally allow `inplace` operations." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Mesh\n", "" ], "text/plain": [ "Mesh(Region(pmin=[-5e-08, -2.5e-08, -1e-08], pmax=[1.5e-07, 7.5e-08, 3.0000000000000004e-08], dims=['x', 'y', 'z'], units=['m', 'm', 'm']), n=[5, 2, 1])" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scaled = mesh.scale(2)\n", "scaled" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Mesh\n", "" ], "text/plain": [ "Mesh(Region(pmin=[-10.0, 0.0, 100.0], pmax=[-9.9999999, 5e-08, 100.00000002], dims=['x', 'y', 'z'], units=['m', 'm', 'm']), n=[5, 2, 1])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "translated = mesh.translate((-10, 0, 100))\n", "translated" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can scale with different factors along different directions and pass in a `reference_point` that does not move when scaling. If not passed the `center` of the mesh is kept fixed." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Mesh\n", "" ], "text/plain": [ "Mesh(Region(pmin=[0.0, -2.5, 0.0], pmax=[1e-07, -2.4999999625, 4e-08], dims=['x', 'y', 'z'], units=['m', 'm', 'm']), n=[5, 2, 1])" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scaled = mesh.scale((1, 0.75, 2), reference_point=(-10, -10, 0))\n", "scaled" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is also possible to scale the mesh `inplace`." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Mesh\n", "" ], "text/plain": [ "Mesh(Region(pmin=[0.0, 6.25e-09, -1e-08], pmax=[1e-07, 4.375e-08, 3.0000000000000004e-08], dims=['x', 'y', 'z'], units=['m', 'm', 'm']), n=[5, 2, 1])" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh.scale((1, 0.75, 2), inplace=True)\n", "mesh" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mesh alignment\n", "We can also check if two meshes are aligned using the `is_aligned` function.\n", "\n", "Two meshes are considered to be aligned if and only if:\n", "\n", " 1. They have same discretisation cell size.\n", "\n", " 2. They have common cell coordinates.\n", "\n", "for a given tolerance value.\n", "\n", "Let's create two meshes which have two different regions which overlap but the cells centers are located in the same places.\n", "\n", "These two meshes are aligned." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p1 = (-50, -25, 0)\n", "p2 = (50, 25, 5)\n", "cell = (5, 5, 5)\n", "region1 = df.Region(p1=p1, p2=p2)\n", "mesh1 = df.Mesh(region=region1, cell=cell)\n", "\n", "p1 = (-45, -20, 0)\n", "p2 = (10, 20, 5)\n", "cell = (5, 5, 5)\n", "region2 = df.Region(p1=p1, p2=p2)\n", "mesh2 = df.Mesh(region=region2, cell=cell)\n", "\n", "mesh1.is_aligned(mesh2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we create a third mesh where the cell centers are not in the same places then this mesh is not aligned." ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p1 = (-42, -20, 0)\n", "p2 = (13, 20, 5)\n", "cell = (5, 5, 5)\n", "region3 = df.Region(p1=p1, p2=p2)\n", "mesh3 = df.Mesh(region=region3, cell=cell)\n", "\n", "mesh1.is_aligned(mesh3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, we can provide the tolerance as an argument so dictate how strict the `is_aligned` check is." ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mesh1.is_aligned(mesh3, tolerance=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Other dimensions\n", "\n", "Similar to regions, a mesh can have any dimension >= 1. Plotting is currently only supported for 3d regions.\n", "\n", "For example, a 1d mesh" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Mesh\n", "" ], "text/plain": [ "Mesh(Region(pmin=[0], pmax=[10], dims=['x'], units=['m']), n=[10])" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "region = df.Region(p1=0, p2=10)\n", "mesh = df.Mesh(region=region, n=10)\n", "mesh" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or a 4d mesh" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "Mesh\n", "" ], "text/plain": [ "Mesh(Region(pmin=[0, 0, 0, 0], pmax=[10, 20, 30, 40], dims=['x0', 'x1', 'x2', 'x3'], units=['m', 'm', 'm', 'm']), n=[2, 3, 4, 5])" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "region = df.Region(p1=(0, 0, 0, 0), p2=(10, 20, 30, 40))\n", "mesh = df.Mesh(region=region, n=(2, 3, 4, 5))\n", "mesh" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.17" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }