# -*- coding: utf-8 -*-
"""
@author: Rochlitz.R
"""
from custEM.meshgen import meshgen_utils as mu
from custEM.meshgen.dem_interpolator import DEM
from custEM.meshgen.bathy_tools import Bathy
from custEM.misc.misc import read_paths
from custEM.misc.misc import make_directories
import pygimli as pg
from pygimli import meshtools
import custEM as ce
import numpy as np
import json
import os
import sys
import shutil
import inspect
[docs]
class BlankWorld:
"""
Toolbox for mesh generation using pyGIMLi functionalities to generate a
world **Omega**, add structures such as interfaces, anomalies, transmitters
etc. to this world and finally call TetGen to built a mesh based on the
exported polyfile.
Important notice
----------------
This class was devloped over four years and provides useful functionalities
for the mesh generation, but there are many parts which need an overhaul
in the future. Maybe an updated version of these tools will be moved to
pyGIMLi in future and updated to work with TetGen 1.6
"""
def __init__(self, **raw_kwargs):
"""
Initialize **BlankWorld** class to build a *fullspace*,
*halfspace* or *layered_earth* mesh.
Keyword arguments
-----------------
**General''
- name = 'NoNameMesh', type str
name of the mesh
- x_dim, y_dim, z_dim = [-1e4, 1e4], type list of two float
specify dimensions of world (without boundary mesh)
- backup_script = True, type bool
store a backup of the used mesh-generation script in the *p_dir*,
useful to see the configuration used to build a mesh afterwards
- triangle_q = 34., type float
quality (see *triangle* documentation) of the triangle-surface mesh
**Surface Mesh**
- inner_area = None, type str
can be used to specify an "inner" area on the surface mesh with
other trinagle sizes than the surrounding "outer" area, useful
to refine an area of the surface mesh without explicitly defining
structurs such as lines or refined receiver postions; possible
options are *'box'* (or *'rectangle'*), *'circle'* and *'ellipse'*
- inner_area_size = [], type list of one or two floats
specify size of inner area:
- half of x- and y-dimensions for *rectangle*
- radius for *'circle'*
- radii a and b for *'ellipse'*
- inner_area_shift = [0., 0.], type list of two float
shift center of inner area about [x, y] - offsets
- inner_area_cell_size = None, type float
set or owerwrite the maximum area of cells on the surface fo
the "inner area"
- outer_area_cell_size = None, type float
set or owerwrite the maximum area of cells on the surface fo
the "outer area"
**General cell sizes**
- interface_cell_sizes = None, type list of floats
specify the maximum area of triangle-cells on subsurface interfaces
- layer_cell_sizes = None, type list of floats
specify cell size for *n_layers* in layered-earth meshes
- subsurface_cell_size = 0., type float
use non-zero value to restrict cell_sizes of any subsurface domain;
overwritten by *layer_cell_sizes*
- airspace_cell_size = 0., type float
use non-zero value to restrict cell_sizes of airspace
- anom_cell_size = 0., type float
general contraint for cell sizes of anomalies; usually it makes
more sense to specify a *cell_size* for each anomaly while adding
it
- boundary_mesh_cell_size = 0., type float
use non-zero value to restrict cell_sizes of boundary_mesh
**Topography**
- topo = None, type str or function
specify topography file to be used to adjust heights of nodes in
the surface mesh; alternatively, define a functions z = f(x, y)
- zeroing = False
shift the complete mesh about the mean height value of the
topography
- subsurface_topos = None, type list of topography information
Modify to incorporate varying subsurface layer thicknesses with
using either synthetic topography definitions or real DEM data
- centering = False, type bool
set **True**, if the given DEM should be centered relative to the
computational domain coordinate frame
- easting_shift, northing_shift = None, type float
if centering is **False**, a custom offset of a given DEM relaive
to the computational domain in x- and y- direction can be set
- rotation (RAD) = None, type float ()
rotation of the DEM around the center coordinate
(set via the *centering* flag or *easting/northing_shift*)
relative to the computational domain, given in rad, not degree
**Bathymetry**
- bathymetry = False, type bool
if **True**, automatically include water domains for any height
values smaller than zero, specified by the topography file
- coast_refinement = None, type bool
deprecated, see *bathy_tools*
- min_dist = 50., type float
minimum distance between nodes on coast_line
- min_points = 10, type int
minimum number of identified nodes on coast line to create
a water domain; used to avoid very tiny water domains
- water_cell_size = 0., type float
use non-zero value to restrict water cell size if bathymetry=True
- water_surface_cell_size = 1e6, type float
constrain triangle sizes on water surface
**Miscellanous**
- save_geometry = False, type bool
write seperate *npz* file containing rx, tx and specific
refinement and geometry information for re-usage to build
inversion meshes
- intersections = False, type bool
so far only for internal usage to remove identical polygones of
"touching" faces
- split_subsurface_anomalies = False, type bool
can be used to split multiple "touching" blocks with different
heights automatically to avoid intersections
- a_tol = 1e-2, type float
absolut tolerance for a couple of geometric calculations
- m_dir = 'meshes', type str
main mesh directory
- r_dir = 'results', type str
main results directory
- p_dir = '/para', type str
parameter directory, relative to *m_dir*
- t_dir = '/para/topo', type str
directory containing topography data files, relative to *m_dir*
- s_dir = '/_mesh', type str
storage directory for *.mesh*-field output, relative to *m_dir*
- w_dir = '/mesh_create', type str
working directory for mesh generation, relative to *m_dir*
"""
# the following parameters are initialized here for internal usage
# do NOT overwrite these parameters with other values!
self.world_type = 'box'
self.Omega = pg.Mesh(3)
self.Omega.n_polys = []
self.region_info = []
self.marker_counter = 0
self.n_layers = 1
self.layer_depths = []
self.edge_node_ids = []
self.outcrop_edges = []
self.preserve_edges = True
self.preserve_nodes = False
self.edge_marker = []
self.poly_marker = []
self.export_surface_topo = True
self.export_omega = True
self.rxs = []
self.ground_tri = []
self.air_tri = []
self.txs = []
self.grounding = []
self.id_offsets = []
self.interfaces = []
self.dem_dict = dict()
# name and default world dimensions, everything is defined in SI units
self.name = 'NoNameMesh'
self.x_dim = [-1e4, 1e4]
self.y_dim = [-1e4, 1e4]
self.z_dim = [-1e4, 1e4]
self.backup_script = True
# adjust these parameters to define a fine-discretized area at the
# surface of subsurface interfaces
self.inner_area_cell_size = 0.
self.outer_area_cell_size = 0.
self.inner_area_size = []
self.inner_area = None
self.inner_area_shift = [0., 0.]
self.interface_cell_sizes = None
self.triangle_q = 32.
self.intersections = False
self.split_subsurface_anomalies = None
self.anom_cell_size = 0.
# volumetric cell constraints for airspace and subsurface
self.subsurface_cell_size = 0.
self.airspace_cell_size = 0.
self.boundary_mesh_cell_size = 0.
self.water_cell_size = 0.
self.layer_cell_sizes = None
self.a_tol = 1e-2
self.bathymetry = False
self.coast_refinement = None
self.water_surface_cell_size = 1e6
self.min_dist = 50.
self.min_points = 10
self.z_approx = 0.
# parameters regarding documentation of mesh generation and topography
self.save_geometry = False
self.topo = None
self.zeroing = False
self.exaggeration = None
self.subsurface_topos = None
self.easting_shift = None
self.northing_shift = None
self.centering = False
self.rotation = None
self.t_dir, self.p_dir, self.s_dir, self.w_dir = None, None, None, None
self.r_dir, self.m_dir = read_paths(os.path.dirname(ce.__file__) +
'/misc/paths.dat')
for key in raw_kwargs:
if key not in self.__dict__:
print('Warning! Unknown BlankWorld attribute set:', key)
print("Aborting before something unexpected happens!")
raise SystemExit
self.update_kwargs(raw_kwargs)
if 'p_dir' not in raw_kwargs:
self.p_dir = self.m_dir + '/para'
if 't_dir' not in raw_kwargs:
self.t_dir = '.'
if 's_dir' not in raw_kwargs:
self.s_dir = self.m_dir + '/_mesh'
if 'w_dir' not in raw_kwargs:
self.w_dir = self.m_dir + '/mesh_create'
make_directories(self.r_dir, self.m_dir, 'E_t', m_dir_only=True)
if self.preserve_nodes:
self.preserve_nodes = []
if not os.path.isdir(self.w_dir):
os.makedirs(self.w_dir)
if not os.path.isdir(self.p_dir):
os.makedirs(self.p_dir)
if not os.path.isdir(self.t_dir):
os.makedirs(self.t_dir)
if not os.path.isdir(self.s_dir):
os.makedirs(self.s_dir)
if self.bathymetry:
# if os.path.isfile(self.w_dir + '/coast_line.npy'):
# self.bathy = Bathy(self.t_dir + '/' + self.topo,
# min_dist=self.min_dist, do_nothing=True)
# self.bathy.marine_marker = np.load(self.w_dir + '/marmar.npy')
# self.bathy.coast_line = np.load(self.w_dir + '/coast_line.npy')
# else:
self.bathy = Bathy(self.t_dir + '/' + self.topo,
min_dist=self.min_dist,
min_points=self.min_points,
coast_refinement=self.coast_refinement,
centering=self.centering,
easting_shift=self.easting_shift,
northing_shift=self.northing_shift)
# for line in self.bathy.coast_line:
# print('coast line #n: ', line)
if self.world_type == 'box':
self.create_box()
print('=' * 80)
print('Domain initialization:')
print('--> x limit = ' + str(self.x_dim))
print('--> y limit = ' + str(self.y_dim))
print('--> z limit = ' + str(self.z_dim))
else:
print('Fatal Error! Currently, only "box"-worlds are allowed!')
raise SystemExit
if self.backup_script:
try:
shutil.copy2(inspect.stack()[1][1], self.p_dir + '/' +
self.name + '_script_backup.py')
except FileNotFoundError:
print('Warning! Failed to save backup script. Continuing ...')
pass
[docs]
def create_box(self):
"""
Create bounding box from given limits and initialize **inner_area** and
**outer_area** properties, if a separation into a fine-discretized
inner (observation) area at the surface and outer area (boundary-area)
was specified during intitializing the **BlankWorld** class.
"""
self.box = np.array([[self.x_dim[0], self.y_dim[0], self.z_dim[1]],
[self.x_dim[1], self.y_dim[0], self.z_dim[1]],
[self.x_dim[1], self.y_dim[1], self.z_dim[1]],
[self.x_dim[0], self.y_dim[1], self.z_dim[1]],
[self.x_dim[0], self.y_dim[0], self.z_dim[0]],
[self.x_dim[1], self.y_dim[0], self.z_dim[0]],
[self.x_dim[1], self.y_dim[1], self.z_dim[0]],
[self.x_dim[0], self.y_dim[1], self.z_dim[0]]])
for j in range(len(self.box)):
self.Omega.createNode(self.box[j, 0], self.box[j, 1],
self.box[j, 2])
self.top_id_mapping = [[1, 0], [2, 1], [3, 2], [0, 3]]
self.bottom_id_mapping = [[5, 4], [6, 5], [7, 6], [4, 7]]
self.box_node_count = self.Omega.nodeCount()
if self.inner_area is None:
self.inner_area_frame = []
elif self.inner_area == 'circle':
print(' - using circular shape for inner area refinement.')
if len(self.inner_area_size) != 1:
print('Warning! Given "inner_area" size should be [radius] '
'for a circle!')
self.inner_area_frame = mu.create_circle(
self.inner_area_size[0], 32, 0.)
elif self.inner_area == 'ellipse':
print(' - using ellipsoidal shape for inner area refinement.')
if len(self.inner_area_size) != 2:
print('Warning! Given "inner_area" size should be '
'[heigth / 2, width / 2] for an ellipse!')
if self.inner_area_size[0] == self.inner_area_size[1]:
print('Notice!, "a" and "b" have the same size. Creating '
'"circle" instead.')
self.inner_area_frame = mu.create_circle(
self.inner_area_size[0], 32, 0.)
else:
self.inner_area_frame = mu.create_ellipse(
self.inner_area_size[0],
self.inner_area_size[1], 32, 0.)
elif self.inner_area == 'rectangle' or 'box':
print(' - using rectangular shape for inner area refinement.')
if len(self.inner_area_size) != 2:
print('Warning! Given "inner_area" size should be '
'[heigth / 2, width / 2] for a rectangle!')
self.inner_area_frame = np.array(
[[-self.inner_area_size[0], -self.inner_area_size[1], 0.],
[-self.inner_area_size[0], self.inner_area_size[1], 0.],
[self.inner_area_size[0], self.inner_area_size[1], 0.],
[self.inner_area_size[0], -self.inner_area_size[1], 0.]])
if len(self.inner_area_frame) > 0:
self.inner_area_shift.extend([0.])
self.inner_area_frame += self.inner_area_shift
[docs]
def build_surface(self, insert_line_tx=[], insert_loop_tx=[],
insert_points=[], insert_lines=[],
insert_loops=[], insert_paths=[],
insert_ground_tri=[],
closed_path=False, subsurface_interface=False,
line_marker=None, loop_marker=None, path_marker=None,
line_tx_marker=None, loop_tx_marker=None,
point_marker=None, add_to_existing=False):
"""
Create pyGIMLi surface polygon for 2D meshing with "triangle" later on.
Note! Structures at the surface such as CSEM transmitters or
observation lines / points have to incorporated here.
Keyword arguments
-----------------
- insert_line_tx = [], type list
list of grounded transmitter paths, numbered successively
- insert_loop_tx = [], type list
list of ungrounded transmitter paths; note, loop_tx will be always
added and numbered after all line tx
- insert_points = [], type list
list of meshgen_utils "pointSet" definitions
- insert_lines = [], type list
list of meshgen_utils "line" definitions
- insert_loops = [], type list
list of meshgen_utils "loop" definitions
- insert_paths = [], type list
list of paths descriped by coordinates
- closed_path = False, type bool or list of len == len(insert_paths)
close the segsment between the first and last point of a given
"PATH" to either generate a grounded (line) or ungrounded (loop)
CSEM source or in general, open or closed polygons
- subsurface_interface = False, type bool
for internal usage only, specifies if the 2D polygon will be on the
surface or a subsurface interface
"""
if isinstance(closed_path, bool): # if True or False, used for all
closed_path = [closed_path] * len(insert_paths)
else:
if len(closed_path) != len(insert_paths):
raise Exception(
'"closed_path" must be a list with boolean '
'values with the same length as the insert_paths list or a'
' single boolean value used for all insert_paths.')
if len(insert_line_tx) != 0:
[self.txs.append(line_tx) for line_tx in insert_line_tx]
self.grounding.extend([True] * len(insert_line_tx))
if len(insert_loop_tx) != 0:
[self.txs.append(loop_tx) for loop_tx in insert_loop_tx]
self.grounding.extend([False] * len(insert_loop_tx))
if point_marker is None:
point_marker = [99] * len(insert_points)
if line_marker is None:
line_marker = [99] * len(insert_lines)
if loop_marker is None:
loop_marker = [99] * len(insert_loops)
if path_marker is None:
path_marker = [99] * len(insert_paths)
if line_tx_marker is None:
line_tx_marker = [49] * len(insert_line_tx)
if loop_tx_marker is None:
loop_tx_marker = [49] * len(insert_loop_tx)
if add_to_existing:
surface = self.surface
else:
if not self.bathymetry:
surface = meshtools.polytools.createPolygon(self.box[:4, 0:2],
isClosed=True)
else:
tmp = []
for jj, line in enumerate(self.bathy.coast_line):
for j in range(len(line)):
if line[j][0] in [self.x_dim[0], self.x_dim[1]] or \
line[j][1] in [self.y_dim[0], self.y_dim[1]]:
if not (line[j][0] in [self.x_dim[0],
self.x_dim[1]] and
line[j][1] in [self.y_dim[0],
self.x_dim[1]]):
print(j)
tmp.append([jj, j])
finalbox = self.box[0, 0:2].reshape(1, 2)
searchbox = np.array([[[self.x_dim[0] + 2 * self.a_tol,
self.y_dim[0]],
[self.x_dim[1] - 2 * self.a_tol,
self.y_dim[0]]],
[[self.x_dim[1],
self.y_dim[0] + 2 * self.a_tol],
[self.x_dim[1],
self.y_dim[1] - 2 * self.a_tol]],
[[self.x_dim[0] + 2 * self.a_tol,
self.y_dim[1]],
[self.x_dim[1] - 2 * self.a_tol,
self.y_dim[1]]],
[[self.x_dim[0],
self.y_dim[0] + 2 * self.a_tol],
[self.x_dim[0],
self.y_dim[1] - 2 * self.a_tol]]])
self.merge_ids = []
for k in range(4):
tmp2 = []
if k != 0:
finalbox = np.vstack((finalbox, self.box[k, 0:2]))
for ids in tmp:
if mu.is_between_2D(
searchbox[k, 0, :],
searchbox[k, 1, :],
self.bathy.coast_line[ids[0]][ids[1], :2]):
tmp2.append(
self.bathy.coast_line[ids[0]][ids[1], :2])
points = np.array(tmp2)
if len(points) > 0:
if k == 0:
to_add = points[points[:, 0].argsort()]
if k == 1:
to_add = points[points[:, 1].argsort()]
if k == 2:
to_add = points[points[:, 0].argsort()][::-1]
if k == 3:
to_add = points[points[:, 1].argsort()][::-1]
else:
to_add = []
for point in to_add:
self.merge_ids.extend([len(finalbox)])
finalbox = np.vstack((finalbox, point))
surface = meshtools.polytools.createPolygon(
finalbox, isClosed=True)
# if len(self.merge_ids) != len(tmp):
# print('Fatal error! with merge ids in line 535')
# print(self.merge_ids, tmp)
# raise SystemExit
edges = []
for j in range(len(insert_points)):
n = surface.createNodeWithCheck([insert_points[j][0],
insert_points[j][1],
insert_points[j][2]])
n.setMarker(point_marker[j])
if type(self.preserve_nodes) is list:
self.preserve_nodes.append(n)
for j in range(len(insert_lines)):
nn = [surface.createNodeWithCheck([o[0], o[1], o[2]]) for
o in insert_lines[j]]
[n.setMarker(line_marker[j]) for n in nn]
if type(self.preserve_nodes) is list:
[self.preserve_nodes.append(n) for n in nn]
surface, edges = mu.add_edges(
surface, edges, nn, len(insert_lines[j]), line_marker[j])
self.edge_marker.extend(
[line_marker[j]] * (len(insert_lines[j]) - 1))
for j in range(len(insert_loops)):
nn = [surface.createNodeWithCheck([o[0], o[1], o[2]]) for
o in insert_loops[j]]
[n.setMarker(loop_marker[j]) for n in nn]
if type(self.preserve_nodes) is list:
[self.preserve_nodes.append(n) for n in nn]
surface, edges = mu.add_edges(
surface, edges, nn, len(insert_loops[j]), loop_marker[j], True)
self.edge_marker.extend([loop_marker[j]] * len(insert_loops[j]))
array = np.array(insert_loops[j])
self.add_area_marker(surface, [np.mean(array[:, 0]),
np.mean(array[:, 1]), 0.],
self.inner_area_cell_size, 395)
for j in range(len(insert_line_tx)):
nn = [surface.createNodeWithCheck([o[0], o[1], o[2]]) for
o in insert_line_tx[j]]
[n.setMarker(line_tx_marker[j]) for n in nn]
if type(self.preserve_nodes) is list:
[self.preserve_nodes.append(n) for n in nn]
surface, edges = mu.add_edges(
surface, edges, nn, len(insert_line_tx[j]), line_tx_marker[j])
self.edge_marker.extend(
[line_tx_marker[j]] * (len(insert_line_tx[j]) - 1))
for j in range(len(insert_loop_tx)):
nn = [surface.createNodeWithCheck([o[0], o[1], o[2]]) for
o in insert_loop_tx[j]]
[n.setMarker(loop_tx_marker[j]) for n in nn]
if type(self.preserve_nodes) is list:
[self.preserve_nodes.append(n) for n in nn]
surface, edges = mu.add_edges(
surface, edges, nn, len(insert_loop_tx[j]),
loop_tx_marker[j], True)
self.edge_marker.extend(
[loop_tx_marker[j]] * len(insert_loop_tx[j]))
array = np.array(insert_loop_tx[j])
self.add_area_marker(surface, [np.mean(array[:, 0]),
np.mean(array[:, 1]), 0.],
self.inner_area_cell_size, 395)
for j in range(len(insert_ground_tri)):
if j == 0:
print('... inserting refinement triangles in surface ... ')
self.ground_tri.append(insert_ground_tri[j])
nn = [surface.createNodeWithCheck([o[0], o[1], o[2]]) for
o in insert_rx_tri[j]]
[n.setMarker(99) for n in nn]
if type(self.preserve_nodes) is list:
[self.preserve_nodes.append(n) for n in nn]
surface, edges = mu.add_edges(
surface, edges, nn, len(insert_paths[j]), 99, False)
self.edge_marker.extend(
[path_marker[j]] * (len(insert_paths[j]) - 1))
for j in range(len(insert_paths)):
nn = [surface.createNodeWithCheck([o[0], o[1], o[2]]) for
o in insert_paths[j]]
[n.setMarker(path_marker[j]) for n in nn]
if type(self.preserve_nodes) is list:
[self.preserve_nodes.append(n) for n in nn]
surface, edges = mu.add_edges(
surface, edges, nn, len(insert_paths[j]), path_marker[j],
closed_path[j])
self.edge_marker.extend(
[path_marker[j]] * (len(insert_paths[j]) - 1))
if closed_path[j]:
array = np.array(insert_paths[j])
self.add_area_marker(surface, [np.mean(array[:, 0]),
np.mean(array[:, 1]), 0.],
self.inner_area_cell_size, 395)
self.edge_marker.extend([path_marker[j]])
if self.bathymetry:
for j in range(len(self.bathy.coast_line)):
nn = [surface.createNodeWithCheck([o[0], o[1], o[2]]) for
o in self.bathy.coast_line[j]]
for k in range(len(self.bathy.coast_line[j]) - 1):
surface.createEdge(nn[k], nn[k + 1], marker=48)
surface.createEdge(nn[-1], nn[0], marker=48)
self.add_area_marker(surface, [self.bathy.marine_marker[j][0],
self.bathy.marine_marker[j][1],
0.],
self.water_surface_cell_size, 801)
if not subsurface_interface:
self.surface = surface
else:
self.subsurface_interfaces.append(surface)
[docs]
def add_inv_domains(self, depth, poly=None, rxs=None, cell_size=0.,
x_frame=1e3, y_frame=1e3, z_frame=1e3, d=100.,
marker_positions=None, overwrite_markers=None,
split_depths=None):
if poly is None and rxs is not None:
allrx = rxs[0]
for ri in range(1, len(rxs)):
allrx = np.vstack((allrx, rxs[ri]))
xmin, xmax = np.min(allrx[:, 0]), np.max(allrx[:, 0])
ymin, ymax = np.min(allrx[:, 1]), np.max(allrx[:, 1])
xamin, xamax = np.argmin(allrx[:, 0]), np.argmax(allrx[:, 0])
yamin, yamax = np.argmin(allrx[:, 1]), np.argmax(allrx[:, 1])
poly = np.array([[allrx[yamin, 0], ymin - d, 0.],
[xmax + d, allrx[xamax, 1], 0.],
[allrx[yamax, 0], ymax + d, 0.],
[xmin - d, allrx[xamin, 1], 0.]])
if marker_positions is None:
mpos1 = np.array(poly[np.argmin(poly[:, 0])])
mpos1[0] = np.mean(poly[:, 0])
mpos1[1] = np.mean(poly[:, 1])
mpos1[2] = depth
mpos2 = np.array(poly[np.argmin(poly[:, 0])])
mpos2[0] -= 2 * self.a_tol
mpos2[1] -= self.a_tol
mpos2[2] = depth - 2 * self.a_tol
mpos = [mpos1, mpos2]
else:
mpos = marker_positions
x = [np.min(poly[:, 0]) - x_frame, np.max(poly[:, 0]) + x_frame]
y = [np.min(poly[:, 1]) - y_frame, np.max(poly[:, 1]) + y_frame]
poly2 = np.array([[x[0], x[1], x[1], x[0]],
[y[0], y[0], y[1], y[1]],
[0., 0., 0., 0.]]).T
self.add_surface_anomaly(insert_paths=[poly, poly2],
depths=[depth, depth - z_frame],
cell_sizes=[cell_size, 0.],
marker_positions=mpos,
split_depths=split_depths)
if overwrite_markers is None:
self.overwrite_markers = [2, 3]
else:
self.overwrite_markers = overwrite_markers
[docs]
def add_surface_anomaly(self, insert_paths=[],
depths=None, cell_sizes=None, marker=None,
dips=None, dip_azimuths=None, split_depths=None,
marker_positions=None, split_coords=None):
self.outcrop_marker_positions = marker_positions
self.build_surface(
insert_paths=insert_paths,
closed_path=True, add_to_existing=True)
for jj in range(len(insert_paths)):
self.outcrop_edges.append(insert_paths[jj])
self.surf_anom_depths = depths
self.surf_anom_cell_sizes = cell_sizes
self.surf_anom_dips = dips
self.surf_anom_dip_azimuths = dip_azimuths
self.surf_anom_marker = marker
# self.split_depths = []
# if all(depth == depths[0] for depth in depths):
# for val in depths:
# self.split_depths.append([val])
# else:
# unique_depths = np.sort(np.unique(np.array(depths)))
# if all(depth > 0. for depth in unique_depths):
# unique_depths = unique_depths[::-1]
# for jj in range(len(insert_paths)):
# tmp = []
# for val in unique_depths[::-1]:
# tmp.append(val)
# if depths[jj] == val:
# break
# self.split_depths.append(tmp)
if split_depths is None:
self.split_depths = [[dep] for dep in self.surf_anom_depths]
else:
self.split_depths = split_depths
if split_coords is None:
self.split_coords = [None] * len(insert_paths)
else:
self.split_coords = split_coords
[docs]
def add_intersecting_anomaly(self,
surface_intersection=False,
surface_path=[],
intersection_paths=[],
intersecting_layers=None,
cell_size=None,
bottom=None,
top=None,
marker=None,
marker_position=None):
if surface_intersection:
self.build_surface(insert_paths=[surface_path],
closed_path=True,
add_to_existing=True,
marker=49)
if not self.intersections:
self.intersections = True
self.il_shift = []
self.intersection_polys = []
self.intersection_paths = []
self.intersecting_layers = []
self.anom_top = []
self.anom_bottom = []
self.ianom_cell_sizes = []
self.anom_marker = []
self.anom_marker_position = []
self.il_shift.append(int(np.min(intersecting_layers)))
self.intersection_polys.append([])
self.intersection_paths.append(intersection_paths)
self.intersection_paths[-1] = np.array(self.intersection_paths[-1])
self.intersecting_layers.append(intersecting_layers)
self.anom_top.append(top)
self.anom_bottom.append(bottom)
self.ianom_cell_sizes.append(cell_size)
self.anom_marker.append(marker)
if marker_position is None:
print('Error! Please provide marker position explicitly for now as'
' automatic detection is not implemented yet! Aborting ...')
raise SystemExit
self.anom_marker_position.append(marker_position)
[docs]
def build_fullspace_mesh(self, cell_size=None, boundary_marker=99):
"""
Create a fullspace mesh.
Keyword arguments:
------------------
- cell_size = 0., type float
create a fullspace mesh with cells of maximum volume ([m^3]) of
**max_cell_size**
"""
print('... building fullspace ... ')
box = np.array([[self.x_dim[0], self.y_dim[0], self.z_dim[1]],
[self.x_dim[1], self.y_dim[0], self.z_dim[1]],
[self.x_dim[1], self.y_dim[1], self.z_dim[1]],
[self.x_dim[0], self.y_dim[1], self.z_dim[1]],
[self.x_dim[0], self.y_dim[0], self.z_dim[0]],
[self.x_dim[1], self.y_dim[0], self.z_dim[0]],
[self.x_dim[1], self.y_dim[1], self.z_dim[0]],
[self.x_dim[0], self.y_dim[1], self.z_dim[0]]])
nn = [self.Omega.createNodeWithCheck([o[0], o[1], o[2]]) for o in box]
f1 = self.Omega.createQuadrangleFace(nn[0], nn[1], nn[2], nn[3])
f2 = self.Omega.createQuadrangleFace(nn[4], nn[5], nn[6], nn[7])
f3 = self.Omega.createQuadrangleFace(nn[0], nn[1], nn[5], nn[4])
f4 = self.Omega.createQuadrangleFace(nn[1], nn[2], nn[6], nn[5])
f5 = self.Omega.createQuadrangleFace(nn[2], nn[3], nn[7], nn[6])
f6 = self.Omega.createQuadrangleFace(nn[3], nn[0], nn[4], nn[7])
f1.setMarker(boundary_marker)
f2.setMarker(boundary_marker)
f3.setMarker(boundary_marker)
f4.setMarker(boundary_marker)
f5.setMarker(boundary_marker)
f6.setMarker(boundary_marker)
if cell_size is None:
cell_size = self.subsurface_cell_size
self.add_marker('fullspace ', [self.x_dim[0] + self.a_tol,
self.y_dim[0] + self.a_tol,
self.z_dim[0] + self.a_tol],
cell_size=cell_size)
[docs]
def build_halfspace_mesh(self, boundary_marker=99, **mesh_build_kwargs):
"""
Create halfspace mesh.
Keyword arguments:
------------------
- q = 34., type float
quality (see *triangle* documentation) of the surface mesh
- outer_area_cell_size = None, type float
if not specified when initializing the **BlankWorld** instance,
set or owerwrite the maximum area of cells on the surface for
the "outer area" here
- inner_area_cell_size = None, type float
if not specified when initializing the **BlankWorld** instance,
set or owerwrite the maximum area of cells on the surface for
the "inner area" here
"""
if self.zeroing and self.bathymetry:
print('Fatal error! Do not know what to to if both, *zeroing* and '
'*bathymetry* options are set *True*.')
raise SystemExit
print('... building halfspace ... ')
x0 = (self.x_dim[0] + self.x_dim[1]) / 2
y0 = (self.y_dim[0] + self.y_dim[1]) / 2
self.update_kwargs(mesh_build_kwargs)
self.add_marker('airspace ',
[x0, y0, self.z_dim[1] - 2. * self.a_tol],
cell_size=self.airspace_cell_size)
self.add_interface()
for k in range(4): # for each of the 4 sides faces
nodes = [int(id_k) for id_k in self.edge_node_ids[0][k]]
self.Omega.n_polys.append(nodes + self.bottom_id_mapping[k])
if 'upper_frame' not in self.__dict__:
self.Omega.n_polys.append(nodes + self.top_id_mapping[k])
else:
nodes = [int(id_k) for id_k in self.upper_frame[k]]
self.Omega.n_polys.append(nodes + self.top_id_mapping[k])
self.poly_marker.extend([boundary_marker] * 8)
top = self.add_quadrangle_face(0, 1, 2, 3)
bottom = self.add_quadrangle_face(4, 5, 6, 7)
top.setMarker(boundary_marker)
bottom.setMarker(boundary_marker)
self.add_marker('subsurface ',
[x0, y0, self.z_dim[0]],
cell_size=self.subsurface_cell_size)
self.extend_anomaly_outcrops()
[docs]
def build_layered_earth_mesh(self, n_layers, layer_depths,
insert_struct=None, closed_struct=None,
tx_struct=None, struct_interface=0,
boundary_marker=99, **le_kwargs):
"""
Create layered earth mesh.
Required arguments:
-------------------
- n_layers, type int
number of subsurface layers
- layer_depths, type list of floats
list with len == n_layers - 1, specifying the depths of 1D
subsurface layer interfaces. the z-axis is pointing upwards, so
negative values correspond to increasing depths. This argument
will be overwritten, if subsurface topography is applied
"""
if self.zeroing and self.bathymetry:
print('Fatal error! Do not know what to to if both, *zeroing* and '
'*bathymetry* options are set *True*.')
raise SystemExit
print('... building layered earth ... ')
self.update_kwargs(le_kwargs)
self.n_layers = n_layers
self.layer_depths = [float(elem) for elem in layer_depths]
self.subsurface_interfaces = []
if n_layers - 1 != len(layer_depths) and self.subsurface_topos is None:
print('Fatal error! Number of layers does not match the number of '
'depths of layer boundaries (n_layers - 1)')
raise SystemExit
if self.layer_cell_sizes is None:
self.layer_cell_sizes = [0. for ll in range(n_layers - 1)]
if self.interface_cell_sizes is None:
self.interface_cell_sizes = [self.outer_area_cell_size for
ll in range(n_layers - 1)]
if self.subsurface_topos is None:
self.subsurface_topos = [None for ll in range(n_layers - 1)]
if insert_struct is not None:
if closed_struct is None:
closed_struct = [False for ii in range(len(insert_struct))]
if tx_struct is None:
tx_struct = [False for ii in range(len(insert_struct))]
# airspace:
self.add_interface()
for k in range(4): # for each of the 4 sides faces
if 'upper_frame' not in self.__dict__:
nodes = [int(id_k) for id_k in self.edge_node_ids[0][k]]
self.Omega.n_polys.append(nodes + self.top_id_mapping[k])
else:
nodes = [int(id_k) for id_k in self.upper_frame[k]]
self.Omega.n_polys.append(nodes + self.top_id_mapping[k])
self.poly_marker.extend([boundary_marker] * 4)
top = self.add_quadrangle_face(0, 1, 2, 3) # add top quadrangle face
top.setMarker(boundary_marker)
self.add_marker('airspace ', [self.x_dim[0], self.y_dim[0],
self.z_dim[1] - 2. * self.a_tol],
cell_size=self.airspace_cell_size)
# n_layers - 1 from surface to depth:
for j in range(n_layers - 1):
nc = self.Omega.nodeCount()
frame = meshtools.polytools.createPolygon(self.box[:4, 0:2],
isClosed=True)
if insert_struct is not None and struct_interface == j:
for jj in range(len(insert_struct)):
print(insert_struct[jj])
nn = [frame.createNodeWithCheck([o[0], o[1], o[2]]) for
o in insert_struct[jj]]
for k in range(len(insert_struct[jj]) - 1):
frame.createEdge(nn[k], nn[k + 1], marker=46)
if closed_struct[jj]:
frame.createEdge(nn[-1], nn[0], marker=46)
if tx_struct[jj]:
self.txs.append(insert_struct[jj])
if closed_struct[jj]:
self.grounding.append(False)
else:
self.grounding.append(True)
if self.intersections:
for ii, il in enumerate(self.intersecting_layers):
if j in il:
nn = [frame.createNodeWithCheck([o[0], o[1], o[2]]) for
o in self.intersection_paths[ii][
j - self.il_shift[ii]]]
for k in range(len(
self.intersection_paths[ii][
j - self.il_shift[ii]])-1):
frame.createEdge(nn[k], nn[k + 1], marker=46)
frame.createEdge(nn[-1], nn[0], marker=46)
self.intersection_polys[ii].append(
np.array([n.id() for n in nn], dtype=int))
self.add_area_marker(frame, [self.x_dim[0] + self.a_tol,
self.y_dim[0] + self.a_tol, 0.],
self.interface_cell_sizes[j], marker=555)
if self.subsurface_topos[j] is None:
self.add_interface(frame, self.layer_depths[j], nc)
m_depth = self.layer_depths[j]
else:
self.add_interface(frame, self.subsurface_topos[j], nc)
m_depth = self.add_topo(np.array([self.x_dim[0],
self.y_dim[0],
0.]).reshape(1, 3),
self.subsurface_topos[j])[0, 2]
for k in range(4): # for each of the 4 sides faces
nodes = [int(id_k) for id_k in self.edge_node_ids[j][k]]
nodes2 = [int(id_k) for id_k in self.edge_node_ids[j + 1][k]]
self.Omega.n_polys.append(nodes + nodes2[::-1])
self.poly_marker.extend([boundary_marker])
self.add_marker('layer #' + str(j) + ' ',
[self.x_dim[0], self.y_dim[0], m_depth],
cell_size=self.layer_cell_sizes[j])
# last layer
for k in range(4): # for each of the 4 sides faces
nodes = [int(id_k) for id_k in self.edge_node_ids[-1][k]]
self.Omega.n_polys.append(nodes + self.bottom_id_mapping[k])
self.poly_marker.extend([boundary_marker] * 4)
bottom = self.add_quadrangle_face(4, 5, 6, 7) # add bottom quad. face
bottom.setMarker(boundary_marker)
self.add_marker('subsurface ', [self.x_dim[0], self.y_dim[0],
self.z_dim[0]],
cell_size=self.subsurface_cell_size)
if self.intersections:
self.extend_intersecting_anomaly()
self.extend_anomaly_outcrops()
[docs]
def add_interface(self, frame=None, topo=None, id_offset=None,
eval_zl=False):
"""
Add more or less horizontal 2D interface (layer interfaces) with or
without topography to Omega.
"""
surf = False
if frame is None:
frame = self.surface
topo = self.topo
surf = True
self.preserve_dummy = False
if self.zeroing:
eval_zl = True
if self.inner_area is not None:
if len(self.inner_area_frame) > 0:
for j in range(len(self.inner_area_frame)):
nn = [self.surface.createNodeWithCheck([
o[0], o[1], o[2]]) for o in self.inner_area_frame]
for k in range(len(self.inner_area_frame) - 1):
self.surface.createEdge(nn[k], nn[k + 1])
self.surface.createEdge(nn[-1], nn[0])
self.add_area_marker(self.surface,
[self.inner_area_shift[0] -
self.inner_area_size[0] + 2. * self.a_tol,
self.inner_area_shift[1], 0.],
self.inner_area_cell_size,
marker=555)
if not self.bathymetry:
self.add_area_marker(self.surface, [self.x_dim[0] + self.a_tol,
self.y_dim[0] + self.a_tol, 0.],
self.outer_area_cell_size,
marker=799)
else:
self.add_area_marker(self.surface, [self.x_dim[0] + self.a_tol,
self.y_dim[0] + self.a_tol, 0.],
self.outer_area_cell_size,
marker=801)
if self.bathymetry:
print('Warning! Hacked area marker for bouillante, should not hurt'
' otherwise, but maybe not ...')
self.add_area_marker(self.surface, [self.x_dim[1] - self.a_tol,
self.y_dim[1] - self.a_tol, 0.],
self.outer_area_cell_size,
marker=802)
interface_mesh = pg.meshtools.createMesh(frame,
quality=self.triangle_q)
# preserve edges of transmitters or observation lines (alternative way)
if self.preserve_edges or self.preserve_dummy:
if surf:
self.preserve_dummy = True
self.preserve_edges = []
for j, edge in enumerate(interface_mesh.boundaries()):
if edge.marker() in [47, 48, 49]:
self.preserve_edges.append(
[edge.node(0).id() + self.box_node_count,
edge.node(1).id() + self.box_node_count])
self.edge_marker.append(edge.marker())
# hack for insert_struct and intersections markers 46
elif edge.marker() in [46]:
self.preserve_edges.append(
[edge.node(0).id() + id_offset,
edge.node(1).id() + id_offset])
self.edge_marker.append(edge.marker())
tmp_marker = [node.marker() for node in interface_mesh.nodes()]
# apply nodes of 2D surface mesh to Omega
interface_mesh.setDimension(3)
surface_node_positions = interface_mesh.positions().array()
surface_node_positions = self.add_topo(surface_node_positions, topo,
eval_zl)
# interface_mesh.exportVTK('smesh.vtk')
if self.export_surface_topo:
if surf:
np.savetxt(self.p_dir + '/' + self.name + '_surf.xyz',
surface_node_positions)
else:
if type(topo) is str:
np.savetxt(self.p_dir + '/' + self.name + '_' +
topo[:-4] + '_triangulated.xyz',
surface_node_positions)
if id_offset is None:
id_offset = self.box_node_count
self.outcrop_surface_ids = []
for jj in range(len(self.outcrop_edges)):
self.outcrop_surface_ids.append(mu.find_edge_node_ids(
self.outcrop_edges[jj][:, :2], interface_mesh,
id_offset=id_offset))
self.id_offsets.extend([id_offset])
self.interfaces.append(interface_mesh)
self.edge_node_ids.append(mu.find_edge_node_ids(
self.box[:4, :2], interface_mesh,
id_offset=id_offset))
if self.bathymetry and surf:
all_bb_ids, ext_bb_ids, int_bb_ids = [], [], []
for j, edge in enumerate(interface_mesh.boundaries()):
if edge.marker() == 48:
if edge.node(0).id() not in all_bb_ids:
all_bb_ids.extend([int(edge.node(0).id())])
if edge.node(1).id() not in all_bb_ids:
all_bb_ids.extend([int(edge.node(1).id())])
all_domain_bb_ids = [item for sl in self.edge_node_ids[0]
for item in sl]
dummy = 0
for idx in all_bb_ids:
if idx + self.box_node_count not in all_domain_bb_ids:
int_bb_ids.extend([idx])
elif idx in self.merge_ids:
int_bb_ids.extend([idx])
else:
ext_bb_ids.extend([idx])
if all(surface_node_positions[
np.array(all_domain_bb_ids, dtype=int) -
self.box_node_count, 2] < 0.):
ext_bb_ids = [item for sl in self.edge_node_ids[0]
for item in sl]
ext_bb_ids[:] = [x - self.box_node_count for x in ext_bb_ids]
for j in int_bb_ids:
surface_node_positions[j] = [surface_node_positions[j, 0],
surface_node_positions[j, 1],
0.0]
bathy_cells = []
bathy_nodes = []
dummy = 0
for jj in range(len(interface_mesh.cells())):
if interface_mesh.cell(jj).marker() == 801:
bathy_cells.extend([jj])
dummy += 1
if interface_mesh.cell(jj).node(0).id() not in bathy_nodes:
bathy_nodes.extend(
[interface_mesh.cell(jj).node(0).id()])
if interface_mesh.cell(jj).node(1).id() not in bathy_nodes:
bathy_nodes.extend(
[interface_mesh.cell(jj).node(1).id()])
if interface_mesh.cell(jj).node(2).id() not in bathy_nodes:
bathy_nodes.extend(
[interface_mesh.cell(jj).node(2).id()])
bathy_ids = np.array(bathy_nodes)
nb = []
for val in bathy_ids:
if surface_node_positions[val, 2] > -2. and \
val not in all_bb_ids:
print('overlapping land or shallow water set to -2.: ',
surface_node_positions[val])
surface_node_positions[val, 2] = -2.
for j in range(len(surface_node_positions)):
node = self.Omega.createNodeWithCheck(
[surface_node_positions[j, 0],
surface_node_positions[j, 1],
surface_node_positions[j, 2]])
if type(self.preserve_nodes) is list:
node.setMarker(tmp_marker[j])
for j in range(len(interface_mesh.cells())):
self.Omega.createTriangleFace(
self.Omega.node(interface_mesh.cell(j).node(
0).id() + id_offset),
self.Omega.node(interface_mesh.cell(j).node(
1).id() + id_offset),
self.Omega.node(interface_mesh.cell(j).node(
2).id() + id_offset),
marker=interface_mesh.cell(j).marker())
if self.bathymetry and surf:
for val in bathy_ids:
node = self.Omega.createNodeWithCheck(
[surface_node_positions[val, 0],
surface_node_positions[val, 1],
0.])
nb.extend([node.id()])
nb = np.array(nb)
if len(ext_bb_ids) > 0:
self.air_frame = [[], [], [], []]
self.upper_frame = [[], [], [], []]
for j in range(4):
mirror_ids = []
for nn in self.edge_node_ids[0][j]:
if nn - self.box_node_count in ext_bb_ids:
tmp = nb[int(np.where(
nn - self.box_node_count == bathy_ids)[0])]
self.air_frame[j].extend([int(tmp)])
mirror_ids.extend([int(nn)])
self.upper_frame[j].extend([int(tmp)])
else:
self.upper_frame[j].extend([int(nn)])
self.air_frame[j].extend([None])
mirror_ids.extend([None])
if None in self.air_frame[j] and not \
all([item is None for item in self.air_frame[j]]):
top_row, bottom_row = [], []
for jj, bbid in enumerate(self.air_frame[j]):
if jj == len(self.air_frame[j]) - 1:
if bbid is None:
if self.air_frame[j][jj - 1] is not None:
top_row.extend([
self.upper_frame[j][jj]])
else:
top_row.extend([bbid])
bottom_row.extend([mirror_ids[jj]])
if len(top_row + bottom_row[::-1]) > 0:
self.Omega.n_polys.append(
top_row + bottom_row[::-1])
self.poly_marker.extend([23])
break
elif (jj != len(self.air_frame[j]) - 1 and
bbid is None):
if self.air_frame[j][jj + 1] is None:
if self.air_frame[j][jj - 1] is not None:
top_row.extend([
self.upper_frame[j][jj]])
self.Omega.n_polys.append(
top_row + bottom_row[::-1])
self.poly_marker.extend([23])
top_row, bottom_row = [], []
else:
top_row.extend([self.upper_frame[j][jj]])
else:
top_row.extend([bbid])
bottom_row.extend([mirror_ids[jj]])
elif (None not in self.air_frame[j] and not
all([item is None for item in self.air_frame[j]])):
if len(self.air_frame[j] + mirror_ids[::-1]) != 0:
self.Omega.n_polys.append(self.air_frame[j] +
mirror_ids[::-1])
self.poly_marker.extend([23])
for jj in bathy_cells:
a = int(nb[np.where(bathy_ids == interface_mesh.cell(
jj).node(0).id())][0])
b = int(nb[np.where(bathy_ids == interface_mesh.cell(
jj).node(1).id())][0])
c = int(nb[np.where(bathy_ids == interface_mesh.cell(
jj).node(2).id())][0])
self.Omega.createTriangleFace(
self.Omega.node(a), self.Omega.node(b),
self.Omega.node(c), marker=222)
# deprecated alternative approach, might be interesting again in future
# new_bathy_frame = mu.find_on_frame(
# self.bathy.coast_line[0], self.Omega.positions().array()[nb])
# dummy = [self.Omega.node(nb[kk]).id() for kk in new_bathy_frame]
# arr = self.Omega.positions().array()[dummy]
# # self.Omega.n_polys.append(dummy)
# # self.poly_marker.extend([23])
# w_surf = meshtools.polytools.createPolygon(
# arr,
# isClosed=True)
# dummy_mesh = pg.meshtools.createMesh(
# w_surf, switches='pzeAY' + 'a' +
# '{:.20f}'.format(self.water_surface_cell_size) +
# 'q' + str(self.triangle_q))
# dummy_mesh.exportVTK('dmesh.vtk')
# dummy_node_positions = dummy_mesh.positions().array()
# dummy_nodes = []
# for jj in range(len(dummy_node_positions)):
# node = self.Omega.createNodeWithCheck(
# [dummy_node_positions[jj, 0],
# dummy_node_positions[jj, 1],
# 0.])
# dummy_nodes.append(node.id())
#
# for kk in range(len(dummy_mesh.cells())):
# self.Omega.createTriangleFace(
# self.Omega.node(
# int(dummy_nodes[dummy_mesh.cell(kk).node(
# 1).id()])),
# self.Omega.node(
# int(dummy_nodes[dummy_mesh.cell(kk).node(
# 2).id()])),
# self.Omega.node(
# int(dummy_nodes[dummy_mesh.cell(kk).node(
# 0).id()])),
# marker=222)
for jj in range(len(self.bathy.marine_marker)):
self.add_marker('marine ', [self.bathy.marine_marker[jj][0],
self.bathy.marine_marker[jj][1],
-2. * self.a_tol],
cell_size=self.water_cell_size,
marker=1)
self.marker_counter += 1
[docs]
def add_quadrangle_face(self, n1, n2, n3, n4):
"""
Add quadrangle facet to Omega defined by nodes n1 - n4.
Required arguments:
-------------------
- n1, n2, n3 n4, type int
four node id's to build a quandrangle.
"""
face = self.Omega.createQuadrangleFace(
self.Omega.node(int(n1)), self.Omega.node(int(n2)),
self.Omega.node(int(n3)), self.Omega.node(int(n4)))
return(face)
[docs]
def add_topo(self, node_positions, topo=None, eval_zl=False):
"""
Apply topography, either DEM information or synthetic descriptions
as z = f(x, y) to the z-values of an interface mesh.
Required arguments:
-------------------
- node_positions, type array of shape (n, 3)
positions of all nodes of the interface mesh
Keyword arguments:
------------------
topo = None, type str or function
if None, **self.topo** will be used for the earth's surface, for
further interfaces, either additional DEM information can be used
or a synthetic definition via a function f(x, y)
"""
if topo is None:
topo = self.topo
if topo is None:
return(node_positions)
if type(topo) is float:
node_positions[:, 2] = topo
elif type(topo) is str:
if topo not in self.dem_dict.keys():
self.dem_dict[topo] = self.init_dem(node_positions, topo)
node_positions[:, 2] = self.dem_dict[topo](node_positions[:, 0],
node_positions[:, 1],
rotation=self.rotation)
else:
node_positions[:, 2] = topo(node_positions[:, 0],
node_positions[:, 1])
if self.zeroing:
if eval_zl:
self.zero_level = np.mean(node_positions[:, 2])
print(' - estimated shift for zero leveling: --> ' +
str(self.zero_level))
node_positions[:, 2] = node_positions[:, 2] - self.zero_level
if self.exaggeration is not None:
node_positions[:, 2] *= self.exaggeration
return(node_positions)
[docs]
def get_topo_vals(self, node_positions, topo=None, z=0.):
"""
Apply topography, either DEM information or synthetic descriptions
as z = f(x, y) to the z-values of an interface mesh.
Required arguments:
-------------------
- node_positions, type array of shape (n, 3)
positions of all nodes of the interface mesh
Keyword arguments:
------------------
topo = None, type str or function
if None, **self.topo** will be used for the earth's surface, for
further interfaces, either additional DEM information can be used
or a synthetic definition via a function f(x, y)
"""
if topo is None:
topo = self.topo
if topo is None:
return(node_positions[:, 2] + z)
if type(topo) is float:
coords = np.ones(len(node_positions)) * topo
elif type(topo) is str:
if topo not in self.dem_dict.keys():
self.dem_dict[topo] = self.init_dem(node_positions, topo)
coords = self.dem_dict[topo](node_positions[:, 0],
node_positions[:, 1],
rotation=self.rotation)
else:
coords = topo(node_positions[:, 0], node_positions[:, 1])
return(coords + z)
[docs]
def extend_world(self, x_fact=10., y_fact=10., z_fact=10., marker=[0, 1],
cell_sizes=None, boundary_marker=99, initial=True):
"""
alias for *there_is_always_a_bigger_world()* method
"""
self.there_is_always_a_bigger_world(x_fact, y_fact, z_fact, marker,
cell_sizes, boundary_marker,
initial=initial)
[docs]
def there_is_always_a_bigger_world(self, x_fact=10., y_fact=10.,
z_fact=10., marker=[0, 1],
cell_sizes=None,
boundary_marker=99, initial=True):
"""
Appends a bigger world (halfspace) to the existing world *Omega*.
This bigger world is synonymous to a "tetrahedron boundary".
Keyword arguments:
------------------
- x_fact, y_fact, z_fact = 10., type float
factor in x-, y- and z-direction multiplied to the original domain
size specified by the **x/y/zlim* class attributes.
- marker = [0, 1], list of two int
specify your own marker values for the bounding halfspace mesh.
"""
print('... appending bigger world ... ')
if initial:
self.old_nodes = []
x = np.mean(self.x_dim)
y = np.mean(self.y_dim)
z = np.mean(self.z_dim)
self.x_dim_big = [x - np.diff(self.x_dim)[0] * x_fact / 2.,
x + np.diff(self.x_dim)[0] * x_fact / 2.]
self.y_dim_big = [y - np.diff(self.y_dim)[0] * y_fact / 2.,
y + np.diff(self.y_dim)[0] * y_fact / 2.]
self.z_dim_big = [z - np.diff(self.z_dim)[0] * z_fact / 2.,
z + np.diff(self.z_dim)[0] * z_fact / 2.]
else:
x = np.mean(self.x_dim_big)
y = np.mean(self.y_dim_big)
z = np.mean(self.z_dim_big)
self.x_dim_big = [x - np.diff(self.x_dim_big)[0] * x_fact / 2.,
x + np.diff(self.x_dim_big)[0] * x_fact / 2.]
self.y_dim_big = [y - np.diff(self.y_dim_big)[0] * y_fact / 2.,
y + np.diff(self.y_dim_big)[0] * y_fact / 2.]
self.z_dim_big = [z - np.diff(self.z_dim_big)[0] * z_fact / 2.,
z + np.diff(self.z_dim_big)[0] * z_fact / 2.]
nc = self.Omega.nodeCount()
id_shift = [[1, 0], [2, 1], [3, 2], [0, 3]]
ids = [idx for e_ids in self.edge_node_ids[0] for idx in e_ids[:-1]]
b_heights = self.Omega.positions().array()[ids]
if self.bathymetry and any([None in ll for ll in self.air_frame]):
if marker == [0, 1]:
marker = [0, 2]
print('... setting original boundary heights to 0 to be able '
'\n to append a conform bigger world PLC ...')
b_heights[:, 2] = 0.
ids2 = [idx for uids in self.upper_frame for idx in uids[:-1]]
for jj, elem in enumerate(ids):
if elem == ids2[jj]:
self.Omega.node(int(elem)).setPos(b_heights[jj])
else:
if not all(z == b_heights[0, 2] for z in b_heights[:, 2]):
print('... averaging original surface boundary heights to be '
'\n able to append a conform bigger world PLC ...')
b_heights[:, 2] = np.mean(b_heights[:, 2])
for jj, elem in enumerate(ids):
self.Omega.node(int(elem)).setPos(b_heights[jj])
self.old_nodes.append(self.Omega.createNode(
[self.x_dim_big[0], self.y_dim_big[0], b_heights[0, 2]]))
self.old_nodes.append(self.Omega.createNode(
[self.x_dim_big[1], self.y_dim_big[0], b_heights[0, 2]]))
self.old_nodes.append(self.Omega.createNode(
[self.x_dim_big[1], self.y_dim_big[1], b_heights[0, 2]]))
self.old_nodes.append(self.Omega.createNode(
[self.x_dim_big[0], self.y_dim_big[1], b_heights[0, 2]]))
self.Omega.createNode([self.x_dim_big[0], self.y_dim_big[0],
self.z_dim_big[0]])
self.Omega.createNode([self.x_dim_big[1], self.y_dim_big[0],
self.z_dim_big[0]])
self.Omega.createNode([self.x_dim_big[1], self.y_dim_big[1],
self.z_dim_big[0]])
self.Omega.createNode([self.x_dim_big[0], self.y_dim_big[1],
self.z_dim_big[0]])
self.Omega.createNode([self.x_dim_big[0], self.y_dim_big[0],
self.z_dim_big[1]])
self.Omega.createNode([self.x_dim_big[1], self.y_dim_big[0],
self.z_dim_big[1]])
self.Omega.createNode([self.x_dim_big[1], self.y_dim_big[1],
self.z_dim_big[1]])
self.Omega.createNode([self.x_dim_big[0], self.y_dim_big[1],
self.z_dim_big[1]])
self.old_nodes.append(self.Omega.createNode(
[np.mean(self.x_dim_big), self.y_dim_big[0], b_heights[0, 2]]))
self.old_nodes.append(self.Omega.createNode(
[self.x_dim_big[1], np.mean(self.y_dim_big), b_heights[0, 2]]))
self.old_nodes.append(self.Omega.createNode(
[np.mean(self.x_dim_big), self.y_dim_big[1], b_heights[0, 2]]))
self.old_nodes.append(self.Omega.createNode(
[self.x_dim_big[0], np.mean(self.y_dim_big), b_heights[0, 2]]))
if initial:
for k in range(4): # for each of the 4 trapezoids
if self.bathymetry and (None in self.air_frame[k] or
self.air_frame[k] == self.upper_frame[k]):
nodes = [int(id_k) for id_k in self.upper_frame[k]]
else:
nodes = [int(id_k) for id_k in self.edge_node_ids[0][k]]
ll = int(len(nodes) / 2.)
self.Omega.n_polys.append(
nodes[:ll+1] + [nc + 12 + int(k)] +
[nc + id_shift[k][1]])
self.Omega.n_polys.append( # add this triangle
nodes[ll:] + [nc + id_shift[k][0]] +
[nc + 12 + int(k)])
self.poly_marker.extend([23, 23])
else:
self.Omega.n_polys.append([nc + 1, nc + 12, nc,
self.old_nodes[-8 -8].id(),
self.old_nodes[-4 -8].id(),
self.old_nodes[-7 -8].id()])
self.Omega.n_polys.append([nc + 2, nc + 13, nc + 1,
self.old_nodes[-7 -8].id(),
self.old_nodes[-3 -8].id(),
self.old_nodes[-6 -8].id()])
self.Omega.n_polys.append([nc + 3, nc + 14, nc + 2,
self.old_nodes[-6 -8].id(),
self.old_nodes[-2 -8].id(),
self.old_nodes[-5 -8].id()])
self.Omega.n_polys.append([nc, nc + 15, nc + 3,
self.old_nodes[-5 -8].id(),
self.old_nodes[-1 -8].id(),
self.old_nodes[-8 -8].id()])
self.poly_marker.extend([23, 23, 23, 23])
for k in range(3):
self.Omega.n_polys.append([nc + int(k)] + [nc + 12 + int(k)] +
[nc + 1 + int(k)] + [nc + 5 + int(k)] +
[nc + 4 + int(k)])
self.Omega.n_polys.append([nc + 3] + [nc + 15] +
[nc] + [nc + 4] + [nc + 7])
self.poly_marker.extend([boundary_marker] * 4)
for k in range(3):
self.Omega.n_polys.append([nc + int(k)] + [nc + 12 + int(k)] +
[nc + 1 + int(k)] + [nc + 9 + int(k)] +
[nc + 8 + int(k)])
self.Omega.n_polys.append([nc + 3] + [nc + 15] +
[nc] + [nc + 8] + [nc + 11])
self.poly_marker.extend([boundary_marker] * 4)
top = self.add_quadrangle_face(nc + 4, nc + 5, nc + 6, nc + 7)
bottom = self.add_quadrangle_face(nc + 8, nc + 9, nc + 10, nc + 11)
top.setMarker(boundary_marker)
bottom.setMarker(boundary_marker)
if cell_sizes is None:
cell_sizes = [self.boundary_mesh_cell_size,
self.boundary_mesh_cell_size]
self.add_marker('subsurface bigger world',
[self.x_dim_big[0], self.y_dim_big[0],
self.z_dim_big[0]],
marker=marker[1],
cell_size=cell_sizes[1])
self.add_marker('airspace bigger world',
[self.x_dim_big[0], self.y_dim_big[0],
self.z_dim_big[1] - 2. * self.a_tol],
marker=marker[0],
cell_size=cell_sizes[0])
[docs]
def extend_anomaly_outcrops(self):
"""
Function for extending anomaly-bodies into a layer. Completely
automated for internal usage. For more details, it is referred to the
**add_surface_anomaly()** method, where users can define the anomalies!
"""
# Define default values if none or provided
if self.outcrop_edges == []:
return
n_sa = len(self.outcrop_edges)
if self.surf_anom_depths is None:
self.surf_anom_depths = [self.z_dim[1] / 20. for j in range(n_sa)]
if self.surf_anom_marker is None:
self.surf_anom_marker = [None for j in range(n_sa)]
if self.surf_anom_dips is None:
self.surf_anom_dips = [0. for j in range(n_sa)]
if self.surf_anom_cell_sizes is None:
self.surf_anom_cell_sizes = [0. for j in range(n_sa)]
if self.surf_anom_dip_azimuths is None:
self.surf_anom_dip_azimuths = [0. for j in range(n_sa)]
if self.split_subsurface_anomalies is None:
self.split_subsurface_anomalies = [False for j in range(n_sa)]
if len(self.surf_anom_depths) != len(self.surf_anom_dips) != \
len(self.surf_anom_dip_azimuths) != len(self.outcrop_edges) != \
len(self.surf_anom_marker) != len(self.surf_anom_cell_sizes):
print('Fatal Error! Provided lists of *depth*, *dip*, and '
'*dip_azimuth* values and given outcrops must be the same!')
raise SystemExit
# iterate over all surface anomalies
lens = []
anom_polys = []
for jj in range(n_sa):
print('... working on extending anomly: ' + str(jj))
n_segs = len(self.outcrop_surface_ids[jj])
surface_frame = []
# iterate over all segments of a poly at the surface and extend
# to 3D facets using the "n_poly" method
old_nn = []
for k in range(len(self.split_depths[jj])):
tmp = []
for ll in range(n_segs):
poly = self.Omega.positions().array()[
self.outcrop_surface_ids[jj][ll]]
# if k == 0 and ll == 0:
# z_approx = np.mean(poly[:, 2])
poly = np.vstack((poly[0, :], poly[-1, :]))
if k == 0:
surface_frame.extend(self.Omega.positions().array()[
self.outcrop_surface_ids[jj][ll]])
poly[:, 2] = self.z_approx + self.split_depths[jj][k]
nn = [self.Omega.createNodeWithCheck(
[o[0], o[1], o[2]]) for o in poly]
tmp.append(nn)
if k == 0:
upper = self.outcrop_surface_ids[jj][ll].tolist()
else:
upper = [old_nn[ll][0].id(), old_nn[ll][1].id()]
lower = [nn[1].id(), nn[0].id()]
anom_polys.append(upper + lower)
lens.append(len(upper + lower))
if k == 0:
frm = np.array(surface_frame)
old_nn = np.array(tmp)
if self.split_subsurface_anomalies[jj]:
self.Omega.n_polys.append([ele[0].id() for ele in tmp])
self.poly_marker.append(23)
if len(self.split_depths) != 1 and k > 0:
if not self.split_subsurface_anomalies[jj]:
self.marker_counter -= 1
if self.split_subsurface_anomalies[jj]:
self.add_marker('surface brick',
[np.mean(frm[:, 0]), np.mean(frm[:, 1]),
self.z_approx + self.split_depths[jj][k] +
self.a_tol],
cell_size=self.surf_anom_cell_sizes[jj])
# add closing poly
if not self.split_subsurface_anomalies[jj]:
self.Omega.n_polys.append([ele[0].id() for ele in tmp])
self.poly_marker.append(23)
if self.surf_anom_dips[jj] != 0.:
for kk, elem in enumerate([ele[0].id() for ele in tmp]):
if kk == 0:
bottom = self.z_approx - \
np.abs(self.surf_anom_depths[jj]) * \
np.cos(np.deg2rad(90 - self.surf_anom_dips[jj]))
hdist = np.abs(self.surf_anom_depths[jj]) * \
np.sin(np.deg2rad(90 - self.surf_anom_dips[jj]))
xx = hdist * \
np.cos(np.deg2rad(self.surf_anom_dip_azimuths[jj]))
yy = hdist * \
np.sin(np.deg2rad(self.surf_anom_dip_azimuths[jj]))
pos = self.Omega.positions().array()[int(elem)]
self.Omega.node(int(elem)).setPos(
[pos[0] + xx, pos[1] + yy, bottom])
# add marker
if not self.split_subsurface_anomalies[jj]:
if self.outcrop_marker_positions is None:
m_x_pos, m_y_pos = np.mean(frm[:, 0]), np.mean(frm[:, 1])
print('... add marker 10 m below interpoalted DEM ...')
m_depth = self.add_topo(
np.array([m_x_pos, m_y_pos, 0.]).reshape(
1, 3))[0, 2] - 10.
else:
m_x_pos, m_y_pos, m_depth =\
self.outcrop_marker_positions[jj]
self.add_marker('surface brick', [m_x_pos, m_y_pos,
m_depth],
cell_size=self.surf_anom_cell_sizes[jj],
marker=self.surf_anom_marker[jj])
print('... incorporation of anomaly ' + str(jj) + ' done ...')
ids = mu.get_unique_poly_ids(anom_polys, np.max(lens))
for idx in ids:
self.Omega.n_polys.append(anom_polys[idx])
self.poly_marker.append(23)
[docs]
def extend_intersecting_anomaly(self):
"""
Function for extending anomaly-bodies which interesect layers.
"""
for ii in range(len(self.intersection_paths)):
self.intersection_ids = []
for k in range(len(self.intersection_polys[ii])):
self.intersection_polys[ii][k] += self.id_offsets[
k + self.il_shift[ii] + 1]
self.intersection_ids.append(mu.find_edge_node_ids(
self.Omega.positions().array()[
self.intersection_polys[ii][k]][:, :2],
self.interfaces[k + self.il_shift[ii] + 1],
id_offset=self.id_offsets[
k + self.il_shift[ii] + 1]))
# close all faces between layers
if len(self.intersection_polys[ii]) > 1:
for k in range(len(self.intersection_polys[ii]) - 1):
for ll in range(len(self.intersection_polys[ii][k])):
self.Omega.n_polys.append(
self.intersection_ids[k][ll].tolist() +
[self.intersection_ids[k+1][ll][0]])
self.Omega.n_polys.append(
self.intersection_ids[k+1][ll].tolist() +
[self.intersection_ids[k][ll][-1]])
self.poly_marker.extend([23])
# close all faces between top and uppermost intersecting layer
if self.anom_top[ii] is not None:
nn = [self.Omega.createNodeWithCheck([o[0], o[1],
self.anom_top[ii]])
for o in self.intersection_paths[ii][0]]
top_ids = [n.id() for n in nn] + [nn[0].id()]
for ll in range(len(self.intersection_polys[ii][0])):
self.Omega.n_polys.append(
self.intersection_ids[0][ll].tolist() +
[top_ids[ll+1], top_ids[ll]])
self.poly_marker.extend([23])
self.Omega.n_polys.append([n.id() for n in nn])
self.poly_marker.extend([23])
# close all faces between bottom and lowermost intersecting layer
if self.anom_bottom[ii] is not None:
nn = [self.Omega.createNodeWithCheck([o[0], o[1],
self.anom_bottom[ii]])
for o in self.intersection_paths[ii][-1]]
bottom_ids = [n.id() for n in nn] + [nn[0].id()]
for ll in range(len(self.intersection_polys[ii][0])):
self.Omega.n_polys.append(
self.intersection_ids[-1][ll].tolist() +
[bottom_ids[ll+1], bottom_ids[ll]])
self.poly_marker.extend([23, 23])
self.Omega.n_polys.append([n.id() for n in nn])
self.poly_marker.extend([23])
self.add_marker('surface brick', self.anom_marker_position[ii],
cell_size=self.ianom_cell_sizes[ii],
marker=self.anom_marker[ii])
[docs]
def add_brick(self, start, stop, cell_size=0., marker=None):
"""
Add rectangular *brick* anomaly to **Omega**. Note that this is a
special case of the **add_plate()** method, which can be used
alternatively.
Required arguments:
-------------------
- start, type list
list of back lower left corner [x_min,, y_min, z_min].
- stop, type list
list of front upper right corner [x_max,, y_max, z_max].
Keyword arguments:
------------------
- cell_size = 0., type float
maximum cell-volume contstraint, default (0.) means no constraint.
"""
print('... adding brick ...')
x_lim = [start[0], stop[0]]
y_lim = [start[1], stop[1]]
z_lim = [start[2], stop[2]]
brick = np.array([[x_lim[0], y_lim[0], z_lim[1]],
[x_lim[1], y_lim[0], z_lim[1]],
[x_lim[1], y_lim[1], z_lim[1]],
[x_lim[0], y_lim[1], z_lim[1]],
[x_lim[0], y_lim[0], z_lim[0]],
[x_lim[1], y_lim[0], z_lim[0]],
[x_lim[1], y_lim[1], z_lim[0]],
[x_lim[0], y_lim[1], z_lim[0]]])
nn = [self.Omega.createNodeWithCheck([o[0], o[1], o[2]])
for o in brick]
# self.Omega.createQuadrangleFace(nn[0], nn[1], nn[2], nn[3])
# self.Omega.createQuadrangleFace(nn[4], nn[5], nn[6], nn[7])
# self.Omega.createQuadrangleFace(nn[0], nn[1], nn[5], nn[4])
# self.Omega.createQuadrangleFace(nn[1], nn[2], nn[6], nn[5])
# self.Omega.createQuadrangleFace(nn[2], nn[3], nn[7], nn[6])
# self.Omega.createQuadrangleFace(nn[3], nn[0], nn[4], nn[7])
# Alternative implementation via n_polys
self.Omega.n_polys.append([nn[0].id(), nn[1].id(),
nn[2].id(), nn[3].id()])
self.Omega.n_polys.append([nn[4].id(), nn[5].id(),
nn[6].id(), nn[7].id()])
self.Omega.n_polys.append([nn[0].id(), nn[1].id(),
nn[5].id(), nn[4].id()])
self.Omega.n_polys.append([nn[1].id(), nn[2].id(),
nn[6].id(), nn[5].id()])
self.Omega.n_polys.append([nn[2].id(), nn[3].id(),
nn[7].id(), nn[6].id()])
self.Omega.n_polys.append([nn[3].id(), nn[0].id(),
nn[4].id(), nn[7].id()])
self.poly_marker.extend([23] * 6)
self.add_marker('brick', [np.mean(x_lim), np.mean(y_lim),
np.mean(z_lim)],
cell_size=cell_size, marker=marker)
[docs]
def add_prism(self, poly, top, bottom, cell_size=0., marker=None,
marker_pos=None):
"""
Add rectangular *brick* anomaly to **Omega**. Note that this is a
special case of the **add_plate()** method, which can be used
alternatively.
Required arguments:
-------------------
- start, type list
list of back lower left corner [x_min,, y_min, z_min].
- stop, type list
list of front upper right corner [x_max,, y_max, z_max].
Keyword arguments:
------------------
- cell_size = 0., type float
maximum cell-volume contstraint, default (0.) means no constraint.
"""
print('... adding prism ...')
nn1 = [self.Omega.createNode([o[0], o[1], top]) for o in poly]
nn2 = [self.Omega.createNode([o[0], o[1], bottom]) for o in poly]
for ll in range(len(nn1) - 1):
self.Omega.createQuadrangleFace(nn1[ll], nn1[ll+1],
nn2[ll+1], nn2[ll])
self.Omega.createQuadrangleFace(nn1[-1], nn1[0], nn2[0], nn2[-1])
self.Omega.n_polys.append([nn.id() for nn in nn1])
self.Omega.n_polys.append([nn.id() for nn in nn2])
self.poly_marker.extend([23, 23])
if marker_pos is None:
print('Oo, need manually defined marker position for now')
raise SystemExit
self.add_marker('brick', [marker_pos[0] + self.a_tol,
marker_pos[1] + self.a_tol,
marker_pos[2] + self.a_tol],
cell_size=cell_size, marker=marker)
[docs]
def add_hollow_cylinder(self, r1, r2, h=100., origin=[0., 0., 0.],
rotation=None, cell_size=0., n_segs=16,
preserve=[], close=False, mark=False,
close1st=False):
"""
Add hollow cylinder-anomaly to **Omega**. In progress!
Required arguments:
-------------------
- r1, type float
inner radius of cylinder.
- r2, type float
outer radius of cylinder.
Keyword arguments:
------------------
- cell_size = 0., type float
maximum cell-volume contstraint, default (0.) means no constraint.
"""
print('... adding cylinder ...')
loops = []
radii = [r1, r2, r1, r2]
h = [origin[2] + h/2., origin[2] + h/2.,
origin[2] - h/2., origin[2] - h/2.]
for j, r in enumerate(radii):
loops.append(mu.loop_c(origin=origin, r=r, n_segs=n_segs, z=h[j]))
loops_rot = []
if rotation is not None:
for circ in loops:
loops_rot.append(mu.rotate_around_point(np.array(circ),
np.deg2rad(rotation), 'x'))
loops = loops_rot
nns = []
for j in range(4):
nns.append([self.Omega.createNodeWithCheck([o[0], o[1], o[2]]) for
o in loops[j]])
if j in preserve:
if type(self.preserve_edges) is not list:
self.preserve_edges = []
for k in range(len(nns[j]) - 1):
self.Omega.createEdge(nns[j][k], nns[j][k + 1])
self.preserve_edges.append(
[nns[j][k].id(), nns[j][k+1].id()])
self.Omega.createEdge(nns[j][-1], nns[j][0])
self.preserve_edges.append(
[nns[j][-1].id(), nns[j][0].id()])
if close and j in [0, 2]:
self.Omega.n_polys.append([int(n.id()) for n in nns[j]])
self.poly_marker.extend([23])
if close1st and j in [0]:
self.Omega.n_polys.append([int(n.id()) for n in nns[j]])
self.poly_marker.extend([23])
for j in range(n_segs - 1):
self.Omega.createQuadrangleFace(nns[0][j], nns[1][j],
nns[1][j+1], nns[0][j+1])
self.Omega.createQuadrangleFace(nns[2][j], nns[3][j],
nns[3][j+1], nns[2][j+1])
self.Omega.createQuadrangleFace(nns[0][j], nns[2][j],
nns[2][j+1], nns[0][j+1])
self.Omega.createQuadrangleFace(nns[1][j], nns[3][j],
nns[3][j+1], nns[1][j+1])
self.Omega.createQuadrangleFace(nns[0][-1], nns[1][-1],
nns[1][0], nns[0][0])
self.Omega.createQuadrangleFace(nns[2][-1], nns[3][-1],
nns[3][0], nns[2][0])
self.Omega.createQuadrangleFace(nns[0][-1], nns[2][-1],
nns[2][0], nns[0][0])
self.Omega.createQuadrangleFace(nns[1][-1], nns[3][-1],
nns[3][0], nns[1][0])
if mark:
self.add_marker('hollow_cylinder',
[origin[0] + r1 + 0.001, origin[1],
origin[2]], cell_size=cell_size)
[docs]
def add_cylinder(self, r, h, origin=[0., 0., -200.], n_segs=16,
dip=0., dip_azimuth=0., cell_size=0., marker=None,
marker_position=None):
"""
Add cylinder-anomaly to **Omega**.
Required arguments
------------------
- r, type float
radius of the cylinder
- h, type float
height of the cylinder
Keyword arguments
-----------------
- origin = [0., 0., -200.], type list/array of three floats
center (of mass) coordinate of the cylinder; it is rotated in case
of non-zero dip or dip azimuth values around this point
- n_segs, type int
number of segments to approximate the circumference
- dip = 0., type float,
dip angle of the plate in degree
- dip_azimuth = 0., type float
dip azimuth angle in degree, zero points towards positive
x-direction
- cell_size = 0., type float
maximum cell-volume contstraint, default (0.) means no constraint.
- marker = None, type int
specify custom marker for plate if required
"""
print('... adding cylinder ...')
loops=[mu.loop_c(r=r, n_segs=n_segs, z=-h/2.),
mu.loop_c(r=r, n_segs=n_segs, z=h/2.)]
nns = []
for j in range(2):
if dip != 0.:
loops[j] = mu.rotate(loops[j], np.deg2rad(dip), axis='y')
if dip_azimuth != 0.:
loops[j] = mu.rotate(loops[j], np.deg2rad(dip_azimuth),
axis='z')
for jj in range(n_segs):
loops[j][jj, :] += origin
nc = self.Omega.nodeCount()
nns.append([self.Omega.createNodeWithCheck([o[0], o[1], o[2]]) for
o in loops[j]])
self.Omega.n_polys.append(np.arange(nc, nc + n_segs).tolist())
self.poly_marker.extend([23])
for j in range(n_segs - 1):
self.Omega.createQuadrangleFace(nns[0][j], nns[1][j],
nns[1][j+1], nns[0][j+1])
self.Omega.createQuadrangleFace(nns[0][-1], nns[1][-1],
nns[1][0], nns[0][0])
if marker_position is None:
marker_position = [origin[0], origin[1], origin[2]]
self.add_marker('cylinder', marker_position,
cell_size=cell_size, marker=marker)
[docs]
def add_plate(self, dx=200., dy=200., dz=50., origin=[0.0, 0.0, -200.0],
dip=0., dip_azimuth=0., cell_size=0., marker=None):
"""
Add rectangular *plate* anomaly to **Omega**. Rectangular bodys can be
added by using the method *add_brick()*.
Keyword arguments
-----------------
- dx, dy = 200., type float
dimensions in x/y-directions (before rotation)
- dz = 50., type float
thickness of the plate
- origin = [0., 0., -200.], type list/array of three floats
center (of mass) coordinate of the plate; plate is rotated in case
of non-zero dip or dip azimuth values around this point
- dip = 0., type float,
dip angle of the plate in degree
- dip_azimuth = 0., type float
dip azimuth angle in degree, zero points towards positive
x-direction
- cell_size = 0., type float
maximum cell-volume contstraint, default (0.) means no constraint.
- marker = None, type int
specify custom marker for plate if required
"""
print('... adding plate ... ')
x_lim = [-dx/2.0, dx/2.0]
y_lim = [-dy/2.0, dy/2.0]
zlim = [-dz/2.0, dz/2.0]
corners = np.array([[x_lim[0], y_lim[0], zlim[1]],
[x_lim[1], y_lim[0], zlim[1]],
[x_lim[1], y_lim[1], zlim[1]],
[x_lim[0], y_lim[1], zlim[1]],
[x_lim[0], y_lim[0], zlim[0]],
[x_lim[1], y_lim[0], zlim[0]],
[x_lim[1], y_lim[1], zlim[0]],
[x_lim[0], y_lim[1], zlim[0]]])
if dip != 0.:
corners = mu.rotate(corners, np.deg2rad(dip), axis='y')
if dip_azimuth != 0.:
corners = mu.rotate(corners, np.deg2rad(dip_azimuth), axis='z')
for j in range(len(corners)):
corners[j, :] += origin
nn = [self.Omega.createNodeWithCheck([o[0], o[1], o[2]])
for o in corners]
self.Omega.createQuadrangleFace(nn[0], nn[1], nn[2], nn[3])
self.Omega.createQuadrangleFace(nn[4], nn[5], nn[6], nn[7])
self.Omega.createQuadrangleFace(nn[0], nn[1], nn[5], nn[4])
self.Omega.createQuadrangleFace(nn[1], nn[2], nn[6], nn[5])
self.Omega.createQuadrangleFace(nn[2], nn[3], nn[7], nn[6])
self.Omega.createQuadrangleFace(nn[3], nn[0], nn[4], nn[7])
self.add_marker('plate', [origin[0], origin[1], origin[2]],
cell_size=cell_size, marker=marker)
[docs]
def add_points(self, points, marker=None):
"""
Add points to **Omega**.
Required arguments
------------------
- points, type array with shape (n, 3)
arrqay of points to be added
Keyword arguments
-----------------
- marker = None, type int
specify custom marker for points if required
"""
print('... adding points ... ')
if marker is None:
marker = [99] * len(points)
if type(marker) is int:
marker = [marker] * len(points)
nn = [self.Omega.createNodeWithCheck(
[o[0], o[1], o[2]]) for o in points]
if type(self.preserve_nodes) is list:
[n.setMarker(marker[j]) for j, n in enumerate(nn)]
[self.preserve_nodes.append(n) for n in nn]
[docs]
def add_air_tri(self, paths):
"""
Add refinement trianles to **Omega**.
Required arguments:
-------------------
- paths, type list
*paths* is a list of pointsets (arrays with shape (n, 3)) which
will be connected from point to point.
"""
print('... adding refinement triangles in airspace ... ')
for j in range(len(paths)):
self.air_tri.append(paths[j])
nn = [self.Omega.createNodeWithCheck([o[0], o[1], o[2]]) for
o in paths[j]]
if type(self.preserve_nodes) is list:
[n.setMarker(77) for n in nn]
[self.preserve_nodes.append(n) for n in nn]
if type(self.preserve_edges) is list:
for k in range(len(paths[j]) - 1):
self.Omega.createEdge(nn[k], nn[k + 1])
self.preserve_edges.append([nn[k].id(), nn[k + 1].id()])
self.edge_marker.extend(
[77] * (len(paths[j]) - 1))
[docs]
def add_paths(self, paths, closed_paths=None, marker=None):
"""
Add (observation) lines to **Omega**.
Required arguments:
-------------------
- paths, type list
*paths* is a list of pointsets (arrays with shape (n, 3)) which
will be connected from point to point.
Keyword arguments:
------------------
- closed_paths = None, type list
a list of same length as *paths* is required with *True* or *False*
values defining which path should be closed and which not.
"""
if closed_paths is None:
closed_paths = [False for elem in range(len(paths))]
if marker is None:
marker = [77] * len(paths)
if type(marker) is int:
marker = [marker] * len(paths)
print('... adding paths ... ')
for j in range(len(paths)):
nn = [self.Omega.createNodeWithCheck([o[0], o[1], o[2]]) for
o in paths[j]]
if type(self.preserve_nodes) is list:
[n.setMarker(marker[j]) for n in nn]
[self.preserve_nodes.append(n) for n in nn]
if type(self.preserve_edges) is list:
for k in range(len(paths[j]) - 1):
self.Omega.createEdge(nn[k], nn[k + 1])
self.preserve_edges.append([nn[k].id(), nn[k + 1].id()])
self.edge_marker.extend(
[marker[j]] * (len(paths[j]) - 1))
if closed_paths[j] is not False:
self.Omega.createEdge(nn[-1], nn[0])
self.preserve_edges.append([nn[-1].id(), nn[0].id()])
self.edge_marker.extend([marker[j]])
[docs]
def add_tx(self, tx, grounded=True, marker=None):
"""
Add transmitter paths to **Omega**.
Required arguments:
-------------------
- tx, type list
list of pointsets (arrays with shape (n, 3)) which
will be connected from point to point.
Keyword arguments:
------------------
- grounded = False, type bool
close first and last point of path to build a grounded loop or not.
"""
if marker is None:
marker = 77
nn = [self.Omega.createNodeWithCheck([o[0], o[1], o[2]]) for o in tx]
if type(self.preserve_nodes) is list:
[n.setMarker(marker) for n in nn]
[self.preserve_nodes.append(n) for n in nn]
if type(self.preserve_edges) is list:
for k in range(len(tx) - 1):
self.Omega.createEdge(nn[k], nn[k + 1])
self.preserve_edges.append([nn[k].id(), nn[k + 1].id()])
self.edge_marker.extend(
[marker] * (len(tx) - 1))
if not grounded:
self.Omega.createEdge(nn[-1], nn[0])
self.preserve_edges.append([nn[-1].id(), nn[0].id()])
self.edge_marker.extend([marker])
self.txs.append(tx)
self.grounding.append(grounded)
[docs]
def add_rx(self, points, marker=None):
"""
Add recevier points for automated interpolation to **Omega**.
Required arguments
------------------
- points, type list or array with shape (n, 3)
list of receiver points to add
"""
print('... adding points ... ')
if marker is None:
marker = [98] * len(points)
if type(marker) is int:
marker = [marker] * len(points)
self.rxs.append(points)
[docs]
def add_marker(self, label, pos, marker=None, cell_size=0.):
"""
Add region marker to domains. This function is called internally for
each layer or anomaly added. By default, domains are numbered
consequtively. It is usually not recommended not to maunaly set domain
marker, even though the possibility is given by overwriting the
keyword argument *marker=None* when adding, e.g., anomalies.
Required arguments:
-------------------
- label, type str
label for defining the domain information print output after
calling TetGen, which is enabled by default.
- pos, type list or array of length (3,)
x, y, z coordinates where the marker should be placed, note that
this point must be located inside the correct domain!
Keyword arguments:
------------------
- marker = None, type int
set manually defined region marker for this domain.
- cell_size = 0., type float
maximum cell-volume contstraint, default (0.) means no constraint.
"""
if marker is None:
marker = self.marker_counter
self.marker_counter += 1
self.region_info.append([label + ' marker:', marker,
'marker position:', [pos[0] + self.a_tol,
pos[1] + self.a_tol,
pos[2] + self.a_tol],
'volume_constraint:', cell_size])
self.Omega.addRegionMarker(pg.RVector3(pos[0] + self.a_tol,
pos[1] + self.a_tol,
pos[2] + self.a_tol),
marker, area=cell_size)
[docs]
def add_refinement_area(self, frame, quality, area, z=0., topo=False):
poly = pg.meshtools.createPolygon(
[point[:2] for point in frame],
isClosed=True)
refinement_mesh = pg.meshtools.createMesh(
poly, quality=quality, area=area)
refinement_mesh.setDimension(3)
coords = refinement_mesh.positions().array()
if not topo:
coords[:, 2] = z
else:
coords = self.add_topo(coords)
coords[:, 2] += z
for ci in range(len(coords)):
nn = [self.Omega.createNodeWithCheck(
[coords[ci, 0], coords[ci, 1], coords[ci, 2]])]
if type(self.preserve_nodes) is list:
[n.setMarker(marker[j]) for j, n in enumerate(nn)]
[self.preserve_nodes.append(n) for n in nn]
[docs]
def add_area_marker(self, poly, pos, face_size=0., marker=777):
"""
Add region marker to the surface or interfaces for 2D meshing, solely
used internally to define area-constraints for cooresponding meshes.
Required arguments:
-------------------
- poly, type pyGIMLi Polygon
polygon to which the marker is defined.
- pos, type list or array of length (3,)
x, y, z coordinates where the marker should be placed, note that
this point must be located inside the correct domain!
Keyword arguments:
------------------
- face_size = 0., type int
maximum facet-area contstraint, default (0.) means no constraint.
- marker = 777, type int
region marker to define sub-areas within the *poly*.
"""
poly.addRegionMarker(pg.RVector3(pos[0], pos[1], pos[2]),
marker, area=face_size)
[docs]
def update_kwargs(self, kwargs):
self.__dict__.update(kwargs)
for key in kwargs:
if key not in self.__dict__:
print('Warning! Unknown keyword argument set:', key)
print('... Continuing ...')
[docs]
def call_tetgen(self, tet_param='default', ignore_tx_z=False,
export_before=True, export_format=' -k',
name=None, copy=True, print_infos=True, suppress='NEF'):
"""
Wrapper-function to call TetGen version 1.5 and copy the resulting mesh
into the default *"m_dir"/_mesh* directory.
Keyword arguments:
------------------
- tet_param = 'default', type str
specify your parameters which should be passed to the TetGen call,
e.g., '**'-pq1.4aYAO2/3'** without including flags for the export.
Alternative keyword arguments are 'raw' ('**'-p'**) or 'faces'
('**'-pYA'**). It is also possible to directly pass the TetGen-
options (aside from export), e.g., *tet_param='-pq1.2aMAO2/3'*.
- export_before = True, type bool
automatically export the pyGIMLi mesh *Omega* as ".poly"-file in
the working directory **self.w_dir** before calling TetGen.
- export_vtk = True, type bool
set False for not exporting *.vtk*-file to be viewed in Paraview.
- name = None, type str
specify custom export name for the mesh, defualt is **self.name**.
- copy = True, type bool
automatically copy the mesh in *.mesh* (Medit)-format to the save-
directory **self.s_dir** for automated conversion on the first
computation run.
- print_infos = True, type bool
set *False* if no domain information should be printed at the end.
"""
if name is not None:
self.name = name
self.eval_unique_markers()
self.write_mesh_parameters(ignore_tx_z)
if os.path.isfile(self.w_dir + '/' + self.name + '.edge'):
os.remove(self.w_dir + '/' + self.name + '.edge')
if os.path.isfile(self.w_dir + '/' + self.name + '.node'):
os.remove(self.w_dir + '/' + self.name + '.node')
if export_before is True:
mu.export_tetgen_poly_file(
self.Omega, self.w_dir + '/' + self.name + '.poly',
poly_marker=self.poly_marker)
if self.export_omega:
self.Omega.exportVTK(self.w_dir + '/' + self.name +
'_omega.vtk')
if self.preserve_edges:
mu.export_tetgen_edge_file(
self.preserve_edges, self.edge_marker,
self.w_dir + '/' + self.name + '.edge')
if type(self.preserve_nodes) is list:
mu.export_tetgen_node_file(
self.preserve_nodes,
self.w_dir + '/' + self.name + '.node')
if tet_param == 'default':
tetgen_parameters = '-pq1.2aA'
elif tet_param == 'raw':
tetgen_parameters = "-p"
elif tet_param == 'faces':
tetgen_parameters = "-pMA"
else:
tetgen_parameters = tet_param
os.system("tetgen " + tetgen_parameters + export_format + suppress +
" " + self.w_dir + '/' + self.name + '.poly')
if copy:
# shutil.copy2(self.w_dir + '/' + self.name + '.1.mesh',
# self.s_dir + '/' + self.name + '.mesh')
shutil.copy2(self.w_dir + '/' + self.name + '.1.vtk',
self.m_dir + '/_vtk/' + self.name + '.vtk')
if print_infos is True:
self.print_region_info()
[docs]
def write_mesh_parameters(self, ignore_tx_z):
"""
Write topography and layered-earth structure information to a para-
meter file to preserve these information for following computations.
Also, tx and rx information is stored.
"""
RX = pg.Mesh(3)
for j, rx in enumerate(self.rxs):
if type(rx) is not list:
self.rxs[j] = rx.tolist()
for pos in self.rxs[j]:
RX.createNode(pos, marker=j)
RX.exportVTK(self.w_dir + '/' + self.name + '_rx.vtk')
TX = pg.Mesh(3)
txs = [tx for tx in self.txs]
for j, tx in enumerate(self.txs):
if type(tx) is not list:
self.txs[j] = tx.tolist()
for j, tx in enumerate(txs):
if not ignore_tx_z:
txs[j][:, 2] = self.get_topo_vals(txs[j])
nn = []
for ni, node in enumerate(txs[j]):
nn.append(TX.createNode(node))
if ni != 0:
TX.createEdge(nn[ni-1], nn[ni], marker=j+1)
if not self.grounding[j]:
TX.createEdge(nn[-1], nn[0], marker=j+1)
TX.exportVTK(self.w_dir + '/' + self.name + '_tx.vtk')
if callable(self.topo):
self.topo = self.topo.__name__
if self.bathymetry:
self.n_layers += 1
self.layer_depths = [0.] + self.layer_depths
container = dict()
container['topo'] = self.topo
container['centering'] = self.centering
container['easting_shift'] = self.easting_shift
container['northing_shift'] = self.northing_shift
container['rotation'] = self.rotation
container['n_layers'] = self.n_layers
container['layer_depths'] = self.layer_depths
container['zeroing'] = self.zeroing
container['bathymetry'] = self.bathymetry
container['rx'] = self.rxs
container['tx'] = self.txs
container['grounding'] = self.grounding
container['domains'] = self.domains
container['n_domains'] = self.n_domains
if hasattr(self, 'post_rotate'):
container['post_rotate'] = self.post_rotate
if hasattr(self, 'overwrite_markers'):
container['overwrite_markers'] = self.overwrite_markers
with open(self.p_dir + '/' + self.name + '_parameters.json',
'w') as outfile:
json.dump(container, outfile, indent=0)
if self.save_geometry:
np.savez(self.p_dir + '/' + self.name + '_geometry.npz',
txs=np.array(self.txs, dtype=object),
rxs=np.array(self.rxs, dtype=object),
ground_tri=self.ground_tri,
air_tri=self.air_tri,
topo = self.topo,
easting = self.easting_shift,
northing = self.northing_shift,
centering = self.centering,
rotation = self.rotation)
[docs]
def print_region_info(self):
"""
Print the domain information of the final mesh after calling TetGen.
"""
print('===========================================================')
print('### region information ###')
print('===========================================================')
for j in self.region_info:
print(j)
print('===========================================================')
print('### Done! ###')
print('===========================================================')
[docs]
def eval_unique_markers(self):
"""
extract unique markers specified in world *Omega*
"""
markers = []
for marker in self.Omega.regionMarkers():
markers.append(marker.marker())
self.domains = np.unique(markers).tolist()
self.n_domains = len(self.domains)
[docs]
def init_dem(self, surface_mesh_coords, topo=None):
"""
Initialize DEM instance for interpolating real-world topography from
digital-elevation-model-files and check if the extent of the dem-file
coordinates is large enough for the specified mesh dimensions.
Required arguments:
-------------------
- surface_mesh_coords, type array of shape (n, 3)
array containing all coordinates of the surface mesh for
comparing the spatial extent in x- and y-direction.
"""
if topo is None:
topo = self.topo
dem = DEM(self.t_dir + '/' + topo,
centering=self.centering,
easting_shift=self.easting_shift,
northing_shift=self.northing_shift)
if np.min(surface_mesh_coords[:, 0]) < np.min(dem.x) or \
np.max(surface_mesh_coords[:, 0]) > np.max(dem.x) or \
np.min(surface_mesh_coords[:, 1]) < np.min(dem.y) or \
np.max(surface_mesh_coords[:, 1]) > np.max(dem.y):
print('Fatal Error! The provided DEM file coordinates do not '
'cover the desired extension of the model domain '
'surface! Maybe this happend due to rotation or wrong '
'shifting of easting and northing? \nIf not, '
'either the model domain size must be decreased or the '
'area of the DEM file must be extended as extrapolation '
'is not allowed.')
print(np.min(surface_mesh_coords[:, 0]), np.min(dem.x))
print(np.max(surface_mesh_coords[:, 0]), np.max(dem.x))
print(np.min(surface_mesh_coords[:, 1]), np.min(dem.y))
print(np.max(surface_mesh_coords[:, 1]), np.max(dem.y))
raise SystemExit
else:
return(dem)
[docs]
def remove_duplicate_poly_faces(self):
"""
Check which *self.Omega.n_poly* faces were added multiple times if
subsurface structures touch each other and remove the duplicates.
"""
tmp, tmp_marker = [], []
max_poly_length = 0
for poly in self.Omega.n_polys:
if len(poly) > max_poly_length:
max_poly_length = len(poly)
ids = mu.get_unique_poly_ids(self.Omega.n_polys, max_poly_length)
for idx in ids:
tmp.append(self.Omega.n_polys[idx])
tmp_marker.append(self.poly_marker[idx])
self.Omega.n_polys = tmp
self.poly_marker = tmp_marker