diff --git a/README.md b/README.md index b04037f..090bb22 100644 --- a/README.md +++ b/README.md @@ -1,234 +1,236 @@ UVW - Universal VTK Writer ========================== [![Build Status](https://travis-ci.com/prs513rosewood/uvw.svg?branch=master)](https://travis-ci.com/prs513rosewood/uvw) [![Coverage Status](https://coveralls.io/repos/github/prs513rosewood/uvw/badge.svg?branch=master)](https://coveralls.io/github/prs513rosewood/uvw?branch=master) UVW is a small utility library to write VTK files from data contained in Numpy arrays. It handles fully-fledged `ndarrays` defined over {1, 2, 3}-d domains, with arbitrary number of components. There are no constraints on the particular order of components, although copy of data can be avoided if the array is Fortran contiguous, as VTK files are written in Fortran order. UVW supports multi-process writing of VTK files, so that it can be used in an MPI environment. ## Getting Started Here is how to install and use `uvw`. ### Prerequisites * Python 3. It may work with python 2, but it hasn't been tested. * [Numpy](http://www.numpy.org/). This code has been tested with Numpy version 1.14.3. * [mpi4py](https://mpi4py.readthedocs.io/en/stable/) only if you wish to use the parallel classes of UVW (i.e. the submodule `uvw.parallel`) ### Installing This library can be installed with `pip`: ``` pip install --user uvw ``` If you want to activate parallel capabilities, run: ``` pip install --user uvw[mpi] ``` which will automatically pull `mpi4py` as a dependency. ### Writing Numpy arrays As a first example, let us write a multi-component numpy array into a rectilinear grid: ```python import numpy as np from uvw import RectilinearGrid, DataArray # Creating coordinates x = np.linspace(-0.5, 0.5, 10) y = np.linspace(-0.5, 0.5, 20) z = np.linspace(-0.9, 0.9, 30) # Creating the file (with possible data compression) grid = RectilinearGrid('grid.vtr', (x, y, z), compression=True) # A centered ball x, y, z = np.meshgrid(x, y, z, indexing='ij') r = np.sqrt(x**2 + y**2 + z**2) ball = r < 0.3 # Some multi-component multi-dimensional data data = np.zeros([10, 20, 30, 3, 3]) data[ball, ...] = np.array([[0, 1, 0], [1, 0, 0], [0, 1, 1]]) # Some cell data cell_data = np.zeros([9, 19, 29]) cell_data[0::2, 0::2, 0::2] = 1 # Adding the point data (see help(DataArray) for more info) grid.addPointData(DataArray(data, range(3), 'ball')) # Adding the cell data grid.addCellData(DataArray(cell_data, range(3), 'checkers')) grid.write() ``` UVW also supports writing data on 2D and 1D physical domains, for example: ```python import sys import numpy as np from uvw import RectilinearGrid, DataArray # Creating coordinates x = np.linspace(-0.5, 0.5, 10) y = np.linspace(-0.5, 0.5, 20) # A centered disk xx, yy = np.meshgrid(x, y, indexing='ij') r = np.sqrt(xx**2 + yy**2) R = 0.3 disk = r < R data = np.zeros([10, 20]) data[disk] = np.sqrt(1-(r[disk]/R)**2) # File object can be used as a context manager # and you can write to stdout! with RectilinearGrid(sys.stdout, (x, y)) as grid: grid.addPointData(DataArray(data, range(2), 'data')) ``` ## Writing in parallel with `mpi4py` The classes contained in the `uvw.parallel` submodule support multi-process writing using `mpi4py`. Here is a code example: ```python import numpy as np from mpi4py import MPI from uvw.parallel import PRectilinearGrid from uvw import DataArray comm = MPI.COMM_WORLD rank = comm.Get_rank() N = 20 # Domain bounds per rank bounds = [ {'x': (-2, 0), 'y': (-2, 0)}, {'x': (-2, 0), 'y': (0, 2)}, {'x': (0, 2), 'y': (-2, 2)}, ] # Domain sizes per rank sizes = [ {'x': N, 'y': N}, {'x': N, 'y': N}, {'x': N, 'y': 2*N-1}, # account for overlap ] # Size offsets per rank offsets = [ [0, 0], [0, N], [N, 0], ] x = np.linspace(*bounds[rank]['x'], sizes[rank]['x']) y = np.linspace(*bounds[rank]['y'], sizes[rank]['y']) xx, yy = np.meshgrid(x, y, indexing='ij', sparse=True) r = np.sqrt(xx**2 + yy**2) data = np.exp(-r**2) # Indicating rank info with a cell array proc = np.ones((x.size-1, y.size-1)) * rank with PRectilinearGrid('pgrid.pvtr', (x, y), offsets[rank]) as rect: rect.addPointData(DataArray(data, range(2), 'gaussian')) rect.addCellData(DataArray(proc, range(2), 'proc')) ``` As you can see, using `PRectilinearGrid` feels just like using `RectilinearGrid`, except that you need to supply the position of the local grid in the global grid numbering (the `offsets[rank]` in the above example). Note that RecilinearGrid VTK files need an overlap in point data, hence why the global grid size ends up being `(2*N-1, 2*N-1)`. If you forget that overlap, Paraview (or another VTK-based software) may complain that some parts in the global grid (aka "extents" in VTK) are missing data. ## List of features Here is a list of what is available in UVW: ### VTK file formats - Image data (`.vti`) - Rectilinear grid (`.vtr`) - Structured grid (`.vts`) - Parallel Rectilinear grid (`.pvtr`) ### Data representation - ASCII - Base64 (raw and compressed: the `compression` argument of file constructors can be `True`, `False`, or an integer in `[-1, 9]` for compression levels) ### Planned developments Here is a list of future developments: - [x] Image data - [ ] Unstructured grid - [x] Structured grid - [x] Parallel writing (`mpi4py`-enabled `PRectilinearGrid` *is now available!*) - [ ] Benchmarking + performance comparison with [pyevtk](https://bitbucket.org/pauloh/pyevtk) ## Developing These instructions will get you a copy of the project up and running on your local machine for development and testing purposes. ### Git repository First clone the git repository: ``` git clone https://github.com/prs513rosewood/uvw.git ``` Then you can use pip in development mode (possibly in [virtualenv](https://virtualenv.pypa.io/en/stable/)): ``` pip install --user -e .[mpi,tests] ``` ## Running the tests The tests can be run using [pytest](https://docs.pytest.org/en/latest/): ``` -cd tests; mpiexec -n 2 pytest --with-mpi +pytest +# or for tests with mpi +mpiexec -n 2 pytest --with-mpi ``` ## License This project is licensed under the MIT License - see the LICENSE.md file for details. ## Acknowledgments * [@PurpleBooth](https://github.com/PurpleBooth)'s [README-Template](https://gist.github.com/PurpleBooth/109311bb0361f32d87a2) diff --git a/tests/conftest.py b/tests/conftest.py index fabe522..127ce47 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,68 +1,71 @@ import pytest import numpy as np from vtk.util.numpy_support import vtk_to_numpy as v2n @pytest.fixture(params=[1, 2, 3]) def field_data(request): N = 4 dim = request.param if dim == 1: x = np.linspace(0, 1, N*2) coords = [x] r = np.abs(x) if dim == 2: x = np.linspace(0, 1, N*2) y = np.linspace(0, 1, N+2) coords = [x, y] xx, yy = np.meshgrid(x, y, indexing='ij', sparse=True) r = np.sqrt(xx**2 + yy**2) if dim == 3: x = np.linspace(0, 1, N*2) y = np.linspace(0, 1, N+2) z = np.linspace(0, 1, N) coords = [x, y, z] xx, yy, zz = np.meshgrid(x, y, z, indexing='ij', sparse=True) r = np.sqrt(xx**2 + yy**2 + zz**2) e_r = np.zeros([s - 1 for s in r.shape] + [3, 3]) e_r[..., :, :] = np.array([[0, 1, 0], [1, 0, 0], [0, 1, 1]]) - e_r[..., :, :] = np.eye(3) return coords, r, e_r @pytest.fixture(params=[False, True]) def compression_fixture(request): return request @pytest.fixture(params=['ascii', 'binary', 'append']) def format_fixture(request): return request def get_vtk_data(reader, sstream): if type(sstream) is str: reader.SetFileName(sstream) else: reader.SetReadFromInputString(True) reader.SetInputString(sstream.getvalue()) reader.Update() output = reader.GetOutput() return v2n(output.GetPointData().GetArray('point')), \ v2n(output.GetCellData().GetArray('cell')) -def transp(dim): - trans = list(range(dim+2)) - trans[-2], trans[-1] = trans[-1], trans[-2] - return trans - +@pytest.fixture(params=['C', 'F']) +def ordering_fixture(request): + def transp(dim): + trans = list(range(dim+2)) + if request.param == 'C': + trans[-2], trans[-1] = trans[-1], trans[-2] + return trans + request.transp = transp + return request diff --git a/tests/test_parallel.py b/tests/test_parallel.py index 5437ba9..d8fb47b 100644 --- a/tests/test_parallel.py +++ b/tests/test_parallel.py @@ -1,106 +1,123 @@ import numpy as np import pytest import vtk import os from mpi4py import MPI from vtk.util.numpy_support import vtk_to_numpy from vtk import vtkXMLPRectilinearGridReader -from conftest import get_vtk_data, transp +from conftest import get_vtk_data from numpy import all from uvw.parallel import PRectilinearGrid from uvw import DataArray -def test_prectilinear_grid(field_data, compression_fixture, format_fixture): +def clean(f): + try: + os.remove(f.pfilename) + os.remove(f.filename) + except FileNotFoundError: + pass + + +@pytest.mark.mpi_skip +def test_prectilinear_grid(field_data, + compression_fixture, + format_fixture, + ordering_fixture): coords, r, e_r = field_data dim = r.ndim out_name = 'test_prectilinear_grid.pvtr' compress = compression_fixture.param format = format_fixture.param - with PRectilinearGrid(out_name, - coords, dim * [0], compression=compress) as rect: - rect.init_master(None) # useless here: for coverage only - rect.addPointData(DataArray(r, range(dim), 'point'), vtk_format=format) - rect.addCellData(DataArray(e_r, range(dim), 'cell'), vtk_format=format) + rect = PRectilinearGrid(out_name, + coords, dim * [0], compression=compress) + rect.init_master(None) # useless here: for coverage only + rect.addPointData( + DataArray( + r, range(dim), 'point', components_order=ordering_fixture.param + ), + vtk_format=format, + ) + rect.addCellData( + DataArray( + e_r, range(dim), 'cell', components_order=ordering_fixture.param + ), + vtk_format=format, + ) + rect.write() reader = vtkXMLPRectilinearGridReader() vtk_r, vtk_e_r = get_vtk_data(reader, out_name) vtk_r = vtk_r.reshape(r.shape, order='F') - vtk_e_r = vtk_e_r.reshape(e_r.shape, order='F').transpose(transp(dim)) + vtk_e_r = vtk_e_r.reshape(e_r.shape, order='F') \ + .transpose(ordering_fixture.transp(dim)) assert all(vtk_r == r) assert all(vtk_e_r == e_r) - os.remove(out_name) - os.remove('test_prectilinear_grid_rank0.vtr') + clean(rect) @pytest.mark.mpi(min_size=2) def test_prectilinear_grid_mpi(compression_fixture, format_fixture): comm = MPI.COMM_WORLD rank = comm.Get_rank() N = 3 bounds = [ (0, 1), (1, 2) ] offsets = [ [0, 0, 0], [2*N, 0, 0], ] x = np.linspace(*bounds[rank], N*2) y = np.linspace(0, 1, N+2) z = np.linspace(0, 1, N) out_name = 'test_prectilinear_grid_mpi.pvtr' xx, yy, zz = np.meshgrid(x, y, z, indexing='ij', sparse=True) r = np.sqrt(xx**2 + yy**2 + zz**2) compress = compression_fixture.param format = format_fixture.param rect = PRectilinearGrid(out_name, (x, y, z), offsets[rank], compression=compress) rect.addPointData(DataArray(r, range(3), 'R'), vtk_format=format) rect.write() if rank == 0: reader = vtk.vtkXMLPRectilinearGridReader() reader.SetFileName(out_name) reader.Update() output = reader.GetOutput() vtk_r = vtk_to_numpy(output.GetPointData().GetArray('R')) else: vtk_r = None vtk_r = comm.bcast(vtk_r, root=0) vtk_r = vtk_r.reshape([4*N-1, N+2, N], order='F') i, j, k = [int(i) for i in offsets[rank]] # Adjusting for overlap if offsets[rank][0] != 0: i -= 1 sub_vtk = vtk_r[i:i+x.size, j:j+y.size, k:k+z.size] assert np.all(sub_vtk == r) - if rank == 0: - try: - os.remove(out_name) - os.remove('test_prectilinear_grid_mpi_rank0.vtr') - os.remove('test_prectilinear_grid_mpi_rank1.vtr') - except FileNotFoundError: - pass + clean(rect) if __name__ == '__main__': test_prectilinear_grid() diff --git a/tests/test_vtk_files.py b/tests/test_vtk_files.py index 61e0912..3e7526a 100644 --- a/tests/test_vtk_files.py +++ b/tests/test_vtk_files.py @@ -1,107 +1,127 @@ import io import numpy as np from numpy import all, min, max from vtk import ( vtkXMLRectilinearGridReader, vtkXMLImageDataReader, vtkXMLStructuredGridReader, ) from vtk.util.numpy_support import vtk_to_numpy -from conftest import get_vtk_data, transp +from conftest import get_vtk_data from uvw import ( ImageData, RectilinearGrid, StructuredGrid, DataArray, ) -def test_rectilinear_grid(field_data, compression_fixture, format_fixture): +def test_rectilinear_grid(field_data, + compression_fixture, + format_fixture, + ordering_fixture): coords, r, e_r = field_data dim = r.ndim f = io.StringIO() compress = compression_fixture.param format = format_fixture.param rect = RectilinearGrid(f, coords, compression=compress) - rect.addPointData(DataArray(r, range(dim), 'point'), vtk_format=format) - rect.addCellData(DataArray(e_r, range(dim), 'cell'), vtk_format=format) + rect.addPointData( + DataArray(r, range(dim), 'point', ordering_fixture.param), + vtk_format=format + ) + rect.addCellData( + DataArray(e_r, range(dim), 'cell', ordering_fixture.param), + vtk_format=format + ) rect.write() reader = vtkXMLRectilinearGridReader() # Testing the xml pretty print output as well pretty_sstream = io.StringIO(str(rect.writer)) for ss in [f, pretty_sstream]: vtk_r, vtk_e_r = get_vtk_data(reader, ss) vtk_r = vtk_r.reshape(r.shape, order='F') - vtk_e_r = vtk_e_r.reshape(e_r.shape, order='F').transpose(transp(dim)) + vtk_e_r = vtk_e_r.reshape(e_r.shape, order='F') \ + .transpose(ordering_fixture.transp(dim)) assert all(vtk_r == r) assert all(vtk_e_r == e_r) -def test_image_data(field_data, compression_fixture, format_fixture): +def test_image_data(field_data, + compression_fixture, + format_fixture, + ordering_fixture): coords, r, e_r = field_data dim = r.ndim f = io.StringIO() compress = compression_fixture.param format = format_fixture.param with ImageData( f, [(min(x), max(x)) for x in coords], [x.size for x in coords], compression=compress) as fh: - fh.addPointData(DataArray(r, range(dim), 'point'), vtk_format=format) - fh.addCellData(DataArray(e_r, range(dim), 'cell'), vtk_format=format) + fh.addPointData( + DataArray(r, range(dim), 'point', ordering_fixture.param), + vtk_format=format + ) + fh.addCellData( + DataArray(e_r, range(dim), 'cell', ordering_fixture.param), + vtk_format=format + ) reader = vtkXMLImageDataReader() vtk_r, vtk_e_r = get_vtk_data(reader, f) vtk_r = vtk_r.reshape(r.shape, order='F') - vtk_e_r = vtk_e_r.reshape(e_r.shape, order='F').transpose(transp(dim)) + vtk_e_r = vtk_e_r.reshape(e_r.shape, order='F') \ + .transpose(ordering_fixture.transp(dim)) assert all(vtk_r == r) assert all(vtk_e_r == e_r) def test_structured_grid(compression_fixture, format_fixture): f = io.StringIO() N = 5 r = np.linspace(0, 1, N) theta = np.linspace(0, 2*np.pi, 5*N) theta, r = np.meshgrid(theta, r, indexing='ij') x = r*np.cos(theta) y = r*np.sin(theta) points = np.vstack([x.ravel(), y.ravel()]).T compress = compression_fixture.param format = format_fixture.param grid = StructuredGrid(f, points, (N, 5*N), compression=compress) data = np.exp(-4*r**2) grid.addPointData(DataArray(data, reversed(range(2)), 'data'), vtk_format=format) grid.write() reader = vtkXMLStructuredGridReader() reader.SetReadFromInputString(True) reader.SetInputString(f.getvalue()) reader.Update() output = reader.GetOutput() vtk_data = vtk_to_numpy(output.GetPointData().GetArray('data')) vtk_data = vtk_data.reshape(data.shape, order='C') assert all(vtk_data == data) diff --git a/uvw/data_array.py b/uvw/data_array.py index 67d25ca..5487847 100644 --- a/uvw/data_array.py +++ b/uvw/data_array.py @@ -1,50 +1,50 @@ import functools class DataArray: """Class holding information on ndarray""" def __init__(self, data, spatial_axes, name='', components_order='C'): """ Data array constructor :param data: the numpy array containing the data (possibly a view) :param spatial_axes: a container of ints that indicate which axes of the array correspond to space dimensions (in order) :param name: the name of the data :param components_order: the order of the non-spatial axes of the array """ self.data = data axes = list(range(data.ndim)) spatial_axes = list(spatial_axes) if data.ndim < len(spatial_axes): raise ValueError( 'Dimensions of data smaller than space dimensions') for ax in spatial_axes: axes.remove(ax) nb_components = functools.reduce( lambda x, y: x * data.shape[y], axes, 1) if components_order == 'C': axes.reverse() - elif components_order == 'F': # pragma: no cover + elif components_order == 'F': pass else: raise ValueError('Unrecognized components order') axes += spatial_axes # Hopefully this is a view self.flat_data = self.data.transpose(*axes).reshape(-1, order='F') self.attributes = { "Name": name, "type": str(self.flat_data.dtype).capitalize(), "NumberOfComponents": str(nb_components) } def __str__(self): # pragma: no cover return self.attributes.__str__()