PythonOCC and SMESH

PythonOCC and SMESH

The SMESH module that is included in the PythonOCC distribution has a lot to offer in terms of advanced CAD meshing techniques and therefore deserves some special attention. With this blog post I will try to shed some light on  SMESH’s capabilities and the merits it has to offer for the PythonOCC CAD platform. I will mainly focus on prism meshing techniques,  simply because this is my current field of interest. But also because it is a fruitful topic for demonstrating  SMESH’s capabilities.

PythonOCC’s SMESH is a standalone version (created by Fotios Sioutis) of the SMESH module that comes shipped with the  Salome Platform. Currently PythonOCC is being updated to the latest versions of the SMESH and GEOM modules (Version 6.5.0) and the latest Opencascade Version 6.5.3 in it’s  Opencascade Community Edition version 0.10 incarnation.

Together with these updates fully functional Netgen (version 4.9.13) and Tetgen (version 1.4.3) meshing plugins will be available to enjoy high quality unstructured grid generation of complex CAD models directly from within PythonOCC.

SMESH mesh framework: batteries included

SMESH is a complete meshing framework for the discretization of Opencascade CAD models that comes shipped with a considerable number of ready-to-use mesh algorithms and mesh hypotheses (the latter are used for fine-grained mesh size and quality control).  Here’s an overview of the most important algorithms and hypotheses:

SMESH Algorithms and Hypothesis

Most important SMESH Mesh algorithms and hypothesis. (*) New in upcoming PythonOCC release; ($$) Plugins that require commercial software.

A (*) indicates algorithms that will be newly included in the upcoming PythonOCC version and ($$) indicates mesh plugins that rely on commercial software (the GHS3D tetrahedral mesher and the BLSurf quad/tri surface mesher) and therefore are not included in PythonOCC.

The Netgen mesher (developed by Joachim Schöberl) provides surface triangulations and volume discretisations with usually high mesh qualities thanks to it’s automatic mesh optimization techniques. Netgen surface grids are usually of significantly better quality  compared to the results that can be obtained with Mefisto2. In some situations, though, Mefisto2 seems to be more stable due to the simplicity of it’s mesh algorithm. In several occasion I saw Netgen version 4.5 fail on spherical surfaces (they are tricky to mesh because of poles in their surface parametrization). However, this seems to have improved with Netgen version 4.9.13. Of course, because SMESH offers the flexibility to arbitrarily mix mesh algorithms , Netgen and Mefisto2 can suitably be used together for discretizing complicated models.

The new Tetgen plugin that will be available in PythonOCC  is capable of performing tetrahedral volume discretizations with a given lower bound on element qualities.  This threshold, that can be defined in addition to a maximum element volume constraint, limits the aspect ratio of all tetrahedra in the final mesh. The quality constraint sets a limit on the ratio q=R/L, where R is the radius of the circumsphere of a tetrahedron and L is the tetrahedrons shortest edge length. An upper bound on q is usually required for FE methods to achieve high solution accuracies. Also q can be used to control the growth rate of tetrahedra e.g. with increasing distance away from a densely discretized surface. Of course within a compound shape different quality constraints can be set on different regions for accurate control of mesh densities and mesh qualities.

SMESH can easily be extended with plugins to provide mesh functionalities that are neither available through any of SMESH’s plugins nor through any of SMESH’s standard mesher (StdMesher* classes). These plugins can be implementing by

  1. deriving from the class SMESH_3D_Algo, (or _2D_ / _1D_) for developing the mesh algorithm,
  2. implementing mesh hypotheses for providing necessary constraints/parameters to the new algorithm (in case none of the existing hypotheses fit) deriving from the SMESH_Hypothesis class and finally
  3. creating python wrappers for the new SMESH plugin.

Besides it’s meshing functionalities SMESH provides classes for the import/export of meshes (currently supported: DAT, STL, UNV file formats), a large set of mesh modification techniques (like mesh smoothing, element merging, splitting of elements, etc.), extensive mesh quality controls  and  techniques for filtering and selecting elements. Documentations for SMESH can be found on the Salome homepage (here and here). Sometimes it is quite instructive to see SMESH in action by experimenting with the meshing functionalities available in the Salome project (if you experience issues installing Salome on your computer, you might want to give CAELinux a try).

Viscous Layers: Prismatic boundary layers

One thing worth pointing out about SMESH version 6.5.0 is the newly available “Viscous Layers” plugin for creating 3D prismatic boundary layers. It allows for the automatic construction of prismatic boundary layers that grow inwards from the bounding faces of a solid shape. The algorithm requires the number of desired layers, a stretching factor to define  cell layer thicknesses and a total thickness value to be specified. Optionally faces can be specified that should not be treated. Though it does sometimes fail on solids with strongly bended concave faces, it seems to work well in many cases. Here’s an example for a cylinder from which a box with filleted edges was carved out:

New "Viscous Layers" algorithm in SMESH version 6.5.0

Example mesh created using the new “Viscous Layers” algorithm (in conjunction with Netgen) available in SMESH version 6.5.0.

For finite element (FE) computations this feature is very useful for refining boundaries between solids that represent strong material contrasts.

Currently the “Viscous Layers” algorithm works together only with the 3D Netgen mesher, i.e. it has not been integrated with the Tetgen mesher yet. The reason for this is a rather technical. The class StdMeshers_ViscousLayers in SMESH is implement as a hypothesis (i.e. it derives from SMESH_Hypothesis) that can be added to a NETGEN_3D algorithm in addition to the StdMeshers_MaxElementVolume Hypothesis. In the Compute()  method of the Netgen plugin class (file NETGENPlugin_NETGEN_3D.cpp) the surface mesh of the solid to be meshed is handed over to the Netgen mesher for performing the volume meshing. To facilitate the boundary layer meshing a proxy surface mesh (class SMESH_ProxyMesh ) is used:

1
2
3
4
5
6
7
SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
    if ( _viscousLayersHyp )
    {
      proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape );
      if ( !proxyMesh )
        return false;
    }
SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
    if ( _viscousLayersHyp )
    {
      proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape );
      if ( !proxyMesh )
        return false;
    }

In case boundary layer meshing should not be performed, SMESH_ProxyMesh is initialized with the input mesh and does nothing else but simply pass the required surface mesh through. Otherwise the proxy mesh is computed by the “Viscous Layers” algorithm: It computes all prismatic boundary layers and then returns the internal triangulated surface of those boundary layers. This new triangulated surface is subsequently used to compute the tetrahedral mesh of the remaining volume.

The problem with integrating the “Viscous Layers” algorithm with Tetgen is that the Tetgen mesher requires region and hole markers to be specified in conjunction with the discretized surfaces that define the geometry to be meshed. These region markers need to be specified in the form of 3D points that unambiguously identify hole and solid regions in the input geometry. Currently I use an algorithm for automatically finding those points that considers the original CAD geometry. Once this algorithm is modified to work on the discretized geometry instead, integrating the “Viscous Layers” algorithm into the Tetgen mesher plugin should be pretty simple.

Defining algorithms and hypotheses

SMESH provides a very general way of defining how each part of a CAD geometry should be discretized. This is done by adding mesh algorithms and hypotheses on subparts of a geometry (solids, faces, edges …) as necessary. Thereby algorithms and hypotheses automatically apply to all descendants of the shape they are applied on. For example, if we set a 1D mesh constraint on a solid shape then all edges of the solid will be discretized respecting that 1D constraint.

It often makes sense to apply 1D/2D constraints on the main shape to be meshed (thereby making those constraints global) and overwrite them locally on edges and faces that need a different treatment. The same is true for algorithms, of course: It might make sense to let Netgen mesh some solids of a compound geometry (local algorithms) and let Tetgen do the job on the remaining solids (global algorithm; this is just an example of course).

Let’s take a look at a simple yet somewhat more advanced example: Meshing a box with prismatic cells:

Meshing a box with prism cells

2D/3D mesh algorithms and hypotheses required for meshing a box with prism cells.

Here we apply a Netgen_2D algorithm on the bottom face of the box and a Projection_2D algorithm together with a ProjectionSource2D hypothesis (specifying the bottom face as the projection source) on it’s top face. All the remaining “wall faces” have Quadrangle_2D algorithms applied on. Note that in this example 1D algorithms/hypothesis are not shown.

To define the number of prism layers here we have defined Regular_1D algorithms in conjunction with NumberOfSegments hypotheses on all edges connecting the bottom with the top face. Of course, we could also have used Arithmetic_1D, StartEndLength or FixedPoints1D hypotheses instead for irregularly spaced prism layers. Note that because SMESH’s Projection_2D algorithm projects a given surface discretisation into the parametric space of the target face, we could also have tilted the boxes’ “wall faces” inwards (using BRepOffsetAPI_DraftAngle) prior to meshing it – the prism meshing would have worked nonetheless.

Here’s another example: Meshing a thick spherical shell with prisms.

Another prism mesh example

SMESH algorithms and hypotheses used to mesh a thick spherical shell with radial prisms.

The benefits of the surface mesh projection technique becomes more evident in this case (the inner triangular mesh created with MEFISTO_2D is projected into the parameter space of the outer spherical face). Note that in this case a RadialQuadrangle_1D2D algorithm is applied to the face that connects the two spherical faces. It takes care of projecting the 1D discretization of the faces inner edge to it’s outer edge and the subsequent discretization of the face with quadrilaterals. Because in this example no edge exists that connects the faces inner edge to it’s outer edge, we need to apply a NumberofLayers2D hypothesis instead of defining the quadrangle layer distribution by specifying 1D edge constraints.

The following example illustrates how to create a simple hybrid (with mixed structured/unstructured cells) mesh:

Hybrid structured/unstructured mesh example

Example on how to create a hybrid structured/unstructured mesh with SMESH.

Such a  geometry can be created using GEOM’s Splitter algorithm (I’ll talk a little bit more about that later on). The important question here is: After we have created our final geometry, how do we find the right solids/faces to apply the mesh algorithms/hypotheses on ? It’s not a simple question because most of Opencascade’s (or GEOM’s) BRep algorithms recreate at least part of the geometry (as they are supposed to) and thereby invalidate the shape or subshape references we might have had before. In this case, we could search for those shapes geometrically:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def HybridMeshExample():
    Box = BRepPrimAPI_MakeBox(gp_Pnt(0,0,0), 100,100,70).Shape()
    Sphere = BRepPrimAPI_MakeSphere(gp_Ax2(gp_Pnt(50,50,35), gp_Dir(0,-1,0)), 30).Shape()
    Box2 = BRepPrimAPI_MakeBox(gp_Pnt(0,0,-10), 100,100,10).Shape()
    # create compound solid with shared faces
    Splitter = GEOMAlgo_Splitter()
    Splitter.AddShape(Box)
    Splitter.AddShape(Box2)
    Splitter.AddShape(Sphere)
    Splitter.Perform()
    result=Splitter.Shape()
    # Locate the bottom face and the face shared between the two boxes
    finder=ShapeFinder()
    faces=finder.FindShapesOnBoxFace(result, gp_Pnt(0,0,-10), gp_Pnt(100,100,100), "Bottom", "ON", TopAbs_FACE)
    BottomFace=TopTools_ListIteratorOfListOfShape(faces)
    faces=finder.FindShapesOnBoxFace(result, gp_Pnt(0,0,0), gp_Pnt(100,100,70), "Bottom", "ON", TopAbs_FACE)
    MiddleFace=TopTools_ListIteratorOfListOfShape(faces[0]).Value()
    # Locate the individual solids
    solids=finder.FindShapesWithPoint(result, TopAbs_SOLID, gp_Pnt(50,50,-5), TopAbs_IN)
    BottomSolid=solids.Value(1)
    solids=finder.FindShapesWithPoint(result, TopAbs_SOLID, gp_Pnt(50,50,3.), TopAbs_IN)
    TopSolid=solids.Value(1)
    solids=finder.FindShapesWithPoint(result, TopAbs_SOLID, gp_Pnt(50,50,35.), TopAbs_IN)
    SphereSolid=solids.Value(1)
    # Create a hybrid mesh
    meshgen = MT.MeshGen()
    surfmesh = meshgen.Mesh(result, "Mesh")
    # Regular 1D discretization of all edges by default
    surfmesh.Segment(MT.REGULAR, result).AutomaticLength(0.1)
    surfmesh.Triangle(MT.MEFISTO, TopSolid); surfmesh.Triangle(MT.MEFISTO, SphereSolid)
    # Triangulation of BottomFace; projection BottomFace -> MiddleFace
    proj2D = surfmesh.Projection2D( MiddleFace )
    proj2D.SourceFace( BottomFace )
    surfmesh.Segment(MT.REGULAR, BottomFace).AutomaticLength(0.1)
    surfmesh.Triangle(MT.MEFISTO, BottomFace)
    # Regular 1D discretization on faces connecting top,bottom faces of BottomSolid
    surfmesh.Segment(MT.REGULAR, BottomSolid).NumberOfSegments(2)
    # Quadrangle algorithm on all faces of bottom solid not yet treated
    surfmesh.Quadrangle(BottomSolid)
    surfmesh.Prism(BottomSolid) # 3D extrusion algorithm on BottomSolid
    surfmesh.Compute()
    surfmesh.Display(display)
def HybridMeshExample():
    Box = BRepPrimAPI_MakeBox(gp_Pnt(0,0,0), 100,100,70).Shape()
    Sphere = BRepPrimAPI_MakeSphere(gp_Ax2(gp_Pnt(50,50,35), gp_Dir(0,-1,0)), 30).Shape()
    Box2 = BRepPrimAPI_MakeBox(gp_Pnt(0,0,-10), 100,100,10).Shape()
    # create compound solid with shared faces
    Splitter = GEOMAlgo_Splitter()
    Splitter.AddShape(Box)
    Splitter.AddShape(Box2)
    Splitter.AddShape(Sphere)
    Splitter.Perform()
    result=Splitter.Shape()
    # Locate the bottom face and the face shared between the two boxes
    finder=ShapeFinder()
    faces=finder.FindShapesOnBoxFace(result, gp_Pnt(0,0,-10), gp_Pnt(100,100,100), "Bottom", "ON", TopAbs_FACE)
    BottomFace=TopTools_ListIteratorOfListOfShape(faces)
    faces=finder.FindShapesOnBoxFace(result, gp_Pnt(0,0,0), gp_Pnt(100,100,70), "Bottom", "ON", TopAbs_FACE)
    MiddleFace=TopTools_ListIteratorOfListOfShape(faces[0]).Value()
    # Locate the individual solids
    solids=finder.FindShapesWithPoint(result, TopAbs_SOLID, gp_Pnt(50,50,-5), TopAbs_IN)
    BottomSolid=solids.Value(1)
    solids=finder.FindShapesWithPoint(result, TopAbs_SOLID, gp_Pnt(50,50,3.), TopAbs_IN)
    TopSolid=solids.Value(1)
    solids=finder.FindShapesWithPoint(result, TopAbs_SOLID, gp_Pnt(50,50,35.), TopAbs_IN)
    SphereSolid=solids.Value(1)
    # Create a hybrid mesh
    meshgen = MT.MeshGen()
    surfmesh = meshgen.Mesh(result, "Mesh")
    # Regular 1D discretization of all edges by default
    surfmesh.Segment(MT.REGULAR, result).AutomaticLength(0.1)
    surfmesh.Triangle(MT.MEFISTO, TopSolid); surfmesh.Triangle(MT.MEFISTO, SphereSolid)
    # Triangulation of BottomFace; projection BottomFace -> MiddleFace
    proj2D = surfmesh.Projection2D( MiddleFace )
    proj2D.SourceFace( BottomFace )
    surfmesh.Segment(MT.REGULAR, BottomFace).AutomaticLength(0.1)
    surfmesh.Triangle(MT.MEFISTO, BottomFace)
    # Regular 1D discretization on faces connecting top,bottom faces of BottomSolid
    surfmesh.Segment(MT.REGULAR, BottomSolid).NumberOfSegments(2)
    # Quadrangle algorithm on all faces of bottom solid not yet treated
    surfmesh.Quadrangle(BottomSolid)
    surfmesh.Prism(BottomSolid) # 3D extrusion algorithm on BottomSolid
    surfmesh.Compute()
    surfmesh.Display(display)

Here the class ShapeFinder, that simply employs BRepClass3d_SolidClassifier and GEOMAlgo_FinderShapeOn1, could be implemented like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
class ShapeFinder(object):
    def __init__(self):
       pass
 
    def FindShapesWithPoint(self, aShape, SubShapeType, P, location, aTol=0.0001):
        # sensible combination of values for SubShapeType and location are:
        # SubShapeType = TopAbs_EDGE:  --> location = TopAbs_ON
        #              = TopAbs_FACE:  -->          = TopAbs_ON
        #              = TopAbs_SHELL: -->          = TopAbs_IN, TopAbs_OUT or TopAbs_ON
        #              = TopAbs_SOLID: -->          = TopAbs_IN or TopAbs_OUT
        if not (location in [TopAbs_IN, TopAbs_ON, TopAbs_OUT, TopAbs_INTERNAL]):
            raise ValueError, "Unexpected Value for SubShapeType"
        Ex=TopExp_Explorer(aShape, SubShapeType); DS=TopoDS.TopoDS()
        result = TopTools_SequenceOfShape()
        while Ex.More():
            subshape=Ex.Current()
            B=BRepClass3d_SolidClassifier(subshape, P, aTol)
            if B.State()==location:
                result.Append(subshape)
            Ex.Next()
        return result
 
    def FindShapesOnBoxFace(self, aShape, P0, P1, BoxSide, SearchLocation, SearchForShapesWithType, aTol = 0.0001):
        BoxBuilder=BRepPrimAPI_MakeBox(P0, P1)
        if BoxSide=="Bottom":
            BoxFace=BoxBuilder.BottomFace()
        elif BoxSide=="Back":
            BoxFace=BoxBuilder.BackFace()
        elif BoxSide=="Front":
            BoxFace=BoxBuilder.FrontFace()
        elif BoxSide=="Left":
            BoxFace=BoxBuilder.LeftFace()
        elif BoxSide=="Right":
            BoxFace=BoxBuilder.RightFace()
        elif BoxSide=="Top":
            BoxFace=BoxBuilder.TopFace()
        else:
            raise ValueError, "BoxSide must be Bottom, Back, Front, Left, Right or Top"
        surface=BRep_Tool().Surface(BoxFace) # a Handle_Geom_Surface
        return self.FindShapesOnSurface(aShape, surface, SearchLocation, SearchForShapesWithType, aTol)
 
    def FindShapesOnSurface(self, aShape, SearchSurface, SearchLocation, SearchForShapesWithType, aTol = 0.0001):
        if not (isinstance(SearchSurface, Handle_Geom_Surface)):
            raise TypeError, "SearchSurface must be an instance of Handle_Geom_Surface"
        if not (SearchForShapesWithType in [TopAbs_VERTEX, TopAbs_EDGE, TopAbs_WIRE, TopAbs_FACE, TopAbs_SHELL, TopAbs_SOLID]):
            raise ValueError, "SearchShapeType must be an instance of TopAbs_ShapeEnum"
        location = {TopAbs_IN:GEOMAlgo_ST_IN, TopAbs_ON:GEOMAlgo_ST_ON, TopAbs_ON:GEOMAlgo_ST_OUT}[SearchLocation]
        aFinder=GEOMAlgo_FinderShapeOn1()
        aFinder.SetShape(aShape)
        aFinder.SetTolerance(aTol)
        aFinder.SetSurface(SearchSurface)
        aFinder.SetShapeType( SearchForShapesWithType )
        aFinder.SetState(location)
        aFinder.SetNbPntsMin(3); aFinder.SetNbPntsMax(100);
        aFinder.Perform()
        if aFinder.ErrorStatus():
            return None
        else:
            return aFinder.Shapes() # returns TopTools_ListOfShape
class ShapeFinder(object):
    def __init__(self):
       pass

    def FindShapesWithPoint(self, aShape, SubShapeType, P, location, aTol=0.0001):
        # sensible combination of values for SubShapeType and location are:
        # SubShapeType = TopAbs_EDGE:  --> location = TopAbs_ON
        #              = TopAbs_FACE:  -->          = TopAbs_ON
        #              = TopAbs_SHELL: -->          = TopAbs_IN, TopAbs_OUT or TopAbs_ON
        #              = TopAbs_SOLID: -->          = TopAbs_IN or TopAbs_OUT
        if not (location in [TopAbs_IN, TopAbs_ON, TopAbs_OUT, TopAbs_INTERNAL]):
            raise ValueError, "Unexpected Value for SubShapeType"
        Ex=TopExp_Explorer(aShape, SubShapeType); DS=TopoDS.TopoDS()
        result = TopTools_SequenceOfShape()
        while Ex.More():
            subshape=Ex.Current()
            B=BRepClass3d_SolidClassifier(subshape, P, aTol)
            if B.State()==location:
                result.Append(subshape)
            Ex.Next()
        return result

    def FindShapesOnBoxFace(self, aShape, P0, P1, BoxSide, SearchLocation, SearchForShapesWithType, aTol = 0.0001):
        BoxBuilder=BRepPrimAPI_MakeBox(P0, P1)
        if BoxSide=="Bottom":
            BoxFace=BoxBuilder.BottomFace()
        elif BoxSide=="Back":
            BoxFace=BoxBuilder.BackFace()
        elif BoxSide=="Front":
            BoxFace=BoxBuilder.FrontFace()
        elif BoxSide=="Left":
            BoxFace=BoxBuilder.LeftFace()
        elif BoxSide=="Right":
            BoxFace=BoxBuilder.RightFace()
        elif BoxSide=="Top":
            BoxFace=BoxBuilder.TopFace()
        else:
            raise ValueError, "BoxSide must be Bottom, Back, Front, Left, Right or Top"
        surface=BRep_Tool().Surface(BoxFace) # a Handle_Geom_Surface
        return self.FindShapesOnSurface(aShape, surface, SearchLocation, SearchForShapesWithType, aTol)

    def FindShapesOnSurface(self, aShape, SearchSurface, SearchLocation, SearchForShapesWithType, aTol = 0.0001):
        if not (isinstance(SearchSurface, Handle_Geom_Surface)):
            raise TypeError, "SearchSurface must be an instance of Handle_Geom_Surface"
        if not (SearchForShapesWithType in [TopAbs_VERTEX, TopAbs_EDGE, TopAbs_WIRE, TopAbs_FACE, TopAbs_SHELL, TopAbs_SOLID]):
            raise ValueError, "SearchShapeType must be an instance of TopAbs_ShapeEnum"
        location = {TopAbs_IN:GEOMAlgo_ST_IN, TopAbs_ON:GEOMAlgo_ST_ON, TopAbs_ON:GEOMAlgo_ST_OUT}[SearchLocation]
        aFinder=GEOMAlgo_FinderShapeOn1()
        aFinder.SetShape(aShape)
        aFinder.SetTolerance(aTol)
        aFinder.SetSurface(SearchSurface)
        aFinder.SetShapeType( SearchForShapesWithType )
        aFinder.SetState(location)
        aFinder.SetNbPntsMin(3); aFinder.SetNbPntsMax(100);
        aFinder.Perform()
        if aFinder.ErrorStatus():
            return None
        else:
            return aFinder.Shapes() # returns TopTools_ListOfShape

Note that for the example above (HybridMeshExample) I employed an interface to SMESH that I adapted from the python interface that is available in Salome. Using this interface makes working with SMESH through python feel a lot more natural and less tedious. IMHO adapting such an interface for PythonOCC would be very helpful for efficiently working with SMESH.

The solution with the ShapeFinder above is not very elegant. Wouldn’t it be nice if we could attach the information as to how Box2 (from the example above) can be meshed with prisms to the shape Box2 itself as a property ? And then find a way to make sure those properties stay valid along the history of modifications that happen to box2 ? For example, we could mark the top and bottom faces of Box2 as “top” and “bottom” (defining them as source and target face for a 2D projection algorithm) and all of it’s remaining faces as “wall faces”.

Attaching such properties to shapes is akin to what is called “topological naming” in Opencascade. It is supported by the Opencascade Application Framework (OCAF), which provides a complete document infrastructure with undo/redo mechanisms.
However, at least for the purposes discussed here, OCAF is IMHO overly complicated to work with. Consequently I implemented my own “lightweight” version of a “topological naming” scheme.

Suppose we have a class ShapeHistory that is able to track shape modifications and a class ShapeDocument (whose singleton instance is linked to a singleton ShapeHistory instance) that is able to store metadata on shapes and their modifications. Then we could accompany Box2 with additional data:

1
2
3
4
5
6
7
8
9
BB = BRepPrimAPI_MakeBox(0,0,-10,100,100,10);
Box2 = BB.Shape()
shapehistory.AddPrimitive(BB.Shape(), "Box2")
Botshape = shapehistory.GetShapeReference(BB.BottomFace())
Topshape = shapehistory.GetShapeReference(BB.TopFace())
shapedocument.SetProperties(BB.Shape(), {'Name':self.GetName(),
                                         'Material-ID':MatID,
                                         'TopShape':Topshape,
                     'BottomShape':Botshape})
BB = BRepPrimAPI_MakeBox(0,0,-10,100,100,10);
Box2 = BB.Shape()
shapehistory.AddPrimitive(BB.Shape(), "Box2")
Botshape = shapehistory.GetShapeReference(BB.BottomFace())
Topshape = shapehistory.GetShapeReference(BB.TopFace())
shapedocument.SetProperties(BB.Shape(), {'Name':self.GetName(),
                                         'Material-ID':MatID,
                                         'TopShape':Topshape,
	   				 'BottomShape':Botshape})

The tricky thing with this approach is to make sure our metadata stays up-to-date with the evolution of Box2. For this to work we need to tell shapehistory about any solids, faces and edges  that get deleted, modified or are being newly generated by any of Opencascades BRep algorithms.
Fortunately, any class in Opencascade that derives from BRepBuilderAPI_MakeShape (and accordingly any class in GEOM that derives from GEOMAlgo_BuilderShape) provides shape evolution support through it’s methods Generated(fromshape), Modified(shape) and IsDeleted(shape).  These functions can be used to build a “shape evolution tree” that consists of knots that hold shapes (or shape references) and edges that are marked as either “modified”, “generated” or “deleted”.
Now whenever we modify a shape by any of Opencascades algorithms, we use the information provided by that algorithm to update the tree structure. Implementing these update mechanisms is a somewhat tedious thing to do because many of Opencacades algorithms  provide their own shape history support methods that behave slightly different to the ones defined in BRepBuilderAPI_MakeShape. Note also that OCAF does not relief the developer from that coding part (at least not that I am aware of; there might be a commercial add-on available for that though).

Having implemented the shape history update mechanism, we can  now write:

1
2
3
4
5
6
7
8
Splitter = GEOMAlgo_Splitter()
Splitter.AddShape(Box)
Splitter.AddShape(Box2)
Splitter.AddShape(Sphere)
Splitter.Perform()
shapehistory.LoadShapeEvolutionFromBuilder(Splitter, 
                                           BuilderType.BuilderSplitter,
                                           "Splitter")
Splitter = GEOMAlgo_Splitter()
Splitter.AddShape(Box)
Splitter.AddShape(Box2)
Splitter.AddShape(Sphere)
Splitter.Perform()
shapehistory.LoadShapeEvolutionFromBuilder(Splitter, 
                                           BuilderType.BuilderSplitter,
                                           "Splitter")

Finally we can now retrieve all required information about how to mesh “Box2″ with prisms after it has been modified by the GEOMAlgo_Splitter algorithm:

1
2
3
4
5
6
7
8
9
10
11
# Get the list of shapes resulting from "Box2"
Box2FinalShapes = shapehistory.CurrentShape(Box2)
BottomSolid = Box2FinalShapes[0] # only 1 in this case
# All shapes resulting from "Box2" automatically
# inherit all the properties of "Box2"
TopShapeRef = shapedocument.GetProperty(BottomSolid, "TopShape")
# Shapes are always stored in the form of references
TopShape = shapehistory.GetShape(TopShapeRef)
# This is the original top face of Box2, most likely it was
# modified as well, get it's image
TopShapeFinalShapes = shapehistory.CurrentShape(TopShape)
# Get the list of shapes resulting from "Box2"
Box2FinalShapes = shapehistory.CurrentShape(Box2)
BottomSolid = Box2FinalShapes[0] # only 1 in this case
# All shapes resulting from "Box2" automatically
# inherit all the properties of "Box2"
TopShapeRef = shapedocument.GetProperty(BottomSolid, "TopShape")
# Shapes are always stored in the form of references
TopShape = shapehistory.GetShape(TopShapeRef)
# This is the original top face of Box2, most likely it was
# modified as well, get it's image
TopShapeFinalShapes = shapehistory.CurrentShape(TopShape)

With these techniques at hand we can now define mesh constraints, material IDs, boundary IDs, prism meshing metadata and so forth on primitive shapes right at the time they come to life. During subsequent CAD operations (like boolean operations, chamfer, fillet, draft, etc) these properties will automatically be updated and can finally be used during the mesh creation and mesh export stages.
This approach is  a lot more intuitive compared to having to define all those properties on the final geometry that might consist of dozens of solids and hundreds of faces and edges. Here’s an illustrative example of the shape evolution of a cylinder that was fused with a torus and subsequently overlain with a box using the Splitter:

Example shape evolution

Example shape evolution for a cylinder fused with a torus and subsequently overlain by a box using GEOMAlgo_Splitter.

Mesh export

Mesh elements in SMESH are stored in a tree of sub-meshes that mirrors the topological tree structure of the shape that was being meshed.  Or, put simply, SMESH allows accessing sub-mesh data structures for each solid, face, edge or vertex in a shape that was meshed.  This is perfect in case we have metadata attached to the subshapes of our CAD model.  Those metadata might be required for further mesh post-processing (e.g. for mesh export, mesh visualization or mesh-based numerical computations).

Let’s assume we have generated a compound shape consisting of several solids  with shared faces (resultshape), that we like to mesh and export for finite element computations. Here’s how we could do it if we have a  shapedocument that provides shape metadata and a MeshExporter that does the low-level mesh file writing available:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
# apply algorithms/hypotheses
meshgen = MT.MeshGen()
mesh = meshgen.Mesh(resultshape, "Mesh")
# only global 1D/2D/3D constraints
mesh.Segment(MT.REGULAR, resultshape).AutomaticLength(0.1)
mesh.Triangle(MT.NETGEN_2D, resultshape).LengthFromEdges()
mesh.Tetrahedron(MT.NETGEN_3D, resultshape).MaxElementVolume(1000)
# compute the mesh
mesh.Compute()
# export the mesh
# 1) mesh nodes
mainmeshds = mesh.GetMeshDS()
for inode in range(meshds.MinNodeID(), meshds.MaxNodeID()+1):
    nd = mainmeshds.FindNode(inode)
    MeshExporter.AddNode(inode, nd.X(), nd.Y(), nd.Z())
# 2) boundary triangles
FSMap = TopTools_IndexedDataMapOfShapeListOfShape()
TopExp().MapShapesAndAncestors(resultshape, TopAbs_FACE, TopAbs_SOLID, FSMap)
for face in Topo(resultshape).faces():
    nbsolidancestors = FSMap.FindFromIndex(FSMap.FindIndex(face)).Extent()
    if (nbsolidancestors>=2):
        continue # not a boundary face
    submesh = mesh.GetSubMesh(face)
    submeshds = submesh.GetSubMeshDS()
    faceboundaryid = shapedocument.GetProperty(face, "Boundary-ID")
    for ielem in range(submeshds.NbElements()):
    elem = submeshds.elemValue(ielem) # SMDS_MeshElement
        if (elem.GetEntityType()==SMDSEntity_Triangle):
            n1 = elem.GetNode(0).GetID(); n2 = elem.GetNode(1).GetID();
            n3 = elem.GetNode(2).GetID()
            MeshExporter.AddBoundaryTriangle(n1, n2, n3, faceboundaryid)
# 3) volume elements
for solid in Topo(resultshape).solids():
    volmesh = mesh.GetSubMesh(solid)
    volmeshds = submesh.GetSubMeshDS()
    materialid = shapedocument.GetProperty(solid, "Material-ID")
    for ielem in range(volmeshds.NbElements()):
        elem = volmeshds.elemValue(ielem) # SMDS_MeshElement
        if (elem.GetEntityType()==SMDSEntity_Tetra):
            n1 = elem.GetNode(0).GetID(); n2 = elem.GetNode(1).GetID()
            n3 = elem.GetNode(2).GetID(); n4 = elem.GetNode(3).GetID()
            MeshExporter.AddTetra(n1, n2, n3, n4, materialid)
MeshExporter.Write("meshfile.dat")
# apply algorithms/hypotheses
meshgen = MT.MeshGen()
mesh = meshgen.Mesh(resultshape, "Mesh")
# only global 1D/2D/3D constraints
mesh.Segment(MT.REGULAR, resultshape).AutomaticLength(0.1)
mesh.Triangle(MT.NETGEN_2D, resultshape).LengthFromEdges()
mesh.Tetrahedron(MT.NETGEN_3D, resultshape).MaxElementVolume(1000)
# compute the mesh
mesh.Compute()
# export the mesh
# 1) mesh nodes
mainmeshds = mesh.GetMeshDS()
for inode in range(meshds.MinNodeID(), meshds.MaxNodeID()+1):
    nd = mainmeshds.FindNode(inode)
    MeshExporter.AddNode(inode, nd.X(), nd.Y(), nd.Z())
# 2) boundary triangles
FSMap = TopTools_IndexedDataMapOfShapeListOfShape()
TopExp().MapShapesAndAncestors(resultshape, TopAbs_FACE, TopAbs_SOLID, FSMap)
for face in Topo(resultshape).faces():
    nbsolidancestors = FSMap.FindFromIndex(FSMap.FindIndex(face)).Extent()
    if (nbsolidancestors>=2):
        continue # not a boundary face
    submesh = mesh.GetSubMesh(face)
    submeshds = submesh.GetSubMeshDS()
    faceboundaryid = shapedocument.GetProperty(face, "Boundary-ID")
    for ielem in range(submeshds.NbElements()):
	elem = submeshds.elemValue(ielem) # SMDS_MeshElement
        if (elem.GetEntityType()==SMDSEntity_Triangle):
            n1 = elem.GetNode(0).GetID(); n2 = elem.GetNode(1).GetID();
            n3 = elem.GetNode(2).GetID()
            MeshExporter.AddBoundaryTriangle(n1, n2, n3, faceboundaryid)
# 3) volume elements
for solid in Topo(resultshape).solids():
    volmesh = mesh.GetSubMesh(solid)
    volmeshds = submesh.GetSubMeshDS()
    materialid = shapedocument.GetProperty(solid, "Material-ID")
    for ielem in range(volmeshds.NbElements()):
        elem = volmeshds.elemValue(ielem) # SMDS_MeshElement
        if (elem.GetEntityType()==SMDSEntity_Tetra):
            n1 = elem.GetNode(0).GetID(); n2 = elem.GetNode(1).GetID()
            n3 = elem.GetNode(2).GetID(); n4 = elem.GetNode(3).GetID()
            MeshExporter.AddTetra(n1, n2, n3, n4, materialid)
MeshExporter.Write("meshfile.dat")

Note that in the current version of PythonOCC the “elemValue” method for SMESHDS_SubMesh is not available, but can easily be added by extending the corresponding SMESH SWIG interface file SMESHDS.i:

1
2
3
4
5
6
7
8
%extend SMESHDS_SubMesh {
    const SMDS_MeshElement * elemValue(int index) {
    int i;
    SMDS_ElemIteratorPtr iterator=$self->GetElements();
    for (i=0; i<=index-1; i++) {iterator->next();}
    return iterator->next();
    }
};
%extend SMESHDS_SubMesh {
    const SMDS_MeshElement * elemValue(int index) {
	int i;
	SMDS_ElemIteratorPtr iterator=$self->GetElements();
	for (i=0; i<=index-1; i++) {iterator->next();}
	return iterator->next();
    }
};

This is certainly not a good solution in terms of efficiency – it should better be done e.g. by wrapping SMDS_ElemIteratorPtr. Here it just serves as a quick-fix solution.

Automatic prism meshing

Above I have shown how we can simplify prism meshing of solids by attaching metada to solid shapes (e.g. the information about “top”, “bottom” faces) and keep that information up-to-date  during subsequent CAD operations.
There are some situations where such information is readily available taking the information into account about how the shape considered was created. Box and cylinder shapes are simple examples of such situations:  BRepPrimAPI_MakeBox for example provides BottomFace() and TopFace() methods that can be used to distinguish bottom/top faces from wall faces for prism meshing.  Extrusion and thick solid shapes are two other examples I’d like to highlight in the following.
Again, here I’ll focus once more on prism meshing techniques, however, the general approach of extracting metadata from Opencascade BRep algorithms, attaching that data to the shapes created and maintaining their validity through the history of shape evolutions can be applied in a lot of different settings.

Automatic prism meshing of extrusion shapes

“Top”/”Bottom” shape metadata for automatic prism meshing is readily available for shapes created by a general (BRepOffsetAPI_MakePipe) or a linear (BRepPrimAPI_MakePrism) sweep assuming we use a face or a shell shape as the generator shape. For a general sweep we could add the required metadata like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def PerformGeneralSweep(path, generator, ShapeName, MaterialID):
    piper = BRepOffsetAPI_MakePipe(TopoDS.TopoDS().wire(path), generator)
    piper.Build(); result = piper.Shape()
    shapehistory.AddPrimitive(result, "General Sweep")
    # result is
    # - a solid in case generator is a face
    # - a compsolid in case generator is a shell
    result_solids = allSolids(result)
    generator_faces = allFaces(generator)
    path_edges = allEdges(path)
    for i, solid in enumerate(result_solids):
        shapehistory.AddPrimitive(solid, "General Sweep solid #"+str(i))
        shapedocument.SetProperties(solid, {'Name':ShapeName, 'MaterialID':MaterialID})
    firstshapefaces = allFaces(piper.FirstShape()) # same as allFaces(generator)
    lastshapefaces = allFaces(piper.LastShape())
    # top/bot shape metadata for prism meshing
    for i, solid in enumerate(result_solids):
        shapedocument.SetProperties(solid,
              {'BottomShape':shapehistory.GetShapeReference(firstshapefaces[i])})
        shapedocument.SetProperties(solid,
              {'TopShape':shapehistory.GetShapeReference(lastshapefaces[i])})
def PerformGeneralSweep(path, generator, ShapeName, MaterialID):
    piper = BRepOffsetAPI_MakePipe(TopoDS.TopoDS().wire(path), generator)
    piper.Build(); result = piper.Shape()
    shapehistory.AddPrimitive(result, "General Sweep")
    # result is
    # - a solid in case generator is a face
    # - a compsolid in case generator is a shell
    result_solids = allSolids(result)
    generator_faces = allFaces(generator)
    path_edges = allEdges(path)
    for i, solid in enumerate(result_solids):
        shapehistory.AddPrimitive(solid, "General Sweep solid #"+str(i))
        shapedocument.SetProperties(solid, {'Name':ShapeName, 'MaterialID':MaterialID})
    firstshapefaces = allFaces(piper.FirstShape()) # same as allFaces(generator)
    lastshapefaces = allFaces(piper.LastShape())
    # top/bot shape metadata for prism meshing
    for i, solid in enumerate(result_solids):
        shapedocument.SetProperties(solid,
              {'BottomShape':shapehistory.GetShapeReference(firstshapefaces[i])})
        shapedocument.SetProperties(solid,
              {'TopShape':shapehistory.GetShapeReference(lastshapefaces[i])})

Here allSolids, allFaces and allEdges are just my own convenient functions around TopExp_Explorer - PythonOCC offers the class Topo instead  for accessing the topology of shapes.
Basically the same approach as show here works for shapes generated by the class BRepPrimAPI_MakePrism, that also provides FirstShape() and LastShape() methods. Both BRepPrimAPI_MakePrism and BRepOffsetAPI_MakePipe create the resulting shapes in a way such that the sequence of faces and edges is the same in the subshapes returned by FirstShape() and LastShape(). Furthermore face number i in the subshape returned by FirstShape() always corresponds to solid number i in the result shape (returned by  Shape()). This makes automatic prism meshing of such geometries particularly simple.

Automatic prism meshing of thick solid shapes

Hybrid mesh with tetrahedral and prismatic boundary layer mesh.

Hybrid mesh with tetrahedral and prismatic boundary layer mesh. The boundary layer was created using BRepOffset_MakeOffset.

As shown above we can employ SMESH’s new “viscous layers” algorithm to add prismatic  boundary meshes on solid shapes meshed with the Netgen tetrahedral mesher. In contrast to this mesh-based solution we could also take another approach for realizing a boundary layer mesh:  Create a boundary layer geometry, add it to the solid considered and fill it with a prism mesh. This approach could be used in cases where the “viscous layers” algorithm fails, or in cases where the boundary layer and the inner volume mesh represent  different material regions.
We can employ BRepOffset_MakeOffset to compute the offset shape of a certain solid. From the outer shell of that offset shape and the outer shell of the original solid we can then create a hollow solid using BRepBuilderAPI_MakeSolid. This hollow solid is the boundary layer shape that we seek to mesh with prisms. Unfortunately, the sequence of the faces in the two shells of the boundary layer shape is not the same (the same is true for the edges). But we can ask the BRepOffset_MakeOffset class  for that information and create corresponding association maps (EdgeAssoc and FaceAssoc) for later automatic prism meshing:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
def AddBoundaryLayer(sourceshape, thickness, materialID):
    mode = BRepOffset_Skin; Intersection = 0; SelfInter = 0; Join = GeomAbs_Intersection
    offsetbuilder = BRepOffset_MakeOffset()
    offsetbuilder.Initialize(sourceshape,thickness,1.e-5,mode,Intersection,SelfInter,Join)
    offsetbuilder.MakeOffsetShape()
    offsetshape = builder.Shape()
    # extract shell from sourceshape
    Ex = TopExp_Explorer(sourceshape,TopAbs_SHELL)
    innershell = TopoDS().shell(Ex.Current())
    # extract shell from offsetshape
    Ex2 = TopExp_Explorer(offsetshape,TopAbs_SHELL)
    outershell = TopoDS.TopoDS().shell(Ex2.Current())
    # build a hollow solid from the two shells
    BMS = BRepBuilderAPI_MakeSolid(innershell, outershell)
    boundarylayer = BMS.Shape()
    shapehistory.AddPrimitive(boundarylayer, "Boundary Layer")
    shapedocument.SetProperties(boundarylayer,
        {'Name':"Boundary Layer", 'MatID':materialID})
    shapedocument.SetProperties(boundarylayer,
        {'BottomShape':shapehistory.GetShapeReference(innershell)})
    shapedocument.SetProperties(boundarylayer,
        {'TopShape':shapehistory.GetShapeReference(outershell)})
    # find edge and face association asking the offset builder
    is_faces = allFaces(innershell); is_edges = allEdges(innershell)
    os_faces_map = TopTools_IndexedMapOfShape()
    os_edges_map = TopTools_IndexedMapOfShape()
    TopExp().MapShapes(outershell, TopAbs_FACE, os_faces_map)
    TopExp().MapShapes(outershell, TopAbs_EDGE, os_edges_map)
    faceimages = offsetbuilder.OffsetFacesFromShapes() # BRepAlgo_Image
    edgeimages = offsetbuilder.OffsetEdgesFromShapes()
    FaceAssoc = dict(); EdgeAssoc = dict()
    # find association inner shell faces --> outer shell faces
    for is_face_id, is_face in enumerate(is_faces):
        generatedfaces = TopTools_ListOfShape()
        faceimages.LastImage(is_face, generatedfaces)
        if generatedfaces.Extent()!=1:
            raise RuntimeError, "Expected exactly 1 face generated from original face"
        generatedface = generatedfaces.Value()
        if (os_faces_map.Contains(generatedface)!=1):
            raise RuntimeError, "This generated face does not belong to the outer shell."
        FaceAssoc[is_face_id] = os_faces_map.FindIndex(generatedface)-1
    # find association inner shell edges --> outer shell edges
    for is_edge_id, is_edge in enumerate(is_edges):
        generatededges = TopTools_ListOfShape()
        edgeimages.LastImage(is_edge, generatededges)
        if generatededges.Extent()!=1:
            raise RuntimeError, "Expected exactly 1 edge generated from original edge"
        generatededge = generatededge.Value()
        if (os_edges_map.Contains(generatededge)!=1):
            raise RuntimeError, "This generated edge does not belong to the outer shell."
        EdgeAssoc[is_edge_id] = os_edges_map.FindIndex(generatededge)-1
    self.shapedocument.SetProperties(boundarylayer,
        {'FaceAssoc':FaceAssoc, 'EdgeAssoc':EdgeAssoc})
    # build a compound of sourceshape and boundarylayer to return as result shape
    result = TopoDS_Compound()
    BB = BRep_Builder()
    BB.MakeCompound(result)
    BB.Add(result, sourceshape)
    BB.Add(result, boundarylayer)
    return result
def AddBoundaryLayer(sourceshape, thickness, materialID):
    mode = BRepOffset_Skin; Intersection = 0; SelfInter = 0; Join = GeomAbs_Intersection
    offsetbuilder = BRepOffset_MakeOffset()
    offsetbuilder.Initialize(sourceshape,thickness,1.e-5,mode,Intersection,SelfInter,Join)
    offsetbuilder.MakeOffsetShape()
    offsetshape = builder.Shape()
    # extract shell from sourceshape
    Ex = TopExp_Explorer(sourceshape,TopAbs_SHELL)
    innershell = TopoDS().shell(Ex.Current())
    # extract shell from offsetshape
    Ex2 = TopExp_Explorer(offsetshape,TopAbs_SHELL)
    outershell = TopoDS.TopoDS().shell(Ex2.Current())
    # build a hollow solid from the two shells
    BMS = BRepBuilderAPI_MakeSolid(innershell, outershell)
    boundarylayer = BMS.Shape()
    shapehistory.AddPrimitive(boundarylayer, "Boundary Layer")
    shapedocument.SetProperties(boundarylayer,
        {'Name':"Boundary Layer", 'MatID':materialID})
    shapedocument.SetProperties(boundarylayer,
        {'BottomShape':shapehistory.GetShapeReference(innershell)})
    shapedocument.SetProperties(boundarylayer,
        {'TopShape':shapehistory.GetShapeReference(outershell)})
    # find edge and face association asking the offset builder
    is_faces = allFaces(innershell); is_edges = allEdges(innershell)
    os_faces_map = TopTools_IndexedMapOfShape()
    os_edges_map = TopTools_IndexedMapOfShape()
    TopExp().MapShapes(outershell, TopAbs_FACE, os_faces_map)
    TopExp().MapShapes(outershell, TopAbs_EDGE, os_edges_map)
    faceimages = offsetbuilder.OffsetFacesFromShapes() # BRepAlgo_Image
    edgeimages = offsetbuilder.OffsetEdgesFromShapes()
    FaceAssoc = dict(); EdgeAssoc = dict()
    # find association inner shell faces --> outer shell faces
    for is_face_id, is_face in enumerate(is_faces):
        generatedfaces = TopTools_ListOfShape()
        faceimages.LastImage(is_face, generatedfaces)
        if generatedfaces.Extent()!=1:
            raise RuntimeError, "Expected exactly 1 face generated from original face"
        generatedface = generatedfaces.Value()
        if (os_faces_map.Contains(generatedface)!=1):
            raise RuntimeError, "This generated face does not belong to the outer shell."
        FaceAssoc[is_face_id] = os_faces_map.FindIndex(generatedface)-1
    # find association inner shell edges --> outer shell edges
    for is_edge_id, is_edge in enumerate(is_edges):
        generatededges = TopTools_ListOfShape()
        edgeimages.LastImage(is_edge, generatededges)
        if generatededges.Extent()!=1:
            raise RuntimeError, "Expected exactly 1 edge generated from original edge"
        generatededge = generatededge.Value()
        if (os_edges_map.Contains(generatededge)!=1):
            raise RuntimeError, "This generated edge does not belong to the outer shell."
        EdgeAssoc[is_edge_id] = os_edges_map.FindIndex(generatededge)-1
    self.shapedocument.SetProperties(boundarylayer,
        {'FaceAssoc':FaceAssoc, 'EdgeAssoc':EdgeAssoc})
    # build a compound of sourceshape and boundarylayer to return as result shape
    result = TopoDS_Compound()
    BB = BRep_Builder()
    BB.MakeCompound(result)
    BB.Add(result, sourceshape)
    BB.Add(result, boundarylayer)
    return result

Note that this approach does not work for concave geometries for which faces would vanish during the offset operation.

Mixed-Material regions in CAD models

For finite element simulations we often need to have volume meshes with several separated regions that we can attach different material properties to. Such volume meshes can be created by discretizing compound shapes that consist of several solids (each representing a different material region) with  shared faces in between.  Faces (and edges) between the solids must be shared to ensure we are able obtain a regular surface mesh, which is required to compute the final volume mesh.

Here’s an example of a geometry that consists of three solids with shared faces:

Hybrid structured/unstructured mesh example on a compound solid with shared faces

Example of a compound geometry consisting of three solids. Two solids are discretized with tetrahedra, one with prismatic cells using arithmetically increasing cell layer thicknesses.

Unfortunately, there is no general method available in Opencascade to create such geoemetries. There are some special cases, where creating them is relatively simple:  One example are the two fused boxes covered by a boundary layer from the prism boundary mesh example above.
But constructing such geometries for the general case is not directly supported by Opencascade because they do not follow the assumption that each edge in a model is only  connected to at most two faces. Models that fulfill this assumption exhibit a “two-manifold” topology (loosely spoken). The spherical shell example above has several edges that are connected to three faces.
If we think about a sphere that is overlain by a box as a simple example, we can easily see that the construction of such a geometry can not be carried out using  any of Opencascade’s boolean operations (or a combination of them):

About the construction of compound solid shapes with shared faces

The construction of compound shapes with shared faces can not be easily done in OCC using any of (or a combination of ) boolean operations.

Luckily, PythonOCC provides some support for the construction of geometries with a “non-manifold” topology through it’s python interface to a standalone version of  Salome’s GEOM module0. The GEOM module is a complete OpenCascade – OCAF based CAD framework. Among many other features it provides an advanced partition/glueing algorithm (GEOMAlgo_Splitter), that supports the creation of compound solids with shared faces (3D problems) as well as compound faces with shared edges (2D problems). Using GEOMAlgo_Splitter the problem of overlaying a sphere with a box can be solved like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
def OverlaySpherBox():
    Box = BRepPrimAPI_MakeBox(gp_Pnt(0.0,0.0,0.0), 1.5,1.5,1.5).Shape()
    Sphere = BRepPrimAPI_MakeSphere(gp_Ax2(gp_Pnt(0,0,0), gp_Dir(0,-1,0)), 1.0).Shape()
    # Perform splitting
    Splitter = GEOMAlgo_Splitter()
    Splitter.AddShape(Sphere)
    Splitter.AddShape(Box)
    Splitter.Perform()
    result = Splitter.Shape()
    # recover resulting sphere solids: all solids that are in Image(Sphere)
    #                                  but not in Image(Box)
    spheresolids = TopTools_ListOfShape()
    ImagesSphere=TopTools_ListOfShape(); Splitter.Images().LastImage(Sphere, ImagesSphere)
    ImagesBox=TopTools_ListOfShape(); Splitter.Images().LastImage(Box, ImagesBox)
    ImagesBoxSet = TopTools_MapOfShape()
    it = TopTools_ListIteratorOfListOfShape(ImagesBox)
    while it.More():
        ImagesBoxSet.Add(it.Value()); it.Next()
    it2 = TopTools_ListIteratorOfListOfShape(ImagesSphere)
    while it2.More():
        solid = it2.Value();
        if not (ImagesBoxSet.Contains(solid)==1):
            spheresolids.Append(solid)
        it2.Next()
    newSphere = TopTools_ListIteratorOfListOfShape(spheresolids).Value() # only one in this case
    # recover new box that has faces shared  with the sphere
    #   collect all faces of all shapes in Image(Box) that are not in "InBoxShapeSet"
    SB  = BRepBuilderAPI_MakeSolid(); abuilder = BRep_Builder()
    shell=TopoDS_Shell(); abuilder.MakeShell(shell)
    InBoxShapeSet = TopTools_MapOfShape()
    InBoxShapeList = Splitter.InParts(Box)
    it3 = TopTools_ListIteratorOfListOfShape(InBoxShapeList)
    while it3.More():
        InBoxShapeSet.Add(it3.Value()); it3.Next()
    it4 = TopTools_ListIteratorOfListOfShape(ImagesBox)
    while it4.More():
        shape = it4.Value()
        Ex = TopExp_Explorer(shape, TopAbs_FACE)
        while Ex.More():
            face = TopoDS().face(Ex.Current())
            if not (InBoxShapeSet.Contains(face)==1):
                abuilder.Add(shell, face)
            Ex.Next()
        it4.Next()
    SB.Add(shell)
    newBox = SB.Solid()
    return newBox, newSphere
def OverlaySpherBox():
    Box = BRepPrimAPI_MakeBox(gp_Pnt(0.0,0.0,0.0), 1.5,1.5,1.5).Shape()
    Sphere = BRepPrimAPI_MakeSphere(gp_Ax2(gp_Pnt(0,0,0), gp_Dir(0,-1,0)), 1.0).Shape()
    # Perform splitting
    Splitter = GEOMAlgo_Splitter()
    Splitter.AddShape(Sphere)
    Splitter.AddShape(Box)
    Splitter.Perform()
    result = Splitter.Shape()
    # recover resulting sphere solids: all solids that are in Image(Sphere)
    #                                  but not in Image(Box)
    spheresolids = TopTools_ListOfShape()
    ImagesSphere=TopTools_ListOfShape(); Splitter.Images().LastImage(Sphere, ImagesSphere)
    ImagesBox=TopTools_ListOfShape(); Splitter.Images().LastImage(Box, ImagesBox)
    ImagesBoxSet = TopTools_MapOfShape()
    it = TopTools_ListIteratorOfListOfShape(ImagesBox)
    while it.More():
        ImagesBoxSet.Add(it.Value()); it.Next()
    it2 = TopTools_ListIteratorOfListOfShape(ImagesSphere)
    while it2.More():
        solid = it2.Value();
        if not (ImagesBoxSet.Contains(solid)==1):
            spheresolids.Append(solid)
        it2.Next()
    newSphere = TopTools_ListIteratorOfListOfShape(spheresolids).Value() # only one in this case
    # recover new box that has faces shared  with the sphere
    #   collect all faces of all shapes in Image(Box) that are not in "InBoxShapeSet"
    SB  = BRepBuilderAPI_MakeSolid(); abuilder = BRep_Builder()
    shell=TopoDS_Shell(); abuilder.MakeShell(shell)
    InBoxShapeSet = TopTools_MapOfShape()
    InBoxShapeList = Splitter.InParts(Box)
    it3 = TopTools_ListIteratorOfListOfShape(InBoxShapeList)
    while it3.More():
        InBoxShapeSet.Add(it3.Value()); it3.Next()
    it4 = TopTools_ListIteratorOfListOfShape(ImagesBox)
    while it4.More():
        shape = it4.Value()
        Ex = TopExp_Explorer(shape, TopAbs_FACE)
        while Ex.More():
            face = TopoDS().face(Ex.Current())
            if not (InBoxShapeSet.Contains(face)==1):
                abuilder.Add(shell, face)
            Ex.Next()
        it4.Next()
    SB.Add(shell)
    newBox = SB.Solid()
    return newBox, newSphere

Using GEOMAlgo_Splitter fairly general 2D/3D algorithms can be developed for the creation of shapes that support mixed-material regions.

Shape/shape contacts

In addition to mixed-material regions we often need to establish face or edge contacts between solids within a model that we like to treat with finite element or other mesh-based numerical simulation techniques. Think of the simple scenario of a spherical or a cylindrical body placed on a substrate layer:

Shape/shape contact situations often required for FE analysis

Different possibilities for solid shapes to be in contact. Solid/solid and face/face contacts are the most important geometric situations occurring in FE models. Point/face and edge/face contacts are more difficult to realize.

The “solid/solid contact” on the left side corresponds to the situation where the top body slightly dips into the layer material. It can be realized in the same way as the example with the sphere and the box above (and technically of course it is a face/face contact situation). The face/face contact on the right side is the most important scenario and can be realized using the GEOMAlgo_Splitter:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def CylinderOnPlate():
    Box = BRepPrimAPI_MakeBox(gp_Pnt(0,0,0), 100,100,20).Shape()
    Cyl = BRepPrimAPI_MakeCylinder(gp_Ax2(gp_Pnt(50,50,20), gp_Dir(0,0,1)), 25, 50).Shape()
    Splitter = GEOMAlgo_Splitter()
    Splitter.AddShape(Box)
    Splitter.AddShape(Cyl)
    Splitter.Perform()
    result=Splitter.Shape()
    # find the shared face
    sharedface = None
    FSMap = TopTools_IndexedDataMapOfShapeListOfShape()
    TopExp().MapShapesAndAncestors(result, TopAbs_FACE, TopAbs_SOLID, FSMap)
    Ex = TopExp_Explorer(result, TopAbs_FACE)
    while Ex.More():
        face = TopoDS().face(Ex.Current())
        index = FSMap.FindIndex(face)
        if (index!=0) and (FSMap.FindFromIndex(index).Extent()==2):
            sharedface = face; break
        Ex.Next()
    return result, sharedface
def CylinderOnPlate():
    Box = BRepPrimAPI_MakeBox(gp_Pnt(0,0,0), 100,100,20).Shape()
    Cyl = BRepPrimAPI_MakeCylinder(gp_Ax2(gp_Pnt(50,50,20), gp_Dir(0,0,1)), 25, 50).Shape()
    Splitter = GEOMAlgo_Splitter()
    Splitter.AddShape(Box)
    Splitter.AddShape(Cyl)
    Splitter.Perform()
    result=Splitter.Shape()
    # find the shared face
    sharedface = None
    FSMap = TopTools_IndexedDataMapOfShapeListOfShape()
    TopExp().MapShapesAndAncestors(result, TopAbs_FACE, TopAbs_SOLID, FSMap)
    Ex = TopExp_Explorer(result, TopAbs_FACE)
    while Ex.More():
        face = TopoDS().face(Ex.Current())
        index = FSMap.FindIndex(face)
        if (index!=0) and (FSMap.FindFromIndex(index).Extent()==2):
            sharedface = face; break
        Ex.Next()
    return result, sharedface

Point/face contacts are difficult to enforce on the CAD model and rather should be realized during the discretization stage. In the example above (sphere on substrate layer) it was only possible to realize by cutting “a cross” made of two edges into the top face of the substrate box.

Face/edge contact example

Face/edge contacts between solid shapes can be established if the contact edges form closed loops on the faces they cut.

Edge/edge contacts, like the one shown above,  can be realized with GEOMAlgo_Splitter as well. However, meshing will fail for cases where the edges shared between the objects under consideration do not form closed loops. This is due to the inability of the SMESH surface mesher (Netgen and Mefisto2) to tread faces with internal edges. Edges that form closed loops will cut the face they interfere with into several faces that each can be surface discretized without any issue.