doc/salome/examples/use_existing_faces.py

   1 # Usage of "Use Faces to be Created Manually" algorithm
   2
   3
   4 import salome
   5 salome.salome_init()
   6 import GEOM
   7 from salome.geom import geomBuilder
   8 geompy = geomBuilder.New(salome.myStudy)
   9
  10 import SMESH, SALOMEDS
  11 from salome.smesh import smeshBuilder
  12 smesh =  smeshBuilder.New(salome.myStudy)
  13 import salome_notebook
  14
  15 import numpy as np
  16
  17 # define my 2D algorithm
  18 def my2DMeshing( geomFace ):
  19
  20     # find gravity center of geomFace
  21     gcXYZ = geompy.PointCoordinates( geompy.MakeCDG( geomFace ))
  22
  23     # define order and orientation of edges
  24     sortedEdges = []
  25     geomEdges = geompy.SubShapeAll( geomFace, geompy.ShapeType["EDGE"])
  26     sortedEdges.append(( geomEdges.pop(0), True ))
  27     while geomEdges:
  28         prevEdge_rev = sortedEdges[ -1 ]
  29         prevVV = geompy.SubShapeAll( prevEdge_rev[0], geompy.ShapeType["VERTEX"])
  30         prevV2 = prevVV[ prevEdge_rev[1] ]
  31         found = False
  32         for iE in range( len( geomEdges )):
  33             v1,v2 = geompy.SubShapeAll( geomEdges[ iE ], geompy.ShapeType["VERTEX"])
  34             same1,same2 = [( geompy.MinDistance( prevV2, v ) < 1e-7 ) for v in [v1,v2] ]
  35             if not same1 and not same2: continue
  36             sortedEdges.append(( geomEdges.pop( iE ), same1 ))
  37             found = True
  38             break
  39         assert found
  40     sortedEdges.reverse()
  41
  42     # put nodes on edges in a right order
  43     nodes = []
  44     for edge, isForward in sortedEdges:
  45         v1,v2 = geompy.SubShapeAll( edge, geompy.ShapeType["VERTEX"])
  46         edgeNodes = mesh.GetSubMeshNodesId( v2,   all=False ) + \
  47                     mesh.GetSubMeshNodesId( edge, all=False ) + \
  48                     mesh.GetSubMeshNodesId( v1,   all=False )
  49         if not isForward: edgeNodes.reverse()
  50         nodes.extend( edgeNodes[:-1] )
  51
  52     # create nodes inside the geomFace
  53     r1 = 0.6
  54     r2 = 1 - r1
  55     nodesInside = []
  56     for n in nodes:
  57         nXYZ = mesh.GetNodeXYZ( n )
  58         newXYZ = np.add( np.multiply( r1, gcXYZ ), np.multiply( r2, nXYZ ))
  59         nodesInside.append( mesh.AddNode( newXYZ[0], newXYZ[1], newXYZ[2] ))
  60         mesh.SetNodeOnFace( nodesInside[-1], geomFace, 0, 0 )
  61
  62     # find out orientation of faces to create
  63     #    geomFace normal
  64     faceNorm = geompy.GetNormal( geomFace )
  65     v1,v2 = [ geompy.PointCoordinates( v ) \
  66               for v in geompy.SubShapeAll( faceNorm, geompy.ShapeType["VERTEX"]) ]
  67     faceNormXYZ = np.subtract( v2, v1 )
  68     outDirXYZ   = np.subtract( v1, [ 50, 50, 50 ] )
  69     if np.dot( faceNormXYZ, outDirXYZ ) < 0: # reversed face
  70         faceNormXYZ = np.multiply( -1., faceNormXYZ )
  71     #   mesh face normal
  72     e1 = np.subtract( mesh.GetNodeXYZ( nodes[0] ), mesh.GetNodeXYZ( nodes[1] ))
  73     e2 = np.subtract( mesh.GetNodeXYZ( nodes[0] ), mesh.GetNodeXYZ( nodesInside[0] ))
  74     meshNorm = np.cross( e1, e2 )
  75     #   faces orientation
  76     reverse = ( np.dot( faceNormXYZ, meshNorm ) < 0 )
  77
  78     # create mesh faces
  79     iN = len( nodes )
  80     while iN:
  81         n1, n2, n3, n4 = nodes[iN-1], nodes[iN-2], nodesInside[iN-2], nodesInside[iN-1]
  82         iN -= 1
  83         if reverse:
  84             f = mesh.AddFace( [n1, n2, n3, n4] )
  85         else:
  86             f = mesh.AddFace( [n4, n3, n2, n1] )
  87         # new faces must be assigned to geometry to allow 3D algorithm finding them
  88         mesh.SetMeshElementOnShape( f, geomFace )
  89
  90     if reverse:
  91         nodesInside.reverse()
  92     polygon = mesh.AddPolygonalFace( nodesInside )
  93     mesh.SetMeshElementOnShape( polygon, geomFace )
  94
  95     return
  96
  97 # create geometry and get faces to mesh with my2DMeshing()
  98 box = geompy.MakeBoxDXDYDZ( 100, 100, 100 )
  99 f1 = geompy.SubShapeAll( box, geompy.ShapeType["FACE"])[0]
 100 f2 = geompy.GetOppositeFace( box, f1 )
 101 geompy.addToStudy( box, "box" )
 102 geompy.addToStudy( f1, "f1" )
 103 geompy.addToStudy( f2, "f2" )
 104
 105 # compute 1D mesh
 106 mesh = smesh.Mesh( box )
 107 mesh.Segment().NumberOfSegments( 5 )
 108 mesh.Compute()
 109
 110 # compute 2D mesh
 111 mesh.Quadrangle()
 112 mesh.UseExistingFaces(f1) # UseExistingFaces() allows using my2DMeshing();
 113 mesh.UseExistingFaces(f2) # assign UseExistingFaces() BEFORE calling my2DMeshing()!
 114 my2DMeshing( f1 )
 115 my2DMeshing( f2 )
 116 assert mesh.Compute()
 117
 118 # compute 3D mesh
 119 mesh.Prism()
 120 assert mesh.Compute()