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   if mesh.NbVolumes()>0:
 183     case2D=False
 184     mesh.Reorient2DBy3D( [ mesh ], mesh, 1 )
 185   else:
 186     case2D=True
 187
 188   mesh=cleanGroups(mesh)
 189
 190   # Node color status
 191   nodeList=mesh.GetNodesId()
 192   volElemList=mesh.GetElementsByType(SMESH.VOLUME)
 193   surfElemList=mesh.GetElementsByType(SMESH.FACE)
 194   edgeElemList=mesh.GetElementsByType(SMESH.EDGE)
 195   colorList=[-1]*len(nodeList)
 196
 197   case2D=True
 198   for group in mesh.GetGroups():
 199     if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides' : case2D=False
 200
 201   sides=[]
 202   for group in mesh.GetGroups():
 203     if case2D:
 204       if group.GetType()==SMESH.FACE and group.GetName()[:5]=='sides':
 205         sides.append(group)
 206     else:
 207       if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides':
 208         sides.append(group)
 209
 210   sortedSides=[None]*len(sides)
 211   for group in sides:
 212     N=group.GetName().replace('sides','').replace('_bset','').replace(' ','')
 213     N=int(N)
 214     sortedSides[N]=group
 215
 216     elems=group.GetIDs()
 217     for elemId in elems:
 218       for elemNodeId in mesh.GetElemNodes(elemId) :
 219         colorList[elemNodeId-1]=N
 220   #print colorList
 221
 222   crackOnly=True
 223   for iN in range(len(sides)/2):
 224     side0=sortedSides[2*iN]
 225     side1=sortedSides[2*iN+1]
 226     elemsOfside0=side0.GetIDs()
 227     elemsOfside1=side1.GetIDs()
 228     NodesOfside0=[]
 229     NodesOfside1=[]
 230     for elem in elemsOfside0: NodesOfside0+=mesh.GetElemNodes(elem)
 231     for elem in elemsOfside1: NodesOfside1+=mesh.GetElemNodes(elem)
 232     front=set(NodesOfside0).intersection(set(NodesOfside1))
 233     if len(front)==0: crackOnly=False
 234
 235   if crackOnly:
 236     mesh.ExportMED(outFile, 0, SMESH.MED_V2_2, 1, None ,1)
 237     return('crack')
 238
 239   # Propagates color using elem connectivity
 240   # Always propagates max color
 241
 242   #elemToTreat=volElemList
 243
 244   #while len(elemToTreat)>0 :
 245     #print len(elemToTreat)
 246     #for elemId in elemToTreat:
 247       #minColor=sys.maxint
 248       #maxColor=-sys.maxint
 249       #for elemNodeId in mesh.GetElemNodes(elemId) :
 250         #nodeColor=colorList[elemNodeId-1]
 251         #if nodeColor<minColor : minColor=nodeColor
 252         #if nodeColor>maxColor : maxColor=nodeColor
 253       #if minColor!=maxColor :
 254         #elemToTreat.remove(elemId)
 255         #for elemNodeId in mesh.GetElemNodes(elemId) :
 256           #colorList[elemNodeId-1]=maxColor
 257
 258   ifChanged=True
 259   if case2D:
 260     elemList=[surfElemList,edgeElemList]
 261     grElemList=[[],[]]
 262   else:
 263     elemList=[volElemList,surfElemList,edgeElemList]
 264     grElemList=[[],[],[]]
 265
 266   while ifChanged :
 267     ifChanged=False
 268     for elemId in elemList[0]:
 269       minColor=sys.maxint
 270       maxColor=-sys.maxint
 271       for elemNodeId in mesh.GetElemNodes(elemId) :
 272         nodeColor=colorList[elemNodeId-1]
 273         if nodeColor<minColor : minColor=nodeColor
 274         if nodeColor>maxColor : maxColor=nodeColor
 275       if minColor!=maxColor :
 276         ifChanged = True
 277         for elemNodeId in mesh.GetElemNodes(elemId) :
 278           colorList[elemNodeId-1]=maxColor
 279
 280   for l in grElemList:
 281     for x in range(len(sides)):
 282       l.append([])
 283
 284   for N, el in enumerate(elemList):
 285     for elemId in el:
 286       elemNodesId=mesh.GetElemNodes(elemId)
 287       elemColor=colorList[elemNodesId[0]-1]
 288       if elemColor>=0:
 289         grElemList[N][elemColor].append(elemId)
 290
 291   #for elemId in surfElemList:
 292     #elemNodesId=mesh.GetElemNodes(elemId)
 293     #elemColor=colorList[elemNodesId[0]-1]
 294     #if elemColor>=0:
 295       #selem[elemColor].append(elemId)
 296
 297   for n in range(len(sides)):
 298     if case2D:
 299       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.FACE,grElemList[0][n])
 300       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.EDGE,grElemList[1][n])
 301     else:
 302       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.VOLUME,grElemList[0][n])
 303       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.FACE,grElemList[1][n])
 304       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.EDGE,grElemList[2][n])
 305
 306   if outFile==None: outFile=meshFile
 307   mesh.ExportMED(outFile, 0, SMESH.MED_V2_2, 1, None ,1)
 308   return(True)
 309
 310
 311 def cleanGroups(mesh):
 312   import SMESH
 313   for group in mesh.GetGroups():
 314     if '_bset' in group.GetName():
 315       group.SetName(group.GetName().replace('_bset',''))
 316
 317     if group.GetType()==SMESH.NODE:
 318       if group.GetName() in ['SURFACE','lip','SFRONT_NODES','FRONT']: mesh.RemoveGroup(group)
 319
 320     #elif group.GetType()==SMESH.EDGE:
 321
 322     elif group.GetType()==SMESH.FACE:
 323       if group.GetName() in ['SURFACE','Nlip']:
 324         mesh.RemoveGroup(group)
 325
 326     elif group.GetType()==SMESH.VOLUME:
 327       if (group.GetName() in ['ELSET0','AUTO']) or (group.GetName()[:4] in ['SIDE']) :
 328         mesh.RemoveGroup(group)
 329
 330   return(mesh)
 331
 332
 333 def getMaxAspectRatio(tmpdir):
 334   logFile=path.join(tmpdir,'MESHING_OUTPUT')
 335   print logFile
 336   if not path.isfile(logFile): return(-1)
 337
 338   import re
 339   f = open(logFile, "r")
 340   for line in f:
 341     if re.search("WORST ELEMENT QUALITY", line):  maxAR=line
 342
 343   f.close()
 344   for r in [' ','WORSTELEMENTQUALITY','\n']: maxAR=maxAR.replace(r,'')
 345   return(float(maxAR))
 346
 347
 348
 349
 350 def removeFromSessionPath(envVar, patern):
 351   if type(patern) is not list: patern=[patern]
 352   if type(envVar) is not list: envVar=[envVar]
 353
 354   for env in envVar:
 355     path=environ[env]
 356     listPath=path.split(':')
 357     for p in listPath:
 358       for pat in patern:
 359         if pat in p:
 360           path=path.replace(p,'')
 361           path.replace('::',':')
 362     environ[env]=path
 363
 364
 365 #def isPlane(geomObject, eps=1.e-9):
 366   #import salome
 367   #salome.salome_init()
 368   #theStudy = salome.myStudy
 369
 370   #import salome_notebook
 371   #notebook = salome_notebook.NoteBook(theStudy)
 372
 373   #import GEOM
 374   #from salome.geom import geomBuilder
 375   #geompy = geomBuilder.New(theStudy)
 376
 377   #Vs=geompy.SubShapeAll(geomObject, geompy.ShapeType["VERTEX"])
 378   #if len(Vs)<=3:
 379     #return(True)
 380   #elif len(Vs)>3:
 381     #P0=numpy.array(geompy.GetPosition(Vs[0])[:3])
 382     #P1=numpy.array(geompy.GetPosition(Vs[1])[:3])
 383     #P2=numpy.array(geompy.GetPosition(Vs[2])[:3])
 384     #V01=P1-P0
 385     #V02=P2-P0
 386     #V12=P2-P1
 387     #norm01=numpy.linalg.norm(V01)
 388     #norm02=numpy.linalg.norm(V02)
 389     #norm12=numpy.linalg.norm(V12)
 390     #if (norm01<eps) or (norm02<eps) or (norm12<eps):
 391       #print 'error'
 392       #return(False)
 393     #else:
 394       #N=numpy.cross(V01,V02)
 395       #N=N/numpy.linalg.norm(N)
 396       #maxDist=0.
 397       #for P in Vs[3:]:
 398         #Pi=numpy.array(geompy.GetPosition(P)[:3])
 399         #V=Pi-P0
 400         #d=numpy.dot(V,N)
 401         #maxDist=numpy.max([maxDist,numpy.abs(d)])
 402   #else:
 403     #print 'error'
 404     #return(False)
 405
 406   #return(maxDist<eps)