Hi all, I am using dolfin 2019.1.0 installed via pip. In this particular example I’m solving the Poisson equation (as described in the classic example https://fenicsproject.org/olddocs/dolfin/1.5.0/python/demo/documented/poisson/python/documentation.html). I am trying to figure out the correct way to evaluate the residual (specifically the difference between the left hand side and the ride hand side of the equation that I am trying to solve). What I find is that if I don’t enforce the boundary conditions correctly, I get extremely high values for the residual.

Here is the minimal code:

```
import dolfin as d
import matplotlib.pyplot as plt
# Creating the mesh and function space:
mesh = d.UnitSquareMesh(32,32)
V = d.FunctionSpace(mesh, "Lagrange", 1)
# Defining the boundary conditions:
def boundary(x):
return x[0] < d.DOLFIN_EPS or x[0] > 1.0 - d.DOLFIN_EPS
u0 = d.Constant(0.0)
bc = d.DirichletBC(V, u0, boundary)
# Defining the Variational problem:
u_sol = d.TrialFunction(V)
v = d.TestFunction(V)
f = d.Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree = 1)
g = d.Expression("sin(5*x[0])", degree = 1)
a = d.inner(d.grad(u_sol), d.grad(v))*d.dx
L = f*v*d.dx + g*v*d.ds
u_sol = d.Function(V)
# Compute the solution:
d.solve(a == L, u_sol, bc)
# Evaluate the residual (i.e. the difference between the LHS and the RHS):
RES = d.inner(d.grad(u_sol), d.grad(v))*d.dx - f*v*d.dx - g*v*d.ds
res1 = d.assemble(RES)
res2 = d.assemble(RES)
# Enforcing the boundary conditions only for res1 and not for res2:
bc.apply(res1, u_sol.vector())
# The max and min values of the residual for the two cases differ greatly:
print("Residual max with BC applied:", res1.max())
print("Residual min with BC applied:", res1.min())
print("Residual max without BC applied:", res2.max())
print("Residual min without BC applied:", res2.min())
# Plotting the value of the residual of the entire domain:
resF1 = d.Function(V,res1)
plt.figure()
col = d.plot(resF1)
plt.colorbar(col)
plt.title('Residual with BC applied')
plt.show()
# The value of the residual with no boundary conditions applied:
resF2 = d.Function(V,res2)
plt.figure()
col = d.plot(resF2)
plt.colorbar(col)
plt.title('Residual without BC applied')
plt.show()
```

So my questions are as follows:

- What does bc.apply() do exactly?
- Am I calculating the residual (i.e. res1) correctly? And why is res2 so large?
- Is there an easier way to evaluate the residual at any given point of the domain?

Thanks!