Skip to content

Commit b34995c

Browse files
committed
[ENH] Better implementations of octree/dense grid interaction
1 parent ff0f126 commit b34995c

File tree

5 files changed

+133
-62
lines changed

5 files changed

+133
-62
lines changed
Lines changed: 17 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,22 @@
11
# Refactoring grid to add transform
22

33
## TODO:
4-
- [ ] Rotation of pyvista grid
5-
- [ ] Rotation of pyvista meshes
4+
- [x] Rotation of pyvista grid
5+
- [x] Rotation of pyvista meshes
6+
- [ ] Run something with octrees again
7+
- [ ] **We do need settings presets**
8+
69

10+
## Notes:
11+
### Optimizing interpolation June 2024
12+
13+
- Create_grad_kernel is call 27 times
14+
15+
### Other:
16+
{'_cached_pivot': [200, 200, 0],
17+
'_is_default_transform': False,
18+
'position': array([0., 0., 0.]),
19+
'rotation': array([ 0., -0., -45.]),
20+
'scale': array([1., 1., 1.])}
21+
722
## Questions:

gempy/API/compute_API.py

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@
1010
from .grid_API import set_custom_grid
1111
from ..core.data.gempy_engine_config import GemPyEngineConfig
1212
from ..core.data.geo_model import GeoModel
13+
from ..modules.data_manipulation.engine_factory import interpolation_input_from_structural_frame
1314
from ..optional_dependencies import require_gempy_legacy
1415

1516

@@ -39,7 +40,7 @@ def compute_model(gempy_model: GeoModel, engine_config: Optional[GemPyEngineConf
3940
)
4041

4142
# TODO: To decide what to do with this.
42-
interpolation_input = gempy_model.interpolation_input_copy
43+
interpolation_input = interpolation_input_from_structural_frame(gempy_model)
4344
gempy_model.taped_interpolation_input = interpolation_input # * This is used for gradient tape
4445

4546
gempy_model.solutions = gempy_engine.compute_model(

gempy/API/initialization_API.py

Lines changed: 11 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66

77
from gempy.API.io_API import read_surface_points, read_orientations
88
from gempy_engine.core.data import InterpolationOptions
9+
from ..core.data.grid_modules import RegularGrid
910
from ..optional_dependencies import require_subsurface
1011
from ..core.data import StructuralElement
1112
from ..core.data.geo_model import GeoModel
@@ -44,10 +45,12 @@ def create_geomodel(
4445
"""
4546
# init resolutions well
4647
if resolution is None:
47-
grid: Grid = Grid(extent=extent)
48-
grid.octree_levels = refinement
48+
grid: Grid = Grid.init_octree_grid(
49+
extent=extent,
50+
octree_levels=refinement
51+
)
4952
else:
50-
grid = Grid(
53+
grid: Grid = Grid.init_dense_grid(
5154
extent=extent,
5255
resolution=resolution
5356
)
@@ -92,18 +95,18 @@ def structural_elements_from_borehole_set(
9295
Returns:
9396
list[StructuralElement]: A list of StructuralElement instances.
9497
"""
95-
98+
9699
ss = require_subsurface()
97100
borehole_set: ss.core.geological_formats.BoreholeSet
98-
101+
99102
elements = []
100103
component_lith: dict[Hashable, np.ndarray] = borehole_set.get_bottom_coords_for_each_lith()
101-
104+
102105
for name, properties in elements_dict.items():
103106
top_coordinates = component_lith.get(properties['id'])
104107
if top_coordinates is None:
105108
raise ValueError(f"Top lithology {properties['id']} not found in borehole set.")
106-
109+
107110
element = StructuralElement(
108111
name=name,
109112
id=properties['id'],
@@ -119,7 +122,7 @@ def structural_elements_from_borehole_set(
119122
)
120123
elements.append(element)
121124
# Reverse the list to have the oldest rocks at the bottom
122-
125+
123126
return elements
124127

125128

gempy/core/data/grid.py

Lines changed: 52 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
from typing import Optional
77

88
from gempy_engine.core.data.centered_grid import CenteredGrid
9+
from gempy_engine.core.data.options import EvaluationOptions
910
from gempy_engine.core.data.transforms import Transform
1011
from .grid_modules import RegularGrid, CustomGrid, Sections
1112
from .grid_modules.topography import Topography
@@ -23,7 +24,6 @@ class GridTypes(enum.Flag):
2324
NONE = 2 ** 10
2425

2526
# ? What should we do with the extent?
26-
_extent: Optional[np.ndarray] # * Model extent should be cross grid
2727

2828
_octree_grid: Optional[RegularGrid] = None
2929
_dense_grid: Optional[RegularGrid] = None
@@ -41,11 +41,25 @@ class GridTypes(enum.Flag):
4141
_octree_levels: int = -1
4242

4343
def __init__(self, extent=None, resolution=None):
44-
self.extent = extent
4544
# Init basic grid empty
4645
if extent is not None and resolution is not None:
4746
self.dense_grid = RegularGrid(extent, resolution)
4847

48+
@classmethod
49+
def init_octree_grid(cls, extent, octree_levels):
50+
grid = cls()
51+
grid._octree_grid = RegularGrid(
52+
extent=extent,
53+
resolution=np.array([2 ** octree_levels] * 3),
54+
)
55+
grid.active_grids |= grid.GridTypes.OCTREE
56+
grid._update_values()
57+
return grid
58+
59+
@classmethod
60+
def init_dense_grid(cls, extent, resolution):
61+
return cls(extent, resolution)
62+
4963
def __str__(self):
5064
active_grid_types_str = [g_type for g_type in self.GridTypes if self.active_grids & g_type]
5165

@@ -54,44 +68,32 @@ def __str__(self):
5468
grid_summary_str = "\n".join(grid_summary)
5569
return f"Grid Object:\n{grid_summary_str}"
5670

57-
58-
5971
@property
6072
def transform(self) -> Transform:
61-
if self._transform is None:
62-
if self.dense_grid is not None:
63-
return self.dense_grid.transform
64-
elif self.octree_grid is not None:
65-
return self.octree_grid.transform
66-
else:
67-
return Transform.init_neutral()
68-
return self._transform
73+
if self.dense_grid is not None:
74+
return self.dense_grid.transform
75+
elif self.octree_grid is not None:
76+
return self.octree_grid.transform
77+
else:
78+
return Transform.init_neutral()
6979

7080
@transform.setter
7181
def transform(self, value: Transform):
7282
self._transform = value
73-
83+
7484
@property
7585
def extent(self):
76-
if self._extent is None:
77-
# Try to get the extent from the dense or octree grid if those are also none raise an error
78-
if self.dense_grid is not None:
79-
return self.dense_grid.extent
80-
elif self.octree_grid is not None:
81-
return self.octree_grid.extent
82-
else:
83-
raise AttributeError('Extent is not defined')
86+
if self.dense_grid is not None:
87+
return self.dense_grid.extent
88+
elif self.octree_grid is not None:
89+
return self.octree_grid.extent
8490
else:
85-
return self._extent
91+
raise AttributeError('Extent is not defined')
8692

87-
@extent.setter
88-
def extent(self, value):
89-
self._extent = value
90-
9193
@property
9294
def corner_min(self):
9395
return self.extent[::2]
94-
96+
9597
@property
9698
def corner_max(self):
9799
return self.extent[1::2]
@@ -135,6 +137,18 @@ def octree_grid(self):
135137

136138
@octree_grid.setter
137139
def octree_grid(self, value):
140+
raise AttributeError('Octree grid is not allowed to be set directly. Use init_octree_grid instead')
141+
142+
def set_octree_grid(self, value: RegularGrid, evaluation_options: EvaluationOptions):
143+
regular_grid_resolution = value.resolution
144+
# Check all directions has the same res
145+
if not np.all(regular_grid_resolution == regular_grid_resolution[0]):
146+
raise AttributeError('Octree resolution must be isotropic')
147+
octree_levels = int(np.log2(regular_grid_resolution[0]))
148+
# Check if octrree levels are the same
149+
if octree_levels != evaluation_options.number_octree_levels:
150+
raise AttributeError('Regular grid resolution does not match octree levels. Resolution must be 2^n')
151+
138152
self._octree_grid = value
139153
self.active_grids |= self.GridTypes.OCTREE
140154
self._update_values()
@@ -181,8 +195,10 @@ def centered_grid(self, value):
181195

182196
@property
183197
def regular_grid(self):
184-
warnings.warn('This property is deprecated. Use the dense_grid or octree_grid instead', DeprecationWarning)
185-
if self.dense_grid is not None and self.octree_grid is not None:
198+
dense_grid_exists_and_active = self.dense_grid is not None and self.GridTypes.DENSE in self.active_grids
199+
octree_grid_exists_and_active = self.octree_grid is not None and self.GridTypes.OCTREE in self.active_grids
200+
201+
if dense_grid_exists_and_active and octree_grid_exists_and_active:
186202
raise AttributeError('Both dense_grid and octree_grid are active. This is not possible.')
187203
elif self.dense_grid is not None:
188204
return self.dense_grid
@@ -191,41 +207,27 @@ def regular_grid(self):
191207
else:
192208
return None
193209

194-
@regular_grid.setter
195-
def regular_grid(self, value):
196-
raise AttributeError('This property is deprecated. Use the dense_grid or octree_grid instead')
197-
198-
@property
199-
def octree_levels(self):
200-
return self._octree_levels
201-
202-
@property
203-
def octree_levels(self):
204-
return self._octree_levels
205-
206-
@octree_levels.setter
207-
def octree_levels(self, value):
208-
self._octree_levels = value
209-
self.octree_grid = RegularGrid(
210-
extent=self.extent,
211-
resolution=np.array([2 ** value] * 3),
212-
)
213-
self.active_grids |= self.GridTypes.OCTREE
214-
210+
# noinspection t
215211
def _update_values(self):
216212
values = []
217213

218214
if self.GridTypes.OCTREE in self.active_grids:
215+
if self.octree_grid is None: raise AttributeError('Octree grid is active but not defined')
219216
values.append(self.octree_grid.values)
220217
if self.GridTypes.DENSE in self.active_grids:
218+
if self.dense_grid is None: raise AttributeError('Dense grid is active but not defined')
221219
values.append(self.dense_grid.values)
222220
if self.GridTypes.CUSTOM in self.active_grids:
221+
if self.custom_grid is None: raise AttributeError('Custom grid is active but not defined')
223222
values.append(self.custom_grid.values)
224223
if self.GridTypes.TOPOGRAPHY in self.active_grids:
224+
if self.topography is None: raise AttributeError('Topography grid is active but not defined')
225225
values.append(self.topography.values)
226226
if self.GridTypes.SECTIONS in self.active_grids:
227+
if self.sections is None: raise AttributeError('Sections grid is active but not defined')
227228
values.append(self.sections.values)
228229
if self.GridTypes.CENTERED in self.active_grids:
230+
if self.centered_grid is None: raise AttributeError('Centered grid is active but not defined')
229231
values.append(self.centered_grid.values)
230232

231233
self.values = np.concatenate(values)

test/test_modules/.benchmarks/test_transformed_space.py

Lines changed: 51 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -55,7 +55,7 @@ def test_plot_transformed_data_including_grid_transform():
5555
model.grid.dense_grid = regular_grid
5656

5757
gp.compute_model(model)
58-
58+
5959
if PLOT:
6060
gpv = require_gempy_viewer()
6161

@@ -80,3 +80,53 @@ def test_plot_transformed_data_including_grid_transform():
8080
'arrow_size': .01
8181
}
8282
)
83+
84+
85+
def test_plot_transformed_data_including_grid_transform_octree():
86+
model = gp.generate_example_model(ExampleModel.ANTICLINE, compute_model=False)
87+
88+
i = 4
89+
# Calculate point_y_axis
90+
regular_grid = gp.data.grid.RegularGrid.from_corners_box(
91+
pivot=(200, 200),
92+
point_x_axis=(800, 800),
93+
distance_point3=1000,
94+
zmin=model.grid.extent[4],
95+
zmax=model.grid.extent[5],
96+
resolution=np.array([2 ** i] * 3),
97+
plot=True
98+
)
99+
100+
# model.grid.active_grids ^= gp.data.grid.Grid.GridTypes.DENSE
101+
102+
options = model.interpolation_options
103+
options.number_octree_levels = 4
104+
105+
model.grid.set_octree_grid(regular_grid, model.interpolation_options.evaluation_options)
106+
107+
gp.compute_model(model)
108+
109+
if PLOT:
110+
gpv = require_gempy_viewer()
111+
112+
gpv.plot_3d(
113+
model,
114+
image=True,
115+
transformed_data=True,
116+
show_boundaries=True,
117+
show_lith=True,
118+
kwargs_plot_data={
119+
'arrow_size': .01
120+
}
121+
)
122+
123+
gpv.plot_3d(
124+
model,
125+
image=False,
126+
transformed_data=False,
127+
show_boundaries=True,
128+
show_lith=True,
129+
kwargs_plot_data={
130+
'arrow_size': .01
131+
}
132+
)

0 commit comments

Comments
 (0)