Skip to content

Commit 19342da

Browse files
benk-mirabenk-mira
committed
refactor
1 parent 7a36da5 commit 19342da

12 files changed

Lines changed: 356 additions & 382 deletions

File tree

‎simpeg_drivers/plate_simulation/driver.py‎

Lines changed: 30 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -16,11 +16,11 @@
1616
import numpy as np
1717
from dask.distributed import Client
1818
from geoapps_utils.base import Driver, get_logger
19+
from geoapps_utils.modelling.plates import Plate
1920
from geoapps_utils.utils.transformations import azimuth_to_unit_vector
2021
from geoh5py.data import FloatData, ReferencedData
2122
from geoh5py.objects import Octree, Points, Surface
2223
from geoh5py.shared.utils import fetch_active_workspace
23-
from grid_apps.octree_creation.driver import OctreeDriver
2424

2525
from simpeg_drivers.driver import (
2626
InversionDriver,
@@ -30,9 +30,9 @@
3030
)
3131
from simpeg_drivers.options import BaseForwardOptions, ModelTypeEnum
3232
from simpeg_drivers.plate_simulation.models.events import Anomaly, Erosion, Overburden
33-
from simpeg_drivers.plate_simulation.models.parametric import Plate
3433
from simpeg_drivers.plate_simulation.models.series import DikeSwarm, Geology
3534
from simpeg_drivers.plate_simulation.options import PlateSimulationOptions
35+
from simpeg_drivers.utils.synthetics.meshes import get_octree_mesh
3636
from simpeg_drivers.utils.utils import validate_out_group
3737

3838

@@ -141,7 +141,7 @@ def plates(self) -> list[Plate]:
141141
depth_offset=-1 * offset,
142142
)
143143
plate = Plate(
144-
self.params.model.plate.model_copy(
144+
self.params.model.plate.geometry.model_copy(
145145
update={
146146
"easting": center[0],
147147
"northing": center[1],
@@ -185,13 +185,21 @@ def make_mesh(self) -> Octree:
185185
"""
186186

187187
logger.info("making the mesh...")
188-
octree_params = self.params.mesh.octree_params(
189-
self.survey,
190-
self.simulation_parameters.active_cells.topography_object,
191-
[p.surface.copy(parent=self._out_group) for p in self.plates],
192-
)
193-
octree_driver = OctreeDriver(octree_params)
194-
mesh = octree_driver.run()
188+
with fetch_active_workspace(self.params.geoh5, mode="r+") as geoh5:
189+
surfaces = [p.surface(geoh5) for p in self.plates]
190+
mesh = get_octree_mesh(
191+
opts=self.params.mesh,
192+
survey=self.survey,
193+
topography=self.simulation_parameters.active_cells.topography_object,
194+
plates=surfaces,
195+
)
196+
# octree_params = self.params.mesh.octree_params(
197+
# self.survey,
198+
# self.simulation_parameters.active_cells.topography_object,
199+
# [p.surface.copy(parent=self._out_group) for p in self.plates],
200+
# )
201+
# octree_driver = OctreeDriver(octree_params)
202+
# mesh = octree_driver.run()
195203
mesh.parent = self._out_group
196204

197205
return mesh
@@ -208,7 +216,10 @@ def make_model(self) -> FloatData:
208216
)
209217

210218
dikes = DikeSwarm(
211-
[Anomaly(plate, plate.params.plate_property) for plate in self.plates],
219+
[
220+
Anomaly(plate, self.params.model.plate.plate_property)
221+
for plate in self.plates
222+
],
212223
name="plates",
213224
)
214225

@@ -278,21 +289,19 @@ def replicate(
278289
plates = []
279290
for i in range(number):
280291
center = (
281-
np.r_[plate.params.geometry.origin]
292+
np.r_[plate.params.origin]
282293
+ azimuth_to_unit_vector(azimuth) * offsets[i]
283294
)
284-
new_geometry = plate.params.geometry.model_copy(
285-
update={
286-
"easting": center[0],
287-
"northing": center[1],
288-
"elevation": center[2],
289-
}
290-
)
291295
new_plate = Plate(
292-
plate.params.model_copy(update={"geometry": new_geometry})
296+
plate.params.model_copy(
297+
update={
298+
"easting": center[0],
299+
"northing": center[1],
300+
"elevation": center[2],
301+
}
302+
)
293303
)
294304

295-
new_plate.params.name = f"{plate.params.name} offset {i + 1}"
296305
plates.append(new_plate)
297306

298307
return plates

‎simpeg_drivers/plate_simulation/models/parametric.py‎

Lines changed: 115 additions & 125 deletions
Original file line numberDiff line numberDiff line change
@@ -13,20 +13,10 @@
1313
from abc import ABC, abstractmethod
1414

1515
import numpy as np
16-
from geoapps_utils.modelling.plates import PlateModel, inside_plate
17-
from geoapps_utils.utils.transformations import (
18-
rotate_points,
19-
rotate_xyz,
20-
x_rotation_matrix,
21-
z_rotation_matrix,
22-
)
2316
from geoh5py.objects import Octree, Surface
24-
from geoh5py.shared.utils import fetch_active_workspace
25-
from geoh5py.workspace import Workspace
2617
from trimesh import Trimesh
2718
from trimesh.proximity import ProximityQuery
2819

29-
from simpeg_drivers.plate_simulation.models.options import PlateOptions
3020
from simpeg_drivers.utils.utils import active_from_xyz
3121

3222

@@ -57,121 +47,121 @@ def mask(self, mesh: Octree) -> np.ndarray:
5747
"""
5848

5949

60-
class Plate(Parametric):
61-
"""
62-
Define a rotated rectangular block in 3D space
63-
64-
:param params: Parameters describing the plate.
65-
:param surface: Surface object representing the plate.
66-
"""
67-
68-
def __init__(
69-
self,
70-
params: PlateOptions,
71-
workspace: Workspace | None = None,
72-
):
73-
self.params = params
74-
self._workspace = workspace
75-
super().__init__(self._create_surface())
76-
77-
def _create_surface(self) -> Surface:
78-
"""
79-
Create a surface object from a plate object.
80-
81-
:param workspace: Workspace object to create the surface in.
82-
:param out_group: Output group to store the surface.
83-
"""
84-
with fetch_active_workspace(self.workspace) as ws:
85-
surface = Surface.create(
86-
ws,
87-
vertices=self.vertices,
88-
cells=self.triangles,
89-
name=self.params.name,
90-
)
91-
92-
return surface
93-
94-
@property
95-
def surface(self):
96-
"""
97-
A surface object representing the plate.
98-
"""
99-
return self._surface
100-
101-
@property
102-
def triangles(self) -> np.ndarray:
103-
"""Triangulation of the block."""
104-
return np.vstack(
105-
[
106-
[0, 2, 1],
107-
[1, 2, 3],
108-
[0, 1, 4],
109-
[4, 1, 5],
110-
[1, 3, 5],
111-
[5, 3, 7],
112-
[2, 6, 3],
113-
[3, 6, 7],
114-
[0, 4, 2],
115-
[2, 4, 6],
116-
[4, 5, 6],
117-
[6, 5, 7],
118-
]
119-
)
120-
121-
@property
122-
def vertices(self) -> np.ndarray:
123-
"""Vertices for triangulation of a rectangular prism in 3D space."""
124-
125-
u_1 = self.params.geometry.origin[0] - (
126-
self.params.geometry.strike_length / 2.0
127-
)
128-
u_2 = self.params.geometry.origin[0] + (
129-
self.params.geometry.strike_length / 2.0
130-
)
131-
v_1 = self.params.geometry.origin[1] - (self.params.geometry.dip_length / 2.0)
132-
v_2 = self.params.geometry.origin[1] + (self.params.geometry.dip_length / 2.0)
133-
w_1 = self.params.geometry.origin[2] - (self.params.geometry.width / 2.0)
134-
w_2 = self.params.geometry.origin[2] + (self.params.geometry.width / 2.0)
135-
136-
vertices = np.array(
137-
[
138-
[u_1, v_1, w_1],
139-
[u_2, v_1, w_1],
140-
[u_1, v_2, w_1],
141-
[u_2, v_2, w_1],
142-
[u_1, v_1, w_2],
143-
[u_2, v_1, w_2],
144-
[u_1, v_2, w_2],
145-
[u_2, v_2, w_2],
146-
]
147-
)
148-
149-
return self._rotate(vertices)
150-
151-
@property
152-
def workspace(self) -> Workspace:
153-
if self._workspace is None:
154-
self._workspace = Workspace()
155-
156-
return self._workspace
157-
158-
def _rotate(self, vertices: np.ndarray) -> np.ndarray:
159-
"""Rotate vertices and adjust for reference point."""
160-
theta = -1 * self.params.geometry.direction
161-
phi = -1 * self.params.geometry.dip
162-
rotated_vertices = rotate_xyz(vertices, self.params.geometry.origin, theta, phi)
163-
164-
return rotated_vertices
165-
166-
def mask(self, mesh: Octree) -> np.ndarray:
167-
rotations = [
168-
z_rotation_matrix(np.deg2rad(self.params.geometry.direction)),
169-
x_rotation_matrix(np.deg2rad(self.params.geometry.dip)),
170-
]
171-
rotated_centers = rotate_points(
172-
mesh.centroids, origin=self.params.geometry.origin, rotations=rotations
173-
)
174-
return inside_plate(rotated_centers, self.params.geometry)
50+
# class Plate(Parametric):
51+
# """
52+
# Define a rotated rectangular block in 3D space
53+
#
54+
# :param params: Parameters describing the plate.
55+
# :param surface: Surface object representing the plate.
56+
# """
57+
#
58+
# def __init__(
59+
# self,
60+
# params: PlateOptions,
61+
# workspace: Workspace | None = None,
62+
# ):
63+
# self.params = params
64+
# self._workspace = workspace
65+
# super().__init__(self._create_surface())
66+
#
67+
# def _create_surface(self) -> Surface:
68+
# """
69+
# Create a surface object from a plate object.
70+
#
71+
# :param workspace: Workspace object to create the surface in.
72+
# :param out_group: Output group to store the surface.
73+
# """
74+
# with fetch_active_workspace(self.workspace) as ws:
75+
# surface = Surface.create(
76+
# ws,
77+
# vertices=self.vertices,
78+
# cells=self.triangles,
79+
# name=self.params.name,
80+
# )
81+
#
82+
# return surface
83+
#
84+
# @property
85+
# def surface(self):
86+
# """
87+
# A surface object representing the plate.
88+
# """
89+
# return self._surface
90+
#
91+
# @property
92+
# def triangles(self) -> np.ndarray:
93+
# """Triangulation of the block."""
94+
# return np.vstack(
95+
# [
96+
# [0, 2, 1],
97+
# [1, 2, 3],
98+
# [0, 1, 4],
99+
# [4, 1, 5],
100+
# [1, 3, 5],
101+
# [5, 3, 7],
102+
# [2, 6, 3],
103+
# [3, 6, 7],
104+
# [0, 4, 2],
105+
# [2, 4, 6],
106+
# [4, 5, 6],
107+
# [6, 5, 7],
108+
# ]
109+
# )
110+
#
111+
# @property
112+
# def vertices(self) -> np.ndarray:
113+
# """Vertices for triangulation of a rectangular prism in 3D space."""
114+
#
115+
# u_1 = self.params.geometry.origin[0] - (
116+
# self.params.geometry.strike_length / 2.0
117+
# )
118+
# u_2 = self.params.geometry.origin[0] + (
119+
# self.params.geometry.strike_length / 2.0
120+
# )
121+
# v_1 = self.params.geometry.origin[1] - (self.params.geometry.dip_length / 2.0)
122+
# v_2 = self.params.geometry.origin[1] + (self.params.geometry.dip_length / 2.0)
123+
# w_1 = self.params.geometry.origin[2] - (self.params.geometry.width / 2.0)
124+
# w_2 = self.params.geometry.origin[2] + (self.params.geometry.width / 2.0)
125+
#
126+
# vertices = np.array(
127+
# [
128+
# [u_1, v_1, w_1],
129+
# [u_2, v_1, w_1],
130+
# [u_1, v_2, w_1],
131+
# [u_2, v_2, w_1],
132+
# [u_1, v_1, w_2],
133+
# [u_2, v_1, w_2],
134+
# [u_1, v_2, w_2],
135+
# [u_2, v_2, w_2],
136+
# ]
137+
# )
138+
#
139+
# return self._rotate(vertices)
140+
#
141+
# @property
142+
# def workspace(self) -> Workspace:
143+
# if self._workspace is None:
144+
# self._workspace = Workspace()
145+
#
146+
# return self._workspace
147+
#
148+
# def _rotate(self, vertices: np.ndarray) -> np.ndarray:
149+
# """Rotate vertices and adjust for reference point."""
150+
# theta = -1 * self.params.geometry.direction
151+
# phi = -1 * self.params.geometry.dip
152+
# rotated_vertices = rotate_xyz(vertices, self.params.geometry.origin, theta, phi)
153+
#
154+
# return rotated_vertices
155+
#
156+
# def mask(self, mesh: Octree) -> np.ndarray:
157+
# rotations = [
158+
# z_rotation_matrix(np.deg2rad(self.params.geometry.direction)),
159+
# x_rotation_matrix(np.deg2rad(self.params.geometry.dip)),
160+
# ]
161+
# rotated_centers = rotate_points(
162+
# mesh.centroids, origin=self.params.geometry.origin, rotations=rotations
163+
# )
164+
# return inside_plate(rotated_centers, self.params.geometry)
175165

176166

177167
class Body(Parametric):

0 commit comments

Comments
 (0)