Extracting fields as 2d arrays on rectangular grid #133

lokgar · 2024-03-05T15:18:03Z

lokgar
Mar 5, 2024

Hi there!

I've been playing with this nice package and need some assistance. My goal is to extract calculated modal field values and represent them as 2D arrays on a rectangular grid. This format is essential for the manipulation and analysis I need to do for my research.

I've delved into the package's source code and the scikit-fem documentation, but I'm still struggling to make it work. I grasp the concept of interpolating and re-mapping values from a triangular FEM mesh. However, the limited documentation makes it challenging to fully understand the interplay between methods related to the different bases, mesh, etc.

If anyone could offer some guidance with code snippets (which, I think, should be relatively simple for those familiar with skfem), I'd be incredibly thankful!

Answered by HelgeGehring

Mar 7, 2024

# Trying interpolator
basis = modes[0].basis
basis_fix = basis.with_element(ElementVector(ElementDG(ElementTriP1())))

(et, et_basis), (ez, ez_basis) = basis.split(modes[0].E)
(et_x, et_x_basis), (et_y, et_y_basis) = basis_fix.split(basis_fix.project(et_basis.interpolate(et)))

node_coordinates = mesh.p
xmin, xmax = node_coordinates[0, :].min() / 2, node_coordinates[0, :].max() / 2
ymin, ymax = node_coordinates[1, :].min() / 2, node_coordinates[1, :].max() / 2

grid_x, grid_y = np.mgrid[xmin:xmax:200j, ymin:ymax:200j]

coordinates = np.array([grid_x.flatten(), grid_y.flatten()])

grid_data_1 = np.array((et_x_basis.interpolator(et_x)(coordinates), et_y_basis.interpolator(et_y)(coordinates), …

View full answer

HelgeGehring · 2024-03-05T17:23:48Z

HelgeGehring
Mar 5, 2024
Maintainer

Yay, thanks!

grid_data = mode.basis.interpolator(mode.E)(array_with_coordinates)

or if used more often

interpolator = mode.basis.interpolator(mode.E)

grid_data_1 = interpolator(array_with_coordinates)
grid_data_2 = interpolator(array_with_other_coordinates)

That should do the job :)

But be aware that you take a slight hit in precision if you do this. It depends a lot on the exact problem if it's problematic.

What do you plan to do with the interpolated data? Maybe there's a way how we could get it done without a rectangular grid? :)

6 replies

lokgar Mar 6, 2024
Author

It doesn't work for me...

interpolator = modes[0].basis.interpolator(modes[0].E)

node_coordinates = mesh.p
xmin, xmax = node_coordinates[0, :].min()/2, node_coordinates[0, :].max()/2
ymin, ymax = node_coordinates[1, :].min()/2, node_coordinates[1, :].max()/2

grid_x, grid_y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]

array_with_coordinates = np.array([grid_x.flatten(), grid_y.flatten()])

grid_data_1 = interpolator(array_with_coordinates)

Output:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[89], [line 11](vscode-notebook-cell:?execution_count=89&line=11)
      [7](vscode-notebook-cell:?execution_count=89&line=7) grid_x, grid_y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
      [9](vscode-notebook-cell:?execution_count=89&line=9) array_with_coordinates = np.array([grid_x.flatten(), grid_y.flatten()])
---> [11](vscode-notebook-cell:?execution_count=89&line=11) grid_data_1 = interpolator(array_with_coordinates)

File [c:\Users\...\miniconda3\envs\femwell\Lib\site-packages\skfem\assembly\basis\cell_basis.py:222](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:222), in CellBasis.interpolator.<locals>.interpfun(x)
    [220](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:220)     shape = x.shape
    [221](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:221)     x = x.reshape(shape[0], -1)
--> [222](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:222) out = self.probes(x) @ y
    [223](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:223) # reshape output back to original shape
    [224](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:224) if shape is not None:

File [c:\Users\...\miniconda3\envs\femwell\Lib\site-packages\skfem\assembly\basis\cell_basis.py:187](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:187), in CellBasis.probes(self, x)
    [180](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:180) pts = self.mapping.invF(x[:, :, np.newaxis], tind=cells)
    [181](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:181) phis = np.array(
    [182](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:182)     [
    [183](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:183)         self.elem.gbasis(self.mapping, pts, k, tind=cells)[0]
    [184](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:184)         for k in range(self.Nbfun)
    [185](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:185)     ]
    [186](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:186) ).flatten()
--> [187](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:187) return coo_matrix(
    [188](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:188)     (
    [189](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:189)         phis,
...
    [151](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/scipy/sparse/_coo.py:151) if self.data.ndim != 1 or self.row.ndim != 1 or \
    [152](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/scipy/sparse/_coo.py:152)         self.col.ndim != 1:
    [153](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/scipy/sparse/_coo.py:153)     raise ValueError('row, column, and data arrays must be 1-D')

ValueError: row, column, and data array must all be the same length

Or when I pass single coordinate:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[90], [line 1](vscode-notebook-cell:?execution_count=90&line=1)
----> [1](vscode-notebook-cell:?execution_count=90&line=1) interpolator(array_with_coordinates[:, 0])

File [c:\Users\...\miniconda3\envs\femwell\Lib\site-packages\skfem\assembly\basis\cell_basis.py:222](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:222), in CellBasis.interpolator.<locals>.interpfun(x)
    [220](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:220)     shape = x.shape
    [221](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:221)     x = x.reshape(shape[0], -1)
--> [222](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:222) out = self.probes(x) @ y
    [223](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:223) # reshape output back to original shape
    [224](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:224) if shape is not None:

File [c:\Users\...\miniconda3\envs\femwell\Lib\site-packages\skfem\assembly\basis\cell_basis.py:179](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:179), in CellBasis.probes(self, x)
    [168](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:168) def probes(self, x: ndarray):
    [169](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:169)     """Return matrix which acts on a solution vector to find its values
    [170](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:170)     on points `x`.
    [171](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:171) 
   (...)
    [177](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:177) 
    [178](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:178)     """
--> [179](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:179)     cells = self.mesh.element_finder(mapping=self.mapping)(*x)
    [180](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:180)     pts = self.mapping.invF(x[:, :, np.newaxis], tind=cells)
    [181](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:181)     phis = np.array(
    [182](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:182)         [
    [183](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/skfem/assembly/basis/cell_basis.py:183)             self.elem.gbasis(self.mapping, pts, k, tind=cells)[0]
...
   [1373](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/numpy/core/einsumfunc.py:1373) # Check the kwargs to avoid a more cryptic error later, without having to
   [1374](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/numpy/core/einsumfunc.py:1374) # repeat default values here
   [1375](file:///C:/Users/.../miniconda3/envs/femwell/Lib/site-packages/numpy/core/einsumfunc.py:1375) valid_einsum_kwargs = ['dtype', 'order', 'casting']

ValueError: einstein sum subscripts string contains too many subscripts for operand 1

Thanks in advance if you can help :)

HelgeGehring Mar 7, 2024
Maintainer

Could you try array_with_coordinates.T?

And maybe post a short working script which I can just execute?

lokgar Mar 7, 2024
Author

I've tried different ways including transpose, but it didn't work.
Check the script below.

# Simulation setup from example
from collections import OrderedDict

import matplotlib.pyplot as plt
import numpy as np
import shapely
import shapely.affinity
from scipy.constants import epsilon_0, speed_of_light
from shapely.ops import clip_by_rect
from skfem import Basis, ElementTriP0, ElementVector, ElementDG, ElementTriP1
from skfem.io.meshio import from_meshio

from femwell.maxwell.waveguide import compute_modes, plot_mode
from femwell.mesh import mesh_from_OrderedDict
from femwell.visualization import plot_domains

wg_width = 0.6
wg_thickness = 0.22
core = shapely.geometry.box(-wg_width / 2, 0, +wg_width / 2, wg_thickness)
env = shapely.affinity.scale(core.buffer(5, resolution=8), xfact=0.5, yfact=0.5)

polygons = OrderedDict(
    core=core,
    box=clip_by_rect(env, -np.inf, -np.inf, np.inf, 0),
    clad=clip_by_rect(env, -np.inf, 0, np.inf, np.inf),
)

resolutions = dict(core={"resolution": 0.05, "distance": 180})

mesh = from_meshio(
    mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=20)
)
mesh.draw().show()

plot_domains(mesh)
plt.show()

basis0 = Basis(mesh, ElementTriP0())
epsilon = basis0.zeros()
for subdomain, n in {"core": 3.4757, "box": 1.444, "clad": 1.444}.items():
    epsilon[basis0.get_dofs(elements=subdomain)] = n**2
basis0.plot(epsilon, colorbar=True).show()

wavelength = 1.55

modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=3, order=2)

for mode in modes:
    print(f"Effective refractive index: {mode.n_eff:.4f}")
    print(f"TE fraction: {mode.te_fraction:.4f}")
    mode.show(mode.E.real, colorbar=True, direction="x")
    mode.show(mode.E.imag, colorbar=True, direction="x")

# Trying interpolator
interpolator = modes[0].basis.interpolator(modes[0].E)

node_coordinates = mesh.p
xmin, xmax = node_coordinates[0, :].min() / 2, node_coordinates[0, :].max() / 2
ymin, ymax = node_coordinates[1, :].min() / 2, node_coordinates[1, :].max() / 2

grid_x, grid_y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]

array_with_coordinates = np.array([grid_x.flatten(), grid_y.flatten()])

grid_data_1 = interpolator(array_with_coordinates)

HelgeGehring Mar 7, 2024
Maintainer

# Trying interpolator
basis = modes[0].basis
basis_fix = basis.with_element(ElementVector(ElementDG(ElementTriP1())))

(et, et_basis), (ez, ez_basis) = basis.split(modes[0].E)
(et_x, et_x_basis), (et_y, et_y_basis) = basis_fix.split(basis_fix.project(et_basis.interpolate(et)))

node_coordinates = mesh.p
xmin, xmax = node_coordinates[0, :].min() / 2, node_coordinates[0, :].max() / 2
ymin, ymax = node_coordinates[1, :].min() / 2, node_coordinates[1, :].max() / 2

grid_x, grid_y = np.mgrid[xmin:xmax:200j, ymin:ymax:200j]

coordinates = np.array([grid_x.flatten(), grid_y.flatten()])

grid_data_1 = np.array((et_x_basis.interpolator(et_x)(coordinates), et_y_basis.interpolator(et_y)(coordinates), ez_basis.interpolator(ez)(coordinates))).T

plt.imshow(grid_data_1[:,0].reshape(grid_x.shape).real)
plt.show()

does the job :)

It's a workaround for kinnala/scikit-fem#973 (forgot in my first answer that I somewhen had exactly that problem and didn't yet find time to see if I can fix that bug :/ )

Answer selected by lokgar

lokgar Mar 7, 2024
Author

Ah, I've even used this fix for plotting separate components of the field, but didn't try here. Thanks a lot!

HelgeGehring Mar 7, 2024
Maintainer

Happy to help :)
For plotting separate components, @lucasgrjn recently made PR #132, to have a function for that :)

lokgar · 2024-03-06T12:40:10Z

lokgar
Mar 6, 2024
Author

Hi @HelgeGehring, thank you for the quick response!

The example use case for me is performing mode decomposition for a field I got after simulating complex structures in FDTD. So, I need to re-mesh eigenmode fields to overlap with the resulting field from FDTD (or vice-versa, which is probably more complicated). It will decrease precision indeed, however we can make a denser mesh in the mode-solver.

Besides this use case, I've encountered issues accessing field components (Ex/Ey/Ez) and plotting them using non-built-in functions.

For instance: I do not fully comprehend those splits and re-basing in overlap calculation, but I would like to take the longitudinal component into account, as I experienced it to be crucial in some problems:

basis_j_fix = basis_j.with_element(ElementVector(ElementTriP1()))

(et, et_basis), (ez, ez_basis) = basis_j.split(E_j)
E_j = basis_j_fix.project(et_basis.interpolate(et), dtype=np.cfloat)
(et_x, et_x_basis), (et_y, et_y_basis) = basis_j_fix.split(E_j)

(et, et_basis), (ez, ez_basis) = basis_j.split(H_j)
H_j = basis_j_fix.project(et_basis.interpolate(et), dtype=np.cfloat)
(ht_x, ht_x_basis), (ht_y, ht_y_basis) = basis_j_fix.split(H_j)

def overlap(w):
        return cross(
            np.conj(w["E_i"][0]),
            np.array((ht_x_basis.interpolator(ht_x)(w.x), ht_y_basis.interpolator(ht_y)(w.x))),
        ) + cross(
            np.array((et_x_basis.interpolator(et_x)(w.x), et_y_basis.interpolator(et_y)(w.x))),
            np.conj(w["H_i"][0]),
        )

2 replies

lokgar Mar 6, 2024
Author

It is also relevant (for me) in full intensity calculation, but et_basis and ez_basis are different (why?) and I'm not sure how to align bases correctly to add longitudinal Ez/Hz components to analysis.

HelgeGehring Mar 6, 2024
Maintainer

Okay, sounds good :)

The reason for the bases of the transversal and normal component being different is the different requirements on the continuity, the normal component is continuous over the domain while only the normal component of the transversal components is continuous across interfaces. Thus we need to chose different elements for the two components
As long as we do the FEM, we don't see the difference too much as we just need to choose the right basis to interpolate/project.

For using things on a finite difference grid, you don't need to use a Form to assemble something. But it's important to either do things in FEM with interpolate and a Functional/... or to use the interpolator to get on a grid. (The example you posted interpolates between two different FEM grids. We implemented that for taking overlap integrals between waveguides of different width for example)

I think there was a problem to use the interpolator for Vector fields, therefore I projected first on scalar linear elements to interpolate from there.

coordinates = some2Dgrid

basis_j_fix = basis_j.with_element(ElementVector(ElementTriP1()))

(et, et_basis), (ez, ez_basis) = basis_j.split(E_j)
E_j = basis_j_fix.project(et_basis.interpolate(et), dtype=np.cfloat)
(et_x, et_x_basis), (et_y, et_y_basis) = basis_j_fix.split(E_j)

np.array((et_x_basis.interpolator(et_x)(coordinates), et_x_basis.interpolator(et_y)(coordinates)), ez_basis.interpolator(ez)(coordinates)))

will probably do the job to interpolate all three coordinates on a 2D grid 🤔 (not yet tested)

From there the overlap integral can be calculated by just summing up the products as usual in finite difference

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Extracting fields as 2d arrays on rectangular grid #133

{{title}}

Replies: 2 comments 8 replies

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

Select a reply

Extracting fields as 2d arrays on rectangular grid #133

lokgar Mar 5, 2024

Replies: 2 comments · 8 replies

HelgeGehring Mar 5, 2024 Maintainer

lokgar Mar 6, 2024 Author

HelgeGehring Mar 7, 2024 Maintainer

lokgar Mar 7, 2024 Author

HelgeGehring Mar 7, 2024 Maintainer

lokgar Mar 7, 2024 Author

HelgeGehring Mar 7, 2024 Maintainer

lokgar Mar 6, 2024 Author

lokgar Mar 6, 2024 Author

HelgeGehring Mar 6, 2024 Maintainer

lokgar
Mar 5, 2024

Replies: 2 comments 8 replies

HelgeGehring
Mar 5, 2024
Maintainer

lokgar Mar 6, 2024
Author

HelgeGehring Mar 7, 2024
Maintainer

lokgar Mar 7, 2024
Author

HelgeGehring Mar 7, 2024
Maintainer

lokgar Mar 7, 2024
Author

HelgeGehring Mar 7, 2024
Maintainer

lokgar
Mar 6, 2024
Author

lokgar Mar 6, 2024
Author

HelgeGehring Mar 6, 2024
Maintainer