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 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

\[\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:

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:

1. The number of discretisation cells in all geometric directions, or

2. The size of a single discretisation cell.

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_{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:

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:

mesh.cell

array([2.0e-08, 2.5e-08, 2.0e-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.

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:

mesh.n

array([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:

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:

mesh.region

Region

• pmin = [0.0, 0.0, 0.0]
• pmax = [1e-07, 5e-08, 2e-08]
• dims = ['x', 'y', 'z']
• units = ['m', 'm', 'm']

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

mesh.region.pmin  # minimum point

array([0., 0., 0.])

mesh.region.edges  # edge lenghts

array([1.e-07, 5.e-08, 2.e-08])

mesh.region.centre  # centre point

array([5.0e-08, 2.5e-08, 1.0e-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})\):

mesh.n

array([5, 2, 1])

and the size of a single discretisation cell:

mesh.cell

array([2.0e-08, 2.5e-08, 2.0e-08])

The total number of discretisation cells is:

len(mesh)

10

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:

# NBVAL_IGNORE_OUTPUT
mesh.indices

<generator object Mesh.indices at 0x7f187843d3c0>

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.

list(mesh.indices)

[(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:

index = (2, 1, 0)

mesh.index2point(index)

array([5.00e-08, 3.75e-08, 1.00e-08])

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:

list(mesh)

[array([1.00e-08, 1.25e-08, 1.00e-08]),
 array([3.00e-08, 1.25e-08, 1.00e-08]),
 array([5.00e-08, 1.25e-08, 1.00e-08]),
 array([7.00e-08, 1.25e-08, 1.00e-08]),
 array([9.00e-08, 1.25e-08, 1.00e-08]),
 array([1.00e-08, 3.75e-08, 1.00e-08]),
 array([3.00e-08, 3.75e-08, 1.00e-08]),
 array([5.00e-08, 3.75e-08, 1.00e-08]),
 array([7.00e-08, 3.75e-08, 1.00e-08]),
 array([9.00e-08, 3.75e-08, 1.00e-08])]

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

for point in mesh:
    print(point)

[1.00e-08 1.25e-08 1.00e-08]
[3.00e-08 1.25e-08 1.00e-08]
[5.00e-08 1.25e-08 1.00e-08]
[7.00e-08 1.25e-08 1.00e-08]
[9.00e-08 1.25e-08 1.00e-08]
[1.00e-08 3.75e-08 1.00e-08]
[3.00e-08 3.75e-08 1.00e-08]
[5.00e-08 3.75e-08 1.00e-08]
[7.00e-08 3.75e-08 1.00e-08]
[9.00e-08 3.75e-08 1.00e-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:

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

mesh.point2index(point)

(2, 1, 0)

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

list(mesh.cells.x)

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

list(mesh.cells.y)

[1.25e-08, 3.75e-08]

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

list(mesh.vertices.x)

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

list(mesh.vertices.y)

[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:

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

mesh == mesh_same

True

mesh == mesh_different

False

mesh != mesh_different

True

Finally, mesh has its representation string:

mesh

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]

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.

Scale and translate

Similarly to regions, a mesh has two methods scale and translate. Both optionally allow inplace operations.

scaled = mesh.scale(2)
scaled

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]

translated = mesh.translate((-10, 0, 100))
translated

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]

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.

scaled = mesh.scale((1, 0.75, 2), reference_point=(-10, -10, 0))
scaled

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]

It is also possible to scale the mesh inplace.

mesh.scale((1, 0.75, 2), inplace=True)
mesh

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]

Mesh alignment

We can also check if two meshes are aligned using the is_aligned function.

Two meshes are considered to be aligned if and only if:

1. They have same discretisation cell size.

2. They have common cell coordinates.

for a given tolerance value.

Let’s create two meshes which have two different regions which overlap but the cells centers are located in the same places.

These two meshes are aligned.

p1 = (-50, -25, 0)
p2 = (50, 25, 5)
cell = (5, 5, 5)
region1 = df.Region(p1=p1, p2=p2)
mesh1 = df.Mesh(region=region1, cell=cell)

p1 = (-45, -20, 0)
p2 = (10, 20, 5)
cell = (5, 5, 5)
region2 = df.Region(p1=p1, p2=p2)
mesh2 = df.Mesh(region=region2, cell=cell)

mesh1.is_aligned(mesh2)

True

If we create a third mesh where the cell centers are not in the same places then this mesh is not aligned.

p1 = (-42, -20, 0)
p2 = (13, 20, 5)
cell = (5, 5, 5)
region3 = df.Region(p1=p1, p2=p2)
mesh3 = df.Mesh(region=region3, cell=cell)

mesh1.is_aligned(mesh3)

False

However, we can provide the tolerance as an argument so dictate how strict the is_aligned check is.

mesh1.is_aligned(mesh3, tolerance=2)

True

Other dimensions

Similar to regions, a mesh can have any dimension >= 1. Plotting is currently only supported for 3d regions.

For example, a 1d mesh

region = df.Region(p1=0, p2=10)
mesh = df.Mesh(region=region, n=10)
mesh

Mesh

• Region
  • pmin = [0]
  • pmax = [10]
  • dims = ['x']
  • units = ['m']
• n = [10]

Or a 4d mesh

region = df.Region(p1=(0, 0, 0, 0), p2=(10, 20, 30, 40))
mesh = df.Mesh(region=region, n=(2, 3, 4, 5))
mesh

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]