src/Tools/ZCracksPlug/utilityFunctions.py

   1 #import sys
   2 #sys.path.append('/home/I60976/00_PROJETS/2015_INTEGRATION_ZCRACKS/zcracks_salome/zcracks')
   3
   4 import numpy, subprocess, sys
   5 from os import remove, getpid, path, environ
   6 from output import message
   7
   8 def calcCoordVectors(normalIN, directionIN):
   9   V3TEMP=numpy.cross(normalIN,directionIN)
  10   directionTEMP=numpy.cross(V3TEMP,normalIN)
  11
  12   normalOUT=numpy.array(normalIN)/numpy.linalg.norm(normalIN)
  13   directionOUT=numpy.array(directionTEMP)/numpy.linalg.norm(directionTEMP)
  14   V3OUT=numpy.array(V3TEMP)/numpy.linalg.norm(V3TEMP)
  15   return(normalOUT, directionOUT, V3OUT)
  16
  17
  18 def testStrictRange(x, inf=0.0, sup=False):
  19   test=False
  20   c1=(type(x)==list)
  21   if c1:
  22     c2=(len(x)==1)
  23     if c2:
  24       c3=(type(x[0])==type(inf))
  25       if c3:
  26         c4=(x[0]>inf)
  27         c5=True
  28         if sup!=False:
  29           c5=(x[0]<sup)
  30         if c4 and c5:
  31           test=True
  32   return(test)
  33
  34 def test3dVector(x):
  35   test=False
  36   c1=(type(x)==list)
  37   if c1:
  38     c2=(len(x)==3)
  39     if c2:
  40       c3=(type(x[0])==float)
  41       c4=(type(x[1])==float)
  42       c5=(type(x[2])==float)
  43       if c3 and c4 and c5:
  44         test=True
  45   return(test)
  46
  47 def testRange(x, inf=0.0, sup=False):
  48   test=False
  49   c1=(type(x)==list)
  50   if c1:
  51     c2=(len(x)==1)
  52     if c2:
  53       c3=(type(x[0])==type(inf))
  54       if c3:
  55         c4=(x[0]>=inf)
  56         c5=True
  57         if sup!=False:
  58           c5=(x[0]<=sup)
  59         if c4 and c5:
  60           test=True
  61   return(test)
  62
  63 def check(data):
  64   OK=True
  65
  66   test=False
  67   c1=(data['crackedName']!='')
  68   if c1:
  69     test=True
  70   if not test:
  71     message('E','Invalid Cracked name',goOn=True)
  72     OK=False
  73
  74   test=False
  75   c1=path.isfile(data['saneName'])
  76   if c1:
  77     c2=(data['saneName']!=data['crackedName'])
  78     if c2:
  79       test=True
  80     else:
  81       message('E','sane mesh and cracked mesh are identical',goOn=True)
  82       OK=False
  83   if not test:
  84     message('E','Bad sane mesh file',goOn=True)
  85     OK=False
  86
  87   test=testStrictRange(data['minSize'])
  88   if not test:
  89     message('E','invalid min size',goOn=True)
  90     OK=False
  91
  92   test=testStrictRange(data['maxSize'])
  93   if not test:
  94     message('E','invalid max size',goOn=True)
  95     OK=False
  96
  97   if OK:
  98     test=(data['maxSize'][0]>=data['minSize'][0])
  99     if not test:
 100       message('E','min size greater than max size',goOn=True)
 101       OK=False
 102
 103   test=testStrictRange(data['extractLength'])
 104   if not test:
 105     message('E','invalid extract length',goOn=True)
 106     OK=False
 107
 108   test=testRange(data['gradation'], inf=1.0)
 109   if not test:
 110     message('E','invalid Gradation',goOn=True)
 111     OK=False
 112
 113   test=testRange(data['layers'], inf=1)
 114   if not test:
 115     message('E','invalid layers',goOn=True)
 116     OK=False
 117
 118   test=testRange(data['iterations'], inf=1)
 119   if not test:
 120     message('E','invalid iterations',goOn=True)
 121     OK=False
 122   return(OK)
 123
 124
 125 def calcElemSize(A, R):
 126   x=R*(1.-numpy.cos(A/2.))
 127   h=R*numpy.sin(A/2.)
 128   return(x, h)
 129
 130 def meshCrack(geomObject, minSize, maxSize, chordal, dim):
 131   import salome
 132
 133   salome.salome_init()
 134   theStudy = salome.myStudy
 135
 136   import  SMESH, SALOMEDS
 137   from salome.smesh import smeshBuilder
 138   smesh = smeshBuilder.New(theStudy)
 139   Maillage = smesh.Mesh(geomObject)
 140
 141   if dim==3:
 142     MG_CADSurf = Maillage.Triangle(algo=smeshBuilder.MG_CADSurf)
 143     MG_CADSurf_Parameters = MG_CADSurf.Parameters()
 144     MG_CADSurf_Parameters.SetPhysicalMesh( 0 )
 145     MG_CADSurf_Parameters.SetGeometricMesh( 1 )
 146     MG_CADSurf_Parameters.SetMinSize( minSize )
 147     MG_CADSurf_Parameters.SetMaxSize( maxSize )
 148     MG_CADSurf_Parameters.SetChordalError( chordal )
 149
 150   elif dim==2:
 151     Regular_1D = Maillage.Segment()
 152     Adaptive = Regular_1D.Adaptive(minSize,maxSize,chordal)
 153     NETGEN_2D_ONLY = Maillage.Triangle(algo=smeshBuilder.NETGEN_2D)
 154   else:
 155     message('E',"error in mesh dimension",goOn=True)
 156
 157   isDone = Maillage.Compute()
 158
 159   #crack1 = Maillage.CreateEmptyGroup( SMESH.NODE, 'crack' )
 160   #nbAdd = crack1.AddFrom( Maillage.GetMesh() )
 161   #crack2 = Maillage.CreateEmptyGroup( SMESH.NODE, 'surface' )
 162   #nbAdd = crack2.AddFrom( Maillage.GetMesh() )
 163
 164   return(Maillage)
 165
 166 def extendElsets(meshFile, outFile=None):
 167
 168   if outFile==None: outFile=meshFile
 169
 170   if not path.isfile(meshFile):
 171     message('E','Mesh med file is not valid')
 172     return('error')
 173
 174   import SMESH, salome
 175   #salome.salome_init()
 176   theStudy = salome.myStudy
 177   from salome.smesh import smeshBuilder
 178   smesh = smeshBuilder.New(theStudy)
 179
 180   ([mesh], status) = smesh.CreateMeshesFromMED(meshFile)
 181
 182   mesh=cleanGroups(mesh)
 183
 184   # Node color status
 185   nodeList=mesh.GetNodesId()
 186   volElemList=mesh.GetElementsByType(SMESH.VOLUME)
 187   surfElemList=mesh.GetElementsByType(SMESH.FACE)
 188   edgeElemList=mesh.GetElementsByType(SMESH.EDGE)
 189   colorList=[-1]*len(nodeList)
 190
 191   case2D=True
 192   for group in mesh.GetGroups():
 193     if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides' : case2D=False
 194
 195   sides=[]
 196   for group in mesh.GetGroups():
 197     if case2D:
 198       if group.GetType()==SMESH.FACE and group.GetName()[:5]=='sides':
 199         sides.append(group)
 200     else:
 201       if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides':
 202         sides.append(group)
 203
 204   sortedSides=[None]*len(sides)
 205   for group in sides:
 206     N=group.GetName().replace('sides','').replace('_bset','').replace(' ','')
 207     N=int(N)
 208     sortedSides[N]=group
 209
 210     elems=group.GetIDs()
 211     for elemId in elems:
 212       for elemNodeId in mesh.GetElemNodes(elemId) :
 213         colorList[elemNodeId-1]=N
 214   #print colorList
 215
 216   crackOnly=True
 217   for iN in range(len(sides)/2):
 218     side0=sortedSides[2*iN]
 219     side1=sortedSides[2*iN+1]
 220     elemsOfside0=side0.GetIDs()
 221     elemsOfside1=side1.GetIDs()
 222     NodesOfside0=[]
 223     NodesOfside1=[]
 224     for elem in elemsOfside0: NodesOfside0+=mesh.GetElemNodes(elem)
 225     for elem in elemsOfside1: NodesOfside1+=mesh.GetElemNodes(elem)
 226     front=set(NodesOfside0).intersection(set(NodesOfside1))
 227     if len(front)==0: crackOnly=False
 228
 229   if crackOnly:
 230     mesh.ExportMED(outFile, 0, SMESH.MED_V2_2, 1, None ,1)
 231     return('crack')
 232
 233   # Propagates color using elem connectivity
 234   # Always propagates max color
 235
 236   #elemToTreat=volElemList
 237
 238   #while len(elemToTreat)>0 :
 239     #print len(elemToTreat)
 240     #for elemId in elemToTreat:
 241       #minColor=sys.maxint
 242       #maxColor=-sys.maxint
 243       #for elemNodeId in mesh.GetElemNodes(elemId) :
 244         #nodeColor=colorList[elemNodeId-1]
 245         #if nodeColor<minColor : minColor=nodeColor
 246         #if nodeColor>maxColor : maxColor=nodeColor
 247       #if minColor!=maxColor :
 248         #elemToTreat.remove(elemId)
 249         #for elemNodeId in mesh.GetElemNodes(elemId) :
 250           #colorList[elemNodeId-1]=maxColor
 251
 252   ifChanged=True
 253   if case2D:
 254     elemList=[surfElemList,edgeElemList]
 255     grElemList=[[],[]]
 256   else:
 257     elemList=[volElemList,surfElemList,edgeElemList]
 258     grElemList=[[],[],[]]
 259
 260   while ifChanged :
 261     ifChanged=False
 262     for elemId in elemList[0]:
 263       minColor=sys.maxint
 264       maxColor=-sys.maxint
 265       for elemNodeId in mesh.GetElemNodes(elemId) :
 266         nodeColor=colorList[elemNodeId-1]
 267         if nodeColor<minColor : minColor=nodeColor
 268         if nodeColor>maxColor : maxColor=nodeColor
 269       if minColor!=maxColor :
 270         ifChanged = True
 271         for elemNodeId in mesh.GetElemNodes(elemId) :
 272           colorList[elemNodeId-1]=maxColor
 273
 274   for l in grElemList:
 275     for x in range(len(sides)):
 276       l.append([])
 277
 278   for N, el in enumerate(elemList):
 279     for elemId in el:
 280       elemNodesId=mesh.GetElemNodes(elemId)
 281       elemColor=colorList[elemNodesId[0]-1]
 282       if elemColor>=0:
 283         grElemList[N][elemColor].append(elemId)
 284
 285   #for elemId in surfElemList:
 286     #elemNodesId=mesh.GetElemNodes(elemId)
 287     #elemColor=colorList[elemNodesId[0]-1]
 288     #if elemColor>=0:
 289       #selem[elemColor].append(elemId)
 290
 291   for n in range(len(sides)):
 292     if case2D:
 293       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.FACE,grElemList[0][n])
 294       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.EDGE,grElemList[1][n])
 295     else:
 296       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.VOLUME,grElemList[0][n])
 297       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.FACE,grElemList[1][n])
 298       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.EDGE,grElemList[2][n])
 299
 300   if outFile==None: outFile=meshFile
 301   mesh.ExportMED(outFile, 0, SMESH.MED_V2_2, 1, None ,1)
 302   return(True)
 303
 304
 305 def cleanGroups(mesh):
 306   import SMESH
 307   for group in mesh.GetGroups():
 308     if '_bset' in group.GetName():
 309       group.SetName(group.GetName().replace('_bset',''))
 310
 311     if group.GetType()==SMESH.NODE:
 312       if group.GetName() in ['SURFACE','lip','SFRONT_NODES','FRONT']: mesh.RemoveGroup(group)
 313
 314     #elif group.GetType()==SMESH.EDGE:
 315
 316     elif group.GetType()==SMESH.FACE:
 317       if group.GetName() in ['SURFACE','Nlip']:
 318         mesh.RemoveGroup(group)
 319
 320     elif group.GetType()==SMESH.VOLUME:
 321       if (group.GetName() in ['ELSET0','AUTO']) or (group.GetName()[:4] in ['SIDE']) :
 322         mesh.RemoveGroup(group)
 323
 324   return(mesh)
 325
 326
 327 def getMaxAspectRatio(tmpdir):
 328   logFile=path.join(tmpdir,'MESHING_OUTPUT')
 329   print logFile
 330   if not path.isfile(logFile): return(-1)
 331
 332   import re
 333   f = open(logFile, "r")
 334   for line in f:
 335     if re.search("WORST ELEMENT QUALITY", line):  maxAR=line
 336
 337   f.close()
 338   for r in [' ','WORSTELEMENTQUALITY','\n']: maxAR=maxAR.replace(r,'')
 339   return(float(maxAR))
 340
 341
 342 #def extendElsets(meshFile):
 343   #if not path.isfile(meshFile):
 344     #message('E','Mesh med file is not valid')
 345     #return(-1)
 346
 347   #import SMESH, salome
 348   ##salome.salome_init()
 349   #theStudy = salome.myStudy
 350   #from salome.smesh import smeshBuilder
 351   #smesh = smeshBuilder.New(theStudy)
 352
 353   #([mesh], status) = smesh.CreateMeshesFromMED(meshFile)
 354
 355   ## Node color status
 356   #nodeList=mesh.GetNodesId()
 357   #colorList=[0]*len(nodeList)
 358
 359   ## Init using SIDE0 SIDE1
 360   #for group in mesh.GetGroups():
 361     #if group.GetType()==SMESH.FACE :
 362       #color=0
 363       #if group.GetName()[0:4]=='SIDE0' :
 364         #color=1
 365       #elif group.GetName()[0:4]=='SIDE1' :
 366         #color=2
 367       #else : continue
 368       ## Get faces
 369       #faces=group.GetIDs()
 370       ## Set faces nodes to given color
 371       #for face_id in faces :
 372         #for face_node_id in mesh.GetElemNodes(face_id) :
 373           #colorList[face_node_id-1]=color
 374
 375   ## Propagates color using elem connectivity
 376   ## Always propagates max color
 377   #volElemList=mesh.GetElementsByType(SMESH.VOLUME)
 378   #ifChanged=True
 379   #while ifChanged :
 380     #ifChanged=False
 381     #minColor=100
 382     #maxColor=0
 383     #for elemId in volElemList:
 384       #for elemNodeId in mesh.GetElemNodes(elemId) :
 385         #nodeColor=colorList[elemNodeId-1]
 386         #if nodeColor<minColor : minColor=nodeColor
 387         #if nodeColor>maxColor : maxColor=nodeColor
 388       #if minColor!=maxColor :
 389         #ifChanged = True
 390         #for elemNodeId in mesh.GetElemNodes(elemId) :
 391           #colorList[elemNodeId-1]=maxColor
 392
 393   #velem0 = []
 394   #velem1 = []
 395   #for elemId in volElemList:
 396     #elemNodesId=mesh.GetElemNodes(elemId)
 397     #elemColor=colorList[elemNodesId[0]-1]
 398     #if(elemColor==1) : velem0.append(elemId)
 399     #if(elemColor==2) : velem1.append(elemId)
 400
 401   #mesh.MakeGroupByIds('SIDE_co',SMESH.VOLUME,velem0)
 402   #mesh.MakeGroupByIds('SIDE_ext',SMESH.VOLUME,velem1)
 403
 404   #surfElemList=mesh.GetElementsByType(SMESH.FACE)
 405   #selem0 = []
 406   #selem1 = []
 407   #nbelem0=0
 408   #nbelem1=0
 409
 410   #for elemId in surfElemList:
 411     #elemNodesId=mesh.GetElemNodes(elemId)
 412     #elemColor=colorList[elemNodesId[0]-1]
 413     #if(elemColor==1) : selem0.append(elemId)
 414     #if(elemColor==2) : selem1.append(elemId)
 415
 416   #mesh.MakeGroupByIds('SIDE_co',SMESH.FACE,selem0)
 417   #mesh.MakeGroupByIds('SIDE_ext',SMESH.FACE,selem1)
 418
 419   #maxAR=0.
 420   #for elem in volElemList:
 421     #maxAR=max(mesh.GetAspectRatio(elem),maxAR)
 422   #for elem in surfElemList:
 423     #maxAR=max(mesh.GetAspectRatio(elem),maxAR)
 424
 425   #mesh.ExportMED(meshFile, 0, SMESH.MED_V2_2, 1, None ,1)
 426   #return(maxAR)
 427
 428
 429 def removeFromSessionPath(envVar, patern):
 430   if type(patern) is not list: patern=[patern]
 431   if type(envVar) is not list: envVar=[envVar]
 432
 433   for env in envVar:
 434     path=environ[env]
 435     listPath=path.split(':')
 436     for p in listPath:
 437       for pat in patern:
 438         if pat in p:
 439           path=path.replace(p,'')
 440           path.replace('::',':')
 441     environ[env]=path
 442
 443
 444 #def isPlane(geomObject, eps=1.e-9):
 445   #import salome
 446   #salome.salome_init()
 447   #theStudy = salome.myStudy
 448
 449   #import salome_notebook
 450   #notebook = salome_notebook.NoteBook(theStudy)
 451
 452   #import GEOM
 453   #from salome.geom import geomBuilder
 454   #geompy = geomBuilder.New(theStudy)
 455
 456   #Vs=geompy.SubShapeAll(geomObject, geompy.ShapeType["VERTEX"])
 457   #if len(Vs)<=3:
 458     #return(True)
 459   #elif len(Vs)>3:
 460     #P0=numpy.array(geompy.GetPosition(Vs[0])[:3])
 461     #P1=numpy.array(geompy.GetPosition(Vs[1])[:3])
 462     #P2=numpy.array(geompy.GetPosition(Vs[2])[:3])
 463     #V01=P1-P0
 464     #V02=P2-P0
 465     #V12=P2-P1
 466     #norm01=numpy.linalg.norm(V01)
 467     #norm02=numpy.linalg.norm(V02)
 468     #norm12=numpy.linalg.norm(V12)
 469     #if (norm01<eps) or (norm02<eps) or (norm12<eps):
 470       #print 'error'
 471       #return(False)
 472     #else:
 473       #N=numpy.cross(V01,V02)
 474       #N=N/numpy.linalg.norm(N)
 475       #maxDist=0.
 476       #for P in Vs[3:]:
 477         #Pi=numpy.array(geompy.GetPosition(P)[:3])
 478         #V=Pi-P0
 479         #d=numpy.dot(V,N)
 480         #maxDist=numpy.max([maxDist,numpy.abs(d)])
 481   #else:
 482     #print 'error'
 483     #return(False)
 484
 485   #return(maxDist<eps)