Commit 5518e2cd authored by Daniel Boeckenhoff's avatar Daniel Boeckenhoff
Browse files

grosser umbau in mesh3D cut

parent 5069715e
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment