UFL - field norm calculation negative

Question

Hi,

I am trying to map an electric field and its norm through out my volume using the electrical potential previously calculated. I applied the example in the tutorial and I came up with this UFL script:

cell          = tetrahedron
scalarElement = FiniteElement("Lagrange", cell, 1)
vectorElement = VectorElement("Lagrange", cell, 1)
#
E = TrialFunction(vectorElement)
w = TestFunction(vectorElement)

# Coefficients
u = Coefficient(scalarElement)  # previously calculated

# Weak variational form for the electric field
# E = -grad(u)
# Bilinear form
a = inner(E, w) * dx
# Linear form
L = - inner( grad(u), w ) * dx

# Weak variational form for the norm of electric field
# ||E||^2 = dot( grad(u), grad(u) )
E_norm_2 = TrialFunction(scalarElement)
w_norm_2 = TestFunction(scalarElement)
#
# Bilinear form for norm
a_norm_2 = E_norm_2 * w_norm_2 * dx
# Linear form for norm
L_norm_2 = dot( grad(u), grad(u) ) * w_norm_2 * dx

When I run my example with a pretty good precision, the electric field seems okay. On the other hand I have negative values for the electric field's norm. Did I write something not okay?

Best regards,
Yann

Comments

Hi, I think it makes more sense to represent the norm in finite element space of elementwise constant functions. Try

element = FiniteElement('Discontinuous Lagrange', cell, 0)
E_norm_2 = TrialFunction(element)
w_norm_2 = TestFunction(element)

---

Hi MiroK,

Thanks for your response. May I ask you why it makes more sense to represent the norm in finite element space of elementwise constant functions?

I tried your solution, FEniCS is stuck at the system building stage, i.e. it does not start calculating Ax=b ...

Thanks again for your help,
Yann

Answer 1

FFC uses the optimized tensor representation by default when compiling your form.

Unfortunately this representation can result in floating point rounding errors for some forms due to expansion of terms such as (a-b)2 into a2 + b**2 - 2ab.

To get around this problem, either ask ffc to use the "quadrature" representation on the commandline:

ffc -r quadrature MyFile.ufl

or equivalently pass the option as metadata to your form by adding this line to the top of your .ufl file

dxq = Measure("dx", metadata={"representation": "quadrature"})

and use dxq instead of dx in the norm forms. That should fix the negative norms.

Besides that, MiroK's answer is sound advice, which will get you the norm piecewise constant for each cell. You may also want to scale the norm with 1/CellVolume(cell), depending on what you're looking for.

Comments

Hi Martinal,

Your explanation sounds convincing. Unfortunately, my results for the now norm span between ~ [-10, 10] V2/m2.
Any other suggestions?

Best regards,
Yann

---

What does your C++ code look like? Are you actually using a_norm_2 and L_norm_2?