From 4f664364127b1adb5b8d97f15c84694a5cbc1fae Mon Sep 17 00:00:00 2001 From: Sigurd Pettersen Date: Tue, 27 Sep 2022 10:55:42 +0200 Subject: [PATCH] ENH: Add grid geometry extraction for VTK's vtkExplicitStructuredGrid (#796) Adds method Grid.get_vtk_esg_geometry_data() that builds and returns grid geometry data in a format optimized for use with VTK's explicit structured grid class. --- .../xtg/grdcp3d_get_vtk_esg_geometry_data.c | 290 ++++++++++++++++++ src/clib/xtg/libxtg.h | 15 + src/xtgeo/grid3d/_grid_etc1.py | 42 +++ src/xtgeo/grid3d/grid.py | 46 +++ tests/test_grid3d/test_grid.py | 49 +++ 5 files changed, 442 insertions(+) create mode 100644 src/clib/xtg/grdcp3d_get_vtk_esg_geometry_data.c diff --git a/src/clib/xtg/grdcp3d_get_vtk_esg_geometry_data.c b/src/clib/xtg/grdcp3d_get_vtk_esg_geometry_data.c new file mode 100644 index 000000000..55d182fcd --- /dev/null +++ b/src/clib/xtg/grdcp3d_get_vtk_esg_geometry_data.c @@ -0,0 +1,290 @@ +/* +**************************************************************************************** +* +* NAME: +* grdcp3d_get_vtk_esg_geometry_data.c +* +* DESCRIPTION: +* Get geometry data that is suitable as input to VTK's vtkExplicitStructuredGrid. +* See: https://vtk.org/doc/nightly/html/classvtkExplicitStructuredGrid.html#details +* +* Basically what we have to do is build an unstructured grid representation where +* all grid cells are represented as hexahedrons with explicit connectivities. +* The connectivity table refers into the accompanying vertex table. +* +* In VTK, cells order increases in I fastest, then J, then K. +* +* This function also tries to remove/weld duplicate vertices, but this is a WIP. +* +* ARGUMENTS: +* ncol, nrow, nlay i Dimensions +* coordsv i Coordinate vector (xtgeo fmt) +* zcornsv i ZCORN vector (xtgeo fmt) +* vert_arr o Will be populated with vertex coordinates (XYZ, XYZ, etc) +* conn_arr o Will be populated with the cell connectivities +* +* RETURNS: +* The actual number of XYZ vertices written into vert_arr. +* Also updates output arrays vert_arr and conn_arr. +* +* TODO/ISSUES/BUGS: +* - verify/enhance removal of duplicate nodes +* - optimize the way cell corners are accessed +* - reduce memory usage in GeoMaker helper class +* +* LICENCE: +* CF XTGeo's LICENSE +*************************************************************************************** +*/ + +#include "libxtg.h" +#include "libxtg_.h" + + +// Forward declarations of static helpers +static struct GeoMaker* gm_create(long cellCountI, long cellCountJ, long cellCountK, double* targetVertexBuffer, long* targetConnBuffer); +static void gm_destroy(struct GeoMaker* gm); +static void gm_addVertex(struct GeoMaker* gm, long i, long j, long k, double v[3]); +static long gm_vertexCount(const struct GeoMaker* gm); +static long gm_connectivityCount(const struct GeoMaker* gm); + + +long grdcp3d_get_vtk_esg_geometry_data(long ncol, + long nrow, + long nlay, + + double *coordsv, + long ncoordin, + float *zcornsv, + long nlaycornin, + + double *vert_arr, + long n_vert_arr, + long *conn_arr, + long n_conn_arr) +{ + const long cellCount = ncol*nrow*nlay; + const long maxVertexCount = 8*cellCount; + const long expectedConnCount = 8*cellCount; + + if (3*maxVertexCount > n_vert_arr) + { + throw_exception("Allocated size of vertex array is too small"); + return 0; + } + if (expectedConnCount != n_conn_arr) + { + throw_exception("Allocated size of connectivity array is too small"); + return 0; + } + + // 7---------6 6---------7 + // /| /| /| /| + // / | / | / | / | + // 4---------5 | 4---------5 | + // | 3------|--2 | 2------|--3 + // | / | / | / | / + // |/ |/ |/ |/ + // 0---------1 0---------1 + // VTK XTGeo + // + // VTK hex cell definition from: + // https://kitware.github.io/vtk-examples/site/VTKBook/05Chapter5/#54-cell-types + + struct GeoMaker* gm = gm_create(ncol, nrow, nlay, vert_arr, conn_arr); + double crs[24]; + for (long k = 0; k < nlay; k++) + { + for (long j = 0; j < nrow; j++) + { + for (long i = 0; i < ncol; i++) + { + grdcp3d_corners(i, j, k, ncol, nrow, nlay, coordsv, 0, zcornsv, 0, crs); + + gm_addVertex(gm, i, j, k, &crs[0*3]); + gm_addVertex(gm, i+1, j, k, &crs[1*3]); + gm_addVertex(gm, i+1, j+1, k, &crs[3*3]); + gm_addVertex(gm, i, j+1, k, &crs[2*3]); + gm_addVertex(gm, i, j, k+1, &crs[4*3]); + gm_addVertex(gm, i+1, j, k+1, &crs[5*3]); + gm_addVertex(gm, i+1, j+1, k+1, &crs[7*3]); + gm_addVertex(gm, i, j+1, k+1, &crs[6*3]); + } + } + } + + const long finalVertexCount = gm_vertexCount(gm); + const long finalConnCount = gm_connectivityCount(gm); + + gm_destroy(gm); + gm = NULL; + + return finalVertexCount; +} + + +// General helpers +// ==================================================================================== + +// ------------------------------------------------------------------------------------ +static long coordsEqual(double a[3], double b[3]) +{ + return (a[0] == b[0] && a[1] == b[1] && a[2] == b[2]); +} + + +// The Geomaker "class" +// ==================================================================================== +typedef struct GeoMaker +{ + long vertexCountI; + long vertexCountJ; + long vertexCountK; + + long* globToOutputVertIdx; + + double* vertexArr; + long numVerticesAdded; + long* connArr; + long numConnAdded; +} GeoMaker; + + +// ------------------------------------------------------------------------------------ +static GeoMaker* gm_create(long cellCountI, long cellCountJ, long cellCountK, double* targetVertexBuffer, long* targetConnBuffer) +{ + GeoMaker* gm = (GeoMaker*)malloc(sizeof(GeoMaker)); + + gm->vertexCountI = cellCountI + 1; + gm->vertexCountJ = cellCountJ + 1; + gm->vertexCountK = cellCountK + 1; + + const long maxGlobVertexIdx = gm->vertexCountI*gm->vertexCountJ*gm->vertexCountK; + gm->globToOutputVertIdx = malloc(8*maxGlobVertexIdx*sizeof(long)); + for (int i = 0; i < 8*maxGlobVertexIdx; i++) + { + gm->globToOutputVertIdx[i] = -1; + } + + gm->vertexArr = targetVertexBuffer; + gm->numVerticesAdded = 0; + gm->connArr = targetConnBuffer; + gm->numConnAdded = 0; + + return gm; +} + +// ------------------------------------------------------------------------------------ +static void gm_destroy(GeoMaker* gm) +{ + if (gm) + { + free(gm->globToOutputVertIdx); + free(gm); + } +} + +// ------------------------------------------------------------------------------------ +static long gm_vertexIJKToIdx(const GeoMaker* gm, long i, long j, long k) +{ + return k*(gm->vertexCountI*gm->vertexCountJ) + j*gm->vertexCountI + i; +} + +// ------------------------------------------------------------------------------------ +static long gm_findVertex(const GeoMaker* gm, long i, long j, long k, double v[3]) +{ + const long globalVertexIdx = gm_vertexIJKToIdx(gm, i, j, k); + + for (int offset = 0; offset < 8; offset++) + { + const long mappedIdx = gm->globToOutputVertIdx[8*globalVertexIdx + offset]; + if (mappedIdx < 0) + { + return -1; + } + + if (coordsEqual(v, &gm->vertexArr[3*mappedIdx])) + { + return mappedIdx; + } + } + + return -1; +} + +// ------------------------------------------------------------------------------------ +static long gm_findVertexAlongIJ(const GeoMaker* gm, long i, long j, long k, double v[3]) +{ + if (k > 0) + { + const long mappedIdx = gm_findVertex(gm, i, j, k - 1, v); + if (mappedIdx >= 0) + { + return mappedIdx; + } + } + + if (k < gm->vertexCountK - 1) + { + const long mappedIdx = gm_findVertex(gm, i, j, k + 1, v); + if (mappedIdx >= 0) + { + return mappedIdx; + } + } + + return -1; +} + +// ------------------------------------------------------------------------------------ +static void gm_addVertex(GeoMaker* gm, long i, long j, long k, double v[3]) +{ + const long globalVertexIdx = gm_vertexIJKToIdx(gm, i, j, k); + + long mappedIdx = -1; + for (int offset = 0; offset < 8; offset++) + { + mappedIdx = gm->globToOutputVertIdx[8*globalVertexIdx + offset]; + if (mappedIdx < 0) + { + // Take a peek above and below + mappedIdx = gm_findVertexAlongIJ(gm, i, j, k, v); + if (mappedIdx >= 0) + { + gm->globToOutputVertIdx[8*globalVertexIdx + offset] = mappedIdx; + gm->connArr[gm->numConnAdded] = mappedIdx; + gm->numConnAdded++; + return; + } + + mappedIdx = gm->numVerticesAdded; + + memcpy(&gm->vertexArr[3*mappedIdx], v, 3*sizeof(double)); + gm->numVerticesAdded++; + + gm->globToOutputVertIdx[8*globalVertexIdx + offset] = mappedIdx; + + break; + } + + if (coordsEqual(v, &gm->vertexArr[3*mappedIdx])) + { + break; + } + } + + gm->connArr[gm->numConnAdded] = mappedIdx; + gm->numConnAdded++; +} + +// ------------------------------------------------------------------------------------ +static long gm_vertexCount(const GeoMaker* gm) +{ + return gm->numVerticesAdded; +} + +// ------------------------------------------------------------------------------------ +static long gm_connectivityCount(const GeoMaker* gm) +{ + return gm->numConnAdded; +} diff --git a/src/clib/xtg/libxtg.h b/src/clib/xtg/libxtg.h index 76e8af696..549244a8f 100644 --- a/src/clib/xtg/libxtg.h +++ b/src/clib/xtg/libxtg.h @@ -1864,6 +1864,21 @@ grdcp3d_corners(long ic, long n_swig_np_flt_inplaceflat_v1, double corners[]); +long +grdcp3d_get_vtk_esg_geometry_data(long ncol, + long nrow, + long nlay, + + double *swig_np_dbl_inplaceflat_v1, // coordsv + long n_swig_np_dbl_inplaceflat_v1, + float *swig_np_flt_inplaceflat_v1, // zcornsv + long n_swig_np_flt_inplaceflat_v1, + + double *swig_np_dbl_aout_v1, // output vertex arr + long n_swig_np_dbl_aout_v1, // allocated length of vertex arr + long *swig_np_lng_aout_v1, // hex connectivity array (out) + long n_swig_np_lng_aout_v1); // allocated length of hex conn arr + void grdcp3d_get_vtk_grid_arrays(long ncol, long nrow, diff --git a/src/xtgeo/grid3d/_grid_etc1.py b/src/xtgeo/grid3d/_grid_etc1.py index ef6cf85b2..7936bdc05 100644 --- a/src/xtgeo/grid3d/_grid_etc1.py +++ b/src/xtgeo/grid3d/_grid_etc1.py @@ -576,6 +576,48 @@ def get_xyz_corners(self, names=("X_UTME", "Y_UTMN", "Z_TVDSS")): return tuple(grid_props) +def get_vtk_esg_geometry_data( + self, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Get geometry data consisting of vertices and cell connectivities suitable for + use with VTK's vtkExplicitStructuredGrid. + + Returned tuple contains: + - numpy array with dimensions in terms of points (not cells) + - vertex array, numpy array with vertex coordinates + - connectivity array for all the cells, numpy array with integer indices + - inactive cell indices, numpy array with integer indices + """ + + self._xtgformat2() + + # Number of elements to allocate in the vertex and connectivity arrays + num_cells = self.ncol * self.nrow * self.nlay + n_vertex_arr = 3 * 8 * num_cells + n_conn_arr = 8 * num_cells + + # Note first value in return tuple which is the actual number of vertices that + # was written into vertex_arr and which we'll use to shrink the array. + vertex_count, vertex_arr, conn_arr = _cxtgeo.grdcp3d_get_vtk_esg_geometry_data( + self.ncol, + self.nrow, + self.nlay, + self._coordsv, + self._zcornsv, + n_vertex_arr, + n_conn_arr, + ) + + # Need to shrink the vertex array + vertex_arr = np.resize(vertex_arr, 3 * vertex_count) + vertex_arr = vertex_arr.reshape(-1, 3) + + point_dims = np.asarray((self.ncol, self.nrow, self.nlay)) + 1 + inact_indices = self.get_actnum_indices(order="F", inverse=True) + + return point_dims, vertex_arr, conn_arr, inact_indices + + def get_vtk_geometries(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """Return actnum, corners and dims arrays for VTK ExplicitStructuredGrid usage.""" self._xtgformat2() diff --git a/src/xtgeo/grid3d/grid.py b/src/xtgeo/grid3d/grid.py index 32ae6651f..effa90956 100644 --- a/src/xtgeo/grid3d/grid.py +++ b/src/xtgeo/grid3d/grid.py @@ -1196,6 +1196,52 @@ def get_dataframe(self, activeonly=True, ijk=True, xyz=True, doubleformat=False) def dataframe(self, *args, **kwargs): return self.get_dataframe(*args, **kwargs) + def get_vtk_esg_geometry_data(self): + """Get grid geometry data suitable for use with VTK's vtkExplicitStructuredGrid. + + Builds and returns grid geometry data in a format tailored for use with VTK's + explicit structured grid (ESG). Essentially this entails building an + unstructured grid representation where all the grid cells are represented as + hexahedrons with explicit connectivities. The cell connectivity array refers + into the accompanying vertex array. + + In VTK, cell order increases in I fastest, then J, then K. + + The returned tuple contains: + - numpy array with dimensions in terms of points (not cells) + - vertex array, numpy array with vertex coordinates + - connectivity array for all the cells, numpy array with integer indices + - inactive cell indices, numpy array with integer indices + + This function also tries to remove/weld duplicate vertices, but this is still + a work in progress. + + Example usage with VTK:: + + dims, vert_arr, conn_arr, inact_arr = xtg_grid.get_vtk_esg_geometry_data() + + vert_arr = vert_arr.reshape(-1, 3) + vtk_points = vtkPoints() + vtk_points.SetData(numpy_to_vtk(vert_arr, deep=1)) + + vtk_cell_array = vtkCellArray() + vtk_cell_array.SetData(8, numpy_to_vtkIdTypeArray(conn_arr, deep=1)) + + vtk_esgrid = vtkExplicitStructuredGrid() + vtk_esgrid.SetDimensions(dims) + vtk_esgrid.SetPoints(vtk_points) + vtk_esgrid.SetCells(vtk_cell_array) + + vtk_esgrid.ComputeFacesConnectivityFlagsArray() + + ghost_arr_vtk = vtk_esgrid.AllocateCellGhostArray() + ghost_arr_np = vtk_to_numpy(ghost_arr_vtk) + ghost_arr_np[inact_arr] = vtkDataSetAttributes.HIDDENCELL + + .. versionadded:: 2.20 + """ + return _grid_etc1.get_vtk_esg_geometry_data(self) + def get_vtk_geometries(self): """Get necessary arrays on correct layout for VTK ExplicitStructuredGrid usage. diff --git a/tests/test_grid3d/test_grid.py b/tests/test_grid3d/test_grid.py index 8701a6e5f..6128bda11 100644 --- a/tests/test_grid3d/test_grid.py +++ b/tests/test_grid3d/test_grid.py @@ -798,3 +798,52 @@ def test_get_vtk_geometries_emerald(show_plot): plt.set_scale(1, 1, 5) plt.show_axes() plt.show() + + +def test_get_vtk_esg_geometry_data_four_cells(): + """Test that extracted VTK ESG connectivity and vertex coordinates are correct + for simple grid with only four cells in a single layer""" + grd = xtgeo.create_box_grid((2, 2, 1), increment=(1, 3, 10)) + + # Note that returned dims is in terms of points + dims, vertex_arr, conn_arr, inactive_cell_indices = grd.get_vtk_esg_geometry_data() + + assert list(dims) == [3, 3, 2] + assert len(vertex_arr) == 3 * 3 * 2 + assert len(conn_arr) == 8 * 4 + assert len(inactive_cell_indices) == 0 + + # Bounding box of first cell should be min/max x=[0,1] y=[0,3], z[0,10] + assert vertex_arr[conn_arr[0]] == pytest.approx([0, 0, 0]) + assert vertex_arr[conn_arr[1]] == pytest.approx([1, 0, 0]) + assert vertex_arr[conn_arr[2]] == pytest.approx([1, 3, 0]) + assert vertex_arr[conn_arr[3]] == pytest.approx([0, 3, 0]) + assert vertex_arr[conn_arr[4]] == pytest.approx([0, 0, 10]) + assert vertex_arr[conn_arr[5]] == pytest.approx([1, 0, 10]) + assert vertex_arr[conn_arr[6]] == pytest.approx([1, 3, 10]) + assert vertex_arr[conn_arr[7]] == pytest.approx([0, 3, 10]) + + # Bounding box of fourth cell should be min/max x=[1,2] y=[3,6], z[0,10] + assert vertex_arr[conn_arr[24]] == pytest.approx([1, 3, 0]) + assert vertex_arr[conn_arr[25]] == pytest.approx([2, 3, 0]) + assert vertex_arr[conn_arr[26]] == pytest.approx([2, 6, 0]) + assert vertex_arr[conn_arr[27]] == pytest.approx([1, 6, 0]) + assert vertex_arr[conn_arr[28]] == pytest.approx([1, 3, 10]) + assert vertex_arr[conn_arr[29]] == pytest.approx([2, 3, 10]) + assert vertex_arr[conn_arr[30]] == pytest.approx([2, 6, 10]) + assert vertex_arr[conn_arr[31]] == pytest.approx([1, 6, 10]) + + +def test_get_vtk_esg_geometry_data_box(): + """Test that extracted VTK ESG geometry looks sensible for small test grid with + incative cells""" + grd = xtgeo.create_box_grid((2, 6, 4), increment=(20, 10, 7)) + grd._actnumsv[0, 0, 0] = 0 + grd._actnumsv[1, 0, 0] = 0 + + dims, vertex_arr, conn_arr, inactive_cell_indices = grd.get_vtk_esg_geometry_data() + + assert list(dims) == [3, 7, 5] + assert len(vertex_arr) == 3 * 7 * 5 + assert len(conn_arr) == 8 * 2 * 6 * 4 + assert len(inactive_cell_indices) == 2