3D magnetic simulations wrong assembling

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

First of all, the example is not reproducible, i.e. no one can copy paste this code onto their system and run it, reproducing your results.

Secondly, create some dolfinx.mesh.MeshTags for the various cells and boundaries (facets) that you are trying to mark, and export them so that you can visualize them in Paraview (with dolfinx.io.XDMFFile.

Please store pts_ferro to a meshtag object and visually inspect it to see that you tag the correct cells.