Dear all,
I am trying to calculate the direction of a static electric field.
The start point is Poisson's equation to determine the potential phi on a 2D domain.
Then I take -grad(phi) to obtain the electric field.
I then calculate the magnitude of the electric field using a projection and project back to plot the scaled field. You'll get the idea from the screenshots.
My problem is the funny artifact at the top middle.
MWE should be this:
import os
from dolfin import *
def convertMesh(geofile, dim=3, noConvert=False):
if geofile[-3:] != "geo":
print("Error: Not a geo file.")
return -1
else:
geofile = geofile[:-4]
if not noConvert:
os.system("gmsh %s.geo -%d -smooth 3 -v 0" % (geofile, dim))
try:
os.system("dolfin-convert %s.msh %s.xml" % (geofile, geofile))
except:
raise
return (("%s.xml" % geofile), ("%s_physical_region.xml" % geofile), ("%s_facet_region.xml" % geofile))
(meshfile, regions, walls) = convertMesh("lcps-250ghz-cross-section.geo", dim=2, noConvert=False)
mesh = Mesh(meshfile)
materials = MeshFunction('size_t', mesh, regions)
boundaries = MeshFunction('size_t', mesh, walls)
nodal_space = FunctionSpace(mesh, 'Lagrange', 3)
L_j = TestFunction(nodal_space)
L_i = TrialFunction(nodal_space)
V1 = 1.0
V2 = -1.0*V1
bc_ground = DirichletBC(nodal_space, Constant(0.0), boundaries, 1)
bc_source_1 = DirichletBC(nodal_space, Constant(V1), boundaries, 2)
bc_source_2 = DirichletBC(nodal_space, Constant((V1+V2)/2), boundaries, 3)
bc_source_3 = DirichletBC(nodal_space, Constant(V2), boundaries, 4)
A = inner(grad(L_i), grad(L_j))*dx
b = Constant(0.0)*L_j*dx
phi = Function( nodal_space )
solve( A == b, phi, [bc_ground, bc_source_1, bc_source_2, bc_source_3] )
plot(phi)
vector_space = FunctionSpace(mesh, 'N1curl', 3)
E = Function(vector_space)
E = -grad(phi)
E_abs = Function(nodal_space)
E_abs = project( sqrt( dot(E, E) ), nodal_space )
plot(project( E/E_abs, vector_space) )
interactive()
The geo-file is available here: https://dl.dropboxusercontent.com/u/2465459/lcps-250ghz-cross-section.geo
PS: This may be considered a follow-up question to: http://fenicsproject.org/qa/3733/strange-behaviour-of-dot-product-of-vector-min-example-incl