I am computing the magnetic field with the following code. The plot of the solution is just a screen with blue overlapped on the 3d grid. Is there a way to investigate what is actually being computed? What the values are like?
Code:
from fenics import *
"""
Simulating Magnetic field from a square loop of wire.
Box is 5x5x1. The wire with lower right hand corner on the origin and at
height z = 0.5.
Equation:
curl (curl A/mu) = 0
curl A = B (what we're looking for)
At the boundary (on the wire), B is defined.
"""
#create mesh and define function space
# box is 5x5x1
length=width= 5.0
height = 0.1
mesh = BoxMesh(Point(0., 0., 0.), Point(length, width, height), 20, 20, 20)
V = VectorFunctionSpace(mesh, 'P', 1)
loop_side = 5.0 #length of square loop
C = 0.1 #constant (~Remanance/4pi) for boundary eqn component
tol = 1E-1#tolerance for boundary definition
loop_height_z = height/2.0 #the sq loop will be at this height
mu = 1.32E-6; #material mu for neodymium magnet
#-----define boundary components-----#
# x component at boundary
b_D_x_str = "2*K*x[2]*(1/pow(pow(l-x[0], 2) + pow(x[2], 2), 2) - 1/pow(pow(x[0], 2) + pow(x[2], 2), 2))"
# y component at boundary
b_D_y_str = "2*K*x[2]*(1/pow(pow(l-x[1], 2) + pow(x[2], 2), 2) - 1/pow(pow(x[1], 2) + pow(x[2], 2), 2))"
# z component at boundary
b_D_z_str = """-2*K*((x[0]-l)/pow(pow(l-x[0],2) + pow(x[2], 2), 2) + x[1]/(pow(pow(x[1], 2) + pow(x[2], 2), 2)) +
x[0]/(pow(pow(x[0], 2) + pow(x[2], 2), 2)) -
(x[1]-l)/pow(pow(l-x[1],2) + pow(x[2], 2), 2))"""
#-----END define boundary components-----#
b_D = Expression((b_D_x_str, b_D_y_str, b_D_z_str), degree=3, l=loop_side, K=C)
""" Checks if vector x is on the boundary: a square loop with lower left hand
corner on the origin. Then shifted up by loop_height_z in the z direction"""
def on_boundary(x):
on_leftright= near(x[0], 0., tol) or near(x[0], loop_side, tol)
on_topbottom = near(x[1], 0., tol) or near(x[1], loop_side, tol)
on_midheight = near(x[2], loop_height_z, tol)
if (near(x[0], 0., tol) or near(x[0], loop_side, tol)) or \
(near(x[1], 0., tol) or near(x[1], loop_side, tol)) and \
(near(x[2], loop_height_z, tol)):
return True
else:
return False
#boundary condition
bc = DirichletBC(V, b_D, on_boundary)
#defining variational problem
# integral over domain (B/mu) * curl(s) = 0
B = TrialFunction(V); #B field
s = TestFunction(V); #test function
f = Constant((0., 0., 0.)) # zero on RHS
a = (1/mu)*dot(B, curl(s))*dx # Left hand side
#F = PETScVector(3) #right hand side
F = inner(f, s)*dx #right hand side
#compute solution
B = Function(V)
solve(a==F, B, bc)
# Plot solution
print(B)
#plot(mesh)
plot(B)
# Save solution to file in VTK formatvtk
vtkfile = File("./output/Bfield/Sqloop.pvd")
vtkfile << B
# Hold plot
from matplotlib import pyplot as plt
from matplotlib import *
axes = plt.gca()
axes.set_zlim(-3, 3)
axes.set_ylim(-1, 7)
axes.set_xlim(-1, 7)
plt.show()