custEM.meshgen package

Submodules

custEM.meshgen.bathy_tools module

@author: Rochlitz.R

class custEM.meshgen.bathy_tools.Bathy(bathy_file, x=None, y=None, centering=False, easting_shift=None, northing_shift=None, min_dist=50.0, min_points=10, coast_refinement=None, do_nothing=False)[source]

Bases: object

Coastline identification class for digital elevation models. Please note that the basic methods are working properly for a couple of specific models, but a more general application of this class needs some debugging and further extensions.

add_point_to_coast_line(idx, direction, clockwise, tol=None)[source]

Add a coordinate to a coastline path if the distance to the previous position on the coastline is larger than min_dist.

Required arguments

  • idx, type list of [int, int]

    x and y index of the DEM grid

  • direction, type str

    direction string, either ‘w’, ‘a’, ‘s’, or ‘d’

  • clockwise, type bool

    go in clock- or counterclockwise direction

Keyword arguments

  • tol = None, type float

    custom search tolerance if required, not recommended to be changed

add_point_towards_east(idx)[source]

Add a coordinate to a coastline path in eastern direction.

Required arguments

  • idx, type list of [int, int]

    x and y index of the DEM grid

add_point_towards_north(idx)[source]

Add a coordinate to a coastline path in northern direction.

Required arguments

  • idx, type list of [int, int]

    x and y index of the DEM grid

add_point_towards_south(idx)[source]

Add a coordinate to a coastline path in southern direction.

Required arguments

  • idx, type list of [int, int]

    x and y index of the DEM grid

add_point_towards_west(idx)[source]

Add a coordinate to a coastline path in western direction.

Required arguments

  • idx, type list of [int, int]

    x and y index of the DEM grid

build_closing_frame(idx, clockwise)[source]

Continue coastline path in clockwise or counterclockwise direction if the horizontal mesh boundary is reached.

Required arguments

  • idx, type list of [int, int]

    x and y index of the DEM grid

  • clockwise, type bool

    go in clock- or counterclockwise direction

check_distance(x, y)[source]

Calculate distance between a coordinate [x, y] and the previous position on the coastline.

Required arguments

  • x, y, typefloat

    x and y coordinates for the horizontal distance evaluation

check_neighbours(kk, jj)[source]

Check z-values of the eight neighbors around a position in the DEM. Utility function that needs an update.

Required arguments

  • kk, jj, type int

    x/y indices of the DEM grid

eval_coast_line()[source]

Main method to find intersection of topo and zero-height level.

find_direction(idx, from_boundary=False)[source]

Identifiy clock- or counterclockwise search path direction based on regular grid DEM data automatically. Does not work in all cases so far.

Required arguments

  • idx, type list of [int, int]

    x and y index of the DEM grid to start the direction search from

Keyword arguments

  • from_boundary = False, type bool

    internally used flag to check if a coastline path was following the horizontal mesh boundaries before entering the grid again

find_last_dir(pre_idx, idx)[source]

Evaluate previous direction during the coastline search.

Required arguments

  • pre_idx, type list of [int, int]

    previous x and y index of the DEM grid

  • idx, type list of [int, int]

    x and y index of the DEM grid

find_start_point(jjj=0)[source]

Evaluate start indices for a new coastline search in the DEM.

Keyword arguments

  • jjj, type int

    counter to continute until all indices in the DEM-grid were considered

go_clockwise(side, idx, from_corner=False)[source]

Continue coastline path in clockwise direction if the horizontal mesh boundary is reached.

Required arguments

  • side, type str

    string specifying one of the four sides

  • idx, type list of [int, int]

    x and y index of the DEM grid

  • from_corner = False, type bool

    internally used flag to continue a path after a corner of the horizontal mesh extent was reached

go_counter_clockwise(side, idx, from_corner=False)[source]

Continue coastline path in counterclockwise direction if the horizontal mesh boundary is reached.

Required arguments

  • side, type str

    string specifying one of the four sides

  • idx, type list of [int, int]

    x and y index of the DEM grid

  • from_corner = False, type bool

    internally used flag to continue a path after a corner of the horizontal mesh extent was reached

loadASC(topo_file, centering, easting_shift, northing_shift)[source]

Load asc file (DEM matrix with location header). This function is redudant to an identical function in the DEM class and will be removed in a major update of the bathy-tools.

Required arguments

  • see Bathy class description

loadTXT(topo_file, centering, easting_shift, northing_shift)[source]

Load txt file with column-based list of DEM coordinates. This function is redudant to an identical function in the DEM class and will be removed in a major update of the bathy-tools.

Required arguments

  • see Bathy class description

translate(x, y, centering=False, easting_shift=None, northing_shift=None, back=False)[source]

Function to translate coordinates in horizontal direction.

Required arguments

  • x, y, type list/array of coordinates

    coordinates to be translated horizontally

Keyword arguments

  • see Bathy class description

custEM.meshgen.dem_interpolator module

@authors: Rochlitz.R & Günther.T

class custEM.meshgen.dem_interpolator.DEM(demfile, x=None, y=None, centering=False, easting_shift=None, northing_shift=None, revert_z=False)[source]

Bases: object

Interpolation class for digital elevation models (DEM). Automatically creates interpolation objects during the class initialization, based on DEM information stored as list of coordinates in txt format or as matrix in asc format.

loadASC(ascfile, centering, easting_shift, northing_shift)[source]

Load ASC file (DEM matrix with location header).

Required arguments

  • see DEM class description

loadTXT(demfile, centering, easting_shift, northing_shift)[source]

Load txt file with column-based list of DEM coordinates.

Required arguments

  • see DEM class description

show(cmap='terrain', cbar=True, ax=None, **kwargs)[source]

Show digital elevation model. Needs an update.

Keyword arguments

  • cmap = “terrain”, type str ()

    matplotlib colormap definiton

  • cbar = True, type bool

    add colorbar to the plot or not

  • ax = None, type matplotlib figure axes object

    add the plot to a given axes object or create a new one

  • **kwargs, type keyword arguments

    add additional keyword arguments for the plot style (e.g., lw)

translate(centering=False, easting_shift=None, northing_shift=None, back=False)[source]

Function to translate coordinates in horizontal direction.

Keyword arguments

  • see DEM class description

  • back = False, type bool

    set True to translate back in the opposite direction; useful for rotating coordinates around a specific origin

custEM.meshgen.invmesh_tools module

@author: Rochlitz.R

class custEM.meshgen.invmesh_tools.PrismWorld(name, x_extent, x_reduction=100.0, y_depth=300.0, z_depth=300.0, prism_quality=32.0, prism_area=10000.0, n_prisms=50.0, orthogonal_txs=None, txs=[], surface_rxs=[], **blankworld_kwargs)[source]

Bases: object

Specific class for automated generation of pseudo-2D / 3D inversion meshes discretized by prisms.

add_prisms(back_nodes, front_nodes)[source]

Add nodes and triangle faces to include horizontal prisms as conformal PLC.

add_triangles(back_nodes, front_nodes)[source]

Add front and back triangle faces of the prisms. All the if statements are required to account for potential nodes on the x-directed boundary edges of the rectangular prism outcrops at the surface.

create_2d_mesh(name, prism_quality, prism_area)[source]

Create and store 2D mesh in x-z slice with triangle, which will be the basis for elongating the triangles to prisms in y-direction later on.

create_surface_rectangles()[source]

Build rectangular frames of the prisms at the surface to be included in the surface mesh later on for avoiding intersections. Check also if transmitter paths are intersecting these rectangles, and if yes, find the intersections and add the corresponding points to the frames.

sort_surface_node_ids(ids)[source]

Sort nodes at the surface for conformal triangle insertion.

custEM.meshgen.mesh_convert module

@author: modified by Rochlitz.R

custEM.meshgen.mesh_convert.mesh2xml2h5(mesh_name, m_dir, rot)[source]

Converts TetGen meshes in Medit (.mesh) format to xml and h5 files for the usage in FEniCS. This function is called automatically during the intialization of the MOD class and converts a mesh if it does not already exist as xml and h5 file.

Required arguments

  • mesh_name, type str

    mesh name

  • m_dir, type str

    path to mesh directory

Functionality

Convert between .mesh and .xml, parser implemented as a state machine:

0 = read ‘Dimension’ 1 = read dimension 2 = read ‘Vertices’ 3 = read number of vertices 4 = read next vertex 5 = read ‘Triangles’ or ‘Tetrahedra’ 6 = read number of cells 7 = read next cell 8 = read next domains 9 = done

custEM.meshgen.mesh_convert.write_cell_interval(ofile, cell, n0, n1)[source]

Write interval (1D) information to domain file.

custEM.meshgen.mesh_convert.write_cell_tetrahedron(ofile, cell, n0, n1, n2, n3)[source]

Write tetrahedron 3D information to file.

custEM.meshgen.mesh_convert.write_cell_triangle(ofile, cell, n0, n1, n2)[source]

Write triangle 2D information to domain file.

custEM.meshgen.mesh_convert.write_domain_marker(ofile, cell, marker, dfile=False)[source]

Write domain markers to domain file.

Write footer cells to file.

Write footer to domain file.

Write footer edges to file.

Write footer graph to file.

Write footer of mesh file.

Write footer vertices to file.

custEM.meshgen.mesh_convert.write_graph_edge(ofile, v1, v2, weight=1)[source]

Write edge graph to file.

custEM.meshgen.mesh_convert.write_graph_vertex(ofile, vertex, num_edges, weight=1)[source]

Write vertex graph to file.

custEM.meshgen.mesh_convert.write_header_cells(ofile, num_cells)[source]

Write header cells to file.

custEM.meshgen.mesh_convert.write_header_domains(ofile, num_domains, dfile=False)[source]

Write header to domain file.

custEM.meshgen.mesh_convert.write_header_edges(ofile, num_edges)[source]

Write header edges to file.

custEM.meshgen.mesh_convert.write_header_graph(ofile, graph_type)[source]

Write header graph to file.

custEM.meshgen.mesh_convert.write_header_mesh(ofile, cell_type, dim)[source]

Write header of mesh file.

custEM.meshgen.mesh_convert.write_header_vertices(ofile, num_vertices)[source]

Write header vertices to file.

custEM.meshgen.mesh_convert.write_vertex(ofile, vertex, x, y, z)[source]

Write vertex to file.

custEM.meshgen.mesh_convert.write_xml_domain_file(ofile_name, num_domains, marker)[source]

Write domain marker to mesh function file.

custEM.meshgen.mesh_convert.xml_to_hdf5(xml, h5)[source]

Convert xml to HDF5 (h5) file.

Required arguments

  • xml, type str

    name of xml file

  • h5, type str

    name of h5 file

custEM.meshgen.meshgen_tools module

@author: Rochlitz.R

class custEM.meshgen.meshgen_tools.BlankWorld(**raw_kwargs)[source]

Bases: object

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

add_air_tri(paths)[source]

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.

add_area_marker(poly, pos, face_size=0.0, marker=777)[source]

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.

add_brick(start, stop, cell_size=0.0, marker=None)[source]

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.

add_cylinder(r, h, origin=[0.0, 0.0, -200.0], n_segs=16, dip=0.0, dip_azimuth=0.0, cell_size=0.0, marker=None, marker_position=None)[source]

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

add_hollow_cylinder(r1, r2, h=100.0, origin=[0.0, 0.0, 0.0], rotation=None, cell_size=0.0, n_segs=16, preserve=[], close=False, mark=False, close1st=False)[source]

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.

add_interface(frame=None, topo=None, id_offset=None, eval_zl=False)[source]

Add more or less horizontal 2D interface (layer interfaces) with or without topography to Omega.

add_intersecting_anomaly(surface_intersection=False, surface_path=[], intersection_paths=[], intersecting_layers=None, cell_size=None, bottom=None, top=None, marker=None, marker_position=None)[source]
add_inv_domains(depth, poly=None, rxs=None, cell_size=0.0, x_frame=1000.0, y_frame=1000.0, z_frame=1000.0, d=100.0, marker_positions=None, overwrite_markers=None, split_depths=None)[source]
add_marker(label, pos, marker=None, cell_size=0.0)[source]

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.

add_paths(paths, closed_paths=None, marker=None)[source]

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.

add_plate(dx=200.0, dy=200.0, dz=50.0, origin=[0.0, 0.0, -200.0], dip=0.0, dip_azimuth=0.0, cell_size=0.0, marker=None)[source]

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

add_points(points, marker=None)[source]

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

add_prism(poly, top, bottom, cell_size=0.0, marker=None, marker_pos=None)[source]

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.

add_quadrangle_face(n1, n2, n3, n4)[source]

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.

add_refinement_area(frame, quality, area, z=0.0, topo=False)[source]
add_rx(points, marker=None)[source]

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

add_surface_anomaly(insert_paths=[], depths=None, cell_sizes=None, marker=None, dips=None, dip_azimuths=None, split_depths=None, marker_positions=None, split_coords=None)[source]
add_topo(node_positions, topo=None, eval_zl=False)[source]

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)

add_tx(tx, grounded=True, marker=None)[source]

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.

build_fullspace_mesh(cell_size=None, boundary_marker=99)[source]

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

build_halfspace_mesh(boundary_marker=99, **mesh_build_kwargs)[source]

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

build_layered_earth_mesh(n_layers, layer_depths, insert_struct=None, closed_struct=None, tx_struct=None, struct_interface=0, boundary_marker=99, **le_kwargs)[source]

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

build_surface(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)[source]

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

call_tetgen(tet_param='default', ignore_tx_z=False, export_before=True, export_format=' -k', name=None, copy=True, print_infos=True, suppress='NEF')[source]

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.

create_box()[source]

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.

eval_unique_markers()[source]

extract unique markers specified in world Omega

extend_anomaly_outcrops()[source]

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!

extend_intersecting_anomaly()[source]

Function for extending anomaly-bodies which interesect layers.

extend_world(x_fact=10.0, y_fact=10.0, z_fact=10.0, marker=[0, 1], cell_sizes=None, boundary_marker=99, initial=True)[source]

alias for there_is_always_a_bigger_world() method

get_topo_vals(node_positions, topo=None, z=0.0)[source]

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)

init_dem(surface_mesh_coords, topo=None)[source]

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.

print_region_info()[source]

Print the domain information of the final mesh after calling TetGen.

remove_duplicate_poly_faces()[source]

Check which self.Omega.n_poly faces were added multiple times if subsurface structures touch each other and remove the duplicates.

there_is_always_a_bigger_world(x_fact=10.0, y_fact=10.0, z_fact=10.0, marker=[0, 1], cell_sizes=None, boundary_marker=99, initial=True)[source]

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.

update_kwargs(kwargs)[source]
write_mesh_parameters(ignore_tx_z)[source]

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.

custEM.meshgen.meshgen_utils module

@author: Rochlitz.R

custEM.meshgen.meshgen_utils.add_edges(surface, edges, nn, size, marker, closed=False)[source]

Add edge to polygon.

custEM.meshgen.meshgen_utils.assign_topography(nodes, t_dir=None, topo=None, z=0.0, centering=True, easting_shift=None, northing_shift=None, revert_z=False, rotation=None)[source]

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

custEM.meshgen.meshgen_utils.closest_node(coord, coords)[source]

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

custEM.meshgen.meshgen_utils.create_circle(r, n, h)[source]

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

custEM.meshgen.meshgen_utils.create_ellipse(a, b, n, h)[source]

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

custEM.meshgen.meshgen_utils.export_tetgen_edge_file(preserve_edges, edge_marker, filename)[source]

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

custEM.meshgen.meshgen_utils.export_tetgen_node_file(preserve_nodes, filename)[source]

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

custEM.meshgen.meshgen_utils.export_tetgen_poly_file(poly, filename, float_format='.3f', poly_marker=None, **kwargs)[source]

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

custEM.meshgen.meshgen_utils.find_closest(l1, l2)[source]

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

custEM.meshgen.meshgen_utils.find_edge_node_ids(frame, mesh, id_offset)[source]

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)

custEM.meshgen.meshgen_utils.find_nodes_on_segment(p0, p1, mesh, topo_f)[source]

Find nodes on a segment.

Required arguments

  • mesh

custEM.meshgen.meshgen_utils.find_on_frame(frame_coords, coords)[source]

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

custEM.meshgen.meshgen_utils.get_closest(coord, coords, n=None, max_dist=1.0)[source]

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

custEM.meshgen.meshgen_utils.get_intersect(a1, a2, b1, b2)[source]

Return intersection coordinate of two segments a and b.

custEM.meshgen.meshgen_utils.get_unique_poly_ids(polys, max_poly_length)[source]

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

custEM.meshgen.meshgen_utils.inside_poly(x, y, path)[source]

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

custEM.meshgen.meshgen_utils.is_between(a, b, c, a_tol=0.01, print_points=False)[source]

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

custEM.meshgen.meshgen_utils.is_between_1D(a, b, c, a_tol=0.01, print_points=False)[source]

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

custEM.meshgen.meshgen_utils.is_between_2D(a, b, c, a_tol=0.01, print_points=False)[source]

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

custEM.meshgen.meshgen_utils.is_between_2D_exact(a, b, c, a_tol=0.01, print_points=False)[source]

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

custEM.meshgen.meshgen_utils.is_between_2Dvert(a, b, c, a_tol=0.01, print_points=False)[source]

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

custEM.meshgen.meshgen_utils.line(start=[0.0, 0.0, 0.0], stop=[0.0, 0.0, 0.0], n_segs=100, z=0.0, topo=None, t_dir=None, **topo_kwargs)[source]

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

custEM.meshgen.meshgen_utils.line_x(start=-100.0, stop=100.0, n_segs=100, y=0.0, z=0.0, topo=None, t_dir=None, **topo_kwargs)[source]

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

custEM.meshgen.meshgen_utils.line_y(start=-100.0, stop=100.0, n_segs=100, x=0.0, z=0.0, topo=None, t_dir=None, **topo_kwargs)[source]

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

custEM.meshgen.meshgen_utils.line_z(start=-100.0, stop=100.0, n_segs=100, x=0.0, y=0.0)[source]

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

custEM.meshgen.meshgen_utils.loop_c(origin=[0.0, 0.0, 0.0], r=50.0, n_segs=100, z=0.0, topo=None, t_dir=None, **topo_kwargs)[source]

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

custEM.meshgen.meshgen_utils.loop_e(origin=[0.0, 0.0, 0.0], a=50.0, b=50.0, n_segs=100, z=0.0, topo=None, t_dir=None, **topo_kwargs)[source]

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

custEM.meshgen.meshgen_utils.loop_r(start=[-100.0, -100.0], stop=[100.0, 100.0], n_segs=4, nx=None, ny=None, z=0.0, topo=None, t_dir=None, loop_rotation=None, **topo_kwargs)[source]

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

custEM.meshgen.meshgen_utils.npoly_to_triangles(world, ids)[source]

Add polygones with more than three edges in form of triangles to Omega.

custEM.meshgen.meshgen_utils.order_list(ref, points, tmp)[source]

Sort values in list in increasing order. Deprecated and not used anymore.

custEM.meshgen.meshgen_utils.pointset(pointset)[source]

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)

custEM.meshgen.meshgen_utils.poly_area(path)[source]

Caluclate enclosed area of polygon.

Required arguments

  • path, type array with shape (:, 2) or (:, 3)

    array of coordinates describing the polygon

custEM.meshgen.meshgen_utils.poly_peri(path)[source]

Calculate perimeter of polygon.

Required arguments

  • path, type array with shape (:, 2) or (:, 3)

    array of coordinates describing the polygon

custEM.meshgen.meshgen_utils.refine_adaptive(coords, txs, r1=20.0, r2=5.0, r3=1.0, d3=200.0, d2=500.0, min_tx_dist=100.0, rot=30.0, n_segs=3)[source]

Decription to add

custEM.meshgen.meshgen_utils.refine_path(path, n_segs=None, length=None)[source]

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

custEM.meshgen.meshgen_utils.refine_rx(coords, r=5.0, rot=None, n_segs=3)[source]

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

custEM.meshgen.meshgen_utils.resolve_rx_overlaps(rxs, refinement_size=10.0, ignore_z=True)[source]
custEM.meshgen.meshgen_utils.rotate(array, alpha, axis='z', tensor=False)[source]

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

custEM.meshgen.meshgen_utils.rotate_around_point(array, shift, angles, axes, pre_rotate=False)[source]

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

custEM.meshgen.meshgen_utils.rx_grid(xlinsp, ylinsp, z=0.0)[source]

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.

custEM.meshgen.meshgen_utils.translate(array, shift, back=False)[source]

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

custEM.meshgen.meshgen_utils.write_synth_topo_to_asc(t_dir, file_name, topo, spacing=25.0, x_min=-10000.0, y_min=-10000.0)[source]

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

custEM.meshgen.meshgen_utils.write_synth_topo_to_xyz(t_dir, file_name, topo)[source]

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

custEM.meshgen.vtk_convert module

@author: modified by Rochlitz.R

custEM.meshgen.vtk_convert.vtk2xml2h5(mesh_name, m_dir, rot, overwrite_markers)[source]

Converts TetGen meshes in VTK (.vtk) format to xml and h5 files for the usage in FEniCS. This function is called automatically during the intialization of the MOD class and converts a mesh if it does not already exist as xml and h5 file.

Required arguments

  • mesh_name, type str

    mesh name

  • m_dir, type str

    path to mesh directory

Functionality

Convert between .vtk and .xml, parser implemented as a state machine:

0 = read number of vertices 1 = read vertices 2 = read number of cells 3 = read cells 4 = read number of markers 5 = read cell markers 6 = done

custEM.meshgen.vtk_convert.write_cell_tetrahedron(ofile, cell, n0, n1, n2, n3)[source]

Write tetrahedron 3D information to file.

custEM.meshgen.vtk_convert.write_domain_marker(ofile, cell, marker, dfile=False)[source]

Write domain markers to domain file.

Write footer cells to file.

Write footer to domain file.

Write footer of mesh file.

Write footer vertices to file.

custEM.meshgen.vtk_convert.write_header_cells(ofile, num_cells)[source]

Write header cells to file.

custEM.meshgen.vtk_convert.write_header_domains(ofile, num_domains, dfile=False)[source]

Write header to domain file.

custEM.meshgen.vtk_convert.write_header_mesh(ofile, cell_type, dim)[source]

Write header of mesh file.

custEM.meshgen.vtk_convert.write_header_vertices(ofile, num_vertices)[source]

Write header vertices to file.

custEM.meshgen.vtk_convert.write_vertex(ofile, vertex, x, y, z)[source]

Write vertex to file.

custEM.meshgen.vtk_convert.write_xml_domain_file(ofile_name, num_domains, markers)[source]

Write domain marker to mesh function file.

custEM.meshgen.vtk_convert.xml_to_hdf5(xml, h5)[source]

Convert xml to HDF5 (h5) file.

Required arguments

  • xml, type str

    name of xml file

  • h5, type str

    name of h5 file

Module contents

meshgen

Submodules:

  • meshgen_tools for tetrahedral mesh generation

  • meshgen_utils providing utility functions for mesh generation

  • invmesh_tools providing specific classes for building inversion meshes

  • mesh_convert for automated conversion of .mesh files to .xml and .h5

  • dem_interpolator for interpolating real digital elevation data

  • bathy_tools for identifying coastline paths