You would define a MeshFunction over your polygonal domains with for instance a single tag indicating the boundary between any polygon. Then, you should send this meshfunction into the dS measure (Interior facet mesh, as subdomain_data) and call the integral restriction.
Imagine you have your MeshFunction ff indicating the polygonal domains, then you would call jump(u)*dS(subdomain_data=ff, subdomain_id=id_of_polygonal_facets).
Note that if the order of “+” and minus matters, you would need to do some extra work, i.e.