--- orphan: true jupytext: text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.13.7 kernelspec: display_name: Python 3 language: python name: python3 --- (examples-force-intro)= # Force and Torque Basics This tutorial demonstrates Magpylib's force and torque computation capabilities through step-by-step examples. For comprehensive API details, see the [Documentation](docs-forcecomp). The underlying physics principles and computational methods are explained in the [Physics and Computation](guide-physics-force-computation) section. -------------------------------- ## Minimal Example Magpylib computes magnetic forces and torques via the `getFT(sources, targets)` function, which calculates the force F and torque T acting on `targets` in the magnetic field generated by `sources`. `getFT()` relies on a numerical scheme of discretizing the targets. **Critical Parameters:** - **`sources`**: Any Magpylib object that generates a magnetic field, or a 1D list thereof. - **`targets`**: Any Magpylib magnet, current, or `Dipole` object, or a 1D list thereof. - **`meshing`**: Each target object must have the `meshing` parameter explicitly set, which defines the discretization level. The following example computes the force and torque acting on a current loop positioned near a cubical permanent magnet: ```python import numpy as np import magpylib as magpy # Set numpy print format np.set_printoptions(formatter={"float": "{:.2e}".format}) # Source cube = magpy.magnet.Cuboid( dimension=(0.01, 0.01, 0.01), # 1cm cube polarization=(0.7, 0.7, 0.7), # units (T) ) # Target loop = magpy.current.Circle( diameter=0.02, # 2cm diameter current=1e3, # 1000 A position=(0, 0, 0.001), # 1mm above magnet center meshing=40, # discretization level ) # Compute force and torque F, T = magpy.getFT(cube, loop) print(f"Force: {F} N") # Force: [3.77e-01 3.77e-01 -7.53e-01] N print(f"Torque: {T} N*m") # Torque: [-3.18e-02 3.18e-02 7.45e-20] N*m ``` **Understanding the results:** - **Force**: The loop experiences attraction toward the magnet center (negative z-component) and symmetric forces in x and y directions due to the magnet's diagonal polarization. - **Torque**: The asymmetric positioning creates a rotational moment about the x and y axes. ```{warning} **Units matter!** Forces are computed in (N) and torques in (N*m). Unlike magnetic field calculations, [scaling invariance](guide-docs-io-scale-invariance) does not apply to force computations. Always use consistent SI units for all input parameters. ``` -------------------------------- ## The Pivot Point The total torque acting on an object consists of two components: the **intrinsic magnetic torque** and the **force moment** resulting from off-center forces. The force moment depends on the chosen pivot point, making torque calculations fundamentally dependent on this reference point, see [Physics and Computation](guide-physics-force-computation). ```{figure} ../../../_static/images/examples_force_intro01.png :width: 60% :align: center :alt: Sketch of force and torques acting on magnetic moment. Sketch of force and torque acting on a magnetic moment that is part of a large mass system. ``` - **Freely floating magnets**: The pivot is the center of mass (barycenter). If the magnet has homogeneous mass density the barycenter coincides with the centroid (geometric center). - **Mounted magnets**: The pivot is given by the barycenter of the larger mass system on which the magnet is mounted, or simply by a mechanical attachment point like an axle, bearing, or hinge. **Input for `pivot`:** - `pivot='centroid'`: this is the default value and selects the geometric center of the target as pivot. - `pivot=(x, y, z)`: any position can be specified as pivot point. - `pivot=None`: only the intrinsic magnetic torque contribution is computed. **Example: electric motor** The following example demonstrates how pivot point selection affects torque calculations in an electric motor configuration, where a permanent magnet mounted on a rotor interacts with current-carrying stator coils: ```python import numpy as np import magpylib as magpy # Set numpy print format np.set_printoptions(formatter={"float": "{:.2e}".format}) # Rotor magnet cube = magpy.magnet.Cuboid( dimension=(0.2, 0.2, 0.2), polarization=(0, 0, 1), meshing=1000 ) # Stator coils loop1 = magpy.current.Circle(diameter=0.2, current=1e3, position=(0.1, 0, 0.15)) loop2 = magpy.current.Circle(diameter=0.2, current=-1e3, position=(-0.1, 0, 0.15)) stator = loop1 + loop2 # Collection of currents # Compare torques for different pivot points pivot_points = [None, "centroid", (0, 0, -0.5), (0, 0, -1)] pivot_labels = ["Intrinsic only", "Centroid", "Axle at z=-0.5", "Axle at z=-1.0"] for pivot, label in zip(pivot_points, pivot_labels): F, T = magpy.getFT(stator, cube, pivot=pivot) print(f"{label}:" f" Force: {np.round(F, 1)} N" f" Torque: {np.round(T, 1)} N*m") # Intrinsic only: # Force: [80.6 -0. 0. ] N # Torque: [0. 5.7 0. ] N*m # # Centroid: # Force: [80.6 -0. 0. ] N # Torque: [ 0. 3.4 -0. ] N*m # # Axle at z=-0.5: # Force: [80.6 -0. 0. ] N # Torque: [ 0. 43.7 -0. ] N*m # # Axle at z=-1.0: # Force: [80.6 -0. 0. ] N # Torque: [ 0. 84. -0.] N*m ``` **Key observations:** - The magnetic force is independent of the pivot point selection. It depends only on the relative position and orientation between source and target. - Lever arm effect: Moving the pivot further away from the force application points (each cell center) increases the force moment proportionally ($T = r \times F$). - With `pivot=None` we can compute the intrinsic torque excluding the force moment. This result is unphysical, but useful for comparison with other software. -------------------------------- ## Meshing Force and torque computation is based on numerical discretization of target objects. The `meshing` parameter controls this discretization by specifying the target number of mesh cells, directly affecting both computational accuracy and performance. Higher `meshing` values mean more cells, better accuracy, but longer computation time. ```{hint} Set `meshing` to an integer representing the target number of mesh cells. The algorithms will select optimal cell distributions (homogeneous cells with aspect ratio close to 1) close to (but often not exactly) your target `meshing` value. ``` **Class-specific meshing rules** - `Sphere` and `Dipole`: no meshing required (automatically uses 1 point at center for exact results) - `Cuboid`: accepts also tuple `meshing=(n1, n2, n3)` for explicit subdivision - `Circle`: minimum `meshing=4` (requires at least 4 points around circumference) - `Polyline`: minimum meshing must equal or exceed the number of line segments **Mesh analysis with `meshreport`** Meshing finesse is critical for estimating computation effort and expected precision. By calling `F,T = getFT(..., meshreport=True)` you can print the number of mesh cells of each target. ```python import magpylib as magpy # Source dipole = magpy.misc.Dipole(moment=(0, 0, 1e6)) # Targets c1 = magpy.magnet.Cylinder(dimension=(2, 1), polarization=(0, 0, 1), meshing=1) c2 = magpy.magnet.Cylinder(dimension=(2, 1), polarization=(0, 0, 1), meshing=10) c3 = magpy.magnet.Cylinder(dimension=(2, 1), polarization=(0, 0, 1), meshing=100) # Compute (F,T) AND generate mesh report F, T = magpy.getFT(dipole, [c1, c2, c3], meshreport=True) # Mesh report: # Target Cylinder(id=2589596340240): 2 points # Target Cylinder(id=2589596291216): 20 points # Target Cylinder(id=2589596292176): 124 points ``` **Key insight:** The algorithm selects optimal cell distributions close to (but not exactly) your target `meshing` value. **Mesh analysis with `return_mesh`** Use `meshes = getFT(..., return_mesh=True)` to inspect the generated meshes and verify proper discretization. In this case `getFT()` will return a list of dictionaries with keys `'pts'`, `'moments'` and `'cvecs'` containing the meshes of each target instead of (F, T). ```{code-cell} ipython3 import numpy as np import magpylib as magpy import matplotlib.pyplot as plt # Source dipole = magpy.misc.Dipole(moment=(0, 0, 1e6)) # Targets with increasing mesh density i0 = 1 t1 = magpy.current.Circle(diameter=1, current=i0, meshing=8) t2 = magpy.current.Circle(diameter=2, current=i0, meshing=24) t3 = magpy.current.Circle(diameter=3, current=i0, meshing=64) # Return mesh data instead of forces meshes = magpy.getFT(dipole, [t1, t2, t3], return_mesh=True) # Visualize mesh fig, ax = plt.subplots(subplot_kw={'projection': '3d'}, dpi=130) ax.view_init(elev=90, azim=0) # Top-down view colors = ['red', 'blue', 'forestgreen'] for mesh, color in zip(meshes, colors): # Load mesh points = mesh["pts"] # mesh points cvecs = mesh["cvecs"] / i0 # current path tangent vectors # Draw current segments start = points - cvecs / 2 # starting point of current path ax.quiver( start[:, 0], start[:, 1], start[:, 2], cvecs[:, 0], cvecs[:, 1], cvecs[:, 2], colors=np.random.rand(len(cvecs),3)*.7, # Random dark color for each segment arrow_length_ratio=0 # No arrow head ) # Draw mesh points ax.scatter(*points.T, s=8, c=color) ax.set_title('Mesh Visualization') plt.tight_layout() plt.show() ``` **Key insights:** `return_mesh=True` enables you to check the meshes in detail. Notice the uniform mesh point distribution and the finer finesse with increasing `meshing` parameter. -------------------------------- ## Finite Difference Gradient Magnet force computation requires the magnetic field gradient $\nabla\vec{B}$, which Magpylib calculates using a finite difference scheme. The finite difference step size (absolute value in units (m)) is controlled by the `eps` parameter in `getFT()`. It is `eps=1e-5` by default. **Good value for `eps`** - as large as possible to avoid numerical cancellation - as small as possible for good gradient approximation A good value can be 0.1 % of the smallest length in your system, but it's best to simply try it out: ```python import numpy as np import magpylib as magpy # Set numpy print format np.set_printoptions(formatter={"float": "{:.5e}".format}) # Example: magnet force computation with different eps values cube1 = magpy.magnet.Cuboid(dimension=(1, 1, 1), polarization=(1, 0, 0)) cube2 = cube1.copy(position=(0.6, 0.7, 0.6), meshing=(4, 4, 4)) # Test different eps values eps_values = [1e-2, 1e-3, 1e-4, 1e-5, 1e-20] for eps in eps_values: F, T = magpy.getFT(cube1, cube2, eps=eps) print(F) # [-1.32116e+04 1.96446e+04 4.72297e+04] - 1 correct digit # [-1.21114e+04 1.94403e+04 4.62385e+04] - 2 correct digits # [-1.20999e+04 1.94382e+04 4.62280e+04] - 4 correct digits # [-1.20998e+04 1.94382e+04 4.62279e+04] - 6 correct digits # [ 0.00000e+00 0.00000e+00 0.00000e+00] - no correct digits ``` **Key observations:** A good choice of `eps` balances good gradient approximation with low numerical noise. ```{note} Current objects are not affected by the `eps` parameter as no gradient field is required for their force computation. ``` -------------------------------- ## Convergence In general, the accuracy of the force computation improves with higher `meshing` values. Excluded from this rule are only `Dipole` and `Sphere` objects which give exact results without meshing. The following example shows a computation that converges quickly with increasing mesh finesse. We employ [reciprocity](examples-force-physics-reciprocity) to compare to the exact solution. ```python import numpy as np import magpylib as magpy dipole = magpy.misc.Dipole(moment=(0, 0, 1e6), position=(1, 1, 1)) loop = magpy.current.Circle(diameter=1, current=1) # Forward (exact) F0, T0 = magpy.getFT(loop, dipole, pivot=(0, 0, 0)) # Backward (numerical), has opposite sign for meshing in [8, 40, 200, 1000]: # Set meshing parameter loop.meshing = meshing # Compute backward force and torque F1, T1 = magpy.getFT(dipole, loop, pivot=(0, 0, 0)) # Compute vectorwise relative error errF = np.linalg.norm(F1 + F0) / np.linalg.norm(F0) errT = np.linalg.norm(T1 + T0) / np.linalg.norm(T0) print(f"Meshing: {meshing:>4}, Force Error: {errF:.2e}, Torque Error: {errT:.2e}") # Meshing: 8, Force Error: 1.02e-02, Torque Error: 3.93e-03 # Meshing: 40, Force Error: 3.82e-04, Torque Error: 8.30e-05 # Meshing: 200, Force Error: 1.53e-05, Torque Error: 3.31e-06 # Meshing: 1000, Force Error: 6.12e-07, Torque Error: 1.33e-07 ``` ```{warning} **Convergence behavior**: Accuracy improvement is not always monotonic with mesh refinement. Some intermediate meshing values may temporarily show worse results before eventually converging to the correct solution. **Numerical sensitivity**: Specific combinations of `meshing` and `eps` parameters can occasionally produce unreliable results due to numerical issues. Always test multiple parameter combinations to verify stability. ``` -------------------------------- ## Summary - Best Practices A general ruleset for working properly with `getFT()` avoiding common pitfalls: **Single call of `getFT()`** ```python # Good: Single vectorized call F, T = magpy.getFT(sources, targets) # Computes all forces at once # Avoid: Multiple separate calls for target in targets: F, T = magpy.getFT(sources, target) # Inefficient! ``` **Unless memory considerations:** Memory usage scales with the total number of mesh points. The bottleneck is creating a large observer array (all mesh points + 6 finite difference points per magnet mesh point), and the subsequent call of `getB()`. For large systems, consider breaking into smaller batches. **Exploit reciprocity:** ```python # Efficient: Use Dipole as target (1 mesh point) F, T = magpy.getFT(cube, dipole) # Less efficient: Use complex object as target (many mesh points) F, T = magpy.getFT(dipole, cube) # Requires meshing the cube ``` **Verify convergence:** Ideally, vary both `meshing` and `eps` parameters for testing convergence. ```python # Test multiple mesh densities meshings = [50, 100, 200] for m in meshings: target.meshing = m F, T = magpy.getFT(source, target) # Compare results... ```