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=(isinstance(x, list))
  21   if c1:
  22     c2=(len(x)==1)
  23     if c2:
  24       c3=(isinstance(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=(isinstance(x, list))
  37   if c1:
  38     c2=(len(x)==3)
  39     if c2:
  40       c3=(isinstance(x[0], float))
  41       c4=(isinstance(x[1], float))
  42       c5=(isinstance(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=(isinstance(x, list))
  50   if c1:
  51     c2=(len(x)==1)
  52     if c2:
  53       c3=(isinstance(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
 135   import  SMESH, SALOMEDS
 136   from salome.smesh import smeshBuilder
 137   smesh = smeshBuilder.New()
 138   Maillage = smesh.Mesh(geomObject)
 139
 140   if dim==3:
 141     MG_CADSurf = Maillage.Triangle(algo=smeshBuilder.MG_CADSurf)
 142     MG_CADSurf_Parameters = MG_CADSurf.Parameters()
 143     MG_CADSurf_Parameters.SetPhysicalMesh( 0 )
 144     MG_CADSurf_Parameters.SetGeometricMesh( 1 )
 145     MG_CADSurf_Parameters.SetMinSize( minSize )
 146     MG_CADSurf_Parameters.SetMaxSize( maxSize )
 147     MG_CADSurf_Parameters.SetChordalError( chordal )
 148
 149   elif dim==2:
 150     Regular_1D = Maillage.Segment()
 151     Adaptive = Regular_1D.Adaptive(minSize,maxSize,chordal)
 152     NETGEN_2D_ONLY = Maillage.Triangle(algo=smeshBuilder.NETGEN_2D)
 153   else:
 154     message('E',"error in mesh dimension",goOn=True)
 155
 156   isDone = Maillage.Compute()
 157
 158   #crack1 = Maillage.CreateEmptyGroup( SMESH.NODE, 'crack' )
 159   #nbAdd = crack1.AddFrom( Maillage.GetMesh() )
 160   #crack2 = Maillage.CreateEmptyGroup( SMESH.NODE, 'surface' )
 161   #nbAdd = crack2.AddFrom( Maillage.GetMesh() )
 162
 163   return(Maillage)
 164
 165 def extendElsets(meshFile, outFile=None):
 166
 167   if outFile==None: outFile=meshFile
 168
 169   if not path.isfile(meshFile):
 170     message('E','Mesh med file is not valid')
 171     return('error')
 172
 173   import SMESH, salome
 174   #salome.salome_init()
 175   from salome.smesh import smeshBuilder
 176   smesh = smeshBuilder.New()
 177
 178   ([mesh], status) = smesh.CreateMeshesFromMED(meshFile)
 179
 180   if mesh.NbVolumes()>0:
 181     case2D=False
 182     mesh.Reorient2DBy3D( [ mesh ], mesh, 1 )
 183   else:
 184     case2D=True
 185
 186   mesh=cleanGroups(mesh)
 187
 188   # Node color status
 189   nodeList=mesh.GetNodesId()
 190   volElemList=mesh.GetElementsByType(SMESH.VOLUME)
 191   surfElemList=mesh.GetElementsByType(SMESH.FACE)
 192   edgeElemList=mesh.GetElementsByType(SMESH.EDGE)
 193   colorList=[-1]*len(nodeList)
 194
 195   case2D=True
 196   for group in mesh.GetGroups():
 197     if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides' : case2D=False
 198
 199   sides=[]
 200   for group in mesh.GetGroups():
 201     if case2D:
 202       if group.GetType()==SMESH.FACE and group.GetName()[:5]=='sides':
 203         sides.append(group)
 204     else:
 205       if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides':
 206         sides.append(group)
 207
 208   sortedSides=[None]*len(sides)
 209   for group in sides:
 210     N=group.GetName().replace('sides','').replace('_bset','').replace(' ','')
 211     N=int(N)
 212     sortedSides[N]=group
 213
 214     elems=group.GetIDs()
 215     for elemId in elems:
 216       for elemNodeId in mesh.GetElemNodes(elemId) :
 217         colorList[elemNodeId-1]=N
 218   #print colorList
 219
 220   crackOnly=True
 221   for iN in range(len(sides)/2):
 222     side0=sortedSides[2*iN]
 223     side1=sortedSides[2*iN+1]
 224     elemsOfside0=side0.GetIDs()
 225     elemsOfside1=side1.GetIDs()
 226     NodesOfside0=[]
 227     NodesOfside1=[]
 228     for elem in elemsOfside0: NodesOfside0+=mesh.GetElemNodes(elem)
 229     for elem in elemsOfside1: NodesOfside1+=mesh.GetElemNodes(elem)
 230     front=set(NodesOfside0).intersection(set(NodesOfside1))
 231     if len(front)==0: crackOnly=False
 232
 233   if crackOnly:
 234     mesh.ExportMED(outFile)
 235     return('crack')
 236
 237   # Propagates color using elem connectivity
 238   # Always propagates max color
 239
 240   #elemToTreat=volElemList
 241
 242   #while len(elemToTreat)>0 :
 243     #print len(elemToTreat)
 244     #for elemId in elemToTreat:
 245       #minColor=sys.maxint
 246       #maxColor=-sys.maxint
 247       #for elemNodeId in mesh.GetElemNodes(elemId) :
 248         #nodeColor=colorList[elemNodeId-1]
 249         #if nodeColor<minColor : minColor=nodeColor
 250         #if nodeColor>maxColor : maxColor=nodeColor
 251       #if minColor!=maxColor :
 252         #elemToTreat.remove(elemId)
 253         #for elemNodeId in mesh.GetElemNodes(elemId) :
 254           #colorList[elemNodeId-1]=maxColor
 255
 256   ifChanged=True
 257   if case2D:
 258     elemList=[surfElemList,edgeElemList]
 259     grElemList=[[],[]]
 260   else:
 261     elemList=[volElemList,surfElemList,edgeElemList]
 262     grElemList=[[],[],[]]
 263
 264   while ifChanged :
 265     ifChanged=False
 266     for elemId in elemList[0]:
 267       minColor=sys.maxsize
 268       maxColor=-sys.maxsize
 269       for elemNodeId in mesh.GetElemNodes(elemId) :
 270         nodeColor=colorList[elemNodeId-1]
 271         if nodeColor<minColor : minColor=nodeColor
 272         if nodeColor>maxColor : maxColor=nodeColor
 273       if minColor!=maxColor :
 274         ifChanged = True
 275         for elemNodeId in mesh.GetElemNodes(elemId) :
 276           colorList[elemNodeId-1]=maxColor
 277
 278   for l in grElemList:
 279     for x in range(len(sides)):
 280       l.append([])
 281
 282   for N, el in enumerate(elemList):
 283     for elemId in el:
 284       elemNodesId=mesh.GetElemNodes(elemId)
 285       elemColor=colorList[elemNodesId[0]-1]
 286       if elemColor>=0:
 287         grElemList[N][elemColor].append(elemId)
 288
 289   #for elemId in surfElemList:
 290     #elemNodesId=mesh.GetElemNodes(elemId)
 291     #elemColor=colorList[elemNodesId[0]-1]
 292     #if elemColor>=0:
 293       #selem[elemColor].append(elemId)
 294
 295   for n in range(len(sides)):
 296     if case2D:
 297       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.FACE,grElemList[0][n])
 298       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.EDGE,grElemList[1][n])
 299     else:
 300       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.VOLUME,grElemList[0][n])
 301       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.FACE,grElemList[1][n])
 302       mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.EDGE,grElemList[2][n])
 303
 304   if outFile==None: outFile=meshFile
 305   mesh.ExportMED(outFile)
 306   return(True)
 307
 308
 309 def cleanGroups(mesh):
 310   import SMESH
 311   for group in mesh.GetGroups():
 312     if '_bset' in group.GetName():
 313       group.SetName(group.GetName().replace('_bset',''))
 314
 315     if group.GetType()==SMESH.NODE:
 316       if group.GetName() in ['SURFACE','lip','SFRONT_NODES','FRONT']: mesh.RemoveGroup(group)
 317
 318     #elif group.GetType()==SMESH.EDGE:
 319
 320     elif group.GetType()==SMESH.FACE:
 321       if group.GetName() in ['SURFACE','Nlip']:
 322         mesh.RemoveGroup(group)
 323
 324     elif group.GetType()==SMESH.VOLUME:
 325       if (group.GetName() in ['ELSET0','AUTO']) or (group.GetName()[:4] in ['SIDE']) :
 326         mesh.RemoveGroup(group)
 327
 328   return(mesh)
 329
 330
 331 def getMaxAspectRatio(tmpdir):
 332   logFile=path.join(tmpdir,'MESHING_OUTPUT')
 333   print(logFile)
 334   if not path.isfile(logFile): return(-1)
 335
 336   import re
 337   f = open(logFile, "r")
 338   for line in f:
 339     if re.search("WORST ELEMENT QUALITY", line):  maxAR=line
 340
 341   f.close()
 342   for r in [' ','WORSTELEMENTQUALITY','\n']: maxAR=maxAR.replace(r,'')
 343   return(float(maxAR))
 344
 345
 346
 347
 348 def removeFromSessionPath(envVar, patern):
 349   if not isinstance(patern, list): patern=[patern]
 350   if not isinstance(envVar, list): envVar=[envVar]
 351
 352   for env in envVar:
 353     path=environ[env]
 354     listPath=path.split(':')
 355     for p in listPath:
 356       for pat in patern:
 357         if pat in p:
 358           path=path.replace(p,'')
 359           path.replace('::',':')
 360     environ[env]=path
 361
 362
 363 #def isPlane(geomObject, eps=1.e-9):
 364   #import salome
 365   #salome.salome_init()
 366   #theStudy = salome.myStudy
 367
 368   #import salome_notebook
 369   #notebook = salome_notebook.NoteBook(theStudy)
 370
 371   #import GEOM
 372   #from salome.geom import geomBuilder
 373   #geompy = geomBuilder.New(theStudy)
 374
 375   #Vs=geompy.SubShapeAll(geomObject, geompy.ShapeType["VERTEX"])
 376   #if len(Vs)<=3:
 377     #return(True)
 378   #elif len(Vs)>3:
 379     #P0=numpy.array(geompy.GetPosition(Vs[0])[:3])
 380     #P1=numpy.array(geompy.GetPosition(Vs[1])[:3])
 381     #P2=numpy.array(geompy.GetPosition(Vs[2])[:3])
 382     #V01=P1-P0
 383     #V02=P2-P0
 384     #V12=P2-P1
 385     #norm01=numpy.linalg.norm(V01)
 386     #norm02=numpy.linalg.norm(V02)
 387     #norm12=numpy.linalg.norm(V12)
 388     #if (norm01<eps) or (norm02<eps) or (norm12<eps):
 389       #print 'error'
 390       #return(False)
 391     #else:
 392       #N=numpy.cross(V01,V02)
 393       #N=N/numpy.linalg.norm(N)
 394       #maxDist=0.
 395       #for P in Vs[3:]:
 396         #Pi=numpy.array(geompy.GetPosition(P)[:3])
 397         #V=Pi-P0
 398         #d=numpy.dot(V,N)
 399         #maxDist=numpy.max([maxDist,numpy.abs(d)])
 400   #else:
 401     #print 'error'
 402     #return(False)
 403
 404   #return(maxDist<eps)