Hey Shahan,
I think it should be possible by using the interpolation points (provided by the integration function). The code would be something like:
IntPoints = []
Integrate(levelset_domain = { "levelset" : lsetp1, "domain_type" : IF}, cf=0, mesh = mesh, ip_container = IntPoints, order=1)
Now, the list IntPoints contains all used integration points are located on the zero level. By increasing the order you can increase the number of returned points.
Best wishes,
Paul