Hi everyone,

I’m solving Poisson’s equation with a source term defined by an image on a structured mesh. This runs fine (even though there could still be mistakes in my approach), even with a random noise as given in this minimal code. If I’m not mistaken, the way I define the variational problem should impose Neumann conditions on the edge. However by projecting my gradient I get something quite far from `du/dn = 0`

. Is there a way to get a “better” gradient? Here is my minimal code :

```
from dolfin import *
import numpy as np
from scipy.interpolate import *
h = 60
w = 50
loss = np.random.rand(h,w)
loss_cont = RectBivariateSpline(np.linspace(0.5,h-0.5,h,dtype = np.float64),np.linspace(0.5,w-0.5,w,dtype = np.float64),loss)
# Create mesh and define function space
my_mesh = RectangleMesh(Point(0,0),Point(h,w),h,w)
V = FunctionSpace(my_mesh, 'Lagrange', 1)
vxmap = vertex_to_dof_map(V)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V)
f = fem.interpolation.interpolate(f, V)
ks = np.array([i % (w+1) for i in range(V.dim())])
ls = np.array([int(i / (h+1)) for i in range(V.dim())])
f.vector()[vxmap[:]] = loss_cont(ls,ks,grid=False)[:] # Define source term at its nodes
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u)
# Gradient
W = VectorFunctionSpace(my_mesh, "CG", 1)
gradu = project(grad(u), W)
wxmap = vertex_to_dof_map(W)
my_grad = np.array(gradu.vector().get_local()[wxmap]).reshape(h+1,w+1,2)
print(my_grad[:,0])
# Setting the borders
my_grad[:,0,1] = 0
my_grad[:,-1,1] = 0
my_grad[0,:,0] = 0
my_grad[-1,:,0] = 0
# Plot gradient
plot(gradu)
```

I’m using dolfin on colab, with the package that fem-on-colab provides.

PS : Not sure if I should ask another question in this topic, but I also encountered issues with the display of the solution that, when presented with a few nodes, gives only two colors in the output.

And sometimes it also leaves blank spaces (wasn’t able to reproduce though).