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.

  • 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
  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • idx, type list of [int, int]
    x and y index of the DEM grid to start the direction search from
  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • x, y, type list/array of coordinates
    coordinates to be translated horizontally
  • 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).

  • see DEM class description
loadTXT(demfile, centering, easting_shift, northing_shift)[source]

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

  • see DEM class description
show(cmap='terrain', cbar=True, ax=None, **kwargs)[source]

Show digital elevation model. Needs an update.

  • 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.

  • 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.

  • mesh_name, type str
    mesh name
  • m_dir, type str
    path to mesh directory

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.

  • 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.

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_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.

  • 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!
  • 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.

  • 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].
  • 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.

  • r, type float
    radius of the cylinder
  • h, type float
    height of the cylinder
  • 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!

  • r1, type float
    inner radius of cylinder.
  • r2, type float
    outer radius of cylinder.
  • 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, h_fact=1.3, v_fact=1.2, d=100.0, marker_positions=None, overwrite_markers=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.

  • 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!
  • 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.

  • paths, type list
    paths is a list of pointsets (arrays with shape (n, 3)) which will be connected from point to point.
  • 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().

  • 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.

  • points, type array with shape (n, 3)
    arrqay of points to be added
  • 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.

  • 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].
  • 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.

  • 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.

  • 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.

  • node_positions, type array of shape (n, 3)
    positions of all nodes of the interface mesh
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.

  • tx, type list
    list of pointsets (arrays with shape (n, 3)) which will be connected from point to point.
  • 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.

  • 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.

  • 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.

  • 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=[], 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.

  • 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', 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.

  • 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.

  • node_positions, type array of shape (n, 3)
    positions of all nodes of the interface mesh
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.

  • 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”.

  • 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()[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.

  • nodes, type list/array with shape (:, 3)
    list or array of 3D coordinates
  • 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(node, nodes)[source]

Find the closest node in a list of 2D/3D corrdinates.

  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • poly, type pyGIMLi mesh
    piecewise linear complex which should be exported to the .poly file
  • filename, type str
    name of output file
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.

  • 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.

  • 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.

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

Check if coordinates are located on a frame.

  • 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_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.

  • 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

  • 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.

  • a, b, c, type list/array of length 3
    three points to be evaluated
  • 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.

  • a, b, c, type list/array of length 3
    three points to be evaluated
  • 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.

  • a, b, c, type list/array of length 3
    three points to be evaluated
  • 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.

  • a, b, c, type list/array of length 3
    three points to be evaluated
  • 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.

  • a, b, c, type list/array of length 3
    three points to be evaluated
  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • 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.

  • path, type array with shape (:, 2) or (:, 3)
    array of coordinates describing the polygon
custEM.meshgen.meshgen_utils.refine_adaptive(coords, txs, r=10.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.

  • 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.

  • coords, type list/array with shape (:, 3)
    list or array of 3D coordinates
  • 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.

  • array, type array of shape (:, 3)
    array of coordinates to be rotated
  • alpha, type float
    angle(RAD) specifying the rotation
  • 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.

  • 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
  • pre_rotate = False, type bool
    conduct an initial roation using the second entry of axes/angles
custEM.meshgen.meshgen_utils.translate(array, shift, back=False)[source]

Translate array about given shift

  • 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
  • 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.

  • 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
  • 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.

  • 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.

  • mesh_name, type str
    mesh name
  • m_dir, type str
    path to mesh directory

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.

  • 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