# -*- coding: utf-8 -*-
"""
@author: Rochlitz.R
"""
from custEM.meshgen import meshgen_utils as mu
from custEM.meshgen.meshgen_tools import BlankWorld
import numpy as np
import pygimli as pg
[docs]
class PrismWorld:
"""
Specific class for automated generation of pseudo-2D / 3D inversion
meshes discretized by prisms.
"""
def __init__(self, name, x_extent, x_reduction=100., y_depth=300.,
z_depth=300., prism_quality=32., prism_area=1e4, n_prisms=50.,
orthogonal_txs=None, txs=[], surface_rxs=[],
**blankworld_kwargs):
"""
Inititialize paramters relevant for prism incorporation and add
everything to a mesh generated with a *BlankWorld* instance. All
*BlankWorld* keyword arguments are passed.
"""
self.x_extent = x_extent
self.x_reduction = x_reduction
self.y_depth = y_depth
self.z_depth = z_depth
self.n_prisms = n_prisms
self.x_vec = np.linspace(x_extent[0], x_extent[1], n_prisms + 1)
self.txs = [tx for tx in txs]
self.surface_rx = surface_rxs
if orthogonal_txs is None:
self.orthogonal_txs = [False] * len(txs)
else:
self.orthogonal_txs = orthogonal_txs
self.topo = None
if 'topo' in blankworld_kwargs:
self.topo = blankworld_kwargs['topo']
self.PrismWorld = BlankWorld(name=name,
preserve_edges=True,
**blankworld_kwargs)
self.create_2d_mesh(name, prism_quality, prism_area)
self.create_surface_rectangles()
self.PrismWorld.build_surface(insert_loops=self.rectangles,
insert_line_tx=self.txs,
insert_paths=self.surface_rx)
self.PrismWorld.build_halfspace_mesh()
back_nodes, front_nodes = [], []
for ci in range(len(self.prism_coords)):
back_nodes.append(self.PrismWorld.Omega.createNodeWithCheck(
[self.prism_coords[ci, 0], -self.y_depth,
self.prism_coords[ci, 2]]))
front_nodes.append(self.PrismWorld.Omega.createNodeWithCheck(
[self.prism_coords[ci, 0], self.y_depth,
self.prism_coords[ci, 2]]))
self.add_triangles(back_nodes, front_nodes)
self.add_prisms(back_nodes, front_nodes)
[docs]
def create_2d_mesh(self, name, prism_quality, prism_area):
"""
Create and store 2D mesh in x-z slice with triangle, which will be
the basis for elongating the triangles to prisms in y-direction later
on.
"""
# shift prism edge at surface on tx if a completely orthogonal tx
# is close to a prism edge, use 20 % of prism distace to identify
dist = np.diff(self.x_vec)[0]
for ti, tx in enumerate(self.txs):
if self.orthogonal_txs[ti]:
for xi, x in enumerate(self.x_vec):
if np.abs(x - tx[0, 0]) < 0.2 * dist:
self.x_vec[xi] = tx[0, 0]
# check for rx intersections and resolve by shifting rx refinement
# edges
for ri, rx in enumerate(self.surface_rx):
for x in self.x_vec:
if (x > np.min(rx[:, 0]) and x < np.max(rx[:, 0])):
self.surface_rx[ri][:, 0] += x - np.min(rx[:, 0])
if self.topo is None:
z = np.zeros(len(self.x_vec))
else:
z = np.zeros((len(self.x_vec), 3))
z[:, 0] = self.x_vec
z = self.PrismWorld.add_topo(z)[:, 2]
self.z_depth -= np.mean(z)
frame = pg.meshtools.createPolygon(
[[self.x_extent[1] - self.x_reduction, -self.z_depth]] +
[[self.x_extent[0] + self.x_reduction, -self.z_depth]] +
[[self.x_vec[xi], z[xi]] for xi in range(self.n_prisms + 1)],
isClosed=True)
[frame.createNodeWithCheck([self.x_vec[xi], z[xi]]) for
xi in range(self.n_prisms)]
self.xzmesh = pg.meshtools.createMesh(
frame, quality=prism_quality, area=prism_area)
self.xzmesh.save(
self.PrismWorld.w_dir + '/' + name + '.bms')
self.xzmesh.setDimension(3)
self.xzmesh.swapCoordinates(1, 2)
self.prism_coords = self.xzmesh.positions().array()
self.surf_nodes = np.arange(2,2+len(z))
[docs]
def create_surface_rectangles(self):
"""
Build rectangular frames of the prisms at the surface to be included
in the surface mesh later on for avoiding intersections. Check also
if transmitter paths are intersecting these rectangles, and if yes,
find the intersections and add the corresponding points to the frames.
"""
tmp_tx = [np.array(tx) for tx in self.txs]
self.rectangles = []
self.rectangles.append(mu.loop_r(
[self.x_vec[0], -self.y_depth],
[self.x_vec[1], self.y_depth], nx=1, ny=1))
for xi in range(self.n_prisms - 1):
self.rectangles.append(mu.loop_r(
[self.x_vec[xi+1], -self.y_depth],
[self.x_vec[xi+2], self.y_depth], nx=1, ny=1))
for ti in range(len(tmp_tx)):
for si in range(len(tmp_tx[ti]) - 1):
isect = mu.get_intersect(tmp_tx[ti][si, :2],
tmp_tx[ti][si + 1, :2],
self.rectangles[xi][1, :2],
self.rectangles[xi][2, :2])
if isect is not None:
self.rectangles[xi] = np.insert(
self.rectangles[xi], 2,
np.append(isect, [0.]), 0)
self.rectangles[xi + 1] = np.insert(
self.rectangles[xi + 1],
len(self.rectangles[xi + 1]),
np.append(isect, [0.]), 0)
self.txs[ti] = np.vstack((self.txs[ti],
np.append(isect, [0.])))
for ti in range(len(self.txs)):
self.txs[ti] = self.txs[ti][np.argsort(self.txs[ti][:, 0])]
[docs]
def add_triangles(self, back_nodes, front_nodes):
"""
Add front and back triangle faces of the prisms. All the if statements
are required to account for potential nodes on the x-directed boundary
edges of the rectangular prism outcrops at the surface.
"""
back_pos = self.PrismWorld.Omega.positions().array()[
np.isclose(self.PrismWorld.Omega.positions().array()[:, 1],
-self.y_depth)]
back_ids = np.arange(self.PrismWorld.Omega.nodeCount())[
np.isclose(self.PrismWorld.Omega.positions().array()[:, 1],
-self.y_depth)]
back_topo = self.PrismWorld.get_topo_vals(back_pos)
front_pos = self.PrismWorld.Omega.positions().array()[
np.isclose(self.PrismWorld.Omega.positions().array()[:, 1],
self.y_depth)]
front_ids = np.arange(self.PrismWorld.Omega.nodeCount())[
np.isclose(self.PrismWorld.Omega.positions().array()[:, 1],
self.y_depth)]
front_topo = self.PrismWorld.get_topo_vals(front_pos)
for ci in range(len(self.xzmesh.cells())):
fulltri1 = True
fulltri2 = True
# resolve intersection of added points on the front and back
# edges of the prisms in the surface mesh
for ijk in [[0, 1, 2], [1, 2, 0], [2, 0, 1]]:
if (self.xzmesh.cell(ci).node(ijk[0]).id() in self.surf_nodes and
self.xzmesh.cell(ci).node(ijk[1]).id() in self.surf_nodes):
isects0 = []
isects1 = []
for idx, pos in enumerate(back_pos):
if np.isclose(pos[2], back_topo[idx]):
if mu.is_between_1D(
[back_nodes[self.xzmesh.cell(
ci).node(ijk[0]).id()].x()],
[back_nodes[self.xzmesh.cell(
ci).node(ijk[1]).id()].x()],
pos[:2]):
fulltri1 = False
isects0.append(back_ids[idx])
for idx, pos in enumerate(front_pos):
if np.isclose(pos[2], front_topo[idx]):
if mu.is_between_1D(
[front_nodes[self.xzmesh.cell(
ci).node(ijk[0]).id()].x()],
[front_nodes[self.xzmesh.cell(
ci).node(ijk[1]).id()].x()],
pos[:2]):
fulltri2 = False
isects1.append(front_ids[idx])
if len(isects0) > 0:
sorted_ids = self.sort_surface_node_ids(
[back_nodes[self.xzmesh.cell(
ci).node(ijk[0]).id()].id()] + isects0 +
[back_nodes[self.xzmesh.cell(
ci).node(ijk[1]).id()].id()])
mu.npoly_to_triangles(
self.PrismWorld.Omega,
sorted_ids + [back_nodes[self.xzmesh.cell(
ci).node(ijk[2]).id()].id()])
else:
fulltri1 = True
if len(isects1) > 0:
sorted_ids = self.sort_surface_node_ids(
[front_nodes[self.xzmesh.cell(
ci).node(ijk[0]).id()].id()] + isects1 +
[front_nodes[self.xzmesh.cell(
ci).node(ijk[1]).id()].id()])
mu.npoly_to_triangles(
self.PrismWorld.Omega,
sorted_ids + [front_nodes[self.xzmesh.cell(
ci).node(ijk[2]).id()].id()])
else:
fulltri2 = True
if fulltri1:
self.PrismWorld.Omega.createTriangleFace(
self.PrismWorld.Omega.node(
back_nodes[self.xzmesh.cell(ci).node(0).id()].id()),
self.PrismWorld.Omega.node(
back_nodes[self.xzmesh.cell(ci).node(1).id()].id()),
self.PrismWorld.Omega.node(
back_nodes[self.xzmesh.cell(ci).node(2).id()].id()))
if fulltri2:
self.PrismWorld.Omega.createTriangleFace(
self.PrismWorld.Omega.node(
front_nodes[self.xzmesh.cell(ci).node(0).id()].id()),
self.PrismWorld.Omega.node(
front_nodes[self.xzmesh.cell(ci).node(1).id()].id()),
self.PrismWorld.Omega.node(
front_nodes[self.xzmesh.cell(ci).node(2).id()].id()))
[docs]
def add_prisms(self, back_nodes, front_nodes):
"""
Add nodes and triangle faces to include horizontal prisms as conformal
PLC.
"""
already_added = []
z = [0., 0., 0.]
for ci in range(len(self.xzmesh.cells())):
if self.topo is not None:
z = self.PrismWorld.add_topo(np.array(
[[self.xzmesh.cell(ci).node(0).x(), 0., 0.],
[self.xzmesh.cell(ci).node(1).x(), 0., 0.],
[self.xzmesh.cell(ci).node(2).x(), 0., 0.]]))[:, 2]
if (self.xzmesh.cell(ci).node(0).z() < z[0] -
self.PrismWorld.a_tol and
self.xzmesh.cell(ci).node(1).z() < z[1] -
self.PrismWorld.a_tol and
self.xzmesh.cell(ci).node(2).z() < z[2] -
self.PrismWorld.a_tol):
for ij in [[0, 1], [1, 2], [2, 0]]:
poly = [back_nodes[self.xzmesh.cell(
ci).node(ij[0]).id()].id(),
back_nodes[self.xzmesh.cell(
ci).node(ij[1]).id()].id(),
front_nodes[self.xzmesh.cell(
ci).node(ij[1]).id()].id(),
front_nodes[self.xzmesh.cell(
ci).node(ij[0]).id()].id()]
if sorted(poly) not in already_added:
already_added.append(sorted(poly))
mu.npoly_to_triangles(self.PrismWorld.Omega, poly)
for ijk in [[0, 1, 2], [1, 2, 0], [2, 0, 1]]:
if (np.isclose(self.xzmesh.cell(ci).node(ijk[0]).z(),
z[ijk[0]], atol=1) and not
np.isclose(self.xzmesh.cell(ci).node(ijk[1]).z(),
z[ijk[1]], atol=1) and not
np.isclose(self.xzmesh.cell(ci).node(ijk[2]).z(),
z[ijk[2]], atol=1)):
ids = mu.find_nodes_on_segment(
[self.prism_coords[
self.xzmesh.cell(ci).node(ijk[0]).id(), 0],
-self.y_depth, 0.],
[self.prism_coords[
self.xzmesh.cell(ci).node(ijk[0]).id(), 0],
self.y_depth, 0.],
self.PrismWorld.Omega, self.topo)
mu.npoly_to_triangles(
self.PrismWorld.Omega,
[back_nodes[self.xzmesh.cell(ci).node(
ijk[1]).id()].id()] +
ids +
[front_nodes[self.xzmesh.cell(ci).node(
ijk[1]).id()].id()])
mu.npoly_to_triangles(
self.PrismWorld.Omega,
[back_nodes[self.xzmesh.cell(ci).node(
ijk[2]).id()].id()] +
ids +
[front_nodes[self.xzmesh.cell(ci).node(
ijk[2]).id()].id()])
if (np.isclose(self.xzmesh.cell(ci).node(ijk[0]).z(),
z[ijk[0]], atol=1) and
(np.isclose(self.xzmesh.cell(ci).node(ijk[0]).x(),
np.min(self.prism_coords[:, 0]), atol=1) or
np.isclose(self.xzmesh.cell(ci).node(ijk[0]).x(),
np.max(self.prism_coords[:, 0]), atol=1)) and
(np.isclose(self.xzmesh.cell(ci).node(ijk[1]).z(),
z[ijk[1]], atol=1) or
np.isclose(self.xzmesh.cell(ci).node(ijk[2]).z(),
z[ijk[2]], atol=1))):
ids = mu.find_nodes_on_segment(
[self.prism_coords[
self.xzmesh.cell(ci).node(ijk[0]).id(), 0],
-self.y_depth, 0.],
[self.prism_coords[
self.xzmesh.cell(ci).node(ijk[0]).id(), 0],
self.y_depth, 0.],
self.PrismWorld.Omega, self.topo)
if np.isclose(self.xzmesh.cell(ci).node(ijk[1]).z(),
z[ijk[1]], atol=1.):
zid = ijk[2]
elif np.isclose(self.xzmesh.cell(ci).node(ijk[2]).z(),
z[ijk[2]], atol=1.):
zid = ijk[1]
mu.npoly_to_triangles(
self.PrismWorld.Omega,
[back_nodes[self.xzmesh.cell(
ci).node(zid).id()].id()] +
ids +
[front_nodes[self.xzmesh.cell(
ci).node(zid).id()].id()])
self.PrismWorld.add_marker('prism ' + str(ci),
[self.xzmesh.cell(ci).center().x(),
self.xzmesh.cell(ci).center().y(),
self.xzmesh.cell(ci).center().z()])
[docs]
def sort_surface_node_ids(self, ids):
"""
Sort nodes at the surface for conformal triangle insertion.
"""
if len(ids) > 3:
return(np.array(ids)[
np.argsort(self.PrismWorld.Omega.positions(
).array()[ids, 0])].tolist())
else:
return(ids)