# Source code for pygimli.meshtools.polytools

# -*- coding: utf-8 -*-
"""Tools to create or manage PLC.

Please note there is currently no collision or intersection check at all.

Volunteers welcome to help creating, adapting or interfacing a basic
geometry system. A lot of things are needed:

* 2D
* 3D
* More geometric primitives
* Boolean operations (union, intersection, difference)
* Collision recognizing
* Cubic spline interpolation for polygons (partly done)
* GUI .. interactive creation
*
"""

import math
import os
from os import system

import numpy as np
import pygimli as pg

def _polyCreateDefaultEdges(poly, boundaryMarker=1, isClosed=True, **kwargs):
"""INTERNAL."""
nEdges = poly.nodeCount() - 1 + isClosed
bm = None
if hasattr(boundaryMarker, '__len__'):
if len(boundaryMarker) == nEdges:
bm = boundaryMarker
else:
raise Exception("marker length != nEdges", len(boundaryMarker),
nEdges)
else:
bm = [boundaryMarker] * nEdges

for i in range(poly.nodeCount() - 1):
poly.createEdge(poly.node(i), poly.node(i + 1), bm[i])

if isClosed:
poly.createEdge(poly.node(poly.nodeCount() - 1), poly.node(0), bm[-1])

def setPolyRegionMarker(poly, marker=1, area=0.0, markerPosition=None,
isHole=False, **kwargs):
"""Internal to set region markers to single elementary geometry.

Internal to set region markers.
If no absolute marker position is given.
The region marker is placed 1mm beside the first node in direction to the
geometry center.

Parameters
----------
poly : :gimliapi:GIMLI::Mesh
The Polygon that will get a marker
marker : int[1]
The region marker, every resulting mesh cell will get this marker.
area : float[0]
The region max cell size, every resulting mesh cell will get a cell
size lower than area in m² or m³ for 3D, respectively.
markerPosition : pg.Pos
Absolute marker position if you don't want the marker in the center of
the geometry.
isHole : bool [False]
Marks the geometry as a hole and will be cut in any merge mesh.

Keyword Arguments
----------------
**kwargs
"""
pos = None
if markerPosition is not None:
pos = markerPosition
else:
center = pg.center(poly.positions())
p0 = poly.node(0).pos()
# region marker near node 0; 1mm in direction to the center
# should be safer than the center itself
pos = p0 + (center-p0).norm() * 0.001

if isHole is True:
else:

[docs]def createRectangle(start=None, end=None, pos=None, size=None, **kwargs):
"""Create rectangle polygon.

Create rectangle with start position and a given size.
Give either start and end OR pos and size.

Parameters
----------
start : [x, y]
Left upper corner. Default [-0.5, 0.5]
end : [x, y]
Right lower corner. Default [0.5, -0.5]
pos : [x, y]
Center position. The rectangle will be moved.
size : [x, y]
Factors for x and y by which the rectangle, defined by **start** and
**width**, are scaled.

Keyword Arguments
-----------------
**kwargs

marker : int [1]
Marker for the resulting triangle cells after mesh generation
markerPosition : floats [x, y] [pos + (end - start) * 0.2]
Absolute position of the marker (works for both regions and
holes).
area : float [0]
Maximum cell size for resulting triangles after mesh generation
isHole : bool [False]
The polygon will become a hole instead of a triangulation
boundaryMarker : int [1]
Marker for the resulting boundary edges
leftDirection : bool [True]
TODO Rotational direction

Returns
-------
poly : :gimliapi:GIMLI::Mesh
The resulting polygon is a :gimliapi:GIMLI::Mesh.

Examples
--------
>>> # no need to import matplotlib, pygimli show does.
>>> import pygimli as pg
>>> import pygimli.meshtools as mt
>>> r1 = mt.createRectangle(pos=[1, -1], size=[4.0, 4.0],
...                      marker=1, area=0.1, markerPosition=[0, -2])
>>> r2 = mt.createRectangle(start=[0.5, -0.5], end=[2, -2],
...                      marker=2, area=1.1)
>>> _ = pg.show(mt.mergePLC([r1, r2]))
>>> pg.wait()
"""
if start is None:
start = [-0.5, 0.5]
if end is None:
end = [0.5, -0.5]

poly = pg.Mesh(dim=2, isGeometry=True)

sPos = pg.Pos(start)
ePos = pg.Pos(end)

verts = [sPos, [sPos[0], ePos[1]], ePos, [ePos[0], sPos[1]]]

# TODO refactor with polyCreatePolygon

if kwargs.pop("leftDirection", False):
for v in verts[::-1]:
poly.createNode(v)
else:
for v in verts:
poly.createNode(v)

# Note that we do not support the usage of start/end AND size/pos. Only one
# of the pairs. Otherwise strange things will happen with the region
# markers!
if size is not None:
poly.scale(size)
if pos is not None:
poly.translate(pos)

_polyCreateDefaultEdges(poly, **kwargs)

kwargs['markerPosition'] = kwargs.pop('markerPosition',
sPos + (ePos - sPos) * 0.2)

setPolyRegionMarker(poly, **kwargs)

return poly

[docs]def createWorld(start, end, marker=1, area=0., layers=None, worldMarker=True,
**kwargs):
"""Create simple rectangular world.

Create simple rectangular [hexagonal] world with appropriate boundary conditions.
Surface boundary is set do pg.core.MARKER_BOUND_HOMOGEN_NEUMANN, i.e, -1
and inner subsurface is set to pg.core.MARKER_BOUND_MIXED, i.e., -2 or
Numbered in ascending order in left direction starting upper left if
worldMarker is set to false.

TODO
----
* 3D with layers

Parameters
----------
start: [x, y, [z]]
Upper/Left/[Front] Corner
end: [x, y, [z]]
Lower/Right/[Back] Corner
marker: int
Marker for the resulting triangle cells after mesh generation.
area: float | list
Maximum cell size for resulting triangles after mesh generation.
If area is a float set it global, if area is a list set it per layer.
layers: [float] [None]
List of depth coordinates for some layers.
worldMarker: bool [True]
Specify kind of preset boundary marker [-1, -2] or ascending order [1, 2, 3, 4 ..]

Other Parameters
----------------
Forwarded to createCube

Returns
-------
poly : :gimliapi:GIMLI::Mesh
The resulting polygon is a :gimliapi:GIMLI::Mesh.

Examples
--------
>>> from pygimli.meshtools import createWorld
>>> from pygimli.viewer.mpl import drawMesh
>>> import matplotlib.pyplot as plt
>>> world = createWorld(start=[-5, 0], end=[5, -5], layers=[-1,-2,-3])
>>>
>>> fig, ax = plt.subplots()
>>> drawMesh(ax, world)
>>> plt.show()
"""
if len(start) == 3 and len(end) == 3:

if layers is not None:
pg.critical("3D with layers is not yet implemented.")

world = createCube(size=pg.Pos(end)-pg.Pos(start),
pos=(pg.Pos(end)-pg.Pos(start))/2.0,
**kwargs)

for i, b in enumerate(world.boundaries()):
if worldMarker is True:
if b.norm()[2] == 1.0:
b.setMarker(pg.core.MARKER_BOUND_HOMOGEN_NEUMANN)
else:
b.setMarker(pg.core.MARKER_BOUND_MIXED)
else:
if b.norm() == [-1, 0, 0]:
b.setMarker(1)
elif b.norm() == [1, 0, 0]:
b.setMarker(2)
elif b.norm() == [0, 0, 1]:
b.setMarker(3)
elif b.norm() == [0, 0, -1]:
b.setMarker(4)
elif b.norm() == [0, -1, 0]:
b.setMarker(5)
elif b.norm() == [0, 1, 0]:
b.setMarker(6)

return world

z = np.array(start[1], dtype=float)
if layers is not None:
z = np.append(z, np.array(layers, dtype=float))
z = np.append(z, end[1])
# ensure - decreasing order if layers are out of bounding box
z = np.sort(z)[::-1]

poly = pg.Mesh(dim=2, isGeometry=True)

if isinstance(area, float) or isinstance(area, int):
area = np.ones(len(z)-1) * float(area)

if len(area) < len(z) -1:
pg.warn('Missing {} area value, padding with zeros'.format(
(len(z) -1) - len(area)))
_area = np.zeros(len(z)-1)
_area[0:len(area)] = area
area = _area

for i, depth in enumerate(z):
n = poly.createNode([start[0], depth])
if i > 0:
if len(z) == 2:
marker=marker, area=area[0])
else:
marker=i, area=area[i - 1])

for i, depth in enumerate(z[::-1]):
poly.createNode([end[0], depth])

_polyCreateDefaultEdges(poly,
boundaryMarker=range(1, poly.nodeCount() + 1))

if worldMarker:
for b in poly.boundaries():
if b.norm()[1] == 1.0:
b.setMarker(pg.core.MARKER_BOUND_HOMOGEN_NEUMANN)
else:
b.setMarker(pg.core.MARKER_BOUND_MIXED)

if layers is not None:
for i in range(len(layers)):
poly.createEdge(poly.node(i + 1),
poly.node(poly.nodeCount() - i - 2),
poly.boundaryCount() + 1)

# pg.warnNonEmptyArgs(kwargs)
return poly

[docs]def createCircle(pos=None, radius=1, segments=12, start=0, end=2.*math.pi,
**kwargs):
"""Create simple circle polygon.

Create simple circle polygon with given attributes.

Parameters
----------
pos : [x, y] [[0.0, 0.0]]
Center position
radius : float | [a,b] [1]
radius or halfaxes of the circle
segments : int
Discrete amount of segments for the circle
start : double [0]
end : double [2*pi]

**kwargs:

marker : int [1]
Marker for the resulting triangle cells after mesh generation
markerPosition : floats [x, y] [0.0, 0.0]
Position of the marker (works for both regions and holes)
area : float [0]
Maximum cell size for resulting triangles after mesh generation
isHole : bool [False]
The polygon will become a hole instead of a triangulation
boundaryMarker : int [1]
Marker for the resulting boundary edges
leftDirection : bool [True]
Rotational direction
isClosed : bool [True]
Add closing edge between last and first node.

Returns
-------
poly : :gimliapi:GIMLI::Mesh
The resulting polygon is a :gimliapi:GIMLI::Mesh.

Examples
--------
>>>  # no need to import matplotlib. pygimli's show does
>>> import math
>>> from pygimli.viewer.mpl import drawMesh
>>> import pygimli.meshtools as mt
>>> c0 = mt.createCircle(pos=(-5.0, 0.0), radius=2, segments=6)
>>> c1 = mt.createCircle(pos=(-2.0, 2.0), radius=1, area=0.01, marker=2)
>>> c2 = mt.createCircle(pos=(0.0, 0.0), segments=5, start=0, end=math.pi)
>>> c3 = mt.createCircle(pos=(5.0, 0.0), segments=3, start=math.pi,
...                      end=1.5*math.pi, isClosed=False)
>>> plc = mt.mergePLC([c0, c1, c2, c3])
>>> fig, ax = pg.plt.subplots()
>>> drawMesh(ax, plc, fillRegion=False)
>>> pg.wait()
"""
# TODO refactor with polyCreatePolygon
if pos is None:
pos = [0.0, 0.0]

poly = pg.Mesh(dim=2, isGeometry=True)

dPhi = (end - start) / (segments)
nPhi = segments + 1

if abs((end % (2. * math.pi) - start)) < 1e-6:
nPhi = segments

for i in range(0, nPhi):
if kwargs.pop('leftDirection', True):
phi = start + i * dPhi
else:
phi = start - i * dPhi

xp = np.cos(phi)
yp = np.sin(phi)
poly.createNode([xp, yp])

else:
poly.translate(pos)

_polyCreateDefaultEdges(poly, **kwargs)

if kwargs.pop('isClosed', True):
setPolyRegionMarker(poly, **kwargs)

# need a better way mess with these or wrong kwargs
# pg.warnNonEmptyArgs(kwargs)

return poly

[docs]def createLine(start, end, segments=1, **kwargs):
"""Create simple line polygon.

Create simple line polygon from start to end.

Parameters
----------
start : [x, y]
start position
end : [x, y]
end position
segments : int
Discrete amount of segments for the line

Keyword Arguments
-----------------
boundaryMarker : int [1]
Marker for the resulting boundary edges
leftDirection : bool [True]
Rotational direction

Returns
-------
poly : :gimliapi:GIMLI::Mesh
The resulting polygon is a :gimliapi:GIMLI::Mesh.

Examples
--------
>>>  # no need to import matplotlib. pygimli's show does
>>> import pygimli as pg
>>> import pygimli.meshtools as mt
>>>
>>> w = mt.createWorld(start=[0, 0], end=[3, 3])
>>> l1 = mt.createLine(start=[1, 1], end=[1, 2], segments=1,
...                    leftDirection=False)
>>> l2 = mt.createLine(start=[1, 1], end=[2, 1], segments=20,
...                    leftDirection=True)
>>>
>>> ax, _ = pg.show(mt.createMesh([w, l1, l2,]))
>>> ax, _ = pg.show([w, l1, l2,], ax=ax, fillRegion=False)
>>> pg.wait()
"""
# TODO refactor with polyCreatePolygon
poly = pg.Mesh(dim=2, isGeometry=True)
startPos = pg.RVector3(start)
endPos = pg.RVector3(end)
a = endPos - startPos

dt = 1. / segments
left = kwargs.pop('leftDirection', True)

for i in range(0, segments + 1):
if left:
p = startPos + a * (dt * i)
else:
p = endPos - a * (dt * i)

poly.createNode(p)

_polyCreateDefaultEdges(poly, isClosed=False, **kwargs)
return poly

**kwargs):
"""Create a polygon from a list of vertices.

All vertices need to be unique and duplicate vertices will be ignored.
If you want the polygon be a closed region you can set the 'isClosed' flag.
Closed region can be attributed by assigning a region marker.
The automatic region marker is placed in the center of all vertices.

Parameters
----------
verts : []
* List of x y pairs [[x0, y0], ... ,[xN, yN]]

isClosed : bool [True]
Add closing edge between last and first node.

interpolate : str ['linear']
Interpolation rule for addnodes. 'linear' or 'spline'. TODO 'harmfit'

**kwargs:

marker : int [None]
Marker for the resulting triangle cells after mesh generation.
markerPosition : floats [x, y] [0.0, 0.0]
Position (absolute) of the marker (works for both regions and
holes)
area : float [0]
Maximum cell size for resulting triangles after mesh generation
isHole : bool [False]
The polygon will become a hole instead of a triangulation
boundaryMarker : int [1]
Marker for the resulting boundary edges
leftDirection : bool [True]
Rotational direction

Returns
-------
poly : :gimliapi:GIMLI::Mesh
The resulting polygon is a :gimliapi:GIMLI::Mesh.

Examples
--------
>>> # no need to import matplotlib, pygimli show does.
>>> import pygimli as pg
>>> import pygimli.meshtools as mt
>>> p1 = mt.createPolygon([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0]],
...                       isClosed=True, marker=3, area=0.1)
>>> p2 = mt.createPolygon([[0.3, 0.15], [0.85, 0.15], [0.85, 0.7]],
...                       isClosed=True, isHole=True)
>>> p3 = mt.createPolygon([[-0.1, 0.2], [-1.1, 0.2], [-1.1, 1.2], [-0.1, 1.2]],
>>> p4 = mt.createPolygon([[-0.1, 0.2], [-1.1, 0.2], [-1.1, 1.2], [-0.1, 1.2]],
...                       marker=4)
>>> ax, _ = pg.show(mt.mergePLC([p1, p2, p3, p4]), showNodes=True)
>>> pg.wait()
"""
poly = pg.Mesh(dim=2, isGeometry=True)

if isClosed:
verts = np.array(verts)
verts = np.vstack([verts, verts[0]])

tV = pg.utils.cumDist(verts)
tI = []

for i, t in enumerate(tV[0:len(tV)-1]):
tI.append(t)
tI.append(tV[i] + dt*(j+1))

if not isClosed:
tI.append(tV[-1])

verts = pg.meshtools.interpolate(verts, tI,
method=interpolate,
periodic=isClosed)

if kwargs.pop("leftDirection", False):
for v in verts[::-1]:
if isinstance(v, float) or isinstance(v, int):
poly.createNodeWithCheck([v, 0], warn=True)
else:
poly.createNodeWithCheck(v, warn=True)
else:
for v in verts:
if isinstance(v, float) or isinstance(v, int):
poly.createNodeWithCheck([v, 0], warn=True)
else:
poly.createNodeWithCheck(v, warn=True)

_polyCreateDefaultEdges(poly, isClosed=isClosed,
boundaryMarker=kwargs.pop('boundaryMarker', 1))

if isClosed:
setPolyRegionMarker(poly, **kwargs)

return poly

[docs]def mergePLC(plcs, tol=1e-3):
"""Merge multiply polygons.

Merge multiply polygons into a single polygon.
Common nodes and common edges will be checked and removed.
When a node touches an edge, the edge will be splited.

3D only OOC with polytools

TODO:
* Crossing or Node/Edge intersections will NOT be
recognized yet.
* Edge on Node touch

Parameters
----------
plcs: [:gimliapi:GIMLI::Mesh]
List of PLC that want to be merged into one new PLC

tol : double
Tolerance to check for duplicated nodes. [1e-3]

Returns
-------
plc : :gimliapi:GIMLI::Mesh
The resulting polygon is a :gimliapi:GIMLI::Mesh.

Examples
--------
>>> import pygimli as pg
>>> import pygimli.meshtools as mt
>>> from pygimli.viewer.mpl import drawMesh
>>> world = mt.createWorld(start=[-10, 0], end=[10, -10], marker=1)
>>> c1 = mt.createCircle([-1, -4], radius=1.5, area=0.1,
...                       marker=2, segments=5)
>>> c2 = mt.createCircle([-6, -5], radius=[1.5, 3.5], isHole=1)
>>> r1 = mt.createRectangle(pos=[3, -5], size=[2, 2], marker=3)
>>> r2 = mt.createRectangle(start=[4, -4], end=[6, -6],
...                          marker=4, area=0.1)
>>> plc = mt.mergePLC([world, c1, c2, r1, r2])
>>> fig, ax = pg.plt.subplots()
>>> drawMesh(ax, plc)
>>> drawMesh(ax, mt.createMesh(plc))
>>> pg.wait()
"""
if plcs[0].dim() == 3:
return mergePLC3D(plcs, tol)

tmp = pg.optImport('tempfile')
names = []
for p in plcs:
_, namePLC = tmp.mkstemp(suffix='.poly')
pg.meshtools.exportPLC(p, namePLC)
names.append(namePLC)

for n in names[1:]:
syscal = 'polyMerge {0} {1} {0}'.format(names[0], n)
pg.debug(syscal)
os.system(syscal)

for n in names:
try:
pg.debug('Remove:', n)
os.remove(n)
except:
print("can't remove:", n)

return plc

# handle 2D geometries
plc = pg.Mesh(dim=2, isGeometry=True)

for p in plcs:
nodes = []
for n in p.nodes():
nn = plc.createNodeWithCheck(n.pos(), tol,
warn=False, edgeCheck=True)
if n.marker() != 0:
nn.setMarker(n.marker())
nodes.append(nn)

for e in p.boundaries():
plc.createEdge(nodes[e.node(0).id()], nodes[e.node(1).id()],
e.marker())

if len(p.regionMarkers()) > 0:
for rm in p.regionMarkers():

if len(p.holeMarker()) > 0:
for hm in p.holeMarker():

return plc

[docs]def mergePLC3D(plcs, tol=1e-3):
"""Experimental replacement for polyMerge. Don't expect too much.
"""
if len(plcs) < 2:
pg.critical("Give at least 2 PLCs.")

if plcs[0].dim() != 3:
pg.warn("2D poly found. redirect to mergePLC")
return mergePLC(plcs, tol)

# first try. merge all into p0 = plcs[0]
#  * will only work if all faces of plcs[1:] does not match any face of p0
#  * or all matching plcs[1:] are lie completely within p0

p0 = pg.Mesh(plcs[0])

for p in plcs[1:]:
for b in p.boundaries():
p0.copyBoundary(b)

if len(p.regionMarkers()) > 0:
for rm in p.regionMarkers():

for hm in p.holeMarker():

return p0

"""API change here .. use createParaMeshPLC instead."""
pg.deprecated("use createParaMeshPLC")
return createParaMeshPLC(*args, **kwargs)

paraMaxCellSize=0.0, boundary=-1, boundaryMaxCellSize=0,
balanceDepth=True,
"""Create a PLC mesh for an inversion parameter mesh.

Create a PLC mesh for an inversion parameter mesh with for a given list of
sensor positions. Sensor positions are assumed to be on the surface and
must be sorted and unique.

You can create a parameter mesh without sensors if you just set [xMin,
xMax] as sensors.

The PLC is a :gimliapi:GIMLI::Mesh and contain nodes, edges and two
region markers, one for the parameters domain (marker=2) and a larger
boundary around the outside (marker=1)

TODO:
* spline interpolation between sensorpoints or addpoints for non closed
* subsurface sensors (partly .. see example)

Parameters
----------
sensors : [RVector3] | DataContainer with sensorPositions() | [xMin, xMax]
Sensor positions. Must be sorted and unique in positive x direction.
Depth need to be y-coordinate.

Relative distance for refinement nodes between two sensors (1=none),
e.g., 0.5 means 1 additional node between two neighboring sensors
e.g., 0.33 means 2 additional equidistant nodes between two sensors

Maximum depth for parametric domain, 0 (default) means 0.4 * maximum
sensor range.
balanceDepth: bool [True]
Equal depth for the parametric domain.

paraBoundary : float, optional
Margin for parameter domain in absolute sensor distances. 2 (default).

paraMaxCellSize: double, optional
Maximum size for parametric size in m*m

boundaryMaxCellSize: double, optional
Maximum cells size in the boundary region in m*m

boundary : float, optional
Boundary width to be appended for domain prolongation in absolute
para domain width.
Values lover 0 force the boundary to be 4 times para domain width.

isClosed : bool [False]
Create a closed geometry from sensor positions.
Region marker is 1. Boundary marker is -1 (homogeneous Neumann)

Returns
-------
poly: :gimliapi:GIMLI::Mesh
piecewise linear complex (PLC) containing nodes and edges

Examples
--------
>>> # no need to import matplotlib, pygimli show does.
>>> import pygimli as pg
>>> import pygimli.meshtools as plc
>>> # Create the simplest paramesh PLC with a para box of 10 m without
>>> # sensors
>>> p = plc.createParaMeshPLC([0,10])
>>> # you can add subsurface sensors now with
>>> for z in range(1,4):
...     n = p.createNode((5,-z), -99)
>>> ax,_ = pg.show(p)
"""
if isClosed:
boundaryMarker=-1, marker=1,
area=paraMaxCellSize, **kwargs)
return plc

noSensors = False
if hasattr(sensors, 'sensorPositions'):  # obviously a DataContainer type
sensors = sensors.sensorPositions()
elif isinstance(sensors, np.ndarray):
if sensors.ndim == 1:
sensors = [pg.RVector3(s, 0) for s in sensors]
else:  # assume 2d array with 2 or 3 values per item
sensors = [pg.RVector3(s) for s in sensors]
elif isinstance(sensors, list):
if len(sensors) == 2:
# guess we have just a desired Pbox with
sensors = [pg.RVector3(sensors[0], 0.0),
pg.RVector3(sensors[1], 0.0)]
noSensors = True
paraBoundary = 0

eSpacing = kwargs.pop('eSpacing', sensors[0].distance(sensors[1]))

iz = 1
xMin, yMin, zMin = sensors[0][0], sensors[0][1], sensors[0][2]
xMax, yMax, zMax = xMin, yMin, zMin
for e in sensors:
xMin = min(xMin, e[0])
xMax = max(xMax, e[0])
yMin = min(yMin, e[1])
yMax = max(yMax, e[1])
zMin = min(zMin, e[2])
zMax = max(zMax, e[2])

if abs(yMin) < 1e-8 and abs(yMax) < 1e-8:
iz = 2

paraBound = eSpacing * paraBoundary

paraDepth = 0.4 * (xMax - xMin)

poly = pg.Mesh(dim=2, isGeometry=True)
# define para domain without surface
n1 = poly.createNode([xMin - paraBound, sensors[0][iz]])
if balanceDepth:
n2 = poly.createNode([xMin - paraBound, bD])
n3 = poly.createNode([xMax + paraBound, bD])
else:
n2 = poly.createNode([xMin - paraBound, sensors[0][iz] - paraDepth])
n3 = poly.createNode([xMax + paraBound, sensors[-1][iz] - paraDepth])
n4 = poly.createNode([xMax + paraBound, sensors[-1][iz]])

if boundary < 0:
boundary = 4

bound = abs(xMax - xMin) * boundary
if bound > paraBound:
# define world without surface
n11 = poly.createNode(n1.pos() - [bound, 0.])
n12 = poly.createNode(n11.pos() - [0., bound + paraDepth])
n14 = poly.createNode(n4.pos() + [bound, 0.])
n13 = poly.createNode(n14.pos() - [0., bound + paraDepth])

poly.createEdge(n1,  n11, pg.core.MARKER_BOUND_HOMOGEN_NEUMANN)
poly.createEdge(n11, n12, pg.core.MARKER_BOUND_MIXED)
poly.createEdge(n12, n13, pg.core.MARKER_BOUND_MIXED)
poly.createEdge(n13, n14, pg.core.MARKER_BOUND_MIXED)
poly.createEdge(n14, n4, pg.core.MARKER_BOUND_HOMOGEN_NEUMANN)
poly.addRegionMarker(n12.pos() + [1e-3, 1e-3], 1, boundaryMaxCellSize)

poly.createEdge(n1, n2, 1)
poly.createEdge(n2, n3, 1)
poly.createEdge(n3, n4, 1)
poly.addRegionMarker(n2.pos() + [1e-3, 1e-3], 2, paraMaxCellSize)

# define surface
nSurface = []
nSurface.append(n1)

if not noSensors:
for i, e in enumerate(sensors):
if iz == 2:
e.rotateX(-math.pi / 2)
nSurface.append(poly.createNode(e, pg.core.MARKER_NODE_SENSOR))
if i < len(sensors) - 1:
e1 = sensors[i + 1]
if iz == 2:
e1.rotateX(-math.pi / 2)
nSurface.append(poly.createNode((e + e1) * 0.5))
# print("Surface add ", e, el, nSurface[-2].pos(),
#        nSurface[-1].pos())
if i > 0:
e1 = sensors[i - 1]
if iz == 2:
e1.rotateX(-math.pi / 2)
nSurface.append(poly.createNode(e - (e - e1) * paraDX))
nSurface.append(poly.createNode(e, pg.core.MARKER_NODE_SENSOR))
if i < len(sensors) - 1:
e1 = sensors[i + 1]
if iz == 2:
e1.rotateX(-math.pi / 2)
nSurface.append(poly.createNode(e + (e1 - e) * paraDX))
# print("Surface add ", nSurface[-3].pos(), nSurface[-2].pos(),
#        nSurface[-1].pos())
nSurface.append(n4)

for i in range(len(nSurface) - 1, 0, -1):
poly.createEdge(nSurface[i], nSurface[i - 1],
pg.core.MARKER_BOUND_HOMOGEN_NEUMANN)

return poly

Read 2D :term:Triangle or 3D :term:Tetgen PLC files.

Parameters
----------
filename: string
Filename *.poly

comment: string ('#')
String containing all characters that define a comment line. Identified
lines will be ignored during import.

Returns
-------
poly :
:gimliapi:GIMLI::Mesh
"""
with open(filename, 'r') as fi:

# Filter comment lines
comment_lines = []
for i, line in enumerate(content):
if line[0] in comment:
comment_lines.append(i)
for j in comment_lines[::-1]:
del(content[j])

fromOne = 0

poly = pg.Mesh(dim=dimension, isGeometry=False)
# isGeometry forces expensive checks: we assume the plc is valid so we set
# this flag in the end

# Nodes section
for i in range(nVerts):
row = content[1 + i].split('\r\n')[0].split()

if len(row) == (1 + dimension + nPointsAttributes + haveNodeMarker):
if i == 0:
fromOne = int(row[0])
if dimension == 2:
n = poly.createNode((float(row[1]), float(row[2])))
elif dimension == 3:
n = poly.createNode((float(row[1]), float(row[2]),
float(row[3])))
if haveNodeMarker:
n.setMarker(int(row[-1]))

else:
print(i, len(row), row,
(1 + dimension + nPointsAttributes + haveNodeMarker))
raise Exception("Poly file seams corrupt: node section line: " +
content[1 + i])

# Segment section
row = content[1 + nVerts].split()

if len(row) != 2:
raise Exception("Format unknown for segment section " + row)

nSegments = int(row[0])
haveBoundaryMarker = int(row[1])

if dimension == 2:
for i in range(nSegments):
row = content[2 + nVerts + i].split()

if len(row) == (3 + haveBoundaryMarker):
marker = 0
if haveBoundaryMarker:
marker = int(row[3])

poly.createEdge(
poly.node(int(row[1]) - fromOne),
poly.node(int(row[2]) - fromOne), marker)
else:
segment_offset = 0
for i in range(nSegments):
row = content[2 + nVerts + i + segment_offset].split()
numBounds = int(row[0])
numHoles = int(row[1])
# if numHoles != '0':
#     pg.error("Can't handle 3D faces with holes yet")
marker = 0
if haveBoundaryMarker:
marker = int(row[2])

face = None
for k in range(numBounds):
boundRow = content[2 + nVerts + i + segment_offset + 1]\
.split()
# nNodes = int(boundRow[0])
nodeIdx = [int(_b) for _b in boundRow[1:]]

if k == 0:
face = poly.createPolygonFace(poly.nodes(nodeIdx),
marker=marker, check=True)
else:
if len(nodeIdx) == 2:
if nodeIdx[0] == nodeIdx[1]:
else:

segment_offset += 1

for k in range(numHoles):
r = content[2 + nVerts + i + segment_offset + 1]\
.split()

segment_offset += 1

nSegments += segment_offset

# Hole section
row = content[2 + nVerts + nSegments].split()

if len(row) != 1:
raise Exception("Format unknown for hole section " + row)

nHoles = int(row[0])
for i in range(nHoles):
row = content[3 + nVerts + nSegments + i].split()
if len(row) == 3:
elif len(row) == 4 and dimension == 3:
else:
raise Exception("Poly file seams corrupt: hole section line (3):" +
row + " : " + str(i) + " " + str(len(row)))

if (3 + nVerts + nSegments + nHoles) < len(content):
# Region section
row = content[3 + nVerts + nSegments + nHoles].split()

if len(row) != 1:
raise Exception("Format unknown for region section " + row)

nRegions = int(row[0])

for i in range(nRegions):
row = content[4 + nVerts + nSegments + nHoles + i].split()
if len(row) == 5:
marker=int(float(row[3])),
area=float(row[4]))
elif len(row) == 6 and dimension == 3:
float(row[3])],
marker=int(float(row[4])),
area=float(row[5]))
else:
raise Exception("Poly file seams corrupt: region section " +
"line (5): " + str(i) + " " + str(len(row)))

poly.setGeometry(True)
return poly

[docs]def exportPLC(poly, fname, **kwargs):
r"""General writer to save piece-wise linear complex (PLC) as poly file.

Choose from poly.dimension() and forward appropriate to
:gimliapi:GIMLI::Mesh::exportAsTetgenPolyFile
and :py:mod:pygimli.meshtools.writeTrianglePoly

Parameters
----------
poly : :gimliapi:GIMLI::Mesh
The polygon to be written.

fname : string
Filename of the file to write (\\*.n, \\*.e).

Examples
--------
>>> import pygimli as pg
>>> import tempfile, os
>>> fname = tempfile.mktemp() # Create temporary string for filename.
>>> world2d = pg.meshtools.createWorld(start=[-10, 0], end=[20, -10])
>>> pg.meshtools.exportPLC(world2d, fname)
Mesh: Nodes: 4 Cells: 0 Boundaries: 4
>>> world3d = pg.createGrid([0, 1], [0, 1], [-1, 0])
>>> pg.meshtools.exportPLC(world3d, fname)
>>> os.remove(fname)
"""
if poly.dimension() == 2:
exportTrianglePoly(poly, fname, **kwargs)
else:
exportTetgenPoly(poly, fname, **kwargs)

[docs]def writePLC(*args, **kwargs):
"""
Backward compatibility.
Please use :py:mod:pygimli.meshtools.exportPLC.
"""
pg.deprecated('use exportPLC')  # 16.08.2019
return exportPLC(*args, **kwargs)

def exportTrianglePoly(poly, fname, float_format='.15e'):
r"""Write :term:Triangle poly.

Write :term:Triangle :cite:Shewchuk96b ASCII file.
See: ://www.cs.cmu.edu/~quake/triangle.html

Parameters
----------
poly : :gimliapi:GIMLI::Mesh
mesh PLC holding nodes, edges, holes & regions

fname : string
Filename of the file to read (\\*.n, \\*.e)

float_format : string
format string for floats according to str.format()

verbose : boolean [False]
Be verbose during import.
"""
if float_format[0] != '{':
pfmt = '{:' + float_format + '}'
else:
pfmt = float_format
with open(fname, 'w') as fid:
fid.write('{:d}\t2\t0\t1\n'.format(poly.nodeCount()))
nm = poly.nodeMarkers()
bm = poly.boundaryMarkers()

fmt = '{:d}' + ('\t' + pfmt) * 2 + '\t{:d}\n'
for i, p in enumerate(poly.positions()):
fid.write(fmt.format(i, p.x(), p.y(), nm[i]))
fid.write('{:d}\t1\n'.format(poly.boundaryCount()))

for i, b in enumerate(poly.boundaries()):
fid.write('{:d}\t{:d}\t{:d}\t{:d}\n'.format(i, b.node(0).id(),
b.node(1).id(), bm[i]))
fid.write('{:d}\n'.format(len(poly.holeMarker())))

fmt = '{:d}' + ('\t' + pfmt) * 2 + '\n'
for i, h in enumerate(poly.holeMarker()):
fid.write(fmt.format(i, h.x(), h.y()))
fid.write('{:d}\n'.format(len(poly.regionMarkers())))

fmt = '{:d}' + ('\t' + pfmt) * 3 + '\t{:.15e}\n'
for i, r in enumerate(poly.regionMarkers()):
fid.write(fmt.format(i, r.x(), r.y(), r.marker(), r.area()))

return

def writeTrianglePoly(*args, **kwargs):
""" Backward compatibility.
Please use :py:mod:pygimli.meshtools.exportTrianglePoly.
"""
return exportTrianglePoly(*args, **kwargs)

def exportTetgenPoly(poly, filename, float_format='.12e', **kwargs):
r"""
Writes a given piecewise linear complex (mesh/poly) into a Ascii file in
:term:Tetgen .poly format.

Parameters
----------
filename: string
Name in which the result will be written. The recommended file
ending is '.poly'.

poly: :gimliapi:GIMLI::Mesh
Piecewise linear complex as :gimliapi:GIMLI::Mesh to be exported.

float_format: format string ('.12e')
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.

kwargs:
* extraBoundaries:

"""
if filename[-5:] != '.poly':
filename = filename + '.poly'
polytxt = ''
sep = '\t'  # standard tab seperated file
linesep = '\n'  # os.linesep does not work in mingwshell, testit!!
assert poly.dim() == 3, 'Exit, only for 3D meshes.'
boundary_marker = 1
attribute_count = 0

# Part 1/4: node list
# intro line
# <nodecount> <dimension (3)> <# of attributes> <boundary markers (0 or 1)>
polytxt += '{0}{5}{1}{5}{2}{5}{3}{4}'.format(poly.nodeCount(), 3,
attribute_count,
boundary_marker,
linesep, sep)
# loop over positions, attributes and marker(node)
# <point idx> <x> <y> <z> [attributes] [boundary marker]
point_str = '{:d}'  # index of the point
for i in range(3):
# coords as float with given precision
point_str += sep + '{:%s}' % (float_format)
point_str += sep + '{:d}' + linesep  # node marker
for j, node in enumerate(poly.nodes()):
fill = [node.id()]
fill.extend([pos for pos in node.pos()])
fill.append(node.marker())
polytxt += point_str.format(*fill)

# Part 2/4: boundary list
# intro line
# <# of facets> <boundary markers (0 or 1)>
nBoundaries = poly.boundaryCount()
# look for extra boundaries present in either the PLC or in kwargs
extraBoundaries = []
if 'extraBoundaries' in kwargs:
extraBoundaries += kwargs.pop('extraBoundaries', [])

if hasattr(poly, 'extraBoundaries'):
extraBoundaries += poly.extraBoundaries

if len(extraBoundaries) > 0:
print("Detected ", len(extraBoundaries), " extra boundaries!")

nBoundaries += len(extraBoundaries)
polytxt += '{0:d}{2}1{1}'.format(nBoundaries, linesep, sep)
# loop over facets, each facet can contain an arbitrary number of holes
# and polygons, in our case, there is always one polygon per facet.

hole_str = '{:d}'
for m in range(3):
hole_str += sep + '{:%s}' % float_format

hole_str += linesep

for bound in poly.boundaries():
# one line per facet
# <# of polygons> [# of holes] [boundary marker]
try:
nSubs = bound.subfaceCount()
except BaseException:
nSubs = 0
try:
nHoles = len(bound.holeMarkers())
except BaseException:
nHoles = 0

npolys = 1 + nSubs + len(bound.secondaryNodes())
polytxt += '{3}{2}{4}{2}{0:d}{1}'.format(bound.marker(), linesep,
sep, npolys, nHoles)
# inner loop over polygons
# <# of corners> <corner 1> <corner 2> ... <corner #>
for l in range(1):
poly_str = '{:d}'.format(bound.nodeCount())
poly_str += sep + sep.join(['{:d}'.format(n) for n in bound.ids()])
polytxt += '{0}{1}'.format(poly_str, linesep)

# loop over subfaces
for l in range(nSubs):
sub = bound.subface(l)
poly_str = '{:d}'.format(len(sub))
poly_str += sep + sep.join(['{:d}'.format(n.id()) for n in sub])
polytxt += '{0}{1}'.format(poly_str, linesep)

# inner loop over holes
if nHoles > 0:
for n, hole in enumerate(bound.holeMarkers()):
polytxt += hole_str.format(n, *hole)

# not necessary yet ?! why is there an extra hole section?
# because this is for 2D holes in facets only

# loop over secondaryNodes add them as single points
for l in range(len(bound.secondaryNodes())):
ind = bound.secondaryNodes()[l].id()
poly_str = '{:d}'.format(2)
poly_str += sep + '{0:d} {0:d}'.format(ind)
polytxt += '{0}{1}'.format(poly_str, linesep)

# part 2b: extra boundaries that cannot be part of mesh class
for nodes in extraBoundaries:
# <# of polygons> [# of holes] [boundary marker]
npolys = 1
polytxt += '1{2}0{2}{0:d}{1}'.format(111, linesep, sep)
# <# of corners> <corner 1> <corner 2> ... <corner #>
poly_str = '{:d}'.format(len(nodes))
for ind in nodes:
poly_str += sep + '{:d}'.format(ind)

polytxt += '{0}{1}'.format(poly_str, linesep)

# part 3/4: hole list
# intro line
# <# of holes>
holes = poly.holeMarker()
polytxt += '{:d}{}'.format(len(holes), linesep)
# loop over hole markers
# <hole #> <x> <y> <z>

for n, hole in enumerate(holes):
polytxt += hole_str.format(n, *hole)

# part 4/4: region attributes and volume constraints (optional)
# intro line
# <# of regions>
regions = poly.regionMarkers()
polytxt += '{:d}{}'.format(len(regions), linesep)
# loop over region markers
# <region #> <x> <y> <z> <region number> <region attribute>
region_str = '{:d}'
for o in range(3):
region_str += sep + '{:%s}' % (float_format)

region_str += sep + '{:d}%s{:%s}' % (sep, float_format) + linesep
for p, region in enumerate(regions):
polytxt += region_str.format(p, region.x(), region.y(), region.z(),
region.marker(),
region.area())

# writing file
with open(filename, 'w') as out:
out.write(polytxt)
out.close()

[docs]def syscallTetgen(filename, quality=1.2, area=0, preserveBoundary=False,
verbose=False, tetgen='tetgen'):
"""Create a mesh with :term:Tetgen from file.

Create a :term:Tetgen :cite:Si2004 mesh from a PLC.

Forwards to system call tetgen, which must be known to your system.

Parameters
----------
filename: str

quality: float [1.2]
Refines mesh (to improve mesh quality). [1.1 ... ]

area: float [0.0]
Maximum cell size (m³)

preserveBoundary: bool [False]
Preserve PLC boundary mesh

verbose: bool [False]
be verbose

tetgen: str | path ['tetgen']
Binary for tetgen. Given as complete path or simple the binary name
if its known in the system path.

Returns
-------
mesh : :gimliapi:GIMLI::Mesh
"""
filebody = filename.replace('.poly', '')

# tetgen -pazVAC -q1.2 \$MESH
# test -O2
syscal = tetgen + ' -pzAC'

if area > 0:
syscal += 'a' + str(area)
else:
syscal += 'a'

syscal += 'q' + str(quality)

if not verbose:
syscal += 'Q'
else:
pass
# syscal += 'V'

if preserveBoundary:
syscal += 'Y'

syscal += ' ' + filebody + '.poly'

if verbose:
print(syscal)
pg.debug(syscal)

system(syscal)

mesh = None
if os.path.isfile(filebody + '.1.node'):
# system('meshconvert -it -BD -o ' + filebody + ' ' + filebody + '.1')
try:
os.remove(filebody + '.1.node')
os.remove(filebody + '.1.ele')
os.remove(filebody + '.1.face')
except BaseException as e:
print(e)
else:
# system('meshconvert -it -BD -o ' + filebody + ' ' + filebody + '-1')
try:
os.remove(filebody + '-1.node')
os.remove(filebody + '-1.ele')
os.remove(filebody + '-1.face')
except BaseException as e:
print(e)

# mesh = pg.Mesh(filebody)
return mesh

def polyCreateWorld(filename, x=None, depth=None, y=None, marker=0,
maxCellSize=0, verbose=True):
"""Create the PLC of a default world.

Create the PLC of a default world.

Out of core wrapper for dcfemlib::polytools::polyCreateWorld

# TODO needs to be converted to the Python-only tools.

Parameters
----------

Returns
-------
"""
if depth is None:
return

if x is None:
return

dimension = 3
z = depth

if y is None:
dimension = 2

syscal = 'polyCreateWorld -d ' + str(dimension) \
+ ' -x ' + str(x) \
+ ' -y ' + str(y) \
+ ' -z ' + str(z) \
+ ' -m ' + str(marker) \

if maxCellSize > 0:
syscal += " -a " + str(maxCellSize)

syscal = syscal + ' ' + filename

if verbose:
print(syscal)

os.system(syscal)

[docs]def createFacet(mesh, boundaryMarker=None, verbose=True):
"""Create a coplanar PLC of a 2d mesh or poly

TODO:
* mesh with cell into plc with boundaries
* poly account for inner edges

"""
if mesh.dimension() != 2:
pg.error("need two dimensional mesh or poly")

if mesh.cellCount() > 0:
pg.critical("Implementme")

poly = pg.Mesh(dim=3, isGeometry=True)

nodes = [poly.createNode(n.pos()).id() for n in mesh.nodes()]

if boundaryMarker is None:
for rm in mesh.regionMarkers():
boundaryMarker = rm.marker()
continue

b = poly.createBoundary(nodes, marker=boundaryMarker or 0)

for h in mesh.holeMarker():

return poly

[docs]def createCube(size=[1.0, 1.0, 1.0],
pos=None, rot=None, boundaryMarker=0, **kwargs):
"""Create cube PLC as geometrie definition.

Parameters
----------
size : [x, y, z]
x, y, and z-size of the cube. Default = [1.0, 1.0, 1.0] in m
pos : pg.Pos [None]
The center position, default is at the origin.
rot : pg.Pos [None]
Rotate on the center.
boundaryMarker : int[0]
Boundary marker for the resulting faces.

** kwargs:
Marker related arguments:
See :py:mod:pygimli.meshtools.polytools.setPolyRegionMarker

Examples
--------
>>> import pygimli.meshtools as mt
>>> cube = mt.createCube()
>>> print(cube)
Mesh: Nodes: 8 Cells: 0 Boundaries: 6
>>> cube = mt.createCube([10, 10, 1])
>>> print(cube.bb())
[[-5.0, -5.0, -0.5], [5.0, 5.0, 0.5]]
>>> cube = mt.createCube([10, 10, 1], pos=[-4.0, 0.0, 0.0])
>>> print(pg.center(cube.positions()))
RVector3: (-4.0, 0.0, 0.0)

Returns
-------
poly : :gimliapi:GIMLI::Mesh
The resulting polygon is a :gimliapi:GIMLI::Mesh.

"""
poly = pg.Mesh(3, isGeometry=True)

for y in [-0.5, 0.5]:
poly.createNode(-0.5, y, -0.5)
poly.createNode(+0.5, y, -0.5)
poly.createNode(+0.5, y, +0.5)
poly.createNode(-0.5, y, +0.5)

faces = [[4, 5, 1, 0],
[5, 6, 2, 1],
[6, 7, 3, 2],
[7, 4, 0, 3],
[0, 1, 2, 3],
[7, 6, 5, 4], ]

if isinstance(boundaryMarker, list):
for i, f in enumerate(faces):
poly.createPolygonFace(poly.nodes(f), marker=boundaryMarker[i])
else:
for f in faces:
poly.createPolygonFace(poly.nodes(f), marker=boundaryMarker)

poly.scale(size)

if rot is not None:
poly.rotate(rot)

if pos is not None:
poly.translate(pos)

setPolyRegionMarker(poly, **kwargs)

return poly

def extrude(p2, z=-1.0, boundaryMarker=0, **kwargs):
"""Create 3D body by extruding a closed 2D poly into z direction

Parameters
----------
p2 : :gimliapi:GIMLI::Mesh
2D geometry

z : float [-1.0]
2D geometry

Keyword Arguments
-----------------
** kwargs:
Marker related arguments:
See :py:mod:pygimli.meshtools.polytools.setPolyRegionMarker

Returns
-------
poly : :gimliapi:GIMLI::Mesh
The resulting polygon is a :gimliapi:GIMLI::Mesh.
"""
if p2.dimension() != 2:
pg.error("need two dimensional mesh or poly")

if p2.cellCount() > 0:
pg.critical("Implementme")

poly = pg.Mesh(3, isGeometry=True)
top = []
for n in p2.nodes():
top.append(poly.createNode(n.pos()).id())

bot = []
for n in p2.nodes():
bot.append(poly.createNode(n.pos() + [0.0, 0.0, z]).id())
N = len(top)

poly.createPolygonFace(poly.nodes(top), marker=boundaryMarker)
poly.createPolygonFace(poly.nodes(bot[::-1]), marker=boundaryMarker)

for i in range(len(top)):
poly.createPolygonFace(poly.nodes([i, N + i,
N + (i + 1) % N,
(i+1) % N]),
marker=boundaryMarker)

setPolyRegionMarker(poly, **kwargs)

return poly

pos=None, rot=None, boundaryMarker=0, **kwargs):
"""Create PLC of a cylinder.

Out of core wrapper for dcfemlib::polytools.

Note, there is a bug in the old polytools which ignores the area settings
for marker == 0.

Parameters
----------

height : float
Height of the cylinder

nSegments : int [8]
Number of segments of the cylinder.

pos : pg.Pos [None]
The center position, default is at the origin.

Keyword Arguments
----------------
** kwargs:
Marker related arguments:
See :py:mod:pygimli.meshtools.polytools.setPolyRegionMarker

Returns
-------
poly : :gimliapi:GIMLI::Mesh
The resulting polygon is a :gimliapi:GIMLI::Mesh.

"""
poly = extrude(circ, z=height, boundaryMarker=boundaryMarker, **kwargs)
# move it to z=0
poly.translate([0.0, 0.0, -height/2])

if rot is not None:
c = pg.center(poly.positions())
poly.translate(-c)
poly.rotate(rot)
poly.translate(c)

if pos is not None:
poly.translate(pos)

return poly

def boundaryPlaneIntersectionLines(boundaries, plane):
"""Create Lines from boundaries that intersect a plane."""
lines = []

for b in boundaries:
ps = []
for i, n in enumerate(b.shape().nodes()):
line = pg.Line(n.pos(), b.shape().node(
(i + 1) % b.shape().nodeCount()).pos())
p = plane.intersect(line, 1e-8, True)
if p.valid():
ps.append(p)

if len(ps) == 2:
lines.append(list(zip([ps[0].x(), ps[1].x()],
[ps[0].z(), ps[1].z()])))
return lines

if __name__ == "__main__":
pass