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:
[tex]
\nabla\cdot\frac{1}{b}\nabla p = f(b)\quad\text{ when }b>b_p\\
p=0\qquad \text{ when }b\leq b_p
[/tex]

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?

Hi
have you had a look at the NGSolve extension ngsxfem?
github.com/ngsxfem/ngsxfem
I think this is what you are looking for.

Best
Christopher

That looks like it may do exactly what I want.
Thanks!