Commit da5c84dc authored by Daniel Boeckenhoff's avatar Daniel Boeckenhoff
Browse files

project functionality added to mesh

parent b6df77dc
......@@ -639,7 +639,9 @@ class Mesh3D(tfields.TensorMaps):
return inst
def project(self, tensor_field,
delta=None, merge_functions=None):
delta=None, merge_functions=None,
point_face_assignment=None,
return_point_face_assignment=False):
"""
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.
......@@ -651,6 +653,11 @@ class Mesh3D(tfields.TensorMaps):
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
point_face_assignment (np.array, dtype=int): array assigning each
point to a face. Each entry position corresponds to a point of the
tensor, each entry value is the index of the assigned face
return_point_face_assignment (bool): if true, return the
point_face_assignment
Examples:
>>> import tfields
......@@ -661,19 +668,23 @@ class Mesh3D(tfields.TensorMaps):
... 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)
>>> points = tfields.Tensors([[0.5, 0.2, 0.0],
... [0.5, 0.02, 0.0],
... [0.5, 0.8, 0.0],
... [0.5, 0.8, 0.1]]) # not contained
>>> m_points = m.project(points, delta=0.01)
>>> 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]]),
>>> fields = [tfields.Tensors([1,3,42, -1]),
... tfields.Tensors([[0,1,2], [2,3,4], [3,4,5], [-1] * 3]),
... tfields.Tensors([[[0, 0]] * 2,
... [[2, 2]] * 2,
... [[3, 3]] * 2])]
... [[3, 3]] * 2,
... [[9, 9]] * 2])]
>>> tf = tfields.TensorFields(points, *fields)
>>> m_tf = m.project(tf)
>>> m_tf = m.project(tf, delta=0.01)
>>> 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],
......@@ -687,13 +698,21 @@ class Mesh3D(tfields.TensorMaps):
... [[np.nan, np.nan]] * 2],
... equal_nan=True)
Returning the calculated point_face_assignment can speed up multiple
results
>>> m_tf, point_face_assignment = m.project(tf, delta=0.01,
... return_point_face_assignment=True)
>>> m_tf_fast = m.project(tf, delta=0.01, point_face_assignment=point_face_assignment)
>>> assert m_tf.equal(m_tf_fast, 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()
# setup empty map fields and collect fields
n_faces = len(self.maps[0])
point_indices = np.arange(len(tensor_field))
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))]
......@@ -710,20 +729,43 @@ class Mesh3D(tfields.TensorMaps):
**kwargs))
if merge_functions is None:
merge_functions = [lambda x: np.mean(x, axis=0)] * len(fields)
# build point_face_assignment if not given.
if point_face_assignment is not None:
if len(point_face_assignment) != len(tensor_field):
raise ValueErrror("Template needs same lenght as tensor_field")
point_face_assignment = point_face_assignment
else:
mask = self.in_faces(tensor_field, delta=delta)
face_indices = np.arange(n_faces)
point_face_assignment = [face_indices[mask[p_index, :]]
for p_index in point_indices]
point_face_assignment = np.array([fa if len(fa) != 0 else [-1]
for fa in point_face_assignment])
point_face_assignment = point_face_assignment.reshape(-1)
point_face_assignment_set = set(point_face_assignment)
# merge the fields according to point_face_assignment
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:
for i, f_index in enumerate(point_face_assignment_set):
if i % 1000 == 0:
print(i, len(point_face_assignment_set))
if f_index == -1:
# point could not be mapped
continue
elif len(res) == 1:
map_field[f] = res
point_in_face_indices = point_indices[point_face_assignment == f_index]
res = field[point_in_face_indices]
if len(res) == 1:
map_field[f_index] = res
else:
map_field[f] = merge_function(res)
map_field[f_index] = merge_function(res)
map_fields.append(map_field)
inst.maps[0].fields = map_fields
if return_point_face_assignment:
return inst, point_face_assignment
return inst
def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False):
......@@ -1121,22 +1163,6 @@ class Mesh3D(tfields.TensorMaps):
if map_index is not None:
if not len(self.maps[0]) == 0:
kwargs['color'] = self.maps[0].fields[map_index]
dim_defined = False
if 'axes' in kwargs:
dim_defined = True
if 'z_index' in kwargs:
if kwargs['z_index'] is not None:
kwargs['dim'] = 3
else:
kwargs['dim'] = 2
dim_defined = True
if 'dim' in kwargs:
dim_defined = True
if not dim_defined:
kwargs['dim'] = 2
return rna.plotting.plot_mesh(self, self.faces, **kwargs)
......
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