# Undesirable behavior of BRepAlgoAPI_Fuse

I'm working on a problem of identifying common faces between solids and re-building solids to share these common faces. Here's the algorithm used (which works very well for the vast majority of cases):

TopoDS_Shell shp1_shell = BRepTools::OuterShell(solid1);

TopoDS_Shell shp2_shell = BRepTools::OuterShell(solid2);

fused_shells = BRepAlgoAPI_Fuse(shp1_shell, shp2_shell).Shape();

TopExp_Explorer Ex;

for (Ex.Init(fused_shells, TopAbs_ShapeEnum.TopAbs_SHELL); Ex.More(); Ex.Next()) {

TopoDS_Shell crt_shell = TopoDS.Shell(Ex.Current());

ShapeFix_Shell FixTopShell(crt_shell);

FixTopShell.Perform();

if (FixTopShell.NbShells() > 1) {

TopoDS_Compound shellComp = OCTopoDS.Compound(FixTopShell.Shape());

TopExp_Explorer ExShls(shellComp, TopAbs_ShapeEnum.TopAbs_SHELL);

for (; ExShls.More(); ExShls.Next()) {

TopoDS_Shell shl1 = TopoDS.Shell(ExShls.Current());

ShapeFix_Shell FixShl = new ShapeFix_Shell(shl1);

FixShl.Perform();

TopoDS_Solid sol_tmp1 = BRepBuilderAPI_MakeSolid(FixShl.Shell()).Solid();

ShapeFix_Solid FixSld(sol_tmp1);

FixSld.Perform();

TopoDS_Solid sol_tmp = TopoDS.Solid(FixSld.Solid());

// add it to collection....

}

} else {

TopoDS_Shell aShell = FixTopShell.Shell();

TopoDS_Solid sol = new BRepBuilderAPI_MakeSolid(aShell).Solid();

// add it to collection....

}

}

Pretty much, the algorith fuses the shells of the two components, then looks at the result and re-builds solids from these resulting shells (after fixing them). This algorithm works very well with two exceptions

1. If the two solids touch each other only on a point or a line, the resulting solids might not be closed (this does not bother me -- yet, at least)

2. If we have a solid that's almost inside other solid, the result is undesirable. For this example:

solid1 = BRepPrimAPI_MakeBox(gp_Pnt(0, 0, 0), 3, 3, 3);

solid2 = BRepPrimAPI_MakeBox(gp_Pnt(1, 1, 0), 1, 1, 1);

solid1 = .... cut solid2 from solid1 ....

Now, when applying the algorithm above, the re-built solid1 has 7 faces (the faces that don't touch solid2 + the bottom ones, one corresponding to the bottom face of solid2 (extending from (1,1,0) to (2,2,0), and the other the surrounding face (extending from (0,0,0) to (3,3,0) without the hole). Therefore, when the solid is re-build, it gets again to the original box. Solid2 is recreated correctly.

What I'd like is to have the BRrepAlgoAPI_Fuse create the shell for solid1 going around the quasi-hole at the bottom (where solid2 is), and not using the bottom face of solid2. Is there some way to get this behavior? Maybe the algorithm should be changed to make it handle this case (though it's not clear how to identify the case)?

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 Jane,

yes, the algorithm I'm employing fuses the shells of the solids, and not the solids themselves. This avoids creating a single solid out of the two initial boxes. The resulting compound of shells is then processed and the solids recreated one by one.

I've done a bit more work on fixing this issue, it's discussed here: http://www.opencascade.org/org/forum/thread_19053/ Unfortunatelly, it seems to fail for cylindrical faces, I'm not sure why; it works for the case of a box contained in another box.

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

Hi Mark,

thanks for posting the code -- I'll have to try it out (unfortunately, now I have to work on another project, so it will take a bit before I can do it -- and that's why it took me a while to answer this thread).

However, I came with a solution to this issue, you might want to check it at http://www.opencascade.org/org/forum/thread_19053/ . It works for this case where there's a box inside another box, but it fails for a cylinder inside another cylinder, I'm not sure why. There might be something related to the orientation of faces that makes BRepTools_ReShape mis-handle the common cylindrical face.

To "recover" the newly created solids (i.e. to figure out which one of the original ones corresponds to a new solid in order to be able to set the materials), I'm storing the properties of the solids obtained from GProp_GProps: as far as I could tell, it's sufficient to store the center of mass (from GProp_GProps::CenterOfMass()) and the "mass" (GProp_GProps::Mass()). It might be a bit faster than trying to find points inside solids.

Hi Mark,

I've tried to apply the Salome Partition algorithm, but i'm still having problems with the case I've mentioned (a small cylinder inside a larger one, both aligned at the bottom face). The result of Partition_Splitter still has the extra round face that should be shared by the two cylinders (the outer face of the small cylinder) duplicated -- and thus the mesher fails (I'm using Gmsh).

One thing I don't understand in your algorithm is the use of 'OverlayShape'. What is that shape? Could you post a full example, maybe for the case of the two cylinders.

## Mark, Thank you very much!

Mark, Thank you very much!