Commit 6c0b3246 authored by dboe's avatar dboe
Browse files

debugging split

parent 0e905d7f
Pipeline #83961 passed with stages
in 56 seconds
......@@ -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()
......@@ -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:
......
Markdown is supported
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