Visualizing streamlines without streamplot about magpylib HOT 24 OPEN

This topic is dicussed in #677 before implementations are considered

from magpylib.

hi @quenot

So I have created a new branch called Docu_example_field_lines.

In that branch I have prepared the file for you:

docs/_pages/gallery/gallery_vis_2D_fieldlines.md

As you can see it is formatted in Markdown. We use some Sphinx flavor, but I will fix the formatting in the end, so dont worry about that.

What you have to do now is:

1. fork the magpylib repo to your account.
2. checkout the "Docu_example_field_lines" branch
3. edit the file that I have prepared for you (best practice: clone your repo to your local machine, edit there, push to your repo)
4. When you are done make a pull request from your branch (Docu_example_field_lines in your repo) to my branch (same name but in the Magpylib repo).

I will see your pull request, and review, edit, discuss, ask for changes and so on.

Problem: if you want to see the result you would have to install Sphinx and a bunch of other packages to compile the documentation locally, which I do not recommend because its a lot of effort. So best check out what the other gallery files look like, and what output this generates. If you are in the end unhappy with the result we can ofc edit it anytime :).

Please let me know if you have questions :).

from magpylib.

Thanks. I will try that as soon as I have some time for it.

from magpylib.

I made a first attempt:
https://github.com/quenot/magpylib/blob/Docu_example_field_lines/docs/_pages/gallery/gallery_vis_2D_fieldlines.md
It is a bit frustrating to not see the final result as preview does not seem to do everything well or completely.
Can you confirm that you see my commit and that it looks good? Thanks.
Georges.

from magpylib.

@OrtnerMichael Hello,
I have improved the cuboid example code in:
https://github.com/quenot/magpylib/blob/Docu_example_field_lines/docs/_pages/gallery/gallery_vis_2D_fieldlines.md
Please let me know if you see it and if it is OK. I still have to complete the accompanying text.
It is easy to change the resolution and see the effect. If I reduce from 500x400 to 250x200, there is very little change but the highest resolution gives something sharper and cleaner at the discontinuities on the vertical boundaries. It is quite fast anyway.
I am also designing a more complex example for my personal use, I might add it too.
Georges.

from magpylib.

hi @quenot

Sorry for the wait. Busy times.

To me this looks good. When you have the feeling we should give it a try to compile, just make a pull request from your repository into the Magpylib repository - doesn't have to be in any way a final version. As soon as your draft is in the magpylib repo, I will set up Read the docs to compile it, so we can have a look at it.

Only when its completely ready I will push into the magpylib main branch and it will appear in the main docs.

Here some feedback concerning the code:

1. I suggest that you focus on the example above, most users will only be interested in that, copy pasting code, and add the explanations below for the users who really want to know.

2. About the grid creation. You can use np.mgrid as shown here in example 2. That avoids initializing another array with np.array as you do in line 20. However, the bottleneck is anyways somewhere else, so the "slow" Python loop is not so relevant.

3. Your lines 35-39 seem a bit complicated to follow for choosing the start points. Why this choice?

4. Your example is quite general. However, by choosing e.g. different grid sizes in different directions (line 6) the code becomes more complicated. But the core of what you want to demonstrate is line 29 and then line 43. I would try to avoid making things fancy but focus on the central message. Most ppl will have more difficulties reading the code like this. And those ppl who want to make different window sizes, grid spacing and so on will anyways know how to do that.

from magpylib.

hi @OrtnerMichael

I agree with your comments. If I finalize it, I will make simpler and I will choose "hard" parameter values. I keep parameters for now as it is easier for investigating the problems and, unfortunately, I found quite serious ones.

While working on my more complex personal example, I discovered potentially fatal problems with the stability of the integration. First, the result depends on the order of the integration and on the location of its origin. This is due to the handling of boundary conditions and can be fixed in principle by changing a bit the use of cumsum on the borders and/or relative to the origin. The solution is just a bit more complex and I think that could provide a clean function for doing the integration optimally.

Second, and much worse, the method, as it currently is, does not work well with magnets not oriented along the axes. Here is what I get for instance by just adding:

cube.rotate_from_angax(40, 'z')

It is not difficult to understand what happens. This is an integration problem linked to the misalignment of the places in the boundary where there is a discontinuity in the B field and the grid used for the integration. The result is wrong and shifted if the boundary does not pass exactly in the middle of adjacent pixel centers. This is amplified when the intensity of the fields decreases after the boundary crossing. There actually also is such an effect for the case in which the magnet boundaries are parallel to the axes but it gets unnoticed. Indeed the proposed method is of very limited use if it cannot handle oriented magnets.

I am thinking of workaround likely to work which I did not try yet but it is not very simple. Basically, it consists in avoiding integration paths crossing the problematic boundaries by placing masks on them and make the integration by "water-shading" around them The process would eventually be completed by removing the mask at the end. I am not sure that it would be worth it though and maybe my original idea of modifying the streamplot code so that it handles closed streamlines is easier (certainly for the final user) and would lead to more accurate streamlines.

Using streamplot multiple times with masks is cumbersome, especially if there are several oriented magnets, but it may currently be the most practical solution for getting complete streamlines with given starting points.

from magpylib.

This is amplified when the intensity of the fields decreases after the boundary crossing

Hi @quenot,

just throwing this in, but would it help to normalize the field beforehand so that the amplitude is 1 everywhere?

from magpylib.

Hi @Alexboiboi
Thanks for the suggestion but that would very likely break the div B = 0 condition which is essential for the method. Also, normalizing the amplitude would not prevent the discontinuity in the direction and therefore in the components.

from magpylib.

hey @quenot

oje, this looks bad. This is what I initially expected from the way it was implemented (same integration start point for everyone, integration crossing charge surfaces).

The problem is that any integration-path based workaround cannot rely anymore on the cumsum trick :(, unless you find some complex criterion when to cumsum from the left and when from the right, but that will fail when magnet assemblies become more complex (multiple magents).

Maybe there is a way to avoid the numerical instability? Why do you think the integration problem is related to the location of the boundary in the cells?

from magpylib.

Hi @quenot , @OrtnerMichael ,

I just checked how good it works with pyvista just for comparison:

import magpylib as magpy
import numpy as np
import plotly.graph_objects as go
import pyvista as pv


def get_projected_cuboid(cuboid):
    # returns only last position/orientation
    p = cuboid._position[-1] / 2
    d = cuboid.dimension
    line = [
        [0, 0, 0],
        [1, 0, 0],
        [1, 1, 0],
        [0, 1, 0],
        [0, 0, 0],
        [0, 0, 1],
        [1, 0, 1],
        [1, 0, 0],
        [1, 0, 1],
        [1, 1, 1],
        [1, 1, 0],
        [1, 1, 1],
        [0, 1, 1],
        [0, 1, 0],
        [0, 1, 1],
        [0, 0, 1],
    ]
    line = (np.array(line) - 0.5) * d + p
    return cube._orientation[-1].apply(line)


def get_concatenated_lines(streamlines):
    # get the lines and points arrays
    lines = streamlines.lines
    points = streamlines.points

    # iterate over the lines array and extract the point coordinates
    line_points = []
    i = 0
    while i < lines.shape[0]:
        l = lines[i]
        inds = lines[i + 1 : i + 1 + l]
        i += l + 1
        line_points.extend([[None, None, None], *points[inds]])
    return np.array(line_points, dtype=float)


# Compute the B-field of a cuboid magnet on the grid
nx, ny, nz = 500, 200, 2  # grid size
wx, wy, wz = 10, 8, 0.01  # window size (scale 1 mm ~ 1 inch)
mx, my, mz = 2.4, 1.8, 10  # magnet size
nisos = 31
cube = magpy.magnet.Cuboid(magnetization=(0, -1000, 0), dimension=(mx, my, mz))
cube.rotate_from_angax(40, "z")

# Create a 3D grid with Pyvista
grid = pv.ImageData(
    dimensions=(nx, ny, nz),
    spacing=(wx / (nx - 1), wy / (ny - 1), wz / (nz - 1)),
    origin=(-wx / 2, -wy / 2, -wz / 2),
)

# compute B-field and add as data to grid
grid["B"] = magpy.getB(cube, grid.points)

# compute field lines
yl = 0
seed = pv.Line((-mx / 2, yl, 0), (mx / 2, yl, 0), nisos)
#seed = pv.Plane(i_size=wx, j_size=wy)
strl = grid.streamlines_from_source(
    seed,
    vectors="B",
    max_time=10,
    initial_step_length=0.01,
    integration_direction="both",
)

fig = go.Figure()
fig.update_layout(height=600, yaxis_scaleanchor="x")
x, y, z = get_concatenated_lines(strl).T
fig.add_scattergl(x=x, y=y, mode="lines", name="streamlines")
x, y, z = get_projected_cuboid(cube).T
fig.add_scattergl(x=x, y=y, mode="lines", name="cuboid")
x, y, z = seed.points.T
fig.add_scattergl(x=x, y=y, mode="markers", name="seeds", marker_size=4)
fig.add_shape(
    type="rect",
    x0=-wx / 2,
    y0=-wy / 2,
    x1=wx / 2,
    y1=wy / 2,
    showlegend=True,
    name="grid window",
)

from magpylib.

@Alexboiboi
Thanks, that is interesting. Unfortunately, it seems that pyvista does not handle the closed streamlines significantly better than matplotlib. I did not rotate the starting points with the magnet from the previous example because that was not necessary for showing the integration problem but if I do so with:

seed = pv.Line((-mx*np.cos(2*np.pi/9)/2, -mx*np.sin(2*np.pi/9)/2, 0), \
               (mx*np.cos(2*np.pi/9)/2, mx*np.sin(2*np.pi/9)/2, 0), nisos)

here is the result (spiraling problem again):

I understand that this is a difficult problem because of the sharp discontinuity in the field direction at some boundaries.

Pyvista does a bit better than matplotlib though as it seems to succeed in closing the "easy" streamlines. Maybe there is some threshold that can be adjusted so that it succeeds everywhere?

from magpylib.

@OrtnerMichael
Crossing "charged" surfaces is not a problem in itself for the method. The problem is with the empirical implementation of the integration and in particular with the fact that when a "charged" boundary is crossed, cumsum does something equivalent to attributing half of the values on the segment to each extremity while the proportions should actually be proportional to the length of the part of the segment which is on each side of the discontinuity, hence the "teeth" pattern (which, I think, confirm this analysis, especially considering where it happens). It does not seem easy to implement the integration in this way but it is likely that this would solve the problem.
Yes, whatever the fix, it is very likely to be much less simple and computationally efficient than a cumsum.

from magpylib.

It is possible to reduce spiraling by trimming the lines that form a loop and prevent them to further drift during integration. There is still the forward and backwards integration that mismatch a bit but maybe we could take only the forward one for closed loops.
The only thing that I found that also helps is to increase the discretization of the window.

import magpylib as magpy
import numpy as np
import plotly.graph_objects as go
import pyvista as pv


def trim_points(points):
    # Calculate vectors between consecutive points
    vec = points[1:] - points[:-1]
    # Calculate angles between consecutive vectors
    n1 = np.linalg.norm(vec[:-1], axis=1, keepdims=True)
    n2 = np.linalg.norm(vec[1:], axis=1, keepdims=True)
    dot = np.einsum("ij,ij->i", vec[:-1] / n1, vec[1:] / n2)
    dot[dot > 1], dot[dot < -1] = 1, -1
    angles = np.arccos(dot)
    # Cumulative sum of angles
    cumulative_angles = np.cumsum(angles)
    # Find the index where cumulative angle exceeds 360 degrees
    m1 = cumulative_angles >= (2 * np.pi)
    m2 = cumulative_angles <= (-2 * np.pi)
    trim_index = np.argmax(np.logical_or(m1, m2))
    if trim_index == 0:
        return points
    # Trim the array accordingly
    return points[: trim_index + 1]


def get_projected_cuboid(cuboid):
    # returns only last position/orientation
    p = cuboid._position[-1] / 2
    d = cuboid.dimension
    line = [
        [0, 0, 0],
        [1, 0, 0],
        [1, 1, 0],
        [0, 1, 0],
        [0, 0, 0],
        [0, 0, 1],
        [1, 0, 1],
        [1, 0, 0],
        [1, 0, 1],
        [1, 1, 1],
        [1, 1, 0],
        [1, 1, 1],
        [0, 1, 1],
        [0, 1, 0],
        [0, 1, 1],
        [0, 0, 1],
    ]
    line = (np.array(line) - 0.5) * d + p
    return cube._orientation[-1].apply(line)


def get_lines(streamlines, concatenate=False):
    # get the lines and points arrays
    lines = streamlines.lines
    points = streamlines.points

    # iterate over the lines array and extract the point coordinates
    line_points = []
    i = 0
    while i < lines.shape[0]:
        l = lines[i]
        inds = lines[i + 1 : i + 1 + l]
        i += l + 1
        pts = points[inds]
        pts = trim_points(pts)
        line_points.append(pts)
    if concatenate:
        return [np.concatenate([[[np.nan] * 3, *l] for l in line_points])]
    return line_points


# Compute the B-field of a cuboid magnet on the grid
nx, ny, nz = 500, 200, 2  # grid size
wx, wy, wz = 10, 8, 0.01  # window size (scale 1 mm ~ 1 inch)
mx, my, mz = 2.4, 1.8, 10  # magnet size
nisos = 31
cube = magpy.magnet.Cuboid(magnetization=(0, -1000, 0), dimension=(mx, my, mz))
cube.rotate_from_angax(40, "z")

# Create a 3D grid with Pyvista
grid = pv.ImageData(
    dimensions=(nx, ny, nz),
    spacing=(wx / (nx - 1), wy / (ny - 1), wz / (nz - 1)),
    origin=(-wx / 2, -wy / 2, -wz / 2),
)

# compute B-field and add as data to grid
grid["B"] = magpy.getB(cube, grid.points)

# compute field lines
yl = 0
line = [(-mx / 2, yl, 0), (mx / 2, yl, 0)]
line = cube._orientation.apply(line)
seed = pv.Line(*line, nisos)
# seed = pv.Plane(i_size=wx, j_size=wy)
strl = grid.streamlines_from_source(
    seed,
    vectors="B",
    max_time=10,
    initial_step_length=0.01,
    integration_direction="both",
)

fig = go.Figure()
fig.update_layout(height=600, yaxis_scaleanchor="x")

lines = get_lines(strl, concatenate=False)
for i, line in enumerate(lines):
    x, y, _ = line.T
    fig.add_scattergl(x=x, y=y, mode="lines", name=f"streamlines_{i:02d}")
x, y, z = get_projected_cuboid(cube).T
fig.add_scattergl(x=x, y=y, mode="lines", name="cuboid")
x, y, z = seed.points.T
fig.add_scattergl(x=x, y=y, mode="markers", name="seeds", marker_size=4)
fig.add_shape(
    type="rect",
    x0=-wx / 2,
    y0=-wy / 2,
    x1=wx / 2,
    y1=wy / 2,
    showlegend=True,
    name="grid window",
)

from magpylib.

Hi @Alexboiboi @OrtnerMichael,
Thanks again. This is almost exactly the kind of fix / enhancement that I was thinking of for matplotlib streamplot (and asked for in the enhancement request #27261 on matplotlib github). In this case, I think that it would have required to edit the code of the streamplot function. I was thinking of it with two minor modifications: (i) use atan2 on the cross product and the dot product for better taking into account large and sharp variations in the direction of the B field as it happens on some boundaries, and (ii) stop the integration after a rotation of (plus or minus) 1*pi only on each direction after the starting point for possibly a better junction.
I also noticed that increasing the resolution led to a better stability and accuracy of the streamlines.

from magpylib.

And this is how it looks like if I remove the redundant backwards integration for closed lines

import magpylib as magpy
import numpy as np
import plotly.graph_objects as go
import pyvista as pv


def trim_points(points):
    # Calculate vectors between consecutive points
    vec = points[1:] - points[:-1]
    # Calculate angles between consecutive vectors
    n1 = np.linalg.norm(vec[:-1], axis=1, keepdims=True)
    n2 = np.linalg.norm(vec[1:], axis=1, keepdims=True)
    dot = np.einsum("ij,ij->i", vec[:-1] / n1, vec[1:] / n2)
    dot[dot > 1], dot[dot < -1] = 1, -1
    angles = np.arccos(dot)
    # Cumulative sum of angles
    cumulative_angles = np.cumsum(angles)
    # Find the index where cumulative angle exceeds 360 degrees
    m1 = cumulative_angles > (2 * np.pi)
    m2 = cumulative_angles < (-2 * np.pi)
    trim_index = np.argmax(np.logical_or(m1, m2))
    if trim_index == 0:
        return points
    # Trim the array accordingly
    return points[: trim_index + 1]


def get_projected_cuboid(cuboid):
    # returns only last position/orientation
    p = cuboid._position[-1] / 2
    d = cuboid.dimension
    line = [
        [0, 0, 0],
        [1, 0, 0],
        [1, 1, 0],
        [0, 1, 0],
        [0, 0, 0],
        [0, 0, 1],
        [1, 0, 1],
        [1, 0, 0],
        [1, 0, 1],
        [1, 1, 1],
        [1, 1, 0],
        [1, 1, 1],
        [0, 1, 1],
        [0, 1, 0],
        [0, 1, 1],
        [0, 0, 1],
    ]
    line = (np.array(line) - 0.5) * d + p
    return cube._orientation[-1].apply(line)


def get_lines(streamlines, concatenate=False):
    # get the lines and points arrays
    lines = streamlines.lines
    points = streamlines.points

    # iterate over the lines array and extract the point coordinates
    lines_list = []
    i = 0
    line_ind = 0
    trimed=[]
    while i < lines.shape[0]:
        l = lines[i]
        inds = lines[i + 1 : i + 1 + l]
        i += l + 1
        pts = points[inds]
        new_pts = trim_points(pts)
        trimed.append(len(new_pts)!= len(pts))
        lines_list.append(new_pts)
        line_ind +=1
    # remove closed lines that iterate from both sides
    # first half is forward integration, second is backwards integration
    lines_list = [l for i,l in enumerate(lines_list) if not trimed[i] or i>=len(lines_list)/2]
    if concatenate:
        return [np.concatenate([[[np.nan] * 3, *l] for l in lines_list])]
    return lines_list


# Compute the B-field of a cuboid magnet on the grid
nx, ny, nz = 500, 200, 2  # grid size
wx, wy, wz = 10, 8, 0.01  # window size (scale 1 mm ~ 1 inch)
mx, my, mz = 2.4, 1.8, 10  # magnet size
nisos = 31
cube = magpy.magnet.Cuboid(magnetization=(0, -1000, 0), dimension=(mx, my, mz))
cube.rotate_from_angax(40, "z")

# Create a 3D grid with Pyvista
grid = pv.ImageData(
    dimensions=(nx, ny, nz),
    spacing=(wx / (nx - 1), wy / (ny - 1), wz / (nz - 1)),
    origin=(-wx / 2, -wy / 2, -wz / 2),
)

# compute B-field and add as data to grid
grid["B"] = magpy.getB(cube, grid.points)

# compute field lines
yl = 0
line = [(-mx / 2, yl, 0), (mx / 2, yl, 0)]
line = cube._orientation.apply(line)
seed = pv.Line(*line, nisos)
# seed = pv.Plane(i_size=wx, j_size=wy)
strl = grid.streamlines_from_source(
    seed,
    vectors="B",
    max_time=10,
    initial_step_length=0.01,
    integration_direction="both",
)

fig = go.Figure()
fig.update_layout(height=600, yaxis_scaleanchor="x")

lines = get_lines(strl, concatenate=False)
for i, line in enumerate(lines):
    x, y, _ = line.T
    fig.add_scattergl(x=x, y=y, mode="lines", name=f"streamlines_{i:02d}")
x, y, z = get_projected_cuboid(cube).T
fig.add_scattergl(x=x, y=y, mode="lines", name="cuboid")
x, y, z = seed.points.T
fig.add_scattergl(x=x, y=y, mode="markers", name="seeds", marker_size=4)
fig.add_shape(
    type="rect",
    x0=-wx / 2,
    y0=-wy / 2,
    x1=wx / 2,
    y1=wy / 2,
    showlegend=True,
    name="grid window",
)

from magpylib.

hmmm ... closed streamlines is a general problem. the way they are computed through a forward iterative algorithm implies that they will not be closed as a result of numerical issues - mainly the grid interpolation i guess.

The solution form @Alexboiboi is quite nice. Should we make this a docu example?

from magpylib.

@Alexboiboi @OrtnerMichael
Thanks Alexandre for the fix.

I am currently working on a further improvement for forcing a closing by slightly distorting the integrated streamlines. If these should be closed for theoretical reasons, the distortion would rather correct the integration drift rather than adding noise. Not really difficult but a bit laborious, I will tell you when it is completed.

About the angle calculation, arccos returns a result between 0 and pi and the accumulation is likely to be incorrect if the angle variation is not monotonous. I have modified the code with:

    # Calculate vectors between consecutive points
    vec = points[1:]-points[:-1]
    # Calculate angles between consecutive vectors
    (vx1, vy1, _), (vx2, vy2, _) = vec[:-1].T, vec[1:].T
    angles = np.arctan2(vx1*vy2-vx2*vy1, vx1*vx2+vy1*vy2)

from magpylib.

@OrtnerMichael @Alexboiboi
Hello,
Here it is (updated on November 14th):

import magpylib as magpy
import numpy as np
import pyvista as pv
import plotly.graph_objects as go

def get_projected_cuboid(cuboid):
    # returns only last position/orientation
    p = cuboid._position[-1] / 2
    d = cuboid.dimension
    line = [
        [0, 0, 0],
        [1, 0, 0],
        [1, 1, 0],
        [0, 1, 0],
        [0, 0, 0],
        [0, 0, 1],
        [1, 0, 1],
        [1, 0, 0],
        [1, 0, 1],
        [1, 1, 1],
        [1, 1, 0],
        [1, 1, 1],
        [0, 1, 1],
        [0, 1, 0],
        [0, 1, 1],
        [0, 0, 1],
    ]
    line = (np.array(line) - 0.5) * d + p
    return cube._orientation[-1].apply(line)

def get_pv_lines(streamlines):
    # get the lines and points arrays
    lines = streamlines.lines
    points = streamlines.points
    # iterate over the lines array and extract the point coordinates
    line_points = []
    i = 0
    while i < lines.shape[0]:
        l = lines[i]
        inds = lines[i+1:i+1+l]
        i += l+1
        pts = points[inds][:,:2]
        line_points.append(pts)
    return line_points

def clean_line(line):
    # Some integrators do produce lines with sequences of multiple identical points
    i, k = 0, line.shape[0]-1
    while i < k:
        if np.array_equal(line[i+1], line[i]):
            line = np.delete(line, i, axis=0)
            k -= 1
        else:
            i += 1
    return(line)

def get_cumulative_angles(line):
    # Calculate vectors between consecutive points
    vec = line[1:]-line[:-1]
    # Calculate angles between consecutive vectors
    (vx1, vy1), (vx2, vy2) = vec[:-1].T, vec[1:].T
    angles = np.arctan2(vx1*vy2-vx2*vy1, vx1*vx2+vy1*vy2)
    # Cumulative sum of angles
    return np.cumsum(angles)

def has_loop(line):
    cumulative_angles = get_cumulative_angles(line)
    # Loops are detected after a + or - 2*pi rotation, that heuristic may fail in some cases
    if abs(cumulative_angles[-1]) < 2*np.pi: return 0  # There is no loop
    if cumulative_angles[-1] > +2*np.pi: return +1 # Trigonometric loop
    if cumulative_angles[-1] < -2*np.pi: return -1 # Clockwise loop

def projection_lambda(line, point, i):
    p0, p1 = line[i:i+2]
    return np.sum((point-p0)*(p1-p0))/np.sum((p1-p0)*(p1-p0))

def get_quasi_loop(line):
    # Quasi-loop is line interrupted at the projection of the seed point after a + or - 2*pi rotation
    cumulative_angles = get_cumulative_angles(line)
    # Start after approximately a + or - 2*pi rotation, that heuristic may fail in some cases
    i = np.argmax(abs(cumulative_angles) > 2*np.pi)
    # Find the projection of the seed
    previous_move = 0  # Avoid infinite back and forth if the projection is on a corner
    while True:
        l = projection_lambda(line, line[0], i)
        if l < 0:
            if previous_move == +1 or i+2 == len(line):
                break  # Projection is on a corner or we would get out of the line
            else:
                i -= 1
                previous_move = -1
        elif l > 1:
            if previous_move == -1 or i == 0:
                break  # Projection is on a corner or we would get out of the line
            else:
                i += 1
                previous_move = +1
        else:
            break  # Projection of teh seed is within the i to i+1 segment, bounds included
    if i == 0 or i+2 >= len(line):
        return None  # We should not be at an end of the line
    l = np.clip(projection_lambda(line, line[0], i), 0, 1)
    quasi_loop = line[0:i+1]
    if l > 0:
        quasi_loop = np.concatenate([quasi_loop, [(1-l)*line[i]+ l*line[i+1]]], axis=0)
    return quasi_loop

def close_loop(line, max_iter=50, verbose=0):
    previous_drift, stop = np.inf, False
    for k in range(max_iter):
        # print(k, previous_drift)
        quasi_loop = get_quasi_loop(line)
        if quasi_loop is None:
            return line, np.inf
        # Quasi-loop normal at the end
        vectors = quasi_loop[1:]-quasi_loop[:-1]
        normal = np.array([vectors[-1,1], -vectors[-1,0]])
        normal /= np.linalg.norm(normal)
        # Quasi-loop drift
        drift = np.sum(normal*(quasi_loop[-1]-quasi_loop[0]))
        if previous_drift == np.inf:
            initial_drift = drift
        if abs(drift) < abs(previous_drift):
            previous_drift = drift
            previous_quasi_loop = quasi_loop
        else:
            break
        # Quasi-loop length
        lengths = np.linalg.norm(vectors, axis=1)
        length = np.sum(lengths)
        # Quasi-loop center
        centers = 0.5*(quasi_loop[:-1]+quasi_loop[1:])
        center = np.sum(centers*np.expand_dims(lengths, axis=1), axis=0)/length
        # From center unit vectors
        nr, nrpoint = line-center, quasi_loop[-1]-center
        nr, nrpoint = nr/np.linalg.norm(nr,axis=1, keepdims=True), nrpoint/np.linalg.norm(nrpoint)
        # Drift correction factor
        dcf = drift/length/np.sum(normal*nrpoint)
        # Correct drift
        padded_cumulative_lengths = np.pad(np.cumsum(np.linalg.norm(line[1:]-line[:-1], axis=1)), (1, 0))
        line = line-dcf*nr*np.expand_dims(padded_cumulative_lengths, axis=1)
    if verbose > 0:
        print(k+1, "iterations, initial drift =", initial_drift, "final_drift =", previous_drift)
    previous_quasi_loop[-1] = previous_quasi_loop[0]  # Actually close the loop after drift is minimized
    return previous_quasi_loop, initial_drift

def merge_or_close_pair(pair, max_iter=50, verbose=0):
    # Assumes this order from streamlines_from_source(), NOT checked
    fline, bline = pair 
    fline, bline = clean_line(fline), clean_line(bline)
    floop, bloop = has_loop(fline), has_loop(bline)
    if floop == 0 and bloop == 0: # No loop
        if verbose > 0:
            print("No loop")
        return np.concatenate([np.flip(bline, axis=0), fline[1:]])
    if floop != 0 and bloop == 0: # Loop on forward line only, infrequent, not tested
        if verbose > 0:
            print("Loop on forward line only")
        return close_loop(fline)[0]
    if floop == 0 and bloop != 0: # Loop on forward line only, infrequent, not tested
        if verbose > 0:
            print("Loop on backward line only")
        return np.flip(close_loop(bline)[0], axis=0)
    # Loop on both forward and backward loops, choose the one with the smallest initial drift
    if verbose > 0:
        print("Loop on both forward and backward lines")
    fline, fdrift = close_loop(fline, max_iter=max_iter, verbose=verbose)
    bline, bdrift = close_loop(bline, max_iter=max_iter, verbose=verbose)
    if fdrift == np.inf and bdrift == np.inf:
        if verbose > 0:
            print("failed to close loop")
        return fline [0:1]
    if abs(fdrift) < abs(bdrift):
        if verbose > 0:
            print("forward line has the smallest initial drift")
        return fline
    else:
        if verbose > 0:
            print("backward line has the smallest initial drift")
        return np.flip(bline, axis=0)

# Compute the B-field of a cuboid magnet on the grid
nx, ny, nz = 500, 400, 1  # grid size
wx, wy, wz = 10, 8, 0.01  # window size (scale 1 mm ~ 1 inch)
mx, my, mz = 2.4, 1.8, 10  # magnet size
nisos = 31
verbose = 1
cube = magpy.magnet.Cuboid(magnetization=(0, -1000, 0), dimension=(mx, my, mz))
cube.rotate_from_angax(40, "z")

# Create a 3D grid with Pyvista
# Grid points are located at the centers of the pixels of a (nx, ny) image
# filling the [-wx/2, wx/2] x [-wx/2, wx/2] area
dx, dy, dz = wx/nx, wy/ny, wz/nz
ox, oy, oz = -wx/2+dx/2, -wy/2+dy/2, 0
grid = pv.ImageData(dimensions=(nx, ny, nz), spacing=(dx, dy, dz), origin=(ox, oy, oz))

# compute B-field and add as data to grid
grid["B"] = magpy.getB(cube, grid.points)

# compute field lines
yl = 0
line = [(-mx/2, yl, 0), (mx/2, yl, 0)]
line = cube._orientation.apply(line)
seed = pv.Line(*line, nisos-1)
strl = grid.streamlines_from_source(
    seed,
    vectors="B",
    max_time=20,
    initial_step_length=0.01,
    integration_direction="both",
    integrator_type=45
)

lines = get_pv_lines(strl)
plines = []
for i, p in enumerate(seed.points):
    pair = [line for line in lines if np.array_equal(line[0], p[0:2])]
    if len(pair) == 2:
        if verbose > 0:
            print("Seed", i, end = " ")
        plines.append(merge_or_close_pair(pair, verbose=verbose))
    else:
        print("Seed", i, "p : forward or backward line is missing,", len(pair), "lines")

fig = go.Figure()
fig.update_layout(height=600, yaxis_scaleanchor="x")

for i, line in enumerate(plines):
    x, y = line.T
    fig.add_scattergl(x=x, y=y, mode="lines", name=f"streamlines_{i:02d}", line_width=1.5)
x, y, _ = get_projected_cuboid(cube).T
fig.add_scattergl(x=x, y=y, mode="lines", name="cuboid")
x, y, _ = seed.points.T
fig.add_scattergl(x=x, y=y, mode="markers", name="seeds", marker_size=3)
fig.add_shape(
    type="rect",
    x0=-wx / 2,
    y0=-wy / 2,
    x1=wx / 2,
    y1=wy / 2,
    # showlegend=True,
    name="grid window",
)

and the result (updated):

The code is quite long. I left the debugging information (verbose). This can be removed and also the parameters might be replaced by their values in the code for a simpler example. The resolution should be set large enough for correct results but this is not due to the line closing.

This method can still be improved for the case of streamlines theoretically including singularities (abrupt changes in directions as on some borders of the magnet). See complements in my next comment.

from magpylib.

@Alexboiboi @OrtnerMichael
As a complement, the method I have just given works but it is not the best I can think of. There actually are several possible ways to distort the lines so that they loop on themselves. All of them are quite equivalent when dealing with "smooth" streamlines but this is not the case with streamlines including singularities (abrupt changes in directions as on some borders of the magnet). I tried so far two of them. The first one consisting in moving the points according to the normal to the line failed because near the singularities for shrinking moves, the distorted streamline crosses itself. The second one consisting in moving the points along rays originating from the line's center of gravity worked and this is the one given above but I think that it is not optimal for several reasons, including the one that it changes the angle of the line at the seed point and at its opposite. Possibly worse and most visible, the correction is significantly exaggerated near the extremities of elongated streamlines.
Considering that the integration is very accurate (very small drift) for smooth streamlines and that therefore most of the drift is due to what happens near the singularities, a better method would concentrate the correction of the distortion in the neighborhood of the singularities. I am then thinking of an alternate method based on an interpolation between aligned versions of the forward and backward lines (this is not trivial but a priori doable) with variations of the interpolation parameter concentrated in the high curvature regions. This method should not need iterations. I will give it a try when I find some time.

After thinking a bit more about it, this is feasible and would almost certainly provide a more accurate result but this will not be a priority for me as the current version works already quite well and increasing the accuracy can alternatively be obtained by increasing the grid resolution.

from magpylib.

hi @quenot

Thanks for sharing this nice, and long piece of code.
As you say, there are multiple ways to make the streamlines closed loops, The one you described last seems also quite feasible. In an actual example I would probably show the method demonstrated by @Alexboiboi above, make it as simple as possible (e.g. by removing the get_projected_cuboid part, and maybe have a second part afterwards where it is shown how to solve the loop-closing problem.

@quenot I propose that you let us know when you feel that you have a "good" line-closing code that you would like to share. I understand that its not a priority.

from magpylib.

@OrtnerMichael @Alexboiboi
I have updated the code in my previous comment. It is now a bit more structured and cleaned, some bugs have been corrected, it is more robust and tested, and no-longer-necessary development printing has been removed. It is not shorter though and still based on deformations by dilation relatively to the loop center. The change in code structure has also been done in order to make easier the implementation of the second solution but I am still not sure that we need it. As an intermediate solution, a modification of the amplitude of the correction in the radial direction might have the same effect. Currently, is is just proportional the the curvilinear abscissa but that can be easily changed for "concentrating" the distortion in the high curvature parts. This is a second-order correction anyway, which is unlikely to be necessary if the grid resolution is high enough and this is necessary for a good initial integration, whose errors cannot be compensated for by the loop closing.

Regarding the practical aspect for the users, maybe the simplest solution is to put all the relevant procedures into a python module which could be imported, making the remainder of the example relatively simple. I would put there everything not related to pyvista so that the module could be used either with pyvista or with matplotlib. Even simpler (for the users) would be to have such a module integrated with matplotlib or pyvista.

from magpylib.

@OrtnerMichael @Alexboiboi

Hello,

I think that I now have a good line-closing code that I am ready to share. I have set up a github repository for it. You are welcome to integrate part of it in magpylib or to refer to it.

I don't think that I will develop something more elaborated. Actually, even better than the last method that I proposed would be to integrate the drift compensation directly inside the integrator (that was my initial idea) while linking the local importance of the compensation to the streamline curvature. I have ideas on how to do that but I don't think that there is a need for it and I don't have time for that either. The current version is now largely sufficient for my needs and it is likely to be so too for those of most users.

from magpylib.

hey @quenot

thanks a lot for sharing this, I will add everything to the discussion that users will find via searches, so that things can be used right away.
I'll check out the codes in detail and think about doku-implementation.
I have to confess that it is quite tempting - you are not the first user who would like to see field lines in 2D.

br michael

from magpylib.