1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226
|
import logging
import numpy as np
# https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html
vtk_to_meshio_type = {
0: "empty",
1: "vertex",
# 2: 'poly_vertex',
3: "line",
# 4: 'poly_line',
5: "triangle",
# 6: 'triangle_strip',
7: "polygon",
# 8: 'pixel',
9: "quad",
10: "tetra",
# 11: 'voxel',
12: "hexahedron",
13: "wedge",
14: "pyramid",
15: "penta_prism",
16: "hexa_prism",
21: "line3",
22: "triangle6",
23: "quad8",
24: "tetra10",
25: "hexahedron20",
26: "wedge15",
27: "pyramid13",
28: "quad9",
29: "hexahedron27",
30: "quad6",
31: "wedge12",
32: "wedge18",
33: "hexahedron24",
34: "triangle7",
35: "line4",
#
# 60: VTK_HIGHER_ORDER_EDGE,
# 61: VTK_HIGHER_ORDER_TRIANGLE,
# 62: VTK_HIGHER_ORDER_QUAD,
# 63: VTK_HIGHER_ORDER_POLYGON,
# 64: VTK_HIGHER_ORDER_TETRAHEDRON,
# 65: VTK_HIGHER_ORDER_WEDGE,
# 66: VTK_HIGHER_ORDER_PYRAMID,
# 67: VTK_HIGHER_ORDER_HEXAHEDRON,
}
def _get_writer(filetype, filename):
import vtk
if filetype in "vtk-ascii":
logging.warning("VTK ASCII files are only meant for debugging.")
writer = vtk.vtkUnstructuredGridWriter()
writer.SetFileTypeToASCII()
elif filetype == "vtk-binary":
writer = vtk.vtkUnstructuredGridWriter()
writer.SetFileTypeToBinary()
elif filetype == "vtu-ascii":
logging.warning("VTU ASCII files are only meant for debugging.")
writer = vtk.vtkXMLUnstructuredGridWriter()
writer.SetDataModeToAscii()
elif filetype == "vtu-binary":
writer = vtk.vtkXMLUnstructuredGridWriter()
writer.SetDataModeToBinary()
elif filetype == "xdmf2":
writer = vtk.vtkXdmfWriter()
elif filetype == "xdmf3":
writer = vtk.vtkXdmf3Writer()
else:
assert filetype == "exodus", f"Unknown file type '{filename}'."
writer = vtk.vtkExodusIIWriter()
# if the mesh contains vtkmodeldata information, make use of it
# and write out all time steps.
writer.WriteAllTimeStepsOn()
return writer
def write(filetype, filename, mesh):
import vtk
def _create_vtkarray(X, name):
array = vtk.util.numpy_support.numpy_to_vtk(X, deep=1)
array.SetName(name)
return array
# assert data integrity
for val in mesh.point_data.values():
assert len(val) == len(mesh.points), "Point data mismatch."
for key, value in mesh.cell_data.items():
assert key in mesh.cells, "Cell data without cell"
for key2 in value:
assert len(value[key2]) == len(mesh.cells[key]), "Cell data mismatch."
vtk_mesh = _generate_vtk_mesh(mesh.points, mesh.cells)
# add point data
pd = vtk_mesh.GetPointData()
for name, X in mesh.point_data.items():
# There is a naming inconsistency in VTK when it comes to multivectors
# in Exodus files:
# If a vector 'v' has two components, they are called 'v_x', 'v_y'
# (note the underscore), if it has three, then they are called 'vx',
# 'vy', 'vz'. See bug <https://vtk.org/Bug/view.php?id=15894>.
# For VT{K,U} files, no underscore is ever added.
pd.AddArray(_create_vtkarray(X, name))
# Add cell data.
# The cell_data is structured like
#
# cell_type ->
# key -> array
# key -> array
# [...]
# cell_type ->
# key -> array
# key -> array
# [...]
# [...]
#
# VTK expects one array for each `key`, so assemble the keys across all
# mesh_types. This requires each key to be present for each mesh_type, of
# course.
all_keys = []
for cell_type in mesh.cell_data:
all_keys += mesh.cell_data[cell_type].keys()
# create unified cell data
for key in all_keys:
for cell_type in mesh.cell_data:
assert key in mesh.cell_data[cell_type]
unified_cell_data = {
key: np.concatenate([value[key] for value in mesh.cell_data.values()])
for key in all_keys
}
# add the array data to the mesh
cd = vtk_mesh.GetCellData()
for name, array in unified_cell_data.items():
cd.AddArray(_create_vtkarray(array, name))
# add field data
fd = vtk_mesh.GetFieldData()
for key, value in mesh.field_data.items():
fd.AddArray(_create_vtkarray(value, key))
writer = _get_writer(filetype, filename)
writer.SetFileName(filename)
try:
writer.SetInput(vtk_mesh)
except AttributeError:
writer.SetInputData(vtk_mesh)
writer.Write()
return
def _generate_vtk_mesh(points, cells):
import vtk
from vtk.util import numpy as numpy_support
mesh = vtk.vtkUnstructuredGrid()
# set points
vtk_points = vtk.vtkPoints()
# Not using a deep copy here results in a segfault.
vtk_array = numpy_support.numpy_to_vtk(points, deep=True)
vtk_points.SetData(vtk_array)
mesh.SetPoints(vtk_points)
# Set cells.
meshio_to_vtk_type = {y: x for x, y in vtk_to_meshio_type.items()}
# create cell_array. It's a one-dimensional vector with
# (num_points2, p0, p1, ... ,pk, numpoints1, p10, p11, ..., p1k, ...
cell_types = []
cell_offsets = []
cell_connectivity = []
len_array = 0
for meshio_type, data in cells.items():
numcells, num_local_nodes = data.shape
vtk_type = meshio_to_vtk_type[meshio_type]
# add cell types
cell_types.append(np.empty(numcells, dtype=np.ubyte))
cell_types[-1].fill(vtk_type)
# add cell offsets
cell_offsets.append(
np.arange(
len_array,
len_array + numcells * (num_local_nodes + 1),
num_local_nodes + 1,
dtype=np.int64,
)
)
cell_connectivity.append(
np.c_[num_local_nodes * np.ones(numcells, dtype=data.dtype), data].flatten()
)
len_array += len(cell_connectivity[-1])
cell_types = np.concatenate(cell_types)
cell_offsets = np.concatenate(cell_offsets)
cell_connectivity = np.concatenate(cell_connectivity)
connectivity = vtk.util.numpy_support.numpy_to_vtkIdTypeArray(
cell_connectivity.astype(np.int64), deep=1
)
# wrap the data into a vtkCellArray
cell_array = vtk.vtkCellArray()
cell_array.SetCells(len(cell_types), connectivity)
# Add cell data to the mesh
mesh.SetCells(
numpy_support.numpy_to_vtk(
cell_types, deep=1, array_type=vtk.vtkUnsignedCharArray().GetDataType()
),
numpy_support.numpy_to_vtk(
cell_offsets, deep=1, array_type=vtk.vtkIdTypeArray().GetDataType()
),
cell_array,
)
return mesh
|