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):
20 c1=(isinstance(x, list))
24 c3=(isinstance(x[0], type(inf)))
36 c1=(isinstance(x, list))
40 c3=(isinstance(x[0], float))
41 c4=(isinstance(x[1], float))
42 c5=(isinstance(x[2], float))
47 def testRange(x, inf=0.0, sup=False):
49 c1=(isinstance(x, list))
53 c3=(isinstance(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):
135 import SMESH, SALOMEDS
136 from salome.smesh import smeshBuilder
137 smesh = smeshBuilder.New()
138 Maillage = smesh.Mesh(geomObject)
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 )
150 Regular_1D = Maillage.Segment()
151 Adaptive = Regular_1D.Adaptive(minSize,maxSize,chordal)
152 NETGEN_2D_ONLY = Maillage.Triangle(algo=smeshBuilder.NETGEN_2D)
154 message('E',"error in mesh dimension",goOn=True)
156 isDone = Maillage.Compute()
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() )
165 def extendElsets(meshFile, outFile=None):
167 if outFile==None: outFile=meshFile
169 if not path.isfile(meshFile):
170 message('E','Mesh med file is not valid')
174 #salome.salome_init()
175 from salome.smesh import smeshBuilder
176 smesh = smeshBuilder.New()
178 ([mesh], status) = smesh.CreateMeshesFromMED(meshFile)
180 if mesh.NbVolumes()>0:
182 mesh.Reorient2DBy3D( [ mesh ], mesh, 1 )
186 mesh=cleanGroups(mesh)
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)
196 for group in mesh.GetGroups():
197 if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides' : case2D=False
200 for group in mesh.GetGroups():
202 if group.GetType()==SMESH.FACE and group.GetName()[:5]=='sides':
205 if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides':
208 sortedSides=[None]*len(sides)
210 N=group.GetName().replace('sides','').replace('_bset','').replace(' ','')
216 for elemNodeId in mesh.GetElemNodes(elemId) :
217 colorList[elemNodeId-1]=N
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()
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
234 mesh.ExportMED(outFile)
237 # Propagates color using elem connectivity
238 # Always propagates max color
240 #elemToTreat=volElemList
242 #while len(elemToTreat)>0 :
243 #print len(elemToTreat)
244 #for elemId in elemToTreat:
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
258 elemList=[surfElemList,edgeElemList]
261 elemList=[volElemList,surfElemList,edgeElemList]
262 grElemList=[[],[],[]]
266 for elemId in elemList[0]:
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 :
275 for elemNodeId in mesh.GetElemNodes(elemId) :
276 colorList[elemNodeId-1]=maxColor
279 for x in range(len(sides)):
282 for N, el in enumerate(elemList):
284 elemNodesId=mesh.GetElemNodes(elemId)
285 elemColor=colorList[elemNodesId[0]-1]
287 grElemList[N][elemColor].append(elemId)
289 #for elemId in surfElemList:
290 #elemNodesId=mesh.GetElemNodes(elemId)
291 #elemColor=colorList[elemNodesId[0]-1]
293 #selem[elemColor].append(elemId)
295 for n in range(len(sides)):
297 mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.FACE,grElemList[0][n])
298 mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.EDGE,grElemList[1][n])
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])
304 if outFile==None: outFile=meshFile
305 mesh.ExportMED(outFile)
309 def cleanGroups(mesh):
311 for group in mesh.GetGroups():
312 if '_bset' in group.GetName():
313 group.SetName(group.GetName().replace('_bset',''))
315 if group.GetType()==SMESH.NODE:
316 if group.GetName() in ['SURFACE','lip','SFRONT_NODES','FRONT']: mesh.RemoveGroup(group)
318 #elif group.GetType()==SMESH.EDGE:
320 elif group.GetType()==SMESH.FACE:
321 if group.GetName() in ['SURFACE','Nlip']:
322 mesh.RemoveGroup(group)
324 elif group.GetType()==SMESH.VOLUME:
325 if (group.GetName() in ['ELSET0','AUTO']) or (group.GetName()[:4] in ['SIDE']) :
326 mesh.RemoveGroup(group)
331 def getMaxAspectRatio(tmpdir):
332 logFile=path.join(tmpdir,'MESHING_OUTPUT')
334 if not path.isfile(logFile): return(-1)
337 f = open(logFile, "r")
339 if re.search("WORST ELEMENT QUALITY", line): maxAR=line
342 for r in [' ','WORSTELEMENTQUALITY','\n']: maxAR=maxAR.replace(r,'')
348 def removeFromSessionPath(envVar, patern):
349 if not isinstance(patern, list): patern=[patern]
350 if not isinstance(envVar, list): envVar=[envVar]
354 listPath=path.split(':')
358 path=path.replace(p,'')
359 path.replace('::',':')
363 #def isPlane(geomObject, eps=1.e-9):
365 #salome.salome_init()
366 #theStudy = salome.myStudy
368 #import salome_notebook
369 #notebook = salome_notebook.NoteBook(theStudy)
372 #from salome.geom import geomBuilder
373 #geompy = geomBuilder.New(theStudy)
375 #Vs=geompy.SubShapeAll(geomObject, geompy.ShapeType["VERTEX"])
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])
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):
392 #N=numpy.cross(V01,V02)
393 #N=N/numpy.linalg.norm(N)
396 #Pi=numpy.array(geompy.GetPosition(P)[:3])
399 #maxDist=numpy.max([maxDist,numpy.abs(d)])