doc/examples/use_existing_faces.py

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