Bug in Dipole and Wire sources

[prisae]

Horrible bug in Dipole and Wire sources

import emg3d
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

h = np.ones(24)*100
grid = emg3d.TensorMesh([h, h, h], 'CCC')

src_coo = np.array([[-1000, +1000, 0], [+1000, -1000, 0]])
source = emg3d.TxElectricDipole(coordinates=src_coo)
frequency = 10
sfield = emg3d.get_source_field(grid, source, frequency)

sfield.grid.plot_3d_slicer(sfield.fx.ravel('F'), view='abs', v_type='Ex', pcolor_opts={'norm': LogNorm()})

Output:

emg3d/fields.py:923: UserWarning: emg3d: Normalizing Source: 58.0000000000.
emg3d/fields.py:923: UserWarning: emg3d: Normalizing Source: 58.0000000000.
emg3d/fields.py:923: UserWarning: emg3d: Normalizing Source: 58.0000000000.

[prisae]

The problem lies in the function emg3d.fields._dipole_vector:
https://github.com/emsig/emg3d/blob/main/emg3d/fields.py#L792C5-L931

It works if:

• dipole (or wire segments) are aligned along principal coordinate axis, or
• dipole (or wire segments) point coordinates are all increasing/decreasing, or
• in all cases, if the points are never further away than one cell.

It does NOT work if:

• dipole (or wire segments) point coordinates are mixed increasing/decreasing and the points are further away than one cell.

NEEDS FIXING

[prisae]

Current Work-Around

Refine your source by adding enough points so that the wire points are never further away than one cell. If emg3d does not raise a Warning any longer, then it is fine.

Here is an example function that you can use to refine your source:

def refine(source):
    pts = source.points
    new_pts = np.zeros((pts.shape[0]*2-1, 3))
    for i in range(pts.shape[0]-1):
        new_pts[2*i, :] = pts[i, :]
        new_pts[2*i+1, :] = (pts[i, :] + pts[i+1, :]) / 2
    new_pts[-1, :] = pts[-1, :]
    return emg3d.TxElectricWire(new_pts)

The above example needs four times refining, after which there is no more warning and the source is correct:

new_src = source
for i in range(4):
    new_src = refine(new_src)
print(new_src.points)
sfield = emg3d.get_source_field(grid, new_src, frequency)
sfield.grid.plot_3d_slicer(sfield.fx.ravel('F'), view='abs', v_type='Ex', pcolor_opts={'norm': LogNorm()})

The output from the printing:

[[-1000.  1000.     0.]
 [ -875.   875.     0.]
 [ -750.   750.     0.]
 [ -625.   625.     0.]
 [ -500.   500.     0.]
 [ -375.   375.     0.]
 [ -250.   250.     0.]
 [ -125.   125.     0.]
 [    0.     0.     0.]
 [  125.  -125.     0.]
 [  250.  -250.     0.]
 [  375.  -375.     0.]
 [  500.  -500.     0.]
 [  625.  -625.     0.]
 [  750.  -750.     0.]
 [  875.  -875.     0.]
 [ 1000. -1000.     0.]]