I try to write a heat equation solver with Neumann boundary condition. There is one on the tutorial A known analytical solution — FEniCSx tutorial (jsdokken.com) but I do not know how to change it into a Neumann boundary condition.

I also try the tutorial code Setting multiple Dirichlet, Neumann, and Robin conditions — FEniCSx tutorial (jsdokken.com) and rewrite it into a heat solver, but it fails , lack of accuracy and speed.

I hope that someone can tell me how to rewrite the first code into a pure zero Neumann boundary condition.

For the second case on tutorial, the following code is the one that I rewrote.

for i in range(num_steps):

u_exact.t += dt

kappa_exact.t += dt

source.t += dt

u_ex.interpolate(u_exact.eval)

```
f=Function(V)
f.interpolate(source.eval)
kappa.interpolate(kappa_exact.eval)
u, v = TrialFunction(V), TestFunction(V)
F = dt *kappa * inner(grad(u), grad(v)) * dx - inner(u_n + dt *f, v) * dx+inner(u, v) * dx
bcs = []
for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F += condition.bc
```

# Solve linear variational problem

```
a = lhs(F)
L = rhs(F)
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
u_n = uh
```

A zero Neumann condition means that there is no additional integral added to your variational form. This follows from deriving the weak form of the problem, where one after integration by parts will get a term like \int_{\partial\Omega} \kappa \frac{\partial u}{\partial n}\cdot v ~\mathrm{d} s which one explicitly removes from the varitional form.