Hello everyone,
I’m trying to set up a simulation for a magnetic field problem, and I’m encountering some issues. Here’s the setup:
I have a 3D magnet, which corresponds to the extrusion of a 2D surface, and as a result, I’m reducing the problem to two dimensions. My mesh is a PyVista mesh, and I’ve expanded it to include a large rectangular domain that represents both the magnet and the surrounding air. The materials (magnet and air) are defined using cell_data in PyMesh.
Now, I’m trying to solve for the magnetic potential in this mesh using DolfinX, but I’ve simplified things by ignoring the permeability for now. My boundary conditions are as follows:
- Dirichlet condition: potential = 0 at infinity (which corresponds to the boundary of my overall mesh).
- Potentially: I’m enforcing potential = 0 at y=0 to impose symmetry along the y-axis (since the magnet is symmetric).
I’ve merged ideas from tutorials and forum posts, but I’m running into issues. The code runs without errors, but the output isn’t correct. It seems like the boundary conditions aren’t being respected, and I’m struggling to figure out why. I’m not sure what’s causing the incorrect output field, and I could use some insights on what I might be doing wrong.
import numpy as np # For vector and matrix manipulation
# Visualization
import pyvista
from mpi4py import MPI
import matplotlib.pyplot as plt
# FEniCS
from ufl import (SpatialCoordinate, TestFunction, TrialFunction, dx, dot, grad, div, diff, inner)
import basix
from petsc4py import PETSc
import ufl
from dolfinx import geometry
from dolfinx.mesh import create_mesh, locate_entities_boundary, exterior_facet_indices
from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, assemble_scalar, form)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.plot import vtk_mesh
from dolfinx import fem
#computational_mesh is a pyvista mesh already computed
# 1️⃣ Convert data into format usable by FEniCSx
pts = computational_mesh.points[:, :2] # Only keep x, y coordinates
cells = computational_mesh.faces # Extract the faces
# 2️⃣ Reshape cells (triangles)
n_cells = cells.shape[0] // 4 # Reshape the cells with a prefix (3, i1, i2, i3)
cells = cells.reshape((int(n_cells), 4))[:, 1:]
# 3️⃣ Convert data to FEniCSx compatible types
cells = cells.astype(np.int32)
pts = pts.astype(np.float64)
# 4️⃣ Define triangular mesh element for FEniCSx
c_el = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(2,)))
# 5️⃣ Create the FEniCSx mesh
domain = create_mesh(MPI.COMM_WORLD, cells, pts, c_el)
ldim = 2 # Lagrange element order
V = functionspace(domain, ("Lagrange", ldim))
tdim = domain.topology.dim
# Dirichlet boundary condition (phi=0)
a0 = fem.Function(V)
a0.x.array[:] = 0
# Find the boundary nodes and define that all of them belong to the Dirichlet boundary (phi=0)
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)
# Localize facets where y = 0 (forcing axisymmetry condition)
dofs_y0 = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[1], 0))
# Apply Dirichlet boundary condition (phi = 0) for y = 0
bc_y0 = dirichletbc(a0, dofs_y0)
del a0
Q = functionspace(domain, ("DG", 0))
mu_r = Function(Q)
cell_labels = computational_mesh.cell_data["source_label"]
cells_ferro = cell_labels == 1 # Cells where magnetization is 1
Hc = 1000
Mc = Constant(domain, PETSc.ScalarType(Hc))
# Defining the permeability values
mu_r.x.array[cells_ferro] = muf
mu_r.x.array[~cells_ferro] = muair
# Magnitude of the magnetization vector
M = Function(Q)
M.x.array[:] = 0
M.x.array[cells_ferro] = Mc
M_ = Function(V)
M_.interpolate(M)
# Trial and test functions
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_y0, bc])
uh = problem.solve()
# Display the solution
u_topology, u_cell_types, u_geometry = vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u"] = uh.x.array.real
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=False)
u_plotter.add_mesh(initial_mesh.extract_feature_edges(boundary_edges=True), color="blue", show_edges=True, line_width=2, label="Initial Mesh")
u_plotter.show()