src/Tools/MacMesh/MacMesh/CutnGroup.py

   1 # Copyright (C) 2014-2016  EDF R&D
   2 #
   3 # This library is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU Lesser General Public
   5 # License as published by the Free Software Foundation; either
   6 # version 2.1 of the License, or (at your option) any later version.
   7 #
   8 # This library is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 # Lesser General Public License for more details.
  12 #
  13 # You should have received a copy of the GNU Lesser General Public
  14 # License along with this library; if not, write to the Free Software
  15 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 #
  17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 #
  19
  20
  21
  22 # This module allows cutting and grouping geometries by defining plane sections, with level of cutting as well as customizable Prefixes.
  23
  24 import math, Config
  25
  26 from salome.geom import geomBuilder
  27 geompy = geomBuilder.New()
  28
  29 def Go(GeoObj, CutPlnLst, OutLvlLst, PrefixLst, Publish):
  30
  31     """
  32     This function cuts any geometry (with infinite trim !) into several subgeometries that are cleanly saved inside the navigation tree. (Check GoTrim for the same functionality with custom trim size)
  33       - GeoObj is the geometrical object to be cut and grouped
  34       - CutPlnLst is a list of plane definitions. Each plane is a 6-tuple (contains 6 values). The first three represent the coordinates of the origin point and the second three represent the coordinates of the normal vector to the plane
  35             Example 1: [(0,0,0,1,0,0)]: cut along a plane passing through the origin and normal to the x-axis
  36             Example 2: [(0,0,0,1,0,0),(50,0,0,0,1,0)]: in addition to the first plane cut, cut through a plane passing by (50,0,0) and normal to the y axis.
  37             Note that the plane size us determined automatically from the size of the geometry in question (using a very big trim size = 100 x length of geometry!)
  38       - OutLvlLst is a list containing integers that represent the inner sectioning level with respect to the original geometry type
  39             A value of 1 means that the section will provide elements of one level lower than the original type. For example a solid sectioned at level 1 will produce faces. A Face sectioned at level 1 will produce edges.
  40             A value of 2 means that a deeper sectioning will be applied. A solid sectioned with level 2 will give faces and edges. A face will give edges and vertices. An edge will give only vertices
  41             The number of elements in this list should be (this is verified in the code) equal to the number of elements in the plane cut list. This is logical.
  42             Example 1: [1]
  43             Example 2: [1, 2], This means that the cut over the second plane will produce two types of elements unlike the first cut which will only output the first level objects.
  44       - PrefixLst is a list of strings that contains the naming Prefixes that are used by the script to generate the subshape names. This is very useful for relating the results to the sectioning requested.
  45             Example 1: ['Entry']
  46             Example 2: ['Entry','Exit']     The resulting groups from the sectioning with plane no.1 will then be saved as "Entry_FACE" and/or "Entry_EDGE" according to the original geometry object type and the cutting level
  47
  48     Imagine that we have a solid called ExampleSolid, an example command will be:
  49     CutnGroup.Go(ExampleSolid,[(0,0,0,1,0,0),(50,0,0,0,1,0)],[1, 2],['Entry','Exit'])
  50     """
  51
  52     NumCuts = CheckInput(CutPlnLst, OutLvlLst, PrefixLst, 1)
  53     OrigType = FindStandType(GeoObj,0)
  54     InvDictionary = {v: k for k, v in geompy.ShapeType.items()}  # Give geometry type name as a function of standard type numbering, ex: 4=FACE, 6=EDGE, 7=VERTEX
  55     TrimSize = geompy.BasicProperties(GeoObj)[0]*100
  56     CutPlane = [] ; Sections = [] ; Parts = []
  57
  58     if NumCuts:
  59         for i in range(0, NumCuts):             # Loop over the cutting planes to create them one by one
  60             CutPlane.append(CreatePlane(CutPlnLst[i],TrimSize))
  61         OutFather = geompy.MakePartition([GeoObj],CutPlane, [], [],FindStandType(GeoObj,1), 0, [], 0)   #Creating the partition object
  62         if Publish: geompy.addToStudy(OutFather,'SectionedObject')
  63         for i in range(0, NumCuts):
  64             for j in range(OrigType+1+2, OrigType+1+2*(OutLvlLst[i]+1),2):
  65                 if j == 8 : j = 7;      # Exception for the vertex case (=7)
  66                 PossSubShapesID = geompy.SubShapeAllIDs(OutFather,j)    # List of subshape IDs than correspond to the required cutting level (section type : face/wire/vertex)
  67                 PossSubShapes = geompy.ExtractShapes(OutFather,j)               # and the corresponding objects
  68                 Accepted = []
  69                 for k in range(0,len(PossSubShapesID)):         # Loop over all the subshapes checking if they belong to the cutting plane! if yes add them to current list
  70                     if  IsOnPlane(PossSubShapes[k], CutPlnLst[i], 1e-7):
  71                         Accepted.append(PossSubShapesID[k])
  72                 if Accepted :                                   # If some element is found, save it as a group with the prescribed Prefix
  73                     dummyObj = geompy.CreateGroup(OutFather, j)
  74                     geompy.UnionIDs(dummyObj, Accepted)
  75                     Sections.append(dummyObj)
  76                     if Publish:geompy.addToStudyInFather(OutFather, dummyObj, PrefixLst[i]+"_"+InvDictionary[j][0:2])
  77                 else :
  78                     print("Warning: For the section no.", i, ", No intersection of type " + InvDictionary[j] + " was found. Hence, no corresponding groups were created")
  79
  80         SubShapesID = geompy.SubShapeAllIDs(OutFather,OrigType+1)               # Saving also the groups corresponding to the sectioned item of the same type: ex. A solid into n sub-solids due to the sections
  81         for i in range(0,len(SubShapesID)):
  82             dummyObj = geompy.CreateGroup(OutFather, OrigType+1)
  83             geompy.UnionIDs(dummyObj, [SubShapesID[i]])
  84             Parts.append(dummyObj)
  85             if Publish: geompy.addToStudyInFather(OutFather, dummyObj, "SB"+"_"+InvDictionary[OrigType+1][0:3]+"_"+str(i+1))
  86
  87         return OutFather, Sections, Parts
  88     else:
  89         print("Fatal error, the routine cannot continue any further, check your input variables")
  90         return -1
  91
  92 def GoTrim(GeoObj, CutPlnLst, OutLvlLst, PrefixLst, Publish):
  93
  94     """
  95     This function cuts any geometry into several subgeometries that are cleanly saved inside the navigation tree with a fully customizable trim size.
  96       - GeoObj is the geometrical object to be cut and grouped
  97       - CutPlnLst is a list of plane definitions. Each plane is a 7-tuple (contains 7 values). The first three represent the coordinates of the origin point and the second three represent the coordinates of the normal vector to the plane, the last value corresponds to the trim size of the planes
  98             Example 1: [(0,0,0,1,0,0,5)]: cut along a plane passing through the origin and normal to the x-axis with a trim size of 5
  99             Example 2: [(0,0,0,1,0,0,5),(50,0,0,0,1,0,10)]: in addition to the first plane cut, cut through a plane passing by (50,0,0) and normal to the y axis with a trim size of 10
 100       - OutLvlLst is a list containing integers that represent the inner sectioning level with respect to the original geometry type
 101             A value of 1 means that the section will provide elements of one level lower than the original type. For example a solid sectioned at level 1 will produce faces. A Face sectioned at level 1 will produce edges.
 102             A value of 2 means that a deeper sectioning will be applied. A solid sectioned with level 2 will give faces and edges. A face will give edges and vertices. An edge will give only vertices
 103             The number of elements in this list should be (this is verified in the code) equal to the number of elements in the plane cut list. This is logical.
 104             Example 1: [1]
 105             Example 2: [1, 2], This means that the cut over the second plane will produce two types of elements unlike the first cut which will only output the first level objects.
 106       - PrefixLst is a list of strings that contains the naming Prefixes that are used by the script to generate the subshape names. This is very useful for relating the results to the sectioning requested.
 107             Example 1: ['Entry']
 108             Example 2: ['Entry','Exit']     The resulting groups from the sectioning with plane no.1 will then be saved as "Entry_FACE" and/or "Entry_EDGE" according to the original geometry object type and the cutting level
 109
 110     Imagine that we have a solid called ExampleSolid, an example command will be:
 111     CutnGroup.Go(ExampleSolid,[(0,0,0,1,0,0,5),(50,0,0,0,1,0,10)],[1, 2],['Entry','Exit'])
 112     """
 113
 114     NumCuts = CheckInput(CutPlnLst, OutLvlLst, PrefixLst, 0)
 115     OrigType = FindStandType(GeoObj,0)
 116     InvDictionary = {v: k for k, v in geompy.ShapeType.items()}  # Give geometry type name as a function of standard type numbering, ex: 4=FACE, 6=EDGE, 7=VERTEX
 117     CutPlane = [] ; Sections = [] ; Parts = []
 118     if NumCuts:
 119         for i in range(0, NumCuts):             # Loop over the cutting planes to create them one by one
 120             CutPlane.append(CreatePlane(CutPlnLst[i][0:6],CutPlnLst[i][6]))
 121         OutFather = geompy.MakePartition([GeoObj],CutPlane, [], [],FindStandType(GeoObj,1), 0, [], 0)   #Creating the partition object
 122         if Publish: geompy.addToStudy(OutFather,'SectionedObject')
 123         for i in range(0, NumCuts):
 124             for j in range(OrigType+1+2, OrigType+1+2*(OutLvlLst[i]+1),2):
 125                 if j == 8 : j = 7;      # Exception for the vertex case (=7)
 126                 PossSubShapesID = geompy.SubShapeAllIDs(OutFather,j)    # List of subshape IDs than correspond to the required cutting level (section type : face/wire/vertex)
 127                 PossSubShapes = geompy.ExtractShapes(OutFather,j)               # and the corresponding objects
 128                 Accepted = []
 129                 for k in range(0,len(PossSubShapesID)):         # Loop over all the subshapes checking if they belong to the cutting plane WITH THE TRIM SIZE CONDITION! if yes add them to current list
 130                     if  IsOnPlane(PossSubShapes[k], CutPlnLst[i], 1e-7) and Distance2Pt(geompy.PointCoordinates(geompy.MakeCDG(PossSubShapes[k])),CutPlnLst[i][0:3])<=CutPlnLst[i][-1]:
 131                         Accepted.append(PossSubShapesID[k])
 132                 if Accepted :                                   # If some element is found, save it as a group with the prescribed Prefix
 133                     dummyObj = geompy.CreateGroup(OutFather, j)
 134                     geompy.UnionIDs(dummyObj, Accepted)
 135                     Sections.append(dummyObj)
 136                     if Publish: geompy.addToStudyInFather(OutFather, dummyObj, PrefixLst[i]+"_"+InvDictionary[j][0:2])
 137                 else :
 138                     print("Warning: For the section no.", i, ", No intersection of type " + InvDictionary[j] + " was found. Hence, no corresponding groups were created")
 139
 140         SubShapesID = geompy.SubShapeAllIDs(OutFather,OrigType+1)               # Saving also the groups corresponding to the sectioned item of the same type: ex. A solid into n sub-solids due to the sections
 141         for i in range(0,len(SubShapesID)):
 142             dummyObj = geompy.CreateGroup(OutFather, OrigType+1)
 143             geompy.UnionIDs(dummyObj, [SubShapesID[i]])
 144             Parts.append(dummyObj)
 145             if Publish: geompy.addToStudyInFather(OutFather, dummyObj, "SB"+"_"+InvDictionary[OrigType+1][0:3]+"_"+str(i+1))
 146
 147         return OutFather, Sections, Parts
 148     else:
 149         print("Fatal error, the routine cannot continue any further, check your input variables")
 150         return -1
 151 def FindStandType(GeoObj, method):
 152     """
 153     Find the standard index for the Geometrical object/compound type input, see dictionary in geompy.ShapeType
 154     """
 155     TopType = GeoObj.GetMaxShapeType().__str__()
 156     UnModType = geompy.ShapeType[TopType]
 157     if method == 0 :
 158         StandType = UnModType-int(not(UnModType%2))             # So that wires and edges and considered the same, faces and shells, and so on
 159     else :
 160         StandType = UnModType
 161
 162     return(StandType)
 163
 164 def CreatePlane(CutPlnVar,Trim):
 165     """
 166     Creates a temporary point and vector in salome in order to build the sectioning planes needed
 167     """
 168     Temp_Vtx = geompy.MakeVertex(CutPlnVar[0], CutPlnVar[1], CutPlnVar[2])
 169     Temp_Vec = geompy.MakeVectorDXDYDZ(CutPlnVar[3], CutPlnVar[4], CutPlnVar[5])
 170     CutPlane = geompy.MakePlane(Temp_Vtx, Temp_Vec, Trim)
 171     return(CutPlane)
 172
 173 def CheckInput(CutPlnLst, OutLvlLst, PrefixLst, AutoTrim):
 174     """
 175     Checks the user input specifically if all needed parameters are provided
 176     """
 177     if not ((len(CutPlnLst) == len(OutLvlLst)) and (len(CutPlnLst) == len(PrefixLst))):
 178         print("Missing information about one or more of the cut planes")
 179         return 0
 180     elif not ((len(CutPlnLst[0]) == 6+int(not AutoTrim))):
 181         print("For each cutting plane you need to specify 6 parameters = 2 x 3 coordinates")
 182         return 0
 183     else:
 184         return len(CutPlnLst)
 185
 186 def IsOnPlane(GeoSubObj, CutPlnVar, tolerance):
 187     """
 188     Checks whether a geometry (vertex, segment, or face) belongs *completely* to the plane defined as a point and a normal vector
 189     """
 190     # lambda function that represents the plane equation, function = 0 <=> Pt defined with Coor belongs to plane
 191     PlaneEq = lambda Coor: CutPlnVar[3]*(Coor[0]-CutPlnVar[0])+CutPlnVar[4]*(Coor[1]-CutPlnVar[1])+CutPlnVar[5]*(Coor[2]-CutPlnVar[2])
 192     OrigType = FindStandType(GeoSubObj,0)
 193     if (OrigType >= 7):             # Vertex
 194         NonTrimDecision = abs(PlaneEq(geompy.PointCoordinates(GeoSubObj))) < tolerance
 195         if len(CutPlnVar) == 6 : return NonTrimDecision # No trim condition used
 196         else : return (NonTrimDecision and Distance2Pt(CutPlnVar[0:3],geompy.PointCoordinates(GeoSubObj))<=CutPlnVar[6]/2)
 197     elif (OrigType >= 5):           # Line, decompose into two points then call recursively IsOnPlane function!
 198         Verdict = True
 199         for i in range(0,2):
 200             Verdict = Verdict and IsOnPlane(geompy.GetVertexByIndex(GeoSubObj,i), CutPlnVar, tolerance)
 201         return Verdict
 202     elif (OrigType >= 3):           # Face, decompose into three points then call recursively IsOnPlane function!
 203         if IsOnPlane(geompy.MakeCDG(GeoSubObj),CutPlnVar, tolerance) : # Center of gravity belongs to plane, check if normal is parallel to plane
 204             NormalP1Coor = geompy.PointCoordinates(geompy.GetVertexByIndex(geompy.GetNormal(GeoSubObj),0))
 205             NormalP2Coor = geompy.PointCoordinates(geompy.GetVertexByIndex(geompy.GetNormal(GeoSubObj),1))
 206             Normal = [NormalP1Coor[0]-NormalP2Coor[0], NormalP1Coor[1]-NormalP2Coor[1], NormalP1Coor[2]-NormalP2Coor[2]]
 207             CrossP = CrossProd(CutPlnVar[3:6],Normal)               # Checks whether normals (of section plane and of face) are parallel or not
 208             if (abs(CrossP[0])<tolerance and abs(CrossP[1])<tolerance and abs(CrossP[2])<tolerance):        # meaning zero cross product => parallel
 209                 return True
 210             else :
 211                 return False
 212         else :
 213             return False
 214
 215
 216 def CrossProd(V1,V2):
 217     """
 218     Determines the cross product of two 3D vectors
 219     """
 220     return ([V1[1]*V2[2]-V1[2]*V2[1], V1[2]*V2[0]-V1[0]*V2[2], V1[0]*V2[1]-V1[1]*V2[0]])
 221
 222 def Distance2Pt(P1,P2):
 223     """
 224     Returns the distance between two points
 225     """
 226     return (math.sqrt((P1[0]-P2[0])**2+(P1[1]-P2[1])**2+(P1[2]-P2[2])**2))