diff --git a/doc/conf.py b/doc/conf.py index f986aa7deb55faf5c786c26ac825b1f32832a8e9..f26f98413f90cf675707dc06a71693eac74a270e 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/examples/vtk-sample-elements.py b/examples/vtk-sample-elements.py new file mode 100644 index 0000000000000000000000000000000000000000..9b8e19545dd41ab66f4d5e5a09329d36b48b17a4 --- /dev/null +++ b/examples/vtk-sample-elements.py @@ -0,0 +1,115 @@ +""" +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 +import numpy.linalg as la + +import matplotlib.pyplot as plt + + +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 # 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: %d" % points.shape[0]) + + ax.grid() + for i, p in enumerate(points.T): + 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 show: + plt.show(block=True) + + +def create_sample_element(cell_type, order=3, visualize=True): + 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() + grid = source.GetOutput() + + 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() + + cell = grid.GetCell(0) + points = vtk_to_numpy(cell.GetPoints().GetData()).T + + dim = cell.GetCellDimension() + 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) + + 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__": + 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/__init__.py b/pyvisfile/vtk/__init__.py index 0abaa022f4fe492e2b6e96091a3e72d38cf11483..6c8ae61173318a996c1a04ec855532dd63ce20fa 100644 --- a/pyvisfile/vtk/__init__.py +++ b/pyvisfile/vtk/__init__.py @@ -38,11 +38,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 ^^^^^^^^^^^^^ @@ -62,6 +62,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 --------------- @@ -87,6 +94,7 @@ Convenience functions .. autofunction:: write_structured_grid """ +# {{{ types VTK_INT8 = "Int8" VTK_UINT8 = "UInt8" @@ -99,6 +107,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 @@ -114,6 +145,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 @@ -129,15 +176,29 @@ 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): def __init__(self): @@ -197,6 +258,13 @@ class XMLRoot(XMLElementBase): # likely a string instance, write it directly file.write(child) +# }}} + + +# {{{ encoded buffers + +_U32CHAR = np.dtype(np.uint32).char + class EncodedBuffer: def encoder(self): @@ -300,6 +368,10 @@ class Base64ZLibEncodedBuffer: return len(self.b64header) + len(self.b64data) +# }}} + + +# {{{ data array class DataArray(object): def __init__(self, name, container, vector_padding=3, @@ -312,48 +384,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: @@ -363,6 +406,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, @@ -374,11 +418,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): @@ -415,6 +463,10 @@ class DataArray(object): def invoke_visitor(self, visitor): return visitor.gen_data_array(self) +# }}} + + +# {{{ grids class UnstructuredGrid(object): """ @@ -450,8 +502,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): @@ -508,8 +559,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" @@ -521,11 +576,27 @@ 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): + """ + :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 @@ -536,13 +607,19 @@ class XMLGenerator(object): else: 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 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) @@ -552,6 +629,7 @@ class XMLGenerator(object): class InlineXMLGenerator(XMLGenerator): """ + .. automethod:: __init__ .. automethod:: __call__ """ @@ -629,8 +707,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 new file mode 100644 index 0000000000000000000000000000000000000000..400fcea2560c75449880abd68e1569b6894fb571 --- /dev/null +++ b/pyvisfile/vtk/vtk_ordering.py @@ -0,0 +1,229 @@ +__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): + return [add_tuples(x, y) for y in ary] + +# }}} + + +# {{{ VTK_LAGRANGE_${SIMPLEX} (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): + """ + :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: + 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_to_index[v] for v in node_tuples] + +# }}} + + +# {{{ VTK_LAGRANGE_${QUAD} (i.e. QUADRILATERAL/HEXAHEDRON) + +# }}}