#!/usr/bin/env python
# coding: utf-8

# # Field operations
# 
# There are several convenience methods that can be used to analyse the field. Let us first define the mesh we are going to work with.

# In[1]:


import numpy as np

import discretisedfield as df

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


# We are going to initialise the vector field (`dim=3`), with
# 
# $$\mathbf{f}(x, y, z) = (xy, 2xy, xyz)$$
# 
# For that, we are going to use the following Python function.

# In[2]:


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


# Finally, our field is

# In[3]:


field = df.Field(mesh, dim=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 3-length iterable and it contains the coordinates of the point.

# In[4]:


point = (0, 0, 0)
field(point)


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

# In[5]:


point = (100, 100, 100)
try:
    field(point)
except ValueError:
    print("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 `x`, `y`, or `z` attribute of the field.

# In[6]:


x_component = field.x
x_component((0, 0, 0))


# In[7]:


field.y


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

# In[8]:


field.components


# It is possible to change the component names:

# In[9]:


field.components = ["mx", "my", "mz"]
field.mx((0, 0, 0))


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

# In[10]:


try:
    field.x
except AttributeError as e:
    print(e)


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

# In[11]:


field.components = ["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 and custom labes have to be specified in order to access individual field components:

# In[12]:


field_4d = df.Field(
    mesh, dim=4, value=[1, 1, 1, 1], components=["c1", "c2", "c3", "c4"]
)
field_4d


# In[13]:


field_4d.c1((0, 0, 0))


# ## 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 `|` operator which checks that. Let us have a look at a few meshes:

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

# In[15]:


mesh1 | mesh2


# In[16]:


mesh1 | mesh3  # discretisation cell is different


# In[17]:


mesh1 | mesh4  # although discretisation cell is the same, mesh4 is shifted in space by 0.5


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

# In[18]:


subregion = df.Region(p1=(1.5, 2.2, 3.9), p2=(6.1, 5.9, 9.9))
field[subregion]


# 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 extact 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`.

# In[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, dim=3, value=field)


# 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 average
# 
# The average of the field can be obtained by calling `discretisedfield.Field.average` property.

# In[20]:


field.average


# Average always return a tuple, independent of the dimension of the field's value.

# In[21]:


field.x.average


# ## Iterating through the field
# 
# The field object itself is an iterable. That means that it can be iterated through. As a result it returns a tuple, where the first element is the coordinate of the mesh point, whereas the second one is its value.

# In[22]:


for coordinate, value in field:
    print(coordinate, value)


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

# In[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.plane` will return a new field object which contains only discretisation cells that belong to that plane. The planes allowed are the planes perpendicular to the axes of the Cartesian coordinate system. For instance, a plane parallel to the $yz$-plane (perpendicular to the $x$-axis) which intesects the $x$-axis at 1, can be written as
# 
# $$x = 1$$

# In[24]:


field.plane(x=1)


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

# In[25]:


field.plane("x")


# ## Cascading the operations
# 
# Let us say we want to compute the average 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.

# In[26]:


field.plane(y=10).x.average


# This gives the same result as for instance

# In[27]:


field.x.plane(y=10).average


# ## Complex fields
# 
# `discretisedfield` supports complex-valued fields.

# In[28]:


cfield = df.Field(mesh, dim=3, value=(1 + 1.5j, 2, 3j))


# We can extract `real` and `imaginary` part.

# In[29]:


cfield.real((0, 0, 0))


# In[30]:


cfield.imag((0, 0, 0))


# Similarly we get `real` and `imaginary` parts of individual components.

# In[31]:


cfield.x.real((0, 0, 0))


# In[32]:


cfield.x.imag((0, 0, 0))


# Complex conjugate.

# In[33]:


cfield.conjugate((0, 0, 0))


# Phase in the complex plane.

# In[34]:


cfield.phase((0, 0, 0))


# ## Algebra operations
# 
# Let us define two fields:

# In[35]:


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, dim=3, value=(1, 1, 0))
f2 = df.Field(mesh, dim=3, value=(2, 1, 3))


# In[36]:


f1.average


# In[37]:


f2.average


# ### `+` operation

# In[38]:


f1 + f2


# In[39]:


(f1 + f2).average


# ### `-` operation

# In[40]:


f1 - f2


# In[41]:


(f1 - f2).average


# ### `*` operation
# 
# Basic multiplication is not defined between vector fields. In that case, we perform either dot or cross product, which we are going to discuss later.

# In[42]:


f1 * f2  # both are vector fields


# Scalar with vector field:

# In[43]:


f1.x * f2


# Scalar with vector field:

# In[44]:


f1.x * f2.y


# ### `/` operation

# In[45]:


f1 / f2  # both are vector fields


# Dividing vector field by a scalar field:

# In[46]:


f1 / f2.x


# Scalar field divided by vector field is not allowed:

# In[47]:


f2.x / f1


# ### `**` operator
# 
# This operator is allowed only on scalar fields:

# In[48]:


f1**2


# In[49]:


f1.x**2


# ### Compund operations

# In[50]:


f1 += f2


# In[51]:


f1 -= f2


# In[52]:


f1 *= f2.x


# In[53]:


f2 /= f2.y


# ## Vector products
# 
# As the title says, these products are applied between vector fields only.
# 
# ### Dot product
# 
# Dot product is implemented through `@` operator:

# In[54]:


f1 @ f2


# ### Cross product
# 
# Cross product between vector fields is performed using `&` operator:

# In[55]:


f1 & f2


# ## Vector calculus
# 
# ### Directional derivative $\left(\frac{\partial}{\partial x_{i}}f\right)$
# 
# Defined on both scalar and vector fields:

# In[56]:


f1.derivative("x")


# ### Gradient $(\nabla f)$
# 
# Defined on scalar fields:

# In[57]:


f1.x.grad


# ### Divergence $(\nabla \cdot f)$
# 
# Defined on vector fields:

# In[58]:


f1.div


# ### Curl  $(\nabla \times f)$
# 
# Defined on vector fields:

# In[59]:


f1.curl


# ### Laplace operator $(\nabla^{2} f)$
# 
# Defined on both vector and scalar fields:

# In[60]:


f1.laplace


# In[61]:


f1.x.laplace


# ## Integrals
# 
# In the most recent version of `discretisedfield`, computing integrals is generalised to accomodate the calculation of different types of integrals. Instead of giving the "theory" behind how it was implemented, we are going to show several examples which hopefully are going to give you an idea how you can compute different integrals.
# 
# Let us first create a field:

# In[62]:


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, dim=3, value=(-3, 0, 4), norm=1e6)


# ### Volume integral
# 
# $$\int_{V}\mathbf{f}(\mathbf{r})\text{d}V$$

# In[63]:


df.integral(f * df.dV)


# Since $\text{d}V = \text{d}x\text{d}y\text{d}z$, we can compute the integral as:

# In[64]:


df.integral(f * df.dx * df.dy * df.dz)


# $$\int_{V}f_{x}(\mathbf{r})\text{d}V$$

# In[65]:


df.integral(f.x * df.dV)


# ### 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}$$
# 
# Like all plane-related operations, the field must be sliced.

# In[66]:


df.integral(f.plane("z") @ df.dS)


# Similarly, we can write $\text{d}S = \text{d}x\text{d}y$ when we cut $z$-plane or we can use $|\text{d}\mathbf{S}|$.
# 
# $$\int_{S}f_{x}(\mathbf{r}) \text{d}x\text{d}y$$

# In[67]:


df.integral(f.x.plane("z") * df.dx * df.dy)


# $$\int_{S}f_{x}(\mathbf{r}) |\text{d}\mathbf{S}|$$

# In[68]:


df.integral(f.x.plane("z") * abs(df.dS))


# ### Line integrals
# 
# $$\int_{0}^{x_\text{max}}\mathbf{f}(\mathbf{r}) \text{d}x$$

# In[69]:


df.integral(f * df.dx, direction="x")


# $$\int_{0}^{y_\text{max}}f_{x}(\mathbf{r}) \text{d}y$$

# In[70]:


df.integral(f.x * df.dy, direction="y")


# ### Improper integrals

# $$\int_{0}^{y}\mathbf{f}(\mathbf{r}) \text{d}y'$$

# In[71]:


df.integral(f.x * df.dy, direction="y", improper=True)


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

# In[72]:


df.integral(f.div * df.dV)


# ## Operation pipelines

# In[73]:


df.integral(f1.x.grad.div.grad.curl.y.grad * df.dV)


# ## Example
# 
# 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$$

# In[74]:


import math

m = field.orientation.plane("z")
S = df.integral(m @ (m.derivative("x") & m.derivative("y")) * df.dx * df.dy) / (
    4 * math.pi
)
S


# Or using Ubermag function:

# In[75]:


import discretisedfield.tools as dft

dft.topological_charge(m)


# Using Berg-Luescher method

# In[76]:


dft.topological_charge(m, method="berg-luescher")


# ## Applying `numpys` universal functions
# All 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](https://numpy.org/doc/stable/reference/ufuncs.html#available-ufuncs).

# In[77]:


import numpy as np


# In[78]:


f1 = df.Field(mesh, dim=1, value=1)
f2 = df.Field(mesh, dim=1, value=np.pi)
f3 = df.Field(mesh, dim=1, value=2)


# In[79]:


np.sin(f1)


# In[80]:


np.sin(f2)((0, 0, 0))


# In[81]:


np.sum((f1, f2, f3))((0, 0, 0))


# In[82]:


np.exp(f1)((0, 0, 0))


# In[83]:


np.power(f3, 2)((0, 0, 0))


# ## Other
# 
# Full description of all existing functionality can be found in the [API Reference](https://ubermag.github.io/api/_autosummary/discretisedfield.Field.html).