# 以xml文本的方式读取vtu文件

Please refresh the page if equations are not rendered correctly.
---------------------------------------------------------------

## 1. 指定单元数和尺寸创建体素网格

import numpy as np
import vtk
import xml.dom.minidom as minidom
from xml.dom import minidom

import timeit

'''

'''
##################################
###  generate points & cells   ###
##################################

starttime1 = timeit.default_timer()

dx, dy, dz = 1, 1, 1

# number of cells in different directions
nrow = 100
ncolumn = 150
ndepth = 60

lx = dx*nrow
ly = dy*ncolumn
lz = dz * ndepth

ncells = nrow * ncolumn * ndepth
npoints = (nrow + 1) * (ncolumn + 1) * (ndepth + 1)

# Coordinates
X = np.arange(0, lx + 0.1*dx, dx, dtype='float64')
Y = np.arange(0, ly + 0.1*dy, dy, dtype='float64')
Z = np.arange(0, lz + 0.1*dz, dz, dtype='float64')

x = np.zeros((nrow + 1, ncolumn + 1, ndepth + 1))
y = np.zeros((nrow + 1, ncolumn + 1, ndepth + 1))
z = np.zeros((nrow + 1, ncolumn + 1, ndepth + 1))

for k in range(ndepth + 1):
for j in range(ncolumn + 1):
for i in range(nrow + 1):
x[i,j,k] = X[i]
y[i,j,k] = Y[j]
z[i,j,k] = Z[k]
# points
x = x.reshape([1,-1])[0]
y = y.reshape([1,-1])[0]
z = z.reshape([1,-1])[0]

points = np.zeros([x.size,4])
points[:,0] = x # points_organized[:,:,:,0]
points[:,1]=y # points_organized[:,:,:,1]
points[:,2]=z # points_organized[:,:,:,2]
points[:,3]=np.arange(x.size) # index of each point, points_organized[:,:,:,3]

##print('-------Points list--------')
##print(points)

# connectivity: Speed up cell allocation
points_organized = points.reshape([nrow+1,ncolumn+1,ndepth+1,4])
con = points_organized[:,:,:,-1]

##print('-------points_organized list--------')
##print(points_organized)

endtime1 = timeit.default_timer()
print("The time for data generation", endtime1 - starttime1)

starttime2 = timeit.default_timer()

connectivity = np.empty(shape=(ncells,8),dtype='uint32')
index_con = 0
for k in range(nrow):
for j in range(ncolumn):
for i in range(ndepth):
arr1 = np.ravel(con[k,:,:][j:j+2,i:i+2], 'F')
arr2 = np.ravel(con[k+1,:,:][j:j+2,i:i+2], 'F')
vertex = np.array(list(zip(arr1,arr2)),dtype='uint32').flatten()
connectivity[index_con,:] = vertex
index_con += 1

if index_con % 6000 == 0:
print(round(index_con/npoints*100,2),'%')

endtime2 = timeit.default_timer()
print("The time for for-loops", endtime2 - starttime2)

print('-------Vertex list--------')
print(connectivity)

starttime3 = timeit.default_timer()

# offsets: only for element type 11.
offsets = np.arange(8,ncells*8+0.5,8, dtype='uint32')

# cellTypes types

cellType = vtk.vtkVoxel().GetCellType() # # get cell type for voxel
cellTypes = np.zeros((nrow, ncolumn, ndepth), dtype='uint32').flatten() + cellType

# dataArray

pressure = np.random.rand(ncells).reshape( (nrow, ncolumn, ndepth))
temp = np.random.rand(npoints).reshape( (nrow + 1, ncolumn + 1, ndepth + 1))

##################################
### create xml-based vtu file  ###
##################################

PolyDMT = minidom.Document()

root = PolyDMT.createElement('VTKFile')
root.setAttribute('type', 'UnstructuredGrid')  # 属性是加在标签名之后、尖括号<>之中的
root.setAttribute('version', '0.1')
root.setAttribute('byte_order', 'LittleEndian')
PolyDMT.appendChild(root)

GridTypeChild = PolyDMT.createElement('UnstructuredGrid')
root.appendChild(GridTypeChild)

piece = PolyDMT.createElement('Piece')
piece.setAttribute('NumberOfPoints', str(npoints))
piece.setAttribute('NumberOfCells', str(ncells))
GridTypeChild.appendChild(piece)

Points = PolyDMT.createElement('Points')
piece.appendChild(Points)

cells = PolyDMT.createElement('Cells')
piece.appendChild(cells)

pointData = PolyDMT.createElement('PointData')
piece.appendChild(pointData)

cellData = PolyDMT.createElement('CellData')
piece.appendChild(cellData)

a = np.arange(0, 10)

def np2str(np_arr):
np_arr = np_arr.flatten()
str_arr = ""
for i in np_arr:
str_arr += str(i) + " "
return str_arr

def dataFun(type, noComponents, dataContent, appendNode, name='', format='ascii'):
'''
创建元素、属性及其dataArray
'''
dataArray = PolyDMT.createElement('DataArray')
if appendNode.tagName != "Points": dataArray.setAttribute('Name', name)
dataArray.setAttribute('type', type)
if appendNode.tagName != "Cells": dataArray.setAttribute('NumberOfComponents', noComponents)
dataArray.setAttribute('format', format)
appendNode.appendChild(dataArray)
dataArray.appendChild(PolyDMT.createTextNode(np2str(dataContent)))

dataFun('Float32', '3', points[:,:3][:npoints], Points)
dataFun('Int32', '1', connectivity[:ncells*8], cells, 'connectivity')
dataFun('Int32', '1', offsets[:ncells], cells, 'offsets')
dataFun('UInt8', '1', cellTypes[:ncells], cells, 'types')

##dataFun('Float32', '1', greylevel, cellData, 'greylevel')
##dataFun('Int32', '3', a, cellData, 'YarnTangent')
##dataFun('Int32', '2', a, cellData, 'Location')
##dataFun('Int32', '1', a, cellData, 'VolumeFraction')
##dataFun('Int32', '1', a, cellData, 'SurfaceDistance')
dataFun('Float32', '1', temp, pointData, 'Temperature')

# write out the xml file.
xml_str = PolyDMT.toprettyxml(indent="\t")
save_path_file = "./createPolyDMT30.vtu"
with open(save_path_file, "w", buffering = 1) as f:
f.write(xml_str)
f.close()
endtime3 = timeit.default_timer()

print("The time for file output", endtime3 - starttime3)


## 2. xml文件树的创建

import xml.dom.minidom as minidom
import numpy as np
import xml.sax

# class VtkHandler(xml.sax.ContentHandler):
#     def startElement(self, name, attrs):
#         print(name, attrs)
# handler = VtkHandler()
# parser = xml.sax.make_parser()
# parser.setContentHandler(handler)
# parser.parse('./texgen.vtu')

# domtree = xml.dom.minidom.parse('./texgen.vtu')
#
# group= domtree.documentElement
#
# for i in domtree.childNodes:
#     print(i.tagName)
# VTKFile = group.getElementsByTagName("VTKFile")
# print(VTKFile)
#
# DataArrays = group.getElementsByTagName('DataArray')
# # print("Name:{} ".format(group.getElementsByTagName('DataArray')[10].childNodes[0].data)) # 输出节点的数值
# # print("%d DataArrays" % DataArrays.length)
# for DataArray in DataArrays:
#     if DataArray.hasAttribute('Name'):
#         print('---'*200)
#         print(DataArray.getAttribute('Name'))
#         print(DataArray.childNodes[0].data)

##
### 修改原有标签的值
##persons[0].setAttribute("id","100")
##persons[1].getElementsByTagName('name')[0].childNodes[0].nodeValue = "New Name"
##persons[2].getElementsByTagName('age')[0].childNodes[0].nodeValue = "-10"
##

from xml.dom import minidom
import os

PolyDMT = minidom.Document()

root = PolyDMT.createElement('VTKFile')
root.setAttribute('type', 'UnstructuredGrid')  # 属性是加在标签名之后、尖括号<>之中的
root.setAttribute('version', '0.1')
root.setAttribute('byte_order', 'LittleEndian')
PolyDMT.appendChild(root)

GridTypeChild = PolyDMT.createElement('UnstructuredGrid')
root.appendChild(GridTypeChild)

piece = PolyDMT.createElement('Piece')
piece.setAttribute('NumberOfPoints', '216')
piece.setAttribute('NumberOfCells', '125')
GridTypeChild.appendChild(piece)

points = PolyDMT.createElement('Points')
piece.appendChild(points)

cells = PolyDMT.createElement('Cells')
piece.appendChild(cells)

pointData = PolyDMT.createElement('PointData')
piece.appendChild(pointData)

cellData = PolyDMT.createElement('CellData')
piece.appendChild(cellData)

a = np.arange(0, 10)

def np2str(np_arr):
np_arr = np_arr.reshape(1, -1)[0]
str_arr = ""
for i in a:
str_arr += str(i) + " "
return str_arr

print(points.tagName=='Points')

def dataFun(type, noComponents, dataContent, appendNode, name='', format='ascii'):
dataArray = PolyDMT.createElement('DataArray')
if appendNode.tagName != "Points": dataArray.setAttribute('Name', name)
dataArray.setAttribute('type', type)
if appendNode.tagName != "Cells": dataArray.setAttribute('NumberOfComponents', noComponents)
dataArray.setAttribute('format', format)
appendNode.appendChild(dataArray)
dataArray.appendChild(PolyDMT.createTextNode(np2str(dataContent)))

dataFun('Float32', '3', a, points)
dataFun('Int32', '1', a, cells, 'connectivity')
dataFun('Int32', '1', a, cells, 'offsets')
dataFun('Int32', '1', a, cells, 'types')
dataFun('Int32', '1', a,  cellData, 'YarnIndex')
dataFun('Int32', '3', a,  cellData, 'YarnTangent')
dataFun('Int32', '2', a,  cellData, 'Location')
dataFun('Int32', '1', a,  cellData, 'VolumeFraction')
dataFun('Int32', '1', a,  cellData, 'SurfaceDistance')
dataFun('Int32', '3', a,  cellData, 'Orientation')

# write out the xml file.
xml_str = PolyDMT.toprettyxml(indent="\t")
save_path_file = "data1.vtu"
with open(save_path_file, "w") as f:
f.write(xml_str)


## 3. 将TexGen生成的体素网格转换为图片序列

import pyvista as pv
import polykriging as pk
mesh_shape = (20, 20, 5) # number of cells in x, y and z direction
# convert the mesh to images
pk.io.voxel2img(mesh, mesh_shape, dataset="YarnIndex", save_path="./img/",
scale=50, img_name="img", format="tif",
scale_algrithm="linear")

import xml.dom.minidom
import numpy as np
from skimage import io
from mayavi import mlab

# 将TexGen生成的模型转换为图片
# 注意像素尺寸

nx, ny, nz = 50, 50, 20  # number of cells in each direction

# load the unstructed grid (voxel model from TexGen with extension of .vtu)
domtree = xml.dom.minidom.parse('C:/Users/palme/Desktop/111.vtu')
rootNode = domtree.documentElement
UnstructuredGrid = rootNode.getElementsByTagName('UnstructuredGrid')

Piece = UnstructuredGrid[0].getElementsByTagName('Piece')
childlist = Piece[0].childNodes

CellData = Piece[0].getElementsByTagName('CellData')
DataArray = CellData[0].getElementsByTagName('DataArray')

###################################################
##io.use_plugin('pil')
##        'D:/04_coding/Python/06_Visualization/Mesh/00_data_acquisition/3D_2_4_5layers/Long_0/Long_0_img/*.tif') #
####io.imshow_collection(ic)
####io.show()
##image_3d = io.concatenate_images(ic) #

#3mlab.contour3d(image_3d[:,70:130,45:130])
#mlab.outline()
#np.save('./greylevel.npy',image_3d[:,70:130,45:130])

def np2str(np_arr):
np_arr = np_arr.reshape(1,-1)[0]
str_arr = ""
for i in np_arr:
str_arr += str(i)+" "
return str_arr

#交换数组的轴至：shape(z,y,x)
image_3d = np.swapaxes(image_3d, 0,1)
image_3d = np.swapaxes(image_3d, 1,2)

imgId= domtree.createElement('DataArray')
# 切片时一定要注意多余数据不影响显示，故ParaView不会报错
imgId.appendChild(domtree.createTextNode(np2str(image_3d[:, :, 0:704:4])))
imgId.setAttribute('Name','image_3d')
imgId.setAttribute('NumberOfComponents','1')
imgId.setAttribute('format','ascii')
imgId.setAttribute('type','Float64')
CellData[0].appendChild(imgId)

# update DataArray
DataArray = CellData[0].getElementsByTagName('DataArray')
# Save vtu as image sequence
for item in DataArray:
if item.getAttribute('Name') == "YarnIndex":
YarnIndex = item.childNodes[0].data
YarnIndex=YarnIndex.split(" ")
img_sequence = np.empty([nz, nx, ny])
for i in range(nz):
start = i*nx*ny
end = (i+1)*nx*ny
img = np.array(YarnIndex[start:end], dtype = float).reshape([nx,ny])
img_sequence[i,:,:]  = img
io.imsave("C:/Users/palme/Desktop/img/" + str(i)+'.jpeg', img)

# Delete unnecessary cell data
for item in DataArray:
if item.getAttribute('Name') in ["SurfaceDistance", "Location", "YarnTangent", "Orientation"]:
CellData[0].removeChild(item)

# write out the xml file.
domtree.writexml(open("C:/Users/palme/Desktop/result.vtu","w"))
### Write the result to a vtu file
##with open("C:/Users/palme/Desktop/result.vtu","w") as fs:
##    fs.write(domtree.toxml())
##    fs.close()


## 4. 相关文件：

Everything not saved will be lost.