Merge branch 'feat_mag_face_voxel'

examples/config/tolerance.yaml

@@ -19,6 +19,11 @@ parallel_sweep:
 #   - if the normalized voxel distance is larger than the threshold, numerical approximations are used
 "integral_simplify": 20.0
 
+# control of the magnetic field is computed for the point cloud
+#   - "face" is using the face currents to compute the magnetic field
+#   - "voxel" is using the voxel currents to compute the magnetic field
+"biot_savart": "face"
+
 # method for dense matrix multiplication
 #   - "fft" for doing the multiplication with circulant tensors and FFT
 #   - "direct" for standard matrix multiplication

pypeec/data/schema_tolerance.yaml

@@ -108,6 +108,7 @@
 "required":
     - "parallel_sweep"
     - "integral_simplify"
+    - "biot_savart"
     - "mult_type"
     - "fft_options"
     - "factorization_options"
@@ -129,6 +130,11 @@
     "integral_simplify":
         "type": "number"
         "minimum": 0
+    "biot_savart":
+        "type": "string"
+        "enum":
+            - "face"
+            - "voxel"
     "mult_type":
         "type": "string"
         "enum":

pypeec/lib_solver/extract_solution.py

@@ -26,20 +26,20 @@
 LOGGER = scilogger.get_logger(__name__, "pypeec")
 
 
-def _get_biot_savart(pts, pts_face, I_face):
+def _get_biot_savart(pts, pts_src, I_src):
     """
     Compute the magnetic field at a specified point.
     The field is created by many currents.
     """
 
     # get the distance between the points and the voxels
-    vec = pts-pts_face
+    vec = pts-pts_src
 
     # get the norm of the distance
     nrm = lna.norm(vec, axis=1, keepdims=True)
 
     # compute the Biot-Savart contributions
-    H_all = (1/(4*np.pi))*(np.cross(I_face, vec, axis=1)/(nrm**3))
+    H_all = (1/(4*np.pi))*(np.cross(I_src, vec, axis=1)/(nrm**3))
 
     # sum the contributions
     H_pts = np.sum(H_all, axis=0)
@@ -71,7 +71,7 @@ def _get_magnetic_charge(pts, pts_src, I_src):
     return H_pts
 
 
-def get_magnetic_field_electric(n, d, idx_fc, A_net_c, I_fc, pts_net_c, pts_cloud):
+def get_magnetic_field_electric(n, d, idx_fc, A_net_c, I_fc, pts_net_c, pts_cloud, biot_savart):
     """
     Compute the magnetic field for the provided points (contributions of the electric domains).
     The Biot-Savart law is used for the electric material contribution.
@@ -89,19 +89,31 @@ def get_magnetic_field_electric(n, d, idx_fc, A_net_c, I_fc, pts_net_c, pts_clou
     idx_fy = np.in1d(idx_fc, np.arange(1*nv, 2*nv, dtype=np.int64))
     idx_fz = np.in1d(idx_fc, np.arange(2*nv, 3*nv, dtype=np.int64))
 
-    # get the face positions
-    pts_face = 0.5*np.abs(A_net_c.transpose())*pts_net_c
-
-    # get the face currents
-    I_face = np.zeros((len(I_fc), 3), dtype=np.complex128)
-    I_face[idx_fx, 0] = dx*I_fc[idx_fx]
-    I_face[idx_fy, 1] = dy*I_fc[idx_fy]
-    I_face[idx_fz, 2] = dz*I_fc[idx_fz]
+    if biot_savart == "voxel":
+        # get the voxel positions
+        pts_src = pts_net_c
+
+        # project the faces current into the voxels
+        I_src = np.zeros((len(pts_net_c), 3), dtype=np.complex128)
+        I_src[:, 0] = dx*0.5*np.abs(A_net_c[:, idx_fx])*I_fc[idx_fx]
+        I_src[:, 1] = dy*0.5*np.abs(A_net_c[:, idx_fy])*I_fc[idx_fy]
+        I_src[:, 2] = dz*0.5*np.abs(A_net_c[:, idx_fz])*I_fc[idx_fz]
+    elif biot_savart == "face":
+        # get the face positions
+        pts_src = 0.5*np.abs(A_net_c.transpose())*pts_net_c
+
+        # get the face currents as vector
+        I_src = np.zeros((len(I_fc), 3), dtype=np.complex128)
+        I_src[idx_fx, 0] = dx*I_fc[idx_fx]
+        I_src[idx_fy, 1] = dy*I_fc[idx_fy]
+        I_src[idx_fz, 2] = dz*I_fc[idx_fz]
+    else:
+        raise ValueError("invalid field computation method")
 
     # for each provided point, compute the magnetic field
     H_pts = np.zeros((len(pts_cloud), 3), dtype=np.complex128)
     for i, pts_tmp in enumerate(pts_cloud):
-        H_pts[i, :] = _get_biot_savart(pts_tmp, pts_face, I_face)
+        H_pts[i, :] = _get_biot_savart(pts_tmp, pts_src, I_src)
 
     return H_pts
 

pypeec/run/solver.py

@@ -182,6 +182,7 @@ def _run_solver_sweep(data_solver, data_internal, data_param, sol_init):
     # extract the data
     n = data_solver["n"]
     d = data_solver["d"]
+    biot_savart = data_solver["biot_savart"]
     condition_options = data_solver["condition_options"]
     solver_options = data_solver["solver_options"]
     pts_cloud = data_solver["pts_cloud"]
@@ -305,7 +306,7 @@ def _run_solver_sweep(data_solver, data_internal, data_param, sol_init):
         integral = extract_solution.get_integral(P_fc, P_fm, W_fc, W_fm, S_total)
 
         # get the cloud point magnetic field (contributions of the electric domains)
-        H_pts_c = extract_solution.get_magnetic_field_electric(n, d, idx_fc, A_net_c, I_fc, pts_net_c, pts_cloud)
+        H_pts_c = extract_solution.get_magnetic_field_electric(n, d, idx_fc, A_net_c, I_fc, pts_net_c, pts_cloud, biot_savart)
 
         # get the cloud point magnetic field (contributions of the magnetic domains)
         H_pts_m = extract_solution.get_magnetic_field_magnetic(A_net_m, I_fm, pts_net_m, pts_cloud)