diff --git a/test/test_mesh3D.py b/test/test_mesh3D.py index 370034d97dd8722fa6184090b8162e18d66aefb8..06aeaf22842de843fb0cbb08c767a207c1a4e68d 100644 --- a/test/test_mesh3D.py +++ b/test/test_mesh3D.py @@ -4,13 +4,28 @@ import unittest import sympy # NOQA: F401 import os import sys -from .test_core import Tensors_Check +from .test_core import TensorFields_Check THIS_DIR = os.path.dirname( os.path.realpath(os.path.join(os.getcwd(), os.path.expanduser(__file__)))) sys.path.append(os.path.normpath(os.path.join(THIS_DIR))) -class Sphere_Test(Tensors_Check, unittest.TestCase): +class Mesh3D_Check(TensorFields_Check, unittest.TestCase): + def test_cut_split(self): + x, y, z = sympy.symbols('x y z') + self._inst.cut(x + 1./100*y > 0, at_intersection='split') + + def test_triangle(self): + tri = self._inst.triangles() + mesh = tri.mesh() + # print(self._inst, self._inst.maps) + # print(tri) + # print(mesh, mesh.maps) + # print(tfields.Tensors(self._inst).equal(mesh)) + self.assertTrue(self._inst.equal(mesh)) + + +class Sphere_Test(Mesh3D_Check): def setUp(self): self._inst = tfields.Mesh3D.grid( (1, 1, 1), @@ -20,15 +35,10 @@ class Sphere_Test(Tensors_Check, unittest.TestCase): self._inst.transform('cartesian') self._inst[:, 1] += 2 - def test_cut_split(self): - x, y, z = sympy.symbols('x y z') - self._inst.cut(x + 1./100*y > 0, at_intersection='split') - def test_triangle(self): - tri = self._inst.triangles() - mesh = tri.mesh() - print(self._inst, mesh) - self.assertTrue(self._inst.equal(mesh)) +class Square_Test(Mesh3D_Check): + def setUp(self): + self._inst = tfields.Mesh3D.plane((0,1,2j), (0,1,2j), (0,0,1j)) class IO_Stl_test(unittest.TestCase): diff --git a/tfields/core.py b/tfields/core.py index b57240f93ad9336990088118d0409451de9ee202..36889ad91083ce6a0bd44625591c9598892cf2ab 100644 --- a/tfields/core.py +++ b/tfields/core.py @@ -1406,27 +1406,40 @@ class Tensors(AbstractNdarray): return inst, template def _cut_template(self, template): - # # Redirect fields - # fields = [] - # if template.fields: - # template_field = np.array(template.fields[0]) - # if len(self) > 0: - # ''' - # if new vertices have been created in the template, it is - # in principle unclear what fields we have to refer to. - # Thus in creating the template, we gave np.nan. - # To make it fast, we replace nan with 0 as a dummy and correct - # the field entries afterwards with np.nan. - # ''' - # nan_mask = np.isnan(template_field) - # template_field[nan_mask] = 0 # dummy reference to index 0. - # template_field = template_field.astype(int) - # for field in self.fields: - # projected_field = field[template_field] - # projected_field[nan_mask] = np.nan # correction for nan - # fields.append(projected_field) - # return type(self)(Tensors(template), *fields) - return self[template.fields[0]] + """ + In principle, what we do is returning + self[template.fields[0]] + + If the templates tensors is given (has no dimension 0), 0))), we switch + to only extruding the field entries according to the indices provided + by template.fields[0]. This allows the template to define additional + points, extending the object it should cut. This becomes relevant for + Mesh3D when adding vertices at the edge of the cut is necessary. + """ + # Redirect fields + fields = [] + if template.fields: + template_field = np.array(template.fields[0]) + if len(self) > 0: + ''' + if new vertices have been created in the template, it is + in principle unclear what fields we have to refer to. + Thus in creating the template, we gave np.nan. + To make it fast, we replace nan with 0 as a dummy and correct + the field entries afterwards with np.nan. + ''' + nan_mask = np.isnan(template_field) + template_field[nan_mask] = 0 # dummy reference to index 0. + template_field = template_field.astype(int) + for field in self.fields: + projected_field = field[template_field] + projected_field[nan_mask] = np.nan # correction for nan + fields.append(projected_field) + if dim(template) == 0: + tensors = self + else: + tensors = template + return type(self)(Tensors(template), *fields) def cut(self, expression, coord_sys=None, return_template=False, **kwargs): """ @@ -2266,16 +2279,16 @@ class TensorMaps(TensorFields): >>> sliced = mesh[2:] >>> assert isinstance(sliced, tfields.TensorMaps) >>> assert isinstance(sliced.fields[0], tfields.Tensors) - >>> assert isinstance(sliced.maps[0], tfields.TensorFields) + >>> assert isinstance(sliced.maps[3], tfields.TensorFields) >>> assert sliced.fields[0].equal([10.5, 1, 1]) - >>> assert sliced.maps[0].equal([[0, 1, 2]]) - >>> assert sliced.maps[0].fields[0].equal([[5, 6]]) + >>> assert sliced.maps[3].equal([[0, 1, 2]]) + >>> assert sliced.maps[3].fields[0].equal([[5, 6]]) Picking >>> picked = mesh[1] >>> assert np.array_equal(picked, [0, 0, 1]) - >>> assert np.array_equal(picked.maps[0], np.empty((0, 3))) + >>> assert np.array_equal(picked.maps[3], np.empty((0, 3))) >>> assert np.array_equal(picked.maps[1], [[0]]) Masking @@ -2284,7 +2297,7 @@ class TensorMaps(TensorFields): >>> assert masked.equal([[0, 0, 0], [0, -1, 0], [1, 1, 1], [-1, -1, -1]]) >>> assert masked.fields[0].equal([42, 10.5, 1, 1]) >>> assert masked.fields[1].equal([1, 3, 3, 3]) - >>> assert masked.maps[0].equal([[1, 2, 3]]) + >>> assert masked.maps[3].equal([[1, 2, 3]]) >>> assert masked.maps[1].equal([[0], [1], [2], [3]]) Iteration @@ -2293,46 +2306,40 @@ class TensorMaps(TensorFields): """ item = super(TensorMaps, self).__getitem__(index) - try: - if issubclass(type(item), TensorMaps): - if isinstance(index, tuple): - index = index[0] - if len(item.maps) == 0: - item.maps = Maps(item.maps) - indices = np.array(range(len(self))) - keep_indices = indices[index] - if isinstance(keep_indices, (int, np.int64, np.int32)): - keep_indices = [keep_indices] - delete_indices = set(indices).difference(set(keep_indices)) - - # correct all maps that contain deleted indices - for map_dim in self.maps.keys(): - # build mask, where the map should be deleted - map_delete_mask = np.full( - (len(self.maps[map_dim]),), False, dtype=bool - ) - for i, mp in enumerate(self.maps[map_dim]): - for index in mp: - if index in delete_indices: - map_delete_mask[i] = True - break - map_mask = ~map_delete_mask - - # build the correction counters - move_up_counter = np.zeros( - self.maps[map_dim].shape, dtype=int - ) - for p in delete_indices: - move_up_counter[self.maps[map_dim] > p] -= 1 - - item.maps[map_dim] = ( - self.maps[map_dim] + move_up_counter - )[map_mask] - except IndexError as err: - warnings.warn( - "Index error occured for field.__getitem__. Error " - "message: {err}".format(**locals()) - ) + if issubclass(type(item), TensorMaps): + if isinstance(index, tuple): + index = index[0] + if len(item.maps) == 0: + item.maps = Maps(item.maps) + indices = np.array(range(len(self))) + keep_indices = indices[index] + if isinstance(keep_indices, (int, np.int64, np.int32)): + keep_indices = [keep_indices] + delete_indices = set(indices).difference(set(keep_indices)) + + # correct all maps that contain deleted indices + for map_dim in self.maps.keys(): + # build mask, where the map should be deleted + map_delete_mask = np.full( + (len(self.maps[map_dim]),), False, dtype=bool + ) + for i, mp in enumerate(self.maps[map_dim]): + for index in mp: + if index in delete_indices: + map_delete_mask[i] = True + break + map_mask = ~map_delete_mask + + # build the correction counters + move_up_counter = np.zeros( + self.maps[map_dim].shape, dtype=int + ) + for p in delete_indices: + move_up_counter[self.maps[map_dim] > p] -= 1 + + item.maps[map_dim] = ( + self.maps[map_dim] + move_up_counter + )[map_mask] return item @@ -2461,7 +2468,7 @@ class TensorMaps(TensorFields): index_lut = np.full(len(self), np.nan) index_lut[template.fields[0]] = np.arange(len(template.fields[0])) for mp_dim, mp in inst.maps.items(): - mp = mp.cut(template.maps[mp_dim]) + mp = mp._cut_template(template.maps[mp_dim]) if template.fields: mp = Maps.to_map(index_lut[mp], *mp.fields) inst.maps[mp_dim] = mp @@ -2603,7 +2610,7 @@ class TensorMaps(TensorFields): """ remove_condition = np.array(remove_condition) - # # built instance that only contains the vaild points + # # build instance that only contains the vaild points # inst = self[~remove_condition].copy() # delete_indices = np.arange(self.shape[0])[remove_condition] # face_keep_masks = self.to_maps_masks(~remove_condition) diff --git a/tfields/triangles3D.py b/tfields/triangles3D.py index 1e707063b83ec894265ca2f87a99fb0b4c6ee504..f81bb128cdfe2efb110df417ce674b9df90fe139 100644 --- a/tfields/triangles3D.py +++ b/tfields/triangles3D.py @@ -204,9 +204,8 @@ class Triangles3D(tfields.TensorFields): mp = tfields.TensorFields( np.arange(len(self)).reshape((-1, 3)), *self.fields) - return tfields.Mesh3D( - self, - maps=[mp]).cleaned() + mesh = tfields.Mesh3D(self, maps=[mp]) + return mesh.cleaned(stale=False) # stale vertices can not occure here @cached_property() def _areas(self):