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).