This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Strange artifact when normalising vector field (with MWE)

0 votes

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.

Potential phi
Result of E/E_abs

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

asked May 13, 2016 by cweickhmann FEniCS Novice (550 points)
...