How to calculate the magnetic force between several fixed magnets and one moving magnet?

[feldnerd]

I would like to calculate how much a moving magnet would move in the presence of two or more fixed magnets.

The first step is to calculate how much magnetic force is experienced by the moving magnet. The below code seems like a reasonable way to do this. Please let me know if this is incorrect and if there is a better way:

B_field = np.zeros(3)
B_field += moving_magnet.getB(fixed_magnet1.position)
B_field += moving_magnet.getB(fixed_magnet2.position)
magnetic_moment = moving_magnet.magnetization
magnetic_force = np.dot(B_field, magnetic_moment)

[OrtnerMichael]

hi @feldnerd

The force acting on a magnet is given by

np.dot(gradB*moment)

so first you would have to compute the gradient. The you have the problem that the field is not homogeneous inside the magnet, so you would have to split up the moving magnet into a lot of small parts, sum up over all the forces from the individual parts - otherwise you are computing the dipole-dipole force.

As it happens we are currently working on force codes. They are still far from being complete and very far from being integrated into Magpylib - this is happening towards the end of the year. However, if you are interested you could make use of our current codes.

br michael

[OrtnerMichael]

def getFTcube(sources, targets, eps=1e-5, anchor=None, squeeze=True):
    """
    Compute force and torque acting on a Cuboid magnet.

    Parameters
    ----------
    sources: source and collection objects or 1D list thereof
        Sources that generate the magnetic field. Can be a single source (or collection)
        or a 1D list of l sources and/or collection objects.

    targets: Cuboid object or 1D list of Cuboid objects
        Force and Torque acting on targets in the magnetic field generated by the sources
        will be computed. A target must have a valid mesh parameter.
    
    eps: float, default=1e-5
        The magnetic field gradient is computed using finite differences (FD). eps is
        the FD step size. A good value is 1e-5 * characteristic system size (magnet size,
        distance between sources and targets, ...).
    
    anchor: array_like, default=None
        The Force adds to the Torque via the anchor point. For a freely floating magnet
        this would be the barycenter. If `anchor=None`, this part of the Torque computation
        is ommitted.

    squeeze: bool, default=True
        The output of the computation has the shape (n,3) where n corresponds to the number
        of targets. By default this is reduced to (3,) when there is only one target.
    """
    anchor = check_input_anchor(anchor)
    targets = check_input_targets_Cuboid(targets)
    # MISSING: allow Collections

    # number of Cuboids
    tgt_number = len(targets)

    # number of instances of each Cuboid
    inst_numbers = np.array([np.prod(tgt.mesh) for tgt in targets])

    # total number of instances
    no_inst = np.sum(inst_numbers)

    # cumsum of number of instances (used for indexing)
    insti = np.r_[0, np.cumsum(inst_numbers)]

    # field computation positions (1xfor B, 6x for gradB)
    POSS = np.zeros((no_inst,7,3))

    # moment of each instance
    MOM = np.zeros((no_inst,3))

    # MISSING: eps should be defined relative to the sizes of the system
    eps_vec = np.array([(0,0,0), (eps,0,0), (-eps,0,0), (0,eps,0), (0,-eps,0), (0,0,eps), (0,0,-eps)])

    for i,tgt in enumerate(targets):
        tgt_vol = np.prod(tgt.dimension)
        inst_mom = tgt.magnetization * tgt_vol / inst_numbers[i]
        MOM[insti[i]:insti[i+1]] = inst_mom

        mesh = mesh_cuboid(tgt)
        for j,ev in enumerate(eps_vec):
            POSS[insti[i]:insti[i+1],j] = mesh + ev + tgt.position

    BB = magpy.getB(sources, POSS, sumup=True)
    gradB = (BB[:,1::2]-BB[:,2::2]) / (2*eps)
    gradB = np.swapaxes(gradB,0,1)

    Fs = np.sum((gradB*MOM),axis=2).T
    #Ts = np.zeros((no_inst,3))
    Ts = np.cross(BB[:,0], MOM)
    if anchor is not None:
        Ts -= np.cross(POSS[:,0]-anchor, Fs)

    T = np.array([np.sum(Ts[insti[i]:insti[i+1]],axis=0) for i in range(tgt_number)])
    F = np.array([np.sum(Fs[insti[i]:insti[i+1]],axis=0) for i in range(tgt_number)])

    if squeeze:
        F = np.squeeze(F)
        T = np.squeeze(T)

    return F, T

[OrtnerMichael]

def mesh_cuboid(object):
    """
    splits cuboid into given mesh
    returns grid positions relative to cuboid position
    """
    a,b,c = object.dimension
    n1,n2,n3 = object.mesh
    xs = np.linspace(-a/2, a/2, n1+1)
    ys = np.linspace(-b/2, b/2, n2+1)
    zs = np.linspace(-c/2, c/2, n3+1)

    dx = xs[1] - xs[0] if len(xs)>1 else a
    dy = ys[1] - ys[0] if len(ys)>1 else b
    dz = zs[1] - zs[0] if len(zs)>1 else c

    xs_cent = xs[:-1] + dx/2 if len(xs)>1 else xs + dx/2
    ys_cent = ys[:-1] + dy/2 if len(ys)>1 else ys + dy/2
    zs_cent = zs[:-1] + dz/2 if len(zs)>1 else zs + dz/2

    if isinstance(a, pint.Quantity):
        permutas = np.array(list(itertools.product(xs_cent.magnitude, ys_cent.magnitude, zs_cent.magnitude)))
        return permutas * xs_cent.units

    return np.array(list(itertools.product(xs_cent, ys_cent, zs_cent)))

[OrtnerMichael]

from magpylib._src.obj_classes.class_magnet_Cuboid import Cuboid
from magpylib._src.obj_classes.class_current_Polyline import Polyline


def check_input_anchor(anchor):
    """
    check and format anchor input
    """
    if anchor is None:
        return None
    if isinstance(anchor, (list, tuple)):
        anchor = np.array(anchor)
    if not isinstance(anchor, np.ndarray):
        raise ValueError("Anchor input must be list tuple or array.")
    if anchor.shape != (3,):
        raise ValueError("Anchor input must have shape (3,).")
    return anchor


def check_input_targets_Cuboid(targets):
    """ check and format targets input """
    if isinstance(targets, Cuboid):
        targets = [targets]
    if not isinstance(targets, list):
        raise ValueError("Bad targets input.")
    return targets

[OrtnerMichael]

let me know if this helps - these codes are working but they are far from being optimized in any way

[feldnerd]

Hi @OrtnerMichael,

Thank you for the swift and detailed response!

One thing that is confusing me, mesh_cuboid(object) takes in an object, but what kind of object. It seems like it should take a magnet, but magnets do not have a .mesh property.

For my use-case, an approximation would work. I do not need to know the exact force, just how the force on the moving magnet changes with different positions of the fixed magnets. Would then just computing the dipole-dipole force work, and if so is that easier to implement?

[OrtnerMichael]

Hi @feldnerd

here is an example how to use the code:

from mforce import getFTcube
import magpylib as magpy

# all inputs and output in SI units

cube1 = magpy.magnet.Cuboid(
    dimension=(1,1,1),
    polarization=(1,1,1),
)
cube2 = cube1.copy(position=(1,0,2))
cube2.mesh = (10,10,10)

F,T = getFTcube(cube1, [cube2, cube2.copy()], anchor=(0,0,0))
print(F)
print(T)

sure it would be faster to use a dipole approximation - which is the same as using cube2.mesh = (1,1,1) ;).

Whats your use case ? If you explain I can maybe tell you which approximation would be justified.

[feldnerd]

Thanks again! I am working on some physics simulations for my research. A vague answer... but I can say more once it is published :)

[OrtnerMichael]

:)

I just published the force codes as a separate package called magpylib-force.
The package development is on github.

pip install magpylib-force will do the trick

[feldnerd]

Amazing! Thanks for this!