Merge branch 'main' of https://github.com/magpylib/magpylib-material-…

…response

.gitignore

@@ -1,3 +1,6 @@
+__temp*
+
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]

CHANGELOG.md

@@ -10,6 +10,14 @@
     - if susceptiblity input is scalar, isotropic susceptibility is assumed, if it is a 3-vector it can be anisotropic
 * Various tests included of interface and computation, isotropic and anisotropic tests confirm computaiton
 
+# Unreleased
+* Improve internals
+* anisotropic susceptibilities are now allowed.
+* Improve suszeptibility input possibilities:
+    - give susceptibility to parent collection
+    - if susceptiblity input is scalar, isotropic susceptibility is assumed, if it is a 3-vector it can be anisotropic
+* Various tests included of interface and computation, isotropic and anisotropic tests confirm computaiton
+
 # 0.2.1a0
 * Fix null polarization for rotated objects ([#7](https://github.com/magpylib/magpylib-material-response/pull/7))
 * Fix docs not building ([#6](https://github.com/magpylib/magpylib-material-response/pull/6))

magpylib_material_response/__init__.py

@@ -2,4 +2,4 @@
 analysis for soft materials and demagnetization of hard magnets. Leveraging the Method
 of Moments, it calculates magnetic material response with high precision."""
 
-__version__ = "0.2.1a0"
+__version__ = "0.2.2"

magpylib_material_response/demag.py

@@ -31,21 +31,71 @@
 logger.configure(**config)
 
 
-def get_susceptibilities(*sources, susceptibility=None):
+def get_susceptibilities(sources, susceptibility):
     """Return a list of length (len(sources)) with susceptibility values
     Priority is given at the source level, hovever if value is not found, it is searched
     up the parent tree, if available. Raises an error if no value is found when reached
     the top level of the tree."""
-    susceptibilities = []
+    n = len(sources)
+
+    # susceptibilities from source attributes
+    if susceptibility is None:
+        susis = []
+        for src in sources:
+            susceptibility = getattr(src, "susceptibility", None)
+            if susceptibility is None:
+                if src.parent is None:
+                    raise ValueError(
+                        "No susceptibility defined in any parent collection"
+                    )
+                susis.extend(get_susceptibilities(src.parent))
+            elif not hasattr(susceptibility, "__len__"):
+                susis.append((susceptibility, susceptibility, susceptibility))
+            elif len(susceptibility) == 3:
+                susis.append(susceptibility)
+            else:
+                raise ValueError("susceptibility is not scalar or array fo length 3")
+    else:
+        # susceptibilities as input to demag function
+        if np.isscalar(susceptibility):
+            susis = np.ones((n, 3)) * susceptibility
+        elif len(susceptibility) == 3:
+            susis = np.tile(susceptibility, (n, 1))
+            if n == 3:
+                raise ValueError(
+                    "Apply_demag input susceptibility is ambiguous - either scalar list or vector single entry. "
+                    "Please choose different means of input or change the number of cells in the Collection."
+                )
+        else:
+            if len(susceptibility) != n:
+                raise ValueError(
+                    "Apply_demag input susceptibility must be scalar, 3-vector, or same length as input Collection."
+                )
+            susis = np.array(susceptibility)
+            if susis.ndim == 1:
+                susis = np.repeat(susis, 3).reshape(n, 3)
+
+    susis = np.reshape(susis, 3 * n, order="F")
+    return np.array(susis)
+
+
+def get_H_ext(*sources, H_ext=None):
+    """Return a list of length (len(sources)) with H_ext values
+    Priority is given at the source level, hovever if value is not found, it is searched up the
+    the parent tree, if available. Sets H_ext to zero if no value is found when reached the top
+    level of the tree"""
+    H_exts = []
     for src in sources:
-        susceptibility = getattr(src, "susceptibility", None)
-        if susceptibility is None:
+        H_ext = getattr(src, "H_ext", None)
+        if H_ext is None:
             if src.parent is None:
-                raise ValueError("No susceptibility defined in any parent collection")
-            susceptibilities.extend(get_susceptibilities(src.parent))
+                # print("Warning: No value for H_ext defined in any parent collection. H_ext set to zero.")
+                H_exts.append((0.0, 0.0, 0.0))
+            else:
+                H_exts.extend(get_H_ext(src.parent))
         else:
-            susceptibilities.append(susceptibility)
-    return susceptibilities
+            H_exts.append(H_ext)
+    return H_exts
 
 
 def demag_tensor(
@@ -341,14 +391,17 @@ def apply_demag(
         )  # shape ii = x1, ... xn, y1, ... yn, z1, ... zn
 
         # set up S
-        if susceptibility is None:
-            susceptibility = get_susceptibilities(*magnets_list)
-        susceptibility = np.array(susceptibility)
-        if len(susceptibility) != n:
+        sus = get_susceptibilities(magnets_list, susceptibility)
+        S = np.diag(sus)  # shape ii, jj
+
+        # set up H_ext
+        H_ext = get_H_ext(*magnets_list)
+        H_ext = np.array(H_ext)
+        if len(H_ext) != n:
             raise ValueError(
-                "Apply_demag input collection and susceptibility must have same length."
+                "Apply_demag input collection and H_ext must have same length."
             )
-        S = np.diag(np.tile(susceptibility, 3))  # shape ii, jj
+        H_ext = np.reshape(H_ext, (3 * n, 1), order="F")
 
         # set up T (3 pol unit, n cells, n positions, 3 Bxyz)
         with timelog("Demagnetization tensor calculation", min_log_time=min_log_time):
@@ -362,7 +415,7 @@ def apply_demag(
             T *= magpy.mu_0
             T = T.swapaxes(2, 3).reshape((3 * n, 3 * n)).T  # shape ii, jj
 
-        pol_tolal = pol_magnets
+        pol_total = pol_magnets
 
         if currents_list:
             with timelog(
@@ -371,14 +424,14 @@ def apply_demag(
                 pos = np.array([src.position for src in magnets_list])
                 pol_currents = magpy.getB(currents_list, pos, sumup=True)
                 pol_currents = np.reshape(pol_currents, (3 * n, 1), order="F")
-                pol_tolal += np.matmul(S, pol_currents)
+                pol_total += np.matmul(S, pol_currents)
 
         # set up Q
         Q = np.eye(3 * n) - np.matmul(S, T)
 
         # determine new polarization vectors
         with timelog("Solving of linear system", min_log_time=1):
-            pol_new = np.linalg.solve(Q, pol_tolal)
+            pol_new = np.linalg.solve(Q, pol_total + np.matmul(S, H_ext))
 
         pol_new = np.reshape(pol_new, (n, 3), order="F")
         # pol_new *= .4*np.pi

tests/test_basic.py

@@ -1,5 +1,48 @@
+import magpylib as magpy
+import numpy as np
+
 import magpylib_material_response
+from magpylib_material_response.demag import apply_demag
+from magpylib_material_response.meshing import mesh_Cuboid
 
 
 def test_version():
     assert isinstance(magpylib_material_response.__version__, str)
+
+
+def test_susceptibility_inputs():
+    """
+    test if different xi inputs give the same result
+    """
+
+    zone = magpy.magnet.Cuboid(
+        dimension=(1, 1, 1),
+        polarization=(0, 0, 1),
+    )
+    mesh = mesh_Cuboid(zone, (2, 2, 2))
+
+    dm1 = apply_demag(mesh, susceptibility=4)
+    dm2 = apply_demag(mesh, susceptibility=(4, 4, 4))
+    dm3 = apply_demag(mesh, susceptibility=[4] * 8)
+    dm4 = apply_demag(mesh, susceptibility=[(4, 4, 4)] * 8)
+
+    zone = magpy.magnet.Cuboid(
+        dimension=(1, 1, 1),
+        polarization=(0, 0, 1),
+    )
+    zone.susceptibility = 4
+    mesh = mesh_Cuboid(zone, (2, 2, 2))
+    dm5 = apply_demag(mesh)
+
+    zone = magpy.magnet.Cuboid(
+        dimension=(1, 1, 1),
+        polarization=(0, 0, 1),
+    )
+    zone.susceptibility = (4, 4, 4)
+    mesh = mesh_Cuboid(zone, (2, 2, 2))
+    dm6 = apply_demag(mesh)
+
+    b1 = dm1.getB((1, 2, 3))
+    for dm in [dm2, dm3, dm4, dm5, dm6]:
+        bb = dm.getB((1, 2, 3))
+        np.testing.assert_allclose(b1, bb)

tests/test_isotropic_anisotropic.py

@@ -0,0 +1,68 @@
+import magpylib as magpy
+import numpy as np
+
+from magpylib_material_response import demag, meshing
+
+
+def test_isotropic_susceptibility():
+
+    cells = 1000  # should be >=1000, otherwise discretization error too large
+
+    magnet = magpy.magnet.Cuboid(dimension=(1e-3, 1e-3, 1e-3), polarization=(0, 0, 1.1))
+    grid = np.loadtxt("tests/testdata/grid_points.pts")
+    field_ansys = np.loadtxt("tests/testdata/isotropic_results_ansys.txt", skiprows=1)
+    field_ansys = field_ansys[:, 3:]
+
+    # isotropic
+    magnet.susceptibility = 0.1
+    magnet_meshed = meshing.mesh_Cuboid(magnet, cells)
+
+    demag.apply_demag(magnet_meshed, inplace=True)
+
+    field_magpylib = magnet_meshed.getB(grid)
+
+    np.testing.assert_allclose(field_ansys, field_magpylib, rtol=0, atol=0.0012)
+
+
+def test_anisotropic_susceptibility():
+
+    cells = 1000  # should be >=1000, otherwise discretization error too large
+
+    magnet = magpy.magnet.Cuboid(dimension=(1e-3, 1e-3, 1e-3), polarization=(0, 0, 1.1))
+    grid = np.loadtxt("tests/testdata/grid_points.pts")
+    field_ansys = np.loadtxt("tests/testdata/anisotropic_results_ansys.txt", skiprows=1)
+    field_ansys = field_ansys[:, 3:]
+
+    # anisotropic
+    magnet.susceptibility = (0.3, 0.2, 0.1)
+    magnet_meshed = meshing.mesh_Cuboid(magnet, cells)
+
+    demag.apply_demag(magnet_meshed, inplace=True)
+
+    field_magpylib = magnet_meshed.getB(grid)
+
+    np.testing.assert_allclose(field_ansys, field_magpylib, rtol=0, atol=0.0012)
+
+
+def test_negative_susceptibility():
+
+    cells = 1000  # should be >=1000, otherwise discretization error too large
+
+    magnet = magpy.magnet.Cuboid(
+        dimension=(1e-3, 1e-3, 1e-3), polarization=(0, 0, -0.1)
+    )
+    grid = np.loadtxt("tests/testdata/grid_points.pts")
+    field_ansys = np.loadtxt(
+        "tests/testdata/negative_susceptibility_ansys.txt", skiprows=1
+    )
+    field_ansys = field_ansys[:, 3:]
+
+    # isotropic
+    magnet.susceptibility = -1.1
+    magnet_meshed = meshing.mesh_Cuboid(magnet, cells)
+
+    demag.apply_demag(magnet_meshed, inplace=True)
+
+    field_magpylib = magnet_meshed.getB(grid)
+
+    np.testing.assert_allclose(field_ansys, field_magpylib, rtol=0, atol=0.0065)