Commit 07604102 authored by Daniel Boeckenhoff's avatar Daniel Boeckenhoff
Browse files

mesh3D.cut - before the idea to manipulate fields in _cut_sympy to append the...

mesh3D.cut - before the idea to manipulate fields in _cut_sympy to append the 'scalarMap' and remove it in the end.
parent 06041647
......@@ -473,9 +473,13 @@ class Tensors(AbstractNdarray):
bases.append(t.coordSys)
except AttributeError:
pass
# get most frequent coordSys
coordSys = sorted(bases, key=Counter(bases).get, reverse=True)[0]
kwargs['coordSys'] = coordSys
if bases:
# get most frequent coordSys
coordSys = sorted(bases, key=Counter(bases).get, reverse=True)[0]
kwargs['coordSys'] = coordSys
else:
default = cls.__slot_defaults__[cls.__slots__.index('coordSys')]
kwargs['coordSys'] = default
''' transform all raw inputs to cls type with correct coordSys. Also
automatically make a copy of those instances that are of the correct
......@@ -1235,6 +1239,25 @@ class TensorFields(Tensors):
for i, field in enumerate(item.fields):
self.fields[i].__setitem__(index, field)
@classmethod
def merged(cls, *objects, **kwargs):
if not all([isinstance(o, cls) for o in objects]):
# TODO: could allow if all faceScalars are none
raise TypeError("Merge constructor only accepts {cls} instances."
.format(**locals()))
inst = super(TensorFields, cls).merged(*objects, **kwargs)
fields = []
if all([len(obj.fields) == len(objects[0].fields)
for obj in objects]):
for fld_idx in range(len(objects[0].fields)):
field = tfields.Tensors.merged(*[obj.fields[fld_idx]
for obj in objects])
fields.append(field)
inst = cls.__new__(cls, inst, *fields)
return inst
def transform(self, coordSys):
super(TensorFields, self).transform(coordSys)
for field in self.fields:
......@@ -1377,10 +1400,10 @@ class TensorMaps(TensorFields):
item = super(TensorMaps, self).__getitem__(index)
try:
if issubclass(type(item), TensorMaps):
item.maps = [mp.copy() for mp in item.maps]
if isinstance(index, tuple):
index = index[0]
if item.maps:
item.maps = [mp.copy() for mp in item.maps]
indices = np.array(range(len(self)))
keep_indices = indices[index]
if isinstance(keep_indices, int):
......
import numpy as np
def ensure_complex(*base_vectors):
# ensure, that the third entry in base_vector of type tuple becomes a complex type
base_vectors = list(base_vectors)
for i, vector in enumerate(base_vectors):
if isinstance(vector, tuple):
if len(vector) == 3:
vector = list(vector)
vector[2] = complex(vector[2])
base_vectors[i] = tuple(vector)
return base_vectors
def to_base_vectors(*base_vectors):
"""
Transform tuples to arrays with np.mgrid
Args:
tuple of lenght 3 with complex third entry -> start, stop, n_steps
Returns:
list if np.array for each base
"""
base_vectors = list(base_vectors)
# transform tuples to arrays with mgrid
for i in range(len(base_vectors)):
if isinstance(base_vectors[i], tuple):
base_vectors[i] = np.mgrid[slice(*base_vectors[i])]
if isinstance(base_vectors[i], list):
base_vectors[i] = np.array(base_vectors[i])
return base_vectors
def igrid(*base_vectors, **kwargs):
"""
Args:
......@@ -55,15 +85,7 @@ def igrid(*base_vectors, **kwargs):
"""
iter_order = kwargs.pop('iter_order', np.arange(len(base_vectors)))
# ensure, that the third entry in base_vector of type tuple becomes a complex type
base_vectors = list(base_vectors)
for i, vector in enumerate(base_vectors):
if isinstance(vector, tuple):
if len(vector) == 3:
vector = list(vector)
vector[2] = complex(vector[2])
base_vectors[i] = tuple(vector)
base_vectors = ensure_complex(*base_vectors)
# create the grid
if all([isinstance(val, tuple) for val in base_vectors]):
......@@ -71,13 +93,7 @@ def igrid(*base_vectors, **kwargs):
obj = np.mgrid[base_vectors]
obj = obj.reshape(obj.shape[0], -1).T
else:
base_vectors = list(base_vectors)
# transform tuples to arrays with mgrid
for i in range(len(base_vectors)):
if isinstance(base_vectors[i], tuple):
base_vectors[i] = np.mgrid[slice(*base_vectors[i])]
if isinstance(base_vectors[i], list):
base_vectors[i] = np.array(base_vectors[i])
base_vectors = to_base_vectors(*base_vectors)
obj = np.empty(shape=(reduce(lambda x, y: x * y,
map(len, base_vectors)),
......@@ -95,11 +111,12 @@ def igrid(*base_vectors, **kwargs):
return i
loop_rec([base_vectors[i] for i in iter_order], len(base_vectors))
swap_indices = compare_permutations(iter_order, np.arange(len(base_vectors)))
swap_columns(obj, *swap_indices)
return obj
def swap_columns(array, *index_tuples):
"""
Args:
......
......@@ -38,6 +38,7 @@ def fields_to_scalars(fields):
def faces_to_maps(faces, *fields):
return [tfields.TensorFields(faces, *fields, dtype=int, dim=3)]
def maps_to_faces(maps):
if len(maps) == 0:
return np.array([])
......@@ -133,9 +134,13 @@ class Mesh3D(tfields.TensorMaps):
raise ValueError("Face dimension should be 3")
return obj
@classmethod
def plane(cls, *base_vectors, **kwargs):
vertices = tfields.Tensors.grid(*base_vectors)
base_vectors = tfields.grid.ensure_complex(*base_vectors)
base_vectors = tfields.grid.to_base_vectors(*base_vectors)
fix_coord = None
for coord in range(3):
if len(base_vectors[coord]) > 1:
......@@ -143,6 +148,8 @@ class Mesh3D(tfields.TensorMaps):
if len(base_vectors[coord]) == 0:
continue
fix_coord = coord
if fix_coord is None:
raise ValueError("Describe a plane with one variable fiexed")
var_coords = list(range(3))
var_coords.pop(var_coords.index(fix_coord))
......@@ -165,6 +172,9 @@ class Mesh3D(tfields.TensorMaps):
if not len(base_vectors) == 3:
raise AttributeError("3 base_vectors vectors required")
base_vectors = tfields.grid.ensure_complex(*base_vectors)
base_vectors = tfields.grid.to_base_vectors(*base_vectors)
indices = [0, -1]
coords = range(3)
baseLengthsAbove1 = [len(b) > 1 for b in base_vectors]
......@@ -226,57 +236,6 @@ class Mesh3D(tfields.TensorMaps):
def nfaces(self):
return self.faces.shape[0]
def toOneSegment(self, mirrorZ=True):
"""
Map the points to the first segment and mirror to positive z
if mirrorZOption is True. Special w7x method
Examples:
Build a mesh in three segments
>>> scalars = np.array([[1,2], [3,4], [5,6], [7,8]], dtype=float)
>>> m = tfields.Mesh3D([[6,-1,1], [6,0,1], [6,1,1],
... [6,-1,0.5], [6,0,0.5], [6,1,0.5]],
... [[0, 3, 4], [0, 1, 4], [1,4,5], [1,2,5]],
... faceScalars=scalars)
>>> c = m.copy()
>>> c.coordinateTransform(c.CYLINDER)
>>> c[:, 1] *= -1
>>> c[:, 2] *= -1
>>> m = tfields.Mesh3D([m, c])
>>> m2 = m.copy()
>>> m2.toSegment(1)
>>> m3 = m.copy()
>>> m3.toSegment(2)
>>> mAll = tfields.Mesh3D([m, m2, m3])
>>> mAll.toOneSegment()
>>> bool((mAll.faceScalars == scalars * 6).all())
True
"""
with self.tempCoordSys(self.CYLINDER):
dropCut = (-2 * np.pi / 10 < y) & (y < 2 * np.pi / 10)
if mirrorZ:
dropCut = dropCut & (z > 0)
dropCutMask = self.evalf(dropCut)
faceKeepMask = self.getFaceMask(dropCutMask)
excludedMesh = self.copy()
self.removeFaces(~faceKeepMask)
excludedMesh.removeFaces(faceKeepMask)
# remove 0 faces to be faster
from sympy.abc import s
zeroMask = tfields.evalf(excludedMesh.faceScalars,
expression=(s == 0),
coords=[s] * excludedMesh.getScalarDepth())
excludedMesh.removeFaces(zeroMask)
# to one segment
super(Mesh3D, self).toOneSegment(mirrorZ=mirrorZ)
centroidSF = tfields.ScalarField3D(excludedMesh.getCentroids(),
excludedMesh.faceScalars)
centroidSF.toOneSegment(mirrorZ=mirrorZ)
self.stackScalars(centroidSF)
def in_faces(self, points, delta, assign_multiple=False):
"""
Check whether points lie within triangles with Barycentric Technique
......@@ -385,16 +344,6 @@ class Mesh3D(tfields.TensorMaps):
return ~faceDeleteMask
def getScalarMap(self, mask):
"""
Return copy of self without vertices where mask is True
Copy because self is immutable
"""
# built instance that only contains the vaild points
faceKeepMask = self.getFaceMask(mask)
scalarMap = np.arange(self.nfaces())[faceKeepMask]
return scalarMap
def removeFaces(self, faceDeleteMask):
"""
Remove faces where faceDeleteMask is True
......@@ -761,18 +710,18 @@ class Mesh3D(tfields.TensorMaps):
def _intersect(self, triangle, plane, facePointsRejected):
nTrue = facePointsRejected.count(True)
lonelyBool = True if nTrue == 1 else False
index = facePointsRejected.index(lonelyBool)
lonely_bool = True if nTrue == 1 else False
index = facePointsRejected.index(lonely_bool)
s0, d0 = self._getSegmentPlaneIntersection(triangle[0], triangle[1], plane)
s1, d1 = self._getSegmentPlaneIntersection(triangle[1], triangle[2], plane)
s2, d2 = self._getSegmentPlaneIntersection(triangle[2], triangle[0], plane)
singlePoint = triangle[index]
couplePoints = [triangle[j] for j in range(3)
if not facePointsRejected[j] == lonelyBool]
if not facePointsRejected[j] == lonely_bool]
# TODO handle special cases. For now triangles with at least two points on plane are excluded
newPoints = None
new_points = None
faces = None
if len(s0) == 2:
......@@ -784,16 +733,16 @@ class Mesh3D(tfields.TensorMaps):
if len(s2) == 2:
# both points on plane
return None, None
if lonelyBool:
if lonely_bool:
# one new triangle
if len(s0) == 1 and len(s1) == 1:
newPoints = np.array([couplePoints[0], couplePoints[1], s0[0], s1[0]])
new_points = np.array([couplePoints[0], couplePoints[1], s0[0], s1[0]])
faces = np.array([[0, 2, 1], [1, 2, 3]])
elif len(s1) == 1 and len(s2) == 1:
newPoints = np.array([couplePoints[0], couplePoints[1], s1[0], s2[0]])
new_points = np.array([couplePoints[0], couplePoints[1], s1[0], s2[0]])
faces = np.array([[0, 1, 2], [0, 2, 3]])
elif len(s0) == 1 and len(s2) == 1:
newPoints = np.array([couplePoints[0], couplePoints[1], s0[0], s2[0]])
new_points = np.array([couplePoints[0], couplePoints[1], s0[0], s2[0]])
faces = np.array([[0, 1, 2], [1, 3, 2]])
else:
return None, None
......@@ -801,23 +750,23 @@ class Mesh3D(tfields.TensorMaps):
else:
# two new triangles
if len(s0) == 1 and len(s1) == 1:
newPoints = np.array([singlePoint, s0[0], s1[0]])
new_points = np.array([singlePoint, s0[0], s1[0]])
faces = np.array([[0, 2, 1]])
elif len(s1) == 1 and len(s2) == 1:
newPoints = np.array([singlePoint, s1[0], s2[0]])
new_points = np.array([singlePoint, s1[0], s2[0]])
faces = np.array([[0, 2, 1]])
elif len(s0) == 1 and len(s2) == 1:
newPoints = np.array([singlePoint, s0[0], s2[0]])
new_points = np.array([singlePoint, s0[0], s2[0]])
faces = np.array([[0, 1, 2]])
else:
return None, None
return newPoints, faces
return new_points, faces
def _cut_sympy(self, expression, at_intersection="remove"):
def _cut_sympy(self, expression, at_intersection="remove", _template=None):
"""
Partition the mesh with the cuts given
Args:
expression (sympy logical expression): see Tensors.cut
expression (sympy logical expression): see Tensors.cut
at_intersection (str): option switch for expression (see
tfields.Mesh3D.cut)
Examples:
......@@ -830,20 +779,20 @@ class Mesh3D(tfields.TensorMaps):
>>> mCut.faceScalars
array([[ 6., 7.]])
Cut with condition will return the manual/instruction on how to
conduct the cut fast (mapMesh)
Cut with 'return_template=True' will return an instruction
on how to conduct the cut fast (template)
>>> from sympy.abc import y
>>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,2,0], [1,2,0]],
... faces=[[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
... faceScalars=[[1],[2],[3],[4]])
>>> mCut, mapMesh = m.cut(y < 1.5, at_intersection='split',
>>> mCut, template = m.cut(y < 1.5, at_intersection='split',
... return_template=True)
>>> mCut.faceScalars.T
array([[ 1., 2., 3., 3., 4.]])
Applying the mapMesh results in the same result as applying the cut
Applying the template results in the same result as applying the cut
with the condition
>>> mCut2 = m.cut(mapMesh)
>>> mCut2 = m.cut(template)
>>> bool((mCut == mCut2).all())
True
>>> bool((mCut.faceScalars == mCut2.faceScalars).all())
......@@ -859,12 +808,13 @@ class Mesh3D(tfields.TensorMaps):
mask = self.evalf(expression)
# remove the points
inst = self[mask].copy()
scalarMap = self.getScalarMap(mask)
faceKeepMask = self.getFaceMask(mask)
scalarMap = np.arange(self.nfaces())[faceKeepMask]
scalarMap = scalarMap.reshape((-1, 1)).astype(float)
if not any(~mask):
# no faces have been removed.
return inst, Mesh3D([inst], faceScalars=scalarMap)
return inst, Mesh3D(inst, faceScalars=scalarMap)
if all(~mask):
# all faces have been removed)
......@@ -872,9 +822,9 @@ class Mesh3D(tfields.TensorMaps):
# add points and faces intersecting with the plane
if at_intersection == 'split' or at_intersection == 'splitRough':
expressionParts = symTools.getExpressionParts(expression)
if len(expressionParts) > 1:
onEdgeMesh = self.copy()
expression_parts = symTools.getExpressionParts(expression) # TODO: move to tfields lib
if len(expression_parts) > 1:
new_mesh = self.copy()
if at_intersection == 'splitRough':
"""
the following is, to speed up the process. Problem is, that
......@@ -886,40 +836,39 @@ class Mesh3D(tfields.TensorMaps):
faceIntersMask = np.full((self.faces.shape[0]), False, dtype=bool)
for i, face in enumerate(self.faces):
facePointsRejected = [-mask[f] for f in face]
faceOnEdge = any(facePointsRejected) and not all(facePointsRejected)
if faceOnEdge:
face_on_edge = any(facePointsRejected) and not all(facePointsRejected)
if face_on_edge:
faceIntersMask[i] = True
onEdgeMesh.removeFaces(-faceIntersMask)
for exprPart in expressionParts:
onEdgeMesh = onEdgeMesh.cut(exprPart, at_intersection='split')
newMesh = onEdgeMesh
elif len(expressionParts) == 1:
new_mesh.removeFaces(-faceIntersMask)
for exprPart in expression_parts:
new_mesh, _template = new_mesh._cut_sympy(exprPart,
at_intersection='split',
_template=_template)
elif len(expression_parts) == 1:
points = [sympy.symbols('x0, y0, z0'),
sympy.symbols('x1, y1, z1'),
sympy.symbols('x2, y2, z2')]
fpr = sympy.symbols('fpr2, fpr1, fpr2')
planeSympy = symTools.getPlane(expression)
n = np.array(planeSympy.normal_vector).astype(float)
d = -n.dot(np.array(planeSympy.p1).astype(float))
plane = {'normal': n, 'd': d}
plane_sympy = symTools.getPlane(expression)
norm_sympy = np.array(plane_sympy.normal_vector).astype(float)
d = -norm_sympy.dot(np.array(plane_sympy.p1).astype(float))
plane = {'normal': norm_sympy, 'd': d}
norm_vectors = self.triangles.norms()
newPoints = np.empty((0, 3))
newFaces = np.empty((0, 3))
newFaceScalars = []
newNormVectors = []
new_points = np.empty((0, 3))
new_faces = np.empty((0, 3))
new_map_fields = [[] for field in self.maps[0].fields]
new_norm_vectors = []
newScalarMap = []
nNew = 0
n_new = 0
for i, face in enumerate(self.faces):
"""
facePointsRejected is a mask for each face that is True, where
a Point is on the rejected side of the plane
"""
facePointsRejected = [~mask[f] for f in face]
faceOnEdge = any(facePointsRejected) and not all(facePointsRejected)
if faceOnEdge:
face_on_edge = any(facePointsRejected) and not all(facePointsRejected)
if face_on_edge:
"""
define the two lines that are intersecting the plane:
one point will be alone on one side of the cutting plane (singlePoint)
......@@ -927,25 +876,27 @@ class Mesh3D(tfields.TensorMaps):
"""
nTrue = facePointsRejected.count(True)
# the bool that occures on
lonelyBool = True if nTrue == 1 else False
lonely_bool = True if nTrue == 1 else False
triangle_points = [self[f] for f in face]
"""
Add the intersection points and faces
"""
# tick = time()
if lonelyBool:
if lonely_bool:
# singlePoint on cut away side
newP, newF = self._intersect(triangle_points, plane, facePointsRejected)
if newP is not None:
newP = np.array(newP)
np.place(newP, np.logical_and(0 < newP, newP <= eps), 0.0)
newPoints = np.concatenate((newPoints, newP))
newFaces = np.concatenate((newFaces, newF + nNew))
newFaceScalars.extend([self.faceScalars[i]] * 2)
newNormVectors.extend([norm_vectors[i]] * 2)
new_points = np.concatenate((new_points, newP))
new_faces = np.concatenate((new_faces, newF +
n_new))
for fld_idx, fld in enumerate(self.maps[0].fields):
new_map_fields[fld_idx].extend([fld[i]] * 2)
new_norm_vectors.extend([norm_vectors[i]] * 2)
newScalarMap.extend([i] * 2)
nNew += 4
n_new += 4
else:
pass
else:
......@@ -954,43 +905,97 @@ class Mesh3D(tfields.TensorMaps):
if newP is not None:
newP = np.array(newP)
np.place(newP, np.logical_and(0 < newP, newP <= eps), 0.0)
newPoints = np.concatenate((newPoints, newP))
newFaces = np.concatenate((newFaces, newF + nNew))
newFaceScalars.extend([self.faceScalars[i]] * 1)
newNormVectors.extend([norm_vectors[i]] * 1)
new_points = np.concatenate((new_points, newP))
new_faces = np.concatenate((new_faces, newF +
n_new))
for fld_idx, fld in enumerate(self.maps[0].fields):
new_map_fields[fld_idx].extend([fld[i]] * 1)
new_norm_vectors.extend([norm_vectors[i]] * 1)
newScalarMap.extend([i])
nNew += 3
n_new += 3
else:
pass
newPoints = [[float(p[0]), float(p[1]), float(p[2])] for p in newPoints]
newMesh = self.__class__(newPoints, faces=newFaces,
coordSys=inst.getCoordSys(),
faceScalars=newFaceScalars)
newMesh.align_norms(newNormVectors)
new_points = [[float(p[0]), float(p[1]), float(p[2])] for p in new_points] # TODO:[map(float, p) for p in new_points]
face_map = tfields.TensorFields(new_faces, *new_map_fields, dtype=int,
coordSys=inst.coordSys)
new_mesh = self.__class__(new_points, maps=[face_map],
coordSys=inst.coordSys,
faceScalars=new_map_fields)
new_mesh.align_norms(new_norm_vectors)
newScalarMap = np.array(newScalarMap).reshape((-1, 1)).astype(float)
scalarMap = np.concatenate((scalarMap, newScalarMap))
inst = self.__class__([inst, newMesh])
else:
raise ValueError("CutExpression is not splitable.")
raise ValueError("Sympy expression is not splitable.")
'''
Merge parts of mesh that are clearly in the cuts (inst) with
parts of mesh on the edge of the cuts, where new vertices and faces
have been defined
'''
inst = self.__class__.merged(inst, new_mesh)
elif at_intersection == 'remove':
pass
else:
raise AttributeError("No at_intersection method called {at_intersection} "
"implemented".format(**locals()))
meshMap = inst.copy()
meshMap.faceScalars = scalarMap
return inst, meshMap
template = inst.copy()
template.faceScalars = scalarMap
return inst, template
def _cut_template(self, template):
if not template:
return Mesh3D([])
inst = Mesh3D(template,
*[field[template.fields[0].astype(int)]
for field in self.fileds],
maps=template.maps)
"""
Args:
template (tfields.Mesh3D)
Examples:
>>> import tfields
>>> import numpy as np
Mesh
>>> mmap = tfields.TensorFields([[0, 1, 2], [0, 3, 4]],
... [[42, 21], [-42, -21]])
>>> m = tfields.Mesh3D([[0]*3, [1]*3, [2]*3, [3]*3, [4]*3],
... fields=[[0.0, 0.1, 0.2, 0.3, 0.4],
... [0.0, -0.1, -0.2, -0.3, -0.4]],
... maps=[mmap])
Template
>>> tmap = tfields.TensorFields([[0, 3, 4], [0, 1, 2]],
... [1, 0])
>>> t = tfields.Mesh3D([[0]*3, [-1]*3, [-2]*3, [-3]*3, [-4]*3],
... fields=[[1, 0, 3, 2, 4]],
... maps=[tmap])
>>> res = m._cut_template(t)
>>> assert np.array_equal(res.fields,
... [[0.1, 0.0, 0.3, 0.2, 0.4],
... [-0.1, 0.0, -0.3, -0.2, -0.4]])
>>> assert np.array_equal(res.maps[0].fields[0],
... [[-42, -21], [42, 21]])
"""
if template.fields:
fields = [field[template.fields[0].astype(int)]
for field in self.fields]
else:
fields = []
maps = []
for mp, template_mp in zip(self.maps, template.maps):
if template_mp.fields:
mp_fields = [field[template_mp.fields[0].astype(int)]
for field in mp.fields]
else:
mp_fields = []
new_mp = tfields.TensorFields(template_mp,
*mp_fields)
maps.append(new_mp)
inst = tfields.Mesh3D(template,
fields=fields,
maps=maps)
return inst
def cut(self, expression, coordSys=None, at_intersection=None,
......@@ -998,7 +1003,16 @@ class Mesh3D(tfields.TensorMaps):
"""
cut method for Mesh3D.
Args:
expression (sympy_logical expression | Mesh3D):
expression (sympy logical expression | Mesh3D):
sympy locical expression: Sympy expression that defines planes
in 3D
Mesh3D: A mesh3D will be interpreted as a template, i.e. a
fast instruction of how to cut the triangles.
It is the second part of the tuple, returned by a previous
cut with a sympy locial expression with 'return_template=True'.
We use the vertices and maps of the Mesh as the sceleton of
the returned mesh. The fields are mapped according to
indices in the template.maps[i].fields.
coordSys (coordinate system to cut in):
at_intersection (str): instruction on what to do, when a cut will intersect a triangle.
Options: "remove" (Default)
......@@ -1007,6 +1021,7 @@ class Mesh3D(tfields.TensorMaps):
to redo the same cut fast
Examples: