# -*- coding: utf-8 -*-
"""
@author: Rochlitz.R
"""
import numpy as np
import os
"""
Auxiliary functions for meshgen_tools.
"""
[docs]
def pointset(pointset):
"""
Return a given set of 3D points, either list or array, as numpy array.
Required arguments
------------------
- pointset, type list
list of coordinates to be converted to an array of shape (:, 3)
"""
return(np.array(pointset))
[docs]
def line(start=[0., 0., 0.], stop=[0., 0., 0.], n_segs=100, z=0.,
topo=None, t_dir=None, **topo_kwargs):
"""
Return 3D coordinates of points on a line in arbitrary direction.
Keyword arguments
-----------------
- start, stop = [0., 0., 0.], type list or array
start and stop points
- n_segs = 100, type int
number of segments on the line
- topo_kwargs
listed below
- z = 0., type float
offset of the line in z-direction
- topo = None, type str or python function
DEM-topography file or synthetic topography (python) function.
- t_dir = None, type str
directory of topography (DEM) files
- 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
- revert_z = False, type bool
revert z-values of the DEM if the original height information was
provided in a donward-oriented coorindate system
- 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
"""
nodes = np.zeros((int(n_segs) + 1, 3))
nodes[:, 0] = np.linspace(start[0], stop[0], n_segs + 1)
nodes[:, 1] = np.linspace(start[1], stop[1], n_segs + 1)
nodes[:, 2] = np.linspace(start[2], stop[2], n_segs + 1)
return(assign_topography(nodes, z=z, topo=topo, t_dir=t_dir,
**topo_kwargs))
[docs]
def line_x(start=-1e2, stop=1e2, n_segs=100, y=0., z=0.,
topo=None, t_dir=None, **topo_kwargs):
"""
Return 3D coordinates of points on a line in x-direction.
Keyword arguments
-----------------
- start, stop = +-1e2, type float
start and stop values in x-direction
- n_segs = 100, type int
number of segments on the line
- y, z = 0., type float
offset of the line in y- and z-direction
- topo_kwargs
see description of *line* method or *meshgen_tools*/*dem_interpolator*
"""
nodes = np.zeros((int(n_segs) + 1, 3))
nodes[:, 0] = np.linspace(start, stop, int(n_segs) + 1)
nodes[:, 1] = y
return(assign_topography(nodes, z=z, topo=topo, t_dir=t_dir,
**topo_kwargs))
[docs]
def line_y(start=-1e2, stop=1e2, n_segs=100, x=0., z=0.,
topo=None, t_dir=None, **topo_kwargs):
"""
Return 3D coordinates of points on a line in y-direction.
Keyword arguments
-----------------
- start, stop = +-1e2, type float
start and stop values in x-direction
- n_segs = 100, type int
number of segments on the line
- x, z = 0., type float
offset of the line in x- and z-direction
- topo_kwargs
see description of *line* method or *meshgen_tools*/*dem_interpolator*
"""
nodes = np.zeros((int(n_segs) + 1, 3))
nodes[:, 0] = x
nodes[:, 1] = np.linspace(start, stop, n_segs + 1)
return(assign_topography(nodes, z=z, topo=topo, t_dir=t_dir,
**topo_kwargs))
[docs]
def line_z(start=-1e2, stop=1e2, n_segs=100, x=0., y=0.):
"""
Return 3D coordinates of points on a line in z-direction.
Keyword arguments
-----------------
- start, stop = +-1e2, type float
start and stop values in x-direction
- n_segs = 100, type int
number of segments on the line
- x, y = 0., type float
offset of the line in x- and y-direction
"""
nodes = np.zeros((int(n_segs) + 1, 3))
nodes[:, 0] = x
nodes[:, 1] = y
nodes[:, 2] = np.linspace(start, stop, n_segs + 1)
return(nodes)
[docs]
def loop_r(start=[-1e2, -1e2], stop=[1e2, 1e2], n_segs=4, nx=None, ny=None,
z=0., topo=None, t_dir=None, loop_rotation=None, **topo_kwargs):
"""
Return 3D coordinates of points on a rectangular loop.
Keyword arguments
-----------------
- start, stop = +-[1e2, 1e2], type list: [float, float]
opposite corner coordinates of rectangular loop
- n_segs = 4, type int
number of overall segments of the loop, this keyword argument is
ignored if *nx* and *ny* are specified
- nx, ny = 100, type int
number of segments on the x- and y-directed edges of the loop
- loop_rotation = None,
rotate loop horizontally around angle in (RAD)
- topo_kwargs
see description of *line* method or *meshgen_tools*/*dem_interpolator*
"""
if nx is None and ny is None:
peri_half = np.abs(start[0] - stop[0]) + np.abs(start[1] - stop[1])
ratio1 = np.abs(start[0] - stop[0]) / peri_half
ratio2 = np.abs(start[1] - stop[1]) / peri_half
nx_segs = int(np.round((n_segs/2.) * ratio1))
ny_segs = int(np.round((n_segs/2.) * ratio2))
if nx_segs == 0:
nx_segs = 1
if ny_segs == 0:
ny_segs = 1
else:
nx_segs, ny_segs, = nx, ny
nodes = line_x(start[0], stop[0], n_segs=nx_segs, y=start[1],
topo=topo, t_dir=t_dir)[:-1, :]
nodes = np.vstack((nodes, line_y(start[1], stop[1], n_segs=ny_segs,
x=stop[0], topo=topo,
t_dir=t_dir)[:-1, :]))
nodes = np.vstack((nodes, line_x(stop[0], start[0], n_segs=nx_segs,
y=stop[1], topo=topo,
t_dir=t_dir)[:-1, :]))
nodes = np.vstack((nodes, line_y(stop[1], start[1], n_segs=ny_segs,
x=start[0], topo=topo,
t_dir=t_dir)[:-1, :]))
if loop_rotation is not None:
shft = [np.mean(nodes[:, 0]), np.mean(nodes[:, 1]), 0.]
nodes = rotate_around_point(nodes, shft,
[np.deg2rad(loop_rotation)], ['z'])
return(assign_topography(nodes, z=z, topo=topo, t_dir=t_dir,
**topo_kwargs))
[docs]
def loop_c(origin=[0., 0., 0.], r=5e1, n_segs=100, z=0.,
topo=None, t_dir=None, **topo_kwargs):
"""
Return 3D coordinates of points on a circular loop.
Keyword arguments
-----------------
- origin = [0., 0., 0.], type list of floats
x-, y-, z- coordinates of the center of the loop
- r = 5e1, type float
radius of the loop
- n_segs = 100, type int
number of segments on the loop circumference
- topo_kwargs
see description of *line* method or *meshgen_tools*/*dem_interpolator*
"""
nodes = create_circle(r, n_segs, z)
for j in range(len(nodes)):
nodes[j, :] += origin
return(assign_topography(nodes, z=z, topo=topo, t_dir=t_dir,
**topo_kwargs))
[docs]
def loop_e(origin=[0., 0., 0.], a=5e1, b=5e1, n_segs=100, z=0.,
topo=None, t_dir=None, **topo_kwargs):
"""
Return 3D coordinates of points on a ellipsoidal loop.
Keyword arguments
-----------------
- origin = [0., 0., 0.], type list of floats
x-, y-, z- coordinates of the center of the loop
- a, b = 5e1, type float
radius of the loop
- n_segs = 100, type int
number of segments on the loop circumference
- topo_kwargs
see description of *line* method or *meshgen_tools*/*dem_interpolator*
"""
nodes = create_ellipse(a, b, n_segs, z)
for j in range(len(nodes)):
nodes[j, :] += origin
return(assign_topography(nodes, z=z, topo=topo, t_dir=t_dir,
**topo_kwargs))
[docs]
def rx_grid(xlinsp, ylinsp, z=0.):
"""
Build Rx grid from two orthogonal vectors.
Required arguments
------------------
- xvec/yvec, type list/array of len 3
start, stop and number of nodes for *linspace* method, building the
vectors in x and y direction.
"""
# define vectors
xr = np.linspace(xlinsp[0], xlinsp[1], int(xlinsp[2]))
yr = np.linspace(ylinsp[0], ylinsp[1], int(ylinsp[2]))
grid = np.zeros((int(xlinsp[2] * ylinsp[2]), 3))
# create grid
a = 0
for x in xr:
for y in yr:
grid[a, 0] = x
grid[a, 1] = y
a += 1
grid[:, 2] = z
return(grid)
[docs]
def rotate(array, alpha, axis='z', tensor=False):
"""
Rotate array about angle alpha(RAD) around given axis
Default rotation axis is z, alternatively x or y.
Required arguments
------------------
- array, type array of shape (:, 3)
array of coordinates to be rotated
- alpha, type float
angle(RAD) specifying the rotation
Keyword arguments
-----------------
- axis = 'z', type str
rotation axis
"""
if axis == 'z':
rotation_matrix = np.array(([np.cos(-alpha), -np.sin(-alpha), 0],
[np.sin(-alpha), np.cos(-alpha), 0],
[0, 0, 1]), dtype=np.float64)
elif axis == 'x':
rotation_matrix = np.array(([1, 0, 0],
[0, np.cos(-alpha), -np.sin(-alpha)],
[0, np.sin(-alpha), np.cos(-alpha)]),
dtype=np.float64)
elif axis == 'y':
rotation_matrix = np.array(([np.cos(-alpha), 0, np.sin(-alpha)],
[0, 1, 0],
[-np.sin(-alpha), 0, np.cos(-alpha)]),
dtype=np.float64)
if not tensor:
for i in range(array.shape[0]):
array[i] = array[i].dot(rotation_matrix)
else:
array = rotation_matrix.dot(array.dot(rotation_matrix.T))
return(array)
[docs]
def translate(array, shift, back=False):
"""
Translate array about given shift
Required arguments
------------------
- array, type array with shape (:, 3)
array of coordinates to ne translated
- shift, type list or array of length 3
x-, y-, and z- values defining the shift
Keyword arguments
-----------------
- back = False, type bool
set "True" is shift in opposite direction is required
"""
if not back:
return(array - np.array(shift).reshape(1, 3))
else:
return(array + np.array(shift).reshape(1, 3))
[docs]
def rotate_around_point(array, shift, angles, axes, pre_rotate=False):
"""
Rotate array around arbitrary point defined by shift, not around centroid.
Required arguments
------------------
- array, type array of shape (:, 3)
array of coordinates to be rotated
- shift, type list or array of length 3
x-, y-, and z- values defining the shift
- angles, type list of floats
list of angles(RAD) specifying the rotation
- axes, type list of str
list of rotation axes
Keyword arguments
-----------------
- pre_rotate = False, type bool
conduct an initial roation using the second entry of axes/angles
"""
if len(angles) != len(axes):
print('Error! Length of lists of *angles* and rotation *axes* must be '
'the same!')
raise SystemExit
array = translate(array, shift)
if pre_rotate:
array = rotate(array, -angles[1], axes[1])
array[:, 1] *= np.cos(angles[0])
for jjj in range(len(angles)):
array = rotate(array, angles[jjj], axes[jjj])
return(translate(array, shift, back=True))
[docs]
def create_circle(r, n, h):
"""
Create a list of **n** (number of) 3D coordinates on the circumference of a
circle with radius **r** and height **h**.
Required arguments
------------------
- r, type float
the radius of the circle
- n, type int
number of nodes on the cirumference of the circle
- h, type float
z-coordinates of the circle
"""
return(np.array([(np.cos(2 * np.pi / n * x)*r,
np.sin(2 * np.pi / n * x)*r,
h) for x in range(0, n)]))
[docs]
def create_ellipse(a, b, n, h):
"""
Create a list of **n** (number of) 3D coordinates on the circumference of
an ellipsum with 'radii' **a** and **b**, and height **h**.
Required arguments
------------------
- a, b, type float
the two radius of the ellipse
- n, type int
number of nodes on the cirumference of the ellipse
- h, type float
z-coordinates of the ellipse
"""
def ellipse(t, a, b):
return a*np.cos(t), b*np.sin(t)
create_ellipseints = [ellipse(t, a, b) for t in np.linspace(0, 2*np.pi, n)]
x, y = [np.array(v) for v in list(zip(*create_ellipseints))]
xy = np.hstack((x.reshape(n, 1), y.reshape(n, 1)))
return(np.hstack((xy, np.ones((n, 1)) * h)))
[docs]
def poly_area(path):
"""
Caluclate enclosed area of polygon.
Required arguments
------------------
- path, type array with shape (:, 2) or (:, 3)
array of coordinates describing the polygon
"""
n = len(path)
area = 0.0
for i in range(n):
j = (i + 1) % n
area += path[i][0] * path[j][1]
area -= path[j][0] * path[i][1]
area = abs(area) / 2.0
return(area)
[docs]
def poly_peri(path):
"""
Calculate perimeter of polygon.
Required arguments
------------------
- path, type array with shape (:, 2) or (:, 3)
array of coordinates describing the polygon
"""
n = len(path)
peri = 0.0
for i in range(n - 1):
peri += np.sqrt((path[i+1, 0] - path[i, 0])**2 +
(path[i+1, 1] - path[i, 1])**2)
peri += np.sqrt((path[-1, 0] - path[0, 0])**2 +
(path[-1, 1] - path[0, 1])**2)
return(peri)
[docs]
def find_edge_node_ids(frame, mesh, id_offset):
"""
Internal utility function for meshgen tools to find the edge node ids of a
horizontal 2D mesh for the extension of the mesh in z-direction.
Required arguments
------------------
- frame, type array with shape (:, 3)
references frame defined by coordinates
- mesh, type pyGIMLi mesh
mesh, whose coordinates should be checked
- id_offset, type int
offset of ids between 2D surface meshes and 3D mesh (should be removed
if the meshgen_tools get an overhaul)
"""
ids = []
frame = np.vstack((frame, frame[0]))
for j in range(len(frame) - 1):
tmp_idx, tmp_pos, dist = [], [], []
for idx, pos in enumerate(mesh.positions().array()):
if is_between_2D(frame[j, :], frame[j+1, :], pos):
tmp_idx.append(idx)
tmp_pos.append(pos[:2])
tmp_pos = np.array(tmp_pos)
for jj in range(len(tmp_pos)):
dist.append(np.linalg.norm(frame[j, :] - tmp_pos[jj]))
tmp_list = [sx for _, sx in sorted(zip(dist, tmp_idx))]
ids.append(np.array(tmp_list) + id_offset)
return(ids)
[docs]
def find_closest(l1, l2):
"""
Find the closest entry in a list of 1D corrdinates in comparison to all
entries of another list.
Required arguments
------------------
- l1, type list
reference list of coordinates
- l2, type list
another list of coordinate, usually with a different length
"""
ll = []
for k in range(len(l1)):
ll.append(min(range(len(l2)), key=lambda i: abs(l2[i]-l1[k])))
return(ll)
[docs]
def resolve_rx_overlaps(rxs, refinement_size=10., ignore_z=True):
allrx = rxs[0]
for ri in range(1, len(rxs)):
allrx = np.vstack((allrx, rxs[ri]))
# go over x-direction
allrx = allrx[np.argsort(allrx[:, 0]), :]
to_keep = []
skip = []
for ii in range(len(allrx)):
if ignore_z:
dists = np.linalg.norm(allrx[ii, :2] - allrx[:, :2], axis=1)
else:
dists = np.linalg.norm(allrx[ii, :] - allrx[:, :], axis=1)
where = np.argwhere(dists < 2.*refinement_size).flatten()
skip.extend(where[where>ii])
if ii not in skip:
to_keep.append(ii)
return(allrx[to_keep])
[docs]
def closest_node(coord, coords):
"""
Find the closest coordinate in a list of 2D/3D corrdinates.
Required arguments
------------------
- node, list/array with length 3
reference node
- nodes, list/array with shape (:, 2) or (:, 3)
list or array of nodes to be checked
"""
coords = np.asarray(coords)
dist_2 = np.sum((coords - coord)**2, axis=1)
return(np.argmin(dist_2))
[docs]
def get_closest(coord, coords, n=None, max_dist=1.):
"""
Find the ids the closest coordinates in a list of 2D/3D corrdinates,
either using a distance threshold or by specifying the number of closests
coordinates.
Required arguments
------------------
- coord, list/array with length 3
reference node
- coords, list/array with shape (:, 2) or (:, 3)
list or array of nodes to be checked
Keyword arguments
-----------------
- n = None, type int
number of closests coordinates, overwrites *max_dist*
- max_dist = 1., type float
return all coordinate ids which are closer than *max_dists*
"""
ids = []
dist = np.sqrt(np.sum((coords - coord)**2, axis=1))
if n is not None:
return(np.argsort(dist)[:n])
else:
return(np.argsort(dist)[dist<=max_dist])
[docs]
def add_edges(surface, edges, nn, size, marker, closed=False):
"""
Add edge to polygon.
"""
for k in range(size - 1):
if sorted([nn[k].id(), nn[k+1].id()]) not in edges:
edges.append(sorted([nn[k].id(), nn[k+1].id()]))
surface.createEdge(nn[k], nn[k + 1], marker=marker)
if closed:
if sorted([nn[-1].id(), nn[0].id()]) not in edges:
edges.append(sorted([nn[-1].id(), nn[0].id()]))
surface.createEdge(nn[-1], nn[0], marker=marker)
return(surface, edges)
[docs]
def find_on_frame(frame_coords, coords):
"""
Check if coordinates are located on a frame.
Required arguments
------------------
- frame_coords, type array with shape (:, 3)
array of frame coordinates
- coords, type array with shape (:, 3)
array of coordinates to be checked
"""
frame = []
frame_coords = np.vstack((frame_coords, frame_coords[0, :]))
for j in range(len(frame_coords)-1):
tmp = []
for k, coord in enumerate(coords):
if is_between_2D(frame_coords[j], frame_coords[j+1], coord):
if k not in frame:
tmp.extend([k])
tmp_ordered = order_list(frame_coords[j], coords[tmp], tmp)
frame.extend(tmp_ordered)
return(frame)
[docs]
def find_nodes_on_segment(p0, p1, mesh, topo_f):
"""
Find nodes on a segment.
Required arguments
------------------
- mesh
"""
ids = []
tmp_idx, tmp_pos, dist = [], [], []
for idx, pos in enumerate(mesh.positions().array()):
if is_between_2D(p0, p1, pos):
tmp_idx.append(idx)
if topo_f is None:
if np.isclose(pos[2], 0., atol=1e-2):
tmp_pos.append(pos)
else:
if np.isclose(pos[2], topo_f(pos[0], pos[1]), atol=1e-2):
tmp_pos.append(pos)
tmp_pos = np.array(tmp_pos)
for jj in range(len(tmp_pos)):
dist.append(np.linalg.norm(p0 - tmp_pos[jj]))
return([sx for _, sx in sorted(zip(dist, tmp_idx))])
[docs]
def order_list(ref, points, tmp):
"""
Sort values in list in increasing order. Deprecated and not used anymore.
"""
dist = np.sqrt((points[:, 0] - ref[0])**2 + (points[:, 1] - ref[1])**2)
return(np.array(np.array(tmp)[dist.argsort()], dtype=int))
[docs]
def get_intersect(a1, a2, b1, b2):
"""
Return intersection coordinate of two segments *a* and *b*.
"""
T = np.array([[0, -1], [1, 0]])
da = np.atleast_2d(a2 - a1)
db = np.atleast_2d(b2 - b1)
dp = np.atleast_2d(a1 - b1)
dap = np.dot(da, T)
denom = np.sum(dap * db, axis=1)
if denom == 0.:
return(None)
else:
num = np.sum(dap * dp, axis=1)
isect = np.atleast_2d(num / denom).T * db + b1
if is_between_2D(a1, a2, isect[0]):
return(isect)
else:
return(None)
[docs]
def is_between(a, b, c, a_tol=1e-2, print_points=False):
"""
Search if a Nedelec-dof is located on the edge between two points in 3D,
needed for crooked sources defined as path.
Search if point c is in between a and b.
Required arguments:
-------------------
- a, b, c, type list/array of length 3
three points to be evaluated
Keyword arguments:
------------------
- a_tol = 1e-2, type float
absolute tolerance to account for numerical inaccuracies
- print_points = False, type bool
set **True** to enable print of coordinates
"""
if print_points:
print(a, b, c)
if (a[0] - a_tol <= c[0] <= b[0] + a_tol or
a[0] + a_tol >= c[0] >= b[0] - a_tol):
if ((a[1] - a_tol <= c[1] <= b[1] + a_tol or
a[1] + a_tol >= c[1] >= b[1] - a_tol)):
if ((a[2] - a_tol <= c[2] <= b[2] + a_tol or
a[2] + a_tol >= c[2] >= b[2] - a_tol)):
if np.isclose((c[1] - a[1]) * (b[0] - a[0]),
(c[0] - a[0]) * (b[1] - a[1]), a_tol):
return(True)
return(False)
[docs]
def is_between_1D(a, b, c, a_tol=1e-2, print_points=False):
"""
Search if a x-coordinate is located in an interval in
in horizontal direction. Search if point c is in between a and b.
Required arguments
------------------
- a, b, c, type list/array of length 3
three points to be evaluated
Keyword arguments
-----------------
- a_tol = 1e-2, type float
absolute tolerance to account for numerical inaccuracies
- print_points = False, type bool
set **True** to enable print of coordinates
"""
if print_points:
print(a, b, c)
if (a[0] < c[0] < b[0] or a[0] > c[0] > b[0]):
return(True)
return(False)
[docs]
def is_between_2D(a, b, c, a_tol=1e-2, print_points=False):
"""
Search if a Nedelec-dof is located on the edge between two points in 2D
in horizontal direction. Search if point c is in between a and b.
Required arguments
------------------
- a, b, c, type list/array of length 3
three points to be evaluated
Keyword arguments
-----------------
- a_tol = 1e-2, type float
absolute tolerance to account for numerical inaccuracies
- print_points = False, type bool
set **True** to enable print of coordinates
"""
if print_points:
print(a, b, c)
if (a[0] - a_tol <= c[0] <= b[0] + a_tol or
a[0] + a_tol >= c[0] >= b[0] - a_tol):
if (a[1] - a_tol <= c[1] <= b[1] + a_tol or
a[1] + a_tol >= c[1] >= b[1] - a_tol):
if np.isclose((c[1] - a[1]) * (b[0] - a[0]),
(c[0] - a[0]) * (b[1] - a[1]), a_tol):
return(True)
return(False)
[docs]
def is_between_2D_exact(a, b, c, a_tol=1e-2, print_points=False):
"""
Search if a Nedelec-dof is located on the edge between two points in 2D
in horizontal direction. Search if point c is in between a and b.
Required arguments
------------------
- a, b, c, type list/array of length 3
three points to be evaluated
Keyword arguments
-----------------
- a_tol = 1e-2, type float
absolute tolerance to account for numerical inaccuracies
- print_points = False, type bool
set **True** to enable print of coordinates
"""
if print_points:
print(a, b, c)
if (a[0] < c[0] < b[0] or a[0] > c[0] > b[0]):
if (a[1] < c[1] < b[1] or a[1] > c[1] > b[1]):
if np.isclose((c[1] - a[1]) * (b[0] - a[0]),
(c[0] - a[0]) * (b[1] - a[1]), a_tol):
return(True)
return(False)
[docs]
def is_between_2Dvert(a, b, c, a_tol=1e-2, print_points=False):
"""
Search if a Nedelec-dof is located on the edge between two points in 2D
in vertical direction. Search if point c is in between a and b.
Required arguments
------------------
- a, b, c, type list/array of length 3
three points to be evaluated
Keyword arguments
-----------------
- a_tol = 1e-2, type float
absolute tolerance to account for numerical inaccuracies
- print_points = False, type bool
set **True** to enable print of coordinates
"""
if print_points:
print(a, b, c)
if (a[0] - a_tol <= c[0] <= b[0] + a_tol or
a[0] + a_tol >= c[0] >= b[0] - a_tol):
if (a[2] - a_tol <= c[2] <= b[2] + a_tol or
a[2] + a_tol >= c[2] >= b[2] - a_tol):
if np.isclose(c[1], a[1], a_tol) and np.isclose(c[1], b[1], a_tol):
if np.isclose((c[2] - a[2]) * (b[0] - a[0]),
(c[0] - a[0]) * (b[2] - a[2]), a_tol):
return(True)
return(False)
[docs]
def write_synth_topo_to_asc(t_dir, file_name, topo, spacing=25.,
x_min=-1e4, y_min=-1e4):
"""
Write an array of topography data to an *.asc* file.
Required arguments
------------------
- t_dir, type str
name of output directory
- file_name, type str
name of exported topography file
- topo, type array with shape (:, 3)
2D array containing only z_coordinates with topography information
Keyword arguments:
------------------
- spacing = 25., type float
horizontal grid spacing in x- and y-direction
- x_min = -1e4, type float
x-coordinate of "lower left" corner
- y_min = -1e4, type float
y-coordiante of "lower_left" corner
"""
with open(t_dir + '/' + file_name + '.asc', "w") as demfile:
demfile.write('ncols ' + str(topo.shape[1]) + '\n'
'nrows ' + str(topo.shape[0]) + '\n'
'xllcorner ' + str(x_min) + '\n'
'yllcorner ' + str(y_min) + '\n'
'cellsize ' + str(spacing) + '\n'
'NODATA_value -9999\n')
with open(t_dir + '/' + file_name + '.asc', "ab") as demfile:
np.savetxt(demfile, topo, '%3.8f')
[docs]
def write_synth_topo_to_xyz(t_dir, file_name, topo):
"""
Write an array of topography data to an *.xyz* file.
Required arguments
------------------
- t_dir, type str
name of output directory
- file_name, type str
name of exported topography file
- topo, type array with shape (:, 3)
array containing coordinates with topography information
"""
with open(t_dir + '/' + file_name + '.xyz', "wb") as demfile:
np.savetxt(demfile, topo)
[docs]
def assign_topography(nodes, t_dir=None, topo=None, z=0., centering=True,
easting_shift=None, northing_shift=None,
revert_z=False, rotation=None):
"""
Assign topography on nodes.
Required arguments
------------------
- nodes, type list/array with shape (:, 3)
list or array of 3D coordinates
Keyword arguments
-----------------
- t_dir = None, type str
directory of topography file
- topo = None, type function or str
topography function or file
- z = 0., type float
add constant offset to topography in z-direction
- centering = True, type bool
Set **True** if extent of topography file should be centered horizontally
around origin (0, 0)
- easting_shift = None, type float
add specified shift to x-coordinates of topography file
- northing_shift = None, type float
add specified shift to y-coordinates of topography file
- rotation = None, type float
rotate complete topography file about specified angle(DEG) horizontally
"""
if topo is None:
nodes[:, 2] = z
else:
if type(topo) is not str:
nodes[:, 2] = topo(nodes[:, 0], nodes[:, 1])
else:
from custEM.meshgen.dem_interpolator import DEM
dem = DEM(t_dir + '/' + topo, centering=bool(centering),
easting_shift=easting_shift,
northing_shift=northing_shift, revert_z=revert_z)
nodes[:, 2] = dem(nodes[:, 0], nodes[:, 1], rotation=rotation)
nodes[:, 2] += z
return(nodes)
[docs]
def inside_poly(x, y, path):
"""
Return True if a coordinate (x, y) is inside a polygon defined by
a list of verticies [(x1, y1), (x2, x2), ... , (xN, yN)].
Reference: http://www.ariel.com.au/a/python-point-int-poly.html
Required arguments
------------------
- x, type float
x coordinate of point to be checked
- y, type float
y coordinate of point to be checked
- path, type list/array with shape (:, 3)
list or array of 3D coordinates
"""
n = len(path)
inside = False
p1x, p1y = path[0]
for i in range(1, n + 1):
p2x, p2y = path[i % n]
if y > min(p1y, p2y):
if y <= max(p1y, p2y):
if x <= max(p1x, p2x):
if p1y != p2y:
xinters = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
if p1x == p2x or x <= xinters:
inside = not inside
p1x, p1y = p2x, p2y
return(inside)
[docs]
def refine_path(path, n_segs=None, length=None):
"""
Refine each section of a 3D path into n_segs equal segments.
Required arguments
------------------
- path, type list/array with shape (:, 3)
list or array of 3D coordinates
- n_segs, type int
number of nodes to be inserted within each segment of the path
- length, type float
maximum length of any segment of the path
"""
if n_segs is not None and length is not None:
print('Error! Specify either *n_segs* (each segment of the path is '
'split into n_segs or a maximum segment *length* for the path.')
raise SystemExit
new_path = np.array(path[0]).reshape(1, 3)
for j in range(len(path) - 1):
if length is not None:
n_segs = int(np.linalg.norm(path[j] - path[j+1]) / length) + 1
refined = np.array(
[np.linspace(path[j, 0], path[j+1, 0], n_segs + 1),
np.linspace(path[j, 1], path[j+1, 1], n_segs + 1),
np.linspace(path[j, 2], path[j+1, 2], n_segs + 1)])
new_path = np.vstack((new_path, refined.T[1:, :]))
return(new_path)
[docs]
def refine_rx(coords, r=5., rot=None, n_segs=3):
"""
Refine rx position by enclosing them with hexagons of radius *r*.
The smaller *r*, the better the refinement. Note that each hexagon will be
subdivided into 6 equilateral triangles, which will support optimum
tetrahedra quality during the meshing process.
Required arguments
------------------
- coords, type list/array with shape (:, 3)
list or array of 3D coordinates
Keyword arguments
-----------------
- r = 5., type float
radius of surrounding equilateral (triangle if *n_segs=3*) shape
- rot = None, type float
horizontal rotation angle(DEG) between subsequent shapes
- n_segs = 3, type int
number of corners of surrounding equilateral polygon shape
"""
rx_refined = []
for rec in coords:
rx_refined.append(loop_c(origin=rec, r=r, z=rec[2], n_segs=n_segs))
if rot is not None:
degs = np.arange(len(coords)) * rot
for j in range(len(coords)):
rx_refined[j] = rotate_around_point(
rx_refined[j], coords[j], [np.deg2rad(degs[j])], ['z'])
return(rx_refined)
[docs]
def refine_adaptive(coords, txs, r1=20., r2=5., r3=1., d3=200., d2=500.,
min_tx_dist=100., rot=30., n_segs=3):
"""
Decription to add
"""
rx_refined = []
new_coords = []
rx_single = [[] for ti in txs]
for pos in coords:
dists = []
for tx in txs:
dists.append(np.min(np.linalg.norm(pos[:2] - tx[:, :2], axis=1)))
dist = min(dists)
if dist < min_tx_dist:
#pass
for ti, tx in enumerate(txs):
if dists[ti] >= min_tx_dist:
rx_single[ti].append(list(pos))
rx_refined.append(loop_c(origin=pos, r=r1, z=pos[2],
n_segs=n_segs))
elif dist < d3:
new_coords.append(list(pos))
rx_refined.append(loop_c(origin=pos, r=r3, z=pos[2],
n_segs=n_segs))
elif dist < d2:
new_coords.append(list(pos))
rx_refined.append(loop_c(origin=pos, r=r2, z=pos[2],
n_segs=n_segs))
else:
new_coords.append(list(pos))
rx_refined.append(loop_c(origin=pos, r=r1, z=pos[2],
n_segs=n_segs))
if rot is not None:
degs = np.arange(len(rx_refined)) * rot
for j in range(len(rx_refined)):
rx_refined[j] = rotate_around_point(
rx_refined[j], np.mean(rx_refined[j], axis=0),
[np.deg2rad(degs[j])], ['z'])
return(new_coords, rx_single, rx_refined)
[docs]
def npoly_to_triangles(world, ids):
"""
Add polygones with more than three edges in form of triangles to *Omega*.
"""
for idx in range(len(ids) - 2):
world.createTriangleFace(
world.node(ids[idx]),
world.node(ids[idx+1]),
world.node(ids[-1]))
[docs]
def get_unique_poly_ids(polys, max_poly_length):
"""
Get unique ids of nodes of multiple polygons.
Required arguments
------------------
- polys, type list
list of polygons
- max_poly_length
maximum number of nodes of within all polygons
"""
arr = np.zeros((len(polys), max_poly_length), dtype=int)
for k, ele in enumerate(polys):
arr[k, :len(ele)] = sorted(ele)
_, ids = np.unique(arr, axis=0, return_index=True)
return(ids)
[docs]
def export_tetgen_poly_file(poly, filename, float_format='.3f',
poly_marker=None, **kwargs):
"""
Write a given piecewise linear complex (mesh/poly) into an Ascii file in
Tetgen *.poly* format.
Required arguments
------------------
- poly, type pyGIMLi mesh
piecewise linear complex which should be exported to the *.poly* file
- filename, type str
name of output file
Keyword arguments
-----------------
float_format = '.12e', type string
format that will be used to write float values in the Ascii file,
default is the exponential float form with a precision of 12 digits
"""
if filename[-5:] != '.poly':
filename = filename + '.poly'
polytxt = ''
sep = '\t' # standard tab seperated file
assert poly.dim() == 3, 'Exit, only for 3D meshes.'
boundary_marker = 1
attribute_count = 0
# Part 1/4: node list
# intro line
# <nodecount> <dimension (3)> <# of attributes> <boundary marker (0 or 1)>
polytxt += '{0}{5}{1}{5}{2}{5}{3}{4}'.format(poly.nodeCount(), 3,
attribute_count,
boundary_marker,
os.linesep, sep)
# loop over positions, attributes and marker(node)
# <point idx> <x> <y> <z> [attributes] [boundary marker]
point_str = '{:d}' # index of the point
for i in range(3):
# coords as float with given precision
point_str += sep + '{:%s}' % (float_format)
point_str += sep + '{:d}' + os.linesep # node marker
for j, node in enumerate(poly.nodes()):
fill = [node.id() + 1]
fill.extend([pos for pos in node.pos()])
fill.append(node.marker())
polytxt += point_str.format(*fill)
# Part 2/4: boundary list
# intro line
# <# of facets> <boundary marker (0 or 1)>
nBoundaries = poly.boundaryCount()
# look for extra boundaries present in either the PLC or in kwargs
n_polys = []
if 'n_polys' in kwargs:
n_polys += kwargs.pop('n_polys', [])
if hasattr(poly, 'n_polys'):
n_polys += poly.n_polys
nBoundaries += len(n_polys)
polytxt += '{0:d}{2}1{1}'.format(nBoundaries, os.linesep, sep)
# loop over facets, each facet can contain an arbitrary number of holes
# and polygons, in our case, there is always one polygon per facet.
for bound in poly.boundaries():
# one line per facet
# <# of polygons> [# of holes] [boundary marker]
npolys = 1
polytxt += '1{2}0{2}{0:d}{1}'.format(bound.marker(), os.linesep, sep)
# inner loop over polygons
# <# of corners> <corner 1> <corner 2> ... <corner #>
for l in range(npolys):
poly_str = '{:d}'.format(bound.nodeCount())
for ind in bound.ids():
poly_str += sep + '{:d}'.format(ind + 1)
polytxt += '{0}{1}'.format(poly_str, os.linesep)
# inner loop over holes
# not necessary yet ?! why is there an extra hole section?
# because this is for 2D holes in facets only
# part 2b: extra boundaries that cannot be part of mesh class
if poly_marker is None:
poly_marker = [111] * len(n_polys)
for nn, nodes in enumerate(n_polys):
# <# of polygons> [# of holes] [boundary marker]
npolys = 1
polytxt += '1{2}0{2}{0:d}{1}'.format(poly_marker[nn], os.linesep, sep)
# <# of corners> <corner 1> <corner 2> ... <corner #>
poly_str = '{:d}'.format(len(nodes))
for ind in nodes:
poly_str += sep + '{:d}'.format(ind + 1)
polytxt += '{0}{1}'.format(poly_str, os.linesep)
# part 3/4: hole list
# intro line
# <# of holes>
holes = poly.holeMarker()
polytxt += '{:d}{}'.format(len(holes), os.linesep)
# loop over hole marker
# <hole #> <x> <y> <z>
hole_str = '{:d}'
for m in range(3):
hole_str += sep + '{:%s}' % float_format
hole_str += os.linesep
for n, hole in enumerate(holes):
polytxt += hole_str.format(n, *hole)
# part 4/4: region attributes and volume constraints (optional)
# intro line
# <# of regions>
try:
regions = poly.regionMarker()
except Exception:
regions = poly.regionMarkers()
polytxt += '{:d}{}'.format(len(regions), os.linesep)
# loop over region marker
# <region #> <x> <y> <z> <region number> <region attribute>
region_str = '{:d}'
for o in range(3):
region_str += sep + '{:%s}' % (float_format)
region_str += sep + '{:d}%s{:%s}' % (sep, float_format) + os.linesep
for p, region in enumerate(regions):
polytxt += region_str.format(p, region.x(), region.y(), region.z(),
region.marker(),
region.area())
# writing file
with open(filename, 'w') as out:
out.write(polytxt)
[docs]
def export_tetgen_edge_file(preserve_edges, edge_marker, filename):
"""
Write a given set of edges which should be preserved into a Ascii file in
Tetgen *.edge* format.
Required arguments
------------------
- preserve_edges, type list
edges which should be exported to the *.edge* file
- edge_marker, type list
list of markers assigned to the edges
- filename, type str
name of output file
"""
edgetxt = ''
# intro line
# <edgecount> <boundary marker (0 or 1)>
edgetxt += '{} {}{}'.format(len(preserve_edges), 1, os.linesep)
for jj, edge in enumerate(preserve_edges):
edgetxt += '{:d} {:d} {:d} {:d}{}'.format(
jj + 1, edge[0] + 1, edge[1] + 1, edge_marker[jj], os.linesep)
with open(filename, 'w') as out:
out.write(edgetxt)
[docs]
def export_tetgen_node_file(preserve_nodes, filename):
"""
Write a given set of nodes which should be preserved into an Ascii file in
Tetgen *.edge* format.
Required arguments
------------------
- preserve_nodes, type list
edges which should be exported to the *.node* file
- filename, type str
name of output file
"""
nodetxt = ''
# intro line
# <nodecount> <boundary markers (0 or 1)>
nodetxt += '{} {}{}'.format(len(preserve_nodes), 0, 1, os.linesep)
for jj, node in enumerate(preserve_nodes):
nodetxt += '{:d} {:.3f} {:.3f} {:.3f} {:d}{}'.format(jj + 1,
node.x(),
node.y(),
node.z(),
node.marker(),
os.linesep)
with open(filename, 'w') as out:
out.write(nodetxt)