812 lines
25 KiB
Python
812 lines
25 KiB
Python
|
# $Id$
|
||
|
|
# This tool converts a msc.marc result file into the vtk format that
|
||
|
|
# can be viewed by Paraview software (Kitware), or MayaVi (needs xml-vtk, or ...
|
||
|
|
#
|
||
|
|
# About the vtk format: http://www.vtk.org/VTK/project/about.html
|
||
|
|
# Some example vtk files: http://people.sc.fsu.edu/~jburkardt/data/vtk/vtk.html
|
||
|
|
# www.paraview.org
|
||
|
|
|
||
|
|
import os,sys,math,time,re
|
||
|
|
# python external
|
||
|
|
try:
|
||
|
|
import numpy as N
|
||
|
|
import numpy
|
||
|
|
except:
|
||
|
|
print('Could not import numpy.')
|
||
|
|
|
||
|
|
import py_post # MSC closed source module to access marc result files
|
||
|
|
|
||
|
|
class MARC_POST():
|
||
|
|
import re
|
||
|
|
def __init__(self):
|
||
|
|
self.projdir='./'
|
||
|
|
|
||
|
|
def opent16(self,incr=None):
|
||
|
|
self.fpath=os.path.join(self.projdir,self.postname)
|
||
|
|
print 'Trying to open ',self.fpath,' ...'
|
||
|
|
self.p=py_post.post_open(self.fpath)
|
||
|
|
if self.p==None:
|
||
|
|
print 'Could not open %s.'%self.postname; #return 'err'#; sys.exit(1)
|
||
|
|
raise Exception('Could not open t16')
|
||
|
|
print 'Postfile %s%s is open ...'%(self.projdir,self.postname)
|
||
|
|
self.maxincr=self.p.increments()
|
||
|
|
print 'and has %i increments'%self.maxincr
|
||
|
|
if incr is None:
|
||
|
|
self.p.moveto(self.maxincr-1)
|
||
|
|
else:
|
||
|
|
self.p.moveto(incr+1)
|
||
|
|
print 'moved to increment ', self.p.increment
|
||
|
|
self.p.extrapolation('translate') # linear, translate, average. query with p.extrapolate
|
||
|
|
print 'extrapolation method is ', self.p.extrapolate
|
||
|
|
print 'collecting model information'
|
||
|
|
self.t16info(printFlag=0)
|
||
|
|
print 't16 is open'
|
||
|
|
self.p
|
||
|
|
|
||
|
|
def t16info(self, printFlag=1
|
||
|
|
):
|
||
|
|
if not self.p:
|
||
|
|
self.p=self.opent16()#post_open(self.postname)
|
||
|
|
print self.p
|
||
|
|
oldincr=self.p.position
|
||
|
|
if oldincr==0:
|
||
|
|
self.p.moveto(1)
|
||
|
|
|
||
|
|
self.nnodes=self.p.nodes() #self.p.node(self.nnodes) crashes; only to nnodes-1 possible
|
||
|
|
self.nodes=range(0,self.nnodes)
|
||
|
|
self.nscals=self.p.node_scalars(); #print 'nscals', nscals
|
||
|
|
self.nscal_list=['']*self.nscals
|
||
|
|
self.nel=self.p.elements()
|
||
|
|
self.elscals=self.p.element_scalars()
|
||
|
|
self.elscal_list=['']*self.elscals
|
||
|
|
self.eltens=self.p.element_tensors()
|
||
|
|
self.elten_list=['']*self.eltens
|
||
|
|
for i in range(0,self.nscals):
|
||
|
|
self.nscal_list[i]=self.p.node_scalar_label(i)
|
||
|
|
for i in range (0,self.elscals):
|
||
|
|
self.elscal_list[i]=self.p.element_scalar_label(i)
|
||
|
|
if printFlag==1: print i, self.elscal_list[i]
|
||
|
|
for i in range (0,self.eltens):
|
||
|
|
self.elten_list[i]=self.p.element_tensor_label(i)
|
||
|
|
if printFlag==1: print i, self.elten_list[i]
|
||
|
|
for i in range(0,self.p.element_tensors()):
|
||
|
|
if printFlag==1: print 'Element Tensor: ', i, self.p.element_tensor_label(i)
|
||
|
|
if printFlag==1: print ''
|
||
|
|
for i in range(0,self.p.element_scalars()):
|
||
|
|
if printFlag==1: print 'Element Scalar: ', i, self.p.element_scalar_label(i)
|
||
|
|
if oldincr==0:
|
||
|
|
self.p.moveto(0)
|
||
|
|
|
||
|
|
def closet16(self):
|
||
|
|
print('Trying to close FEM result file ...')
|
||
|
|
try:
|
||
|
|
if self.p:
|
||
|
|
self.p.close()
|
||
|
|
print('FEM result file closed.')
|
||
|
|
self.p=None
|
||
|
|
else:
|
||
|
|
print('post object not open?')
|
||
|
|
except:
|
||
|
|
print 'ERROR. Could not close FEM result file.'
|
||
|
|
|
||
|
|
def getLabelNr(self, label=None, type='Scalar'):
|
||
|
|
if type[0]=='S' or type[0]=='s': # element scalar
|
||
|
|
labelNr=self.elscal_list.index(label)
|
||
|
|
elif type[0]=='N': # node scalar
|
||
|
|
labelNr=self.nscal_list.index(label)
|
||
|
|
elif type[0]=='T' or type[0]=='t': # tensor
|
||
|
|
labelNr=self.elten_list.index(label)
|
||
|
|
print('Found label %s at index %i'%(label,labelNr))
|
||
|
|
return labelNr
|
||
|
|
|
||
|
|
|
||
|
|
def writeNodes2VTK(self, fobj):
|
||
|
|
self.points=[]
|
||
|
|
self.VTKcnt=200 # number of values per line in vtk file
|
||
|
|
ndCnt=1
|
||
|
|
fobj.write('POINTS %i'%self.p.nodes()+' float\n')
|
||
|
|
self.nodes_dict={} # store the node IDs in case of holes in the numbering
|
||
|
|
for iNd in self.nodes:
|
||
|
|
nd=self.p.node(iNd)
|
||
|
|
disp=self.p.node_displacement(iNd)
|
||
|
|
nd_xyz=[nd.x+disp[0], nd.y+disp[1], nd.z+disp[2]]
|
||
|
|
self.points.append(nd_xyz) # for pyvtk
|
||
|
|
fobj.write('%f %f %f \n'%
|
||
|
|
(nd.x+disp[0], nd.y+disp[1], nd.z+disp[2]))
|
||
|
|
self.nodes_dict[nd.id-1]=iNd
|
||
|
|
fobj.write('\n')
|
||
|
|
print('Nodes written to VTK: %i'%self.p.nodes())
|
||
|
|
|
||
|
|
def writeElements2VTK(self, fobj):
|
||
|
|
fobj.write('\nCELLS %i %i'%(self.p.elements(),self.p.elements()*9)+'\n')
|
||
|
|
self.cells=[] #for pyvtk
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
el=self.p.element(iEl)
|
||
|
|
cell_nodes=[] # for pyvtk
|
||
|
|
ndlist=el.items
|
||
|
|
#for k in [0, 1, 3, 2, 4, 5, 7, 6]: # FOR CELL TPYE VTK_VOXEL
|
||
|
|
#for k in [0, 4, 3, 1, 5, 7, 6, 2]:
|
||
|
|
for k in [0, 1, 2, 3, 4, 5, 6, 7]: # FOR CELL TYPE VTK_HEXAHEDRON
|
||
|
|
node=ndlist[k]-1
|
||
|
|
cell_nodes.append(self.nodes_dict[node])
|
||
|
|
self.cells.append(cell_nodes) # for pyvtk
|
||
|
|
for e in self.cells:
|
||
|
|
fobj.write('8 ')
|
||
|
|
for n in e:
|
||
|
|
fobj.write('%6i '%n)
|
||
|
|
fobj.write('\n')
|
||
|
|
fobj.write('\nCELL_TYPES %i'%self.p.elements()+'\n')
|
||
|
|
cnt=0
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
cnt=cnt+1
|
||
|
|
#fobj.write('11\n') #VTK_VOXEL
|
||
|
|
fobj.write('12 ') #VTK_HEXAHEDRON
|
||
|
|
if cnt>self.VTKcnt:
|
||
|
|
fobj.write('\n');cnt=0
|
||
|
|
fobj.write('\n')
|
||
|
|
print('Elements written to VTK: %i'%self.p.elements())
|
||
|
|
#print('Nr of nodes: ',self.nodes)
|
||
|
|
|
||
|
|
def writeElScalars2NodesVTK(self,fobj):
|
||
|
|
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
|
||
|
|
nScal=12
|
||
|
|
nComponents=1+nScal
|
||
|
|
fobj.write('SCALARS scalars float %i\n'%nComponents)
|
||
|
|
fobj.write('LOOKUP_TABLE default\n')
|
||
|
|
idxScal=self.nscal_list.index('Displacement Z')
|
||
|
|
for iNd in self.nodes:
|
||
|
|
#fobj.write('%f %f '%(self.p.node_scalar(iNd,idxScal), N.random.rand()))
|
||
|
|
fobj.write('%f '%(self.p.node_scalar(iNd,idxScal)))
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
el=self.p.element(iEl)
|
||
|
|
ndlist=el.items
|
||
|
|
if (iNd+1) in ndlist:
|
||
|
|
idx=ndlist.index(iNd+1)
|
||
|
|
for iV in range(0,nScal):
|
||
|
|
elData=self.p.element_scalar(iEl,35+iV)
|
||
|
|
fobj.write('%f '%(elData[idx].value))
|
||
|
|
break
|
||
|
|
fobj.write('\n')
|
||
|
|
fobj.write('\n')
|
||
|
|
|
||
|
|
def writeNodeScalars2VTK(self,fobj):
|
||
|
|
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
|
||
|
|
nNdDat=self.nscals
|
||
|
|
nComponents=1+nNdDat
|
||
|
|
self.pointDataScalars=[]
|
||
|
|
for idxNdScal in range(-3,self.nscals): #now include node x,y,z
|
||
|
|
if idxNdScal>=0:
|
||
|
|
datalabel=self.nscal_list[idxNdScal]
|
||
|
|
datalabel=re.sub("\s",'_',datalabel)
|
||
|
|
else:
|
||
|
|
if idxNdScal==-3: datalabel='node.x'
|
||
|
|
if idxNdScal==-2: datalabel='node.y'
|
||
|
|
if idxNdScal==-1: datalabel='node.z'
|
||
|
|
fobj.write('SCALARS %s float %i\n'%(datalabel,1))#nComponents))
|
||
|
|
fobj.write('LOOKUP_TABLE default\n')
|
||
|
|
cnt=0
|
||
|
|
for iNd in range(0,self.nnodes):
|
||
|
|
cnt=cnt+1
|
||
|
|
if idxNdScal>=0:
|
||
|
|
ndData=self.p.node_scalar(iNd,idxNdScal)
|
||
|
|
else:
|
||
|
|
nd=self.p.node(iNd)
|
||
|
|
if idxNdScal==-3: ndData=nd.x
|
||
|
|
if idxNdScal==-2: ndData=nd.y
|
||
|
|
if idxNdScal==-1: ndData=nd.z
|
||
|
|
fobj.write('%E '%(ndData))
|
||
|
|
if cnt>self.VTKcnt:
|
||
|
|
fobj.write('\n')
|
||
|
|
cnt=0
|
||
|
|
fobj.write('\n')
|
||
|
|
fobj.write('\n')
|
||
|
|
|
||
|
|
def writeElementData2VTK(self,fobj):
|
||
|
|
self.sig_vMises=[]
|
||
|
|
self.sig33=[]
|
||
|
|
idx_sig_vMises=self.getLabelNr('Equivalent Von Mises Stress')
|
||
|
|
idx_sig33=self.getLabelNr('Comp 33 of Cauchy Stress')
|
||
|
|
fobj.write('\nCELL_DATA %i\n'%self.p.elements())
|
||
|
|
nElDat=self.elscals
|
||
|
|
nComponents=1+nElDat
|
||
|
|
for idxElScal in range(0,self.elscals):
|
||
|
|
datalabel=self.elscal_list[idxElScal]
|
||
|
|
datalabel=re.sub("\s",'_',datalabel)
|
||
|
|
fobj.write('\n\nSCALARS %s float %i\n'%(datalabel,1))
|
||
|
|
fobj.write('LOOKUP_TABLE default\n')
|
||
|
|
cnt=0
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
cnt=cnt+1
|
||
|
|
elData=self.p.element_scalar(iEl,idxElScal)
|
||
|
|
avgScal=0.0
|
||
|
|
if datalabel in ['phi1', 'PHI','phi2']: # Euler angles should not be averaged
|
||
|
|
avgScal=avgScal+elData[0].value
|
||
|
|
else:
|
||
|
|
for IP in range(0,8):
|
||
|
|
avgScal=avgScal+elData[IP].value
|
||
|
|
avgScal=avgScal/8.
|
||
|
|
fobj.write('%E '%(avgScal))
|
||
|
|
if idxElScal==idx_sig_vMises:
|
||
|
|
self.sig_vMises.append(avgScal)
|
||
|
|
elif idxElScal==idx_sig33:
|
||
|
|
self.sig33.append(avgScal)
|
||
|
|
if cnt>self.VTKcnt:
|
||
|
|
fobj.write('\n')
|
||
|
|
cnt=0
|
||
|
|
fobj.write('\n')
|
||
|
|
|
||
|
|
def get_avg_el_scal(self,idxElScal):
|
||
|
|
result=[]
|
||
|
|
datalabel=self.elscal_list[idxElScal]
|
||
|
|
print('Collecting %s from all elements'%datalabel)
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
elData=self.p.element_scalar(iEl,idxElScal)
|
||
|
|
avgScal=0.0
|
||
|
|
for IP in range(0,8):
|
||
|
|
avgScal=avgScal+elData[IP].value
|
||
|
|
avgScal=avgScal/8.
|
||
|
|
result.append(avgScal)
|
||
|
|
return result
|
||
|
|
|
||
|
|
def writeUniaxiality2VTK(self,fobj):
|
||
|
|
#fobj.write('\nCELL_DATA %i\n'%self.p.elements())
|
||
|
|
datalabel='uniaxiality_sig_vMises_durch_sig33'
|
||
|
|
fobj.write('SCALARS %s float %i\n'%(datalabel,1))
|
||
|
|
fobj.write('LOOKUP_TABLE default\n')
|
||
|
|
cnt=0
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
cnt=cnt+1
|
||
|
|
#if abs(self.sig33[iEl])<1e-5:
|
||
|
|
if abs(self.sig_vMises[iEl])<1e-5:
|
||
|
|
datum=0.
|
||
|
|
else:
|
||
|
|
#datum=self.sig_vMises[iEl]/self.sig33[iEl]
|
||
|
|
datum=self.sig33[iEl]/self.sig_vMises[iEl]
|
||
|
|
fobj.write('%E '%(datum))
|
||
|
|
if cnt>self.VTKcnt:
|
||
|
|
fobj.write('\n')
|
||
|
|
cnt=0
|
||
|
|
fobj.write('\n')
|
||
|
|
|
||
|
|
def stress_per_element(self):
|
||
|
|
self.stress=[]
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
sig=self.avg_elten(2,elID=iEl)
|
||
|
|
self.stress.append(sig[0])
|
||
|
|
|
||
|
|
def mean_stress_per_element(self):
|
||
|
|
self.mean_stress=[]
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
sig=self.stress[iEl]
|
||
|
|
self.mean_stress.append(self.meanStress(sig))
|
||
|
|
|
||
|
|
def triaxiality_per_element(self):
|
||
|
|
# classical triaxiality
|
||
|
|
# 1/3 : uniax tension
|
||
|
|
self.triaxiality=[]
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
t=self.mean_stress[iEl]/self.sig_vMises[iEl]
|
||
|
|
self.triaxiality.append(t)
|
||
|
|
|
||
|
|
def moreElData2VTK(self,fobj,data=[],label='datalabel'):
|
||
|
|
fobj.write('SCALARS %s float %i\n'%(label,1))
|
||
|
|
fobj.write('LOOKUP_TABLE default\n')
|
||
|
|
cnt=0
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
cnt=cnt+1
|
||
|
|
fobj.write('%E '%(data[iEl]))
|
||
|
|
if cnt>self.VTKcnt:
|
||
|
|
fobj.write('\n')
|
||
|
|
cnt=0
|
||
|
|
fobj.write('\n')
|
||
|
|
|
||
|
|
def calc_lode_parameter(self):
|
||
|
|
# [-1 ... +1] see e.g. Wippler & Boehlke triaxiality measures doi:10.1002/pamm.201010061
|
||
|
|
# +1 : uniax tensile?
|
||
|
|
# 0 : shear
|
||
|
|
# -1 : uniax compr ?
|
||
|
|
self.lode=[]
|
||
|
|
try:
|
||
|
|
self.stress
|
||
|
|
except:
|
||
|
|
self.stress_per_element()
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
sig=self.stress[iEl]
|
||
|
|
lode=self.stress2lode(sig)
|
||
|
|
self.lode.append(lode)
|
||
|
|
|
||
|
|
def stress2lode(self,stress):
|
||
|
|
[pStress,pAxes]=self.princStress(stress)
|
||
|
|
s1=pStress[0]
|
||
|
|
s2=pStress[1]
|
||
|
|
s3=pStress[2]
|
||
|
|
lode=(2*s2-s1-s3) / ( s1 - s3 )
|
||
|
|
return lode
|
||
|
|
|
||
|
|
def princStress(self, stress):
|
||
|
|
"""
|
||
|
|
Function to compute 3D principal stresses and sort them.
|
||
|
|
from: http://geodynamics.org/svn/cig/short/3D/PyLith/trunk/playpen/postproc/vtkcff.py
|
||
|
|
"""
|
||
|
|
stressMat=N.array(stress)
|
||
|
|
(princStress, princAxes) = numpy.linalg.eigh(stressMat)
|
||
|
|
idx = princStress.argsort()
|
||
|
|
princStressOrdered = princStress[idx]
|
||
|
|
princAxesOrdered = princAxes[:,idx]
|
||
|
|
return princStressOrdered, princAxesOrdered
|
||
|
|
|
||
|
|
def avg_elten(self,
|
||
|
|
idxElTen, mat=0, elID=None):
|
||
|
|
tensum=N.zeros((3,3));
|
||
|
|
T=N.zeros((3,3));
|
||
|
|
pts=0;
|
||
|
|
avg=N.zeros((3,3));
|
||
|
|
#print 'Element Scalars'
|
||
|
|
#print self.p.element_scalar_label(elscal2)
|
||
|
|
if elID==None:
|
||
|
|
averaged_elements=range(0,self.nel)
|
||
|
|
else:
|
||
|
|
averaged_elements=[elID]
|
||
|
|
#for i in range (0,self.nel):
|
||
|
|
for i in averaged_elements:
|
||
|
|
if mat==0 or int(self.p.element_scalar(i,4)[0].value)==mat:
|
||
|
|
eldata=self.p.element(i)
|
||
|
|
T=self.p.element_tensor(i,idxElTen)
|
||
|
|
for k in range (0,8):
|
||
|
|
tensum[0][0] = tensum[0][0] + T[k].t11
|
||
|
|
tensum[0][1] = tensum[0][1] + T[k].t12
|
||
|
|
tensum[0][2] = tensum[0][2] + T[k].t13
|
||
|
|
#tensum1[1][0] = tensum1[1][0] + T1[k].t21
|
||
|
|
tensum[1][1] = tensum[1][1] + T[k].t22
|
||
|
|
tensum[1][2] = tensum[1][2] + T[k].t23
|
||
|
|
#tensum1[2][0] = tensum1[2][0] + T1[k].t31
|
||
|
|
#tensum1[2][1] = tensum1[2][1] + T1[k].t32
|
||
|
|
tensum[2][2] = tensum[2][2] + T[k].t33
|
||
|
|
pts=pts+1
|
||
|
|
avg=tensum/pts
|
||
|
|
#print avg
|
||
|
|
avg=self.fillComponents(avg)
|
||
|
|
#print avg
|
||
|
|
del [T]
|
||
|
|
return (avg,tensum,pts)
|
||
|
|
|
||
|
|
def fillComponents(self,
|
||
|
|
halftensor
|
||
|
|
):
|
||
|
|
halftensor[1][0]=halftensor[0][1]
|
||
|
|
halftensor[2][0]=halftensor[0][2]
|
||
|
|
halftensor[2][1]=halftensor[1][2]
|
||
|
|
return halftensor
|
||
|
|
|
||
|
|
def vMises(self,tensor33):
|
||
|
|
t=tensor33
|
||
|
|
s=(t[0,0]-t[1,1])**2+(t[1,1]-t[2,2])**2+(t[0,0]-t[2,2])**2+\
|
||
|
|
6*(t[0,1]**2+t[1,2]**2+t[2,0]**2)
|
||
|
|
vM=N.sqrt(s/2.)
|
||
|
|
return vM
|
||
|
|
|
||
|
|
def meanStress(self,tensor33):
|
||
|
|
t=tensor33
|
||
|
|
s=t[0,0]+t[1,1]+t[2,2]
|
||
|
|
ms=s/3.
|
||
|
|
return ms
|
||
|
|
|
||
|
|
def invariants(self,tensor33):
|
||
|
|
t=tensor33
|
||
|
|
I1=t[0,0]+t[1,1]+t[2,2]
|
||
|
|
I2=t[0,0]*t[1,1]+t[1,1]*t[2,2]+t[0,0]*t[2,2]-\
|
||
|
|
t[0,1]**2-t[1,2]**2-t[0,2]**2
|
||
|
|
# I3 = det(t)
|
||
|
|
I3=t[0,0]*t[1,1]*t[2,2]+\
|
||
|
|
2*t[0,1]*t[1,2]*t[2,0]-\
|
||
|
|
t[2,2]*t[0,1]**2-t[0,0]*t[1,2]**2-t[1,1]*t[0,2]**2
|
||
|
|
return [ I1, I2, I3 ]
|
||
|
|
|
||
|
|
|
||
|
|
class VTK_WRITER():
|
||
|
|
'''
|
||
|
|
The resulting vtk-file can be imported in Paraview 3.12
|
||
|
|
Then use Filters: Cell Data to Point Data + Contour
|
||
|
|
to plot semi-transparent iso-surfaces.
|
||
|
|
'''
|
||
|
|
import re
|
||
|
|
def __init__(self):
|
||
|
|
self.p=MARC_POST() # self.p
|
||
|
|
|
||
|
|
def openFile(self, filename='test.vtp'):
|
||
|
|
#if not self.f:#==None:
|
||
|
|
self.f=open(filename,'w+')
|
||
|
|
self.fname=filename
|
||
|
|
|
||
|
|
def writeFirstLines(self,
|
||
|
|
vtkFile=None,
|
||
|
|
version='2.0',
|
||
|
|
comment='Test',
|
||
|
|
dformat='ASCII', # BINARY | [ASCII]
|
||
|
|
dtype='UNSTRUCTURED_GRID' # UNSTRUCTURED GRID
|
||
|
|
):
|
||
|
|
if vtkFile==None:
|
||
|
|
vtkFile=self.f
|
||
|
|
# First Line contains Data format version
|
||
|
|
self.versionVTK=version
|
||
|
|
vtkFile.write('# vtk DataFile Version %s\n'%self.versionVTK)
|
||
|
|
# Comment goes to 2nd line and has maximum 256 chars
|
||
|
|
vtkFile.write(comment+'\n')
|
||
|
|
vtkFile.write(dformat+'\n')
|
||
|
|
vtkFile.write('DATASET '+dtype+'\n')
|
||
|
|
|
||
|
|
|
||
|
|
def marc2vtkBatch(self):
|
||
|
|
for iori in range(1,63):
|
||
|
|
#self.p=msc_post.MSC_POST()
|
||
|
|
self.p.postname='indent_fric0.3_R2.70_cA146.0_h0.320_ori%03i_OST_h19d.t16'%(iori)
|
||
|
|
if os.path.exists(self.p.postname):
|
||
|
|
self.marc2vtk(mode='fast', batchMode=1)
|
||
|
|
|
||
|
|
def marc2vtk(self, label=None, mode='fast',
|
||
|
|
batchMode=0,
|
||
|
|
incRange=None,
|
||
|
|
incStepMult=1.):
|
||
|
|
if batchMode==0:
|
||
|
|
try:
|
||
|
|
self.p
|
||
|
|
except:
|
||
|
|
self.p=MARC_POST()
|
||
|
|
### ---- CHANGE dir/job/model to process here
|
||
|
|
os.chdir('M:/nicu');
|
||
|
|
jobname='ori001'
|
||
|
|
self.p.postname='indent_fric0.3_R0.25_cA90.0_h0.010_4320els_'+jobname+'.t16'
|
||
|
|
### ----
|
||
|
|
|
||
|
|
self.p.opent16()
|
||
|
|
self.p.t16info()
|
||
|
|
|
||
|
|
incMax=self.p.p.increments();
|
||
|
|
if incRange is None:
|
||
|
|
incStep=5
|
||
|
|
incRange=range(0,incMax+1,incStep)
|
||
|
|
self.vtkPath=os.getcwd()+'/vtk_%s/'%self.p.postname
|
||
|
|
if not os.path.exists(self.vtkPath):
|
||
|
|
os.mkdir(self.vtkPath)
|
||
|
|
for inc in incRange:
|
||
|
|
print('Increment: %i'%inc)
|
||
|
|
self.p.p.moveto(inc)
|
||
|
|
t=self.p.p.time
|
||
|
|
sys.stdout.write('inc:%i, time:%.3f\n'%(self.p.p.increment,t))
|
||
|
|
self.incStr='inc%03i'%(inc*incStepMult)
|
||
|
|
self.openFile(filename=self.vtkPath+self.p.postname[0:-4]+'_'+
|
||
|
|
self.incStr+'.vtk')
|
||
|
|
self.writeFirstLines(comment=self.p.postname,
|
||
|
|
dtype='UNSTRUCTURED_GRID')
|
||
|
|
self.p.writeNodes2VTK(fobj=self.f)
|
||
|
|
self.p.writeElements2VTK(fobj=self.f)
|
||
|
|
self.p.writeNodeScalars2VTK(fobj=self.f)
|
||
|
|
self.p.writeElementData2VTK(fobj=self.f)
|
||
|
|
# insert generation and writing of derived post values
|
||
|
|
# *here*
|
||
|
|
|
||
|
|
self.f.close()
|
||
|
|
print('Increment (self.p.p.increment): %i'%self.p.p.increment)
|
||
|
|
print('Data written.')
|
||
|
|
print(self.p.postname+' ready.')
|
||
|
|
|
||
|
|
|
||
|
|
def scaleBar(self, length=1.0, posXYZ=[0., 0., 0.]):
|
||
|
|
self.fsb=open('micronbar_l%.1f.vtp'%length,'w+')
|
||
|
|
self.writeFirstLines(self.fsb, comment='micronbar')
|
||
|
|
pts=N.array([])
|
||
|
|
width=length*1.
|
||
|
|
height=length*1.
|
||
|
|
wVec=N.array([0., width, 0.])
|
||
|
|
lVec=N.array([length,0.,0.])
|
||
|
|
hVec=N.array([0.,0.,height])
|
||
|
|
posXYZ=posXYZ-0.5*wVec-0.5*lVec#-0.5*hVec # CENTERING Y/N
|
||
|
|
posXYZ=N.array(posXYZ)
|
||
|
|
pts=[posXYZ, posXYZ+lVec,
|
||
|
|
posXYZ+wVec,
|
||
|
|
posXYZ+wVec+lVec]
|
||
|
|
pts.extend([pts[0]+hVec,pts[1]+hVec,pts[2]+hVec,pts[3]+hVec])
|
||
|
|
print len(pts), pts
|
||
|
|
self.fsb.write('POINTS %i float\n'%len(pts))
|
||
|
|
for npts in range(0,len(pts)):
|
||
|
|
self.fsb.write('%f %f %f\n'%(pts[npts][0], pts[npts][1], pts[npts][2]))
|
||
|
|
if 1: #Triad
|
||
|
|
nCells=3
|
||
|
|
#nCells=1 #One Line
|
||
|
|
ptsPerCell=2 # Lines (Type=3)
|
||
|
|
#ptsPerCell=4 # Quads (Type=9)
|
||
|
|
#ptsPerCell=8 # Hexahedron (Type=12)
|
||
|
|
cellSize=(ptsPerCell+1)*nCells
|
||
|
|
self.fsb.write('CELLS %i %i\n'%(nCells,cellSize))
|
||
|
|
self.fsb.write('2 0 1\n') #X-Line
|
||
|
|
self.fsb.write('2 0 2\n') #Y-Line
|
||
|
|
self.fsb.write('2 0 4\n') #Z-Line
|
||
|
|
#self.fsb.write('4 0 1 3 2\n') #Quad
|
||
|
|
#self.fsb.write('%i 0 1 3 2 4 5 7 6\n'%ptsPerCell) #Hexahedron
|
||
|
|
self.fsb.write('CELL_TYPES %i\n'%(nCells))
|
||
|
|
self.fsb.write('3\n3\n3\n')#Line
|
||
|
|
#self.fsb.write('12\n')#Hexahedron
|
||
|
|
else: # Cube, change posXYZ
|
||
|
|
nCells=1
|
||
|
|
ptsPerCell=2 # Lines (Type=3)
|
||
|
|
#ptsPerCell=4 # Quads (Type=9)
|
||
|
|
#ptsPerCell=8 # Hexahedron (Type=12)
|
||
|
|
cellSize=(ptsPerCell+1)*nCells
|
||
|
|
self.fsb.write('CELLS %i %i\n'%(nCells,cellSize))
|
||
|
|
self.fsb.write('2 0 1\n') #Line
|
||
|
|
#self.fsb.write('4 0 1 3 2\n') #Quad
|
||
|
|
#self.fsb.write('%i 0 1 3 2 4 5 7 6\n'%ptsPerCell) #Hexahedron
|
||
|
|
self.fsb.write('CELL_TYPES %i\n'%(nCells))
|
||
|
|
self.fsb.write('3\n')#Line
|
||
|
|
#self.fsb.write('12\n')#Hexahedron
|
||
|
|
|
||
|
|
|
||
|
|
self.fsb.write('\n')
|
||
|
|
self.fsb.close()
|
||
|
|
print self.fsb
|
||
|
|
|
||
|
|
def example_unstructured(self):
|
||
|
|
self.openFile(filename='example_unstructured_grid.vtk')
|
||
|
|
#self.writeFirstLines()
|
||
|
|
self.f.write('''
|
||
|
|
# vtk DataFile Version 2.0
|
||
|
|
example_unstruct_grid
|
||
|
|
ASCII
|
||
|
|
|
||
|
|
POINTS 12 float
|
||
|
|
0 0 0
|
||
|
|
1 0 0
|
||
|
|
1 1 0
|
||
|
|
0 1 0
|
||
|
|
0 0 1
|
||
|
|
1 0 1
|
||
|
|
1 1 1
|
||
|
|
0 1 1
|
||
|
|
0 0 1.9
|
||
|
|
1 0 1.9
|
||
|
|
1 1 1.9
|
||
|
|
0 1 1.9
|
||
|
|
|
||
|
|
|
||
|
|
CELLS 2 18
|
||
|
|
8 0 1 2 3 4 5 6 7
|
||
|
|
8 4 5 6
|
||
|
|
7 8 9 10 11
|
||
|
|
|
||
|
|
CELL_TYPES 2
|
||
|
|
12
|
||
|
|
12
|
||
|
|
|
||
|
|
POINT_DATA 12
|
||
|
|
SCALARS nodex float 1
|
||
|
|
LOOKUP_TABLE default
|
||
|
|
2.34E+12
|
||
|
|
2.00
|
||
|
|
0.00
|
||
|
|
1.62
|
||
|
|
5.03
|
||
|
|
1.02
|
||
|
|
1.50
|
||
|
|
0.00
|
||
|
|
3 5 6 23423423423423423423.23423423''')
|
||
|
|
self.f.close()
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
def writeNodes2VTK(self, fobj):
|
||
|
|
self.VTKcnt=200 # how many numbers per line in vtk file
|
||
|
|
#self.VTKcnt=6
|
||
|
|
ndCnt=1
|
||
|
|
#self.nodes=range(0,10)
|
||
|
|
fobj.write('POINTS %i'%self.p.nodes()+' float\n')
|
||
|
|
for iNd in self.nodes:
|
||
|
|
nd=self.p.node(iNd)
|
||
|
|
disp=self.p.node_displacement(iNd)
|
||
|
|
#contact=self.p.node_scalar(iNd,contactNr)
|
||
|
|
#ndCnt=ndCnt+1
|
||
|
|
fobj.write('%f %f %f \n'%
|
||
|
|
#(nd.x, nd.y, nd.z))
|
||
|
|
(nd.x+disp[0], nd.y+disp[1], nd.z+disp[2]))
|
||
|
|
|
||
|
|
#if ndCnt>6:
|
||
|
|
# fobj.write('\n')
|
||
|
|
# ndCnt=1
|
||
|
|
fobj.write('\n')
|
||
|
|
print('Nodes written to VTK: %i'%self.p.nodes())
|
||
|
|
#print('Nr of nodes: ',self.nodes)
|
||
|
|
|
||
|
|
def writeElements2VTK(self, fobj):
|
||
|
|
fobj.write('\nCELLS %i %i'%(self.p.elements(),self.p.elements()*9)+'\n')
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
el=self.p.element(iEl)
|
||
|
|
#disp=self.p.node_displacement(iNd)
|
||
|
|
#contact=self.p.node_scalar(iNd,contactNr)
|
||
|
|
#ndCnt=ndCnt+1
|
||
|
|
fobj.write('8 ')
|
||
|
|
ndlist=el.items
|
||
|
|
#for k in [0, 1, 3, 2, 4, 5, 7, 6]: # FOR CELL TPYE VTK_VOXEL
|
||
|
|
#for k in [0, 4, 3, 1, 5, 7, 6, 2]:
|
||
|
|
for k in [0, 1, 2, 3, 4, 5, 6, 7]: # FOR CELL TYPE VTK_HEXAHEDRON
|
||
|
|
fobj.write('%6i '%(ndlist[k]-1))
|
||
|
|
fobj.write('\n')
|
||
|
|
#if ndCnt>6:
|
||
|
|
# fobj.write('\n')
|
||
|
|
# ndCnt=1
|
||
|
|
fobj.write('\nCELL_TYPES %i'%self.p.elements()+'\n')
|
||
|
|
cnt=0
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
cnt=cnt+1
|
||
|
|
#fobj.write('11\n') #VTK_VOXEL
|
||
|
|
fobj.write('12 ') #VTK_HEXAHEDRON
|
||
|
|
if cnt>self.VTKcnt:
|
||
|
|
fobj.write('\n');cnt=0
|
||
|
|
fobj.write('\n')
|
||
|
|
print('Elements written to VTK: %i'%self.p.elements())
|
||
|
|
#print('Nr of nodes: ',self.nodes)
|
||
|
|
|
||
|
|
def writeElScalars2NodesVTK(self,fobj):
|
||
|
|
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
|
||
|
|
nScal=12
|
||
|
|
nComponents=1+nScal
|
||
|
|
fobj.write('SCALARS scalars float %i\n'%nComponents)
|
||
|
|
fobj.write('LOOKUP_TABLE default\n')
|
||
|
|
idxScal=self.nscal_list.index('Displacement Z')
|
||
|
|
for iNd in self.nodes:
|
||
|
|
fobj.write('%f '%(self.p.node_scalar(iNd,idxScal)))
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
el=self.p.element(iEl)
|
||
|
|
ndlist=el.items
|
||
|
|
if (iNd+1) in ndlist:
|
||
|
|
idx=ndlist.index(iNd+1)
|
||
|
|
for iV in range(0,nScal):
|
||
|
|
elData=self.p.element_scalar(iEl,35+iV)
|
||
|
|
fobj.write('%f '%(elData[idx].value))
|
||
|
|
break
|
||
|
|
fobj.write('\n')
|
||
|
|
fobj.write('\n')
|
||
|
|
|
||
|
|
def writeNodeScalars2VTK(self,fobj):
|
||
|
|
#print('writeElementData2VTK')
|
||
|
|
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
|
||
|
|
nNdDat=self.nscals
|
||
|
|
nComponents=1+nNdDat
|
||
|
|
for idxNdScal in range(-3,self.nscals): # include node x,y,z
|
||
|
|
if idxNdScal>=0:
|
||
|
|
datalabel=self.nscal_list[idxNdScal]
|
||
|
|
datalabel=re.sub("\s",'_',datalabel)
|
||
|
|
else:
|
||
|
|
if idxNdScal==-3: datalabel='node.x'
|
||
|
|
if idxNdScal==-2: datalabel='node.y'
|
||
|
|
if idxNdScal==-1: datalabel='node.z'
|
||
|
|
fobj.write('SCALARS %s float %i\n'%(datalabel,1))#nComponents))
|
||
|
|
fobj.write('LOOKUP_TABLE default\n')
|
||
|
|
cnt=0
|
||
|
|
for iNd in range(0,self.nnodes):
|
||
|
|
cnt=cnt+1
|
||
|
|
if idxNdScal>=0:
|
||
|
|
ndData=self.p.node_scalar(iNd,idxNdScal)
|
||
|
|
else:
|
||
|
|
nd=self.p.node(iNd)
|
||
|
|
if idxNdScal==-3: ndData=nd.x
|
||
|
|
if idxNdScal==-2: ndData=nd.y
|
||
|
|
if idxNdScal==-1: ndData=nd.z
|
||
|
|
fobj.write('%E '%(ndData))
|
||
|
|
if cnt>self.VTKcnt:
|
||
|
|
fobj.write('\n')
|
||
|
|
cnt=0
|
||
|
|
fobj.write('\n')
|
||
|
|
fobj.write('\n')
|
||
|
|
|
||
|
|
def writeElementData2VTK(self,fobj):
|
||
|
|
#print('writeElementData2VTK')
|
||
|
|
fobj.write('\nCELL_DATA %i\n'%self.p.elements())
|
||
|
|
nElDat=self.elscals
|
||
|
|
nComponents=1+nElDat
|
||
|
|
for idxElScal in range(0,self.elscals):
|
||
|
|
datalabel=self.elscal_list[idxElScal]
|
||
|
|
datalabel=re.sub("\s",'_',datalabel)
|
||
|
|
fobj.write('\n\nSCALARS %s float %i\n'%(datalabel,1))#nComponents))
|
||
|
|
fobj.write('LOOKUP_TABLE default\n')
|
||
|
|
cnt=0
|
||
|
|
for iEl in range(0,self.nel):
|
||
|
|
cnt=cnt+1
|
||
|
|
elData=self.p.element_scalar(iEl,idxElScal)
|
||
|
|
avgScal=0.0
|
||
|
|
if datalabel in ['phi1', 'PHI','phi2']:
|
||
|
|
avgScal=avgScal+elData[0].value
|
||
|
|
else:
|
||
|
|
for IP in range(0,8):
|
||
|
|
avgScal=avgScal+elData[IP].value
|
||
|
|
avgScal=avgScal/8.
|
||
|
|
fobj.write('%E '%(avgScal))
|
||
|
|
if cnt>self.VTKcnt:
|
||
|
|
fobj.write('\n')
|
||
|
|
cnt=0
|
||
|
|
fobj.write('\n')
|
||
|
|
|
||
|
|
|
||
|
|
def example1(self):
|
||
|
|
self.openFile()
|
||
|
|
self.writeFirstLines()
|
||
|
|
self.f.write('''DATASET POLYDATA
|
||
|
|
POINTS 8 float
|
||
|
|
0.0 0.0 0.0
|
||
|
|
1.0 0.0 0.0
|
||
|
|
1.0 1.0 0.0
|
||
|
|
0.0 1.0 0.0
|
||
|
|
0.0 0.0 1.0
|
||
|
|
1.0 0.0 1.0
|
||
|
|
1.0 1.0 1.0
|
||
|
|
0.0 1.0 1.0
|
||
|
|
POLYGONS 6 30
|
||
|
|
4 0 1 2 3
|
||
|
|
4 4 5 6 7
|
||
|
|
4 0 1 5 4
|
||
|
|
4 2 3 7 6
|
||
|
|
4 0 4 7 3
|
||
|
|
4 1 2 6 5
|
||
|
|
|
||
|
|
CELL_DATA 6
|
||
|
|
SCALARS cell_scalars int 1
|
||
|
|
LOOKUP_TABLE default
|
||
|
|
0
|
||
|
|
1
|
||
|
|
2
|
||
|
|
3
|
||
|
|
4
|
||
|
|
5
|
||
|
|
NORMALS cell_normals float
|
||
|
|
0 0 -1
|
||
|
|
0 0 1
|
||
|
|
0 -1 0
|
||
|
|
0 1 0
|
||
|
|
-1 0 0
|
||
|
|
1 0 0
|
||
|
|
FIELD FieldData 2
|
||
|
|
cellIds 1 6 int
|
||
|
|
0 1 2 3 4 5
|
||
|
|
faceAttributes 2 6 float
|
||
|
|
0.0 1.0 1.0 2.0 2.0 3.0 3.0 4.0 4.0 5.0 5.0 6.0
|
||
|
|
|
||
|
|
POINT_DATA 8
|
||
|
|
SCALARS sample_scalars float 1
|
||
|
|
LOOKUP_TABLE my_table
|
||
|
|
0.0
|
||
|
|
1.0
|
||
|
|
2.0
|
||
|
|
3.0
|
||
|
|
4.0
|
||
|
|
5.0
|
||
|
|
6.0
|
||
|
|
7.0
|
||
|
|
LOOKUP_TABLE my_table 8
|
||
|
|
0.0 0.0 0.0 1.0
|
||
|
|
1.0 0.0 0.0 1.0
|
||
|
|
0.0 1.0 0.0 1.0
|
||
|
|
1.0 1.0 0.0 1.0
|
||
|
|
0.0 0.0 1.0 1.0
|
||
|
|
1.0 0.0 1.0 1.0
|
||
|
|
0.0 1.0 1.0 1.0
|
||
|
|
1.0 1.0 1.0 1.0''')
|
||
|
|
self.f.close()
|
||
|
|
|
||
|
|
|
||
|
|
import pyvtk
|
||
|
|
class marc_to_vtk():
|
||
|
|
'''
|
||
|
|
Anybody wants to implement it with pyvtk?
|
||
|
|
The advantage would be that pyvtk can also wirte the
|
||
|
|
<xml>-VTK format and binary.
|
||
|
|
These can be plotted with mayavi.
|
||
|
|
'''
|
||
|
|
def __init__(self):
|
||
|
|
self.p=[]#MARC_POST() # self.p
|
||
|
|
|
||
|
|
def run(self):
|
||
|
|
vtk = pyvtk.VtkData(\
|
||
|
|
pyvtk.UnstructuredGrid(self.p.points,
|
||
|
|
hexahedron=self.p.cells),
|
||
|
|
'm2v output')
|
||
|
|
vtk.tofile('m2v_file')
|
||
|
|
#vtk.tofile('example3b','binary')
|
||
|
|
#VtkData('example3')
|