# Magnetostatics 3D - wrong 3D magnetic vector potential

**URL:** <https://fenicsproject.discourse.group/t/magnetostatics-3d-wrong-3d-magnetic-vector-potential/10643>\
**Category:** General\
**Created:** [March 9, 2023, 10:50pm UTC](https://fenicsproject.discourse.group/t/magnetostatics-3d-wrong-3d-magnetic-vector-potential/10643 "2023-03-09T22:50:30Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![ftrillaudp](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ftrillaudp/32/4769_2.png) [@ftrillaudp](https://fenicsproject.discourse.group/u/ftrillaudp)\
**Post date:** [March 9, 2023, 10:50pm UTC](https://fenicsproject.discourse.group/t/magnetostatics-3d-wrong-3d-magnetic-vector-potential/10643/1 "2023-03-09T22:50:30Z")

</div>

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.

```auto
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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [March 10, 2023, 5:23am UTC](https://fenicsproject.discourse.group/t/magnetostatics-3d-wrong-3d-magnetic-vector-potential/10643/2 "2023-03-10T05:23:50Z")

</div>

> [@ftrillaudp](#):
>
> 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.

You should consider  
[https://jsdokken.com/dolfinx-tutorial/chapter3/em.html](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.

---

<div class="post-metadata">

**Author:** ![ftrillaudp](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ftrillaudp/32/4769_2.png) [@ftrillaudp](https://fenicsproject.discourse.group/u/ftrillaudp)\
**Post date:** [March 16, 2023, 5:17pm UTC](https://fenicsproject.discourse.group/t/magnetostatics-3d-wrong-3d-magnetic-vector-potential/10643/3 "2023-03-16T17:17:34Z")

</div>

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:

 ![BField](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/0/039c804b13b02980d029ec12a7aa99b27d9916bb.jpeg)

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:

```auto
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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [March 16, 2023, 5:42pm UTC](https://fenicsproject.discourse.group/t/magnetostatics-3d-wrong-3d-magnetic-vector-potential/10643/4 "2023-03-16T17:42:48Z")

</div>

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

> **[DefElement](https://defelement.com/families/index.html)**
>
> An encyclopedia of finite element definitions

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`

---

<div class="post-metadata">

**Author:** ![ftrillaudp](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ftrillaudp/32/4769_2.png) [@ftrillaudp](https://fenicsproject.discourse.group/u/ftrillaudp)\
**Post date:** [March 16, 2023, 9:20pm UTC](https://fenicsproject.discourse.group/t/magnetostatics-3d-wrong-3d-magnetic-vector-potential/10643/5 "2023-03-16T21:20:33Z")

</div>

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.

```auto
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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [March 16, 2023, 9:55pm UTC](https://fenicsproject.discourse.group/t/magnetostatics-3d-wrong-3d-magnetic-vector-potential/10643/6 "2023-03-16T21:55:43Z")

</div>

> [@ftrillaudp](#):
>
> 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.

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

> **[GitHub - ornladios/ADIOS2: Next generation of ADIOS developed in the Exascale...](https://github.com/ornladios/ADIOS2)**
>
> Next generation of ADIOS developed in the Exascale Computing Program - GitHub - ornladios/ADIOS2: Next generation of ADIOS developed in the Exascale Computing Program

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](https://github.com/FEniCS/dolfinx/blob/main/docker/Dockerfile#L218)

It is preinstalled in

> **[Build software better, together](https://github.com/FEniCS/dolfinx/pkgs/container/dolfinx%2Fdolfinx)**
>
> GitHub is where people build software. More than 100 million people use GitHub to discover, fork, and contribute to over 420 million projects.

and it is Also part of the conda recipe

> <https://github.com/conda-forge/fenics-dolfinx-feedstock/blob/main/recipe/meta.yaml>

[https://anaconda.org/conda-forge/fenics-dolfinx](https://anaconda.org/conda-forge/fenics-dolfinx)

---

<div class="post-metadata">

**Author:** ![ftrillaudp](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ftrillaudp/32/4769_2.png) [@ftrillaudp](https://fenicsproject.discourse.group/u/ftrillaudp)\
**Post date:** [March 17, 2023, 12:35am UTC](https://fenicsproject.discourse.group/t/magnetostatics-3d-wrong-3d-magnetic-vector-potential/10643/7 "2023-03-17T00:35:28Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ftrillaudp](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ftrillaudp/32/4769_2.png) [@ftrillaudp](https://fenicsproject.discourse.group/u/ftrillaudp)\
**Post date:** [March 17, 2023, 1:36am UTC](https://fenicsproject.discourse.group/t/magnetostatics-3d-wrong-3d-magnetic-vector-potential/10643/8 "2023-03-17T01:36:37Z")

</div>

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:

```auto
# 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