Undesirable behavior of BRepAlgoAPI_Fuse

Hi Tiberiu Chelcea,

Im currently working on implementing a CAD kernel for a maxwell equation finite element solver. For this I need to be able to create models with overlapping regions that represent different materials. I had to learn that this is actually much more complicated to archieve with OCC than I thought. Therefore I was really happy to read all your comments in the OCC forum about this issue. It would be very nice if we could work together on finding a solution for this problem.
I adapted the code you posted (see class below, Im working with pythonocc for now, but will switch to C++ later on) and extended it to work with an arbitrary number of solids. If I check each new solid against all the solids that where added before (not creating a component of them) and use TopExp_Explorer instead of OuterShell to obtain the shell of a solid then it seems to work and shared faces dont get lost.
However, I've got some problems recovering the individual shells after the fuse operation with ShapeFix_Shell. For some reason, it sometimes returns 3 or 1 shells instead of the 2 I expect. For example, overlaying a Box with a Sphere will return 3 shells with shell number 2 being the shared face (which should be part of shell number 2!). Ovelaying the same Sphere with a Box works as expected. Next thing I will try is to re-implement the algorithm in ShapeFix_Shell to recover the shells more reliably. Other solution Im thinking about would be to glue the solids instead of fusing them, but how to find the edges/faces that need to be glued? And how to handle some special cases like two boxes on top of each other?
I found an interesting piece of C++ code in Gmsh (OCC_Connect.h/.OCC_Connect.cpp) which seems to implement an algorithm for "cut and merge" of shells. I didnt had the time yet to digg further into the code, but maybe it could be helpful for you (however, the author states in the Gmesh forum that the algorithm is quite slow and not very stable...)

class ShapeOverlay3D(object):
def __init__(self):
self.solids=[]

def add(self, OverlayShape):
if len(self.solids)==0:
self.solids.append(OverlayShape)
else:
for solidID, solid in enumerate(self.solids):
OverlayShape_modified, solid_modified, SolidsDidOverlap=self.OverlayTwoSolids(solid, OverlayShape)
if SolidsDidOverlap:
OverlayShape=OverlayShape_modified # the OverlayShape will contain additional edges/faces
self.solids[solidID]=solid_modified # update the shape that we did cut
self.solids.append(OverlayShape_modified)

def GetSolids(self):
return self.solids

def GetCompoundShape(self):
OverlayShape = TopoDS_Compound()
CmpBuilder = BRep_Builder()
CmpBuilder.MakeCompound(OverlayShape)
for aSolid in self.solids:
CmpBuilder.Add(OverlayShape, aSolid)
return OverlayShape

def GetShell(self, aSolid):
return TopExp_Explorer(aSolid, TopAbs_SHELL).Current()
#return TopoDS().Shell(BRepTools().OuterShell(TopoDS().Solid(aSolid)))

def MakeCutSolids(self, solid1, solid2):
Cutter = BRepAlgoAPI_Cut(solid1,solid2)
Cutter.Build()
aComp = Cutter.Shape()
# boolean operation in OCC always return a compound,
# which contains the solid, extract and return only the solid
return TopExp_Explorer(aComp, TopAbs_SOLID).Current(), Cutter.HasModified()

def MakeFuseShells(self, shell1, shell2):
fused_shells = BRepAlgoAPI_Fuse(shell1, shell2).Shape()
return fused_shells

def BuildSolidFromShell(self, aShell, fixSolid):
aSolid = BRepBuilderAPI_MakeSolid(aShell).Solid()
if fixSolid:
return aSolid
else:
FixSld = ShapeFix_Solid(aSolid);
FixSld.Perform()
aSolid=TopoDS().Solid(FixSld.Solid())
return aSolid

def RecoverIndividualShells(self, aCompound):
shells=[]
FixTopShell=ShapeFix_Shell(aCompound);
#FixTopShell.FixFaceOrientation(aCompound, False, True)
FixTopShell.Perform()
if (FixTopShell.NbShells() > 1):
shellComp = TopoDS().Compound(FixTopShell.Shape())
ExShls = TopExp_Explorer(shellComp, TopAbs_SHELL)
while ExShls.More():
shl1 = TopoDS().Shell(ExShls.Current())
# fix potential issues with the shell, e.g. it might
# contain faces that are not oriented properly
FixShl = ShapeFix_Shell(shl1); FixShl.Perform();
shells.append(FixShl.Shell())
ExShls.Next();
else: # only one shell recovered, probably the original solids didnt overlap
shells.append(FixTopShell.Shell())
return shells

def OverlayTwoSolids(self, CutSolid, OverlaySolid):
self.solidPairs_1.append(CutSolid)
self.solidPairs_2.append(OverlaySolid)
recovered_solids=[];
# cut solid with the shape we want to "overlay"
CutSolid_, SolidsOverlap = self.MakeCutSolids(CutSolid, OverlaySolid)
if not SolidsOverlap:
return CutSolid, OverlaySolid, SolidsOverlap
CutSolid_shell = self.GetShell(CutSolid_)
OverlaySolid_shell = self.GetShell(OverlaySolid)
# the fuse operation removes dublicate faces and returns a compound
# that contains all faces in a single shell
# the individual shells can be recovered with ShapeFix_Shell
# Note: in case the two shells did not overlap, the compound would
# contain two individual shells (but this case is handled above)
Combined_shell = TopoDS().Shell(TopExp_Explorer(fused_shells, TopAbs_SHELL).Current())
# recover the individual shells ...
recovered_shells = self.RecoverIndividualShells(Combined_shell)
# rebuild solids from shells ...
for aShell in recovered_shells:
aSolid = self.BuildSolidFromShell(aShell, True)
recovered_solids.append(aSolid)
if len(recovered_solids)==2:
return recovered_solids[0], recovered_solids[1], SolidsOverlap
elif len(recovered_solids)==3:
return recovered_solids[0], recovered_solids[2], SolidsOverlap
else:
raise ValueError, "Invalid number of solids recovered: " + str(len(recovered_solids))

Hi, Tiberiu Chelcea

I just read your discussion, and found it helps in terms of alternative ways to solve the boolean Unite (Fuse) operation problems. We know if we directly use BRepAlgoAPI_Fuse of two solids, excluding your 2 exceptions, it generates one Solid but seems the surfaces and edges are not united. And once you get this way, it's pretty hard to fuse the surfaces and edges. A simple example will be two identical cube-box aligned and touching at one face, after fuse, it's one Solid, but 10 surfaces instead of 6 surfaces for a large box. If my understanding is correct, you are using this shell fuse, then fix, then recreate solid method to avoid above situation to happen, thus will create a 6 surface box eventually, right?

For you exception number two, my thought is, it might be hard to distinguish if you need the cap of the solid2 be kept or not when they touch on their bottom surfaces, if you try to make your solid2 as solid2 = BRepPrimAPI_MakeBox(gp_Pnt(1, 1, -1), 1, 1, 1); so when you do cut, it sure will leave a hole on the bottom, and you then can make a solid with a hole. I haven't tried it myself, but thought it worth trying if we remove the ambiguity, see if it works.

Jane

Hi Tiberiu Chelcea,

I think I made some progress in solving the problem you describe. Im using the Partition_Spliter() algorithm in Salomes Geom module (note theres also a standalone version of the module available at http://sourceforge.net/projects/salomegeometry/, so you do not need to use the Salome package).

You can find a stripped-down version of the code I use below. In the attached PDF-file theres an image of an example with 3 spheres, 3 boxes and 1 cylinder that I combined such that shared faces wont get meshed twice and each solid can have its own material parameters. It also works for co-planar faces, e.g. for your example with the two boxes, which is really nice, cause it saves a lot of headaches trying to identify this case.

A problem I am having with this approach is that I can not reliably recover the individual solids after all faces have been splited by the partition algorithm (in the example five solids result for seven original solids that where combined). However, after I created a surface mesh with MEFISTO, it is sufficient for me to locate a point inside each of the original solids. Because this is the information that the volume mesher I use (tetgen) needs to identify the different material regions.

To find a point inside a shape I simply use a center point on one of its faces and move along the faces normal inwards by a certain distance. If it doesnt work for a face I just move to the next and try again. Using BRepClass3d_SolidClassifier I can check if a point is inside a shape.

I hope this helps also for your case. I would be happy to hear about your progress.

Cheers,
Mark

def PointInShape(aShape, aPoint, tol=1e-5):
solidClassifier = BRepClass3d_SolidClassifier(aShape, aPoint, tol)
return solidClassifier.State()==TopAbs_IN

def NormalOnFace(aFace):
DS=TopoDS()
bf = BRepGProp_Face(DS.Face(aFace))
bounds = bf.Bounds()
vec = gp_Vec()
pt = gp_Pnt()
u=min(bounds[0],bounds[1]) + abs(bounds[1]-bounds[0])/2.0
p=min(bounds[2],bounds[3]) + abs(bounds[3]-bounds[2])/2.0
pt=gp_Pnt(); vec=gp_Vec()
bf.Normal(u,p,pt,vec)
return vec, pt

def FindPointInShape(aShape, eps):
# an attempt to reliably identify a point that
# is located inside a given Shape, e.g. useful for tetgen region marker
faces = allFaces(aShape)
assert(len(faces)>0), "FindPointInShape: Invalid shape"
for face in faces:
# if aShape is a valid OCC Solid, then all of its face normals should
# point outward
N, P=NormalOnFace(face)
N.Normalize()
PInside = gp_Pnt(P.X()-eps*N.X(), P.Y()-eps*N.Y(), P.Z()-eps*N.Z())
if PointInShape(aShape, PInside, 1e-8):
return PInside
raise ValueError, "FindPointInShape: Cant find a point inside shape, try with different eps"

class ShapeOverlay3D(object):
def __init__(self, display):
self.Partitioner=Partition_Spliter()
self.solids=[]
self.origsolids=[]
self.display=display

def add(self, OverlayShape):
self.origsolids.append(OverlayShape)
if len(self.solids)==0:
self.solids.append(OverlayShape)
self.Partitioner.AddShape(OverlayShape)
else:
for solid in self.solids:
self.Partitioner.AddShape(solid)
self.Partitioner.AddShape(OverlayShape)
self.Partitioner.Compute()
self.Partitioner.RemoveShapesInside(OverlayShape)
result = self.Partitioner.Shape() # a compound with solids that share faces
self.solids=allSolids(result)
self.Partitioner.Clear()

def GetSolids(self):
return self.solids

def GetCenterPoints(self):
centerpoints=[]
for solid in self.origsolids:
centerpoints.append(FindPointInShape(solid, 3.0))
return centerpoints

def GetCompoundShape(self):
resultShape = TopoDS_Compound()
CmpBuilder = BRep_Builder()
CmpBuilder.MakeCompound(resultShape)
for solid in self.solids:
CmpBuilder.Add(resultShape, solid)
return resultShape