feat: Add function for creating 3D grids (#36)

Adds `form3DGrid` function under the new `gridtools3d` module.
 - Currently supported grids are:
     - Pure Tets
     - Pure Hexs
     - Mixed Hex/Tet/Pyr meshes
 - Also adds example for mixed Hex/Tet/Pyr meshes in `examples/mixedMesh3D`

Author: jrwrigh

example/mixedMesh3D/makeVTK.py

@@ -0,0 +1,31 @@
+#!/usr/bin/env python
+import numpy as np
+import vtk
+from pathlib import Path
+import pyvista as pv
+
+## ---- Set path to this git repository ----
+# (so that you can use a "custom" version easily)
+import sys
+sys.path.insert(0, '../../')
+
+import vtkpytools as vpt
+
+#%% ---- File inputs ----
+    # Existing files
+coordsPath = Path('./mixedMesh3D.crd')
+connecPath = Path('./mixedMesh3D.cnn')
+
+vtkPath = Path('result/mixedmesh3d.vtk')
+
+#%% ---- Building the VTK Grid object ----
+coords = np.loadtxt(coordsPath, skiprows=1)
+coords = coords[:,:]
+
+connec = np.loadtxt(connecPath, dtype=np.int64, skiprows=1)
+connec = connec[:,:] -1
+
+    # Create the 2D grid from the given information
+grid = vpt.form3DGrid(coords, connectivity_array=connec)
+grid.save(vtkPath)
+

example/mixedMesh3D/mixedMesh3D.cnn

@@ -0,0 +1,4 @@
+  3
+1 2 3  4 9 6 7 8
+4 3 7  8 5 5 5 5
+7 5 10 3 3 3 3 3

example/mixedMesh3D/mixedMesh3D.crd

@@ -0,0 +1,13 @@
+  10
+0.0  0.0  0.0
+1.0  0.0  0.0
+1.0  1.0  0.0
+0.0  1.0  0.0
+0.5  2.0 -0.5
+1.0  0.0 -1.0
+1.0  1.0 -1.0
+0.0  1.0 -1.0
+0.0  0.0 -1.0
+1.25  1.75 -0.5
+
+

example/mixedMesh3D/result/mixedmesh3d.vtk

@@ -1,4 +1,5 @@
 from .barfiletools import *
+from .gridtools3d import *
 from .common import (getGeometricSeries, unstructuredToPoly, orderPolyDataLine,
                      vCutter, Profile, globFile, rotateTensor, makeRotationTensor,
                      symmetric2FullTensor, full2SymmetricTensor, calcStrainRate)

vtkpytools/__init__.py

@@ -0,0 +1 @@
+from .core import form3DGrid

vtkpytools/gridtools3d/__init__.py

@@ -0,0 +1,60 @@
+import numpy as np
+import vtk
+from pathlib import Path
+import pyvista as pv
+
+
+def form3DGrid(coords_array, connectivity_array) -> pv.UnstructuredGrid:
+    """Create 3D VTK UnstructuredGrid from coordinates and connectivity
+
+    Currently supports tet grids, mixed tet/pyramid/hex grids.
+
+    Parameters
+    ----------
+    coords_array : numpy.ndarray
+        Coordinates of the points. Shape in (nPoints,2).
+    connectivity_array : numpy.ndarray
+        Connectivity array of the cells. Shape is (nCells,PointsPerCell). The
+        index of points should start from 0. The type of cell will be inferred
+        based on the number of points per cell. Mixed meshes inferred by
+        repeated node IDs in a single element.
+
+    Returns
+    -------
+    pv.UnstructuredGrid
+        Pyvista UnstructuredGrid object with the grid loaded.
+    """
+    nCells = connectivity_array.shape[0]
+    nPnts = connectivity_array.shape[1]
+    cell_type = None
+    if nPnts == 4:
+        cell_type = vtk.VTK_TETRA # ==int(10)
+    elif nPnts == 8:
+        tetCells = connectivity_array[:,3] == connectivity_array[:,4]
+        pyrCells = (connectivity_array[:,4] == connectivity_array[:,5]) & ~tetCells
+        if not np.any(pyrCells + tetCells): # not a mixed mesh
+            cell_type = vtk.VTK_HEXAHEDRON # ==int(12)
+    else:
+        raise ValueError('This connectivity file has the wrong number of points.'
+                            ' Must be either 4 or 8 points per cell, this has {}'.format(nPnts))
+
+    if cell_type:
+        connectivity_array = np.hstack((np.ones((nCells,1), dtype=np.int64)*nPnts, connectivity_array))
+        offsets = np.arange(0, connectivity_array.size+1, nPnts+1, dtype=np.int64)
+        cell_types = np.ones(nCells, dtype=np.int64) * cell_type
+    else: # mixed mesh
+        hexCells = np.invert(tetCells + pyrCells)
+        cell_types = tetCells*vtk.VTK_TETRA + pyrCells*vtk.VTK_PYRAMID + hexCells*vtk.VTK_HEXAHEDRON
+        offset     = tetCells*np.int64(4)   + pyrCells*np.int64(5) + hexCells*np.int64(8)
+
+        connectivity_array[tetCells, 4:] = np.int64(-1)
+        connectivity_array[pyrCells, 5:] = np.int64(-1)
+        connectivity_array = np.hstack((offset[:,None], connectivity_array))
+        connectivity_array = connectivity_array[connectivity_array != -1]
+
+        offsets = np.roll(np.cumsum(offset+1), 1)
+        offsets[0] = 0
+
+    grid = pv.UnstructuredGrid(offsets, connectivity_array, cell_types, coords_array)
+
+    return grid