Field operations#

There are several convenience methods that can be used to analyse the field. Let us first define a mesh with three spatial dimensions (ndim=3) which we can work with.

[1]:
import numpy as np

import discretisedfield as df

p1 = (-50, -50, -50)
p2 = (50, 50, 50)
n = (4, 2, 2)
mesh = df.Mesh(p1=p1, p2=p2, n=n)

We are going to initialise a three dimensional vector field (nvdim=3), with

\[\mathbf{f}(x, y, z) = (xy, 2xy, xyz)\]

For that, we are going to use the following Python function.

[2]:
def value_function(pos):
    x, y, z = pos
    return x * y, 2 * x * y, x * y * z

Finally, our field is

[3]:
field = df.Field(mesh, nvdim=3, value=value_function)

Sampling the field#

As we have shown previously, a field can be sampled by calling it. The argument must be a iterable with the same length as mesh.ndim which in our case is three - as we are working in three spatial dimensions. The coordinates of the point must also be contained inside the region.

[4]:
point = (0, 0, 0)
field(point)
[4]:
array([ 312.5,  625. , 7812.5])

However, if the point is outside the mesh, an exception is raised.

[5]:
point = (100, 100, 100)
try:
    field(point)
except ValueError:
    print("Exception raised.")
Exception raised.

Extracting the component of a vector field#

A three-dimensional vector field can be understood as three separate scalar fields, where each scalar field is a component of a vector field value. A scalar field of a component can be extracted by accessing the associated vdims attribute of the field. In our case this is x, y, or z.

[6]:
x_component = field.x
x_component((0, 0, 0))
[6]:
array([312.5])
[7]:
field.y
[7]:
Field
  • Mesh
    • Region
      • pmin = [-50, -50, -50]
      • pmax = [50, 50, 50]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [4, 2, 2]
  • nvdim = 1

Default names x, y, and z are only available for fields with dimensionality (nvdim) 2 or 3.

[8]:
field.vdims
[8]:
['x', 'y', 'z']

It is possible to change the component names:

[9]:
field.vdims = ["mx", "my", "mz"]
field.mx((0, 0, 0))
[9]:
array([312.5])

This overrides the component labels and the old x, y and z cannot be used anymore:

[10]:
try:
    field.x
except AttributeError as e:
    print(e)
Object has no attribute x.

We change the component labels back to x, y, and z for the rest of this notebook.

[11]:
field.vdims = ["x", "y", "z"]

Custom component names can optionally also be specified during field creation. If not specified, the default values are used for fields with dimensions 2 or 3. Higher-dimensional fields have no defaults v0, v1, v2, v3, etc.:

[12]:
field_4d = df.Field(mesh, nvdim=4, value=[1, 1, 1, 1])
field_4d
[12]:
Field
  • Mesh
    • Region
      • pmin = [-50, -50, -50]
      • pmax = [50, 50, 50]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [4, 2, 2]
  • nvdim = 4
  • vdims:
    • v0
    • v1
    • v2
    • v3
[13]:
field_4d.v1((0, 0, 0))
[13]:
array([1.])

Extracting smaller region#

Let us say we are not interested in the entire field but only in its smaller portion - only some discretisation cells. In that case, we have two options. Before we discuss them, let us first define what we mean by “aligned” meshes:

  • Mesh A is aligned to mesh B if and only if all cell coordinates of mesh A are also the coordinates of (some) cells in mesh B.

There is is_aligned method which checks that. Let us have a look at a few meshes:

[14]:
mesh1 = df.Mesh(region=df.Region(p1=(0, 0, 0), p2=(10, 10, 10)), cell=(1, 1, 1))
mesh2 = df.Mesh(region=df.Region(p1=(3, 3, 3), p2=(6, 6, 6)), cell=(1, 1, 1))
mesh3 = df.Mesh(region=df.Region(p1=(0, 0, 0), p2=(10, 10, 10)), cell=(2, 2, 2))
mesh4 = df.Mesh(
    region=df.Region(p1=(3.5, 3.5, 3.5), p2=(6.5, 6.5, 6.5)), cell=(1, 1, 1)
)

Let us now have a look if those meshes are aligned:

[15]:
mesh1.is_aligned(mesh2)
[15]:
True
[16]:
mesh1.is_aligned(mesh3)  # discretisation cell is different
[16]:
False
[17]:
mesh1.is_aligned(
    mesh4
)  # although discretisation cell is the same, mesh4 is shifted in space by 0.5
[17]:
False

Extracting subfield on aligned mesh#

If we want to get a subfield whose mesh is aligned to the field we want to take part of, we use [] operator. The resulting field is going to have a minimum-sized mesh which contains the region we pass as an argument.

[18]:
subregion = df.Region(p1=(1.5, 2.2, 3.9), p2=(6.1, 5.9, 9.9))
field[subregion]
[18]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [25.0, 50.0, 50.0]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [1, 1, 1]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

We can see that the resulting field’s mesh has the minimum dimesions aligned mesh should have in order to contain the subregion. The resulting field has the same discretisation cell as the original one.

Extracting field on any mesh#

If we want to extract part of the field on any mesh which is contained inside the field, we do that by “resampling”. We create a new field on a submesh and pass the field we want take subfield from as value.

[19]:
subregion = df.Region(p1=(1.5, 2.5, 3.5), p2=(5.5, 5.5, 6.5))
submesh = df.Mesh(region=subregion, cell=(0.5, 0.5, 0.5))
df.Field(submesh, nvdim=3, value=field)
[19]:
Field
  • Mesh
    • Region
      • pmin = [1.5, 2.5, 3.5]
      • pmax = [5.5, 5.5, 6.5]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [8, 6, 6]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

One could ask why don’t we always use resampling because it is a generalised case. The reason is because computing a subfield using [] operator is much faster.

Computing the mean#

The mean of the field can be obtained by calling discretisedfield.Field.mean function.

[20]:
field.mean()
[20]:
array([0., 0., 0.])

Mean always return a numpy array, independent of the dimension of the field’s value.

[21]:
field.x.mean()
[21]:
array([0.])

Iterating through the field#

The field object itself is an iterable. That means that it can be iterated through. As a result it returns the values stored in the field.

[22]:
for val in field:
    print(val)
[   937.5   1875.  -23437.5]
[  312.5   625.  -7812.5]
[-312.5 -625.  7812.5]
[ -937.5 -1875.  23437.5]
[ -937.5 -1875.  23437.5]
[-312.5 -625.  7812.5]
[  312.5   625.  -7812.5]
[   937.5   1875.  -23437.5]
[  937.5  1875.  23437.5]
[ 312.5  625.  7812.5]
[ -312.5  -625.  -7812.5]
[  -937.5  -1875.  -23437.5]
[  -937.5  -1875.  -23437.5]
[ -312.5  -625.  -7812.5]
[ 312.5  625.  7812.5]
[  937.5  1875.  23437.5]

Sampling the field along the line#

To sample the points of the field which are on a certain line, discretisedfield.Field.line method is used. It takes two points p1 and p2 that define the line and an integer n which defines how many mesh coordinates on that line are required. The default value of n is 100.

[23]:
line = field.line(p1=(-10, 0, 0), p2=(10, 0, 0), n=5)

Intersecting the field with a plane#

If we intersect the field with a plane, discretisedfield.Field.sel can be used to return a new field object which contains only discretisation cells that belong to a ndim - 1 plane. The planes allowed are the planes perpendicular to the geometric axes. For instance, for a three-dimensional geometry, a plane parallel to the \(yz\)-plane (perpendicular to the \(x\)-axis) which intersects the \(x\)-axis at 1, can be written as

\[x = 1\]
[24]:
field.sel(x=1)
[24]:
Field
  • Mesh
    • Region
      • pmin = [-50, -50]
      • pmax = [50, 50]
      • dims = ['y', 'z']
      • units = ['m', 'm']
    • n = [2, 2]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

If we want to cut through the middle of the mesh, we do not need to provide a particular value for a coordinate.

[25]:
field.sel("x")
[25]:
Field
  • Mesh
    • Region
      • pmin = [-50, -50]
      • pmax = [50, 50]
      • dims = ['y', 'z']
      • units = ['m', 'm']
    • n = [2, 2]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

Selecting part of the field#

We can also select part of the field using, discretisedfield.Field.sel. This can be used to return a new field object which contains only discretisation cells in the specific region. The dimensionality of the field remains the same. The selections hat are allowed are the geometric axes. For instance, for a three-dimensional geometry, we can select the \(x\)-axis between -10 and 30. This can be written as

\[x = (-10, 30)\]

This will return three cells in the \(x\) direction that between \(-25\) and \(50\). This is as there are originally four cells which have boundaries:

   Cell 0     Cell 1     Cell 2     Cell 3
|----------|----------|----------|----------|
-50       -25         0         25         50

If we select our lower point at \(-10\) then the whole of Cell 1 will be selected and if we select our upper point at \(+10\) then the whole of Cell 3 will be selected. All the cells between these two points will also be selected.

[26]:
field.sel(x=(-10, 30))
[26]:
Field
  • Mesh
    • Region
      • pmin = [-25.0, -50.0, -50.0]
      • pmax = [50.0, 50.0, 50.0]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [3, 2, 2]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

Cascading the operations#

Let us say we want to compute the mean() of an \(x\) component of the field on the plane \(y=10\). In order to do that, we can cascade several operation in a single line.

[27]:
field.sel(y=10).x.mean()
[27]:
array([0.])

This gives the same result as for instance

[28]:
field.x.sel(y=10).mean()
[28]:
array([0.])

Complex fields#

discretisedfield supports complex-valued fields.

[29]:
cfield = df.Field(mesh, nvdim=3, value=(1 + 1.5j, 2, 3j))

We can extract real and imaginary part.

[30]:
cfield.real((0, 0, 0))
[30]:
array([1., 2., 0.])
[31]:
cfield.imag((0, 0, 0))
[31]:
array([1.5, 0. , 3. ])

Similarly we get real and imaginary parts of individual components.

[32]:
cfield.x.real((0, 0, 0))
[32]:
array([1.])
[33]:
cfield.x.imag((0, 0, 0))
[33]:
array([1.5])

Complex conjugate.

[34]:
cfield.conjugate((0, 0, 0))
[34]:
array([1.-1.5j, 2.-0.j , 0.-3.j ])

Phase in the complex plane.

[35]:
cfield.phase((0, 0, 0))
[35]:
array([0.98279372, 0.        , 1.57079633])

Algebra operations#

Let us define two fields:

[36]:
region = df.Region(p1=(0, 0, 0), p2=(10e-9, 10e-9, 10e-9))
mesh = df.Mesh(region=region, n=(10, 10, 10))

f1 = df.Field(mesh, nvdim=3, value=(1, 1, 1))
f2 = df.Field(mesh, nvdim=3, value=(2, 1, 3))
[37]:
f1.mean()
[37]:
array([1., 1., 1.])
[38]:
f2.mean()
[38]:
array([2., 1., 3.])

+ operation#

[39]:
f1 + f2
[39]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z
[40]:
(f1 + f2).mean()
[40]:
array([3., 2., 4.])

- operation#

[41]:
f1 - f2
[41]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z
[42]:
(f1 - f2).mean()
[42]:
array([-1.,  0., -2.])

* operation#

Basic multiplication between vector fields is carried out element wise.

[43]:
f1 * f2  # both are vector fields
[43]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

Scalar with vector field:

[44]:
f1.x * f2
[44]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

Scalar with vector field:

[45]:
f1.x * f2.y
[45]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 1

/ operation#

Dividing vector field by a vector field:

[46]:
f1 / f2
[46]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

Dividing vector field by a scalar field:

[47]:
f1 / f2.x
[47]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

We can also divide a scalar field by a vector field:

[48]:
f2.x / f1
[48]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

** operator#

This operator is applied element wise to all fields

[49]:
f1**2
[49]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z
[50]:
f1.x**2
[50]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 1

Compund operations#

[51]:
f1 += f2
[52]:
f1 -= f2
[53]:
f1 *= f2.x
[54]:
f2 /= f2.y

Vector products#

As the title says, these products are applied between vector fields only.

Dot product#

Dot product is implemented through dot method:

[55]:
f1.dot(f2)
[55]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 1

Cross product#

Cross product between vector fields is performed using cross method:

[56]:
f1.cross(f2)
[56]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

Vector calculus#

Directional derivative \(\left(\frac{\partial}{\partial x_{i}}f\right)\)#

Defined on both scalar and vector fields:

[57]:
f1.diff("x")
[57]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

Gradient \((\nabla f)\)#

Defined on scalar fields:

[58]:
f1.x.grad
[58]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

Divergence \((\nabla \cdot f)\)#

Defined on nvdim=2 or nvdim=3 vector fields:

[59]:
f1.div
[59]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 1

Curl \((\nabla \times f)\)#

Defined on nvdim=3 vector fields:

[60]:
f1.curl
[60]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

Laplace operator \((\nabla^{2} f)\)#

Defined on both vector (nvdim=3) and scalar (nvdim=1) fields:

[61]:
f1.laplace
[61]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z
[62]:
f1.x.laplace
[62]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-08, 1e-08, 1e-08]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 10, 10]
  • nvdim = 1

Integrals#

Computing integrals is performed using integrate method. We are going to show several examples which will give you an idea how you can compute different integrals.

Let us first create a field:

[63]:
import discretisedfield as df

p1 = (0, 0, 0)
p2 = (100e-9, 100e-9, 100e-9)
cell = (2e-9, 2e-9, 2e-9)
region = df.Region(p1=p1, p2=p2)
mesh = df.Mesh(region=region, cell=cell)
f = df.Field(mesh, nvdim=3, value=(-3, 0, 4), norm=1e6)

Volume integral#

\[\int_{V}\mathbf{f}(\mathbf{r})\text{d}V\]
[64]:
f.integrate()
[64]:
array([-6.e-16,  0.e+00,  8.e-16])

Surface integral#

There is disretisedfield.dS value which is a vector field perpendicular to the surface with magnitude equal to the area of \(\text{d}S\).

\[\int_{S}\mathbf{f}(\mathbf{r}) \cdot \text{d}\mathbf{S}\]

There are two common ways to do surface integrals.

[65]:
f.sel("z").z.integrate()
[65]:
array([8.e-09])
[66]:
f.dot((0, 0, 1)).sel("z").integrate()
[66]:
array([8.e-09])

Line integrals#

\[\int_{0}^{x_\text{max}}\mathbf{f}(\mathbf{r}) \text{d}x\]
[67]:
f.integrate("x")
[67]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0]
      • pmax = [1e-07, 1e-07]
      • dims = ['y', 'z']
      • units = ['m', 'm']
    • n = [50, 50]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z
\[\int_{0}^{y_\text{max}}f_{x}(\mathbf{r}) \text{d}y\]
[68]:
f.x.integrate("y")
[68]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0]
      • pmax = [1e-07, 1e-07]
      • dims = ['x', 'z']
      • units = ['m', 'm']
    • n = [50, 50]
  • nvdim = 1

Example#

We have showed how to compute an integral when integrand is just a field. It is important to have in mind that this can be any field after some operations have been applied on it. For instance:

\[\int_{V}\nabla\cdot\mathbf{f}(\mathbf{r})\text{d}V\]
[69]:
f.div.integrate()
[69]:
array([0.])

Here we implement skyrmion number calculations using operations on fields:

\[S = \frac{1}{4\pi} \int \mathbf{m} \cdot \left(\frac{\partial \mathbf{m}}{\partial x} \times \frac{\partial \mathbf{m}}{\partial y}\right) dxdy\]
[70]:
import math

m = field.orientation.sel("z")
S = m.dot(m.diff("x").cross(m.diff("y"))).integrate() / (4 * math.pi)
S
[70]:
array([0.])

Or using Ubermag function:

[71]:
import discretisedfield.tools as dft

dft.topological_charge(m)
[71]:
0.0

Using Berg-Luescher method

[72]:
dft.topological_charge(m, method="berg-luescher")
[72]:
0.0

Resampling the field#

It is possible to resample the field on a new mesh while keeping pmin and pmax stay unchanged using the resample function with the desired number of new number cells passed as an argument.

The values of the new cells are taken from the nearest old cell and no interpolation is performed.

[73]:
f.resample((10, 9, 5))
[73]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-07, 1e-07, 1e-07]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [10, 9, 5]
  • nvdim = 3
  • vdims:
    • x
    • y
    • z

Applying numpys universal functions#

numpy universal functions can be applied to discretisedfield.Field objects. Below we show a different examples. For available functions please refer to the numpy documentation.

[75]:
f1 = df.Field(mesh, nvdim=1, value=1)
f2 = df.Field(mesh, nvdim=1, value=np.pi)
f3 = df.Field(mesh, nvdim=1, value=2)
[76]:
np.sin(f1)
[76]:
Field
  • Mesh
    • Region
      • pmin = [0.0, 0.0, 0.0]
      • pmax = [1e-07, 1e-07, 1e-07]
      • dims = ['x', 'y', 'z']
      • units = ['m', 'm', 'm']
    • n = [50, 50, 50]
  • nvdim = 1
[77]:
np.sin(f2)((0, 0, 0))
[77]:
array([1.2246468e-16])
[78]:
np.sum((f1, f2, f3))((0, 0, 0))
[78]:
array([6.14159265])
[79]:
np.exp(f1)((0, 0, 0))
[79]:
array([2.71828183])
[80]:
np.power(f3, 2)((0, 0, 0))
[80]:
array([4.])

Other#

Full description of all existing functionality can be found in the API Reference.