Hello,
found the error: I mixed up 'region_hcs' and 'region_conductivities' in line 122.
Correct is:
Code:
m = BilinearForm(fes, symmetric=False)
m += SymbolicBFI(CoefficientFunction(region_densities * region_hcs) * u * v)
m.Assemble()
Now the solution fits also the analytic solution.