src/Tools/ZCracksPlug/utilityFunctions.py

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