2 #sys.path.append('/home/I60976/00_PROJETS/2015_INTEGRATION_ZCRACKS/zcracks_salome/zcracks')
4 import numpy, subprocess, sys
5 from os import remove, getpid, path, environ
6 from output import message
8 def calcCoordVectors(normalIN, directionIN):
9 V3TEMP=numpy.cross(normalIN,directionIN)
10 directionTEMP=numpy.cross(V3TEMP,normalIN)
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)
18 def testStrictRange(x, inf=0.0, sup=False):
24 c3=(type(x[0])==type(inf))
40 c3=(type(x[0])==float)
41 c4=(type(x[1])==float)
42 c5=(type(x[2])==float)
47 def testRange(x, inf=0.0, sup=False):
53 c3=(type(x[0])==type(inf))
67 c1=(data['crackedName']!='')
71 message('E','Invalid Cracked name',goOn=True)
75 c1=path.isfile(data['saneName'])
77 c2=(data['saneName']!=data['crackedName'])
81 message('E','sane mesh and cracked mesh are identical',goOn=True)
84 message('E','Bad sane mesh file',goOn=True)
87 test=testStrictRange(data['minSize'])
89 message('E','invalid min size',goOn=True)
92 test=testStrictRange(data['maxSize'])
94 message('E','invalid max size',goOn=True)
98 test=(data['maxSize'][0]>=data['minSize'][0])
100 message('E','min size greater than max size',goOn=True)
103 test=testStrictRange(data['extractLength'])
105 message('E','invalid extract length',goOn=True)
108 test=testRange(data['gradation'], inf=1.0)
110 message('E','invalid Gradation',goOn=True)
113 test=testRange(data['layers'], inf=1)
115 message('E','invalid layers',goOn=True)
118 test=testRange(data['iterations'], inf=1)
120 message('E','invalid iterations',goOn=True)
125 def calcElemSize(A, R):
126 x=R*(1.-numpy.cos(A/2.))
130 def meshCrack(geomObject, minSize, maxSize, chordal, dim):
134 theStudy = salome.myStudy
136 import SMESH, SALOMEDS
137 from salome.smesh import smeshBuilder
138 smesh = smeshBuilder.New(theStudy)
139 Maillage = smesh.Mesh(geomObject)
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 )
151 Regular_1D = Maillage.Segment()
152 Adaptive = Regular_1D.Adaptive(minSize,maxSize,chordal)
153 NETGEN_2D_ONLY = Maillage.Triangle(algo=smeshBuilder.NETGEN_2D)
155 message('E',"error in mesh dimension",goOn=True)
157 isDone = Maillage.Compute()
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() )
166 def extendElsets(meshFile, outFile=None):
168 if outFile==None: outFile=meshFile
170 if not path.isfile(meshFile):
171 message('E','Mesh med file is not valid')
175 #salome.salome_init()
176 theStudy = salome.myStudy
177 from salome.smesh import smeshBuilder
178 smesh = smeshBuilder.New(theStudy)
180 ([mesh], status) = smesh.CreateMeshesFromMED(meshFile)
182 if mesh.NbVolumes()>0:
184 mesh.Reorient2DBy3D( [ mesh ], mesh, 1 )
188 mesh=cleanGroups(mesh)
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)
198 for group in mesh.GetGroups():
199 if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides' : case2D=False
202 for group in mesh.GetGroups():
204 if group.GetType()==SMESH.FACE and group.GetName()[:5]=='sides':
207 if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides':
210 sortedSides=[None]*len(sides)
212 N=group.GetName().replace('sides','').replace('_bset','').replace(' ','')
218 for elemNodeId in mesh.GetElemNodes(elemId) :
219 colorList[elemNodeId-1]=N
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()
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
236 mesh.ExportMED(outFile, 0, SMESH.MED_V2_2, 1, None ,1)
239 # Propagates color using elem connectivity
240 # Always propagates max color
242 #elemToTreat=volElemList
244 #while len(elemToTreat)>0 :
245 #print len(elemToTreat)
246 #for elemId in elemToTreat:
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
260 elemList=[surfElemList,edgeElemList]
263 elemList=[volElemList,surfElemList,edgeElemList]
264 grElemList=[[],[],[]]
268 for elemId in elemList[0]:
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 :
277 for elemNodeId in mesh.GetElemNodes(elemId) :
278 colorList[elemNodeId-1]=maxColor
281 for x in range(len(sides)):
284 for N, el in enumerate(elemList):
286 elemNodesId=mesh.GetElemNodes(elemId)
287 elemColor=colorList[elemNodesId[0]-1]
289 grElemList[N][elemColor].append(elemId)
291 #for elemId in surfElemList:
292 #elemNodesId=mesh.GetElemNodes(elemId)
293 #elemColor=colorList[elemNodesId[0]-1]
295 #selem[elemColor].append(elemId)
297 for n in range(len(sides)):
299 mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.FACE,grElemList[0][n])
300 mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.EDGE,grElemList[1][n])
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])
306 if outFile==None: outFile=meshFile
307 mesh.ExportMED(outFile, 0, SMESH.MED_V2_2, 1, None ,1)
311 def cleanGroups(mesh):
313 for group in mesh.GetGroups():
314 if '_bset' in group.GetName():
315 group.SetName(group.GetName().replace('_bset',''))
317 if group.GetType()==SMESH.NODE:
318 if group.GetName() in ['SURFACE','lip','SFRONT_NODES','FRONT']: mesh.RemoveGroup(group)
320 #elif group.GetType()==SMESH.EDGE:
322 elif group.GetType()==SMESH.FACE:
323 if group.GetName() in ['SURFACE','Nlip']:
324 mesh.RemoveGroup(group)
326 elif group.GetType()==SMESH.VOLUME:
327 if (group.GetName() in ['ELSET0','AUTO']) or (group.GetName()[:4] in ['SIDE']) :
328 mesh.RemoveGroup(group)
333 def getMaxAspectRatio(tmpdir):
334 logFile=path.join(tmpdir,'MESHING_OUTPUT')
336 if not path.isfile(logFile): return(-1)
339 f = open(logFile, "r")
341 if re.search("WORST ELEMENT QUALITY", line): maxAR=line
344 for r in [' ','WORSTELEMENTQUALITY','\n']: maxAR=maxAR.replace(r,'')
350 def removeFromSessionPath(envVar, patern):
351 if type(patern) is not list: patern=[patern]
352 if type(envVar) is not list: envVar=[envVar]
356 listPath=path.split(':')
360 path=path.replace(p,'')
361 path.replace('::',':')
365 #def isPlane(geomObject, eps=1.e-9):
367 #salome.salome_init()
368 #theStudy = salome.myStudy
370 #import salome_notebook
371 #notebook = salome_notebook.NoteBook(theStudy)
374 #from salome.geom import geomBuilder
375 #geompy = geomBuilder.New(theStudy)
377 #Vs=geompy.SubShapeAll(geomObject, geompy.ShapeType["VERTEX"])
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])
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):
394 #N=numpy.cross(V01,V02)
395 #N=N/numpy.linalg.norm(N)
398 #Pi=numpy.array(geompy.GetPosition(P)[:3])
401 #maxDist=numpy.max([maxDist,numpy.abs(d)])