(BUGFIX) supports for discrete interpolation

were being reused, which may have caused unexpected behavior...

Author: lachlangrose

LoopStructural/modelling/core/geological_model.py

@@ -158,6 +158,10 @@ def __init__(self, origin, maximum, rescale=True, nsteps=(40, 40, 40),
         self.bounding_box /= self.scale_factor
         self.support = {}
         self.reuse_supports = reuse_supports
+        if reuse_supports:
+            logger.warning("Supports are shared between geological features \n"
+                                "this may cause unexpected behaviour and should only\n"
+                                "be use by advanced users")
         logger.info('Reusing interpolation supports: {}'.format(self.reuse_supports))
         self.stratigraphic_column = None
         self.parameters = {'features': [], 'model': {'bounding_box': self.origin.tolist() + self.maximum.tolist(),
@@ -480,18 +484,22 @@ def get_interpolator(self, interpolatortype='PLI', nelements=1e5,
             # number of steps is the length of the box / step vector
             nsteps = np.ceil((bb[1, :] - bb[0, :]) / step_vector).astype(int)
             # create a structured grid using the origin and number of steps
-            mesh_id = 'mesh_{}'.format(nelements)
-            mesh = self.support.get(mesh_id,
-                                    TetMesh(origin=bb[0, :], nsteps=nsteps,
-                                            step_vector=step_vector))
-            if mesh_id not in self.support:
-                self.support[mesh_id] = mesh
+            if reuse_supports:
+                mesh_id = 'mesh_{}'.format(nelements)
+                mesh = self.support.get(mesh_id,
+                                        TetMesh(origin=bb[0, :], nsteps=nsteps,
+                                                step_vector=step_vector))
+                if mesh_id not in self.support:
+                    self.support[mesh_id] = mesh
+            else:
+                mesh = TetMesh(origin=bb[0, :], nsteps=nsteps, step_vector=step_vector)
             logger.info("Creating regular tetrahedron mesh with %i elements \n"
                         "for modelling using PLI" % (mesh.ntetra))
 
             return PLI(mesh)
 
         if interpolatortype == 'FDI':
+
             # find the volume of one element
             if element_volume is None:
                 element_volume = box_vol / nelements
@@ -504,12 +512,15 @@ def get_interpolator(self, interpolatortype='PLI', nelements=1e5,
                 logger.error("Cannot create interpolator: number of steps is too small")
                 return None
             # create a structured grid using the origin and number of steps
-            grid_id = 'grid_{}'.format(nelements)
-            grid = self.support.get(grid_id, StructuredGrid(origin=bb[0, :],
+            if self.reuse_supports:
+                grid_id = 'grid_{}'.format(nelements)
+                grid = self.support.get(grid_id, StructuredGrid(origin=bb[0, :],
                                                             nsteps=nsteps,
                                                             step_vector=step_vector))
-            if grid_id not in self.support:
-                self.support[grid_id] = grid
+                if grid_id not in self.support:
+                    self.support[grid_id] = grid
+            else:
+                grid = StructuredGrid(origin=bb[0, :], nsteps=nsteps,step_vector=step_vector)
             logger.info("Creating regular grid with %i elements \n"
                         "for modelling using FDI" % grid.n_elements)
             return FDI(grid)

LoopStructural/utils/bounding_box.py

@@ -0,0 +1,10 @@
+# class BoundingBox:
+#     def __init__(self,origin,maximum):
+#         self.origin = origin
+#         self.maximum = maximum
+#         self.name_map = {'xmin':(0,0),'ymin':(0,1),'zmin':(0,2),'xmax':(1,0),'ymax':(1,1),'zmax':(1,2)
+#                          'lower':(0,2),'upper':(1,2),
+#                          'minx':(0,0),'miny':(0,1),'minz':(0,2),'maxx':(1,0),'maxy':(1,1),'maxz':(1,2)}
+#     def get_value(self,name):
+#         ix,iy = self.name_map.get(name,(-1,-1))            
+