porepy.geometry.intersections module

Module with functions for computing intersections between geometric objects.

line_tessellation(p1, p2, l1, l2)[source]

Compute intersection of two line segment tessellations of a line.

The function will identify partly overlapping line segments between l1 and l2, and compute their common length.

Parameters

p1 (ndarray) –
shape=(3, n_p1)

Points in first tessellation.
p2 (ndarray) –
shape=(3, n_p2)

Points in second tessellation.
l1 (ndarray) –
shape=(2, n_tri_1)

Line segments in first tessellation, referring to indices in p2.
l2 (ndarray) –
shape= (2, n_tri_1)

Line segments in second tessellation, referring to indices in p2.

Raises

AssertionError – If porepy.intersections.segments_3d returns an unknown shape.

Returns

List of tuples with each representing an overlap. The tuple consists of

int:: the index of the segment in the first tesselation.
int:: the index of the segment in the second tesselation.
float:: the common length of the two segments in the two tessellations.

Return type

list[tuple[int, int, float]]

polygons_3d(polys, target_poly=None, tol=1e-8, include_point_contact=True)[source]

Compute the intersection between polygons embedded in 3d.

In addition to intersection points, the function also decides:

Whether intersection points lie in the interior, on a segment or a vertex. If segment or vertex, the index of the segment or vertex is returned.

Whether a pair of intersection points lie on the same boundary segment of a polygon, that is, if the polygon has a T or L-type intersection with another polygon.

Assumptions:

All polygons are convex. Non-convex polygons will simply be treated in a wrong way. To circumvent this, split the non-convex polygon into convex parts.

No polygon contains three points on a line, that is, an angle of pi. This can be included, possibly by temporarily stripping the hanging node from the polygon definition.

If two polygons meet in a vertex, this is not considered an intersection.

If two polygons lie in the same plane, intersection types (vertex, segment, interior) are not classified. This will be clear from the returned values. Inclusion of this should be possible, but it has not been a priority.

Contact between polygons in a single point may not be accurately calculated.

Parameters

polys (list[numpy.ndarray]) –
shape=(3, np)

Each list item represents a polygon, specified by its vertices as a numpy array. There should be at least three vertices in the polygon.
target_poly (Optional[Union[int, ndarray]]) –
default=None

Index in poly of the polygons that should be targeted for intersection findings. These will be compared with the whole set in poly.

If not provided, all polygons are compared with each other.
tol (float) –
default=1e-8

Geometric tolerance for the computations.
include_point_contact (bool) –
default=True

If True, point contacts will be considered an intersection. This is an experimental feature, use with care.

Returns

ndarray: shape=(3, np)

Intersection coordinates.

ndarray:: For each of the polygons, give the index of the intersection points, referring to the columns of the intersection coordinates.
ndarray:: For each polygon, a list telling whether each of the intersections is on the boundary of the polygon or not. For polygon i, the first element in this list tells whether the point formed by point-indices 0 and 1 in the previous return argument is on the boundary.
list:: Each list element is a 2-tuple with the indices of intersecting polygons.
ndarray:: For each polygon, for all intersection points (same order as the second return value), a 2-tuple, where the first value gives an index, the second is a Boolean, True if the intersection is on a segment, False if vertex. The index identifies the vertex, or the first vertex of the segment. If the intersection is in the interior of a polygon, the tuple is replaced by an empty list.
ndarray:: For each polygon, for all intersection points, True if this intersection is formed by a single point.

Return type

Returns a tuple consisting of

segments_2d(start_1, end_1, start_2, end_2, tol=1e-8)[source]

Check if two line segments, defined by their start- and endpoints, intersect.

The lines are assumed to be in 2D.

If the lines are (almost) parallel, i.e. if

\[\begin{split}\begin{vmatrix} start_1 & end_1 \\ start_2 & end_2 \end{vmatrix}\leq tol\times|start_1-end_1|\times|start_2-end_2|,\end{split}\]

a line segment is returned instead of an intersection point.

Todo

This function can be replaced by a call to segments_3d().

Example

>>> segments_2d([0, 0], [1, 1], [0, 1], [1, 0])
array([[ 0.5],
       [ 0.5]])

>>> segments_2d([0, 0], [1, 1], [0, 0], [2, 2])
array([[0., 1.],
       [0., 1.]])

>>> segments_2d([0, 0], [1, 0], [0, 1], [1, 1]) is None
True

Parameters

start_1 (ndarray) –
shape=(2,)

Coordinates of start point for first line.
end_1 (ndarray) –
shape=(2,)

Coordinates of end point for first line.
start_2 (ndarray) –
shape=(2,)

Coordinates of start point for first line.
end_2 (ndarray) –
shape=(2,)

Coordinates of end point for first line.
tol (float) –
default=1e-8

Tolerance for detecting parallel lines.

Raises

ValueError – If the start and endpoints of a line are the same.

Returns

Coordinates of intersection point, or the endpoints of the intersection segments if relevant. In the case of a segment, the first point (column) will be closest to start_1. Shape is (2, np), where np is 1 for a point intersection, or 2 for a segment intersection. If the lines do not intersect, None is returned.

Return type

Optional[ndarray]

segments_3d(start_1, end_1, start_2, end_2, tol=1e-8)[source]

Find intersection points (or segments) of two 3d lines.

Parameters

start_1 (ndarray) –
shape=(3,)

Coordinates of start point for first line.
end_1 (ndarray) –
shape=(3,)

Coordinates of end point for first line.
start_2 (ndarray) –
shape=(3,)

Coordinates of start point for first line.
end_2 (ndarray) –
shape=(3,)

Coordinates of end point for first line.
tol (float) –
default=1e-8

Tolerance for detecting parallel lines.

Returns

Coordinates of intersection points. Shape is (3, np), where np is 1 for a point intersection, or 2 for a segment intersection. If the lines do not intersect, None is returned.

Return type

Optional[ndarray]

segments_polygon(start, end, poly, tol=1e-5)[source]

Compute the internal intersection from line segments to a polygon in 3d. Intersections with the boundary of the polygon are not computed.

Note

It is required that all points lie in a plane. A sanity check will be performed.

Example

>>> import numpy as np
>>> import porepy as pp
>>> start = np.array([0.5, 0.5, -0.5])
>>> end = np.array([0.5, 0.5, 0.5])
>>> poly = np.array([[0, 1, 1, 0], [0, 0, 1, 1], [0, 0, 0, 0]])
>>> is_cross, pt = pp.intersections.segments_polygon(start, end, poly)
>>> print(is_cross)
>>> print(pt)

Parameters

start (ndarray) –
shape=(nd, num_segments)

One endpoint of segments.
end (ndarray) –
shape=(nd, num_segments)

Other endpoint of segments.
poly (ndarray) –
shape=(nd, num_vertices)

Vertices of polygon.
tol (float) –
default=1e-5

Tolerance for the geometric computations.

Returns

A tuple consisting of

ndarray: shape=(num_segments)

boolean array, identifying whether a segment has an intersection with the polygon (useful to filter the second return parameter).

ndarray: shape=(nd, num_segments)

float array containing the intersection points.

Return type

tuple[numpy.ndarray, numpy.ndarray]

segments_polyhedron(start, end, poly, tol=1e-5)[source]

Compute the intersection from line segments to the interior of a convex polyhedron. Intersections with the boundary of the polyhedron are not computed.

Note

There are four possibilities for each segment:

1. the segment is completely inside the polyhedron, meaning that its vertices are both inside the polyhedron; 2. the segment has only one vertex in the polyhedron; 3. the segment is completely outside the polyhedron; 4. the segment has in intersection but both vertices are outside the polyhedron.

Example

>>> import numpy as np
>>> import porepy as pp
>>> s = np.array([0.5, 0.5, 0.25])
>>> e = np.array([0.5, 0.5, 0.75])
>>> p = np.array(
>>>     [
>>>         [[0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 1.0]],
>>>         [[1.0, 1.0, 1.0, 1.0], [0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 1.0]],
>>>         [[0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 1.0]],
>>>         [[0.0, 1.0, 0.0, 1.0], [1.0, 1.0, 1.0, 1.0], [0.0, 0.0, 1.0, 1.0]],
>>>         [[0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 1.0], [0.0, 0.0, 0.0, 0.0]],
>>>         [[0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0]],
>>>     ]
>>> )
>>> pts, s_in, e_in, perc = pp.intersections.segments_polyhedron(s, e, p)

Parameters

start (ndarray) –
shape=(nd, num_segments)

One endpoint of segments.
end (ndarray) –
shape=(nd, num_segments)

Other endpoint of segments.
poly (ndarray) –
shape(nd, num_vertices)

Vertices of polyhedron organised face by face.
tol (float) –
default=1e-5

Tolerance for the geometric computations.

Returns

A tuple consisting of

ndarray: shape=(num_segments,)

Intersection points with the polyhedron, start and end points are not included in this list.

ndarray: shape=(num_segments,)

Boolean array indicating whether the start of a segment is inside the polyhedron.

ndarray: shape=(num_segments,)

Boolean array indicating whether the end of a segment is inside the polyhedron.

ndarray: shape=(num_segments,)

Length percentage of a segment inside the polyhedron.

Return type

tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]

split_intersecting_segments_2d(p, e, tol=1e-8, return_argsort=False)[source]

Process a set of points and connections between them so that the result is an extended point set and new connections that do not intersect.

The function is written for gridding of fractured domains, but may be of use in other cases as well. The geometry is assumed to be 2D.

The connections are defined by their start and endpoints, and can also have tags assigned. If so, the tags are preserved as connections are split. The connections are uniquified, so that no combination of point indices occurs more than once.

Note

For (partly) overlapping segments, only one of the tags will survive the uniquification. The other can be reconstructed by using the third output.

Parameters

p (ndarray) –
shape=(2, n_pt)

Coordinates of points to be processed.
e (ndarray) –
shape=(n, n_con)

Connections between lines. n >= 2, row 0 and 1 are index of start and endpoints, additional rows are tags.
tol (float) –
default=1e-8

Tolerance used for comparing equal points.
return_argsort (bool) –
default=False

Return the mapping between the input segments and the output segments.

Returns

A tuple with three or four entries (last only returned if return_argsort is True). The entries are as follows

ndarray: shape=(2, n_pt)

Points.

ndarray: shape=)2, n_edges)

new, non-intersecting edges.

2-tuple of ndarray:: two arrays with length n_con with the first item being a set of tags, before uniquification of the edges, and the second being a column mapping from the unique edges to all edges. To recover lost tags associated with the points in column i, first find all original columns which maps to i (tuple[1] == i), then recover the tags by the hits.

ndarray: shape=(n_edges,)

mapping of the new edges to the input edges.

Return type

Union[tuple[numpy.ndarray, numpy.ndarray, tuple[numpy.ndarray, numpy.ndarray], numpy.ndarray], tuple[numpy.ndarray, numpy.ndarray, tuple[numpy.ndarray, numpy.ndarray]]]

surface_tessellations(poly_sets, return_simplexes=False)[source]

Intersect a set of surface tessellations to find a finer subdivision that does not intersect with any of the input tessellations.

It is assumed that the polygon sets are 2d.

The implementation relies heavily on shapely’s intersection finders.

Parameters

poly_sets (list[list[numpy.ndarray]]) – Lists of polygons to be intersected.
return_simplexes (bool) –
default=False

If True, the subdivision is further split into a triangulation. The mappings from the original polygons is updated accordingly.

Raises

NotImplementedError – If a triangulation of a non-convex polygon is attempted. Can only happen if return_simplexes is True.

Returns

Tuple consisting of

list of ndarray:: Each element being a polygon so that the list together form a subdivision of the intersection of all polygons in the input sets.
list of spmatrix:: Mappings from each of the input polygons to the intersected polygons. If the mapping’s item[i][j, k] is non-zero, polygon k in set i has a (generally partial) overlap with polygon j in the intersected polygon set. Specifically the value will be 1.

Return type

tuple[list[numpy.ndarray], list[scipy.sparse._csr.csr_matrix]]

triangulations(p_1, p_2, t_1, t_2)[source]

Compute intersection of two triangle tessellations of a surface.

The function will identify partly overlapping triangles between t_1 and t_2, and compute their common area. If parts of domain 1 or 2 are covered by one tessellation only, this will simply be ignored by the function.

Note

The function relies on the intersection algorithm in shapely.geometry.Polygon.