mesh3D.py 46 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
#!/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
11
import rna
12
import tfields
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
13
14

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


20
21
22
23
24
25
def _dist_from_plane(point, plane):
    return plane['normal'].dot(point) + plane['d']


def _segment_plane_intersection(p0, p1, plane):
    """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
26
27
28
29
30
    Get the intersection between the ray spanned by p0 and p1 with the plane.
    Args:
        p0: array of length 3
        p1: array of length 3
        plane: 
31
32
    Returns:
        points, direction
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
33
            points (list of arrays of length 3): 3d points
34
35
36
    """
    distance0 = _dist_from_plane(p0, plane)
    distance1 = _dist_from_plane(p1, plane)
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
37
38
    p0_on_plane = abs(distance0) < np.finfo(float).eps
    p1_on_plane = abs(distance1) < np.finfo(float).eps
39
    points = []
40
    direction = 0
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
41
    if p0_on_plane:
42
        points.append(p0)
43

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
44
    if p1_on_plane:
45
        points.append(p1)
46
    # remove duplicate points
47
48
    if len(points) > 1:
        points = np.unique(points, axis=0)
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
49
    if p0_on_plane and p1_on_plane:
50
        return points, direction
51
52

    if distance0 * distance1 > np.finfo(float).eps:
53
        return points, direction
54
55
56

    direction = np.sign(distance0)
    if abs(distance0) < np.finfo(float).eps:
57
        return points, direction
58
    elif abs(distance1) < np.finfo(float).eps:
59
        return points, direction
60
61
62
    if abs(distance0 - distance1) > np.finfo(float).eps:
        t = distance0 / (distance0 - distance1)
    else:
63
        return points, direction
64

65
    points.append(p0 + t * (p1 - p0))
66
    # remove duplicate points
67
68
69
    if len(points) > 1:
        points = np.unique(points, axis=0)
    return points, direction
70
71
72
73
74
75
76


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:
77
78
79
80
81
82
83
84
85
        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
86
    """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
87
88
    n_true = vertices_rejected.count(True)
    lonely_bool = True if n_true == 1 else False
89
90
91
92
93
    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)

94
95
96
    single_index = index
    couple_indices = [j for j in range(3)
                      if not vertices_rejected[j] == lonely_bool]
97
98
99
100
101
102

    # 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
103
        return new_points
104
105
    if len(s1) == 2:
        # both points on plane
106
        return new_points
107
108
    if len(s2) == 2:
        # both points on plane
109
        return new_points
110
    if lonely_bool:
111
        # two new triangles
112
        if len(s0) == 1 and len(s1) == 1:
113
114
            new_points = [[couple_indices[0], s0[0], couple_indices[1]],
                          [couple_indices[1], complex(1), s1[0]]]
115
        elif len(s1) == 1 and len(s2) == 1:
116
117
            new_points = [[couple_indices[0], couple_indices[1], s1[0]],
                          [couple_indices[0], complex(2), s2[0]]]
118
        elif len(s0) == 1 and len(s2) == 1:
119
120
            new_points = [[couple_indices[0], couple_indices[1], s0[0]],
                          [couple_indices[1], s2[0], complex(2)]]
121
    else:
122
        # one new triangle
123
        if len(s0) == 1 and len(s1) == 1:
124
            new_points = [[single_index, s1[0], s0[0]]]
125
        elif len(s1) == 1 and len(s2) == 1:
126
            new_points = [[single_index, s2[0], s1[0]]]
127
        elif len(s0) == 1 and len(s2) == 1:
128
129
            new_points = [[single_index, s0[0], s2[0]]]
    return new_points
130
131


132
133
134
135
136
137
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
138

139
140
141
def fields_to_scalars(fields):
    return np.array(fields)

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
142

143
def faces_to_maps(faces, *fields):
144
    return [tfields.TensorFields(faces, *fields, dtype=int, dim=3)]
145

146

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


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

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

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

        a list of scalars is assigned to each face
188
189
190
        >>> 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
191
192

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

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

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

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
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
        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))

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
315
316
317
318
319
320
321
322
323
324
325
    @classmethod
    def _load_stl(cls, path):
        """
        Factory method
        Given a path to a stl file, construct the object
        """
        import stl.mesh
        mesh = stl.mesh.Mesh.from_file(path)
        mesh = cls(vertices, faces=faces)
        return mesh

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
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
383
384
385
386
387
388
389
390
391
392
393
394
395
396
    @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

397
398
    @classmethod
    def plane(cls, *base_vectors, **kwargs):
399
400
401
402
403
404
405
        """
        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
406
        vertices = tfields.Tensors.grid(*base_vectors, **kwargs)
407
408
409

        base_vectors = tfields.grid.ensure_complex(*base_vectors)
        base_vectors = tfields.grid.to_base_vectors(*base_vectors)
410
411
412
413
414
415
416
        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
417
418
        if fix_coord is None:
            raise ValueError("Describe a plane with one variable fiexed")
419
420
421
422
423
424
425
426
427
428
429
430
431
432

        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
433
        inst = cls.__new__(cls, vertices, faces=faces)
434
435
436
437
        return inst

    @classmethod
    def grid(cls, *base_vectors, **kwargs):
438
439
        """
        Construct 'cuboid' along base_vectors
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
440
441
442
443
444
445
446
447
448
449
450
451
        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
452
453
454
455
            >>> oktaeder = tfields.Mesh3D.grid((1, 1, 1),
            ...                                (-np.pi, np.pi, 5),
            ...                                (-np.pi / 2, np.pi / 2, 3),
            ...                                coord_sys='spherical')
456
            >>> oktaeder.transform('cartesian')
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
457
458
459
460
461
462
463
464
465
466
467

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

470
        """
471
472
473
        if not len(base_vectors) == 3:
            raise AttributeError("3 base_vectors vectors required")

474
475
476
        base_vectors = tfields.grid.ensure_complex(*base_vectors)
        base_vectors = tfields.grid.to_base_vectors(*base_vectors)

477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
        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
495
                planes.append(cls.plane(*basePart, **kwargs))
496
497
498
        inst = cls.merged(*planes, **kwargs)
        return inst

499
500
501
502
503
504
505
506
507
508
509
510
511
    @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
512
513
    def faceScalars(self, scalars):
        self.maps[0].fields = scalars_to_fields(scalars)
514

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
515
516
    @cached_property()
    def _triangles(self):
517
518
519
520
521
522
        """
        with the decorator, this should be handled like an attribute though it is a function

        """
        if self.faces.size == 0:
            return tfields.Triangles3D([])
523
524
525
526
        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)
527

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
528
    def triangles(self):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
529
530
531
532
533
534
535
536
        """
        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
537
538
539
540
541
542
543
        return self._triangles

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

    @cached_property()
    def _planes(self):
544
545
        if self.faces.size == 0:
            return tfields.Planes3D([])
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
546
547
548
549
        return tfields.Planes3D(self.centroids(), self.triangles().norms())

    def planes(self):
        return self._planes
550

551
    def nfaces(self):
552
553
        return self.faces.shape[0]

554
    def in_faces(self, points, delta, assign_multiple=False):
555
556
        """
        Check whether points lie within triangles with Barycentric Technique
557
        see Triangles3D.in_triangles
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
558
559
560
561
        If multiple requests are done on huge meshes,
        this can be hugely optimized by requesting the property
        self.tree or setting it to self.tree = <saved tree> before
        calling in_faces
562
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
563
564
565
566
567
568
569
570
571
572
        key = 'mesh_tree'
        if hasattr(self, '_cache') and key in self._cache:
            log = logging.getLogger()
            log.info(
                "Using cached decision tree to speed up point - face mapping.")
            masks = self.tree.in_faces(
                points, delta, assign_multiple=assign_multiple)
        else:
            masks = self.triangles().in_triangles(
                points, delta, assign_multiple=assign_multiple)
573
574
        return masks

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
    @property
    def tree(self):
        """
        Cached property to retrieve a bounding_box Searcher. This searcher can
        hugely optimize 'in_faces' searches

        Examples:
            >>> mesh = tfields.Mesh3D.grid((0, 1, 3), (1, 2, 3), (2, 3, 3))
            >>> _ = mesh.tree
            >>> assert hasattr(mesh, '_cache')
            >>> assert 'mesh_tree' in mesh._cache
            >>> mask = mesh.in_faces(tfields.Points3D([[0.2, 1.2, 2.0]]),
            ...                      0.00001)
            >>> assert mask.sum() == 1  # one point in one triangle
        """
        if not hasattr(self, '_cache'):
            self._cache = {}
        key = 'mesh_tree'
        if key in self._cache:
            tree = self._cache[key]
        else:
            tree = tfields.bounding_box.Searcher(self)
            self._cache[key] = tree
        return tree

    @tree.setter
    def tree(self, tree):
        if not hasattr(self, '_cache'):
            self._cache = {}
        key = 'mesh_tree'
        self._cache[key] = tree

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
607
    def removeFaces(self, face_delete_mask):
608
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
609
        Remove faces where face_delete_mask is True
610
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
611
612
613
        face_delete_mask = np.array(face_delete_mask, dtype=bool)
        self.faces = self.faces[~face_delete_mask]
        self.faceScalars = self.faceScalars[~face_delete_mask]
614

615
    def template(self, sub_mesh):
616
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
617
        'Manual' way to build a template that can be used with self.cut
618
        Returns:
619
620
            Mesh3D: template (see cut), can be used as template to retrieve
                sub_mesh from self instance
621
        Examples:
622
623
            >>> import tfields
            >>> from sympy.abc import y
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
624
625
            >>> mp = tfields.TensorFields([[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
            ...                           [1, 2, 3, 4])
626
            >>> 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
627
628
629
630
            ...                     maps=[mp])
            >>> m_cut = m.cut(y < 1.5, at_intersection='split')
            >>> template = m.template(m_cut)
            >>> assert m_cut.equal(m.cut(template))
631

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
632
633
        TODO:
            fields template not yet implemented
634
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
635
636
        face_indices = np.arange(self.maps[0].shape[0])
        cents = tfields.Tensors(sub_mesh.centroids())
637
        mask = self.in_faces(cents, delta=None)
638
        inst = sub_mesh.copy()
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
639
640
641
642
643
644
645
646
647
648
649
        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)
                        ]
650
651
        return inst

652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
    def project(self, tensor_field,
                delta=None, merge_functions=None):
        """
        project the points of the tensor_field to a copy of the mesh
        and set the face values accord to the field to the maps field.
        If no field is present in tensor_field, the number of points in a mesh
        is counted.

        Args:
            tensor_field (Tensors | TensorFields)
            delta (float | None): forwarded to Mesh3D.in_faces
            merge_functions (callable): if multiple Tensors lie in the same face,
                they are mapped with the merge_function to one value

        Examples:
            >>> import tfields
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
668
            >>> import numpy as np
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
            >>> mp = tfields.TensorFields([[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
            ...                           [1, 2, 3, 4])
            >>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,2,0], [1,2,0]],
            ...                     maps=[mp])

            Projecting points onto the mesh gives the count
            >>> points = tfields.Tensors([[0.5, 0.2, 0.0], [0.5, 0.02, 0.0], [0.5, 0.8, 0.0]])
            >>> m_points = m.project(points)
            >>> assert m_points.maps[0].fields[0].equal([2, 1, 0, 0])

            TensorFields with arbitrary size are projected,
            combinging the fields automatically
            >>> fields = [tfields.Tensors([1,3,42]),
            ...           tfields.Tensors([[0,1,2], [2,3,4], [3,4,5]]),
            ...           tfields.Tensors([[[0, 0]] * 2,
            ...                            [[2, 2]] * 2,
            ...                            [[3, 3]] * 2])]
            >>> tf = tfields.TensorFields(points, *fields)
            >>> m_tf = m.project(tf)
            >>> assert m_tf.maps[0].fields[0].equal([2, 42, np.nan, np.nan],
            ...                                     equal_nan=True)
            >>> assert m_tf.maps[0].fields[1].equal([[1, 2, 3],
            ...                                      [3, 4, 5],
            ...                                      [np.nan] * 3,
            ...                                      [np.nan] * 3],
            ...                                     equal_nan=True)
            >>> assert m_tf.maps[0].fields[2].equal([[[1, 1]] * 2,
            ...                                      [[3, 3]] * 2,
            ...                                      [[np.nan, np.nan]] * 2,
            ...                                      [[np.nan, np.nan]] * 2],
            ...                                     equal_nan=True)

        """
        if not issubclass(type(tensor_field), tfields.Tensors):
            tensor_field = tfields.TensorFields(tensor_field)
        mask = self.in_faces(tensor_field, delta=None)
        inst = self.copy()
        
        n_faces = len(self.maps[0])
        if not hasattr(tensor_field, 'fields') or len(tensor_field.fields) == 0:
            fields = [np.full(len(tensor_field), 1)]
            empty_map_fields = [tfields.Tensors(np.full(n_faces, 0))]
            if merge_functions is None:
                merge_functions = [np.sum]
        else:
            fields = tensor_field.fields
            empty_map_fields = []
            for field in fields:
                cls = type(field)
                kwargs = {key: getattr(field, key) for key in cls.__slots__}
                shape = (n_faces,) + field.shape[1:]
                empty_map_fields.append(cls(np.full(shape, np.nan),
                                            **kwargs))
            if merge_functions is None:
                merge_functions = [lambda x: np.mean(x, axis=0)] * len(fields)

        map_fields = []
        for field, map_field, merge_function in \
                zip(fields, empty_map_fields, merge_functions):
            for f in range(len(self.maps[0])):
                res = field[mask[:, f]]
                if len(res) == 0:
                    continue
                elif len(res) == 1:
                    map_field[f] = res
                else:
                    map_field[f] = merge_function(res)
            map_fields.append(map_field)
        inst.maps[0].fields = map_fields
        return inst

740
    def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False):
741
        """
742
        Partition the mesh with the cuts given and return the template
743
        """
744
745
746
747
748
        eps = 0.000000001
        # direct return if self is empty
        if len(self) == 0:
            return self.copy(), self.copy()

749
750
751
752
753
754
755
756
757
758
759
        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))))

760
        # mask for points that do not fulfill the cut expression
761
        mask = inst.evalf(expression)
762
        # remove the points
763

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
764
        if not any(~mask):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
765
766
767
768
            # no vertex is valid
            inst = inst[mask]
        elif all(~mask):
            # all vertices are valid
769
            inst = inst[mask]
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
        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()
797
798
        elif at_intersection == 'split' or at_intersection == 'splitRough':
            '''
799
            add vertices and faces that are at the border of the cuts
800
            '''
801
            expression_parts = tfields.lib.symbolics.split_expression(expression)
802
            if len(expression_parts) > 1:
803
                new_mesh = inst.copy()
804
805
806
807
808
809
810
811
                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.
                    """
812
813
                    faceIntersMask = np.full((inst.faces.shape[0]), False, dtype=bool)
                    for i, face in enumerate(inst.faces):
814
815
                        vertices_rejected = [-mask[f] for f in face]
                        face_on_edge = any(vertices_rejected) and not all(vertices_rejected)
816
                        if face_on_edge:
817
                            faceIntersMask[i] = True
818
                    new_mesh.removeFaces(-faceIntersMask)
819

820
                for exprPart in expression_parts:
821
822
823
                    inst, _ = inst._cut_sympy(exprPart,
                                              at_intersection='split',
                                              _in_recursion=True)
824
            elif len(expression_parts) == 1:
825
                # TODO maps[0] -> smthng like inst.get_map(dim=3)
826
827
828
                points = [sympy.symbols('x0, y0, z0'),
                          sympy.symbols('x1, y1, z1'),
                          sympy.symbols('x2, y2, z2')]
829
                plane_sympy = tfields.lib.symbolics.to_plane(expression)
830
831
832
                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}
833

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
834
                norm_vectors = inst.triangles().norms()
835
836
                new_points = np.empty((0, 3))
                new_faces = np.empty((0, 3))
837
                new_fields = [tfields.Tensors(np.empty((0,) + field.shape[1:]),
838
                                              coord_sys=field.coord_sys)
839
840
                              for field in inst.fields]
                new_map_fields = [[] for field in inst.maps[0].fields]
841
                new_norm_vectors = []
842
                newScalarMap = []
843
                n_new = 0
844

845
846
847
848
849
850
851
                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]):
852
                    """
853
                    vertices_rejected is a mask for each face that is True, where
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
854
                    a point is on the rejected side of the plane
855
                    """
856
                    vertices_rejected = [~mask[f] for f in face]
857
858
859
860
861
                    if any(vertices_rejected):
                        # delete face
                        face_delete_indices.add(i)
                    if any(vertices_rejected) and not all(vertices_rejected):
                        # face on edge
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
862
863
                        n_true = vertices_rejected.count(True)
                        lonely_bool = True if n_true == 1 else False
864

Daniel Boeckenhoff's avatar
dunno    
Daniel Boeckenhoff committed
865
                        triangle_points = [vertices[f] for f in face]
866
                        """
867
                        Add the intersection points and faces
868
                        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
869
                        intersection = _intersect(triangle_points, plane, vertices_rejected)
870
                        last_idx = len(vertices) - 1
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
871
                        for tri_list in intersection:
872
                            new_face = []
873
874
875
                            for item in tri_list:
                                if isinstance(item, int):
                                    # reference to old vertex
876
                                    new_face.append(face[item])
877
878
879
                                elif isinstance(item, complex):
                                    # reference to new vertex that has been
                                    # concatenated already
880
                                    new_face.append(last_idx + int(item.imag))
881
882
                                else:
                                    # new vertex
883
                                    new_face.append(len(vertices))
884
                                    vertices = np.append(vertices,
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
885
                                                         [[float(x) for x in item]],
886
887
888
889
890
                                                         axis=0)
                                    fields = [np.append(field,
                                                        np.full((1,) + field.shape[1:], np.nan),
                                                        axis=0)
                                              for field in fields]
891
                            faces = np.append(faces, [new_face], axis=0)
892
                            faces_fields = [np.append(field,
893
                                                      [field[i]],
894
895
896
897
898
                                                      axis=0)
                                            for field in faces_fields]
                            faces_fields[-1][-1] = i

                face_map = tfields.TensorFields(faces, *faces_fields,
899
                                                dtype=int,
900
                                                coord_sys=inst.maps[0].coord_sys)
901
902
903
                inst = tfields.Mesh3D(vertices,
                                      *fields,
                                      maps=[face_map] + inst.maps[1:],
904
                                      coord_sys=inst.coord_sys)
905
906
907
908
909
                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]
910
            else:
911
                raise ValueError("Sympy expression is not splitable.")
912
            inst = inst.cleaned()
913
        elif at_intersection == 'remove':
914
            inst = inst[mask]
915
        else:
916
917
            raise AttributeError("No at_intersection method called {at_intersection} "
                                 "implemented".format(**locals()))
918
919
920
921
922
923
924
925
926
927
928
929
930

        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)
931
        return inst, template
932
933

    def _cut_template(self, template):
934
935
936
937
938
939
940
941
        """
        Args:
            template (tfields.Mesh3D)

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

942
            Build mesh
943
944
945
            >>> 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],
946
947
            ...                    [0.0, 0.1, 0.2, 0.3, 0.4],
            ...                    [0.0, -0.1, -0.2, -0.3, -0.4],
948
949
            ...                    maps=[mmap])

950
            Build template
951
952
953
            >>> tmap = tfields.TensorFields([[0, 3, 4], [0, 1, 2]],
            ...                             [1, 0])
            >>> t = tfields.Mesh3D([[0]*3, [-1]*3, [-2]*3, [-3]*3, [-4]*3],
954
            ...                    [1, 0, 3, 2, 4],
955
956
            ...                    maps=[tmap])

957
            Use template as instruction to make a fast cut
958
959
960
961
962
963
964
965
966
            >>> 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]])
                                   
        """
967
968
        # Possible Extension (small todo): check: len(field(s)) == len(self/maps)

969
        # Redirect fields
970
        fields = []
971
        if template.fields:
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
            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)
988

989
        # Redirect maps and their fields
990
991
        maps = []
        for mp, template_mp in zip(self.maps, template.maps):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
992
            mp_fields = []
993
994
995
996
            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
997
                    mp_fields.append(field[template_mp.fields[0].astype(int)])
998
            new_mp = tfields.TensorFields(tfields.Tensors(template_mp),
999
1000
1001
                                          *mp_fields)
            maps.append(new_mp)

1002
1003
        inst = tfields.Mesh3D(tfields.Tensors(template),
                              *fields,
1004
                              maps=maps)
1005
1006
        return inst

1007
    def cut(self, expression, coord_sys=None, at_intersection=None,
1008
            return_template=False):
1009
1010
1011
        """
        cut method for Mesh3D.
        Args:
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
            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.
1022
            coord_sys (coordinate system to cut in):
1023
            at_intersection (str): instruction on what to do, when a cut will intersect a triangle.
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1024
1025
1026
                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.
1027
1028
            return_template (bool): If True: return the template
                            to redo the same cut fast
1029
1030
        Examples:
            define the cut
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1031
            >>> import numpy as np
1032
            >>> import tfields
1033
            >>> from sympy.abc import x,y,z
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1034
            >>> cut_expr = x > 1.5
1035

1036
1037
            >>> m = tfields.Mesh3D.grid((0, 3, 4),
            ...                         (0, 3, 4),
1038
            ...                         (0, 0, 1))
1039
1040
1041
1042
1043
1044
            >>> 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
1045
            >>> mNew = m.cut(cut_expr)
1046
            >>> len(mNew)
1047
            8
1048
            >>> mNew.nfaces()
1049
1050
1051
1052
            6
            >>> float(mNew[:, 0].min())
            2.0

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
            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')
1063
            >>> float(m_split[:, 0].min())
1064
            1.5
1065
            >>> len(m_split)
1066
            15
1067
            >>> m_split.nfaces()
1068
1069
            15

1070
1071
            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
1072
            >>> m_split_2, template = m.cut(cut_expr, at_intersection='split',
1073
            ...                                    return_template=True)
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
            >>> 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)
1084

1085
            This seems irrelevant at first but consider, the map field or the
1086
1087
1088
1089
1090
1091
1092
1093
            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:
1094
1095
            >>> 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')
1096
1097
1098

        Returns:
            copy of cut mesh
1099
            * optional: template
1100
1101

        """
1102
        with self.tmp_transform(coord_sys or self.coord_sys):
1103
            if isinstance(expression, Mesh3D):
1104
1105
                template = expression
                obj = self._cut_template(template)
1106
1107
1108
1109
1110
1111
            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
1112

1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
    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
1126
    def plot(self, **kwargs):  # pragma: no cover
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1127
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1128
        Forwarding to plotTools.plot_mesh
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1129
1130
1131
1132
1133
1134
1135
1136
        """
        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
1137
        if 'axes' in kwargs:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1138
            dim_defined = True
1139
1140
        if 'z_index' in kwargs:
            if kwargs['z_index'] is not None:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
                kwargs['dim'] = 3
            else:
                kwargs['dim'] = 2
            dim_defined = True
        if 'dim' in kwargs:
            dim_defined = True

        if not dim_defined:
            kwargs['dim'] = 2

1151
        return rna.plotting.plot_mesh(self, self.faces, **kwargs)
1152

1153

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1154
if __name__ == '__main__':  # pragma: no cover
1155
1156
    import doctest

1157
    doctest.run_docstring_examples(Mesh3D.project, globals())
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1158
    doctest.run_docstring_examples(Mesh3D.tree, globals())
1159
    # doctest.testmod()