mesh3D.py 39.9 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

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

94
95
96
97
    print(triangle, plane, vertices_rejected)
    print(lonely_bool, index)
    print(s0, s1, s2)
    print(d0, d1, d2)
98
99
    if len(s0) == 2:
        # both points on plane
100
        return new_points
101
102
    if len(s1) == 2:
        # both points on plane
103
        return new_points
104
105
    if len(s2) == 2:
        # both points on plane
106
        return new_points
107
    if lonely_bool:
108
        # two new triangles
109
        if len(s0) == 1 and len(s1) == 1:
110
111
            new_points = [[couple_indices[0], s0[0], couple_indices[1]],
                          [couple_indices[1], complex(1), s1[0]]]
112
        elif len(s1) == 1 and len(s2) == 1:
113
114
            new_points = [[couple_indices[0], couple_indices[1], s1[0]],
                          [couple_indices[0], complex(2), s2[0]]]
115
        elif len(s0) == 1 and len(s2) == 1:
116
117
            new_points = [[couple_indices[0], couple_indices[1], s0[0]],
                          [couple_indices[1], s2[0], complex(2)]]
118
    else:
119
        # one new triangle
120
        if len(s0) == 1 and len(s1) == 1:
121
            new_points = [[single_index, s1[0], s0[0]]]
122
        elif len(s1) == 1 and len(s2) == 1:
123
            new_points = [[single_index, s2[0], s1[0]]]
124
        elif len(s0) == 1 and len(s2) == 1:
125
126
            new_points = [[single_index, s0[0], s2[0]]]
    return new_points
127
128


129
130
131
132
133
134
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
135

136
137
138
def fields_to_scalars(fields):
    return np.array(fields)

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
139

140
def faces_to_maps(faces, *fields):
141
    return [tfields.TensorFields(faces, *fields, dtype=int, dim=3)]
142

143

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


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

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

        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
180
        >>> assert m.triangles().equal(tfields.Triangles3D([[ 1.,  0.,  0.],
181
182
        ...                                               [ 0.,  1.,  0.],
        ...                                               [ 0.,  0.,  0.]]))
183
184

        a list of scalars is assigned to each face
185
186
187
        >>> 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
188
189

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

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

    """
213
    def __new__(cls, tensors, *fields, **kwargs):
214
215
        if not issubclass(type(tensors), Mesh3D):
            kwargs['dim'] = 3
216
217
218
        faces = kwargs.pop('faces', None)
        faceScalars = kwargs.pop('faceScalars', [])
        maps = kwargs.pop('maps', None)
219
        if maps is not None and faces is not None:
220
221
222
            raise ValueError("Conflicting options maps and faces")
        if maps is not None:
            kwargs['maps'] = maps
223
        if len(faceScalars) > 0:
224
225
226
227
228
229
230
231
232
233
234
235
236
            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
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
    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):
            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.
Daniel Boeckenhoff's avatar
dunno    
Daniel Boeckenhoff committed
261
262
            import matplotlib.colors as colors
            import matplotlib.pyplot as plt
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
263
264
265
266
267
268
            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

Daniel Boeckenhoff's avatar
dunno    
Daniel Boeckenhoff committed
269
270
271
        if len(kwargs) != 0:
            raise ValueError("Unused arguments.")

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
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
376
377
378
379
380
381
382
        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

383
384
    @classmethod
    def plane(cls, *base_vectors, **kwargs):
385
386
387
388
389
390
391
        """
        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
392
        vertices = tfields.Tensors.grid(*base_vectors, **kwargs)
393
394
395

        base_vectors = tfields.grid.ensure_complex(*base_vectors)
        base_vectors = tfields.grid.to_base_vectors(*base_vectors)
396
397
398
399
400
401
402
        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
403
404
        if fix_coord is None:
            raise ValueError("Describe a plane with one variable fiexed")
405
406
407
408
409
410
411
412
413
414
415
416
417
418

        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
419
        inst = cls.__new__(cls, vertices, faces=faces)
420
421
422
423
        return inst

    @classmethod
    def grid(cls, *base_vectors, **kwargs):
424
425
        """
        Construct 'cuboid' along base_vectors
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
426
427
428
429
430
431
432
433
434
435
436
437
        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
Daniel Boeckenhoff's avatar
dunno    
Daniel Boeckenhoff committed
438
439
440
441
            >>> oktaeder = tfields.Mesh3D.grid((1, 1, 1),
            ...                                (-np.pi, np.pi, 5),
            ...                                (-np.pi / 2, np.pi / 2, 3),
            ...                                coord_sys='spherical')
442
            >>> oktaeder.transform('cartesian')
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
443
444
445
446
447
448
449
450
451
452
453

            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')
454
            >>> cylinder.transform('cartesian')
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
455

456
        """
457
458
459
        if not len(base_vectors) == 3:
            raise AttributeError("3 base_vectors vectors required")

460
461
462
        base_vectors = tfields.grid.ensure_complex(*base_vectors)
        base_vectors = tfields.grid.to_base_vectors(*base_vectors)

463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
        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
481
                planes.append(cls.plane(*basePart, **kwargs))
482
483
484
        inst = cls.merged(*planes, **kwargs)
        return inst

485
486
487
488
489
490
491
492
493
494
495
496
497
    @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
498
499
    def faceScalars(self, scalars):
        self.maps[0].fields = scalars_to_fields(scalars)
500

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
501
502
    @cached_property()
    def _triangles(self):
503
504
505
506
507
508
        """
        with the decorator, this should be handled like an attribute though it is a function

        """
        if self.faces.size == 0:
            return tfields.Triangles3D([])
509
510
511
512
        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)
513

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
514
    def triangles(self):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
515
516
517
518
519
520
521
522
        """
        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
523
524
525
526
527
528
529
        return self._triangles

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

    @cached_property()
    def _planes(self):
530
531
        if self.faces.size == 0:
            return tfields.Planes3D([])
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
532
533
534
535
        return tfields.Planes3D(self.centroids(), self.triangles().norms())

    def planes(self):
        return self._planes
536

537
    def nfaces(self):
538
539
        return self.faces.shape[0]

540
    def in_faces(self, points, delta, assign_multiple=False):
541
542
        """
        Check whether points lie within triangles with Barycentric Technique
543
        see Triangles3D.in_triangles
544
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
545
        masks = self.triangles().in_triangles(points, delta,
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
546
                                              assign_multiple=assign_multiple)
547
548
        return masks

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
549
    def removeFaces(self, face_delete_mask):
550
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
551
        Remove faces where face_delete_mask is True
552
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
553
554
555
        face_delete_mask = np.array(face_delete_mask, dtype=bool)
        self.faces = self.faces[~face_delete_mask]
        self.faceScalars = self.faceScalars[~face_delete_mask]
556

557
    def template(self, sub_mesh):
558
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
559
        'Manual' way to build a template that can be used with self.cut
560
        Returns:
561
562
            Mesh3D: template (see cut), can be used as template to retrieve
                sub_mesh from self instance
563
        Examples:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
564
565
            >>> mp = tfields.TensorFields([[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
            ...                           [1, 2, 3, 4])
566
            >>> 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
567
            ...                     maps=[mp])
568
            >>> from sympy.abc import y
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
569
570
571
            >>> m_cut = m.cut(y < 1.5, at_intersection='split')
            >>> template = m.template(m_cut)
            >>> assert m_cut.equal(m.cut(template))
572

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
573
574
        TODO:
            fields template not yet implemented
575
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
576
577
        face_indices = np.arange(self.maps[0].shape[0])
        cents = tfields.Tensors(sub_mesh.centroids())
578
        mask = self.in_faces(cents, delta=None)
579
        inst = sub_mesh.copy()
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
580
581
582
583
584
585
586
587
588
589
590
        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)
                        ]
591
592
        return inst

593
    def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False):
594
        """
595
        Partition the mesh with the cuts given and return the template
596
        """
597
598
599
600
601
        eps = 0.000000001
        # direct return if self is empty
        if len(self) == 0:
            return self.copy(), self.copy()

602
603
604
605
606
607
608
609
610
611
612
        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))))

613
        # mask for points that do not fulfill the cut expression
614
        mask = inst.evalf(expression)
615
        # remove the points
616

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
617
        if not any(~mask):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
618
619
620
621
            # no vertex is valid
            inst = inst[mask]
        elif all(~mask):
            # all vertices are valid
622
            inst = inst[mask]
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
        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()
650
651
        elif at_intersection == 'split' or at_intersection == 'splitRough':
            '''
652
            add vertices and faces that are at the border of the cuts
653
            '''
654
            expression_parts = tfields.lib.symbolics.split_expression(expression)
655
            if len(expression_parts) > 1:
656
                new_mesh = inst.copy()
657
658
659
660
661
662
663
664
                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.
                    """
665
666
                    faceIntersMask = np.full((inst.faces.shape[0]), False, dtype=bool)
                    for i, face in enumerate(inst.faces):
667
668
                        vertices_rejected = [-mask[f] for f in face]
                        face_on_edge = any(vertices_rejected) and not all(vertices_rejected)
669
                        if face_on_edge:
670
                            faceIntersMask[i] = True
671
                    new_mesh.removeFaces(-faceIntersMask)
672

673
                for exprPart in expression_parts:
674
675
676
                    inst, _ = inst._cut_sympy(exprPart,
                                              at_intersection='split',
                                              _in_recursion=True)
677
            elif len(expression_parts) == 1:
678
                # TODO maps[0] -> smthng like inst.get_map(dim=3)
679
680
681
                points = [sympy.symbols('x0, y0, z0'),
                          sympy.symbols('x1, y1, z1'),
                          sympy.symbols('x2, y2, z2')]
682
                plane_sympy = tfields.lib.symbolics.to_plane(expression)
683
684
685
                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}
686

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
687
                norm_vectors = inst.triangles().norms()
688
689
                new_points = np.empty((0, 3))
                new_faces = np.empty((0, 3))
690
                new_fields = [tfields.Tensors(np.empty((0,) + field.shape[1:]),
691
                                              coord_sys=field.coord_sys)
692
693
                              for field in inst.fields]
                new_map_fields = [[] for field in inst.maps[0].fields]
694
                new_norm_vectors = []
695
                newScalarMap = []
696
                n_new = 0
697

698
699
700
701
702
703
704
                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]):
705
                    """
706
                    vertices_rejected is a mask for each face that is True, where
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
707
                    a point is on the rejected side of the plane
708
                    """
709
                    vertices_rejected = [~mask[f] for f in face]
710
711
712
713
714
                    if any(vertices_rejected):
                        # delete face
                        face_delete_indices.add(i)
                    if any(vertices_rejected) and not all(vertices_rejected):
                        # face on edge
715
                        nTrue = vertices_rejected.count(True)
716
                        lonely_bool = True if nTrue == 1 else False
717

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

                face_map = tfields.TensorFields(faces, *faces_fields,
753
                                                dtype=int,
754
                                                coord_sys=inst.maps[0].coord_sys)
755
756
757
                inst = tfields.Mesh3D(vertices,
                                      *fields,
                                      maps=[face_map] + inst.maps[1:],
758
                                      coord_sys=inst.coord_sys)
759
760
761
762
763
                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]
764
            else:
765
                raise ValueError("Sympy expression is not splitable.")
766
            inst = inst.cleaned()
767
        elif at_intersection == 'remove':
768
            inst = inst[mask]
769
        else:
770
771
            raise AttributeError("No at_intersection method called {at_intersection} "
                                 "implemented".format(**locals()))
772
773
774
775
776
777
778
779
780
781
782
783
784

        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)
785
        return inst, template
786
787

    def _cut_template(self, template):
788
789
790
791
792
793
794
795
        """
        Args:
            template (tfields.Mesh3D)

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

796
            Build mesh
797
798
799
            >>> 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],
800
801
            ...                    [0.0, 0.1, 0.2, 0.3, 0.4],
            ...                    [0.0, -0.1, -0.2, -0.3, -0.4],
802
803
            ...                    maps=[mmap])

804
            Build template
805
806
807
            >>> tmap = tfields.TensorFields([[0, 3, 4], [0, 1, 2]],
            ...                             [1, 0])
            >>> t = tfields.Mesh3D([[0]*3, [-1]*3, [-2]*3, [-3]*3, [-4]*3],
808
            ...                    [1, 0, 3, 2, 4],
809
810
            ...                    maps=[tmap])

811
            Use template as instruction to make a fast cut
812
813
814
815
816
817
818
819
820
            >>> 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]])
                                   
        """
821
822
        # Possible Extension (small todo): check: len(field(s)) == len(self/maps)

823
        # Redirect fields
824
        fields = []
825
        if template.fields:
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
            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)
842

843
        # Redirect maps and their fields
844
845
        maps = []
        for mp, template_mp in zip(self.maps, template.maps):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
846
            mp_fields = []
847
848
849
850
            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
851
                    mp_fields.append(field[template_mp.fields[0].astype(int)])
852
            new_mp = tfields.TensorFields(tfields.Tensors(template_mp),
853
854
855
                                          *mp_fields)
            maps.append(new_mp)

856
857
        inst = tfields.Mesh3D(tfields.Tensors(template),
                              *fields,
858
                              maps=maps)
859
860
        return inst

861
    def cut(self, expression, coord_sys=None, at_intersection=None,
862
            return_template=False):
863
864
865
        """
        cut method for Mesh3D.
        Args:
866
867
868
869
870
871
872
873
874
875
            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.
876
            coord_sys (coordinate system to cut in):
877
            at_intersection (str): instruction on what to do, when a cut will intersect a triangle.
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
878
879
880
                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.
881
882
            return_template (bool): If True: return the template
                            to redo the same cut fast
883
884
        Examples:
            define the cut
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
885
            >>> import numpy as np
886
            >>> import tfields
887
            >>> from sympy.abc import x,y,z
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
888
            >>> cut_expr = x > 1.5
889

890
891
            >>> m = tfields.Mesh3D.grid((0, 3, 4),
            ...                         (0, 3, 4),
892
            ...                         (0, 0, 1))
893
894
895
896
897
898
            >>> 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
899
            >>> mNew = m.cut(cut_expr)
900
            >>> len(mNew)
901
            8
902
            >>> mNew.nfaces()
903
904
905
906
            6
            >>> float(mNew[:, 0].min())
            2.0

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
907
908
909
910
911
912
913
914
915
916
            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')
917
            >>> float(m_split[:, 0].min())
918
            1.5
919
            >>> len(m_split)
920
            15
921
            >>> m_split.nfaces()
922
923
            15

924
925
            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
926
            >>> m_split_2, template = m.cut(cut_expr, at_intersection='split',
927
            ...                                    return_template=True)
928
929
930
931
932
933
934
935
936
937
            >>> 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)
938

939
            This seems irrelevant at first but consider, the map field or the
940
941
942
943
944
945
946
947
            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:
948
949
            >>> 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')
950
951
952

        Returns:
            copy of cut mesh
953
            * optional: template
954
955

        """
956
        with self.tmp_transform(coord_sys or self.coord_sys):
957
            if isinstance(expression, Mesh3D):
958
959
                template = expression
                obj = self._cut_template(template)
960
961
962
963
964
965
            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
966

967
968
969
970
971
972
973
974
975
976
977
978
979
    def disjoint_parts(self, return_template=False):
        mp_description = self.disjoint_map(0)
        parts = self.parts(mp_description)
        if not return_template:
            return parts
        else:
            templates = []
            for i, part in enumerate(parts):
                template = part.copy()
                template.maps[0].fields[0] = tfields.Tensors(mp_description[1][i])
                templates.append(template)
            return parts, templates

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
980
    def plot(self, **kwargs):  # pragma: no cover
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
981
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
982
        Forwarding to plotTools.plot_mesh
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
        """
        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
1005
        return tfields.plotting.plot_mesh(self, self.faces, **kwargs)
1006

1007

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1008
if __name__ == '__main__':  # pragma: no cover
1009
1010
    import doctest

1011
1012
    doctest.run_docstring_examples(Mesh3D.cut, globals())
    # doctest.testmod()