added project functionality to mesh3d

tfields/mesh3D.py

@@ -564,11 +564,12 @@ class Mesh3D(tfields.TensorMaps):
             Mesh3D: template (see cut), can be used as template to retrieve
                 sub_mesh from self instance
         Examples:
+            >>> import tfields
+            >>> from sympy.abc import y
             >>> mp = tfields.TensorFields([[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
             ...                           [1, 2, 3, 4])
             >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,2,0], [1,2,0]],
             ...                     maps=[mp])
-            >>> from sympy.abc import y
             >>> m_cut = m.cut(y < 1.5, at_intersection='split')
             >>> template = m.template(m_cut)
             >>> assert m_cut.equal(m.cut(template))
@@ -593,6 +594,93 @@ class Mesh3D(tfields.TensorMaps):
                         ]
         return inst
 
+    def project(self, tensor_field,
+                delta=None, merge_functions=None):
+        """
+        project the points of the tensor_field to a copy of the mesh
+        and set the face values accord to the field to the maps field.
+        If no field is present in tensor_field, the number of points in a mesh
+        is counted.
+
+        Args:
+            tensor_field (Tensors | TensorFields)
+            delta (float | None): forwarded to Mesh3D.in_faces
+            merge_functions (callable): if multiple Tensors lie in the same face,
+                they are mapped with the merge_function to one value
+
+        Examples:
+            >>> import tfields
+            >>> mp = tfields.TensorFields([[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
+            ...                           [1, 2, 3, 4])
+            >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,2,0], [1,2,0]],
+            ...                     maps=[mp])
+
+            Projecting points onto the mesh gives the count
+            >>> points = tfields.Tensors([[0.5, 0.2, 0.0], [0.5, 0.02, 0.0], [0.5, 0.8, 0.0]])
+            >>> m_points = m.project(points)
+            >>> assert m_points.maps[0].fields[0].equal([2, 1, 0, 0])
+
+            TensorFields with arbitrary size are projected,
+            combinging the fields automatically
+            >>> fields = [tfields.Tensors([1,3,42]),
+            ...           tfields.Tensors([[0,1,2], [2,3,4], [3,4,5]]),
+            ...           tfields.Tensors([[[0, 0]] * 2,
+            ...                            [[2, 2]] * 2,
+            ...                            [[3, 3]] * 2])]
+            >>> tf = tfields.TensorFields(points, *fields)
+            >>> m_tf = m.project(tf)
+            >>> assert m_tf.maps[0].fields[0].equal([2, 42, np.nan, np.nan],
+            ...                                     equal_nan=True)
+            >>> assert m_tf.maps[0].fields[1].equal([[1, 2, 3],
+            ...                                      [3, 4, 5],
+            ...                                      [np.nan] * 3,
+            ...                                      [np.nan] * 3],
+            ...                                     equal_nan=True)
+            >>> assert m_tf.maps[0].fields[2].equal([[[1, 1]] * 2,
+            ...                                      [[3, 3]] * 2,
+            ...                                      [[np.nan, np.nan]] * 2,
+            ...                                      [[np.nan, np.nan]] * 2],
+            ...                                     equal_nan=True)
+
+        """
+        if not issubclass(type(tensor_field), tfields.Tensors):
+            tensor_field = tfields.TensorFields(tensor_field)
+        mask = self.in_faces(tensor_field, delta=None)
+        inst = self.copy()
+        
+        n_faces = len(self.maps[0])
+        if not hasattr(tensor_field, 'fields') or len(tensor_field.fields) == 0:
+            fields = [np.full(len(tensor_field), 1)]
+            empty_map_fields = [tfields.Tensors(np.full(n_faces, 0))]
+            if merge_functions is None:
+                merge_functions = [np.sum]
+        else:
+            fields = tensor_field.fields
+            empty_map_fields = []
+            for field in fields:
+                cls = type(field)
+                kwargs = {key: getattr(field, key) for key in cls.__slots__}
+                shape = (n_faces,) + field.shape[1:]
+                empty_map_fields.append(cls(np.full(shape, np.nan),
+                                            **kwargs))
+            if merge_functions is None:
+                merge_functions = [lambda x: np.mean(x, axis=0)] * len(fields)
+
+        map_fields = []
+        for field, map_field, merge_function in \
+                zip(fields, empty_map_fields, merge_functions):
+            for f in range(len(self.maps[0])):
+                res = field[mask[:, f]]
+                if len(res) == 0:
+                    continue
+                elif len(res) == 1:
+                    map_field[f] = res
+                else:
+                    map_field[f] = merge_function(res)
+            map_fields.append(map_field)
+        inst.maps[0].fields = map_fields
+        return inst
+
     def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False):
         """
         Partition the mesh with the cuts given and return the template
@@ -1010,5 +1098,5 @@ class Mesh3D(tfields.TensorMaps):
 if __name__ == '__main__':  # pragma: no cover
     import doctest
 
-    doctest.run_docstring_examples(Mesh3D.cut, globals())
+    doctest.run_docstring_examples(Mesh3D.project, globals())
     # doctest.testmod()