How to interpolate solution on the mesh _including_ boundary?

[przemarbor]

Hi,
The problem seems trivial but I cannot find a way of solving it
basing only on package documentation.
I want to solve 2d BVP and generate 3d surface plot of the solution.
Let it be given simple Laplace problem on a 2D grid [0,1]x[0,1]:

nabla^2 u = 0

with BCs as follows:

u = 0, for x = 0
u = 1, for x = 1
du/dn = 0, for y = 0 and y = 1

with exact solution u(x,y) = x.

There are some problems I stumbled upon:

1. It seems that there is no method built-in py-pde package to generate
   3D plots...
   OK, that's no problem - I can always use matplotlib's functions.
   In order to do that I need to interpolate obtained solution on a
   new grid but...
2. ... now, how can I interpolate the field basing on
   solution and_including the values on the boundary?
   When I print res.data it gives me field values only for internal points.
   The same happens when I try to interpolate - res.interpolate((0,0))
   returns 0.099999999, not 0, or when the field is plotted with
   built-in plot() function - the plot is from 0.1 to 0.9 while it is expected to be from 0 to 1.

Here's my code:

import numpy as np
from pde import CartesianGrid, solve_laplace_equation

xmin = 0
xmax = 1
ymin = 0
ymax = 1
xgridNo = 5
ygridNo = 5

grid =  CartesianGrid([[xmin, xmax], [ymin, ymax]], [xgridNo, ygridNo])

bc_x = [{"value": 0}, {"value": "x"}]
bc_y = [{"derivative": 0}, {"derivative": 0}]
bcs =[bc_x, bc_y]

res = solve_laplace_equation(grid, bcs)
res.plot()

x = np.linspace(xmin,xmax)
y = np.linspace(ymin,ymax)
X, Y = np.meshgrid(x,y)

Z = res.interpolate((0,0))
print(Z)
print(res.data)

And it's output:

0.09999999999999998
[[0.1 0.1 0.1 0.1 0.1]
 [0.3 0.3 0.3 0.3 0.3]
 [0.5 0.5 0.5 0.5 0.5]
 [0.7 0.7 0.7 0.7 0.7]
 [0.9 0.9 0.9 0.9 0.9]]

[david-zwicker]

This is expected behavior since py-pde keeps track of the field values at the center of the cells; see documentation. The resulting field has no information about the boundary conditions and uses nearest-neighbor interpolation at the border to not distort anything (it uses linear interpolation in the inside of the domain).

To obtain the values of the field on the boundary, you need to supply the boundary conditions to the interpolation function:

res.interpolate((0, 0.5), bc=bcs)

which returns 0, as expected. Note that the interpolation has trouble dealing with the corners, so that res.interpolate((0, 0), bc=bcs) can return non-sensical values. #547 should improve the situation and will hopefully arrive in the main branch soon!