set_style generic

test/test_core.py

@@ -70,6 +70,19 @@ class Base_Test(object):
     def tearDown(self):
         del self._inst
 
+
+
+class Tensor_Fields_Test(object):
+    def test_fields(self):
+        # field is of type list
+        self.assertTrue(isinstance(self._inst.fields, list))
+        self.assertTrue(len(self._inst.fields) == len(self._fields))
+
+        for field, target_field in zip(self._inst.fields, self._fields):
+            self.assertTrue(np.array_equal(field, target_field))
+            # fields are copied not reffered by a pointer
+            self.assertFalse(field is target_field)
+
 """
 EMPTY TESTS
 """
@@ -80,20 +93,10 @@ class Tensors_Empty_Test(Base_Test, unittest.TestCase):
         self._inst = tfields.Tensors([], dim=3)
 
 
-class TensorFields_Empty_Test(Tensors_Empty_Test):
+class TensorFields_Empty_Test(Tensors_Empty_Test, Tensor_Fields_Test):
     def setUp(self):
         self._fields = []
         self._inst = tfields.TensorFields([], dim=3)
-
-    def test_fields(self):
-        # field is of type list
-        self.assertTrue(isinstance(self._inst.fields, list))
-        self.assertTrue(len(self._inst.fields) == len(self._fields))
-
-        for field, target_field in zip(self._inst.fields, self._fields):
-            self.assertTrue(np.array_equal(field, target_field))
-            # fields are copied not reffered by a pointer
-            self.assertFalse(field is target_field)
     
 
 class TensorFields_Copy_Test(TensorFields_Empty_Test):
@@ -104,6 +107,9 @@ class TensorFields_Copy_Test(TensorFields_Empty_Test):
         tensors = tfields.Tensors.grid(*base)
         self._inst = tfields.TensorFields(tensors, *self._fields)
 
+        self.assertTrue(self._fields[0].coord_sys, 'cylinder')
+        self.assertTrue(self._fields[1].coord_sys, 'cartesian')
+
 
 class TensorMaps_Empty_Test(TensorFields_Empty_Test):
     def setUp(self):

tfields/core.py

@@ -650,6 +650,9 @@ class Tensors(AbstractNdarray):
             ...              [1, 4, 6],
             ...              [1, 4, 7]])
             True
+
+            Lists or arrays are accepted also.
+            Furthermore, the iteration order can be changed
             >>> lins = tfields.Tensors.grid(np.linspace(3, 4, 2), np.linspace(0, 1, 2),
             ...                             np.linspace(6, 7, 2), iter_order=[1, 0, 2])
             >>> lins.equal([[3, 0, 6],
@@ -675,8 +678,22 @@ class Tensors(AbstractNdarray):
             ...              [1, 4, 7]])
             True
 
-        """
-        inst = cls.__new__(cls, tfields.lib.grid.igrid(*base_vectors, **kwargs))
+            When given the coord_sys argument, the grid is performed in the
+            given coorinate system:
+            >>> lins3 = tfields.Tensors.grid(np.linspace(4, 9, 2),
+            ...                              np.linspace(np.pi/2, np.pi/2, 1),
+            ...                              np.linspace(4, 4, 1),
+            ...                              iter_order=[2, 0, 1],
+            ...                              coord_sys=tfields.bases.CYLINDER)
+            >>> assert lins3.coord_sys == 'cylinder'
+            >>> lins3.transform('cartesian')
+            >>> assert np.array_equal(lins3[:, 1], [4, 9])
+
+        """
+        cls_kwargs = {attr: kwargs.pop(attr) for attr in list(kwargs.keys()) if attr in cls.__slots__}
+        inst = cls.__new__(cls,
+                           tfields.lib.grid.igrid(*base_vectors, **kwargs),
+                           **cls_kwargs)
         return inst
 
     @property
@@ -1144,9 +1161,10 @@ class Tensors(AbstractNdarray):
 
         """
         indices = np.arange(self.shape[0])
-        dists = self.distances(self)
+        dists = self.distances(self)  # this takes long
         distsInEpsilon = dists <= epsilon
-        return [indices[die] for die in distsInEpsilon]
+        indices = [indices[die] for die in distsInEpsilon]  # this takes long
+        return indices
 
     def _weights(self, weights, rigid=True):
         """

tfields/mesh3D.py

@@ -9,7 +9,11 @@ part of tfields library
 import numpy as np
 import sympy
 import tfields
+
+# obj imports
 from tfields.lib.decorators import cached_property
+import logging
+import os
 
 
 def _dist_from_plane(point, plane):
@@ -227,6 +231,148 @@ class Mesh3D(tfields.TensorMaps):
             raise ValueError("Face dimension should be 3")
         return obj
 
+    def _save_obj(self, path, **kwargs):
+        """
+        Save obj as wavefront/.obj file
+        """
+        obj = kwargs.pop('object', None)
+        group = kwargs.pop('group', None)
+
+        cmap = kwargs.pop('cmap', 'viridis')
+        map_index = kwargs.pop('map_index', None)
+
+        path = path.replace('.obj', '')
+        directory, name = os.path.split(path)
+
+        if not (self.faceScalars.size == 0 or map_index is None):
+            import matplotlib.colors as colors
+            scalars = self.maps[0].fields[map_index]
+            min_scalar = scalars[~np.isnan(scalars)].min()
+            max_scalar = scalars[~np.isnan(scalars)].max()
+            vmin = kwargs.pop('vmin', min_scalar)
+            vmax = kwargs.pop('vmax', max_scalar)
+            if vmin == vmax:
+                if vmin == 0.:
+                    vmax = 1.
+                else:
+                    vmin = 0.
+            norm = colors.Normalize(vmin, vmax)
+            color_map = plt.get_cmap(cmap)
+        else:
+            # switch for not coloring the triangles and thus not producing the materials
+            norm = None
+
+        if norm is not None:
+            mat_name = name + '_frame_{0}.mat'.format(map_index)
+            scalars[np.isnan(scalars)] = min_scalar - 1
+            sorted_scalars = scalars[scalars.argsort()]
+            sorted_scalars[sorted_scalars == min_scalar - 1] = np.nan
+            sorted_faces = self.faces[scalars.argsort()]
+            scalar_set = np.unique(sorted_scalars)
+            scalar_set[scalar_set == min_scalar - 1] = np.nan
+            mat_path = os.path.join(directory, mat_name)
+            with open(mat_path, 'w') as mf:
+                for s in scalar_set:
+                    if np.isnan(s):
+                        mf.write("newmtl nan")
+                        mf.write("Kd 0 0 0\n\n")
+                    else:
+                        mf.write("newmtl mtl_{0}\n".format(s))
+                        mf.write("Kd {c[0]} {c[1]} {c[2]}\n\n".format(c=color_map(norm(s))))
+        else:
+            sorted_faces = self.faces
+
+        # writing of the obj file
+        with open(path + '.obj', 'w') as f:
+            f.write("# File saved with tfields Mesh3D._save_obj method\n\n")
+            if norm is not None:
+                f.write("mtllib ./{0}\n\n".format(mat_name))
+            if obj is not None:
+                f.write("o {0}\n".format(obj))
+            if group is not None:
+                f.write("g {0}\n".format(group))
+            for vertex in self:
+                f.write("v {v[0]} {v[1]} {v[2]}\n".format(v=vertex))
+
+            last_scalar = None
+            for i, face in enumerate(sorted_faces + 1):
+                if norm is not None:
+                    if not last_scalar == sorted_scalars[i]:
+                        last_scalar = sorted_scalars[i]
+                        f.write("usemtl mtl_{0}\n".format(last_scalar))
+                f.write("f {f[0]} {f[1]} {f[2]}\n".format(f=face))
+
+    @classmethod
+    def _load_obj(cls, path, *group_names):
+        """
+        Factory method
+        Given a path to a obj/wavefront file, construct the object
+        """
+        import csv
+        log = logging.getLogger()
+
+        with open(path, mode='r') as f:
+            reader = csv.reader(f, delimiter=' ')
+            groups = []
+            group = None
+            vertex_no = 1
+            for line in reader:
+                if not line:
+                    continue
+                if line[0] == '#':
+                    continue
+                if line[0] == 'g':
+                    if group:
+                        groups.append(group)
+                    group = dict(name=line[1], vertices={}, faces=[])
+                elif line[0] == 'v':
+                    if not group:
+                        log.warning("No group specified. I invent one myself.")
+                        group = dict(name='Group', vertices={}, faces=[])
+                    vertex = list(map(float, line[1:4]))
+                    group['vertices'][vertex_no] = vertex
+                    vertex_no += 1
+                elif line[0] == 'f':
+                    face = []
+                    for v in line[1:]:
+                        w = v.split('/')
+                        face.append(int(w[0]))
+                    group['faces'].append(face)
+
+        vertices = []
+        for g in groups[:]:
+            vertices.extend(g['vertices'].values())
+
+        if len(group_names) != 0:
+            groups = [g for g in groups if g['name'] in group_names]
+
+        faces = []
+        for g in groups:
+            faces.extend(g['faces'])
+        faces = np.add(np.array(faces), -1).tolist()
+
+        """
+        Building the class from retrieved vertices and faces
+        """
+        if len(vertices) == 0:
+            return cls([])
+        faceLenghts = [len(face) for face in faces]
+        for i in reversed(range(len(faceLenghts))):
+            length = faceLenghts[i]
+            if length == 3:
+                continue
+            if length == 4:
+                log.warning("Given a Rectangle. I will split it but "
+                            "sometimes the order is different.")
+                faces.insert(i + 1, faces[i][2:] + faces[i][:1])
+                faces[i] = faces[i][:3]
+            else:
+                raise NotImplementedError()
+        mesh = cls(vertices, faces=faces)
+        if group_names:
+            mesh = mesh.cleaned()
+        return mesh
+
     @classmethod
     def plane(cls, *base_vectors, **kwargs):
         """
@@ -236,7 +382,7 @@ class Mesh3D(tfields.TensorMaps):
                 be one-dimensional
             **kwargs: forwarded to __new__
         """
-        vertices = tfields.Tensors.grid(*base_vectors)
+        vertices = tfields.Tensors.grid(*base_vectors, **kwargs)
 
         base_vectors = tfields.grid.ensure_complex(*base_vectors)
         base_vectors = tfields.grid.to_base_vectors(*base_vectors)
@@ -263,13 +409,41 @@ class Mesh3D(tfields.TensorMaps):
                 idx_bot_right = len(base1) * (i0 + 1) + (i1 + 1)
                 faces.append([idx_top_left, idx_top_right, idx_bot_left])
                 faces.append([idx_top_right, idx_bot_left, idx_bot_right])
-        inst = cls.__new__(cls, vertices, faces=faces, **kwargs)
+        inst = cls.__new__(cls, vertices, faces=faces)
         return inst
 
     @classmethod
     def grid(cls, *base_vectors, **kwargs):
         """
         Construct 'cuboid' along base_vectors
+        Examples:
+            Building symmetric geometries were never as easy:
+
+            Approximated sphere with radius 1, translated in y by 2 units
+            >>> sphere = tfields.Mesh3D.grid((1, 1, 1),
+            ...                              (-np.pi, np.pi, 12),
+            ...                              (-np.pi / 2, np.pi / 2, 12),
+            ...                              coord_sys='spherical')
+            >>> sphere.transform('cartesian')
+            >>> sphere[:, 1] += 2
+
+            Oktaeder
+            >>> oktder = tfields.Mesh3D.grid((1, 1, 1),
+            ...                              (-np.pi, np.pi, 5),
+            ...                              (-np.pi / 2, np.pi / 2, 3),
+            ...                              coord_sys='spherical')
+
+            Cube with edge length of 2 units
+            >>> cube = tfields.Mesh3D.grid((-1, 1, 2),
+            ...                            (-1, 1, 2),
+            ...                            (-5, -3, 2))
+
+            Cylinder 
+            >>> cylinder = tfields.Mesh3D.grid((1, 1, 1),
+            ...                                (-np.pi, np.pi, 12),
+            ...                                (-5, 3, 12),
+            ...                                coord_sys='cylinder')
+
         """
         if not len(base_vectors) == 3:
             raise AttributeError("3 base_vectors vectors required")
@@ -295,8 +469,7 @@ class Mesh3D(tfields.TensorMaps):
                 basePart = base_vectors[:]
                 basePart[coord] = np.array([base_vectors[coord][ind]],
                                            dtype=float)
-
-                planes.append(cls.plane(*basePart))
+                planes.append(cls.plane(*basePart, **kwargs))
         inst = cls.merged(*planes, **kwargs)
         return inst
 
@@ -352,16 +525,16 @@ class Mesh3D(tfields.TensorMaps):
         Check whether points lie within triangles with Barycentric Technique
         """
         masks = self.triangles().in_triangles(points, delta,
-                                            assign_multiple=assign_multiple)
+                                              assign_multiple=assign_multiple)
         return masks
 
-    def removeFaces(self, faceDeleteMask):
+    def removeFaces(self, face_delete_mask):
         """
-        Remove faces where faceDeleteMask is True
+        Remove faces where face_delete_mask is True
         """
-        faceDeleteMask = np.array(faceDeleteMask, dtype=bool)
-        self.faces = self.faces[~faceDeleteMask]
-        self.faceScalars = self.faceScalars[~faceDeleteMask]
+        face_delete_mask = np.array(face_delete_mask, dtype=bool)
+        self.faces = self.faces[~face_delete_mask]
+        self.faceScalars = self.faceScalars[~face_delete_mask]
 
     def template(self, sub_mesh, delta=1e-9):
         """
@@ -384,13 +557,19 @@ class Mesh3D(tfields.TensorMaps):
         """
         face_indices = np.arange(self.maps[0].shape[0])
         cents = tfields.Tensors(sub_mesh.centroids())
-        scalars = []
-        mask = self.in_faces(cents, delta=delta)
-        scalars = []
-        for face_mask in mask:
-            scalars.append(face_indices[face_mask][0])
+        mask = self.in_faces(cents, delta)
         inst = sub_mesh.copy()
-        inst.maps[0].fields = [tfields.Tensors(scalars)]
+        if inst.maps:
+            scalars = []
+            for face_mask in mask:
+                scalars.append(face_indices[face_mask][0])
+            inst.maps[0].fields = [tfields.Tensors(scalars, dim=1)]
+        else:
+            inst.maps = [tfields.TensorFields([],
+                                              tfields.Tensors([], dim=1),
+                                              dim=3,
+                                              dtype=int)
+                        ]
         return inst
 
     def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False):

tfields/plotting/__init__.py

@@ -10,6 +10,7 @@ import matplotlib.pyplot as plt
 import matplotlib as mpl
 import numpy as np
 from .mpl import *
+from six import string_types
 
 
 def set_default(dictionary, attr, value):
@@ -121,7 +122,7 @@ class PlotOptions(object):
         vmax = self.get('vmax', vmaxDefault)
         return cmap, vmin, vmax
 
-    def formatColors(self, colors, fmt='rgba', length=None):
+    def format_colors(self, colors, fmt='rgba', length=None):
         """
         format colors according to fmt argument
         Args:
@@ -134,10 +135,11 @@ class PlotOptions(object):
             colors in fmt
         """
         hasIter = True
-        if not hasattr(colors, '__iter__'):
+        if not hasattr(colors, '__iter__') or isinstance(colors, string_types):
             # colors is just one element
             hasIter = False
             colors = [colors]
+
         if hasattr(colors[0], '__iter__') and fmt == 'norm':
             # rgba given but norm wanted
             cmap, vmin, vmax = self.getNormArgs(cmapDefault='NotSpecified',
@@ -148,9 +150,9 @@ class PlotOptions(object):
             self.plotKwargs['vmax'] = vmax
             self.plotKwargs['cmap'] = cmap
         elif fmt == 'rgba':
-            if isinstance(colors[0], str) or isinstance(colors[0], unicode):
+            if isinstance(colors[0], string_types):
                 # string color defined
-                colors = map(mpl.colors.to_rgba, colors)
+                colors = [mpl.colors.to_rgba(color) for color in colors]
             else:
                 # norm given rgba wanted
                 cmap, vmin, vmax = self.getNormArgs(cmapDefault='NotSpecified',
@@ -218,3 +220,8 @@ class PlotOptions(object):
         if len(args) != 1:
             raise ValueError("Invalid number of args ({0})".format(len(args)))
         return self.retrieve(args[0], default, keep=keep)
+
+
+if __name__ == '__main__':
+    import doctest
+    doctest.testmod()

tfields/plotting/mpl.py

@@ -7,6 +7,7 @@ import numpy as np
 import warnings
 import os
 import matplotlib as mpl
+from matplotlib import style
 import matplotlib.pyplot as plt
 from matplotlib.patches import Circle
 import mpl_toolkits.mplot3d as plt3D
@@ -37,7 +38,7 @@ def gca(dim=None, **kwargs):
     return axis
 
 
-def upgrade_style(style, source, dest="~/.config/matplotlib/"):
+def upgrade_style(style, source, dest=None):
     """
     Copy a style file at <origionalFilePath> to the <dest> which is the foreseen
     local matplotlib rc dir by default
@@ -48,35 +49,42 @@ def upgrade_style(style, source, dest="~/.config/matplotlib/"):
         dest (str): local directory to copy the file to. Matpotlib has to
             search this directory for mplstyle files!
     """
-    styleExtension = 'mplstyle'
-    path = tfields.lib.in_out.resolve(os.path.join(dest, style + '.' + styleExtension))
+    if dest is None:
+        dest = mpl.get_configdir()
+    style_extension = 'mplstyle'
+    path = tfields.lib.in_out.resolve(os.path.join(dest, style + '.' +
+                                                   style_extension))
     source = tfields.lib.in_out.resolve(source)
     tfields.lib.in_out.cp(source, path)
 
 
-def set_style(style='tfields', dest="~/.config/matplotlib/"):
+def set_style(style='tfields', dest=None):
     """
     Set the matplotlib style of name
     Important:
         Either you
     Args:
         style (str)
-        dest (str): local directory to use file from. if None, use default maplotlib styles
+        dest (str): local directory to use file from. if None, use default
+            maplotlib destination
     """
     if dest is None:
-        path = style
-    else:
-        styleExtension = 'mplstyle'
-        path = tfields.lib.in_out.resolve(os.path.join(dest, style + '.' + styleExtension))
-    try:
+        dest = mpl.get_configdir()
+
+    style_extension = 'mplstyle'
+    path = tfields.lib.in_out.resolve(os.path.join(dest, style + '.' +
+                                                   style_extension))
+    if style in mpl.style.available:
+        plt.style.use(style)
+    elif os.path.exists(path):
         plt.style.use(path)
-    except IOError:
+    else:
         log = logging.getLogger()
         if style == 'tfields':
             log.warning("I will copy the default style to {dest}."
                         .format(**locals()))
             source = os.path.join(os.path.dirname(__file__),
-                                  style + '.' + styleExtension)
+                                  style + '.' + style_extension)
             try:
                 upgrade_style(style, source, dest)                
                 set_style(style)
@@ -163,6 +171,7 @@ def plot_array(array, **kwargs):
         Artist or list of Artists (imitating the axis.scatter/plot behaviour).
         Better Artist not list of Artists
     """
+    array = np.array(array)
     tfields.plotting.set_default(kwargs, 'methodName', 'scatter')
     po = tfields.plotting.PlotOptions(kwargs)
 
@@ -194,6 +203,8 @@ 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
@@ -241,9 +252,9 @@ def plot_mesh(vertices, faces, **kwargs):
         """
         sort out color arguments
         """
-        facecolors = po.formatColors(facecolors,
-                                     fmt='norm',
-                                     length=nFacesInitial)
+        facecolors = po.format_colors(facecolors,
+                                      fmt='norm',
+                                      length=nFacesInitial)
         if not full:
             facecolors = facecolors[dotProduct > 0]
         po.plotKwargs['facecolors'] = facecolors
@@ -257,9 +268,9 @@ def plot_mesh(vertices, faces, **kwargs):
         color = po.retrieve_chain('color', 'c', 'facecolors',
                                   default='grey',
                                   keep=False)
-        color = po.formatColors(color,
-                                fmt='rgba',
-                                length=len(faces))
+        color = po.format_colors(color,
+                                 fmt='rgba',
+                                 length=len(faces))
         nanMask = np.isnan(color)
         if nanMask.any():
             warnings.warn("nan found in colors. Removing the corresponding faces!")

tfields/triangles3D.py

@@ -215,21 +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.bulk[aIndices, :] -
-                                                    self.bulk[bIndices, :]).T),
+                                                   (self[aIndices, :] -
+                                                    self[bIndices, :]).T),
                                    axis=0)
                 b = np.linalg.norm(np.linalg.solve(transform.T,
-                                                   (self.bulk[aIndices, :] -
-                                                    self.bulk[cIndices, :]).T),
+                                                   (self[aIndices, :] -
+                                                    self[cIndices, :]).T),
                                    axis=0)
                 c = np.linalg.norm(np.linalg.solve(transform.T,
-                                                   (self.bulk[bIndices, :] -
-                                                    self.bulk[cIndices, :]).T),
+                                                   (self[bIndices, :] -
+                                                    self[cIndices, :]).T),
                                    axis=0)
             else:
-                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)
+                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)
 
             # sort by length for numerical stability
             lengths = np.concatenate((a.reshape(-1, 1), b.reshape(-1, 1), c.reshape(-1, 1)), axis=1)
@@ -459,7 +459,9 @@ class Triangles3D(tfields.TensorFields):
         Barycentric method to optain, wheter a point is in any of the triangles
         Args:
             point (list of len 3)
-            delta (float): acceptance in +- norm vector direction
+            delta (float / None):
+                float: acceptance in +- norm vector direction
+                None: accept the face with the minimum distance to the point
         Returns:
             np.array: boolean mask, True where point in a triangle within delta
         Examples:
@@ -487,6 +489,11 @@ class Triangles3D(tfields.TensorFields):
             ...     m2.triangles()._in_triangles(np.array([0.2, 0.2, 0.1]), 0.2),
             ...     np.array([ True, False], dtype=bool))
 
+            if you set delta to None, the minimal distance point(s) are accepted
+            >>> assert np.array_equal(
+            ...     m2.triangles()._in_triangles(np.array([0.2, 0.2, 0.1]), None),
+            ...     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]],
@@ -515,16 +522,33 @@ class Triangles3D(tfields.TensorFields):
         except:
             raise
 
+        # min_dist_method switch if delta is None
+        if delta is None:
+            delta = 1.        
+            min_dist_method = True
+        else:
+            min_dist_method = False
+
         u, v, w = self._baricentric(point, delta=delta).T
 
         if delta == 0.:
             w[np.isnan(w)] = 0.  # division by 0 in baricentric makes w = 0 nan.
 
+
         with warnings.catch_warnings():
             warnings.filterwarnings('ignore', message="invalid value encountered in less_equal")
-            orthogonalAcceptance = (abs(w) <= 1)
             barycentric_bools = ((v <= 1.) & (v >= 0.)) & ((u <= 1.) & (u >= 0.)) & ((v + u <= 1.0))
-        return np.array(barycentric_bools & orthogonalAcceptance)
+            if min_dist_method:
+                orthogonal_acceptance = np.full(barycentric_bools.shape, False)
+                closest_indices = np.argmin(abs(w)[barycentric_bools])
+                # transform the indices to the whole array, not only the
+                # barycentric_bools selection