Magnetostatics 3D - wrong 3D magnetic vector potential

Dear all,

I am finishing up the 3D case study. I am joining the *.py file generating the 3D mesh and solving the magnetic vector potential A knowing a current source (current density).

I think that I am doing something wrong as the magnetic vector potential is not computed correctly. The second point is that, i am puzzled how to derive from the vector A, the magnetic flux density B in fenicsx, I would need some help on this one which would help to get a better view on the solution.

from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx.io import XDMFFile
from dolfinx.fem import (Expression, Function, FunctionSpace, VectorFunctionSpace)
from dolfinx.fem import (dirichletbc, locate_dofs_topological)
from dolfinx import *

import ufl
from ufl import as_vector, dot, dx, grad, curl, inner

from petsc4py.PETSc import ScalarType
import numpy as np

rank = MPI.COMM_WORLD.rank

import gmsh
import sys
import math

rank = MPI.COMM_WORLD.rank

gmsh.initialize()

gmsh.model.add("magnetostatics3D")

gmsh.option.setColor("Mesh.Color.Lines", 0, 0, 0)

model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 3

R = 0.2
th = 0.05

Ro = 2*(R+th)

lc = th / 4

flag_hexa = 0

if mesh_comm.rank == model_rank:
	cascade = gmsh.model.occ
	meshing = gmsh.model.mesh

	cascade.addPoint(R, 0, -th/2, lc, 1)
	cascade.addPoint(R+th, 0, -th/2, lc, 2)
	cascade.addPoint(R+th, 0, th/2, lc, 3)
	cascade.addPoint(R, 0, th/2, lc, 4)
	cascade.addLine(1, 2, 1)
	cascade.addLine(2, 3, 2)
	cascade.addLine(3, 4, 3)
	cascade.addLine(4, 1, 4)
	cascade.addCurveLoop([1, 2, 3, 4], 1)
	cascade.addPlaneSurface([1], 1)
	
	cascade.synchronize()
	

	if (flag_hexa == 1):
		meshing.setTransfiniteCurve(1, 10)
		meshing.setTransfiniteCurve(2, 10)
		meshing.setTransfiniteCurve(3, 10)
		meshing.setTransfiniteCurve(4, 10)
		meshing.setTransfiniteSurface(1, "Left", [1, 2, 3, 4])
		meshing.setRecombine(2, 1)
	
		# Automatic create a tag, here 1 and 2
		ls_1 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, -math.pi, [20], recombine = True)
		ls_2 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, math.pi, [20], recombine = True)
	else:
		ls_1 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, -math.pi, [20], recombine = False)
		ls_2 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, math.pi, [20], recombine = False)
		
	cascade.fragment(gmsh.model.occ.getEntities(3), [])
	
	cascade.addSphere(0, 0, 0, Ro, 3)
	# Cut out the coil from the sphere keeping a unique interface for both
	cascade.fragment([(3, 3)], [(3, 1) ,(3, 2)])
	
	cascade.removeAllDuplicates
	meshing.removeDuplicateNodes
	meshing.removeDuplicateElements
	
	cascade.synchronize()
	
	up, down = gmsh.model.getAdjacencies(gdim,3)

	gmsh.model.addPhysicalGroup(gdim, [1, 2], 1, name="Coil")
	gmsh.model.setColor(ls_1, 255, 0, 0)
	gmsh.model.setColor(ls_2, 255, 0, 0)
	gmsh.model.addPhysicalGroup(gdim, [3], 2, name="Air")
	gmsh.model.setColor([(gdim, 3)], 0, 0, 255)
	gmsh.model.addPhysicalGroup(gdim-1, [down[0]], 3, name="Boundary")
	gmsh.model.setColor([(gdim-1, down[0])], 0, 0, 255)

	meshing.generate(gdim)

	gmsh.write("magnetostatics3D.msh")

# Read in by dolfinx
mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

# Output DOLFINx meshes to file
with XDMFFile(MPI.COMM_WORLD, "meshCheck.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(facet_tags)
    xdmf.write_meshtags(cell_tags)

# Function space: edge elements
V0 = FunctionSpace(mesh, ("N1curl", 1))

# Create facet to cell connectivity required to determine boundary facets
dofs_facets = locate_dofs_topological(V0, mesh.topology.dim-1, facet_tags.find(3))

# Working boundary condition
u0 = Function(V0)
u0.x.set(0.)
bc = dirichletbc(u0, dofs_facets)

# Trial and test functions
u = ufl.TrialFunction(V0)
v = ufl.TestFunction(V0)

##Define Current Density and Mu0
V1 = VectorFunctionSpace(mesh, ("CG", 1))
J = Function(V1, dtype=ScalarType, name='Current_density')
cells1 = cell_tags.find(1)

J0 = 1e8
def Je(x, J0):
    return np.vstack([-J0*np.sin(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), J0*np.cos(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), np.zeros(x.shape[1])])

J.interpolate(lambda x: Je(x, J0), cells1)
J.x.scatter_forward()

mu = 4*np.pi*10**-7

a = (1 / mu) * inner(curl(u), curl(v)) * dx
l = -inner(J, v) * dx

A = Function(V0)
problem = fem.petsc.LinearProblem(a, l, u = A, bcs = [bc])
problem.solve()


from dolfinx import io
with io.XDMFFile(mesh.comm, "J.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(J)

from dolfinx import io
with io.XDMFFile(mesh.comm, "A.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(A)

Best regards,

Frederic

1 Like

You should consider
https://jsdokken.com/dolfinx-tutorial/chapter3/em.html
and in particular how B is derved from A. In 3D, you can use B=ufl.curl(A) and interpolate that into a suitable FunctionSpace.

Dear Dokken,

I am getting there. I did what you suggested and played around with the linear solver.

I have an issue with the output of the B field as shown below:

There is on top of the expected circular field an additional field in the x direction. I have no idea where it comes from. I am guessing that I am doing something wrong here and likely in the choice of FunctionSpace but it gets a bit too deep in FEM for me at this stage. I am stuck here.

I am joining anew the python code:

from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx.io import XDMFFile
from dolfinx.fem import (Expression, Function, FunctionSpace, VectorFunctionSpace)
from dolfinx.fem import (dirichletbc, locate_dofs_topological)
from dolfinx import *
from dolfinx import io

import ufl
from ufl import as_vector, dot, dx, grad, curl, inner

from petsc4py.PETSc import ScalarType
import numpy as np

rank = MPI.COMM_WORLD.rank

import gmsh
import sys
import math

gmsh.initialize()

gmsh.model.add("magnetostatics3D")

gmsh.option.setColor("Mesh.Color.Lines", 0, 0, 0)

model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 3

Ri = 0.1
th = 0.05

Ra = 0.3

lc = th / 10

flag_hexa = 0

if mesh_comm.rank == model_rank:
	cascade = gmsh.model.occ
	meshing = gmsh.model.mesh

	cascade.addPoint(Ri, 0, -th/2, lc, 1)
	cascade.addPoint(Ri+th, 0, -th/2, lc, 2)
	cascade.addPoint(Ri+th, 0, th/2, lc, 3)
	cascade.addPoint(Ri, 0, th/2, lc, 4)
	cascade.addLine(1, 2, 1)
	cascade.addLine(2, 3, 2)
	cascade.addLine(3, 4, 3)
	cascade.addLine(4, 1, 4)
	cascade.addCurveLoop([1, 2, 3, 4], 1)
	cascade.addPlaneSurface([1], 1)
	
	cascade.synchronize()
	

	if (flag_hexa == 1):
		meshing.setTransfiniteCurve(1, 10)
		meshing.setTransfiniteCurve(2, 10)
		meshing.setTransfiniteCurve(3, 10)
		meshing.setTransfiniteCurve(4, 10)
		meshing.setTransfiniteSurface(1, "Left", [1, 2, 3, 4])
		meshing.setRecombine(2, 1)
	
		# Automatic create a tag, here 1 and 2
		ls_1 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, -math.pi, [40], recombine = True)
		ls_2 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, math.pi, [40], recombine = True)
	else:
		ls_1 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, -math.pi, [40], recombine = False)
		ls_2 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, math.pi, [40], recombine = False)
		
	cascade.fragment(gmsh.model.occ.getEntities(3), [])
	
	cascade.addSphere(0, 0, 0, Ra, 3)
	# Cut out the coil from the sphere keeping a unique interface for both
	cascade.fragment([(3, 3)], [(3, 1) ,(3, 2)])
	
	cascade.removeAllDuplicates
	meshing.removeDuplicateNodes
	meshing.removeDuplicateElements
	
	cascade.synchronize()
	
	up, down = gmsh.model.getAdjacencies(gdim,3)

	gmsh.model.addPhysicalGroup(gdim, [1, 2], 1, name="Coil")
	gmsh.model.setColor(ls_1, 255, 0, 0)
	gmsh.model.setColor(ls_2, 255, 0, 0)
	gmsh.model.addPhysicalGroup(gdim, [3], 2, name="Air")
	gmsh.model.setColor([(gdim, 3)], 0, 0, 255)
	gmsh.model.addPhysicalGroup(gdim-1, [down[0]], 3, name="Boundary")
	gmsh.model.setColor([(gdim-1, down[0])], 0, 0, 255)

	meshing.generate(gdim)

	gmsh.write("magnetostatics3D.msh")

# Read in by dolfinx
mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

# Output DOLFINx meshes to file
with XDMFFile(MPI.COMM_WORLD, "meshCheck.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(facet_tags)
    xdmf.write_meshtags(cell_tags)

# Function space: edge elements
V0 = FunctionSpace(mesh, ("N1curl", 1))

# Create facet to cell connectivity required to determine boundary facets
dofs_facets = locate_dofs_topological(V0, mesh.topology.dim-1, facet_tags.find(3))

# Working boundary condition
u0 = Function(V0)
u0.x.set(0.)
bc = dirichletbc(u0, dofs_facets)

# Trial and test functions
u = ufl.TrialFunction(V0)
v = ufl.TestFunction(V0)

##Define Current Density and Mu0
V1 = VectorFunctionSpace(mesh, ("CG", 1))
J = Function(V1, dtype=ScalarType, name='Current_density')
cells1 = cell_tags.find(1)

J0 = 1e8
def Je(x, J0):
    return np.vstack([-J0*np.sin(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), J0*np.cos(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), np.zeros(x.shape[1])])

J.interpolate(lambda x: Je(x, J0), cells1)
J.x.scatter_forward()

mu = 4*np.pi*1e-7

a = (1 / mu) * inner(curl(u), curl(v)) * dx
l = inner(J, v) * dx

A = Function(V0)
problem = fem.petsc.LinearProblem(a, l, u = A, bcs = [bc], petsc_options = {"ksp_type": "gmres", "ksp_rtol":1e-8, "ksp_atol":1e-12, "ksp_max_it": 1000, "pc_type": "none"})
problem.solve()

B = Function(V1)
B_expr = Expression(ufl.curl(A), V1.element.interpolation_points())
B.interpolate(B_expr)


with io.XDMFFile(mesh.comm, "J.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(J)

with io.XDMFFile(mesh.comm, "A.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(A)
    
with io.XDMFFile(mesh.comm, "B.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(B)

Best,

Frederic

A comes from N1 curl. The curl of something in N1Curl is in Raviart-Thomas, see

Which is a sub-space of DG-1.
Therefore, I would suggest that you create the space for B as VectorFunctionSpace(mesh, ("DG", 1))
and use dolfinx.io.VTXWriter to write the function to file, so you do not get the interpolation on to vertices, which is what happens with XDMFFile

Thanks for the quick answer and the link. I made the change.

I could not get ADIOS2 to compile correctly on my ubuntu 22.04LTS machine to get VTXWriter working. It seems to be an issue of debian systems.

from mpi4py import MPI
from dolfinx.io import gmshio 
from dolfinx.io import XDMFFile
from dolfinx.fem import (Expression, Function, FunctionSpace, VectorFunctionSpace)
from dolfinx.fem import (dirichletbc, locate_dofs_topological)
from dolfinx import *
from dolfinx import io

import ufl
from ufl import as_vector, dot, dx, grad, curl, inner

from petsc4py.PETSc import ScalarType
import numpy as np

rank = MPI.COMM_WORLD.rank

import gmsh
import sys
import math

gmsh.initialize()

gmsh.model.add("magnetostatics3D")

gmsh.option.setColor("Mesh.Color.Lines", 0, 0, 0)

model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 3

Ri = 0.1
th = 0.05
Ra = 0.3
la = 0.5

lc = th / 10
nl = 60

if mesh_comm.rank == model_rank:
	cascade = gmsh.model.occ
	meshing = gmsh.model.mesh

	cascade.addPoint(Ri, 0, -th/2, lc, 1)
	cascade.addPoint(Ri+th, 0, -th/2, lc, 2)
	cascade.addPoint(Ri+th, 0, th/2, lc, 3)
	cascade.addPoint(Ri, 0, th/2, lc, 4)
	cascade.addLine(1, 2, 1)
	cascade.addLine(2, 3, 2)
	cascade.addLine(3, 4, 3)
	cascade.addLine(4, 1, 4)
	cascade.addCurveLoop([1, 2, 3, 4], 1)
	cascade.addPlaneSurface([1], 1)
	
	cascade.addPoint(0., 0., -la/2, lc, 5)
	cascade.addPoint(0., 0., la/2, lc, 6)
	cascade.addLine(5, 6, 5)
	
	cascade.synchronize()

	ls_1 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, -math.pi, [nl], recombine = False)
	ls_2 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, math.pi, [nl], recombine = False)
		
	cascade.fragment(gmsh.model.occ.getEntities(3), [])
	
	cascade.addSphere(0, 0, 0, Ra, 3)
	# Cut out the coil from the sphere keeping a unique interface for both
	cascade.fragment([(3, 3)], [(3, 1) ,(3, 2), (1, 5)])
	
	cascade.removeAllDuplicates
	meshing.removeDuplicateNodes
	meshing.removeDuplicateElements
	
	cascade.synchronize()
	
	up, down = gmsh.model.getAdjacencies(gdim, 3)

	gmsh.model.addPhysicalGroup(gdim, [1, 2], 1, name="Coil")
	gmsh.model.setColor(ls_1, 255, 0, 0)
	gmsh.model.setColor(ls_2, 255, 0, 0)
	gmsh.model.addPhysicalGroup(gdim, [3], 2, name="Air")
	gmsh.model.setColor([(gdim, 3)], 0, 0, 255)
	gmsh.model.addPhysicalGroup(gdim-1, [down[0]], 3, name="Boundary")
	gmsh.model.setColor([(gdim-1, down[0])], 0, 0, 255)
	gmsh.model.addPhysicalGroup(1, [5], 4, name="Axis")

	meshing.generate(gdim)

	gmsh.write("magnetostatics3D.msh")

# Read in by dolfinx
mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

# Output DOLFINx meshes to file
with XDMFFile(MPI.COMM_WORLD, "meshCheck.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(facet_tags)
    xdmf.write_meshtags(cell_tags)

# Function space: edge elements
V0 = FunctionSpace(mesh, ("N1curl", 1))

# Create facet to cell connectivity required to determine boundary facets
dofs_facets = locate_dofs_topological(V0, mesh.topology.dim-1, facet_tags.find(3))

# Working boundary condition
u0 = Function(V0)
u0.x.set(0.)
bc = dirichletbc(u0, dofs_facets)

# Trial and test functions
u = ufl.TrialFunction(V0)
v = ufl.TestFunction(V0)

##Define Current Density and Mu0
V1 = VectorFunctionSpace(mesh, ("CG", 1))
J = Function(V1, dtype=ScalarType, name='Current_density')
cells1 = cell_tags.find(1)

J0 = 1e8
def Je(x, J0):
    return np.vstack([-J0*np.sin(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), J0*np.cos(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), np.zeros(x.shape[1])])

J.interpolate(lambda x: Je(x, J0), cells1)
J.x.scatter_forward()

mu = 4*np.pi*1e-7

a = (1 / mu) * inner(curl(u), curl(v)) * dx
l = inner(J, v) * dx

A = Function(V0)
problem = fem.petsc.LinearProblem(a, l, u = A, bcs = [bc], petsc_options = {"ksp_type": "gmres", "ksp_rtol":1e-8, "ksp_atol":1e-12, "ksp_max_it": 1000, "pc_type": "none"})
problem.solve()

V2 = VectorFunctionSpace(mesh, ("DG", 1))
B = Function(V2)
B_expr = Expression(ufl.curl(A), V2.element.interpolation_points())
B.interpolate(B_expr)

with io.XDMFFile(mesh.comm, "J.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(J)

with io.XDMFFile(mesh.comm, "A.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(A)
    
with io.XDMFFile(mesh.comm, "B.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_tags) 
    xdmf.write_function(B)

#with dolfinx.io.VTXWriter(mesh.comm, "B.bp", [B]) as vtx:
#    vtx.write(0.0)

Best,

Frederic

What do you mean by this? How did you install DOLFINx and how did you try to install adios2?

There are various ways of making adios work on ubuntu22.04.
See

The DOLFINx dockerfile has the instructions on how to install it (and it is based on ubuntu 22.04)
https://github.com/FEniCS/dolfinx/blob/main/docker/Dockerfile#L218

It is preinstalled in

and it is Also part of the conda recipe

https://anaconda.org/conda-forge/fenics-dolfinx

I use the PPA for ubuntu to install fenicsx.

I have tried to compile adios2 from source.

I moved to pyvista which works quite nicely.

Best,

Frederic

I did something with pyvistaqt.

I have only one small issue, I do not know how to plot the norm of the B field on the coil mesh. I tried different ways but I get the full mesh:

# Install gmsh via "pip install gmsh"

from mpi4py import MPI
from dolfinx.io import gmshio 
from dolfinx.io import XDMFFile
from dolfinx.fem import (Expression, Function, FunctionSpace, VectorFunctionSpace)
from dolfinx.fem import (dirichletbc, locate_dofs_topological)
from dolfinx import *
from dolfinx import io

import ufl
from ufl import as_vector, dot, dx, grad, curl, inner, sqrt

from petsc4py.PETSc import ScalarType
import numpy as np

rank = MPI.COMM_WORLD.rank

import gmsh
import sys
import math

gmsh.initialize()

gmsh.model.add("magnetostatics3D")

gmsh.option.setColor("Mesh.Color.Lines", 0, 0, 0)

gmsh.option.setNumber('Mesh.Optimize', 1)
gmsh.option.setNumber('Mesh.OptimizeNetgen', 1)
gmsh.option.setNumber('Mesh.Algorithm3D', 10)

model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 3

Ri = 0.1
th = 0.05
Ra = 0.3
la = 0.5

lc = th / 8
nl = 80

if mesh_comm.rank == model_rank:
	cascade = gmsh.model.occ
	meshing = gmsh.model.mesh

	cascade.addPoint(Ri, 0, -th/2, lc, 1)
	cascade.addPoint(Ri+th, 0, -th/2, lc, 2)
	cascade.addPoint(Ri+th, 0, th/2, lc, 3)
	cascade.addPoint(Ri, 0, th/2, lc, 4)
	cascade.addLine(1, 2, 1)
	cascade.addLine(2, 3, 2)
	cascade.addLine(3, 4, 3)
	cascade.addLine(4, 1, 4)
	cascade.addCurveLoop([1, 2, 3, 4], 1)
	cascade.addPlaneSurface([1], 1)
	
	cascade.addPoint(0., 0., -la/2, lc, 5)
	cascade.addPoint(0., 0., la/2, lc, 6)
	cascade.addLine(5, 6, 5)
	
	cascade.synchronize()

	ls_1 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, -math.pi, [nl], recombine = False)
	ls_2 = cascade.revolve([(2, 1)], 0, 0, 0, 0, 0, 1, math.pi, [nl], recombine = False)
		
	cascade.fragment(gmsh.model.occ.getEntities(3), [])
	
	cascade.addSphere(0, 0, 0, Ra, 3)
	# Cut out the coil from the sphere keeping a unique interface for both
	cascade.fragment([(3, 3)], [(3, 1) ,(3, 2), (1, 5)])
	
	cascade.removeAllDuplicates
	
	cascade.synchronize()
	
	up, down = gmsh.model.getAdjacencies(gdim, 3)

	gmsh.model.addPhysicalGroup(gdim, [1, 2], 1, name="Coil")
	gmsh.model.setColor(ls_1, 255, 0, 0)
	gmsh.model.setColor(ls_2, 255, 0, 0)
	gmsh.model.addPhysicalGroup(gdim, [3], 2, name="Air")
	gmsh.model.setColor([(gdim, 3)], 0, 0, 255)
	gmsh.model.addPhysicalGroup(gdim-1, [down[0]], 3, name="Boundary")
	gmsh.model.setColor([(gdim-1, down[0])], 0, 0, 255)
	gmsh.model.addPhysicalGroup(1, [5], 4, name="Axis")
	
	meshing.removeDuplicateNodes
	meshing.removeDuplicateElements

	meshing.generate(gdim)

#	gmsh.write("magnetostatics3D.msh")

# Read in by dolfinx
mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

print("### MESH: DONE ###")

# Output DOLFINx meshes to file
#with XDMFFile(MPI.COMM_WORLD, "meshCheck.xdmf", "w") as xdmf:
#    xdmf.write_mesh(mesh)
#    xdmf.write_meshtags(facet_tags)
#    xdmf.write_meshtags(cell_tags)

# Function space: edge elements
V0 = FunctionSpace(mesh, ("N1curl", 1))

# Create facet to cell connectivity required to determine boundary facets
dofs_facets = locate_dofs_topological(V0, mesh.topology.dim-1, facet_tags.find(3))

# Working boundary condition
u0 = Function(V0)
u0.x.set(0.)
bc = dirichletbc(u0, dofs_facets)

# Trial and test functions
u = ufl.TrialFunction(V0)
v = ufl.TestFunction(V0)

##Define Current Density and Mu0
V1 = VectorFunctionSpace(mesh, ("CG", 1))
J = Function(V1, dtype=ScalarType, name='Current_density')
cells1 = cell_tags.find(1)

J0 = 1e8
def Je(x, J0):
    return np.vstack([-J0*np.sin(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), J0*np.cos(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), np.zeros(x.shape[1])])

J.interpolate(lambda x: Je(x, J0), cells1)
J.x.scatter_forward()

mu = 4*np.pi*1e-7

a = (1 / mu) * inner(curl(u), curl(v)) * dx
l = inner(J, v) * dx

A = Function(V0)
problem = fem.petsc.LinearProblem(a, l, u = A, bcs = [bc], petsc_options = {"ksp_type": "gmres", "ksp_rtol":1e-8, "ksp_atol":1e-12, "ksp_max_it": 1000, "pc_type": "none"})
problem.solve()

print("### SOLVE: DONE ###")

V2 = VectorFunctionSpace(mesh, ("DG", 1))
B = Function(V2)
B_expr = Expression(ufl.curl(A), V2.element.interpolation_points())
B.interpolate(B_expr)

V3 = FunctionSpace(mesh, ("CG", 1))
Bn = Function(V3)
Bn_expr = Expression(ufl.sqrt(ufl.dot(B, B)), V3.element.interpolation_points())
Bn.interpolate(Bn_expr, cells1)


print("### POSTPROCESS: DONE ###")

#with io.XDMFFile(mesh.comm, "J.xdmf", "w") as xdmf:
#    xdmf.write_mesh(mesh)
#    xdmf.write_meshtags(cell_tags) 
#    xdmf.write_function(J)

#with io.XDMFFile(mesh.comm, "A.xdmf", "w") as xdmf:
#    xdmf.write_mesh(mesh)
#    xdmf.write_meshtags(cell_tags) 
#    xdmf.write_function(A)
#    
#with io.XDMFFile(mesh.comm, "B.xdmf", "w") as xdmf:
#    xdmf.write_mesh(mesh)
#    xdmf.write_meshtags(cell_tags) 
#    xdmf.write_function(B)


import pyvista
from pyvistaqt import BackgroundPlotter
import os
os.environ['PYVISTA_OFF_SCREEN'] = 'false'

# If environment variable PYVISTA_OFF_SCREEN is set to true save a png
# otherwise create interactive plot
if pyvista.OFF_SCREEN:
    pyvista.start_xvfb(wait=0.1)


# Set some global options for all plots
transparent = False
figsize = 800
pyvista.rcParams["background"] = [0.5, 0.5, 0.5]

# Create a plotter with two subplots, and add mesh tag plot to the
# first sub-window
subplotter = BackgroundPlotter(shape=(1, 3))

subplotter.subplot(0, 0)
subplotter.add_text("B field", font_size=14, color="black", position="upper_edge")
# Next, we create a pyvista.UnstructuredGrid based on the mesh
mesh_cells, cell_types, x = plot.create_vtk_mesh(mesh)
gridmesh = pyvista.UnstructuredGrid(mesh_cells, cell_types, x)
# Attach the cells tag data to the pyvita grid
gridmesh.cell_data["Marker"] = cell_tags.values
gridmesh.set_active_scalars("Marker")
# Next, we create a pyvista.UnstructuredGrid based on the mesh
cells, cell_types, x = plot.create_vtk_mesh(V2)
grid = pyvista.UnstructuredGrid(cells, cell_types, x)
# Add this grid (as a wireframe) to the plotter
subplotter.add_mesh(gridmesh, style="wireframe", line_width=2, color="black")
# Create point cloud of vertices, and add the vertex values to the cloud
grid.point_data["B"] = B.x.array.reshape(x.shape[0], V2.dofmap.index_map_bs)
glyphs1 = grid.glyph(orient = "B", factor = 0.04)
subplotter.add_mesh(glyphs1, scalar_bar_args={'title': 'B (T)'})

subplotter.subplot(0, 1)
subplotter.add_text("Coil", font_size=14, color="black", position="upper_edge")
cells, types, x = plot.create_vtk_mesh(mesh, entities=cell_tags.find(1))
sub_mesh = pyvista.UnstructuredGrid(cells, types, x)
subplotter.add_mesh(sub_mesh, style="wireframe", line_width=2, color="black")

cells, cell_types, x = plot.create_vtk_mesh(V3, entities=cell_tags.find(1))
grid2 = pyvista.UnstructuredGrid(cells, cell_types, x)
grid2.point_data["Bn"] = Bn.x.array
grid2.set_active_scalars("Bn")
subplotter.add_mesh(grid2, scalar_bar_args={'title': '|B| (T)'})

subplotter.subplot(0, 2)
subplotter.add_text("Coil", font_size=14, color="black", position="upper_edge")
cells, types, x = plot.create_vtk_mesh(mesh, entities=cell_tags.find(1))
sub_mesh = pyvista.UnstructuredGrid(cells, types, x)
subplotter.add_mesh(sub_mesh, style="wireframe", line_width=2, color="black")

cells, cell_types, x = plot.create_vtk_mesh(V1)
grid3 = pyvista.UnstructuredGrid(cells, cell_types, x)
grid3.point_data["Je"] = J.x.array.reshape(x.shape[0], V1.dofmap.index_map_bs)
glyphs3 = grid3.glyph(orient = "Je", factor = 0.5e-9)
subplotter.add_mesh(glyphs3, scalar_bar_args={'title': 'Je (A/m²)'})

# Save as png if we are using a container with no rendering
if pyvista.OFF_SCREEN:
	subplotter.screenshot("fields.png", transparent_background=transparent, window_size=[figsize, figsize])
else:
	subplotter.show()

subplotter.screenshot("fields.png", transparent_background=transparent, window_size=[figsize, figsize])

subplotter.app.exec_()

The magnitude of the field is similar to what I get with other software. Nice!

Best,

Fredeiic

1 Like