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
