Commit 8b71fffa authored by Daniel Böckenhoff (Laptop)'s avatar Daniel Böckenhoff (Laptop)
Browse files

triangles neat, mesh3D messy

parent 64f4270c
......@@ -7,17 +7,17 @@ Mail: daniel.boeckenhoff@ipp.mpg.de
core of tfields library
contains numpy ndarray derived bases of the tfields package
"""
import tfields.bases
import numpy as np
import warnings
import os
import pathlib
from six import string_types
from contextlib import contextmanager
from collections import Counter
import numpy as np
import sympy
import scipy as sp
import scipy.spatial # NOQA: F401
import os
from six import string_types
import pathlib
import warnings
import tfields.bases
np.seterr(all='warn', over='raise')
......@@ -356,19 +356,23 @@ class Tensors(AbstractNdarray):
__slot_setters__ = [tfields.bases.get_coord_system_name]
def __new__(cls, tensors, **kwargs):
dtype = kwargs.pop('dtype', np.float64)
dtype = kwargs.pop('dtype', None)
order = kwargs.pop('order', None)
dim = kwargs.pop('dim', None)
''' copy constructor extracts the kwargs from tensors'''
if issubclass(type(tensors), Tensors):
dtype = tensors.dtype
if dim is not None:
dim = tensors.dim
coordSys = kwargs.pop('coordSys', tensors.coordSys)
tensors = tensors.copy()
tensors.transform(coordSys)
kwargs['coordSys'] = coordSys
if dtype is None:
dtype = tensors.dtype
else:
if dtype is None:
dtype = np.float64
''' demand iterable structure '''
try:
......@@ -443,7 +447,7 @@ class Tensors(AbstractNdarray):
Merge also shifts the maps to still refer to the same tensors
>>> tm_a = tfields.TensorMaps(merge, maps=[[[0, 1, 2]]])
>>> tm_b = tm_a.copy()
>>> tm_a.coordSys
>>> assert tm_a.coordSys == 'cylinder'
>>> tm_merge = tfields.TensorMaps.merged(tm_a, tm_b)
>>> assert tm_merge.coordSys == 'cylinder'
>>> assert tm_merge.maps[0].equal([[0, 1, 2],
......@@ -853,11 +857,14 @@ class Tensors(AbstractNdarray):
>>> import numpy
>>> import sympy
>>> x, y, z = sympy.symbols('x y z')
>>> p = tfields.Points3D([[1., 2., 3.], [4., 5., 6.], [1, 2, -6],
... [-5, -5, -5], [1,0,-1], [0,1,-1]])
>>> p = tfields.Tensors([[1., 2., 3.], [4., 5., 6.], [1, 2, -6],
... [-5, -5, -5], [1,0,-1], [0,1,-1]])
>>> np.array_equal(p.evalf(x > 0),
... [True, True, True, False, True, False])
True
>>> np.array_equal(p.evalf(x >= 0),
... [True, True, True, False, True, True])
True
And combination
>>> np.array_equal(p.evalf((x > 0) & (y < 3)),
......@@ -1181,14 +1188,18 @@ class TensorFields(Tensors):
"""
item = super(TensorFields, self).__getitem__(index)
if issubclass(type(item), TensorFields):
if isinstance(index, slice):
item.fields = [field.__getitem__(index) for field in item.fields]
elif isinstance(index, tuple):
item.fields = [field.__getitem__(index[0]) for field in item.fields]
else:
if item.fields:
try:
if issubclass(type(item), TensorFields):
if isinstance(index, tuple):
index = index[0]
if isinstance(index, slice):
item.fields = [field.__getitem__(index) for field in item.fields]
else:
if item.fields:
item.fields = [field.__getitem__(index) for field in item.fields]
except IndexError as err:
warnings.warn("Index error occured for field.__getitem__. Error "
"message: {err}".format(**locals()))
return item
......@@ -1527,8 +1538,8 @@ class TensorMaps(TensorFields):
if __name__ == '__main__': # pragma: no cover
import doctest
# doctest.testmod()
doctest.run_docstring_examples(TensorMaps, globals())
doctest.testmod()
# doctest.run_docstring_examples(TensorMaps, globals())
# doctest.run_docstring_examples(TensorFields, globals())
# doctest.run_docstring_examples(AbstractNdarray.copy, globals())
# doctest.run_docstring_examples(TensorMaps.equal, globals())
......@@ -42,9 +42,13 @@ def evalf(array, cutExpression=None, coords=None):
If array of other shape than (?, 3) is given, the coords need to be specified
>>> a0, a1 = sympy.symbols('a0 a1')
>>> assert np.array_equal(tfields.evalf([[0., 1.], [-1, 3]], a1 > 2,
... coords=[a0, a1]),
... coords=[a0, a1]),
... np.array([False, True], dtype=bool))
>= is taken care of
>>> assert np.array_equal(tfields.evalf(a, x >= 0),
... np.array([ True, True, True, False, True, True]))
"""
if isinstance(array, list):
array = np.array(array)
......@@ -57,7 +61,6 @@ def evalf(array, cutExpression=None, coords=None):
coords = sympy.symbols('x y z')
else:
raise ValueError("coords are None and shape is not (?, 3)")
elif len(coords) != array.shape[1]:
raise ValueError("Length of coords is not {0} but {1}".format(array.shape[1], len(coords)))
......@@ -65,7 +68,7 @@ def evalf(array, cutExpression=None, coords=None):
cutExpression,
modules={'&': np.logical_and, '|': np.logical_or})
mask = np.array([preMask(*x) for x in array], dtype=bool)
mask = np.array([preMask(*vals) for vals in array], dtype=bool)
return mask
......
This diff is collapsed.
......@@ -39,6 +39,67 @@ class Triangles3D(tfields.TensorFields):
.format(len(obj)))
return obj
def ntriangles(self):
"""
Returns:
int: number of triangles
"""
return len(self) // 3
def _to_triangles_mask(self, mask):
mask = np.array(mask)
mask = mask.reshape((self.ntriangles(), 3))
mask = mask.all(axis=1)
return mask
def __getitem__(self, index):
"""
In addition to the usual, also slice fields
Examples:
>>> import tfields
>>> import numpy as np
>>> vectors = tfields.Tensors(np.array([range(30)] * 3).T)
>>> triangles = tfields.Triangles3D(vectors, range(10))
>>> assert np.array_equal(triangles[3:6],
... [[3] * 3,
... [4] * 3,
... [5] * 3])
>>> assert triangles[3:6].fields[0][0] == 1
"""
item = super(tfields.TensorFields, self).__getitem__(index)
try: # __iter__ will try except __getitem__(i) until IndexError
if issubclass(type(item), Triangles3D): # block int, float, ...
if len(item) % 3 != 0:
item = tfields.Tensors(item)
elif item.fields:
# build triangle index / indices / mask when possible
tri_index = None
if isinstance(index, tuple):
index = index[0]
if isinstance(index, int):
pass
elif isinstance(index, slice):
start = index.start or 0
stop = index.stop or len(self)
step = index.step
if start % 3 == 0 and (stop - start) % 3 == 0 and step is None:
tri_index = slice(start // 3, stop // 3)
else:
tri_index = self._to_triangles_mask(index)
# apply triangle index to fields
if not tri_index is None:
item.fields = [field.__getitem__(tri_index)
for field in item.fields]
else:
item = tfields.Tensors(item)
except IndexError as err:
warnings.warn("Index error occured for field.__getitem__. Error "
"message: {err}".format(**locals()))
return item
@classmethod
def merged(cls, *objects, **kwargs):
with warnings.catch_warnings():
......@@ -61,20 +122,25 @@ class Triangles3D(tfields.TensorFields):
def triangleScalars(self, scalars):
self.fields = tfields.scalars_to_fields(scalars)
def get_ntriangles(self):
"""
Returns:
int: number of triangles
"""
return len(self) // 3
def evalf(self, expression=None, coordSys=None):
"""
Triangle3D implementation
Examples:
>>> from sympy.abc import x
>>> t = tfields.Triangles3D([[1., 2., 3.], [-4., 5., 6.], [1, 2, -6],
... [5, -5, -5], [1,0,-1], [0,1,-1],
... [-5, -5, -5], [1,0,-1], [0,1,-1]])
>>> mask = t.evalf(x >= 0)
>>> assert np.array_equal(t[mask],
... tfields.Triangles3D([[ 5., -5., -5.],
... [ 1., 0., -1.],
... [ 0., 1., -1.]]))
Returns:
np.array: mask which is True, where expression evaluates True
"""
mask = super(Triangles3D, self).evalf(expression, coordSys=coordSys)
mask = mask.reshape((self.get_ntriangles(), 3))
mask = mask.all(axis=1)
mask = self._to_triangles_mask(mask)
mask = np.array([mask] * 3).T.reshape((len(self)))
return mask
......@@ -86,25 +152,23 @@ class Triangles3D(tfields.TensorFields):
>>> import tfields
>>> x, y, z = sympy.symbols('x y z')
>>> t = tfields.Triangles3D([[1., 2., 3.], [-4., 5., 6.], [1, 2, -6],
... [5, -5, -5], [1,0,-1], [0,1,-1], [-5, -5, -5],
... [1,0,-1], [0,1,-1]])
... [5, -5, -5], [1, 0, -1], [0, 1, -1],
... [-5, -5, -5], [1, 0, -1], [0, 1, -1]])
>>> tc = t.cut(x >= 0)
>>> tc
Triangles3D([[ 5., -5., -5.],
[ 1., 0., -1.],
[ 0., 1., -1.]])
>>> tc.triangleScalars
array([], shape=(1, 0), dtype=float64)
>>> t.appendScalars([1,2,3])
>>> assert tc.equal(tfields.Triangles3D([[ 5., -5., -5.],
... [ 1., 0., -1.],
... [ 0., 1., -1.]]))
>>> assert np.array_equal(tc.triangleScalars, [])
>>> t.fields.append(tfields.Tensors([1,2,3]))
>>> tc2 = t.cut(x >= 0)
>>> tc2.triangleScalars
array([[ 2.]])
>>> assert np.array_equal(tc2.triangleScalars, np.array([[2.]]))
"""
mask = self.evalf(expression, coordSys=coordSys)
inst = self[mask].copy()
inst.triangleScalars = inst.triangleScalars[mask.reshape((self.get_ntriangles(), 3)).T[0]]
# inst.triangleScalars =
# inst.triangleScalars[mask.reshape((self.ntriangles(), 3)).T[0]]
return inst
......@@ -113,7 +177,7 @@ class Triangles3D(tfields.TensorFields):
"""
Cached method to retrieve areas of triangles
"""
transform = np.eye()
transform = np.eye(3)
return self.areas(transform=transform)
def areas(self, transform=None):
......@@ -124,31 +188,30 @@ class Triangles3D(tfields.TensorFields):
The triangle points are transformed with transform if given
before calclulating the area
Examples:
>>> m = tfields.Mesh3D([[1,0,0], [0,0,1], [0,0,0]], [[0, 1, 2]]);
>>> m.triangles.areas()
array([ 0.5])
>>> m = tfields.Mesh3D([[1,0,0], [0,0,1], [0,0,0]],
... faces=[[0, 1, 2]])
>>> assert np.allclose(m.triangles.areas(), np.array([0.5]))
>>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0], [0,0,1]], [[0, 1, 2], [1, 2, 3]]);
>>> m.triangles.areas()
array([ 0.5, 0.5])
>>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0], [0,0,1]],
... faces=[[0, 1, 2], [1, 2, 3]])
>>> assert np.allclose(m.triangles.areas(), np.array([0.5, 0.5]))
>>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [1,1,0], [0,0,1], [1,0,1]],
... [[0, 1, 2], [0, 3, 4]]);
>>> m.triangles.areas()
array([ 0.5, 0.5])
... faces=[[0, 1, 2], [0, 3, 4]])
>>> assert np.allclose(m.triangles.areas(), np.array([0.5, 0.5]))
"""
if transform is None:
return self._areas
else:
indices = range(self.get_ntriangles())
indices = range(self.ntriangles())
aIndices = [i * 3 for i in indices]
bIndices = [i * 3 + 1 for i in indices]
cIndices = [i * 3 + 2 for i in indices]
# define 3 vectors building the triangle, transform it back and take their norm
if not np.array_equal(transform, np.eye):
if not np.array_equal(transform, np.eye(3)):
a = np.linalg.norm(np.linalg.solve(transform.T,
(self[aIndices, :] - self[bIndices, :]).T),
axis=0)
......@@ -175,7 +238,7 @@ class Triangles3D(tfields.TensorFields):
Returns:
three np.arrays with corner points of triangles
"""
indices = range(self.get_ntriangles())
indices = range(self.ntriangles())
aIndices = [i * 3 for i in indices]
bIndices = [i * 3 + 1 for i in indices]
cIndices = [i * 3 + 2 for i in indices]
......@@ -191,12 +254,13 @@ class Triangles3D(tfields.TensorFields):
Examples:
>>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]],
... [[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]);
>>> m.triangles.circumcenters()
Points3D([[ 5.00000000e-01, 5.00000000e-01, 0.00000000e+00],
[ -5.00000000e-01, 5.00000000e-01, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 2.22044605e-16],
[ 3.33333333e-01, 3.33333333e-01, 3.33333333e-01]])
... faces=[[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]);
>>> assert np.allclose(
... m.triangles.circumcenters(),
... [[0.5, 0.5, 0.0],
... [-0.5, 0.5, 0.0],
... [0.0, 0.0, 0.0],
... [1.0 / 3, 1.0 / 3, 1.0 / 3]])
"""
pointA, pointB, pointC = self.corners()
......@@ -221,7 +285,7 @@ class Triangles3D(tfields.TensorFields):
this version is faster but takes much more ram also.
So much that i get memory error with a 32 GB RAM
"""
nT = self.get_ntriangles()
nT = self.ntriangles()
mat = np.ones((1, 3)) / 3.0
# matrix product calculatesq center of all triangles
return tfields.Points3D(np.dot(mat, self.reshape(nT, 3, 3))[0], coordSys=self.coordSys)
......@@ -239,7 +303,7 @@ class Triangles3D(tfields.TensorFields):
"""
Examples:
>>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]],
... [[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]);
... faces=[[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]);
>>> m.triangles.centroids()
Points3D([[ 0.33333333, 0.33333333, 0. ],
[-0.33333333, 0.33333333, 0. ],
......@@ -267,12 +331,13 @@ class Triangles3D(tfields.TensorFields):
"""
Examples:
>>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]],
... [[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]);
>>> m.triangles.norms()
array([[ 0. , 0. , 1. ],
[ 0. , 0. , -1. ],
[ 0. , 1. , 0. ],
[ 0.57735027, 0.57735027, 0.57735027]])
... faces=[[0, 1, 3],[0, 2, 3],[1,2,4], [1, 3, 4]]);
>>> assert np.allclose(m.triangles.norms(),
... [[0.0, 0.0, 1.0],
... [0.0, 0.0, -1.0],
... [0.0, 1.0, 0.0],
... [0.57735027] * 3],
... atol=1e-8)
"""
ab, ac = self.edges()
......@@ -303,21 +368,25 @@ class Triangles3D(tfields.TensorFields):
... [0.1, 0.0, 0.1]])
if point lies in the plane, return np.nan, else inf for w if delta is exactly 0.
>>> assert np.array_equal(m2.triangles._baricentric(np.array([0.2, 0.2, 0]),
... delta=0.),
... [[0.1, 0.1, np.nan],
... [0.1, 0.0, np.inf]])
>>> baric = m2.triangles._baricentric(np.array([0.2, 0.2, 0]),
... delta=0.),
>>> baric_expected = np.array([[0.1, 0.1, np.nan],
... [0.1, 0.0, np.inf]])
>>> assert ((baric == baric_expected) | (np.isnan(baric) &
... np.isnan(baric_expected))).all()
If you define triangles that have colinear side vectors or in general lead to
not invertable matrices [ab, ac, ab x ac] the values will be nan
>>> m3 = tfields.Mesh3D([[0,0,0], [2,0,0], [4,0,0], [0,1,0]], [[0, 1, 2], [0, 1, 3]]);
>>> m3 = tfields.Mesh3D([[0,0,0], [2,0,0], [4,0,0], [0,1,0]], faces=[[0, 1, 2], [0, 1, 3]]);
>>> bc = m3.triangles._baricentric(np.array([0.2, 0.2, 0]), delta=0.3)
>>> assert np.array_equal(bc, [[ np.nan, np.nan, np.nan], [ 0.1, 0.2, 0. ]])
>>> bc_exp = np.array([[ np.nan, np.nan, np.nan], [ 0.1, 0.2, 0. ]])
>>> assert ((bc == bc_exp) | (np.isnan(bc) &
... np.isnan(bc_exp))).all()
Returns:
np.ndarray: barycentric coordinates u, v, w of point with respect to each triangle
"""
if self.get_ntriangles() == 0:
if self.ntriangles() == 0:
return np.array([])
a, _, _ = self.corners()
......@@ -375,9 +444,10 @@ class Triangles3D(tfields.TensorFields):
Returns:
np.array: boolean mask, True where point in a triangle within delta
Examples:
>>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], [[0, 1, 2]]);
>>> m.triangles._in_triangles(tfields.Points3D([[0.2, 0.2, 0]]))
array([ True], dtype=bool)
>>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], faces=[[0, 1, 2]]);
>>> assert np.array_equal(
... m.triangles._in_triangles(tfields.Points3D([[0.2, 0.2, 0]])),
... np.array([True], dtype=bool))
wrong format of point will raise a ValueError
>>> m.triangles._in_triangles(tfields.Points3D([[0.2, 0.2, 0], [0.2, 0.2, 0]]))
......@@ -387,32 +457,44 @@ class Triangles3D(tfields.TensorFields):
All Triangles are tested
>>> m2 = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0], [4,0,0], [4, 4, 0], [8, 0, 0]],
... [[0, 1, 2], [3, 4, 5]]);
... faces=[[0, 1, 2], [3, 4, 5]]);
>>> m2.triangles._in_triangles(np.array([0.2, 0.2, 0]))
array([ True, False], dtype=bool)
>>> m2.triangles._in_triangles(np.array([5, 2, 0]))
array([False, True], dtype=bool)
>>> assert np.array_equal(
... m2.triangles._in_triangles(np.array([0.2, 0.2, 0])),
... np.array([True, False], dtype=bool))
>>> assert np.array_equal(
... m2.triangles._in_triangles(np.array([5, 2, 0])),
... np.array([False, True], dtype=bool))
delta allows to accept points that lie within delta orthogonal to the tringle plain
>>> m2.triangles._in_triangles(np.array([0.2, 0.2, 9000]), 0.0)
array([False, False], dtype=bool)
>>> m2.triangles._in_triangles(np.array([0.2, 0.2, 0.1]), 0.2)
array([ True, False], dtype=bool)
>>> assert np.array_equal(
... m2.triangles._in_triangles(np.array([0.2, 0.2, 9000]), 0.0),
... np.array([False, False], dtype=bool))
>>> assert np.array_equal(
... m2.triangles._in_triangles(np.array([0.2, 0.2, 0.1]), 0.2),
... np.array([ True, False], dtype=bool))
If you define triangles that have colinear side vectors or in general lead to
not invertable matrices the you will always get False
>>> m3 = tfields.Mesh3D([[0,0,0], [2,0,0], [4,0,0], [0,1,0]], [[0, 1, 2], [0, 1, 3]]);
>>> m3 = tfields.Mesh3D([[0,0,0], [2,0,0], [4,0,0], [0,1,0]],
... faces=[[0, 1, 2], [0, 1, 3]]);
>>> import pytest
>>> with pytest.warns(UserWarning) as record:
... mask = m3.triangles._in_triangles(np.array([0.2, 0.2, 0]), delta=0.3)
>>> mask
array([False, True], dtype=bool)
>>> mask = m3.triangles._in_triangles(np.array([0.2, 0.2, 0]), delta=0.3)
>>> assert np.array_equal(mask,
... np.array([False, True], dtype=bool))
"""
if self.get_ntriangles() == 0:
if self.ntriangles() == 0:
return np.array([], dtype=bool)
try:
point = np.reshape(point, 3)
except ValueError:
raise ValueError("point must be castable to shape 3 but is of shape {0}"
.format(point.shape))
except:
raise
u, v, w = self._baricentric(point, delta=delta).T
if delta == 0.:
......@@ -451,10 +533,10 @@ class Triangles3D(tfields.TensorFields):
"""
log = logging.getLogger(__name__)
if self.get_ntriangles() == 0:
if self.ntriangles() == 0:
return np.empty((tensors.shape[0], 0), dtype=bool)
masks = np.zeros((len(tensors), self.get_ntriangles()), dtype=bool)
masks = np.zeros((len(tensors), self.ntriangles()), dtype=bool)
with tensors.tmp_transform(self.coordSys):
for i, point in enumerate(tensors):
if i % 1000 == 0:
......@@ -495,15 +577,17 @@ class Triangles3D(tfields.TensorFields):
of a triangle
Examples:
>>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [-1,0,0], [0,1,0], [0,0,1]],
... [[0, 1, 3],[0, 2, 3],[1,2,4]]);
... faces=[[0, 1, 3],[0, 2, 3],[1,2,4]]);
Corner points are found
>>> m.triangles.on_edges(tfields.Points3D([[0,1,0]]))
array([ True, True, False], dtype=bool)
>>> assert np.array_equal(
... m.triangles._on_edges(tfields.Points3D([[0,1,0]])),
... np.array([ True, True, False], dtype=bool))
Side points are found, too
>>> m.triangles.on_edges(tfields.Points3D([[0.5,0,0.5]]))
array([False, False, True], dtype=bool)
>>> assert np.array_equal(
... m.triangles._on_edges(tfields.Points3D([[0.5,0,0.5]])),
... np.array([False, False, True], dtype=bool))
"""
u, v, w = self._baricentric(point, 1.).T
......@@ -528,7 +612,7 @@ class Triangles3D(tfields.TensorFields):
"""
# set weights to 1.0 if weights is None
if weights is None:
weights = np.ones(self.get_ntriangles())
weights = np.ones(self.ntriangles())
return super(Triangles3D, self)._weights(weights, rigid=rigid)
......
Supports Markdown
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