Magnetostatics example in dolfinx

dolfinx
1

I’ve been spending a lot of time recently trying to use FEniCS for electromagnetics simulation, and have managed to get the magnetostatics demo working, using pygmsh for the mesh generation and dolfinx for solving. I am new to FEniCS and FEM simulation so I’m not sure all of my modifications were valid, but I figured I would post the code in any case if it is useful for others. I haven’t compared the result to the working example in dolfin (to do) but the magnetic vector potential result looks superficially similar.

Mesh generation:

from pygmsh.opencascade import Geometry
from pygmsh import generate_mesh
import meshio
import numpy as np


a = 1.0   # inner radius of iron cylinder
b = 1.2   # outer radius of iron cylinder
c_1 = 0.8 # radius for inner circle of copper wires
c_2 = 1.4 # radius for outer circle of copper wires
r = 0.1   # radius of copper wires
R = 5.0   # radius of domain
n = 10    # number of windings
lc1 = 0.1 #Maximum length of mesh
lc2 = 0.05
lc3 = 0.03

geom = Geometry()

center_point = [0,0,0]

# Define geometry for wires (N = North (up), S = South (down))
angles_N = [i*2*np.pi/n for i in range(n)]
angles_S = [(i + 0.5)*2*np.pi/n for i in range(n)]
wires_N = [geom.add_disk([c_1*np.cos(v), c_1*np.sin(v),0], r, char_length=lc3) for v in angles_N]
wires_S = [geom.add_disk([c_2*np.cos(v), c_2*np.sin(v), 0], r, char_length=lc3) for v in angles_S]
wires_N = geom.boolean_union(wires_N)
wires_S = geom.boolean_union(wires_S)

circle1 = geom.add_disk([0,0,0], b, char_length=lc2)
circle2 = geom.add_disk(center_point, a, char_length=lc2)
cylinder = geom.boolean_difference([circle1], [circle2])
#wires = wires_S+wires_N

domain = geom.add_disk(center_point, R, char_length=lc1)
domain_sub = geom.boolean_difference([domain], [wires_S]+[wires_N]+[cylinder], delete_other=False)
system = geom.boolean_fragments([domain_sub], [cylinder]+[wires_S]+[wires_N])

geom.add_raw_code("Physical Surface(1) = {" + cylinder.id + "};") 
geom.add_raw_code("Physical Surface(2) = {" + wires_N.id + "};")
geom.add_raw_code("Physical Surface(3) = {" + wires_S.id + "};")
geom.add_raw_code("Physical Surface(4) = {" + domain_sub.id + "};")
geom.add_raw_code("Physical Curve(1) = {24};")

mesh = generate_mesh(geom, geo_filename = "test.geo", prune_z_0 = True)

# What follows is this (strange) mesh conversion stuff ... . ;-)
line_cells     = None
line_data     = None
triangle_cells = None
triangle_data  = None

for key in mesh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = mesh.cell_data_dict["gmsh:physical"][key]
        print(set(triangle_data))
    elif key == "line":
        line_data = mesh.cell_data_dict["gmsh:physical"][key]

for cell in mesh.cells:
    if cell.type == "line":
        if line_cells is None:
            line_cells = cell.data
        else:
            line_cells = np.vstack([line_cells, cell.data])
    elif cell.type == "triangle":
        if triangle_cells is None:
            triangle_cells = cell.data
        else:
            triangle_cells = np.vstack([triangle_cells, cell.data])

line_mesh = meshio.Mesh(points=mesh.points,
                         cells={"line": line_cells},
                         cell_data={"name_to_read":[line_data]})
triangle_mesh =meshio.Mesh(points=mesh.points,
                           cells=[("triangle", triangle_cells)],
                           cell_data={"name_to_read":[triangle_data]})
meshio.xdmf.write("line_mesh.xdmf", line_mesh)
meshio.xdmf.write("triangle_mesh.xdmf", triangle_mesh)

Solver:

import dolfinx
import numpy as np

from dolfinx import Mesh, geometry, Function, DirichletBC, Constant, solve, FunctionSpace
from dolfinx.io import XDMFFile
from dolfinx.plotting import plot

import ufl
from ufl import TrialFunction, TestFunction, dot, grad, dx, inner, Measure

from mpi4py import MPI
comm = MPI.COMM_WORLD


with XDMFFile(MPI.COMM_WORLD,"triangle_mesh.xdmf",'r') as infile:
    mesh = infile.read_mesh(dolfinx.cpp.mesh.GhostMode.none, name="Grid")

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(MPI.COMM_WORLD, "triangle_mesh.xdmf", "r") as xdmf_infile:
    mt = xdmf_infile.read_meshtags(mesh, name="Grid")

with XDMFFile(MPI.COMM_WORLD, "line_mesh.xdmf", "r") as xdmf_infile:
    mt_line = xdmf_infile.read_meshtags(mesh, name="Grid")

n=10 #ten wires on each side of the cylinder
# Define function space
CG = ufl.FiniteElement("P", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, CG)

# Define boundary condition
inner_facets = mt_line.indices[mt_line.values == 1]
dofs_outer = dolfinx.fem.locate_dofs_topological(V, 1, inner_facets)



ubcOuter = Function(V)
ubcOuter.vector.set(0)
bc = DirichletBC(ubcOuter, dofs_outer)
#bc = DirichletBC(V, Constant(mesh, 0), 'on_boundary')

# Define subdomain markers and integration measure
dx = Measure('dx', domain=mesh, subdomain_data=mt)

# Define current densities
J_N = Constant(mesh, 1.0)
J_S = Constant(mesh, -1.0)

#Locate wires
facets_N = mt.indices[mt.values == 2]
dofs_N = dolfinx.fem.locate_dofs_topological(V, 2, facets_N)

facets_S = mt.indices[mt.values == 3]
dofs_S = dolfinx.fem.locate_dofs_topological(V, 2, facets_S)

facets_cylinder = mt.indices[mt.values == 1]
dofs_cylinder = dolfinx.fem.locate_dofs_topological(V, 2, facets_cylinder)

#determine if coordinate is in wire:

#First generate the list of coordinates that is within wire_N
dof_array = V.dofmap.list.array()
coordinates = V.tabulate_dof_coordinates()

coords_cylinder = coordinates[dofs_cylinder[:,0]]
coords_N = coordinates[dofs_N[:,0]]
coords_S = coordinates[dofs_S[:,0]]
coords_copper = np.vstack((coords_N, coords_S))

# Define magnetic permeability
def permeability(x):
    #initialize vector to be returned
    perm = np.full(x.shape[1], 0, dtype=np.float64)
    for i in range(x.shape[1]):
        #Wire
        if np.any(np.all(x[:,i]==coords_copper, axis=1)):
            perm[i] = 1.26e-3
        elif np.any(np.all(x[:,i]==coords_cylinder, axis=1)):
            perm[i] = 1e-5 # Iron (but not right)
        else:
            perm[i] = 4*np.pi*1e-7
    print(perm)
    return perm

mu = Function(V)
mu.interpolate(permeability)

#def in_wire_N(x):
#    #Generate boolean array
#    for i in x.shape[1]:
#        if x[:,] in coords_N



# Define variational problem
A_z = TrialFunction(V)
v = TestFunction(V)
a = (1 / mu)*inner(grad(A_z), grad(v))*dx
L_N = inner(J_N,v)*dx(2)
L_S = inner(J_S,v)*dx(3)
L = L_N + L_S

# Solve variational problem
A_z = Function(V)
solve(a == L, A_z, bc)

with XDMFFile(MPI.COMM_WORLD, "A_z.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(A_z)

# Compute magnetic field (B = curl A)
#W = VectorFunctionSpace(mesh, 'P', 1)
#B = project(as_vector((A_z.dx(1), -A_z.dx(0))), W)

The most hacky part of the demo is the interpolation of mu, I couldn’t find an elegant way in dolfinx to have mu be determined topologically rather than by coordinate. So until I find such a utility function, I have a search function that maps each coordinate to a subdomain. Is there an elegant way to do this in dolfinx?

My last question is about the project() function. I couldn’t find it in dolfinx, is there a way to achieve the same functionality currently?

2

For projection, you can consider the approach taken in dolfiny by either copying this approach or installing this convenience library.

Also it is beneficial to have a look at:

Finally, I have made a minimal example highlighting the combination of these operations:

import dolfinx
import dolfinx.io

import ufl
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 5,5, dolfinx.cpp.mesh.CellType.quadrilateral)
num_cells_local = mesh.topology.index_map(mesh.topology.dim).size_local
# Create a meshtag for the first and last local indices of the mesh
indices = np.array([0,1,2,3, num_cells_local-4, num_cells_local-3, num_cells_local-2, num_cells_local-1], dtype=np.int32)
values = np.array([1,1,1,1,2,2,2,2], dtype=np.int32)
mt = dolfinx.MeshTags(mesh, mesh.topology.dim, indices, values)


# Create a DG function for these cells
V = dolfinx.FunctionSpace(mesh, ("DG",0))
v = dolfinx.Function(V)

for i in range(num_cells_local):
    if i in mt.indices:
        local_index = np.flatnonzero(mt.indices == i)[0]
        if mt.values[local_index] == 1:
            v.vector.setValueLocal(i, 1.1)
        elif mt.values[local_index] == 2:
            v.vector.setValueLocal(i, 2.1)
    else:
        v.vector.setValueLocal(i, 0.1)


# Project DG function to CG space
W= dolfinx.FunctionSpace(mesh, ("CG",1))
u, w = ufl.TrialFunction(W), ufl.TestFunction(W)
a= ufl.inner(u, w) * ufl.dx
l = ufl.inner(v, w) * ufl.dx
uh = dolfinx.Function(W)
A = dolfinx.fem.assemble_matrix(a)
A.assemble()
b = dolfinx.fem.assemble_vector(l)

solver = PETSc.KSP().create(MPI.COMM_WORLD)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
solver.setOperators(A)
solver.solve(b, uh.vector)
uh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

# Write mesh tags and projected function to output files
xdmf_file = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "w")
xdmf_file.write_mesh(mesh)
mesh.topology.create_connectivity(mesh.topology.dim-1,mesh.topology.dim)
xdmf_file.write_meshtags(mt)
xdmf_file.close()
dolfinx.io.VTKFile("projectedDG.pvd").write(uh)
3 Likes