can't evaluate gradient at given point in nonlinear poisson demo

Question

I have a feeling there's a simple fix for this problem, but evaluating the gradient function is not working the way I expect in the standard nonlinear-poisson demo. I'm trying to evaluate the gradient at a particular point. I know the solution and the gradient is working fine, since the demo plots it. But I'm trying to access the specific value at (0,0) (I then want to test using that value as a Neumann condition).

So I'm adding

print "u(0,0)=", u((0,0))
gradu = project(grad(u),V)
print "grad u(0,0)=", gradu((0,0))

But projecting the gradient fails with

ufl.log.UFLException: Shape mismatch.

Just taking gradu = grad(u) without projecting also fails, with the call to gradu((0,0)) then returning the error

    component, i = component[:-1], component[-1]
IndexError: tuple index out of range

I would have expected the projection operator to work as a function, even if the direct grad operator does not.

The full code from demo_nonlinear-poisson.py with the grad reference added, is

from dolfin import *

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
File("mesh.pvd") << mesh

V = FunctionSpace(mesh, "CG", 1)

# Define boundary condition
g = Constant(1.0)
bc = DirichletBC(V, g, DirichletBoundary())

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Expression("x[0]*sin(x[1])", degree=2)
F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx

# Compute solution
solve(F == 0, u, bc, solver_parameters={"newton_solver": {"relative_tolerance": 1e-6}})

print "u(0,0)=", u((0,0))
#gradu = grad(u)
gradu = project(grad(u),V)
print "grad u(0,0)=", gradu((0,0))    

# Plot solution and solution gradient
plot(u, title="Solution")
plot(grad(u), title="Solution gradient")
interactive()

What am I doing wrong?

Comments

Hello， I have encountered the same problem. Can you tell me how you solved it?
Best wishes!

Answer 1

You were projecting onto a scalar function space. Instead consider projecting onto a vector function space since $\nabla u \in \mathbb{R}^2$ in this problem:

V_vec = VectorFunctionSpace(mesh, "CG", 1)
gradu = project(grad(u),V_vec)
print "grad u(0,0)=", gradu((0,0))    

Comments

I see it now. Thanks nate.