diff --git a/tests/test_mesh3D.py b/tests/test_mesh3D.py
index bc02bc6f3a16f16579df2635d5eac05572e2db8b..8742b4ed0ea36e6fa7ee7d1103041710b5cb2f95 100644
--- a/tests/test_mesh3D.py
+++ b/tests/test_mesh3D.py
@@ -2,7 +2,7 @@ from tempfile import NamedTemporaryFile
import tfields
import numpy as np
import unittest
-import sympy # NOQA: F401
+import sympy
import os
import sys
from tests.test_core import TensorMaps_Check
@@ -58,6 +58,56 @@ class Square_Dense_Test(Mesh3D_Check, unittest.TestCase):
self.assertEqual(len(cut_inst), 10)
+class Box_Sparse_Test(Mesh3D_Check, unittest.TestCase):
+ def setUp(self):
+ self._inst = tfields.Mesh3D.grid((0, 1, 2), (1, 2, 2), (2, 3, 2))
+
+ def test_cut_split(self):
+ from sympy.abc import x
+
+ cuts = self._inst.cut(x < 0.2, at_intersection="split")
+ print(cuts)
+ self.assertTrue(
+ cuts.equal(
+ [
+ [0.0, 1.0, 2.0],
+ [0.0, 1.0, 3.0],
+ [0.0, 2.0, 2.0],
+ [0.0, 2.0, 3.0],
+ [0.2, 1.0, 2.8],
+ [0.2, 1.0, 2.0],
+ [0.2, 1.0, 3.0],
+ [0.2, 1.8, 2.0],
+ [0.2, 2.0, 2.0],
+ [0.2, 2.0, 2.8],
+ [0.2, 2.0, 3.0],
+ [0.2, 1.8, 3.0],
+ ],
+ rtol=1e-12,
+ )
+ )
+ self.assertTrue(
+ cuts.maps[3].equal(
+ [
+ [0, 1, 2],
+ [1, 2, 3],
+ [0, 1, 4],
+ [0, 4, 5],
+ [1, 4, 6],
+ [0, 2, 7],
+ [0, 7, 5],
+ [2, 7, 8],
+ [2, 3, 9],
+ [2, 9, 8],
+ [3, 9, 10],
+ [1, 3, 11],
+ [1, 11, 6],
+ [3, 11, 10],
+ ]
+ )
+ )
+
+
class Sphere_Test(Mesh3D_Check, unittest.TestCase):
def setUp(self):
basis_points = 4
@@ -81,12 +131,3 @@ class IO_Stl_test(unittest.TestCase): # no Mesh3D_Check for speed
if __name__ == "__main__":
unittest.main()
- # TODO!
- import rna
- from sympy.abc import x
-
- kwargs = dict(dim=3, edgecolor="k")
- mesh = tfields.Mesh3D.grid((0, 1, 2), (1, 2, 2), (2, 3, 2))
- cuts = mesh.cut(x < 0.2, at_intersection="split")
- cuts.plot(**kwargs)
- rna.plotting.show()
diff --git a/tfields/mesh3D.py b/tfields/mesh3D.py
index 603be92afb0ec0837d18ff11692152168b03098e..2ed7f51ebc80e1253cc517179c03a54d5f1b1ec4 100644
--- a/tfields/mesh3D.py
+++ b/tfields/mesh3D.py
@@ -7,7 +7,6 @@ Mail: daniel.boeckenhoff@ipp.mpg.de
part of tfields library
"""
import numpy as np
-import sympy
import rna
import tfields
@@ -95,43 +94,43 @@ def _intersect(triangle, plane, vertices_rejected):
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
+ new_vertices = None
if len(s0) == 2:
# both points on plane
- return new_points
+ return new_vertices
if len(s1) == 2:
# both points on plane
- return new_points
+ return new_vertices
if len(s2) == 2:
# both points on plane
- return new_points
+ return new_vertices
if lonely_bool:
# two new triangles
if len(s0) == 1 and len(s1) == 1:
- new_points = [
+ new_vertices = [
[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 = [
+ new_vertices = [
[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 = [
+ new_vertices = [
[couple_indices[0], couple_indices[1], s0[0]],
[couple_indices[1], s2[0], complex(2)],
]
else:
# one new triangle
if len(s0) == 1 and len(s1) == 1:
- new_points = [[single_index, s1[0], s0[0]]]
+ new_vertices = [[single_index, s1[0], s0[0]]]
elif len(s1) == 1 and len(s2) == 1:
- new_points = [[single_index, s2[0], s1[0]]]
+ new_vertices = [[single_index, s2[0], s1[0]]]
elif len(s0) == 1 and len(s2) == 1:
- new_points = [[single_index, s0[0], s2[0]]]
- return new_points
+ new_vertices = [[single_index, s0[0], s2[0]]]
+ return new_vertices
class Mesh3D(tfields.TensorMaps):
@@ -801,9 +800,9 @@ class Mesh3D(tfields.TensorMaps):
expression_parts = tfields.lib.symbolics.split_expression(expression)
if len(expression_parts) > 1:
new_mesh = inst.copy()
- for exprPart in expression_parts:
+ for expr_part in expression_parts:
inst, _ = inst._cut_sympy(
- exprPart, at_intersection=at_intersection, _in_recursion=True
+ expr_part, at_intersection=at_intersection, _in_recursion=True
)
elif len(expression_parts) == 1:
face_delete_indices = set([])
@@ -849,23 +848,20 @@ class Mesh3D(tfields.TensorMaps):
face_inters_mask[i] = True
new_mesh.remove_faces(-face_inters_mask)
- for exprPart in expression_parts:
+ for expr_part in expression_parts:
inst, _ = inst._cut_sympy(
- exprPart, at_intersection="split", _in_recursion=True
+ expr_part, at_intersection="split", _in_recursion=True
)
elif len(expression_parts) == 1:
- points = [
- sympy.symbols("x0, y0, z0"),
- sympy.symbols("x1, y1, z1"),
- sympy.symbols("x2, y2, z2"),
- ]
+ # build plane from cut expression
plane_sympy = tfields.lib.symbolics.to_plane(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}
+ # initialize empty containers
norm_vectors = inst.triangles().norms()
- new_points = np.empty((0, 3))
+ new_vertices = np.empty((0, 3))
new_faces = np.empty((0, 3))
new_fields = [
tfields.Tensors(
@@ -875,18 +871,17 @@ class Mesh3D(tfields.TensorMaps):
]
new_map_fields = [[] for field in inst.maps[3].fields]
new_norm_vectors = []
- newScalarMap = []
- n_new = 0
+ # copy TODO?
vertices = np.array(inst)
faces = np.array(inst.maps[3])
fields = [np.array(field) for field in inst.fields]
faces_fields = [np.array(field) for field in inst.maps[3].fields]
- face_delete_indices = set([])
+ face_delete_indices = set([]) # indexing not intersected faces
for i, face in enumerate(inst.maps[3]):
"""
- vertices_rejected is a mask for each face that is True,
+ 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]
@@ -905,7 +900,7 @@ class Mesh3D(tfields.TensorMaps):
intersection = _intersect(
triangle_points, plane, vertices_rejected
)
- last_idx = len(vertices) - 1
+ last_idx = len(vertices)
for tri_list in intersection:
new_face = []
for item in tri_list: