Plane test in mesh3d

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):

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)

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):