From e7bec9f4369b50e927a3b662ac9e0324cd2f5fef Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 7 Jun 2020 19:55:54 -0500 Subject: [PATCH 01/13] add VTK_LAGRANGE element ids and some cleanup --- pyvisfile/vtk/__init__.py | 146 ++++++++++++++++++++++++++------------ 1 file changed, 101 insertions(+), 45 deletions(-) diff --git a/pyvisfile/vtk/__init__.py b/pyvisfile/vtk/__init__.py index 5b7ccdb..8050e1f 100644 --- a/pyvisfile/vtk/__init__.py +++ b/pyvisfile/vtk/__init__.py @@ -17,11 +17,11 @@ Vector formats .. data:: VF_LIST_OF_COMPONENTS - ``[[x0,y0,z0], [x1,y1,z1]`` + ``[[x0, y0, z0], [x1, y1, z1]]`` .. data:: VF_LIST_OF_VECTORS - ``[[x0,x1], [y0,y1], [z0,z1]]`` + ``[[x0, x1], [y0, y1], [z0, z1]]`` Element types ^^^^^^^^^^^^^ @@ -41,6 +41,13 @@ Element types .. data:: VTK_WEDGE .. data:: VTK_PYRAMID +.. data:: VTK_LAGRANGE_CURVE +.. data:: VTK_LAGRANGE_TRIANGLE +.. data:: VTK_LAGRANGE_QUADRILATERAL +.. data:: VTK_LAGRANGE_TETRAHEDRON +.. data:: VTK_LAGRANGE_HEXAHEDRON +.. data:: VTK_LAGRANGE_WEDGE + Building blocks --------------- @@ -66,6 +73,7 @@ Convenience functions .. autofunction:: write_structured_grid """ +# {{{ types VTK_INT8 = "Int8" VTK_UINT8 = "UInt8" @@ -78,6 +86,29 @@ VTK_UINT64 = "UInt64" VTK_FLOAT32 = "Float32" VTK_FLOAT64 = "Float64" + +NUMPY_TO_VTK_TYPES = { + np.int8: VTK_INT8, + np.uint8: VTK_UINT8, + np.int16: VTK_INT16, + np.uint16: VTK_UINT16, + np.int32: VTK_INT32, + np.uint32: VTK_UINT32, + np.int64: VTK_INT64, + np.uint64: VTK_UINT64, + np.float32: VTK_FLOAT32, + np.float64: VTK_FLOAT64, + } + +# }}} + + +# {{{ cell types + +# NOTE: should keep in sync with +# https://gitlab.kitware.com/vtk/vtk/-/blob/master/Common/DataModel/vtkCellType.h + +# linear cells VTK_VERTEX = 1 VTK_POLY_VERTEX = 2 VTK_LINE = 3 @@ -93,6 +124,22 @@ VTK_HEXAHEDRON = 12 VTK_WEDGE = 13 VTK_PYRAMID = 14 +# NOTE: these were added in VTK 8.1 as part of the commit +# https://gitlab.kitware.com/vtk/vtk/-/commit/cc5101a805386f205631357bba782b2a7d17531a + +# high-order Lagrange cells +VTK_LAGRANGE_CURVE = 68 +VTK_LAGRANGE_TRIANGLE = 69 +VTK_LAGRANGE_QUADRILATERAL = 70 +VTK_LAGRANGE_TETRAHEDRON = 71 +VTK_LAGRANGE_HEXAHEDRON = 72 +VTK_LAGRANGE_WEDGE = 73 + +# }}} + + +# {{{ cell node counts + CELL_NODE_COUNT = { VTK_VERTEX: 1, # VTK_POLY_VERTEX: no a-priori size @@ -108,14 +155,28 @@ CELL_NODE_COUNT = { VTK_HEXAHEDRON: 8, VTK_WEDGE: 6, VTK_PYRAMID: 5, + # VTK_LAGRANGE_CURVE: no a-priori size + # VTK_LAGRANGE_TRIANGLE: no a-priori size + # VTK_LAGRANGE_QUADRILATERAL: no a-priori size + # VTK_LAGRANGE_TETRAHEDRON: no a-priori size + # VTK_LAGRANGE_HEXAHEDRON: no a-priori size + # VTK_LAGRANGE_WEDGE: no a-priori size } +# }}} -VF_LIST_OF_COMPONENTS = 0 # [[x0,y0,z0], [x1,y1,z1] -VF_LIST_OF_VECTORS = 1 # [[x0,x1], [y0,y1], [z0,z1]] -_U32CHAR = np.dtype(np.uint32).char +# {{{ vector format +# e.g. [[x0, y0, z0], [x1, y1, z1]] +VF_LIST_OF_COMPONENTS = 0 +# e.g. [[x0, x1], [y0, y1], [z0, z1]] +VF_LIST_OF_VECTORS = 1 + +# }}} + + +# {{{ xml # Ah, the joys of home-baked non-compliant XML goodness. class XMLElementBase(object): @@ -176,6 +237,13 @@ class XMLRoot(XMLElementBase): # likely a string instance, write it directly file.write(child) +# }}} + + +# {{{ + +_U32CHAR = np.dtype(np.uint32).char + class EncodedBuffer: def encoder(self): @@ -291,48 +359,19 @@ class DataArray(object): self.encoded_buffer = container.encoded_buffer return - def vec_type(vec): - if vec.dtype == np.int8: - return VTK_INT8 - elif vec.dtype == np.uint8: - return VTK_UINT8 - elif vec.dtype == np.int16: - return VTK_INT16 - elif vec.dtype == np.uint16: - return VTK_UINT16 - elif vec.dtype == np.int32: - return VTK_INT32 - elif vec.dtype == np.uint32: - return VTK_UINT32 - elif vec.dtype == np.int64: - return VTK_INT64 - elif vec.dtype == np.uint64: - return VTK_UINT64 - elif vec.dtype == np.float32: - return VTK_FLOAT32 - elif vec.dtype == np.float64: - return VTK_FLOAT64 - else: - raise TypeError( - "Unsupported vector type '%s' in VTK writer" % (vec.dtype)) - if not isinstance(container, np.ndarray): raise ValueError( "cannot convert object of type `%s' to DataArray" % type(container)) - if not isinstance(container, np.ndarray): - raise TypeError("expected numpy array, got '%s' instead" - % type(container)) - - if container.dtype == object: + if container.dtype.char == "O": for subvec in container: if not isinstance(subvec, np.ndarray): raise TypeError("expected numpy array, got '%s' instead" % type(subvec)) container = np.array(list(container)) - assert container.dtype != object + assert container.dtype.char != "O" if len(container.shape) > 1: if vector_format == VF_LIST_OF_COMPONENTS: @@ -342,6 +381,7 @@ class DataArray(object): "numpy vectors of rank >2 are not supported" assert container.strides[1] == container.itemsize, \ "2D numpy arrays must be row-major" + if vector_padding > container.shape[1]: container = np.asarray(np.hstack(( container, @@ -353,11 +393,15 @@ class DataArray(object): self.components = container.shape[1] else: self.components = 1 - self.type = vec_type(container) + + self.type = NUMPY_TO_VTK_TYPES.get(container.dtype.type, None) + if self.type is None: + raise TypeError("unsupported vector type: '%s'" % (container.dtype)) + if not container.flags.c_contiguous: container = container.copy() - buf = buffer(container) + buf = buffer(container) self.encoded_buffer = BinaryEncodedBuffer(buf) def get_encoded_buffer(self, encoder, compressor): @@ -429,8 +473,7 @@ class UnstructuredGrid(object): def copy(self): return UnstructuredGrid( (self.point_count, self.points), - (self.cell_count, self.cell_connectivity, - self.cell_offsets), + (self.cell_count, self.cell_connectivity, self.cell_offsets), self.cell_types) def vtk_extension(self): @@ -596,8 +639,15 @@ class InlineXMLGenerator(XMLGenerator): return el def gen_data_array(self, data): - el = XMLElement("DataArray", type=data.type, Name=data.name, - NumberOfComponents=data.components, format="binary") + keys = { + "type": data.type, + "Name": data.name, + "format": "binary" + } + if data.components > 1: + keys["NumberOfComponents"] = data.components + + el = XMLElement("DataArray", **keys) data.encode(self.compressor, el) el.add_child("\n") return el @@ -624,10 +674,16 @@ class AppendedDataXMLGenerator(InlineXMLGenerator): return xmlroot def gen_data_array(self, data): - el = XMLElement("DataArray", type=data.type, Name=data.name, - NumberOfComponents=data.components, format="appended", - offset=self.base64_len) + keys = { + "type": data.type, + "Name": data.name, + "format": "appended", + "offset": self.base64_len, + } + if data.components > 1: + keys["NumberOfComponents"] = data.components + el = XMLElement("DataArray", **keys) self.base64_len += data.encode(self.compressor, self.app_data) return el -- GitLab From d5323a17365caeeb28c08c445852a4726853f5d6 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 7 Jun 2020 19:56:16 -0500 Subject: [PATCH 02/13] examples: plot vtk elements using vtk itself to get numbering --- examples/vtk-sample-elements.py | 74 +++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 examples/vtk-sample-elements.py diff --git a/examples/vtk-sample-elements.py b/examples/vtk-sample-elements.py new file mode 100644 index 0000000..62f2be3 --- /dev/null +++ b/examples/vtk-sample-elements.py @@ -0,0 +1,74 @@ +import numpy as np +import matplotlib.pyplot as plt + + +def plot_node_ordering(filename, points): + if points.shape[0] == 2: + fig = plt.figure(figsize=(8, 8), dpi=300) + ax = fig.gca() + ax.plot(points[0], points[1], "o") + elif points.shape[0] == 3: + from mpl_toolkits.mplot3d import art3d + + fig = plt.figure(figsize=(8, 8), dpi=300) + ax = fig.gca(projection="3d") + ax.plot(points[0], points[1], points[2], "o") + else: + raise ValueError("dimension not supported") + + ax.grid() + for i, p in enumerate(points.T): + ax.text(*p, str(i), color="k") + + print("output: %s.png" % filename) + fig.savefig(filename) + if 0: + plt.show(block=True) + + +def create_sample_element(cell_type, order=3): + try: + import vtk + from vtkmodules.util.numpy_support import vtk_to_numpy + except ImportError: + raise ImportError("python bindings for vtk cannot be found") + + cell_type = cell_type.upper() + vtk_cell_type = getattr(vtk, cell_type, None) + if vtk_cell_type is None: + raise ValueError("unknown cell type: '%s'" % cell_type) + + source = vtk.vtkCellTypeSource() + source.SetCellType(vtk_cell_type) + source.SetBlocksDimensions(1, 1, 1) + if "LAGRANGE" in cell_type: + source.SetCellOrder(order) + source.Update() + + filename = "sample_%s.vtu" % cell_type.lower() + print("cell type: %s" % cell_type) + print("output: %s" % filename) + + grid = source.GetOutput() + + # write vtk file + writer = vtk.vtkXMLUnstructuredGridWriter() + writer.SetFileName(filename) + writer.SetCompressorTypeToNone() + writer.SetDataModeToAscii() + writer.SetInputData(grid) + writer.Write() + + # write numbered matplotlib file + cell = grid.GetCell(0) + points = vtk_to_numpy(cell.GetPoints().GetData()).T + + dim = cell.GetCellDimension() + points = points[0:max(dim, 2)] + + filename = "sample_%s" % cell_type.lower() + plot_node_ordering(filename, points) + + +if __name__ == "__main__": + create_sample_element("VTK_LAGRANGE_TRIANGLE", order=5) -- GitLab From a4021ed61eb66bcf28381445dc9c706611e7baff Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 9 Jun 2020 21:43:07 -0500 Subject: [PATCH 03/13] add some ordering maps for lagrange elements --- examples/vtk-sample-elements.py | 79 ++++++++++----- pyvisfile/vtk/vtk_ordering.py | 167 ++++++++++++++++++++++++++++++++ 2 files changed, 223 insertions(+), 23 deletions(-) create mode 100644 pyvisfile/vtk/vtk_ordering.py diff --git a/examples/vtk-sample-elements.py b/examples/vtk-sample-elements.py index 62f2be3..fd99bbb 100644 --- a/examples/vtk-sample-elements.py +++ b/examples/vtk-sample-elements.py @@ -1,32 +1,47 @@ import numpy as np +import numpy.linalg as la + import matplotlib.pyplot as plt -def plot_node_ordering(filename, points): +VTK_LAGRANGE_SIMPLICES = [ + "VTK_LAGRANGE_CURVE", + "VTK_LAGRANGE_TRIANGLE", + "VTK_LAGRANGE_TETRAHEDRON", + ] + + +def plot_node_ordering(filename, points, show=False): + if points.shape[0] == 1: + points = np.hstack([points, np.zeros_like(points)]) + if points.shape[0] == 2: fig = plt.figure(figsize=(8, 8), dpi=300) ax = fig.gca() ax.plot(points[0], points[1], "o") elif points.shape[0] == 3: - from mpl_toolkits.mplot3d import art3d + from mpl_toolkits.mplot3d import art3d # noqa: F401 fig = plt.figure(figsize=(8, 8), dpi=300) ax = fig.gca(projection="3d") ax.plot(points[0], points[1], points[2], "o") + + ax.view_init(0, 90) else: - raise ValueError("dimension not supported") + raise ValueError("dimension not supported: %d" % points.shape[0]) ax.grid() for i, p in enumerate(points.T): - ax.text(*p, str(i), color="k") + if abs(p[1]) < 1.0e-14: + ax.text(*p, str(i), color="k", fontsize=12) print("output: %s.png" % filename) fig.savefig(filename) - if 0: + if show: plt.show(block=True) -def create_sample_element(cell_type, order=3): +def create_sample_element(cell_type, order=3, visualize=True): try: import vtk from vtkmodules.util.numpy_support import vtk_to_numpy @@ -44,31 +59,49 @@ def create_sample_element(cell_type, order=3): if "LAGRANGE" in cell_type: source.SetCellOrder(order) source.Update() - - filename = "sample_%s.vtu" % cell_type.lower() - print("cell type: %s" % cell_type) - print("output: %s" % filename) - grid = source.GetOutput() - # write vtk file - writer = vtk.vtkXMLUnstructuredGridWriter() - writer.SetFileName(filename) - writer.SetCompressorTypeToNone() - writer.SetDataModeToAscii() - writer.SetInputData(grid) - writer.Write() + if visualize: + filename = "sample_%s.vtu" % cell_type.lower() + print("cell type: %s" % cell_type) + print("output: %s" % filename) + + writer = vtk.vtkXMLUnstructuredGridWriter() + writer.SetFileName(filename) + writer.SetCompressorTypeToNone() + writer.SetDataModeToAscii() + writer.SetInputData(grid) + writer.Write() - # write numbered matplotlib file cell = grid.GetCell(0) points = vtk_to_numpy(cell.GetPoints().GetData()).T dim = cell.GetCellDimension() - points = points[0:max(dim, 2)] + points = points[0:dim] + + if cell_type in VTK_LAGRANGE_SIMPLICES: + from pyvisfile.vtk.vtk_ordering import vtk_lagrange_simplex_node_tuples + nodes = np.array(vtk_lagrange_simplex_node_tuples(dim, order)) + + error = la.norm(nodes/order - points.T) + print("error[%d]: %.5e" % (order, error)) + assert error < 5.0e-7 + + if visualize: + filename = "sample_%s" % cell_type.lower() + plot_node_ordering(filename, points, show=False) - filename = "sample_%s" % cell_type.lower() - plot_node_ordering(filename, points) + if cell_type in VTK_LAGRANGE_SIMPLICES: + filename = "sample_%s_new" % cell_type.lower() + plot_node_ordering(filename, nodes.T, show=False) if __name__ == "__main__": - create_sample_element("VTK_LAGRANGE_TRIANGLE", order=5) + if 1: + for cell_type in VTK_LAGRANGE_SIMPLICES: + print("cell_type: ", cell_type) + for order in range(1, 11): + create_sample_element( + cell_type, order=order, visualize=False) + + create_sample_element("VTK_LAGRANGE_TETRAHEDRON", order=4) diff --git a/pyvisfile/vtk/vtk_ordering.py b/pyvisfile/vtk/vtk_ordering.py new file mode 100644 index 0000000..17de1b8 --- /dev/null +++ b/pyvisfile/vtk/vtk_ordering.py @@ -0,0 +1,167 @@ +from pytools import ( + add_tuples, + generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam) + + +# {{{ + +def add_tuple_to_list(ary, x): + return [add_tuples(x, y) for y in ary] + +# }}} + + +# {{{ VTK_LAGRANGE_SIMPLICES (i.e. CURVE/TRIANGLE/TETRAHEDRON) + +def vtk_lagrange_curve_node_tuples(order, is_consistent=False): + if is_consistent: + node_tuples = [(0,), (order,)] + [(i,) for i in range(1, order)] + else: + node_tuples = [(i,) for i in range(order + 1)] + + return node_tuples + + +def vtk_lagrange_triangle_node_tuples(order): + nodes = [] + offset = (0, 0) + + if order < 0: + return nodes + + while order >= 0: + if order == 0: + nodes += [offset] + break + + # y + # ^ + # | + # 2 + # |`\ + # | `\ + # | `\ + # | `\ + # | `\ + # 0----------1 --> x + + # add vertices + vertices = [(0, 0), (order, 0), (0, order)] + nodes += add_tuple_to_list(vertices, offset) + if order == 1: + break + + # add faces + face_ids = range(1, order) + faces = ( + # vertex 0 -> 1 + [(i, 0) for i in face_ids] + # vertex 1 -> 2 + + [(order - i, i) for i in face_ids] + # vertex 2 -> 0 + + [(0, order - i) for i in face_ids]) + nodes += add_tuple_to_list(faces, offset) + + order = order - 3 + offset = add_tuples(offset, (1, 1)) + + return nodes + + +def vtk_lagrange_tetrahedron_node_tuples(order): + nodes = [] + offset = (0, 0, 0) + + while order >= 0: + if order == 0: + nodes += [offset] + break + + # z + # ,/ + # 3 + # ,/|`\ + # ,/ | `\ + # ,/ '. `\ + # ,/ | `\ + # ,/ | `\ + # 0-----------'.--------2---> y + # `\. | ,/ + # `\. | ,/ + # `\. '. ,/ + # `\. |/ + # `1 + # `. + # x + + vertices = [(0, 0, 0), (order, 0, 0), (0, order, 0), (0, 0, order)] + nodes += add_tuple_to_list(vertices, offset) + if order == 1: + break + + # add edges + edge_ids = range(1, order) + edges = ( + # vertex 0 -> 1 + [(i, 0, 0) for i in edge_ids] + # vertex 1 -> 2 + + [(order - i, i, 0) for i in edge_ids] + # vertex 2 -> 0 + + [(0, order - i, 0) for i in edge_ids] + # vertex 0 -> 3 + + [(0, 0, i) for i in edge_ids] + # vertex 1 -> 3 + + [(order - i, 0, i) for i in edge_ids] + # vertex 2 -> 3 + + [(0, order - i, i) for i in edge_ids]) + nodes += add_tuple_to_list(edges, offset) + + # add faces + face_ids = add_tuple_to_list( + vtk_lagrange_triangle_node_tuples(order - 3), (1, 1)) + faces = ( + # face between vertices (0, 2, 3) + [(i, 0, j) for i, j in face_ids] + # face between vertices (1, 2, 3) + + [(j, order - (i + j), i) for i, j in face_ids] + # face between vertices (0, 1, 3) + + [(0, j, i) for i, j in face_ids] + # face between vertices (0, 1, 2) + + [(j, i, 0) for i, j in face_ids]) + nodes += add_tuple_to_list(faces, offset) + + order = order - 4 + offset = add_tuples(offset, (1, 1, 1)) + + return nodes + + +def vtk_lagrange_simplex_node_tuples(dims, order, is_consistent=False): + if dims == 1: + return vtk_lagrange_curve_node_tuples(order, is_consistent=is_consistent) + elif dims == 2: + return vtk_lagrange_triangle_node_tuples(order) + elif dims == 3: + return vtk_lagrange_tetrahedron_node_tuples(order) + else: + raise ValueError("unsupported dimension: %d" % dims) + + +def vtk_lagrange_simplex_node_tuples_to_permutation(node_tuples): + order = max([sum(i) for i in node_tuples]) + dims = len(node_tuples[0]) + + node_to_index = dict( + (node_tuple, i) + for i, node_tuple in enumerate(gnitstam(order, dims)) + ) + + assert len(node_tuples) == len(node_to_index) + return [node_tuples[v] for v in node_tuples] + +# }}} + + +# {{{ VTK_LAGRANGE_QUADS (i.e. QUADRILATERAL/HEXAHEDRON) + +# }}} -- GitLab From e1c9b5552c7994ff30e3eaa4d106b5b53b89464b Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 9 Jun 2020 22:03:59 -0500 Subject: [PATCH 04/13] call xml attributes attributes --- pyvisfile/vtk/__init__.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pyvisfile/vtk/__init__.py b/pyvisfile/vtk/__init__.py index 8050e1f..5fb7791 100644 --- a/pyvisfile/vtk/__init__.py +++ b/pyvisfile/vtk/__init__.py @@ -639,15 +639,15 @@ class InlineXMLGenerator(XMLGenerator): return el def gen_data_array(self, data): - keys = { + attrs = { "type": data.type, "Name": data.name, "format": "binary" } if data.components > 1: - keys["NumberOfComponents"] = data.components + attrs["NumberOfComponents"] = data.components - el = XMLElement("DataArray", **keys) + el = XMLElement("DataArray", **attrs) data.encode(self.compressor, el) el.add_child("\n") return el @@ -674,16 +674,16 @@ class AppendedDataXMLGenerator(InlineXMLGenerator): return xmlroot def gen_data_array(self, data): - keys = { + attrs = { "type": data.type, "Name": data.name, "format": "appended", "offset": self.base64_len, } if data.components > 1: - keys["NumberOfComponents"] = data.components + attrs["NumberOfComponents"] = data.components - el = XMLElement("DataArray", **keys) + el = XMLElement("DataArray", **attrs) self.base64_len += data.encode(self.compressor, self.app_data) return el -- GitLab From 59fa6ce85f241146c6de49073ff7b682f6669c74 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 9 Jun 2020 22:15:04 -0500 Subject: [PATCH 05/13] pass VTKFile version to xml generators --- pyvisfile/vtk/__init__.py | 31 ++++++++++++++++++++++++------- pyvisfile/vtk/vtk_ordering.py | 2 +- 2 files changed, 25 insertions(+), 8 deletions(-) diff --git a/pyvisfile/vtk/__init__.py b/pyvisfile/vtk/__init__.py index 5fb7791..8798f0a 100644 --- a/pyvisfile/vtk/__init__.py +++ b/pyvisfile/vtk/__init__.py @@ -240,7 +240,7 @@ class XMLRoot(XMLElementBase): # }}} -# {{{ +# {{{ encoded buffers _U32CHAR = np.dtype(np.uint32).char @@ -347,6 +347,10 @@ class Base64ZLibEncodedBuffer: return len(self.b64header) + len(self.b64data) +# }}} + + +# {{{ data array class DataArray(object): def __init__(self, name, container, vector_padding=3, @@ -438,6 +442,10 @@ class DataArray(object): def invoke_visitor(self, visitor): return visitor.gen_data_array(self) +# }}} + + +# {{{ grids class UnstructuredGrid(object): """ @@ -530,8 +538,12 @@ class StructuredGrid(object): def add_celldata(self, data_array): self.celldata.append(data_array) +# }}} + -def make_vtkfile(filetype, compressor): +# {{{ vtk xml writers + +def make_vtkfile(filetype, compressor, version="0.1"): import sys if sys.byteorder == "little": bo = "LittleEndian" @@ -543,11 +555,11 @@ def make_vtkfile(filetype, compressor): kwargs["compressor"] = "vtkZLibDataCompressor" return XMLElement("VTKFile", - type=filetype, version="0.1", byte_order=bo, **kwargs) + type=filetype, version=version, byte_order=bo, **kwargs) class XMLGenerator(object): - def __init__(self, compressor=None): + def __init__(self, compressor=None, vtk_file_version=None): if compressor == "zlib": try: import zlib # noqa @@ -558,13 +570,18 @@ class XMLGenerator(object): else: raise ValueError("Invalid compressor name `%s'" % compressor) + if vtk_file_version is None: + vtk_file_version = "0.1" + + self.vtk_file_version = vtk_file_version self.compressor = compressor def __call__(self, vtkobj): """Return an :class:`XMLElement`.""" child = self.rec(vtkobj) - vtkf = make_vtkfile(child.tag, self.compressor) + vtkf = make_vtkfile(child.tag, self.compressor, + version=self.vtk_file_version) vtkf.add_child(child) return XMLRoot(vtkf) @@ -658,8 +675,8 @@ class AppendedDataXMLGenerator(InlineXMLGenerator): .. automethod:: __call__ """ - def __init__(self, compressor=None): - InlineXMLGenerator.__init__(self, compressor) + def __init__(self, compressor=None, vtk_file_version=None): + InlineXMLGenerator.__init__(self, compressor, vtk_file_version) self.base64_len = 0 self.app_data = XMLElement("AppendedData", encoding="base64") diff --git a/pyvisfile/vtk/vtk_ordering.py b/pyvisfile/vtk/vtk_ordering.py index 17de1b8..90bdeb8 100644 --- a/pyvisfile/vtk/vtk_ordering.py +++ b/pyvisfile/vtk/vtk_ordering.py @@ -157,7 +157,7 @@ def vtk_lagrange_simplex_node_tuples_to_permutation(node_tuples): ) assert len(node_tuples) == len(node_to_index) - return [node_tuples[v] for v in node_tuples] + return [node_to_index[v] for v in node_tuples] # }}} -- GitLab From c4eebaac5b9f8a1551cf958af79658a564cb6e6d Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 16 Jun 2020 01:30:11 +0200 Subject: [PATCH 06/13] Apply suggestion to pyvisfile/vtk/__init__.py --- pyvisfile/vtk/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyvisfile/vtk/__init__.py b/pyvisfile/vtk/__init__.py index 8798f0a..e191c2e 100644 --- a/pyvisfile/vtk/__init__.py +++ b/pyvisfile/vtk/__init__.py @@ -571,6 +571,7 @@ class XMLGenerator(object): raise ValueError("Invalid compressor name `%s'" % compressor) if vtk_file_version is None: + # https://www.paraview.org/Wiki/VTK_XML_Formats vtk_file_version = "0.1" self.vtk_file_version = vtk_file_version -- GitLab From ce49388866938088ea7921736f4619de60a20320 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 16 Jun 2020 01:30:18 +0200 Subject: [PATCH 07/13] Apply suggestion to pyvisfile/vtk/vtk_ordering.py --- pyvisfile/vtk/vtk_ordering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyvisfile/vtk/vtk_ordering.py b/pyvisfile/vtk/vtk_ordering.py index 90bdeb8..a031730 100644 --- a/pyvisfile/vtk/vtk_ordering.py +++ b/pyvisfile/vtk/vtk_ordering.py @@ -11,7 +11,7 @@ def add_tuple_to_list(ary, x): # }}} -# {{{ VTK_LAGRANGE_SIMPLICES (i.e. CURVE/TRIANGLE/TETRAHEDRON) +# {{{ VTK_LAGRANGE_${SIMPLEX} (i.e. CURVE/TRIANGLE/TETRAHEDRON) def vtk_lagrange_curve_node_tuples(order, is_consistent=False): if is_consistent: -- GitLab From 68b5ee37bc6f9b54e33daae907e2bdb0350a2c0f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 16 Jun 2020 01:30:25 +0200 Subject: [PATCH 08/13] Apply suggestion to pyvisfile/vtk/vtk_ordering.py --- pyvisfile/vtk/vtk_ordering.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyvisfile/vtk/vtk_ordering.py b/pyvisfile/vtk/vtk_ordering.py index a031730..5a03cc9 100644 --- a/pyvisfile/vtk/vtk_ordering.py +++ b/pyvisfile/vtk/vtk_ordering.py @@ -162,6 +162,6 @@ def vtk_lagrange_simplex_node_tuples_to_permutation(node_tuples): # }}} -# {{{ VTK_LAGRANGE_QUADS (i.e. QUADRILATERAL/HEXAHEDRON) +# {{{ VTK_LAGRANGE_${QUAD} (i.e. QUADRILATERAL/HEXAHEDRON) # }}} -- GitLab From 6a95bad9f77ab2474197603ca2aab880b1ff9332 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 15 Jun 2020 18:41:58 -0500 Subject: [PATCH 09/13] restore NumberOfComponents for scalar data --- pyvisfile/vtk/__init__.py | 23 +++++------------------ 1 file changed, 5 insertions(+), 18 deletions(-) diff --git a/pyvisfile/vtk/__init__.py b/pyvisfile/vtk/__init__.py index c4fe2e9..544b687 100644 --- a/pyvisfile/vtk/__init__.py +++ b/pyvisfile/vtk/__init__.py @@ -678,15 +678,8 @@ class InlineXMLGenerator(XMLGenerator): return el def gen_data_array(self, data): - attrs = { - "type": data.type, - "Name": data.name, - "format": "binary" - } - if data.components > 1: - attrs["NumberOfComponents"] = data.components - - el = XMLElement("DataArray", **attrs) + el = XMLElement("DataArray", type=data.type, Name=data.name, + NumberOfComponents=data.components, format="binary") data.encode(self.compressor, el) el.add_child("\n") return el @@ -713,16 +706,10 @@ class AppendedDataXMLGenerator(InlineXMLGenerator): return xmlroot def gen_data_array(self, data): - attrs = { - "type": data.type, - "Name": data.name, - "format": "appended", - "offset": self.base64_len, - } - if data.components > 1: - attrs["NumberOfComponents"] = data.components + el = XMLElement("DataArray", type=data.type, Name=data.name, + NumberOfComponents=data.components, format="appended", + offset=self.base64_len) - el = XMLElement("DataArray", **attrs) self.base64_len += data.encode(self.compressor, self.app_data) return el -- GitLab From d206d7908146225dc6cd17ae4a3773a1f1077400 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 15 Jun 2020 19:09:43 -0500 Subject: [PATCH 10/13] add some docs to vtk ordering file --- doc/conf.py | 10 +++++- pyvisfile/vtk/vtk_ordering.py | 62 +++++++++++++++++++++++++++++++++++ 2 files changed, 71 insertions(+), 1 deletion(-) diff --git a/doc/conf.py b/doc/conf.py index f986aa7..f26f984 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -23,6 +23,8 @@ # coming with Sphinx (named 'sphinx.ext.*') or your custom ones. extensions = [ "sphinx.ext.autodoc", + "sphinx.ext.intersphinx", + "sphinx.ext.mathjax", ] # Add any paths that contain templates here, relative to this directory. @@ -145,7 +147,7 @@ htmlhelp_basename = 'PyVisfiledoc' # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, author, document class [howto/manual]). -latex_documents = [("index", "pyvisfile.tex", +latex_documents = [("index", "pyvisfile.tex", "PyVisfile documentation", "Andreas Kloeckner", "manual")] @@ -157,3 +159,9 @@ latex_documents = [("index", "pyvisfile.tex", # If false, no module index is generated. #latex_use_modindex = True + +intersphinx_mapping = { + 'http://docs.python.org/': None, + 'http://docs.scipy.org/doc/numpy/': None, + 'https://documen.tician.de/modepy': None, +} diff --git a/pyvisfile/vtk/vtk_ordering.py b/pyvisfile/vtk/vtk_ordering.py index 5a03cc9..400fcea 100644 --- a/pyvisfile/vtk/vtk_ordering.py +++ b/pyvisfile/vtk/vtk_ordering.py @@ -1,8 +1,54 @@ +__copyright__ = "Copyright (C) 2020 Alexandru Fikl" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + from pytools import ( add_tuples, generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam) +__doc__ = """ +VTK High-Order Lagrange Elements +-------------------------------- + +The high-order elements are described in +`this blog post `_. The ordering in the new elements is as +follows: + + 1. vertices of the element (in an order that matches the linear elements, + e.g. :data:`~pyvisfile.vtk.VTK_TRIANGLE`) + 2. Edge (or face in 2D) nodes in sequence. + 3. Face (3D only) nodes. These are only the nodes interior to the face, + i.e. without the edges, and they are reported by the same rules + recursively. + 4. Interior nodes are also defined recursively. + +To a large extent, matches the order used by ``gmsh`` and described +`here `_. + +.. autofunction:: vtk_lagrange_simplex_node_tuples +.. autofunction:: vtk_lagrange_simplex_node_tuples_to_permutation +""" # noqa + + # {{{ def add_tuple_to_list(ary, x): @@ -137,6 +183,22 @@ def vtk_lagrange_tetrahedron_node_tuples(order): def vtk_lagrange_simplex_node_tuples(dims, order, is_consistent=False): + """ + :arg dims: dimension of the simplex, i.e. 1 corresponds to a curve, 2 to + a triangle, etc. + :arg order: order of the polynomial representation, which also defines + the number of nodes on the simplex. + :arg is_consistent: If *True*, 1D curve node ordering will follow the + same rules as the higher-dimensional simplices by putting the + vertices first. This is not the default in VTK as of version 8.1, + when the higher-order elements were introduced. + + :return: a :class:`list` of ``dims``-dimensional tuples of integers + up to ``order`` in the ordering expected by VTK. This list can be + passed to :func:`vtk_lagrange_simplex_node_tuples_to_permutation` + to obtain a permutation from the order used by :mod:`modepy`. + """ + if dims == 1: return vtk_lagrange_curve_node_tuples(order, is_consistent=is_consistent) elif dims == 2: -- GitLab From 21097fa41cc43e1c4f0947c89e7d50afb37567ed Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 15 Jun 2020 19:37:09 -0500 Subject: [PATCH 11/13] add some information about the vtk version --- pyvisfile/vtk/__init__.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/pyvisfile/vtk/__init__.py b/pyvisfile/vtk/__init__.py index 544b687..6c8ae61 100644 --- a/pyvisfile/vtk/__init__.py +++ b/pyvisfile/vtk/__init__.py @@ -581,6 +581,22 @@ def make_vtkfile(filetype, compressor, version="0.1"): class XMLGenerator(object): def __init__(self, compressor=None, vtk_file_version=None): + """ + :arg vtk_file_version: a string ``"x.y"`` with the desired VTK + XML file format version. Relevant versions are as follows: + + * ``"0.1"`` is the original version. + * ``"1.0"`` added support for 64-bit indices and offsets, as + described `here `_. + * ``"2.0"`` added support for ghost array data, as + described `here `_. + * ``"2.1"``: added support for writing additional information + attached to a :class:`DataArray` using + `information keys `_. + * ``"2.2"``: changed the node numbering of the hexahedron, as + described `here `_. + """ # noqa + if compressor == "zlib": try: import zlib # noqa @@ -613,6 +629,7 @@ class XMLGenerator(object): class InlineXMLGenerator(XMLGenerator): """ + .. automethod:: __init__ .. automethod:: __call__ """ -- GitLab From 6c8035e6b520732b304bb38b9fa823e0211946c1 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Mon, 15 Jun 2020 20:28:41 -0500 Subject: [PATCH 12/13] add some note to the example --- examples/vtk-sample-elements.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/examples/vtk-sample-elements.py b/examples/vtk-sample-elements.py index fd99bbb..4ed164a 100644 --- a/examples/vtk-sample-elements.py +++ b/examples/vtk-sample-elements.py @@ -1,3 +1,9 @@ +""" +This example can be used to check if the VTK ordering and the ones implemented +in :mod:`pyvisfil.vtk.vtk_ordering` match. It can be useful for debugging +and implementing additional element types. +""" + import numpy as np import numpy.linalg as la @@ -97,11 +103,10 @@ def create_sample_element(cell_type, order=3, visualize=True): if __name__ == "__main__": - if 1: - for cell_type in VTK_LAGRANGE_SIMPLICES: - print("cell_type: ", cell_type) - for order in range(1, 11): - create_sample_element( + for cell_type in VTK_LAGRANGE_SIMPLICES: + print("cell_type: ", cell_type) + for order in range(1, 11): + create_sample_element( cell_type, order=order, visualize=False) create_sample_element("VTK_LAGRANGE_TETRAHEDRON", order=4) -- GitLab From c30f61f506f4309797d683cdb6d3437c9fdc0f28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Tue, 16 Jun 2020 03:51:24 +0200 Subject: [PATCH 13/13] Apply suggestion to examples/vtk-sample-elements.py --- examples/vtk-sample-elements.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/vtk-sample-elements.py b/examples/vtk-sample-elements.py index 4ed164a..9b8e195 100644 --- a/examples/vtk-sample-elements.py +++ b/examples/vtk-sample-elements.py @@ -2,6 +2,9 @@ This example can be used to check if the VTK ordering and the ones implemented in :mod:`pyvisfil.vtk.vtk_ordering` match. It can be useful for debugging and implementing additional element types. + +To facilitate this comparison, and unlike the rest of the package, this example +makes use of the Vtk Python bindings. """ import numpy as np -- GitLab