I'm trying to model diffusion through a cylindrical domain $D = \{ (x,y,z) : x^2 + y^2 \leq 1, \;\; 0 \leq z \leq 1\}$.
The is an initial concentration of the diffusant at the upper flat surface, say $C(x,y,z=1) = C_{init}$, and initially zero concentration within the domain.
The diffusant can move through both the $z=0$ and $z=1$ interfaces; there is no flux through the vertical walls.
The diffusion equation is $$ \frac{\partial C}{\partial t} = \Delta C. $$
The boundary conditions are the problem: $$ \frac{\partial C}{\partial t} (x, y, z=0, t) = K_0 \int_{z=0 \text{ cap}} \frac{\partial C}{\partial z} \;dS, $$ $$ \frac{\partial C}{\partial t} (x, y, z=1, t) = K_1 \int_{z=1 \text{ cap}} \frac{\partial C}{\partial z} \;dS. $$ These come from assuming that the diffusion through the end caps occurs into finite-volume chambers, and the accumulated diffusant affects the rate of diffusion across the boundary. The $K_i$ are constants.
Trying to set up the variational formulation for Fenics, with $C_N$ equal to the concentration at time step $N$, I get
$$ \int_D \frac{C_{N} - C_{N-1}}{\Delta t} v \; dV + \int_D \nabla C_N \cdot \nabla v \; d V - \int_{\partial D} v \nabla C_N\; d S = 0, $$
or equivalently $$ \int_D \frac{C_{N} - C_{N-1}}{\Delta t} v \; dV + \int_D \nabla C_N \cdot \nabla v \; d V - \int_{z=1 \text{ cap}} v \frac{\partial C_N}{\partial z}\; d S + \int_{z=0 \text{ cap}} v \frac{\partial C_N}{\partial z}\; d S= 0. $$
I don't know where to go from here. The boundary conditions tell about the integral of $\partial C / \partial z$ over the caps, but I don't know how to translate that into information about the integral against $v$ over the caps. Any advice would be welcome.
N.B.: I know that the problem I've stated is essentially one dimensional, and can be solved that way. The problem I'm actually interested in has a more complex domain, so I gave a simpler version that should capture the B.C. issue.