Field distortion by a perfect mu-metal shield

import pyvista as pv
from bfieldtools._example_plots import mesh_actor
from bfieldtools._example_plots import triangle_mesh, vector_glyphs

pv.set_plot_theme("document")


import numpy as np
import trimesh

from bfieldtools.mesh_conductor import MeshConductor, StreamFunction
from bfieldtools.mesh_calculus import gradient
from bfieldtools.utils import load_example_mesh, combine_meshes

# This doesn't matter, the problem is scale-invariant
scaling_factor = 1


planemesh = load_example_mesh("10x10_plane_hires")
planemesh.apply_scale(scaling_factor)

# Specify coil plane geometry
center_offset = np.array([9, 0, 0]) * scaling_factor
standoff = np.array([0, 4, 0]) * scaling_factor

# Create coil plane pairs
coil_plus = trimesh.Trimesh(
    planemesh.vertices + center_offset + standoff, planemesh.faces, process=False
)

coil_minus = trimesh.Trimesh(
    planemesh.vertices + center_offset - standoff, planemesh.faces, process=False
)

joined_planes = combine_meshes((coil_plus, coil_minus))
planemesh = joined_planes

# Create mesh class object
coil = MeshConductor(mesh_obj=joined_planes, fix_normals=True)

# Separate object for shield geometry
shieldmesh = load_example_mesh("closed_cylinder_remeshed")
shieldmesh.apply_scale(15)

shield = MeshConductor(mesh_obj=shieldmesh, process=True, fix_normals=True)

N = 80
points = np.zeros((N * N, 3))
w = 12
X, Y = np.meshgrid(np.linspace(-w, w, N), np.linspace(-w, w, N), indexing="ij")
X += planemesh.vertices.mean(axis=0)[0]
Y += planemesh.vertices.mean(axis=0)[1]
points[:, 0] = X.flatten()
points[:, 1] = Y.flatten()
points[:, 2] += planemesh.vertices.mean(axis=0)[2]


def plot_plane():
    plotter.add_mesh(
        triangle_mesh(
            [X.min(), X.max(), X.max(), X.min()],
            [Y.min(), Y.min(), Y.max(), Y.max()],
            np.zeros(4),
            np.array([[0, 1, 2], [2, 3, 0]]),
        ),
        color=(0.1, 0.1, 0.1),
        opacity=0.6,
        show_scalar_bar=False,
    )

plotter = fig = pv.Plotter()
s0 = plotter.add_mesh(
    triangle_mesh(*shieldmesh.vertices.T, shieldmesh.faces),
    color=(0.5, 0.5, 0.5),
    opacity=0.3,
    show_scalar_bar=False,
)
mesh_actor(s0).prop.backface_culling = False
mesh_actor(s0).prop.ambient = 0.5

I_prim = np.load("../publication_software/Shielded coil/biplanar_streamfunction.npy")
sprim = StreamFunction(coil.vert2inner @ I_prim, coil)
m = max(abs(sprim))
plotter = s1 = sprim.plot(False, 16, vmin=-m, vmax=m)
plotter = s2 = sprim.plot(True, 20)
mesh_actor(s2).mapper.scalar_visibility = False
mesh_actor(s2).prop.line_width = 1.2
#                          colormap='RdBu')
mesh_actor(s2).prop.render_lines_as_tubes = True
plot_plane()
scene = plotter
scene.camera.position = [39.154871143623325, -40.509425675368334, 26.56155776567048]
scene.camera.focal_point = [
    4.239945673839333,
    1.041549923485209,
    -0.0005302515738243585,
]
scene.camera.view_angle = 30.0
scene.camera.up = [-0.28276020498745635, 0.33658483701858727, 0.898196701154387]
scene.camera.clipping_range = [16.139073445910277, 116.31572537292347]
scene.render()

%% Calculate primary potential matrix

Compute slightly inside

d = 1e-3
P_prim = coil.U_coupling(shieldmesh.vertices - d * shieldmesh.vertex_normals)

plotter = pv.Plotter()
s = plotter.add_mesh(
    triangle_mesh(*shieldmesh.vertices.T, shieldmesh.faces),
    scalars=P_prim @ sprim,
    opacity=1.0,
    show_scalar_bar=False,
)
mesh_actor(s).mapper.dataset = mesh_actor(s).mapper.dataset.contour(10)

%% Calculate linear collocation BEM matrix

P_shield = shield.U_coupling(shieldmesh.vertices - d * shieldmesh.vertex_normals)

%% Solve equivalent stream function for the perfect linear mu-metal layer

I_shield = np.linalg.solve(-P_shield, P_prim @ sprim)
# I_shield = P_prim @ I_prim
s_shield = StreamFunction(I_shield, shield)
g = gradient(s_shield, shieldmesh, rotated=True)

plotter = fig = pv.Plotter()
s0 = plotter.add_mesh(
    triangle_mesh(*shieldmesh.vertices.T, shieldmesh.faces),
    color=(0.5, 0.5, 0.5),
    opacity=0.3,
    show_scalar_bar=False,
)
mesh_actor(s0).prop.backface_culling = False
plotter = s1 = s_shield.plot(False, 256)
mesh_actor(s1).prop.backface_culling = False
# s2 = s_shield.plot(True, 10)
plotter.add_mesh(
    vector_glyphs(
        *shieldmesh.triangles_center.T,
        *g,
        scale_factor=8e-07,
        mode="arrow",
        scale_mode="vector"
    ),
    color=(1, 1, 1),
    show_scalar_bar=False,
)

scene = plotter
plot_plane()
scene.camera.position = [39.154871143623325, -40.509425675368334, 26.56155776567048]
scene.camera.focal_point = [
    4.239945673839333,
    1.041549923485209,
    -0.0005302515738243585,
]
scene.camera.view_angle = 30.0
scene.camera.up = [-0.28276020498745635, 0.33658483701858727, 0.898196701154387]
scene.camera.clipping_range = [16.139073445910277, 116.31572537292347]
scene.render()

U2_shield = shield.U_coupling(points)
U2_prim = coil.U_coupling(points)

from bfieldtools.contour import scalar_contour

cc1 = scalar_contour(planemesh, planemesh.vertices[:, 2], contours=[-0.1])
cc1 = np.vstack(cc1)
cc1a = cc1[: cc1.shape[0] // 2]
cc1b = cc1[cc1.shape[0] // 2 :]
cc2 = scalar_contour(shieldmesh, shieldmesh.vertices[:, 2], contours=[-0.1])
cc2 = np.vstack(cc2)
cc2a = cc1[: cc2.shape[0] // 2]
cc2b = cc1[cc2.shape[0] // 2 :]

import matplotlib.pyplot as plt
import matplotlib.colors as colors


def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=256):
    new_cmap = colors.LinearSegmentedColormap.from_list(
        "trunc({n},{a:.2f},{b:.2f})".format(n=cmap.name, a=minval, b=maxval),
        cmap(np.linspace(minval, maxval, n)),
    )
    return new_cmap


cmap = plt.get_cmap("RdBu")
# cmap.set_over((0.95,0.95,0.95))
# cmap.set_under((0.95,0.95,0.95))

u0 = abs(
    np.sum(U2_shield, axis=1).reshape(N, N)
)  # Solid angle of the shield, zero outside
u0 /= u0.max()
u0[u0 < 1e-6] = 0
u1 = (U2_prim @ sprim).reshape(N, N)
u2 = (U2_shield @ I_shield).reshape(N, N) * u0
u3 = (u1 + u2) * u0

vmax = np.max(abs(u3)) * 0.99
levels = np.linspace(-vmax, vmax, 120)
levels = np.hstack((-np.max(abs(u3)), levels, np.max(abs(u3))))
p = plt.contourf(
    X,
    Y,
    u1,
    levels=levels,
    cmap=cmap,
    vmin=-vmax,
    vmax=vmax,
    norm=colors.SymLogNorm(linthresh=0.2 * vmax, linscale=0.8, vmin=-vmax, vmax=vmax),
)
plt.plot(cc1a[:, 0], cc1a[:, 1], linewidth=3, color="gray")
plt.plot(cc1b[:, 0], cc1b[:, 1], linewidth=3, color="gray")
plt.axis("image")
plt.axis("off")
xlims = p.axes.get_xlim()
plt.figure()
p = plt.contourf(
    X,
    Y,
    u2,
    levels=levels,
    cmap=cmap,
    vmin=-vmax,
    vmax=vmax,
    norm=colors.SymLogNorm(linthresh=0.2 * vmax, linscale=0.8, vmin=-vmax, vmax=vmax),
)
plt.plot(cc1a[:, 0], cc1a[:, 1], linewidth=3, color="gray")
plt.plot(cc1b[:, 0], cc1b[:, 1], linewidth=3, color="gray")
plt.plot(cc2[:, 0], cc2[:, 1], linewidth=3, color="gray")
plt.axis("image")
plt.axis("off")
p.axes.set_xlim(xlims)
plt.figure()
p = plt.contourf(
    X,
    Y,
    u3,
    levels=levels,
    cmap=cmap,
    vmin=-vmax,
    vmax=vmax,
    norm=colors.SymLogNorm(linthresh=0.2 * vmax, linscale=0.8, vmin=-vmax, vmax=vmax),
)
plt.plot(cc1a[:, 0], cc1a[:, 1], linewidth=3, color="gray")
plt.plot(cc1b[:, 0], cc1b[:, 1], linewidth=3, color="gray")
plt.plot(cc2[:, 0], cc2[:, 1], linewidth=3, color="gray")
plt.axis("image")
plt.axis("off")
p.axes.set_xlim(xlims)

# Display the assembled scene (Sphinx-Gallery captures it during docs builds).
if "plotter" in globals() and not plotter._closed:
    plotter.show()