Commit 5fa6ef11 authored by Daniel Boeckenhoff's avatar Daniel Boeckenhoff
Browse files

__iter__ for fast field iteration avoiding __getitem__; futhermore allowed...

__iter__ for fast field iteration avoiding __getitem__; futhermore allowed empty fields for empty template maps in cut_template
parent cbaecdea
......@@ -9,7 +9,7 @@ __all__ = [
"__license__",
"__copyright__",
]
__version__ = '0.1.0.dev9'
__version__ = '0.1.0.dev10'
__title__ = 'tfields'
__summary__ = "numpy + sympy implementation of tensor fields with attached coordinate systems"
__keywords__ = "tensors coordinate system trafo sympy numpy"
......
......@@ -126,12 +126,12 @@ class AbstractNdarray(np.ndarray):
>>> from tfields import Tensors, TensorFields
>>> scalars = Tensors([0, 1, 2])
>>> vectors = Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
>>> scalarField = TensorFields(vectors, scalars, coord_sys='cylinder')
>>> scalar_field = TensorFields(vectors, scalars, coord_sys='cylinder')
Save it and restore it
>>> out_file = NamedTemporaryFile(suffix='.pickle')
>>> pickle.dump(scalarField,
>>> pickle.dump(scalar_field,
... out_file)
>>> _ = out_file.seek(0)
......@@ -301,7 +301,7 @@ class AbstractNdarray(np.ndarray):
Recursively walk trough all __slots__ and describe all elements
"""
d = {}
d['bulk'] = np.array(self)
d['bulk'] = self.bulk
d['bulk_type'] = self.__class__.__name__
for attr in self._iter_slots():
value = getattr(self, attr)
......@@ -517,6 +517,21 @@ class Tensors(AbstractNdarray):
return obj
def __iter__(self):
"""
Forwarding iterations to the bulk array. Otherwise __getitem__ would
kick in and slow down imensely.
Examples:
>>> import tfields
>>> vectors = tfields.Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
>>> scalar_field = tfields.TensorFields(vectors, [42, 21, 10.5], [1, 2, 3])
>>> [(point.rank, point.dim) for point in scalar_field]
[(0, 1), (0, 1), (0, 1)]
"""
for index in range(len(self)):
yield super(Tensors, self).__getitem__(index).view(Tensors)
@classmethod
def merged(cls, *objects, **kwargs):
"""
......@@ -664,6 +679,14 @@ class Tensors(AbstractNdarray):
inst = cls.__new__(cls, tfields.lib.grid.igrid(*base_vectors, **kwargs))
return inst
@property
def bulk(self):
"""
The pure ndarray version of the actual state
-> nothing attached
"""
return np.array(self)
@property
def rank(self):
"""
......@@ -1203,10 +1226,10 @@ class TensorFields(Tensors):
>>> from tfields import Tensors, TensorFields
>>> scalars = Tensors([0, 1, 2])
>>> vectors = Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
>>> scalarField = TensorFields(vectors, scalars)
>>> scalarField.rank
>>> scalar_field = TensorFields(vectors, scalars)
>>> scalar_field.rank
1
>>> scalarField.fields[0].rank
>>> scalar_field.fields[0].rank
0
>>> vectorField = TensorFields(vectors, vectors)
>>> vectorField.fields[0].rank
......@@ -1288,26 +1311,26 @@ class TensorFields(Tensors):
>>> import tfields
>>> import numpy as np
>>> vectors = tfields.Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]])
>>> scalarField = tfields.TensorFields(vectors, [42, 21, 10.5], [1, 2, 3])
>>> scalar_field = tfields.TensorFields(vectors, [42, 21, 10.5], [1, 2, 3])
Slicing
>>> sliced = scalarField[2:]
>>> sliced = scalar_field[2:]
>>> assert isinstance(sliced, tfields.TensorFields)
>>> assert isinstance(sliced.fields[0], tfields.Tensors)
>>> assert sliced.fields[0].equal([10.5])
Picking
>>> picked = scalarField[1]
>>> picked = scalar_field[1]
>>> assert np.array_equal(picked, [0, 0, 1])
Masking
>>> masked = scalarField[[True, False, True]]
>>> masked = scalar_field[[True, False, True]]
>>> assert masked.equal([[0, 0, 0], [0, -1, 0]])
>>> assert masked.fields[0].equal([42, 10.5])
>>> assert masked.fields[1].equal([1, 3])
Iteration
>>> _ = [point for point in scalarField]
>>> _ = [point for point in scalar_field]
"""
item = super(TensorFields, self).__getitem__(index)
......
......@@ -590,6 +590,8 @@ class Mesh3D(tfields.TensorMaps):
... [[-42, -21], [42, 21]])
"""
# Possible Extension (small todo): check: len(field(s)) == len(self/maps)
# Redirect fields
fields = []
if template.fields:
......@@ -610,12 +612,14 @@ class Mesh3D(tfields.TensorMaps):
projected_field[nan_mask] = np.nan # correction for nan
fields.append(projected_field)
# Redirect maps fields
# Redirect maps and their fields
maps = []
for mp, template_mp in zip(self.maps, template.maps):
mp_fields = []
if len(template_mp.fields) > 0:
for field in mp.fields:
for field in mp.fields:
if len(template_mp) == 0 and len(template_mp.fields) == 0:
mp_fields.append(field[0:0]) # np.empty
else:
mp_fields.append(field[template_mp.fields[0].astype(int)])
new_mp = tfields.TensorFields(tfields.Tensors(template_mp),
*mp_fields)
......
......@@ -194,8 +194,6 @@ def plot_mesh(vertices, faces, **kwargs):
vmin
vmax
"""
vertices = np.array(vertices)
faces = np.array(faces)
if faces.shape[0] == 0:
warnings.warn("No faces to plot")
return None
......
......@@ -215,18 +215,21 @@ class Triangles3D(tfields.TensorFields):
if not np.array_equal(transform, np.eye(3)):
a = np.linalg.norm(np.linalg.solve(transform.T,
(self[aIndices, :] - self[bIndices, :]).T),
(self.bulk[aIndices, :] -
self.bulk[bIndices, :]).T),
axis=0)
b = np.linalg.norm(np.linalg.solve(transform.T,
(self[aIndices, :] - self[cIndices, :]).T),
(self.bulk[aIndices, :] -
self.bulk[cIndices, :]).T),
axis=0)
c = np.linalg.norm(np.linalg.solve(transform.T,
(self[bIndices, :] - self[cIndices, :]).T),
(self.bulk[bIndices, :] -
self.bulk[cIndices, :]).T),
axis=0)
else:
a = np.linalg.norm(self[aIndices, :] - self[bIndices, :], axis=1)
b = np.linalg.norm(self[aIndices, :] - self[cIndices, :], axis=1)
c = np.linalg.norm(self[bIndices, :] - self[cIndices, :], axis=1)
a = np.linalg.norm(self.bulk[aIndices, :] - self.bulk[bIndices, :], axis=1)
b = np.linalg.norm(self.bulk[aIndices, :] - self.bulk[cIndices, :], axis=1)
c = np.linalg.norm(self.bulk[bIndices, :] - self.bulk[cIndices, :], axis=1)
# sort by length for numerical stability
lengths = np.concatenate((a.reshape(-1, 1), b.reshape(-1, 1), c.reshape(-1, 1)), axis=1)
......@@ -245,9 +248,9 @@ class Triangles3D(tfields.TensorFields):
bIndices = [i * 3 + 1 for i in indices]
cIndices = [i * 3 + 2 for i in indices]
a = self[aIndices, :]
b = self[bIndices, :]
c = self[cIndices, :]
a = self.bulk[aIndices, :]
b = self.bulk[bIndices, :]
c = self.bulk[cIndices, :]
return a, b, c
def circumcenters(self):
......@@ -500,6 +503,7 @@ class Triangles3D(tfields.TensorFields):
ValueError: point must be castable to shape 3 but is of shape (2, 3)
"""
if self.ntriangles() == 0:
return np.array([], dtype=bool)
......@@ -537,9 +541,9 @@ class Triangles3D(tfields.TensorFields):
Returns:
np.ndarray: 2-d mask which is True, where a point is in a triangle
triangle index ->
triangle index (axis = 1)->
points index
1 0 0
(axis = 0) 1 0 0
| 1 0 0
V 0 0 1
......@@ -548,7 +552,7 @@ class Triangles3D(tfields.TensorFields):
triangle indices can have multiple true values
For Example, if you want to know the number of points in one
face, just do:
face, just sum over all points:
>> tris.in_triangles(poits).sum(axis=0)[face_index]
"""
......@@ -557,9 +561,10 @@ class Triangles3D(tfields.TensorFields):
masks = np.zeros((len(tensors), self.ntriangles()), dtype=bool)
with tensors.tmp_transform(self.coord_sys):
for i, point in enumerate(tensors):
for i, point in enumerate(iter(tensors)):
# print("Status: {i} points processed. Yet {nP} are found "
# "to be in triangles.".format(nP=masks.sum() - 1, **locals()))
# "to be in triangles (before handling assign_multiple)."
# .format(nP=masks.sum(), **locals())) # SLOW!
masks[i] = self._in_triangles(point, delta)
# print("Status: All points processed. Yet {nP} are found "
# "to be in triangles.".format(nP=masks.sum() - 1, **locals()))
......
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