HexaShrink, an exact scalable framework for hexahedral meshes with attributes and discontinuities: multiresolution rendering and storage of geoscience models
Jean-Luc Peyrot, Laurent Duval, Frédéric Payan, Lauriane Bouard, Lénaïc Chizat, Sébastien Schneider, Marc Antonini
CComputational Geosciences
HexaShrink, an exact scalable framework for hexahedralmeshes with attributes and discontinuities: multiresolutionrendering and storage of geoscience models
Jean-Luc Peyrot · Laurent Duval · Fr´ed´eric Payan · Lauriane Bouard · L´ena¨ıc Chizat · S´ebastien Schneider · Marc Antonini
Last compiled: Tuesday 7 th May, 2019, 01:18
Abstract
With huge data acquisition progresses real-ized in the past decades and acquisition systems nowable to produce high resolution grids and point clouds,the digitization of physical terrains becomes increas-ingly more precise. Such extreme quantities of gener-ated and modeled data greatly impact computationalperformances on many levels of high-performance com-puting (HPC): storage media, memory requirements,transfer capability, and finally simulation interactiv-ity, necessary to exploit this instance of big data. Ef-ficient representations and storage are thus becoming
This work was partly presented in [43]. This is a post-peer-review, pre-copyedit version of an article published in Compu-tational Geosciences. The final authenticated version is avail-able online at: http://dx.doi.org/10.1007/s10596-019-9816-2Laurent Duval and L´ena¨ıc ChizatIFP Energies nouvelles1 et 4 avenue de Bois-Pr´eauF-92852 Rueil-MalmaisonTel.: +33 (0) 1 47 52 61 02E-mail: [email protected],[email protected] DuvalUniversity Paris-Est, LIGMESIEE ParisF-93162 Noisy-le-GrandJean-Luc Peyrot, S´ebastien Schneider and Lauriane BouardIFP Energies nouvellesRond-Point de l’´echangeur de Solaize, BP3F-69360 SolaizeE-mail: [email protected],[email protected]´ed´eric Payan and Marc AntoniniUniversit´e Cˆote d’Azur, CNRS, I3S2000, route des LuciolesF-06900 Sophia AntipolisTel.: +33(0) 4 89 15 43 22 / +33(0) 4 92 94 27 18E-mail: [email protected],[email protected]. Chizat is now with INRIA, ENS, PSL Research UniversityParis, France. S´ebastien Schneider is now with HoloMake. “enabling technologies” in HPC experimental and simu-lation science [17]. We propose HexaShrink, an originaldecomposition scheme for structured hexahedral vol-ume meshes. The latter are used for instance in biomed-ical engineering, materials science, or geosciences. Hexa-Shrink provides a comprehensive framework allowingefficient mesh visualization and storage. Its exactly re-versible multiresolution decomposition yields a hierar-chy of meshes of increasing levels of details, in terms ofeither geometry, continuous or categorical properties ofcells.Starting with an overview of volume meshes com-pression techniques, our contribution blends coherentlydifferent multiresolution wavelet schemes in different di-mensions. It results in a global framework preservingdiscontinuities (faults) across scales, implemented as afully reversible upscaling at different resolutions. Ex-perimental results are provided on meshes of varyingsize and complexity. They emphasize the consistencyof the proposed representation, in terms of visualiza-tion, attribute downsampling and distribution at dif-ferent resolutions. Finally, HexaShrink yields gains instorage space when combined to lossless compressiontechniques.
Keywords
Compression · Corner point grid · Discretewavelet transform · Geometrical discontinuities · Hexahedral volume meshes · High-PerformanceComputing · Multiscale methods · Simulation · Upscaling
Simulation sciences and scientific modelling in high-performance computing employ meshes with increas-ing precision and dynamics. Among them, hexahedral a r X i v : . [ c s . G R ] M a y Peyrot, Duval, Payan, Bouard, Chizat, Schneider, Antonini meshes are commonly handled in biomedical engineer-ing [28], computational materials science [39] and ingeosciences. They are for instance used by geologists tostudy flow simulations for reservoir modelling [7], andbenefit from an increasing interest for geologic modelbuilding [8,34]. Huge progresses in data acquisition pro-duce increasingly more accurate digitization of physicalterrains. The tremendous quantity of data thus gener-ated prominently impacts computational resources andperformances: memory size required for their storageand visualization, but also their transmission and trans-fer, and ultimately their processing. Consequently, itgreatly affects the overall simulation interactivity. Thistrend affects the oil and gas sector at large [41].We propose HexaShrink, an efficient multiscale rep-resentation dedicated to hexahedral meshes with at-tributes and discontinuities. HexaShrink combines fourwavelet-like decompositions to adapt to the heteroge-neous nature of geoscience meshes. Geometrical, contin-uous and categorical properties are consistently down-sampled (upscaled in geoscience terms [7, pp. 181–204,Simulation and Upscaling]) in an exactly reversible man-ner. In addition to regular structures, HexaShrink alsostrives to manage mesh externalities like boundariesand borders tagged as inactive cells for simulation pur-poses. It produces a hierarchy of meshes at dyadic res-olutions maintaining geometrical coherency over scales,consistent with geomodeler/simulator upscaling opera-tions. It finally lends itself to efficient lossless storage,in combination with state-of-the-art compression algo-rithms.The paper is structured as follows: Section 2 presentsthe specificities of structured hexahedral meshes andthe discontinuities they may contain. Section 3 reviewsprior methods for volume mesh representation or com-pression [16]. We introduce the HexaShrink structuredmesh representation in Section 4. We detail the fourmain multiscale wavelet-like schemes that entail an ex-actly invertible hierarchy of downsampled meshes, con-sistently with respect to geometry and properties. Aspecial care is taken on the accurate representation offaults and the management of mesh borders. Section5 presents visual results, and evaluates the quality ofthe HexaShrink with respect to categorical property co-herency across scales, and dyadic upscaling by a geo-modeler. An exhaustive evaluation obtained from com-bining HexaShrink with different lossless compressionalgorithms, at different resolution levels, shows the in-terest of the proposed representation in terms of losslesscompression for storage. Finally, Section 6 summarizesour contributions and proposes future works. D objects. They partition their innerspace with a set of three-dimensional elements namedcells (or zones). While pyramid and triangular prismpartitions exist, most of the existing VMs are composedof tetrahedral (4 faces) or hexahedral (6 faces) elements.They are called tets or hexes (sometimes bricks), re-spectively. A VM composed of different kinds of cells,tetrahedra and hexahedra for instance, is termed hy-brid. A VM is described by the location of vertices in3 D space ( geometry ) and the incidence information be-tween cells, edges, and vertices ( connectivity ). In func-tion of the application domain, VMs also contain phys-ical properties associated to vertices, edges, or cells. Ingeosciences, properties can be scalar (like a single poros-ity value, Figure 1) or vectorial (a vector of composi-tional proportions in different rocks within a cell. . . ).We also distinguish categorical or nominal variables(symbolic and discrete values describing the composi-tion of rocks: sandstone (0), limestone (1), shale (2). . . )from continuous properties (saturation, porosity, per-meability, or a temperature taking values in a givenrange [ − T (cid:91) , T (cid:93) ]).Fig. 1: Example of a VM used in geosciences (left); samemesh with the associated porosity property (right).Note that this mesh is hybrid and unstructured, withboth hexahedral and tetrahedral elements.A non-degenerate hex has 6 faces named quads, 12 edges , and 8 vertices . Depending on incidence infor-mation between cells, edges and vertices, hexahedralmeshes are either unstructured or structured. The de-gree of an edge is the number of adjacent faces. Anhexahedral mesh is unstructured if cells are placed ir-regularly in the volume domain, i.e. , if degrees are notthe same for all edges of the same nature. Unstruc-tured meshes have an important memory footprint, asall the connectivity information must be described ex-plicitly. However, they are well-suited to model complexvolumes, Computer-aided design (CAD) models for in-stance, as shown in Figure 2.An hexahedral mesh is structured if cells are regu-larly organized in the volume domain, i.e. , if the degree exaShrink: multiresolution reversible compression of hexahedral meshes for geoscience simulation models 3 Fig. 2: CAD model defined by an unstructured VM.is equal to four for interior edges (inside the volume),and equal to two for border edges (on a border of thevolume). In that case, the set of hexahedral cells is topo-logically aligned on a 3 D Cartesian grid (see Figure 3).Each vertex of the mesh can be associated to a nodeof the grid. Hence, each cell can be indexed by onlyone triplet ( i, j, k ), and the connectivity informationbecomes implicit: only the position of the vertices isneeded to model the mesh.
Cell ( elementary volume ) ( hexahedron , t e trah e dr on ) Cell face
Edge N ode : grid vertices Gri d Fig. 3: Structured hexahedral mesh composed of(5 × ×
3) cells.2.2 Hexahedral meshes & geometrical discontinuitiesHexahedral meshes in geosciences are generally struc-tured, and thus based on a Cartesian grid. But thesemeshes may contain geometrical discontinuities . Theycorrespond, for instance, to geological faults. It inducesa vertex disparity in space at the same node. The as-sociation of one node of the Cartesian grid with 8 ver-tices (one for each adjacent cell) handles this specificity.Figure 4-(a) provides an illustration of a fault-free vol-ume. On Figure 4-(b), we see that this structure allowsto describe for instance a vertical fault (by positioningvertices differently about the node), while preservingthe Cartesian grid. node with vertex node with equal vertices1 node with equal vertices node with equal vertices (a) Free-fault area. node with 4 vertices at position P and vertices at position P’ (b) Area with a verticalfault. Fig. 4: A fault-free and a fault area. The most popular data structure for structured hex-ahedral meshes with geometrical discontinuities is the
Corner Point Grid [33, 46] tessellation of an Euclidean3 D volume. This structure is often termed pillar grid . Itis based on a set of vertical or inclined pillars runningfrom the top to the bottom of the geological model.A cell is defined by its 8 adjacent vertices (2 on eachadjacent pillar, see Figure 5), and the vertices of theadjacent cells are described independently, in order tomodel faults and gaps. Across the associated Cartesiangrid, each cell can be indexed by a triplet ( i, j, k ). i j k N ode (0,1,1) Hexahedral cell (0,0,0)
Lower extremity of pillar (1,1) P illar (0,1) Gri d (i,j) Fig. 5: An hexahedron, according to the pillar grid structure.This pillar grid also allows to model geological col-lapses (or erosion surfaces), by using degenerate cells,i.e., cells with (at least) two vertices on one pillar lo-cated at the same position (see Figure 6 for differentdegenerate configurations). (a) Single. (b) Two opposed. (c) Two adjacent.
Fig. 6: Degenerate hexahedral cells due to a single col-lapsed pillar (a) or two different collapsed pillar loca-tions (b), (c).
We deal in this section mostly with the ontological de-scription of volume meshes, leaving aside the specifici-ties related to actual data format.
Peyrot, Duval, Payan, Bouard, Chizat, Schneider, Antonini quantization ofvertex coordinates. It consists of constraining the ver-tex coordinates to a discrete and finite set of values.Hence, it becomes possible to encode each coordinatewith an integer index, instead of a 32-bit floating-pointvalue. It is common to quantize the coordinates with12 or 16 bits, reducing the geometry information bya compression factor of 2 .
60 or 2 respectively. Quanti-zation inevitably introduces an irreversible loss in ac-curacy. Visualization typically tolerates precision loss(as long as visual distortion remains negligible), un-like some numerical simulations requiring more precisecomputations.
Prediction (as well as related interpolation meth-ods) further improves the geometry compactness. Pre-dictive coding resorts to estimating the position of avertex from already encoded neighbor vertices.
Pre-diction errors (differences between predicted and ac-tual positions) are generally smaller in amplitude andsparser, which makes their entropy coding (which codesdifferently frequently occurring patterns) efficient [48,p. 63 sq.].Regarding connectivity, when meshes are unstruc-tured, the most frequent technique performs a traversalof mesh elements, and describes the incidence configu-rations with a reduced list of symbols. These symbolsare then entropy coded. When meshes are structured,the connectivity is implicit, reducing its cost to zero.For such meshes, the only additional information toencode are geometrical discontinuities describing faultsand gaps.3.2 Volume mesh compression: prior worksThe basic tools previously presented can be implementedon the ontological structure of meshes, and improved inmany different ways. Their combination, with the as-sistance of advanced compression techniques, permits more efficient tetrahedral or hexahedral mesh coding.Previously proposed algorithms are presented below,classified into two categories: single-rate and progres-sive/multiresolution . They lead to a compact mesh representation, most ofthe time driven by efficient connectivity encoding. Thefirst method for tetrahedral meshes, Grow & Fold, waspresented by Szymczak and Rossignac [51] at the endof the nineties. It is an extension of
EdgeBreaker [47]developed for triangle meshes. The method consists inbuilding a tetrahedral spanning tree from a root tetra-hedron. The traversal is arbitrary among the three neigh-boring tets (Section 2.1) of the cell currently processed,and 3 bits are needed to encode each cell. The resultingspanning tree does not retain the same topology as theoriginal mesh, because some vertices are replicated dur-ing the traversal. “Fold” and “glue” techniques are thusneeded during encoding to restore the original meshfrom the tetrahedron tree. The additional cost is 4 bits,leading to a total cost of 7 bits per tetrahedron.The cut-border initiated in [19] was adapted to tetra-hedral meshes [18]. It denotes the frontier between tetra-hedra already encoded and those to encode. At eachiteration, either a triangle or an adjacent tetrahedronis added to the cut-border. In this case, if the addedvertex is not already in the cut-border, this latter isincluded by a connect operation, and is given a localindex. As the indexing is done locally, the integers to en-code are very small, leading to a compact connectivityrepresentation. In addition, two methods are proposedto encode geometry and associated properties, basedon prediction and entropy coding. This method yieldsgood bit rates (2 .
40 bits per tetrahedron for connec-tivity) for usual meshes, handles non-manifold borders,but worst-cases severely impact bitrates and runtimes(which tend to be quadratic).Isenburg and Alliez [22] are the first to deal withhexahedral VMs. The connectivity is encoded as a se-quence of edge degrees — in a way similar to [52] fortriangular meshes — via a region-growing process of aconvex hull called hull surface . It relies on the assump-tion that hexahedral meshes are often highly regular,which implies that the majority of vertices are sharedby 8 cells. It involves an almost constant edge degree allover the mesh, which significantly decreases the entropyof the connectivity information. A context-based arith-metic coder [53] is then used to encode the connectivityat very low bit rates, between 0 .
18 and 1 .
55 bits per hex-ahedron. Regarding geometry, a user-defined quantiza-tion first restricts the number of bits for coordinates, exaShrink: multiresolution reversible compression of hexahedral meshes for geoscience simulation models 5 and then a predictive scheme based on the parallelo-gram rule encodes the position of vertices added duringthe region-growing process.Krivograd et al. [30] proposed a variant to [22] thatencodes the vertex degrees — number of non-compressedhexahedra around a given vertex — instead of the edgedegrees. On the one hand, this variant achieves bettercompression performances than [22] for dense meshes.On the other hand, it only deals with manifold meshes,and the algorithm is complex as interior cells are en-coded after border cells (it involves many specific casesto process when encoding the connectivity).Lindstrom and Isenburg proposed an original algo-rithm for unstructured meshes called Hexzip [35]. Thisalgorithm is considered as fully lossless, because theinitial indexing of vertices and hexahedra is preserved.For this purpose, connectivity is encoded directly inits indexed structure, by predicting the eight indices ofan hexahedron from preceding ones. This technique issuitable because hexahedral meshes generally have reg-ular strides between indices of subsequent hexahedra. Ahash-table is then used to transform the index structureinto a very redundant and byte-aligned list of symbols,that can be compressed efficiently with gzip (discussedin Section 5.4). Concerning the geometry, spectral pre-diction [21] is used. This algorithm is faster and lessmemory intensive than [22] as the connectivity is notmodified. It handles non-manifold meshes and degen-erate elements. On the other hand, bitrates are higherbecause of the lossless constraints. Unlike methods pre-sented above, Chen et al. [11] focus on geometry com-pression for tetrahedral meshes. The authors proposed a flipping approach based on an extension of the parallel-ogram rule (initially proposed for triangle meshes [52])to tetrahedra. It consists in predicting the position ofan outer vertex of two face-adjacent tetrahedra, withrespect to the other vertices. To globally optimize thegeometry compression, a Minimum Spanning Tree min-imizing the global prediction error for the whole meshis computed. This method is more efficient than priorflipping approaches whose traversal does not depend onthe geometry, but solely on the connectivity.Streaming compression is a subcategory of singlerate compression, dedicated to huge data that cannotfit entirely in the core memory. A particular attention toI/O efficiency is thus required, to enable the encodingof huge meshes with a small memory footprint. Isen-burg and coworkers are the first to propose streamingcompression for VMs (extended from his method for tri-angular meshes [25]): for tetrahedral meshes [24], andthen for hexahedral meshes [14]. In the latter, for in-stance, the compressor does not require the knowledgeof the full list of vertices and cells before encoding. The compressor starts encoding the mesh as soon as the firsthexahedron and its eight adjacent vertices have beenread. For a given hexahedron: i) its face-adjacency isfirst encoded in function of its configuration with hex-ahedra already processed; ii) positions of vertices thatare referenced for the first time are predicted (spectralprediction from adjacent cells); iii) prediction errors areencoded; iv) data structures relative to vertices, becom-ing useless (because their incidence has been entirelydescribed) are finally removed from memory. Comparedto other single rate techniques, streaming tends to achievesimilar compression performances for geometry, but poorerperformances for connectivity. algorithms (also called scalable or multires-olution , see Section 4.1 for details) enable the originalmeshes to be represented and compressed at successiveLODs (levels of details). The main advantage is that itis not necessary to decompress a mesh entirely beforevizualising it. A coarse approximation of the mesh (alsoknown as its lowest resolution) is first decompressedand displayed. Then this coarse mesh is updated withthe successive LOD (termed higher resolutions) that aredecompressed progressively. While they cannot achieveyet compression performance of single-rate algorithms,progressive algorithms are popular because they enableLOD, and also adaptive transmission and displaying,in function of user constraints (network, bandwidth,screen resolution. . . ).Pajarola et al. [40] are the first to propose in 1999a progressive algorithm dedicated to VM compression.This work is inspired by a simplification technique fortetrahedral meshes [49]. It simplifies a given tetrahe-dral mesh progressively, by using successive edge col-lapses [20]. Each time an edge is collapsed, its adjacentcells are removed, and all the information required toreverse this operation is stored: index of the vertex tosplit, and the set of incident faces to “cut”. Thus, duringdecompression, the LODs can be also recovered itera-tively, by using the stored data describing vertex splits .During coding, an edge is selected such as its collapseleads to the minimal error, with respect to specific costfunctions. This algorithm gives a bitrate inferior to 6bits per tetrahedron (for connectivity only).In 2003, Danovaro et al. [15] proposed two progres-sive representations based on a decomposition of a fielddomain into tetrahedral cells. The first is based on ver-tex splits , as the previous method, the second is basedon tetrahedron bisections . This operation consists insubdividing a tetrahedron into two tetrahedra by addinga vertex in the middle of its longest edge. Unlike with
Peyrot, Duval, Payan, Bouard, Chizat, Schneider, Antonini vertex splits, the representation based on tetrahedronbisections is obtained by following a coarse-to-fine ap-proach, i.e. , by applying successive bisections to an ini-tial coarse mesh. Also, this representation only needsto encode the difference vectors between the verticesadded by bisections and theirs real positions. This rep-resentation is thus more compact, as the mesh topol-ogy does not need to be encoded, but it only deals withstructured meshes.VMs multiresolution decomposition based on wavelets[27] was proposed by Boscard´ın et al. [4] for tetrahe-dral meshes. It is based on the tetrahedron subdivi-sion scheme [3] that transforms a tetrahedron into 8sub-tetrahedra, by introducing 6 new vertices on eachedge. After analysis, the input mesh is replaced by abase tetrahedral mesh, and several sets of detail waveletcoefficients. Although coefficients corresponding to dif-ferences between two resolutions could be encoded formesh synthesis, this work does not provide an actualcompression scheme.In [12], Chizat proposed a prototype for a multireso-lution decomposition of geoscientific hexahedral mesheswith the pillar grid structure (Section 2.2). His maincontribution resides in a multiresolution analysis (MRA)that partially manages geometrical discontinuities rep-resenting the faults. It can be achieved by using a mor-phological wavelet transform (Section 4.2.2). This non-separable transform enables the preservation of somefault shapes at different resolutions, as depicted in Fig-ure 7.Fig. 7: Dyadic non-separable multiresolution renderingon a simple geologic mesh [12].The latter work is a seed for the upcoming descrip-tion of HexaShrink. approximation can be interpreted asa decomposition of data at different resolutions, LOD or scales, through a recursive analysis process. It iscalled exact, reversible or invertible when a synthesis scheme can retrieve the original data. Inter-scale rela-tionships [9, 10] often yield sparsification or increased compressibility on sufficiently regular datasets. In dis-crete domains, each analysis stage transforms a set ofvalues (continuous or categorical, in one or several di-mensions), denoted by S . The resulting representationconsists in one subset S − of approximation coefficientsat a lower resolution identified by a negative index S − that approximates the original signal, plus one sub-set of details D − or a combination thereof. The lat-ter represents information missing in the approxima-tion S − . Depending on the MRA scheme, the lower resolution S − may represent a coarsening or “low fre-quencies” of the original samples, or an upscaling ingeosciences (cf. Section 4). The subset D − representsrefinements, fast variations or “high-frequency” detailsremoved from S . We consider here exact systems, al-lowing the perfect recovery of S from a combination ofsubsets S − and D − . Hence, a similar analysis stagecan be applied iteratively, and perfectly again, to thelower resolution S − , in a so-called pyramid scheme.Thus, with the non-positive extremum decompositionlevel L , and indices 0 ≥ l ≥ L , after an | L | -level mul-tiresolution decomposition, the input set S is now de-composed and represented by the subset S L — a (very)coarse approximation of S — and | L | subsets of details D L , . . . , D l , . . . , D − , representing information missingbetween each two consecutive approximations.We consider in the following four different MRA fla-vors, all called wavelets for simplicity. They stem fromiterated, (rounded) linear or non-linear combinations ofcoefficients, as well as separable (applied separately in1 D on each direction) or non-separable ones. Withoutgoing into technicalities here (cf. Section 4.2.3), com-putations are performed using the lifting scheme . Itsuffices to mention that lifting uses complementary in-terleaved grids of values, often indexed with odd andeven indices. Values on one grid are usually predicted(approximations) and updated (details) from the oth-ers. The main interests reside in reduced computationalload, in-place computations and the possibility to main-tain exact integer precision, using for instance only dyadic-rational coefficients (written as m/ n , ( m, n ) ∈ Z × N )and rounding. We refer to [27, sections 2.3, 3.2 and4.3] for a concise account on both non-separable andnon-linear wavelet MRAs, and to [5, 29, 44, 50] for morecomprehensive visions of wavelets and their lifting im-plementations. A recent use in geological model upscal-ing is given in [45].More simply put, for our hexahedral VMs, the dyadicanalysis stage transforms each cell block C l of values exaShrink: multiresolution reversible compression of hexahedral meshes for geoscience simulation models 7 around 2 = 8 contiguous cells (possibly borrowingvalues from a limited cell neighborhood) at resolution l . They are turned into one approximating cell (lowerresolution ( S l − ), and a subset of 2 − D l − , as depicted in Figure 8 along with the re-verse synthesis stage. Hence, if a VM at resolution l Analysis
Original group of (2,2,2) coefficients
Detail coefficients Approximation coefficient
Synthesis
Detail coefficients Approximation coefficient Original group of (2,2,2) coefficients
Fig. 8: Analysis and synthesis stages for VMs.is composed of (cid:0) C li × C lj × C lk (cid:1) cells, C li being the num-ber of cells in direction i , the VM of lower resolutionwill be of dimension (cid:108) C li (cid:109) × (cid:24) C lj (cid:25) × (cid:108) C lk (cid:109) , to take intoaccount non-power-of-two sized grids. As several dig-ital attributes are associated to each cell (geometry,continuous or categorical properties), different types ofMRA are performed separately on the different vari-ables defining these properties, as explained in the fol-lowing sections.4.2 Multiresolution scheme for geometryStandard linear MRA schemes rely on smoothing oraveraging and difference filters for approximations anddetails, respectively. To preserve coherency of represen-tation of geometrical discontinuities — whatever theresolution — a special care is taken to avoid excessivesmoothing, while at the same time allowing the reversesynthesis. As the pillar grid format is used (see Section2.2), vertices are inevitably positioned along pillars. So,our multiresolution scheme for geometry data only fo-cuses on: – the z coordinates of the 8 vertices associated toeach node. According to the naming convention pre-sented in Figure 9, those 8 vertices can be differ-entiated according to their relative positions [Back(B)/Front (F), Bottom (B)/Top (T), Left (L)/Right(R)]; – the x and y coordinates of the nodes describingthe low (bottom) and high (top) extremities of allthe pillars (the x and y coordinates of intermediarynodes being implicit). The nodes are called here-inafter the floor and ceil nodes, respectively. (a) A node and its 8 sur-rounding cells. FTLBTL BB R FBL FTR BTR BB L FBR (b) Splitting view of thenode into its 8 vertices.
Fig. 9: Vertex naming with the Back (B)/Front (F),Bottom (B)/Top (T), Left (L)/Right (R) convention.Actual 3 D meshes can exhibit very irregular bound-aries. Hence, a Boolean field called Actnum may beassociated to each cell to inactivate its display (andits influence during simulations as well). It enables thedescription of either mesh boundaries (Figure 10), orcaves/overhangs. Resultantly, this
Actnum field mustbe carefully considered during the multiresolution anal-ysis of the geometry data, to avoid artifacts at lower res-olutions on frontiers between active and inactive cells(see Section 4.2.4 and Figure 16).Fig. 10: Mesh
Actnum field.By construction, most geological VMs have no hori-zontal fault, as there is no vertical gap between any twoadjacent layers of cells. For every node, each of the fourtop vertices has the same z coordinate as its counterpartbottom vertex. Therefore, from now on, our geometrymultiscale representation method only deals with the z coordinates of the bottom vertices BBL, BBR, FBL, Peyrot, Duval, Payan, Bouard, Chizat, Schneider, Antonini and FBR of each node.An instance of the decomposition shown in Figure 8can be implemented with the proposed two-step tech-nique depicted by arrows in Figure 11: – A non-linear and non-separable 2 D morphologi-cal wavelet transform applied on the nodes, inorder to detect the faults in the input VM, andthen to preserve their coherency in the lower res-olutions. This step relies on a fault segmentation within the input VM obtained by studying all pos-sible fault configurations for the top view of the VM(see Figure 12); – A non-linear 1 D wavelet transform applied onthe output of the above first step to analyze the z coordinates of the vertices along each pillar.The same 1 D wavelet transform is also applied onthe sets of x and y coordinates of the floor and ceil nodes, to complete the “horizontal” decomposition. Group of (2x2x4) cells
Morphological wavelet
Group of (4x4x4) cells(resolution l ) Group of (2x2x2) cells(resolution l-1 ) constrained 1D wavelet i j k Top view for fault segmentation
Fig. 11: HexaShrink multiresolution scheme for geom-etry: (left) input grid and its top view; (middle) out-put from the non-separable, non-linear 2 D morpholog-ical wavelet based on a fault segmentation (based onthe top view); (right) non-linear 1 D wavelet transformalong pillars (orange lines). This stage detects the faults in the original mesh, in or-der to preserve them during the morphological waveletanalysis. For each node, a dozen of fault configurations,depending on BBL, BBR, FBL, and FBR, is possible:fault-free (1), straight (2), corner (4), T-oriented (4) orcross (1), as illustrated in Figure 12.
Fault -free
BBL FBL BBL FBL Horizontal Vertical Top-right Bottom-right Top-left Bottom-left T-south T-east T-west Cross T-north
BBR BBL
FBR FBL BBR BBL
FBR FBL BBR BBL
FBR FBL BBR BBL
FBR FBL BBR BBL FBR FBL BBR BBL FBR FBL
BBR BBL
FBR FBL
BBR BBL
FBR FBL BBR BBL
FBR FBL
BBR BBL
FBR FBL
BBR BBL
FBR FBL
BBR BBL
FBR FBL
Fig. 12: The 12 possible fault configurations (in blacklines) at a given node.Each configuration depends on the four orientationsof the cardinal axes (north, south, east and west), whichare either active or inactive. For instance, the T-northconfiguration has its south axis inactive, while the threeremaining ones are active. Assuming that a fault con-figuration is z -invariant, meaning that the nodes be-longing to the same pillar present the same fault con-figuration, a single 2 D configuration map is sufficientto represent the fault configuration of the whole mesh,as illustrated by Figure 13. Horizontal
BBR
Top-left Vertical Bottom-right D T o p v i e w BBL FBRFBL BBRBBL FBRFBL BBRBBL FBRFBL BBRBBL FBRFBL
Fig. 13: Fault segmentation within the original mesh. D morphological wavelet transform The fault segmentation guides the multiresolution anal-ysis to preserve faults, as much as possible, all over thedecomposition process. The fault configuration of 4 as-sociated nodes at resolution l is used to predict theextension of the downsampled fault structure at reso-lution l − exaShrink: multiresolution reversible compression of hexahedral meshes for geoscience simulation models 9 This horizontal prediction is based on the logicalfunction OR ( ∨ ), computed on each side of each groupof 4 nodes. For instance, a resulting fault node config-uration contains a west axis if the fault configurationsof the 2 left nodes contain at least 1 west axis, as il-lustrated in Figure 14. By repeating the procedure forthe north, south and east axes of each resulting node,fault node configurations at lower resolutions are fullypredicted. This non-linear and peculiar choice is meantto maintain a directional flavor of orientated faults forflows; other choices could be devised, depending onphysical rules and geological intuitions.Finally, from this prediction, the node whose config-uration minimizes its distance with the predicted one,corresponds to the aforementioned approximation co-efficient, which will be part of the novel Z matrix atlower resolution l −
1. The same procedure can be ap-plied recursively until the wanted resolution.Fig. 14: Prediction of a fault node at resolution l − l ,orange ovals denoting ∨ operands. D wavelet transform This 1 D wavelet transform is applied on the output ofthe above horizontal 2 D morphological wavelet, to an-alyze the z coordinates of the 4 sets of vertices BDR,FDR, BDL and FDL separately, along each selected pil-lar. Hence, HexaShrink here decomposes at each scale z l into a subsampled pillar coordinate z l − and its associ-ated detail d l − . By geomodel construction, coordinatebehavior along the pillars is expected to be relativelysmooth. This entails the use of a modified, longer splinewavelet. The latter can be termed LeGall [31], or CDF5 / . 𝑧 𝑙−1 [1] 𝑧 𝑙 [2]𝑧 𝑙 [1]𝑧 𝑙 [0]𝑧 𝑙 [4]𝑧 𝑙 [3] 𝑧 𝑙−1 [0]𝑧 𝑙−1 [2]𝑑 𝑙−1 [0] P 𝑑 𝑙−1 [1] P UU 𝑑 𝑙−1 [2] P Details Approximation coefficients (lower resolution)
Fig. 15: Principle of the lifting scheme (Prediction andUpdate) for the rounded linear 1 D wavelet from (1)-(2), to analyze z coordinates of vertices BDR, FDR,BDL and FDL along each pillar, as well as x and y coordinates of the floor and ceil nodes.The lifting analysis operations Prediction and
Up-date are depicted by Figure 15. To retrieve respectivelythe sets of details d l and the z l coordinates at resolu-tion l − l , the following equations are used( ∀ n ∈ N ): d l − [ n ] = z l [2 n + 1] − (cid:22) z l [2 n ] + z l [2 n + 2]2 (cid:23) , (1) z l − [ n ] = z l [2 n + 0] + (cid:22) d l − [ n −
1] + d l − [ n ]4 (cid:23) , (2)where both dyadic integers and rounding are evident(see Section 4.1). With rounding, lifting schemes canthus manage integer-to-integer transformations [6]. Forsynthesis, to reconstruct resolution l from resolution l −
1, we only have to reverse the order and the signof the equations: z l [2 n ] = z l − [ n ] − (cid:22) d l − [ n −
1] + d l − [ n ]4 (cid:23) , (3) z l [2 n + 1] = d l − [ n ] + (cid:22) z l [2 n ] + z l [2 n + 2]2 (cid:23) . (4) A pertinent multiresolution on complex meshes requiresto cope with externalities that may hamper their han-dling: floor and ceil borders and outer boundaries (Fig-ure 10). First, to keep borders unchanged from the orig-inal mesh, throughout all resolutions, the following con-straints must be met: z l − [0] = z l [0] , (5) z l − [ n l − k −
1] = z l [ n lk − . (6)Both constraints can be fulfilled if one satisfies thefollowing conditions: – Floor border condition to meet (5): d l − [ −
1] = − d l − [0] , (7) – Ceil border condition to meet (6): d l − [ n l − k −
1] = − d l − [ n l − k − , ( n lk odd) (8) d l − [ n l − k −
1] = − d l − [ n l − k −
2] ( n lk even)+4 z l [ n lk − − z l [ n lk − . (9)To complete the MRA of the geometry, the samerounded 1 D wavelet as in Section 4.2.3 is finally appliedto the sets of floor and ceil nodes of the initial VM, toget the x and y coordinates of the extremities of theremaining pillars at the lower resolution.Second, the Actnum field should also be consideredto lessen mesh boundary artifacts. Indeed, severe dis-turbances may appear at lower resolutions if not wiselyprocessed during analysis, as shown in Figure 16. (a) Boundary artifactswithout proper
Actnum management. (b) Lower mesh res-olution with efficient
Actnum
Boolean valuescare-taking.
Fig. 16: Inadequate
Actnum fields management duringanalysis may lead to severe boundary artifacts (left)that can be dealt with (right) as examplified withmesh l . During our study,we find that one vertex at resolution l − l (Section4.2.2) are active. So, a cell at resolution l − × l .4.3 Multiresolution scheme for continuous propertiesOnce the geometry is coded, one can focus on associatedcontinuous properties. For scalar ones, a value p i ∈ R is associated to each cell i in the mesh. Consistently with the handling of cell blocks C of 2 × × Haar wavelet. The resolution l − l . The approximation coef-ficient p l − is thus the average value of the related eightproperty coefficients { p l , p l , . . . , p l } . The seven detailsrequired for synthesis are differences with respect to theapproximation coefficient: p l − = 18 (cid:88) n =1 p ln ; d l − n = p ln − p l − , ∀ n (cid:54) = 1 . To deal with real-valued (floating-point) properties,and avoid accuracy imprecision due to the divide op-erator, we introduce the following modifications. First,reals are mapped into integers up to a user-defined pre-cision, here with a 10 factor. Second, we disable thedivision by using a sum. The analysis system thus be-comes: p l − = (cid:88) n =1 p ln ; d l − n = 8 p ln − p l − , ∀ n (cid:54) = 1 , and the synthesis system turns into: p ln = 18 ( d l − n + p l − ) , ∀ n (cid:54) = 1 ; p l = p l − − (cid:88) n =2 p ln . Approximation and coefficients are stored as is. Torecover the accurately scaled values, the division oper-ator should however be applied as a simple linear post-processing.4.4 Multiresolution scheme for categorical propertiesWe finally complete the global mesh multiresolution de-composition with an original categorical-valued schemecalled modelet [2]. We assume that a mesh cell categorybelongs to a set of classes Ω = { ω , ω , ..., ω W } , tak-ing discrete values. The cell block C l { p l , p l , . . . , p l } thus contains, at resolution l , integers indexing cat-egories from Ω l . They take values in a subset of Ω .The multiresolution scheme is expected to produce, atlower resolutions, discrete values in embedded subsets: Ω ⊃ Ω − ⊃ · · · ⊃ Ω l ⊃ · · · . In other words, a cell cate-gory can only belong to an existing category at an upperresolution. We choose here the modal value (mode) i.e. ,the most frequently represented in C l . If | ω w | denotesthe cardinal of this class, then (cid:80) Ww =1 | ω lw | = |C l | = 8.We choose for the modelet: p l − = arg max {| ω lw | , ω lw ∈ Ω l } . It may happen that the above definition does not yielda unique maximum. If two or more categories dominate exaShrink: multiresolution reversible compression of hexahedral meshes for geoscience simulation models 11 a cell block, a generic approach consists in taking intoaccount its first block cell neighborhood (the surround-ing 26 cells, except at mesh borders and boundaries).We affect the dominant value in the first neighborhoodto p l − . In case of a draw again, the second-order sur-rounding can be used, iteratively. In practice for thepresented version of HexaShrink, we limit to the first-order neighborhood, and choose the lowest indexed cat-egory when the maximum is not unique. Equipped withthis unique lower resolution representative value, weproceed similarly to Section 4.3 for details, by using dif-ferences between original categories and the mode. Asclasses are often indexed by positive integers, a slightmotivation allows to get only non-negative indices. Byavoiding negative values, one expects a decrease in dataentropy of around 5 %, which benefits to compression.We thus change the sign of a detail coefficient ifand only if it generates a value out of the range of { ω , ω , . . . , ω W } , and then control this condition dur-ing reconstruction. So, all details { d l − n } for a cell block C are determined by: d l − n = ( − ( p ln − p l − < ∧ ((2 p l − − p ln ) / ∈ Ω ) × ( p ln − p l − ) . During synthesis, coefficients { p li } are obtained thanksto the closed-form equation: p ln = p l − + ( − ( p l − + d l − n ) / ∈ Ω ) × d l − n . ™ . The latter also en-compasses structural information required to exchangemodels, generating information overhead. Finally, de-tail simplification through multiscale decompositionsdoes not possess well-established quality metrics. Allof the above hampers exhaustive objective evaluations such as possible in image processing, where metrics andbenchmarks have been evaluated for decades.To evaluate the performance of HexaShrink, we baseour analysis on a set of seven geological meshes, withgeometries ranging from smooth to fractured, and di-verse categorical and continuous properties. Their maincharacteristics and properties are summarized in Table1. As will be seen, they appear representative enoughto allow one to derive consistent observations and con-clusions for different data handling purposes.They are initially stored in the GRDECL (“GRiDECLipse”) file format [42]. Originally complex geomod-els are thus described with details rendering their ge-ometry explicit and structured, an important featurefor geomodelers or flow simulation software (Petrel ™ ,SKUA-GOCAD ™ , Eclipse ™ . . . ).As observed in the state-of-the-art (Section 3.2), toour knowledge, reversible multiscale representations ofgeometry and properties — together with discontinu-ity preservation — of hexahedral meshes do not exist.Even if the standardized 3 D extension to the JPEG2000image compression format (termed J2K-3D) could pro-cess the properties as volumetric images, it is not perse suited to volume meshes, especially with categoricalproperties. We thus focus on visualizations, compar-isons with geomodeler upscaling capabilities, and theembedding of our multiscale decompositions into sev-eral all-purpose compression algorithms.The evaluation methodology is twofold. First, weexemplify the outcome of HexaShrink on meshes on ei-ther their geometry with a continuous and a categori-cal property at different dyadic scales. This reversibleframework is put into perspective with similar down-scaling processes in a reference geomodeler. Second, acomprehensive evaluation of lossless compression per-formance is provided, using state-of-the-art coders oneither the raw files or their multiscale decomposed coun-terparts.5.2 Reversible multiscale mesh representationFigures 17 and 18 present the reversible multiscale de-compositions generated by HexaShrink for mesh Actnum
Continuous Categorical1 93,600 80 × ×
26 No 4 .
62 MB 100 % Porosity Rock type2 1,000,000 100 × ×
100 No 42 .
46 MB 100 % — —3 36,816 59 × ×
16 Yes 1 .
46 MB 100 % — —4 210,000 100 × ×
21 Yes 7 .
88 MB 20 % — —5 450,576 149 × ×
16 Yes 22 .
73 MB 46 % Porosity, Permeability —6 5,577,325 227 × ×
305 Yes 274 .
57 MB 97 % Porosity Rock type7 13,947,600 240 × ×
197 Yes 580 .
94 MB 100 % Porosity Rock type
Table 1: Meshes chosen for evaluation: a compendium of their ontological characteristics and geological properties.to an anticlinal on the mesh at original resolution re-mains perceptible on the final “Lego brick” resolution.Looking at the porosity property (middle column), oneobserves how the values are progressively homogenizedon coarser hexes. Concerning the rock type (last col-umn), one observes that the modelet scheme tends tolocally maintain predominant categories resolution af-ter resolution, which is very satisfactory.On Figure 18, the much larger mesh × × ,these two structural discontinuities are still present,while keeping a good shape fidelity, globally. Concern-ing the attributes, the decompositions are also ade-quate.To emphasize further the capacity of our schemeto maintain coherency across the resolutions for theproperties, Figure 19 shows the evolution of the RockType distribution until the third resolution for mesh / ™ vsHexaShrinkSKUA-GOCAD ™ or PETREL ™ are frequently used ingeosciences to handle geological objects and to gener-ate meshes for flow simulation. These specific meshesdescribe structural discontinuities whose impact is sig-nificant on simulation. To obtain such meshes, the geo-model — a surface description of horizons or faults — isfitted into a grid at the desirable resolution. Pillar orien-tation is influenced by fault dip and cell layer thicknessis adapted to the distance between horizons. Additionalproperties can then be assigned to mesh cells: porosity,saturation, rock type. . . from well data or geological in-terpretation. Would one wish to lower the resolution,the process described above should be reiterated.A simpler alternative proposed by geomodelers con-sists in upscaling meshes. Such methods are usuallyflexible yet often ad-hoc , converting geometry and prop-erties in a non-reversible manner. Figure 20 confrontsmeshes ™ and HexaShrink. Hexa-Shrink tends to better preserve faults (colored in red),as compared to SKUA-GOCAD ™ . Figures emphasizean improved preservation of mesh borders, with an effi-cient management of Actnum throughout resolutions.Some artifacts may appear with SKUA-GOCAD ™ ’s up-scaling, which are automatically averted by HexaShrink,leading to nicer meshes at low resolution. As a sum-mary, HexaShrink, while being fully reversible at dyadic exaShrink: multiresolution reversible compression of hexahedral meshes for geoscience simulation models 13Geometry. Porosity. Rock type.Original mesh − − − − Fig. 17: Original mesh − − − − Fig. 18: Original mesh exaShrink: multiresolution reversible compression of hexahedral meshes for geoscience simulation models 15
Fig. 19: Evolution of the distribution of the
Rock type categories across resolutions for mesh / Mesh Level gzip bzip2 LZMA1 none 3 .
73 4 .
98 6 .
431 5 .
62 6 .
07 7 . .
67 6 . .
13 7 . .
442 none 3 .
23 8 .
41 10 .
121 6 .
49 10 .
82 11 . . .
58 12 . .
03 13 .
353 none 2 .
67 2 .
99 3 .
631 3 .
88 4 .
70 5 . . .
05 4 . .
93 5 . .
484 none 1 .
83 1 .
89 2 .
211 2 .
64 3 .
06 3 . .
76 3 . .
23 3 . .
655 none 2 .
46 2 .
55 3 .
331 3 .
14 2 .
83 3 . . .
26 2 . .
92 3 . .
816 none 2 .
31 2 .
25 3 .
041 3 .
31 3 .
53 4 . . .
24 4 . .
68 5 . .
737 none 3 .
20 5 .
98 12 .
521 5 .
42 7 .
07 8 . . .
72 7 . .
12 9 . . Table 2: Comparative lossless coding performanceswith compression ratios at different HexaShrink resolu-tion levels combining HexaShrink with gzip, bzip2 andLZMA.We exhaustively compare compression performancesin Table 2. For the sake of clarity, recall that differentcomputational methodologies exist: a compression ra-tio of 4 is given by the fraction between the sizes ofthe original file and the compressed one (the larger theratio, the better the compression). The latter can alsobe related to its inverse, the smallest file representing25 % of the raw data ( × (cid:0) − (cid:1) × i.e. by directly compressing the binary formats, gains in filesize are already observed, from 62 .
50 % ( (cid:0) − . (cid:1) × .
50 % for LZMA. We first remark thatimprovements sensitively increase as we use more recententropic coders, with only one exception for mesh
Fig. 20: After dyadic downsampling/upscaling, HexaShrink (bottom) better preserves faults, and manages non-active cells ( i.e. , with null
Actnum values) across scales, yielding nicer borders at each resolution, contrary toGOCAD. From left to right: resolution − −
2, and −
3, respectively. exaShrink: multiresolution reversible compression of hexahedral meshes for geoscience simulation models 17 nation of a 1-level HexaShrink associated with LZMAyields a compressed mesh twice as small (81 .
90 %) asa direct gzip compaction on the original mesh. This isadvantageous, as the proposed method either providesaccess to a two-fold downsampled mesh together witha smaller overall size.One can wish to have access to further levels of ap-proximation. The third line of each table block speci-fies the range of additional available levels (dependingon mesh size), here from 2–4, with the minimum andmaximum compression ratios attained. While we stillobserve a marginal improvement over a 1-level Hexa-Shrink (and again a slightly anomalous behavior formesh . .
35 fold compression.In rare cases, HexaShrink does not result in clearcompression improvement. For mesh . I of small integers, with offset o and scale s : s × I + o .LZMA’s superior capability owes to its capacity to cap-ture complex models of byte patterns. By contrast, witha wavelet decomposition, the affine relationship in sucha case is poorly captured throughout approximations,due to the rounding in wavelet lifting (Section 4.2.3).Hence, multiscale decompositions may slightly reducethe raw compression performance for meshes presentinginitial “numerical format” artifacts or illusory floating-point precision. This however does not hamper the us-ability of HexaShrink for storage and visualization, asthe direct access to a hierarchy of resolutions respect-ing discontinuities is granted, while already providingimpressive compression rates of about 8–9, superior tomost results for the others meshes. Speed performance is highly dependent on mesh size,discontinuity complexity, levels of details. For a baselineevaluation, a Java implementation was run on a laptopwith Intel Core i7-6820HQ CPU @ 2 .
70 GHz processorand 16 GB RAM. Each mesh (from our dataset in Table2) was compressed to the maximum level, and decom-pressed, twelve times. As the outcomes were relativelystable, they were averaged. Timings for analysis or syn-thesis alone, and cumulated with lossless encoding anddecoding, are summarized in Table 3. [A]nalysis [A]+gzip [A]+bzip2 [A]+LZMAMin. 2 .
80 3 .
06 5 .
69 6 . .
790 1 .
03 3 .
35 3 . Table 3: Cumulative lossless compression and decom-pression durations of HexaShrink with gzip, bzip2 andLZMA, in seconds.Analysis is slightly slower than synthesis, both tak-ing from a couple of seconds to a couple of minutes.Concerning the coders, gzip is the fastest, adding lit-tle overhead to HexaShrink speed. However, bzip2 orLZMA durations can reveal more expensive than anal-ysis alone. The largest mesh is compressed in a littleless than 13 min. What is more appealing in applica-tions is that only a few seconds are necessary for smallmeshes. Lossless decompression is very fast, even for thelargest mesh (13 s–20 s), as synthesis takes most of thetime. And LZMA is the fastest here, adding only 5 %overhead to synthesis, for recovering the whole origi-nal mesh. Decompressing lower resolutions only is ofcourse even faster. This is interesting, as here we obtaina beneficial asymmetrical compression-decompressionscheme, termed “compress once, decompress many”.Once a model is built, one can afford to compress it onlyonce, whatever the time it takes. Then, after transfer-ring, handling it with the benefit of the smaller sizes, de-compressing it many times is less expensive. The aboveperformance could be greatly improved with more in-volved acceleration techniques.As a result, on all tested examples, we demonstratethe possibility of storing independently multiscale meshproperties as approximations and details, while preserv-ing geometry (hence faults). This method additionallyhas an important benefit in data handling, visualiza-tion and compression, all at a reasonable computationalcost.
HexaShrink offers a comprehensive and efficient frame-work for a scalable representation of hexahedral mesheswith continuous and categorical properties, at dyadicresolutions. It is first dedicated to the visualization ofmassive structured meshes, as used in geosciences, thatcan contain geometrical discontinuities, to describe faultsfor instance. Four adapted multiresolution representa-tions are matched to the underlying nature of each datafield. They permit to decompose such specific meshesprogressively. In particular, this framework includes amorphological transform that takes into account geo-metrical discontinuities relative to any fault configura-tion, and preserve their rendering across resolutions,while maintaining structural coherency.HexaShrink can process any mesh size, thanks to aGPU-based out-of-core algorithm. This is crucial, giventhe constant evolution of data acquisition density thatyields increasingly massive and accurate dataset, as-sociated to more demanding simulations. The Hexa-Shrink decomposition is consistent with respect to meshrescaling in geosciences, and provides an option for a re-versible upscaling; [37] recently proposed such a wavelet-inspired scheme. Finally, it lends itself to an efficientlossless compression, which can be used for storage andtransfer.Perspectives can deploy into many directions. Moti-vated by preliminary progressive lossless compressionresults, we aim at developing a more versatile mul-tiresolution compression scheme, to manage the po-tential evolution of mesh geometry or properties overtime, with a special care for simulation-related qualitymetrics. As multiresolution analysis and synthesis werecomputed independently from compression, we also en-vision a better matched coding of approximations anddetails, for additional performance, toward lossy com-pression [32, 36].
Acknowledgements
The authors would like to thank C. Dawson and M. F.Wheeler for their help, as well as the reviewers whosecomments helped improve the compression performanceassessment and comparison.
References
1. Abdelfattah, M.S., Hagiescu, A., Singh, D.: Gzip on achip: high performance lossless data compression on FP-GAs using OpenCL. In: Proc. Int. Workshop OpenCL(2014). DOI 10.1145/2664666.2664670 17 2. Antonini, M., Payan, F., Schneider, S., Duval, L., Pey-rot, J.-L.: Method of exploitation of hydrocarbons of anunderground formation by means of optimized scaling(2017). Patent 103. Bey, J.: Tetrahedral grid refinement. Computing (4),355–378 (1995). DOI 10.1007/BF02238487. 64. Boscard´ın, L.B., Castro, L.R., Castro, S.M., Giusti, A.D.:Wavelets bases defined over tetrahedra. INSTEC J. Com-puter Science and Technology (1), 46–52 (2006) 65. Bruekers, F.A.M.L., van den Enden, A.W.M.: New net-works for perfect inversion and perfect reconstruction.IEEE J. Sel. Areas Comm. (1), 129–137 (1992). DOI10.1109/49.124464 66. Calderbank, A.R., Daubechies, I., Sweldens, W., Yeo,B.L.: Wavelet transforms that map integers to integers.Appl. Comput. Harmon. Analysis (3), 332–369 (1998).DOI 10.1006/acha.1997.0238. 97. Cannon, S.: Reservoir Modelling: A Practical Guide.John Wiley & Sons (2018) 28. Caumon, G., Gray, G., Antoine, C., Titeux, M.O.: Three-dimensional implicit stratigraphic model building fromremote sensing data on tetrahedral meshes: Theory andapplication to a regional model of La Popa basin, NEMexico. IEEE Trans. Geosci. Remote Sens. (3), 1613–1621 (2013). DOI 10.1109/tgrs.2012.2207727. 29. Chaux, C., Duval, L., Benazza-Benyahia, A., Pesquet,J.C.: A nonlinear Stein based estimator for multichan-nel image denoising. IEEE Trans. Signal Process. (8),3855–3870 (2008). DOI 10.1109/TSP.2008.921757 610. Chaux, C., Pesquet, J.C., Duval, L.: Noise covarianceproperties in dual-tree wavelet decompositions. IEEETrans. Inform. Theory (12), 4680–4700 (2007). DOI10.1109/TIT.2007.909104 611. Chen, D., Chiang, Y.J., Memon, N., Wu, X.: Geometrycompression of tetrahedral meshes using optimized pre-diction. In: Proc. Eur. Sig. Image Proc. Conf. (2005).512. Chizat, L.: Multiresolution signal compression: Explo-ration and application. Master’s thesis, ENS Cachan(2014) 6, 1213. Cohen, A., Daubechies, I., Feauveau, J.C.: Biorthogonalbases of compactly supported wavelets. Commun. ACM (5), 485–560 (1992). DOI 10.1002/cpa.3160450502 914. Courbet, C., Isenburg, M.: Streaming compression of hex-ahedral meshes. Vis. Comput. (6-8), 1113–1122 (2010).DOI 10.1007/s00371-010-0481-7. 515. Danovaro, E., De Floriani, L., Lee, M.T., Samet, H.: Mul-tiresolution tetrahedral meshes: An analysis and a com-parison. In: Proc. Shape Modeling International, pp. 83–91 (2002). DOI 10.1109/SMI.2002.1003532. 516. Dupont, F., Lavou´e, G., Antonini, M.: 3D mesh compres-sion. In: L. Lucas, C. Loscos, Y. Remion (eds.) 3D Videofrom Capture to Diffusion. Wiley-ISTE (2013) 217. Foster, I., Ainsworth, M., Allen, B., Bessac, J., Cappello,F., Choi, J.Y., Constantinescu, E., Davis, P.E., Di, S.,Di, W., Guo, H., Klasky, S., Dam, K.K.V., Kurc, T.,Liu, Q., Malik, A., Mehta, K., Mueller, K., Munson, T.,Ostouchov, G., Parashar, M., Peterka, T., Pouchard, L.,Tao, D., Tugluk, O., Wild, S., Wolf, M., Wozniak, J.M.,Xu, W., Yoo, S.: Computing just what you need: On-line data analysis and reduction at extreme scales. In:Lecture Notes in Computer Science, pp. 3–19. SpringerInternational Publishing (2017). 118. Gumhold, S., Guthe, S., Straßer, W.: Tetrahedral meshcompression with the cut-border machine. In: Proc. IEEEVisualization Conf., pp. 51–58 (1999). 4exaShrink: multiresolution reversible compression of hexahedral meshes for geoscience simulation models 1919. Gumhold, S., Straßer, W.: Real time compression of tri-angle mesh connectivity. In: Proc. ACM SIGGRAPHComput. Graph., pp. 133–140 (1998). DOI 10.1145/280814.280836. 420. Hoppe, H., DeRose, T., Duchamp, T., McDonald, J.,Stuetzle, W.: Mesh optimization. In: Proc. ACM SIG-GRAPH Comput. Graph., pp. 19–26 (1993). DOI10.1145/166117.166119. 521. Ibarria, L., Lindstrom, P., Rossignac, J.: Spectral pre-dictors. In: Proc. Data Compression Conf., pp. 163–172(2007). DOI 10.1109/DCC.2007.72. 522. Isenburg, M., Alliez, P.: Compressing hexahedral volumemeshes. Graph. Model. (4), 239–257 (2003). DOI10.1016/s1524-0703(03)00044-4. 4, 523. Isenburg, M., Gumhold, S.: Out-of-core compression forgigantic polygon meshes. In: Proc. SIGGRAPH Int.Conf. Comput. Graph. Interactive Tech., pp. 935–942(2003). DOI 10.1145/1201775.882366. 1224. Isenburg, M., Lindstrom, P., Gumhold, S., Shewchuk, J.:Streaming compression of tetrahedral volume meshes. In:Proc. Graphics Interface, pp. 115–121 (2006). 525. Isenburg, M., Lindstrom, P., Snoeyink, J.: Streamingcompression of triangle meshes. In: Proc. EurographicsSymp. Geom. Process., vol. 255, pp. 111–118 (2005). 526. ITU-T T.809: JPEG2000 image coding system: Exten-sions for three-dimensional data (2011). ISO/IEC 15444-10:2011 1227. Jacques, L., Duval, L., Chaux, C., Peyr´e, G.: A panoramaon multiscale geometric representations, intertwiningspatial, directional and frequency selectivity. Signal Pro-cess. (12), 2699–2730 (2011). DOI 10.1016/j.sigpro.2011.04.025. 628. Kober, C., M¨uller-Hannemann, M.: A case study inhexahedral mesh generation: Simulation of the humanmandible. Eng. Comput. (3), 249–260 (2001). DOI10.1007/pl00013389 229. Kovaˇcevi´c, J., Goyal, V., Vetterli, M.: Signal ProcessingFourier and Wavelet Representations (2012) 630. Krivograd, S., Trlep, M., ˇZalik, B.: A hexahedral meshconnectivity compression with vertex degrees. Comput.Aided Des. (12), 1105–1112 (2008). DOI 10.1016/j.cad.2008.10.013. 531. Le Gall, D., Tabatabai, A.: Sub-band coding of digital im-ages using symmetric short kernel filters and arithmeticcoding techniques. In: Proc. Int. Conf. Acoust. SpeechSignal Process. (1988). DOI 10.1109/icassp.1988.196696932. Liang, X., Di, S., Tao, D., Li, S., Li, S., Guo, H.,Chen, Z., Cappello, F.: Error-controlled lossy compres-sion optimized for high compression ratios of scientificdatasets. In: IEEE Int. Conf. Big Data (2018). DOI10.1109/bigdata.2018.8622520 1833. Lie, K.A.: An Introduction to Reservoir Simulation UsingMATLAB. User Guide for the Matlab Reservoir Simula-tion Toolbox (MRST) (2016). SINTEF ICT, Departe-ment of Applied Mathematics 334. Lie, K.A., Møyner, O., Natvig, J.R., Kozlova, A.,Bratvedt, K., Watanabe, S., Li, Z.: Successful applicationof multiscale methods in a real reservoir simulator envi-ronment. Computat. Geosci. (5-6), 981–998 (2017).DOI 10.1007/s10596-017-9627-2 235. Lindstrom, P., Isenburg, M.: Lossless compression of hex-ahedral meshes. In: Proc. Data Compression Conf., pp.192–201 (2008). DOI 10.1109/dcc.2008.12. 536. Lu, T., Liu, Q., He, X., Luo, H., Suchyta, E., Choi, J.,Podhorszki, N., Klasky, S., Wolf, M., Liu, T., Qiao, Z.:Understanding and modeling lossy compression schemes on HPC scientific data. In: IEEE International Paralleland Distributed Processing Symposium (2018). DOI10.1109/ipdps.2018.00044. 1837. Misaghian, N., Assareh, M., Sadeghi, M.: An upscal-ing approach using adaptive multi-resolution upgriddingand automated relative permeability adjustment. Com-putat. Geosci. (1), 261–282 (2018). DOI 10.1007/s10596-017-9688-2 1838. Nelson, M., Gailly, J.L.: The Data Compression Book.Wiley (1995) 1539. Owen, S.J., Brown, J.A., Ernst, C.D., Lim, H., Long,K.N.: Hexahedral mesh generation for computational ma-terials modeling. Procedia Eng. , 167–179 (2017).DOI 10.1016/j.proeng.2017.09.803 240. Pajarola, R., Rossignac, J., Szymczak, A.: Implantsprays: Compression of progressive tetrahedral mesh con-nectivity. In: Proc. IEEE Visualization Conf., pp. 299–305 (1999). DOI 10.1109/VISUAL.1999.809901. 541. Perrons, R.K., Jensen, J.W.: Data as an asset: What theoil and gas sector can learn from other industries about”big data”. Energy Pol. , 117–121 (2015). DOI 10.1016/j.enpol.2015.02.020. 242. Pettersen, Ø.: Basics of reservoir simulation with theEclipse reservoir simulator. Department of Mathemat-ics, University of Bergen, Norway (2006). Lecture Notes1143. Peyrot, J.-L., Duval, L., Schneider, S., Payan, F., An-tonini, M.: (H)exashrink: Multiresolution compression oflarge structured hexahedral meshes with discontinuitiesin geosciences. In: Proc. Int. Conf. Image Process., pp.1101–1105 (2016). DOI 10.1109/ICIP.2016.7532528 144. Rao, R.M., Bopardikar, A.S.: Wavelet Transforms: Intro-duction to Theory and Applications. Prentice Hall (1998)645. Rezapour, A., Ortega, A., Sahimi, M.: Upscaling of geo-logical models of oil reservoirs with unstructured gridsusing lifting-based graph wavelet transforms. Transp.Porous Med. , 661–684 (2019). DOI 10.1007/s11242-018-1219-7 646. Røe, P., Hauge, R.: A volume-conserving representationof cell faces in corner point grids. Computat. Geosci. (3), 453–460 (2016). DOI 10.1007/s10596-015-9500-0.347. Rossignac, J.: Edgebreaker: Connectivity compression fortriangle meshes. IEEE Trans. Visual Comput. Graph. (1), 47–61 (1999). DOI 10.1109/2945.764870. 448. Salomon, D., Motta, G.: Handbook of Data Compression.Springer (2009) 4, 1549. Staadt, O.G., Gross, M.H.: Progressive tetrahedraliza-tions. In: Proc. IEEE Visualization Conf., pp. 397–402(1998). DOI 10.1109/VISUAL.1998.745329. 550. Sweldens, W.: The lifting scheme: a custom-design con-struction of biorthogonal wavelets. Appl. Comput. Har-mon. Analysis (2), 186–200 (1996). DOI 10.1006/acha.1996.0015. 651. Szymczak, A., Rossignac, J.: Grow & fold: Compress-ing the connectivity of tetrahedral meshes. Comput.Aided Des. (8-9), 527–537 (2000). DOI 10.1016/S0010-4485(00)00040-3. 452. Touma, C., Gotsman, C.: Triangle mesh compression. In:Proc. Graphics Interface Conf., pp. 26–34. Vancouver,Canada (1998). DOI 10.20380/GI1998.04 4, 553. Witten, I.H., Neal, R.M., Cleary, J.G.: Arithmetic codingfor data compression. Commun. ACM30