PySubdiv 1.0: open-source geological modelling and reconstruction by non-manifold subdivision surfaces

. Sealed geological models are commonly used as an input to process simulations, for example, in hydrogeological or 10 geomechanical studies. Creating these meshes often requires tedious manual work – and it is, therefore, difficult to adjust a once-created model. In this work, we propose a flexible framework to create and interact with geological models using explicit surface representations. The essence of the work lies in the determination of the control mesh and the definition of semi-sharp creases values, which, in combination, enable the representation of complex structural settings with a low number of control points. We achieve this flexibility through the adaptation of recent algorithms from the field of computer graphics to the 15 specific requirements of geological modelling, specifically the representation of non-manifold topologies and sharp features. We combine the method with a particle swarm optimisation (PSO) approach to enable the automatic optimisation of vertex position and crease sharpness values. The result of this work is implemented in an open-source software (PySubdiv) for reconstructing geological structures while resulting in a model which is (1) sealed/watertight, (2) controllable with a control mesh, and (3) topologically similar to the input geological structure. Also, the reconstructed model may include a fewer number 20 of vertices compared to the input geological structure, which results in reducing the cost of modelling and simulation. In

addition to enabling a manual adjustment of sealed geological models, the algorithm also provides a method for the integration of explicit surface representations in inverse frameworks and the consideration of uncertainties.

Introduction
Parametric surface-based representation is one of the major approaches to the surface representation of geological objects (De Paor 1996, Farin and Hamann 1997, De Kemp 1999).Several studies considered a surface-based approach in geological modelling since the outstanding features of the structure, e.g.heterogeneity, are explicitly demonstrated by the surfaces of the boundary without depending on grid cells (Jacquemyn et al. 2019, Moulaeifard et al. 2023).It is worth mentioning that Jacquemyn et al. (2019) comprehensively investigate the advantages of surface-based modelling compared to grid-based representation.Botsch et al. (2010) consider NURBS (Non-Uniform Rational B-Splines) and subdivision surfaces as two major methods of parametric surface-based representation for generating controllable free-form surfaces.Subdivision surfaces convert an initial mesh (control mesh) to a smooth mesh by successive refinements until the final mesh is sufficiently fine.The final smooth mesh is a controllable free-form surface which can be controlled easily by the control points (vertices of the initial mesh).Although both NURBS and Subdivision surfaces produce controllable free-form surfaces, subdivision surfaces give freedom from topological limitations, whereas NURBS underscore the smooth manipulation of the model (Cashman 2010).
Previous studies of geological modelling have dealt with NURBS (Paluszny et al. 2007, Börner et al. 2015, Jacquemyn et al. 2019).However, NURBS surfaces suffer from two limitations: (1) for generating a model with complex topology, several NURBS patches have to be stitched together, which leads to difficulty in modelling, and (2) modification of NURBS surfaces is challenging in the complex cases since NURBS surfaces are based on grid structures (DeRose et al. 1998, Sederberg et al. 2008, Botsch et al. 2010) Although extensive research has been carried out on exploiting NURBS for geological and reservoir modelling, researchers have not considered subdivision surfaces in much detail.Lévy and Mallet (1999) and Mallet (2002) investigate the subdivision surface method, which exploits discrete smooth interpolation (D.S.I) as the refinement scheme at each iteration.In their approach, all vertices (except the control vertices) are free to move in each iteration since D.S.I is a discrete fairing method.
Although their work deserves appreciation, they mention that the formal proof of the convergence of the series of recursively refined meshes is an open question; however, it can be observed experimentally.Also, they do not investigate the compatibility of subdivision surfaces for the structures with non-manifold topologies, which is broadly observed in complex geological and reservoir modelling, e.g.connection between faults and horizons, channel intersection, and hierarchical layered structures.In this paper, we focus on the non-manifold subdivision surface method based on the Loop algorithm, in which the convergence of the series of refined meshes is formally proved (Loop 1987).
One of the typical examples of non-manifold surfaces in geological modelling is where several faces of the mesh share one edge (Fig. 2).Representations of non-manifold structures require complex algorithms (Rossignac & Cardoze, 1999).The lack of supporting non-manifold surfaces is one of the limitations of classical subdivision surfaces.Ying and Zorin (2001) modified the subdivision surfaces algorithm and made it compatible with the non-manifold topology which is implemented in the core of PySubdiv.Classic subdivision surface algorithms tend to generate smooth surfaces.However, complex geological structures inherently consist of creases and corners (e.g., sharp fault bend), which require modified subdivision algorithms.Subdivision surfaces with the semi-sharp creases approach generate creases on the mesh by increasing the resistance of the vertices to the smoothing procedure.In practice, the resistance of the vertices can be created by applying the crease sharpness value to the edges consisting of the respective vertices.The higher the crease sharpness value of the edge is, the sharper the creases on the structure.It is worth mentioning that PySubdiv exploits the semi-sharp creases for a better representation of the creases of complex geological structures.Furthermore, the subdivision surfaces are used for reconstruction by fitting smooth surfaces to mesh or dense point cloud data (Ma et al. 2015).The goal of reconstruction in this paper is to generate the control mesh that gives control over the reconstructed mesh with a few numbers of control points (Fig. 1).In other words, the reconstructed mesh not only is fitted to the input mesh but also can be controlled by vertices of the control mesh.Moulaeifard et al. (2023) investigate the reconstruction of geological structures by non-manifold subdivision surfaces.However, they manually fit the smooth surfaces to the input data, which is a time-consuming process in complex geological structures.This paper presents the automatic reconstruction of geological structures which is implemented in the PySubdiv.
The structure of the control mesh is based on two critical variables: (1) the position of the control points and (2) the crease sharpness value of each edge.Therefore, for fitting the smooth mesh to input data, the position of the control points and crease sharpness values have to be set.Several researchers have studied the simplification method to estimate the control mesh from the input mesh (Hoppe et al. 1994, Suzuki et al. 1999, Wu et al. 2017).However, Ma et al. (2015) advise exploiting the distinguished features of the input mesh instead of the simplification method since simplification would be time-consuming for structures with complex topology.It is worth mentioning that extensive research has been carried out on exploiting the key features of the input mesh for generating the control mesh; Ma et al. (2015) use the umbilics and ridges.Also, Marinov and Kobbelt (2005) and Kälberer et al. (2007) consider the curvatures for generating the control mesh.
After estimating the control mesh, we have to optimise the control mesh to find the best fit to input mesh.Wu et al. (2017) investigate automatic fitting to solve the optimisation problem by using the augmented Lagrangian method.They show that their method provides significant gains over previous works, e.g.Marinov and Kobbelt (2005), by generating the reconstructed mesh with fewer control points while consisting of comparable errors.Therefore, for the sake of efficiency, PySubdiv makes use of the automatic reconstruction method proposed by Wu et al. (2017), which is explained in section 2.4.The final reconstructed structure is controllable by control points, sealed (or watertight), and topologically similar to the input model.
To date, there is no practical software to generate and manipulate sealed geological models using subdivision surface approaches.We attempt to close this gap with the approach described in this paper and the implementation in the accompanying software package PySubdiv.This study investigates the advantages and limitations of PySubdiv for the modelling and reconstruction of geological and reservoir structures with a case study.Also, PySubdiv can export the final files as 3D objects based on common object formats, e.g.obj, which can be read by most computer graphics and meshing software.

Methods
The core functionality of PySubdiv consists of four fundamental parts, which are investigated in the following section: (1) subdivision surface algorithm, (2) modelling with semi-sharp creases, (3) supporting non-manifold topology, and (4) automatic reconstruction.

Subdivision surface algorithm
The subdivision surface algorithm refines the input (control) mesh to generate the final desirable mesh based on mathematical rules (Halstead et al. 1993, Reif 1995, Stam 1998, Peters 2015).Subdivision surfaces follow two steps at each refinement: (1) splitting step, which includes implementing the new vertices on the surface, and (2) averaging step for updating the location of the vertices.There are several subdivision surface schemes based on different criteria, e.g.type of input mesh (triangular or quad) and the approach for refinement.Loop Scheme (Loop 1987) is one of the common subdivision schemes for triangular meshes which is already implemented in PySubdiv.In the following, the Loop algorithm is explained.Loop (1987) defines the loop algorithm to generate smooth surfaces for triangular meshes by using splitting and averaging steps in each refinement stage.In the splitting step, a new vertex is inserted on the midpoint of each edge (blue vertex in Fig. 3), which results in the splitting of each triangle of the control mesh into four triangles (Fig. 3).(2) Also, to compute the new location of the midpoint of an edge (ℎ) enclosed by four existing vertices ( 1 ,  2 ,  3 ,  4 ) the Loop algorithm use (Fig. 4b):

Piecewise Smooth Subdivision Surfaces
Representation of the objects consisting of sharp features by smooth subdivision algorithms leads to unsatisfactory results.Hoppe et al. (1994) propose to add new rules for the representation of the creases and corners to the Loop subdivision scheme in which the new location of the vertex  depends on the number of connected sharp edges as represented in Table 1 (DeRose et al. 1998).

Table 1
The rules for updating the positions of vertices connected to sharp edges based on piecewise smooth subdivision surfaces (DeRose et al.The piecewise subdivision surface method presents an acceptable solution for generating the sharp regions of the model.
However, in complex geological modelling, it is vital to model the semi-sharp regions which are not quite sharp.

Modeling with semi-sharp creases
The creases and corners can be made by considering the crease sharpness values for the edges of the control mesh (Fig. 5).
This value can be between zero to one, which indicates zero and infinite sharpness, respectively.Modelling different geometric objects become flexible by regulating crease sharpness.As mentioned in section 2.1.2,Hoppe et al. (1994) define the new rules for generating sharp regions during the subdivision procedure.DeRose et al. ( 1998) generalize the method of Hoppe et al. (1994) for updating the position of the vertices based on semi-sharp subdivision surfaces (Table 2).

Table 2
The rules of semi-sharp creases scheme for updating the positions of vertices (DeRose et al. 1998).

Number of adjacent sharp edges
The rule for updating the position of the vertex 0 or 1 smooth subdivision rules (section 2.1.1) > 2 corner rule of Table 1 2 crease rule of Table 1 Fig. 5a represents the control mesh with eight control points (red vertices).Also, Fig. 5b shows the final smooth mesh after applying three times subdivision surfaces when all edges of the control mesh have no crease sharpness values.However, Fig.

Supporting non-manifold topology
Non-manifold topologies are extensively observed in complex geological and reservoir modelling.The classic subdivision surfaces can not support non-manifold topology since there is at least one irregular vertex or/and edge.Ying and Zorin (2001) propose the non-manifold subdivision surface algorithm, which supports a wide variety of non-manifold topology problems.A full explanation of this method is beyond the scope of this paper, and in the following, we mention part of their work which is practical in geological modelling.However, Moulaeifard et al. (2023) extensively investigated nonmanifold subdivision surfaces for geological and reservoir modelling with multiple examples.Ying and Zorin (2001) categorize the vertices of the mesh into three types and define the rule for updating the position of vertices (Table 3); (1) regular vertices, (2) simple singular non-manifold vertices that are connected only to two neighbourhood vertices by two shared non-manifold edges (Fig. 7) and, (3) complex singular non-manifold vertices which are the other nonmanifold vertices.PySubdiv exploits the non-manifold subdivision surfaces for geological modelling.The rules for updating the position of the vertices based on the non-manifold subdivision surface algorithm (Ying and Zorin 2001).

Vertex type
The rule for updating the position of the vertices, including non-manifold

Automatic reconstruction
The geological model can be constructed by various methods, e.g.marching cube from the implicit model.As mentioned in Section 1, the goal of surface reconstruction in this paper is to make the geological model manageable with a few numbers of control points (of the control mesh).In order to reconstruct the input mesh, first, the control mesh should be estimated based on the salient features of the input mesh (e.g., the location of minima, maxima, and saddle points).Exploiting the salient features of the input mesh leads to preserving the critical features of the input mesh (Ma et al. 2015).Then, the location of the control points and crease sharpness values of the edges are optimized to find the best fit of the reconstructed and input meshes.
It is worth mentioning that applying the subdivision surface algorithm to the watertight control mesh results in a watertight reconstructed mesh.Fig. 8 represents the general workflow of the reconstruction process of an anticlinal structure.The input mesh contains two triangulated surfaces which are generated by GemPy software (de la Varga et al. 2019).In order to obtain the initial control mesh, vertices on the input mesh are extracted (purple spheres) to generate the edges, faces and, finally, a watertight initial control mesh.The initial control mesh runs through the optimisation process to be optimised.The reconstructed mesh is generated by applying the subdivision surface algorithm to the optimised control mesh.Also, the following conceptualized approach can be mentioned for the reconstruction process (Fig. 9): A noisy input mesh ℳ (black curve), e.g. a rough geological interface, is generated from laser scanning or geophysical measurements.A coarse control mesh К (orange curve) is generated based on the features of the input mesh (e.g., minima and maxima), which contains the control points vi.An unfitted subdivision surface SК is generated based on the control mesh.On the subdivision surface, s(vi) indicates the position of vertices vi after subdivision.Hence, the distances di between s(vi) and the input surface can be calculated by projecting the normal from the input mesh or finding the closest vertices (Marinov and to the edges.and automatic methods, the first step is the estimation of the control points () with outstanding features of the input geological mesh, e.g.local maxima and minima.PySubdiv proposes some critical points of the input mesh as suitable candidates for the generation of the control mesh.However, it is always possible to arbitrarily select a set of control points.Then PySubdiv used the Delaunay triangulation approach to generate the initial control mesh.The second step is the reconstruction through optimising the location of control points and crease sharpness value for each edge of the control mesh to find the best fit of the reconstructed mesh to the input mesh.
Given the input geological mesh consisting of  vertices  = { 1 ,  2 ,  3 , . .  }, the control mesh consisting of  vertices  = { 1 ,  2 ,  3 , . .  } and  edges with the crease sharpness values ℎ = {ℎ 1 , ℎ 2 , ℎ 3 , . .ℎ  }.The goal is to minimize the vector , which indicates the difference between  and the projection of  onto the reconstructed (subdivided) mesh.Hoppe et al. (1994) mention that each vertex of the reconstructed mesh can be approximated as an affine combination of control points ().

CORE of PySubdiv
As mentioned in section 2, the core functionality of PySubdiv includes (1) subdivision surface algorithm, (2) modelling with semi-sharp creases, (3) supporting non-manifold topology, and (4) automatic reconstruction.PySubdiv is written with an object-oriented approach using Python programming language (Rossignac and Cardoze 1999).Also, PySubdiv exploits different varieties of open-source external libraries which are integrated into the core.Table 4 represents the main external libraries implemented in PySubdiv.

Table 4
The main external libraries implemented in PySubdiv.PySubdiv accepts the structures with triangle mesh, while the input mesh can be generated by any arbitrary method (e.g.implicit or explicit).Also, the input mesh can be either watertight or non-watertight.If the input mesh is not triangulated, PySubdiv can convert the input mesh to the triangle mesh by using the triangulate function of the PyVista library (Sullivan and Kaszynski 2019).

Library
The estimation of the position and number of control points is the key stage in reconstruction which leads to the generation of the watertight modelling by providing the control points at surface intersections.PySubdiv offers some candidates to the user based on the features of the input geological structure, e.g.minima, maxima and boundaries.However, it is always possible for the user to select or generate the control points based on the requirements for the interpretation.For example, the geological models generated based on real data may be associated with uncertainties (Wellmann and Caumon 2018).The user can consider some control points at the suspicious locations regardless if the locations are on the boundary or body part of the layer.
In order to show the general workflow of generating control mesh in PySubdiv, an anticlinal structure model is investigated.First, the top mesh is selected (the colour changes from green to red after selection).Then, the control points should be chosen among the vertices of the selected mesh (Fig. 11a).The individual points can be selected by pressing the "P" button on the keyboard while the mouse cursor hovers above them.In this example, the maximal points on the anticline axis have been chosen, as well as an additional point in the centre.Three points on the flanks have been sampled to support the control mesh at the turning point.Further important points lie in the corner of the mesh.The next step is to triangulate the selected points.
"boundary" during the triangulation.If declined, the Delaunay algorithm is executed directly on the automatically detected boundaries.However, unwanted orthogonal faces might appear.It is worth noting that the polygon boundary for the Delaunay algorithm can be defined by the user.Finally, the triangulated control mesh is generated (Fig. 11b).Repeating the mentioned steps for the second input surface yields the second simplified surface for the control mesh.In order to create a watertight mesh, four additional sides should be created for the anticline model.Therefore, the boundary vertices of the two simplified surfaces are selected and triangulated.The vertices should be sampled from the simplified surfaces to generate the watertight control mesh (Fig. 12a).At the moment, the four sides must be meshed individually and then merged together.When all four sides are generated, the different meshes can be stitched together by merging.Finally, the sealed control mesh is generated and shown in (Fig. 12b).The next step is to assign the crease sharpness value for each edge of the initial control mesh (Fig. 12b).Hoppe et al. (1994), Wu et al. (2017) offer to consider a threshold angle (θ0) as a criterion for tagging the edges of the initial control mesh.For each edge, it is tagged as sharp (e.i., crease sharpness value equal to one) if the angle between the normal of two adjacent faces (θe) is more than the threshold angle (θe > θ0).It is worth reminding that this value is just an initial estimation for crease sharpness values.
(III) Reconstruction (Optimisation) of Control Mesh (automatically done by software) In the first step of optimisation, the generated (initial) control mesh and the original input meshes are imported as the input data for the reconstruction algorithm (Algorithm 1)."MeshOptimizer" class instantiates an object based on the input data, which can start the fitting process with the "optimize" method.It is worth mentioning that the final subdivided (reconstructed) mesh can be exported by applying the subdivision surface algorithm to the optimized control mesh.

Case study
As a case study, a part of Upper Rhine Graben (URG) is reconstructed by PySubdiv.URG is a long geological structure in the central part of the European Cenozoic rift system that contains geothermal energy resources.The data set of the URG consists of several different geological units (grid nodes) published by Freymark et al. (2020), which consists of 616464 individual nodes (Fig. 13).The dimensions of the original model are 292 km in the x-direction, 525 km in the y-direction and 130 km deep (z-direction).
As mentioned in section 3, the input data of PySubdiv should consist of the triangular mesh, which can be generated by any arbitrary method or software.In this case study, PyVista generates the triangular mesh from the individual grid nodes, resulting in 18 different surfaces and ten volumes.

Estimation of the control mesh
The initial watertight control mesh is prepared based on the prominent features of the input mesh, such as (a) minimal and maximal points and (b) surface intersections of different layers to ensure that the final mesh is watertight.Fig. 14 represents the control mesh consisting of 832 control points distributed over 18 individual surfaces.
Also, the edges of the control mesh, which consist of a threshold angle of 80° (section 3.1), are considered sharp and given the crease sharpness values equal to one.Other edges are considered smooth and assigned crease sharpness values equal to zero.
From our experience, exploiting the angle of 80° as the threshold angle in this case study leads to acceptable results.However, this angle can be different depending on the complexity of the model.

Optimisation of the control mesh
The reconstructed mesh is generated after applying two subdivisions to the control mesh (Fig. 15).It consists of roughly 15000 vertices.In order to evaluate the reconstructed structure, the distance metric between the points of the original and the reconstruction model is computed.Red-coloured areas imply high deviations of the subdivided surface from the original model.
The mean error for the whole domain is approximately 496 m ± 362 m which is around 0.6% of the total elevation height.
Areas of high error concentrate mainly in regions where two individual geological units are connected and further on the lower boundary of the Lithospheric Mantle (lower boundaries of Fig. 13).On the Lithospheric Mantle, the highest errors are up to 5700 m which is around 4.75% of the total elevation (approximately 120 km).Most parts of the model consist of small deviations indicated by the dark blue colour, especially the higher elevated layers are well-fitted.

External libraries limitation
From the computational point of view, the major goal of PySubdiv is the calculation of the new location for the vertices of the control mesh based on the subdivision algorithm.Therefore, NumPy (Van Der Walt et al. 2011) and Scipy (Virtanen et al. 2020) libraries play key roles in the core of the software, which can be the source of problems when the input mesh is large, e.g.limited memory.
Although the structure of the PySubdiv is planned to avoid the generation of unnecessary big matrixes, the initial input data can also help to tackle this problem.For example, the suitable number of subdivisions is one of the key input data in reconstruction (the number of subdivisions controls the smoothness of the final mesh, which is different from the number of iterations of the optimisation process).Applying a small number of subdivisions can not guarantee the acceptable generation of smooth and semi-sharp parts of a reconstructed mesh.However, applying a large number of subdivisions not only does not make significant differences in the smoothness of the reconstructed mesh but also remarkably increases the unnecessary calculation costs since more subdivisions mean dealing with more vertices.

Conclusion
This study illustrated the framework (PySubdiv) to generate suitable control meshes and fitted reconstructed meshes for complex geological structures and reservoir models based on the non-manifold subdivision surface algorithm.The reconstructed mesh is watertight and topologically similar to the input mesh.Also, the control mesh consists of those control points which play a key role in the flexibility and management of the reconstructed mesh.Subdivision surfaces, unlike spline surfaces, support arbitrary topology, which gives more freedom to the user during generating the control mesh.
(18) Merges two surface meshes together to form a watertight mesh.In order to form a watertight mesh, vertices on the intersection must be the same on both meshes.
(19) Starts the decimation of a selected mesh.The reduction factor can be set by the user.

Figure 2 :
Figure 2: A common example of a non-manifold topology in geological modelling; multiple faces share one edge.Non-manifold vertices (green) and edge (red).

Figure 3 :
Figure 3: Splitting step of Loop scheme; (a) The control mesh; (b) Splitting each triangle into four by inserting new vertices (yellow vertices) on the mid of each edge.

Figure 4 :
Figure 4: Averaging step of generating the smooth surface by Loop subdivision scheme; (a) Updating the position of the existing vertices; (b) Updating the position of the new midpoint vertices.
The rule for updating the position of the   one edge dart smooth subdivision rules (section 2.1.1)Two edges (   and   ) where the   and   are two adjacent vertices connected to  by corner   =  (corner rule) 5c represents the effect of the resistance of three sharp edges to the smoothing procedure (black edges, with crease sharpness values equal to one).

Figure 5 :
Figure 5: Generating creases on a mesh by applying three times subdivision surface algorithm; (a) Control mesh (b); All edges of the control mesh are smooth edges (red edges); (c) Four edges are crease edges (blue edges), and nine edges are smooth (red edges).
Fig.6represents the surface intersections as a common example of non-manifold topology in geological modelling where multiple faces are shared by one edge ( e.g., the intersection of different faults or intersections between a horizon and a fault).

Figure 6 :
Figure 6: (a) A schematic representation of the intersection between a fault and a layer; (b) Mesh structure of the surface intersection includes the simple singular vertices (brown circles).

Figure 8 :
Figure 8: General workflow of subdivision surface fitting with an example of an anticlinal structure.

Figure 9 :
Figure 9: Synthesized and conceptualized approach for surface fitting with subdivision surfaces.

Figure 10
Figure 10 represents the GUI window of PySubdiv containing the anticline model (two non-watertight surfaces) which is generated by Gempy (de la Varga et al. 2019).It is worth mentioning that watertight or non-watertight meshes can be imported as the input mesh.All of the elements inside the GUI are explained in the appendix of this paper.

Figure 10 :
Figure 10: The GUI of PySubdiv allows the user to construct the control mesh for the (geological) model.The important GUI elements are labelled and explained in the following.

Figure 11 :
Figure 11: (a) Control points are generated based on the vertices of the input mesh; (b) Triangulated control mesh is constructed based on the control points and the Delaunay triangulation.

Figure 12 :
Figure 12: (a) Stitching and merging the different parts of the control mesh to make it watertight ; (b) Final watertight control mesh

Figure 13 :
Figure 13: Grid-based representation of part of URG model which contains 616464 individual nodes and 10 different volumes generated based on the data of Freymark et al. (2020).

Figure 14 :
Figure 14: (a) Initial and unfitted control mesh of the URG model.Surfaces are coloured concerning the different geological units.The surfaces of the Lower Crystalline Crust and the Lithosphere mantel are hidden by the boundary surface.The 832 control points are coloured in red; (b) Approximately representation of the 3D side view of the model along the profile  ′ .

Figure 15 :
Figure 15: (a) Distance map of the subdivided surfaces compared to the input mesh (Fig. 13).Red areas indicate high deviation from the original input mesh while blue areas indicate low deviation; (b) Exploded view on the distance map of the four extracted sub-horizontal layers.Distances are scaled for each layer individually to emphasize the areas of high and low errors independent of the maximal global error.

Figure 16 :
Figure 16: Distinguishing different parts of the surface by adding control points to avoid failure in the reconstruction process.