【Python】將網格資料寫入到VTK檔案

FE-有限元鹰發表於2024-08-29

1. vtk 檔案格式

根據官網進行總結

vtk檔案組成:5個部分.

第一部分,第一行:表明檔案版本.寫"# vtk DataFile Version 2.0"就行

第二部分,第二行:表明標題(title).隨便寫.

第三部分,第三行:ASCII或者 BINARY

第四部分,開始定義dataset structure.這部分用於描述資料集的幾何和拓撲結構.以關鍵字"DATASET type"的一行開始.後續部分根據不同的dataset type,使用其他關鍵字/資料組合定義實際資料.

支援以下type:
STRUCTURED_POINTS
STRUCTURED_GRID
UNSTRUCTURED_GRID
POLYDATA
STRUCTURED_POINTS
RECTILINEAR_GRID
FIELD
關鍵字不區分大小寫,可以用空格分隔.

第五部分,定義點資料或者實體資料等.稱為:dataset attributes 如:
POINT_DATA type
...
CELL_DATA type
...

注意:Indices are 0-offset. Thus the first point is point id 0.;Cell types and indices are of type int; The geometry/topology description must occur prior to the data attribute description.

1.1 Dataset Format _ 定義幾何/拓撲資料和關係

VTK支援5種資料集格式:structured points, structured grid, rectilinear grid, unstructured grid, and polygonal data.

這裡只描述非結構化網格,其他不說,用不到.

1.2 Unstructured Grid

非結構化網格資料集由任何可用的單元型別的任意組合組成(比如四六面體混合).非結構化網格由點、單元格和單元格型別定義.

CELLS keyword 需要兩個引數:

  • n : the number of cells
  • size : the size of the cell list. The cell list size is CELLS關鍵字下所有整數的個數 (i.e., sum of numPoints and connectivity indices over each cell).

例如:60代表的是從第一行以下所有整數的個數

CELLS 11 60
8 0 1 4 3 6 7 10 9
8 1 2 4 5 7 8 10 11
4 6 10 9 12
4 11 14 10 13
6 15 16 17 14 13 12
6 18 15 19 16 20 17
4 22 23 20 19
3 21 22 18
3 22 19 18
2 26 25
1 24

CELL_TYPES keyword需要一個引數: the number of cells n.這個引數值對應於CELLS的n引數.這個關鍵字的資料行,是指定每個單元格的單元型別的單個整數值

案例:定義非結構化網格的拓撲

DATASET UNSTRUCTURED_GRID
POINTS n dataType
p0x p0y p0z
p1x p1y p1z
…
p(n-1)x p(n-1)y p(n-1)z

CELLS n size
numPoints0, i0, j0, k0, …
numPoints1, i1, j1, k1, …
numPoints2, i2, j2, k2, …
…
numPointsn-1, in-1, jn-1, kn-1, …

CELL_TYPES n
type0
type1
type2
…
typen-1

1.3 Dataset Attribute Format --- 定義附著在節點或者單元上的資料

支援以下資料集屬性:

  • 標量(一到四個分量)
  • 向量
  • 法線
  • 紋理座標(1D, 2D和3D)
  • 張量
  • field data

此外,可以使用RGBA顏色規範定義lookup table,與標量資料關聯.

定義dataset Attribute 要給定一個獨特的dataName,不能有空格.一個檔案中可以包含多個相同型別的dataset Attribute. For example, two different scalar fields defined on the dataset points, pressure and temperature, can be contained in the same file.

1.3.1 Scalars.

SCALARS dataName dataType numComp
LOOKUP_TABLE tableName
s0
s1
…
sn-1

標量定義包括lookup table的指定.lookup table 是可選的,如果不指定,vtk會使用預設表(tableName指定為default);numComp選項也是可選的,如果不給定,預設為1.引數範圍:1~4

color scalars( 直接對映到顏色的無符號字元值 )的定義varies depending upon the number of values (nValues) per scalar.** If the file format is ASCII, the color scalars are defined using nValues float values between (0,1). **

COLOR_SCALARS dataName nValues
c00 c01 … c0(nValues-1)
c10 c11 … c1(nValues-1)
…
c(n-1)0 c(n-1)1 … c(n-1)(nValues-1)

1.3.2 Lookup Table.

LOOKUP_TABLE tableName size
r0 g0 b0 a0
r1 g1 b1 a1
…
rsize-1 gsize-1 bsize-1 asize-1

tableName 是表的標識,不能有空格.每一行表示rgba陣列.a表示透明度,a=0時顏色透明.如果檔案是ASCII格式,每個值範圍為(0-1的浮點數;

1.3.3 Vectors.

VECTORS dataName dataType
v0x v0y v0z
v1x v1y v1z
…
v(n-1)x v(n-1)y v(n-1)z

1.3.4 Normals.

Normals are assumed normalized |n| = 1.

NORMALS dataName dataType
n0x n0y n0z
n1x n1y n1z
…
n(n-1)x n(n-1)y n(n-1)z

1.3.5 Tensors.

Currently only real-valued, symmetric tensors are supported.

TENSORS dataName dataType
t000 t001 t002
t010 t011 t012
t020 t021 t022

t100 t101 t102
t110 t111 t112
t120 t121 t122

tn - 100 tn - 101 tn - 102
tn - 110 tn - 111 tn - 112
tn - 120 tn - 121 tn - 122

1.4 example

# vtk DataFile Version 2.0
Unstructured Grid Example
ASCII
DATASET UNSTRUCTURED_GRID

POINTS 27 float
0 0 0  1 0 0  2 0 0  0 1 0  1 1 0  2 1 0
0 0 1  1 0 1  2 0 1  0 1 1  1 1 1  2 1 1
0 1 2  1 1 2  2 1 2  0 1 3  1 1 3  2 1 3
0 1 4  1 1 4  2 1 4  0 1 5  1 1 5  2 1 5
0 1 6  1 1 6  2 1 6

CELLS 11 60
8 0 1 4 3 6 7 10 9
8 1 2 4 5 7 8 10 11
4 6 10 9 12
4 11 14 10 13
6 15 16 17 14 13 12
6 18 15 19 16 20 17
4 22 23 20 19
3 21 22 18
3 22 19 18
2 26 25
1 24

CELL_TYPES 11
12
11
10
8
7
6
9
5
4
3
1

POINT_DATA 27
SCALARS scalars float 1
LOOKUP_TABLE default
0.0 1.0 2.0 3.0 4.0 5.0
6.0 7.0 8.0 9.0 10.0 11.0
12.0 13.0 14.0 15.0 16.0 17.0
18.0 19.0 20.0 21.0 22.0 23.0
24.0 25.0 26.0

VECTORS vectors float
1 0 0  1 1 0  0 2 0  1 0 0  1 1 0  0 2 0
1 0 0  1 1 0  0 2 0  1 0 0  1 1 0  0 2 0
0 0 1  0 0 1  0 0 1  0 0 1  0 0 1  0 0 1
0 0 1  0 0 1  0 0 1  0 0 1  0 0 1  0 0 1
0 0 1  0 0 1  0 0 1

CELL_DATA 11
SCALARS scalars float 1
LOOKUP_TABLE CellColors
0.0 1.0 2.0 3.0 4.0 5.0
6.0 7.0 8.0 9.0 10.0

LOOKUP_TABLE CellColors 11
.4 .4 1 1
.4 1 .4 1
.4 1 1 1
1 .4 .4 1
1 .4 1 1
1 1 .4 1
1 1 1 1
1 .5 .5 1
.5 1 .5 1
.5 .5 .5 1
1 .5 .4 1

2.使用python將網格和資料寫入vtk檔案

def export_vtk(output_file:str,
                nodes:np.ndarray,
                element_Blocks:List[np.ndarray],
                element_types:List[str],
                point_data: Dict[str:Dict[str:np.ndarray]] = {},
                cell_data:  Dict[str:Dict[str:np.ndarray]] = {},
                **kwargs)->bool:
    """匯出vtk格式的網格資料檔案
    inputs:
    output_file: 輸出檔名
    nodes: 節點座標矩陣,shape=(n,3)
    element_Blocks: 單元列表,
    element_types: 單元型別列表,
    point_data: 節點資料字典.e.g. {'SCALARS':{'disp':disp_array,'stress':stress_array},'VECTORS':{'velocity':velocity_array}}
    cell_data: 單後設資料字典.e.g. {'SCALARS':{'disp':disp_array,'stress':stress_array},'VECTORS':{'velocity':velocity_array}}
    """
    def cal_cell_num(element_Blocks:List[np.ndarray])->tuple:
        cell_num,cell_piont_num=0,0
        for _,block in enumerate(element_Blocks):
            cell_num+=block.shape[0]
            cell_piont_num+=block.shape[0]*(block.shape[1]+1)
        return cell_num,cell_piont_num
    
    SEP="  " if 'sep'not in kwargs.keys()  else kwargs['sep']
    
    type_dict={"beam": 3, # vtk_line, 梁單元
                "tri": 5, # vtk_triangle, 三角形單元
                "quad":9, # vtk_quad, 四邊形單元
                "tet":10, # vtk_tetra, 四面體單元
                "hex":12, # vtk_hexahedron, 六面體單元
                "wedge":13, # vtk_wedge, 三稜柱單元
                "pyramid":14, # vtk_pyramid, 四稜錐單元
                }
    
    f=open(output_file,'w')
    
    # 寫入檔案頭
    f.write('# vtk DataFile Version 3.0\n')
    f.write(f'{kwargs["title"]}\n') if 'title' in kwargs.keys() else f.write('Finite Element Results\n')
    f.write('ASCII\n')
    f.write('DATASET UNSTRUCTURED_GRID\n')
    
    # 寫入點座標
    f.write(f'\nPOINTS {nodes.shape[0]} float\n')
    for i in range(nodes.shape[0]):
        f.write(SEP.join(map(str,nodes[i,:]))+'\n')
        
    # 寫入單後設資料
    n,size=cal_cell_num(element_Blocks)
    
    # 構造單後設資料
    f.write(f'\nCELLS {n} {size}\n')
    cells_lines,cell_type_lines="",""
    for block ,typeStr in zip(element_Blocks,element_types):
        for i in range(block.shape[0]):
            cells_lines+=f"{block.shape[1]} "+SEP.join(map(str,block[i,:]))+'\n'
            cell_type_lines+=f"{type_dict[typeStr]}\n"
    f.write(cells_lines)
    
    # 寫入單元型別資料
    f.write(f'\nCELL_TYPES {n}\n')
    f.write(cell_type_lines)
    
    # 寫入節點資料,節點字典資料為空則跳過
    if point_data:
        f.write(f'\nPOINT_DATA {nodes.shape[0]}\n') 
        if 'SCALARS' in point_data.keys():
            for var_name,array in point_data['SCALARS'].items():
                if array.ndim==1:
                    # 如果是一維陣列,則擴充套件成二維陣列
                    array=array.reshape((-1,1))
                # 標量的Ncomp=1~4
                f.write(f'SCALARS {var_name} float {array.shape[1]}\n')
                f.write(f'LOOKUP_TABLE default\n')
                for i in range(array.shape[0]):
                    f.write(SEP.join(map(str,array[i,:]))+'\n')
        
        # VECTORS dataName dataType
        # v0x v0y v0z
        # v1x v1y v1z
        # …
        # v(n-1)x v(n-1)y v(n-1)z
        if 'VECTORS' in point_data.keys():
            for var_name,array in point_data['VECTORS'].items():
                if array.ndim==1:
                    # 如果是一維陣列,則擴充套件成二維陣列,向量必須是(-1,3)
                    array=array.reshape((-1,3))
                assert array.shape[1]==3, "VECTORS data must be (-1,3) shape"
                
                f.write(f'VECTORS {var_name} float\n')
                for i in range(array.shape[0]):
                    f.write(SEP.join(map(str,array[i,:]))+'\n')
    else:
        pass
    
    # 寫入單後設資料,單元字典資料為空則跳過
    if cell_data:
        f.write(f'\nCELL_DATA {n}\n')
        if 'SCALARS' in cell_data.keys():
            for var_name,array in cell_data['SCALARS'].items():
                if array.ndim==1:
                    # 如果是一維陣列,則擴充套件成二維陣列
                    array=array.reshape((-1,1))
                # 標量的Ncomp=1~4
                f.write(f'SCALARS {var_name} float {array.shape[1]}\n')
                f.write(f'LOOKUP_TABLE default\n')
                for i in range(array.shape[0]):
                    f.write(SEP.join(map(str,array[i,:]))+'\n')
                    
        if 'VECTORS' in cell_data.keys():
            for var_name,array in cell_data['VECTORS'].items():
                if array.ndim==1:
                    # 如果是一維陣列,則擴充套件成二維陣列,向量必須是(-1,3)
                    array=array.reshape((-1,3))
                assert array.shape[1]==3, "VECTORS data must be (-1,3) shape"
                
                f.write(f'VECTORS {var_name} float\n')
                for i in range(array.shape[0]):
                    f.write(SEP.join(map(str,array[i,:]))+'\n')
    else:
        pass

相關文章