Commit d023ce70 authored by Franjo's avatar Franjo

Added contributions from Ivor Clifford

parent a6bc8074
#!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
# =============================================================================
# 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
mergeSurfacePatches.C
EXE = $(FOAM_USER_APPBIN)/mergeSurfacePatches
EXE_INC = \
-I../../meshLibrary/lnInclude \
-I$(FOAM_SRC)/meshTools/lnInclude \
-I$(FOAM_SRC)/triSurface/lnInclude
EXE_LIBS = \
-ltriSurface \
-L$(FOAM_USER_LIBBIN) \
-lmeshLibrary \
-lmeshTools
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
(