# Neumann boundary conditions on arbitrary surface for finite difference diffusion

I am facing the following problem, formulated in practical terms: I have a region $$\Omega$$ in two or three dimensions, represented as a binary mask, and an initial density $$u_0$$ within that region that I'd like to evolve following a diffusion equation:

$$\partial_t u(\vec{x},t)=\nabla(D\nabla u)+f(u)$$

where D can be a quadratic matrix. I'm struggling with how to ensure Neumann boundary conditions of the form $$(D\nabla u)\cdot\vec{n}|_{\partial\Omega}=0$$, i.e. the direction of diffusion should be parallel to the surface along the entire surface. All I could find (and actually understand...) was in 1D, which I could translate to rectangular regions and scalar $$D$$, but now I'm trying to move to arbitrary regions. My first guess was to get the surface normals (e.g. from the gradient of a distance transform) and then project diffusion vectors along the surface onto a plane orthogonal to the surface normal, but I have no idea whether that's a good way to do it.

Are Neumann boundary conditions supposed to imply conservation of mass? Because I can't see how that's going to work in a discrete scheme. If it's in any way relevant, I'm using a simple forward model in time and central differences in space. Do I need an implicit time scheme instead? Any help is much appreciated!

• What do you mean by quadratic matrix? – nicoguaro May 22 '19 at 15:31
• What's your spacial discretization scheme? How you go about applying a Neumann BC is going to be specific to how you're representing your solution (finite volume, finite element, etc). – John May 22 '19 at 21:59
• Btw, a mixed formulation could be useful, since then your BC will be posed as Dirichlet condition for the heat flux and normal components are usually included into the dofs (e.g., for Raviart-Thomas finite elements) – VorKir May 23 '19 at 8:46
• Re matrix I only meant to say that D is not scalar. In fact, it can also be a function of $\vec{x}$. My spatial discretization scheme is finite differences, I guess... I'm given maps of D and the mask as pixel arrays (from MR images, because I want to model the diffusion of tumor cells in the brain) – Jens Petersen May 27 '19 at 9:07