Solving a PDE on a submesh

I’m trying to solve a PDE system on a 2D domain. The PDE has the form [tex] \partial_tb + \Delta p=\lambda b [/tex] where p is a function that satisfies different equations on different subdomains:
\nabla\cdot\frac{1}{b}\nabla p = f(b)\quad\text{ when }b>b_p\
p=0\qquad \text{ when }b\leq b_p

At each time-step I want to find the region(s) where b>b_p (b_p is a constant, but the region defined by this threshold will change over time), and then solve the PDE for p on that subdomain/submesh.

Is there a way to do this kind of thing (i.e. solving a PDE on a subset of the domain) using NGSolve?

have you had a look at the NGSolve extension ngsxfem?

I think this is what you are looking for.


That looks like it may do exactly what I want.