Source code for custEM.meshgen.invmesh_tools

# -*- 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)