Source code for meshpy.geometry

```__doc__ = """

Geometry builder
----------------

.. autoclass:: GeometryBuilder

Geometries
----------

These functions are designed so that their output can be splat-passed to

builder = GeometryBuilder()

.. autoclass:: Marker
:members:
:undoc-members:

.. autofunction:: make_box
.. autofunction:: make_circle
.. autofunction:: make_ball
.. autofunction:: make_cylinder

Extrusions and surfaces of revolution
-------------------------------------

.. data:: EXT_OPEN
.. data:: EXT_CLOSED_IN_RZ

.. autofunction:: generate_extrusion
.. autofunction:: generate_surface_of_revolution
"""

import numpy as np

# {{{ geometry building

def bounding_box(points):
return (
np.asarray(np.min(points, axis=0), dtype=np.float64),
np.asarray(np.max(points, axis=0), dtype=np.float64))

def is_multi_polygon(facets):
if not len(facets):
return False

try:
facets[0][0][0]  # facet 0, poly 0, point 0
except TypeError:
# pure python raises this
return False
except IndexError:
# numpy raises this
return False
else:
return True

def offset_point_indices(facets, offset):
if is_multi_polygon(facets):
return [[tuple(p_i+offset for p_i in poly)
for poly in facet]
for facet in facets]
else:
return [tuple(p_i+offset for p_i in facet) for facet in facets]

[docs]
class GeometryBuilder:
"""
.. automethod:: set
.. automethod:: wrap_in_box
.. automethod:: bounding_box
.. automethod:: center
.. automethod:: apply_transform
"""
def __init__(self):
self.points = []
self.facets = []
self.facet_hole_starts = None
self.facet_markers = None
self.point_markers = None

[docs]
facet_markers=None, point_markers=None):
if isinstance(facet_markers, int):
facet_markers = len(facets) * [facet_markers]

if facet_hole_starts and not self.facet_hole_starts:
self.facet_hole_starts = len(self.facets) * []
if facet_markers and not self.facet_markers:
self.facet_markers = len(self.facets) * [0]
if point_markers and not self.point_markers:
self.point_markers = len(self.points) * [0]

if not facet_hole_starts and self.facet_hole_starts:
facet_hole_starts = len(facets) * [[]]
if not facet_markers and self.facet_markers:
facet_markers = len(facets) * [0]
if not point_markers and self.point_markers:
point_markers = len(points) * [0]

if is_multi_polygon(facets) and not is_multi_polygon(self.facets):
self.facets = [[facet] for facet in self.facets]

if not is_multi_polygon(facets) and is_multi_polygon(self.facets):
facets = [[facet] for facet in facets]

self.facets.extend(offset_point_indices(facets, len(self.points)))
self.points.extend(points)

if facet_markers:
self.facet_markers.extend(facet_markers)
assert len(facets) == len(facet_markers)
if facet_hole_starts:
self.facet_hole_starts.extend(facet_hole_starts)
assert len(facets) == len(facet_hole_starts)
if point_markers:
self.point_markers.extend(point_markers)
assert len(points) == len(point_markers)

def make_facets():
end = len(points)-1
for i in range(end):
yield i, i+1
yield end, 0

facet_markers=facet_markers,
point_markers=point_markers)

def dimensions(self):
return len(self.points[0])

[docs]
def set(self, mesh_info):
"""Transfer the built geometry into a :class:`meshpy.triangle.MeshInfo`
or a :class:`meshpy.tet.MeshInfo`.
"""

mesh_info.set_points(self.points, self.point_markers)
if self.facet_hole_starts or is_multi_polygon(self.facets):
mesh_info.set_facets_ex(self.facets,
self.facet_hole_starts, self.facet_markers)
else:
mesh_info.set_facets(self.facets, self.facet_markers)

def mesher_module(self):
dim = self.dimensions()
if dim == 2:
import meshpy.triangle
return meshpy.triangle
elif dim == 3:
import meshpy.tet
return meshpy.tet
else:
raise ValueError("unsupported dimensionality %d" % dim)

[docs]
def bounding_box(self):
return bounding_box(self.points)

[docs]
def center(self):
a, b = bounding_box(self.points)
return (a+b)/2

[docs]
def wrap_in_box(self, distance, subdivisions=None):
"""
:param subdivisions: is a tuple of integers specifying
the number of subdivisions along each axis.
"""

a, b = bounding_box(self.points)
points, facets, _, facet_markers = \
make_box(a-distance, b+distance, subdivisions)

[docs]
def apply_transform(self, f):
self.points = [f(x) for x in self.points]

# }}}

# {{{ actual geometries

[docs]
class Marker:
MINUS_X = 1
PLUS_X = 2
MINUS_Y = 3
PLUS_Y = 4
MINUS_Z = 5
PLUS_Z = 6
SHELL = 100

FIRST_USER_MARKER = 1000

[docs]
def make_box(a, b, subdivisions=None):
"""
:param subdivisions: is a tuple of integers specifying
the number of subdivisions along each axis.
"""

a = [float(ai) for ai in a]
b = [float(bi) for bi in b]

assert len(a) == len(b)

dimensions = len(a)
if dimensions == 2:
# CAUTION: Do not change point or facet order here.
# Other code depends on this staying the way it is.

points = [
(a[0], a[1]),
(b[0], a[1]),
(b[0], b[1]),
(a[0], b[1]),
]

facets = [(0, 1), (1, 2), (2, 3), (3, 0)]

facet_markers = [
Marker.MINUS_Y, Marker.PLUS_X,
Marker.PLUS_Y, Marker.MINUS_X]

elif dimensions == 3:
#    7--------6
#   /|       /|
#  4--------5 |  z
#  | |      | |  ^
#  | 3------|-2  | y
#  |/       |/   |/
#  0--------1    +--->x

points = [
(a[0], a[1], a[2]),
(b[0], a[1], a[2]),
(b[0], b[1], a[2]),
(a[0], b[1], a[2]),
(a[0], a[1], b[2]),
(b[0], a[1], b[2]),
(b[0], b[1], b[2]),
(a[0], b[1], b[2]),
]

facets = [
(0, 1, 2, 3),
(0, 1, 5, 4),
(1, 2, 6, 5),
(7, 6, 2, 3),
(7, 3, 0, 4),
(4, 5, 6, 7)
]

facet_markers = [Marker.MINUS_Z, Marker.MINUS_Y, Marker.PLUS_X,
Marker.PLUS_Y, Marker.MINUS_X, Marker.PLUS_Z]
else:
raise ValueError("unsupported dimension count: %d" % len(a))

if subdivisions is not None:
if dimensions != 2:
raise NotImplementedError(
"subdivision not implemented for any "
"dimension count other than 2")

from meshpy.triangle import subdivide_facets
points, facets, facet_markers = subdivide_facets(
[subdivisions[0], subdivisions[1],
subdivisions[0], subdivisions[1]],
points, facets, facet_markers)

return points, facets, None, facet_markers

[docs]
def make_circle(r, center=(0, 0), subdivisions=40, marker=Marker.SHELL):
def round_trip_connect(seq):
result = []
for i in range(len(seq)):
result.append((i, (i+1) % len(seq)))
return result

phi = np.linspace(0, 2*np.pi, num=subdivisions, endpoint=False)
cx, cy = center
x = r*np.cos(phi) + cx
y = r*np.sin(phi) + cy

return ([np.array(pt) for pt in zip(x, y)],
round_trip_connect(list(range(subdivisions))),
None,
subdivisions*[marker])

[docs]
def make_ball(r, subdivisions=10):
from math import pi, cos, sin

dphi = pi/subdivisions

def truncate(my_r):
if abs(my_r) < 1e-9*r:
return 0
else:
return my_r

rz = [(truncate(r*sin(i*dphi)), r*cos(i*dphi)) for i in range(subdivisions+1)]

return generate_surface_of_revolution(

[docs]
height_subdivisions=1):
dz = height/height_subdivisions
rz = [(0, 0)] \
+ [(radius, i*dz) for i in range(height_subdivisions+1)] \
+ [(0, height)]
ring_markers = [Marker.MINUS_Z] \
+ ((height_subdivisions)*[Marker.SHELL]) \
+ [Marker.PLUS_Z]

return generate_surface_of_revolution(rz,
ring_markers=ring_markers)

# }}}

# {{{ extrusions

def _is_same_float(a, b, threshold=1e-10):
if abs(a) > abs(b):
a, b = b, a

# now abs(a) <= abs(b) always
return abs(b) < threshold or abs(a-b) < threshold*abs(b)

EXT_OPEN = 0
EXT_CLOSED_IN_RZ = 1

[docs]
def generate_extrusion(rz_points, base_shape, closure=EXT_OPEN,
point_idx_offset=0, ring_point_indices=None,
ring_markers=None, rz_closure_marker=0):
"""Extrude a given connected *base_shape* (a list of (x,y) points)
along the z axis. For each step in the extrusion, the base shape
is multiplied by a radius and shifted in the z direction. Radius
and z offset are given by *rz_points*, which is a list of
(r, z) tuples.

Returns ``(points, facets, facet_holestarts, markers)``, where *points* is
a list of (3D) points and facets is a list of polygons. Each polygon is, in
turn, represented by a tuple of indices into *points*. If
*point_idx_offset* is not zero, these indices start at that number.
*markers* is a list equal in length to *facets*, each specifying the facet
marker of that facet.  *facet_holestarts* is also equal in length to
*facets*, each element is a list of hole starting points for the
corresponding facet.

Use :meth:`~meshpy.tet.MeshInfo.set_facets_ex` to add the extrusion to a
:class:`~meshpy.tet.MeshInfo` structure.

The extrusion proceeds by generating quadrilaterals connecting each
ring.  If any given radius in *rz_points* is 0, triangle fans are

If *closure* is :data:`EXT_OPEN`, no efforts are made to put end caps on the
extrusion.

If *closure* is :data:`EXT_CLOSED_IN_RZ`, then a torus-like structure
is assumed and the last ring is just connected to the first.

If *ring_markers* is not None, it is an list of markers added to each
ring. There should be len(rz_points)-1 entries in this list.
corresponding *XXX_closure_marker*.  If *facet_markers* is given, this function
returns (points, facets, markers), where markers is is a list containing
a marker for each generated facet. Unspecified markers generally
default to 0.

If *ring_point_indices* is given, it must be a list of the same
length as *rz_points*. Each entry in the list may either be None,
or a list of point indices. This list must contain the same number
of points as the *base_shape*; it is taken as the indices of
pre-existing points that are to be used for the given ring, instead
of generating new points.
"""

assert len(rz_points) > 0

if ring_markers is not None:
assert len(rz_points) == len(ring_markers)+1

def get_ring(ring_idx):
try:
return rings[ring_idx]
except KeyError:
# need to generate fresh ring, continue
pass

p_indices = None
if ring_point_indices is not None:
p_indices = ring_point_indices[ring_idx]

first_idx = point_idx_offset+len(points)

r, z = rz_points[ring_idx]

if r == 0:
p_indices = (first_idx,)
points.append((0, 0, z))
else:
p_indices = tuple(range(first_idx, first_idx+len(base_shape)))
points.extend([(x*r, y*r, z) for (x, y) in base_shape])

rings[ring_idx] = p_indices
return p_indices

def pair_with_successor(ln):
n = len(ln)
return [(ln[i], ln[(i+1) % n]) for i in range(n)]

"""Add several new facets, each polygon in new_polys corresponding
to a new facet.
"""
facets.extend([poly] for poly in new_polys)
markers.extend(len(new_polys)*[marker])
holelists.extend(len(new_polys)*[[]])

"""Add a single facet, with each polygon in *facet_polygons*
belonging to a single facet.
"""
facets.append(facet_polygons)
markers.append(marker)
holelists.append(holestarts)

def connect_ring(ring1_idx, ring2_idx, marker):
r1, z1 = rz_points[ring1_idx]
r2, z2 = rz_points[ring2_idx]

if _is_same_float(z2, z1):
assert not _is_same_float(r1, r2)
# we're moving purely outward--this doesn't need fans, only plane
# surfaces. Special casing this leads to more freedom for TetGen
# and hence better meshes.

if r1 == 0:
# make opening surface
if r2 != 0:
elif r2 == 0:
# make closing surface
else:
# make single-surface interface with hole
get_ring(ring1_idx),
get_ring(ring2_idx),
],
holestarts=[(0, 0, z1)], marker=marker)
else:
ring1 = get_ring(ring1_idx)
ring2 = get_ring(ring2_idx)
if r1 == 0:
# make opening fan
assert len(ring1) == 1
start_pt = ring1[0]

if r2 != 0:
[(start_pt, succ, pt)
for pt, succ in pair_with_successor(ring2)],
marker=marker)
elif r2 == 0:
# make closing fan
assert len(ring2) == 1
end_pt = ring2[0]
[(pt, succ, end_pt)
for pt, succ in pair_with_successor(ring1)],
marker=marker)
else:
pairs1 = pair_with_successor(ring1)
pairs2 = pair_with_successor(ring2)
[(a, b, c, d) for ((a, b), (d, c)) in zip(pairs1, pairs2)],
marker=marker)

points = []
facets = []
markers = []
holelists = []

rings = {}

# pre-populate ring dict with ring_point_indices
if ring_point_indices is not None:
for i, ring_points in enumerate(ring_point_indices):
if ring_points is not None:
assert isinstance(ring_points, tuple)

if rz_points[i][0] == 0:
assert len(ring_points) == 1
else:
assert len(ring_points) == len(base_shape), \
("Ring points length (%d) does not "
"match base shape length (%d)"
% (len(ring_points), len(base_shape)))

rings[i] = ring_points

for i in range(len(rz_points)-1):
if ring_markers is not None:
ring_marker = ring_markers[i]
else:
ring_marker = 0

connect_ring(i, i+1, ring_marker)

if closure == EXT_CLOSED_IN_RZ:
connect_ring(len(rz_points)-1, 0, rz_closure_marker)

return points, facets, holelists, markers

[docs]
def generate_surface_of_revolution(rz_points,
point_idx_offset=0, ring_point_indices=None,
ring_markers=None, rz_closure_marker=0):
from math import sin, cos, pi