I am trying to use the multi-union solid (g4 version 10.6.2, macOS), but I am finding some likely-buggy behaviour in the overlap checking functionality. Overlap checking seems to work fine for conventional Union (finishing in a very reasonable length of time), but hangs badly when using multi-unions (even with a completely equivalent solid), with occasional warnings of the following style:
Checking overlaps for volume R3_lv_PV (G4MultiUnion) ...
-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : GeomSolids1001
issued by : G4BooleanSolid::GetPointOnSurface()
Solid - Sub3_B49_e_R3_zone4
All 100k attempts to generate a point on the surface have failed!
The solid created may be an invalid Boolean construct!
*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------
Now to convince you that there is indeed a problem I attach two GDML files describing two equivalent geometries (with extension .txt as I am forbidden from uploading .gdml by this forum’s configuration) that demonstrate this difference in behaviour. If you do a diff of these two files you can clearly see that they are identical except for the use of multiunion in one and conventional binary union in the other (similarly if you visualise they are indeed same solid). Running overlap check on r3-union.txt terminates and passes quickly, and r3-multi-union.txt does not terminate, and occasionally emits the kind of warnings noted above. The diff:
So my question to the experts here is why is there a difference in the behaviour, where overlap checking with multi-union hanging indefinitely and binary union is OK for the attached file? Is this is a bug or is there some undocumented behaviour relating to multi-unions that affects the way the algorithm works? Naively this should not matter, if I understand how it works then it will just pick a node of the multi union at random and then pick a point on the surface of that node. So why the problem?
The warning message on a failure to generate a random point on the surface of a boolean solid may appear if the result of a boolean operation (subtraction or intersection) is a null object, or the size of the resulting solid is much smaller than the size of the constituent solids.
In your case the problem is not in the union but in the solid Sub3_B49_e_R3_zone4. I could not find the description of Sub3_B49_e_R3_zone4 in the attached files, but I guess it is a boolean subtraction where at least one of the constituent solids is much bigger than the resulting solid.
OK I have spent some more time trying to understand this issue and I think it’s clear to me now that this is not related to the primitives being too small. In fact the warning message that I posted is a red herring, I don’t think it has anything to do with the problem. In G4MultiUnion::GetPointOnSurface, where a random node of the multi-union is selected, solid.GetPointOnSurface(); has no problem generating a point on the surface of the node. The problem is that the points on the surfaces of these nodes are apparently never on the surface of the multi-union as a whole. To me this is geometrically, I want to say “impossible”, I cannot imagine how it is possible to create the union of nodes where none of the surfaces of the nodes correspond to the surface of the final multi-union.
So to me this suggests that the issue is in G4MultiUnion::InsideWithExclusion. It seems that in G4MultiUnion::InsideWithExclusion , at the line G4int limit = fVoxels.GetCandidatesVoxelArray(aPoint, candidates, exclusion);, fVoxels.GetCandidatesVoxelArray is always returning 0 for the multi-union solid I have created. So the next two loops are never run and the EInside is alwayskOutside. The reason for this is in G4Voxelizer:
G4int G4Voxelizer::GetCandidatesVoxelArray(const G4ThreeVector &point,
std::vector<G4int> &list, G4SurfBits *crossed) const
{
// Method returning the candidates corresponding to the passed point
list.clear();
for (auto i = 0; i <= 2; ++i)
{
if(point[i] < fBoundaries[i].front() || point[i] >= fBoundaries[i].back())
return 0; // THIS METHOD ALWAYS RETURNS HERE!!!!
}
...
So to be clear every single randomly generated point is always failing this first fBoundaries check in GetCandidatesVoxelArray (which I do not understand, sorry this code is hard for me to interpret as the datatypes involve are not abstracted). So what does this mean? Is this a bug in the multi-union implementation or is there something subtly wrong about my geometry such that GetCandidatesVoxelArray is always failing? Is there a known failure-case for this?
I am running these tests with the example in extended/persistency/gdml/G01.
The algorithm for selecting a random point on the surface of a Boolean solid is as follows: Geant4 generates a random point on the surface of one of the components of the Boolean operation. If the generated point is not on the surface of the resulting object, Geant4 retries. After 100,000 unsuccessful attempts, Geant4 gives up and simply returns the last generated point. Obviously, this point will not be on the surface of the resulting object, and, therefore, it will not be on the surface of other Boolean solids where the solid will be used.
So, the problem is in Sub3_B49_e_R3_zone4, which I believe is G4SubtractionSolid. I suggest you compare the volume of Sub3_B49_e_R3_zone4 and the volume of the solids that were used to construct it.
BTW, to use G4MultyUnion only for a couple of components is not very rational. It will not give any advantage in comparison with G4UnionSolid.
Yes I understand, and when it fails, it emits the warning I posted above. However there is a subtlety here in that the 100k limit is only for the binary Boolean solids—there is no 100k limit (and associated warning) for multi-union. What is happening is that it is getting stuck in the do/while loop of MultiUnion::GetPointOnSurface. I have added comments to the loop to try to explain:
// MultiUnion::GetPointOnSurface
do
{
G4long rnd = G4RandFlat::shootInt(G4long(0), size);
G4VSolid& solid = *fSolids[rnd];
point = solid.GetPointOnSurface(); // This part is fine, generally successfully generating points on the surfaces of a node.
const G4Transform3D& transform = fTransformObjs[rnd];
point = GetGlobalPoint(transform, point);
}
while (Inside(point) != EInside::kSurface); // This part fails repeatedly!
IT is the Inside(point) in G4MultiUnion::GetPointOnSurface that is failing repeatedly. The occasional warning of the type I posted above is just a rare occurrence that will occasionally happen as the above do/while loop is looping many times. That is why I said that it is a red herring, it is just a side-effect of the fact that the do/while loop is stuck in an infinite (or at least very very very long) loop.
To be clear: if the point = solid.GetPointOnSurface(); is successfully getting a point on the surface of a node, how can it be that it is never on the surface of the multi-union itself? I cannot imagine such a shape, and the example I uploaded is actually very simple geometrically anyway. If Geant4 can generate a point on the surface of a node then that point should also correspond to the surface of the multi-union itself a good proportion of the time unless the geometry is obscenely pathological.
Multi-union has one advantage even for only 2 components over binary union in that its nodes can be disconnected (whereas G4Union the two components must be overlapping somehow), admittedly in this instance it is not relevant as the two nodes are overlapping. Otherwise I understand yes it will offer no speedup and may even be slower, but to demonstrate the issue that I am facing, I have reduced the multi-union to feature fewer nods. In my real application this multi-union has 5 nodes, but I am simplifying the problem as much as I can. Perhaps I can go down to a single node and see if it persists. I can try this tomorrow.
What I am trying to say is that I think there’s either a bug in G4MultiUnion or G4Voxelizer or there are undocumented constraints (or at least ones that I have not read) about the type of CSG trees that multi-unions can be used with that I wish to know.
I am dependent on being able to make arbitrary multi-unions (in code that I have written elsewhere for automatic writing of geometry) and if Geant4 is unable to handle arbitrary multi-unions then I need to understand why, and if it turns out to be a bug in Geant4 then I would like to fix it. Simply using G4Union instead is not really an option for me due to reasons that are not important and are unrelated to this problem.
I’ve done some investigation and can confirm that it is a bug. G4MultiUnion object defined in r3-multi-union.txt is unable to generate a random point on its surface. Please create a report in the Geant4 Problem Tracking System:
As a temporary fix you can replace in G4MultiUnion::GetPointOnSurface() while (Inside(point) != EInside::kSurface);
with while (InsideNoVoxels(point) != EInside::kSurface);
I thank you for the workaround you have posted, this works well for me for now to progress in the mean time. What I actually did was change the G4MultiUnion::Inside to forward its call straight to InsideNoVoxels, as editing GetPointOnSurface will result in a bugged G4PVPlacement::CheckOverlaps as the if (solid->Inside(ptmp) != kSurface) will fail otherwise. As follows:
{
// Classify point location with respect to solid:
// o eInside - inside the solid
// o eSurface - close to surface within tolerance
// o eOutside - outside the solid
// Hitherto, it is considered that only parallelepipedic nodes can be
// added to the container
// Implementation using voxelisation techniques:
// ---------------------------------------------
// return InsideIterator(aPoint);
return InsideNoVoxels(aPoint);
// EInside location = InsideWithExclusion(aPoint);
// return location;
}
