DAMASK_EICMD/lib/damask/result/marc2vtk.py

812 lines
25 KiB
Python
Raw Normal View History

# $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')