Commit c87f626a by Daniel Böckenhoff (Laptop)

### found problem with copy constructor with child

parent ce425b9d
 ... ... @@ -15,6 +15,7 @@ from .lib import * # NOQA # classes: from .points3D import Points3D # NOQA from .mesh3D import Mesh3D # NOQA from .mesh3D import fields_to_scalars, scalars_to_fields from .triangles3D import Triangles3D # NOQA from .scalarField3D import ScalarField3D # NOQA from .vectorField3D import VectorField3D # NOQA ... ...
 import tfields import sympy import numpy as np import warnings CARTESIAN = 'cartesian' CYLINDER = 'cylinder' ... ... @@ -24,6 +25,18 @@ cylinder.connect_to(cartesian, Z], inverse=False, fill_in_gaps=False) def cylinder_to_cartesian(array): rPoints = np.copy(array[:, 0]) phiPoints = np.copy(array[:, 1]) array[:, 0] = rPoints * np.cos(phiPoints) array[:, 1] = rPoints * np.sin(phiPoints) cylinder.to_cartesian = cylinder_to_cartesian cartesian.connect_to(cylinder, [x, y, z], [sympy.sqrt(x**2 + y**2), ... ... @@ -36,27 +49,19 @@ cartesian.connect_to(cylinder, fill_in_gaps=False) def cylinder_to_cartesian(array): rPoints = np.copy(array[:, 0]) phiPoints = np.copy(array[:, 1]) array[:, 0] = rPoints * np.cos(phiPoints) array[:, 1] = rPoints * np.sin(phiPoints) def cartesian_to_cylinder(array): x_vals = np.copy(array[:, 0]) y_vals = np.copy(array[:, 1]) problemPhiIndices = np.where((x_vals < 0) & (y_vals == 0)) problem_phi_indices = np.where((x_vals < 0) & (y_vals == 0)) array[:, 0] = np.sqrt(x_vals**2 + y_vals**2) # r np.seterr(divide='ignore', invalid='ignore') # phi array[:, 1] = np.sign(y_vals) * np.arccos(x_vals / array[:, 0]) np.seterr(divide='warn', invalid='warn') array[:, 1][problemPhiIndices] = np.pi array[:, 1][problem_phi_indices] = np.pi tfields.convert_nan(array, 0.) # for cases like cartesian 0, 0, 1 cylinder.to_cartesian = cylinder_to_cartesian cartesian.to_cylinder = cartesian_to_cylinder # spherical ... ... @@ -64,20 +69,81 @@ r, phi, theta = sympy.symbols('r, phi, theta') spherical = sympy.diffgeom.CoordSystem(SPHERICAL, patch, ['r', 'phi', 'theta']) spherical.connect_to(cartesian, [r, phi, theta], [r * sympy.cos(phi) * sympy.sin(theta), r * sympy.sin(phi) * sympy.sin(theta), r * sympy.cos(theta)], [r * sympy.cos(phi - np.pi) * sympy.sin(theta - np.pi / 2), r * sympy.sin(phi - np.pi) * sympy.sin(theta - np.pi / 2), r * sympy.cos(theta - np.pi / 2)], # [r * sympy.cos(phi) * sympy.sin(theta), # r * sympy.sin(phi) * sympy.sin(theta), # r * sympy.cos(theta)], inverse=False, fill_in_gaps=False) def spherical_to_cartesian(array): r = np.copy(array[:, 0]) phi = np.copy(array[:, 1]) theta = np.copy(array[:, 2]) array[:, 0] = r * np.sin(theta - np.pi / 2) * np.cos(phi - np.pi) array[:, 1] = r * np.sin(theta - np.pi / 2) * np.sin(phi - np.pi) array[:, 2] = r * np.cos(theta - np.pi / 2) # classical # array[:, 0] = r * np.sin(theta) * np.cos(phi) # array[:, 1] = r * np.sin(theta) * np.sin(phi) # array[:, 2] = r * np.cos(theta) spherical.to_cartesian = spherical_to_cartesian cartesian.connect_to(spherical, [x, y, z], [sympy.sqrt(x**2 + y**2 + z**2), sympy.Piecewise( (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0))), (sympy.atan(y / x), True)), (sympy.pi, ((x < 0) & sympy.Eq(y, 0))), (0., sympy.Eq(sympy.sqrt(x**2 + y**2), 0)), (sympy.sign(y) * sympy.acos(x / sympy.sqrt(x**2 + y**2)), True)), sympy.Piecewise( (0., sympy.Eq(x**2 + y**2 + z**2, 0)), # (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0) & sympy.Eq(z, 0))), (sympy.acos(z / sympy.sqrt(x**2 + y**2 + z**2)), True))], (sympy.atan2(z, sympy.sqrt(x**2 + y**2)), True))], # sympy.Piecewise( # (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0))), # (sympy.atan(y / x), True)), # sympy.Piecewise( # (0., sympy.Eq(x**2 + y**2 + z**2, 0)), # # (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0) & sympy.Eq(z, 0))), # (sympy.acos(z / sympy.sqrt(x**2 + y**2 + z**2)), True))], inverse=False, fill_in_gaps=False) def cartesian_to_spherical(array): """ convert array to r, phi, theta r in [0, oo] phi in [-pi, pi] theta element [0, pi]: elevation angle defined from - Z-axis up """ xy = array[:, 0]**2 + array[:, 1]**2 # phi for phi between -pi, pi problemPhiIndices = np.where((array[:, 0] < 0) & (array[:, 1] == 0)) with warnings.catch_warnings(): warnings.filterwarnings('ignore', message="invalid value encountered in divide") array[:, 1] = np.sign(array[:, 1]) * np.arccos(array[:, 0] / np.sqrt(xy)) array[:, 1][problemPhiIndices] = np.pi tfields.convert_nan(array, 0.) # for cases like cartesian 0, 0, 1 # array[:,1] = np.arctan2(array[:,1], array[:,0]) # phi for phi between 0, 2pi # r array[:, 0] = np.sqrt(xy + array[:, 2]**2) # theta # for elevation angle defined from - Z-axis up: # array[:, 2] = np.arctan2(np.sqrt(xy), array[:, 2]) # for elevation angle defined from XY-plane up: array[:, 2] = np.arctan2(array[:, 2], np.sqrt(xy)) cartesian.to_spherical = cartesian_to_spherical
 ... ... @@ -360,7 +360,8 @@ class Tensors(AbstractNdarray): order = kwargs.pop('order', None) ''' copy constructor ''' if issubclass(type(tensors), cls): print type(tensors), cls, issubclass(type(tensors), cls) if issubclass(type(tensors), cls) or issubclass(cls, type(tensors)): obj = tensors.copy() if dtype != obj.dtype or order is not None: obj = obj.astype(dtype, order=order) ... ... @@ -386,7 +387,9 @@ class Tensors(AbstractNdarray): ''' process empty inputs ''' if len(tensors) == 0: if dim is not None: if issubclass(type(tensors), tfields.Tensors): tensors = np.empty(tensors.shape, dtype=tensors.dtype) elif dim is not None: tensors = np.empty((0, dim)) else: raise ValueError("Empty tensors need dimension " ... ... @@ -430,16 +433,41 @@ class Tensors(AbstractNdarray): >>> import numpy as np >>> import tfields >>> import tfields.bases >>> vecA = tfields.Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]]) >>> vecB = tfields.Tensors([[5, 4, 1]], coordSys=tfields.bases.cylinder) >>> vecC = tfields.Tensors([[5, 4, 1]], coordSys=tfields.bases.cylinder) >>> merge = tfields.Tensors.merged(vecA, vecB, vecC, [[2, 0, 1]]) >>> assert merge.coordSys == 'cylinder' Use of most frequent coordinate system >>> vec_a = tfields.Tensors([[0, 0, 0], [0, 0, 1], [0, -1, 0]]) >>> vec_b = tfields.Tensors([[5, 4, 1]], coordSys=tfields.bases.cylinder) >>> vec_c = tfields.Tensors([[4, 2, 3]], coordSys=tfields.bases.cylinder) >>> merge = tfields.Tensors.merged(vec_a, vec_b, vec_c, [[2, 0, 1]]) >>> assert merge.coordSys == 'cylinder' >>> assert merge.equal([[0, 0, 0], ... [0, 0, 1], ... [1, -np.pi / 2, 0], ... [5, 4, 1], ... [4, 2, 3], ... [2, 0, 1]]) >>> tm_a = tfields.TensorMaps(merge, maps=[[[0, 1, 2]]]) >>> tm_b = tm_a.copy() >>> tm_a.coordSys >>> tm_merge = tfields.TensorMaps.merged(tm_a, tm_b) >>> assert tm_merge.coordSys == 'cylinder' >>> assert tm_merge.maps[0].equal([[0, 1, 2], ... list(range(len(merge), ... len(merge) + 3, ... 1))]) >>> obj_list = [tfields.Tensors([[1, 2, 3]], coordSys=tfields.bases.CYLINDER), ... tfields.Tensors([[3] * 3]), ... tfields.Tensors([[5, 1, 3]])] >>> merge2 = tfields.Tensors.merged(*obj_list, coordSys=tfields.bases.CARTESIAN) >>> assert merge2.equal([[-0.41614684, 0.90929743, 3.], ... [3, 3, 3], [5, 1, 3]], atol=1e-8) """ ''' get most frequent coordSys or predefined coordSys ''' coordSys = kwargs.get('coordSys', None) dimension = kwargs.get('dim', None) if coordSys is None: bases = [] for t in objects: ... ... @@ -468,6 +496,11 @@ class Tensors(AbstractNdarray): for i, obj in enumerate(remainingObjects): tensors = np.append(tensors, obj, axis=0) if len(tensors) == 0 and dimension is None: for obj in objects: kwargs['dim'] = dim(obj) return cls.__new__(cls, tensors, **kwargs) @classmethod ... ... @@ -565,10 +598,10 @@ class Tensors(AbstractNdarray): >>> assert t[1, 1] == 0. >>> assert t[2, 1] == 0. theta is 0 at (0, 0, 1) and pi at (0, 0, -1) >>> assert round(t[1, 2], 10) == round(np.pi / 2, 10) >>> assert t[2, 2] == np.pi >>> assert t[3, 2] == 0. theta is 0 at (0, 0, 1) and pi / 2 at (0, 0, -1) >>> assert round(t[1, 2], 10) == round(0, 10) >>> assert t[2, 2] == -np.pi / 2 >>> assert t[3, 2] == np.pi / 2 theta is defined 0 for R == 0 >>> assert t[4, 0] == 0. ... ... @@ -979,6 +1012,14 @@ class TensorFields(Tensors): """ Discrete Tensor Field Args: tensors (array): base tensors *fields (array): multiple fields assigned to one base tensor. Fields themself are also of type tensor **kwargs: rigid (bool): demand equal field and tensor lenght ... : see tfields.Tensors Examples: >>> from tfields import Tensors, TensorFields >>> scalars = Tensors([0, 1, 2]) ... ... @@ -1013,15 +1054,34 @@ class TensorFields(Tensors): >>> cp = TensorFields(vectorField) >>> assert vectorField.equal(cp) Copying takes care of coordSys >>> cp.transform(tfields.bases.CYLINDER) >>> cp_cyl = TensorFields(cp) >>> assert cp_cyl.coordSys == tfields.bases.CYLINDER Copying with changing type >>> tcp = TensorFields(vectorField, dtype=int) >>> assert vectorField.equal(tcp) >>> assert tcp.dtype == int Raises: TypeError: >>> import tfields >>> tfields.TensorFields([1, 2, 3], [3]) # doctest: +ELLIPSIS Traceback (most recent call last): ... ValueError: Length of base (3) should be the same as the length of all fields ([1]). This error can be suppressed by setting rigid=False >>> loose = tfields.TensorFields([1, 2, 3], [3], rigid=False) """ __slots__ = ['coordSys', 'fields'] def __new__(cls, tensors, *fields, **kwargs): rigid = kwargs.pop('rigid', True) obj = super(TensorFields, cls).__new__(cls, tensors, **kwargs) if issubclass(type(tensors), TensorFields): if tensors.fields is None: ... ... @@ -1033,12 +1093,13 @@ class TensorFields(Tensors): # (over)write fields obj.fields = [Tensors(field) for field in fields] if not all([len(f) == len(obj) for f in obj.fields]): if rigid: olen = len(obj) field_lengths = [len(f) for f in obj.fields] raise ValueError("Length of base ({olen}) should be the same as" " the length of all fields ({field_lengths})." .format(**locals())) if not all([flen == olen for flen in field_lengths]): raise ValueError("Length of base ({olen}) should be the same as" " the length of all fields ({field_lengths})." .format(**locals())) return obj def __getitem__(self, index): ... ... @@ -1162,6 +1223,16 @@ class TensorMaps(TensorFields): >>> assert mesh.equal(tfields.TensorFields(vectors, scalars)) >>> assert mesh.maps[0].fields[0].equal(maps[0].fields[0]) Copy constructor >>> mesh_copy = tfields.TensorMaps(mesh) Copying takes care of coordSys >>> mesh_copy.transform(tfields.bases.CYLINDER) >>> mesh_copy.coordSys >>> mesh_cp_cyl = TensorMaps(mesh_copy) >>> mesh_cp_cyl.coordSys >>> assert mesh_cp_cyl.coordSys == tfields.bases.CYLINDER Raises: >>> import tfields >>> tfields.TensorMaps([1] * 4, dim=3, maps=[[1, 2, 3]]) # +doctest: ELLIPSIS ... ... @@ -1181,10 +1252,37 @@ class TensorMaps(TensorFields): raise ValueError("Incorrect map rank {mp.rank}" .format(**locals())) maps_cp.append(mp) kwargs['maps'] = maps_cp maps = maps_cp obj = super(TensorMaps, cls).__new__(cls, tensors, *fields, **kwargs) obj.maps = maps return obj @classmethod def merged(cls, *objects, **kwargs): if not all([isinstance(o, cls) for o in objects]): # TODO: could allow if all faceScalars are none raise TypeError("Merge constructor only accepts {cls} instances." .format(**locals())) tensor_lengths = [len(o) for o in objects] cum_tensor_lengths = [sum(tensor_lengths[:i]) for i in range(len(objects))] maps = [] dims = [] for i, o in enumerate(objects): for map_field in o.maps: map_field = map_field + cum_tensor_lengths[i] try: map_index = dims.index(map_field.dim) except ValueError: maps.append(map_field) dims.append(map_field.dim) else: maps[map_index] = TensorFields.merged(maps[map_index], map_field) kwargs['maps'] = maps obj = super(TensorMaps, cls).merged(*objects, **kwargs) return cls.__new__(cls, obj, maps=maps) def equal(self, other, **kwargs): """ Test, whether the instance has the same content as other. ... ... @@ -1367,8 +1465,8 @@ class TensorMaps(TensorFields): if __name__ == '__main__': # pragma: no cover import doctest doctest.testmod() # doctest.run_docstring_examples(TensorMaps.cleaned, 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())
 ... ... @@ -39,7 +39,11 @@ def faces_to_maps(faces, *fields): return [tfields.TensorFields(faces, *fields, dtype=int)] def maps_to_faces(maps): return maps[0] if len(maps) == 0: return np.array([]) elif len(maps) > 1: raise NotImplementedError("Multiple maps") return np.array(maps[0]) class Mesh3D(tfields.TensorMaps): ... ... @@ -47,13 +51,15 @@ class Mesh3D(tfields.TensorMaps): """ Points3D child used as vertices combined with faces to build a geometrical mesh of triangles Examples: >>> import tfields >>> import numpy as np >>> m = tfields.Mesh3D([[1,2,3], [3,3,3], [0,0,0], [5,6,7]], faces=[[0, 1, 2], [1, 2, 3]]) >>> m.equal([[1, 2, 3], ... [3, 3, 3], ... [0, 0, 0], ... [5, 6, 7]]) True >>> m.faces.equal([[0, 1, 2], [1, 2, 3]]) >>> np.array_equal(m.faces, [[0, 1, 2], [1, 2, 3]]) True conversion to points only ... ... @@ -79,15 +85,16 @@ class Mesh3D(tfields.TensorMaps): True adding together two meshes: >>> m2 = tfields.Mesh3D([[1,0,0],[2,0,0],[0,3,0]], faces=[[0,1,2]], faceScalars=[.7]) >>> msum = tfields.Mesh3D([mScalar, m2]) >>> msum Mesh3D([[ 1., 0., 0.], [ 0., 1., 0.], [ 0., 0., 0.], [ 1., 0., 0.], [ 2., 0., 0.], [ 0., 3., 0.]]) >>> m2 = tfields.Mesh3D([[1,0,0],[2,0,0],[0,3,0]], ... faces=[[0,1,2]], faceScalars=[.7]) >>> msum = tfields.Mesh3D.merged(mScalar, m2) >>> msum.equal([[ 1., 0., 0.], ... [ 0., 1., 0.], ... [ 0., 0., 0.], ... [ 1., 0., 0.], ... [ 2., 0., 0.], ... [ 0., 3., 0.]]) True >>> msum.faces array([[0, 1, 2], [3, 4, 5]]) ... ... @@ -95,9 +102,9 @@ class Mesh3D(tfields.TensorMaps): Saving and reading >>> from tempfile import NamedTemporaryFile >>> outFile = NamedTemporaryFile(suffix='.npz') >>> m.savez(outFile.name) >>> m.save(outFile.name) >>> _ = outFile.seek(0) >>> m1 = tfields.Mesh3D.createFromFile(outFile.name) >>> m1 = tfields.Mesh3D.load(outFile.name) >>> bool(np.all(m == m1)) True >>> m1.faces ... ... @@ -111,7 +118,7 @@ class Mesh3D(tfields.TensorMaps): faces = kwargs.pop('faces', None) faceScalars = kwargs.pop('faceScalars', []) maps = kwargs.pop('maps', None) if maps and faces: if maps is not None and faces is not None: raise ValueError("Conflicting options maps and faces") if maps is not None: kwargs['maps'] = maps ... ... @@ -142,8 +149,8 @@ class Mesh3D(tfields.TensorMaps): return fields_to_scalars(self.maps[0].fields) @faceScalars.setter def faceScalars(self, faceScalars): self.maps[0].faces = scalars_to_fields(faceScalars) def faceScalars(self, scalars): self.maps[0].fields = scalars_to_fields(scalars) @classmethod def _updateSlotKwargs(cls, kwargs, skipCache=True): ... ...
 ... ... @@ -22,17 +22,15 @@ class Triangles3D(tfields.TensorFields): """ Points3D child restricted to n * 3 Points. 3 Points always group together to one triangle. Examples: >>> t = tfields.Triangles3D([[1,2,3], [3,3,3], [0,0,0]]); t Triangles3D([[ 1., 2., 3.], [ 3., 3., 3.], [ 0., 0., 0.]]) >>> t = tfields.Triangles3D([[1,2,3], [3,3,3], [0,0,0]]) >>> t.triangleScalars array([], shape=(1, 0), dtype=float64) >>> _ = t.coordinateTransform(t.CYLINDER) >>> a = [tfields.Points3D([[1],[2],[3]], coordSys=tfields.Points3D.CYLINDER), ... tfields.Points3D([[3],[3],[3]]), tfields.Points3D([[5],[1],[3]])] >>> tfields.Triangles3D(a, coordSys=tfields.Points3D.CARTESIAN) >>> _ = t.transform(tfields.bases.CYLINDER) >>> a = [tfields.Points3D([[1, 2, 3]], coordSys=tfields.bases.CYLINDER), ... tfields.Points3D([[3] * 3]), ... tfields.Points3D([[5, 1, 3]])] >>> tfields.Triangles3D.merged(*a, coordSys=tfields.bases.CARTESIAN) Triangles3D([[-0.41614684, 0.90929743, 3. ], [ 3. , 3. , 3. ], [ 5. , 1. , 3. ]]) ... ... @@ -41,74 +39,40 @@ class Triangles3D(tfields.TensorFields): >>> t.appendScalars([1]) """ __slots__ = ['coordSys', 'triangleScalars', '_cache'] def __new__(cls, inputArray, triangleScalars=None, **kwargs): if any([isinstance(a, Triangles3D) for a in inputArray]) and triangleScalars is None: if not all([isinstance(a, Triangles3D) for a in inputArray]): # TODO: could allow if all triangleScalars are none raise TypeError("Copy constructor of Triangles3D only" "accepts Triangles3D instances.") if not len(set([o.triangleScalars.shape for o in inputArray])) == 1: raise TypeError("Shape of triangleScalars that " "should be concatenated must be the same.") triangleScalars = np.concatenate([o.triangleScalars for o in inputArray], axis=0) obj = tfields.Points3D.__new__(cls, inputArray, **kwargs) if not obj.getNPoints() % 3 == 0: raise TypeError("Points 3D does not describe triangles. " "Size of obj is {0}, has no divider 3".format(len(obj))) # Add triangleScalars if triangleScalars is None: triangleScalars = np.empty((len(obj) // 3, 0)) triangleScalars = np.array(triangleScalars) if not len(triangleScalars) == obj.getNPoints() / 3: raise ValueError("Length of triangleScalars has to be the same as" "number of triangles ({0}) but is {1}." .format(obj.getNPoints() / 3, len(triangleScalars))) if len(triangleScalars.shape) == 1: triangleScalars = triangleScalars.reshape(-1, 1) obj.triangleScalars = triangleScalars # __slots__ = ['coordSys', 'triangleScalars', '_cache'] def __new__(cls, tensors, *fields, **kwargs): if not issubclass(type(tensors), Triangles3D): kwargs['dim'] = 3 kwargs['rigid'] = False obj = super(Triangles3D, cls).__new__(cls, tensors, *fields, **kwargs) if not len(obj) % 3 == 0: raise TypeError("Input object of size({0}) has no divider 3 and" " does not describe triangles." .format(len(obj))) return obj @classmethod def merged(cls, *objects, **kwargs): obj = tfields.TensorFields.merged(*objects, **kwargs) return cls.__new__(cls, obj) @property def triangleScalars(self): """ Deprecated scalars implementation """ return tfields.fields_to_scalars(self.fields) @triangleScalars.setter def triangleScalars(self, scalars): self.fields = tfields.scalars_to_fields(scalars) def saveTxt(self, filePath): self.getMesh3D().save(filePath) @classmethod def createFromTxtFile(cls, filePath): with ioTools.TextFile(filePath, 'r') as f: matrix = f.process(seperator=' ', lineBreak='\n') parsedMatrix = [] for line in matrix: if line[0].startswith('#'): continue if '' in line: pyTools.array.removeValues(line, '') if not line: continue elif len(line) == 3: parsedMatrix.append(map(float, line)) elif len(line) == 1: parsedMatrix.append([eval(line[0])]) else: raise ValueError("Length of line is not 1 or 3 but %s" % len(line)) infoPeriodicity = 4 scalars = [x[0] for i, x in enumerate(parsedMatrix) if i % infoPeriodicity == 0] pyTools.array.removePositions(parsedMatrix, [i for i in range(len(parsedMatrix)) if i % infoPeriodicity == 0]) return cls(parsedMatrix, triangleScalars=scalars) def convertNaN(self, value=0.): super(Triangles3D, self).convertNaN(value) nanIndicesScalars = np.isnan(self.triangleScalars) self.triangleScalars[nanIndicesScalars] = value def appendScalars(self, scalars): """ Examples: ... ...
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!