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

part of tfields library
"""
import numpy as np
import sympy
import tfields
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
12
from tfields.lib.decorators import cached_property
13
14


15
16
17
18
19
20
21
22
23
24
25
26
27
def _dist_from_plane(point, plane):
    return plane['normal'].dot(point) + plane['d']


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

    if p1OnPlane:
34
        points.append(p1)
35
    # remove duplicate points
36
37
    if len(points) > 1:
        points = np.unique(points, axis=0)
38
    if p0OnPlane and p1OnPlane:
39
        return points, direction
40
41

    if distance0 * distance1 > np.finfo(float).eps:
42
        return points, direction
43
44
45

    direction = np.sign(distance0)
    if abs(distance0) < np.finfo(float).eps:
46
        return points, direction
47
    elif abs(distance1) < np.finfo(float).eps:
48
        return points, direction
49
50
51
    if abs(distance0 - distance1) > np.finfo(float).eps:
        t = distance0 / (distance0 - distance1)
    else:
52
        return points, direction
53

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


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:
66
67
68
69
70
71
72
73
74
        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
75
76
77
78
79
80
81
82
    """
    nTrue = vertices_rejected.count(True)
    lonely_bool = True if nTrue == 1 else False
    index = vertices_rejected.index(lonely_bool)
    s0, d0 = _segment_plane_intersection(triangle[0], triangle[1], plane)
    s1, d1 = _segment_plane_intersection(triangle[1], triangle[2], plane)
    s2, d2 = _segment_plane_intersection(triangle[2], triangle[0], plane)

83
84
85
    single_index = index
    couple_indices = [j for j in range(3)
                      if not vertices_rejected[j] == lonely_bool]
86
87
88
89
90
91

    # 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
92
        return new_points
93
94
    if len(s1) == 2:
        # both points on plane
95
        return new_points
96
97
    if len(s2) == 2:
        # both points on plane
98
        return new_points
99
    if lonely_bool:
100
        # two new triangles
101
        if len(s0) == 1 and len(s1) == 1:
102
103
            new_points = [[couple_indices[0], s0[0], couple_indices[1]],
                          [couple_indices[1], complex(1), s1[0]]]
104
        elif len(s1) == 1 and len(s2) == 1:
105
106
            new_points = [[couple_indices[0], couple_indices[1], s1[0]],
                          [couple_indices[0], complex(2), s2[0]]]
107
        elif len(s0) == 1 and len(s2) == 1:
108
109
            new_points = [[couple_indices[0], couple_indices[1], s0[0]],
                          [couple_indices[1], s2[0], complex(2)]]
110
111

    else:
112
        # one new triangle
113
        if len(s0) == 1 and len(s1) == 1:
114
            new_points = [[single_index, s1[0], s0[0]]]
115
        elif len(s1) == 1 and len(s2) == 1:
116
            new_points = [[single_index, s2[0], s1[0]]]
117
        elif len(s0) == 1 and len(s2) == 1:
118
119
            new_points = [[single_index, s0[0], s2[0]]]
    return new_points
120
121


122
123
124
125
126
127
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
128

129
130
131
def fields_to_scalars(fields):
    return np.array(fields)

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
132

133
def faces_to_maps(faces, *fields):
134
    return [tfields.TensorFields(faces, *fields, dtype=int, dim=3)]
135

136

137
def maps_to_faces(maps):
138
139
140
141
142
    if len(maps) == 0:
        return np.array([])
    elif len(maps) > 1:
        raise NotImplementedError("Multiple maps")
    return np.array(maps[0])
143
144


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

        conversion to points only
162
163
164
165
166
        >>> tfields.Points3D(m).equal([[1, 2, 3],
        ...                            [3, 3, 3],
        ...                            [0, 0, 0],
        ...                            [5, 6, 7]])
        True
167
168
169
170
171
172

        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
173
        >>> assert m.triangles().equal(tfields.Triangles3D([[ 1.,  0.,  0.],
174
175
        ...                                               [ 0.,  1.,  0.],
        ...                                               [ 0.,  0.,  0.]]))
176
177

        a list of scalars is assigned to each face
178
179
180
        >>> 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
181
182

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

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

    """
206
    def __new__(cls, tensors, *fields, **kwargs):
207
208
        if not issubclass(type(tensors), Mesh3D):
            kwargs['dim'] = 3
209
210
211
        faces = kwargs.pop('faces', None)
        faceScalars = kwargs.pop('faceScalars', [])
        maps = kwargs.pop('maps', None)
212
        if maps is not None and faces is not None:
213
214
215
            raise ValueError("Conflicting options maps and faces")
        if maps is not None:
            kwargs['maps'] = maps
216
        if len(faceScalars) > 0:
217
218
219
220
221
222
223
224
225
226
227
228
229
            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

230
231
    @classmethod
    def plane(cls, *base_vectors, **kwargs):
232
233
234
235
236
237
238
        """
        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__
        """
239
        vertices = tfields.Tensors.grid(*base_vectors)
240
241
242

        base_vectors = tfields.grid.ensure_complex(*base_vectors)
        base_vectors = tfields.grid.to_base_vectors(*base_vectors)
243
244
245
246
247
248
249
        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
250
251
        if fix_coord is None:
            raise ValueError("Describe a plane with one variable fiexed")
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270

        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])
        inst = cls.__new__(cls, vertices, faces=faces, **kwargs)
        return inst

    @classmethod
    def grid(cls, *base_vectors, **kwargs):
271
272
273
        """
        Construct 'cuboid' along base_vectors
        """
274
275
276
        if not len(base_vectors) == 3:
            raise AttributeError("3 base_vectors vectors required")

277
278
279
        base_vectors = tfields.grid.ensure_complex(*base_vectors)
        base_vectors = tfields.grid.to_base_vectors(*base_vectors)

280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
        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)

                planes.append(cls.plane(*basePart))
        inst = cls.merged(*planes, **kwargs)
        return inst

303
304
305
306
307
308
309
310
311
312
313
314
315
    @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
316
317
    def faceScalars(self, scalars):
        self.maps[0].fields = scalars_to_fields(scalars)
318

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
319
320
    @cached_property()
    def _triangles(self):
321
322
323
324
325
326
        """
        with the decorator, this should be handled like an attribute though it is a function

        """
        if self.faces.size == 0:
            return tfields.Triangles3D([])
327
328
329
330
        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)
331

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
332
333
334
335
336
337
338
339
    def triangles(self):
        return self._triangles

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

    @cached_property()
    def _planes(self):
340
341
        if self.faces.size == 0:
            return tfields.Planes3D([])
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
342
343
344
345
        return tfields.Planes3D(self.centroids(), self.triangles().norms())

    def planes(self):
        return self._planes
346

347
    def nfaces(self):
348
349
        return self.faces.shape[0]

350
    def in_faces(self, points, delta, assign_multiple=False):
351
352
353
        """
        Check whether points lie within triangles with Barycentric Technique
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
354
        masks = self.triangles().in_triangles(points, delta,
355
                                            assign_multiple=assign_multiple)
356
357
358
359
360
361
362
363
364
365
        return masks

    def removeFaces(self, faceDeleteMask):
        """
        Remove faces where faceDeleteMask is True
        """
        faceDeleteMask = np.array(faceDeleteMask, dtype=bool)
        self.faces = self.faces[~faceDeleteMask]
        self.faceScalars = self.faceScalars[~faceDeleteMask]

366
    def template(self, sub_mesh, delta=1e-9):
367
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
368
        'Manual' way to build a template that can be used with self.cut
369
        Returns:
370
371
            Mesh3D: template (see cut), can be used as template to retrieve
                sub_mesh from self instance
372
        Examples:
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
373
374
            >>> mp = tfields.TensorFields([[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
            ...                           [1, 2, 3, 4])
375
            >>> 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
376
            ...                     maps=[mp])
377
            >>> from sympy.abc import y
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
378
379
380
            >>> m_cut = m.cut(y < 1.5, at_intersection='split')
            >>> template = m.template(m_cut)
            >>> assert m_cut.equal(m.cut(template))
381

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
382
383
        TODO:
            fields template not yet implemented
384
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
385
386
        face_indices = np.arange(self.maps[0].shape[0])
        cents = tfields.Tensors(sub_mesh.centroids())
387
        scalars = []
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
388
389
390
391
        mask = self.in_faces(cents, delta=delta)
        scalars = []
        for face_mask in mask:
            scalars.append(face_indices[face_mask][0])
392
        inst = sub_mesh.copy()
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
393
        inst.maps[0].fields = [tfields.Tensors(scalars)]
394
395
        return inst

396
    def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False):
397
        """
398
        Partition the mesh with the cuts given and return the template
399
400

        """
401
402
403
404
405
        eps = 0.000000001
        # direct return if self is empty
        if len(self) == 0:
            return self.copy(), self.copy()

406
407
408
409
410
411
412
413
414
415
416
        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))))

417
        # mask for points that do not fulfill the cut expression
418
        mask = inst.evalf(expression)
419
        # remove the points
420
421

        if not any(~mask) or all(~mask):
422
            inst = inst[mask]
423
424
        elif at_intersection == 'split' or at_intersection == 'splitRough':
            '''
425
            add vertices and faces that are at the border of the cuts
426
            '''
427
            expression_parts = tfields.lib.symbolics.split_expression(expression)
428
            if len(expression_parts) > 1:
429
                new_mesh = inst.copy()
430
431
432
433
434
435
436
437
                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.
                    """
438
439
                    faceIntersMask = np.full((inst.faces.shape[0]), False, dtype=bool)
                    for i, face in enumerate(inst.faces):
440
441
                        vertices_rejected = [-mask[f] for f in face]
                        face_on_edge = any(vertices_rejected) and not all(vertices_rejected)
442
                        if face_on_edge:
443
                            faceIntersMask[i] = True
444
                    new_mesh.removeFaces(-faceIntersMask)
445

446
                for exprPart in expression_parts:
447
448
449
                    inst, _ = inst._cut_sympy(exprPart,
                                              at_intersection='split',
                                              _in_recursion=True)
450
            elif len(expression_parts) == 1:
451
                # TODO maps[0] -> smthng like inst.get_map(dim=3)
452
453
454
                points = [sympy.symbols('x0, y0, z0'),
                          sympy.symbols('x1, y1, z1'),
                          sympy.symbols('x2, y2, z2')]
455
                plane_sympy = tfields.lib.symbolics.to_plane(expression)
456
457
458
                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}
459

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
460
                norm_vectors = inst.triangles().norms()
461
462
                new_points = np.empty((0, 3))
                new_faces = np.empty((0, 3))
463
                new_fields = [tfields.Tensors(np.empty((0,) + field.shape[1:]),
464
                                              coord_sys=field.coord_sys)
465
466
                              for field in inst.fields]
                new_map_fields = [[] for field in inst.maps[0].fields]
467
                new_norm_vectors = []
468
                newScalarMap = []
469
                n_new = 0
470

471
472
473
474
475
476
477
                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]):
478
                    """
479
                    vertices_rejected is a mask for each face that is True, where
480
481
                    a Point is on the rejected side of the plane
                    """
482
                    vertices_rejected = [~mask[f] for f in face]
483
484
485
486
487
                    if any(vertices_rejected):
                        # delete face
                        face_delete_indices.add(i)
                    if any(vertices_rejected) and not all(vertices_rejected):
                        # face on edge
488
                        nTrue = vertices_rejected.count(True)
489
                        lonely_bool = True if nTrue == 1 else False
490

491
                        triangle_points = [inst[f] for f in face]
492
                        """
493
                        Add the intersection points and faces
494
                        """
495
496
497
                        newP = _intersect(triangle_points, plane, vertices_rejected)
                        last_idx = len(vertices) - 1
                        for tri_list in newP:
498
                            new_face = []
499
500
501
                            for item in tri_list:
                                if isinstance(item, int):
                                    # reference to old vertex
502
                                    new_face.append(face[item])
503
504
505
                                elif isinstance(item, complex):
                                    # reference to new vertex that has been
                                    # concatenated already
506
                                    new_face.append(last_idx + int(item.imag))
507
508
                                else:
                                    # new vertex
509
                                    new_face.append(len(vertices))
510
511
512
513
514
515
516
                                    vertices = np.append(vertices,
                                                         [map(float, item)],
                                                         axis=0)
                                    fields = [np.append(field,
                                                        np.full((1,) + field.shape[1:], np.nan),
                                                        axis=0)
                                              for field in fields]
517
                            faces = np.append(faces, [new_face], axis=0)
518
                            faces_fields = [np.append(field,
519
                                                      [field[i]],
520
521
522
523
524
                                                      axis=0)
                                            for field in faces_fields]
                            faces_fields[-1][-1] = i

                face_map = tfields.TensorFields(faces, *faces_fields,
525
                                                dtype=int,
526
                                                coord_sys=inst.maps[0].coord_sys)
527
528
529
                inst = tfields.Mesh3D(vertices,
                                      *fields,
                                      maps=[face_map] + inst.maps[1:],
530
                                      coord_sys=inst.coord_sys)
531
532
533
534
535
                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]
536
            else:
537
                raise ValueError("Sympy expression is not splitable.")
538
            inst = inst.cleaned()
539
        elif at_intersection == 'remove':
540
            inst = inst[mask]
541
        else:
542
543
            raise AttributeError("No at_intersection method called {at_intersection} "
                                 "implemented".format(**locals()))
544
545
546
547
548
549
550
551
552
553
554
555
556

        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)
557
        return inst, template
558
559

    def _cut_template(self, template):
560
561
562
563
564
565
566
567
        """
        Args:
            template (tfields.Mesh3D)

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

568
            Build mesh
569
570
571
            >>> 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],
572
573
            ...                    [0.0, 0.1, 0.2, 0.3, 0.4],
            ...                    [0.0, -0.1, -0.2, -0.3, -0.4],
574
575
            ...                    maps=[mmap])

576
            Build template
577
578
579
            >>> tmap = tfields.TensorFields([[0, 3, 4], [0, 1, 2]],
            ...                             [1, 0])
            >>> t = tfields.Mesh3D([[0]*3, [-1]*3, [-2]*3, [-3]*3, [-4]*3],
580
            ...                    [1, 0, 3, 2, 4],
581
582
            ...                    maps=[tmap])

583
            Use template as instruction to make a fast cut
584
585
586
587
588
589
590
591
592
            >>> 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]])
                                   
        """
593
        # Redirect fields
594
        fields = []
595
        if template.fields:
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
            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)
612
613

        # Redirect maps fields
614
615
        maps = []
        for mp, template_mp in zip(self.maps, template.maps):
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
616
617
618
619
            mp_fields = []
            if len(template_mp.fields) > 0:
                for field in mp.fields:
                    mp_fields.append(field[template_mp.fields[0].astype(int)])
620
            new_mp = tfields.TensorFields(tfields.Tensors(template_mp),
621
622
623
                                          *mp_fields)
            maps.append(new_mp)

624
625
        inst = tfields.Mesh3D(tfields.Tensors(template),
                              *fields,
626
                              maps=maps)
627
628
        return inst

629
    def cut(self, expression, coord_sys=None, at_intersection=None,
630
            return_template=False):
631
632
633
        """
        cut method for Mesh3D.
        Args:
634
635
636
637
638
639
640
641
642
643
            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.
644
            coord_sys (coordinate system to cut in):
645
            at_intersection (str): instruction on what to do, when a cut will intersect a triangle.
646
647
                Options:    "remove" (Default)
                            "split" - Create new triangles that make up the old one.
648
649
            return_template (bool): If True: return the template
                            to redo the same cut fast
650
651
        Examples:
            define the cut
652
            >>> import tfields
653
654
655
            >>> from sympy.abc import x,y,z
            >>> cutExpr = x > 1.5

656
657
            >>> m = tfields.Mesh3D.grid((0, 3, 4),
            ...                         (0, 3, 4),
658
            ...                         (0, 0, 1))
659
660
661
662
663
664
            >>> 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]))))
665
            >>> mNew = m.cut(cutExpr)
666
            >>> len(mNew)
667
            8
668
            >>> mNew.nfaces()
669
670
671
672
673
            6
            >>> float(mNew[:, 0].min())
            2.0

            Cutting with the split option will create new triangles on the edge:
674
675
            >>> m_split = m.cut(cutExpr, at_intersection='split')
            >>> float(m_split[:, 0].min())
676
            1.5
677
            >>> len(m_split)
678
            15
679
            >>> m_split.nfaces()
680
681
            15

682
683
684
685
            Cut with 'return_template=True' will return the exact same mesh but
            additionally an instruction to conduct the exact same cut fast (template)
            >>> m_split_2, template = m.cut(cutExpr, at_intersection='split',
            ...                                    return_template=True)
686
687
688
689
690
691
692
693
694
695
            >>> 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)
696
697
698
699
700
701
702
703
704
705
706

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


            The cut expression may be a sympy.BooleanFunction:
707
708
            >>> 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')
709
710
711

        Returns:
            copy of cut mesh
712
            * optional: template
713
714

        """
715
        with self.tmp_transform(coord_sys or self.coord_sys):
716
717
718
719
720
721
722
723
            if isinstance(expression, Mesh3D):
                obj = self._cut_template(expression)
            else:
                at_intersection = at_intersection or "remove"
                obj, template = self._cut_sympy(expression, at_intersection=at_intersection)
        if return_template:
            return obj, template
        return obj
724

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
725
    def plot(self, **kwargs):  # pragma: no cover
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
726
        """
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
727
        Forwarding to plotTools.plot_mesh
Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
        """
        scalars_demanded = any([v in kwargs for v in ['vmin', 'vmax', 'cmap']])
        map_index = kwargs.pop('map_index', None if not scalars_demanded else 0)
        if map_index is not None:
            if not len(self.maps[0]) == 0:
                kwargs['color'] = self.maps[0].fields[map_index]

        dim_defined = False
        if 'axis' in kwargs:
            dim_defined = True
        if 'zAxis' in kwargs:
            if kwargs['zAxis'] is not None:
                kwargs['dim'] = 3
            else:
                kwargs['dim'] = 2
            dim_defined = True
        if 'dim' in kwargs:
            dim_defined = True

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

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

752

Daniel Boeckenhoff's avatar
Daniel Boeckenhoff committed
753
if __name__ == '__main__':  # pragma: no cover
754
755
    import doctest

756
    # doctest.run_docstring_examples(Mesh3D.cut, globals())
757
    doctest.testmod()