diff --git a/src/Tools/MacMesh/CMakeLists.txt b/src/Tools/MacMesh/CMakeLists.txt new file mode 100644 index 000000000..c50d20707 --- /dev/null +++ b/src/Tools/MacMesh/CMakeLists.txt @@ -0,0 +1,68 @@ +# Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE + +# Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +# CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. + +# This library 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 +# Lesser General Public License for more details. + +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com + +SET(MACMESH_INSTALL_PY ${SALOME_SMESH_INSTALL_PLUGINS}/MacMesh) +SET(MACMESH_INSTALL_DOC ${SALOME_INSTALL_DOC}/gui/SMESH/MacMesh) + +SALOME_CONFIGURE_FILE(Doc/index.html.in index.html) +SALOME_CONFIGURE_FILE(Example/PressureValve.py.in PressureValve.py) + +SET(plugin_DOC_FILES + Doc/MacMesh.v.10avr.pdf + ${CMAKE_CURRENT_BINARY_DIR}/index.html + Doc/snap.jpg + ) + +IF(SALOME_BUILD_DOC) + INSTALL(FILES ${plugin_DOC_FILES} DESTINATION ${MACMESH_INSTALL_DOC}) +ENDIF(SALOME_BUILD_DOC) + +# --- scripts --- + +# scripts / static +SET(plugin_SCRIPTS + MacMesh/Alarms.py + MacMesh/CentralUnrefine.py + MacMesh/CompositeBox.py + MacMesh/CompositeBoxF.py + MacMesh/Config.py + MacMesh/CutnGroup.py + MacMesh/Cylinder.py + MacMesh/GenFunctions.py + MacMesh/MacObject.py + MacMesh/PublishGroups.py + MacMesh/SharpAngle.py +) +SET(sample_SCRIPT + ${CMAKE_CURRENT_BINARY_DIR}/PressureValve.py +) +# --- rules --- + +SALOME_INSTALL_SCRIPTS("${plugin_SCRIPTS}" ${MACMESH_INSTALL_PY}) +SALOME_INSTALL_SCRIPTS("${sample_SCRIPT}" ${SALOME_INSTALL_SCRIPT_PYTHON} DEF_PERMS) + +SET(testname MacMesh_Example_PressureValve) +SALOME_GENERATE_TESTS_ENVIRONMENT(tests_env) +ADD_TEST( + NAME ${testname} + COMMAND ${PYTHON_EXECUTABLE} -B ${CMAKE_SOURCE_DIR}/doc/salome/examples/testme.py ${sample_SCRIPT}) +SET_TESTS_PROPERTIES(${testname} PROPERTIES ENVIRONMENT "${tests_env}") + diff --git a/src/Tools/MacMesh/Doc/MacMesh.v.10avr.pdf b/src/Tools/MacMesh/Doc/MacMesh.v.10avr.pdf new file mode 100755 index 000000000..0c59d56a6 Binary files /dev/null and b/src/Tools/MacMesh/Doc/MacMesh.v.10avr.pdf differ diff --git a/src/Tools/MacMesh/Doc/index.html.in b/src/Tools/MacMesh/Doc/index.html.in new file mode 100644 index 000000000..fc87fe214 --- /dev/null +++ b/src/Tools/MacMesh/Doc/index.html.in @@ -0,0 +1,20 @@ + + + + + MacMesh for Salome + +

The multi-purpose Salome plug-in for regular 2D meshing

+ + + + diff --git a/src/Tools/MacMesh/Doc/snap.jpg b/src/Tools/MacMesh/Doc/snap.jpg new file mode 100644 index 000000000..a9765fc59 Binary files /dev/null and b/src/Tools/MacMesh/Doc/snap.jpg differ diff --git a/src/Tools/MacMesh/Example/PressureValve.py.in b/src/Tools/MacMesh/Example/PressureValve.py.in new file mode 100755 index 000000000..e489faf57 --- /dev/null +++ b/src/Tools/MacMesh/Example/PressureValve.py.in @@ -0,0 +1,94 @@ +################################################################## +# Header for salome initialization ############################### + +import sys, salome, geompy, smesh, SMESH, math, os +sys.path.append( os.path.join( os.getenv('SMESH_ROOT_DIR'), '@MACMESH_INSTALL_PY@')) + +from MacObject import * +from SharpAngle import * +from CentralUnrefine import * +from PublishGroups import * +from CompositeBox import * +from CompositeBoxF import * + + +import Config, GenFunctions + +Config.theStudy = salome.myStudy; +geompy.init_geom(Config.theStudy) + +################################################################## +# Mesh name ###################################################### + +Config.StudyName = "SRV_X." + +################################################################## +# Definition of geometrical parameters ########################### + +X = 1.0 # Valve initial opening + +Config.StudyName += str(X)+"mm" + +R = 7.5 # Upstream pipe radius +T = 5.0 # Upstream pipe thickness +H = 20.0 # Upstream pipe height +J = 6.0 # Jet distance +E = 60.0 # Exit extent + +################################################################## +# Definition of meshing parameters ############################### + +d = 0.1 # Meshing element size at the inner corner +Nl = 1 # Number of levels in the local refinement +################################################################## +Bloc = [] + +# Object No. 1 # +Bloc.append( SharpAngleOut(0.,0.,X,1.5*X,X,d,'NE',Nl, + groups=['PH','PV_IN','VH',None,None,None]) ) + +# Object No. 2 # +Bloc.append( CompositeBox(X/2.+0.5*(R-X/2.),0.5*(X+X/2.)-X/2.,R-X/2.,X+X/2., + groups=[None,'VH',None,'AXIS'] ) ) + +# Object No. 3 # +Bloc.append( CompositeBoxF((0.,-X/2.),(R,-X/2.),(R,-H),(0.,-H), + groups=['IN',None,'PV_IN','AXIS'] ) ) + +# Object No. 4 # +Bloc.append( SharpAngleOut(-T,0.,X,X,X,d,'NW',Nl, + groups=['PH','PV_OUT',None,None,None,None]) ) + +# Object No. 5 # +Bloc.append( SharpAngleOut(-T,X,X,X,X,d,'SW',Nl, + groups=['VH','VV',None,None,None,None]) ) + +if X < T : + gap = T-X + Bloc.append( MacObject('CompBoxF',[(-X/2.-gap/2.,X/2.),(gap,X)], + ['auto'],groups=['PH','VH',None,None] ) ) + +# Object No. 6 # +Bloc.append( MacObject('CompBoxF',[(-T-X/2.-(J-X/2.)/2.,X/2.),(J-X/2.,2.*X)], + ['auto'],groups=[None,None,None,None] ) ) + +# Object No. 7 # +Bloc.append( CentralUnrefine(-T-J,X/2.,2.*E-J,E,'EW', + groups=[None,None,None,'OUT_V','OUT_H_HI','OUT_H_LO'])) + +# Object No. 8 # +Bloc.append( CompositeBox(-T-J/2.,-X/2.-0.5*((E-X)/2.-X/2.),J,(E-X)/2.-X/2., + groups=['OUT_H_LO',None,None,'PV_OUT'] ) ) + +# Object No. 9 # +Bloc.append( CompositeBox(-T-J/2.,X+X/2.+0.5*((E-X)/2.-X/2.),J,(E-X)/2.-X/2., + groups=[None,'OUT_H_HI',None,'VV'] ) ) + +SRVMesh = PublishGroups() + +RealLocalMeshing = Bloc[0][0].GeoPar[1][0]/Bloc[0][0].DirectionalMeshParams[0] +ExtrusionAngle = 2. * math.asin(RealLocalMeshing/(2*R))*180./math.pi +print "\nThe mesh will be revolved with an angle of :",ExtrusionAngle + +RevolveMesh(SRVMesh, Center=[R+0.01,0,0], Direction=[0,1,0], AngleDeg=ExtrusionAngle, Scale=0.001) + diff --git a/src/Tools/MacMesh/MacMesh/Alarms.py b/src/Tools/MacMesh/MacMesh/Alarms.py new file mode 100644 index 000000000..939b40e7c --- /dev/null +++ b/src/Tools/MacMesh/MacMesh/Alarms.py @@ -0,0 +1,10 @@ +def Message (code) : + import sys + MessageString = { 1 : lambda x: "Successfully created \n", + 2 : lambda x: "Fatal: Incorrect input \n", + 3 : lambda x: "Fatal: Overlapping objects detected \n", + 4 : lambda x: "Fatal: Incompatible object type with neighbouring objects" }[code](str(code)) + print MessageString + #if code > 1 : sys.exit() + return 1 + diff --git a/src/Tools/MacMesh/MacMesh/CentralUnrefine.py b/src/Tools/MacMesh/MacMesh/CentralUnrefine.py new file mode 100644 index 000000000..1daeb4dbe --- /dev/null +++ b/src/Tools/MacMesh/MacMesh/CentralUnrefine.py @@ -0,0 +1,166 @@ +# This object allows unrefining from a central point (actually, a line) to the exterior +# X0 and Y0 are the center points of the origin point and not those of the center of the generated block + + + +import sys, salome, geompy, smesh, SMESH, math, commands +CWD = commands.getoutput('pwd') +sys.path.append(CWD) + +from MacObject import * +import Config, GenFunctions + +def CentralUnrefine (X0 , Y0 , DX , DY , Orientation, **args ) : + + DirPar = {'SN' : lambda : ['NW', 'NE', 'EW', 'NW', 'SN', 'SN', 'NE', 'WE'], + 'NS' : lambda : ['SE', 'SW', 'WE', 'SE', 'NS', 'NS', 'SW', 'EW'], + 'EW' : lambda : ['NW', 'SW', 'SN', 'NW', 'EW', 'EW', 'SW', 'NS'], + 'WE' : lambda : ['SE', 'NE', 'NS', 'SE', 'WE', 'WE', 'NE', 'SN'], }[Orientation]() + + CoefVer = {'SN' : lambda : 1., + 'NS' : lambda : -1., + 'EW' : lambda : 0., + 'WE' : lambda : 0., }[Orientation]() + + CoefHor = {'SN' : lambda : 0., + 'NS' : lambda : 0., + 'EW' : lambda : -1., + 'WE' : lambda : 1., }[Orientation]() + + + MacObject('CompBoxF',[(X0+CoefHor*DX/2,Y0+CoefVer*DY/2),(DX,DY)],['auto'],publish=0) + ToLook1 = {'SN' : lambda : 2, + 'NS' : lambda : 3, + 'EW' : lambda : 1, + 'WE' : lambda : 0, }[Orientation]() + + ToLook2 = {'SN' : lambda : 0, + 'NS' : lambda : 0, + 'EW' : lambda : 2, + 'WE' : lambda : 2, }[Orientation]() + + ToLook3 = {'SN' : lambda : [0,1,2,3], + 'NS' : lambda : [1,0,3,2], + 'EW' : lambda : [3,2,1,0], + 'WE' : lambda : [2,3,0,1], }[Orientation]() + + if args.__contains__('groups') : + GroupNames = args['groups'] + else : GroupNames = [None, None, None, None, None, None] + + ExistingSegments = Config.ListObj[-1].DirectionalMeshParams[ToLook1] + ObjIDs = Config.Connections[-1][ToLook1] + RemoveLastObj() + + ExtensionSegments = math.ceil(ExistingSegments/12.)*12. + Dmin = 1.E50 + Dmax = -1.E50 + for ObjID in ObjIDs : + Boundaries = Config.ListObj[ObjID].Boundaries() + if Boundaries[ToLook2] < Dmin : Dmin = Boundaries[ToLook2] + if Boundaries[ToLook2+1] > Dmax : Dmax = Boundaries[ToLook2+1] + dx = 0 + if ExtensionSegments > ExistingSegments : + dn = (ExtensionSegments-ExistingSegments)/2. + dx = dn*(Dmax-Dmin)/ExistingSegments + #MacObject('CompBoxF',[(X0-CoefHor*dx/2+CoefVer*(-X0+Dmin-dx/2),Y0-CoefVer*dx/2+CoefHor*(-Y0+Dmin-dx/2)),(dx,dx)],[(dn,dn)],publish=0) + #MacObject('CompBoxF',[(X0-CoefHor*dx/2+CoefVer*(-X0+Dmax+dx/2),Y0-CoefVer*dx/2+CoefHor*(-Y0+Dmax+dx/2)),(dx,dx)],[(dn,dn)],publish=0) + + BoxSide = (Dmax-Dmin+2*dx)/2. + + Obj = [] + Obj.append(MacObject('BoxAng32',[(X0+CoefHor*(BoxSide/2)+CoefVer*(-BoxSide/2),Y0+CoefVer*(BoxSide/2)+CoefHor*(-BoxSide/2)),(BoxSide,BoxSide)],[int(ExtensionSegments/6),DirPar[0]],groups=GroupArray(ToLook3[0],GroupNames[0]))) + Obj.append(MacObject('BoxAng32',[(X0+CoefHor*(BoxSide/2)+CoefVer*(BoxSide/2),Y0+CoefVer*(BoxSide/2)+CoefHor*(BoxSide/2)),(BoxSide,BoxSide)],[int(ExtensionSegments/6),DirPar[1]],groups=GroupArray(ToLook3[0],GroupNames[0]))) + + NLevOpt = 0 + for NLevels in range (1,100) : + DX1 = abs(CoefVer)*BoxSide*2.**(NLevels+1)+abs(CoefHor)*BoxSide*2.**(NLevels) + DY1 = abs(CoefHor)*BoxSide*2.**(NLevels+1)+abs(CoefVer)*BoxSide*2.**(NLevels) + if DX1 > DX or DY1 > DY : + NLevOpt = NLevels-1 + DXinner = DX1/2. + DYinner = DY1/2. + break + + dummyArray = [DXinner,DYinner,DYinner,DXinner] + D1inner = dummyArray[ToLook2] # = DXinner for SN and NS orientations + D2inner = dummyArray[ToLook2+1] # = DYinner for SN and NS orientations + + dummyArray = [DX,DY,DY,DX] + D1 = dummyArray[ToLook2] # = DX for SN and NS orientations + D2 = dummyArray[ToLook2+1] # = DY for SN and NS orientations + + if D1inner < D1 : + GN0a = GroupArray(ToLook3[0],GroupNames[1]) + GN0b = GroupArray(ToLook3[0],GroupNames[2]) + GN01 = GroupArray(ToLook3[0],GroupNames[1]) + GN02 = GroupArray(ToLook3[0],GroupNames[2]) + if D2inner < D2 : + GN10 = [None,None,None,None] + GN11 = [None,None,None,None] + GN20 = [None,None,None,None] + else : + GN10 = GroupArray(ToLook3[1],GroupNames[3]) + GN11 = GroupArray(ToLook3[1],GroupNames[3]) + GN20 = GroupArray(ToLook3[1],GroupNames[3]) + else : + GN0a = GroupArray(ToLook3[0],GroupNames[1]) + GN0b = GroupArray(ToLook3[0],GroupNames[2]) + GN01 = GroupArray([ToLook3[0],ToLook3[2]],[GroupNames[1],GroupNames[4]]) + GN02 = GroupArray([ToLook3[0],ToLook3[3]],[GroupNames[2],GroupNames[5]]) + if D2inner < D2 : + GN10 = GroupArray(ToLook3[2],GroupNames[4]) + GN11 = GroupArray(ToLook3[3],GroupNames[5]) + GN20 = [None,None,None,None] + else : + GN10 = GroupArray([ToLook3[1],ToLook3[2]],[GroupNames[3],GroupNames[4]]) + GN11 = GroupArray([ToLook3[1],ToLook3[3]],[GroupNames[3],GroupNames[5]]) + GN20 = GroupArray(ToLook3[1],GroupNames[3]) + + for N in range (1,NLevOpt+1): + n=N-1 + D = BoxSide*(2.**n) + if N < NLevOpt : + Obj.append(MacObject('Box42' ,[(X0+D*(CoefHor*1/2-CoefVer*3/2) , Y0+D*(CoefVer*1/2-CoefHor*3/2) ) , (D,D)],['auto',DirPar[2]], groups=GN0a)) + Obj.append(MacObject('BoxAng32',[(X0+D*(CoefHor*3/2-CoefVer*3/2) , Y0+D*(CoefVer*3/2-CoefHor*3/2) ) , (D,D)],['auto',DirPar[3]])) + Obj.append(MacObject('Box42' ,[(X0+D*(CoefHor*3/2-CoefVer*1/2) , Y0+D*(CoefVer*3/2-CoefHor*1/2) ) , (D,D)],['auto',DirPar[4]])) + Obj.append(MacObject('Box42' ,[(X0+D*(CoefHor*3/2+CoefVer*1/2) , Y0+D*(CoefHor*1/2+CoefVer*3/2) ) , (D,D)],['auto',DirPar[5]])) + Obj.append(MacObject('BoxAng32',[(X0+D*(CoefVer*3/2+CoefHor*3/2) , Y0+D*(CoefVer*3/2+CoefHor*3/2) ) , (D,D)],['auto',DirPar[6]])) + Obj.append(MacObject('Box42' ,[(X0+D*(CoefVer*3/2+CoefHor*1/2) , Y0+D*(CoefHor*3/2+CoefVer*1/2) ) , (D,D)],['auto',DirPar[7]], groups=GN0b)) + else : + Obj.append(MacObject('Box42' ,[(X0+D*(CoefHor*1/2-CoefVer*3/2) , Y0+D*(CoefVer*1/2-CoefHor*3/2) ) , (D,D)],['auto',DirPar[2]], groups=GN01)) + Obj.append(MacObject('BoxAng32',[(X0+D*(CoefHor*3/2-CoefVer*3/2) , Y0+D*(CoefVer*3/2-CoefHor*3/2) ) , (D,D)],['auto',DirPar[3]], groups=GN10)) + Obj.append(MacObject('Box42' ,[(X0+D*(CoefHor*3/2-CoefVer*1/2) , Y0+D*(CoefVer*3/2-CoefHor*1/2) ) , (D,D)],['auto',DirPar[4]], groups=GN20)) + Obj.append(MacObject('Box42' ,[(X0+D*(CoefHor*3/2+CoefVer*1/2) , Y0+D*(CoefHor*1/2+CoefVer*3/2) ) , (D,D)],['auto',DirPar[5]], groups=GN20)) + Obj.append(MacObject('BoxAng32',[(X0+D*(CoefVer*3/2+CoefHor*3/2) , Y0+D*(CoefVer*3/2+CoefHor*3/2) ) , (D,D)],['auto',DirPar[6]], groups=GN11)) + Obj.append(MacObject('Box42' ,[(X0+D*(CoefVer*3/2+CoefHor*1/2) , Y0+D*(CoefHor*3/2+CoefVer*1/2) ) , (D,D)],['auto',DirPar[7]], groups=GN02)) + + + if CoefVer and DX>DXinner : + Obj.append(MacObject('CompBoxF',[(X0-CoefVer*0.25*(DX+DXinner),Y0+CoefVer*DYinner/2),((DX-DXinner)/2,DYinner)],['auto'], groups = GroupArray([ToLook3[0],ToLook3[2]],[GroupNames[1],GroupNames[4]]))) + Obj.append(MacObject('CompBoxF',[(X0+CoefVer*0.25*(DX+DXinner),Y0+CoefVer*DYinner/2),((DX-DXinner)/2,DYinner)],['auto'], groups = GroupArray([ToLook3[0],ToLook3[3]],[GroupNames[2],GroupNames[5]]))) + if DY>DYinner : + Obj.append(MacObject('CompBoxF',[(X0-CoefVer*0.25*(DX+DXinner),Y0+CoefVer*(DY+DYinner)/2.),((DX-DXinner)/2,DY-DYinner)],['auto'], groups = GroupArray([ToLook3[1],ToLook3[2]],[GroupNames[3],GroupNames[4]]))) + Obj.append(MacObject('CompBoxF',[(X0+CoefVer*0.25*(DX+DXinner),Y0+CoefVer*(DY+DYinner)/2.),((DX-DXinner)/2,DY-DYinner)],['auto'], groups = GroupArray([ToLook3[1],ToLook3[3]],[GroupNames[3],GroupNames[5]]))) + Obj.append(MacObject('CompBoxF',[(X0,Y0+CoefVer*(DY+DYinner)/2.),(DXinner,DY-DYinner)],['auto'], groups = GroupArray(ToLook3[1],GroupNames[3]))) + elif CoefHor and DY>DYinner : + Obj.append(MacObject('CompBoxF',[(X0+CoefHor*DXinner/2,Y0-CoefHor*0.25*(DY+DYinner)),(DXinner,(DY-DYinner)/2)],['auto'], groups = GroupArray([ToLook3[0],ToLook3[2]],[GroupNames[1],GroupNames[4]]))) + Obj.append(MacObject('CompBoxF',[(X0+CoefHor*DXinner/2,Y0+CoefHor*0.25*(DY+DYinner)),(DXinner,(DY-DYinner)/2)],['auto'], groups = GroupArray([ToLook3[0],ToLook3[3]],[GroupNames[2],GroupNames[5]]))) + if DX>DXinner : + Obj.append(MacObject('CompBoxF',[(X0+CoefHor*(DX+DXinner)/2.,Y0-CoefHor*0.25*(DY+DYinner)),(DX-DXinner,(DY-DYinner)/2)],['auto'], groups = GroupArray([ToLook3[1],ToLook3[2]],[GroupNames[3],GroupNames[4]]))) + Obj.append(MacObject('CompBoxF',[(X0+CoefHor*(DX+DXinner)/2.,Y0+CoefHor*0.25*(DY+DYinner)),(DX-DXinner,(DY-DYinner)/2)],['auto'], groups = GroupArray([ToLook3[1],ToLook3[3]],[GroupNames[3],GroupNames[5]]))) + Obj.append(MacObject('CompBoxF',[(X0+CoefHor*(DX+DXinner)/2.,Y0),(DX-DXinner,DYinner)],['auto'], groups = GroupArray(ToLook3[1],GroupNames[3]))) + return Obj + +def RemoveLastObj() : + Config.ListObj = Config.ListObj[:-1] + Config.Connections = Config.Connections[:-1] + +def GroupArray(indices, GroupNames) : + if type(indices) is int : + indices = [indices] + GroupNames = [GroupNames] + Output = [None,None,None,None] + for i, ind in enumerate(indices) : + Output[ind] = GroupNames[i] + return Output diff --git a/src/Tools/MacMesh/MacMesh/CompositeBox.py b/src/Tools/MacMesh/MacMesh/CompositeBox.py new file mode 100644 index 000000000..042f69c9c --- /dev/null +++ b/src/Tools/MacMesh/MacMesh/CompositeBox.py @@ -0,0 +1,165 @@ +# INTRODUCTION HERE + +import sys, salome, geompy, smesh, SMESH, math, copy, commands +CWD = commands.getoutput('pwd') +sys.path.append(CWD) + +from MacObject import * +import Config, GenFunctions + +def CompositeBox (X0 , Y0 , DX , DY , **args ) : + + if args.__contains__('groups') : + GroupNames = args['groups'] + else : GroupNames = [None, None, None, None] + # Create a full Box just to inherit, globally, the mesh parameters of bounding objects + MacObject('CompBoxF',[(X0,Y0),(DX,DY)],['auto'],publish=0) + + # Save the existing number of segments on each direction + ExistingSegments = Config.ListObj[-1].DirectionalMeshParams + + # Sort the connection list for the full Box + ObjIDLists = SortObjLists(Config.Connections[-1],X0 , Y0 , DX , DY ) + RemoveLastObj() + + print "ObjIDLists: ", ObjIDLists + + RealSegments = [] + Direction = [] + flag = 0 + if not(args.__contains__('recursive')) : Config.Count = 0 + print "Config.Count : ", Config.Count + Config.Criterion = GetCriterion(ObjIDLists) + for index, ObjList in enumerate(ObjIDLists) : + if not (ObjList[0] == -1 or Config.Count >= Config.Criterion): + if len(ObjList)>1 : flag = 1 + else : flag = 0 + for ObjID in ObjList: + ToLook0 = [2,2,0,0][index] + ToLook1 = [3,2,1,0][index] + CommonSide = FindCommonSide(Config.ListObj[ObjID].DirBoundaries(ToLook1),[X0-DX/2.,X0+DX/2.,Y0-DY/2.,Y0+DY/2.][ToLook0:ToLook0+2]) + ToLook2 = [1,0,3,2][index] + RealSegments.append(Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]*IntLen(CommonSide)/IntLen(Config.ListObj[ObjID].DirBoundaries(ToLook1))) + Direction.append(ToLook0/2) + + if flag and Config.Count < Config.Criterion: + if index < 2 : + if abs(CommonSide[0] - (Y0-DY/2.))<1e-7 : SouthGR = GroupNames[0] + else : SouthGR = None + if abs(CommonSide[1] - (Y0+DY/2.))<1e-7 : NorthGR = GroupNames[1] + else : NorthGR = None + CompositeBox (X0, CommonSide[0]+IntLen(CommonSide)/2., DX,IntLen(CommonSide), recursive=1, groups = [SouthGR,NorthGR]+GroupNames[2:4]) + else : + if abs(CommonSide[0] - (X0-DX/2.))<1e-7 : EastGR = GroupNames[2] + else : EastGR = None + if abs(CommonSide[1] - (X0+DX/2.))<1e-7 : WestGR = GroupNames[3] + else : WestGR = None + CompositeBox (CommonSide[0]+IntLen(CommonSide)/2., Y0, IntLen(CommonSide),DY, recursive=1, groups = GroupNames[0:2]+[EastGR,WestGR]) + + if Config.Count >= Config.Criterion : + break + if flag == 0 and Config.Count < Config.Criterion: + #print "Dir : ", Direction + #print "RealSegments : ", RealSegments + + #Xind = Direction.index(0) + #Yind = Direction.index(1) + #MacObject('CompBoxF',[(X0,Y0),(DX,DY)] ,[(RealSegments[Xind],RealSegments[Yind])], groups = GroupNames) + MacObject('CompBoxF',[(X0,Y0),(DX,DY)] ,['auto'], groups = GroupNames) + + Config.Count += 1 + + +def FindCommonSide (Int1, Int2) : + if abs(min(Int1[1],Int2[1])-max(Int1[0],Int2[0])) < 1e-5: return [0,0] + else : return [max(Int1[0],Int2[0]), min(Int1[1],Int2[1])] + +def IntLen (Interval) : + return abs(Interval[1]-Interval[0]) + +def RemoveLastObj() : + Config.ListObj = Config.ListObj[:-1] + Config.Connections = Config.Connections[:-1] + +def GetCriterion (ObjListIDs): + return max(Config.Criterion, max(len(ObjListIDs[0]),len(ObjListIDs[1]))*max(len(ObjListIDs[2]),len(ObjListIDs[3]))) + +def SortObjLists (List,X0,Y0,DX,DY) : + """ + This function sorts the list of neighbouring objects on each side, according to their intersection + with the object being created. From South to North and from East to West + """ + Output = List + # First find the directions where no neighbour exists + # Important : Here we assume that exactly two directions have no neighbours !!! + # Should we change this to allow a more general case ???? + dummy = IndexMultiOcc(List,(-1,)) + + # dummy[0] is either 0, meaning there is no neighbour on X- (West) + # or 1, meaning there is no neighbour on X+ (East) + # Similarly dummy[1] can be either 2 or 3 (South and North respectively) + # In order to get back to the formalism of groups (SNWE) + # => we do the following to define Sense of no neighbours and then the Direction list + # is calculated as to include uniquely the directions where we DO have neighbours + if len(dummy) == 1 : + # This adds a second direction where neighbours are not regarded, it is either 0 or 2 + dummy.append(2*(dummy[0]+2<4)) + print("Careful, you have neighbours on 3 or more sides of the box, we will not check if on two parallel sides the boxes are compatible !!!") + if len(dummy) == 2 or len(dummy) == 1 : + # Sense contains : Vertical then Horizontal + Sense = [dummy[1]%2,dummy[0]] + DirList = [[1,0][dummy[0]],[3,2][dummy[1]%2]] + for index,Direction in enumerate(DirList) : + ObjList = List[Direction] + RankMin = [] + ToLook0 = [2,2,0,0][Direction] + ToLook1 = [3,2,1,0][Direction] + for index1,ObjID in enumerate(ObjList) : + RankMin.append([-1.,1.][Sense[index]] * FindCommonSide(Config.ListObj[ObjID].DirBoundaries(ToLook1),[X0-DX/2.,X0+DX/2.,Y0-DY/2.,Y0+DY/2.][ToLook0:ToLook0+2])[Sense[index]]) + Output[Direction] = SortList(ObjList,RankMin) + + elif len(dummy) == 3 : + # We find the direction where we do have neighbours and then we sort the object list along it + Sense = dummy[0]%2 + Direction = [ i not in dummy for i in range(4) ].index(True) + ObjList = List[Direction] + RankMin = [] + ToLook0 = [2,2,0,0][Direction] + ToLook1 = [3,2,1,0][Direction] + for index1,ObjID in enumerate(ObjList) : + RankMin.append([-1.,1.][Sense] * FindCommonSide(Config.ListObj[ObjID].DirBoundaries(ToLook1),[X0-DX/2.,X0+DX/2.,Y0-DY/2.,Y0+DY/2.][ToLook0:ToLook0+2])[Sense]) + Output[Direction] = SortList(ObjList,RankMin) + else : + print ("Error : the composite box being created has no neighbours, how on earth do you want us to inherit its mesh parameters!!!") + + + return Output + +def IndexMultiOcc (Array,Element) : + """ + This functions returns the occurrences indices of Element in Array. + As opposed to Array.index(Element) method, this allows determining + multiple entries rather than just the first one! + """ + Output = [] + try : Array.index(Element) + except ValueError : print "No more occurrences" + else : Output.append(Array.index(Element)) + + if not(Output == []) and len(Array) > 1 : + for index, ArrElem in enumerate(Array[Output[0]+1:]) : + if ArrElem == Element : Output.append(index+Output[0]+1) + + return Output + +def SortList (ValList, CritList): + Output = [] + SortedCritList = copy.copy(CritList) + SortedCritList.sort() + for i in range(0,len(ValList)): + index = CritList.index(SortedCritList[i]) + Output.append(ValList[index]) + return Output + + + diff --git a/src/Tools/MacMesh/MacMesh/CompositeBoxF.py b/src/Tools/MacMesh/MacMesh/CompositeBoxF.py new file mode 100644 index 000000000..62c1414fa --- /dev/null +++ b/src/Tools/MacMesh/MacMesh/CompositeBoxF.py @@ -0,0 +1,230 @@ +# INTRODUCTION HERE + +import sys, salome, geompy, smesh, SMESH, math, copy, commands +CWD = commands.getoutput('pwd') +sys.path.append(CWD) + +from MacObject import * +import Config, GenFunctions + +def CompositeBoxF (Pt1 , Pt2 , Pt3 , Pt4 , **args ) : + [Pt1 , Pt2 , Pt3 , Pt4] = GenFunctions.SortPoints([Pt1 , Pt2 , Pt3 , Pt4]) + if args.__contains__('groups') : + GroupNames = args['groups'] + else : GroupNames = [None, None, None, None] + # Create a full NonOrtho box just to inherit, globally, the mesh parameters of bounding objects + dummy = MacObject('NonOrtho',[Pt1,Pt2,Pt3,Pt4],['auto'],publish=0) + # Save the existing number of segments on each direction + ExistingSeg0 = Config.ListObj[-1].DirectionalMeshParams + Convention = [2,3,0,1] + ExistingSegments = [ExistingSeg0[Convention[i]] for i in range(4)] + # Save Boundary lengths on each direction + BoundaryLengths = [IntLen(dummy.DirBoundaries(i)) for i in range(4) ] + # Calculate global mesh element size on each direction + GlobalDelta = [1.*BoundaryLengths[i]/ExistingSegments[i] for i in range(4) ] + print "GlobalDelta :",GlobalDelta + # Sort the connection list for the full Box + [(X0,Y0),(DX,DY)] = dummy.GeoPar + ObjIDLists = SortObjLists(Config.Connections[-1],X0 , Y0 , DX , DY ) + [Xmin,Xmax,Ymin,Ymax] = dummy.Boundaries() # Used for groups determination + RemoveLastObj() + + RealSegments = [] + Direction = [] + flag = 0 + if not(args.__contains__('recursive')) : + Config.Count = 0 + + Config.Criterion = GetCriterion(ObjIDLists) + for index, ObjList in enumerate(ObjIDLists) : + if not (ObjList[0] == -1 or Config.Count >= Config.Criterion): + if not(args.__contains__('recursive')) : + Config.DirIndex = index + if index > 1 : Config.RefPts = [Pt2, Pt3] + elif index == 0 : Config.RefPts = [Pt1, Pt2] + else : Config.RefPts = [Pt4, Pt3] + + if len(ObjList)>1 : flag = 1 + else : flag = 0 + for ObjID in ObjList: + ToLook0 = [2,3,0,1][index] + ToLook1 = [3,2,1,0][index] + CommonSide = FindCommonSide(Config.ListObj[ObjID].DirBoundaries(ToLook1),dummy.DirBoundaries(ToLook0)) + ToLook2 = [1,0,3,2][index] + RealSegments = Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]*IntLen(CommonSide)/IntLen(Config.ListObj[ObjID].DirBoundaries(ToLook1)) + LocalDelta = 1.*IntLen(CommonSide)/RealSegments + print "Direction:", ["West","East","South","North"][ToLook2] + print "IntLen(CommonSide):",IntLen(CommonSide) + print "RealSegments:",RealSegments + print "LocalDelta:",LocalDelta + if flag and Config.Count < Config.Criterion: + if index ==0 : + if abs(CommonSide[0] - Ymin)<1e-7 : SouthGR = GroupNames[0] + else : SouthGR = None + if abs(CommonSide[1] - Ymax)<1e-7 : NorthGR = GroupNames[1] + else : NorthGR = None + + NDelta = Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]* (LocalDelta-GlobalDelta[Convention[index]]) + [Pt1,Pt2] = Config.RefPts + Coef = [1.,-1.][index] + Vref1 = [Coef*(Pt2[0]-Pt1[0]),Coef*(Pt2[1]-Pt1[1])] + Vref2 = NormalizeVector([Pt2[0]-Pt3[0],Pt2[1]-Pt3[1]]) + Ptref = Config.ListObj[ObjID].PtCoor[[2,3][index]] + NewPt = ExtrapPoint (Ptref,Vref1,Vref2,NDelta) + CompositeBoxF (Pt1, Pt2, NewPt, Ptref, recursive=1, groups = [SouthGR,NorthGR]+GroupNames[2:4]) + elif index == 1: + if abs(CommonSide[0] - Ymin)<1e-7 : SouthGR = GroupNames[0] + else : SouthGR = None + if abs(CommonSide[1] - Ymax)<1e-7 : NorthGR = GroupNames[1] + else : NorthGR = None + + NDelta = Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]* (LocalDelta-GlobalDelta[Convention[index]]) + [Pt4,Pt3] = Config.RefPts + Coef = 1. + Vref1 = [Coef*(Pt4[0]-Pt3[0]),Coef*(Pt4[1]-Pt3[1])] + Vref2 = NormalizeVector([Pt1[0]-Pt4[0],Pt1[1]-Pt4[1]]) + Ptref = Config.ListObj[ObjID].PtCoor[0] + NewPt = ExtrapPoint (Ptref,Vref1,Vref2,NDelta) + CompositeBoxF (NewPt, Ptref, Pt3, Pt4, recursive=1, groups = [SouthGR,NorthGR]+GroupNames[2:4]) + else : + if abs(CommonSide[0] - Xmin)<1e-7 : WestGR = GroupNames[2] + else : WestGR = None + if abs(CommonSide[1] - Xmax)<1e-7 : EastGR = GroupNames[3] + else : EastGR = None + + NDelta = Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]* (LocalDelta-GlobalDelta[Convention[index]]) + [Pt2,Pt3] = Config.RefPts + Coef = [1.,-1.][index-2] + Vref1 = [Coef*(Pt3[0]-Pt2[0]),Coef*(Pt3[1]-Pt2[1])] + Vref2 = NormalizeVector([Pt3[0]-Pt4[0],Pt3[1]-Pt4[1]]) + Ptref = Config.ListObj[ObjID].PtCoor[[3,0][index-2]] + NewPt = ExtrapPoint (Ptref,Vref1,Vref2,NDelta) + CompositeBoxF (Ptref, Pt2, Pt3, NewPt, recursive=1, groups = GroupNames[0:2] + [WestGR,EastGR]) + + if Config.Count >= Config.Criterion : + break + if flag == 0 and Config.Count < Config.Criterion: + print "Creating NonOrtho object with the points:", Pt1,Pt2,Pt3,Pt4 + MacObject('NonOrtho',[Pt1,Pt2,Pt3,Pt4] ,['auto'], groups = GroupNames) + + Config.Count += 1 + if Config.DirIndex > 1 : Config.RefPts = [Pt1, Pt4] + elif Config.DirIndex==0 : Config.RefPts = [Pt4, Pt3] + else : Config.RefPts = [Pt1, Pt2] + +def FindCommonSide (Int1, Int2) : + if max(Int1[0],Int2[0]) we do the following to define Sense of no neighbours and then the Direction list + # is calculated as to include uniquely the directions where we DO have neighbours + if len(dummy) == 1 : + # This adds a second direction where neighbours are not regarded, it is either 0 or 2 + dummy.append(2*(dummy[0]+2<4)) + print("Careful, you have neighbours on 3 or more sides of the box, we will not check if on two parallel sides the boxes are compatible !!!") + if len(dummy) == 2 or len(dummy) == 1 : + # Sense contains : Vertical then Horizontal + Sense = [dummy[1]%2,dummy[0]] + DirList = [[1,0][dummy[0]],[3,2][dummy[1]%2]] + for index,Direction in enumerate(DirList) : + ObjList = List[Direction] + RankMin = [] + ToLook0 = [2,2,0,0][Direction] + ToLook1 = [3,2,1,0][Direction] + for index1,ObjID in enumerate(ObjList) : + RankMin.append([-1.,1.][Sense[index]] * FindCommonSide(Config.ListObj[ObjID].DirBoundaries(ToLook1),[X0-DX/2.,X0+DX/2.,Y0-DY/2.,Y0+DY/2.][ToLook0:ToLook0+2])[Sense[index]]) + Output[Direction] = SortList(ObjList,RankMin) + + elif len(dummy) == 3 : + # We find the direction where we do have neighbours and then we sort the object list along it + Sense = dummy[0]%2 + Direction = [ i not in dummy for i in range(4) ].index(True) + ObjList = List[Direction] + RankMin = [] + ToLook0 = [2,2,0,0][Direction] + ToLook1 = [3,2,1,0][Direction] + for index1,ObjID in enumerate(ObjList) : + RankMin.append([-1.,1.][Sense] * FindCommonSide(Config.ListObj[ObjID].DirBoundaries(ToLook1),[X0-DX/2.,X0+DX/2.,Y0-DY/2.,Y0+DY/2.][ToLook0:ToLook0+2])[Sense]) + Output[Direction] = SortList(ObjList,RankMin) + else : + print ("Error : the composite box being created has no neighbours, how on earth do you want us to inherit its mesh parameters!!!") + + + return Output + +def IndexMultiOcc (Array,Element) : + """ + This functions returns the occurrences indices of Element in Array. + As opposed to Array.index(Element) method, this allows determining + multiple entries rather than just the first one! + """ + Output = [] + try : Array.index(Element) + except ValueError : print "No more occurrences" + else : Output.append(Array.index(Element)) + + if not(Output == []) and len(Array) > 1 : + for index, ArrElem in enumerate(Array[Output[0]+1:]) : + if ArrElem == Element : Output.append(index+Output[0]+1) + + return Output + +def SortList (ValList, CritList): + Output = [] + SortedCritList = copy.copy(CritList) + SortedCritList.sort() + for i in range(0,len(ValList)): + index = CritList.index(SortedCritList[i]) + Output.append(ValList[index]) + return Output + +def ExtrapPoint (Ptref,Vref1,Vref2,Delta): + """ + This function allows determining the absolute coordinates of an extrapolation point + as shown in the following : + + + ExtrapPoint x---Vref2->--------o + / delta_glob |Vref1 + / | + Ptref x---------------------+ + delta_loc * Nseg + + Delta = (delta_loc - delta_glob) * Nseg + """ + + X = Ptref[0] + Vref1[0] + Delta*Vref2[0] + Y = Ptref[1] + Vref1[1] + Delta*Vref2[1] + return (X,Y,) + diff --git a/src/Tools/MacMesh/MacMesh/Config.py b/src/Tools/MacMesh/MacMesh/Config.py new file mode 100644 index 000000000..a86220e79 --- /dev/null +++ b/src/Tools/MacMesh/MacMesh/Config.py @@ -0,0 +1,29 @@ +dz = 1 +debug = 1 + +try : ListObj +except NameError : ListObj = [] + +try: Connections +except NameError: Connections = [] + +try : publish +except NameError : publish = 1 + +try : Groups +except NameError : Groups = [] + +try : StudyName +except NameError : StudyName = "Compound" + +try : Criterion +except NameError : Criterion = 1 + +try : Count +except NameError : Count = 0 + +try : RefPts +except NameError : RefPts = [] + +try : DirIndex +except NameError : DirIndex = 0 diff --git a/src/Tools/MacMesh/MacMesh/CutnGroup.py b/src/Tools/MacMesh/MacMesh/CutnGroup.py new file mode 100644 index 000000000..207862eb6 --- /dev/null +++ b/src/Tools/MacMesh/MacMesh/CutnGroup.py @@ -0,0 +1,201 @@ +# This module allows cutting and grouping geometries by defining plane sections, with level of cutting as well as customizable Prefixes. +import geompy, math + +def Go(GeoObj, CutPlnLst, OutLvlLst, PrefixLst, Publish): + + """ + 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) + - GeoObj is the geometrical object to be cut and grouped + - 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 + Example 1: [(0,0,0,1,0,0)]: cut along a plane passing through the origin and normal to the x-axis + 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. + 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!) + - OutLvlLst is a list containing integers that represent the inner sectioning level with respect to the original geometry type + A value of 1 means that the section will provide elements of one level lower than the original type. For exemple a solid sectioned at level 1 will produce faces. A Face sectioned at level 1 will produce edges. + 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 + 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. + Example 1: [1] + 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. + - 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. + Example 1: ['Entry'] + 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 + + Imagine that we have a solid called ExampleSolid, an example command will be: + CutnGroup.Go(ExampleSolid,[(0,0,0,1,0,0),(50,0,0,0,1,0)],[1, 2],['Entry','Exit']) + """ + + NumCuts = CheckInput(CutPlnLst, OutLvlLst, PrefixLst, 1) + OrigType = FindStandType(GeoObj,0) + InvDictionary = dict((v,k) for k, v in geompy.ShapeType.iteritems()) # Give geometry type name as a function of standard type numbering, ex: 4=FACE, 6=EDGE, 7=VERTEX + TrimSize = geompy.BasicProperties(GeoObj)[0]*100 + CutPlane = [] ; Sections = [] ; Parts = [] + + if NumCuts: + for i in range(0, NumCuts): # Loop over the cutting planes to create them one by one + CutPlane.append(CreatePlane(CutPlnLst[i],TrimSize)) + OutFather = geompy.MakePartition([GeoObj],CutPlane, [], [],FindStandType(GeoObj,1), 0, [], 0) #Creating the partition object + if Publish: geompy.addToStudy(OutFather,'SectionedObject') + for i in range(0, NumCuts): + for j in range(OrigType+1+2, OrigType+1+2*(OutLvlLst[i]+1),2): + if j == 8 : j = 7; # Exception for the vertex case (=7) + PossSubShapesID = geompy.SubShapeAllIDs(OutFather,j) # List of subshape IDs than correspond to the required cutting level (section type : face/wire/vertex) + PossSubShapes = geompy.ExtractShapes(OutFather,j) # and the corresponding objects + Accepted = [] + 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 + if IsOnPlane(PossSubShapes[k], CutPlnLst[i], 1e-7): + Accepted.append(PossSubShapesID[k]) + if Accepted : # If some element is found, save it as a group with the prescribed Prefix + dummyObj = geompy.CreateGroup(OutFather, j) + geompy.UnionIDs(dummyObj, Accepted) + Sections.append(dummyObj) + if Publish:geompy.addToStudyInFather(OutFather, dummyObj, PrefixLst[i]+"_"+InvDictionary[j][0:2]) + else : + print "Warning: For the section no.", i, ", No intersection of type " + InvDictionary[j] + " was found. Hence, no corresponding groups were created" + + 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 + for i in range(0,len(SubShapesID)): + dummyObj = geompy.CreateGroup(OutFather, OrigType+1) + geompy.UnionIDs(dummyObj, [SubShapesID[i]]) + Parts.append(dummyObj) + if Publish: geompy.addToStudyInFather(OutFather, dummyObj, "SB"+"_"+InvDictionary[OrigType+1][0:3]+"_"+str(i+1)) + + return OutFather, Sections, Parts + else: + print("Fatal error, the routine cannot continue any further, check your input variables") + return -1 + +def GoTrim(GeoObj, CutPlnLst, OutLvlLst, PrefixLst, Publish): + + """ + This function cuts any geometry into several subgeometries that are cleanly saved inside the navigation tree with a fully customizable trim size. + - GeoObj is the geometrical object to be cut and grouped + - 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 + 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 + 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 + - OutLvlLst is a list containing integers that represent the inner sectioning level with respect to the original geometry type + A value of 1 means that the section will provide elements of one level lower than the original type. For exemple a solid sectioned at level 1 will produce faces. A Face sectioned at level 1 will produce edges. + 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 + 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. + Example 1: [1] + 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. + - 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. + Example 1: ['Entry'] + 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 + + Imagine that we have a solid called ExampleSolid, an example command will be: + CutnGroup.Go(ExampleSolid,[(0,0,0,1,0,0,5),(50,0,0,0,1,0,10)],[1, 2],['Entry','Exit']) + """ + + NumCuts = CheckInput(CutPlnLst, OutLvlLst, PrefixLst, 0) + OrigType = FindStandType(GeoObj,0) + InvDictionary = dict((v,k) for k, v in geompy.ShapeType.iteritems()) # Give geometry type name as a function of standard type numbering, ex: 4=FACE, 6=EDGE, 7=VERTEX + CutPlane = [] ; Sections = [] ; Parts = [] + if NumCuts: + for i in range(0, NumCuts): # Loop over the cutting planes to create them one by one + CutPlane.append(CreatePlane(CutPlnLst[i][0:6],CutPlnLst[i][6])) + OutFather = geompy.MakePartition([GeoObj],CutPlane, [], [],FindStandType(GeoObj,1), 0, [], 0) #Creating the partition object + if Publish: geompy.addToStudy(OutFather,'SectionedObject') + for i in range(0, NumCuts): + for j in range(OrigType+1+2, OrigType+1+2*(OutLvlLst[i]+1),2): + if j == 8 : j = 7; # Exception for the vertex case (=7) + PossSubShapesID = geompy.SubShapeAllIDs(OutFather,j) # List of subshape IDs than correspond to the required cutting level (section type : face/wire/vertex) + PossSubShapes = geompy.ExtractShapes(OutFather,j) # and the corresponding objects + Accepted = [] + 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 + if IsOnPlane(PossSubShapes[k], CutPlnLst[i], 1e-7) and Distance2Pt(geompy.PointCoordinates(geompy.MakeCDG(PossSubShapes[k])),CutPlnLst[i][0:3])<=CutPlnLst[i][-1]: + Accepted.append(PossSubShapesID[k]) + if Accepted : # If some element is found, save it as a group with the prescribed Prefix + dummyObj = geompy.CreateGroup(OutFather, j) + geompy.UnionIDs(dummyObj, Accepted) + Sections.append(dummyObj) + if Publish: geompy.addToStudyInFather(OutFather, dummyObj, PrefixLst[i]+"_"+InvDictionary[j][0:2]) + else : + print "Warning: For the section no.", i, ", No intersection of type " + InvDictionary[j] + " was found. Hence, no corresponding groups were created" + + 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 + for i in range(0,len(SubShapesID)): + dummyObj = geompy.CreateGroup(OutFather, OrigType+1) + geompy.UnionIDs(dummyObj, [SubShapesID[i]]) + Parts.append(dummyObj) + if Publish: geompy.addToStudyInFather(OutFather, dummyObj, "SB"+"_"+InvDictionary[OrigType+1][0:3]+"_"+str(i+1)) + + return OutFather, Sections, Parts + else: + print("Fatal error, the routine cannot continue any further, check your input variables") + return -1 +def FindStandType(GeoObj, method): + """ + Find the standard index for the Geometrical object/compound type input, see dictionary in geompy.ShapeType + """ + TopType = GeoObj.GetMaxShapeType().__str__() + UnModType = geompy.ShapeType[TopType] + if method == 0 : + StandType = UnModType-int(not(UnModType%2)) # So that wires and edges and considered the same, faces and shells, and so on + else : + StandType = UnModType + + return(StandType) + +def CreatePlane(CutPlnVar,Trim): + """ + Creates a temporary point and vector in salome in order to build the sectioning planes needed + """ + Temp_Vtx = geompy.MakeVertex(CutPlnVar[0], CutPlnVar[1], CutPlnVar[2]) + Temp_Vec = geompy.MakeVectorDXDYDZ(CutPlnVar[3], CutPlnVar[4], CutPlnVar[5]) + CutPlane = geompy.MakePlane(Temp_Vtx, Temp_Vec, Trim) + return(CutPlane) + +def CheckInput(CutPlnLst, OutLvlLst, PrefixLst, AutoTrim): + """ + Checks the user input specifically if all needed parameters are provided + """ + if not ((len(CutPlnLst) == len(OutLvlLst)) and (len(CutPlnLst) == len(PrefixLst))): + print("Missing information about one or more of the cut planes") + return 0 + elif not ((len(CutPlnLst[0]) == 6+int(not AutoTrim))): + print("For each cutting plane you need to specify 6 parameters = 2 x 3 coordinates") + return 0 + else: + return len(CutPlnLst) + +def IsOnPlane(GeoSubObj, CutPlnVar, tolerance): + """ + Checks whether a geometry (vertex, segment, or face) belongs *completely* to the plane defined as a point and a normal vector + """ + # lambda function that represents the plane equation, function = 0 <=> Pt defined with Coor belongs to plane + PlaneEq = lambda Coor: CutPlnVar[3]*(Coor[0]-CutPlnVar[0])+CutPlnVar[4]*(Coor[1]-CutPlnVar[1])+CutPlnVar[5]*(Coor[2]-CutPlnVar[2]) + OrigType = FindStandType(GeoSubObj,0) + if (OrigType >= 7): # Vertex + NonTrimDecision = abs(PlaneEq(geompy.PointCoordinates(GeoSubObj))) < tolerance + if len(CutPlnVar) == 6 : return NonTrimDecision # No trim condition used + else : return (NonTrimDecision and Distance2Pt(CutPlnVar[0:3],geompy.PointCoordinates(GeoSubObj))<=CutPlnVar[6]/2) + elif (OrigType >= 5): # Line, decompose into two points then call recursively IsOnPlane function! + Verdict = True + for i in range(0,2): + Verdict = Verdict and IsOnPlane(geompy.GetVertexByIndex(GeoSubObj,i), CutPlnVar, tolerance) + return Verdict + elif (OrigType >= 3): # Face, decompose into three points then call recursively IsOnPlane function! + if IsOnPlane(geompy.MakeCDG(GeoSubObj),CutPlnVar, tolerance) : # Center of gravity belongs to plane, check if normal is parallel to plane + NormalP1Coor = geompy.PointCoordinates(geompy.GetVertexByIndex(geompy.GetNormal(GeoSubObj),0)) + NormalP2Coor = geompy.PointCoordinates(geompy.GetVertexByIndex(geompy.GetNormal(GeoSubObj),1)) + Normal = [NormalP1Coor[0]-NormalP2Coor[0], NormalP1Coor[1]-NormalP2Coor[1], NormalP1Coor[2]-NormalP2Coor[2]] + CrossP = CrossProd(CutPlnVar[3:6],Normal) # Checks whether normals (of section plane and of face) are parallel or not + if (abs(CrossP[0]) parallel + return True + else : + return False + else : + return False + + +def CrossProd(V1,V2): + """ + Determines the cross product of two 3D vectors + """ + 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]]) + +def Distance2Pt(P1,P2): + """ + Returns the distance between two points + """ + return (math.sqrt((P1[0]-P2[0])**2+(P1[1]-P2[1])**2+(P1[2]-P2[2])**2)) diff --git a/src/Tools/MacMesh/MacMesh/Cylinder.py b/src/Tools/MacMesh/MacMesh/Cylinder.py new file mode 100644 index 000000000..201cd6c7a --- /dev/null +++ b/src/Tools/MacMesh/MacMesh/Cylinder.py @@ -0,0 +1,94 @@ +# This is an automation of the cylinder-box object, defined with the coordinates of its center, its radius, and the box's +# boundary size. +# The pitch ratio is calculated automatically from the minimum of the box dimensions on x and y. +# This functions can take a groups input containing the group names of 4 sides in addition to the internal circular boundary +# in the following order : [South,North,West,East,Internal]. + +import sys, salome, geompy, smesh, SMESH, math, commands +CWD = commands.getoutput('pwd') +sys.path.append(CWD) + + +from MacObject import * +import Config, GenFunctions + +def Cylinder (X0 , Y0 , D , DX , DY , LocalMeshing , **args) : + if args.__contains__('DLocal') : DLocal = float(args['DLocal']) + else : DLocal = float(min(DX,DY)) + + # K is the pitch ratio + K = float(D)/(DLocal-D) + print "A local pitch ratio of K =", K ," will be used. " + NumCuts = 2*GenFunctions.QuarCylParam(K) + InternalMeshing = int(math.ceil(math.pi*D/(4*NumCuts*LocalMeshing))) + if InternalMeshing == 0 : InternalMeshing = 1 # This sets a minimum meshing condition in order to avoid an error. The user is notified of the value considered for the local meshing + print "Possible Local meshing is :", math.pi*D/(4*NumCuts*InternalMeshing), "\nThis value is returned by this function for your convenience.\n" + if args.__contains__('groups') : + GroupNames = args['groups'] + else : GroupNames = [None, None, None, None, None] + + if DY == DLocal : + if DX == DLocal: + GN1 = [None,GroupNames[1],None,GroupNames[3],GroupNames[4]] + GN2 = [None,GroupNames[1],GroupNames[2],None,GroupNames[4]] + GN3 = [GroupNames[0],None,GroupNames[2],None,GroupNames[4]] + GN4 = [GroupNames[0],None,None,GroupNames[3],GroupNames[4]] + else : + GN1 = [None,GroupNames[1],None,None,GroupNames[4]] + GN2 = [None,GroupNames[1],None,None,GroupNames[4]] + GN3 = [GroupNames[0],None,None,None,GroupNames[4]] + GN4 = [GroupNames[0],None,None,None,GroupNames[4]] + + GN5 = [GroupNames[0],GroupNames[1],None,GroupNames[3]] + GN6 = [GroupNames[0],GroupNames[1],GroupNames[2],None] + else : + if DX == DLocal: + GN1 = [None,None,None,GroupNames[3],GroupNames[4]] + GN2 = [None,None,GroupNames[2],None,GroupNames[4]] + GN3 = [None,None,GroupNames[2],None,GroupNames[4]] + GN4 = [None,None,None,GroupNames[3],GroupNames[4]] + GN7 = [GroupNames[0],None,GroupNames[2],GroupNames[3]] + GN8 = [None,GroupNames[1],GroupNames[2],GroupNames[3]] + else : + GN1 = [None,None,None,None,GroupNames[4]] + GN2 = [None,None,None,None,GroupNames[4]] + GN3 = [None,None,None,None,GroupNames[4]] + GN4 = [None,None,None,None,GroupNames[4]] + + GN5 = [None,None,None,GroupNames[3]] + GN6 = [None,None,GroupNames[2],None] + + GN9 = [GroupNames[0],None,None,GroupNames[3]] + GN10 = [GroupNames[0],None,None,None] + GN11 = [GroupNames[0],None,GroupNames[2],None] + + GN12 = [None,GroupNames[1],None,GroupNames[3]] + GN13 = [None,GroupNames[1],None,None] + GN14 = [None,GroupNames[1],GroupNames[2],None] + + Obj = [] + + Obj.append(MacObject('QuartCyl',[(X0+DLocal/4.,Y0+DLocal/4.),(DLocal/2.,DLocal/2.)],[InternalMeshing,'NE',K], groups = GN1)) + Obj.append(MacObject('QuartCyl',[(X0-DLocal/4.,Y0+DLocal/4.),(DLocal/2.,DLocal/2.)],['auto','NW',K], groups = GN2)) + Obj.append(MacObject('QuartCyl',[(X0-DLocal/4.,Y0-DLocal/4.),(DLocal/2.,DLocal/2.)],['auto','SW',K], groups = GN3)) + Obj.append(MacObject('QuartCyl',[(X0+DLocal/4.,Y0-DLocal/4.),(DLocal/2.,DLocal/2.)],['auto','SE',K], groups = GN4)) + + if DX > DLocal : + dX = (DX - DLocal)/2. + Obj.append(MacObject('CompBoxF',[(X0+DLocal/2.+dX/2.,Y0),(dX,DLocal)],['auto'], groups = GN5)) + Obj.append(MacObject('CompBoxF',[(X0-DLocal/2.-dX/2.,Y0),(dX,DLocal)],['auto'], groups = GN6)) + + if DY > DLocal : + dY = (DY - DLocal)/2. + if DX > DLocal : + Obj.append(MacObject('CompBoxF',[(X0+DLocal/2.+dX/2.,Y0-DLocal/2.-dY/2.),(dX,dY)],['auto'], groups = GN9)) + Obj.append(MacObject('CompBoxF',[(X0,Y0-DLocal/2.-dY/2.),(DLocal,dY)],['auto'], groups = GN10)) + Obj.append(MacObject('CompBoxF',[(X0-DLocal/2.-dX/2.,Y0-DLocal/2.-dY/2.),(dX,dY)],['auto'], groups = GN11)) + Obj.append(MacObject('CompBoxF',[(X0+DLocal/2.+dX/2.,Y0+DLocal/2.+dY/2.),(dX,dY)],['auto'], groups = GN12)) + Obj.append(MacObject('CompBoxF',[(X0,Y0+DLocal/2.+dY/2.),(DLocal,dY)],['auto'], groups = GN13)) + Obj.append(MacObject('CompBoxF',[(X0-DLocal/2.-dX/2.,Y0+DLocal/2.+dY/2.),(dX,dY)],['auto'], groups = GN14)) + else: + Obj.append(MacObject('CompBoxF',[(X0,Y0-DLocal/2.-dY/2.),(DLocal,dY)],['auto'], groups = GN7)) + Obj.append(MacObject('CompBoxF',[(X0,Y0+DLocal/2.+dY/2.),(DLocal,dY)],['auto'], groups = GN8)) + + return Obj diff --git a/src/Tools/MacMesh/MacMesh/GenFunctions.py b/src/Tools/MacMesh/MacMesh/GenFunctions.py new file mode 100644 index 000000000..12610d927 --- /dev/null +++ b/src/Tools/MacMesh/MacMesh/GenFunctions.py @@ -0,0 +1,844 @@ +# In this file are all the generation functions for manipulating the different created macro-objects +import geompy, smesh +import math, copy +import Config +import CutnGroup +import CompositeBox + +########################################################################################################## + +def Box11 (MacObject): + if Config.debug : print "Generating regular box" + + dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1) + RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0) + + MacObject.GeoChildren.append(RectFace) + MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1)) + + if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames) + + if Config.publish : + MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh + Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis + + EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary + Reg1D = MacObject.Mesh[0].Segment() + Reg1D.NumberOfSegments(MacObject.MeshPar[0]) + + MacObject.Mesh[0].Compute() # Generates the mesh + + MacObject.DirectionalMeshParams = [MacObject.MeshPar[0],MacObject.MeshPar[0],MacObject.MeshPar[0],MacObject.MeshPar[0]] + + MacObject.status = 1 + Config.ListObj.append(MacObject) + return MacObject + +########################################################################################################## + +def Box42 (MacObject): + if Config.debug : print "Generating box 4-2 reducer" + + Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.) + RotAngle = {'SN' : lambda : 0, + 'NS' : lambda : math.pi, + 'EW' : lambda : math.pi/2, + 'WE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]() + dummy0 = geompy.MakeRotation( ElemBox42 () , Z_Axis, RotAngle ) + dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1) + RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0) + + MacObject.GeoChildren.append(RectFace) + MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1)) + + if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames) + + if Config.publish : + MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh + Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis + + EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary + Reg1D = MacObject.Mesh[0].Segment() + Reg1D.NumberOfSegments(MacObject.MeshPar[0]) + + MacObject.Mesh[0].Compute() # Generates the mesh + + MacObject.status = 1 + + x = MacObject.MeshPar[0] + MacObject.DirectionalMeshParams = {'SN' : lambda : [3*x, 3*x, 4*x, 2*x], + 'NS' : lambda : [3*x, 3*x, 2*x, 4*x], + 'EW' : lambda : [2*x, 4*x, 3*x, 3*x], + 'WE' : lambda : [4*x, 2*x, 3*x, 3*x], }[MacObject.MeshPar[1]]() + + Config.ListObj.append(MacObject) + return MacObject + + +########################################################################################################## + +def BoxAng32 (MacObject): + if Config.debug : print "Generating sharp angle" + Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.) + RotAngle = {'NE' : lambda : 0, + 'NW' : lambda : math.pi/2, + 'SW' : lambda : math.pi, + 'SE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]() + dummy0 = geompy.MakeRotation( ElemEdge32 () , Z_Axis, RotAngle ) + dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1) + RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0) + + MacObject.GeoChildren.append(RectFace) + MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1)) + + if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames) + + if Config.publish : + MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh + Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis + + EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary + Reg1D = MacObject.Mesh[0].Segment() + Reg1D.NumberOfSegments(MacObject.MeshPar[0]) + + MacObject.Mesh[0].Compute() # Generates the mesh + + MacObject.status = 1 + + x = MacObject.MeshPar[0] + MacObject.DirectionalMeshParams = {'NE' : lambda : [3*x, 2*x, 3*x, 2*x], + 'NW' : lambda : [2*x, 3*x, 3*x, 2*x], + 'SW' : lambda : [2*x, 3*x, 2*x, 3*x], + 'SE' : lambda : [3*x, 2*x, 2*x, 3*x], }[MacObject.MeshPar[1]]() + + Config.ListObj.append(MacObject) + return MacObject +########################################################################################################## +def CompBox (MacObject): + if Config.debug : print "Generating composite box" + + dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1) + RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0) + + MacObject.GeoChildren.append(RectFace) + MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1)) + + if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames) + + if Config.publish : + MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh + Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis + + EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary + + ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1]) + + Reference = [0,0,0] + Vec = [(1,0,0),(0,1,0)] + for Edge in EdgeIDs: + for i in range(0,2): + if IsParallel(Edge,Vec[i]): + if not Reference[i]: # If this is the first found edge to be parallel to this direction, apply user preferences for meshing + Reference[i] = Edge + ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(ReducedRatio[i]*MacObject.MeshPar[0]))) + break + else: # If there already exists an edge parallel to this direction, then use a 1D projection + Apply1DProjMesh(MacObject.Mesh[0],Edge,Reference[i]) + break + + MacObject.Mesh[0].Compute() # Generates the mesh + + MacObject.DirectionalMeshParams = [MacObject.MeshPar[0]*ReducedRatio[1],MacObject.MeshPar[0]*ReducedRatio[1],MacObject.MeshPar[0]*ReducedRatio[0],MacObject.MeshPar[0]*ReducedRatio[0]] + + MacObject.status = 1 + Config.ListObj.append(MacObject) + return MacObject + +########################################################################################################## + +def CompBoxF (MacObject): + if Config.debug : print "Generating composite box" + + dummy1 = geompy.MakeScaleAlongAxes( ElemBox11 (), None , MacObject.GeoPar[1][0], MacObject.GeoPar[1][1], 1) + RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0) + + MacObject.GeoChildren.append(RectFace) + MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1)) + + if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames) + + if Config.publish : + MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh + Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis + + EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary + + #ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1]) + + Reference = [0,0,0] + Vec = [(1,0,0),(0,1,0)] + for Edge in EdgeIDs: + for i in range(0,2): + if IsParallel(Edge,Vec[i]): + if not Reference[i]: # If this is the first found edge to be parallel to this direction, apply user preferences for meshing + Reference[i] = Edge + ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(MacObject.MeshPar[0][i]))) + break + else: # If there already exists an edge parallel to this direction, then use a 1D projection + Apply1DProjMesh(MacObject.Mesh[0],Edge,Reference[i]) + break + + MacObject.Mesh[0].Compute() # Generates the mesh + + MacObject.DirectionalMeshParams = [MacObject.MeshPar[0][1],MacObject.MeshPar[0][1],MacObject.MeshPar[0][0],MacObject.MeshPar[0][0]] + + MacObject.status = 1 + Config.ListObj.append(MacObject) + return MacObject +########################################################################################################## + +def NonOrtho (MacObject): + if Config.debug : print "Generating Non-orthogonal quadrangle" + + RectFace = Quadrangler (MacObject.PtCoor) + + MacObject.GeoChildren.append(RectFace) + MacObject.GeoChildrenNames.append("Quad_"+ str(len(Config.ListObj)+1)) + + + if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames) + + if Config.publish : + MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh + Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis + + EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary + + #ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0],MacObject.GeoPar[1][1]) + + Vec = [MacObject.DirVectors(i) for i in range(4)] + for Edge in EdgeIDs: + Dir = [IsParallel(Edge,Vec[j]) for j in range(4)].index(True) + DirConv = [0,0,1,1][Dir] + ApplyConstant1DMesh(MacObject.Mesh[0],Edge,int(round(MacObject.MeshPar[0][DirConv]))) + + MacObject.Mesh[0].Compute() # Generates the mesh + + MacObject.DirectionalMeshParams = [MacObject.MeshPar[0][1],MacObject.MeshPar[0][1],MacObject.MeshPar[0][0],MacObject.MeshPar[0][0]] + + MacObject.status = 1 + Config.ListObj.append(MacObject) + return MacObject + +########################################################################################################## + +def QuartCyl (MacObject): + if Config.debug : print "Generating quarter cylinder" + Z_Axis = geompy.MakeVectorDXDYDZ(0., 0., 1.) + RotAngle = {'NE' : lambda : 0, + 'NW' : lambda : math.pi/2, + 'SW' : lambda : math.pi, + 'SE' : lambda : -math.pi/2, }[MacObject.MeshPar[1]]() + dummy0 = geompy.MakeRotation( ElemQuartCyl(MacObject.MeshPar[2]) , Z_Axis, RotAngle ) + dummy1 = geompy.MakeScaleAlongAxes( dummy0, None , MacObject.GeoPar[1][0]/10., MacObject.GeoPar[1][1]/10., 1) + RectFace = geompy.MakeTranslation(dummy1, MacObject.GeoPar[0][0], MacObject.GeoPar[0][1], 0) + + MacObject.GeoChildren.append(RectFace) + MacObject.GeoChildrenNames.append("Box_"+ str(len(Config.ListObj)+1)) + + if Config.debug : Publish (MacObject.GeoChildren,MacObject.GeoChildrenNames) + + if Config.publish : + MacObject.Mesh.append(smesh.Mesh(RectFace)) # Creation of a new mesh + Quad2D = MacObject.Mesh[0].Quadrangle() # Applying a quadrangle hypothesis + + EdgeIDs = geompy.SubShapeAllSorted(RectFace,6) # List of Edge IDs belonging to RectFace, 6 = Edge in salome dictionary + Reg1D = MacObject.Mesh[0].Segment() + + #if MacObject.MeshPar[0] == 2 and MacObject.MeshPar[2] <= 2.: + # print("Due to a bug in Salome 6.3, we are forced to either increase or decrease the local refinement by 50%, we choose in this case to increase the model's refinement.") + # MacObject.MeshPar[0] = 3 + + Reg1D.NumberOfSegments(MacObject.MeshPar[0]) + + MacObject.Mesh[0].Compute() # Generates the mesh + + MacObject.status = 1 + + x = MacObject.MeshPar[0] + N = QuarCylParam(MacObject.MeshPar[2])+1 + + MacObject.DirectionalMeshParams = {'NE' : lambda : [2*x, N*x, 2*x, N*x], + 'NW' : lambda : [N*x, 2*x, 2*x, N*x], + 'SW' : lambda : [N*x, 2*x, N*x, 2*x], + 'SE' : lambda : [2*x, N*x, N*x, 2*x], }[MacObject.MeshPar[1]]() + + Config.ListObj.append(MacObject) + return MacObject + +########################################################################################################## +# Below this are the elementary calculation/visualization functions +########################################################################################################## + +def Publish (ObjToPublish,NamesToPublish): + i = 0 + for GeoObj in ObjToPublish : + geompy.addToStudy(GeoObj,NamesToPublish[i]) + i = i+1 + +def IsParallel (Edge, Vector): + """ + Function checks whether a given edge object is parallel to a reference vector. + Output can be 0 (not parallel) or 1 (parallel and same sense) or 2 (parallel and opposite sense). + If the reference vector is null, the function returns 0 + """ + if Vector == (0,0,0) : return 0 + else : + P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0)) + P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1)) + V0 = [ P1[0] - P2[0], P1[1] - P2[1], P1[2] - P2[2] ] + if Distance2Pt((0,0,0),CrossProd(V0,Vector))<1e-7 and DotProd(V0,Vector) > 0 : return 1 + elif Distance2Pt((0,0,0),CrossProd(V0,Vector))<1e-7 and DotProd(V0,Vector) < 0 : return 2 + else : return 0 + +def IsOnCircle (Edge, Center, Radius): + """ + Function checks whether a given edge object belong to the periphery of a circle defined by its + center and radius. + Output can be 0 (does not belong) or 1 (belongs). + If the reference Radius is null, the function returns 0 + Note that this function is basic in the sense that it only checks if the two border points of a + given edge belong to the arc of reference. + """ + if Radius == 0 : return 0 + else : + P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0)) + P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1)) + if abs(Distance2Pt(Center,P1)-Radius) < 1e-6 and abs(Distance2Pt(Center,P2)-Radius) < 1e-6: + return 1 + else : + return 0 + +def CrossProd(V1,V2): + """ + Determines the cross product of two 3D vectors + """ + 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]]) + +def QuarCylParam(PitchRatio): + R = float(PitchRatio)/(PitchRatio+1) + Eps = 1. - R + X = (R+Eps/2.)*math.sin(math.pi/4)+Eps/2. + N = int(math.floor((math.pi*R/4.)/(Eps/2.))) + return N + +def DotProd(V1,V2): + """ + Determines the dot product of two 3D vectors + """ + if len(V1)==2 : V1.append(0) + if len(V2)==2 : V2.append(0) + + return (V1[0]*V2[0]+V1[1]*V2[1]+V1[2]*V2[2]) + +def Distance2Pt(P1,P2): + """ + Returns the distance between two points + """ + return (math.sqrt((P1[0]-P2[0])**2+(P1[1]-P2[1])**2+(P1[2]-P2[2])**2)) + +def ApplyConstant1DMesh (ParentMsh, Edge, Nseg): + Reg1D = ParentMsh.Segment(geom=Edge) + Len = Reg1D.NumberOfSegments(Nseg) + +def Apply1DProjMesh (ParentMsh, Edge, Ref): + Proj1D = ParentMsh.Projection1D(geom=Edge) + SrcEdge = Proj1D.SourceEdge(Ref,None,None,None) + +def EdgeLength (Edge): + """ + This function returns the edge object length. + """ + P1 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0)) + P2 = geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1)) + return Distance2Pt(P1,P2) + + +def D2R (Angle): + return Angle*math.pi/180 + +def R2D (Angle): + return Angle*180/math.pi + +def F2D (FloatNumber): + return round(FloatNumber*100.)/100. + +def BezierGen (PointA, PointB, AngleA, AngleB): + + if AngleA == 0 and AngleB == 0 : return (geompy.MakeEdge(PointA, PointB)) + else : + A = geompy.PointCoordinates(PointA) + B = geompy.PointCoordinates(PointB) + dAB = Distance2Pt(A,B) + dAC = dAB * (math.tan(AngleA)*math.tan(AngleB)) / (math.sin(AngleA) * ( math.tan(AngleA)+math.tan(AngleB) ) ) + AngleOX_AB = math.acos((B[0]-A[0])/dAB) + PointC = geompy.MakeVertex(A[0]+math.cos(AngleA+AngleOX_AB)*dAC,A[1]+math.sin(AngleA+AngleOX_AB)*dAC,0) + CurveACB = geompy.MakeBezier([PointA,PointC,PointB]) + return CurveACB + +def GetSideAngleForBezier (PointA , PointB): + """ + This function takes for input two points A and B where the bezier line is needed. It calculates the incident + angle needed at point A so that the final curve is either at 0 or 90 degrees from the x'Ox axis + """ + A = geompy.PointCoordinates(PointA) + B = geompy.PointCoordinates(PointB) + ABx = B[0]-A[0] + dAB = Distance2Pt(A,B) + Alpha = math.acos(ABx/dAB) + #print "New angle request" + #print ABx, dAB, R2D(Alpha) + if Alpha < math.pi/4 : + #print "returning", R2D(-Alpha) + return -Alpha + elif Alpha < 3*math.pi/4 : + #print "returning", R2D(-(Alpha-math.pi/2)) + return -(Alpha-math.pi/2) + else : + #print "returning", R2D(-(Alpha-math.pi)) + return -(Alpha-math.pi) + +def VecDivRatio (Vec1, Vec2): + """ + This function tries to find the ratio of Vec1 on Vec2 while neglecting any zero term in Vec1. This is used afterwards + for determining the global mesh parameter from automatically detected directional mesh params. If no compatibility is + possible, the function returns -1 + """ + Vec3 = [] + for i in range(len(Vec1)) : + Vec3.append(float(Vec1[i])/Vec2[i]) + Ratio=[] + for i in Vec3 : + if not (abs(i)<1e-7) : Ratio.append(i) + if Ratio : + if min(Ratio) == max(Ratio) and min(Ratio)==int(min(Ratio)) : return(min(Ratio)) + else : return -1 + else : + return -2 + + +def ReduceRatio (dx, dy): + """ + This function transforms a decimal ratio into a scale between two integers, for example : [0.2,0.05] --> [4,1] ; + """ + Output = [0,0] + ratio = float(dy)/dx + if isinteger(ratio) : return [1,ratio] + elif dx == 1 : # when this function is called recursively! + for i in range(1,20) : # searches over 20 decimals + if isinteger(ratio * (10**i) ) : + Output = GetScale((10**i),int(round(ratio * (10**i) ) ) ) + break + else : + for n in range(0,i) : + if isinteger(ratio * ( 10**(i)-10**(n) )) : + Output = GetScale( 10**(i)-10**(n) , int(round(ratio * ( 10**(i)-10**(n) ) ) ) ) + break + if not (Output==[0,0]) : break + return Output + else : + for i in range(1,10) : # searches over 10 decimals + if isinteger(ratio * (10**i) ) : + Output = GetScale((10**i),int(round(ratio * (10**i) ) ) ) + break + else : + for n in range(0,i) : + if isinteger(ratio * ( 10**(i)-10**(n) )) : + Output = GetScale( 10**(i)-10**(n) , int(round(ratio * ( 10**(i)-10**(n) ) ) ) ) + break + if not (Output==[0,0]) : break + + if Output == [0,0] : + print "We are having some trouble while interpreting the following ratio: ",ratio, "\nWe will try a recursive method which may in some cases take some time..." + if dy > dx : + A = ReduceRatio (dx, dy-dx) + return ([A[0],A[1]+A[0]]) + else : + A = ReduceRatio (dy, dx-dy) + return ([A[1]+A[0],A[0]]) + + else : return Output + +def GetScale (X,Y): + """ + This function is called within ReduceRatio and aims to reduce down two integers X and Y by dividing them with their common divisors; + Example: 25 and 5 ---> 5 and 1 / 63 and 12 ---> 21 and 4 + """ + MaxDiv = max(X,Y) + Divisor = 2 # Initializing the divisor + while MaxDiv >= Divisor : + X0 = 0 + Y0 = 0 + if not(X%Divisor) : + X0 = X/Divisor + MaxDiv = max(MaxDiv,X0) + if not(Y%Divisor) : + Y0 = Y/Divisor + MaxDiv = max(MaxDiv,Y0) + if (X0*Y0) : + X = X0 + Y = Y0 + else : + Divisor = Divisor + 1 + return [X,Y] + +def isinteger (x) : + """ + This functions applies a simple check if the entered value is an integer + """ + x = float('%.5f' % (x)) #Truncate x to 5 digits after the decimal point + if math.ceil(x) == math.floor(x) : return True + else : return False +########################################################################################## +# Below this are the functions that create the elementary forms for the macro objects +########################################################################################## + +def ElemBox11 (): + """ + This function returns a simple square face of 1 side length + """ + RectFace = geompy.MakeFaceHW(1, 1, 1) + return RectFace + +def ElemBox42 (): + """ + This function returns a square face of 1 side length, partitioned + according to the elementary 4 to 2 reductor method + """ + OrigRectFace = geompy.MakeFaceHW(1, 1, 1) + + SouthPt1 = geompy.MakeVertex (-.25, -.5, 0) + SouthPt2 = geompy.MakeVertex (0, -.5, 0) + SouthPt3 = geompy.MakeVertex (.25, -.5, 0) + WestPt1 = geompy.MakeVertex (-.5, -.5+1./3, 0) + WestPt2 = geompy.MakeVertex (-.5, -.5+2./3, 0) + EastPt1 = geompy.MakeVertex (.5, -.5+1./3, 0) + EastPt2 = geompy.MakeVertex (.5, -.5+2./3, 0) + NorthPt = geompy.MakeVertex (0, .5, 0) + MidPt1 = geompy.MakeVertex (0, .05, 0) + MidPt2 = geompy.MakeVertex (.2, -.18, 0) + MidPt3 = geompy.MakeVertex (0, -.28, 0) + MidPt4 = geompy.MakeVertex (-.2, -.18, 0) + + Cutter = [] + Cutter.append(geompy.MakeEdge(SouthPt2, MidPt3)) + Cutter.append(geompy.MakeEdge(MidPt1, NorthPt)) + Cutter.append(BezierGen(SouthPt1, MidPt4, GetSideAngleForBezier(SouthPt1,MidPt4), D2R(15))) + Cutter.append(BezierGen(SouthPt3, MidPt2, GetSideAngleForBezier(SouthPt3,MidPt2), D2R(-15))) + Cutter.append(BezierGen(WestPt1, MidPt4, GetSideAngleForBezier(WestPt1,MidPt4), D2R(-10))) + Cutter.append(BezierGen(EastPt1, MidPt2, GetSideAngleForBezier(EastPt1,MidPt2), D2R(10))) + Cutter.append(BezierGen(WestPt2, MidPt1, GetSideAngleForBezier(WestPt2,MidPt1), D2R(-10))) + Cutter.append(BezierGen(EastPt2, MidPt1, GetSideAngleForBezier(EastPt2,MidPt1), D2R(10))) + Cutter.append(BezierGen(MidPt2, MidPt1, D2R(-15), D2R(-15))) + Cutter.append(BezierGen(MidPt3, MidPt2, D2R(10), D2R(15))) + Cutter.append(BezierGen(MidPt3, MidPt4, D2R(-10), D2R(-15))) + Cutter.append(BezierGen(MidPt4, MidPt1, D2R(15), D2R(15))) + + RectFace = geompy.MakePartition([OrigRectFace],Cutter, [], [],4, 0, [], 0) #Creating the partition object + #i=1 + #for SingleCut in Cutter : + # geompy.addToStudy(SingleCut,'Cutter'+str(i)) + # i = i+1 + #geompy.addToStudy(RectFace,'RectFace') + return RectFace + +def ElemEdge32 (): + """ + This function returns a square face of 1 side length, partitioned + according to the elementary edge with 3 to 2 reductor + """ + OrigRectFace = geompy.MakeFaceHW(1., 1., 1) + + SouthPt1 = geompy.MakeVertex (-1./6, -0.5, 0.) + SouthPt2 = geompy.MakeVertex ( 1./6, -0.5, 0.) + WestPt1 = geompy.MakeVertex (-0.5, -1./6, 0.) + WestPt2 = geompy.MakeVertex (-0.5, 1./6, 0.) + EastPt = geompy.MakeVertex ( 0.5, 0., 0.) + NorthPt = geompy.MakeVertex (0., 0.5, 0.) + + MidPt1 = geompy.MakeVertex (-0.2, -0.2, 0.) + MidPt2 = geompy.MakeVertex ( -0.02, -0.02, 0.) + + Cutter = [] + Cutter.append(BezierGen(SouthPt1, MidPt1, GetSideAngleForBezier(SouthPt1,MidPt1) , D2R(-5))) + Cutter.append(BezierGen( WestPt1, MidPt1, GetSideAngleForBezier(WestPt1 ,MidPt1) , D2R(-5))) + Cutter.append(BezierGen(SouthPt2, MidPt2, GetSideAngleForBezier(SouthPt2,MidPt2) , D2R(-10))) + Cutter.append(BezierGen( EastPt, MidPt2, GetSideAngleForBezier(EastPt ,MidPt2) , D2R(5))) + Cutter.append(BezierGen( WestPt2, MidPt2, GetSideAngleForBezier(WestPt2 ,MidPt2) , D2R(-10))) + Cutter.append(BezierGen( MidPt2, NorthPt, GetSideAngleForBezier(NorthPt ,MidPt2) , D2R(-5))) + + Cutter.append(geompy.MakeEdge(MidPt1, MidPt2)) + + RectFace = geompy.MakePartition([OrigRectFace],Cutter, [], [],4, 0, [], 0) #Creating the partition object + #i=1 + #for SingleCut in Cutter : + # geompy.addToStudy(SingleCut,'Cutter'+str(i)) + # i = i+1 + #geompy.addToStudy(RectFace,'RectFace') + return RectFace + +def Quadrangler (Points): + """ + This function returns a quadranglar face based on four points, non of which 3 are non-colinear. + The points are defined by their 2D [(x1,y1),(x2,y2)..] coordinates. + Note that the list of points is already arranged upon the creation in MacObject + """ + Pt = [] + for Point in Points: Pt.append(geompy.MakeVertex(Point[0], Point[1], 0)) + # The first point is added at the end of the list in order to facilitate the line creation + Pt.append(Pt[0]) + #Draw the lines in order to form the 4 side polygon + Ln=[] + for i in range(4) : Ln.append(geompy.MakeLineTwoPnt(Pt[i],Pt[i+1])) + RectFace = geompy.MakeQuad (Ln[0],Ln[1],Ln[2],Ln[3]) + return RectFace + +def ElemQuartCyl(K): + """ + This function returns a quarter cylinder to box relay of 1 side length, partitioned + with a pitch ratio of K, In other words the side of the box is R*(1+(1/K)) + """ + R = 10.*float(K)/(K+1) + Eps = 10.- R + + Config.theStudy.SetReal("R" , R) + Config.theStudy.SetReal("minusR" , -R) + Config.theStudy.SetReal("Eps", Eps) + + CylWire = geompy.MakeSketcher("Sketcher:F 'R' 0:R 0:L 'Eps':TT 10. 10.0:R 90:L 10.0:R 90:L 'Eps':R 90:C 'minusR' 90.0:WW", [0, 0, 0, 0, 0, 1, 1, 0, -0]) + CylFace = geompy.MakeFace(CylWire, 1) + + SouthPt = geompy.MakeVertex (R+Eps/2., 0., 0) + SouthWestPt = geompy.MakeVertex ( 0.,0., 0) #The origin can be used for practical partionning objectifs + WestPt = geompy.MakeVertex (0., R+Eps/2., 0) + + N = int(math.floor((math.pi*R/4.)/(Eps/2.))) + X = 10.*(1.-1./(N+1)) + + + EastPt = geompy.MakeVertex (10.0, X, 0.) + NorthPt = geompy.MakeVertex ( X, 10.0, 0.) + + DivFactor = 8./(F2D(math.log(K))-0.223) + #MidPt = geompy.MakeVertex ((R+Eps)*math.cos(math.pi/4), (R+Eps)*math.sin(math.pi/4), 0.) + MidPt = geompy.MakeVertex (X-Eps/DivFactor, X-Eps/DivFactor, 0.) + + Cutter = [] + Cutter.append(BezierGen(SouthWestPt, MidPt, GetSideAngleForBezier(SouthWestPt,MidPt) , D2R(-5))) + Cutter.append(BezierGen( EastPt, MidPt, GetSideAngleForBezier(EastPt,MidPt) , D2R(5))) + Cutter.append(BezierGen( MidPt, NorthPt, (-1)**((K<1.25)*1)*D2R(-5), GetSideAngleForBezier(NorthPt,MidPt))) + SMBezier = BezierGen( SouthPt, MidPt, GetSideAngleForBezier(SouthPt ,MidPt) , D2R((K<1.25)*180-5)) + WMBezier = BezierGen( WestPt, MidPt, GetSideAngleForBezier(WestPt, MidPt) , D2R(-5)) + Cutter.append(WMBezier) + Cutter.append(SMBezier) + + for i in range(1,N) : + # Determining intermediate points on the bezier lines and then performing additional cuts + + TempAnglePlus = (math.pi/4)*(1+float(i)/N) + SectionResult = CutnGroup.Go(WMBezier, [(0,0,0,math.sin(TempAnglePlus),-math.cos(TempAnglePlus),0)], [1], ['Dummy'], 0) + TempPt1 = SectionResult[1][0] + TempPt11 = geompy.MakeVertex ((N-i)*X/N, 10., 0) + + TempAngleMinus = (math.pi/4)*(1-float(i)/N) + SectionResult = CutnGroup.Go(SMBezier, [(0,0,0,math.sin(TempAngleMinus),-math.cos(TempAngleMinus),0)], [1], ['Dummy'], 0) + TempPt2 = SectionResult[1][0] + TempPt21 = geompy.MakeVertex (10., (N-i)*X/N, 0) + + Cutter.append(geompy.MakeEdge(SouthWestPt, TempPt1)) + Cutter.append(geompy.MakeEdge(SouthWestPt, TempPt2)) + Cutter.append(geompy.MakeEdge(TempPt1, TempPt11)) + Cutter.append(geompy.MakeEdge(TempPt2, TempPt21)) + + CylFace = geompy.MakePartition([CylFace],Cutter, [], [],4, 0, [], 0) #Creating the partition object + CylFace = geompy.MakeTranslation(CylFace, -5., -5., 0.0) + + return CylFace + +def CompatibilityTest(MacObject): + Type = MacObject.Type + if Type == 'Box11' : + BaseDirPar = [1,1,1,1] + return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar)) + elif Type == 'Box42' : + BaseDirPar = {'SN' : lambda : [3, 3, 4, 2], + 'NS' : lambda : [3, 3, 2, 4], + 'EW' : lambda : [2, 4, 3, 3], + 'WE' : lambda : [4, 2, 3, 3], }[MacObject.MeshPar[1]]() + return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar)) + elif Type == 'BoxAng32' : + BaseDirPar = {'NE' : lambda : [3, 2, 3, 2], + 'NW' : lambda : [2, 3, 3, 2], + 'SW' : lambda : [2, 3, 2, 3], + 'SE' : lambda : [3, 2, 2, 3], }[MacObject.MeshPar[1]]() + return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar)) + elif Type == 'CompBox' : + #print "dx is: ", MacObject.GeoPar[1][1], ". dy is: ",MacObject.GeoPar[1][0] + ReducedRatio = ReduceRatio(MacObject.GeoPar[1][0], MacObject.GeoPar[1][1]) + #print ReducedRatio + BaseDirPar = [ReducedRatio[1], ReducedRatio[1], ReducedRatio[0], ReducedRatio[0]] + return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar)) + + elif Type == 'QuartCyl' : + N = QuarCylParam(MacObject.MeshPar[2])+1 + BaseDirPar = {'NE' : lambda : [2, N, 2, N], + 'NW' : lambda : [N, 2, 2, N], + 'SW' : lambda : [N, 2, N, 2], + 'SE' : lambda : [2, N, N, 2], }[MacObject.MeshPar[1]]() + return int(VecDivRatio(MacObject.DirectionalMeshParams, BaseDirPar)) + elif Type == 'CompBoxF' : + RealRatio = MacObject.GeoPar[1][1]/MacObject.GeoPar[1][0] + Xd = 0 + Yd = 0 + if MacObject.DirectionalMeshParams[2]+MacObject.DirectionalMeshParams[3] : + A = int(max(MacObject.DirectionalMeshParams[2:4])) + Xd = int(VecDivRatio([A,0,0,0], [1,1,1,1])) + if MacObject.DirectionalMeshParams[0]+MacObject.DirectionalMeshParams[1] : + A = int(max(MacObject.DirectionalMeshParams[0:2])) + Yd = int(VecDivRatio([0,0,A,0], [1,1,1,1])) + + if Xd == 0 and Yd : Xd = int(round(Yd/RealRatio)) + elif Yd == 0 : Yd = int(round(RealRatio*Xd)) + + return [Xd,Yd] + elif Type == 'NonOrtho' : + MeanDX = 0.5*(IntLen(MacObject.DirBoundaries(0))+IntLen(MacObject.DirBoundaries(1))) + MeanDY = 0.5*(IntLen(MacObject.DirBoundaries(2))+IntLen(MacObject.DirBoundaries(3))) + RealRatio = MeanDY/MeanDX + Xd = 0 + Yd = 0 + if MacObject.DirectionalMeshParams[2]+MacObject.DirectionalMeshParams[3] : + A = int(max(MacObject.DirectionalMeshParams[2:4])) + Xd = int(VecDivRatio([A,0,0,0], [1,1,1,1])) + if MacObject.DirectionalMeshParams[0]+MacObject.DirectionalMeshParams[1] : + A = int(max(MacObject.DirectionalMeshParams[0:2])) + Yd = int(VecDivRatio([0,0,A,0], [1,1,1,1])) + + if Xd == 0 and Yd : Xd = int(round(Yd/RealRatio)) + elif Yd == 0 : Yd = int(round(RealRatio*Xd)) + + return [Xd,Yd] + +def IntLen (Interval) : + """ + This function returns the length of a given interval even if the latter is not sorted correctly. + """ + return abs(Interval[1]-Interval[0]) + +def NextTo (RefBox, Direction, Extension): + """ + This functions returns geometrical parameters for easy positioning of neighbouring objects. + The input (RefBox) and output are in the form : [(X0,Y0),(DX,DY)] + """ + X0_0 = RefBox[0][0] + Y0_0 = RefBox[0][1] + DX_0 = RefBox[1][0] + DY_0 = RefBox[1][1] + + DirectionalCoef = {'Above' : lambda : [ 0, 1], + 'Below' : lambda : [ 0,-1], + 'Right' : lambda : [ 1, 0], + 'Left ' : lambda : [-1, 0], }[Direction]() + + X0_1 = X0_0+ DirectionalCoef[0] * (DX_0/2.+Extension/2.) + DX_1 = abs(DirectionalCoef[0]) * (Extension) + abs(DirectionalCoef[1])*DX_0 + Y0_1 = Y0_0+ DirectionalCoef[1] * (DY_0/2.+Extension/2.) + DY_1 = abs(DirectionalCoef[1]) * (Extension) + abs(DirectionalCoef[0])*DY_0 + + return [(X0_1,Y0_1),(DX_1,DY_1)] + +def GeomMinMax (PtA, PtB): + """ + This function returns geometrical parameters in the format [(X0,Y0),(DX,DY)]. The input being + the coordinates of two points (Xa,Ya), (Xb,Yb). + """ + # First test that the vector relying the two points is oblique + AB = [PtB[0]- PtA[0],PtB[1]- PtA[1]] + if 0 in AB : + print ("Error: the two points are not correctly defined. In the orthonormal system XOY, it is impossible to define a rectangle with these two points") + return -1 + else: + X0 = 0.5*(PtA[0]+PtB[0]) + Y0 = 0.5*(PtA[1]+PtB[1]) + DX = abs(AB[0]) + DY = abs(AB[1]) + return [(X0,Y0),(DX,DY)] + +def AddIfDifferent (List, Element): + if not(Element in List): + List = List+(Element,) + return List + +def IndexMultiOcc (Array,Element) : + """ + This functions returns the occurrences indices of Element in Array. + As opposed to Array.index(Element) method, this allows determining + multiple entries rather than just the first one! + """ + Output = [] + try : Array.index(Element) + except ValueError : print "No more occurrences" + else : Output.append(Array.index(Element)) + + if not(Output == []) and len(Array) > 1 : + for index, ArrElem in enumerate(Array[Output[0]+1:]) : + if ArrElem == Element : Output.append(index+Output[0]+1) + + return Output + +def SortList (ValList, CritList): + Output = [] + SortedCritList = copy.copy(CritList) + SortedCritList.sort() + for i in range(0,len(ValList)): + if i > 0 : + if not(SortedCritList[i]==SortedCritList[i-1]): + index = IndexMultiOcc(CritList,SortedCritList[i]) + Output= Output + [ValList[j] for j in index] + else : + index = IndexMultiOcc(CritList,SortedCritList[i]) + Output= Output + [ValList[j] for j in index] + + return Output + +def SortPoints(Points): + """ + This function sorts a list of the coordinates of N points as to start at + an origin that represents Xmin and Xmax and then proceed in a counter + clock-wise sense + """ + NbPts = len(Points) + Xmin = min([Points[i][0] for i in range(NbPts)]) + Ymin = min([Points[i][1] for i in range(NbPts)]) + Xmax = max([Points[i][0] for i in range(NbPts)]) + Ymax = max([Points[i][1] for i in range(NbPts)]) + Crit = [(abs(Point[0]-Xmin)+0.1*(Xmax-Xmin))*(abs(Point[1]-Ymin)+0.1*(Ymax-Ymin)) for Point in Points] + #print "Input Points : ", Points + #print "Sorting Criterion : ", Crit + Order = SortList (range(NbPts), Crit) + #print "Sorted Results : ", Order + Output = [] + Output.append(Points[Order[0]]) + + Point0 = Points[Order[0]] + #print "Reference point :", Point0 + + V = [[Point1[0]-Point0[0],Point1[1]-Point0[1]] for Point1 in Points] + Cosines = [-(vec[0]-1E-10)/(math.sqrt(DotProd(vec,vec)+1e-25)) for vec in V] + #print "Cosines criterion :", Cosines + Order = SortList(range(NbPts),Cosines) + #print "Ordered points:", Order + for PtIndex in Order[:-1]: Output.append(Points[PtIndex]) + + return Output + diff --git a/src/Tools/MacMesh/MacMesh/MacObject.py b/src/Tools/MacMesh/MacMesh/MacObject.py new file mode 100644 index 000000000..445cd5d73 --- /dev/null +++ b/src/Tools/MacMesh/MacMesh/MacObject.py @@ -0,0 +1,276 @@ +class MacObject: + """ + This represents a python class definition which contains + all necessary information about the macro object being created + in Salome + """ + + def __init__( self, ObjectType, GeoParameters, MeshParameters, **args ): + """ + Initializes the macro object to be created, saves parameters inside of it, checks for neighboring objects, + determines meshing parameters if necessary and finally launches the generation process. + """ + import Config,GenFunctions + if Config.debug : print "Initializing object No. " + str(len(Config.ListObj)+1) + + if 'publish' in args : + if args['publish']==0 : Config.publish = 0 + else : Config.publish = 1 + else : Config.publish = 1 + + if 'groups' in args : + self.GroupNames = args['groups'] + for group in args['groups'] : + if not(group in Config.Groups) and group : Config.Groups.append(group) + else : self.GroupNames = [None, None, None, None] + + if ObjectType == 'NonOrtho': + if not(len(GeoParameters)==4): print "Error: trying to construct a non-ortho object but the 4 constitutive vertices are not given!" + else : + Xmin = min([GeoParameters[i][0] for i in range(4)]) + Xmax = max([GeoParameters[i][0] for i in range(4)]) + Ymin = min([GeoParameters[i][1] for i in range(4)]) + Ymax = max([GeoParameters[i][1] for i in range(4)]) + self.GeoPar = [(0.5*(Xmin+Xmax),0.5*(Ymin+Ymax)),(Xmax-Xmin,Ymax-Ymin)] + self.PtCoor = GenFunctions.SortPoints(GeoParameters) + else: + self.GeoPar = GeoParameters + [Xmin,Ymin,Xmax,Ymax] = [ self.GeoPar[0][0]-0.5*self.GeoPar[1][0], self.GeoPar[0][1]-0.5*self.GeoPar[1][1] ] + [ self.GeoPar[0][0]+0.5*self.GeoPar[1][0], self.GeoPar[0][1]+0.5*self.GeoPar[1][1] ] + self.PtCoor = [(Xmin,Ymin),(Xmax,Ymin),(Xmax,Ymax),(Xmin,Ymax)] + + self.Type = ObjectType + self.LowBound = [ self.GeoPar[0][0]-0.5*self.GeoPar[1][0], self.GeoPar[0][1]-0.5*self.GeoPar[1][1] ] + self.UpperBound = [ self.GeoPar[0][0]+0.5*self.GeoPar[1][0], self.GeoPar[0][1]+0.5*self.GeoPar[1][1] ] + self.MeshPar = MeshParameters + self.GeoChildren = [] + self.GeoChildrenNames = [] + self.Mesh = [] + self.MeshGroups = [] + self.CheckInterfaces() + if 'auto' in MeshParameters : self.AutoParam() + if not(self.MeshPar[0]<0): self.Generate() + else : + Config.ListObj.append(self) + print("Aborting object creation\n ") + + def Generate(self) : + """ + This method generates the geometrical object with the corresponding mesh once all verifications (CheckInterfaces and AutoParam) + have been accomplished + """ + import GenFunctions, Alarms, Config + self = {'Box11' : lambda : GenFunctions.Box11(self), + 'Box42' : lambda : GenFunctions.Box42(self), + 'BoxAng32' : lambda : GenFunctions.BoxAng32(self), + 'CompBox' : lambda : GenFunctions.CompBox(self), + 'CompBoxF' : lambda : GenFunctions.CompBoxF(self), + 'NonOrtho' : lambda : GenFunctions.NonOrtho(self), + 'QuartCyl' : lambda : GenFunctions.QuartCyl(self) }[self.Type]() + + if Config.debug : Alarms.Message(self.status) # notification on the result of the generation algorithm + + + def CheckInterfaces(self): + """ + This method searches for neighbours for the object being created and saves them inside the Config.Connections + array. This array contains 4 entries per object corresponding to West, East, South, and North neighbours. + Note that an object may have more than one neighbour for a given direction. + """ + import Alarms, Config + from GenFunctions import AddIfDifferent + from CompositeBox import FindCommonSide + + Config.Connections.append([(-1,),(-1,),(-1,),(-1,)]) + itemID = len(Config.ListObj) + # In all cases except non ortho, PrincipleBoxes is unitary and contains the box in question + # In the non-ortho case it contains all possible combinations of boxes with 3 vertices + PrincipleBoxes = self.PrincipleBoxes() + for i, TestObj in enumerate(Config.ListObj): + SecondaryBoxes = TestObj.PrincipleBoxes() + ConnX = 0 + ConnY = 0 + for Box0 in PrincipleBoxes: + for Box1 in SecondaryBoxes: + # Along X + CenterDis = abs(Box1[0][0]-Box0[0][0]) + Extension = 0.5*(Box1[1][0]+Box0[1][0]) + if CenterDis - Extension < -1e-7 : + ConnX = -1 + elif CenterDis - Extension < 1e-7 : + if not(FindCommonSide(self.DirBoundaries(2),TestObj.DirBoundaries(3))==[0,0]) and Box1[0][0] < Box0[0][0] : ConnX = 1 + elif not(FindCommonSide(self.DirBoundaries(3),TestObj.DirBoundaries(2))==[0,0]) and Box1[0][0] >= Box0[0][0]: ConnX = 2 + else : ConnX = 0 + + # Along Y + CenterDis = abs(Box1[0][1]-Box0[0][1]) + Extension = 0.5*(Box1[1][1]+Box0[1][1]) + if CenterDis - Extension < -1e-7 : + ConnY = -1 + elif CenterDis - Extension < 1e-7 : + if not(FindCommonSide(self.DirBoundaries(0),TestObj.DirBoundaries(1))==[0,0]) and Box1[0][1] < Box0[0][1] : ConnY = 1 + elif not(FindCommonSide(self.DirBoundaries(1),TestObj.DirBoundaries(0))==[0,0]) and Box1[0][1] >= Box0[0][1]: ConnY = 2 + else : ConnY = 0 + + if not (ConnX*ConnY == 0) : + if max(ConnX,ConnY) == -1 and not('NonOrtho' in [self.Type,TestObj.Type]) : Alarms.Message(3) + else: + if ConnX == 1 and ConnY == -1: + if Config.Connections[i][1] == (-1,) : Config.Connections[i][1] = (itemID,) + else : Config.Connections[i][1] = AddIfDifferent(Config.Connections[i][1],itemID) + if Config.Connections[itemID][0] == (-1,) : Config.Connections[itemID][0] = (i,) + else : Config.Connections[itemID][0] = AddIfDifferent(Config.Connections[itemID][0],i) + elif ConnX == 2 and ConnY == -1: + if Config.Connections[i][0] == (-1,) : Config.Connections[i][0] = (itemID,) + else : Config.Connections[i][0] = AddIfDifferent(Config.Connections[i][0],itemID) + if Config.Connections[itemID][1] == (-1,) : Config.Connections[itemID][1] = (i,) + else : Config.Connections[itemID][1] = AddIfDifferent(Config.Connections[itemID][1],i) + elif ConnY == 1 and ConnX == -1: + if Config.Connections[i][3] == (-1,) : Config.Connections[i][3] = (itemID,) + else : Config.Connections[i][3] = AddIfDifferent(Config.Connections[i][3],itemID) + if Config.Connections[itemID][2] == (-1,) : Config.Connections[itemID][2] = (i,) + else : Config.Connections[itemID][2] = AddIfDifferent(Config.Connections[itemID][2],i) + elif ConnY ==2 and ConnX == -1: + if Config.Connections[i][2] == (-1,) : Config.Connections[i][2] = (itemID,) + else : Config.Connections[i][2] = AddIfDifferent(Config.Connections[i][2],itemID) + if Config.Connections[itemID][3] == (-1,) : Config.Connections[itemID][3] = (i,) + else : Config.Connections[itemID][3] = AddIfDifferent(Config.Connections[itemID][3],i) + + def AutoParam (self): + """ + This method is called only if the 'auto' keyword is used inside the meshing algorithm. It is based on the + connection results per object and tries to find the correct parameters for obtaining a final compatible mesh + between the objects already present and the one being created. If this is not possible, the method gives an error + message. + """ + import Alarms, Config, GenFunctions, CompositeBox + MeshPar = [0,0,0,0] # initialize the mesh parameter value to be used to -1 + [(X0,Y0),(DX,DY)] = self.GeoPar + ObjectsInvolved = [] + for i, Conn in enumerate(Config.Connections[-1]): + if not ( Conn == (-1,) ): # Meaning that there is one or more neighbors on this direction + for ObjID in Conn : + ToLook0 = [2,3,0,1][i] + ToLook1 = [3,2,1,0][i] + CommonSide = CompositeBox.FindCommonSide(Config.ListObj[ObjID].DirBoundaries(ToLook1),self.DirBoundaries(ToLook0)) + #print "Common Side is:", CommonSide + ToLook2 = [1,0,3,2][i] + #print "Full Side is:", CompositeBox.IntLen(Config.ListObj[ObjID].DirBoundaries(ToLook1)) + #print "Full Segments on this direction are:", Config.ListObj[ObjID].DirectionalMeshParams[ToLook2] + RealSegments = round(Config.ListObj[ObjID].DirectionalMeshParams[ToLook2]*CompositeBox.IntLen(CommonSide)/CompositeBox.IntLen(Config.ListObj[ObjID].DirBoundaries(ToLook1))) + #print "RealSegments :", RealSegments + + MeshPar[i] = MeshPar[i] + RealSegments + ObjectsInvolved.append(ObjID+1) + self.DirectionalMeshParams = MeshPar + self.MeshPar[0] = GenFunctions.CompatibilityTest(self) + + if self.MeshPar[0] < 0 : + Alarms.Message(4) + if self.MeshPar[0] == -1 : print ("Problem encountered with object(s) no. "+str(ObjectsInvolved)) + elif self.MeshPar[0] == -2 : print ("This object has no neighbours !!!") + + def Boundaries (self): + """ + This method returns the global boundaries of the MacObject. [Xmin,Xmax,Ymin,Ymax] + """ + Xmin = min([self.DirBoundaries(i)[0] for i in [0,1]]) + Xmax = max([self.DirBoundaries(i)[1] for i in [0,1]]) + Ymin = min([self.DirBoundaries(i)[0] for i in [2,3]]) + Ymax = max([self.DirBoundaries(i)[1] for i in [2,3]]) + + return [Xmin,Xmax,Ymin,Ymax] + + def DirBoundaries (self, Direction): + """ + This method returns a single interval giving [Xmin,Xmax] or [Ymin,Ymax] according to the required direction. + This works particularly well for nonorthogonal objects. + Direction : [0,1,2,3] <=> [South, North, West, East] + """ + PtCoor = self.PtCoor + PtCoor.append(self.PtCoor[0]) + if type(Direction) is str : + Dir = { 'South' : lambda : 0, + 'North' : lambda : 1, + 'West' : lambda : 2, + 'East' : lambda : 3,}[Direction]() + else : Dir = int(Direction) + + PtIndex = [0,2,3,1][Dir] + DirIndex = [0,0,1,1][Dir] + + return sorted([PtCoor[PtIndex][DirIndex],PtCoor[PtIndex+1][DirIndex]]) + def DirVectors (self, Direction): + """ + This method returns for a given object, the real vectors which define a given direction + The interest in using this method is for non-orthogonal objects where the sides can be + deviated from the orthogonal basis vectors + """ + if type(Direction) is str : + Dir = { 'South' : lambda : 0, + 'North' : lambda : 1, + 'West' : lambda : 2, + 'East' : lambda : 3,}[Direction]() + else : Dir = int(Direction) + PtCoor = self.PtCoor + PtCoor.append(self.PtCoor[0]) + PtIndex = [0,2,3,1][Dir] + return [PtCoor[PtIndex+1][0]-PtCoor[PtIndex][0],PtCoor[PtIndex+1][1]-PtCoor[PtIndex][1],0.] + + def GetBorder (self, Criterion): + import geompy, GenFunctions + + if type(Criterion) is str : + Crit = {'South' : lambda : 0, + 'North' : lambda : 1, + 'West' : lambda : 2, + 'East' : lambda : 3,}[Criterion]() + else : Crit = int(Criterion) + + AcceptedObj = [] + if Crit < 4 : + Boundaries = self.Boundaries() + Research = {0 : lambda : [self.DirVectors(0),1,Boundaries[2]], + 1 : lambda : [self.DirVectors(1),1,Boundaries[3]], + 2 : lambda : [self.DirVectors(2),0,Boundaries[0]], + 3 : lambda : [self.DirVectors(3),0,Boundaries[1]], }[Crit]() + + for i,ElemObj in enumerate(self.GeoChildren): + EdgeIDs = geompy.ExtractShapes(ElemObj,6)# List of Edge IDs belonging to ElemObj + for Edge in EdgeIDs: + if GenFunctions.IsParallel(Edge,Research[0]): + if abs( geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,0))[Research[1]] - Research[2] )< 1e-6 or abs( geompy.PointCoordinates(geompy.GetVertexByIndex(Edge,1))[Research[1]] - Research[2] )< 1e-6 : + AcceptedObj.append(Edge) + else : + CenterSrchPar = {'NE' : lambda : [-1., -1.], + 'NW' : lambda : [ 1., -1.], + 'SW' : lambda : [ 1., 1.], + 'SE' : lambda : [-1., 1.], }[self.MeshPar[1]]() + Radius = self.GeoPar[1][1]*float(self.MeshPar[2])/(self.MeshPar[2]+1) + Center = (self.GeoPar[0][0]+CenterSrchPar[0]*self.GeoPar[1][0]/2.,self.GeoPar[0][1]+CenterSrchPar[1]*self.GeoPar[1][1]/2.,0.) + for i,ElemObj in enumerate(self.GeoChildren): + EdgeIDs = geompy.ExtractShapes(ElemObj,6)# List of Edge IDs belonging to ElemObj + for Edge in EdgeIDs: + if GenFunctions.IsOnCircle(Edge,Center,Radius): + AcceptedObj.append(Edge) + return AcceptedObj + + def PrincipleBoxes (self): + """ + This function returns all possible combination rectangular shape objects that can contain at least 3 of the principle vertices + constituting the MacObject. This is indispensible for the Non-ortho types and shall return a number of 24 possible combinations + """ + from itertools import combinations + Boxes = [] + if self.Type == 'NonOrtho': + for combi in combinations(range(4),3): + Xmin = min([self.PtCoor[i][0] for i in combi]) + Xmax = max([self.PtCoor[i][0] for i in combi]) + Ymin = min([self.PtCoor[i][1] for i in combi]) + Ymax = max([self.PtCoor[i][1] for i in combi]) + Boxes.append([(0.5*(Xmin+Xmax),0.5*(Ymin+Ymax)),(Xmax-Xmin,Ymax-Ymin)]) + else : + Boxes = [self.GeoPar] + + return Boxes + + diff --git a/src/Tools/MacMesh/MacMesh/PublishGroups.py b/src/Tools/MacMesh/MacMesh/PublishGroups.py new file mode 100644 index 000000000..dc8dcd899 --- /dev/null +++ b/src/Tools/MacMesh/MacMesh/PublishGroups.py @@ -0,0 +1,215 @@ +# +import geompy, smesh, SMESH +import math +import Config +########################################################################################################## + +def PublishGroups (): + aFilterManager = smesh.CreateFilterManager() + + # Building geometric and mesh compounds and groups ############################################## + if Config.debug : print "Searching for geometric groups and publishing final compound" + + TempGEOList = [] + TempMESHList = [] + + for MacroObj in Config.ListObj : + TempGEOList += MacroObj.GeoChildren + TempMESHList += MacroObj.Mesh + + FinalCompound = geompy.MakeCompound(TempGEOList) + geompy.addToStudy (FinalCompound,Config.StudyName) + MeshCompound = smesh.Concatenate(TempMESHList, 1, 1, 1e-5) + MeshCompound.SetName(Config.StudyName) + + GroupGEO = [] + for group in Config.Groups : + + # Geometric groups definition + TempGEOList = [] + TempNames = [] + for MacroObj in Config.ListObj : + if group in MacroObj.GroupNames : + Occurences = IndexMultiOcc(MacroObj.GroupNames, group) + for Occ in Occurences : + TempGEOList += MacroObj.GetBorder(Occ) + GroupGEO.append(geompy.MakeCompound(TempGEOList)) + geompy.addToStudyInFather(FinalCompound,GroupGEO[-1],'GR_'+group) + + # Mesh groups definition + Criterion = smesh.Filter.Criterion(18,39,0,'GR_'+group,'GR_'+group,39,39,1e-06,smesh.EDGE,7) + MeshCompound.MakeGroupByCriterion(group,Criterion) + + StudyBuilder = Config.theStudy.NewBuilder() + for MeshObj in TempMESHList: + SO = Config.theStudy.FindObjectIOR(Config.theStudy.ConvertObjectToIOR(MeshObj)) + if SO is not None: StudyBuilder.RemoveObjectWithChildren(SO) + + return MeshCompound + + +def IndexMultiOcc (Array,Element) : + """ + This function returns the occurrences indices of Element in Array. + As opposed to Array.index(Element) method, this allows determining + multiple entries rather than just the first one! + """ + Output = [] + try : Array.index(Element) + except ValueError : print "No more occurrences" + else : Output.append(Array.index(Element)) + + if not(Output == [-1]) and len(Array) > 1 : + for index, ArrElem in enumerate(Array[Output[0]+1:]) : + if ArrElem is Element : Output.append(index+Output[0]+1) + + return Output + +def Publish (ObjToPublish): + for i,GeoObj in enumerate(ObjToPublish) : geompy.addToStudy(GeoObj,"Sub_"+str(i)) + +def RevolveMesh(MainMesh,**args): + """ + This function premits to revolute and scale a 2D mesh while transforming the edge + groups into face groups. Moreover, the function automatically creates the face groups + corresponding to the symmetry lower and upper faces + Facultatif arguments are : + - Center [X,Y,Z], origin being the default + - Direction [VX,VY,VZ], x-axis being the default + - AngleDeg or AngleRad : ALPHA, 10 degrees being the default + - Scale : BETA, no scaling being default + """ + ################################################################################ + # Reading input arguments and proceeding to defaults if necessary + ################################################################################ + if 'Center' in args : CenterCoor = [float(Coor) for Coor in args['Center']] + else : + print "\nThe coordinates of the center of revolution were not given\nThe origin is used by default." + CenterCoor = [0.,0.,0.] + + if 'Direction' in args : Direction = [float(Dir) for Dir in args['Direction']] + else : + print "\nThe axis vector of revolution was not given\nThe x-axis is used by default." + Direction = [1.,0.,0.] + + if 'AngleDeg' in args : Angle = float(args['AngleDeg'])*math.pi/180. + elif 'AngleRad' in args : Angle = float(args['AngleRad']) + else : + print "\nThe revolution angle was not given\nAn angle of 10 degrees is used by default." + Angle = 10.*math.pi/180. + + if 'Scale' in args : Scale = float(args['Scale']) + else : Scale = 1. + + + # Creating the lower face group LOFAC + LOFAC = MainMesh.CreateEmptyGroup( SMESH.FACE, 'LOFAC' ) + LOFAC.AddFrom(MainMesh.GetMesh()) + + GR_Names = MainMesh.GetGroupNames() + GRs = MainMesh.GetGroups() + Rev3DMeshGroups = MainMesh.RotationSweepObject2D( MainMesh, SMESH.AxisStruct( CenterCoor[0], CenterCoor[1], CenterCoor[2], Direction[0], Direction[1], Direction[2] ), Angle, 1, 1e-05 ,True) + + # Adding an EDGE suffix to the edge groups (to be deleted eventually by the user...) + for GR in GRs: + CurrentName = GR.GetName() + if CurrentName in GR_Names and not(CurrentName=='LOFAC'): # Meaning that this is an old edge group + GR.SetName(CurrentName+'_EDGE') + + # Removing the _rotated prefix from the rotated FACE groups + for GR in Rev3DMeshGroups: + CurrentName = GR.GetName() + if CurrentName=='LOFAC_rotated' : + GR.SetName('VOL') + else : + #Index = [ GR_Names[i] in CurrentName for i in range(0,len(GR_Names)) ].index(True) + #GR.SetName(GR_Names[Index]) + GR.SetName(CurrentName[:-8]) + + # Creating the upper face group HIFAC + ALLFAC = MainMesh.CreateEmptyGroup( SMESH.FACE, 'ALLFAC' ) + ALLFAC.AddFrom(MainMesh.GetMesh()) + + HIFAC = MainMesh.GetMesh().CutListOfGroups( [ ALLFAC ], [LOFAC] + [ MeshGroup for MeshGroup in Rev3DMeshGroups if not(MeshGroup.GetName()=='VOL') ], 'HIFAC' ) + + # Scaling down the mesh to meter units + if not(Scale==1.): + MeshEditor = MainMesh.GetMeshEditor() + MeshEditor.Scale( MainMesh.GetMesh(), SMESH.PointStruct( 0, 0, 0 ) ,[ Scale, Scale, Scale ], 0 ) + + +def ExtrudeMesh(MainMesh,**args): + """ + This function premits to extrude and scale a 2D mesh while transforming the edge + groups into face groups. Moreover, the function automatically creates the face groups + corresponding to the symmetry lower and upper faces + Facultatif arguments are : + - Direction [VX,VY,VZ], z-axis being default + - Distance : D, default is 1 + - NSteps : the object will be extruded by NSteps*Distance, default is Nsteps = 1 + - Scale : BETA, no scaling being default + """ + ################################################################################ + # Reading input arguments and proceeding to defaults if necessary + ################################################################################ + if 'Distance' in args : Distance = float(args['Distance']) + else : + print "\nThe extrusion distance was not given\nA default value of 1 is used." + Distance = 1. + + if 'Direction' in args : Direction = NormalizeVector([float(Dir) for Dir in args['Direction']],Distance) + else : + print "\nThe extrusion vector of revolution was not given\nThe z-axis is used by default." + Direction = NormalizeVector([0.,0.,1.],Distance) + + if 'Scale' in args : Scale = float(args['Scale']) + else : Scale = 1. + + if 'NSteps' in args : NSteps = int(args['NSteps']) + else : NSteps = 1 + + # Creating the lower face group LOFAC + LOFAC = MainMesh.CreateEmptyGroup( SMESH.FACE, 'LOFAC' ) + LOFAC.AddFrom(MainMesh.GetMesh()) + + GR_Names = MainMesh.GetGroupNames() + GRs = MainMesh.GetGroups() + Ext3DMeshGroups = MainMesh.ExtrusionSweepObject2D(MainMesh,SMESH.DirStruct(SMESH.PointStruct(Direction[0],Direction[1],Direction[2])), NSteps, True) + + # Adding an EDGE suffix to the edge groups (to be deleted eventually by the user...) + for GR in GRs: + CurrentName = GR.GetName() + if CurrentName in GR_Names and not(CurrentName=='LOFAC'): # Meaning that this is an old edge group + GR.SetName(CurrentName+'_EDGE') + + # Removing the _rotated prefix from the rotated FACE groups + for GR in Ext3DMeshGroups: + CurrentName = GR.GetName() + if CurrentName=='LOFAC_extruded' : + GR.SetName('VOL') + else : + #Index = [ GR_Names[i] in CurrentName for i in range(0,len(GR_Names)) ].index(True) + #GR.SetName(GR_Names[Index]) + GR.SetName(CurrentName[:-9]) + + # Creating the upper face group HIFAC + ALLFAC = MainMesh.CreateEmptyGroup( SMESH.FACE, 'ALLFAC' ) + ALLFAC.AddFrom(MainMesh.GetMesh()) + + HIFAC = MainMesh.GetMesh().CutListOfGroups( [ ALLFAC ], [LOFAC] + [ MeshGroup for MeshGroup in Ext3DMeshGroups if not(MeshGroup.GetName()=='VOL') ], 'HIFAC' ) + + # Scaling down the mesh to meter units + if not(Scale==1.): + MeshEditor = MainMesh.GetMeshEditor() + MeshEditor.Scale( MainMesh.GetMesh(), SMESH.PointStruct( 0, 0, 0 ) ,[ Scale, Scale, Scale ], 0 ) + + +def NormalizeVector (V,Norm): + """ + This function returns a normalized vector (magnitude = Norm), parallel to the entered one + """ + V = [float(Coor) for Coor in V] + Norm = float(Norm) + MagV = math.sqrt(V[0]*V[0]+V[1]*V[1]+V[2]*V[2]) + return [Coor*Norm/MagV for Coor in V] + diff --git a/src/Tools/MacMesh/MacMesh/SharpAngle.py b/src/Tools/MacMesh/MacMesh/SharpAngle.py new file mode 100644 index 000000000..3ddc89582 --- /dev/null +++ b/src/Tools/MacMesh/MacMesh/SharpAngle.py @@ -0,0 +1,226 @@ +# This is an automation of the sharp angle object, with a corner at (X0,Y0), side length : Extension and a fine local meshing : LocalMeshing +# The corner orientation is defined as NE (North-East) , NW (North-West), SE, or SW. The object's "arm" is 8/14 of Extension +# | | 8 6 +# ------- --------- +# ----> | | <---- +# | NW NE | oo +# _____| |_____ + +import sys, salome, geompy, smesh, SMESH, math, commands +CWD = commands.getoutput('pwd') +sys.path.append(CWD) + +from MacObject import * +from CompositeBox import * +import Config, GenFunctions + +def SharpAngleOut (X0 , Y0 , DX , DY , DLocal, LocalMeshing , CornerOrientation , NLevels, **args) : + if DLocal == 'auto' : DLocal = float(min(DX,DY)) + + BoxSide = DLocal/(2.**(NLevels+1)) + InternalMeshing = int(math.ceil(BoxSide/(3*LocalMeshing))) + InternalMeshing = InternalMeshing+InternalMeshing%2 # An even number is needed, otherwise the objects would not be compatible once created + if InternalMeshing == 0 : InternalMeshing = 2 # This sets a minimum meshing condition in order to avoid an error. The user is notified of the value considered for the local meshing + print "Possible Local meshing is :", BoxSide/(3*InternalMeshing), "\nThis value is returned by this function for your convenience" + + DirPar = {'NE' : lambda : ['NE', 'NW', 'SE', 'EW', 'NW', 'SN', 'SN', 'NE', 'WE', 'WE', 'SE', 'NS'], + 'NW' : lambda : ['NW', 'NE', 'SW', 'WE', 'NE', 'SN', 'SN', 'NW', 'EW', 'EW', 'SW', 'NS'], + 'SE' : lambda : ['SE', 'SW', 'NE', 'EW', 'SW', 'NS', 'NS', 'SE', 'WE', 'WE', 'NE', 'SN'], + 'SW' : lambda : ['SW', 'SE', 'NW', 'WE', 'SE', 'NS', 'NS', 'SW', 'EW', 'EW', 'NW', 'SN'], }[CornerOrientation]() + + CoefVer = {'NE' : lambda : 1, + 'NW' : lambda : 1, + 'SE' : lambda : -1, + 'SW' : lambda : -1, }[CornerOrientation]() + + CoefHor = {'NE' : lambda : 1, + 'NW' : lambda : -1, + 'SE' : lambda : 1, + 'SW' : lambda : -1, }[CornerOrientation]() + + ToLook = {'NE' : lambda : [0,2,1,3], + 'NW' : lambda : [0,3,1,2], + 'SE' : lambda : [1,2,0,3], + 'SW' : lambda : [1,3,0,2], }[CornerOrientation]() + + if args.__contains__('groups') : + GroupNames = args['groups'] + else : GroupNames = [None, None, None, None, None, None] + + GN00 = GroupArray(ToLook[0],GroupNames[0]) + GN01 = GroupArray(ToLook[1],GroupNames[1]) + + GN1 = GroupArray([ToLook[0],ToLook[1]],[GroupNames[0],GroupNames[5]]) + GN7 = GroupArray([ToLook[0],ToLook[1]],[GroupNames[4],GroupNames[1]]) + + if DY == DLocal : + GN2 = GroupArray([ToLook[1],ToLook[2]],[GroupNames[5],GroupNames[2]]) + GN3 = GroupArray(ToLook[2],GroupNames[2]) + if DX == DLocal: + GN4 = GroupArray([ToLook[2],ToLook[3]],[GroupNames[2],GroupNames[3]]) + GN5 = GroupArray(ToLook[3],GroupNames[3]) + GN6 = GroupArray([ToLook[3],ToLook[0]],[GroupNames[3],GroupNames[4]]) + else : + GN4 = GroupArray(ToLook[2],GroupNames[2]) + GN5 = [None,None,None,None] + GN6 = GroupArray(ToLook[0],GroupNames[4]) + GN21 = GroupArray([ToLook[3],ToLook[0],ToLook[2]],[GroupNames[3],GroupNames[4],GroupNames[2]]) + else : + GN2 = GroupArray(ToLook[1],GroupNames[5]) + GN3 = [None,None,None,None] + if DX == DLocal: + GN4 = GroupArray(ToLook[3],GroupNames[3]) + GN5 = GroupArray(ToLook[3],GroupNames[3]) + GN6 = GroupArray([ToLook[3],ToLook[0]],[GroupNames[3],GroupNames[4]]) + GN22 = GroupArray([ToLook[1],ToLook[2],ToLook[3]],[GroupNames[5],GroupNames[2],GroupNames[3]]) + else : + GN4 = [None,None,None,None] + GN5 = [None,None,None,None] + GN6 = GroupArray(ToLook[0],GroupNames[4]) + GN21 = GroupArray([ToLook[3],ToLook[0]],[GroupNames[3],GroupNames[4]]) + GN22 = GroupArray([ToLook[1],ToLook[2]],[GroupNames[5],GroupNames[2]]) + GN23 = GroupArray([ToLook[2],ToLook[3]],[GroupNames[2],GroupNames[3]]) + + Obj = [] + + Obj.append(MacObject('BoxAng32',[(X0+CoefHor*BoxSide/2,Y0+CoefVer*BoxSide/2),(BoxSide,BoxSide)],[InternalMeshing,DirPar[0]])) + Obj.append(MacObject('BoxAng32',[(X0-CoefHor*BoxSide/2,Y0+CoefVer*BoxSide/2),(BoxSide,BoxSide)],['auto',DirPar[1]], groups = GroupArray(ToLook[0],GroupNames[0]))) + Obj.append(MacObject('BoxAng32',[(X0+CoefHor*BoxSide/2,Y0-CoefVer*BoxSide/2),(BoxSide,BoxSide)],['auto',DirPar[2]], groups = GroupArray(ToLook[1],GroupNames[1]))) + + for N in range (1,NLevels+1): + n = N-1 + if N < NLevels : + Obj.append(MacObject('Box42',[(X0-CoefHor*BoxSide*(2**n)*3/2,Y0+CoefVer*(2**n)*BoxSide/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[3]] , groups = GN00)) + Obj.append(MacObject('BoxAng32',[(X0-CoefHor*(2**n)*BoxSide*3/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[4]] )) + Obj.append(MacObject('Box42',[(X0-CoefHor*(2**n)*BoxSide/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[5]] )) + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[6]] )) + Obj.append(MacObject('BoxAng32',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[7]] )) + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0+CoefVer*(2**n)*BoxSide/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[8]] )) + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0-CoefVer*(2**n)*BoxSide/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[9]] )) + Obj.append(MacObject('BoxAng32',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0-CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[10]] )) + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide/2,Y0-CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[11]] , groups = GN01)) + else : + Obj.append(MacObject('Box42',[(X0-CoefHor*BoxSide*(2**n)*3/2,Y0+CoefVer*(2**n)*BoxSide/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[3]] , groups = GN1)) + Obj.append(MacObject('BoxAng32',[(X0-CoefHor*(2**n)*BoxSide*3/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[4]] , groups = GN2)) + Obj.append(MacObject('Box42',[(X0-CoefHor*(2**n)*BoxSide/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[5]] , groups = GN3)) + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[6]] , groups = GN3)) + Obj.append(MacObject('BoxAng32',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[7]] , groups = GN4)) + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0+CoefVer*(2**n)*BoxSide/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[8]] , groups = GN5)) + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0-CoefVer*(2**n)*BoxSide/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[9]] , groups = GN5)) + Obj.append(MacObject('BoxAng32',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0-CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[10]], groups = GN6)) + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide/2,Y0-CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[11]] , groups = GN7)) + + OuterMeshing = (3/2)*InternalMeshing*2**(NLevels-1) + OuterSegLength = (DLocal/OuterMeshing) + + if DX > DLocal : + dX = DX - DLocal + Obj.append(MacObject('CompBoxF',[(X0+CoefHor*(DX)/2.,Y0),(dX,DLocal)],['auto'], groups = GN21)) + if DY > DLocal : + dY = DY - DLocal + if DX > DLocal : + Obj.append(MacObject('CompBoxF',[(X0+CoefHor*DX/2.,Y0+CoefVer*(DY)/2.),(DX-DLocal,dY)],['auto'], groups = GN23)) + + Obj.append(MacObject('CompBoxF',[(X0,Y0+CoefVer*(DY)/2.),(DLocal,dY)],['auto'], groups = GN22)) + + return Obj + +def SharpAngleIn (X0 , Y0 , DX , DY , DLocal, LocalMeshing , CornerOrientation , NLevels, **args) : + if DLocal == 'auto' : DLocal = float(min(DX,DY)) + + BoxSide = DLocal/(2.**(NLevels)) + InternalMeshing = int(math.ceil(BoxSide/(3*LocalMeshing))) + InternalMeshing = InternalMeshing+InternalMeshing%2 # An even number is needed, otherwise the objects would not be compatible once created + if InternalMeshing == 0 : InternalMeshing = 2 # This sets a minimum meshing condition in order to avoid an error. The user is notified of the value considered for the local meshing + print "Possible Local meshing is :", BoxSide/(3*InternalMeshing), "\nThis value is returned by this function for your convenience..." + + DirPar = {'NE' : lambda : ['NE', 'SN', 'NE', 'WE'], + 'NW' : lambda : ['NW', 'SN', 'NW', 'EW'], + 'SE' : lambda : ['SE', 'NS', 'SE', 'WE'], + 'SW' : lambda : ['SW', 'NS', 'SW', 'EW'], }[CornerOrientation]() + + CoefVer = {'NE' : lambda : 1, + 'NW' : lambda : 1, + 'SE' : lambda : -1, + 'SW' : lambda : -1, }[CornerOrientation]() + + CoefHor = {'NE' : lambda : 1, + 'NW' : lambda : -1, + 'SE' : lambda : 1, + 'SW' : lambda : -1, }[CornerOrientation]() + + ToLook = {'NE' : lambda : [0,2,1,3], + 'NW' : lambda : [0,3,1,2], + 'SE' : lambda : [1,2,0,3], + 'SW' : lambda : [1,3,0,2], }[CornerOrientation]() + + if args.__contains__('groups') : + GroupNames = args['groups'] + else : GroupNames = [None, None, None, None] + + GN01 = GroupArray([ToLook[0],ToLook[1]],[GroupNames[ToLook[0]],GroupNames[ToLook[1]]]) + GN02 = GroupArray(ToLook[1],GroupNames[ToLook[1]]) + GN03 = [None, None, None, None] + GN04 = GroupArray(ToLook[0],GroupNames[ToLook[0]]) + + if DY == DLocal : + GN05 = GroupArray([ToLook[1],ToLook[2]],[GroupNames[ToLook[1]],GroupNames[ToLook[2]]]) + GN08 = GroupArray([ToLook[0],ToLook[2],ToLook[3]],[GroupNames[ToLook[0]],GroupNames[ToLook[2]],GroupNames[ToLook[3]]]) + if DX == DLocal: + GN06 = GroupArray([ToLook[2],ToLook[3]],[GroupNames[ToLook[2]],GroupNames[ToLook[3]]]) + GN07 = GroupArray([ToLook[0],ToLook[3]],[GroupNames[ToLook[0]],GroupNames[ToLook[3]]]) + else : + GN06 = GroupArray(ToLook[2],GroupNames[ToLook[2]]) + GN07 = GroupArray(ToLook[0],GroupNames[ToLook[0]]) + else : + GN05 = GroupArray(ToLook[1],GroupNames[ToLook[1]]) + if DX == DLocal : + GN06 = GroupArray(ToLook[3],GroupNames[ToLook[3]]) + GN07 = GroupArray([ToLook[0],ToLook[3]],[GroupNames[ToLook[0]],GroupNames[ToLook[3]]]) + GN10 = GroupArray([ToLook[1],ToLook[2],ToLook[3]],[GroupNames[ToLook[1]],GroupNames[ToLook[2]],GroupNames[ToLook[3]]]) + else : + GN06 = [None, None, None, None] + GN07 = GroupArray(ToLook[0],GroupNames[ToLook[0]]) + GN08 = GroupArray([ToLook[0],ToLook[3]],[GroupNames[ToLook[0]],GroupNames[ToLook[3]]]) + GN09 = GroupArray([ToLook[2],ToLook[3]],[GroupNames[ToLook[2]],GroupNames[ToLook[3]]]) + GN10 = GroupArray([ToLook[1],ToLook[2]],[GroupNames[ToLook[1]],GroupNames[ToLook[2]]]) + + Obj = [] + + Obj.append(MacObject('BoxAng32',[(X0+CoefHor*BoxSide/2,Y0+CoefVer*BoxSide/2),(BoxSide,BoxSide)],[InternalMeshing,DirPar[0]],groups = GN01)) + + for N in range (1,NLevels+1): + n = N-1 + if N < NLevels : + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[1]],groups = GN02)) + Obj.append(MacObject('BoxAng32',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[2]],groups = GN03)) + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0+CoefVer*(2**n)*BoxSide/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[3]],groups = GN04)) + else : + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[1]],groups = GN05)) + Obj.append(MacObject('BoxAng32',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0+CoefVer*(2**n)*BoxSide*3/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[2]],groups = GN06)) + Obj.append(MacObject('Box42',[(X0+CoefHor*(2**n)*BoxSide*3/2,Y0+CoefVer*(2**n)*BoxSide/2),((2**n)*BoxSide,(2**n)*BoxSide)],['auto',DirPar[3]],groups = GN07)) + + OuterMeshing = (3/2)*InternalMeshing*2**(NLevels-1) + OuterSegLength = (DLocal/OuterMeshing) + + if DX > DLocal : + dX = DX - DLocal + Obj = Obj + CompositeBox(X0+CoefHor*(DLocal+dX/2.),Y0+CoefVer*(DLocal)/2.,dX,DLocal, groups = GN08) + if DY > DLocal : + dY = DY - DLocal + + if DX > DLocal : + Obj = Obj + CompositeBox(X0+CoefHor*(DLocal+(DX-DLocal)/2.),Y0+CoefVer*(DLocal+dY/2.),DX-DLocal,dY, groups = GN09) + + Obj = Obj + CompositeBox(X0+CoefHor*DLocal/2,Y0+CoefVer*(DLocal+dY/2.),DLocal,dY,groups = GN10) + + return Obj + +def GroupArray(indices, GroupNames) : + if type(indices) is int : + indices = [indices] + GroupNames = [GroupNames] + Output = [None,None,None,None] + for i, ind in enumerate(indices) : + Output[ind] = GroupNames[i] + return Output