mesh3D.py 39 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#!/usr/bin/env
# encoding: utf-8
"""
Author:     Daniel Boeckenhoff
Mail:       daniel.boeckenhoff@ipp.mpg.de

part of tfields library
"""
import numpy as np
import sympy
import tfields
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
12
13

# obj imports
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
14
from tfields.lib.decorators import cached_property
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
15
16
import logging
import os
17
18


19
20
21
22
23
24
25
26
27
28
29
30
31
def _dist_from_plane(point, plane):
    return plane['normal'].dot(point) + plane['d']


def _segment_plane_intersection(p0, p1, plane):
    """
    Returns:
        points, direction
    """
    distance0 = _dist_from_plane(p0, plane)
    distance1 = _dist_from_plane(p1, plane)
    p0OnPlane = abs(distance0) < np.finfo(float).eps
    p1OnPlane = abs(distance1) < np.finfo(float).eps
32
    points = []
33
34
    direction = 0
    if p0OnPlane:
35
        points.append(p0)
36
37

    if p1OnPlane:
38
        points.append(p1)
39
    # remove duplicate points
40
41
    if len(points) > 1:
        points = np.unique(points, axis=0)
42
    if p0OnPlane and p1OnPlane:
43
        return points, direction
44
45

    if distance0 * distance1 > np.finfo(float).eps:
46
        return points, direction
47
48
49

    direction = np.sign(distance0)
    if abs(distance0) < np.finfo(float).eps:
50
        return points, direction
51
    elif abs(distance1) < np.finfo(float).eps:
52
        return points, direction
53
54
55
    if abs(distance0 - distance1) > np.finfo(float).eps:
        t = distance0 / (distance0 - distance1)
    else:
56
        return points, direction
57

58
    points.append(p0 + t * (p1 - p0))
59
    # remove duplicate points
60
61
62
    if len(points) > 1:
        points = np.unique(points, axis=0)
    return points, direction
63
64
65
66
67
68
69


def _intersect(triangle, plane, vertices_rejected):
    """
    Intersect a triangle with a plane. Give the info, which side of the
    triangle is rejected by passing the mask vertices_rejected
    Returns:
70
71
72
73
74
75
76
77
78
        list of list. The inner list is of length 3 and refers to the points of
        new triangles. The reference is done with varying types:
            int: reference to triangle index
            complex: reference to duplicate point. This only happens in case
                two triangles are returned. Then only in the second triangle
            iterable: new vertex

    TODO:
        align norm vectors with previous face
79
80
81
82
83
84
85
86
    """
    nTrue = vertices_rejected.count(True)
    lonely_bool = True if nTrue == 1 else False
    index = vertices_rejected.index(lonely_bool)
    s0, d0 = _segment_plane_intersection(triangle[0], triangle[1], plane)
    s1, d1 = _segment_plane_intersection(triangle[1], triangle[2], plane)
    s2, d2 = _segment_plane_intersection(triangle[2], triangle[0], plane)

87
88
89
    single_index = index
    couple_indices = [j for j in range(3)
                      if not vertices_rejected[j] == lonely_bool]
90
91
92
93
94
95

    # TODO handle special cases. For now triangles with at least two points on plane are excluded
    new_points = None

    if len(s0) == 2:
        # both points on plane
96
        return new_points
97
98
    if len(s1) == 2:
        # both points on plane
99
        return new_points
100
101
    if len(s2) == 2:
        # both points on plane
102
        return new_points
103
    if lonely_bool:
104
        # two new triangles
105
        if len(s0) == 1 and len(s1) == 1:
106
107
            new_points = [[couple_indices[0], s0[0], couple_indices[1]],
                          [couple_indices[1], complex(1), s1[0]]]
108
        elif len(s1) == 1 and len(s2) == 1:
109
110
            new_points = [[couple_indices[0], couple_indices[1], s1[0]],
                          [couple_indices[0], complex(2), s2[0]]]
111
        elif len(s0) == 1 and len(s2) == 1:
112
113
            new_points = [[couple_indices[0], couple_indices[1], s0[0]],
                          [couple_indices[1], s2[0], complex(2)]]
114
115

    else:
116
        # one new triangle
117
        if len(s0) == 1 and len(s1) == 1:
118
            new_points = [[single_index, s1[0], s0[0]]]
119
        elif len(s1) == 1 and len(s2) == 1:
120
            new_points = [[single_index, s2[0], s1[0]]]
121
        elif len(s0) == 1 and len(s2) == 1:
122
123
            new_points = [[single_index, s0[0], s2[0]]]
    return new_points
124
125


126
127
128
129
130
131
def scalars_to_fields(scalars):
    scalars = np.array(scalars)
    if len(scalars.shape) == 1:
        return [tfields.Tensors(scalars)]
    return [tfields.Tensors(fs) for fs in scalars]

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
132

133
134
135
def fields_to_scalars(fields):
    return np.array(fields)

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
136

137
def faces_to_maps(faces, *fields):
138
    return [tfields.TensorFields(faces, *fields, dtype=int, dim=3)]
139

140

141
def maps_to_faces(maps):
142
143
144
145
146
    if len(maps) == 0:
        return np.array([])
    elif len(maps) > 1:
        raise NotImplementedError("Multiple maps")
    return np.array(maps[0])
147
148


149
150
151
152
153
class Mesh3D(tfields.TensorMaps):
    # pylint: disable=R0904
    """
    Points3D child used as vertices combined with faces to build a geometrical mesh of triangles
    Examples:
154
155
        >>> import tfields
        >>> import numpy as np
156
        >>> m = tfields.Mesh3D([[1,2,3], [3,3,3], [0,0,0], [5,6,7]], faces=[[0, 1, 2], [1, 2, 3]])
157
158
159
160
161
        >>> m.equal([[1, 2, 3],
        ...          [3, 3, 3],
        ...          [0, 0, 0],
        ...          [5, 6, 7]])
        True
162
        >>> np.array_equal(m.faces, [[0, 1, 2], [1, 2, 3]])
163
        True
164
165

        conversion to points only
166
167
168
169
170
        >>> tfields.Points3D(m).equal([[1, 2, 3],
        ...                            [3, 3, 3],
        ...                            [0, 0, 0],
        ...                            [5, 6, 7]])
        True
171
172
173
174
175
176

        Empty instances
        >>> m = tfields.Mesh3D([]);

        going from Mesh3D to Triangles3D instance is easy and will be cached.
        >>> m = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], faces=[[0, 1, 2]]);
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
177
        >>> assert m.triangles().equal(tfields.Triangles3D([[ 1.,  0.,  0.],
178
179
        ...                                               [ 0.,  1.,  0.],
        ...                                               [ 0.,  0.,  0.]]))
180
181

        a list of scalars is assigned to each face
182
183
184
        >>> mScalar = tfields.Mesh3D([[1,0,0], [0,1,0], [0,0,0]], faces=[[0, 1, 2]], faceScalars=[.5]);
        >>> np.array_equal(mScalar.faceScalars, [[ 0.5]])
        True
185
186

        adding together two meshes:
187
188
189
190
191
192
193
194
195
196
        >>> 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
197
        >>> assert np.array_equal(msum.faces, [[0, 1, 2], [3, 4, 5]])
198
199
200
201

        Saving and reading
        >>> from tempfile import NamedTemporaryFile
        >>> outFile = NamedTemporaryFile(suffix='.npz')
202
        >>> m.save(outFile.name)
203
        >>> _ = outFile.seek(0)
204
        >>> m1 = tfields.Mesh3D.load(outFile.name)
205
206
        >>> bool(np.all(m == m1))
        True
207
        >>> assert np.array_equal(m1.faces, np.array([[0, 1, 2]]))
208
209

    """
210
    def __new__(cls, tensors, *fields, **kwargs):
211
212
        if not issubclass(type(tensors), Mesh3D):
            kwargs['dim'] = 3
213
214
215
        faces = kwargs.pop('faces', None)
        faceScalars = kwargs.pop('faceScalars', [])
        maps = kwargs.pop('maps', None)
216
        if maps is not None and faces is not None:
217
218
219
            raise ValueError("Conflicting options maps and faces")
        if maps is not None:
            kwargs['maps'] = maps
220
        if len(faceScalars) > 0:
221
222
223
224
225
226
227
228
229
230
231
232
233
            map_fields = scalars_to_fields(faceScalars)
        else:
            map_fields = []
        if faces is not None:
            kwargs['maps'] = faces_to_maps(faces,
                                           *map_fields)
        obj = super(Mesh3D, cls).__new__(cls, tensors, *fields, **kwargs)
        if len(obj.maps) > 1:
            raise ValueError("Mesh3D only allows one map")
        if obj.maps and obj.maps[0].dim != 3:
            raise ValueError("Face dimension should be 3")
        return obj

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
    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

376
377
    @classmethod
    def plane(cls, *base_vectors, **kwargs):
378
379
380
381
382
383
384
        """
        Alternative constructor for creating a plane from
        Args:
            *base_vectors: see grid constructors in core. One base_vector has to
                be one-dimensional
            **kwargs: forwarded to __new__
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
385
        vertices = tfields.Tensors.grid(*base_vectors, **kwargs)
386
387
388

        base_vectors = tfields.grid.ensure_complex(*base_vectors)
        base_vectors = tfields.grid.to_base_vectors(*base_vectors)
389
390
391
392
393
394
395
        fix_coord = None
        for coord in range(3):
            if len(base_vectors[coord]) > 1:
                continue
            if len(base_vectors[coord]) == 0:
                continue
            fix_coord = coord
396
397
        if fix_coord is None:
            raise ValueError("Describe a plane with one variable fiexed")
398
399
400
401
402
403
404
405
406
407
408
409
410
411

        var_coords = list(range(3))
        var_coords.pop(var_coords.index(fix_coord))

        faces = []
        base0, base1 = base_vectors[var_coords[0]], base_vectors[var_coords[1]]
        for i1 in range(len(base1) - 1):
            for i0 in range(len(base0) - 1):
                idx_top_left = len(base1) * (i0 + 0) + (i1 + 0)
                idx_top_right = len(base1) * (i0 + 0) + (i1 + 1)
                idx_bot_left = len(base1) * (i0 + 1) + (i1 + 0)
                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])
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
412
        inst = cls.__new__(cls, vertices, faces=faces)
413
414
415
416
        return inst

    @classmethod
    def grid(cls, *base_vectors, **kwargs):
417
418
        """
        Construct 'cuboid' along base_vectors
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
        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')

447
        """
448
449
450
        if not len(base_vectors) == 3:
            raise AttributeError("3 base_vectors vectors required")

451
452
453
        base_vectors = tfields.grid.ensure_complex(*base_vectors)
        base_vectors = tfields.grid.to_base_vectors(*base_vectors)

454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
        indices = [0, -1]
        coords = range(3)
        baseLengthsAbove1 = [len(b) > 1 for b in base_vectors]
        # if one plane is given: rearrange indices and coords
        if not all(baseLengthsAbove1):
            indices = [0]
            for i, b in enumerate(baseLengthsAbove1):
                if not b:
                    coords = [i]
                    break

        base_vectors = list(base_vectors)
        planes = []
        for ind in indices:
            for coord in coords:
                basePart = base_vectors[:]
                basePart[coord] = np.array([base_vectors[coord][ind]],
                                           dtype=float)
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
472
                planes.append(cls.plane(*basePart, **kwargs))
473
474
475
        inst = cls.merged(*planes, **kwargs)
        return inst

476
477
478
479
480
481
482
483
484
485
486
487
488
    @property
    def faces(self):
        return maps_to_faces(self.maps)

    @faces.setter
    def faces(self, faces):
        self.maps = faces_to_maps(faces)

    @property
    def faceScalars(self):
        return fields_to_scalars(self.maps[0].fields)

    @faceScalars.setter
489
490
    def faceScalars(self, scalars):
        self.maps[0].fields = scalars_to_fields(scalars)
491

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
492
493
    @cached_property()
    def _triangles(self):
494
495
496
497
498
499
        """
        with the decorator, this should be handled like an attribute though it is a function

        """
        if self.faces.size == 0:
            return tfields.Triangles3D([])
500
501
502
503
        tris = tfields.Tensors.merged(*[self[mp.flatten()] for mp in self.maps])
        map_fields = [mp.fields for mp in self.maps]
        fields = [tfields.Tensors.merged(*fields) for fields in zip(*map_fields)]
        return tfields.Triangles3D(tris, *fields)
504

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
505
    def triangles(self):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
506
507
508
509
510
511
512
513
        """
        Cached method to retrieve the triangles, belonging to this mesh
        Examples:
            >>> import tfields
            >>> mesh = tfields.Mesh3D.grid((0, 1, 3), (1, 2, 3), (2, 3, 3))
            >>> assert mesh.triangles() is mesh.triangles()

        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
514
515
516
517
518
519
520
        return self._triangles

    def centroids(self):
        return self.triangles().centroids()

    @cached_property()
    def _planes(self):
521
522
        if self.faces.size == 0:
            return tfields.Planes3D([])
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
523
524
525
526
        return tfields.Planes3D(self.centroids(), self.triangles().norms())

    def planes(self):
        return self._planes
527

528
    def nfaces(self):
529
530
        return self.faces.shape[0]

531
    def in_faces(self, points, delta, assign_multiple=False):
532
533
        """
        Check whether points lie within triangles with Barycentric Technique
534
        see Triangles3D.in_triangles
535
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
536
        masks = self.triangles().in_triangles(points, delta,
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
537
                                              assign_multiple=assign_multiple)
538
539
        return masks

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
540
    def removeFaces(self, face_delete_mask):
541
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
542
        Remove faces where face_delete_mask is True
543
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
544
545
546
        face_delete_mask = np.array(face_delete_mask, dtype=bool)
        self.faces = self.faces[~face_delete_mask]
        self.faceScalars = self.faceScalars[~face_delete_mask]
547

548
    def template(self, sub_mesh):
549
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
550
        'Manual' way to build a template that can be used with self.cut
551
        Returns:
552
553
            Mesh3D: template (see cut), can be used as template to retrieve
                sub_mesh from self instance
554
        Examples:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
555
556
            >>> mp = tfields.TensorFields([[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
            ...                           [1, 2, 3, 4])
557
            >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,2,0], [1,2,0]],
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
558
            ...                     maps=[mp])
559
            >>> from sympy.abc import y
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
560
561
562
            >>> m_cut = m.cut(y < 1.5, at_intersection='split')
            >>> template = m.template(m_cut)
            >>> assert m_cut.equal(m.cut(template))
563

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
564
565
        TODO:
            fields template not yet implemented
566
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
567
568
        face_indices = np.arange(self.maps[0].shape[0])
        cents = tfields.Tensors(sub_mesh.centroids())
569
        mask = self.in_faces(cents, delta=None)
570
        inst = sub_mesh.copy()
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
571
572
573
574
575
576
577
578
579
580
581
        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)
                        ]
582
583
        return inst

584
    def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False):
585
        """
586
        Partition the mesh with the cuts given and return the template
587
588

        """
589
590
591
592
593
        eps = 0.000000001
        # direct return if self is empty
        if len(self) == 0:
            return self.copy(), self.copy()

594
595
596
597
598
599
600
601
602
603
604
        inst = self.copy()

        '''
        add the indices of the vertices and maps to the fields. They will be
        removed afterwards
        '''
        if not _in_recursion:
            inst.fields.append(tfields.Tensors(np.arange(len(inst))))
            for mp in inst.maps:
                mp.fields.append(tfields.Tensors(np.arange(len(mp))))

605
        # mask for points that do not fulfill the cut expression
606
        mask = inst.evalf(expression)
607
        # remove the points
608

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
609
        if not any(~mask):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
610
611
612
613
            # no vertex is valid
            inst = inst[mask]
        elif all(~mask):
            # all vertices are valid
614
            inst = inst[mask]
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
        elif at_intersection == 'keep':
            expression_parts = tfields.lib.symbolics.split_expression(expression)
            if len(expression_parts) > 1:
                new_mesh = inst.copy()
                for exprPart in expression_parts:
                    inst, _ = inst._cut_sympy(exprPart,
                                              at_intersection=at_intersection,
                                              _in_recursion=True)
            elif len(expression_parts) == 1:
                face_delete_indices = set([])
                for i, face in enumerate(inst.maps[0]):
                    """
                    vertices_rejected is a mask for each face that is True, where
                    a Point is on the rejected side of the plane
                    """
                    vertices_rejected = [~mask[f] for f in face]
                    if all(vertices_rejected):
                        # delete face
                        face_delete_indices.add(i)
                mask = np.full(len(inst.maps[0]), True, dtype=bool)
                for face_idx in range(len(inst.maps[0])):
                    if face_idx in face_delete_indices:
                        mask[face_idx] = False
                inst.maps[0] = inst.maps[0][mask]
            else:
                raise ValueError("Sympy expression is not splitable.")
            inst = inst.cleaned()
642
643
        elif at_intersection == 'split' or at_intersection == 'splitRough':
            '''
644
            add vertices and faces that are at the border of the cuts
645
            '''
646
            expression_parts = tfields.lib.symbolics.split_expression(expression)
647
            if len(expression_parts) > 1:
648
                new_mesh = inst.copy()
649
650
651
652
653
654
655
656
                if at_intersection == 'splitRough':
                    """
                    the following is, to speed up the process. Problem is, that
                    triangles can exist, where all points lie outside the cut,
                    but part of the area
                    still overlaps with the cut.
                    These are at the intersection line between two cuts.
                    """
657
658
                    faceIntersMask = np.full((inst.faces.shape[0]), False, dtype=bool)
                    for i, face in enumerate(inst.faces):
659
660
                        vertices_rejected = [-mask[f] for f in face]
                        face_on_edge = any(vertices_rejected) and not all(vertices_rejected)
661
                        if face_on_edge:
662
                            faceIntersMask[i] = True
663
                    new_mesh.removeFaces(-faceIntersMask)
664

665
                for exprPart in expression_parts:
666
667
668
                    inst, _ = inst._cut_sympy(exprPart,
                                              at_intersection='split',
                                              _in_recursion=True)
669
            elif len(expression_parts) == 1:
670
                # TODO maps[0] -> smthng like inst.get_map(dim=3)
671
672
673
                points = [sympy.symbols('x0, y0, z0'),
                          sympy.symbols('x1, y1, z1'),
                          sympy.symbols('x2, y2, z2')]
674
                plane_sympy = tfields.lib.symbolics.to_plane(expression)
675
676
677
                norm_sympy = np.array(plane_sympy.normal_vector).astype(float)
                d = -norm_sympy.dot(np.array(plane_sympy.p1).astype(float))
                plane = {'normal': norm_sympy, 'd': d}
678

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
679
                norm_vectors = inst.triangles().norms()
680
681
                new_points = np.empty((0, 3))
                new_faces = np.empty((0, 3))
682
                new_fields = [tfields.Tensors(np.empty((0,) + field.shape[1:]),
683
                                              coord_sys=field.coord_sys)
684
685
                              for field in inst.fields]
                new_map_fields = [[] for field in inst.maps[0].fields]
686
                new_norm_vectors = []
687
                newScalarMap = []
688
                n_new = 0
689

690
691
692
693
694
695
696
                vertices = np.array(inst)
                faces = np.array(inst.maps[0])
                fields = [np.array(field) for field in inst.fields]
                faces_fields = [np.array(field) for field in inst.maps[0].fields]

                face_delete_indices = set([])
                for i, face in enumerate(inst.maps[0]):
697
                    """
698
                    vertices_rejected is a mask for each face that is True, where
699
700
                    a Point is on the rejected side of the plane
                    """
701
                    vertices_rejected = [~mask[f] for f in face]
702
703
704
705
706
                    if any(vertices_rejected):
                        # delete face
                        face_delete_indices.add(i)
                    if any(vertices_rejected) and not all(vertices_rejected):
                        # face on edge
707
                        nTrue = vertices_rejected.count(True)
708
                        lonely_bool = True if nTrue == 1 else False
709

710
                        triangle_points = [inst[f] for f in face]
711
                        """
712
                        Add the intersection points and faces
713
                        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
714
                        intersection = _intersect(triangle_points, plane, vertices_rejected)
715
                        last_idx = len(vertices) - 1
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
716
                        for tri_list in intersection:
717
                            new_face = []
718
719
720
                            for item in tri_list:
                                if isinstance(item, int):
                                    # reference to old vertex
721
                                    new_face.append(face[item])
722
723
724
                                elif isinstance(item, complex):
                                    # reference to new vertex that has been
                                    # concatenated already
725
                                    new_face.append(last_idx + int(item.imag))
726
727
                                else:
                                    # new vertex
728
                                    new_face.append(len(vertices))
729
                                    vertices = np.append(vertices,
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
730
                                                         [[float(x) for x in item]],
731
732
733
734
735
                                                         axis=0)
                                    fields = [np.append(field,
                                                        np.full((1,) + field.shape[1:], np.nan),
                                                        axis=0)
                                              for field in fields]
736
                            faces = np.append(faces, [new_face], axis=0)
737
                            faces_fields = [np.append(field,
738
                                                      [field[i]],
739
740
741
742
743
                                                      axis=0)
                                            for field in faces_fields]
                            faces_fields[-1][-1] = i

                face_map = tfields.TensorFields(faces, *faces_fields,
744
                                                dtype=int,
745
                                                coord_sys=inst.maps[0].coord_sys)
746
747
748
                inst = tfields.Mesh3D(vertices,
                                      *fields,
                                      maps=[face_map] + inst.maps[1:],
749
                                      coord_sys=inst.coord_sys)
750
751
752
753
754
                mask = np.full(len(inst.maps[0]), True, dtype=bool)
                for face_idx in range(len(inst.maps[0])):
                    if face_idx in face_delete_indices:
                        mask[face_idx] = False
                inst.maps[0] = inst.maps[0][mask]
755
            else:
756
                raise ValueError("Sympy expression is not splitable.")
757
            inst = inst.cleaned()
758
        elif at_intersection == 'remove':
759
            inst = inst[mask]
760
        else:
761
762
            raise AttributeError("No at_intersection method called {at_intersection} "
                                 "implemented".format(**locals()))
763
764
765
766
767
768
769
770
771
772
773
774
775

        if _in_recursion:
            template = None
        else:
            template_field = inst.fields.pop(-1)
            template_maps = []
            for mp in inst.maps:
                t_mp = tfields.TensorFields(tfields.Tensors(mp),
                                            mp.fields.pop(-1))
                template_maps.append(t_mp)
            template = tfields.Mesh3D(tfields.Tensors(inst),
                                      template_field,
                                      maps=template_maps)
776
        return inst, template
777
778

    def _cut_template(self, template):
779
780
781
782
783
784
785
786
        """
        Args:
            template (tfields.Mesh3D)

        Examples:
            >>> import tfields
            >>> import numpy as np

787
            Build mesh
788
789
790
            >>> mmap = tfields.TensorFields([[0, 1, 2], [0, 3, 4]],
            ...                             [[42, 21], [-42, -21]])
            >>> m = tfields.Mesh3D([[0]*3, [1]*3, [2]*3, [3]*3, [4]*3],
791
792
            ...                    [0.0, 0.1, 0.2, 0.3, 0.4],
            ...                    [0.0, -0.1, -0.2, -0.3, -0.4],
793
794
            ...                    maps=[mmap])

795
            Build template
796
797
798
            >>> tmap = tfields.TensorFields([[0, 3, 4], [0, 1, 2]],
            ...                             [1, 0])
            >>> t = tfields.Mesh3D([[0]*3, [-1]*3, [-2]*3, [-3]*3, [-4]*3],
799
            ...                    [1, 0, 3, 2, 4],
800
801
            ...                    maps=[tmap])

802
            Use template as instruction to make a fast cut
803
804
805
806
807
808
809
810
811
            >>> res = m._cut_template(t)
            >>> assert np.array_equal(res.fields,
            ...                       [[0.1, 0.0, 0.3, 0.2, 0.4],
            ...                        [-0.1, 0.0, -0.3, -0.2, -0.4]])

            >>> assert np.array_equal(res.maps[0].fields[0],
            ...                       [[-42, -21], [42, 21]])
                                   
        """
812
813
        # Possible Extension (small todo): check: len(field(s)) == len(self/maps)

814
        # Redirect fields
815
        fields = []
816
        if template.fields:
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
            template_field = np.array(template.fields[0])
            if len(self) > 0:
                '''
                if new vertices have been created in the template, it is
                in principle unclear what fields we have to refer to.
                Thus in creating the template, we gave np.nan.
                To make it fast, we replace nan with 0 as a dummy and correct
                the field entries afterwards with np.nan.
                '''
                nan_mask = np.isnan(template_field)
                template_field[nan_mask] = 0  # dummy reference to index 0.
                template_field = template_field.astype(int)
                for field in self.fields:
                    projected_field = field[template_field]
                    projected_field[nan_mask] = np.nan  # correction for nan
                    fields.append(projected_field)
833

834
        # Redirect maps and their fields
835
836
        maps = []
        for mp, template_mp in zip(self.maps, template.maps):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
837
            mp_fields = []
838
839
840
841
            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:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
842
                    mp_fields.append(field[template_mp.fields[0].astype(int)])
843
            new_mp = tfields.TensorFields(tfields.Tensors(template_mp),
844
845
846
                                          *mp_fields)
            maps.append(new_mp)

847
848
        inst = tfields.Mesh3D(tfields.Tensors(template),
                              *fields,
849
                              maps=maps)
850
851
        return inst

852
    def cut(self, expression, coord_sys=None, at_intersection=None,
853
            return_template=False):
854
855
856
        """
        cut method for Mesh3D.
        Args:
857
858
859
860
861
862
863
864
865
866
            expression (sympy logical expression | Mesh3D):
                sympy locical expression: Sympy expression that defines planes
                    in 3D
                Mesh3D: A mesh3D will be interpreted as a template, i.e. a
                    fast instruction of how to cut the triangles.
                    It is the second part of the tuple, returned by a previous
                    cut with a sympy locial expression with 'return_template=True'.
                    We use the vertices and maps of the Mesh as the sceleton of
                    the returned mesh. The fields are mapped according to
                    indices in the template.maps[i].fields.
867
            coord_sys (coordinate system to cut in):
868
            at_intersection (str): instruction on what to do, when a cut will intersect a triangle.
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
869
870
871
                Options:    'remove' (Default) - remove the faces that are on the edge
                            'keep' - keep the faces that are on the edge
                            'split' - Create new triangles that make up the old one.
872
873
            return_template (bool): If True: return the template
                            to redo the same cut fast
874
875
        Examples:
            define the cut
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
876
            >>> import numpy as np
877
            >>> import tfields
878
            >>> from sympy.abc import x,y,z
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
879
            >>> cut_expr = x > 1.5
880

881
882
            >>> m = tfields.Mesh3D.grid((0, 3, 4),
            ...                         (0, 3, 4),
883
            ...                         (0, 0, 1))
884
885
886
887
888
889
            >>> m.fields.append(tfields.Tensors(np.linspace(0, len(m) - 1,
            ...                                             len(m))))
            >>> m.maps[0].fields.append(
            ...     tfields.Tensors(np.linspace(0,
            ...                                 len(m.maps[0]) - 1,
            ...                                 len(m.maps[0]))))
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
890
            >>> mNew = m.cut(cut_expr)
891
            >>> len(mNew)
892
            8
893
            >>> mNew.nfaces()
894
895
896
897
            6
            >>> float(mNew[:, 0].min())
            2.0

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
898
899
900
901
902
903
904
905
906
907
            Cutting with the 'keep' option will leave triangles on the edge
            untouched:
            >>> m_keep = m.cut(cut_expr, at_intersection='keep')
            >>> float(m_keep[:, 0].min())
            1.0
            >>> m_keep.nfaces()
            12

            Cutting with the 'split' option will create new triangles on the edge:
            >>> m_split = m.cut(cut_expr, at_intersection='split')
908
            >>> float(m_split[:, 0].min())
909
            1.5
910
            >>> len(m_split)
911
            15
912
            >>> m_split.nfaces()
913
914
            15

915
916
            Cut with 'return_template=True' will return the exact same mesh but
            additionally an instruction to conduct the exact same cut fast (template)
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
917
            >>> m_split_2, template = m.cut(cut_expr, at_intersection='split',
918
            ...                                    return_template=True)
919
920
921
922
923
924
925
926
927
928
            >>> m_split_template = m.cut(template)
            >>> assert m_split.equal(m_split_2, equal_nan=True)
            >>> assert m_split.equal(m_split_template, equal_nan=True)
            >>> assert len(template.fields) == 1
            >>> assert len(m_split.fields) == 1
            >>> assert len(m_split_template.fields) == 1
            >>> assert m_split.fields[0].equal(
            ...     list(range(8, 16)) + [np.nan] * 7, equal_nan=True)
            >>> assert m_split_template.fields[0].equal(
            ...     list(range(8, 16)) + [np.nan] * 7, equal_nan=True)
929
930
931
932
933
934
935
936
937
938
939

            This seems irrelevant at first but Consider, the map field or the
            tensor field changes:
            >>> m_altered_fields = m.copy()
            >>> m_altered_fields[0] += 42
            >>> assert not m_split.equal(m_altered_fields.cut(template))
            >>> assert tfields.Tensors(m_split).equal(m_altered_fields.cut(template))
            >>> assert tfields.Tensors(m_split.maps[0]).equal(m_altered_fields.cut(template).maps[0])


            The cut expression may be a sympy.BooleanFunction:
940
941
            >>> cut_expr_bool_fun = (x > 1.5) & (y < 1.5) & (y >0.2) & (z > -0.5)
            >>> m_split_bool = m.cut(cut_expr_bool_fun, at_intersection='split')
942
943
944

        Returns:
            copy of cut mesh
945
            * optional: template
946
947

        """
948
        with self.tmp_transform(coord_sys or self.coord_sys):
949
950
951
952
953
954
955
956
            if isinstance(expression, Mesh3D):
                obj = self._cut_template(expression)
            else:
                at_intersection = at_intersection or "remove"
                obj, template = self._cut_sympy(expression, at_intersection=at_intersection)
        if return_template:
            return obj, template
        return obj
957

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
958
    def plot(self, **kwargs):  # pragma: no cover
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
959
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
960
        Forwarding to plotTools.plot_mesh
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
        """
        scalars_demanded = any([v in kwargs for v in ['vmin', 'vmax', 'cmap']])
        map_index = kwargs.pop('map_index', None if not scalars_demanded else 0)
        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 'axis' in kwargs:
            dim_defined = True
        if 'zAxis' in kwargs:
            if kwargs['zAxis'] 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

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
983
        return tfields.plotting.plot_mesh(self, self.faces, **kwargs)
984

985

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
986
if __name__ == '__main__':  # pragma: no cover
987
988
    import doctest

989
    # doctest.run_docstring_examples(Mesh3D.cut, globals())
990
    doctest.testmod()