mesh3D.py 47.7 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
    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
dboe's avatar
dboe committed
30
        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)
dboe's avatar
dboe committed
210
        >>> m1 = tfields.Mesh3D.load(outFile.name, allow_pickle=True)
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

            Cube with edge length of 2 units
            >>> cube = tfields.Mesh3D.grid((-1, 1, 2),
            ...                            (-1, 1, 2),
            ...                            (-5, -3, 2))

dboe's avatar
dboe committed
463
            Cylinder
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
464
465
466
467
            >>> 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
    def project(self, tensor_field,
653
654
655
                delta=None, merge_functions=None,
                point_face_assignment=None,
                return_point_face_assignment=False):
656
657
658
659
660
661
662
663
664
665
666
        """
        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
667
668
669
670
671
            point_face_assignment (np.array, dtype=int): array assigning each
                point to a face. Each entry position corresponds to a point of the
                tensor, each entry value is the index of the assigned face
            return_point_face_assignment (bool): if true, return the
                point_face_assignment
672
673
674

        Examples:
            >>> import tfields
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
675
            >>> import numpy as np
676
677
678
679
680
681
            >>> 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
682
683
684
685
686
            >>> points = tfields.Tensors([[0.5, 0.2, 0.0],
            ...                           [0.5, 0.02, 0.0],
            ...                           [0.5, 0.8, 0.0],
            ...                           [0.5, 0.8, 0.1]])  # not contained
            >>> m_points = m.project(points, delta=0.01)
687
688
689
690
            >>> assert m_points.maps[0].fields[0].equal([2, 1, 0, 0])

            TensorFields with arbitrary size are projected,
            combinging the fields automatically
691
692
            >>> fields = [tfields.Tensors([1,3,42, -1]),
            ...           tfields.Tensors([[0,1,2], [2,3,4], [3,4,5], [-1] * 3]),
693
694
            ...           tfields.Tensors([[[0, 0]] * 2,
            ...                            [[2, 2]] * 2,
695
696
            ...                            [[3, 3]] * 2,
            ...                            [[9, 9]] * 2])]
697
            >>> tf = tfields.TensorFields(points, *fields)
698
            >>> m_tf = m.project(tf, delta=0.01)
699
700
701
702
703
704
705
706
707
708
709
710
711
            >>> 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)

712
713
714
715
716
717
718
            Returning the calculated point_face_assignment can speed up multiple
            results
            >>> m_tf, point_face_assignment = m.project(tf, delta=0.01,
            ...     return_point_face_assignment=True)
            >>> m_tf_fast = m.project(tf, delta=0.01, point_face_assignment=point_face_assignment)
            >>> assert m_tf.equal(m_tf_fast, equal_nan=True)

719
720
721
722
        """
        if not issubclass(type(tensor_field), tfields.Tensors):
            tensor_field = tfields.TensorFields(tensor_field)
        inst = self.copy()
723
724

        # setup empty map fields and collect fields
725
        n_faces = len(self.maps[0])
726
        point_indices = np.arange(len(tensor_field))
727
        if not hasattr(tensor_field, 'fields') or len(tensor_field.fields) == 0:
728
729
730
731
            # if not fields is existing use int type fields and empty_map_fields
            # in order to generate a sum
            fields = [np.full(len(tensor_field), 1, dtype=int)]
            empty_map_fields = [tfields.Tensors(np.full(n_faces, 0, dtype=int))]
732
733
734
735
736
737
738
739
740
741
742
743
744
            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)
745

746
747
748
749
750
751
752
        # build point_face_assignment if not given.
        if point_face_assignment is not None:
            if len(point_face_assignment) != len(tensor_field):
                raise ValueErrror("Template needs same lenght as tensor_field")
            point_face_assignment = point_face_assignment
        else:
            mask = self.in_faces(tensor_field, delta=delta)
753

754
755
756
757
758
759
760
761
762
            face_indices = np.arange(n_faces)
            point_face_assignment = [face_indices[mask[p_index, :]]
                                      for p_index in point_indices]
            point_face_assignment = np.array([fa if len(fa) != 0 else [-1]
                                               for fa in point_face_assignment])
            point_face_assignment = point_face_assignment.reshape(-1)
        point_face_assignment_set = set(point_face_assignment)

        # merge the fields according to point_face_assignment
763
764
765
        map_fields = []
        for field, map_field, merge_function in \
                zip(fields, empty_map_fields, merge_functions):
766
            for i, f_index in enumerate(point_face_assignment_set):
767
                # if i % 1000 == 0:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
768
                #     print(i, len(point_face_assignment_set))
769
770
                if f_index == -1:
                    # point could not be mapped
771
                    continue
772
773
774
775
                point_in_face_indices = point_indices[point_face_assignment == f_index]
                res = field[point_in_face_indices]
                if len(res) == 1:
                    map_field[f_index] = res
776
                else:
777
                    map_field[f_index] = merge_function(res)
778
779
            map_fields.append(map_field)
        inst.maps[0].fields = map_fields
780
781
        if return_point_face_assignment:
            return inst, point_face_assignment
782
783
        return inst

784
    def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False):
785
        """
786
        Partition the mesh with the cuts given and return the template
787
        """
788
789
790
791
792
        eps = 0.000000001
        # direct return if self is empty
        if len(self) == 0:
            return self.copy(), self.copy()

793
794
795
796
797
798
799
800
801
802
803
        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))))

804
        # mask for points that do not fulfill the cut expression
805
        mask = inst.evalf(expression)
806
        # remove the points
807

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
808
        if not any(~mask):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
809
810
811
812
            # no vertex is valid
            inst = inst[mask]
        elif all(~mask):
            # all vertices are valid
813
            inst = inst[mask]
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
        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()
841
842
        elif at_intersection == 'split' or at_intersection == 'splitRough':
            '''
843
            add vertices and faces that are at the border of the cuts
844
            '''
845
            expression_parts = tfields.lib.symbolics.split_expression(expression)
846
            if len(expression_parts) > 1:
847
                new_mesh = inst.copy()
848
849
850
851
852
853
854
855
                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.
                    """
856
857
                    faceIntersMask = np.full((inst.faces.shape[0]), False, dtype=bool)
                    for i, face in enumerate(inst.faces):
858
859
                        vertices_rejected = [-mask[f] for f in face]
                        face_on_edge = any(vertices_rejected) and not all(vertices_rejected)
860
                        if face_on_edge:
861
                            faceIntersMask[i] = True
862
                    new_mesh.removeFaces(-faceIntersMask)
863

864
                for exprPart in expression_parts:
865
866
867
                    inst, _ = inst._cut_sympy(exprPart,
                                              at_intersection='split',
                                              _in_recursion=True)
868
            elif len(expression_parts) == 1:
869
                # TODO maps[0] -> smthng like inst.get_map(dim=3)
870
871
872
                points = [sympy.symbols('x0, y0, z0'),
                          sympy.symbols('x1, y1, z1'),
                          sympy.symbols('x2, y2, z2')]
873
                plane_sympy = tfields.lib.symbolics.to_plane(expression)
874
875
876
                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}
877

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
878
                norm_vectors = inst.triangles().norms()
879
880
                new_points = np.empty((0, 3))
                new_faces = np.empty((0, 3))
881
                new_fields = [tfields.Tensors(np.empty((0,) + field.shape[1:]),
882
                                              coord_sys=field.coord_sys)
883
884
                              for field in inst.fields]
                new_map_fields = [[] for field in inst.maps[0].fields]
885
                new_norm_vectors = []
886
                newScalarMap = []
887
                n_new = 0
888

889
890
891
892
893
894
895
                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]):
896
                    """
897
                    vertices_rejected is a mask for each face that is True, where
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
898
                    a point is on the rejected side of the plane
899
                    """
900
                    vertices_rejected = [~mask[f] for f in face]
901
902
903
904
905
                    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
906
907
                        n_true = vertices_rejected.count(True)
                        lonely_bool = True if n_true == 1 else False
908

Daniel Boeckenhoff's avatar
dunno    
Daniel Boeckenhoff committed
909
                        triangle_points = [vertices[f] for f in face]
910
                        """
911
                        Add the intersection points and faces
912
                        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
913
                        intersection = _intersect(triangle_points, plane, vertices_rejected)
914
                        last_idx = len(vertices) - 1
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
915
                        for tri_list in intersection:
916
                            new_face = []
917
918
919
                            for item in tri_list:
                                if isinstance(item, int):
                                    # reference to old vertex
920
                                    new_face.append(face[item])
921
922
923
                                elif isinstance(item, complex):
                                    # reference to new vertex that has been
                                    # concatenated already
924
                                    new_face.append(last_idx + int(item.imag))
925
926
                                else:
                                    # new vertex
927
                                    new_face.append(len(vertices))
928
                                    vertices = np.append(vertices,
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
929
                                                         [[float(x) for x in item]],
930
931
932
933
934
                                                         axis=0)
                                    fields = [np.append(field,
                                                        np.full((1,) + field.shape[1:], np.nan),
                                                        axis=0)
                                              for field in fields]
935
                            faces = np.append(faces, [new_face], axis=0)
936
                            faces_fields = [np.append(field,
937
                                                      [field[i]],
938
939
940
941
942
                                                      axis=0)
                                            for field in faces_fields]
                            faces_fields[-1][-1] = i

                face_map = tfields.TensorFields(faces, *faces_fields,
943
                                                dtype=int,
944
                                                coord_sys=inst.maps[0].coord_sys)
945
946
947
                inst = tfields.Mesh3D(vertices,
                                      *fields,
                                      maps=[face_map] + inst.maps[1:],
948
                                      coord_sys=inst.coord_sys)
949
950
951
952
953
                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]
954
            else:
955
                raise ValueError("Sympy expression is not splitable.")
956
            inst = inst.cleaned()
957
        elif at_intersection == 'remove':
958
            inst = inst[mask]
959
        else:
960
961
            raise AttributeError("No at_intersection method called {at_intersection} "
                                 "implemented".format(**locals()))
962
963
964
965
966
967
968
969
970
971
972
973
974

        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)
975
        return inst, template
976
977

    def _cut_template(self, template):
978
979
980
981
982
983
984
985
        """
        Args:
            template (tfields.Mesh3D)

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

986
            Build mesh
987
988
989
            >>> 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],
990
991
            ...                    [0.0, 0.1, 0.2, 0.3, 0.4],
            ...                    [0.0, -0.1, -0.2, -0.3, -0.4],
992
993
            ...                    maps=[mmap])

994
            Build template
995
996
997
            >>> tmap = tfields.TensorFields([[0, 3, 4], [0, 1, 2]],
            ...                             [1, 0])
            >>> t = tfields.Mesh3D([[0]*3, [-1]*3, [-2]*3, [-3]*3, [-4]*3],
998
            ...                    [1, 0, 3, 2, 4],
999
1000
            ...                    maps=[tmap])

1001
            Use template as instruction to make a fast cut
1002
1003
1004
1005
1006
1007
1008
            >>> 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]])
dboe's avatar
dboe committed
1009

1010
        """
1011
1012
        # Possible Extension (small todo): check: len(field(s)) == len(self/maps)

1013
        # Redirect fields
1014
        fields = []
1015
        if template.fields:
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
            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)
1032

1033
        # Redirect maps and their fields
1034
1035
        maps = []
        for mp, template_mp in zip(self.maps, template.maps):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1036
            mp_fields = []
1037
1038
1039
1040
            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
1041
                    mp_fields.append(field[template_mp.fields[0].astype(int)])
1042
            new_mp = tfields.TensorFields(tfields.Tensors(template_mp),
1043
1044
1045
                                          *mp_fields)
            maps.append(new_mp)

1046
1047
        inst = tfields.Mesh3D(tfields.Tensors(template),
                              *fields,
1048
                              maps=maps)
1049
1050
        return inst

1051
    def cut(self, *args, **kwargs):
1052
1053
1054
        """
        cut method for Mesh3D.
        Args:
1055
1056
1057
1058
1059
1060
1061
            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'.
dboe's avatar
dboe committed
1062
                    We use the vertices and maps of the Mesh as the skelleton of
1063
1064
                    the returned mesh. The fields are mapped according to
                    indices in the template.maps[i].fields.
1065
            coord_sys (coordinate system to cut in):
1066
            at_intersection (str): instruction on what to do, when a cut will intersect a triangle.
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1067
1068
1069
                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.
1070
1071
            return_template (bool): If True: return the template
                            to redo the same cut fast
1072
1073
        Examples:
            define the cut
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1074
            >>> import numpy as np
1075
            >>> import tfields
1076
            >>> from sympy.abc import x,y,z
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1077
            >>> cut_expr = x > 1.5
1078

1079
1080
            >>> m = tfields.Mesh3D.grid((0, 3, 4),
            ...                         (0, 3, 4),
1081
            ...                         (0, 0, 1))
1082
1083
1084
1085
1086
1087
            >>> 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
1088
            >>> mNew = m.cut(cut_expr)
1089
            >>> len(mNew)
1090
            8
1091
            >>> mNew.nfaces()
1092
1093
1094
1095
            6
            >>> float(mNew[:, 0].min())
            2.0

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
            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')
1106
            >>> float(m_split[:, 0].min())
1107
            1.5
1108
            >>> len(m_split)
1109
            15
1110
            >>> m_split.nfaces()
1111
1112
            15

1113
1114
            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
1115
            >>> m_split_2, template = m.cut(cut_expr, at_intersection='split',
1116
            ...                                    return_template=True)
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
            >>> 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)
1127

1128
            This seems irrelevant at first but consider, the map field or the
1129
1130
1131
1132
1133
1134
1135
1136
            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:
1137
1138
            >>> 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')
1139
1140
1141

        Returns:
            copy of cut mesh
1142
            * optional: template
1143
1144

        """
1145
        return super().cut(*args, **kwargs)
1146

1147
    def disjoint_parts(self, return_template=False):
dboe's avatar
dboe committed
1148
1149
1150
1151
        """
        Returns:
            disjoint_parts(List(cls)), templates(List(cls))
        """
1152
1153
1154
1155
1156
1157
1158
1159
        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()
dboe's avatar
dboe committed
1160
1161
1162
                template.maps[0].fields = [
                    tfields.Tensors(mp_description[1][i])
                ]
1163
1164
1165
                templates.append(template)
            return parts, templates