Added contributions from Ivor Clifford

python/Salome/extractFeatureEdges.py

+#!python
+# =============================================================================
+# Salome GEOM script to extract the feature edges from a body and add them
+# to the group named 'featureEdges'.
+# Tested on Salome 7.4.0 and python 2.7 on 64-bit Windows
+#
+# Author: Ivor Clifford <ivor.clifford@psi.ch>
+#
+# =============================================================================
+
+def extractFeatureEdges(body, minFeatureAngle = 5):
+    '''
+    Find all feature edges on the supplied body and return them as a list
+    of edge ids.
+
+        body - A Salome solid, compound, shell or face object to find all
+                feature edges on.
+        minFeatureAngle - the angle (in degrees) between adjacent surfaces
+                above which the edge will be considered a feature angle.
+    '''
+    import salome
+    salome.salome_init()
+
+    import GEOM
+    from salome.geom import geomBuilder
+    geompy = geomBuilder.New(salome.myStudy)
+    
+    # Check the body type
+    if not (body.GetShapeType() in [GEOM.SHELL, GEOM.SOLID, GEOM.FACE, GEOM.COMPOUND]):
+       raise RuntimeError('Supplied object is not a solid, shell or face.')
+    
+    print 'Extracting edges of %s with feature angle > %g.' % (body.GetName(), minFeatureAngle)
+    
+    # Extract basic info
+    faces = geompy.SubShapeAll(body, geompy.ShapeType["FACE"])
+    curves = geompy.SubShapeAll(body, geompy.ShapeType["EDGE"])
+    points = geompy.SubShapeAll(body, geompy.ShapeType["VERTEX"])
+    
+    faceIds = geompy.GetSubShapesIDs(body, faces)
+    curveIds = geompy.GetSubShapesIDs(body, curves)
+    nodeIds = geompy.GetSubShapesIDs(body, points)
+    
+    maxFaceId = max(faceIds)
+    maxCurveId = max(curveIds)
+    maxNodeId = max(nodeIds)
+    
+    # Reverse mapping from curve id to local curve arrays
+    faceMap = [-1 for i in xrange(maxFaceId+1)]
+    for localId, id in enumerate(faceIds):
+        faceMap[id] = localId
+    
+    curveMap = [-1 for i in xrange(maxCurveId+1)]
+    for localId, id in enumerate(curveIds):
+        curveMap[id] = localId
+    
+    nodeMap = [-1 for i in xrange(maxNodeId+1)]
+    for localId, id in enumerate(nodeIds):
+        nodeMap[id] = localId
+    
+    
+    # Get curves on each face
+    faceCurveIds = [[curveMap[id] for id in geompy.GetSubShapesIDs(
+            body,
+            geompy.SubShapeAll(face, geompy.ShapeType["EDGE"])
+        )] for face in faces]
+    
+    # Get faces attached to each curve
+    curveFaceIds = [[] for id in curveIds]
+    
+    for faceI, ids in enumerate(faceCurveIds):
+        for id in ids:
+            curveFaceIds[id].append(faceI)
+    
+    # Now that we have the connectivity for curves and faces find the
+    # feature edges
+    featureEdgeIds = []
+    for curveId, curve, adjFaceIds in zip(curveIds, curves, curveFaceIds):
+        if len(adjFaceIds) == 2:
+            # Curve with 2 adjacent faces - Test feature angle
+            # Determine break angle at each curve
+            # If greater than the feature edge angle, add the curve to group featureEdges
+            face1 = faces[adjFaceIds[0]]
+            face2 = faces[adjFaceIds[1]]
+            point = geompy.GetFirstVertex(curve)    # Test at the first vertex
+            n1 = geompy.GetNormal(face1, point)
+            n2 = geompy.GetNormal(face2, point)
+            angle = geompy.GetAngle(n1, n2)
+            if angle > minFeatureAngle:
+                featureEdgeIds.append(curveId)
+        
+        elif len(adjFaceIds) == 1:
+            # Curve on standalone face - Add by default
+            featureEdgeIds.append(curveId)
+        
+        elif len(adjFaceIds) == 0:
+            # Standalone curve - Ignore
+            None
+        
+        else:
+            raise RuntimeError('Curve found sharing %d faces. This is unexpected for fully enclosed bodies.' % len(adjFaceIds))
+    
+    # Done
+    print "%d feature edges found" % len(featureEdgeIds)
+    
+    return featureEdgeIds
+
+
+# If run as a standalone script, use the current Salome GUI selection
+# and add the feature edges to group named 'featureEdges'
+if __name__ == '__main__':
+    import salome
+    salome.salome_init()
+    
+    import GEOM
+    from salome.geom import geomBuilder
+    geompy = geomBuilder.New(salome.myStudy)
+    
+    # Get current GUI selection
+    selected = salome.sg.getAllSelected()
+    if len(selected) != 1:
+        raise RuntimeError('A single solid, shell or face object must be selected.')
+    
+    body = salome.myStudy.FindObjectID(selected[0]).GetObject()
+    
+    # Get feature edges and add to the group 'featureEdges'
+    featureEdges = geompy.CreateGroup(body, geompy.ShapeType["EDGE"])
+    geompy.UnionIDs(featureEdges, extractFeatureEdges(body))
+    geompy.addToStudyInFather(body, featureEdges, 'featureEdges')
+    
+    if salome.sg.hasDesktop():
+        salome.sg.updateObjBrowser(1)

python/Salome/salomeTriSurf.py

+#!python
+# =============================================================================
+# Python module for writing OpenFOAM feature edge and triSurface files from
+# within Salome platform.
+# Tested on Salome 7.4.0 and python 2.7 on 64-bit Windows
+#
+# Author: Ivor Clifford <ivor.clifford@psi.ch>
+#
+# =============================================================================
+
+def foamHeader(className, objectName):
+    '''
+    Return the OpenFOAM file header block as a string.
+    '''
+    return '''/*--------------------------------*- C++ -*----------------------------------*\\
+| =========                 |                                                 |
+| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\\\    /   O peration     | Version:  2.2.1                                 |
+|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\\\/     M anipulation  |                                                 |
+\\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       %s;
+    object      %s;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+''' % (className, objectName)
+
+
+class triSurf:
+    def __init__(self, object = None, allEdges = False):
+        '''
+        Construct from the supplied Salome mesh.
+
+            object - the mesh object (must be a triangular surface mesh). If no
+                    object is supplied, the current Salome selection is used.
+            allEdges - If true, all edges on the mesh are included as feature
+                    edges, otherwise only edges contained in groups are included
+
+        NOTE: All face groups are assumed to represent patches. No face subsets
+        are written. All edge groups are added as feature edge subsets. All point
+        groups are added as point subsets.
+        '''
+        # =============================================================================
+        # Initialize salome
+        import salome
+        import SMESH, SALOMEDS
+        from salome.smesh import smeshBuilder
+        from operator import itemgetter
+        from collections import OrderedDict
+        
+        import SMESH, SALOMEDS
+        
+        # =============================================================================
+        # Get the Salome mesh object
+        if object is None:
+            selected = salome.sg.getAllSelected()
+            if len(selected) != 1:
+                raise RuntimeError('A single Salome mesh object must be selected.')
+            
+            object = salome.myStudy.FindObjectID(selected[0]).GetObject()
+        
+        try:
+            object.GetMesh()
+        except:
+            raise RuntimeError('Supplied object is not a Salome SMESH mesh.')
+        
+        smesh = smeshBuilder.New(salome.myStudy)
+        mesh = smesh.Mesh(object)
+        
+        print "Converting SMESH Mesh '%s'" % mesh.GetName()
+        
+        # =============================================================================
+        # Get basic mesh info
+        nNodes = mesh.NbNodes()
+        nFaces = mesh.NbFaces()
+        nTris = mesh.NbTriangles()
+        nEdges = mesh.NbEdges()
+        nodeIds = mesh.GetNodesId()
+        faceIds = mesh.GetElementsByType(SMESH.FACE)
+        edgeIds = mesh.GetElementsByType(SMESH.EDGE)
+        
+        # Check that mesh is strictly triangular
+        if nFaces != nTris:
+            raise RuntimeError('Mesh is not strictly triangular')
+        
+        # Get patch and subset names & ids
+        # All SMESH.FACE groups are assumed to be patches
+        # All SMESH.EDGE groups are assumed to be feature subsets
+        # All SMESH.NODE groups are assumed to be point subsets
+        patches = OrderedDict()
+        pointSubsets = OrderedDict()
+        featureEdgeSubsets = OrderedDict()
+        
+        for group in mesh.GetGroups():
+            if group.GetType() == SMESH.FACE:
+                patches[group.GetName()] = group.GetIDs()
+            elif group.GetType() == SMESH.EDGE:
+                featureEdgeSubsets[group.GetName()] = group.GetIDs()
+            elif group.GetType() == SMESH.NODE:
+                pointSubsets[group.GetName()] = group.GetIDs()
+        
+        # =============================================================================
+        # Process faces and patches
+        # Get patchId for each face
+        lastPatchId = len(patches)
+        patchIds = [lastPatchId] * max(faceIds)
+        patchId = 0
+        for name, ids in patches.iteritems():
+            for faceId in ids:
+                if patchIds[faceId-1] == lastPatchId:
+                    patchIds[faceId-1] = patchId
+                else:
+                    print "Face %d is assigned to both groups %s and %s" % (faceId, name, patches.keys()[patchIds[faceId-1]])
+                    raise RuntimeError('Groups of faces are not unique, i.e. they overlap.')
+                
+            patchId += 1
+        
+        # Compact and reorder patchIds to match faceIds
+        patchIds = [patchIds[faceId-1] for faceId in faceIds]
+        
+        # Reorder faces by increasing group id
+        faceAndpatchIds = sorted(zip(faceIds, patchIds), key=itemgetter(1))
+        faceIds, patchIds = zip(*faceAndpatchIds)
+        
+        # Add unused faces to the default patch
+        defaultFaces = [faceId for faceId, patchId in faceAndpatchIds if patchId == lastPatchId]
+        if len(defaultFaces) > 0:
+            patches['defaultFaces'] = defaultFaces
+        
+        defaultFaces = None
+        
+        # =============================================================================
+        # Process feature edges
+        if not allEdges:
+            edgeIds = []
+            for name, ids in featureEdgeSubsets.iteritems():
+                edgeIds += ids
+            
+            edgeIds = list(set(edgeIds))
+            nEdges = len(edgeIds)
+        
+        # Reverse mapping of edge ids since they aren't necessarily numbered 1..nEdges
+        if len(edgeIds):
+            edgeMap = [-1] * max(edgeIds)
+        else:
+            edgeMap = []
+        
+        i=0
+        for edgeId in edgeIds:
+            edgeMap[edgeId-1] = i
+            i += 1
+        
+        # =============================================================================
+        # Process nodes
+        # Reverse mapping of node ids since nodes aren't necessarily numbered 1..nNodes
+        nodeMap = [-1] * max(nodeIds)
+        i=0
+        for nodeId in nodeIds:
+            nodeMap[nodeId-1] = i
+            i += 1
+        
+        # =============================================================================
+
+        self._mesh = mesh
+        
+        self._nodeIds = nodeIds
+        self._edgeIds = edgeIds
+        self._faceIds = faceIds
+        
+        self._nodeMap = nodeMap
+        self._edgeMap = edgeMap
+        self._faceMap = []
+        
+        self._patches = patches
+        self._pointSubsets = pointSubsets
+        self._featureEdgeSubsets = featureEdgeSubsets
+        self._faceSubsets = {}
+        
+        print 'Done'
+    
+    def nNodes(self):
+        '''
+        Return the number of nodes
+        '''
+        return len(self._nodeIds)
+    
+    def nEdges(self):
+        '''
+        Return the number of edges
+        '''
+        return len(self._edgeIds)
+    
+    def nFacets(self):
+        '''
+        Return the number of triangular facets
+        '''
+        return len(self._faceIds)
+    
+    def nPatches(self):
+        '''
+        Return the number of patches
+        '''
+        return len(self._patches)
+    
+    def _writePatchDefs(self, f, typeName = 'wall'):
+        '''
+        Write the patch definitions to file as an OpenFOAM geometricSurfacePatchList.
+        NOTE: All patches are assumed to be walls.
+        '''
+        patches = self._patches
+        
+        f.write('%d\n(\n' % len(patches))
+        for name in patches.iterkeys():
+            f.write('%s\t%s\n' % (name, typeName))
+        
+        f.write(')\n')
+    
+    def _writeNodes(self, f):
+        '''
+        Write the nodes to file as an OpenFOAM pointField.
+        '''
+        mesh = self._mesh
+        nodeIds = self._nodeIds
+        
+        f.write('%d\n(\n' % len(nodeIds))
+        
+        for x, y, z in [mesh.GetNodeXYZ(nodeId) for nodeId in nodeIds]:
+            f.write( '( %g %g %g )\n' % (x, y, z))
+        
+        f.write(')\n')
+    
+    def _writeFeatureEdges(self, f):
+        '''
+        Write the feature edges to file as an OpenFOAM edgeList.
+        '''
+        mesh = self._mesh
+        nodeMap = self._nodeMap
+        edgeIds = self._edgeIds
+        
+        f.write('%d\n(\n' % len(edgeIds))
+        
+        for edgeId in edgeIds:
+            nodes = mesh.GetElemNodes(edgeId)
+            f.write( '(' + ' '.join([str(nodeMap[nodeId-1]) for nodeId in nodes]) + ')\n')
+        
+        f.write(')\n')
+    
+    def _writeFacets(self, f):
+        '''
+        Write the facets to file as an OpenFOAM List of labelledTri.
+        '''
+        from itertools import chain
+        
+        mesh = self._mesh
+        nodeMap = self._nodeMap
+        patches = self._patches
+        
+        f.write('%d\n(\n' % sum([len(patch) for patch in patches.itervalues()]))
+        
+        patchId = 0
+        for patchId, (patchName, faceIds) in enumerate(patches.iteritems()):
+            for faceId in faceIds:
+                nodes = mesh.GetElemNodes(faceId)
+                f.write( '((' + ' '.join([str(nodeMap[nodeId-1]) for nodeId in nodes]) + ') %d)\n' % patchId)
+        
+        f.write(')\n')
+    
+    def _writeSubsets(self, f, subsets, map, typeId):
+        '''
+        General function to write a subset to file as an OpenFOAM Map<meshSubset>.
+        '''
+        f.write('%d\n(\n' % len(subsets))
+        for name, ids in subsets.iteritems():
+            f.write('%s %s %d ( %s )\n' % (name, typeId, len(ids), ' '.join([str(map[id-1]) for id in ids])))
+        
+        f.write(')\n')
+    
+    def _writePointSubsets(self, f):
+        '''
+        Write the point subsets to file as and OpenFOAM Map<meshSubset>.
+        '''
+        self._writeSubsets(f, self._pointSubsets, self._nodeMap, '2')
+    
+    def _writeFaceSubsets(self, f):
+        '''
+        Write the face subsets to file as and OpenFOAM Map<meshSubset>.
+        '''
+        self._writeSubsets(f, self._faceSubsets, self._faceMap, '4')
+        
+    def _writeFeatureEdgeSubsets(self, f):
+        '''
+        Write the feature edge subsets to file as and OpenFOAM Map<meshSubset>.
+        '''
+        self._writeSubsets(f, self._featureEdgeSubsets, self._edgeMap, '8')
+    
+    def writeEdgeMesh(self, fileName):
+        '''
+        Write to file as an OpenFOAM edgeMesh
+        
+            fileName - The file name to write
+        '''
+        # Create file
+        f = open(fileName, 'wb') # NOTE: file opened as binary to ensure unix-style line breaks
+        
+        # Write header
+        f.write(foamHeader('edgeMesh', self._mesh.GetName()))
+        
+        self._writeNodes(f)
+        self._writeFeatureEdges(f)
+        
+        f.close()
+        
+        print 'edgeMesh written to %s' % fileName
+    
+    def writeFtr(self, fileName):
+        '''
+        Write to file as an OpenFOAM cfMesh FTR file
+        
+            fileName - the file name to write
+        '''
+        # Create file
+        f = open(fileName, 'wb') # NOTE: file opened as binary to ensure unix-style line breaks
+        
+        self._writePatchDefs(f)
+        self._writeNodes(f)
+        self._writeFacets(f)
+        
+        f.close()
+        
+        print 'triSurf written to %s' % fileName
+        
+    def writeFms(self, fileName):
+        '''
+        Write to file as an OpenFOAM cfMesh FMS file
+        
+            fileName - the file name to write
+        '''
+        # Create file
+        f = open(fileName, 'wb') # NOTE: file opened as binary to ensure unix-style line breaks
+        
+        self._writePatchDefs(f)
+        self._writeNodes(f)
+        self._writeFacets(f)
+        self._writeFeatureEdges(f)
+        self._writePointSubsets(f)
+        self._writeFaceSubsets(f)
+        self._writeFeatureEdgeSubsets(f)
+        
+        f.close()
+        
+        print 'triSurf written to %s' % fileName
+
+
+if __name__ == '__main__':
+    import salome
+    salome.salome_init()
+    
+    import SMESH, SALOMEDS

utilities/mergeSurfacePatches/Make/files

+mergeSurfacePatches.C
+
+EXE = $(FOAM_USER_APPBIN)/mergeSurfacePatches

utilities/mergeSurfacePatches/Make/options

+EXE_INC = \
+    -I../../meshLibrary/lnInclude \
+    -I$(FOAM_SRC)/meshTools/lnInclude \
+    -I$(FOAM_SRC)/triSurface/lnInclude
+
+EXE_LIBS = \
+    -ltriSurface \
+    -L$(FOAM_USER_LIBBIN) \
+    -lmeshLibrary \
+    -lmeshTools

utilities/mergeSurfacePatches/mergeSurfacePatches.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | cfMesh: A library for mesh generation
+   \\    /   O peration     |
+    \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
+     \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of cfMesh.
+
+    cfMesh is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    cfMesh is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with cfMesh.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    cfMesh utility to merge the supplied list of patches onto a single
+    patch.
+
+Author
+    Ivor Clifford <ivor.clifford@psi.ch>
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "autoPtr.H"
+#include "Time.H"
+#include "triSurf.H"
+#include "triSurfModifier.H"
+#include "demandDrivenData.H"
+#include "Pair.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Find the supplied list of patch names and return a list of patch Ids
+void getPatchIds
+(
+    const triSurf& origSurf,
+    const wordList& patchNames,
+    DynamicList<label>& patchIds
+)
+{
+    const geometricSurfacePatchList& origPatches = origSurf.patches();
+    
+    // Create patch name map
+    HashSet<word> patchNameHash(patchNames);
+        
+    // Find selected patches
+    label nFound = 0;
+    forAll(origPatches, patchI)
+    {
+        if (patchNameHash.found(origPatches[patchI].name()))
+        {
+            patchIds.append(patchI);
+            nFound++;
+        }
+    }
+    
+    if (nFound != patchNames.size())
+    {
+        WarningIn("getPatchIds")
+            << "Not all supplied patch names were found on the surface mesh" << endl;
+    }
+}
+
+
+// Copy all face subsets from one triSurf to another
+void copyFaceSubsets
+(
+    const triSurf& origSurf,
+    triSurf& newSurf
+)
+{
+    DynList<label> subsetIds;
+    origSurf.facetSubsetIndices(subsetIds);
+    
+    forAll(subsetIds, subsetI)
+    {
+        label newSubsetId = newSurf.addFacetSubset
+        (
+            origSurf.facetSubsetName(subsetI)
+        );
+        
+        labelList origFaces;        
+        origSurf.facetsInSubset(subsetI, origFaces);
+        
+        forAll(origFaces, faceI)
+        {
+            newSurf.addFacetToSubset
+            (