I really don't like these magic numbers.  How about making it a tuple of names?

Also, I think it is wrong.  Imagine I have two boundary conditions on the top (each applying to half of the top).  Then I will end up apply bcs on the bottom as well (because the first time through assemble_bcs is 0, so I set it to 2, the next time it is 2 so I set it to 3).

I think you want:
````
bottom = any(bc.subdomain == "bottom" for bc in bcs)
top = any(bc.subdomain == "top" for bc in bcs)

assemble_bcs = (bottom, top)
````
then in Pyop2 assemble_bcs can be a tuple as well.