diff --git a/tfields/mesh3D.py b/tfields/mesh3D.py
index a9a6e5bd52fc90a6b19dd403aa117088bdfc5545..2849fa477e398ac7fa44fcf9310a554919af08e2 100644
--- a/tfields/mesh3D.py
+++ b/tfields/mesh3D.py
@@ -38,37 +38,37 @@ def _segment_plane_intersection(p0, p1, plane):
distance1 = _dist_from_plane(p1, plane)
p0OnPlane = abs(distance0) < np.finfo(float).eps
p1OnPlane = abs(distance1) < np.finfo(float).eps
- outPoints = []
+ points = []
direction = 0
if p0OnPlane:
- outPoints.append(p0)
+ points.append(p0)
if p1OnPlane:
- outPoints.append(p1)
+ points.append(p1)
# remove duplicate points
- if len(outPoints) > 1:
- outPoints = tfields.Points3D(np.unique(outPoints, axis=0))
+ if len(points) > 1:
+ points = np.unique(points, axis=0)
if p0OnPlane and p1OnPlane:
- return outPoints, direction
+ return points, direction
if distance0 * distance1 > np.finfo(float).eps:
- return outPoints, direction
+ return points, direction
direction = np.sign(distance0)
if abs(distance0) < np.finfo(float).eps:
- return outPoints, direction
+ return points, direction
elif abs(distance1) < np.finfo(float).eps:
- return outPoints, direction
+ return points, direction
if abs(distance0 - distance1) > np.finfo(float).eps:
t = distance0 / (distance0 - distance1)
else:
- return outPoints, direction
+ return points, direction
- outPoints.append(p0 + t * (p1 - p0))
+ points.append(p0 + t * (p1 - p0))
# remove duplicate points
- if len(outPoints) > 1:
- outPoints = tfields.Points3D(np.unique(outPoints, axis=0))
- return outPoints, direction
+ if len(points) > 1:
+ points = np.unique(points, axis=0)
+ return points, direction
def _intersect(triangle, plane, vertices_rejected):
@@ -76,7 +76,15 @@ def _intersect(triangle, plane, vertices_rejected):
Intersect a triangle with a plane. Give the info, which side of the
triangle is rejected by passing the mask vertices_rejected
Returns:
- tuple: points, faces
+ list of list. The inner list is of length 3 and refers to the points of
+ new triangles. The reference is done with varying types:
+ int: reference to triangle index
+ complex: reference to duplicate point. This only happens in case
+ two triangles are returned. Then only in the second triangle
+ iterable: new vertex
+
+ TODO:
+ align norm vectors with previous face
"""
nTrue = vertices_rejected.count(True)
lonely_bool = True if nTrue == 1 else False
@@ -85,56 +93,43 @@ def _intersect(triangle, plane, vertices_rejected):
s1, d1 = _segment_plane_intersection(triangle[1], triangle[2], plane)
s2, d2 = _segment_plane_intersection(triangle[2], triangle[0], plane)
- singlePoint = triangle[index]
- couplePoints = [triangle[j] for j in range(3)
- if not vertices_rejected[j] == lonely_bool]
+ single_index = index
+ couple_indices = [j for j in range(3)
+ if not vertices_rejected[j] == lonely_bool]
# TODO handle special cases. For now triangles with at least two points on plane are excluded
new_points = None
- faces = None
if len(s0) == 2:
# both points on plane
- return None, None
+ return new_points
if len(s1) == 2:
# both points on plane
- return None, None
+ return new_points
if len(s2) == 2:
# both points on plane
- return None, None
+ return new_points
if lonely_bool:
- # one new triangle
+ # two new triangles
if len(s0) == 1 and len(s1) == 1:
- new_points = np.array([couplePoints[0], couplePoints[1], s0[0], s1[0]])
- faces = np.array([[0, 2, 1], [1, 2, 3]])
+ new_points = [[couple_indices[0], s0[0], couple_indices[1]],
+ [couple_indices[1], complex(1), s1[0]]]
elif len(s1) == 1 and len(s2) == 1:
- new_points = np.array([couplePoints[0], couplePoints[1], s1[0], s2[0]])
- faces = np.array([[0, 1, 2], [0, 2, 3]])
+ new_points = [[couple_indices[0], couple_indices[1], s1[0]],
+ [couple_indices[0], complex(2), s2[0]]]
elif len(s0) == 1 and len(s2) == 1:
- 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
+ new_points = [[couple_indices[0], couple_indices[1], s0[0]],
+ [couple_indices[1], s2[0], complex(2)]]
else:
- # two new triangles
+ # one new triangle
if len(s0) == 1 and len(s1) == 1:
- new_points = np.array([singlePoint, s0[0], s1[0]])
- faces = np.array([[0, 2, 1]])
+ new_points = [[single_index, s1[0], s0[0]]]
elif len(s1) == 1 and len(s2) == 1:
- new_points = np.array([singlePoint, s1[0], s2[0]])
- faces = np.array([[0, 2, 1]])
+ new_points = [[single_index, s2[0], s1[0]]]
elif len(s0) == 1 and len(s2) == 1:
- new_points = np.array([singlePoint, s0[0], s2[0]])
- faces = np.array([[0, 1, 2]])
- else:
- return None, None
- return new_points, faces
-
-
-
-
-
+ new_points = [[single_index, s0[0], s2[0]]]
+ return new_points
def scalars_to_fields(scalars):
@@ -790,11 +785,9 @@ class Mesh3D(tfields.TensorMaps):
# mask for points that do not fulfill the cut expression
mask = inst.evalf(expression)
# remove the points
- inst_pre_mask = inst.copy()
- inst = inst[mask]
if not any(~mask) or all(~mask):
- pass
+ inst = inst[mask]
elif at_intersection == 'split' or at_intersection == 'splitRough':
# add points and faces intersecting with the plane
expression_parts = tfields.lib.symbolics.split_expression(expression)
@@ -804,7 +797,7 @@ class Mesh3D(tfields.TensorMaps):
cuts
'''
if len(expression_parts) > 1:
- new_mesh = inst_pre_mask.copy()
+ new_mesh = inst.copy()
if at_intersection == 'splitRough':
"""
the following is, to speed up the process. Problem is, that
@@ -813,8 +806,8 @@ class Mesh3D(tfields.TensorMaps):
still overlaps with the cut.
These are at the intersection line between two cuts.
"""
- faceIntersMask = np.full((self.faces.shape[0]), False, dtype=bool)
- for i, face in enumerate(self.faces):
+ faceIntersMask = np.full((inst.faces.shape[0]), False, dtype=bool)
+ for i, face in enumerate(inst.faces):
vertices_rejected = [-mask[f] for f in face]
face_on_edge = any(vertices_rejected) and not all(vertices_rejected)
if face_on_edge:
@@ -822,11 +815,11 @@ class Mesh3D(tfields.TensorMaps):
new_mesh.removeFaces(-faceIntersMask)
for exprPart in expression_parts:
- new_mesh, _ = new_mesh._cut_sympy(exprPart,
- at_intersection='split',
- _in_recursion=True)
+ inst, _ = inst._cut_sympy(exprPart,
+ at_intersection='split',
+ _in_recursion=True)
elif len(expression_parts) == 1:
- # TODO maps[0] -> smthng like self.get_map(dim=3)
+ # TODO maps[0] -> smthng like inst.get_map(dim=3)
points = [sympy.symbols('x0, y0, z0'),
sympy.symbols('x1, y1, z1'),
sympy.symbols('x2, y2, z2')]
@@ -835,7 +828,7 @@ class Mesh3D(tfields.TensorMaps):
d = -norm_sympy.dot(np.array(plane_sympy.p1).astype(float))
plane = {'normal': norm_sympy, 'd': d}
- norm_vectors = self.triangles.norms()
+ norm_vectors = inst.triangles.norms()
new_points = np.empty((0, 3))
new_faces = np.empty((0, 3))
new_fields = [tfields.Tensors(np.empty((0,) + field.shape[1:]),
@@ -846,94 +839,80 @@ class Mesh3D(tfields.TensorMaps):
newScalarMap = []
n_new = 0
- for i, face in enumerate(self.maps[0]):
+ vertices = np.array(inst)
+ faces = np.array(inst.maps[0])
+ fields = [np.array(field) for field in inst.fields]
+ faces_fields = [np.array(field) for field in inst.maps[0].fields]
+
+ face_delete_indices = set([])
+ for i, face in enumerate(inst.maps[0]):
"""
vertices_rejected is a mask for each face that is True, where
a Point is on the rejected side of the plane
"""
vertices_rejected = [~mask[f] for f in face]
- face_on_edge = any(vertices_rejected) and not all(vertices_rejected)
- 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)
- the other two will be on ther other side (couplePoints)
- """
+ if any(vertices_rejected):
+ # delete face
+ face_delete_indices.add(i)
+ if any(vertices_rejected) and not all(vertices_rejected):
+ # face on edge
nTrue = vertices_rejected.count(True)
- # the bool that occures on
lonely_bool = True if nTrue == 1 else False
- triangle_points = [self[f] for f in face]
+ triangle_points = [inst[f] for f in face]
"""
Add the intersection points and faces
"""
print "______________"
- print triangle_points, vertices_rejected
- newP, newF = _intersect(triangle_points, plane, vertices_rejected)
- print newP, newF
- if newP is not None:
- if lonely_bool:
- '''
- singlePoint on cut away side
- build two new triangles with two known and two new
- points
- '''
- newP = np.array(newP)
- np.place(newP, np.logical_and(0 < newP, newP <= eps), 0.0)
- new_points = np.concatenate((new_points, newP))
- new_fields = [tfields.Tensors.merged(field,
- np.full((4,) + field.shape[1:], np.nan))
- for field in new_fields]
- new_faces = np.concatenate((new_faces, newF +
- n_new))
- for field_idx, field in enumerate(self.maps[0].fields):
- new_map_fields[field_idx].extend([field[i]] * 2)
- if not _in_recursion:
- new_map_fields[-1].extend([i] * 2)
- new_norm_vectors.extend([norm_vectors[i]] * 2)
- newScalarMap.extend([i] * 2)
- n_new += 4
- else:
- '''
- singlePoint on keep side
- build one new triangle with one known and two new points
- '''
- newP = np.array(newP)
- np.place(newP, np.logical_and(0 < newP, newP <= eps), 0.0)
- new_points = np.concatenate((new_points, newP))
- new_fields = [tfields.Tensors.merged(field,
- np.full((3,) + field.shape[1:], np.nan))
- for field in new_fields]
- new_faces = np.concatenate((new_faces, newF +
- n_new))
- for field_idx, field in enumerate(self.maps[0].fields):
- new_map_fields[field_idx].extend([field[i]] * 1)
- if not _in_recursion:
- new_map_fields[-1].extend([i] * 1)
- new_norm_vectors.extend([norm_vectors[i]] * 1)
- newScalarMap.extend([i])
- n_new += 3
-
- new_points = [map(float, p) for p in new_points]
- face_map = tfields.TensorFields(new_faces, *new_map_fields,
+ newP = _intersect(triangle_points, plane, vertices_rejected)
+ last_idx = len(vertices) - 1
+ for tri_list in newP:
+ print vertices
+ face = []
+ for item in tri_list:
+ if isinstance(item, int):
+ # reference to old vertex
+ face.append(item)
+ elif isinstance(item, complex):
+ # reference to new vertex that has been
+ # concatenated already
+ face.append(last_idx + int(item.imag))
+ else:
+ # new vertex
+ face.append(len(vertices))
+ vertices = np.append(vertices,
+ [map(float, item)],
+ axis=0)
+ fields = [np.append(field,
+ np.full((1,) + field.shape[1:], np.nan),
+ axis=0)
+ for field in fields]
+ fields[-1][-1] = len(vertices)
+ print vertices
+ print newP
+ print face
+ faces = np.append(faces, [face], axis=0)
+ faces_fields = [np.append(field,
+ np.full((1,) + field.shape[1:], np.nan),
+ axis=0)
+ for field in faces_fields]
+ faces_fields[-1][-1] = i
+
+ face_map = tfields.TensorFields(faces, *faces_fields,
dtype=int,
- coordSys=inst.coordSys)
- new_mesh = tfields.Mesh3D(new_points,
- *new_fields,
- maps=[face_map],
- coordSys=inst.coordSys)
- new_mesh.align_norms(new_norm_vectors)
+ coordSys=inst.maps[0].coordSys)
+ inst = tfields.Mesh3D(vertices,
+ *fields,
+ maps=[face_map] + inst.maps[1:],
+ coordSys=inst.coordSys)
+ mask = np.full(len(inst.maps[0]), True, dtype=bool)
+ for face_idx in range(len(inst.maps[0])):
+ if face_idx in face_delete_indices:
+ mask[face_idx] = False
+ inst.maps[0] = inst.maps[0][mask]
else:
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 = tfields.Mesh3D.merged(inst, new_mesh)
-
elif at_intersection == 'remove':
pass
else:
@@ -1132,6 +1111,12 @@ class Mesh3D(tfields.TensorMaps):
self.faces[mask, 1] = self.faces[mask, 2]
self.faces[mask, 2] = temp
+ def plot(self):
+ import mplTools as mpt
+ mpt.plotMesh(self, self.faces, color=self.maps[0].fields[0], vmin=0,
+ vmax=20, axis=mpt.gca(3))
+ mpt.plt.show()
+
if __name__ == '__main__':
import doctest