Lets say that one has a scalar field defined in 3D space for whose gradient he wants to find the Morse-Smale Complex for later performing an integration of several hexa-dimensional functions over attractor domains. The functions to be integrated are mostly concentrated around the attractor with relatively small contributions at the boundary. As well, the boundary can have a spiky shape. What approach (or algorithm or even software) would be appropriate to tackle this problem?
Currently the problem is divided in two steps: first the boundary of each attractor is found by testing in many points in space if following the gradient drives the search to other attractor or not. The initial points are chosen along rays starting at the attractor. Those rays point in the direction of the points generated by a Lebedev grid on a sphere of unit radius. So later the integration is separated in radial and angular variables. For each point in the radial quadrature (that goes no further than the point of the boundary found furthest from the attractor) the corresponding sphere is integrated checking for each angular point if it is inside or outside of the region of integration.
I am wondering if there is a more efficient approach around it that has enough resolution to treat the spiky boundaries and the maintains the benefits of the separation in angular and radial integration.