1 # Copyright (C) 2016-2024 CEA, EDF
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.
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.
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
17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 #sys.path.append('/home/I60976/00_PROJETS/2015_INTEGRATION_ZCRACKS/zcracks_salome/zcracks')
23 import numpy, subprocess, sys
24 from os import remove, getpid, path, environ
25 from .output import message
27 def calcCoordVectors(normalIN, directionIN):
28 V3TEMP=numpy.cross(normalIN,directionIN)
29 directionTEMP=numpy.cross(V3TEMP,normalIN)
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)
37 def testStrictRange(x, inf=0.0, sup=False):
39 c1=(isinstance(x, list))
43 c3=(isinstance(x[0], type(inf)))
55 c1=(isinstance(x, list))
59 c3=(isinstance(x[0], float))
60 c4=(isinstance(x[1], float))
61 c5=(isinstance(x[2], float))
66 def testRange(x, inf=0.0, sup=False):
68 c1=(isinstance(x, list))
72 c3=(isinstance(x[0], type(inf)))
86 c1=(data['crackedName']!='')
90 message('E','Invalid Cracked name',goOn=True)
94 c1=path.isfile(data['saneName'])
96 c2=(data['saneName']!=data['crackedName'])
100 message('E','sane mesh and cracked mesh are identical',goOn=True)
103 message('E','Bad sane mesh file',goOn=True)
106 test=testStrictRange(data['minSize'])
108 message('E','invalid min size',goOn=True)
111 test=testStrictRange(data['maxSize'])
113 message('E','invalid max size',goOn=True)
117 test=(data['maxSize'][0]>=data['minSize'][0])
119 message('E','min size greater than max size',goOn=True)
122 test=testStrictRange(data['extractLength'])
124 message('E','invalid extract length',goOn=True)
127 test=testRange(data['gradation'], inf=1.0)
129 message('E','invalid Gradation',goOn=True)
132 test=testRange(data['layers'], inf=1)
134 message('E','invalid layers',goOn=True)
137 test=testRange(data['iterations'], inf=1)
139 message('E','invalid iterations',goOn=True)
144 def calcElemSize(A, R):
145 x=R*(1.-numpy.cos(A/2.))
149 def meshCrack(geomObject, minSize, maxSize, chordal, dim):
154 import SMESH, SALOMEDS
155 from salome.smesh import smeshBuilder
156 smesh = smeshBuilder.New()
157 Maillage = smesh.Mesh(geomObject)
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 )
169 Regular_1D = Maillage.Segment()
170 Adaptive = Regular_1D.Adaptive(minSize,maxSize,chordal)
171 NETGEN_2D_ONLY = Maillage.Triangle(algo=smeshBuilder.NETGEN_2D)
173 message('E',"error in mesh dimension",goOn=True)
175 isDone = Maillage.Compute()
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() )
184 def extendElsets(meshFile, outFile=None):
186 if outFile==None: outFile=meshFile
188 if not path.isfile(meshFile):
189 message('E','Mesh med file is not valid')
193 #salome.salome_init()
194 from salome.smesh import smeshBuilder
195 smesh = smeshBuilder.New()
197 ([mesh], status) = smesh.CreateMeshesFromMED(meshFile)
199 if mesh.NbVolumes()>0:
201 mesh.Reorient2DBy3D( [ mesh ], mesh, 1 )
205 mesh=cleanGroups(mesh)
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)
215 for group in mesh.GetGroups():
216 if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides' : case2D=False
219 for group in mesh.GetGroups():
221 if group.GetType()==SMESH.FACE and group.GetName()[:5]=='sides':
224 if group.GetType()==SMESH.VOLUME and group.GetName()[:5]=='sides':
227 sortedSides=[None]*len(sides)
229 N=group.GetName().replace('sides','').replace('_bset','').replace(' ','')
235 for elemNodeId in mesh.GetElemNodes(elemId) :
236 colorList[elemNodeId-1]=N
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()
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
253 mesh.ExportMED(outFile)
256 # Propagates color using elem connectivity
257 # Always propagates max color
259 #elemToTreat=volElemList
261 #while len(elemToTreat)>0 :
262 #print len(elemToTreat)
263 #for elemId in elemToTreat:
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
277 elemList=[surfElemList,edgeElemList]
280 elemList=[volElemList,surfElemList,edgeElemList]
281 grElemList=[[],[],[]]
285 for elemId in elemList[0]:
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 :
294 for elemNodeId in mesh.GetElemNodes(elemId) :
295 colorList[elemNodeId-1]=maxColor
298 for x in range(len(sides)):
301 for N, el in enumerate(elemList):
303 elemNodesId=mesh.GetElemNodes(elemId)
304 elemColor=colorList[elemNodesId[0]-1]
306 grElemList[N][elemColor].append(elemId)
308 #for elemId in surfElemList:
309 #elemNodesId=mesh.GetElemNodes(elemId)
310 #elemColor=colorList[elemNodesId[0]-1]
312 #selem[elemColor].append(elemId)
314 for n in range(len(sides)):
316 mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.FACE,grElemList[0][n])
317 mesh.MakeGroupByIds('Extended_side%d' %n ,SMESH.EDGE,grElemList[1][n])
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])
323 if outFile==None: outFile=meshFile
324 mesh.ExportMED(outFile)
328 def cleanGroups(mesh):
330 for group in mesh.GetGroups():
331 if '_bset' in group.GetName():
332 group.SetName(group.GetName().replace('_bset',''))
334 if group.GetType()==SMESH.NODE:
335 if group.GetName() in ['SURFACE','lip','SFRONT_NODES','FRONT']: mesh.RemoveGroup(group)
337 #elif group.GetType()==SMESH.EDGE:
339 elif group.GetType()==SMESH.FACE:
340 if group.GetName() in ['SURFACE','Nlip']:
341 mesh.RemoveGroup(group)
343 elif group.GetType()==SMESH.VOLUME:
344 if (group.GetName() in ['ELSET0','AUTO']) or (group.GetName()[:4] in ['SIDE']) :
345 mesh.RemoveGroup(group)
350 def getMaxAspectRatio(tmpdir):
351 logFile=path.join(tmpdir,'MESHING_OUTPUT')
353 if not path.isfile(logFile): return(-1)
356 f = open(logFile, "r")
358 if re.search("WORST ELEMENT QUALITY", line): maxAR=line
361 for r in [' ','WORSTELEMENTQUALITY','\n']: maxAR=maxAR.replace(r,'')
367 def removeFromSessionPath(envVar, patern):
368 if not isinstance(patern, list): patern=[patern]
369 if not isinstance(envVar, list): envVar=[envVar]
373 listPath=path.split(':')
377 path=path.replace(p,'')
378 path.replace('::',':')
382 #def isPlane(geomObject, eps=1.e-9):
384 #salome.salome_init()
385 #theStudy = salome.myStudy
387 #import salome_notebook
388 #notebook = salome_notebook.NoteBook(theStudy)
391 #from salome.geom import geomBuilder
392 #geompy = geomBuilder.New(theStudy)
394 #Vs=geompy.SubShapeAll(geomObject, geompy.ShapeType["VERTEX"])
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])
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):
411 #N=numpy.cross(V01,V02)
412 #N=N/numpy.linalg.norm(N)
415 #Pi=numpy.array(geompy.GetPosition(P)[:3])
418 #maxDist=numpy.max([maxDist,numpy.abs(d)])