Hello everyone,
I’m reaching out for some help or suggestions regarding a magnetic simulation I’ve been working on.
Here’s the context:
I’ve implemented a full pipeline in a class that inherits from a custom DatasetTransform. The input is a sample containing a 2D mesh representing a magnet, along with a specified thickness. The goal is to perform a magnetic simulation in 3D based on this data. tetra creates a system magnet + air and tetrahedralize it. Then I build the fenics simulation, run it and try to visualise the results.
The simulation runs without errors, but when I compute the assembled quantities, the results aren’t physically consistent. For instance, the integrated global magnetization should match magnetization × volume, which is clearly not the case here.
I suspect the issue lies in how I identify or tag the relevant points or regions in the domain — but it’s difficult to verify since PyVista’s visualization is a bit unstable in this case.
Do you have any ideas on what could be going wrong, or suggestions for improving this process?
Thanks a lot in advance! Here is my code:
def tetra(s, thickness, size_box = [100,20], target_edge_length_rect=5):
mesh = s.mesh.to_pyvista()
mesh_extrude = mesh.extrude([0, 0, thickness], capping=True).triangulate()
nm = mesh_extrude.subdivide_adaptive(max_edge_len=0.)
nm = nm.compute_normals(cell_normals=False, point_normals=True, inplace=False)
sxy = size_box[0]//2
sz = size_box[1]//2
pointa = [-sxy, -sxy, sz]
pointb = [-sxy, sxy, sz]
pointc = [sxy, sxy, sz]
rectangle = pv.Rectangle([pointa, pointb, pointc])
rectangle_extrude = rectangle.extrude([0, 0, -2*sz], capping=True).triangulate()
s["extruded_rect"] = read_surfmesh(rectangle_extrude)
remesher_mesh = RemeshSurfaceAdvanced(apply_to="extruded_rect", target_edge_length=target_edge_length_rect)
rect = remesher_mesh(s).extruded_rect.to_pyvista()
tot = rect + nm
tet1 = tetgen.TetGen(tot)
tet1.tetrahedralize(
opt_iterations= 20,
addsteiner_algo=2.0,
minratio=1.1,
nobisect=True,
steinerleft=10000,
opt_max_asp_ratio=2,
mindihedral=20.0,
optmaxdihedral=150.0,
supsteiner_level=4.0,
collinear_ang_tol=170,
verbose=0,
)
grid = tet1.grid # final 3D volumetric mesh
# Step 1: Get the cell centers
cell_centers = grid.cell_centers()
# Step 2: Check which centers are inside the extruded magnet
enclosed = cell_centers.select_enclosed_points(nm, check_surface=True, tolerance=0)
# Step 3: Binary mask per cell (1 = inside, 0 = outside)
mask = enclosed.point_data["SelectedPoints"].astype(int)
# Step 4: Attach the mask as cell data
grid.cell_data["inside_geometry"] = mask
return grid
def locate_ferro(grid, domain):
def pt_in_mesh(x):
xT = x.transpose()
points = np.array([grid.find_closest_point(xx) for xx in xT])
cells = grid.cells
n_cells = cells.shape[0] // 5
cells = cells.reshape((int(n_cells), 5))[:, 1:]
point_flags = np.zeros(points.shape[0], dtype=int)
for i, cell in enumerate(cells):
if grid.cell_data["inside_geometry"][i] == 1:
point_flags[cell] = 1
return point_flags[points].astype(bool)
pts_ferro = locate_entities(domain, domain.topology.dim, pt_in_mesh)
print(pts_ferro.shape)
return pts_ferro
def build_func(Q, V, cells_ferro, val_out, val_in):
# Create a DG-0 function with binary values (0 outside, 1 inside ferro)
func = Function(Q)
func.x.array[:] = val_out
func.x.array[cells_ferro] = val_in
func_ = Function(V)
func_.interpolate(func)
return func_
def compute_for_fenics(computational_mesh):
# 1️⃣ Convert data to FEniCSx format
pts = computational_mesh.points[:, :3]
cells = computational_mesh.cells
# 2️⃣ Reshape cells (tetrahedra)
n_cells = cells.shape[0] // 5
cells = cells.reshape((int(n_cells), 5))[:, 1:]
# 3️⃣ Cast to compatible types
cells = cells.astype(np.int32)
pts = pts.astype(np.float64)
# 4️⃣ Define mesh element for FEniCSx
c_el = ufl.Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3,)))
# 5️⃣ Create the FEniCSx mesh
domain = create_mesh(MPI.COMM_WORLD, cells, pts, c_el)
ldim = 1 # Lagrange element order
V = functionspace(domain, ("Lagrange", ldim))
Q = functionspace(domain, ("DG", 0))
tdim = domain.topology.dim
# Dirichlet BC (phi = 0)
a0 = fem.Function(V)
a0.x.array[:] = 0
# Locate all boundary facets and apply the Dirichlet BC
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = dirichletbc(a0, boundary_dofs)
return Q, V, domain, bc
class MagneticSolverDT_3D(DatasetTransform):
def __init__(self, thickness, magnetization, mu_r, **kwargs):
super().__init__(apply_to=kwargs.pop("apply_to", ""), **kwargs)
self.thickness = thickness
self.mu_r = mu_r
self.magnetization = magnetization
def apply(self, sample):
"""Applies remeshing and magnetic simulation in 3D."""
sample = copy.deepcopy(sample)
grid = tetra(sample, self.thickness)
Q, V, domain, bc = compute_for_fenics(grid)
cells_ferro = locate_ferro(grid, domain)
indicator_ferro = build_func(Q, V, cells_ferro, 0, 1)
M = indicator_ferro * self.magnetization
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v)) * dx
L = M.dx(0) * v * dx
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# Compute H = grad(u)
H = grad(uh)
Hx_f = Function(V)
Hx_f.interpolate(Expression(H[0], V.element.interpolation_points()))
Hx_f.x.array[:] *= indicator_ferro.x.array[:]
# Export to PyVista
u_topology, u_cell_types, u_geometry = vtk_mesh(V)
u_grid = pv.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["i"] = indicator_ferro.x.array.real
u_grid.point_data["u"] = uh.x.array.real
u_grid.point_data["Hx"] = Hx_f.x.array.real
def assemble(f):
compiled_form = form(f)
local_form = assemble_scalar(compiled_form)
return compiled_form.mesh.comm.allreduce(local_form, op=MPI.SUM)
# Physical constants
g = 9.81 # gravity in m/s²
B_earth = 50e-6 # Earth's magnetic field in Tesla
# Compute integrals
integral_Mx = assemble(M * dx)
integral_Hx = assemble(Hx_f * dx)
# Compute effective magnetization and energy
M_eff = np.abs(self.mu_r * integral_Hx - integral_Mx)
mB = -np.abs(M_eff * B_earth) # Effective energy in A·mm²
# Store results in the sample dictionary
sample["scalars.integral_Mx_3d (Amm2)"] = np.array([integral_Mx])
sample["scalars.integral_Hx_3d (Amm2)"] = np.array([integral_Hx])
sample["scalars.mB_3d (µNm)"] = np.array([mB])
sample["scalars.demagnetization_factor_3d (%)"] = np.array([abs(integral_Hx / (integral_Mx + 1e-8)) * 100])
# Extract point cloud inside the magnet
mask = u_grid["i"] > 0.5
selected_points = u_grid.points[mask]
point_cloud = pv.PolyData(selected_points)
point_cloud["u"] = u_grid["u"][mask]
sample["sim_results_3d"] = read_pointcloud(point_cloud)
return sample