# 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:

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

- deriving from the class
*SMESH_3D_Algo*, (or*_2D_*/*_1D_*) for developing the mesh algorithm, - 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 - 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:

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:

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.

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:

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:

### 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

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:

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):

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:

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.

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.

*amazing* contribution Mark, thanks so much!

I think the MeshExporter you refer to is a custom file format for your in-house solver, right? For most cases, mesh.ExportMED / mesh.ExportUNV would be the call to make, for instance for exporting to Code-Aster.

I’m curious if you have an idea how to achieve this: I’d love to map the scalar values [ say Von Mises strain ] onto the BRep geometry, such to have a close coupling from the FEM mesh to the original geometry. Are you aware of a method to interpolate a point on the BRep to that of the resulting FEM values?

Many, many thanks for such a detailed and interesting article. I hope to be able to ask you a few more things on the subject of GEOM’s partitioning algorithms one of these days.

Thanks a lot for your comments Jelle! The MeshExporter class I refer to can be any custom file format writer. The intention of the mesh export example was just to show how BRep related metadata (boundary IDs for example) can be taken into account during mesh export (in my case as you pointed out, MeshExporter writes a format suitable for our in-house FEM solver).

About interpolating FEM results to points on the BRep (not sure I understood your question correctly, please correct me if not): As far as I know SMESH keeps the (u,v) coordinates of all mesh nodes for discretized faces and the (u) coordinates of all mesh nodes for discretized edges. And using Opencascade we can project a 3D point on an edge or a face to find out about it’s (u) or (u,v) coordinates.

Given these informations the task of interpolating a FEM result to a point on a BRep (an edge or a face) is as simple as doing such an interpolation in a 2D plane. So I think we only need to find the triangle the point is located in, get it’s barycentric coordinates within that triangle and use the shape functions to interpolate the solution value. Would be interesting to think about how to implement such a link in an efficient manner.

That’s really an interesting question… we are currently thinking about how to implement isoparametric elements, for which we would need that link as well. I can keep you up-to-date on this if I find out more about it.