Treatment of non-linear coefficient

1

Hi all,
I have a question regarding the implementaiton of non-linear coefficient such as in the non-linear Poisson equation example:

def q(u):
“Return nonlinear coefficient”
return 1 + u**2

My question is how FEniCS evaluates “q(u)” if it is in a form expression such as

a = inner(q(u), v)*dx

My understanding is, given
u = \sum_{j=1}^{N}U_j \phi_j and q(u) = 1 + u^2,
q(u) = 1 + (\sum_{j=1}^{N}U_j \phi_j)^2.

Does FEniCS evaluate the q(u) in the above way?
Or does FEniCS do simpler calculation like
q(u) = 1 + \sum_{j=1}^{N}(1 + U_j^2) \phi_j
?

2

FEniCS will evaluate the first one at quadrature points within the element. You can verify this with the following quick test

from dolfin import *

mesh = UnitIntervalMesh(1)
V = FunctionSpace(mesh,"CG",1)
u = project(SpatialCoordinate(mesh)[0],V)

# One quadrature point in the middle of the element:
dx = dx(metadata={"quadrature_degree":0})

print(assemble((u**2)*dx))

which prints the midpoint evaluation of 0.25 when using one quadrature point at x = 0.5. The exact integral, without restricting quadrature, is 1/3, while interpolating at nodes first would give 0.5.

Note that nodal interpolation is used for Expressions in variational forms, e.g., continuing the example above with

expr = Expression("u*u",u=u,degree=1)
print(assemble(expr*dx(domain=mesh)))

prints 0.5, as expected from nodal interpolation of u^2 into a space of degree 1 (as indicated by the Expression constructor).

2 Likes
3

Hi Kamensky,

Thank you for your detailed explaination. It now makes sense! I appreciate your help!