Mesh basics¶

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:

$$\text{MESH} = \text{REGION} + \text{DISCRETISATION}$$

In this tutorial, we show how to define it, as well as some basic operations we can perform with meshes.

As we showed in previous tutorials, region is always cubic 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

$$\mathbf{p}_{1} = (0, 0, 0)$$$$\mathbf{p}_{2} = (l_{x}, l_{y}, l_{z})$$

with $l_{x} = 100 \,\text{nm}$, $l_{y} = 50 \,\text{nm}$, and $l_{z} = 20 \,\text{nm}$.

So, let us start by defining the region:

In [1]:

import discretisedfield as df

p1 = (0, 0, 0)
p2 = (100e-9, 50e-9, 20e-9)

region = df.Region(p1=p1, p2=p2)

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:

The number of discretisation cells in all 3 ($x$, $y$, and $z$) directions, or
The size of a single discretisation cell.

Let us start with the first case. The number of discretisation cells in all three directions can be passed using n argument, which is a length-3 tuple:

$$n = (n_{x}, n_{y}, n_{z})$$

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:

In [2]:

n = (5, 2, 1)

mesh = df.Mesh(region=region, n=n)

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:

In [3]:

mesh.cell

Out[3]:

(2e-08, 2.5e-08, 2e-08)

Knowing this value, we could have defined the mesh passing this value using cell argument, and we would have got exactly the same mesh.

In [4]:

cell = (20e-9, 25e-9, 20e-9)

mesh = df.Mesh(region=region, cell=cell)

If we now ask our new mesh about the number of discretisation cells:

In [5]:

mesh.n

Out[5]:

(5, 2, 1)

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:

In [6]:

try:
    mesh = df.Mesh(region=region, cell=(3e-9, 3e-9, 3e-9))
except ValueError:
    print("Exception raised.")

Exception raised.

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:

In [7]:

mesh.region

Out[7]:

Region

p1 = (0, 0, 0)
p2 = (1e-07, 5e-08, 2e-08)

Therefore, we can perform all the operations on the region we saw previously, but now through the mesh object (mesh.region). For instance:

In [8]:

mesh.region.pmin  # minimum point

Out[8]:

(0.0, 0.0, 0.0)

In [9]:

mesh.region.edges  # edge lenghts

Out[9]:

(1e-07, 5e-08, 2e-08)

In [10]:

mesh.region.centre  # centre point

Out[10]:

(5e-08, 2.5e-08, 1e-08)

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})$:

In [11]:

mesh.n

Out[11]:

(5, 2, 1)

and the size of a single discretisation cell:

In [12]:

mesh.cell

Out[12]:

(2e-08, 2.5e-08, 2e-08)

The total number of discretisation cells is:

In [13]:

len(mesh)

Out[13]:

This number is simply $n_{x}n_{y}n_{z}$. We can conclude that the entire region is now divided into 100 small cubes (discretisation cells). Each cell in the mesh has its index and its coordinate. We can get indices of all discretisation cells:

In [14]:

# NBVAL_IGNORE_OUTPUT
mesh.indices

Out[14]:

<generator object Mesh.indices at 0x7f8d19a7df90>

This gives us a generator object we can use as an iterable in different Pyhton contexts. For instance, we can give it to the list.

In [15]:

list(mesh.indices)

Out[15]:

[(0, 0, 0),
 (1, 0, 0),
 (2, 0, 0),
 (3, 0, 0),
 (4, 0, 0),
 (0, 1, 0),
 (1, 1, 0),
 (2, 1, 0),
 (3, 1, 0),
 (4, 1, 0)]

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.

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:

In [16]:

index = (2, 1, 0)

mesh.index2point(index)

Out[16]:

(5e-08, 3.75e-08, 1e-08)

We can also get coordinates of neighbouring cells.

In [17]:

mesh.neighbours((0, 0, 0))

Out[17]:

[(1, 0, 0), (0, 1, 0)]

In [18]:

mesh.neighbours((1, 1, 0))

Out[18]:

[(0, 1, 0), (2, 1, 0), (1, 0, 0)]

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:

In [19]:

list(mesh)

Out[19]:

[(1e-08, 1.25e-08, 1e-08),
 (3.0000000000000004e-08, 1.25e-08, 1e-08),
 (5e-08, 1.25e-08, 1e-08),
 (7e-08, 1.25e-08, 1e-08),
 (9e-08, 1.25e-08, 1e-08),
 (1e-08, 3.75e-08, 1e-08),
 (3.0000000000000004e-08, 3.75e-08, 1e-08),
 (5e-08, 3.75e-08, 1e-08),
 (7e-08, 3.75e-08, 1e-08),
 (9e-08, 3.75e-08, 1e-08)]

Since mesh object is an iterator itself, we can use it, for example, in for loops:

In [20]:

for point in mesh:
    print(point)

(1e-08, 1.25e-08, 1e-08)
(3.0000000000000004e-08, 1.25e-08, 1e-08)
(5e-08, 1.25e-08, 1e-08)
(7e-08, 1.25e-08, 1e-08)
(9e-08, 1.25e-08, 1e-08)
(1e-08, 3.75e-08, 1e-08)
(3.0000000000000004e-08, 3.75e-08, 1e-08)
(5e-08, 3.75e-08, 1e-08)
(7e-08, 3.75e-08, 1e-08)
(9e-08, 3.75e-08, 1e-08)

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:

In [21]:

point = (41.6e-9, 35.2e-9, 4.71e-9)

mesh.point2index(point)

Out[21]:

(2, 1, 0)

We can also ask the mesh to give us the midpoints along a certain axis:

In [22]:

list(mesh.midpoints.x)

Out[22]:

[1e-08, 3.0000000000000004e-08, 5e-08, 7e-08, 9e-08]

In [23]:

list(mesh.midpoints.y)

Out[23]:

[1.25e-08, 3.75e-08]

Similarly, we can get the vertices of the cells along a certain axis:

In [24]:

list(mesh.vertices.x)

Out[24]:

[0.0, 2e-08, 4e-08, 6.000000000000001e-08, 8e-08, 1e-07]

In [25]:

list(mesh.vertices.y)

Out[25]:

[0.0, 2.5e-08, 5e-08]

We can compare meshes using == and != relational operators. Let us define two meshes and compare them to the one we have already defined:

In [26]:

mesh_same = df.Mesh(region=region, n=(5, 2, 1))
mesh_different = df.Mesh(region=region, n=(10, 5, 7))

mesh == mesh_same

Out[26]:

True

In [27]:

mesh == mesh_different

Out[27]:

False

In [28]:

mesh != mesh_different

Out[28]:

True

Finally, mesh has its representation string:

In [29]:

mesh

Out[29]:

Mesh

Region
- p1 = (0, 0, 0)
- p2 = (1e-07, 5e-08, 2e-08)
n = (5, 2, 1)
attributes:
- unit: m
- fourierspace: False
- isplane: False

In the representation string, we see p1, p2, 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.