Calculation of derivatives of a eigenmode

[HelgeGehring]

How would I apply an operator like

on the efield I calculate in https://github.com/HelgeGehring/waveguidemodes/blob/main/waveguidemodes/mode_solver.py ?

I've tried to calculate it like

def calculate_hfield(basis, xs, beta):
    xs = xs.astype(complex)

    @BilinearForm(dtype=np.complex64)
    def aform(e_t, e_z, v_t, v_z, w):
        return (-1j * beta * e_t[1] + e_z.grad[1]) * v_t[0] \
               + (1j * beta * e_t[0] - e_z.grad[0]) * v_t[1] \
               + e_t.curl * v_z

    a_operator = aform.assemble(basis)

    @BilinearForm(dtype=np.complex64)
    def bform(e_t, e_z, v_t, v_z, w):
        return dot(e_t, v_t) + e_z * v_z

    b_operator = bform.assemble(basis)

    return scipy.sparse.linalg.spsolve(b_operator, a_operator @ xs)

When I use this, the result looks like

The x-component of the result looks fine, the y-component seems wrong and should look more like

and I really don't know what is going on in z-direction.
Is this approach to apply the operator the right one? Or how would I apply it?

[kinnala]

What are you plotting there? The magnitude of the complex number?

[HelgeGehring]

If I use instead of the nedelec element a vector of linear elements, the calculation of e_z seems to work (except at the outer borders of the domain):

import numpy as np
import matplotlib.pyplot as plt

from skfem import *
from skfem.helpers import *

mesh = MeshTri.init_tensor(x=np.linspace(-1, 1, 50), y=np.linspace(-1, 1, 50))

basis = Basis(mesh, ElementVector(ElementTriP1()) * ElementTriP1())


# e_t(x,y)=[y,-x], e_z = 2
@LinearForm
def from_form(v_t, v_z, w):
    return w.x[1] * v_t[0] - w.x[0] * v_t[1] + 2 * v_z


@BilinearForm
def to_form(e_t, e_z, v_t, v_z, w):
    return inner(e_t, v_t) + inner(e_z, v_z)


x_0 = solve(to_form.assemble(basis), from_form.assemble(basis))
(et, et_basis), (ez, ez_basis) = basis.split(x_0)

et_basis.plot(et)
plt.show()
ez_basis.plot(ez, colorbar=True)
plt.show()


# Calculate

@BilinearForm
def from_form(e_t, e_z, v_t, v_z, w):
    return (-e_t.grad[0, 1] + e_t.grad[1, 0]) * v_z + e_z * (-v_t.grad[0, 1] + v_t.grad[1, 0])


@BilinearForm
def to_form(e_t, e_z, v_t, v_z, w):
    return inner(e_t, v_t) + inner(e_z, v_z)


x_1 = solve(*condense(to_form.assemble(basis), from_form.assemble(basis) @ x_0, D=basis.get_dofs()))
(et, et_basis), (ez, ez_basis) = basis.split(x_1)

et_basis.plot(et)
plt.show()
ez_basis.plot(ez, colorbar=True)
plt.show()

[kinnala]

This should be what you want to do in one line: ez = ez_basis.project(basis_n0.interpolate(x_0)[0].curl)

However, now I see that it gives the wrong result. I think we added ElementTriN0 quite hastily without adding too many tests so perhaps there is a bug lurking somewhere in its implementation.

[kinnala]

I might've found the bug. I'll add this as a test to ElementHcurl 2D implementation.

[kinnala]

I started fixing it in #958. Unfortunately I had to change the rotation direction of 2D curl to be consistent with what we had defined before in skfem.helpers.curl for curl of scalar field.

[HelgeGehring]

Amazing, thanks! Also nice that the curl will be consistent with the helper function!

Thanks also for the short form, nice to see that it can be written like that :)
Would there be a way to write the projection of the vector-field in the beginning also that short?

x_0 = basis.project(lambda x: [[x[1], -x[0]], 2])

would be what I've guessed (np.array doesn't work here as it's a nested list), but lists are not accepted there